From d54253d5732e5581a52fbe644e0a36f3c577fd99 Mon Sep 17 00:00:00 2001
From: Dimitris Tsapetis <dimtsap@hotmail.com>
Date: Wed, 6 Apr 2022 13:49:55 -0400
Subject: [PATCH] Replaces RunModel in examples

---
 .../plot_bayes_parameter_IS_regression.py     |  7 +++---
 .../plot_bayes_parameter_MCMC_regression.py   |  7 +++---
 .../plot_selection_regression_model.py        | 12 +++++-----
 .../inference/mle/plot_regression_model.py    |  9 +++----
 .../plot_subset_resonance.py                  | 13 +++++-----
 .../plot_subset_rosenbrock.py                 | 24 +++++++++----------
 .../plot_adaptive_kriging_branin_hoo.py       |  5 ++--
 .../plot_adaptive_kriging_normal.py           | 14 +++++------
 docs/code/surrogates/gpr/plot_gpr_custom2D.py | 11 +++++----
 docs/code/surrogates/gpr/plot_gpr_sine.py     |  9 +++----
 docs/code/transformations/nataf/plot_nataf.py |  4 ++--
 11 files changed, 61 insertions(+), 54 deletions(-)

diff --git a/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_IS_regression.py b/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_IS_regression.py
index d03b06c34..a659865d6 100644
--- a/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_IS_regression.py
+++ b/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_IS_regression.py
@@ -22,6 +22,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
+from UQpy import PythonModel
 from UQpy.sampling.ImportanceSampling import ImportanceSampling
 from UQpy.inference import BayesParameterEstimation, ComputationalModel
 from UQpy.run_model.RunModel_New import RunModel_New  # required to run the quadratic model
@@ -43,8 +44,9 @@ def pdf_from_kde(domain, samples1d):
 param_true = np.array([1.0, 2.0]).reshape((1, -1))
 print('Shape of true parameter vector: {}'.format(param_true.shape))
 
-h_func = RunModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', vec=False,
-                  var_names=['theta_0', 'theta_1'])
+model = PythonModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', delete_files=True,
+                    var_names=['theta_0', 'theta_1'])
+h_func = RunModel_New(model=model)
 h_func.run(samples=param_true)
 
 # Add noise
@@ -80,4 +82,3 @@ def pdf_from_kde(domain, samples1d):
 
 plt.show()
 
-shutil.rmtree(h_func.model_dir)
diff --git a/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_MCMC_regression.py b/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_MCMC_regression.py
index fcbb2285d..aad83b44c 100644
--- a/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_MCMC_regression.py
+++ b/docs/code/inference/bayes_parameter_estimation/plot_bayes_parameter_MCMC_regression.py
@@ -22,6 +22,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
+from UQpy import PythonModel
 from UQpy.sampling.mcmc.MetropolisHastings import MetropolisHastings
 from UQpy.inference.inference_models.ComputationalModel import ComputationalModel
 from UQpy.run_model.RunModel_New import RunModel_New
@@ -45,8 +46,9 @@ def pdf_from_kde(domain, samples1d):
 param_true = np.array([1.0, 2.0]).reshape((1, -1))
 print('Shape of true parameter vector: {}'.format(param_true.shape))
 
-h_func = RunModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', vec=False,
-                  var_names=['theta_0', 'theta_1'])
+model = PythonModel(model_script='local_pfn_models.py', model_object_name='model_quadratic',
+                    var_names=['theta_0', 'theta_1'])
+h_func = RunModel_New(model=model)
 h_func.run(samples=param_true)
 data_clean = np.array(h_func.qoi_list[0])
 print(data_clean.shape)
@@ -99,4 +101,3 @@ def pdf_from_kde(domain, samples1d):
 
 print(bayes_estimator.sampler.samples[:4])
 
-shutil.rmtree(h_func.model_dir)
diff --git a/docs/code/inference/info_model_selection/plot_selection_regression_model.py b/docs/code/inference/info_model_selection/plot_selection_regression_model.py
index c7350543d..28d937e31 100644
--- a/docs/code/inference/info_model_selection/plot_selection_regression_model.py
+++ b/docs/code/inference/info_model_selection/plot_selection_regression_model.py
@@ -25,6 +25,7 @@
 #%%
 import shutil
 
+from UQpy import PythonModel
 from UQpy.inference import InformationModelSelection, MLE
 from UQpy.run_model.RunModel_New import RunModel_New
 import numpy as np
@@ -43,8 +44,8 @@
 param_true = np.array([1.0, 2.0]).reshape((1, -1))
 print('Shape of true parameter vector: {}'.format(param_true.shape))
 
-h_func = RunModel(model_script='pfn_models.py', model_object_name='model_quadratic', vec=False,
-                  var_names=['theta_0', 'theta_1'])
+model = PythonModel(model_script='pfn_models.py', model_object_name='model_quadratic', var_names=['theta_0', 'theta_1'])
+h_func = RunModel_New(model=model)
 h_func.run(samples=param_true)
 
 # Add noise
@@ -54,7 +55,6 @@
 data_1 = data_clean + noise
 print('Shape of data: {}'.format(data_1.shape))
 
-shutil.rmtree(h_func.model_dir)
 
 #%% md
 #
@@ -66,8 +66,9 @@
 estimators = []
 
 for i in range(3):
-    h_func = RunModel(model_script='pfn_models.py', model_object_name=names[i], vec=False,
-                      var_names=['theta_{}'.format(j) for j in range(i + 1)])
+    model = PythonModel(model_script='pfn_models.py', model_object_name=names[i],
+                        var_names=['theta_{}'.format(j) for j in range(i + 1)])
+    h_func = RunModel_New(model=model)
     M = ComputationalModel(runmodel_object=h_func, n_parameters=i + 1,
                            name=names[i], error_covariance=error_covariance)
     estimators.append(MLE(inference_model=M, data=data_1))
@@ -104,7 +105,6 @@
     model.runmodel_object.run(samples=estim.mle.reshape((1, -1)), append_samples=False)
     y = model.runmodel_object.qoi_list[-1].reshape((-1,))
     ax.plot(domain, y, label=selector.candidate_models[i].name)
-    shutil.rmtree(model.runmodel_object.model_dir)
 
 plt.plot(domain, data_1, linestyle='none', marker='.', label='data')
 plt.xlabel('x')
diff --git a/docs/code/inference/mle/plot_regression_model.py b/docs/code/inference/mle/plot_regression_model.py
index 5b513c499..fdfa1a77d 100644
--- a/docs/code/inference/mle/plot_regression_model.py
+++ b/docs/code/inference/mle/plot_regression_model.py
@@ -22,6 +22,8 @@
 import shutil
 
 import numpy as np
+
+from UQpy import PythonModel
 from UQpy.inference import ComputationalModel, MLE
 from UQpy.distributions import Normal
 from UQpy.inference import MinimizeOptimizer
@@ -38,8 +40,9 @@
 param_true = np.array([1.0, 2.0]).reshape((1, -1))
 print('Shape of true parameter vector: {}'.format(param_true.shape))
 
-h_func = RunModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', vec=False,
-                  var_names=['theta_0', 'theta_1'])
+model = PythonModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', delete_files=True,
+                    var_names=['theta_0', 'theta_1'])
+h_func = RunModel_New(model=model)
 h_func.run(samples=param_true)
 
 # Add noise
@@ -63,5 +66,3 @@
 ml_estimator = MLE(inference_model=candidate_model, data=data_3, n_optimizations=1)
 print('fitted parameters: theta_0={0:.3f} (true=1.), and theta_1={1:.3f} (true=2.)'.format(
     ml_estimator.mle[0], ml_estimator.mle[1]))
-
-shutil.rmtree(h_func.model_dir)
\ No newline at end of file
diff --git a/docs/code/reliability/subset_simulation/plot_subset_resonance.py b/docs/code/reliability/subset_simulation/plot_subset_resonance.py
index 742356704..46ec655a4 100644
--- a/docs/code/reliability/subset_simulation/plot_subset_resonance.py
+++ b/docs/code/reliability/subset_simulation/plot_subset_resonance.py
@@ -29,6 +29,7 @@
 # %%
 import shutil
 
+from UQpy import PythonModel
 from UQpy.reliability import SubsetSimulation
 from UQpy.run_model.RunModel_New import RunModel_New
 from UQpy.sampling import Stretch, ModifiedMetropolisHastings, MonteCarloSampling
@@ -111,7 +112,8 @@
 plt.tight_layout()
 plt.show()
 
-model = RunModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel", ntasks=1)
+m = PythonModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel")
+model = RunModel_New(model=m)
 
 # %% md
 #
@@ -124,7 +126,6 @@
 
 model.run(samples=x_mcs.samples)
 
-shutil.rmtree(model.model_dir)
 A = np.asarray(model.qoi_list) < 0
 pf = np.shape(np.asarray(model.qoi_list)[np.asarray(model.qoi_list) < 0])[0] / 1000000
 print(pf)
@@ -141,7 +142,8 @@
 C[1, 1] = 20 ** 2
 
 for i in range(ntrials):
-    model = RunModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel", ntasks=1)
+    m1 = PythonModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel")
+    model = RunModel_New(model=m1)
     dist = MultivariateNormal(mean=m, cov=C)
     xx = dist.rvs(nsamples=1000, random_state=123)
     xx1 = dist.rvs(nsamples=100, random_state=123)
@@ -149,7 +151,6 @@
     sampling=Stretch(dimension=2, n_chains=100, log_pdf_target=dist.log_pdf)
     x_ss_stretch = SubsetSimulation(sampling=sampling, runmodel_object=model, conditional_probability=0.1,
                                     nsamples_per_subset=1000, samples_init=xx, )
-    shutil.rmtree(model.model_dir)
     pf_stretch[i] = x_ss_stretch.failure_probability
     cov1_stretch[i] = x_ss_stretch.independent_chains_CoV
     cov2_stretch[i] = x_ss_stretch.dependent_chains_CoV
@@ -185,7 +186,8 @@
 C[1, 1] = 20 ** 2
 
 for i in range(ntrials):
-    model = RunModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel", ntasks=1)
+    m1 = PythonModel(model_script='local_Resonance_pfn.py', model_object_name="RunPythonModel")
+    model = RunModel_New(model=m1)
     dist = MultivariateNormal(mean=m, cov=C)
     xx = dist.rvs(nsamples=1000, random_state=123)
     xx1 = dist.rvs(nsamples=100, random_state=123)
@@ -194,7 +196,6 @@
     x_ss_mmh = SubsetSimulation(sampling=sampling, runmodel_object=model, conditional_probability=0.1,
                                 nsamples_per_subset=1000, samples_init=xx)
 
-    shutil.rmtree(model.model_dir)
     pf_mmh[i] = x_ss_mmh.failure_probability
     cov1_mmh[i] = x_ss_mmh.independent_chains_CoV
     cov2_mmh[i] = x_ss_mmh.dependent_chains_CoV
diff --git a/docs/code/reliability/subset_simulation/plot_subset_rosenbrock.py b/docs/code/reliability/subset_simulation/plot_subset_rosenbrock.py
index 08c1568a1..5fdbcbeae 100644
--- a/docs/code/reliability/subset_simulation/plot_subset_rosenbrock.py
+++ b/docs/code/reliability/subset_simulation/plot_subset_rosenbrock.py
@@ -12,18 +12,18 @@
 
 # %%
 
-import shutil
-
-from UQpy.reliability import SubsetSimulation
 import matplotlib.pyplot as plt
-from UQpy.sampling import ModifiedMetropolisHastings, Stretch
-# First import the file that contains the newly defined Rosenbrock distribution
-from local_Rosenbrock import Rosenbrock
 import numpy as np
+import scipy.stats as stats
+
+from UQpy import PythonModel
 # Import this newly defined Rosenbrock distribution into the Distributions module
 from UQpy.distributions import Normal
+from UQpy.reliability import SubsetSimulation
 from UQpy.run_model.RunModel_New import RunModel_New
-import scipy.stats as stats
+from UQpy.sampling import ModifiedMetropolisHastings, Stretch
+# First import the file that contains the newly defined Rosenbrock distribution
+from local_Rosenbrock import Rosenbrock
 
 # %% md
 #
@@ -31,7 +31,8 @@
 
 # %%
 
-model = RunModel(model_script='local_Rosenbrock_pfn.py', model_object_name="RunPythonModel", ntasks=1)
+m = PythonModel(model_script='local_Rosenbrock_pfn.py', model_object_name="RunPythonModel")
+model = RunModel_New(model=m)
 dist = Rosenbrock(p=100.)
 dist_prop1 = Normal(loc=0, scale=1)
 dist_prop2 = Normal(loc=0, scale=10)
@@ -57,7 +58,6 @@
 plt.tight_layout()
 plt.show()
 
-shutil.rmtree(model.model_dir)
 
 print(x_ss_MMH.failure_probability)
 
@@ -68,7 +68,8 @@
 
 # %%
 
-model = RunModel(model_script='local_Rosenbrock_pfn.py', model_object_name="RunPythonModel", ntasks=1)
+m = PythonModel(model_script='local_Rosenbrock_pfn.py', model_object_name="RunPythonModel")
+model = RunModel_New(model=m)
 dist = Rosenbrock(p=100.)
 
 x = stats.norm.rvs(loc=0, scale=1, size=(100, 2), random_state=83276)
@@ -103,5 +104,4 @@
 plt.yticks(np.arange(-20, 180, step=20))
 plt.xlim((-10, 15))
 plt.tight_layout()
-plt.show()
-shutil.rmtree(model.model_dir)
\ No newline at end of file
+plt.show()
\ No newline at end of file
diff --git a/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_branin_hoo.py b/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_branin_hoo.py
index 855de8060..710247463 100644
--- a/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_branin_hoo.py
+++ b/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_branin_hoo.py
@@ -40,6 +40,7 @@
 import numpy as np
 from matplotlib import pyplot as plt
 
+from UQpy import PythonModel
 from UQpy.surrogates import GaussianProcessRegression
 from UQpy.sampling import MonteCarloSampling, AdaptiveKriging
 from UQpy.run_model.RunModel_New import RunModel_New
@@ -64,7 +65,8 @@
 
 # %%
 
-rmodel = RunModel(model_script='local_BraninHoo.py', vec=False)
+model = PythonModel(model_script='local_BraninHoo.py', model_object_name='function')
+rmodel = RunModel_New(model=model)
 
 # %% md
 #
@@ -119,7 +121,6 @@
         tem = np.array([[X[i, j], Y[i, j]]])
         Z[i, j] = function(tem)
 
-shutil.rmtree(rmodel.model_dir)
 
 fig, ax = plt.subplots(1, 1)
 cp = ax.contourf(X, Y, Z, 10)
diff --git a/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_normal.py b/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_normal.py
index 66559afa2..495a6e3f7 100644
--- a/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_normal.py
+++ b/docs/code/sampling/adaptive_kriging/plot_adaptive_kriging_normal.py
@@ -16,6 +16,7 @@
 # %%
 import shutil
 
+from UQpy import PythonModel
 from UQpy.surrogates.gaussian_process import GaussianProcessRegression
 from UQpy.sampling import MonteCarloSampling, AdaptiveKriging
 from UQpy.run_model.RunModel_New import RunModel_New
@@ -42,7 +43,8 @@
 
 # %%
 
-rmodel = RunModel(model_script='local_series.py', vec=False)
+model = PythonModel(model_script='local_series.py', model_object_name='series')
+rmodel = RunModel_New(model=model)
 
 
 # %% md
@@ -153,7 +155,8 @@ def evaluate_function(self, distributions, n_add, surrogate, population, qoi=Non
 optimizer = MinimizeOptimizer(method="L-BFGS-B", bounds=bounds)
 K1 = GaussianProcessRegression(regression_model=LinearRegression(), kernel=RBF(), optimizer=optimizer,
                                hyperparameters=[1, 1, 0.1], optimizations_number=1)
-rmodel1 = RunModel(model_script='local_series.py', vec=False)
+model = PythonModel(model_script='local_series.py', model_object_name='series')
+rmodel1 = RunModel_New(model=model)
 
 # %% md
 #
@@ -205,7 +208,8 @@ def evaluate_function(self, distributions, n_add, surrogate, population, qoi=Non
 
 # Code
 b = MonteCarloSampling(distributions=marginals, nsamples=10 ** 4, random_state=4)
-r1model = RunModel(model_script='local_series.py', vec=False)
+model = PythonModel(model_script='local_series.py', model_object_name='series')
+r1model = RunModel_New(model=model)
 r1model.run(samples=b.samples)
 
 
@@ -215,10 +219,6 @@ def evaluate_function(self, distributions, n_add, surrogate, population, qoi=Non
 elapsed_time = time.time() - start_time
 time.strftime("%H:%M:%S", time.gmtime(elapsed_time))
 
-shutil.rmtree(r1model.model_dir)
-shutil.rmtree(rmodel1.model_dir)
-shutil.rmtree(rmodel.model_dir)
-
 # %% md
 #
 # Results from Monte Carlo Simulation.
diff --git a/docs/code/surrogates/gpr/plot_gpr_custom2D.py b/docs/code/surrogates/gpr/plot_gpr_custom2D.py
index f68fb681d..7ab09298d 100644
--- a/docs/code/surrogates/gpr/plot_gpr_custom2D.py
+++ b/docs/code/surrogates/gpr/plot_gpr_custom2D.py
@@ -22,6 +22,7 @@
 
 import shutil
 
+from UQpy import PythonModel
 from UQpy.surrogates.gaussian_process.regression_models import ConstantRegression
 from UQpy.sampling import RectangularStrata
 from UQpy.sampling import TrueStratifiedSampling
@@ -66,9 +67,10 @@
 
 # %%
 
-rmodel = RunModel(model_script='local_python_model_function.py', vec=False)
+model = PythonModel(model_script='local_python_model_function.py', model_object_name="y_func")
+rmodel = RunModel_New(model=model)
+
 rmodel.run(samples=x.samples)
-shutil.rmtree(rmodel.model_dir)
 
 # %% md
 #
@@ -104,12 +106,11 @@
 x1gv, x2gv = x1g.reshape(x1g.size, 1), x2g.reshape(x2g.size, 1)
 
 y2 = K.predict(np.concatenate([x1gv, x2gv], 1)).reshape(x1g.shape[0], x1g.shape[1])
-r2model = RunModel(model_script='local_python_model_function.py')
+model = PythonModel(model_script='local_python_model_function.py', model_object_name="y_func")
+r2model = RunModel_New(model=model)
 r2model.run(samples=np.concatenate([x1gv, x2gv], 1))
 y_act = np.array(r2model.qoi_list).reshape(x1g.shape[0], x1g.shape[1])
 
-shutil.rmtree(r2model.model_dir)
-
 fig1 = plt.figure()
 ax = fig1.gca(projection='3d')
 surf = ax.plot_surface(x1g, x2g, y_act, cmap=cm.coolwarm, linewidth=0, antialiased=False)
diff --git a/docs/code/surrogates/gpr/plot_gpr_sine.py b/docs/code/surrogates/gpr/plot_gpr_sine.py
index 601ebaaeb..7a9227e04 100644
--- a/docs/code/surrogates/gpr/plot_gpr_sine.py
+++ b/docs/code/surrogates/gpr/plot_gpr_sine.py
@@ -22,6 +22,7 @@
 
 import shutil
 
+from UQpy import PythonModel
 from UQpy.sampling.stratified_sampling.strata import RectangularStrata
 from UQpy.sampling import TrueStratifiedSampling
 from UQpy.run_model.RunModel_New import RunModel_New
@@ -63,8 +64,9 @@
 
 # %%
 
-
-rmodel = RunModel(model_script='local_python_model_1Dfunction.py', delete_files=True)
+model = PythonModel(model_script='local_python_model_1Dfunction.py', model_object_name='y_func',
+                    delete_files=True)
+rmodel = RunModel_New(model=model)
 rmodel.run(samples=x.samples)
 
 from UQpy.surrogates.gaussian_process.regression_models import LinearRegression
@@ -99,7 +101,6 @@
 # %%
 
 rmodel.run(samples=x1, append_samples=False)
-shutil.rmtree(rmodel.model_dir)
 
 # %% md
 #
@@ -110,7 +111,7 @@
 
 fig = plt.figure()
 ax = plt.subplot(111)
-plt.plot(x1, rmodel.qoi_list, label='Sine')
+plt.plot(np.squeeze(x1), np.squeeze(rmodel.qoi_list), label='Sine')
 plt.plot(x1, y, label='Surrogate')
 # plt.plot(x1, y_grad, label='Gradient')
 plt.scatter(K.samples, K.values, label='Data')
diff --git a/docs/code/transformations/nataf/plot_nataf.py b/docs/code/transformations/nataf/plot_nataf.py
index b1a2df368..8f5d0ecc6 100644
--- a/docs/code/transformations/nataf/plot_nataf.py
+++ b/docs/code/transformations/nataf/plot_nataf.py
@@ -158,8 +158,8 @@
     w3[i] = ww[0, 1]
 
 plt.plot(rho, w3)
-plt.xlabel('$\rho_X$')
-plt.ylabel('$\rho_Z$')
+plt.xlabel('rho_X')
+plt.ylabel('rho_Z')
 plt.show()