Repository: damar-wicaksono/uqtestfuns (language: python)

Issue:
Add the series of test functions from Genz (1984)
The paper by Genz (1984)[^1] contains a series of test functions used in numerical integration exercises. Some of the functions were often revisited in the context of sensitivity analysis (see, for instance, Zhang and Pandey (2014) [^2]).

Implement and include the functions in UQTestFuns:

- [x] Oscillatory (resolved via PR #416)
- [x] Product Peak (resolved via PR #413)
- [x] Corner Peak (resolved via PR #412)
- [x] Gaussian (resolved via PR #415)
- [x] Continuous (resolved via PR #414)
- [x] Discontinuous (resolved via PR #417)

[^1]: A. Genz, “Testing Multidimensional Integration Routines,” in Proc. of International Conference on Tools, Methods and Languages for Scientific and Engineering Computation, USA: Elsevier North-Holland, Inc., 1984, pp. 81–94.
[^2]: X. Zhang and M. D. Pandey, “An effective approximation for variance-based global sensitivity analysis,” Reliability Engineering & System Safety, vol. 121, pp. 164–174, Jan. 2014, doi: [10.1016/j.ress.2013.07.010](https://doi.org/10.1016/j.ress.2013.07.010).


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0b5a771..1816e3c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- The M-dimensional Gaussian function from Genz (1984) for integration
+  exercises; one parameter set from the literature is included.
 - The M-dimensional continuous function from Genz (1984) for integration
   exercises; one parameter set from the literature is included.
 - The M-dimensional product peak function from Genz (1984) for integration
diff --git a/docs/_toc.yml b/docs/_toc.yml
index b5ff945..558807f 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -98,6 +98,8 @@ parts:
             title: Genz (Continuous)
           - file: test-functions/genz-corner-peak
             title: Genz (Corner Peak)
+          - file: test-functions/genz-gaussian
+            title: Genz (Gaussian)
           - file: test-functions/genz-product-peak
             title: Genz (Product Peak)
           - file: test-functions/gramacy-1d-sine
diff --git a/docs/fundamentals/integration.md b/docs/fundamentals/integration.md
index 3cbbc0c..b938145 100644
--- a/docs/fundamentals/integration.md
+++ b/docs/fundamentals/integration.md
@@ -18,17 +18,18 @@ kernelspec:
 The table below listed the available test functions typically used
 in the testing and comparison of numerical integration method.
 
-|                                    Name                                     | Input Dimension |     Constructor     |
-|:---------------------------------------------------------------------------:|:---------------:|:-------------------:|
-|        {ref}`Bratley et al. (1992) A <test-functions:bratley1992a>`         |        M        |  `Bratley1992a()`   |
-|        {ref}`Bratley et al. (1992) B <test-functions:bratley1992b>`         |        M        |  `Bratley1992b()`   |
-|        {ref}`Bratley et al. (1992) C <test-functions:bratley1992c>`         |        M        |  `Bratley1992c()`   |
-|        {ref}`Bratley et al. (1992) D <test-functions:bratley1992d>`         |        M        |  `Bratley1992d()`   |
-|          {ref}`Genz (Continuous) <test-functions:genz-continuous>`          |        M        | `GenzContinuous()`  |
-|         {ref}`Genz (Corner Peak) <test-functions:genz-corner-peak>`         |        M        | `GenzCornerPeak()`  |
-|        {ref}`Genz (Product Peak) <test-functions:genz-product-peak>`        |        M        | `GenzProductPeak()` |
-|                  {ref}`Sobol'-G <test-functions:sobol-g>`                   |        M        |     `SobolG()`      |
-|            {ref}`Welch et al. (1992) <test-functions:welch1992>`            |       20        |    `Welch1992()`    |
+|                                   Name                                    | Input Dimension |     Constructor     |
+|:-------------------------------------------------------------------------:|:---------------:|:-------------------:|
+|       {ref}`Bratley et al. (1992) A <test-functions:bratley1992a>`        |        M        |  `Bratley1992a()`   |
+|       {ref}`Bratley et al. (1992) B <test-functions:bratley1992b>`        |        M        |  `Bratley1992b()`   |
+|       {ref}`Bratley et al. (1992) C <test-functions:bratley1992c>`        |        M        |  `Bratley1992c()`   |
+|       {ref}`Bratley et al. (1992) D <test-functions:bratley1992d>`        |        M        |  `Bratley1992d()`   |
+|         {ref}`Genz (Continuous) <test-functions:genz-continuous>`         |        M        | `GenzContinuous()`  |
+|        {ref}`Genz (Corner Peak) <test-functions:genz-corner-peak>`        |        M        | `GenzCornerPeak()`  |
+|           {ref}`Genz (Gaussian) <test-functions:genz-gaussian>`           |        M        |  `GenzGaussian()`   |
+|       {ref}`Genz (Product Peak) <test-functions:genz-product-peak>`       |        M        | `GenzProductPeak()` |
+|                 {ref}`Sobol'-G <test-functions:sobol-g>`                  |        M        |     `SobolG()`      |
+|           {ref}`Welch et al. (1992) <test-functions:welch1992>`           |       20        |    `Welch1992()`    |
 
 In a Python terminal, you can list all the available functions relevant
 for metamodeling applications using ``list_functions()`` and filter the results


Test patch:
diff --git a/docs/test-functions/available.md b/docs/test-functions/available.md
index 959fc38..dd38934 100644
--- a/docs/test-functions/available.md
+++ b/docs/test-functions/available.md
@@ -52,6 +52,7 @@ regardless of their typical applications.
 |                    {ref}`Gayton Hat <test-functions:gayton-hat>`                    |        2        |          `GaytonHat()`          |
 |              {ref}`Genz (Continuous) <test-functions:genz-continuous>`              |        M        |       `GenzContinuous()`        |
 |             {ref}`Genz (Corner Peak) <test-functions:genz-corner-peak>`             |        M        |       `GenzCornerPeak()`        |
+|                {ref}`Genz (Gaussian) <test-functions:genz-gaussian>`                |        M        |        `GenzGaussian()`         |
 |            {ref}`Genz (Product Peak) <test-functions:genz-product-peak>`            |        M        |       `GenzProductPeak()`       |
 |           {ref}`Gramacy (2007) 1D Sine <test-functions:gramacy-1d-sine>`            |        1        |        `Gramacy1DSine()`        |
 |               {ref}`Hyper-sphere Bound <test-functions:hyper-sphere>`               |        2        |         `HyperSphere()`         |
diff --git a/docs/test-functions/genz-continuous.md b/docs/test-functions/genz-continuous.md
index ee3caaa..4abaf61 100644
--- a/docs/test-functions/genz-continuous.md
+++ b/docs/test-functions/genz-continuous.md
@@ -102,6 +102,8 @@ designed to test the performance of numerical integration routines:
   a prominent peak at the center of the multidimensional space.
 - {ref}`Corner peak <test-functions:genz-corner-peak>` function features
   a prominent peak in one corner of the multidimensional space.
+- {ref}`Gaussian <test-functions:genz-gaussian>` function features
+  a bell-shaped peak at the center of the multidimensional space.
 - {ref}`Continuous <test-functions:genz-continuous>` function features
   an exponential decay from the center of the multidimensional space.
   The function is continuous everywhere, but non-differentiable at the center.
diff --git a/docs/test-functions/genz-corner-peak.md b/docs/test-functions/genz-corner-peak.md
index c2ff01a..1a6e56d 100644
--- a/docs/test-functions/genz-corner-peak.md
+++ b/docs/test-functions/genz-corner-peak.md
@@ -106,6 +106,8 @@ designed to test the performance of numerical integration routines:
 - {ref}`Corner peak <test-functions:genz-corner-peak>` function features
   a prominent peak in one corner of the multidimensional space.
   (_this function_)
+- {ref}`Gaussian <test-functions:genz-gaussian>` function features
+  a bell-shaped peak at the center of the multidimensional space.
 - {ref}`Continuous <test-functions:genz-continuous>` function features
   an exponential decay from the center of the multidimensional space.
   The function is continuous everywhere, but non-differentiable at the center.
diff --git a/docs/test-functions/genz-gaussian.md b/docs/test-functions/genz-gaussian.md
new file mode 100644
index 0000000..5a856b5
--- /dev/null
+++ b/docs/test-functions/genz-gaussian.md
@@ -0,0 +1,217 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+(test-functions:genz-gaussian)=
+# Genz Gaussian Function
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import uqtestfuns as uqtf
+```
+
+The Genz Gaussian function is an $M$-dimensional scalar-valued
+function commonly used to assess the accuracy of numerical
+integration routines.
+It is one of six functions introduced by Genz {cite}`Genz1984`;
+see the box below.
+
+The function features a multidimensional bell-shaped peak at the center
+of the multidimensional space.
+The plots for one-dimensional and two-dimensional Genz Gaussian function
+with the default parameters can be seen below.
+
+```{code-cell} ipython3
+:tags: [remove-input]
+
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# --- Create 1D data from Genz Gaussian
+my_fun_1d = uqtf.GenzGaussian(input_dimension=1)
+xx_1d = np.linspace(0, 1, 1000)[:, np.newaxis]
+yy_1d = my_fun_1d(xx_1d)
+
+# --- Create 2D data from Genz Gaussian
+my_fun_2d = uqtf.GenzGaussian(input_dimension=2)
+mesh_2d = np.meshgrid(xx_1d, xx_1d)
+xx_2d = np.array(mesh_2d).T.reshape(-1, 2)
+yy_2d = my_fun_2d(xx_2d)
+
+# --- Create two-dimensional plots
+fig = plt.figure(figsize=(15, 5))
+
+# 1D
+axs_1 = plt.subplot(131)
+axs_1.plot(xx_1d, yy_1d, color="#8da0cb")
+axs_1.grid()
+axs_1.set_xlabel("$x$", fontsize=14)
+axs_1.set_ylabel("$\mathcal{M}(x)$", fontsize=14)
+axs_1.set_title("1D Genz Gaussian")
+
+# Surface
+axs_2 = plt.subplot(132, projection='3d')
+axs_2.plot_surface(
+    mesh_2d[0],
+    mesh_2d[1],
+    yy_2d.reshape(1000, 1000).T,
+    cmap="plasma",
+    linewidth=0,
+    antialiased=False,
+    alpha=0.5
+)
+axs_2.set_xlabel("$x_1$", fontsize=14)
+axs_2.set_ylabel("$x_2$", fontsize=14)
+axs_2.set_zlabel("$\mathcal{M}(x_1, x_2)$", fontsize=14)
+axs_2.set_title("Surface plot of 2D Genz Gaussian", fontsize=14)
+
+# Contour
+axs_3 = plt.subplot(133)
+cf = axs_3.contourf(
+    mesh_2d[0], mesh_2d[1], yy_2d.reshape(1000, 1000).T, cmap="plasma"
+)
+axs_3.set_xlabel("$x_1$", fontsize=14)
+axs_3.set_ylabel("$x_2$", fontsize=14)
+axs_3.set_title("Contour plot of 2D Genz Gaussian", fontsize=14)
+divider = make_axes_locatable(axs_3)
+cax = divider.append_axes('right', size='5%', pad=0.05)
+fig.colorbar(cf, cax=cax, orientation='vertical')
+axs_3.axis('scaled')
+
+fig.tight_layout(pad=3.0)
+plt.gcf().set_dpi(150);
+```
+
+```{note}
+Genz {cite}`Genz1984` introduced six challenging
+parameterized $M$-dimensional functions
+designed to test the performance of numerical integration routines:
+
+- {ref}`Product peak <test-functions:genz-product-peak>` function features
+  a prominent peak at the center of the multidimensional space.
+- {ref}`Corner peak <test-functions:genz-corner-peak>` function features
+  a prominent peak in one corner of the multidimensional space.
+- {ref}`Gaussian <test-functions:genz-gaussian>` function features
+  a bell-shaped peak at the center of the multidimensional space.
+  (_this function_)
+- {ref}`Continuous <test-functions:genz-continuous>` function features
+  an exponential decay from the center of the multidimensional space.
+  The function is continuous everywhere, but non-differentiable at the center.
+
+The functions are further characterized by offset (shift) and shape parameters.
+While the offset parameter has minimal impact on the integral's value,
+the shape parameter significantly affects
+the difficulty of the integration problem.
+```
+
+## Test function instance
+
+To create a default instance of the test function:
+
+```{code-cell} ipython3
+my_testfun = uqtf.GenzGaussian()
+```
+
+Check if it has been correctly instantiated:
+
+```{code-cell} ipython3
+print(my_testfun)
+```
+
+By default, the input dimension is set to $2$[^default_dimension].
+To create an instance with another value of input dimension,
+pass an integer to the parameter `input_dimension` (the first parameter).
+For example, to create an instance of the test function in six dimensions,
+type:
+
+```python
+my_testfun = uqtf.GenzGaussian(input_dimension=6)
+```
+
+## Description
+
+The Genz continuous is defined as:
+
+$$
+\mathcal{M}(\boldsymbol{x}; \boldsymbol{a}, \boldsymbol{b}) = \exp{\left[ - \sum_{i = 1}^M a^2_i \left( x_i - b_i \right)^2 \right]},
+$$
+
+where $\boldsymbol{x} = \left( x_1, \ldots, x_M \right)$
+is the $M$-dimensional vector of input variables;
+and $\boldsymbol{a} = \left( a_1, \ldots, a_M \right)$ 
+$\boldsymbol{b} = \left( b_1, \ldots, b_M \right)$ are $M$-dimensional
+vectors corresponding to the (fixed) shape and offset parameters, respectively.
+Further details about these parameters are provided below.
+
+## Probabilistic input
+
+The input specification for the Genz Gaussian function is shown below.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+print(my_testfun.prob_input)
+```
+
+## Parameters
+
+The parameters of the Genz Gaussian function consists
+of the vector of shape (scale) parameters $\boldsymbol{a}$
+and the vector of offset parameters $\boldsymbol{b}$. 
+
+The shape parameters determines the extent of the product peaking;
+Larger values of $\boldsymbol{a}$
+increase the prominence of the peak,
+making the integration problem more challenging.
+The offset parameters, on the other hand, does not affect significantly
+the difficulty of the problem and can be chosen randomly.
+
+The default parameter is shown below.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+print(my_testfun.parameters)
+```
+
+## Reference results
+
+This section provides several reference results of typical UQ analyses involving
+the test function.
+
+### Sample histogram
+
+Shown below is the histogram of the output based on $100'000$ random points:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+my_testfun.prob_input.reset_rng(42)
+xx_test = my_testfun.prob_input.get_sample(100000)
+yy_test = my_testfun(xx_test)
+
+plt.hist(yy_test, bins="auto", color="#8da0cb");
+plt.grid();
+plt.ylabel("Counts [-]");
+plt.xlabel("$\mathcal{M}(\mathbf{X})$");
+plt.gcf().set_dpi(150);
+```
+
+## References
+
+```{bibliography}
+:filter: docname in docnames
+```
+
+[^default_dimension]: This default dimension applies to all variable dimension
+test functions. It will be used if the `input_dimension` argument is not given.
\ No newline at end of file
diff --git a/docs/test-functions/genz-product-peak.md b/docs/test-functions/genz-product-peak.md
index ec176e7..d4d574c 100644
--- a/docs/test-functions/genz-product-peak.md
+++ b/docs/test-functions/genz-product-peak.md
@@ -101,6 +101,8 @@ designed to test the performance of numerical integration routines:
   (_this function_)
 - {ref}`Corner peak <test-functions:genz-corner-peak>` function features
   a prominent peak in one corner of the multidimensional space.
+- {ref}`Gaussian <test-functions:genz-gaussian>` function features
+  a bell-shaped peak at the center of the multidimensional space.
 - {ref}`Continuous <test-functions:genz-continuous>` function features
   an exponential decay from the center of the multidimensional space.
   The function is continuous everywhere, but non-differentiable at the center.
diff --git a/src/uqtestfuns/test_functions/__init__.py b/src/uqtestfuns/test_functions/__init__.py
index d35370e..7b8cbda 100644
--- a/src/uqtestfuns/test_functions/__init__.py
+++ b/src/uqtestfuns/test_functions/__init__.py
@@ -20,7 +20,7 @@
 from .franke import Franke1, Franke2, Franke3, Franke4, Franke5, Franke6
 from .friedman import Friedman6D, Friedman10D
 from .gayton_hat import GaytonHat
-from .genz import GenzContinuous, GenzCornerPeak, GenzProductPeak
+from .genz import GenzContinuous, GenzCornerPeak, GenzGaussian, GenzProductPeak
 from .gramacy2007 import Gramacy1DSine
 from .hyper_sphere import HyperSphere
 from .ishigami import Ishigami
@@ -82,6 +82,7 @@
     "GaytonHat",
     "GenzContinuous",
     "GenzCornerPeak",
+    "GenzGaussian",
     "GenzProductPeak",
     "Gramacy1DSine",
     "HyperSphere",
diff --git a/src/uqtestfuns/test_functions/genz.py b/src/uqtestfuns/test_functions/genz.py
index a03c13f..3adedf6 100644
--- a/src/uqtestfuns/test_functions/genz.py
+++ b/src/uqtestfuns/test_functions/genz.py
@@ -39,7 +39,12 @@
 )
 from uqtestfuns.core.uqtestfun_abc import UQTestFunVarDimABC
 
-__all__ = ["GenzContinuous", "GenzCornerPeak", "GenzProductPeak"]
+__all__ = [
+    "GenzContinuous",
+    "GenzCornerPeak",
+    "GenzGaussian",
+    "GenzProductPeak",
+]
 
 
 MARGINALS_GENZ1984: MarginalSpecs = [
@@ -260,6 +265,80 @@ class GenzProductPeak(UQTestFunVarDimABC):
     evaluate = staticmethod(evaluate_product_peak)  # type: ignore
 
 
+def evaluate_gaussian(
+    xx: np.ndarray, aa: np.ndarray, bb: np.ndarray
+) -> np.ndarray:
+    """Compute the Gaussian function on a set of input values.
+
+    Parameters
+    ----------
+    xx : np.ndarray
+        M-Dimensional input values given by an N-by-M array where
+        N is the number of input values.
+    aa : np.ndarray
+        The vector of shape parameters of the Genz family of integrand
+        functions; the larger the value the more difficult the integrations.
+        The length of the vector must M.
+    bb : np.ndarray
+        The vector of offset parameters of the Genz family of integrand
+        functions; the values do not alter significantly the difficulty of
+        the integration problem.
+        The length of the vector must be M.
+
+    Returns
+    -------
+    np.ndarray
+        The output of the test function evaluated on the input values.
+        The output is a 1-dimensional array of length N.
+    """
+    yy = np.exp(-1 * np.sum((xx - bb) ** 2 * aa**2, axis=1))
+
+    return yy
+
+
+class GenzGaussian(UQTestFunVarDimABC):
+    """A concrete implementation of the Genz Gaussian function."""
+
+    _tags = ["integration"]
+    _description = "Gaussian integrand from Genz (1984)"
+    _available_inputs: ProbInputSpecs = {
+        "Genz1984": {
+            "function_id": "GenzGaussian",
+            "description": (
+                "Input specification for the Genz family of integrand"
+            ),
+            "marginals": MARGINALS_GENZ1984,
+            "copulas": None,
+        }
+    }
+    _available_parameters: FunParamSpecs = {
+        "Genz1984": {
+            "function_id": "GenzGaussian",
+            "description": (
+                "Parameter set for the continuous function from Genz (1984);"
+                " constant shape and offset values"
+            ),
+            "declared_parameters": [
+                {
+                    "keyword": "aa",
+                    "value": _get_aa_genz_1984,
+                    "type": np.ndarray,
+                    "description": "Shape parameter",
+                },
+                {
+                    "keyword": "bb",
+                    "value": _get_bb_genz_1984,
+                    "type": np.ndarray,
+                    "description": "Offset parameter",
+                },
+            ],
+        },
+    }
+    _default_parameters_id = "Genz1984"
+
+    evaluate = staticmethod(evaluate_gaussian)  # type: ignore
+
+
 def evaluate_continuous(
     xx: np.ndarray, aa: np.ndarray, bb: np.ndarray
 ) -> np.ndarray: