Reproducibility of the PathSolver results

[llacinak]

I've gone through multiple recent existing discussions regarding this topic (sorry if I missed something):

• #851
• #937
• #1071

which decribe the situations when the PathSolver outputs might be non-deterministic even if the seed is specified.

In fact, in Sionna RT 2.0.0, a valid seed is always specified since the default value is 42.

I've tried to prepare a simple scenario which attemps to avoid the configurations that were provided as explanantion for the non-deterministic behavior and test it with only the specular_reflection enabled:

import mitsuba as mi
import numpy as np
import sionna.rt as rt
from sionna.rt import PathSolver, PlanarArray, Receiver, Scene, Transmitter, load_scene


def compute_paths(solver: PathSolver, scene: Scene, seed: int) -> dict:
    # Does this have any impact on the solver (or is the behavior supposed to be fully controlled by the solver seed)?
    np.random.seed(seed)
    # Using the default values of samples_per_src and max_num_paths_per_src which are identical (1000000) in an effort
    # to avoid issues with dropping the paths due to running out of entries in the buffer when
    # max_num_paths_per_src << samples_per_src reported in
    # https://github.com/NVlabs/sionna/discussions/851#discussioncomment-12876333
    paths = solver(scene=scene,
                   los=False,
                   specular_reflection=True,
                   diffuse_reflection=False,
                   refraction=False,
                   diffraction=False,
                   edge_diffraction=False,
                   diffraction_lit_region=False,
                   seed=seed)

    # Sort the rays based on the propagation delays to suppress the issue with non-deterministic ordering of paths
    # in the solver output mentioned in
    # https://github.com/NVlabs/sionna/discussions/851#discussioncomment-12876333
    cir_a, cir_tau = paths.cir(normalize_delays=False, out_type='numpy')
    delay_sorted_path_indices = np.argsort(cir_tau, axis=None)
    return {
        "cir_a": cir_a[:, :, :, :, delay_sorted_path_indices, :],
        "cir_tau": cir_tau[:, :, delay_sorted_path_indices],
        "doppler": np.array(paths.doppler)[:, :, delay_sorted_path_indices],
        "phi_r": np.array(paths.phi_r)[:, :, delay_sorted_path_indices],
        "theta_r": np.array(paths.theta_r)[:, :, delay_sorted_path_indices]
    }


def test_seeds():
    munich_scene = load_scene(rt.scene.munich)

    # Configure minimal antenna arrays for TX/RX
    munich_scene.tx_array = PlanarArray(num_rows=1, num_cols=1, pattern="iso", polarization="V")
    munich_scene.rx_array = PlanarArray(num_rows=1, num_cols=1, pattern="iso", polarization="V")

    # Create simple single transmitter and single receiver scenario that should not have the issue
    # with non-deterministic behavior due to hash collisions when simulating many receivers as described in
    # https://github.com/NVlabs/sionna/discussions/937#discussioncomment-13520118
    tx = Transmitter(name="tx", position=mi.Point3f(85.0, 15.0, 18.0))
    munich_scene.add(tx)
    rx = Receiver(name="rx", position=mi.Point3f(50.0, 20.0, 0.25))
    munich_scene.add(rx)

    solver = PathSolver()

    # Reference solver run
    ref_seed = 1981
    ref_paths = compute_paths(solver=solver, scene=munich_scene, seed=ref_seed)
    print(f"Sorted path delays with reference seed {ref_seed}:")
    print(f"ref_cir_tau: {ref_paths['cir_tau']}")

    # Test solver run for the same seed as the reference
    test1_seed = ref_seed
    test1_paths = compute_paths(solver=solver, scene=munich_scene, seed=test1_seed)
    print(f"\nComparison of the results with reference seed {ref_seed} and test seed {test1_seed}:")
    print(f"test1_cir_a - ref_cir_a: {test1_paths['cir_a'] - ref_paths['cir_a']}")
    print(f"test1_cir_tau - ref_cir_tau: {test1_paths['cir_tau'] - ref_paths['cir_tau']}")
    print(f"test1_doppler - ref_doppler: {test1_paths['doppler'] - ref_paths['doppler']}")
    print(f"test1_phi_r - ref_phi_r: {test1_paths['phi_r'] - ref_paths['phi_r']}")
    print(f"test1_theta_r - ref_theta_r: {test1_paths['theta_r'] - ref_paths['theta_r']}")

    # Test solver run for different seed comparing to the reference
    test2_seed = ref_seed + 1
    test2_paths = compute_paths(solver=solver, scene=munich_scene, seed=test2_seed)
    print(f"\nComparison of the results with reference seed {ref_seed} and test seed {test2_seed}:")
    print(f"test2_cir_a - ref_cir_a: {test2_paths['cir_a'] - ref_paths['cir_a']}")
    print(f"test2_cir_tau - ref_cir_tau: {test2_paths['cir_tau'] - ref_paths['cir_tau']}")
    print(f"test2_doppler - ref_doppler: {test2_paths['doppler'] - ref_paths['doppler']}")
    print(f"test2_phi_r - ref_phi_r: {test2_paths['phi_r'] - ref_paths['phi_r']}")
    print(f"test2_theta_r - ref_theta_r: {test2_paths['theta_r'] - ref_paths['theta_r']}")


if __name__ == "__main__":
    test_seeds()

The paths returned by the solver when simulating on a CPU still show sporadic numeric differences and the behavior seems to be independent of the used seed. In fact the following printout shows matching results for different seeds but differences for the same seed:

Sorted path delays with reference seed 1981:
ref_cir_tau: [[[6.8951073e-07 6.8965591e-07 7.1405566e-07 7.1419583e-07 7.1623856e-07
   7.1637822e-07 7.6630971e-07 7.6644034e-07 7.9020339e-07 7.9188050e-07
   9.6530630e-07 1.1211164e-06 1.1455212e-06 1.1474718e-06 1.4558577e-06]]]

Comparison of the results with reference seed 1981 and test seed 1981:
test1_cir_a - ref_cir_a: [[[[[[0.+0.0000000e+00j]
     [0.+1.8189894e-12j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]
     [0.+0.0000000e+00j]]]]]]
test1_cir_tau - ref_cir_tau: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]
test1_doppler - ref_doppler: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]
test1_phi_r - ref_phi_r: [[[ 0.0000000e+00 -1.1920929e-07  0.0000000e+00  0.0000000e+00
    0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00
    0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00
    0.0000000e+00  0.0000000e+00  0.0000000e+00]]]
test1_theta_r - ref_theta_r: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]

Comparison of the results with reference seed 1981 and test seed 1982:
test2_cir_a - ref_cir_a: [[[[[[0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]
     [0.+0.j]]]]]]
test2_cir_tau - ref_cir_tau: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]
test2_doppler - ref_doppler: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]
test2_phi_r - ref_phi_r: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]
test2_theta_r - ref_theta_r: [[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]

Is there still some source of randomness that I haven't controlled properly?

Is there actually a scenario with which one would get bit-exactly reproducible results between runs (scenario that would use the random number generators in the solver)? I am interested in the exact reproducibility for the purpose of regression testing...

Any feedback is highly appreciated.

[merlinND]

Hello @llacinak,

Thank you for checking the previous discussions. As this is a feature that is often requested, we are working on a path candidates generation algorithm that always produces the same paths regardless of output buffer constraints (max_num_paths_per_src) or nondeterministic thread ordering. One downside will be that it uses much more memory while running.

Note however that the goal is only to produce the same set of paths, not to generate exactly the same numerical results each time (bitwise). The reason is that any quantity computed with a parallel reduction (e.g. dr.sum()) introduces small nondeterministic numerical errors.
For regression testing, it should be possible to add numerical tolerance with e.g. dr.allclose() / np.allclose(). In your test above, are the differences really significant, or could they be considered "close enough"?

[llacinak]

Thank you for a quick reply. Few further clarification questions:

1. Does setting max_num_paths_per_src and samples_per_src to the same value as done in the example above work as workaround for the output buffer constraints issues mentioned above? Or is it not enough?
2. Is the solver seed value itself enough to achieve as much reproducibility as currently possible? Or does one need to control other seeds (e.g. np.random.seed()) as well?

When it comes to the "close enough" question, I see differences in the affected channel coefficients with the relative error greater then 0.1% (e.g 5.0264234e-06+2.62638764e-06j vs. 5.0212852e-06+2.63619950e-06j). That seems a bit excessive to be just a numerical error.

[merlinND]

Hello @llacinak,

1. That should be almost enough, except that the LoS paths are not included in samples_per_src, so you should increase max_num_paths_per_src by the number of receivers.
2. Yes, the seed parameter controls all randomness, there should be no need to set np.random.seed().

If you observe these differences in channel coefficients despite the paths being exactly the same, I think the only source of variation is really the nondeterministic accumulation I described above.

[merlinND]

Closing as the new deterministic PathSolver option in the upcoming release will hopefully fix this.