Update 13_02

models/cDiazepam_parameters.json

@@ -1,8 +1,8 @@
 {
      "parameters":{
         "commentsInput":"X",
-        "_venousInput":{"value":0.010,"commentUnit":"sameUnitAsVolumes"},
-        "venousInput":{"function":"venousInput"},
+        "venousInput":{"value":0.010,"commentUnit":"sameUnitAsVolumes"},
+        "_venousInput":{"function":"venousInput"},
         "inputA":{"value":0.25e6,"commentUnit":"AU;reported concentrations are in AU/volumeUnit","dist":"normal","cv":0},
         "inputTau":{"value":1,"commentUnit":"sameAsTimeUnit","dist":"normal","cv":0},
         "commentPC":"X",

models/cDiazepam_parameters1.json

@@ -1,8 +1,8 @@
 {
      "parameters":{
         "commentsInput":"X",
-        "venousInput":{"value":0.3e6,"commentUnit":"sameUnitAsVolumes"},
-        "_venousInput":{"function":"venousInput"},
+        "_venousInput":{"value":0.3e6,"commentUnit":"sameUnitAsVolumes"},
+        "venousInput":{"function":"venousInput"},
         "inputA":{"value":0.01,"commentUnit":"sameAsVolumes"},
         "inputTau":{"value":1,"commentUnit":"sameAsTimeUnit"},
         "commentPC":"X",

pythonScripts/runSolver.py

@@ -97,35 +97,30 @@ def getFiles(setup, getAll=False):
       inFiles.append({x:'{}{}.txt'.format(x,l) for x in readable})
    return inFiles
 
-
-def main(parFiles,jobDir,startDir='NONE'):
-   model=cModel.model()
-   setupFile=parFiles[0]
-   modelFile=parFiles[1]
-   parameterFile=parFiles[2]
-   model.parse(modelFile,parameterFile)
-   setup=parseSetup(setupFile)
-   scale=getScale(setup)
-   tmax=float(setup['tmax'])*scale
-   print('Using {}'.format(jobDir))
+def findStartPoint(startDir,setup):
    if startDir!='NONE':
       print('Using solution from {}'.format(startDir))
       #t0,y0,s0,lut,lutSE=getStartPointFromDir(startDir)
       startPoint=getStartPointFromDir(startDir)
       print('t0={}'.format(startPoint['t0']))
-   else:
-      try:
-         print('Using solution from {}'.format(setup['startFromRef']))
-      except KeyError:
-         pass
-      startPoint=getStartPoint(setup)
-      print('t0={}'.format(startPoint['t0']))
+      return startPoint
+   try:
+      print('Using solution from {}'.format(setup['startFromRef']))
+   except KeyError:
+      pass
+   startPoint=getStartPoint(setup)
+   print('t0={}'.format(startPoint['t0']))
+   return startPoint
 
+
+
+def solve(model,setup,startPoint,jobDir):
+   #only uses values of parameters from memory
    #we should reorder S1 to match model.lutSE
-   #S1 is nxm, so we should reorder columns
    t0=startPoint['t0']
    y0=startPoint['y0']
    s0=startPoint['s0']
+   #S1 is nxm, so we should reorder columns
    sIn=numpy.zeros(s0.shape)
    if s0.size>0:
       for x in model.lutSE:
@@ -143,6 +138,8 @@ def main(parFiles,jobDir,startDir='NONE'):
       pass
    setup['t0']=t0 
    strides=getStrides(setup)
+   full_t=[]
+   full_sol=[]
    for step in strides:
 
       #check if solution exists
@@ -154,12 +151,12 @@ def main(parFiles,jobDir,startDir='NONE'):
          #set t0,y0,sIn
          latestFiles={x:'{}{}.txt'.format(x,l) for x in readable}
          data=loadSolutionFromFiles(setup,jobDir,[latestFiles])
-         startPoint=startPointObject(data)
-         t0=startPoint['t0']
-         y0=startPoint['y0']
-         s0=startPoint['s0']
+         _startPoint=startPointObject(data)
+         t0=_startPoint['t0']
+         y0=_startPoint['y0']
+         s0=_startPoint['s0']
          for x in model.lutSE:
-            sIn[:,model.lutSE[x]]=s0[:,startPoint['lutSE'][x]]
+            sIn[:,model.lutSE[x]]=s0[:,_startPoint['lutSE'][x]]
          print('Completed step {}: {}'.format(l,filesPresent))
          continue
       
@@ -210,8 +207,36 @@ def main(parFiles,jobDir,startDir='NONE'):
             write3D(x,writeable[x],model.lut,model.lutSE)
          else:
             numpy.savetxt(x,writeable[x])
+
+      full_t.extend(t)
+      full_sol.extend(sol)
       print('Completed step {}'.format(l))
+   
    print('Completed run')
+   return t,sol,se,qt,sOut,full_t,full_sol
+
+
+
+def main(parFiles,jobDir,startDir='NONE'):
+   print(parFiles)
+   model=cModel.model()
+   setupFile=parFiles[0]
+   modelFile=parFiles[1]
+   parameterFile=parFiles[2]
+   model.parse(modelFile,parameterFile)
+   setup=parseSetup(setupFile)
+   scale=getScale(setup)
+   tmax=float(setup['tmax'])*scale
+   print('Using {}'.format(jobDir))
+
+   # set start point
+   startPoint=findStartPoint(startDir,setup)
+   t,sol,se,qt,sOut,full_t,full_sol=solve(model,setup,startPoint,jobDir)
+
+   #new solution
+   #model.updateParameters(...)
+   #t,sol,se,qt,sOut=solve(model,setup,startPoint,jobDir)
+
    #update setup with new content (particularly t0)
    setupOut=os.path.join(jobDir,'setup.json')
    with open(setupOut,'w+') as f:

pythonScripts/uporablja_run_solver.py

@@ -0,0 +1,369 @@
+import ivp
+import solveMatrix
+import os
+import  cModel
+import sys
+import json
+from contextlib import contextmanager
+import numpy as np
+import time
+import scipy.interpolate
+import shutil
+import importlib
+from scipy.optimize import least_squares
+importlib.reload(solveMatrix)
+
+import runSolver as runSolver
+importlib.reload(runSolver)
+
+def updateSetup(job):
+    setupFileSrc = job['setupFile']
+    setupFile = os.path.join(job['jobDir'], 'setupInput.json')
+    with open(setupFileSrc, 'r') as f:
+        setup = json.load(f)
+    try:
+        setup.update(job['setupUpdates'])
+    except KeyError:
+        pass
+    with open(setupFile, 'w+') as f:
+        f.write(json.dumps(setup))
+    return setupFile
+
+
+def generate_noisy_data(full_sol, noise_std=20.5):
+    full_sol = np.asarray(full_sol)  # Ensure input is a NumPy array
+    random_state = np.random.RandomState(seed=None)
+    
+    # Create highly volatile noise
+    noise = random_state.normal(loc=0, scale=noise_std, size=full_sol.shape) ** 2  # Squared to amplify large deviations
+    noise2 = np.random.laplace(loc=0, scale=noise_std * 10, size=full_sol.shape)  # Increased scale for even more extreme spikes
+    
+    # Introduce highly chaotic scaling factors
+    random_scaling = random_state.uniform(-200, 400, size=full_sol.shape)  # Further extended range for more variation
+    random_exponent = random_state.uniform(4, 12, size=full_sol.shape)  # Increased exponentiation for more instability
+    
+    # Apply extreme, chaotic transformations
+    noisy_data = np.zeros_like(full_sol)
+    for i in range(len(full_sol)):
+        perturbation = random_scaling[i] * (abs(noise[i]) ** random_exponent[i]) - noise2[i] * np.sin(full_sol[i])
+        noisy_data[i] = (1 + np.tanh(perturbation)) * full_sol[i] * np.cosh(noise[i])
+        
+        # Introduce frequent large perturbations to increase instability
+        if random_state.rand() < 0.75:  # 30% chance of a drastic shift
+            noisy_data[i] *= random_state.uniform(2, 5)
+    
+    return noisy_data
+
+
+def calculate_residuals(full_sol, noisy_data):
+      """
+      Calculates residuals based on the solution and noisy data.
+      """
+      residuals = np.ravel((noisy_data - full_sol))
+      #print(residuals)
+      return residuals
+
+
+# Updated cost function
+def cost_function(params,parnames,noisy_data,model,startPoint,setup,jobDir):
+    """
+         Cost function that calculates residuals based on the model.
+         """     
+    #print(params)
+    model.setValues(parnames,params)
+    clean_dir(jobDir)
+    t,sol,se,qt,sOut,full_t,full_sol=runSolver.solve(model,setup,startPoint,jobDir)
+    residuals = calculate_residuals(np.array(full_sol), noisy_data)
+    return residuals
+
+
+def create_name_value_dict(original_dict):
+        return {
+            val['name']: val['value']
+            for key, val in original_dict.items()
+            if isinstance(val, dict) and 'name' in val and 'value' in val
+        }
+    
+    # Create parameter dictionary
+def clean_dir(jobDir):
+      
+    shutil.rmtree(jobDir)
+    os.mkdir(jobDir)
+
+
+
+#PRENOS ZAČETNIH PODATKOV 
+##############################################################################################################################################
+
+fh = os.path.expanduser("~")
+jobDir = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen',"rezultati")
+setup_file = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'setup', 'setupMinute.json')
+model_file = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'models', 'cDiazepam.json')
+parameter_file1 = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'models', 'cDiazepam_parameters.json')
+
+
+# Prvi model s prvim parameter file
+model1 = cModel.model()
+model1.parse(model_file, parameter_file1)
+setup1 = runSolver.parseSetup(setup_file)
+name_value_dict1 = create_name_value_dict(model1.parSetup['parameters'])
+
+
+# Lahko zdaj uporabiš name_value_dict1 in name_value_dict2 ločeno
+
+
+ # Excluded parameters
+excluded_params1 = {'venousInput',"heartPC","kidneyPC","remainderPC","skinPC","splanchnicPC","stomachPC","testesPC","one","venousVolume","heartVolume","brainVolume","adiposeVolume","remainderFlow"}
+
+    # Create initial guesses without excluded parameters
+initial_guesses = np.array(
+        [val for par, val in name_value_dict1.items() if par not in excluded_params1],
+        dtype=np.float64
+    )
+
+    # Create list of parameter names without excluded
+param_names = [par for par in name_value_dict1.keys() if par not in excluded_params1]
+initial_pars = np.array([x for x in initial_guesses])
+#initial_pars[3]*=1.1
+clean_dir(jobDir)
+#print("kolicnik")
+#print(initial_guesses/initial_pars)
+#print(param_names)
+
+
+
+# TO NAREDI REŠITVE Z ZAČETNIMI PARAMETRI
+start_point= runSolver.findStartPoint("NONE",setup1)
+model1.setValues(param_names,initial_pars)
+t,sol,se,qt,sOut,full_t,full_sol = runSolver.solve(model1,setup1,start_point,jobDir)
+clean_dir(jobDir)
+
+
+# Step 4: Create a new dictionary for the optimized parameters
+# Map the optimized parameters back to the correct order in the original parameter dictionary
+optimized_params_dict1 = {par: initial_guesses[i] for i, par in enumerate(param_names)}
+
+# Step 5: Include the excluded parameters with their original values
+for par in excluded_params1:
+    optimized_params_dict1[par] = name_value_dict1[par]
+
+# Step 6: Ensure that everything stays in order in the final dictionary
+# This will maintain the order of parameters as they appear in the original `name_value_dict2`
+final_dict1 = {par: optimized_params_dict1.get(par, name_value_dict1[par]) for par in name_value_dict1}
+
+
+#################################################################################################################333
+
+
+
+
+
+
+# Z DRUGIMI PARAMETRI
+##########################################################################################################################
+
+# Drugi model z drugim parameter file
+parameter_file2 = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'models', 'cDiazepam_parameters1.json')
+model2 = cModel.model()
+model2.parse(model_file, parameter_file2)
+setup2 = runSolver.parseSetup(setup_file)
+name_value_dict2 = create_name_value_dict(model2.parSetup['parameters'])
+
+ # Excluded parameters
+excluded_params = {'venousInput',"heartPC","kidneyPC","remainderPC","skinPC","splanchnicPC","stomachPC","testesPC","one","venousVolume","heartVolume","brainVolume","adiposeVolume","remainderFlow"}
+
+
+# Step 1: Filter out excluded parameters for initial guesses
+initial_guesses2 = np.array(
+    [val for par, val in name_value_dict2.items() if par not in excluded_params],
+    dtype=np.float64
+)
+
+param_names2 = [par for par in name_value_dict2.keys() if par not in excluded_params]
+
+
+# Step 4: Create a new dictionary for the optimized parameters
+# Map the optimized parameters back to the correct order in the original parameter dictionary
+optimized_params_dict22 = {par: initial_guesses2[i] for i, par in enumerate(param_names2)}
+
+# Step 5: Include the excluded parameters with their original values
+for par in excluded_params1:
+    optimized_params_dict22[par] = name_value_dict2[par]
+
+# Step 6: Ensure that everything stays in order in the final dictionary
+# This will maintain the order of parameters as they appear in the original `name_value_dict2`
+final_dict22 = {par: optimized_params_dict22.get(par, name_value_dict1[par]) for par in name_value_dict2}
+
+# Step 2: Run optimization
+result = least_squares(
+    cost_function, initial_guesses2,
+    args=(param_names2, full_sol, model2, start_point, setup2, jobDir)
+)
+
+# Step 3: Retrieve optimized parameters
+optimized_params2 = result.x
+
+# Step 4: Create a new dictionary for the optimized parameters
+# Map the optimized parameters back to the correct order in the original parameter dictionary
+optimized_params_dict2 = {par: optimized_params2[i] for i, par in enumerate(param_names2)}
+
+# Step 5: Include the excluded parameters with their original values
+for par in excluded_params:
+    optimized_params_dict2[par] = name_value_dict2[par]
+
+# Step 6: Ensure that everything stays in order in the final dictionary
+# This will maintain the order of parameters as they appear in the original `name_value_dict2`
+final_optimized_dict = {par: optimized_params_dict2.get(par, name_value_dict2[par]) for par in name_value_dict2}
+
+# Now, final_optimized_dict contains the optimized values in the same order as in `name_value_dict2`
+print(name_value_dict2)
+print(final_optimized_dict)
+
+vrednosti3 = np.array(
+        [val for par, val in final_optimized_dict.items()],
+        dtype=np.float64
+    )
+
+param_names3 = [par for par in final_optimized_dict.keys()]
+#model2.setValues(param_names3,vrednosti3)
+#################################################################################################################################
+# PREVRBA ČE SO 3 različne rešitve
+clean_dir(jobDir)
+t,sol,se,qt,sOut,full_t,full_sol1 = runSolver.solve(model1,setup1,start_point,jobDir)
+clean_dir(jobDir)
+t,sol,se,qt,sOut,full_t,full_sol2 = runSolver.solve(model2,setup2,start_point,jobDir)
+
+
+model2.setValues(param_names3,vrednosti3)
+clean_dir(jobDir)
+t,sol,se,qt,sOut,full_t,full_sol3 = runSolver.solve(model2,setup2,start_point,jobDir)
+
+
+print(np.array(full_sol2)/np.array(full_sol1))
+print(np.array(full_sol3)/np.array(full_sol1))
+print(np.array(full_sol3)/np.array(full_sol2))
+
+
+print(initial_guesses2/initial_guesses)
+print(optimized_params2/initial_guesses2)
+print(optimized_params2/initial_guesses)
+##################################################################################################################################
+
+
+
+
+# SHRANIMO REŠITVE ZA PARAMETRE
+import os
+import json
+import shutil
+
+def create_or_clean_dir(directory):
+    """Create the directory if it doesn't exist, otherwise clean it."""
+    if os.path.exists(directory):
+        shutil.rmtree(directory)  # Remove the existing directory and its contents
+    os.makedirs(directory)  # Create a fresh directory
+
+def save_parameters(file_path, data):
+    """Save dictionary data to a JSON file."""
+    with open(file_path, 'w') as f:
+        json.dump(data, f, indent=4)
+    return file_path  # Return the path for immediate use
+
+# Define unique subdirectories
+jobDir1 = os.path.join(jobDir, "run1")
+jobDir2 = os.path.join(jobDir, "run2")
+jobDir3 = os.path.join(jobDir, "run3")
+
+# Ensure directories exist and are clean
+create_or_clean_dir(jobDir1)
+create_or_clean_dir(jobDir2)
+create_or_clean_dir(jobDir3)
+
+# Solve and save results in respective directories
+t, sol, se, qt, sOut, full_t, full_sol1 = runSolver.solve(model1, setup1, start_point, jobDir1)
+save_parameters(os.path.join(jobDir1, "parameters.json"), final_dict1)
+
+t, sol, se, qt, sOut, full_t, full_sol2 = runSolver.solve(model2, setup2, start_point, jobDir2)
+save_parameters(os.path.join(jobDir2, "parameters.json"), final_dict22)
+
+model2.setValues(param_names3, vrednosti3)
+t, sol, se, qt, sOut, full_t, full_sol3 = runSolver.solve(model2, setup2, start_point, jobDir3)
+
+# Save parameters and immediately store the path for further use
+param_file3 = save_parameters(os.path.join(jobDir3, "parameters.json"), final_optimized_dict)
+
+# Now you can directly use param_file3 in further processing
+print(f"Parameters for run3 saved at: {param_file3}")
+
+
+########################################################################################################################
+
+
+
+
+
+
+
+
+
+
+
+#ZGRADIMO 3 MODEL (da bo lažje pol risat)
+##############################################################################################################################################
+# param_file3 already contains the full path to "parameters.json"
+directory = os.path.dirname(param_file3)
+model3 = cModel.model()
+parameter_file3 = os.path.join(directory,"parameters.json")  # Use it as a path
+model3.parse(model_file, parameter_file3)
+# Load the model using the saved parameter file
+
+#model3.parse(model_file, parameter_file3)  # Use param_file3 here
+# Parse the setup
+setup3 = runSolver.parseSetup(setup_file)
+# Create the name-value dictionary
+name_value_dict3 = create_name_value_dict(model3.parSetup['parameters'])
+
+########################################################################################################3333
+
+
+
+
+
+runSolver.main([setup_file,model3,parameter_file3,jobDir])
+
+
+
+
+
+
+
+
+#RISANJE
+########################################################################################################################
+
+job_configs = {
+    'cDiazepam_IVP': {
+        'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'rezultati','run1')
+    },
+    'cDiazepam1_IVP': {
+        'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'rezultati','run2')
+    },
+    'cDiazepamF_IVP': {
+        'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_ociscen', 'rezultati','run3')
+    }
+}
+
+
+
+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
+
+