Merge pull request #190 from SURGroup/Development

Development

README.rst

@@ -1,7 +1,6 @@
-|AzureDevops| |PyPIdownloads| |PyPI| |CondaSURG| |CondaPlatforms| |GithubRelease| |Binder| |Docs| |bear-ified|
+|AzureDevops| |PyPIdownloads| |PyPI| |CondaPlatforms| |GithubRelease| |Binder| |Docs| |bear-ified|
 
 .. |Docs| image:: https://img.shields.io/readthedocs/uqpy?style=plastic  :alt: Read the Docs
-.. |CondaSURG| image:: https://img.shields.io/conda/vn/SURG_JHU/uqpy?style=plastic   :alt: Conda (channel only)
 .. |CondaPlatforms| image:: https://img.shields.io/conda/pn/SURG_JHU/uqpy?style=plastic   :alt: Conda
 .. |GithubRelease| image:: https://img.shields.io/github/v/release/SURGroup/UQpy?style=plastic   :alt: GitHub release (latest by date)
 .. |AzureDevops| image:: https://img.shields.io/azure-devops/build/UQpy/5ce1851f-e51f-4e18-9eca-91c3ad9f9900/1?style=plastic   :alt: Azure DevOps builds
@@ -35,7 +34,7 @@ Uncertainty Quantification with python (UQpy)
 +                       +                                                                  +
 |                       | Promit Chakroborty, Lukáš Novák, Andrew Solanto                  |
 +-----------------------+------------------------------------------------------------------+
-| **Contributors:**     | Michael Gardner                                                  |
+| **Contributors:**     | Michael Gardner, Prateek Bhustali, Julius Schultz, Ulrich Römer  |
 +-----------------------+------------------------------------------------------------------+
 
 Contact
@@ -91,7 +90,6 @@ Using Conda
             * ::
 
                         conda install -c conda-forge uqpy
-                        conda install -c surg_jhu uqpy (latest version)
 
 Clone your fork of the UQpy repo from your GitHub account to your local disk (to get the latest version): 
 

azure-pipelines.yml

@@ -152,18 +152,18 @@ jobs:
       displayName: Install Anaconda packages
       condition: eq(variables['Build.SourceBranch'], 'refs/heads/master')
 
-    - bash: |
-        source activate myEnvironment
-        conda build . recipe --variants "{'version': ['$(GitVersion.SemVer)']}"
-      displayName: Build Noarch conda packages
-      condition: eq(variables['Build.SourceBranch'], 'refs/heads/master')
-
-    - bash: |
-        source activate myEnvironment
-        anaconda login --username $(ANACONDAUSER) --password $(ANACONDAPW)
-        anaconda upload /usr/local/miniconda/envs/myEnvironment/conda-bld/noarch/*.tar.bz2
-      displayName: Upload conda packages
-      condition: eq(variables['Build.SourceBranch'], 'refs/heads/master')
+#    - bash: |
+#        source activate myEnvironment
+#        conda build . recipe --variants "{'version': ['$(GitVersion.SemVer)']}"
+#      displayName: Build Noarch conda packages
+#      condition: eq(variables['Build.SourceBranch'], 'refs/heads/master')
+#
+#    - bash: |
+#        source activate myEnvironment
+#        anaconda login --username $(ANACONDAUSER) --password $(ANACONDAPW)
+#        anaconda upload /usr/local/miniconda/envs/myEnvironment/conda-bld/noarch/*.tar.bz2
+#      displayName: Upload conda packages
+#      condition: eq(variables['Build.SourceBranch'], 'refs/heads/master')
 
 - job: "Create_Docker_images"
   dependsOn: Build_UQpy_and_run_tests

docs/code/inference/bayes_model_selection/bayes_model_selection.py

@@ -46,8 +46,9 @@
 error_covariance = var_n * np.eye(50)
 print(param_true.shape)
 
-z = RunModel(samples=param_true, model_script='local_pfn_models.py', model_object_name='model_quadratic', vec=False,
-             var_names=['theta_1', 'theta_2'])
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+m=PythonModel(model_script='local_pfn_models.py', model_object_name='model_quadratic', var_names=['theta_1', 'theta_2'])
+z = RunModel(samples=param_true, model=m)
 data_clean = z.qoi_list[0].reshape((-1,))
 data = data_clean + Normal(scale=np.sqrt(var_n)).rvs(nsamples=data_clean.size, random_state=456).reshape((-1,))
 
@@ -66,55 +67,10 @@
 model_prior_stds = [[10.], [1., 1.], [1., 2., 0.25]]
 
 
-evidences = []
-model_posterior_means = []
-model_posterior_stds = []
-for n, model in enumerate(model_names):
-    # compute matrix X
-    X = np.linspace(0, 10, 50).reshape((-1, 1))
-    if n == 1:  # quadratic model
-        X = np.concatenate([X, X ** 2], axis=1)
-    if n == 2:  # cubic model
-        X = np.concatenate([X, X ** 2, X ** 3], axis=1)
-
-    # compute posterior pdf
-    m_prior = np.array(model_prior_means[n]).reshape((-1, 1))
-    S_prior = np.diag(np.array(model_prior_stds[n]) ** 2)
-    S_posterior = np.linalg.inv(1 / var_n * np.matmul(X.T, X) + np.linalg.inv(S_prior))
-    m_posterior = np.matmul(S_posterior,
-                            1 / var_n * np.matmul(X.T, data.reshape((-1, 1))) + np.matmul(np.linalg.inv(S_prior),
-                                                                                          m_prior))
-    m_prior = m_prior.reshape((-1,))
-    m_posterior = m_posterior.reshape((-1,))
-    model_posterior_means.append(list(m_posterior))
-    model_posterior_stds.append(list(np.sqrt(np.diag(S_posterior))))
-    print('posterior mean and covariance for ' + model)
-    print(m_posterior, S_posterior)
-
-    # compute evidence, evaluate the formula at the posterior mean
-    like_theta = multivariate_normal.pdf(data, mean=np.matmul(X, m_posterior).reshape((-1,)), cov=error_covariance)
-    prior_theta = multivariate_normal.pdf(m_posterior, mean=m_prior, cov=S_prior)
-    posterior_theta = multivariate_normal.pdf(m_posterior, mean=m_posterior, cov=S_posterior)
-    evidence = like_theta * prior_theta / posterior_theta
-    evidences.append(evidence)
-    print('evidence for ' + model + '= {}\n'.format(evidence))
-
-# compute the posterior probability of each model
-tmp = [1 / 3 * evidence for evidence in evidences]
-model_posterior_probas = [p / sum(tmp) for p in tmp]
-
-print('posterior probabilities of all three models')
-print(model_posterior_probas)
-
-#%% md
-#
-# Define the models for use in UQpy
-
-#%%
-
 candidate_models = []
 for n, model_name in enumerate(model_names):
-    run_model = RunModel(model_script='local_pfn_models.py', model_object_name=model_name, vec=False)
+    m=PythonModel(model_script='local_pfn_models.py', model_object_name=model_name,)
+    run_model = RunModel(model=m)
     prior = JointIndependent([Normal(loc=m, scale=std) for m, std in
                               zip(model_prior_means[n], model_prior_stds[n])])
     model = ComputationalModel(n_parameters=model_n_params[n],
@@ -123,60 +79,18 @@
                                name=model_name)
     candidate_models.append(model)
 
-
-#%% md
-#
-# Run MCMC for one model
-
-#%%
-
-# Quadratic model
-sampling = MetropolisHastings(args_target=(data, ),
-                              log_pdf_target=candidate_models[1].evaluate_log_posterior,
-                              jump=10, burn_length=100,
-                              proposal=JointIndependent([Normal(scale=0.1), ] * 2),
-                              seed=[0., 0.], random_state=123)
-
-bayesMCMC = BayesParameterEstimation(sampling_class=sampling,
-                                     inference_model=candidate_models[1],
-                                     data=data,
-                                     nsamples=3500)
-
-# plot prior, true posterior and estimated posterior
-fig, ax = plt.subplots(1, 2, figsize=(16, 5))
-for n_p in range(2):
-    domain_plot = np.linspace(-0.5, 3, 200)
-    ax[n_p].plot(domain_plot, norm.pdf(domain_plot, loc=model_prior_means[1][n_p], scale=model_prior_stds[1][n_p]),
-                 label='prior', color='green', linestyle='--')
-    ax[n_p].plot(domain_plot, norm.pdf(domain_plot, loc=model_posterior_means[1][n_p],
-                                       scale=model_posterior_stds[1][n_p]),
-                 label='true posterior', color='red', linestyle='-')
-    ax[n_p].hist(bayesMCMC.sampler.samples[:, n_p], density=True, bins=30, label='estimated posterior MCMC')
-    ax[n_p].legend()
-    ax[n_p].set_title('MCMC for quadratic model')
-plt.show()
-
-
-#%% md
-#
-# Run Bayesian Model Selection for all three models
-
-#%%
-
 proposals = [Normal(scale=0.1),
              JointIndependent([Normal(scale=0.1), Normal(scale=0.1)]),
              JointIndependent([Normal(scale=0.15), Normal(scale=0.1), Normal(scale=0.05)])]
-nsamples = [2000, 6000, 14000]
-nburn = [500, 2000, 4000]
-jump = [5, 10, 25]
+nsamples = [2000, 2000, 2000]
+nburn = [1000, 1000, 1000]
+jump = [2, 2, 2]
 
 
 sampling_inputs=list()
 estimators = []
 for i in range(3):
-    sampling = MetropolisHastings(args_target=(data, ),
-                                  log_pdf_target=candidate_models[i].evaluate_log_posterior,
-                                  jump=jump[i],
+    sampling = MetropolisHastings(jump=jump[i],
                                   burn_length=nburn[i],
                                   proposal=proposals[i],
                                   seed=model_prior_means[i],
@@ -185,7 +99,6 @@
                                                   sampling_class=sampling))
 
 selection = BayesModelSelection(parameter_estimators=estimators,
-                                data=data,
                                 prior_probabilities=[1. / 3., 1. / 3., 1. / 3.],
                                 nsamples=nsamples)
 

docs/code/sampling/mcmc/plot_mcmc_metropolis_hastings.py

@@ -58,7 +58,7 @@ def log_rosenbrock_with_param(x, p):
 plt.show()
 
 plt.figure()
-x = MetropolisHastings(dimension=2, pdf_target=log_rosenbrock_with_param, burn_length=500,
+x = MetropolisHastings(dimension=2, log_pdf_target=log_rosenbrock_with_param, burn_length=500,
                        jump=50, n_chains=1, args_target=(20,),
                        nsamples=500)
 plt.plot(x.samples[:, 0], x.samples[:, 1], 'o')

docs/code/sensitivity/chatterjee/README.rst

@@ -0,0 +1,15 @@
+Chatterjee indices
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+These examples serve as a guide for using the Chatterjee sensitivity module. They have been taken from various papers to enable validation of the implementation and have been referenced accordingly.
+
+1. **Ishigami function**
+
+    In addition to the Pick and Freeze scheme, the Sobol indices can be estimated using the rank statistics approach :cite:`gamboa2020global`. We demonstrate this estimation of the Sobol indices using the Ishigami function.
+
+2. **Exponential function**
+
+    For the Exponential model, analytical Cramér-von Mises indices are available :cite:`CVM` and since they are equivalent to the Chatterjee indices in the sample limit, they are shown here.
+
+3. **Sobol function**
+
+    This example was considered in :cite:`gamboa2020global` (page 18) to compare the Pick and Freeze scheme with the rank statistics approach for estimating the Sobol indices.

docs/code/sensitivity/chatterjee/chatterjee_exponential.py

@@ -0,0 +1,76 @@
+"""
+
+Exponential function
+==============================================
+
+The exponential function was used in [1]_ to demonstrate the 
+Cramér-von Mises indices. Chattererjee indices approach the Cramér-von Mises 
+indices in the sample limit and will be demonstrated via this example.
+
+.. math::
+    f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
+
+.. [1] Gamboa, F., Klein, T., & Lagnoux, A. (2018). Sensitivity Analysis Based on \
+Cramér-von Mises Distance. SIAM/ASA Journal on Uncertainty Quantification, 6(2), \
+522-548. doi:10.1137/15M1025621. (`Link <https://doi.org/10.1137/15M1025621>`_)
+
+"""
+
+# %%
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Normal
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.ChatterjeeSensitivity import ChatterjeeSensitivity
+from UQpy.sensitivity.PostProcess import *
+
+np.random.seed(123)
+
+# %% [markdown]
+# **Define the model and input distributions**
+
+# Create Model object
+model = PythonModel(
+    model_script="local_exponential.py",
+    model_object_name="evaluate",
+    var_names=[
+        "X_1",
+        "X_2",
+    ],
+    delete_files=True,
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Normal(0, 1)] * 2)
+
+# %% [markdown]
+# **Compute Chatterjee indices**
+
+# %% [markdown]
+SA = ChatterjeeSensitivity(runmodel_obj, dist_object)
+
+# Compute Chatterjee indices using the pick and freeze algorithm
+SA.run(n_samples=1_000_000)
+
+# %% [markdown]
+# **Chattererjee indices**
+#
+# Chattererjee indices approach the Cramér-von Mises indices in the sample limit.
+#
+# Expected value of the sensitivity indices:
+#
+# :math:`S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145`
+#
+# :math:`S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693`
+
+# %%
+SA.first_order_chatterjee_indices
+
+# **Plot the Chatterjee indices**
+fig1, ax1 = plot_sensitivity_index(
+    SA.first_order_chatterjee_indices[:, 0],
+    plot_title="Chatterjee indices",
+    color="C2",
+)