diff --git a/src/bayesvalidrox/bayes_inference/bayes_inference.py b/src/bayesvalidrox/bayes_inference/bayes_inference.py
index c8073ba76ebeb3249646e17b7f395f8a2472e2f6..3246f405815383cdc26e2502f2e4e6e92292f404 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_inference.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_inference.py
@@ -331,7 +331,7 @@ class BayesInference:
for itr_idx, data in tqdm(
enumerate(self.perturbed_data),
total=self.n_bootstrap_itrs,
- desc="Boostraping the BME calculations", ascii=True
+ desc="Bootstrapping the BME calculations", ascii=True
):
# ---------------- Likelihood calculation ----------------
@@ -352,6 +352,16 @@ class BayesInference:
output_names[i]: data[j:k] for i, (j, k) in
enumerate(indices)
}
+ #print(output_names)
+ #print(indices)
+ #print(numbers)
+ #print(nobs)
+ #print(self.measured_data)
+ #for i, (j, k) in enumerate(indices):
+ # print(i,j,k)
+ #print(data)
+ #print(data_dict)
+ #stop
# Unknown sigma2
if opt_sigma == 'C' or hasattr(self, 'sigma2s'):
@@ -370,7 +380,7 @@ class BayesInference:
# BME (log)
log_BME[itr_idx] = np.log(
np.nanmean(np.exp(logLikelihoods[:, itr_idx],
- dtype=np.float128))
+ dtype=np.longdouble))#float128))
)
# BME correction when using Emulator
@@ -1042,9 +1052,9 @@ class BayesInference:
# Take the first column of Likelihoods (Observation data without noise)
if self.just_analysis or self.bayes_loocv:
index = self.n_tot_measurement-1
- likelihoods = np.exp(self.log_likes[:, index], dtype=np.float128)
+ likelihoods = np.exp(self.log_likes[:, index], dtype=np.longdouble)#np.float128)
else:
- likelihoods = np.exp(self.log_likes[:, 0], dtype=np.float128)
+ likelihoods = np.exp(self.log_likes[:, 0], dtype=np.longdouble)#np.float128)
n_samples = len(likelihoods)
norm_ikelihoods = likelihoods / np.max(likelihoods)
diff --git a/src/bayesvalidrox/bayes_inference/mcmc.py b/src/bayesvalidrox/bayes_inference/mcmc.py
index 6b9ca94823fb7064baa2f0588d0f97fb4c9d1d44..cc632c2de57140ab0ba4bd4c81a49316406dacd1 100755
--- a/src/bayesvalidrox/bayes_inference/mcmc.py
+++ b/src/bayesvalidrox/bayes_inference/mcmc.py
@@ -251,7 +251,7 @@ class MCMC:
# output current autocorrelation estimate
if self.verbose:
- print(f"Mean autocorr time estimate: {np.nanmean(tau):.3f}")
+ print(f"Mean autocorr. time estimate: {np.nanmean(tau):.3f}")
list_gr = np.round(self.gelman_rubin(sampler.chain), 3)
print("Gelman-Rubin Test*: ", list_gr)
@@ -282,15 +282,15 @@ class MCMC:
# Print summary
print('\n')
print('-'*15 + 'Posterior diagnostics' + '-'*15)
- print(f"mean auto-correlation time: {np.nanmean(tau):.3f}")
- print(f"thin: {thin}")
- print(f"burn-in: {burnin}")
- print(f"flat chain shape: {finalsamples.shape}")
- print(f"Mean acceptance fraction: {acc_fr:.3f}")
- print("Gelman-Rubin Test*: ", list_gr)
+ print(f"Mean auto-correlation time: {np.nanmean(tau):.3f}")
+ print(f"Thin: {thin}")
+ print(f"Burn-in: {burnin}")
+ print(f"Flat chain shape: {finalsamples.shape}")
+ print(f"Mean acceptance fraction*: {acc_fr:.3f}")
+ print("Gelman-Rubin Test**: ", list_gr)
print("\n* This value must lay between 0.234 and 0.5.")
- print("* These values must be smaller than 1.1.")
+ print("** These values must be smaller than 1.1.")
print('-'*50)
print(f"\n>>>> Bayesian inference with MCMC for {self.BayesOpts.name} "
@@ -777,7 +777,7 @@ class MCMC:
)
if ~np.isfinite(tmp['logml']):
warnings.warn(
- "logml could not be estimated within maxiter, rerunning with "
+ "Logml could not be estimated within maxiter, rerunning with "
"adjusted starting value. Estimate might be more variable than"
" usual.")
# use geometric mean as starting value
diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index 4db354bbd541affdc50a8dce37eb293188091f73..e1795a74ce7a7b7a7c9519b30b53e1ea11d4eab8 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -237,8 +237,8 @@ class PostProcessing:
elif samples is not None:
self.n_samples = samples.shape[0]
else:
- raise Exception("Please provide either samples or pass number of "
- "samples!")
+ raise Exception("Please provide either samples or pass the number"
+ " of samples!")
# Generate random samples if necessary
Samples = self._get_sample() if samples is None else samples
@@ -922,7 +922,7 @@ class PostProcessing:
ax = plt.gca()
_, p = stats.shapiro(residuals)
if p < 0.01:
- annText = "The residuals seem to come from Gaussian process."
+ annText = "The residuals seem to come from a Gaussian Process."
else:
annText = "The normality assumption may not hold."
at = AnchoredText(annText, prop=dict(size=30), frameon=True,
diff --git a/src/bayesvalidrox/pylink/pylink.py b/src/bayesvalidrox/pylink/pylink.py
index 7af4a14b8f8acb62c21862d5df4a527dc7fe57b9..84695119a2631ab5d37e7c5c2d78b0d5a25b5e09 100644
--- a/src/bayesvalidrox/pylink/pylink.py
+++ b/src/bayesvalidrox/pylink/pylink.py
@@ -170,9 +170,9 @@ class PyLinkForwardModel(object):
obs = self.observations
else:
raise Exception("Please provide the observation data as a "
- "dictionary via observations attribute or pass"
- " the csv-file path to MeasurementFile "
- "attribute")
+ "dictionary via the observations attribute or "
+ "pass the csv-file path to the MeasurementFile"
+ " attribute")
elif case.lower() == 'valid':
if isinstance(self.observations_valid, dict) and \
bool(self.observations_valid):
@@ -184,9 +184,9 @@ class PyLinkForwardModel(object):
obs = self.observations_valid
else:
raise Exception("Please provide the observation data as a "
- "dictionary via Observations attribute or pass"
- " the csv-file path to MeasurementFile "
- "attribute")
+ "dictionary via the observations attribute or "
+ "pass the csv-file path to the MeasurementFile"
+ " attribute")
# Compute the number of observation
n_obs = obs[self.Output.names].notnull().sum().values.sum()
@@ -316,9 +316,8 @@ class PyLinkForwardModel(object):
Process = os.system(f'./../{command}')
if Process != 0:
print('\nMessage 1:')
- print(f'\tIf value of \'{Process}\' is a non-zero value, '
- 'then compilation problems \n' % Process)
-
+ print(f'\tIf the value of \'{Process}\' is a non-zero value'
+ ', then compilation problems occur \n' % Process)
os.chdir("..")
# Read the output
@@ -438,6 +437,12 @@ class PyLinkForwardModel(object):
# Initilization
n_c_points = len(c_points)
all_outputs = {}
+
+ # TODO: if umbridge, just run, no parallel from here
+ # If the link type is UM-Bridge, then no parallel needs to be started from here
+ if self.link_type.lower() == 'umbridge':
+ # For now just do same as usual
+ self.link_type = 'function'
# Extract the function
if self.link_type.lower() == 'function':
@@ -499,7 +504,7 @@ class PyLinkForwardModel(object):
if len(NaN_idx) != 0:
print('\n')
print('*'*20)
- print("\nThe following parametersets have been removed:\n",
+ print("\nThe following parameter sets have been removed:\n",
c_points[NaN_idx])
print("\n")
print('*'*20)
@@ -656,6 +661,6 @@ class PyLinkForwardModel(object):
shutil.rmtree(path)
print("\n")
- print(f'{dir_name}.zip file has been created successfully!\n')
+ print(f'{dir_name}.zip has been created successfully!\n')
return
diff --git a/src/bayesvalidrox/surrogate_models/bayes_linear.py b/src/bayesvalidrox/surrogate_models/bayes_linear.py
index a7d6b5929a83fc89d15d7ab8f369187d0542923c..767e6fff73df0d524dcf94e620b0a67ae0742d83 100644
--- a/src/bayesvalidrox/surrogate_models/bayes_linear.py
+++ b/src/bayesvalidrox/surrogate_models/bayes_linear.py
@@ -138,7 +138,7 @@ class EBLinearRegression(BayesianLinearRegression):
normalize=True, perfect_fit_tol = 1e-6, alpha = 1, copy_X = True, verbose = False):
super(EBLinearRegression,self).__init__(n_iter, tol, fit_intercept, copy_X, verbose)
if optimizer not in ['em','fp']:
- raise ValueError('Optimizer can be either "em" of "fp" ')
+ raise ValueError('Optimizer can be either "em" or "fp" ')
self.optimizer = optimizer
self.alpha = alpha
self.perfect_fit = False
diff --git a/src/bayesvalidrox/surrogate_models/meta_model_engine.py b/src/bayesvalidrox/surrogate_models/meta_model_engine.py
index 2df2dee5390ae4e6dc7eb88343c2469dbd88aad6..c2e9ec409a7d17ceb697c694d9845ede185a0d0c 100644
--- a/src/bayesvalidrox/surrogate_models/meta_model_engine.py
+++ b/src/bayesvalidrox/surrogate_models/meta_model_engine.py
@@ -1067,7 +1067,7 @@ class MetaModelEngine():
maxfun=max_func_itr)
if verbose:
- print(f"global minimum: xmin = {Res_Global.x}, "
+ print(f"Global minimum: xmin = {Res_Global.x}, "
f"f(xmin) = {Res_Global.fun:.6f}, nfev = {Res_Global.nfev}")
return (Run_No, Res_Global.x)
@@ -1237,7 +1237,7 @@ class MetaModelEngine():
elapsed_time = time.time() - start_time
print("\n")
- print(f"elapsed_time: {round(elapsed_time,2)} sec.")
+ print(f"Elapsed_time: {round(elapsed_time,2)} sec.")
print('-'*20)
elif explore_method == 'LOOCV':
diff --git a/src/bayesvalidrox/surrogate_models/reg_fast_ard.py b/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
index 44073da8e78642ba3b3914f6ce55a2d01986b1f1..e6883a3edd6d247c219b8be328f5206b75780fbb 100755
--- a/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
+++ b/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
@@ -286,9 +286,9 @@ class RegressionFastARD(LinearModel, RegressorMixin):
# raise warning in case cholesky fails
if warning_flag == 1:
warnings.warn(("Cholesky decomposition failed! Algorithm uses "
- "pinvh, which is significantly slower, if you "
+ "pinvh, which is significantly slower. If you "
"use RVR it is advised to change parameters of "
- "kernel"))
+ "the kernel!"))
# compute quality & sparsity parameters
s, q, S, Q = self._sparsity_quality(XX, XXd, XY, XYa, Aa, Ri,
@@ -319,7 +319,7 @@ class RegressionFastARD(LinearModel, RegressorMixin):
if converged or i == self.n_iter - 1:
if converged and self.verbose:
- print('Algorithm converged !')
+ print('Algorithm converged!')
break
# after last update of alpha & beta update parameters
diff --git a/src/bayesvalidrox/surrogate_models/reg_fast_laplace.py b/src/bayesvalidrox/surrogate_models/reg_fast_laplace.py
index bdff324ede818a42d226e9aa55aaf01666ca8fc8..84faa6026e18f9013ed04b1bb08e2740671d3af3 100644
--- a/src/bayesvalidrox/surrogate_models/reg_fast_laplace.py
+++ b/src/bayesvalidrox/surrogate_models/reg_fast_laplace.py
@@ -307,7 +307,7 @@ class RegressionFastLaplace():
# recomputation
# zero if regressor has not been chosen yet
if not ind_global_to_local[ind_L_max]:
- raise Exception('cannot recompute index{0} -- not yet\
+ raise Exception('Cannot recompute index{0} -- not yet\
part of the model!'.format(ind_L_max))
Sigma = np.atleast_2d(Sigma)
mu = np.atleast_2d(mu)
diff --git a/src/bayesvalidrox/surrogate_models/sequential_design.py b/src/bayesvalidrox/surrogate_models/sequential_design.py
index fc81dcd4529ca0708dfba47385aef4415992eb3e..1e0935deed32ad5c398a82879d65723b22f03080 100644
--- a/src/bayesvalidrox/surrogate_models/sequential_design.py
+++ b/src/bayesvalidrox/surrogate_models/sequential_design.py
@@ -1099,7 +1099,7 @@ class SeqDesign():
maxfun=max_func_itr)
if verbose:
- print(f"global minimum: xmin = {Res_Global.x}, "
+ print(f"Global minimum: xmin = {Res_Global.x}, "
f"f(xmin) = {Res_Global.fun:.6f}, nfev = {Res_Global.nfev}")
return (Run_No, Res_Global.x)
@@ -1269,7 +1269,7 @@ class SeqDesign():
elapsed_time = time.time() - start_time
print("\n")
- print(f"elapsed_time: {round(elapsed_time,2)} sec.")
+ print(f"Elapsed_time: {round(elapsed_time,2)} sec.")
print('-'*20)
elif explore_method == 'LOOCV':
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index c51ea8c6679aaeec685da069154e8f460c7c4450..96cb1d8c54a40911951b70b76b09ad15a4d6908a 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -136,7 +136,7 @@ class MetaModel():
if self.ExpDesign.method == 'sequential':
raise Exception(
- "Please use MetaModelEngine class for the sequential design!"
+ "Please use the MetaModelEngine class for sequential design!"
)
elif self.ExpDesign.method == 'normal':
@@ -238,7 +238,7 @@ class MetaModel():
if verbose and self.n_bootstrap_itrs > 1:
enum_obj = tqdm(range(self.n_bootstrap_itrs),
total=self.n_bootstrap_itrs,
- desc="Boostraping the metamodel",
+ desc="Bootstrapping the metamodel",
ascii=True)
else:
enum_obj = range(self.n_bootstrap_itrs)
@@ -516,7 +516,7 @@ class MetaModel():
self.ModelOutputDict = ExpDesign.Y
else:
raise Exception('Please provide either a dictionary or a hdf5'
- 'file to ExpDesign.hdf5_file argument.')
+ 'file as the ExpDesign.hdf5_file argument.')
return ED_X_tr, self.ModelOutputDict
@@ -882,9 +882,9 @@ class MetaModel():
print(' '*10 + ' Summary of results ')
print('='*50)
- print("scores:\n", scores)
- print("Best score's degree:", self.deg_array[best_deg-1])
- print("NO. of terms:", len(basis_indices))
+ print("Scores:\n", scores)
+ print("Degree of best score:", self.deg_array[best_deg-1])
+ print("No. of terms:", len(basis_indices))
print("Sparsity index:", round(len(basis_indices)/P, 3))
print("Best Indices:\n", basis_indices)
@@ -980,7 +980,7 @@ class MetaModel():
# only the trace is, to save memory)
PsiM = np.dot(psi_sparse, M)
- h = np.sum(np.multiply(PsiM, psi_sparse), axis=1, dtype=np.float128)
+ h = np.sum(np.multiply(PsiM, psi_sparse), axis=1, dtype=np.longdouble)#float128)
# ------ Calculate Error Loocv for each measurement point ----
# Residuals
diff --git a/src/test_umbridge.py b/src/test_umbridge.py
new file mode 100644
index 0000000000000000000000000000000000000000..8e7c32cf4c868445d641ee834ca850c1b279f4da
--- /dev/null
+++ b/src/test_umbridge.py
@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Dec 11 14:03:23 2023
+
+Example and test for using umbridge with bayesvalidrox.
+Example model is the tsunami model by
+docker run -it -p 4242:4242 -v ~/tsunami_output:/output linusseelinger/model-exahype-tsunami
+
+@author: rkohlhaas
+"""
+
+import joblib
+import numpy as np
+import umbridge
+
+import bayesvalidrox as bv
+
+
+if __name__ == '__main__':
+
+ # This model has 2 inputs and four outputs
+ n_prior_sample = 1000
+ priors = bv.Input()
+ priors.add_marginals()
+ priors.Marginals[0].name = 'x'
+ priors.Marginals[0].input_data = np.random.uniform(50,150,n_prior_sample)
+ priors.add_marginals()
+ priors.Marginals[1].name = 'y'
+ priors.Marginals[1].input_data = np.random.uniform(50,150,n_prior_sample)
+
+ # Define the model
+ model = bv.PyLinkForwardModel()
+ #model.link_type = 'Function'
+ model.link_type = 'umbridge'
+ model.py_file = 'tsunami_model'
+ model.name = 'tsunami_model'
+ model.Output.names = ['T1', 'T2', 'H1', 'H2']
+ #model.observations = data_dict1_times
+
+ # Create the surrogate
+ surrogate_opts1 = bv.MetaModel(priors, model)
+
+ # Select your metamodel method
+ surrogate_opts1.meta_model_type = 'aPCE'
+ surrogate_opts1.pce_reg_method = 'FastARD'
+ surrogate_opts1.pce_deg = np.arange(1, 5)
+ surrogate_opts1.pce_q_norm = 0.4#1.0
+
+ # Select your metamodel method
+ surrogate_opts1.add_ExpDesign()
+ surrogate_opts1.ExpDesign.method = 'normal'
+ surrogate_opts1.ExpDesign.n_init_samples = 20
+ surrogate_opts1.ExpDesign.sampling_method = 'latin-hypercube'
+ surrogate1 = surrogate_opts1.create_metamodel()
+ print('Surrogate completed')
+ print('')
+
+ # Save surrogate
+ with open('tsunami.pk1', 'wb') as output:
+ joblib.dump(surrogate1, output, 2)
+
+ # Post processing
+ L2_PostPCE = bv.PostProcessing(surrogate1)
+ L2_PostPCE.plot_moments(plot_type='line')
+ # Plot to check validation visually.
+ L2_PostPCE.valid_metamodel(n_samples=1)
+ # Compute and print RMSE error
+ L2_PostPCE.check_accuracy(n_samples=300)
+ total_sobol = L2_PostPCE.sobol_indices()
+
+ # Test BayesInference
+ import copy
+ from tsunami_model import tsunami_model
+ true_data = tsunami_model(np.array([[90.0,60.0]]))
+ model.observations = true_data
+ true_data_nox = copy.deepcopy(true_data)
+ true_data_nox.pop('x_values')
+
+ # Test surrogate output shape
+ mean, std = surrogate1.eval_metamodel(np.array([[90.0,60.0]]))
+
+ # Set uncertainty
+ import pandas as pd
+ obsData = pd.DataFrame(true_data_nox, columns=model.Output.names)
+ DiscrepancyOpts = bv.Discrepancy('')
+ DiscrepancyOpts.type = 'Gaussian'
+ DiscrepancyOpts.parameters = (obsData*0.15)**2
+
+ # Parameter estimation / single model validation via TOM
+ BayesOpts = bv.BayesInference(surrogate1)
+ BayesOpts.emulator= True
+ BayesOpts.plot_post_pred = True
+ #BayesOpts.inference_method = 'rejection'
+ BayesOpts.Discrepancy = DiscrepancyOpts
+ BayesOpts.bootstrap = True
+ BayesOpts.n_bootstrap_itrs = 10
+ BayesOpts.bootstrap_noise = 0.05
+ BayesOpts.out_dir = ''
+
+ # Set the MCMC parameters
+ import emcee
+ BayesOpts.inference_method = "MCMC"
+ BayesOpts.mcmc_params = {
+ 'n_steps': 1e3,
+ 'n_walkers': 30,
+ 'moves': emcee.moves.KDEMove(),
+ 'multiprocessing': False,
+ 'verbose': False
+ }
+
+ Bayes = BayesOpts.create_inference()
+ # Pkl the inference results
+ with open(f'Bayes_{model.name}.pkl', 'wb') as output:
+ joblib.dump(Bayes, output, 2)
\ No newline at end of file
diff --git a/src/tsunami.pk1 b/src/tsunami.pk1
new file mode 100644
index 0000000000000000000000000000000000000000..3f5e6991ba13f525c9e5d63d9e044796c36de9e0
Binary files /dev/null and b/src/tsunami.pk1 differ
diff --git a/src/tsunami_model.py b/src/tsunami_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..6e9b0e50bed2ea2d0e01eb9125ff6c2d8cb12a0d
--- /dev/null
+++ b/src/tsunami_model.py
@@ -0,0 +1,16 @@
+# -*- coding: utf-8 -*-
+"""
+Runs the umbridge command for the tsunami model
+
+@author: Rebecca Kohlhaas
+"""
+import umbridge
+import numpy as np
+
+def tsunami_model(params):
+ model = umbridge.HTTPModel('http://localhost:4242', 'forward')
+ out = np.array(model(np.ndarray.tolist(params)))
+
+ return {'T1': out[:,0], 'T2': out[:,1], 'H1': out[:,2], 'H2': out[:,3], 'x_values': [0]}
+
+
\ No newline at end of file