"""
Calculate probability holes using Monte Carlo integration.

For straight legs, this computes the probability mass from the lateral distribution
that drifts into structures/depths when drifting in each of 8 compass directions.

PERFORMANCE: Uses Monte Carlo integration (10-100x faster than dblquad for indicator-type
integrands) with parallel processing across multiple CPU cores.
"""
from typing import Any, Callable, Optional
import geopandas as gpd
import numpy as np
from shapely.geometry import LineString
import time
import logging
import sys
from concurrent.futures import ThreadPoolExecutor, as_completed
from multiprocessing import cpu_count
import os
from functools import wraps
from dataclasses import dataclass, field
from threading import Lock

# Set up logger
logger = logging.getLogger(__name__)

# =============================================================================
# PROFILING / TIMING INFRASTRUCTURE
# =============================================================================

# Threshold in seconds for logging slow operations
SLOW_THRESHOLD = 2.0

# Global profiling statistics (thread-safe)
_profiling_lock = Lock()


@dataclass
class ProfilingStats:
    """Container for profiling statistics."""
    monte_carlo_calls: int = 0
    monte_carlo_total_time: float = 0.0
    monte_carlo_max_time: float = 0.0
    monte_carlo_slow_count: int = 0  # calls > SLOW_THRESHOLD
    monte_carlo_times: list = field(default_factory=list)  # individual call times

    geometry_ops_calls: int = 0
    geometry_ops_total_time: float = 0.0
    geometry_ops_max_time: float = 0.0

    ray_intersection_calls: int = 0
    ray_intersection_total_time: float = 0.0

    polygon_extraction_calls: int = 0
    polygon_extraction_total_time: float = 0.0

    task_times: list = field(default_factory=list)  # (leg_idx, dir_idx, time)
    task_slow_count: int = 0

    def reset(self):
        """Reset all statistics."""
        self.monte_carlo_calls = 0
        self.monte_carlo_total_time = 0.0
        self.monte_carlo_max_time = 0.0
        self.monte_carlo_slow_count = 0
        self.monte_carlo_times = []
        self.geometry_ops_calls = 0
        self.geometry_ops_total_time = 0.0
        self.geometry_ops_max_time = 0.0
        self.ray_intersection_calls = 0
        self.ray_intersection_total_time = 0.0
        self.polygon_extraction_calls = 0
        self.polygon_extraction_total_time = 0.0
        self.task_times = []
        self.task_slow_count = 0


# Global profiling stats instance
_profiling_stats = ProfilingStats()


def get_profiling_stats() -> ProfilingStats:
    """Get the global profiling statistics."""
    return _profiling_stats


def reset_profiling_stats() -> None:
    """Reset profiling statistics."""
    with _profiling_lock:
        _profiling_stats.reset()


def _record_monte_carlo_time(elapsed: float) -> None:
    """Thread-safe recording of Monte Carlo integration timing."""
    with _profiling_lock:
        _profiling_stats.monte_carlo_calls += 1
        _profiling_stats.monte_carlo_total_time += elapsed
        _profiling_stats.monte_carlo_times.append(elapsed)
        if elapsed > _profiling_stats.monte_carlo_max_time:
            _profiling_stats.monte_carlo_max_time = elapsed
        if elapsed > SLOW_THRESHOLD:
            _profiling_stats.monte_carlo_slow_count += 1
            # Note: logger.warning removed - causes errors in QGIS threads


def _record_task_time(leg_idx: int, dir_idx: int, elapsed: float) -> None:
    """Thread-safe recording of task timing."""
    with _profiling_lock:
        _profiling_stats.task_times.append((leg_idx, dir_idx, elapsed))
        if elapsed > SLOW_THRESHOLD:
            _profiling_stats.task_slow_count += 1
            # Note: logger.warning removed - causes errors in QGIS threads


def _record_geometry_op(elapsed: float) -> None:
    """Thread-safe recording of geometry operation timing."""
    with _profiling_lock:
        _profiling_stats.geometry_ops_calls += 1
        _profiling_stats.geometry_ops_total_time += elapsed
        if elapsed > _profiling_stats.geometry_ops_max_time:
            _profiling_stats.geometry_ops_max_time = elapsed


def _record_ray_intersection(elapsed: float) -> None:
    """Thread-safe recording of ray intersection timing."""
    with _profiling_lock:
        _profiling_stats.ray_intersection_calls += 1
        _profiling_stats.ray_intersection_total_time += elapsed


def _record_polygon_extraction(elapsed: float) -> None:
    """Thread-safe recording of polygon extraction timing."""
    with _profiling_lock:
        _profiling_stats.polygon_extraction_calls += 1
        _profiling_stats.polygon_extraction_total_time += elapsed


def timed_function(func):
    """Decorator to time function execution and log slow calls."""
    @wraps(func)
    def wrapper(*args, **kwargs):
        start = time.perf_counter()
        try:
            result = func(*args, **kwargs)
            return result
        finally:
            elapsed = time.perf_counter() - start
            # Note: logger.warning removed - causes errors in QGIS threads
    return wrapper


def _summary_lines_mc(stats: ProfilingStats) -> list[str]:
    out = ["\n[MONTE CARLO INTEGRATION - Primary computation]",
           f"  Total calls: {stats.monte_carlo_calls}",
           f"  Total time: {stats.monte_carlo_total_time:.2f}s"]
    if stats.monte_carlo_calls > 0:
        avg = stats.monte_carlo_total_time / stats.monte_carlo_calls
        out += [f"  Average time per call: {avg:.4f}s",
                f"  Max single call time: {stats.monte_carlo_max_time:.4f}s",
                f"  Slow calls (>{SLOW_THRESHOLD}s): {stats.monte_carlo_slow_count}"]
        if stats.monte_carlo_times:
            t = np.array(stats.monte_carlo_times)
            out += [f"  Time distribution:",
                    f"    Min: {t.min():.4f}s",
                    f"    Median: {np.median(t):.4f}s",
                    f"    95th percentile: {np.percentile(t, 95):.4f}s",
                    f"    99th percentile: {np.percentile(t, 99):.4f}s",
                    f"    Max: {t.max():.4f}s"]
    return out


def _summary_lines_geometry(stats: ProfilingStats) -> list[str]:
    out = ["\n[GEOMETRY OPERATIONS]",
           f"  Distance checks: {stats.geometry_ops_calls}",
           f"  Total time: {stats.geometry_ops_total_time:.2f}s"]
    if stats.geometry_ops_calls > 0:
        out += [f"  Average time: {stats.geometry_ops_total_time/stats.geometry_ops_calls:.6f}s",
                f"  Max time: {stats.geometry_ops_max_time:.6f}s"]
    out += ["\n[RAY-POLYGON INTERSECTION]",
            f"  Total calls: {stats.ray_intersection_calls}",
            f"  Total time: {stats.ray_intersection_total_time:.2f}s"]
    if stats.ray_intersection_calls > 0:
        out.append(f"  Average time: {stats.ray_intersection_total_time/stats.ray_intersection_calls:.8f}s")
    out += ["\n[POLYGON EXTRACTION]",
            f"  Total calls: {stats.polygon_extraction_calls}",
            f"  Total time: {stats.polygon_extraction_total_time:.2f}s"]
    if stats.polygon_extraction_calls > 0:
        out.append(f"  Average time: {stats.polygon_extraction_total_time/stats.polygon_extraction_calls:.6f}s")
    return out


def _summary_lines_tasks(stats: ProfilingStats) -> list[str]:
    out = ["\n[TASK-LEVEL TIMING]",
           f"  Total tasks: {len(stats.task_times)}",
           f"  Slow tasks (>{SLOW_THRESHOLD}s): {stats.task_slow_count}"]
    if stats.task_times:
        durs = [t[2] for t in stats.task_times]
        total = sum(durs)
        out += [f"  Total task time: {total:.2f}s",
                f"  Average task time: {total/len(durs):.4f}s",
                f"  Max task time: {max(durs):.4f}s",
                f"  Top 5 slowest tasks:"]
        for i, (li, di, dur) in enumerate(sorted(stats.task_times, key=lambda x: x[2], reverse=True)[:5]):
            out.append(f"    {i+1}. Leg {li}, Dir {di}: {dur:.2f}s")
    return out


def _summary_lines_breakdown(stats: ProfilingStats) -> list[str]:
    total = (stats.monte_carlo_total_time + stats.geometry_ops_total_time
             + stats.ray_intersection_total_time + stats.polygon_extraction_total_time)
    out = ["\n[TIME BREAKDOWN]"]
    if total > 0:
        out += [f"  monte carlo: {stats.monte_carlo_total_time:.2f}s ({100*stats.monte_carlo_total_time/total:.1f}%)",
                f"  geometry ops: {stats.geometry_ops_total_time:.2f}s ({100*stats.geometry_ops_total_time/total:.1f}%)",
                f"  ray intersection: {stats.ray_intersection_total_time:.2f}s ({100*stats.ray_intersection_total_time/total:.1f}%)",
                f"  polygon extraction: {stats.polygon_extraction_total_time:.2f}s ({100*stats.polygon_extraction_total_time/total:.1f}%)"]
    return out


def print_profiling_summary() -> str:
    """Generate a human-readable profiling summary."""
    stats = _profiling_stats
    lines = ["\n" + "=" * 80, "PROFILING SUMMARY - BOTTLENECK IDENTIFICATION", "=" * 80]
    lines.extend(_summary_lines_mc(stats))
    lines.extend(_summary_lines_geometry(stats))
    lines.extend(_summary_lines_tasks(stats))
    lines.extend(_summary_lines_breakdown(stats))
    lines.append("=" * 80 + "\n")
    summary = "\n".join(lines)
    logger.info(summary)
    return summary


def _log(message: str) -> None:
    """Log message to logger only (minimal console spam)."""
    logger.info(message)


def _extract_polygon_rings(geom) -> list[np.ndarray]:
    """
    Extract all rings from a Polygon or MultiPolygon as NumPy coordinate arrays.

    Args:
        geom: Shapely Polygon or MultiPolygon

    Returns:
        List of NumPy arrays, each of shape (n_points, 2) representing a ring
    """
    start_time = time.perf_counter()

    from shapely.geometry import Polygon, MultiPolygon

    rings = []

    if isinstance(geom, Polygon):
        # Extract exterior ring
        rings.append(np.array(geom.exterior.coords))
        # Extract interior rings (holes)
        for interior in geom.interiors:
            rings.append(np.array(interior.coords))
    elif isinstance(geom, MultiPolygon):
        # Handle each polygon in the MultiPolygon
        for poly in geom.geoms:
            rings.append(np.array(poly.exterior.coords))
            for interior in poly.interiors:
                rings.append(np.array(interior.coords))
    else:
        # Fallback for other geometry types
        logger.warning(f"Unexpected geometry type: {type(geom)}")

    elapsed = time.perf_counter() - start_time
    _record_polygon_extraction(elapsed)

    return rings


def _compute_probability_monte_carlo(
    leg_start: np.ndarray,
    leg_vec: np.ndarray,
    leg_length: float,
    perp_dir: np.ndarray,
    drift_vec: np.ndarray,
    distance: float,
    lateral_range: float,
    polygon_rings: list[np.ndarray],
    dists: list,
    weights: np.ndarray,
    n_samples: int = 500,
) -> float:
    """
    Monte Carlo integration for probability hole calculation.
    10-100x faster than dblquad for this type of integrand.

    Args:
        leg_start: Start point of the leg [x, y]
        leg_vec: Vector from leg start to leg end
        leg_length: Length of the leg in meters
        perp_dir: Perpendicular direction for lateral offsets
        drift_vec: Unit vector for drift direction
        distance: Maximum drift distance in meters
        lateral_range: Lateral range in meters (typically ±5σ)
        polygon_rings: List of polygon rings as coordinate arrays
        dists: List of scipy.stats distributions for lateral offset
        weights: Normalized weights for each distribution
        n_samples: Number of Monte Carlo samples

    Returns:
        Probability value between 0 and 1
    """
    rng = np.random.default_rng()

    # Generate random samples in the integration domain
    # s in [0, 1], y in [-lateral_range, lateral_range]
    s_samples = rng.random(n_samples)
    y_samples = rng.uniform(-lateral_range, lateral_range, n_samples)

    # Compute positions along leg with lateral offset
    positions = (leg_start
                 + np.outer(s_samples, leg_vec)
                 + np.outer(y_samples, perp_dir))

    # Compute end positions after drifting
    end_positions = positions + distance * drift_vec

    # Check which rays hit the polygon (batch vectorized)
    hits = _batch_ray_intersects_polygon(positions, end_positions, polygon_rings)

    if not np.any(hits):
        return 0.0

    # Compute PDF values only for hits
    y_hits = y_samples[hits]
    pdf_values = np.zeros(len(y_hits))
    for weight, dist in zip(weights, dists):
        pdf_values += weight * dist.pdf(y_hits)

    # Monte Carlo estimate
    # Integral = mean(f(x)) * domain_area
    # domain_area = 1 * 2*lateral_range
    domain_area = 2.0 * lateral_range
    probability = np.sum(pdf_values) / n_samples * domain_area

    return float(np.clip(probability, 0.0, 1.0))


def _ray_intersects_polygon(ray_start: np.ndarray, ray_end: np.ndarray,
                           polygon_rings: list[np.ndarray]) -> bool:
    """
    Fast ray-polygon intersection test using ray-crossing algorithm.

    This is much faster than creating LineString objects and using Shapely's intersects().
    Uses the ray-casting algorithm: count how many times the ray crosses polygon edges.

    Args:
        ray_start: Start point of ray [x, y]
        ray_end: End point of ray [x, y]
        polygon_rings: List of polygon rings (exterior + holes), each as Nx2 array

    Returns:
        True if ray intersects any part of the polygon
    """
    start_time = time.perf_counter()

    # Quick bounding box check first
    ray_min = np.minimum(ray_start, ray_end)
    ray_max = np.maximum(ray_start, ray_end)

    result = False
    for ring in polygon_rings:
        if len(ring) < 2:
            continue

        # Check if ray bounding box overlaps with ring bounding box
        ring_min = ring.min(axis=0)
        ring_max = ring.max(axis=0)

        if (ray_max[0] < ring_min[0] or ray_min[0] > ring_max[0] or
            ray_max[1] < ring_min[1] or ray_min[1] > ring_max[1]):
            continue  # No overlap, skip this ring

        # Test ray-segment intersection for each edge in the ring
        # Ray parameterized as: P(t) = ray_start + t * (ray_end - ray_start), t ∈ [0,1]
        # Edge parameterized as: Q(s) = p1 + s * (p2 - p1), s ∈ [0,1]

        ray_dir = ray_end - ray_start
        ray_len_sq = np.dot(ray_dir, ray_dir)

        if ray_len_sq == 0:
            continue  # Degenerate ray

        # Vectorized edge intersection testing
        p1 = ring[:-1]  # All points except last
        p2 = ring[1:]   # All points except first

        # For each edge, solve for intersection
        edge_dir = p2 - p1  # Shape: (n_edges, 2)
        edge_to_ray = ray_start - p1  # Shape: (n_edges, 2)

        # Cross products for 2D
        # Standard segment-segment intersection (Wikipedia / Ericson):
        #   P(t) = ray_start + t * ray_dir   (t ∈ [0,1])
        #   Q(s) = p1 + s * edge_dir         (s ∈ [0,1])
        #
        # edge_to_ray = ray_start - p1, so (p1 - ray_start) = -edge_to_ray.
        # Correct formulas:
        #   t = -( edge_to_ray × edge_dir ) / ( ray_dir × edge_dir )
        #   s = -( edge_to_ray × ray_dir  ) / ( ray_dir × edge_dir )
        cross_ray_edge = ray_dir[0] * edge_dir[:, 1] - ray_dir[1] * edge_dir[:, 0]

        # Avoid division by zero (parallel lines)
        valid = np.abs(cross_ray_edge) > 1e-10

        if not np.any(valid):
            continue

        # Calculate intersection parameters
        cross_to_edge = edge_to_ray[:, 0] * edge_dir[:, 1] - edge_to_ray[:, 1] * edge_dir[:, 0]
        cross_to_ray = edge_to_ray[:, 0] * ray_dir[1] - edge_to_ray[:, 1] * ray_dir[0]

        t = np.where(valid, -cross_to_edge / cross_ray_edge, -1)
        s = np.where(valid, -cross_to_ray / cross_ray_edge, -1)

        # Check if intersection occurs within both segments
        intersects = valid & (t >= 0) & (t <= 1) & (s >= 0) & (s <= 1)

        if np.any(intersects):
            result = True
            break

    elapsed = time.perf_counter() - start_time
    _record_ray_intersection(elapsed)

    return result


def _batch_ray_intersects_polygon(
    ray_starts: np.ndarray,  # Shape: (N, 2)
    ray_ends: np.ndarray,    # Shape: (N, 2)
    polygon_rings: list[np.ndarray],
) -> np.ndarray:
    """
    Batch ray-polygon intersection test - tests multiple rays simultaneously.

    Returns:
        Boolean array of shape (N,) indicating which rays intersect
    """
    n_rays = len(ray_starts)
    hits = np.zeros(n_rays, dtype=bool)

    ray_dirs = ray_ends - ray_starts

    for ring in polygon_rings:
        if len(ring) < 2:
            continue

        # Quick bounding box rejection
        ring_min = ring.min(axis=0)
        ring_max = ring.max(axis=0)

        ray_mins = np.minimum(ray_starts, ray_ends)
        ray_maxs = np.maximum(ray_starts, ray_ends)

        # Rays that might intersect this ring
        possible = ((ray_maxs[:, 0] >= ring_min[0]) &
                    (ray_mins[:, 0] <= ring_max[0]) &
                    (ray_maxs[:, 1] >= ring_min[1]) &
                    (ray_mins[:, 1] <= ring_max[1]) &
                    ~hits)  # Skip already-hit rays

        if not np.any(possible):
            continue

        # Test against all edges
        p1s = ring[:-1]
        p2s = ring[1:]
        edge_dirs = p2s - p1s

        for ray_idx in np.where(possible)[0]:
            ray_start = ray_starts[ray_idx]
            ray_dir = ray_dirs[ray_idx]

            edge_to_ray = ray_start - p1s
            cross_ray_edge = ray_dir[0] * edge_dirs[:, 1] - ray_dir[1] * edge_dirs[:, 0]

            valid = np.abs(cross_ray_edge) > 1e-10
            if not np.any(valid):
                continue

            cross_to_edge = edge_to_ray[:, 0] * edge_dirs[:, 1] - edge_to_ray[:, 1] * edge_dirs[:, 0]
            cross_to_ray = edge_to_ray[:, 0] * ray_dir[1] - edge_to_ray[:, 1] * ray_dir[0]

            t = np.where(valid, -cross_to_edge / cross_ray_edge, -1)
            s = np.where(valid, -cross_to_ray / cross_ray_edge, -1)

            if np.any(valid & (t >= 0) & (t <= 1) & (s >= 0) & (s <= 1)):
                hits[ray_idx] = True

    return hits


def _setup_leg_vars(
    line: LineString,
) -> tuple[np.ndarray, np.ndarray, float, np.ndarray]:
    """Return (leg_start, leg_vec, leg_len, perp_dir) for a leg geometry."""
    coords = np.array(line.coords)
    leg_start, leg_end = coords[0], coords[-1]
    leg_vec = leg_end - leg_start
    leg_len = line.length
    leg_dir = leg_vec / leg_len if leg_len > 0 else np.array([1.0, 0.0])
    perp_dir = np.array([-leg_dir[1], leg_dir[0]])
    return leg_start, leg_vec, leg_len, perp_dir


def _compute_lateral_ranges(
    distributions: list, weights: list,
) -> list[float]:
    result = []
    for dists, wgts in zip(distributions, weights):
        w = np.array(wgts)
        if w.sum() == 0:
            w = np.ones_like(w)
        w = w / w.sum()
        weighted_std = float(np.sqrt(sum(wt * (d.std() ** 2) for d, wt in zip(dists, w))))
        result.append(5.0 * weighted_std)
    return result


def _build_tasks(
    lines: list, distributions: list, weights: list,
    drift_angles: list, lateral_ranges: list,
    objs_gdf_list: list, obj_index_map: list, distance: float,
) -> list:
    tasks = []
    for leg_idx, (line, dists, wgts) in enumerate(zip(lines, distributions, weights)):
        coords = np.array(line.coords)
        if len(coords) < 2:
            continue
        for dir_idx, angle_deg in enumerate(drift_angles):
            tasks.append((
                leg_idx, dir_idx, angle_deg, line, dists, wgts,
                objs_gdf_list, obj_index_map, distance, lateral_ranges[leg_idx],
            ))
    return tasks


def _log_computation_start(n_lines: int, n_dirs: int, n_objs: int, n_total: int, max_workers: int) -> None:
    _log(f"\n{'='*80}")
    _log("COMPUTING PROBABILITY HOLES (PARALLEL OPTIMIZED)")
    _log(f"{'='*80}")
    _log(f"Legs: {n_lines} | Drift directions: {n_dirs} | Objects: {n_objs}")
    _log(f"Total hole calculations: {n_total} | Worker threads: {max_workers}")
    _log(f"{'='*80}")


def _log_computation_end(
    total_holes: int, skipped_full: int, skipped_far: int, start_time: float,
) -> None:
    total_time = time.time() - start_time
    actually_computed = total_holes - skipped_full - skipped_far
    _log(f"\n{'='*80}")
    _log("HOLE CALCULATION COMPLETE")
    _log(f"{'='*80}")
    _log(f"Total: {total_holes} | Computed: {actually_computed} | Skipped-coverage: {skipped_full} | Skipped-far: {skipped_far}")
    _log(f"Total time: {int(total_time/60)}m {int(total_time%60)}s ({total_time:.1f}s)")
    if actually_computed > 0:
        _log(f"Avg time/computed hole: {total_time/actually_computed:.2f}s")
    _log(f"{'='*80}\n")


def _build_obj_index_map(objs_gdf_list: list) -> list[tuple[int, int]]:
    """Flatten list of GeoDataFrames into (gdf_idx, row_idx) pairs."""
    obj_index_map: list[tuple[int, int]] = []
    for gi, gdf in enumerate(objs_gdf_list):
        for ri in range(len(gdf)):
            obj_index_map.append((gi, ri))
    return obj_index_map


def _check_progress_cancel(
    progress_callback: Callable,
    completed_tasks: int,
    total_tasks: int,
    start_time: float,
    last_pct: int,
    per_leg_dir_obj: list,
    futures: dict,
) -> tuple[int, bool]:
    """Invoke progress callback and cancel futures if user requests stop. Returns (new_last_pct, should_cancel)."""
    current = int((completed_tasks / total_tasks) * 100)
    if current <= last_pct:
        return last_pct, False
    elapsed = time.time() - start_time
    avg = elapsed / completed_tasks if completed_tasks > 0 else 0
    eta_s = avg * (total_tasks - completed_tasks)
    msg = f"Progress: {completed_tasks}/{total_tasks} tasks | ETA: {int(eta_s/60)}m {int(eta_s%60)}s"
    if not progress_callback(completed_tasks, total_tasks, msg):
        logger.info("Calculation cancelled by user")
        for f in futures:
            f.cancel()
        return current, True
    return current, False


def _compute_single_direction(
    leg_idx: int,
    dir_idx: int,
    angle_deg: float,
    line: LineString,
    dists: list[Any],
    wgts: list[float],
    objs_gdf_list: list[gpd.GeoDataFrame],
    obj_index_map: list[tuple[int, int]],
    distance: float,
    lateral_range: float,
) -> tuple[int, int, list[float], int, int]:
    """
    Worker function to compute probability holes for one (leg, direction) combination.

    This function is designed to be called by parallel workers.

    Args:
        leg_idx: Index of the leg
        dir_idx: Index of the direction
        angle_deg: Drift angle in degrees
        line: LineString geometry of the leg
        dists: Distribution list for this leg
        wgts: Weight list for this leg
        objs_gdf_list: List of GeoDataFrames containing objects
        obj_index_map: Flattened object indexing
        distance: Maximum drift distance in meters
        lateral_range: Lateral range in meters (±5σ)

    Returns:
        Tuple of (leg_idx, dir_idx, dir_holes, skipped_full_coverage, skipped_too_far)
    """
    task_start = time.perf_counter()
    w = np.array(wgts)
    if w.sum() == 0: w = np.ones_like(w)
    w = w / w.sum()
    leg_start, leg_vec, leg_len, perp_dir = _setup_leg_vars(line)
    angle_rad = np.radians(angle_deg)
    drift_vec = np.array([np.cos(angle_rad), np.sin(angle_rad)])
    polygon_coords_dict = {
        (gi, ri): _extract_polygon_rings(objs_gdf_list[gi].geometry.iloc[ri])
        for gi, ri in obj_index_map}
    dir_holes: list[float] = []
    skipped_too_far = 0
    for obj_idx, (gi, ri) in enumerate(obj_index_map):
        obj = objs_gdf_list[gi].geometry.iloc[ri]
        polygon_rings = polygon_coords_dict[(gi, ri)]
        geom_t = time.perf_counter()
        min_dist = line.distance(obj)
        _record_geometry_op(time.perf_counter() - geom_t)
        if min_dist > distance + lateral_range:
            dir_holes.append(0.0)
            skipped_too_far += 1
            continue
        try:
            mc_t = time.perf_counter()
            p = _compute_probability_monte_carlo(
                leg_start, leg_vec, leg_len, perp_dir, drift_vec,
                distance, lateral_range, polygon_rings, dists, w, n_samples=2000,
            )
            _record_monte_carlo_time(time.perf_counter() - mc_t)
            dir_holes.append(max(0.0, min(1.0, p)))
        except Exception as e:
            logger.debug(f"Monte Carlo failed for leg {leg_idx}, dir {dir_idx}, obj {obj_idx}: {e}")
            dir_holes.append(0.0)
    _record_task_time(leg_idx, dir_idx, time.perf_counter() - task_start)
    return (leg_idx, dir_idx, dir_holes, 0, skipped_too_far)


def compute_probability_holes(
    lines: list[LineString],
    distributions: list[list[Any]],
    weights: list[list[float]],
    objs_gdf_list: list[gpd.GeoDataFrame],
    distance: float,
    progress_callback: Optional[Callable[[int, int, str], bool]] = None,
) -> list[list[list[float]]]:
    """
    Calculate probability holes for all legs, drift directions, and objects.

    Uses Monte Carlo integration with geometric ray intersection to compute
    the probability that vessels drifting from the leg hit each object.
    Monte Carlo is 10-100x faster than dblquad for indicator-type integrands.

    Args:
        lines: List of LineString geometries representing traffic legs
        distributions: List of distribution lists (one per leg, each containing scipy.stats distributions)
        weights: List of weight lists (one per leg, weights for each distribution)
        objs_gdf_list: List of GeoDataFrames containing object geometries (structures or depths)
        distance: Maximum drift distance in meters
        progress_callback: Optional callback function(completed, total, message) -> bool.
                          Returns False to cancel, True to continue.

    Returns:
        3-level nested list: [leg_idx][direction_idx][object_idx] = probability (0-1)
        where direction_idx corresponds to: 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°
    """
    reset_profiling_stats()
    obj_index_map = _build_obj_index_map(objs_gdf_list)
    drift_angles = [0, 45, 90, 135, 180, 225, 270, 315]
    n_objs = len(obj_index_map)
    total_holes = len(lines) * len(drift_angles) * n_objs
    skipped_full_coverage, skipped_too_far = 0, 0
    start_time = time.time()
    max_workers = max(1, cpu_count() - 1)
    _log_computation_start(len(lines), len(drift_angles), n_objs, total_holes, max_workers)
    lateral_ranges = _compute_lateral_ranges(distributions, weights)
    per_leg_dir_obj = [[[0.0] * n_objs for _ in drift_angles] for _ in lines]
    tasks = _build_tasks(
        lines, distributions, weights, drift_angles, lateral_ranges,
        objs_gdf_list, obj_index_map, distance,
    )
    total_tasks, completed_tasks, last_pct = len(tasks), 0, -1
    with ThreadPoolExecutor(max_workers=max_workers) as executor:
        futures = {executor.submit(_compute_single_direction, *t): t for t in tasks}
        for future in as_completed(futures):
            try:
                leg_idx, dir_idx, dir_holes, skip_cov, skip_far = future.result()
                per_leg_dir_obj[leg_idx][dir_idx] = dir_holes
                skipped_full_coverage += skip_cov
                skipped_too_far += skip_far
                completed_tasks += 1
                if progress_callback:
                    last_pct, cancel = _check_progress_cancel(
                        progress_callback, completed_tasks, total_tasks,
                        start_time, last_pct, per_leg_dir_obj, futures,
                    )
                    if cancel:
                        return per_leg_dir_obj
            except Exception as e:
                logger.error(f"Task failed with error: {e}")
    _log_computation_end(total_holes, skipped_full_coverage, skipped_too_far, start_time)
    print_profiling_summary()
    return per_leg_dir_obj
