diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 80dbdd46..c7310eed 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,10 +1,24 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit - # Ruff version. - rev: v0.1.8 + rev: v0.15.12 hooks: - # Run the linter. + # Scope matches CI's `ruff check src`. - id: ruff args: [--fix] - # Run the formatter. + files: ^src/ - id: ruff-format + files: ^src/ + + - repo: local + hooks: + - id: pytest + name: pytest (full suite, with coverage) + # CI also passes --cov-report=xml --junitxml ... -o junit_family=legacy; + # those exist only for Codecov upload and don't affect pass/fail. + # language: system uses the active venv's pytest so the editable + # install and dev extras resolve. + entry: pytest --cov --cov-branch + language: system + pass_filenames: false + always_run: true + stages: [pre-push] diff --git a/examples/plots/darestani_failure_vs_age_by_class.png b/examples/plots/darestani_failure_vs_age_by_class.png new file mode 100644 index 00000000..8ba7d164 Binary files /dev/null and b/examples/plots/darestani_failure_vs_age_by_class.png differ diff --git a/examples/plots/darestani_failure_vs_conductor_area_by_class.png b/examples/plots/darestani_failure_vs_conductor_area_by_class.png new file mode 100644 index 00000000..63081198 Binary files /dev/null and b/examples/plots/darestani_failure_vs_conductor_area_by_class.png differ diff --git a/examples/plots/darestani_fragility_curves.png b/examples/plots/darestani_fragility_curves.png new file mode 100644 index 00000000..2cd61607 Binary files /dev/null and b/examples/plots/darestani_fragility_curves.png differ diff --git a/examples/plots/darestani_mu_vs_age_by_class.png b/examples/plots/darestani_mu_vs_age_by_class.png new file mode 100644 index 00000000..eadbad1c Binary files /dev/null and b/examples/plots/darestani_mu_vs_age_by_class.png differ diff --git a/examples/plots/darestani_sigma_vs_age_by_class.png b/examples/plots/darestani_sigma_vs_age_by_class.png new file mode 100644 index 00000000..e2a7d64a Binary files /dev/null and b/examples/plots/darestani_sigma_vs_age_by_class.png differ diff --git a/examples/plots/plot_darestani_fragility.py b/examples/plots/plot_darestani_fragility.py new file mode 100644 index 00000000..4fbf0fd3 --- /dev/null +++ b/examples/plots/plot_darestani_fragility.py @@ -0,0 +1,268 @@ +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +matplotlib.rcParams.update({ + "font.size": 14, + "axes.titlesize": 16, + "axes.labelsize": 14, + "xtick.labelsize": 12, + "ytick.labelsize": 12, + "legend.fontsize": 12, + "legend.title_fontsize": 13, +}) + +from erad.models.custom_distributions import Darestani2019 +from erad.enums import PoleClass, PoleConstructionMaterial +from erad.quantities import Speed, WindAngle, ConductorArea, PoleAge +from erad.models.asset import DistributionPole +from erad.probability_builder import ProbabilityFunctionBuilder + + +# Baseline inputs. +POLE_MATERIAL = PoleConstructionMaterial.WOOD +WIND_ANGLE = WindAngle(90, "degree") +CONDUCTOR_AREA = ConductorArea(2, "m^2") +POLE_AGE = PoleAge(50, "year") + +# Sensitivity settings. +AGE_VALUES = list(range(0, 101, 5)) # years +CONDUCTOR_AREA_VALUES = [i / 2 for i in range(0, 17)] # m^2 from 0 to 8 +SUMMARY_WIND_SPEED_MPH = 80 + +WIND_SPEEDS_MPH = list(range(0, 250, 2)) + + +def _failure_probability_curve(asset: DistributionPole) -> list[float]: + dist_instance = Darestani2019(asset) + prob_builder = ProbabilityFunctionBuilder("Darestani2019", [Speed(0, "m/s")], dist_instance) + return [prob_builder.probability(Speed(ws, "mph")) for ws in WIND_SPEEDS_MPH] + + +def _failure_probability_at_speed(asset: DistributionPole, wind_speed_mph: float) -> float: + dist_instance = Darestani2019(asset) + prob_builder = ProbabilityFunctionBuilder("Darestani2019", [Speed(0, "m/s")], dist_instance) + return prob_builder.probability(Speed(wind_speed_mph, "mph")) + + +def _darestani_params(asset: DistributionPole) -> tuple[float, float]: + dist_instance = Darestani2019(asset) + return dist_instance.mu, dist_instance.sigma + + +def _add_fixed_params_label(ax, text: str) -> None: + ax.text( + 0.02, + 0.98, + text, + transform=ax.transAxes, + verticalalignment="top", + fontsize=10, + bbox={"boxstyle": "round,pad=0.3", "facecolor": "white", "alpha": 0.85}, + ) + + +def _wind_speed_at_probability( + wind_speeds: list[float], probabilities: list[float], target_probability: float = 0.5 +) -> float | None: + if not wind_speeds or not probabilities or len(wind_speeds) != len(probabilities): + return None + + if target_probability < probabilities[0] or target_probability > probabilities[-1]: + return None + + for i in range(1, len(probabilities)): + p_prev, p_curr = probabilities[i - 1], probabilities[i] + if p_prev <= target_probability <= p_curr: + x_prev, x_curr = wind_speeds[i - 1], wind_speeds[i] + if p_curr == p_prev: + return x_curr + frac = (target_probability - p_prev) / (p_curr - p_prev) + return x_prev + frac * (x_curr - x_prev) + + return None + + +def _plot_baseline_by_pole_class(asset: DistributionPole, out_dir: Path) -> None: + fig, ax1 = plt.subplots(figsize=(10, 6)) + class_labels: list[str] = [] + class_indices: list[int] = [] + ws50_values: list[float] = [] + + for idx, pole_class in enumerate(PoleClass, start=1): + asset.pole_class = pole_class + failure_probability = _failure_probability_curve(asset) + ax1.plot(WIND_SPEEDS_MPH, failure_probability, label=pole_class.name, linewidth=2) + + ws50 = _wind_speed_at_probability(WIND_SPEEDS_MPH, failure_probability, 0.5) + if ws50 is not None: + class_labels.append(pole_class.name) + class_indices.append(idx) + ws50_values.append(ws50) + + ax1.set_xlabel("Wind Speed (mph)") + ax1.set_ylabel("Failure Probability") + ax1.set_title("Darestani2019 Fragility Curves by Pole Class") + ax1.grid(True, alpha=0.3) + ax1.legend(title="Pole Class", loc="lower left") + + ax2 = ax1.twinx() + ax2.set_ylabel("Pole Class") + if ws50_values: + ax2.plot( + ws50_values, + class_indices, + color="black", + linestyle="--", + marker="o", + linewidth=2, + label="Wind Speed at 50% Failure", + ) + ax2.set_yticks(class_indices) + ax2.set_yticklabels(class_labels) + ax2.legend(loc="lower right") + + fixed_text = ( + f"Fixed: Material={asset.pole_material.name}, " + f"Wind Angle={asset.wind_angle.magnitude:g} {asset.wind_angle.units}, " + f"Conductor Area={asset.conductor_area.magnitude:g} {asset.conductor_area.units}, " + f"Pole Age={asset.pole_age.magnitude:g} {asset.pole_age.units}" + ) + _add_fixed_params_label(ax1, fixed_text) + + fig.tight_layout() + out_path = out_dir / "darestani_fragility_curves.png" + fig.savefig(out_path, dpi=180) + plt.close(fig) + print(out_path) + + +def _plot_mu_vs_age_by_class(asset: DistributionPole, out_dir: Path) -> None: + plt.figure(figsize=(10, 6)) + + import math + + for pole_class in PoleClass: + asset.pole_class = pole_class + exp_mu_values = [] + for age in AGE_VALUES: + asset.pole_age = PoleAge(age, "year") + mu, _ = _darestani_params(asset) + exp_mu_values.append(math.exp(mu)) + plt.plot(AGE_VALUES, exp_mu_values, label=pole_class.name, linewidth=2) + + plt.xlabel("Pole Age (years)") + plt.ylabel("exp(μ) (mph)") + plt.title("Median Failure Wind Speed vs. Pole Age by Pole Class") + plt.legend(title="Pole Class") + plt.grid(True, alpha=0.3) + + fixed_text = ( + f"Fixed: Material={asset.pole_material.name}, " + f"Wind Angle={asset.wind_angle.magnitude:g} {asset.wind_angle.units}, " + f"Conductor Area={asset.conductor_area.magnitude:g} {asset.conductor_area.units}" + ) + _add_fixed_params_label(plt.gca(), fixed_text) + + plt.tight_layout() + out_path = out_dir / "darestani_mu_vs_age_by_class.png" + plt.savefig(out_path, dpi=180) + plt.close() + print(out_path) + + +def _plot_sigma_vs_age_by_class(asset: DistributionPole, out_dir: Path) -> None: + plt.figure(figsize=(10, 6)) + + for pole_class in PoleClass: + asset.pole_class = pole_class + sigma_values = [] + for age in AGE_VALUES: + asset.pole_age = PoleAge(age, "year") + _, sigma = _darestani_params(asset) + sigma_values.append(sigma) + plt.plot(AGE_VALUES, sigma_values, label=pole_class.name, linewidth=2) + + plt.xlabel("Pole Age (years)") + plt.ylabel("σ") + plt.ylim(0, 0.8) + plt.title("Failure Wind Speed Variability (σ) vs. Pole Age") + plt.legend(title="Pole Class") + plt.grid(True, alpha=0.3) + + fixed_text = ( + f"Fixed: Material={asset.pole_material.name}, " + f"Wind Angle={asset.wind_angle.magnitude:g} {asset.wind_angle.units}, " + f"Conductor Area={asset.conductor_area.magnitude:g} {asset.conductor_area.units}" + ) + _add_fixed_params_label(plt.gca(), fixed_text) + + plt.tight_layout() + out_path = out_dir / "darestani_sigma_vs_age_by_class.png" + plt.savefig(out_path, dpi=180) + plt.close() + print(out_path) + + +def _plot_failure_vs_conductor_area_by_class(asset: DistributionPole, out_dir: Path) -> None: + plt.figure(figsize=(10, 6)) + + for pole_class in PoleClass: + asset.pole_class = pole_class + y_values = [] + for area in CONDUCTOR_AREA_VALUES: + asset.conductor_area = ConductorArea(area, "m^2") + y_values.append(_failure_probability_at_speed(asset, SUMMARY_WIND_SPEED_MPH)) + plt.plot(CONDUCTOR_AREA_VALUES, y_values, label=pole_class.name, linewidth=2) + + plt.xlabel("Conductor Area (m^2)") + plt.ylabel("Failure Probability") + plt.title(f"Failure Probability vs Conductor Area at {SUMMARY_WIND_SPEED_MPH} mph") + plt.legend(title="Pole Class") + plt.grid(True, alpha=0.3) + + fixed_text = ( + f"Fixed: Material={asset.pole_material.name}, " + f"Wind Angle={asset.wind_angle.magnitude:g} {asset.wind_angle.units}, " + f"Pole Age={asset.pole_age.magnitude:g} {asset.pole_age.units}, " + f"Wind Speed={SUMMARY_WIND_SPEED_MPH} mph" + ) + _add_fixed_params_label(plt.gca(), fixed_text) + + plt.tight_layout() + out_path = out_dir / "darestani_failure_vs_conductor_area_by_class.png" + plt.savefig(out_path, dpi=180) + plt.close() + print(out_path) + + +def main() -> None: + asset = DistributionPole.example() + asset.pole_material = POLE_MATERIAL + asset.wind_angle = WIND_ANGLE + asset.conductor_area = CONDUCTOR_AREA + asset.pole_age = POLE_AGE + + out_dir = Path("examples/plots") + out_dir.mkdir(parents=True, exist_ok=True) + + _plot_baseline_by_pole_class(asset, out_dir) + + asset.conductor_area = CONDUCTOR_AREA + asset.pole_age = POLE_AGE + _plot_mu_vs_age_by_class(asset, out_dir) + + asset.conductor_area = CONDUCTOR_AREA + asset.pole_age = POLE_AGE + _plot_sigma_vs_age_by_class(asset, out_dir) + + asset.conductor_area = CONDUCTOR_AREA + asset.pole_age = POLE_AGE + _plot_failure_vs_conductor_area_by_class(asset, out_dir) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 14efd6a9..5b2bb119 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,8 +6,8 @@ build-backend = "setuptools.build_meta" name = "NREL-erad" dynamic = ["version"] description = "Graph based scalable tool for computing energy resilience metrics for distribution systems." -readme = "README.md" requires-python = ">=3.10" +readme = "README.md" authors = [ { name = "Aadil Latif", email = "aadil.altif@nrel.gov" }, { name = "Kapil Duwadi", email = "kapil.duwadi@nrel.gov" }, @@ -67,6 +67,7 @@ dev = [ "mkdocs-material", "mkdocstrings[python]", "pylint", + "pre-commit", "pytest", "pytest-asyncio", "pytest-cov", @@ -84,6 +85,8 @@ doc = [ [project.urls] Homepage = "https://github.com/NLR-Distribution-Suite/erad" +[tool.poetry] +packages = [{ include = "erad", from = "src" }] [project.scripts] erad = "erad.cli:main" erad-mcp = "erad.mcp:main" @@ -100,7 +103,6 @@ exclude = [ "env", "venv", ] -ignore-init-module-imports = true line-length = 99 indent-width = 4 target-version = "py311" diff --git a/src/erad/data/pole_coefficients_data.py b/src/erad/data/pole_coefficients_data.py new file mode 100644 index 00000000..83b451de --- /dev/null +++ b/src/erad/data/pole_coefficients_data.py @@ -0,0 +1,176 @@ +from erad.enums import PoleClass, PoleConstructionMaterial + +""" This module contains pole coefficients data from Darestani et al. (2019) Tables 14-20""" +Darestani2019_PoleCoefficients = { + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_1): { + "mu_coefficients": [ + 5.581, + -4.306e-3, + -2.408e-2, + 5.986e-3, + 2.569e-5, + -1.231e-3, + 3.524e-3, + 0, + -1.324e-5, + -1.327e-4, + ], + "sigma_coefficients": [ + 2.974e-01, + -1.411e-03, + -9.462e-03, + -2.678e-03, + 7.782e-06, + -5.088e-05, + 6.199e-04, + 6.805e-06, + 5.630e-05, + 6.037e-05, + ], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_2): { + "mu_coefficients": [ + 5.770, + -6.730e-03, + -3.904e-02, + 6.090e-03, + 4.342e-05, + -1.383e-03, + 4.972e-03, + 0, + 0, + -1.328e-04, + ], + "sigma_coefficients": [ + 2.080e-01, + -5.938e-04, + -3.661e-03, + -2.254e-03, + 3.320e-06, + 0, + 2.357e-04, + 3.039e-06, + 1.932e-05, + 6.262e-05, + ], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_3): { + "mu_coefficients": [ + 5.874, + -8.946e-03, + -5.238e-02, + 6.291e-03, + 5.994e-05, + -1.468e-03, + 6.167e-03, + 0, + 0, + -1.342e-04, + ], + "sigma_coefficients": [ + 1.698e-01, + -1.165e-04, + 0, + -2.016e-03, + 0, + 0, + 0, + 9.217e-07, + 0, + 6.326e-05, + ], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_4): { + "mu_coefficients": [ + 5.930e0, + -1.091e-02, + -6.420e-02, + 6.387e-03, + 7.453e-05, + -1.521e-03, + 7.169e-03, + 0, + 0, + -1.352e-04, + ], + "sigma_coefficients": [1.568e-01, 0, 0, -1.990e-03, 0, 0, 0, 0, 0, 6.428e-05], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_5): { + "mu_coefficients": [ + 6.029, + -1.378e-02, + -8.249e-02, + 6.381e-03, + 9.546e-05, + -1.561e-03, + 8.640e-03, + 0, + 0, + -1.358e-04, + ], + "sigma_coefficients": [ + 1.743e-01, + -1.431e-04, + 0, + -2.042e-03, + 1.085e-06, + 0, + 0, + 0, + 0, + 6.477e-05, + ], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_6): { + "mu_coefficients": [ + 6.111, + -1.683e-02, + -9.992e-02, + 6.409e-03, + 1.182e-04, + -1.597e-03, + 1.002e-02, + 0, + 0, + -1.365e-04, + ], + "sigma_coefficients": [ + 1.535e-01, + 1.930e-04, + 2.421e-03, + -1.790e-03, + 0, + 0, + 0, + -3.331e-06, + -3.473e-05, + 6.519e-05, + ], + }, + (PoleConstructionMaterial.WOOD, PoleClass.CLASS_7): { + "mu_coefficients": [ + 6.130, + -1.873e-02, + -1.130e-01, + 6.291e-03, + 1.314e-04, + -1.595e-03, + 1.101e-02, + 0, + 0, + -1.371e-04, + ], + "sigma_coefficients": [ + 1.460e-01, + 3.313e-04, + 2.857e-03, + -1.609e-03, + 0, + -1.538e-05, + 0, + -6.106e-06, + -5.067e-05, + 6.576e-05, + ], + }, +} diff --git a/src/erad/enums.py b/src/erad/enums.py index 2b4eed89..8bf80e90 100644 --- a/src/erad/enums.py +++ b/src/erad/enums.py @@ -1,4 +1,4 @@ -from enum import IntEnum +from enum import IntEnum, StrEnum class ScenarioTypes(IntEnum): @@ -38,3 +38,18 @@ def has_value(cls, value): @classmethod def has_asset(cls, asset): return asset in cls.__members__ + + +class PoleClass(StrEnum): + CLASS_1 = "Class 1" + CLASS_2 = "Class 2" + CLASS_3 = "Class 3" + CLASS_4 = "Class 4" + CLASS_5 = "Class 5" + CLASS_6 = "Class 6" + CLASS_7 = "Class 7" + + +class PoleConstructionMaterial(StrEnum): + WOOD = "wood" or "Wood" or "WOOD" + STEEL = "steel" or "Steel" or "STEEL" diff --git a/src/erad/gdm_mapping.py b/src/erad/gdm_mapping.py index d2a1f5dc..09579dcf 100644 --- a/src/erad/gdm_mapping.py +++ b/src/erad/gdm_mapping.py @@ -34,50 +34,51 @@ AssetTypes.distribution_underground_cables: [ ComponentFilterModel( component_type=gdc.GeometryBranch, - component_filter=lambda x: x.equipment.conductors[0].__class__ - == gde.ConcentricCableEquipment - and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0, + component_filter=lambda x: ( + x.equipment.conductors[0].__class__ == gde.ConcentricCableEquipment + and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0 + ), ), ComponentFilterModel( component_type=gdc.MatrixImpedanceBranch, - component_filter=lambda x: x.equipment.c_matrix[0, 0] - .to("microfarad/kilometer") - .magnitude - > 0.05 - and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0, + component_filter=lambda x: ( + x.equipment.c_matrix[0, 0].to("microfarad/kilometer").magnitude > 0.05 + and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0 + ), ), ], AssetTypes.distribution_overhead_lines: [ ComponentFilterModel( component_type=gdc.GeometryBranch, component_filter=lambda x: ( - x.equipment.conductors[0].__class__ == gde.BareConductorEquipment - ) - and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0, + (x.equipment.conductors[0].__class__ == gde.BareConductorEquipment) + and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0 + ), ), ComponentFilterModel( component_type=gdc.MatrixImpedanceBranch, - component_filter=lambda x: x.equipment.c_matrix[0, 0] - .to("microfarad/kilometer") - .magnitude - < 0.05 - and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0, + component_filter=lambda x: ( + x.equipment.c_matrix[0, 0].to("microfarad/kilometer").magnitude < 0.05 + and x.buses[0].rated_voltage.to("kilovolt").magnitude < 35.0 + ), ), ], AssetTypes.transmission_underground_cables: [ ComponentFilterModel( component_type=gdc.GeometryBranch, - component_filter=lambda x: x.equipment.conductors[0].__class__ - == gde.ConcentricCableEquipment - and x.buses[0].rated_voltage.to("kilovolt").magnitude > 35.0, + component_filter=lambda x: ( + x.equipment.conductors[0].__class__ == gde.ConcentricCableEquipment + and x.buses[0].rated_voltage.to("kilovolt").magnitude > 35.0 + ), ), ], AssetTypes.transmission_overhead_lines: [ ComponentFilterModel( component_type=gdc.GeometryBranch, - component_filter=lambda x: x.equipment.conductors[0].__class__ - == gde.BareConductorEquipment - and x.buses[0].rated_voltage.to("kilovolt").magnitude > 35.0, + component_filter=lambda x: ( + x.equipment.conductors[0].__class__ == gde.BareConductorEquipment + and x.buses[0].rated_voltage.to("kilovolt").magnitude > 35.0 + ), ), ], } diff --git a/src/erad/models/asset.py b/src/erad/models/asset.py index 12cb6bc2..2a8bf553 100644 --- a/src/erad/models/asset.py +++ b/src/erad/models/asset.py @@ -1,16 +1,18 @@ from datetime import datetime from uuid import UUID import math - -from pydantic import computed_field, Field +from pydantic import computed_field, Field, field_validator from infrasys.quantities import Distance from geopy.distance import geodesic from shapely.geometry import Point from pyhigh import get_elevation from infrasys import Component from loguru import logger +from typing import Literal +from types import SimpleNamespace -# from erad.constants import RASTER_DOWNLOAD_PATH +from erad.probability_builder import ProbabilityFunctionBuilder +from erad.models.custom_distributions import CUSTOM_DISTRIBUTIONS from erad.quantities import Acceleration, Speed from erad.models.probability import ( AccelerationProbability, @@ -19,7 +21,8 @@ RatioProbability, SpeedProbability, ) -from erad.enums import AssetTypes +from erad.enums import AssetTypes, PoleClass, PoleConstructionMaterial +from erad.quantities import WindAngle, ConductorArea, PoleAge import erad.models.hazard as hz @@ -290,7 +293,7 @@ def update_survival_probability( asset_location, hazard_model, self.elevation + self.height ) else: - raise (f"Unsupported hazard type {hazard_model.__class__.__name__}") + raise Exception(f"Unsupported hazard type {hazard_model.__class__.__name__}") self.calculate_probabilities(asset_state, frag_curves) if asset_state not in self.asset_state: @@ -306,7 +309,7 @@ def calculate_probabilities(self, asset_state: AssetState, frag_curves): for field in fields: prob_model = getattr(asset_state, field) if prob_model and prob_model.survival_probability == 1: - curve = self.get_vadid_curve(frag_curves, field) + curve = self.get_valid_curve(frag_curves, field) if curve is None: raise Exception( f"No fragility curve found for field - {field} for asset type {self.asset_type.name}" @@ -320,13 +323,12 @@ def calculate_probabilities(self, asset_state: AssetState, frag_curves): quantity = getattr(prob_model, quantity_name) prob_model.survival_probability = 1 - prob_inst.probability(quantity) - def get_vadid_curve(self, frag_curves, field: SyntaxError) -> float: + def get_valid_curve(self, frag_curves, field: str): for frag_curve in frag_curves: if frag_curve.asset_state_param == field: for curve in frag_curve.curves: if curve.asset_type == self.asset_type: return curve.prob_function - return None @classmethod @@ -340,3 +342,58 @@ def example(cls) -> "Asset": longitude=-122.4194, asset_state=[AssetState.example()], ) + + +class DistributionPole(Asset): + pole_class: PoleClass | None = Field( + description="Class of the pole (only for distribution poles)" + ) + pole_material: PoleConstructionMaterial | None = Field( + description="Construction material of the pole (only for distribution poles)" + ) + wind_angle: WindAngle | None = Field( + description="Angle between wind direction and pole" + ) # Hazard property? + conductor_area: ConductorArea | None = Field(description="Conductor area in m^2") + pole_age: PoleAge | None = Field(description="Age of the pole in years") + probability_dist: Literal[*CUSTOM_DISTRIBUTIONS] | None = Field( + description="Custom distribution function to model pole failure probability. Run list_custom_distributions() for available options." + ) + + @field_validator("asset_type") + @classmethod + def validate_asset_type(cls, v): + if v != AssetTypes.distribution_poles: + raise ValueError("Invalid asset type for DistributionPole, must be distribution_poles") + return v + + def get_valid_curve(self, frag_curves, field: str): + if self.probability_dist and field == "wind_speed": + custom_dist_instance = CUSTOM_DISTRIBUTIONS[self.probability_dist](self) + prob_builder = ProbabilityFunctionBuilder( + dist=self.probability_dist, + params=[Speed(0, "m/s")], # Only to set output type + custom_dist_instance=custom_dist_instance, + ) + wrapper = SimpleNamespace() + wrapper.prob_model = prob_builder + return wrapper + return super().get_valid_curve(frag_curves, field) + + @classmethod + def example(cls) -> "DistributionPole": + return DistributionPole( + name="Asset 1", + asset_type=AssetTypes.distribution_poles, + distribution_asset=UUID("123e4567-e89b-12d3-a456-426614174000"), + height=Distance(100, "m"), + latitude=37.7749, + longitude=-122.4194, + asset_state=[AssetState.example()], + pole_class=PoleClass.CLASS_5, + pole_material=PoleConstructionMaterial.WOOD, + wind_angle=WindAngle(90, "degree"), + conductor_area=ConductorArea(0.0005, "m**2"), + pole_age=PoleAge(30, "year"), + probability_dist="Darestani2019", + ) diff --git a/src/erad/models/custom_distributions.py b/src/erad/models/custom_distributions.py new file mode 100644 index 00000000..b0717cb0 --- /dev/null +++ b/src/erad/models/custom_distributions.py @@ -0,0 +1,94 @@ +import numpy as np +from scipy.stats import norm, rv_continuous +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from erad.models.asset import DistributionPole + +from erad.data import pole_coefficients_data + + +def validate_range(param: float, min_value: float, max_value: float, param_name: str): + """Common validation function to check if a parameter is within a specified range.""" + if not (min_value <= param <= max_value): + raise ValueError(f"{param_name} out of range [{min_value}, {max_value}]") + + +class Darestani2019(rv_continuous): + """Custom distribution based on Darestani et al. (2019) model.""" + + def __init__(self, asset: "DistributionPole"): + super().__init__(name="Darestani2019") # Initialize the base class + for attr in ["wind_angle", "conductor_area", "pole_age", "pole_class", "pole_material"]: + if getattr(asset, attr) is None: + raise ValueError(f"{attr} must be provided for Darestani2019 distribution") + + self.wind_angle = asset.wind_angle + self.conductor_area = asset.conductor_area + self.pole_age = asset.pole_age + self.pole_class = asset.pole_class + self.pole_material = asset.pole_material + + validate_range(self.wind_angle.magnitude, 0, 90, "Wind angle") + validate_range(self.conductor_area.magnitude, 0, 8, "Conductor area") + validate_range(self.pole_age.magnitude, 0, 100, "Pole age") + + self.mu, self.sigma = self._calculate_params() + + def _cdf(self, x): + """Cumulative distribution function.""" + if x < 2: + return 0 + if x > 112: + return 1 + log_x = np.log(2.23694 * x) + return norm.cdf((log_x - self.mu) / self.sigma) + + def _rvs(self, size=1, random_state=None): + """Random variates.""" + rng = np.random.default_rng(random_state) + log_transformed_speed = rng.normal(self.mu, self.sigma, size) + sampled_wind_speed = np.exp(log_transformed_speed) / 2.23694 + sampled_wind_speed = np.clip(sampled_wind_speed, 2, 112) + return sampled_wind_speed + + def _calculate_params(self) -> tuple[float, float]: + """Calculate mu and sigma params for Darestani et al. (2019) distribution.""" + + key = (self.pole_material, self.pole_class) + if key not in pole_coefficients_data.Darestani2019_PoleCoefficients: + raise ValueError("No coefficients found") + + coeff_data = pole_coefficients_data.Darestani2019_PoleCoefficients[key] + + params = [ + 1, + self.wind_angle.magnitude, + self.conductor_area.magnitude, + self.pole_age.magnitude, + self.wind_angle.magnitude**2, + self.wind_angle.magnitude * self.conductor_area.magnitude, + self.conductor_area.magnitude**2, + self.wind_angle.magnitude * self.pole_age.magnitude, + self.conductor_area.magnitude * self.pole_age.magnitude, + self.pole_age.magnitude**2, + ] + + mu = np.dot(coeff_data["mu_coefficients"], params) + sigma = np.dot(coeff_data["sigma_coefficients"], params) + + if sigma <= 0: + raise ValueError("Sigma must be positive") + + return mu, sigma + + +CUSTOM_DISTRIBUTIONS = { + "Darestani2019": Darestani2019, # Add other custom distributions here as needed +} + + +# Helper functions +def list_custom_distributions(): + """List available custom functions.""" + return list(CUSTOM_DISTRIBUTIONS.keys()) diff --git a/src/erad/models/fragility_curve.py b/src/erad/models/fragility_curve.py index 3d4a1170..2eb80270 100644 --- a/src/erad/models/fragility_curve.py +++ b/src/erad/models/fragility_curve.py @@ -2,7 +2,6 @@ from typing import Literal from pathlib import Path - from infrasys import Component, BaseQuantity from scipy.stats import _continuous_distns from pydantic import field_validator @@ -13,6 +12,7 @@ from erad.models.asset import AssetState from erad.enums import AssetTypes from erad.quantities import Speed +from erad.models.custom_distributions import CUSTOM_DISTRIBUTIONS FRAGILITY_CURVE_TYPES = [ @@ -23,17 +23,26 @@ class ProbabilityFunction(Component): name: str = "" - distribution: Literal[*SUPPORTED_CONT_DIST] + distribution: str parameters: list[float | BaseQuantity] + @field_validator("distribution") + def validate_distribution(cls, value): + if isinstance(value, str): + if value not in SUPPORTED_CONT_DIST and value not in CUSTOM_DISTRIBUTIONS: + raise ValueError(f"Unsupported distribution: {value}") + return value + @field_validator("parameters") - def validate_parameters(cls, value): + def validate_parameters(cls, value, info): if not any(isinstance(v, BaseQuantity) for v in value): raise ValueError("There should be atleast one BaseQuantity in the parameters") - units = set([v.units for v in value if isinstance(v, BaseQuantity)]) - if not len(units) == 1: - raise ValueError("All BaseQuantities should have the same units") + distribution = info.data.get("distribution") + if distribution in SUPPORTED_CONT_DIST: + units = set([v.units for v in value if isinstance(v, BaseQuantity)]) + if len(units) > 1: + raise ValueError("All BaseQuantities should have the same units") return value @cached_property @@ -42,6 +51,7 @@ def prob_model(self) -> ProbabilityFunctionBuilder: @classmethod def example(cls) -> "ProbabilityFunction": + """Example for a scipy distribution.""" return ProbabilityFunction( distribution="norm", parameters=[Speed(1.5, "m/s"), 2], @@ -55,6 +65,7 @@ class FragilityCurve(Component): @classmethod def example(cls) -> "FragilityCurve": + """Example for a fragility curve with a scipy distribution.""" return FragilityCurve( asset_type=AssetTypes.substation, prob_function=ProbabilityFunction.example(), diff --git a/src/erad/plotting.py b/src/erad/plotting.py new file mode 100644 index 00000000..4e86e0b2 --- /dev/null +++ b/src/erad/plotting.py @@ -0,0 +1,73 @@ +"""Plotting helpers for ERAD. Kept separate from model modules so that +plotly is only imported when plotting code is actually used.""" + +from pathlib import Path + +import numpy as np +import plotly.express as px + +from erad.models.asset import Asset + + +def plot_custom_fragility_curves( + assets: list[Asset], + asset_state_param: str, + file_path: Path | None = None, + x_min: float = 0, + x_max: float = 60, + number_of_points: int = 1000, +): + """Plot fragility curves derived from each asset's custom probability + distribution, one line per asset. + + Each asset must have its custom-distribution field set (e.g. + DistributionPole.probability_dist = "Darestani2019") and the + attributes that distribution requires. The curve is obtained via + Asset.get_valid_curve(frag_curves=[], field=asset_state_param), which + routes to the existing asset-side custom-dist path. + """ + if not assets: + raise ValueError("No assets to plot") + if file_path is not None: + file_path = Path(file_path) + assert file_path.suffix.lower() == ".html", "File path should be an HTML file" + + sample_curve = assets[0].get_valid_curve(frag_curves=[], field=asset_state_param) + if sample_curve is None: + raise ValueError( + f"Asset {assets[0].name!r} has no curve for '{asset_state_param}'. " + "Ensure its custom-distribution field (e.g. probability_dist) is set." + ) + quantity_cls = sample_curve.prob_model.quantity + units = sample_curve.prob_model.units + + x = np.linspace(x_min, x_max, number_of_points) + plot_data = {"x": [], "y": [], "Asset": []} + + for asset in assets: + curve = asset.get_valid_curve(frag_curves=[], field=asset_state_param) + if curve is None: + continue + # Iterate scalar-wise: custom rv_continuous subclasses (e.g. Darestani2019) + # implement _cdf with Python conditionals that aren't array-safe. + y = [curve.prob_model.probability(quantity_cls(xi, units)) for xi in x] + label = asset.name or type(asset).__name__ + plot_data["x"].extend(x) + plot_data["y"].extend(y) + plot_data["Asset"].extend([label] * len(x)) + + x_label = asset_state_param.replace("_", " ").title() + fig = px.line( + plot_data, + x="x", + y="y", + color="Asset", + labels={"x": f"{x_label} [{units}]", "y": "Probability of Failure"}, + title=f"Custom Fragility Curves for {x_label}", + ) + + fig.show() + if file_path: + fig.write_html(file_path) + + return fig diff --git a/src/erad/probability_builder.py b/src/erad/probability_builder.py index 5a7165ab..9864465c 100644 --- a/src/erad/probability_builder.py +++ b/src/erad/probability_builder.py @@ -3,36 +3,42 @@ class ProbabilityFunctionBuilder: - """Class containing utility fuctions for sceario definations.""" + """Class containing utility fuctions for scenario definations.""" - def __init__(self, dist, params: list[float | BaseQuantity]): + def __init__(self, dist: str, params: list[float | BaseQuantity], custom_dist_instance=None): """Constructor for BaseScenario class. Args: - dist (str): Name of teh distribution. Should follow Scipy naming convention - params (list): A list of parameters for the chosen distribution function. See Scipy.stats documentation + dist (str or custom distribution instance): + - String name for scipy.stats distributions (See Scipy.stats documentation) + - String name for custom distributions (Call list_custom_distributions() to see available functions) + params (list): Parameters for the distribution or custom function. """ - base_quantity = [p for p in params if isinstance(p, BaseQuantity)][0] - self.quantity = base_quantity.__class__ - self.units = base_quantity.units - self.dist = getattr(stats, dist) + base_quantities = [p for p in params if isinstance(p, BaseQuantity)] + self.quantity = base_quantities[0].__class__ + self.units = base_quantities[0].units self.params = [p.magnitude if isinstance(p, BaseQuantity) else p for p in params] + + if hasattr(stats, dist): + dist_class = getattr(stats, dist) + self.dist = dist_class(*self.params) # Create a frozen scipy distribution instance + elif custom_dist_instance: + self.dist = custom_dist_instance # Use the provided custom distribution instance + else: + raise ValueError(f"Unsupported distribution: {dist}.") + return def sample(self): """Sample the distribution""" - return self.quantity(self.dist.rvs(*self.params, size=1)[0], self.units) + return self.quantity(self.dist.rvs(size=1)[0], self.units) def probability(self, value: BaseQuantity) -> float: """Calculates survival probability of a given asset. Args: - value (float): value for vetor of interest. Will change with scenarions + value (float): value for vector of interest. Will change with scenarios """ assert isinstance(value, BaseQuantity), "Value must be a BaseQuantity" - - cdf = self.dist.cdf - try: - return cdf(value.to(self.units).magnitude, *self.params) - except Exception: - return cdf(value, *self.params) + val = value.to(self.units).magnitude + return self.dist.cdf(val) diff --git a/src/erad/quantities.py b/src/erad/quantities.py index 72b2dec9..816661dd 100644 --- a/src/erad/quantities.py +++ b/src/erad/quantities.py @@ -31,6 +31,22 @@ class Flow(BaseQuantity): __base_unit__ = "feet**3/second" +class WindAngle(BaseQuantity): + """Quantity representing angle between wind and conductor.""" + + __base_unit__ = "degree" + + +class ConductorArea(BaseQuantity): + """Quantity representing cross-sectional area of a conductor.""" + + __base_unit__ = "m**2" + + +class PoleAge(BaseQuantity): + """Quantity representing age of a pole.""" + + __base_unit__ = "year" class Ratio(BaseQuantity): """Quantity representing a dimensionless ratio (e.g. soil saturation fraction).""" diff --git a/src/erad/systems/asset_system.py b/src/erad/systems/asset_system.py index 7f54182f..b1fd1b04 100644 --- a/src/erad/systems/asset_system.py +++ b/src/erad/systems/asset_system.py @@ -29,9 +29,9 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) def add_component(self, component, **kwargs): - assert isinstance( - component, ASSET_TYPES - ), f"Unsupported model type {component.__class__.__name__}" + assert isinstance(component, ASSET_TYPES), ( + f"Unsupported model type {component.__class__.__name__}" + ) return super().add_component(component, **kwargs) def add_components(self, *components, **kwargs): diff --git a/src/erad/systems/hazard_system.py b/src/erad/systems/hazard_system.py index f406db46..7fdd901b 100644 --- a/src/erad/systems/hazard_system.py +++ b/src/erad/systems/hazard_system.py @@ -14,9 +14,9 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) def add_component(self, component, **kwargs): - assert isinstance( - component, HAZARD_MODELS - ), f"Unsupported model type {component.__class__.__name__}" + assert isinstance(component, HAZARD_MODELS), ( + f"Unsupported model type {component.__class__.__name__}" + ) return super().add_component(component, **kwargs) def add_components(self, *components, **kwargs): diff --git a/tests/test_asset_system.py b/tests/test_asset_system.py index 59ef5155..1b601f90 100644 --- a/tests/test_asset_system.py +++ b/tests/test_asset_system.py @@ -1,6 +1,6 @@ from infrasys import Component -from gdm.distribution import DistributionSystem +from gdm.distribution.distribution_system import DistributionSystem import pytest from erad.systems.asset_system import AssetSystem diff --git a/tests/test_custom_distributions.py b/tests/test_custom_distributions.py new file mode 100644 index 00000000..9c530a44 --- /dev/null +++ b/tests/test_custom_distributions.py @@ -0,0 +1,107 @@ +import pytest +from erad.models.custom_distributions import ( + Darestani2019, + CUSTOM_DISTRIBUTIONS, + list_custom_distributions, +) +from erad.models.asset import DistributionPole +from erad.enums import PoleConstructionMaterial, PoleClass +import re +import tempfile +import json + + +@pytest.mark.parametrize( + "distribution_name,distribution_class", + [ + ("Darestani2019", Darestani2019), + ], +) +def test_custom_distributions_registry(distribution_name, distribution_class): + assert distribution_name in CUSTOM_DISTRIBUTIONS + assert CUSTOM_DISTRIBUTIONS[distribution_name] is distribution_class + + +@pytest.mark.parametrize("distribution_name", ["Darestani2019"]) +def test_list_custom_distributions(distribution_name): + available_custom_distributions = list_custom_distributions() + assert distribution_name in available_custom_distributions + assert isinstance(available_custom_distributions, list) + assert ( + len(available_custom_distributions) > 0 + ), "Should have at least one custom distribution listed" + assert all( + isinstance(name, str) for name in available_custom_distributions + ), "All names should be strings" + + +@pytest.mark.parametrize( + "attr,value,err_msg", + [ + ("wind_angle", None, "wind_angle must be provided for Darestani2019 distribution"), + ("pole_age", None, "pole_age must be provided for Darestani2019 distribution"), + ], +) +def test_darestani2019_invalid_single_parameter(attr, value, err_msg): + asset = DistributionPole.example() + setattr(asset, attr, value) + with pytest.raises(ValueError, match=err_msg): + Darestani2019(asset) + + +def test_darestani2019_invalid_material_and_class(): + asset = DistributionPole.example() + asset.pole_material = "steel" + asset.pole_class = "Class 3" + with pytest.raises(ValueError, match="No coefficients found"): + Darestani2019(asset) + + +@pytest.mark.parametrize( + "attr,value,err_msg", + [ + ("wind_angle", -10, "Wind angle out of range [0, 90]"), + ("conductor_area", 10, "Conductor area out of range [0, 8]"), + ("pole_age", 150, "Pole age out of range [0, 100]"), + ], +) +def test_invalid_range_parameters(attr, value, err_msg): + asset = DistributionPole.example() + setattr(asset, attr, value) + with pytest.raises(ValueError, match=re.escape(err_msg)): + Darestani2019(asset) + + +def test_darestani2019_sampling(): + asset = DistributionPole.example() + custom_dist = Darestani2019(asset) + + # Test sampling and CDF + samples = custom_dist.rvs(size=1000) + assert len(samples) == 1000 + assert all(2 <= s <= 112 for s in samples) + cdf_points = [1, 2, 50, 112, 113] + cdf_values = [custom_dist.cdf(x) for x in cdf_points] + assert cdf_values[0] == 0 + assert cdf_values[4] == 1 + assert all(0 <= cdf <= 1 for cdf in cdf_values) + assert all( + cdf_values[i] <= cdf_values[i + 1] for i in range(len(cdf_values) - 1) + ) # non-decreasing + +def test_darestani2019_calculate_params(): + asset = DistributionPole.example() + custom_dist = Darestani2019(asset) + mu, sigma = custom_dist._calculate_params() + assert isinstance(mu, float) + assert isinstance(sigma, float) + assert sigma > 0 + + # Test with different valid attributes + asset.pole_material = PoleConstructionMaterial.WOOD + asset.pole_class = PoleClass.CLASS_2 + mu2, sigma2 = Darestani2019(asset)._calculate_params() + assert isinstance(mu2, float) + assert isinstance(sigma2, float) + assert sigma2 > 0 + assert (mu, sigma) != (mu2, sigma2) diff --git a/tests/test_fragility_curves.py b/tests/test_fragility_curves.py index 2441c76b..b272483a 100644 --- a/tests/test_fragility_curves.py +++ b/tests/test_fragility_curves.py @@ -31,7 +31,7 @@ def test_invalid_hazard_fragility_model(): FragilityCurve( asset_type=AssetTypes.battery_storage, prob_function=ProbabilityFunction( - distribution="lognorm", parameters=[Speed(35, "cm/s"), 0.5] + distribution="not_valid_dist", parameters=[Speed(35, "cm/s"), 0.5] ), ), ], diff --git a/tests/test_gdm_model_interface.py b/tests/test_gdm_model_interface.py index 42c68450..3a25f599 100644 --- a/tests/test_gdm_model_interface.py +++ b/tests/test_gdm_model_interface.py @@ -1,7 +1,7 @@ from datetime import datetime -from gdm.distribution.components import DistributionBus -from gdm.distribution import DistributionSystem +from gdm.distribution.components.distribution_bus import DistributionBus +from gdm.distribution.distribution_system import DistributionSystem from gdm.quantities import Distance from shapely.geometry import Point diff --git a/tests/test_plot_probability_fuctions.py b/tests/test_plot_probability_fuctions.py index 31ad9a91..e414f83a 100644 --- a/tests/test_plot_probability_fuctions.py +++ b/tests/test_plot_probability_fuctions.py @@ -1,4 +1,8 @@ from erad.default_fragility_curves import DEFAULT_FRAGILTY_CURVES +from erad.plotting import plot_custom_fragility_curves +from erad.models.asset import DistributionPole +from erad.enums import PoleClass, PoleConstructionMaterial +from erad.quantities import WindAngle, ConductorArea, PoleAge def test_plotting(tmp_path): @@ -9,3 +13,33 @@ def test_plotting(tmp_path): hazard_curves.plot(img, 0, 80, 1000) assert img.exists(), "Plotting failed, image not created." + + +def test_plot_custom_fragility_curves(tmp_path): + """Test plotting of per-asset custom fragility curves (Darestani2019).""" + base = DistributionPole.example() + assets = [ + base.model_copy( + update={ + "name": pole_class.name, + "pole_class": pole_class, + "pole_material": PoleConstructionMaterial.WOOD, + "wind_angle": WindAngle(90, "degree"), + "conductor_area": ConductorArea(2, "m**2"), + "pole_age": PoleAge(50, "year"), + "probability_dist": "Darestani2019", + } + ) + for pole_class in PoleClass + ] + + img = tmp_path / "test_plot_custom_fragility.html" + plot_custom_fragility_curves( + assets, + asset_state_param="wind_speed", + file_path=img, + x_min=0, + x_max=60, + number_of_points=200, + ) + assert img.exists(), "Plotting failed, image not created." diff --git a/tests/test_probability_builder.py b/tests/test_probability_builder.py new file mode 100644 index 00000000..50ddbfb0 --- /dev/null +++ b/tests/test_probability_builder.py @@ -0,0 +1,69 @@ +import pytest +from erad.probability_builder import ProbabilityFunctionBuilder +from erad.quantities import Speed +from erad.models.custom_distributions import CUSTOM_DISTRIBUTIONS +from erad.models.asset import DistributionPole +from scipy.stats._distn_infrastructure import rv_continuous_frozen + +# Example scipy distributions to test +SCIPY_DISTS = [ + ("norm", [Speed(1.5, "m/s"), 2], Speed, "meter/second"), + ("lognorm", [Speed(3, "m/s"), 1], Speed, "meter/second"), +] +CUSTOM_DISTS = [ + ("Darestani2019", [Speed(0, "m/s")], Speed, "meter/second"), + # Add other custom distributions here when defined +] + + +@pytest.mark.parametrize("dist_name,params,quantity_cls,units", SCIPY_DISTS) +class TestScipyDistributions: + def test_init_with_scipy_distribution(self, dist_name, params, quantity_cls, units): + b = ProbabilityFunctionBuilder(dist_name, params) + assert isinstance(b.dist, rv_continuous_frozen) + assert b.params == [p.magnitude if hasattr(p, "magnitude") else p for p in params] + assert b.quantity == quantity_cls + assert b.units == units + + def test_probability_for_scipy_distributions(self, dist_name, params, quantity_cls, units): + b = ProbabilityFunctionBuilder(dist_name, params) + p = b.probability(params[0]) + assert isinstance(p, float) and 0 <= p <= 1 + + def test_sample_returns_quantity_instance(self, dist_name, params, quantity_cls, units): + b = ProbabilityFunctionBuilder(dist_name, params) + sample = b.sample() + assert isinstance(sample, b.quantity) and sample.units == b.units + + +@pytest.mark.parametrize("custom_dist_name,params,quantity_cls,units", CUSTOM_DISTS) +class TestCustomDistributions: + def test_init_with_custom_distribution(self, custom_dist_name, params, quantity_cls, units): + asset = DistributionPole.example() + dist_instance = CUSTOM_DISTRIBUTIONS[custom_dist_name](asset) + b = ProbabilityFunctionBuilder(custom_dist_name, params, dist_instance) + assert b.dist is not None + assert b.params == [p.magnitude if hasattr(p, "magnitude") else p for p in params] + assert b.quantity == quantity_cls + assert b.units == units + + def test_probability_for_custom_distributions( + self, custom_dist_name, params, quantity_cls, units + ): + asset = DistributionPole.example() + dist_instance = CUSTOM_DISTRIBUTIONS[custom_dist_name](asset) + b = ProbabilityFunctionBuilder(custom_dist_name, params, dist_instance) + p = b.probability(params[0]) + assert isinstance(p, float) and 0 <= p <= 1 + + def test_custom_distribution_sample(self, custom_dist_name, params, quantity_cls, units): + asset = DistributionPole.example() + dist_instance = CUSTOM_DISTRIBUTIONS[custom_dist_name](asset) + b = ProbabilityFunctionBuilder(custom_dist_name, params, dist_instance) + sample = b.sample() + assert isinstance(sample, b.quantity) and sample.units == b.units + + +def test_unsupported_distribution_raises(): + with pytest.raises(ValueError, match="Unsupported distribution"): + ProbabilityFunctionBuilder("not_a_real_dist", [Speed(1.0, "m/s")])