Merge pull request #192 from damar-wicaksono/dev-189

Add an implementation of the McLain S5 function

CHANGELOG.md

@@ -18,6 +18,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   exercises, is added as a UQ test function.
 - The two-dimensional (third) Franke function, relevant for metamodeling
   exercises, is added as a UQ test function.
+- The two-dimensional McLain S5 function, relevant for metamodeling exercises,
+  is added as a UQ test function
 
 ## [0.1.1] - 2023-07-03
 

docs/_toc.yml

@@ -41,6 +41,8 @@ parts:
             title: Oakley-O'Hagan 1D
           - file: test-functions/otl-circuit
             title: OTL Circuit
+          - file: test-functions/mclain-s5
+            title: McLain S5
           - file: test-functions/piston
             title: Piston Simulation
           - file: test-functions/sobol-g

docs/fundamentals/metamodeling.md

@@ -28,6 +28,7 @@ in the comparison of metamodeling approaches.
 |        {ref}`(3rd) Franke <test-functions:franke-3>`        |         2         |     `Franke3()`      |
 | {ref}`Oakley-O'Hagan 1D <test-functions:oakley-ohagan-1d>`  |         1         |  `OakleyOHagan1D()`  |
 |       {ref}`OTL Circuit <test-functions:otl-circuit>`       |      6 / 20       |    `OTLCircuit()`    |
+|         {ref}`McLain S5 <test-functions:mclain-s5>`         |         2         |     `McLainS5()`     |
 |      {ref}`Piston Simulation <test-functions:piston>`       |      7 / 20       |      `Piston()`      |
 |            {ref}`Sulfur <test-functions:sulfur>`            |         9         |      `Sulfur()`      |
 |       {ref}`Wing Weight <test-functions:wing-weight>`       |        10         |    `WingWeight()`    |

docs/test-functions/available.md

@@ -30,6 +30,7 @@ available in the current UQTestFuns, regardless of their typical applications.
 |          {ref}`Ishigami <test-functions:ishigami>`          |         3         |     `Ishigami()`     |
 | {ref}`Oakley-O'Hagan 1D <test-functions:oakley-ohagan-1d>`  |         1         |  `OakleyOHagan1D()`  |
 |       {ref}`OTL Circuit <test-functions:otl-circuit>`       |      6 / 20       |    `OTLCircuit()`    |
+|         {ref}`McLain S5 <test-functions:mclain-s5>`         |         2         |     `McLainS5()`     |
 |      {ref}`Piston Simulation <test-functions:piston>`       |      7 / 20       |      `Piston()`      |
 |          {ref}`Sobol'-G <test-functions:sobol-g>`           |         M         |      `SobolG()`      |
 |            {ref}`Sulfur <test-functions:sulfur>`            |         9         |      `Sulfur()`      |

docs/test-functions/mclain-s5.md

@@ -0,0 +1,159 @@
+---
+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:mclain-s5)=
+# McLain S5 Function
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import uqtestfuns as uqtf
+```
+
+The McLain S5 function is a two-dimensional scalar-valued function.
+The function was introduced in {cite}`McLain1974` in the context of a procedure 
+for drawing contours from a given set of points.
+
+```{note}
+The McLain's test functions are a set of five two-dimensional functions 
+that mathematically defines a surface. The functions are:
+
+- {ref}`S5 <test-functions:mclain-s5>`: A plateau and plain separated by a steep cliff (this function)
+```
+
+```{code-cell} ipython3
+:tags: [remove-input]
+
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+my_fun = uqtf.McLainS5()
+
+# --- Create 2D data
+xx_1d = np.linspace(1.0, 10.0, 1000)[:, np.newaxis]
+mesh_2d = np.meshgrid(xx_1d, xx_1d)
+xx_2d = np.array(mesh_2d).T.reshape(-1, 2)
+yy_2d = my_fun(xx_2d)
+
+# --- Create two-dimensional plots
+fig = plt.figure(figsize=(10, 5))
+
+# Surface
+axs_1 = plt.subplot(121, projection='3d')
+axs_1.plot_surface(
+    mesh_2d[0],
+    mesh_2d[1],
+    yy_2d.reshape(1000,1000),
+    linewidth=0,
+    cmap="plasma",
+    antialiased=False,
+    alpha=0.5
+)
+axs_1.set_xlabel("$x_1$", fontsize=14)
+axs_1.set_ylabel("$x_2$", fontsize=14)
+axs_1.set_zlabel("$\mathcal{M}(x_1, x_2)$", fontsize=14)
+axs_1.set_title("Surface plot of McLain S5", fontsize=14)
+
+# Contour
+axs_2 = plt.subplot(122)
+cf = axs_2.contourf(
+    mesh_2d[0], mesh_2d[1], yy_2d.reshape(1000, 1000), cmap="plasma"
+)
+axs_2.set_xlabel("$x_1$", fontsize=14)
+axs_2.set_ylabel("$x_2$", fontsize=14)
+axs_2.set_title("Contour plot of McLain S5", fontsize=14)
+divider = make_axes_locatable(axs_2)
+cax = divider.append_axes('right', size='5%', pad=0.05)
+fig.colorbar(cf, cax=cax, orientation='vertical')
+axs_2.axis('scaled')
+
+fig.tight_layout(pad=4.0)
+plt.gcf().set_dpi(75);
+```
+
+As shown in the plots above, the function features two nearly flat regions 
+of height $0.0$ and (approximately) $1.0$.
+The two regions are joined by a sharp rise that runs diagonally from
+$(0.0, 10.0)$ to $(10.0, 0.0)$.
+
+```{note}
+The McLain S5 function appeared in a modified form in the report 
+of Franke {cite}`Franke1979`
+(specifically the {ref}`(2nd) Franke function <test-functions:franke-2>`).
+
+In fact, four of Franke's test functions are 
+slight modifications of McLain's, including the translation of the input domain
+from $[1.0, 10.0]$ to $[0.0, 1.0]$.
+```
+
+## Test function instance
+
+To create a default instance of the McLain S5 function:
+
+```{code-cell} ipython3
+my_testfun = uqtf.McLainS5()
+```
+
+Check if it has been correctly instantiated:
+
+```{code-cell} ipython3
+print(my_testfun)
+```
+
+## Description
+
+The (2nd) Franke function is defined as follows:
+
+$$
+\mathcal{M}(\boldsymbol{x}) = \tanh{\left( x_1 + x_2 - 11 \right)}
+$$
+where $\boldsymbol{x} = \{ x_1, x_2 \}$
+is the two-dimensional vector of input variables further defined below.
+
+## Probabilistic input
+
+Based on {cite}`McLain1974`, the probabilistic input model
+for the function consists of two independent random variables as shown below.
+
+```{code-cell} ipython3
+my_testfun.prob_input
+```
+
+## 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]
+
+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}
+:style: plain
+:filter: docname in docnames
+```

src/uqtestfuns/test_functions/__init__.py

@@ -11,6 +11,7 @@
 from .ishigami import Ishigami
 from .oakley_ohagan_1d import OakleyOHagan1D
 from .otl_circuit import OTLCircuit
+from .mclain import McLainS5
 from .piston import Piston
 from .sobol_g import SobolG
 from .sulfur import Sulfur
@@ -30,6 +31,7 @@
     "Ishigami",
     "OakleyOHagan1D",
     "OTLCircuit",
+    "McLainS5",
     "Piston",
     "SobolG",
     "Sulfur",

src/uqtestfuns/test_functions/mclain.py

@@ -0,0 +1,115 @@
+"""
+Module with an implementation of the McLain's test functions.
+
+The McLain's test functions consists of five two-dimensional scalar-valued
+functions. The functions were introduced in [1] in the context of drawing
+contours from a given set of points.
+
+There are five test functions in McLain's paper each models a mathematically
+defined surface:
+
+- S5: A plateau and plain separated by a steep cliff
+
+Four of the functions (S2-S5) appeared in modified forms in [2].
+
+References
+----------
+
+1. D. H. McLain, "Drawing contours from arbitrary data points," The Computer
+   Journal, vol. 17, no. 4, pp. 318-324, 1974.
+   DOI: 10.1093/comjnl/17.4.318
+2. Richard Franke, "A critical comparison of some methods for interpolation
+   of scattered data," Naval Postgraduate School, Monterey, Canada,
+   Technical Report No. NPS53-79-003, 1979.
+   URL: https://core.ac.uk/reader/36727660
+"""
+import numpy as np
+
+from typing import Optional
+
+from ..core.prob_input.univariate_distribution import UnivDist
+from ..core.uqtestfun_abc import UQTestFunABC
+from .available import create_prob_input_from_available
+
+__all__ = ["McLainS5"]
+
+INPUT_MARGINALS_MCLAIN1974 = [  # From Ref. [1]
+    UnivDist(
+        name="X1",
+        distribution="uniform",
+        parameters=[1.0, 10.0],
+    ),
+    UnivDist(
+        name="X2",
+        distribution="uniform",
+        parameters=[1.0, 10.0],
+    ),
+]
+
+AVAILABLE_INPUT_SPECS = {
+    "McLain1974": {
+        "name": "McLain-1974",
+        "description": (
+            "Input specification for the McLain's test functions "
+            "from McLain (1974)."
+        ),
+        "marginals": INPUT_MARGINALS_MCLAIN1974,
+        "copulas": None,
+    }
+}
+
+DEFAULT_INPUT_SELECTION = "McLain1974"
+
+
+class McLainS5(UQTestFunABC):
+    """A concrete implementation of the McLain S5 function."""
+
+    _TAGS = ["metamodeling"]
+
+    _AVAILABLE_INPUTS = tuple(AVAILABLE_INPUT_SPECS.keys())
+
+    _AVAILABLE_PARAMETERS = None
+
+    _DEFAULT_SPATIAL_DIMENSION = 2
+
+    _DESCRIPTION = "McLain S5 function from McLain (1974)"
+
+    def __init__(
+        self,
+        *,
+        prob_input_selection: Optional[str] = DEFAULT_INPUT_SELECTION,
+        name: Optional[str] = None,
+        rng_seed_prob_input: Optional[int] = None,
+    ):
+        # --- Arguments processing
+        prob_input = create_prob_input_from_available(
+            prob_input_selection,
+            AVAILABLE_INPUT_SPECS,
+            rng_seed=rng_seed_prob_input,
+        )
+        # Process the default name
+        if name is None:
+            name = McLainS5.__name__
+
+        super().__init__(prob_input=prob_input, name=name)
+
+    def evaluate(self, xx: np.ndarray):
+        """Evaluate the McLain S5 function on a set of input values.
+
+        Parameters
+        ----------
+        xx : np.ndarray
+            Two-Dimensional input values given by N-by-2 arrays where
+            N is the number of input values.
+
+        Returns
+        -------
+        np.ndarray
+            The output of the McLain S5 function evaluated
+            on the input values.
+            The output is a 1-dimensional array of length N.
+        """
+        # Compute the (second) Franke function
+        yy = np.tanh(xx[:, 0] + xx[:, 1] - 11)
+
+        return yy