import pandas as pd import numpy as np from lxml import etree as et from lxml import objectify import math import LogPostProcessing import os import subprocess from scipy import optimize import sys def change_in_situ_strain(ligament_name, febio_spec_root, in_situ_strain): LoadData = febio_spec_root.find("LoadData") ligament_lc = None for lc in LoadData: try: if ligament_name in lc.attrib["name"]: ligament_lc = lc break except KeyError: # just in case a name is not defined for one of the load curves pass # change the in situ strain value for all times greater than 1 for point in ligament_lc: time = int(point.text.split(',')[0]) if time >= 1.0: point.text = str(time)+ ',' + str(in_situ_strain) def get_experiment_data(experiment_kinematics_csv, loading_or_all, ligament_name=None): # # get the experiment kinematics we are comparing with # dominant_axes = {"ACL":"Knee JCS Anterior Translation [mm]", "PCL":"Knee JCS Anterior Translation [mm]", # "MCL":"Knee JCS Adduction Rotation [deg]", "LCL": "Knee JCS Adduction Rotation [deg]"} df = pd.read_csv(experiment_kinematics_csv, encoding='utf7') if loading_or_all == "loading": # find the name of the data along the dominant loading axis depending on the ligament dominant_axes = {} column_labels = df.keys() for n in column_labels: if "Adduction" in n or "Abduction" in n: dominant_axes["MCL"] = n dominant_axes["LCL"] = n elif "Anterior" in n or "Posterior" in n: dominant_axes["ACL"] = n dominant_axes["PCL"] = n kinematic_data = list(map(float, df[dominant_axes[ligament_name]].values)) kinematic_data = np.asarray(kinematic_data) elif loading_or_all == "all": # get data from all kinematics channels kinematic_data = [None, None, None, None, None, None] column_labels = df.keys() for n in column_labels: i = None if "Medial" in n or "Lateral" in n: i=0 elif "Anterior" in n or "Posterior" in n: i=1 elif "Superior" in n or "Inferior" in n: i=2 elif "Flexion" in n or "Extension" in n: i = 3 elif "Adduction" in n or "Abduction" in n: i= 4 elif "External" in n or "Internal" in n: i=5 if i is not None: kinematic_data[i] = list(map(float, df[n].values)) kinematic_data = np.asarray(kinematic_data) else: print("unknown rms_error type. options are 'loading' or 'all'") return None return kinematic_data def get_model_data(febio_file, febio_spec_root, loading_or_all, ligament_name=None, relative = False): directory = os.path.dirname(febio_file) results_directory = os.path.join(directory, 'Processed_Results' ) model_kinematics_csv = os.path.join(results_directory, 'Tibiofemoral_Kinematics.csv') try: df = pd.read_csv(model_kinematics_csv, encoding='utf7') except: df = pd.read_csv(model_kinematics_csv, encoding='utf8') time_data = list(map(float, df['# time_steps'].values)) time_data = np.asarray(time_data) if relative: t = 1.0 else: t = 2.0 print(t) # want to extract only the data from time point t onwards as this is the loading time we wish to compare to experimental data idx_2 = np.where(time_data >= t)[0] time_data = time_data[idx_2] # pull out only the time points that were set in the load curves of the febio file LoadData = febio_spec_root.find("LoadData") relevant_time_points = [] for lc in LoadData: if "Load " in lc.attrib["name"]: for point in lc: time_point = float(point.text.split(',')[0]) if time_point >= t: relevant_time_points.append(time_point) break if relative: relevant_time_points = relevant_time_points[1:] # hack solution, # time 1.0 is requested first for prestrain, and then again fro loading. ingnore the first if loading_or_all == "loading": # find the name of the data along the dominant loading axis depending on the ligament dominant_axes = {} column_labels = df.keys() for n in column_labels: if "Adduction" in n: if "Translation" in n: dominant_axes["ACL"] = n dominant_axes["PCL"] = n elif "Rotation" in n: dominant_axes["MCL"] = n dominant_axes["LCL"] = n kinematic_data = list(map(float, df[dominant_axes[ligament_name]].values)) kinematic_data = np.asarray(kinematic_data) kinematic_data = kinematic_data[idx_2] model_data = [] # extract only the data at the relevant time points, add nan if data is missing at that time point for i in relevant_time_points: try: relevant_idx = np.where(time_data == i)[0][0] model_data.append(kinematic_data[relevant_idx]) except IndexError: # if the model did not solve at that time model_data.append(np.nan) model_data = np.asarray(model_data) elif loading_or_all == "all": kinematic_data = [None, None, None, None, None, None] for name in df.columns: i = None if "force" in name.lower() or "load" in name.lower() or "translation" in name.lower(): if "extension" in name.lower(): i = 0 elif "adduction" in name.lower(): i = 1 elif "internal" in name.lower(): i = 2 elif "moment" in name.lower() or "rotation" in name.lower(): if "extension" in name.lower(): i = 3 elif "adduction" in name.lower(): i = 4 elif "internal" in name.lower(): i = 5 if i is not None: kinematic_data[i] = np.asarray(list(map(float, df[name].values))) kinematic_data = np.asarray(kinematic_data) kinematic_data = kinematic_data[:,idx_2] model_data = [] # extract only the data at the relevant time points, add nan if data is missing at that time point for i in relevant_time_points: try: relevant_idx = np.where(time_data == i)[0][0] model_data.append(kinematic_data[:,relevant_idx]) except IndexError: # if the model did not solve at that time model_data.append([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]) model_data = np.asarray(model_data) model_data = model_data.T else: print("unknown rms_error type. options are 'loading' or 'all'") return None return model_data def write_file(xml_tree, new_filename): # Write the New File xml_tree.write(new_filename, xml_declaration=True, pretty_print=True) # make sure that it's pretty printing properly, by parsing and re-writing parser = et.XMLParser(remove_blank_text=True) new_feb_tree = et.parse(new_filename, parser) new_feb_tree.write(new_filename, xml_declaration=True, pretty_print=True) def Find_InSituStrain(febio_file, experiment_kinematics_csv, model_properties_xml, ligament_name, run, initial_guess, rms_type, febioCommand, relative =False): """given an febio knee model, the ligament we want to optimize, find the in situ strain value to minize the difference in kinematics along the dominant loading axis """ # create a text file to print useful info to file_path = os.path.dirname(febio_file) optimization_file = os.path.join(file_path, 'Optimization{}.txt'.format(run) ) # file_path = ('/').join(febio_file.split('/')[:-1]) # optimization_file = file_path + '/Optimization{}.txt'.format(run) # get the kinematics data to compare with experiment_data = get_experiment_data(experiment_kinematics_csv, rms_type, ligament_name) # remove nans from experiment data if rms_type == "loading": experiment_data = experiment_data[~np.isnan(experiment_data)] elif rms_type == "all": new_experiment_data = [] for row in experiment_data: new_row = row[~np.isnan(row)] new_experiment_data.append(new_row) experiment_data = np.asarray(new_experiment_data) def get_rms(in_situ_strain, febio_file, ligament_name, experiment_data): info_file = open(optimization_file, "a") info_file.write('In Situ Strain = {} \n'.format(in_situ_strain)) # check that in situ strain is not negative, if yes, return large rms if in_situ_strain < 0 or in_situ_strain > 2.0: info_file.write('warning: negative in situ strain, returning large rms') info_file.close() return 100000 info_file.close() # close file so we can open and look at it while the models are running Febio_tree = et.parse(febio_file) febio_spec_root = Febio_tree.getroot() # make sure the full path to geometry file is given Geometry = febio_spec_root.find("Geometry") Geometry.attrib["from"] = os.path.join(file_path , 'Geometry_custom.feb') # Geometry.attrib["from"] = file_path + '/Geometry_custom.feb' # change the in situ strain of the ligament of interest change_in_situ_strain(ligament_name, febio_spec_root, in_situ_strain) # rewrite the file - overwrite existing febio file write_file(Febio_tree, febio_file) # run the model subprocess.call([febioCommand, febio_file]) # run post-processing script logfile = febio_file.replace('.feb', '.log') try: LogPostProcessing.MakeGraphs(logfile, model_properties_xml) # get the kinematics data from the model model_data = get_model_data(febio_file, febio_spec_root, rms_type, ligament_name, relative=relative) # any missing data is nan except: # post processing will fail if model did not converge even one time step. In that case, set model data to nans model_data = np.empty(np.shape(experiment_data)) model_data[:] = np.nan # note - nans have been removed from experiment data already, as loads were not specified for those # time points in the model. any nans in model data means the model did not solve at that requested time point # find the rms error dif = model_data - experiment_data dif2 = np.square(dif) # if any model data is nan, sums will return nan, remove any nan values # - this converts dif2 to a 1d array if all square differences (if rms_type is all this included all kinematics axes) dif2 = dif2[~np.isnan(dif2)] rms = np.sqrt((np.sum(dif2))/len(dif2)) # if the rms error is nan ie the model did not solve at any of the requested time points (likely failed to converge) # return large rms error if math.isnan(rms): rms = 100000 # store all the results in the text file info_file = open(optimization_file, "a") info_file.write('Model data for rms calculation: {} \n'.format(model_data)) info_file.write('RMS error = {} \n'.format(rms)) info_file.write('\n') info_file.write('=====================================================================') info_file.write('\n') info_file.close() return rms # res = optimize.minimize_scalar(get_rms, bounds = (1.0, 1.25), args=(febio_file, ligament_name, experiment_data), method='Bounded', tol=0.001, options={'xatol': 0.001, 'maxiter': 30}) # this is done to allow for quicker convergence when the solution is very close to the guess. # has little effect on speed of convergence when solution is not close to guess a= initial_guess c= initial_guess - 0.001 res = optimize.minimize_scalar(get_rms, bracket = (a, c), args=(febio_file, ligament_name, experiment_data), method='brent', tol=0.001, options={'xatol': 0.001, 'maxiter': 30}) opt_in_situ_strain = res.x opt_in_situ_strain = round(opt_in_situ_strain, 3) # round the result to 3 decimal places # write the final optimized in situ strain result to a text file info_file = open(optimization_file, "a") info_file.write('\n') info_file.write('The optimized in situ strain is {}'.format(opt_in_situ_strain)) try: num_iters= res.nit success = res.success term_message = res.message info_file.write('termination message: ' + term_message) info_file.write('success: {}'.format(success)) info_file.write('# iters: {}'.format(num_iters)) except: pass info_file.close() return opt_in_situ_strain def Find_All_InSituStrains(febio_files, experiment_kinematics_files, model_properties_xml, strains_dict, rms_type, febioCommand, relative=False): # create a text file to print useful info to file_path = os.path.dirname(febio_files["ACL"]) optimization_file = os.path.join(file_path, 'Optimization_MultiVar.txt') # get the kinematics data to compare with all_experiment_data = {} for ligament_name in experiment_kinematics_files.keys(): experiment_data = get_experiment_data(experiment_kinematics_files[ligament_name], rms_type, ligament_name) # remove nans from experiment data if rms_type == "loading": experiment_data = experiment_data[~np.isnan(experiment_data)] elif rms_type == "all": new_experiment_data = [] for row in experiment_data: new_row = row[~np.isnan(row)] new_experiment_data.append(new_row) experiment_data = np.asarray(new_experiment_data) all_experiment_data[ligament_name] = experiment_data # we want to run multivariable optimization where the input is all 4 in situ strains, and the output is the rms error # (either along the loading direction or along all directions, depending on rms_type) # input for multi variable optimization function must be as a list, we will give all inputs as lists in the same order all_ligaments = strains_dict.keys() lig_names = list(all_ligaments) strains_list = [] febio_files_list = [] all_experiment_data_list = [] for nm in lig_names: strains_list.append(strains_dict[nm]) febio_files_list.append(febio_files[nm]) all_experiment_data_list.append(all_experiment_data[nm]) def get_rms_multi(strains_list, febio_files_list, all_experiment_data_list, lig_names): info_file = open(optimization_file, "a") info_file.write('In Situ Strain Values = {} \n'.format(strains_list)) # check if any in situ strains are less than 0 or larger than 2 if any(t<0 for t in strains_list) or any(b>2 for b in strains_list): info_file.write('warning: invalid in situ strain guess, returning large rms \n') info_file.close() return 10000 info_file.close() # close file so we can open and look at it while the models are running # for each febio file, change in situ strains of all ligaments to those in the inputs, run the file, # get model kinematics, calculate calculate the squared differences between model and experiment # kinematics on all axes. all_sqaured_difs = [] for a in range(len(lig_names)): febio_file = febio_files_list[a] ligament_name= lig_names[a] experiment_data = all_experiment_data_list[a] Febio_tree = et.parse(febio_file) febio_spec_root = Febio_tree.getroot() # make sure the full path to geometry file is given Geometry = febio_spec_root.find("Geometry") Geometry.attrib["from"] = os.path.join(os.path.dirname(febio_file), 'Geometry_custom.feb') # change the in situ strain of ALL ligaments for number, name in enumerate(lig_names): change_in_situ_strain(name, febio_spec_root, strains_list[number]) # rewrite the file - overwrite existing febio file write_file(Febio_tree, febio_file) # run the model subprocess.call([febioCommand, febio_file]) # run post-processing script logfile = febio_file.replace('.feb', '.log') try: LogPostProcessing.MakeGraphs(logfile, model_properties_xml) # get the kinematics data from the model model_data = get_model_data(febio_file, febio_spec_root, rms_type, ligament_name, relative=relative) # any missing data is nan except: # post processing will fail when model didnt converge one time step. set all model data to nans model_data = np.empty(np.shape(experiment_data)) model_data[:] = np.nan info_file = open(optimization_file, "a") info_file.write('Model kinematics for {} model: {} \n'.format(ligament_name, model_data)) # note - nans have been removed from experiment data already, as loads were not specified for those # time points in the model. any nans in model data means the model did not solve at that requested time point # find the rms error - will be 1d or 2d array depending on rms_type dif = model_data - experiment_data dif2 = np.square(dif) # if any model data is nan, sums will return nan, remove any nan values - this will turn it into a 1d array dif2 = dif2[~np.isnan(dif2)] # if dif2 is all nans - model didnt converge at all - return, penalize rms instead of continuing to run models if len(dif2) == 0: info_file.write("{} model didn't converge. exiting iteration and penalizing output \n".format(ligament_name)) info_file.write("RMS error = 10000 \n") info_file.write('\n') info_file.write('=====================================================================') info_file.write('\n') info_file.close() return 10000 info_file.close() all_sqaured_difs.append(dif2) # calculate rms error for ALL kinematics channels for ALL 4 models that were run all_sqaured_difs = [item for sublist in all_sqaured_difs for item in sublist] # flatten the list of squared differences all_sqaured_difs = np.asarray(all_sqaured_difs) rms = np.sqrt((np.sum(all_sqaured_difs)) / len(all_sqaured_difs)) # store all the results in the text file info_file = open(optimization_file, "a") info_file.write('RMS error = {} \n'.format(rms)) info_file.write('\n') info_file.write('=====================================================================') info_file.write('\n') info_file.close() return rms # res = optimize.minimize(get_rms_multi, x0=np.array(strains_list), args=(febio_files_list,all_experiment_data_list,lig_names), method='L-BFGS-B', bounds=[(0, 2), (0, 2), (0, 2), (0, 2)], tol=0.001, # options={'xatol': 0.001}) # res = optimize.minimize(get_rms_multi, x0=np.array(strains_list), args=(febio_files_list,all_experiment_data_list,lig_names), method='SLSQP', options={'xatol': 0.001}) opt_in_situ_strain = res.x # write the final optimized in situ strain result to a text file info_file = open(optimization_file, "a") info_file.write('\n') info_file.write('The optimized in situ strains are {}'.format(opt_in_situ_strain)) def update_in_situ_strains(febio_file, strains_dict): # open the file, Febio_tree = et.parse(febio_file) febio_spec_root = Febio_tree.getroot() # change the in situ strains for ligament_name, in_situ_strain in iter(strains_dict.items()): change_in_situ_strain(ligament_name, febio_spec_root, in_situ_strain) # rewrite the file - overwrite existing febio file write_file(Febio_tree, febio_file) def run_from_xml(xml_file, febio_command = None, relative=False): """ read the xml file containing the information about the ligaments and files, and run the optimization""" if not febio_command: febio_command = '/home/schwara2/Programs/FEBio/FEBio2.8.2/bin/febio2.lnx64' # # in situ strain intial guesses # strains_dict = {"ACL": 1.016, "PCL": 1.0, "LCL": 1.027, "MCL": 1.034} - there were the initial guesses from lit. Optimization_tree = et.parse(xml_file) Optimization = Optimization_tree.getroot() # get the model properties file name files = Optimization.find("Files") model_properties_file = files.find('model_properties').text # get the options for rms error - loading or all, and opt- single or multi try: options = Optimization.find("Options") try: rms_error = options.find("rms_error") rms_type = rms_error.attrib["type"] except: rms_type = "loading" try: opt = options.find("opt") opt_type = opt.attrib["type"] except: opt_type = "single" except: print("no Options element found in xml, default options will be used") rms_type = "loading" opt_type = "single" ligaments = Optimization.find("Ligaments") strains_dict = {} # if we already know the iss for a ligament, can exclude from the xml and put the known iss here febio_files = {} kinematics_files = {} ligaments_to_optimize = [] # go through each ligament in the file and collect the data for lig in ligaments: if lig.tag is et.Comment: continue ligament_name = lig.tag ligaments_to_optimize.append(ligament_name) febio_file = lig.find('febio_file') febio_files[ligament_name] = febio_file.text experiment_kinematics_csv = lig.find('experiment_kinematics') kinematics_files[ligament_name] = experiment_kinematics_csv.text intitial_iss = lig.find('in_situ_strain') strains_dict[ligament_name] = float(intitial_iss.text) # for single varibale optimization - #run the optimization for each ligament specified in the xml, in the order they appear in the document, # continue looping until ligmanet prestrains dont change by more than 0.001 if opt_type == "single": run = 0 all_converged = False while not all_converged: converged_list = [] for l in ligaments_to_optimize: update_in_situ_strains(febio_files[l], strains_dict) prev_iss = strains_dict[l] lig_iss = Find_InSituStrain(febio_files[l], kinematics_files[l], model_properties_file, ligament_name=l, run=run, initial_guess=prev_iss, rms_type=rms_type, febioCommand=febio_command, relative=relative) # make sure the lig iss was outpt as a single variable not as a list strains_dict[l] = lig_iss change = abs(prev_iss - lig_iss) change = round(change, 5) # to fix issues with python rounding counting 0.001 as not <= 0.001 if change <= 0.001: converged_list.append(l) # check if all ligaments have converged if len(converged_list) == len(ligaments_to_optimize): all_converged = True run += 1 elif opt_type == "multi": # For multi valriable optimization, we will want to optimize all in situ strains simultaneously # so that the rms error for the given models are minimized - can still do rms type loading or all? Find_All_InSituStrains(febio_files, kinematics_files, model_properties_file, strains_dict, rms_type, febioCommand=febio_command, relative=relative) else: print("exiting. opt type not found. use 'single' or 'multi' ") def test(xml_file): Optimization_tree = et.parse(xml_file) Optimization = Optimization_tree.getroot() # get the model properties file name files = Optimization.find("Files") model_properties_file = files.find('model_properties').text # get the options for rms error - loading or all, and opt- single or multi try: options = Optimization.find("Options") try: rms_error = options.find("rms_error") rms_type = rms_error.attrib["type"] except: rms_type = "loading" try: opt = options.find("opt") opt_type = opt.attrib["type"] except: opt_type = "single" except: print("no Options element found in xml, default options will be used") rms_type = "loading" opt_type = "single" ligaments = Optimization.find("Ligaments") strains_dict = {} # if we already know the iss for a ligament, can exclude from the xml and put the known iss here febio_files = {} kinematics_files = {} ligaments_to_optimize = [] # go through each ligament in the file and collect the data for lig in ligaments: if lig.tag is et.Comment: continue ligament_name = lig.tag ligaments_to_optimize.append(ligament_name) febio_file = lig.find('febio_file') febio_files[ligament_name] = febio_file.text experiment_kinematics_csv = lig.find('experiment_kinematics') kinematics_files[ligament_name] = experiment_kinematics_csv.text intitial_iss = lig.find('in_situ_strain') strains_dict[ligament_name] = float(intitial_iss.text) # check reading of experiment and model data, and finding rms error # for single opt type - checked ligament = 'MCL' feb_file = febio_files[ligament] kinematics_csv = kinematics_files[ligament] # open the file, Febio_tree = et.parse(feb_file) febio_spec_root = Febio_tree.getroot() exp= get_experiment_data(kinematics_csv, rms_type, ligament) # remove nans from experiment data if rms_type == "loading": exp = exp[~np.isnan(exp)] elif rms_type == "all": new_experiment_data = [] for row in exp: new_row = row[~np.isnan(row)] new_experiment_data.append(new_row) exp = np.asarray(new_experiment_data) try: mod = get_model_data(feb_file, febio_spec_root, rms_type, ligament, relative=False) # any missing data is nan print("successfully loaded model data") except: # post processing will fail if model did not converge even one time step. In that case, set model data to nans print("loading model data unsuccessful, creating nans same shape as experiment data") mod = np.empty(np.shape(exp)) mod[:] = np.nan print(mod) # find the rms error dif = mod - exp dif2 = np.square(dif) # if any model data is nan, sums will return nan, remove any nan values # - this converts dif2 to a 1d array if all square differences (if rms_type is all this included all kinematics axes) dif2 = dif2[~np.isnan(dif2)] rms = np.sqrt((np.sum(dif2)) / len(dif2)) print(rms) # # for multi opt type - checked # # all_experiment_data = {} # # for ligament_name in kinematics_files.keys(): # experiment_data = get_experiment_data(kinematics_files[ligament_name], rms_type, ligament_name) # # # remove nans from experiment data # if rms_type == "loading": # experiment_data = experiment_data[~np.isnan(experiment_data)] # elif rms_type == "all": # new_experiment_data = [] # for row in experiment_data: # new_row = row[~np.isnan(row)] # new_experiment_data.append(new_row) # experiment_data = np.asarray(new_experiment_data) # # all_experiment_data[ligament_name] = experiment_data # # all_sqaured_difs= [] # # all_ligaments = strains_dict.keys() # lig_names = list(all_ligaments) # strains_list = [] # febio_files_list = [] # all_experiment_data_list = [] # # for nm in lig_names: # strains_list.append(strains_dict[nm]) # febio_files_list.append(febio_files[nm]) # all_experiment_data_list.append(all_experiment_data[nm]) # # for a in range(len(lig_names)): # # febio_file = febio_files_list[a] # ligament_name = lig_names[a] # experiment_data = all_experiment_data_list[a] # # Febio_tree = et.parse(febio_file) # febio_spec_root = Febio_tree.getroot() # # # get the kinematics data from the model # model_data = get_model_data(febio_file, # febio_spec_root, rms_type, ligament_name) # any missing data is nan # # # find the rms error - will be 1d or 2d array depending on rms_type # dif = model_data - experiment_data # dif2 = np.square(dif) # # # if any model data is nan, sums will return nan, remove any nan values # # this converts dif2 to a 1d array if all square differences (if rms_type is all this includes all kinematics axes) # # dif2 = dif2[~np.isnan(dif2)] # # dif2 = dif2[:, ~np.isnan(dif2).any(axis=0)] # # # if dif2 is all nans - model didnt converge at all - return, penalize rms instead of continuing to run models # if len(dif2) == 0: # rms= 1000000 # # all_sqaured_difs.append(dif2) # # # calculate rms error for ALL kinematics channels for ALL 4 models that were run # all_sqaured_difs = [item for sublist in all_sqaured_difs for item in # sublist] # flatten the list of squared differences # all_sqaured_difs = np.asarray(all_sqaured_difs) # rms = np.sqrt((np.sum(all_sqaured_difs)) / len(all_sqaured_difs)) print(rms) if __name__ == '__main__': # on ariels ccf laptop # febio command = "C:\Users\schwara2\Programs\Febio\bin\FEBio2.exe" # run_from_xml(*sys.argv[1:]) # xml_file = "C:\\Users\schwara2\Documents\Open_Knees\Optimization\InSituOpt.xml" # test(*sys.argv[1:])