import math
import os
import tempfile
import xml.etree.ElementTree as ET
from collections.abc import Iterator
from contextlib import contextmanager
from importlib.metadata import PackageNotFoundError, version
from importlib.util import find_spec
from pathlib import Path
from openworld_radio_twin.environment import (
RF_MATERIAL_PROFILE_VERSION,
RF_MATERIAL_PROFILES,
EnvironmentData,
SurfaceClass,
)
from openworld_radio_twin.models import (
SIONNA_PROPAGATION_PARAMETERS,
BuildingFeature,
SimulationRequest,
SionnaEngineConfig,
)
from openworld_radio_twin.simulation.radio_map import RadioMapData
from openworld_radio_twin.simulation.scene_builder import (
MeasurementSurfaceAssets,
SceneAssets,
build_scene_assets,
measurement_surface_path,
transmitter_scene_position,
write_measurement_surface,
)
SIONNA_RADIO_MAP_CONTRACT_VERSION = 2
[docs]
def availability() -> tuple[bool, str]:
if find_spec("sionna") is not None or find_spec("sionna_rt") is not None:
return True, "Sionna RT runtime detected"
return False, "Install the optional 'rt' dependency to enable Sionna RT"
[docs]
def runtime_details() -> dict[str, str | int | float | bool]:
ready, message = availability()
details: dict[str, str | int | float | bool] = {"available": ready, "message": message}
if ready:
try:
details["version"] = version("sionna-rt")
except PackageNotFoundError:
details["version"] = "unknown"
return details
def _select_mitsuba_variant(mi, device: str) -> str:
"""Select a Sionna-compatible backend before importing ``sionna.rt``."""
variants = {
"cuda": "cuda_ad_mono_polarized",
"cpu": "llvm_ad_mono_polarized",
}
if device not in {"auto", *variants}:
raise ValueError("device must be 'auto', 'cuda', or 'cpu'")
current = mi.variant()
compatible = set(variants.values())
if current in compatible and device == "auto":
return current
requested = variants.get(device)
if current is not None:
if current == requested:
return current
requirement = requested or "a Sionna-compatible CUDA or LLVM variant"
raise RuntimeError(
f"Mitsuba variant is already '{current}' on this thread; device='{device}' "
f"requires {requirement}. Select the device before importing sionna.rt, or "
"use a separate process."
)
try:
if device == "auto":
mi.set_variant("cuda_ad_mono_polarized", "llvm_ad_mono_polarized")
else:
assert requested is not None
mi.set_variant(requested)
except (ImportError, RuntimeError) as exc:
if device == "cuda":
raise RuntimeError(
"Sionna was configured for CUDA, but Mitsuba/Dr.Jit cannot initialize a "
"visible NVIDIA GPU. Check the driver, CUDA allocation, container GPU access, "
"and CUDA_VISIBLE_DEVICES, or use device='auto'/'cpu'."
) from exc
raise RuntimeError(f"Sionna could not initialize its {device} backend: {exc}") from exc
selected = mi.variant()
if selected not in compatible:
raise RuntimeError("Mitsuba did not select a Sionna-compatible execution variant")
return selected
def _shape_is_enabled(shape_id: str, config: SionnaEngineConfig) -> bool:
"""Return whether an XML scene shape participates in this solver run."""
name = shape_id.removeprefix("mesh-")
enabled = {
"terrain-ground": config.include_ground,
"terrain-vegetation": config.include_vegetation,
"terrain-paved": config.include_paved,
"terrain-water": config.include_water,
}
if name in enabled:
return enabled[name]
if name.startswith("building_") and name.endswith("_rooftop"):
return config.include_building_roofs
if name.startswith("building_") and name.endswith("_wall"):
return config.include_building_walls
return True
def _geometry_class(shape_name: str) -> str:
if shape_name.startswith("terrain-"):
return shape_name.removeprefix("terrain-")
if shape_name.startswith("building_") and shape_name.endswith("_rooftop"):
return "building_roofs"
if shape_name.startswith("building_") and shape_name.endswith("_wall"):
return "building_walls"
return "other"
@contextmanager
def _masked_scene_path(
source: Path,
config: SionnaEngineConfig,
) -> Iterator[tuple[Path, tuple[str, ...]]]:
"""Yield an XML scene with disabled shape classes removed before Sionna loads it."""
tree = ET.parse(source)
root = tree.getroot()
removed: list[str] = []
for shape in list(root.findall("shape")):
shape_id = shape.get("id", "")
if not _shape_is_enabled(shape_id, config):
root.remove(shape)
removed.append(shape_id.removeprefix("mesh-"))
if not removed:
yield source, ()
return
handle = tempfile.NamedTemporaryFile(
prefix=".owrt-sionna-mask-",
suffix=".xml",
dir=source.parent,
delete=False,
)
target = Path(handle.name)
handle.close()
try:
tree.write(target, encoding="utf-8", xml_declaration=True)
yield target, tuple(removed)
finally:
target.unlink(missing_ok=True)
[docs]
def simulate_sionna(
request: SimulationRequest,
buildings: list[BuildingFeature],
environment: EnvironmentData | None = None,
) -> tuple[RadioMapData, dict[str, str | int | float | bool], list[str]]:
"""Build an ephemeral local scene and compute a georeferenced Sionna radio map."""
if not availability()[0]:
raise RuntimeError("Sionna RT is not installed")
cache_dir = Path(tempfile.gettempdir()) / "openworld-radio-twin" / "matplotlib"
cache_dir.mkdir(parents=True, exist_ok=True)
os.environ.setdefault("MPLCONFIGDIR", str(cache_dir))
origin = request.scene_origin
with tempfile.TemporaryDirectory(prefix="owrt-scene-") as temp_directory:
assets = build_scene_assets(
Path(temp_directory),
buildings,
origin.longitude,
origin.latitude,
request.radius_m,
environment,
True,
request.material_profile,
)
return simulate_sionna_scene(request, assets, len(buildings))
[docs]
def simulate_sionna_scene(
request: SimulationRequest,
assets: SceneAssets,
source_building_count: int | None = None,
) -> tuple[RadioMapData, dict[str, str | int | float | bool], list[str]]:
"""Compute a radio map from persistent scene assets shared by batch cases."""
if not availability()[0]:
raise RuntimeError("Sionna RT is not installed")
import mitsuba as mi
config = SionnaEngineConfig.model_validate(request.engine_config)
variant = _select_mitsuba_variant(mi, config.device)
import numpy as np
from sionna.rt import PlanarArray, RadioMapSolver, Transmitter, load_mesh, load_scene
reference = request.reference_transmitter
geometry_masked = not all(
(
config.include_ground,
config.include_vegetation,
config.include_paved,
config.include_water,
config.include_building_roofs,
config.include_building_walls,
)
)
with _masked_scene_path(assets.scene_xml_path, config) as (
solver_scene_path,
removed_geometry,
):
scene = load_scene(str(solver_scene_path))
scene.frequency = reference.frequency_ghz * 1e9
pattern = "iso" if reference.antenna_pattern == "isotropic" else "tr38901"
scene.tx_array = PlanarArray(
num_rows=1 if pattern == "iso" else 8,
num_cols=1 if pattern == "iso" else 2,
pattern=pattern,
polarization="V",
)
scene.rx_array = PlanarArray(
num_rows=1,
num_cols=1,
pattern="iso",
polarization="V",
)
for index, tx_config in enumerate(request.transmitters):
tx_position = transmitter_scene_position(
assets,
tx_config.position.longitude,
tx_config.position.latitude,
tx_config.position.altitude_m,
)
azimuth = math.radians(tx_config.azimuth_deg)
downtilt = math.radians(tx_config.downtilt_deg)
horizontal = 100 * math.cos(downtilt)
look_at = [
tx_position[0] + horizontal * math.sin(azimuth),
tx_position[1] + horizontal * math.cos(azimuth),
tx_position[2] - 100 * math.sin(downtilt),
]
scene.add(
Transmitter(
name=f"tx-{index}",
position=list(tx_position),
look_at=look_at,
power_dbm=tx_config.power_dbm,
)
)
measurement = write_measurement_surface(
measurement_surface_path(assets, request.resolution_m, request.receiver_height_m),
assets,
request.resolution_m,
request.receiver_height_m,
)
measurement_mesh = load_mesh(str(measurement.path), flip_normals=False)
propagation = config.model_dump(include=set(SIONNA_PROPAGATION_PARAMETERS))
radio_map = RadioMapSolver()(
scene,
measurement_surface=measurement_mesh,
samples_per_tx=request.samples_per_tx,
max_depth=request.max_depth,
seed=request.seed,
**propagation,
)
thermal_noise_power_w = float(scene.thermal_noise_power[0])
data = _aggregate_mesh_radio_map(
radio_map,
measurement,
request.association_metric,
thermal_noise_power_w=thermal_noise_power_w,
)
if not np.any(data.rss > 0):
raise RuntimeError("Sionna returned no positive RSS samples")
served_cells = np.any(data.rss > 0, axis=0)
receiver_cells = data.rows * data.columns
transmitter_count = len(request.transmitters)
excluded_classes = sorted({_geometry_class(name) for name in removed_geometry})
details: dict[str, str | int | float | bool] = {
"version": version("sionna-rt"),
"variant": variant,
"device": "GPU" if variant.startswith("cuda") else "CPU",
"device_policy": config.device,
"geometry_masked": geometry_masked,
"excluded_geometry_count": len(removed_geometry),
"excluded_geometry_classes": ",".join(excluded_classes) or "none",
"ground_geometry": (
"excluded"
if not config.include_ground
else "terrain-following"
if assets.environment is not None
else "flat"
),
"buildings_in_scene": (
assets.building_count
if config.include_building_roofs or config.include_building_walls
else 0
),
"source_buildings": assets.building_count,
"triangles": assets.triangle_count,
"triangle_count_scope": "source assets before solver geometry mask",
"samples_per_tx": request.samples_per_tx,
"frequency_hz": float(reference.frequency_ghz * 1e9),
"noise_bandwidth_hz": float(scene.bandwidth[0]),
"temperature_k": float(scene.temperature[0]),
"thermal_noise_power_w": thermal_noise_power_w,
"transmitter_count": transmitter_count,
"tx_array": "1x1_iso" if pattern == "iso" else "8x2_tr38901",
"array_contract": "linear_float32_(n_tx,rows,columns)",
"cell_size_m": request.resolution_m,
"receiver_cells": receiver_cells,
"measurement_area_m2": (2 * request.radius_m) ** 2,
"initial_ray_samples_total": request.samples_per_tx * transmitter_count,
"sample_budget_scope": "per_transmitter",
"samples_to_cells_ratio_per_tx": request.samples_per_tx / receiver_cells,
"samples_to_cells_ratio_note": (
"samples_per_tx / receiver_cells; diagnostic only, not measured cell hits"
),
"served_cell_fraction": float(np.mean(served_cells)),
"association_metric": request.association_metric,
"max_depth": request.max_depth,
"solver_seed": request.seed,
"propagation": "Sionna RadioMapSolver keyword arguments",
**propagation,
"rf_material_profile_version": RF_MATERIAL_PROFILE_VERSION,
"materials": ",".join(
sorted(
{profile.sionna_material for profile in RF_MATERIAL_PROFILES.values()}
| {"itu_metal", "itu_concrete"}
)
),
"materials_scope": "source scene declarations; excluded_geometry_classes apply",
"coordinate_frame": "WGS84 topocentric ENU in meters",
"terrain_enabled": assets.environment is not None,
"measurement_surface": (
"terrain-following triangular mesh"
if assets.environment is not None
else "planar triangular mesh"
),
"measurement_surface_triangles": int(measurement.triangle_areas_m2.size),
"triangle_to_cell_reduction": "surface-area-weighted linear mean",
"triangle_to_cell_reduction_metrics": "path_gain,rss",
"sinr_definition": "cell_rss / (sum_other_transmitter_cell_rss + thermal_noise_power_w)",
"radio_map_contract_version": SIONNA_RADIO_MAP_CONTRACT_VERSION,
}
if assets.environment is not None:
water_model = next(
(
operation
for operation in assets.environment.terrain.model_operations
if operation.get("operation") == "water_surface_hydroflattening"
),
None,
)
details.update(
{
"terrain_provider": assets.environment.terrain.source.provider,
"terrain_mesh_spacing_m": assets.environment.terrain.spacing_m,
"terrain_vertical_datum": assets.environment.terrain.source.vertical_datum,
"measurement_surface_path": str(measurement.path.relative_to(assets.directory)),
"surface_material_profiles": "semantic classes with documented RF surrogates",
"water_material": RF_MATERIAL_PROFILES[SurfaceClass.WATER].sionna_material,
"marine_water_surface": (
"horizontal vertical-datum-zero interface"
if water_model is not None and water_model["marine_cell_count"]
else "not present"
),
"marine_water_cells": (
int(water_model["marine_cell_count"]) if water_model is not None else 0
),
"water_surface_components": (
int(water_model["component_count"]) if water_model is not None else 0
),
"dem_estimated_water_components": (
int(water_model["dem_estimated_component_count"])
if water_model is not None
else 0
),
"water_anchored_buildings": sum(
building.anchor_surface_class == SurfaceClass.WATER
for building in assets.local_buildings
),
}
)
warnings = []
if source_building_count is not None and assets.building_count < source_building_count:
warnings.append(
f"Skipped {source_building_count - assets.building_count} invalid building footprints."
)
return data, details, warnings
def _aggregate_mesh_radio_map(
radio_map,
measurement: MeasurementSurfaceAssets,
association_metric: str,
*,
thermal_noise_power_w: float,
) -> RadioMapData:
"""Area-average triangle powers, then compute SINR on the exported cells."""
import numpy as np
if not math.isfinite(thermal_noise_power_w) or thermal_noise_power_w < 0:
raise ValueError("Thermal noise power must be finite and non-negative")
rows, columns = measurement.cell_shape
weights = measurement.triangle_areas_m2
weight_sum = np.sum(weights, axis=-1)
def aggregate(values) -> np.ndarray:
array = np.asarray(values.numpy(), dtype=np.float32)
if array.ndim == 1:
array = array[np.newaxis, :]
expected = rows * columns * 2
if array.shape[-1] != expected:
raise RuntimeError(
f"Sionna mesh radio map returned {array.shape[-1]} triangles; expected {expected}"
)
triangles = array.reshape(array.shape[0], rows, columns, 2)
return np.asarray(
np.sum(triangles * weights[np.newaxis, ...], axis=-1) / weight_sum[np.newaxis, ...],
dtype=np.float32,
)
path_gain = aggregate(radio_map.path_gain)
rss = aggregate(radio_map.rss)
sinr = np.empty_like(rss)
for tx_index in range(rss.shape[0]):
# Sum other Tx powers directly: total - desired can erase weak
# interference when the desired signal is many orders stronger.
interference = np.sum(rss[:tx_index], axis=0, dtype=np.float64)
interference += np.sum(rss[tx_index + 1 :], axis=0, dtype=np.float64)
with np.errstate(divide="ignore", invalid="ignore"):
sinr[tx_index] = rss[tx_index] / (interference + thermal_noise_power_w)
selected = {"path_gain": path_gain, "rss": rss, "sinr": sinr}[association_metric]
association = np.argmax(selected, axis=0).astype(np.int32)
association[~np.any(selected > 0, axis=0)] = -1
return RadioMapData(
path_gain=path_gain,
rss=rss,
sinr=sinr,
association=association,
cell_centers=measurement.cell_centers,
)
