# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import pymc as pm import scipy.stats as st import chaospy as cp import itertools import pandas as pd import json import pickle from skopt import gp_minimize from skopt.space import Real from skopt.utils import use_named_args import arviz as az import glob import matplotlib.patches as mpatches def generate_points(num_points=80, PowerLaw_params=None, sigma=0.15, seed=42): """ Generate points for the given number of points, PowerLaw_params , sigma, and seed. Parameters: - num_points: Number of points to generate. - PowerLaw_params: Parameters for PowerLaw model. - sigma: Standard deviation for noise. - seed: Random seed for reproducibility. Returns: - mesh_data: Tuple containing mesh grid data. - data_extended: Dictionary containing extended data. """ if PowerLaw_params is None: PowerLaw_params = {'C': 1.228e-5, 'alpha': 0.6606, 'beta': 1.9918} np.random.seed(seed) points_per_side = int(np.sqrt(num_points)) tau_values = np.linspace(50, 320, points_per_side) time_values = np.linspace(0.039, 1.48, points_per_side) tau_mesh, time_mesh = np.meshgrid(tau_values, time_values) hi_values = PowerLaw_params['C'] * tau_mesh**PowerLaw_params['beta'] * time_mesh**PowerLaw_params['alpha'] tau_flat = tau_mesh.flatten() time_flat = time_mesh.flatten() hi_flat = hi_values.flatten() tau_extended = np.tile(tau_flat, 3) time_extended = np.tile(time_flat, 3) hi_extended = np.concatenate([hi_flat, np.random.normal(hi_flat, sigma, len(hi_flat)), np.random.normal(hi_flat, sigma, len(hi_flat))]) data_extended = { 'Shear Stress in Pa': tau_extended, 'Exposure Time in s': time_extended, 'HI in %': hi_extended } mesh_data = (tau_mesh, time_mesh, hi_values) return mesh_data, data_extended def plot_points(mesh_data, data_extended): """ Plot points using mesh data and extended data. Parameters: - mesh_data: Tuple containing mesh grid data. - data_extended: Dictionary containing extended data. """ tau_mesh, time_mesh, hi_values = mesh_data tau_extended = data_extended['Shear Stress in Pa'] time_extended = data_extended['Exposure Time in s'] hi_extended = data_extended['HI in %'] fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(111, projection='3d') ax.plot_surface(tau_mesh, time_mesh, hi_values, color='grey', alpha=0.2) ax.scatter(tau_extended, time_extended, hi_extended, color='black', s=10, edgecolors='none') ax.invert_xaxis() ax.set_xlabel('Shear Stress (Pa)', fontsize=10) ax.set_ylabel('Exposure Time (s)', fontsize=10) ax.set_zlabel('HI (%)', fontsize=10) plt.show() def run_mcmc(df, distribution='student', C_fixed=None, beta_fixed=None, mcmc_settings=None): """ Run MCMC with flexible settings for distribution type, fixed parameters, and MCMC configurations. Parameters: - df: DataFrame containing the data. - distribution: 'normal' or 'student'. Determines the likelihood distribution. - C_fixed: Fixed value for C if not None. Otherwise, C is treated as a Uniform distribution. - beta_fixed: Fixed value for beta if not None. Otherwise, beta is treated as a Normal distribution. - mcmc_settings: Dictionary with MCMC settings (e.g., samples, tune, target_accept, chains). If None, defaults are used. Returns: - trace: Trace object containing MCMC samples. """ if mcmc_settings is None: mcmc_settings = {'samples': 10000, 'tune': 1000, 'target_accept': 0.95, 'chains': 2} with pm.Model() as model: if C_fixed is not None: C = C_fixed else: C = pm.Uniform('C', lower=0, upper=1) alpha = pm.Normal('alpha', mu=0, sigma=1) sigma = pm.Normal('sigma', mu=0, sigma=1) if beta_fixed is not None: beta = beta_fixed else: beta = pm.Normal('beta', mu=0, sigma=1) hemolysis_est = C * np.array(df["Shear Stress in Pa"])**beta * np.array(df["Exposure Time in s"])**alpha if distribution == 'normal': likelihood = pm.Normal('likelihood', mu=hemolysis_est, sigma=sigma, observed=np.array(df['HI in %'])) elif distribution == 'student': nu = pm.Exponential('nu', 1/30) likelihood = pm.StudentT('likelihood', mu=hemolysis_est, sigma=sigma, nu=nu, observed=np.array(df['HI in %'])) trace = pm.sample( mcmc_settings['samples'], tune=mcmc_settings['tune'], target_accept=mcmc_settings['target_accept'], chains=mcmc_settings['chains'], return_inferencedata=True ) return trace def create_and_fit_distributions(trace, parameter_names, confidence_level=0.99): """ Create and fit distributions for MCMC trace parameters. Parameters: - trace: Trace object containing MCMC samples. - parameter_names: List of parameter names to fit distributions. - confidence_level: Confidence level for calculating confidence intervals. Returns: - distributions: Dictionary containing fitted distributions. - CI_bounds: Dictionary containing confidence interval bounds. """ distributions = {} CI_bounds = {} for param in parameter_names: samples = trace.posterior[param].values.flatten() dist_fit = st.lognorm.fit(samples, floc=0) CI_bounds[param] = st.lognorm.interval(confidence_level, *dist_fit) distribution_cp = cp.LogNormal(np.log(dist_fit[-1]), dist_fit[0]) distributions[param] = { "fit": dist_fit, "distribution": cp.LogNormal(np.log(dist_fit[-1]), dist_fit[0]), "truncated_distribution": cp.Trunc(distribution_cp, lower=CI_bounds[param][0], upper=CI_bounds[param][1]) } return distributions, CI_bounds def save_distributions_and_CI_bounds(distributions, CI_bounds, output_dir): """ Save distribution parameters and confidence interval bounds to JSON files. Parameters: - distributions: Dictionary of fitted distribution details. - CI_bounds: Dictionary of confidence interval bounds. - output_dir: Directory where the output files will be saved. """ dist_params = {param: {"fit": dist_info["fit"]} for param, dist_info in distributions.items()} ci_data = {param: list(bounds) for param, bounds in CI_bounds.items()} dist_filename = output_dir + "_distribution_params.json" with open(dist_filename, 'w') as f: json.dump(dist_params, f, indent=4) ci_filename = output_dir + "_CI_bounds.json" with open(ci_filename, 'w') as f: json.dump(ci_data, f, indent=4) def generate_samples(distributions, CI_bounds, num_samples, add_CI_combinations=False, sampling_method="equidistant", seed=None): """ Generate samples from truncated distributions with specified sampling method. Parameters: - distributions: Dictionary of truncated distributions for each parameter. - CI_bounds: Dictionary of confidence interval bounds for each parameter. - num_samples: Total number of samples to generate (for equidistant, defines the total grid size). - add_CI_combinations: Whether to add all combinations of CI bounds to the samples. - sampling_method: Specifies the sampling method ('equidistant' or 'latin_hypercube'). - seed: Optional seed for reproducibility. Returns: - samples_final: Final set of samples, including optional CI combinations. """ samples_final = [] relevant_parameters = [param for param in CI_bounds if param in ['C', 'alpha', 'beta']] if sampling_method == "latin_hypercube": truncated_distributions = [distributions[param]["truncated_distribution"] for param in relevant_parameters] joint_distribution_trunc = cp.J(*truncated_distributions) samples_final = joint_distribution_trunc.sample(num_samples, rule="latin_hypercube", seed=seed).T elif sampling_method == "random": truncated_distributions = [distributions[param]["truncated_distribution"] for param in relevant_parameters] joint_distribution_trunc = cp.J(*truncated_distributions) samples_final = joint_distribution_trunc.sample(num_samples, rule="random", seed=seed).T elif sampling_method == "equidistant": samples_final = generate_equidistant_samples_in_CI_bounds(CI_bounds, num_samples) if add_CI_combinations: combinations = list(itertools.product(*[CI_bounds[param] for param in relevant_parameters])) samples_final = np.concatenate((samples_final, np.array(combinations))) return samples_final def generate_equidistant_samples_in_CI_bounds(CI_bounds, num_samples): """ Generate a grid of equidistant samples within the confidence interval bounds. Parameters: - CI_bounds: Dictionary of confidence interval bounds for each parameter. - num_samples: Total number of samples desired. Returns: - samples: Array of equidistant samples within the CI bounds. """ relevant_parameters = [param for param in CI_bounds if param in ['C', 'alpha', 'beta']] num_dimensions = len(relevant_parameters) samples_per_dim = int(np.round(num_samples ** (1 / num_dimensions))) equidistant_samples = [np.linspace(CI_bounds[param][0], CI_bounds[param][1], samples_per_dim) for param in relevant_parameters] mesh = np.array(np.meshgrid(*equidistant_samples)).reshape(num_dimensions, -1).T return mesh def process_trace_and_generate_Training_Data_Samples(trace, output_dir, num_samples=200): """ Process MCMC trace and generate samples. Parameters: - trace: Trace object containing MCMC samples. - output_dir: Directory where the output files will be saved. - num_samples: Number of samples to generate. """ output_dir = output_dir.split('.n')[0] parameter_names = list(trace.posterior.data_vars) relevant_parameters = [param for param in parameter_names if param in ['C', 'alpha', 'beta', 'sigma']] relevant_parameters_sample_data = [param for param in parameter_names if param in ['C', 'alpha', 'beta']] distributions, CI_bounds = create_and_fit_distributions(trace, relevant_parameters) # save_distributions_and_CI_bounds(distributions, CI_bounds, output_dir) latin_hypercube_samples = generate_samples(distributions, CI_bounds, num_samples, add_CI_combinations=True, sampling_method="latin_hypercube") pd.DataFrame(latin_hypercube_samples, columns=relevant_parameters_sample_data).to_csv(output_dir + '_latin_hypercube_samples.csv', index=False) equidistant_samples = generate_samples(distributions, CI_bounds, num_samples, add_CI_combinations=True, sampling_method="equidistant") pd.DataFrame(equidistant_samples, columns=relevant_parameters_sample_data).to_csv(output_dir + '_equidistant_samples.csv', index=False) EvaluationSamples = generate_samples(distributions, CI_bounds, 10000, add_CI_combinations=False,sampling_method="random",seed=42) pd.DataFrame(EvaluationSamples, columns=relevant_parameters_sample_data).to_csv(output_dir + '_EvaluationSamples.csv', index=False) def find_optimal_C(data_extended, TracePath): """ Find the optimal value of C using Bayesian Optimization. Parameters: - data_extended: Extended dataset to be used in MCMC sampling. - TracePath: Path where the result object (res_gp) will be saved. Returns: - res_gp: The result object from the Bayesian Optimization, containing the optimal C and other details. """ space = [Real(1e-6, 1e-3, "log-uniform", name='C')] C_values = [] objective_values = [] @use_named_args(space) def objective(**params): """ Objective function for Bayesian Optimization. Runs MCMC sampling and extracts the median of sigma. Parameters: - params: Dictionary containing the value of C. Returns: - sigma_median: The median value of sigma to be minimized. """ C_value = params['C'] trace = run_mcmc(data_extended, C_fixed=C_value) # Extract the median of sigma from the trace sigma_median = np.median(trace.posterior["sigma"].values.flatten()) C_values.append(C_value) objective_values.append(sigma_median) return sigma_median # Minimizing this objective # Perform Bayesian Optimization res_gp = gp_minimize(objective, space, n_calls=10, random_state=42) print(f"Optimal C: {res_gp.x[0]}") # Save the res_gp object to a file with open(TracePath + "res_gp.pkl", 'wb') as file: pickle.dump(res_gp, file) return res_gp def plot_optimization_process(res_gp): """ Plot the optimization process, showing the values of C and corresponding objective function values. Parameters: - res_gp: The result object from the Bayesian Optimization, containing the optimization history. """ # Extracting C values and corresponding objective function values C_values = [x[0] for x in res_gp.x_iters] objective_values = res_gp.func_vals # Find the index of the minimum objective function value min_obj_index = objective_values.argmin() # Use the index to find the corresponding C value min_obj_C_value = C_values[min_obj_index] # Creating the plot fig, ax = plt.subplots(figsize=(10, 6)) ax.scatter(C_values, objective_values, color='black') ax.plot(C_values, objective_values, color='grey', alpha=0.5) # Optional: connect points to visualize the sequence ax.set_xlabel('C') ax.set_ylabel('Objective Function Value \n'+r' (Median of $\sigma$)') ax.axvline(x=min_obj_C_value, color='red', linestyle='--', label=f'Optimal C= {min_obj_C_value:.3e}') ax.set_xscale('log') # Use log scale if C values vary over orders of magnitude plt.legend(fontsize=6) plt.show() ################################################################################################################################################# def load_json(filepath): """Utility function to load data from a JSON file.""" with open(filepath, 'r') as file: return json.load(file) def recreate_distributions_and_CIs(distributions_filepath, CI_bounds_filepath): """ Recreate distributions and confidence intervals (CIs) from saved parameters and CI bounds. Parameters: - distributions_filepath: Path to the JSON file containing the distribution parameters. - CI_bounds_filepath: Path to the JSON file containing the CI bounds. Returns: - recreated_distributions: Dictionary containing the recreated distributions, their CI bounds, and truncated distributions based on the CI bounds. """ dist_params = load_json(distributions_filepath) CI_bounds = load_json(CI_bounds_filepath) # Recreate the distributions using chaospy recreated_distributions = {} for param, params in dist_params.items(): # Extract parameters for LogNormal distribution s, loc, scale = params["fit"] # Recreate the chaospy distribution distribution_cp = cp.LogNormal(np.log(scale), s) recreated_distributions[param] = { "distribution_cp": distribution_cp, "CI_bounds": CI_bounds[param], "truncated_distribution":cp.Trunc(cp.LogNormal(np.log(scale), s),lower=CI_bounds[param][0],upper=CI_bounds[param][1])} return recreated_distributions def load_single_data(PathToFolder, NameKey): """ Load MCMC data and associated files for a single dataset. Parameters: - PathToFolder: Path to the folder containing the data files. - NameKey: Key name to identify the dataset. Returns: - MCMC_dict: Dictionary containing loaded MCMC trace, CI bounds, distribution parameters, equidistant samples, Latin hypercube samples, and recreated distributions. """ MCMC_dict = {} MCMC_dict[NameKey] = {} # Construct file paths nc_path = f"{PathToFolder}{NameKey}.nc" ci_bounds_path = f"{PathToFolder}{NameKey}_CI_bounds.json" distribution_params_path = f"{PathToFolder}{NameKey}_distribution_params.json" equidistant_samples_path = f"{PathToFolder}{NameKey}_equidistant_samples.csv" latin_hypercube_samples_path = f"{PathToFolder}{NameKey}_latin_hypercube_samples.csv" # Load data MCMC_dict[NameKey]['trace'] = az.from_netcdf(nc_path) MCMC_dict[NameKey]['ci_bounds'] = load_json(ci_bounds_path) MCMC_dict[NameKey]['distribution_params'] = load_json(distribution_params_path) MCMC_dict[NameKey]['equidistant_samples'] = pd.read_csv(equidistant_samples_path, delimiter=',', skiprows=1, names=["alpha", "beta"]) MCMC_dict[NameKey]['latin_hypercube_samples'] = pd.read_csv(latin_hypercube_samples_path, delimiter=',', skiprows=1, names=["alpha", "beta"]) MCMC_dict[NameKey]['CP_distributions'] = recreate_distributions_and_CIs(distribution_params_path,ci_bounds_path) return MCMC_dict def load_NICPE(PathToCFDPost, MCMC_dict, polyorder = 4): """ Load NICPE data from CSV files and fit polynomial chaos expansions. Parameters: - PathToCFDPost: Path to the folder containing CFD post-processing data files. - MCMC_dict: Dictionary containing MCMC data and distributions. - polyorder: Polynomial order for the chaos expansion (default is 4). Returns: - dataframes: List of dataframes containing the approximated solvers for each data file. - joint_distribution_trunc: Truncated joint distribution created from the truncated individual distributions. """ file_pattern = f"{PathToCFDPost}/*CFD_EvaluationsHe*.csv" # List to hold the dataframes dataframes = [] # Loop through the matched files and read them for file_path in glob.glob(file_pattern): df = pd.read_csv(file_path, delimiter=';') truncated_distributions = [MCMC_dict["CP_distributions"][param]['truncated_distribution'] for param in list(MCMC_dict["equidistant_samples"].columns)] joint_distribution_trunc = cp.J(*truncated_distributions) polynomial_expansion = cp.generate_expansion(polyorder, joint_distribution_trunc) approx_solver = cp.fit_regression(polynomial_expansion, np.array(df.iloc[:,:-1]).T, np.array(df.iloc[:,-1])) dataframes.append(approx_solver) return dataframes,joint_distribution_trunc # return dataframes def create_violin_plot(normalized_samples): fig, ax = plt.subplots(figsize=(2.5, 5)) parts = ax.violinplot(normalized_samples, showmeans=False, showmedians=False, showextrema=True) for pc in parts['bodies']: pc.set_facecolor('#fd8d3c') pc.set_edgecolor('#fd8d3c') pc.set_alpha(1) # Set alpha to 1 for full opacity for component in ['cmaxes', 'cmins', 'cbars']: pc = parts[component] pc.set_edgecolor('#fd8d3c') pc.set_linewidth(1) violin_legend = mpatches.Patch(color='#fd8d3c', label='MCMC') # Remove existing x-tick labels ax.set_xticks([]) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.set_ylim(0, 2) handles, labels = ax.get_legend_handles_labels() handles.append(violin_legend) ax.legend(handles=handles,loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=1,fontsize=15) ################################################################################################################################################# ########################################################## + + + Pre processing + + + ########################################################### ################################################################################################################################################# ''' To avoid repeatedly performing the computationally intensive MCMC calculations, it makes sense to split the script into a pre-processing part and a post-processing part. In the pre-processing part, computationally intensive MCMC calculations are performed and the results are saved at the desired location. Between the two parts, the CFD evaluation of the hemolysis model for the exported Trainingdatasamples must be performed within the CFD software. In the subsequent post-processing part, these calculations are then loaded. Pre-processing steps: + Generation of synthetic experimental data (only if you do not have your own experimental data) + Execution of MCMC (with flexible configuration options) + Creation of training data based on the MCMC results + Finding an optimal C if C was fixed in the MCMC ''' ############################################################################### ########################## Set paths to save the data ######################### ############################################################################### TracePath = "PathToFolder/Tutorial/" ############################################################################### ################ Generate and plot synthetic experimental data ################ ############################################################################### mesh_data, data_extended = generate_points() plot_points(mesh_data, data_extended) ############################################################################### ######################## Generate and save MCMC trace ######################### ############################################################################### # depending on the settings this will take a while, for this tutorial the number of samples and chains are reduced # for the publication these settings are used: {'samples': 50000, 'tune': 1000, 'target_accept': 0.95, 'chains': 4} trace_student_fix_C_1_228e5 = run_mcmc(data_extended, C_fixed=1.228e-5) ### export trace trace_student_fix_C_1_228e5.to_netcdf(TracePath+"data_extended.nc") ### create samples and export process_trace_and_generate_Training_Data_Samples(trace_student_fix_C_1_228e5,TracePath+"data_extended.nc") ############################################################################### ################################ Find optimal C ############################### ############################################################################### # For this tutorial the number of evaluations (n_calls) are reduced from 100 to 10 # res_gp = find_optimal_C(data_extended, TracePath) # Plot the optimization process # plot_optimization_process(res_gp) ### if optimal C is found, a new MCMC can be run and saved with run_mcmc() and process_trace_and_generate_Training_Data_Samples() ################################################################################################################################################# ########################################################## + + + Post processing + + + ########################################################## ################################################################################################################################################# ''' In this part, the CFD evaluation are loaded and plotted. If NIPCE should be used to accelerate the CFD evaluations, NIPCE is trained and evaluated. Post-processing steps: + Load CFD evaluations + (Train NIPCE) + Plot Hemolysis Distribution ''' ############################################################################### ########################## Set paths to save the data ######################### ############################################################################### TracePath = "PathToFolder/Tutorial/" ############################################################################### ################################# Loading data ################################ ############################################################################### '''Assuming that the CFD_Evaluations are stored as "CFD_EvaluationsHem.csv" file with first parameters and then corresponding hemolysis value''' CFD_Evaluations = pd.read_csv(TracePath+"CFD_EvaluationsHem.csv", delimiter=";") ############################################################################### ################################# Train NIPCE ################################# ############################################################################### data_extended_loaded = load_single_data(TracePath, "data_extended") data_extended_loaded_NIPCE = load_NICPE(TracePath, data_extended_loaded["data_extended"], polyorder = 4) NIPCE_Output = data_extended_loaded_NIPCE[0][0](*data_extended_loaded_NIPCE[1].sample(10000, seed=42)).T normalized_samples = NIPCE_Output / np.median(NIPCE_Output) ############################################################################### ################################### Plotting ################################## ############################################################################### ### For NIPCE create_violin_plot(normalized_samples) ### For CFD_Evaluations create_violin_plot(CFD_Evaluations.iloc[:,-1]/np.median(CFD_Evaluations.iloc[:,-1]))