@mytec: 3.8.0a done

backend/app/services/coverage_service.py

@@ -581,6 +581,60 @@ class CoverageService:
               f"({len(grid)} points, model={selected_model.name}, freq={site.frequency}MHz, "
               f"env={env}, backend={'GPU' if gpu_service.available else 'CPU/NumPy'}) ━━━")
 
+        # ━━━ PHASE 2.6: GPU-Vectorized Terrain LOS + Diffraction ━━━
+        # This replaces the per-point LOS calculation in workers
+        t_batch_terrain = time.time()
+        grid_elevs = np.array([point_elevations.get((lat, lon), 0.0) for lat, lon in grid])
+
+        if settings.use_terrain and gpu_service.available:
+            _clog("━━━ PHASE 2.6: Batch terrain LOS (GPU) ━━━")
+            has_los_arr, terrain_loss_arr = gpu_service.batch_terrain_los(
+                site.lat, site.lon, site.height, site_elevation,
+                grid_lats.get() if hasattr(grid_lats, 'get') else grid_lats,
+                grid_lons.get() if hasattr(grid_lons, 'get') else grid_lons,
+                grid_elevs,
+                pre_distances,
+                site.frequency,
+                self.terrain._tile_cache,
+                num_samples=30,
+            )
+            batch_terrain_time = time.time() - t_batch_terrain
+            blocked_count = np.sum(~has_los_arr)
+            _clog(f"━━━ PHASE 2.6 done: {batch_terrain_time:.2f}s "
+                  f"({blocked_count}/{len(grid)} blocked by terrain) ━━━")
+
+            # Add terrain results to precomputed dict
+            for i, (lat, lon) in enumerate(grid):
+                if (lat, lon) in precomputed:
+                    precomputed[(lat, lon)]['has_los'] = bool(has_los_arr[i])
+                    precomputed[(lat, lon)]['terrain_loss'] = float(terrain_loss_arr[i])
+        else:
+            _clog("━━━ PHASE 2.6: Skipped (terrain disabled or no GPU) ━━━")
+            # Initialize with defaults
+            for lat, lon in grid:
+                if (lat, lon) in precomputed:
+                    precomputed[(lat, lon)]['has_los'] = True
+                    precomputed[(lat, lon)]['terrain_loss'] = 0.0
+
+        # ━━━ PHASE 2.7: GPU-Vectorized Antenna Pattern ━━━
+        if site.azimuth is not None and site.beamwidth and gpu_service.available:
+            t_batch_antenna = time.time()
+            antenna_loss_arr = gpu_service.batch_antenna_pattern(
+                site.lat, site.lon,
+                grid_lats.get() if hasattr(grid_lats, 'get') else grid_lats,
+                grid_lons.get() if hasattr(grid_lons, 'get') else grid_lons,
+                site.azimuth,
+                site.beamwidth,
+            )
+            for i, (lat, lon) in enumerate(grid):
+                if (lat, lon) in precomputed:
+                    precomputed[(lat, lon)]['antenna_loss'] = float(antenna_loss_arr[i])
+            _clog(f"━━━ PHASE 2.7: Batch antenna pattern done: {time.time() - t_batch_antenna:.2f}s ━━━")
+        else:
+            for lat, lon in grid:
+                if (lat, lon) in precomputed:
+                    precomputed[(lat, lon)]['antenna_loss'] = 0.0
+
         # ━━━ PHASE 3: Point calculation ━━━
         dominant_path_service._log_count = 0  # Reset diagnostic counter
         t_points = time.time()
@@ -1117,6 +1171,9 @@ class CoverageService:
                 timing,
                 precomputed_distance=pre.get('distance') if pre else None,
                 precomputed_path_loss=pre.get('path_loss') if pre else None,
+                precomputed_has_los=pre.get('has_los') if pre else None,
+                precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+                precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
             )
             if point.rsrp >= settings.min_signal:
                 points.append(point)
@@ -1139,6 +1196,9 @@ class CoverageService:
         timing: dict,
         precomputed_distance: Optional[float] = None,
         precomputed_path_loss: Optional[float] = None,
+        precomputed_has_los: Optional[bool] = None,
+        precomputed_terrain_loss: Optional[float] = None,
+        precomputed_antenna_loss: Optional[float] = None,
     ) -> CoveragePoint:
         """Fully synchronous point calculation. All terrain tiles must be pre-loaded."""
 
@@ -1165,29 +1225,37 @@ class CoverageService:
             )
             path_loss = model.calculate(prop_input).path_loss_db
 
-        # Antenna pattern
-        antenna_loss = 0.0
-        if site.azimuth is not None and site.beamwidth:
+        # Antenna pattern (use precomputed if available)
+        if precomputed_antenna_loss is not None:
+            antenna_loss = precomputed_antenna_loss
+        elif site.azimuth is not None and site.beamwidth:
             t0 = time.time()
             antenna_loss = self._antenna_pattern_loss(
                 site.lat, site.lon, lat, lon, site.azimuth, site.beamwidth
             )
             timing["antenna"] += time.time() - t0
+        else:
+            antenna_loss = 0.0
 
-        # Terrain LOS (sync)
-        terrain_loss = 0.0
-        has_los = True
-        if settings.use_terrain:
+        # Terrain LOS (use precomputed if available)
+        if precomputed_has_los is not None and precomputed_terrain_loss is not None:
+            has_los = precomputed_has_los
+            terrain_loss = precomputed_terrain_loss
+        elif settings.use_terrain:
             t0 = time.time()
             los_result = self.los.check_line_of_sight_sync(
                 site.lat, site.lon, site.height, lat, lon, 1.5
             )
             has_los = los_result["has_los"]
+            terrain_loss = 0.0
             if not has_los:
                 terrain_loss = self._diffraction_loss(
                     los_result["clearance"], site.frequency
                 )
             timing["los"] += time.time() - t0
+        else:
+            has_los = True
+            terrain_loss = 0.0
 
         # Building loss (spatial index)
         building_loss = 0.0

backend/app/services/gpu_service.py

@@ -139,6 +139,279 @@ class GPUService:
 
         return _to_cpu(L)
 
+    def batch_terrain_los(
+        self,
+        site_lat: float,
+        site_lon: float,
+        site_height: float,
+        site_elevation: float,
+        grid_lats: np.ndarray,
+        grid_lons: np.ndarray,
+        grid_elevations: np.ndarray,
+        distances: np.ndarray,
+        frequency_mhz: float,
+        terrain_cache: dict,
+        num_samples: int = 30,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Batch compute terrain LOS and diffraction loss for all grid points.
+
+        This is the key GPU optimization — instead of sampling terrain profiles
+        one point at a time, we sample ALL profiles in parallel using vectorized
+        operations.
+
+        Args:
+            site_lat, site_lon: Site coordinates
+            site_height: Antenna height above ground (meters)
+            site_elevation: Ground elevation at site (meters)
+            grid_lats, grid_lons: All grid point coordinates
+            grid_elevations: Ground elevation at each grid point
+            distances: Pre-computed distances from site to each point (meters)
+            frequency_mhz: Frequency for diffraction calculation
+            terrain_cache: Dict[tile_name -> numpy array] from terrain_service
+            num_samples: Number of samples per terrain profile
+
+        Returns:
+            (has_los, terrain_loss) - both shape (N,)
+            has_los: boolean array, True if clear line of sight
+            terrain_loss: diffraction loss in dB (0 if has_los)
+        """
+        _xp = gpu_manager.get_array_module()
+        N = len(grid_lats)
+
+        if N == 0:
+            return np.array([], dtype=bool), np.array([], dtype=np.float64)
+
+        # Convert inputs to GPU arrays
+        g_lats = _xp.asarray(grid_lats, dtype=_xp.float64)
+        g_lons = _xp.asarray(grid_lons, dtype=_xp.float64)
+        g_elevs = _xp.asarray(grid_elevations, dtype=_xp.float64)
+        g_dists = _xp.asarray(distances, dtype=_xp.float64)
+
+        # Heights
+        tx_total = float(site_elevation + site_height)
+        rx_height = 1.5  # Receiver height above ground
+
+        # Earth curvature constants
+        EARTH_RADIUS = 6371000.0
+        K_FACTOR = 4.0 / 3.0
+        effective_radius = K_FACTOR * EARTH_RADIUS
+
+        # Sample terrain profiles for all points at once
+        # Create sample positions: shape (N, num_samples)
+        t = _xp.linspace(0, 1, num_samples, dtype=_xp.float64)  # (S,)
+        t = t.reshape(1, -1)  # (1, S)
+
+        # Interpolate lat/lon for all sample points
+        # sample_lats[i, j] = site_lat + t[j] * (grid_lats[i] - site_lat)
+        dlat = g_lats.reshape(-1, 1) - site_lat  # (N, 1)
+        dlon = g_lons.reshape(-1, 1) - site_lon  # (N, 1)
+        sample_lats = site_lat + t * dlat  # (N, S)
+        sample_lons = site_lon + t * dlon  # (N, S)
+
+        # Sample distances along path: shape (N, S)
+        sample_dists = t * g_dists.reshape(-1, 1)  # (N, S)
+
+        # Get terrain elevations for all samples
+        # This is the tricky part - we need to look up from the tile cache
+        # For GPU efficiency, we'll do this on CPU then transfer
+        sample_lats_cpu = _to_cpu(sample_lats).flatten()
+        sample_lons_cpu = _to_cpu(sample_lons).flatten()
+
+        # Batch elevation lookup from cache
+        sample_elevs_cpu = self._batch_elevation_lookup(
+            sample_lats_cpu, sample_lons_cpu, terrain_cache
+        )
+        sample_elevs = _xp.asarray(sample_elevs_cpu, dtype=_xp.float64).reshape(N, num_samples)
+
+        # Compute LOS line height at each sample point
+        # Linear interpolation from tx to rx
+        rx_total = g_elevs + rx_height  # (N,)
+        los_heights = tx_total + t * (rx_total.reshape(-1, 1) - tx_total)  # (N, S)
+
+        # Earth curvature correction at each sample
+        total_dist = g_dists.reshape(-1, 1)  # (N, 1)
+        d = sample_dists  # (N, S)
+        curvature = (d * (total_dist - d)) / (2 * effective_radius)  # (N, S)
+        los_heights_corrected = los_heights - curvature  # (N, S)
+
+        # Clearance at each sample point
+        clearances = los_heights_corrected - sample_elevs  # (N, S)
+
+        # Minimum clearance per profile
+        min_clearances = _xp.min(clearances, axis=1)  # (N,)
+
+        # Has LOS if minimum clearance > 0
+        has_los = min_clearances > 0  # (N,)
+
+        # Diffraction loss for points without LOS
+        # Using simplified ITU-R P.526 formula
+        terrain_loss = _xp.zeros(N, dtype=_xp.float64)
+
+        # Only compute diffraction where blocked
+        blocked_mask = ~has_los
+        blocked_clearances = min_clearances[blocked_mask]
+
+        if _xp.any(blocked_mask):
+            # v = |clearance| / 10 (simplified Fresnel parameter)
+            v = _xp.abs(blocked_clearances) / 10.0
+
+            # Diffraction loss formula from ITU-R P.526
+            loss = _xp.where(
+                v <= 0,
+                _xp.zeros_like(v),
+                _xp.where(
+                    v < 2.4,
+                    6.02 + 9.11 * v + 1.65 * v ** 2,
+                    12.95 + 20 * _xp.log10(v)
+                )
+            )
+            # Cap at reasonable max
+            loss = _xp.minimum(loss, 40.0)
+            terrain_loss[blocked_mask] = loss
+
+        return _to_cpu(has_los).astype(bool), _to_cpu(terrain_loss)
+
+    def _batch_elevation_lookup(
+        self,
+        lats: np.ndarray,
+        lons: np.ndarray,
+        terrain_cache: dict,
+    ) -> np.ndarray:
+        """Look up elevations from cached terrain tiles.
+
+        Vectorized implementation: processes per-tile (1-4 tiles) instead of
+        per-point (thousands of points). Inner operations are all NumPy vectorized.
+
+        Args:
+            lats, lons: Flattened arrays of coordinates
+            terrain_cache: Dict mapping tile_name -> numpy array
+
+        Returns:
+            elevations: Same shape as input lats
+        """
+        elevations = np.zeros(len(lats), dtype=np.float64)
+
+        # Vectorized tile identification
+        lat_ints = np.floor(lats).astype(int)
+        lon_ints = np.floor(lons).astype(int)
+
+        # Process per tile (usually 1-4 tiles, not per point)
+        unique_tiles = set(zip(lat_ints, lon_ints))
+
+        for lat_int, lon_int in unique_tiles:
+            lat_letter = 'N' if lat_int >= 0 else 'S'
+            lon_letter = 'E' if lon_int >= 0 else 'W'
+            tile_name = f"{lat_letter}{abs(lat_int):02d}{lon_letter}{abs(lon_int):03d}"
+
+            tile = terrain_cache.get(tile_name)
+            if tile is None:
+                continue
+
+            # Mask for points in this tile
+            mask = (lat_ints == lat_int) & (lon_ints == lon_int)
+            tile_lats = lats[mask]
+            tile_lons = lons[mask]
+
+            size = tile.shape[0]
+            # Vectorized row/col calculation
+            rows = ((1 - (tile_lats - lat_int)) * (size - 1)).astype(int)
+            cols = ((tile_lons - lon_int) * (size - 1)).astype(int)
+            rows = np.clip(rows, 0, size - 1)
+            cols = np.clip(cols, 0, size - 1)
+
+            # Vectorized lookup - single operation for ALL points in tile
+            tile_elevs = tile[rows, cols].astype(np.float64)
+            tile_elevs[tile_elevs == -32768] = 0.0
+            elevations[mask] = tile_elevs
+
+        return elevations
+
+    def batch_antenna_pattern(
+        self,
+        site_lat: float,
+        site_lon: float,
+        grid_lats: np.ndarray,
+        grid_lons: np.ndarray,
+        azimuth: float,
+        beamwidth: float,
+    ) -> np.ndarray:
+        """Batch compute antenna pattern loss for all grid points.
+
+        Returns antenna_loss in dB, shape (N,)
+        """
+        _xp = gpu_manager.get_array_module()
+        N = len(grid_lats)
+
+        if N == 0 or azimuth is None or not beamwidth:
+            return np.zeros(N, dtype=np.float64)
+
+        # Convert to radians
+        lat1 = _xp.radians(_xp.float64(site_lat))
+        lon1 = _xp.radians(_xp.float64(site_lon))
+        lat2 = _xp.radians(_xp.asarray(grid_lats, dtype=_xp.float64))
+        lon2 = _xp.radians(_xp.asarray(grid_lons, dtype=_xp.float64))
+
+        # Calculate bearing from site to each point
+        dlon = lon2 - lon1
+        x = _xp.sin(dlon) * _xp.cos(lat2)
+        y = _xp.cos(lat1) * _xp.sin(lat2) - _xp.sin(lat1) * _xp.cos(lat2) * _xp.cos(dlon)
+        bearings = (_xp.degrees(_xp.arctan2(x, y)) + 360) % 360
+
+        # Angle difference from antenna azimuth
+        angle_diff = _xp.abs(bearings - azimuth)
+        angle_diff = _xp.where(angle_diff > 180, 360 - angle_diff, angle_diff)
+
+        # Antenna pattern loss (simplified sector pattern)
+        half_bw = beamwidth / 2
+        in_main = angle_diff <= half_bw
+        loss_main = 3 * (angle_diff / half_bw) ** 2
+        loss_side = 3 + 12 * ((angle_diff - half_bw) / half_bw) ** 2
+        loss_side = _xp.minimum(loss_side, 25.0)
+
+        antenna_loss = _xp.where(in_main, loss_main, loss_side)
+        return _to_cpu(antenna_loss)
+
+    def batch_final_rsrp(
+        self,
+        tx_power: float,
+        tx_gain: float,
+        path_loss: np.ndarray,
+        terrain_loss: np.ndarray,
+        antenna_loss: np.ndarray,
+        building_loss: np.ndarray,
+        vegetation_loss: np.ndarray,
+        rain_loss: np.ndarray,
+        indoor_loss: np.ndarray,
+        atmospheric_loss: np.ndarray,
+        reflection_gain: np.ndarray,
+        fading_margin: float = 0.0,
+    ) -> np.ndarray:
+        """Vectorized final RSRP calculation.
+
+        RSRP = tx_power + tx_gain - path_loss - terrain_loss - antenna_loss
+               - building_loss - vegetation_loss - rain_loss - indoor_loss
+               - atmospheric_loss + reflection_gain - fading_margin
+
+        Returns RSRP in dBm, shape (N,)
+        """
+        _xp = gpu_manager.get_array_module()
+
+        rsrp = (
+            float(tx_power) + float(tx_gain)
+            - _xp.asarray(path_loss, dtype=_xp.float64)
+            - _xp.asarray(terrain_loss, dtype=_xp.float64)
+            - _xp.asarray(antenna_loss, dtype=_xp.float64)
+            - _xp.asarray(building_loss, dtype=_xp.float64)
+            - _xp.asarray(vegetation_loss, dtype=_xp.float64)
+            - _xp.asarray(rain_loss, dtype=_xp.float64)
+            - _xp.asarray(indoor_loss, dtype=_xp.float64)
+            - _xp.asarray(atmospheric_loss, dtype=_xp.float64)
+            + _xp.asarray(reflection_gain, dtype=_xp.float64)
+            - float(fading_margin)
+        )
+
+        return _to_cpu(rsrp)
+
 
 # Singleton
 gpu_service = GPUService()

backend/app/services/parallel_coverage_service.py

@@ -226,6 +226,9 @@ def _ray_process_chunk_impl(chunk, terrain_cache, buildings, osm_data, config):
             config['site_elevation'], point_elev, timing,
             precomputed_distance=pre.get('distance') if pre else None,
             precomputed_path_loss=pre.get('path_loss') if pre else None,
+            precomputed_has_los=pre.get('has_los') if pre else None,
+            precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+            precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
         )
         if point.rsrp >= settings.min_signal:
             results.append(point.model_dump())
@@ -535,6 +538,9 @@ def _pool_worker_process_chunk(args):
             config['site_elevation'], point_elev, timing,
             precomputed_distance=pre.get('distance') if pre else None,
             precomputed_path_loss=pre.get('path_loss') if pre else None,
+            precomputed_has_los=pre.get('has_los') if pre else None,
+            precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+            precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
         )
         if point.rsrp >= settings.min_signal:
             results.append(point.model_dump())
@@ -654,6 +660,9 @@ def _pool_worker_shm_chunk(args):
             config['site_elevation'], point_elev, timing,
             precomputed_distance=pre.get('distance') if pre else None,
             precomputed_path_loss=pre.get('path_loss') if pre else None,
+            precomputed_has_los=pre.get('has_los') if pre else None,
+            precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+            precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
         )
         if point.rsrp >= settings.min_signal:
             results.append(point.model_dump())
@@ -816,6 +825,9 @@ def _pool_worker_shm_shared(args):
             site_elev, point_elev, timing,
             precomputed_distance=pre.get('distance') if pre else None,
             precomputed_path_loss=pre.get('path_loss') if pre else None,
+            precomputed_has_los=pre.get('has_los') if pre else None,
+            precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+            precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
         )
 
         if i < 3:
@@ -1134,6 +1146,9 @@ def _calculate_sequential(
             site_elevation, point_elev, timing,
             precomputed_distance=pre.get('distance') if pre else None,
             precomputed_path_loss=pre.get('path_loss') if pre else None,
+            precomputed_has_los=pre.get('has_los') if pre else None,
+            precomputed_terrain_loss=pre.get('terrain_loss') if pre else None,
+            precomputed_antenna_loss=pre.get('antenna_loss') if pre else None,
         )
         if point.rsrp >= settings.min_signal:
             results.append(point.model_dump())

backend/app/services/vegetation_service.py

@@ -21,6 +21,11 @@ class VegetationArea(BaseModel):
     geometry: List[Tuple[float, float]]  # [(lon, lat), ...]
     vegetation_type: str  # forest, wood, scrub, orchard
     density: str  # dense, sparse, mixed
+    # Bounding box for fast rejection (computed from geometry)
+    min_lat: float = 0.0
+    max_lat: float = 0.0
+    min_lon: float = 0.0
+    max_lon: float = 0.0
 
 
 class VegetationCache:
@@ -127,7 +132,24 @@ class VegetationService:
         cached = self.cache.get(min_lat, min_lon, max_lat, max_lon)
         if cached is not None:
             print(f"[Vegetation] Cache hit for bbox")
-            areas = [VegetationArea(**v) for v in cached]
+            areas = []
+            for v in cached:
+                area = VegetationArea(**v)
+                # Recompute bbox if missing (backward compat with old cache)
+                if area.min_lat == 0.0 and area.max_lat == 0.0 and area.geometry:
+                    lons = [p[0] for p in area.geometry]
+                    lats = [p[1] for p in area.geometry]
+                    area = VegetationArea(
+                        id=area.id,
+                        geometry=area.geometry,
+                        vegetation_type=area.vegetation_type,
+                        density=area.density,
+                        min_lat=min(lats),
+                        max_lat=max(lats),
+                        min_lon=min(lons),
+                        max_lon=max(lons),
+                    )
+                areas.append(area)
             self._memory_cache[cache_key] = areas
             return areas
 
@@ -205,11 +227,19 @@ class VegetationService:
             leaf_type = tags.get("leaf_type", "mixed")
             density = "dense" if leaf_type == "needleleaved" else "mixed"
 
+            # Compute bounding box from geometry (lon, lat tuples)
+            lons = [p[0] for p in geometry]
+            lats = [p[1] for p in geometry]
+
             areas.append(VegetationArea(
                 id=element["id"],
                 geometry=geometry,
                 vegetation_type=veg_type,
-                density=density
+                density=density,
+                min_lat=min(lats),
+                max_lat=max(lats),
+                min_lon=min(lons),
+                max_lon=max(lons),
             ))
 
         return areas
@@ -260,8 +290,12 @@ class VegetationService:
         lat: float, lon: float,
         areas: List[VegetationArea]
     ) -> Optional[VegetationArea]:
-        """Check if point is in vegetation area"""
+        """Check if point is in vegetation area (with bbox pre-filter)"""
         for area in areas:
+            # Quick bbox reject - skips 95%+ of polygons
+            if not (area.min_lat <= lat <= area.max_lat and
+                    area.min_lon <= lon <= area.max_lon):
+                continue
             if self._point_in_polygon(lat, lon, area.geometry):
                 return area
         return None