Add pert distribution to gen_kw

Author: larsevj

docs/ert/reference/configuration/data_types.rst

@@ -394,6 +394,42 @@ To define a triangular distribution with a minimum of 1, mode (peak) of 3, and m
 
 .. image:: fig/triangular.png
 
+PERT: PERT Distribution
+^^^^^^^^^^^^^^^^^^^^^^^
+
+The ``PERT`` keyword defines a beta-PERT distribution from a minimum, mode,
+and maximum value. An optional scale controls how strongly the distribution
+is concentrated around the mode.
+
+Syntax
+~~~~~~
+::
+
+    VAR PERT MIN MODE MAX [SCALE]
+
+Parameters
+~~~~~~~~~~
+- **MIN**: The minimum value of the distribution.
+- **MODE**: The most likely value. It must be strictly between ``MIN`` and
+  ``MAX``.
+- **MAX**: The maximum value of the distribution.
+- **SCALE**: An optional value greater than zero. The default is ``4``.
+  Larger values concentrate the distribution more strongly around ``MODE``.
+
+Examples
+~~~~~~~~
+To use the default scale:
+
+::
+
+    VAR_PERT PERT 1 3 5
+
+To use an explicit scale:
+
+::
+
+    VAR_PERT_SCALED PERT 1 3 5 2
+
 3D field parameters: ``FIELD``
 ------------------------------
 

src/ert/config/distribution.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 from pydantic import BaseModel, Field, field_validator, model_validator
-from scipy.special import ndtr
+from scipy.special import betaincinv, ndtr
 
 from .parsing import ConfigValidationError, ConfigWarning, ErrorInfo
 
@@ -28,6 +28,10 @@ def get_param_names(cls) -> list[str]:
             if field.is_required() or (field.default is not None and name != "name")
         ]
 
+    @classmethod
+    def get_min_param_count(cls) -> int:
+        return len(cls.get_param_names())
+
 
 class UnifSettings(TransSettingsValidation):
     name: Literal["uniform"] = "uniform"
@@ -380,6 +384,56 @@ def transform_numpy(self, x: np.ndarray) -> np.ndarray:
         return np.clip(result, self.min, self.max)
 
 
+class PertSettings(TransSettingsValidation):
+    model_config = {"allow_inf_nan": False}
+
+    name: Literal["pert"] = "pert"
+    min: float = 0.0
+    mode: float = 0.5
+    max: float = 1.0
+    scale: float = 4.0
+
+    @classmethod
+    def get_min_param_count(cls) -> int:
+        return 3
+
+    @model_validator(mode="after")
+    def valid_pert_params(self) -> Self:
+        errors = []
+        if not self.min < self.max:
+            errors.append(
+                ErrorInfo(
+                    message=f"Minimum {self.min} must be strictly less than"
+                    f" the maximum {self.max} for pert distribution"
+                )
+            )
+        if not (self.min < self.mode < self.max):
+            errors.append(
+                ErrorInfo(
+                    message=f"The mode {self.mode} must be strictly between"
+                    f" the minimum {self.min} and maximum {self.max}"
+                    " for pert distribution"
+                )
+            )
+        if self.scale <= 0:
+            errors.append(
+                ErrorInfo(
+                    message=f"Scale {self.scale} must be strictly greater than 0"
+                    " for pert distribution"
+                )
+            )
+        if errors:
+            raise ConfigValidationError.from_collected(errors)
+        return self
+
+    def transform_numpy(self, x: np.ndarray) -> np.ndarray:
+        span = self.max - self.min
+        alpha = 1 + self.scale * (self.mode - self.min) / span
+        beta = 1 + self.scale * (self.max - self.mode) / span
+        result = self.min + span * betaincinv(alpha, beta, ndtr(x))
+        return np.clip(result, self.min, self.max)
+
+
 DistributionSettings = Annotated[
     UnifSettings
     | LogNormalSettings
@@ -391,7 +445,8 @@ def transform_numpy(self, x: np.ndarray) -> np.ndarray:
     | TruncNormalSettings
     | ErrfSettings
     | DerrfSettings
-    | TriangularSettings,
+    | TriangularSettings
+    | PertSettings,
     Field(discriminator="name"),
 ]
 
@@ -407,6 +462,7 @@ def transform_numpy(self, x: np.ndarray) -> np.ndarray:
     "TRIANGULAR": TriangularSettings,
     "ERRF": ErrfSettings,
     "DERRF": DerrfSettings,
+    "PERT": PertSettings,
 }
 
 

src/ert/config/gen_kw_config.py

@@ -371,7 +371,12 @@ def _parse_distribution(
             )
         dist_cls = DISTRIBUTION_CLASSES[dist_name]
 
-        if len(values) != len(dist_cls.get_param_names()):
+        param_count = len(values)
+        if not (
+            dist_cls.get_min_param_count()
+            <= param_count
+            <= len(dist_cls.get_param_names())
+        ):
             raise ConfigValidationError.with_context(
                 f"Incorrect number of values: {values}, provided for variable "
                 f"{param_name} with distribution {dist_name}.",

tests/ert/ui_tests/cli/test_update.py

@@ -31,6 +31,7 @@
     "LOGNORMAL": ["MEAN", "STD"],
     "TRUNCATED_NORMAL": ["MEAN", "STD", "MIN", "MAX"],
     "TRIANGULAR": ["MIN", "MODE", "MAX"],
+    "PERT": ["MIN", "MODE", "MAX", "SCALE"],
     "UNIFORM": ["MIN", "MAX"],
     "DUNIF": ["STEPS", "MIN", "MAX"],
     "ERRF": ["MIN", "MAX", "SKEWNESS", "WIDTH"],

tests/ert/unit_tests/config/test_gen_kw_config.py

@@ -142,8 +142,25 @@ def test_short_definition_raises_config_error(tmp_path):
             {"name": "KEY10", "distribution": {"name": "const", "value": 10}},
             {"key": "KEY10", "function": "CONST", "parameters": {"VALUE": 10}},
         ),
+        (
+            {
+                "name": "KEY11",
+                "distribution": {
+                    "name": "pert",
+                    "min": 0,
+                    "mode": 1,
+                    "max": 2,
+                    "scale": 4,
+                },
+            },
+            {
+                "key": "KEY11",
+                "function": "PERT",
+                "parameters": {"MIN": 0, "MODE": 1, "MAX": 2, "SCALE": 4},
+            },
+        ),
     ],
-    ids=[f"KEY{i}" for i in range(1, 11)],
+    ids=[f"KEY{i}" for i in range(1, 12)],
 )
 def test_gen_kw_config_get_priors(spec, expected):
     cfg = GenKwConfig(**spec)
@@ -187,6 +204,7 @@ def test_gen_kw_config_get_priors(spec, expected):
         ("ERRF 1 2 0.1 0.1", False, r"KW_NAME:MY_KEYWORD " + number_regex),
         ("DERRF 10 1 2 0.1 0.1", False, r"KW_NAME:MY_KEYWORD " + number_regex),
         ("TRIANGULAR 0 0.5 1", False, r"KW_NAME:MY_KEYWORD " + number_regex),
+        ("PERT 0 0.5 1", False, r"KW_NAME:MY_KEYWORD " + number_regex),
     ],
 )
 async def test_gen_kw_is_log_or_not(
@@ -403,6 +421,74 @@ def test_gen_kw_params_parsing(tmpdir, params, error):
             )
 
 
+@pytest.mark.parametrize(
+    ("values", "expected_scale"),
+    [
+        (["0", "0.5", "1"], 4.0),
+        (["0", "0.5", "1", "2"], 2.0),
+    ],
+)
+def test_that_pert_uses_default_or_explicit_scale(values, expected_scale):
+    distribution = GenKwConfig._parse_distribution("MYNAME", "PERT", values)
+
+    assert distribution.model_dump() == {
+        "name": "pert",
+        "min": 0.0,
+        "mode": 0.5,
+        "max": 1.0,
+        "scale": expected_scale,
+    }
+
+
+@pytest.mark.parametrize(
+    "values",
+    [
+        ["0", "0.5"],
+        ["0", "0.5", "1", "4", "5"],
+    ],
+)
+def test_that_pert_rejects_parameter_counts_other_than_three_or_four(values):
+    with pytest.raises(ConfigValidationError, match="Incorrect number of values"):
+        GenKwConfig._parse_distribution("MYNAME", "PERT", values)
+
+
+def test_that_pert_requires_minimum_strictly_less_than_maximum():
+    with pytest.raises(
+        ConfigValidationError,
+        match=r"Minimum .* must be strictly less than the maximum",
+    ):
+        GenKwConfig._parse_distribution("MYNAME", "PERT", ["1", "1", "1"])
+
+
+@pytest.mark.parametrize("mode", ["-1", "0", "1", "2"])
+def test_that_pert_requires_mode_strictly_between_bounds(mode):
+    with pytest.raises(ConfigValidationError, match="must be strictly between"):
+        GenKwConfig._parse_distribution("MYNAME", "PERT", ["0", mode, "1"])
+
+
+@pytest.mark.parametrize(
+    "parameter_index",
+    range(4),
+    ids=["minimum", "mode", "maximum", "scale"],
+)
+@pytest.mark.parametrize("value", ["nan", "inf", "-inf"])
+def test_that_pert_requires_finite_parameters(parameter_index, value):
+    values = ["0", "0.5", "1", "4"]
+    values[parameter_index] = value
+
+    with pytest.raises(ConfigValidationError, match="finite"):
+        GenKwConfig._parse_distribution("MYNAME", "PERT", values)
+
+
+@pytest.mark.parametrize("scale", ["0", "-1"])
+def test_that_pert_requires_positive_scale(scale):
+    with pytest.raises(
+        ConfigValidationError,
+        match=r"strictly greater than 0",
+    ):
+        GenKwConfig._parse_distribution("MYNAME", "PERT", ["0", "0.5", "1", scale])
+
+
 @pytest.mark.parametrize(
     ("params", "xinput", "expected"),
     [

tests/ert/unit_tests/config/test_transfer_functions.py

@@ -4,6 +4,7 @@
 import pytest
 from hypothesis import assume, given
 from hypothesis import strategies as st
+from scipy.stats import beta as beta_distribution
 from scipy.stats import norm
 
 from ert.config.distribution import TransSettingsValidation, get_distribution
@@ -245,3 +246,64 @@ def test_that_triangular_is_monotonic(args):
             assert y1 >= y2
         else:
             assert y1 <= y2
+
+
+def valid_pert_params():
+    return nice_floats(min_value=-1e6, max_value=1e6).flatmap(
+        lambda mode: st.tuples(
+            st.floats(
+                min_value=mode - 2,
+                max_value=mode - 1,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            st.just(mode),
+            st.floats(
+                min_value=mode + 1,
+                max_value=mode + 2,
+                allow_nan=False,
+                allow_infinity=False,
+            ).filter(lambda maximum: maximum > mode),
+            nice_floats(min_value=0.1, max_value=20),
+        )
+    )
+
+
+@given(nice_floats(), valid_pert_params())
+def test_that_pert_stays_within_bounds(x, args):
+    minimum, _, maximum, _ = args
+    dist = get_distribution("PERT", args)
+
+    assert minimum <= transform_scalar(dist, x) <= maximum
+
+
+def test_that_pert_transform_does_not_exceed_max_from_roundoff():
+    dist = get_distribution("PERT", [-1.003, 0, 1, 1])
+
+    assert transform_scalar(dist, 9) == 1
+
+
+@given(
+    st.tuples(nice_floats(), nice_floats()).map(sorted),
+    valid_pert_params(),
+)
+def test_that_pert_is_non_strictly_monotonic(x_values, args):
+    x1, x2 = x_values
+    dist = get_distribution("PERT", args)
+
+    assert transform_scalar(dist, x1) <= transform_scalar(dist, x2)
+
+
+@given(
+    nice_floats(min_value=-8, max_value=8),
+    valid_pert_params(),
+)
+def test_that_pert_matches_scaled_beta_quantiles(x, args):
+    minimum, mode, maximum, scale = args
+    span = maximum - minimum
+    alpha = 1 + scale * (mode - minimum) / span
+    beta = 1 + scale * (maximum - mode) / span
+    expected = beta_distribution.ppf(norm.cdf(x), alpha, beta, loc=minimum, scale=span)
+    dist = get_distribution("PERT", args)
+
+    assert np.isclose(transform_scalar(dist, x), expected)