"""
/***************************************************************************
 Intersection
                                 A QGIS plugin
 Topaze
                             -------------------
        begin                : 2025-08-28
        git sha              : $Format:%H$
        copyright            : (C) 2025 by Jean-Marie ARSAC
        email                : jmarsac@arsac.wf
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""

import json
import math
import os
from dataclasses import dataclass
from itertools import combinations
from math import atan2, cos, pi, sin, tan

import numpy as np
from qgis.core import QgsProject
from scipy.optimize import least_squares

from topaze.calc.calc_helpers import (
    az_gr_north_cw,
    line_intersection,
    wrap_gr_0_400,
    wrap_gr_m200_200,
)
from topaze.calc.elevation_utils import compute_target_z_from_verticals_intersection
from topaze.file_utils import FileUtils
from topaze.report_helpers import fmt
from topaze.report_utils import ReportUtils  # uses the auto Markdown→HTML/PDF pipeline
from topaze.toolbelt import PlgLogger
from topaze.toolbelt.i18n import tr
from topaze.topaze_calculator import TopazeCalculator
from topaze.topaze_utils import TopazeUtils
from topaze.topo_sight import TopoSight

# -------------------------- data structures --------------------------


@dataclass
class LOS:  # Line Of Sight used for the target
    station_id: str
    xs: float
    ys: float
    ah_gr: float  # measured horizontal angle (gr)
    v0_gr: float  # station orientation (gr)
    dir_abs_gr: float  # absolute direction to target: v0 + ah (gr)


# -------------------------- core computations --------------------------


def _compute_v0_all_stations_weighted(
    stations_list, measures, target_id, verbose: bool = True
):
    """
    Compute per-station V0 as a distance-weighted circular mean of (bearing - AH),
    using ALL references with known coordinates (no special treatment for AH=0).

    Returns {station_id: {"v0":..., "sigma0_gr":..., "n":..., "details":[...]}}
    If verbose=True, prints per-station diagnostics (one block per station).
    """
    per_station_refs = {}
    for m in measures:
        # keep only references with known coords and not the (unknown) intersection target
        if (
            m.id in stations_list
            and m.id_target != target_id
            and m.x_target is not None
            and m.y_target is not None
            and m.ah is not None
        ):
            per_station_refs.setdefault(m.id, []).append(m)

    v0_by = {}
    g2r = pi / 200.0

    for st_id, refs in per_station_refs.items():
        if not refs:
            continue

        if verbose:
            print(
                tr("V0 computation for station {st}: using {n} reference(s)").format(
                    st=st_id, n=len(refs)
                )
            )

        C = 0.0  # sum w*cos()
        S = 0.0  # sum w*sin()
        diffs = []  # (bearing - AH) in grads
        weights = []  # distances

        # 1) Build (diff, weight) per reference and print inputs
        for r in refs:
            bea = TopazeCalculator.bearing_gon(
                r.x, r.y, r.x_target, r.y_target
            )  # grads 0..400
            diff = wrap_gr_m200_200(bea - r.ah)  # (-200,200]
            d = float(np.hypot(r.x_target - r.x, r.y_target - r.y))
            w = d if d > 0.0 else 1.0

            diffs.append(diff)
            weights.append(w)
            C += w * cos(diff * g2r)
            S += w * sin(diff * g2r)

            if verbose:
                print(
                    tr(
                        "  ref {to}: bearing={b:.5f} gr, AH={ah:.5f} gr, diff={d:.5f} gr, "
                        "distance={R:.3f} m, weight={w:.3f}"
                    ).format(to=r.id_target, b=bea, ah=r.ah, d=diff, R=d, w=w)
                )

        # 2) Weighted circular mean
        v0 = wrap_gr_0_400(atan2(S, C) * 200.0 / pi)
        if verbose:
            print(tr("  V0 (weighted mean) = {v0:.5f} gr").format(v0=v0))

        # 3) Residuals around mean and sigma0
        res = [wrap_gr_m200_200(d - v0) for d in diffs]
        dof = max(len(res) - 1, 0)
        sigma0 = None if dof == 0 else float(np.sqrt(np.sum(np.square(res)) / dof))
        if verbose:
            for r, rres in zip(refs, res):
                print(
                    tr("    residual to mean for ref {to}: {rv:.5f} gr").format(
                        to=r.id_target, rv=rres
                    )
                )
            if sigma0 is None:
                print(tr("  σ₀ (gr): —  (DOF=0)"))
            else:
                print(
                    tr("  σ₀ (gr): {s:.5f}  (n={n}, DOF={d})").format(
                        s=sigma0, n=len(res), d=dof
                    )
                )

        # 4) Store results
        v0_by[st_id] = {
            "v0": v0,
            "sigma0_gr": sigma0,
            "n": len(diffs),
            "details": [
                {
                    "ref_to": r.id_target,
                    "bearing_gr": TopazeCalculator.bearing_gon(
                        r.x, r.y, r.x_target, r.y_target
                    ),
                    "ah_gr": r.ah,
                    "residual_gr": wrap_gr_m200_200(
                        (
                            wrap_gr_m200_200(
                                TopazeCalculator.bearing_gon(
                                    r.x, r.y, r.x_target, r.y_target
                                )
                                - r.ah
                            )
                        )
                        - v0
                    ),
                }
                for r in refs
            ],
        }

    if not v0_by and verbose:
        PlgLogger.log(
            message=tr(
                "No station orientations (V0) could be computed from known references."
            ),
            log_level=2,
            push=True,
        )

    return v0_by


def _build_los_list(target_id, stations_dic, measures, v0_by_station):
    """
    Build the list of lines-of-sight to the selected target from stations that have a computed V0.
    """
    los_list = []
    for m in measures:
        if m.id_target != target_id:
            continue
        if m.id not in v0_by_station:
            continue
        st = v0_by_station[m.id]
        v0 = st["v0"]
        dir_abs = wrap_gr_0_400(v0 + (m.ah if m.ah is not None else 0.0))
        los_list.append(
            LOS(m.id, stations_dic[m.id][0], stations_dic[m.id][1], m.ah, v0, dir_abs)
        )
    return los_list


def _raw_intersection_covadis(los_list):
    """
    Covadis-style raw intersection:
      - sort LOS by absolute direction (0..400 g),
      - build cyclic adjacent pairs (i, i+1 mod n) → n intersections,
      - drop the WORST pair (smallest |sin Δθ|) to keep n-1 intersections,
      - return the simple mean of the remaining intersection points.
    """
    n = len(los_list)
    if n < 2:
        raise ValueError(tr("At least 2 sightings are required for an intersection."))
    if n == 2:
        # only one pair; its intersection is both raw and compensated if consistent
        a, b = los_list[0], los_list[1]
        azr = a.dir_abs_gr * pi / 200.0
        bzr = b.dir_abs_gr * pi / 200.0
        P, Q = (a.xs, a.ys), (b.xs, b.ys)
        pt = line_intersection(P, azr, Q, bzr)
        if pt is None:
            raise RuntimeError(
                tr("The two sight lines are parallel or nearly parallel.")
            )
        return pt[0], pt[1]

    los_sorted = sorted(los_list, key=lambda l: l.dir_abs_gr)
    candidates = []  # (cond, x, y, i, j) with cond = |sin(Δθ)|

    for i in range(n):
        j = (i + 1) % n
        A, B = los_sorted[i], los_sorted[j]
        ai = A.dir_abs_gr * pi / 200.0
        bj = B.dir_abs_gr * pi / 200.0
        P, Q = (A.xs, A.ys), (B.xs, B.ys)
        pt = line_intersection(P, ai, Q, bj)
        if pt is None:
            continue
        cond = abs(sin(ai - bj))  # small => nearly parallel (worst)
        candidates.append((cond, pt[0], pt[1], i, j))

    if not candidates:
        raise RuntimeError(tr("All adjacent sight pairs are nearly parallel."))

    # Keep the best k = min(n-1, len(candidates)) intersections: drop the single worst-conditioned
    candidates.sort(key=lambda t: t[0])  # ascending by conditioning
    k = min(n - 1, len(candidates))
    selected = candidates[-k:]  # largest |sin Δθ|

    x_raw = float(sum(x for _, x, _, _, _ in selected) / len(selected))
    y_raw = float(sum(y for _, _, y, _, _ in selected) / len(selected))
    return x_raw, y_raw


def _lsq_intersection(
    los_list,
    x0,
    y0,
    robust_first=True,
    loss="soft_l1",
    f_scale_gr=0.0005,
    reject_threshold_gr=0.5,
):
    """
    2D absolute-directions LSQ with optional robust pre-fit and automatic outlier rejection.
    - Step 1 (optional): robust fit (soft-L1) → residuals
    - Step 2: reject sights with |residual| > reject_threshold_gr
    - Step 3: classic LSQ (LM) on the kept sights
    Returns ((x_hat, y_hat), stats) where stats contains per-sight residuals and EMQs.
    """
    xs = np.array([l.xs for l in los_list], float)
    ys = np.array([l.ys for l in los_list], float)
    z = np.array([l.dir_abs_gr for l in los_list], float)
    n = len(los_list)

    def resid(u):
        x, y = u[0], u[1]
        az = np.arctan2(x - xs, y - ys) * (200.0 / np.pi)
        az = np.mod(az, 400.0)
        v = z - az
        v = (v + 200.0) % 400.0 - 200.0
        return v

    u0 = np.array([x0, y0], float)

    # Step 1: robust pre-fit (optional)
    if robust_first:
        pre = least_squares(
            resid,
            u0,
            method="trf",
            loss=loss,
            f_scale=f_scale_gr,
            xtol=1e-12,
            ftol=1e-12,
            gtol=1e-12,
            max_nfev=300,
        )
        v_pre = resid(pre.x)
        keep_idx = [i for i, rv in enumerate(v_pre) if abs(rv) <= reject_threshold_gr]
        if len(keep_idx) < 2:
            # keep the two best anyway
            keep_idx = sorted(range(n), key=lambda i: abs(v_pre[i]))[:2]
        xs_k = xs[keep_idx]
        ys_k = ys[keep_idx]
        z_k = z[keep_idx]
    else:
        keep_idx = list(range(n))
        xs_k, ys_k, z_k = xs, ys, z

    # Step 3: classic LSQ on kept
    def resid_k(u):
        x, y = u[0], u[1]
        az = np.arctan2(x - xs_k, y - ys_k) * (200.0 / np.pi)
        az = np.mod(az, 400.0)
        v = z_k - az
        v = (v + 200.0) % 400.0 - 200.0
        return v

    lm = least_squares(
        resid_k,
        u0 if not robust_first else pre.x,
        method="lm",
        xtol=1e-12,
        ftol=1e-12,
        gtol=1e-12,
    )
    x_hat, y_hat = float(lm.x[0]), float(lm.x[1])

    # Residuals on kept and full lists (for reporting)
    v_k = resid_k(np.array([x_hat, y_hat], float))
    # propagate back to the full list (mark rejected as None)
    v_full = [None] * n
    for j, i in enumerate(keep_idx):
        v_full[i] = float(v_k[j])

    dof = max(len(keep_idx) - 2, 0)
    rmse = float(np.sqrt(np.mean(v_k**2))) if len(v_k) else 0.0
    sigma0 = None if dof == 0 else float(np.sqrt(np.sum(v_k**2) / dof))

    # Covariance / EMQ
    std_x = std_y = None
    cov = None
    if sigma0 is not None and hasattr(lm, "jac") and lm.jac is not None:
        J = lm.jac
        JTJ = J.T @ J
        try:
            JTJ_inv = np.linalg.inv(JTJ)
            cov = (sigma0**2) * JTJ_inv
            std_x = float(np.sqrt(cov[0, 0]))
            std_y = float(np.sqrt(cov[1, 1]))
        except np.linalg.LinAlgError:
            pass

    # Per-sight metric equivalents (for kept only; others are marked rejected)
    per_sight = []
    for i, l in enumerate(los_list):
        if v_full[i] is None:
            per_sight.append({"station": l.station_id, "rejected": True})
            continue
        az_gr = az_gr_north_cw(x_hat - l.xs, y_hat - l.ys)
        eps_gr = v_full[i]
        R = float(np.hypot(x_hat - l.xs, y_hat - l.ys))
        eps_rad = eps_gr * np.pi / 200.0
        trans_m = R * np.tan(eps_rad)
        az_rad = az_gr * np.pi / 200.0
        dE = trans_m * np.cos(az_rad)
        dN = -trans_m * np.sin(az_rad)
        per_sight.append(
            {
                "station": l.station_id,
                "v0_gr": l.v0_gr,
                "ah_gr": l.ah_gr,
                "dir_abs_gr": l.dir_abs_gr,
                "computed_gr": az_gr,
                "residual_gr": eps_gr,
                "range_m": R,
                "transverse_m": trans_m,
                "transverse_mm": trans_m * 1000.0,
                "dE_m": dE,
                "dN_m": dN,
                "rejected": False,
            }
        )

    stats = {
        "n": n,
        "kept": [los_list[i].station_id for i in keep_idx],
        "rejected": [los_list[i].station_id for i in range(n) if i not in keep_idx],
        "dof": dof,
        "rmse_gr": rmse,
        "sigma0_gr": sigma0,
        "std_x_m": std_x,
        "std_y_m": std_y,
        "cov": cov.tolist() if cov is not None else None,
        "residuals": per_sight,
    }
    return (x_hat, y_hat), stats, keep_idx


# -------------------------- reporting --------------------------
# --- inside export_markdown_report(...) ---


def add_z_block_intersection(L, z_raw, z_cmp, zstats):
    """Append the 'Elevation (Z)' block for an intersection target."""
    L.append("## " + tr("Elevation (Z)"))
    L.append("")

    if z_raw is None and z_cmp is None:
        L.append(tr("No vertical angles available to compute target elevation."))
        L.append("")
        return

    # Summary line
    L.append(
        "- {}: {} / {} {}".format(
            tr("Target elevation (raw / compensated)"),
            fmt(z_raw, "{:.4f}"),
            fmt(z_cmp, "{:.4f}"),
            tr("m"),
        )
    )
    if zstats:
        n = zstats.get("n")
        dof = zstats.get("dof")
        emq = zstats.get("EMQ_Z_m")
        rm = zstats.get("rmse_gr")
        s0 = zstats.get("sigma0_gr")

        # NEW: apparent level line
        ap = zstats.get("apparent_level", {}) or {}
        enabled = bool(ap.get("enabled", False))
        k_ref = ap.get("k_refraction", None)
        earth_R = ap.get("earth_radius_m", None)
        onoff = tr("ON") if enabled else tr("OFF")
        line = "- {}: {}".format(tr("Apparent level correction"), onoff)
        if k_ref is not None and earth_R is not None:
            line += f" (k={fmt(k_ref, '{:.3f}')}, R={fmt(earth_R, '{:.0f}')} {tr('m')})"
        L.append(line)

        L.append(
            "- {}: {}  —  {}: {}".format(
                tr("Observations"),
                fmt(n, "{:d}"),
                tr("DOF"),
                fmt(dof, "{:d}"),
            )
        )
        L.append(
            "- {} (m): {}  —  RMSE (gr): {}  —  σ₀ (gr): {}".format(
                tr("EMQ Z (1σ)"),
                fmt(emq, "{:.4f}"),
                fmt(rm),
                fmt(s0),
            )
        )
        L.append("")

        residuals = zstats.get("residuals") or []
        if residuals:
            # NEW: add C (m) column if present
            show_C = any("C_m" in r for r in residuals)
            L.append("### " + tr("Per-sight verticals"))
            L.append("")
            if show_C:
                L.append(
                    "| "
                    + tr("Station")
                    + " | "
                    + tr("Distance (m)")
                    + " | C (m) | "
                    + tr("AV measured (gr)")
                    + " | "
                    + tr("AV computed (gr)")
                    + " | "
                    + tr("Residual (gr)")
                    + " | "
                    + tr("Z estimate (m)")
                    + " | "
                    + tr("Z residual (m)")
                    + " |"
                )
                L.append("|---|---:|---:|---:|---:|---:|---:|---:|")
            else:
                L.append(
                    "| "
                    + tr("Station")
                    + " | "
                    + tr("Distance (m)")
                    + " | "
                    + tr("AV measured (gr)")
                    + " | "
                    + tr("AV computed (gr)")
                    + " | "
                    + tr("Residual (gr)")
                    + " | "
                    + tr("Z estimate (m)")
                    + " | "
                    + tr("Z residual (m)")
                    + " |"
                )
                L.append("|---|---:|---:|---:|---:|---:|---:|")

            for r in residuals:
                if show_C:
                    L.append(
                        f"| {r.get('station','')} | "
                        f"{fmt(r.get('r_m'), '{:.3f}')} | {fmt(r.get('C_m'), '{:.4f}')} | "
                        f"{fmt(r.get('av_meas_gr'))} | {fmt(r.get('av_comp_gr'))} | "
                        f"{fmt(r.get('residual_gr'))} | "
                        f"{fmt(r.get('z_est_m'), '{:.4f}')} | {fmt(r.get('residual_z_m'), '{:.4f}')} |"
                    )
                else:
                    L.append(
                        f"| {r.get('station','')} | "
                        f"{fmt(r.get('r_m'), '{:.3f}')} | "
                        f"{fmt(r.get('av_meas_gr'))} | {fmt(r.get('av_comp_gr'))} | "
                        f"{fmt(r.get('residual_gr'))} | "
                        f"{fmt(r.get('z_est_m'), '{:.4f}')} | {fmt(r.get('residual_z_m'), '{:.4f}')} |"
                    )
            L.append("")


def _export_markdown_report(
    target_id,
    stations_dic,
    v0_by_station,
    los_list,
    raw_xy,
    lsq_xy,
    lsq_stats,
    z_raw=None,
    z_cmp=None,
    zstats=None,
):
    L = []
    L.append("# " + tr("Intersection computing"))
    L.append("")
    L.append("**{k}:** {v}".format(k=tr("Target ID"), v=target_id))
    L.append("")
    L.append("## " + tr("Sight lines used"))
    L.append("")
    L.append(
        "| "
        + tr("Station")
        + " | X | Y | "
        + tr("V0 (gr)")
        + " | "
        + tr("AH (gr)")
        + " | "
        + tr("Absolute dir (gr)")
        + " |"
    )
    L.append("|---|---:|---:|---:|---:|---:|")
    for l in los_list:
        L.append(
            f"| {l.station_id} | {l.xs:.3f} | {l.ys:.3f} | {l.v0_gr:.5f} | {l.ah_gr:.5f} | {l.dir_abs_gr:.5f} |"
        )
    L.append("")

    L.append("## " + tr("Station orientations (V0)"))
    L.append("")
    L.append(
        "| "
        + tr("Station")
        + " | "
        + tr("V0 (gr)")
        + " | "
        + tr("n refs")
        + " | σ₀ (gr) |"
    )
    L.append("|---|---:|---:|---:|")
    for st, info in v0_by_station.items():
        s0 = "—" if info["sigma0_gr"] is None else f"{info['sigma0_gr']:.5f}"
        L.append(f"| {st} | {info['v0']:.5f} | {info['n']} | {s0} |")
    L.append("")

    # Results
    L.append("## " + tr("Results"))
    L.append("")
    L.append("| " + tr("Method") + " | X | Y |")
    L.append("|---|---:|---:|")
    L.append(
        "| {} | {:.4f} | {:.4f} |".format(
            tr("Raw (pairwise)"),
            raw_xy[0],
            raw_xy[1],
        )
    )
    if lsq_xy is not None:
        L.append(
            "| {} | {:.4f} | {:.4f} |".format(
                tr("Compensated (Least Squares)"),
                lsq_xy[0],
                lsq_xy[1],
            )
        )
    L.append("")

    # LSQ block
    if lsq_stats is not None:
        s0 = (
            "—"
            if lsq_stats["sigma0_gr"] is None
            else "{:.5f}".format(lsq_stats["sigma0_gr"])
        )
    L.append("## " + tr("Least Squares summary"))
    L.append("")
    L.append(
        "- {}: {}  —  {}: {}".format(
            tr("Observations"),
            lsq_stats["n"],
            tr("DOF"),
            lsq_stats["dof"],
        )
    )
    L.append(
        "- RMSE (gr): {:.5f}  —  σ₀ (gr): {}".format(
            lsq_stats["rmse_gr"],
            s0,
        )
    )
    if lsq_stats["std_x_m"] is not None and lsq_stats["std_y_m"] is not None:
        L.append(
            "- {} (m): {:.4f}  —  {} (m): {:.4f}".format(
                tr("EMQ X (1σ)"),
                lsq_stats["std_x_m"],
                tr("EMQ Y (1σ)"),
                lsq_stats["std_y_m"],
            )
        )
        L.append("")
        # per-sight residuals (angular and metric)
        L.append("### " + tr("Per-sight residuals"))
        L.append("")
        L.append(
            "| "
            + tr("Station")
            + " | "
            + tr("Abs dir (gr)")
            + " | "
            + tr("Computed (gr)")
            + " | "
            + tr("Residual (gr)")
            + " | "
            + tr("Distance (m)")
            + " | "
            + tr("Transverse (mm)")
            + " |"
        )
        L.append("|---|---:|---:|---:|---:|---:|")
        for r in lsq_stats["residuals"]:
            L.append(
                f"| {r['station']} | {r['dir_abs_gr']:.5f} | {r['computed_gr']:.5f} | {r['residual_gr']:.5f} | {r['range_m']:.3f} | {r['transverse_mm']:.1f} |"
            )
        L.append("")
        add_z_block_intersection(L, z_raw, z_cmp, zstats)

    md = "\n".join(L)

    return md


# -------------------------- public API --------------------------


class Intersection:
    @staticmethod
    def compute_intersection():
        """
        Root function.
        Reads 'intersection.json', filters useful observations, computes:
          - station V0 (compensated if multiple refs),
          - raw XY from pairwise LOS,
          - compensated XY by LSQ (if >= 3 LOS),
        builds & exports a Markdown/HTML/PDF report,
        and returns ((x_raw, y_raw), (x_cmp, y_cmp), stats_dict).
        """
        data = FileUtils.load_temp_file("intersection.json")
        if not data:
            PlgLogger.log(
                message=tr("No data from intersection.json"),
                log_level=2,
                push=True,
            )
            out_dir = QgsProject.instance().homePath() + "/rapport/"
            md_name = "resection_report.md"
            FileUtils.remove_temp_file(md_name, out_dir)
            return (None, None, None), (None, None, None), None

        obs_data = json.loads(data)

        corrections = obs_data.get("corrections", {})  # may be empty

        target_id = obs_data["calcul"][0]
        # Build stations dictionary
        stations_dic = {}
        for st in obs_data["stations"]:
            stations_dic[st["matricule"]] = (st["x"], st["y"])

        # Build TopoSight list (we keep only H angles that have both station & target coords known)
        measures: list[TopoSight] = []
        for o in obs_data["obs"]:
            m = None
            if (
                o.get("ah", None) is not None
                or o.get("av", None) is not None
                or o.get("di", None) is not None
            ):
                m = TopoSight(
                    id=o["origine"],
                    id_target=o["cible"],
                    hi=o["hi"],
                    ah=o["ah"],
                    av=o["av"],
                    di=o["di"],
                    hp=o["hp"],
                )
                measures.append(m)
        for m in measures:
            for station in obs_data["stations"]:
                if station["matricule"] == m.id:
                    m.x = station["x"]
                    m.y = station["y"]
                    m.z = station["z"]
                if station["matricule"] == m.id_target:
                    m.x_target = station["x"]
                    m.y_target = station["y"]
                    m.z_target = station["z"]
            """
            print(
                f"obs is id: {m.id} x: {m.x} y: {m.y} hi: {m.hi} "
                f"id_target: {m.id_target} x_target: {m.x_target} "
                f"y_target: {m.y_target} ah: {m.ah} hp: {m.hp}"
            )
            """

        # 1) V0 per station (based on references != target_id)
        stations_list = []
        for st in obs_data["stations"]:
            stations_list.append(st["matricule"])
        v0_by_station = _compute_v0_all_stations_weighted(
            stations_list, measures, target_id
        )

        # 2) Build LOS to target from stations with V0
        los_list = _build_los_list(target_id, stations_dic, measures, v0_by_station)
        if len(los_list) < 2:
            PlgLogger.log(
                message=tr(
                    "Not enough sight lines to compute intersection (need at least 2)."
                ),
                log_level=2,
                push=True,
            )
            return (None, None, None), (None, None, None), None

        # 3) Raw (Covadis-style n-1 adjacent mean)
        x_raw, y_raw = _raw_intersection_covadis(los_list)
        # print(tr("Raw intersection: ({x}, {y})").format(x=x_raw, y=y_raw))

        # 4) LSQ: directions-only adjustment (if >=3)
        x_cmp = y_cmp = None
        stats = None
        if len(los_list) >= 3:
            (x_cmp, y_cmp), stats, keep_idx = _lsq_intersection(
                los_list, x_raw, y_raw, robust_first=False
            )
            # print(tr("Compensated intersection: ({x}, {y})").format(x=x_cmp, y=y_cmp))
        else:
            # With exactly 2 lines, LSQ is exactly determined: use raw as compensated
            x_cmp, y_cmp = x_raw, y_raw
            stats = {
                "n": 2,
                "dof": 0,
                "rmse_gr": 0.0,
                "sigma0_gr": None,
                "std_x_m": None,
                "std_y_m": None,
                "cov": None,
                "residuals": [],
            }

        # 4-bis ) Z computation from verticals (if any)
        # Build station dict with Z for verticals
        stations_by_id = {
            st["matricule"]: {"x": st["x"], "y": st["y"], "z": st.get("z")}
            for st in obs_data["stations"]
        }

        # Choose XY for Z (prefer compensated)
        xt, yt = (
            (x_cmp, y_cmp)
            if (x_cmp is not None and y_cmp is not None)
            else (x_raw, y_raw)
        )

        z_raw = z_cmp = None
        zstats = None
        if xt is not None and yt is not None:
            z_raw, z_cmp, zstats = compute_target_z_from_verticals_intersection(
                target_id,
                xt,
                yt,
                measures,
                stations_by_id,
                verbose=True,
                corrections=corrections,
            )
            if z_cmp is not None:
                print(
                    tr("Target elevation (raw/comp): {zr:.4f} / {zc:.4f} m").format(
                        zr=z_raw, zc=z_cmp
                    )
                )

        # 5) Report
        md_report = _export_markdown_report(
            target_id,
            stations_dic,
            v0_by_station,
            los_list,
            (x_raw, y_raw),
            (x_cmp, y_cmp),
            stats,
            z_raw=z_raw,
            z_cmp=z_cmp,
            zstats=zstats,
        )
        out_dir = QgsProject.instance().homePath() + "/rapport/"
        os.makedirs(out_dir, exist_ok=True)
        md_name = "intersection_report.md"
        FileUtils.save_temp_file(md_name, md_report, out_dir)
        md_path = os.path.join(out_dir, md_name)
        try:
            ReportUtils.markdown_to_html(md_path, title=tr("Intersection Report"))
        except Exception as ex:
            PlgLogger.log(
                message=tr("HTML export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )
        try:
            ReportUtils.markdown_to_pdf(
                md_path, title=tr("Intersection Report"), margins_mm=(7, 15, 7, 15)
            )
        except Exception as ex:
            PlgLogger.log(
                message=tr("PDF export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )
        # Return ((raw), (cmp), stats) as asked
        stats["zstats"] = zstats
        return (x_raw, y_raw, z_raw), (x_cmp, y_cmp, z_cmp), stats
