Compare radio map solvers

Run the deterministic, Monte Carlo Basic and BDPT solvers on one radio map, align their result layouts, and see where and why their values differ.

This tutorial computes one small radio map with all three radio-map solvers and compares them cell by cell. It shows how to bring their results into one layout and explains the differences you will see. Use it to decide which solver to trust for your own maps.

The scene #

A concrete-like wall (eps_r = 4.0, sigma_e = 0.02 S/m) stands in the plane x = 2.5. A transmitter at (0, -1, 0.5) covers a 4 by 8 grid with 0.5 m spacing, 0.5 m above the transmitter. Paths are line of sight and single reflections.

Note

The grid is deliberately not at the transmitter height. Monte Carlo Basic counts reflected rays where they cross the grid plane, and with only vertical walls no reflected ray crosses a plane at the transmitter height. See Radio maps.

Run all solvers #

import time
 
import torch
 
from witwin.core import AntennaState, Mesh, PhysicalMaterial, ReceiverGrid, Scene, Structure
from witwin.channel import deterministic
from witwin.channel.montecarlo import basic, bdpt
 
FREQUENCY_HZ = 3.0e9
COMPONENTS = {"los", "reflection"}
 
vertices = torch.tensor([[2.5, -3.0, -1.0], [2.5, 3.0, -1.0], [2.5, -3.0, 2.0], [2.5, 3.0, 2.0]])
faces = torch.tensor([[0, 1, 2], [1, 3, 2]])
wall = Structure(Mesh(vertices, faces, recenter=False, fill_mode="surface"), PhysicalMaterial(eps_r=4.0, sigma_e=0.02), name="wall")
tx = AntennaState(1, "tx", torch.tensor([0.0, -1.0, 0.5]), power_w=1.0)
# 4 rows along +x and 8 columns along -y, 0.5 m apart, 0.5 m above the transmitter.
grid = ReceiverGrid(
    3,
    origin=torch.tensor([-2.0, 2.0, 1.0]),
    x_axis=torch.tensor([1.0, 0.0, 0.0]),
    y_axis=torch.tensor([0.0, -1.0, 0.0]),
    shape=(4, 8),
    spacing=(0.5, 0.5),
)
scene = Scene(structures=[wall], endpoints=[tx, grid])
 
 
def timed(label, run):
    run()  # warm-up: the first call compiles the scene and loads native kernels
    torch.cuda.synchronize()
    start = time.perf_counter()
    result = run()
    torch.cuda.synchronize()
    print(f"{label:24s} shape {tuple(result.path_gain.shape)}  {1e3 * (time.perf_counter() - start):7.1f} ms")
    return result
 
 
coherent = timed("deterministic coherent", lambda: deterministic.solve(
    scene, deterministic.Config(max_depth=1, components=COMPONENTS), reference_frequency_hz=FREQUENCY_HZ))
incoherent = timed("deterministic incoherent", lambda: deterministic.solve(
    scene, deterministic.Config(max_depth=1, components=COMPONENTS, coherent=False), reference_frequency_hz=FREQUENCY_HZ))
mc = timed("montecarlo.basic", lambda: basic.solve(
    scene, basic.Config(samples=1 << 20, max_depth=1, components=COMPONENTS, seed=1), reference_frequency_hz=FREQUENCY_HZ))
bd = timed("montecarlo.bdpt", lambda: bdpt.solve(
    scene, bdpt.Config(samples=1 << 20, max_depth=1, components=COMPONENTS, seed=1), reference_frequency_hz=FREQUENCY_HZ))
 
# Bring every map to the deterministic layout (T, columns, rows) and compare with the incoherent sum.
reference = incoherent.path_gain
maps = {
    "deterministic coherent": coherent.path_gain,
    "montecarlo.basic": mc.path_gain.reshape(reference.shape),
    "montecarlo.bdpt": bd.path_gain,
}
for label, values in maps.items():
    ratio = (values.sum() / reference.sum()).item()
    error_db = (10 * torch.log10(values / reference)).abs()
    print(f"{label:24s} total / incoherent {ratio:.4f}  median |dB| {error_db.median().item():.2f}  max |dB| {error_db.max().item():.2f}")

Output (RTX 5080, times after warm-up):

deterministic coherent   shape (1, 8, 4)     11.6 ms
deterministic incoherent shape (1, 8, 4)     10.8 ms
montecarlo.basic         shape (1, 32)        3.7 ms
montecarlo.bdpt          shape (1, 8, 4)     12.1 ms
deterministic coherent   total / incoherent 0.9845  median |dB| 0.21  max |dB| 1.01
montecarlo.basic         total / incoherent 1.2215  median |dB| 0.25  max |dB| 3.01
montecarlo.bdpt          total / incoherent 1.0000  median |dB| 0.00  max |dB| 0.00

Align the layouts #

The grid has shape=(4, 8), that is 4 rows and 8 columns. The deterministic and BDPT solvers return (T, columns, rows), here (1, 8, 4). Monte Carlo Basic returns a flat (T, rows * columns) map, here (1, 32); reshape(T, columns, rows) puts it into the same layout. Its component_maps entries already have that layout. See Radio maps.

Read the differences #

The reference is the incoherent deterministic map: every traced path, powers added.

Coherent vs incoherent deterministic. The coherent map adds complex fields, so line of sight and reflection interfere. Cells differ from the power sum by up to 1 dB in either direction, and the totals agree within 2 %. The coherent map is what a narrowband receiver at that exact point sees; the incoherent one is its average over small movements.

BDPT matches the incoherent deterministic map in this scene to four digits, polarization included. Do not expect this in general: in scenes with deeper paths, diffraction or transmission it is an estimate like any Monte Carlo result.

Monte Carlo Basic is 22 % higher in total, with the largest difference, 3.01 dB, at the cell closest to the transmitter. The cause is polarization, not noise. capabilities() lists no polarization support for Monte Carlo Basic (see Solvers and path types). For line of sight in this scene its value is the Friis gain times sin(theta)**2, where theta is the angle from the vertical. The deterministic value is the Friis gain times sin(theta)**4: one factor for the vertical transmitter's pattern and one for the vertical receiver's projection. The two agree for receivers in the transmitter's horizontal plane and differ most at steep angles. The nearest cell, at (-0.5, -1.0, 1.0), sits 45 degrees above the transmitter, where sin(theta)**2 = 0.5: a factor of 2, which is exactly 3.01 dB.

How many samples #

Monte Carlo Basic computes line of sight per cell, so samples only affects the sampled components. Add this to the script to watch the reflection map converge:

reference_reflection = incoherent.component_power["reflection"]
for samples in (1 << 12, 1 << 16, 1 << 20, 1 << 22):
    estimate = basic.solve(
        scene, basic.Config(samples=samples, max_depth=1, components=COMPONENTS, seed=1), reference_frequency_hz=FREQUENCY_HZ
    ).component_maps["reflection"]
    empty = int((estimate == 0).sum())
    ratio = (estimate.sum() / reference_reflection.sum()).item()
    print(f"samples {samples:8d}  reflection total / deterministic {ratio:.4f}  empty cells {empty}")

samples     4096  reflection total / deterministic 0.9545  empty cells 27
samples    65536  reflection total / deterministic 1.0465  empty cells 0
samples  1048576  reflection total / deterministic 1.0204  empty cells 0
samples  4194304  reflection total / deterministic 1.0174  empty cells 0

With 4,096 rays (the default), 27 of the 32 cells receive no reflected ray at all. From 65,536 rays on, every cell gets reflected rays; from about one million rays the total stays within 2 % of the deterministic value. As a rule, raise samples until no cell is empty and the map stops changing; finer grids need more rays, because each cell is smaller.

Which solver to use #

Need	Solver
Exact, repeatable values with interference, polarization and gradients	deterministic (`coherent=True`)
A smooth coverage map without fringes	deterministic (`coherent=False`) or BDPT
A fast power estimate over a large grid	Monte Carlo Basic, with enough `samples`; mind polarization
Paths that are hard to reach from the transmitter side	BDPT

The scene#

Run all solvers#

Align the layouts#

Read the differences#

How many samples#

Which solver to use#