Sionna RT 1.0.2 Seeds / Reproducibility

[rajb245]

The script below loads a standard scene, puts down a 2D grid of receivers, and then instantiates a new solver and traces from transmitter to the receivers twice (forward() function is called twice and that instantiates the solver) with the seed set, and then tries to subtract the two lowpass filtered sets of CIR data from each other to get zero. This doesn't work and the loss is non-zero. What am I missing to get deterministic / reproducible ray traces through my scene is Sionna RT 1.0.2?

import drjit as dr
import sionna.rt
import matplotlib.pyplot as plt
import numpy as np
import math
# Import relevant components from Sionna RT
from sionna.rt import load_scene, PlanarArray, Transmitter, Receiver, Camera, PathSolver

print(dr.__version__)
print(sionna.rt.__version__)

scene = load_scene(sionna.rt.scene.munich, merge_shapes=True) # Merge shapes to speed-up computations
my_cam = Camera(position=[-250,250,150], look_at=[-15,30,28])

# Configure antenna array for all transmitters
scene.tx_array = PlanarArray(num_rows=1,
                             num_cols=1,
                             vertical_spacing=0.5,
                             horizontal_spacing=0.5,
                             pattern="tr38901",
                             polarization="V")

# Configure antenna array for all receivers
scene.rx_array = PlanarArray(num_rows=1,
                             num_cols=1,
                             vertical_spacing=0.5,
                             horizontal_spacing=0.5,
                             pattern="dipole",
                             polarization="V") #cross is cross pol two dipoles

# Create transmitter
tx = Transmitter(name="tx",
                 position=[8.5,21,27],
                 display_radius=2)

# Add transmitter instance to scene
scene.add(tx)

rx_center = [45,75,1.5]
# a rectangular equipsaced grid
dx = 3 #what are these units
dy = 3
x_extent = 50
y_extent = 110

nx = math.floor(x_extent/dx)
ny = math.floor(y_extent/dy)

# Create the initial meshgrid
rx_x = rx_center[0] + x_extent*np.arange(nx)/nx - x_extent/2
rx_y = rx_center[1] + y_extent*np.arange(ny)/ny - y_extent/2
rx_x, rx_y = np.meshgrid(rx_x, rx_y)

# Define the angle of rotation in radians
theta = math.radians(-24)

# Define the center of rotation
center_x = rx_center[0]
center_y = rx_center[1]

# Translate the grid so the center of rotation is at the origin
# We use temporary variables to store the translated coordinates
temp_x = rx_x - center_x
temp_y = rx_y - center_y

# Perform the rotation and store directly into rx_x and rx_y
# It's important to calculate both new coordinates before assigning them
# to avoid using an already partially updated rx_x or rx_y in the second calculation.
new_rx_x = temp_x * math.cos(theta) - temp_y * math.sin(theta)
new_rx_y = temp_x * math.sin(theta) + temp_y * math.cos(theta)

# Translate the grid back to its original position and assign to original variables
rx_x = new_rx_x + center_x
rx_y = new_rx_y + center_y

rx_z = rx_center[2]*np.ones_like(rx_x) #put them all at the same height

rx_positions = np.stack((rx_x,rx_y,rx_z)).reshape(3,-1)
print(rx_positions.shape)
rx = [None,]*rx_positions.shape[-1]
for i, position in enumerate(rx_positions.T):
    rx[i] = Receiver(name="rx%05d"%i,
              position=position,
              display_radius=5)
    scene.add(rx[i])

tx.look_at(rx_center) # Transmitter points towards receiver

# OFDM system parameters
# these are for mu=0, 40 MHz 5G
Nrb = 106
num_subcarriers = 2**math.ceil(math.log2(Nrb*12))
subcarrier_spacing=15e3

#forward trace through the scene that returns lowpass filtered sampled CIRs
def forward():
    p_solver = PathSolver()
    p_solver.loop_mode = 'evaluated'
    paths = p_solver(scene=scene,
                 max_depth=5,
                 los=True,
                 specular_reflection=True,
                 diffuse_reflection=False,
                 refraction=True,
                 synthetic_array=False,
                 seed=41)
    return paths.taps(bandwidth=num_subcarriers*subcarrier_spacing, # Bandwidth to which the channel is low-pass filtered
                  l_min=-6,        # Smallest time lag
                  l_max=100,       # Largest time lag
                  sampling_frequency=None, # Sampling at Nyquist rate, i.e., 1/bandwidth
                  normalize=True,  # Normalize energy
                  normalize_delays=True,
                  )

# this is the ground truth from the unperturbed scene
Y_r, Y_i = forward() # drjit tensors are real valued so this is re and im
#we will perturb the scene here, but the initial check is to see zero loss
# get the output of the model and calculate the loss
Yhat_r, Yhat_i = forward()          # each is a Dr.Jit tensor
loss = dr.sum(dr.squared_norm(Yhat_r - Y_r) + dr.squared_norm(Yhat_i - Y_i))
print(loss)

[merlinND]

Hello @rajb245,

I am closing as duplicate of #851, please see this answer in particular.
Feel free to post again here if you think it is due to a different issue though.

[rajb245]

No that makes sense; I should have checked that the order wasn't different. So if I manage to sort things after the raytrace into some canonical order, say by receiver, this snippet should give zero error between the two calls to the solver right? Will try to implement that soon.

[jhoydis]

Hi @rajb245,

I think that the order does not matter here because you compute the channel taps which is a low-pass filtered and discretely sampled version of the channel impulse response. We'll look into it.

[jhoydis]

I have looked into this. The problem in your code is that you simulate a very large number of receivers at the same time together with a large maximum depth. This can lead to non-deterministic hash collisions which can impact the propagation paths that are found.

You can increase the number of samples per source and/or reduce the number of receivers to make the results increasingly deterministic. There are also differences due to certain mathematical operations on large tensors which are non-deterministic on GPUs with floating-point arithmetic.

We are aware that this is a current limitation and might consider improving it in a future release.

[rajb245]

Increasing samples_per_src seems to do the trick; at least I can get numerically small values for the difference between runs. Thansk!