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/docs/_toc.yml b/docs/_toc.yml
index a21924b..8d35e34 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -100,6 +100,8 @@ parts:
             title: Genz (Corner Peak)
           - file: test-functions/genz-gaussian
             title: Genz (Gaussian)
+          - file: test-functions/genz-discontinuous
+            title: Genz (Discontinuous)
           - file: test-functions/genz-oscillatory
             title: Genz (Oscillatory)
           - file: test-functions/genz-product-peak
diff --git a/docs/fundamentals/integration.md b/docs/fundamentals/integration.md
index 4025582..0733aed 100644
--- a/docs/fundamentals/integration.md
+++ b/docs/fundamentals/integration.md
@@ -18,19 +18,20 @@ 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 (Gaussian) <test-functions:genz-gaussian>`          |        M        |   `GenzGaussian()`   |
-|       {ref}`Genz (Oscillatory) <test-functions:genz-oscillatory>`       |        M        | `GenzOscillatory()`  |
-|      {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 (Discontinuous) <test-functions:genz-discontinuous>`  |        M        | `GenzDiscontinuous()`  |
+|      {ref}`Genz (Gaussian) <test-functions:genz-gaussian>`       |        M        |    `GenzGaussian()`    |
+|   {ref}`Genz (Oscillatory) <test-functions:genz-oscillatory>`    |        M        |  `GenzOscillatory()`   |
+|  {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
diff --git a/docs/fundamentals/sensitivity.md b/docs/fundamentals/sensitivity.md
index 84c11c8..7c67e2b 100644
--- a/docs/fundamentals/sensitivity.md
+++ b/docs/fundamentals/sensitivity.md
@@ -18,31 +18,32 @@ kernelspec:
 The table below listed the available test functions typically used
 in the comparison of sensitivity analysis methods.
 
-|                             Name                             | Input Dimension |     Constructor      |
-|:------------------------------------------------------------:|:---------------:|:--------------------:|
-|          {ref}`Borehole <test-functions:borehole>`           |        8        |     `Borehole()`     |
-| {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}`Damped Oscillator <test-functions:damped-oscillator>`  |        7        | `DampedOscillator()` |
-|             {ref}`Flood <test-functions:flood>`              |        8        |      `Flood()`       |
-|      {ref}`Friedman (6D) <test-functions:friedman-6d>`       |        6        |    `Friedman6D()`    |
-| {ref}`Genz (Corner Peak) <test-functions:genz-corner-peak>`  |        M        |  `GenzCornerPeak()`  |
-|          {ref}`Ishigami <test-functions:ishigami>`           |        3        |     `Ishigami()`     |
-|        {ref}`Moon (2010) 3D <test-functions:moon3d>`         |        3        |      `Moon3D()`      |
-|   {ref}`Morris et al. (2006) <test-functions:morris2006>`    |        M        |    `Morris2006()`    |
-|       {ref}`OTL Circuit <test-functions:otl-circuit>`        |     6 / 20      |    `OTLCircuit()`    |
-|       {ref}`Piston Simulation <test-functions:piston>`       |     7 / 20      |      `Piston()`      |
-| {ref}`Simple Portfolio Model <test-functions:portfolio-3d>`  |        3        |   `Portfolio3D()`    |
-|    {ref}`SaltelliLinear <test-functions:saltelli-linear>`    |        M        |  `SaltelliLinear()`  |
-|           {ref}`Sobol'-G <test-functions:sobol-g>`           |        M        |      `SobolG()`      |
-|        {ref}`Sobol'-G* <test-functions:sobol-g-star>`        |        M        |    `SobolGStar()`    |
-|     {ref}`Sobol'-Levitan <test-functions:sobol-levitan>`     |        M        |   `SobolLevitan()`   |
-|     {ref}`Solar Cell Model <test-functions:solar-cell>`      |        5        |    `SolarCell()`     |
-|            {ref}`Sulfur <test-functions:sulfur>`             |        9        |      `Sulfur()`      |
-|    {ref}`Welch et al. (1992) <test-functions:welch1992>`     |       20        |    `Welch1992()`     |
-|       {ref}`Wing Weight <test-functions:wing-weight>`        |       10        |    `WingWeight()`    |
+|                               Name                               | Input Dimension |      Constructor       |
+|:----------------------------------------------------------------:|:---------------:|:----------------------:|
+|            {ref}`Borehole <test-functions:borehole>`             |        8        |      `Borehole()`      |
+|   {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}`Damped Oscillator <test-functions:damped-oscillator>`    |        7        |  `DampedOscillator()`  |
+|               {ref}`Flood <test-functions:flood>`                |        8        |       `Flood()`        |
+|        {ref}`Friedman (6D) <test-functions:friedman-6d>`         |        6        |     `Friedman6D()`     |
+|   {ref}`Genz (Corner Peak) <test-functions:genz-corner-peak>`    |        M        |   `GenzCornerPeak()`   |
+| {ref}`Genz (Discontinuous) <test-functions:genz-discontinuous>`  |        M        | `GenzDiscontinuous()`  |
+|            {ref}`Ishigami <test-functions:ishigami>`             |        3        |      `Ishigami()`      |
+|          {ref}`Moon (2010) 3D <test-functions:moon3d>`           |        3        |       `Moon3D()`       |
+|     {ref}`Morris et al. (2006) <test-functions:morris2006>`      |        M        |     `Morris2006()`     |
+|         {ref}`OTL Circuit <test-functions:otl-circuit>`          |     6 / 20      |     `OTLCircuit()`     |
+|         {ref}`Piston Simulation <test-functions:piston>`         |     7 / 20      |       `Piston()`       |
+|   {ref}`Simple Portfolio Model <test-functions:portfolio-3d>`    |        3        |    `Portfolio3D()`     |
+|      {ref}`SaltelliLinear <test-functions:saltelli-linear>`      |        M        |   `SaltelliLinear()`   |
+|             {ref}`Sobol'-G <test-functions:sobol-g>`             |        M        |       `SobolG()`       |
+|          {ref}`Sobol'-G* <test-functions:sobol-g-star>`          |        M        |     `SobolGStar()`     |
+|       {ref}`Sobol'-Levitan <test-functions:sobol-levitan>`       |        M        |    `SobolLevitan()`    |
+|       {ref}`Solar Cell Model <test-functions:solar-cell>`        |        5        |     `SolarCell()`      |
+|              {ref}`Sulfur <test-functions:sulfur>`               |        9        |       `Sulfur()`       |
+|      {ref}`Welch et al. (1992) <test-functions:welch1992>`       |       20        |     `Welch1992()`      |
+|         {ref}`Wing Weight <test-functions:wing-weight>`          |       10        |     `WingWeight()`     |
 
 In a Python terminal, you can list all the available functions relevant
 for metamodeling applications using ``list_functions()``


Test patch:
diff --git a/docs/test-functions/available.md b/docs/test-functions/available.md
index 683179e..d49ebe1 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 (Discontinuous) <test-functions:genz-discontinuous>`           |        M        |      `GenzDiscontinuous()`      |
 |                {ref}`Genz (Gaussian) <test-functions:genz-gaussian>`                |        M        |        `GenzGaussian()`         |
 |             {ref}`Genz (Oscillatory) <test-functions:genz-oscillatory>`             |        M        |       `GenzOscillatory()`       |
 |            {ref}`Genz (Product Peak) <test-functions:genz-product-peak>`            |        M        |       `GenzProductPeak()`       |
diff --git a/docs/test-functions/genz-continuous.md b/docs/test-functions/genz-continuous.md
index ab2d7c2..07ea7c0 100644
--- a/docs/test-functions/genz-continuous.md
+++ b/docs/test-functions/genz-continuous.md
@@ -110,6 +110,10 @@ designed to test the performance of numerical integration routines:
   an exponential decay from the center of the multidimensional space.
   The function is continuous everywhere, but non-differentiable at the center.
   (_this function_)
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity.
 
 The functions are further characterized by offset (shift) and shape parameters.
 While the offset parameter has minimal impact on the integral's value,
diff --git a/docs/test-functions/genz-corner-peak.md b/docs/test-functions/genz-corner-peak.md
index c18c2bf..ecbe6d5 100644
--- a/docs/test-functions/genz-corner-peak.md
+++ b/docs/test-functions/genz-corner-peak.md
@@ -113,6 +113,10 @@ designed to test the performance of numerical integration routines:
 - {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.
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity.
 
 The functions are further characterized by offset (shift) and shape parameters.
 While the offset parameter has minimal impact on the integral's value,
diff --git a/docs/test-functions/genz-discontinuous.md b/docs/test-functions/genz-discontinuous.md
new file mode 100644
index 0000000..4399f15
--- /dev/null
+++ b/docs/test-functions/genz-discontinuous.md
@@ -0,0 +1,238 @@
+---
+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-discontinuous)=
+# Genz Discontinuous Function
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import uqtestfuns as uqtf
+```
+
+The Genz discontinuous 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 was later featured in {cite}`Zhang2014` as a test function
+for a global sensitivity analysis method albeit in a non-standard form[^non-standard].
+
+The function features an exponential function on one particular corner
+of the multidimensional space.
+Outside the boundary, the function values drop to zero creating discontinuity.
+The plots for one-dimensional and two-dimensional Genz discontinuous 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 Discontinuous
+my_fun_1d = uqtf.GenzDiscontinuous(input_dimension=1)
+xx_1d = np.linspace(0, 1, 1000)[:, np.newaxis]
+yy_1d = my_fun_1d(xx_1d)
+
+# --- Create 2D data from Genz Discontinuous
+my_fun_2d = uqtf.GenzDiscontinuous(input_dimension=2)
+my_fun_2d.parameters["aa"] = np.array([2, 2])
+my_fun_2d.parameters["bb"] = np.array([0.8, 0.3])
+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 discontinuous")
+
+# 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 discontinuous", 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 discontinuous", 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}`Oscillatory <test-functions:genz-oscillatory>` function features
+  an oscillating shape in the multidimensional space.
+- {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.
+- {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.
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity. (_this function_)
+
+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.GenzDiscontinuous()
+```
+
+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.GenzDiscontinuous(input_dimension=6)
+```
+
+## Description
+
+The Genz continuous is defined as:
+
+$$
+\mathcal{M}(\boldsymbol{x}; \boldsymbol{a}, \boldsymbol{b}) = 
+\begin{cases}
+0, & x_i > b_i, i = 1, \ldots, M,\\
+\exp{\left( \sum_{i = 1}^M a_i x_i \right)}, & \text{otherwise}
+\end{cases}
+$$
+
+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 continuous function is shown below.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+print(my_testfun.prob_input)
+```
+
+## Parameters
+
+The parameters of the Genz corner peak 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, do not affect significantly
+the difficulty of the problem and can be chosen randomly.
+However, larger values of the offset parameters result in a larger integral,
+as the non-zero region of the function expands.
+
+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.
+
+[^non-standard]: referring to the formula, the one appeared in {cite}`Zhang2014`
+have zero-value bounds flip, e.g., $x_1 < b_1, x_2 > b_2$ for two-dimensional
+case.
\ No newline at end of file
diff --git a/docs/test-functions/genz-gaussian.md b/docs/test-functions/genz-gaussian.md
index 246c57f..2fd5231 100644
--- a/docs/test-functions/genz-gaussian.md
+++ b/docs/test-functions/genz-gaussian.md
@@ -109,6 +109,10 @@ designed to test the performance of numerical integration routines:
 - {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.
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity.
 
 The functions are further characterized by offset (shift) and shape parameters.
 While the offset parameter has minimal impact on the integral's value,
diff --git a/docs/test-functions/genz-oscillatory.md b/docs/test-functions/genz-oscillatory.md
index 563ac74..753e19f 100644
--- a/docs/test-functions/genz-oscillatory.md
+++ b/docs/test-functions/genz-oscillatory.md
@@ -108,6 +108,10 @@ designed to test the performance of numerical integration routines:
 - {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.
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity.
 
 The functions are further characterized by offset (shift) and shape parameters.
 While the offset parameter has minimal impact on the integral's value,
diff --git a/docs/test-functions/genz-product-peak.md b/docs/test-functions/genz-product-peak.md
index 319ea9d..011d67b 100644
--- a/docs/test-functions/genz-product-peak.md
+++ b/docs/test-functions/genz-product-peak.md
@@ -108,7 +108,10 @@ designed to test the performance of numerical integration routines:
 - {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.
-  
+- {ref}`Discontinuous <test-functions:genz-discontinuous>` function features
+  an exponential rise from corner of the multidimensional space up to the
+  offset parameter value, after which the function value drops to zero
+  everywhere, creating discontinuity.
 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
diff --git a/src/uqtestfuns/test_functions/__init__.py b/src/uqtestfuns/test_functions/__init__.py
index 31689c2..d32e3c1 100644
--- a/src/uqtestfuns/test_functions/__init__.py
+++ b/src/uqtestfuns/test_functions/__init__.py
@@ -21,6 +21,7 @@
 from .friedman import Friedman6D, Friedman10D
 from .gayton_hat import GaytonHat
 from .genz import (
+    GenzDiscontinuous,
     GenzContinuous,
     GenzCornerPeak,
     GenzGaussian,
@@ -88,6 +89,7 @@
     "GaytonHat",
     "GenzContinuous",
     "GenzCornerPeak",
+    "GenzDiscontinuous",
     "GenzGaussian",
     "GenzOscillatory",
     "GenzProductPeak",
diff --git a/src/uqtestfuns/test_functions/genz.py b/src/uqtestfuns/test_functions/genz.py
index 2b2453e..ef57e5a 100644
--- a/src/uqtestfuns/test_functions/genz.py
+++ b/src/uqtestfuns/test_functions/genz.py
@@ -4,9 +4,23 @@
 Genz [1] introduced six challenging parameterized $M$-dimensional functions
 designed to test the performance of numerical integration routines:
 
+- Oscillatory function features an oscillating shape in the multidimensional
+  space.
+- Product peak function features a prominent peak at the center
+  of the multidimensional space.
 - Corner peak function features a prominent peak in one corner
   of the multidimensional space. This function was featured as a test function
-  in sensitivity analysis ([2]) and metamodeling [3] exercises.
+  in sensitivity analysis ([2]) and metamodeling ([3]) exercises.
+- Gaussian function features a bell-shaped peak at the center
+  of the multidimensional space.
+- Continuous function features an exponential decay from the center
+  of the multidimensional space. The function is continuous everywhere,
+  but non-differentiable at the center.
+- Discontinuous function features an exponential rise from corner
+  of the multidimensional space up to the offset parameter value,
+  after which the function value drops to zero everywhere,
+  creating discontinuity. This function was featured as a test function
+  in sensitivity analysis ([2])
 
 The functions are further characterized by shift and scale parameters.
 While the shift parameter has minimal impact on the integral's value,
@@ -42,6 +56,7 @@
 __all__ = [
     "GenzContinuous",
     "GenzCornerPeak",
+    "GenzDiscontinuous",
     "GenzGaussian",
     "GenzOscillatory",
     "GenzProductPeak",
@@ -484,3 +499,83 @@ class GenzContinuous(UQTestFunVarDimABC):
     _default_parameters_id = "Genz1984"
 
     evaluate = staticmethod(evaluate_continuous)  # type: ignore
+
+
+def evaluate_discontinuous(
+    xx: np.ndarray, aa: np.ndarray, bb: np.ndarray
+) -> np.ndarray:
+    """Compute the discontinuous 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(np.sum(xx * aa, axis=1))
+
+    # Filter the output to make the function discontinuous
+    mask = np.zeros(xx.shape[0], dtype=bool)
+    for i in range(xx.shape[1]):
+        mask |= xx[:, i] > bb[i]
+    yy[mask] = 0
+
+    return yy
+
+
+class GenzDiscontinuous(UQTestFunVarDimABC):
+    """A concrete implementation of the Genz discontinuous function."""
+
+    _tags = ["integration", "sensitivity"]
+    _description = "Discontinuous integrand from Genz (1984)"
+    _available_inputs: ProbInputSpecs = {
+        "Genz1984": {
+            "function_id": "GenzDiscontinuous",
+            "description": (
+                "Input specification for the Genz family of integrand"
+            ),
+            "marginals": MARGINALS_GENZ1984,
+            "copulas": None,
+        }
+    }
+    _available_parameters: FunParamSpecs = {
+        "Genz1984": {
+            "function_id": "GenzDiscontinuous",
+            "description": (
+                "Parameter set for the discontinuous 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_discontinuous)  # type: ignore