@mytec: WebGL works

backend/app/services/gpu_service.py

@@ -436,6 +436,139 @@ class GPUService:
 
         return _to_cpu(rsrp)
 
+    def calculate_interference(
+        self,
+        rsrp_grids: list,
+        frequencies: list,
+    ) -> tuple:
+        """Calculate C/I (carrier-to-interference) ratio for multi-site scenarios.
+
+        For each grid point:
+        - C = signal strength from strongest (serving) cell
+        - I = sum of signal strengths from all other co-frequency cells
+        - C/I = C(dBm) - 10*log10(sum of linear interference powers)
+
+        Args:
+            rsrp_grids: List of RSRP arrays, one per site, shape (N,) each
+            frequencies: List of frequencies (MHz) for each site
+
+        Returns:
+            (ci_ratio, best_server_idx, best_rsrp)
+            ci_ratio: C/I in dB, shape (N,)
+            best_server_idx: Index of serving cell per point, shape (N,)
+            best_rsrp: RSRP of serving cell per point, shape (N,)
+        """
+        _xp = gpu_manager.get_array_module()
+
+        if len(rsrp_grids) < 2:
+            # Single site - no interference, return infinity C/I
+            if rsrp_grids:
+                n_points = len(rsrp_grids[0])
+                return (
+                    np.full(n_points, 50.0, dtype=np.float64),  # 50 dB = effectively no interference
+                    np.zeros(n_points, dtype=np.int32),
+                    np.array(rsrp_grids[0], dtype=np.float64),
+                )
+            return np.array([]), np.array([]), np.array([])
+
+        # Stack RSRP grids: shape (num_sites, num_points)
+        rsrp_stack = _xp.stack([_xp.asarray(g, dtype=_xp.float64) for g in rsrp_grids], axis=0)
+        num_sites, num_points = rsrp_stack.shape
+
+        # Convert to linear power (mW)
+        rsrp_linear = _xp.power(10.0, rsrp_stack / 10.0)
+
+        # Best server per point
+        best_server_idx = _xp.argmax(rsrp_stack, axis=0)
+        best_rsrp = _xp.take_along_axis(rsrp_stack, best_server_idx[_xp.newaxis, :], axis=0)[0]
+        best_rsrp_linear = _xp.take_along_axis(rsrp_linear, best_server_idx[_xp.newaxis, :], axis=0)[0]
+
+        # Group sites by frequency for co-channel interference
+        freq_array = _xp.asarray(frequencies, dtype=_xp.float64)
+
+        # Calculate interference only from co-frequency sites
+        interference_linear = _xp.zeros(num_points, dtype=_xp.float64)
+
+        for point_idx in range(num_points):
+            serving_site = int(_to_cpu(best_server_idx[point_idx]))
+            serving_freq = frequencies[serving_site]
+
+            # Sum power from all other sites on same frequency
+            for site_idx in range(num_sites):
+                if site_idx != serving_site and frequencies[site_idx] == serving_freq:
+                    interference_linear[point_idx] += rsrp_linear[site_idx, point_idx]
+
+        # C/I ratio in dB
+        # Avoid log10(0) with small epsilon
+        epsilon = 1e-30
+        ci_ratio = 10 * _xp.log10(best_rsrp_linear / (interference_linear + epsilon))
+
+        # Clip to reasonable range (-20 to 50 dB)
+        ci_ratio = _xp.clip(ci_ratio, -20, 50)
+
+        return (
+            _to_cpu(ci_ratio),
+            _to_cpu(best_server_idx).astype(np.int32),
+            _to_cpu(best_rsrp),
+        )
+
+    def calculate_interference_vectorized(
+        self,
+        rsrp_grids: list,
+        frequencies: list,
+    ) -> tuple:
+        """Fully vectorized C/I calculation (faster for GPU).
+
+        Same as calculate_interference but avoids Python loops.
+        """
+        _xp = gpu_manager.get_array_module()
+
+        if len(rsrp_grids) < 2:
+            if rsrp_grids:
+                n_points = len(rsrp_grids[0])
+                return (
+                    np.full(n_points, 50.0, dtype=np.float64),
+                    np.zeros(n_points, dtype=np.int32),
+                    np.array(rsrp_grids[0], dtype=np.float64),
+                )
+            return np.array([]), np.array([]), np.array([])
+
+        # Stack RSRP grids: shape (num_sites, num_points)
+        rsrp_stack = _xp.stack([_xp.asarray(g, dtype=_xp.float64) for g in rsrp_grids], axis=0)
+        num_sites, num_points = rsrp_stack.shape
+
+        # Convert to linear power (mW)
+        rsrp_linear = _xp.power(10.0, rsrp_stack / 10.0)
+
+        # Best server per point
+        best_server_idx = _xp.argmax(rsrp_stack, axis=0)
+        best_rsrp = _xp.take_along_axis(rsrp_stack, best_server_idx[_xp.newaxis, :], axis=0)[0]
+        best_rsrp_linear = _xp.take_along_axis(rsrp_linear, best_server_idx[_xp.newaxis, :], axis=0)[0]
+
+        # Create frequency match matrix: (num_sites, num_sites)
+        freq_array = _xp.asarray(frequencies, dtype=_xp.float64)
+        freq_match = freq_array[:, _xp.newaxis] == freq_array[_xp.newaxis, :]
+
+        # Total power from all sites
+        total_power = _xp.sum(rsrp_linear, axis=0)
+
+        # For simplified calculation (all sites same frequency):
+        # Interference = total - serving
+        interference_linear = total_power - best_rsrp_linear
+
+        # C/I ratio in dB
+        epsilon = 1e-30
+        ci_ratio = 10 * _xp.log10(best_rsrp_linear / (interference_linear + epsilon))
+
+        # Clip to reasonable range
+        ci_ratio = _xp.clip(ci_ratio, -20, 50)
+
+        return (
+            _to_cpu(ci_ratio),
+            _to_cpu(best_server_idx).astype(np.int32),
+            _to_cpu(best_rsrp),
+        )
+
 
 # Singleton
 gpu_service = GPUService()