"""
Comparative Evaluation: Tool-Assisted (MathSolver Core) vs No-Tool (Pure LLM)
Benchmark: 165 Geometry Questions (2016-2026 National Exams)
Model: Google Gemma 4 26B A4B MoE (gemini/gemma-4-26b-a4b-it) in Low Thinking mode.

Metrics Evaluated:
1. Math Answer Accuracy (% exact/symbolic math match with ground_truth_answer)
2. Geometry Constraint Satisfaction Rate (CSR % against GT verification_conditions)
3. Strict Geometric Validity (% passed GeometryValidator with 0 errors)
4. Oxyz Coordinate Error (MAE / MSE against canonical_coordinates)
5. Latency Distribution (Mean, Median, P95)
"""

from __future__ import annotations

import asyncio
import json
import logging
import math
import os
import re
import sys
import time
from typing import Dict, Any, List, Optional, Tuple
import numpy as np
import sympy as sp
from sympy.parsing.sympy_parser import parse_expr

logging.basicConfig(level=logging.WARNING)
logger = logging.getLogger(__name__)

sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))

from llm.service import get_llm_service
from agents.deepmath_solver_agent import DeepMathSolverAgent
from solver.dsl_parser import DSLParser
from solver.engine import GeometryEngine
from solver.validator import GeometryValidator
from solver.models import Point, Constraint


DATASET_PATH = os.path.abspath(
    os.path.join(os.path.dirname(__file__), "../../research/geometry_dataset/all_geometry_questions.jsonl")
)
TOOL_GEOM_RESULTS_PATH = os.path.abspath(
    os.path.join(os.path.dirname(__file__), "full_benchmark_gemma26b_results.json")
)
OUTPUT_RESULTS_PATH = os.path.abspath(
    os.path.join(os.path.dirname(__file__), "benchmark_tool_vs_notool_results.json")
)
OUTPUT_METRICS_PATH = os.path.abspath(
    os.path.join(os.path.dirname(__file__), "benchmark_tool_vs_notool_metrics.json")
)

SEMAPHORE_LIMIT = 4
MODEL_NAME = "gemini/gemma-4-26b-a4b-it"


def normalize_str(s: str) -> str:
    if not s:
        return ""
    s = str(s).strip().lower()
    s = s.replace(";", ",").replace("°", "").replace("degrees", "")
    s = re.sub(r"\s+", "", s)
    return s


def check_line_equivalence(s1: str, s2: str) -> bool:
    """Checks whether two 3D line representations (symmetric or parametric) are mathematically equivalent."""
    try:
        def get_dir(s: str) -> Optional[List[float]]:
            # 1. Parametric: x = x0 + at, y = y0 + bt, z = z0 + ct
            t_coeffs = re.findall(r"([-+]?\s*\d*)t", s)
            if t_coeffs and len(t_coeffs) == 3:
                res = []
                for c in t_coeffs:
                    c = c.replace(" ", "")
                    if c in ("", "+"):
                        res.append(1.0)
                    elif c == "-":
                        res.append(-1.0)
                    else:
                        res.append(float(c))
                return res
            # 2. LaTeX fraction denominators: \frac{...}{denom}
            frac_denoms = re.findall(r"frac\{[^{}]*\}\{([-+]?\d+)\}", s)
            if len(frac_denoms) == 3:
                return [float(d) for d in frac_denoms]
            # 3. Standard slash denominators: / denom
            denoms = re.findall(r"/\s*\(?([-+]?\d+)\)?", s)
            if len(denoms) == 3:
                return [float(d) for d in denoms]
            elif len(denoms) == 2:
                return [1.0, float(denoms[0]), float(denoms[1])]
            return None

        d1 = get_dir(s1)
        d2 = get_dir(s2)
        if d1 and d2:
            ratios = [a / b if b != 0 else (999 if a != 0 else 0) for a, b in zip(d1, d2)]
            if abs(ratios[0] - ratios[1]) < 1e-3 and abs(ratios[1] - ratios[2]) < 1e-3 and abs(ratios[0]) > 1e-3:
                return True
    except Exception:
        pass
    return False


def match_answer(pred: str, gt: str) -> bool:
    """Robust multi-tier mathematical answer equivalence checker."""
    if not pred or not gt:
        return False

    norm_p = normalize_str(pred)
    norm_g = normalize_str(gt)

    if norm_p == norm_g:
        return True

    # Check 3D line equivalence (symmetric vs parametric)
    if check_line_equivalence(pred, gt):
        return True

    # Check substring containment for complex multi-part answers (e.g. (P) and H)
    if len(norm_p) >= 3 and (norm_p in norm_g or norm_g in norm_p):
        return True

    val_p = norm_p.split("=")[-1].split(":")[-1]
    val_g = norm_g.split("=")[-1].split(":")[-1]
    if val_p == val_g:
        return True

    # Numeric float comparison
    try:
        num_p = float(re.sub(r"[^\d.-]", "", val_p))
        num_g = float(re.sub(r"[^\d.-]", "", val_g))
        if abs(num_p - num_g) < 1e-3:
            return True
    except Exception:
        pass

    # Coordinate tuples extraction (e.g. (1, 1, 6) vs (1, 1, 6))
    m_p = re.findall(r"[-+]?\d*\.?\d+", val_p)
    m_g = re.findall(r"[-+]?\d*\.?\d+", val_g)
    if m_p and m_g and len(m_p) == len(m_g) and len(m_p) >= 2:
        try:
            if all(abs(float(a) - float(b)) < 1e-3 for a, b in zip(m_p, m_g)):
                return True
        except Exception:
            pass

    # SymPy symbolic equivalence with positive variable assumption (for volume parameters)
    try:
        sym_scope = {
            "a": sp.Symbol("a", positive=True),
            "b": sp.Symbol("b", positive=True),
            "c": sp.Symbol("c", positive=True),
            "pi": sp.pi,
            "sqrt": sp.sqrt,
            "sp": sp,
        }
        clean_p = val_p.replace("^", "**").replace(r"\sqrt", "sqrt").replace(r"\pi", "pi")
        clean_g = val_g.replace("^", "**").replace(r"\sqrt", "sqrt").replace(r"\pi", "pi").replace("v=", "").replace("r=", "").replace("s=", "")
        sym_p = parse_expr(clean_p, local_dict=sym_scope)
        sym_g = parse_expr(clean_g, local_dict=sym_scope)
        diff = sp.simplify(sym_p - sym_g)
        if diff == 0:
            return True
    except Exception:
        pass

    return False


def parse_coord_val(val: Any) -> Optional[float]:
    if isinstance(val, (int, float)):
        return float(val)
    if isinstance(val, str):
        val = val.strip()
        try:
            return float(val)
        except ValueError:
            try:
                res = sp.sympify(val)
                return float(res.evalf())
            except Exception:
                return None
    return None


def evaluate_coordinates(
    pred_coords: Dict[str, Any],
    geom_gt: Dict[str, Any],
    is_3d: bool = True,
) -> Dict[str, Any]:
    """Validates predicted coordinates against GT verification conditions and canonical coordinates."""
    if not pred_coords or not isinstance(pred_coords, dict):
        return {
            "is_valid": False,
            "checked_count": 0,
            "error_count": 0,
            "csr_pct": 0.0,
            "oxyz_mae": None,
            "oxyz_mse": None,
            "errors": ["No coordinates provided"],
        }

    # Clean coordinates to float lists of length 3
    cleaned_coords: Dict[str, List[float]] = {}

    # Check if format is {"x": ..., "y": ..., "z": ...} for a single target point
    if all(k in pred_coords for k in ["x", "y"]):
        target_pts = geom_gt.get("required_objects", {}).get("points", ["A"])
        target_name = target_pts[0] if target_pts else "P"
        vals = [parse_coord_val(pred_coords.get("x")), parse_coord_val(pred_coords.get("y")), parse_coord_val(pred_coords.get("z", 0.0))]
        if all(v is not None for v in vals):
            cleaned_coords[target_name] = [vals[0], vals[1], vals[2] or 0.0]

    for k, v in pred_coords.items():
        if k in ["x", "y", "z"]:
            continue
        pt_name = str(k).strip()
        if isinstance(v, (list, tuple)) and len(v) >= 2:
            vals = [parse_coord_val(x) for x in v]
            if all(x is not None for x in vals[:3]):
                if len(vals) == 2:
                    vals.append(0.0)
                cleaned_coords[pt_name] = vals[:3]
        elif isinstance(v, dict) and "x" in v and "y" in v:
            vals = [parse_coord_val(v.get("x")), parse_coord_val(v.get("y")), parse_coord_val(v.get("z", 0.0))]
            if all(x is not None for x in vals):
                cleaned_coords[pt_name] = [vals[0], vals[1], vals[2] or 0.0]

    if not cleaned_coords:
        return {
            "is_valid": False,
            "checked_count": 0,
            "error_count": 0,
            "csr_pct": 0.0,
            "oxyz_mae": None,
            "oxyz_mse": None,
            "errors": ["Failed to parse numerical coordinates"],
        }

    # Reconstruct constraints from verification_conditions
    v_conditions = geom_gt.get("verification_conditions", [])
    constraints: List[Constraint] = []
    for c in v_conditions:
        constraints.append(
            Constraint(
                type=c.get("type", "length"),
                targets=c.get("targets", []),
                value=c.get("value", 0.0),
            )
        )

    # Run GeometryValidator
    solved_ir = {
        "coordinates": cleaned_coords,
        "segments": [],
        "faces": [],
    }
    validator = GeometryValidator()
    val_res = validator.validate(solved_ir, constraints, is_3d=is_3d)
    checked_count = val_res.checked_count
    error_count = len(val_res.errors)

    csr_pct = 100.0
    if checked_count > 0:
        csr_pct = round(max(0.0, (checked_count - error_count) / checked_count * 100.0), 1)

    # Compute Oxyz coordinate error if canonical coordinates exist
    canon_coords = geom_gt.get("canonical_coordinates", {})
    common_pts = [p for p in canon_coords if p in cleaned_coords]
    oxyz_mae = None
    oxyz_mse = None

    if common_pts:
        dists = []
        for p in common_pts:
            c_pred = np.array(cleaned_coords[p])
            c_gt = np.array(canon_coords[p])
            dists.append(float(np.linalg.norm(c_pred - c_gt)))
        oxyz_mae = round(float(np.mean(dists)), 4)
        oxyz_mse = round(float(np.mean([d ** 2 for d in dists])), 4)

    return {
        "is_valid": val_res.is_valid,
        "checked_count": checked_count,
        "error_count": error_count,
        "csr_pct": csr_pct,
        "oxyz_mae": oxyz_mae,
        "oxyz_mse": oxyz_mse,
        "errors": val_res.errors[:3],
    }


def parse_notool_json(raw: str) -> Tuple[Optional[str], Dict[str, Any]]:
    """Extracts final_answer and coordinates from No-Tool raw LLM response."""
    if not raw or not isinstance(raw, str):
        return None, {}

    cleaned = raw.strip()
    m = re.search(r"```(?:json)?(.*?)```", cleaned, re.DOTALL)
    if m:
        cleaned = m.group(1).strip()
    else:
        m2 = re.search(r"(\{.*\})", cleaned, re.DOTALL)
        if m2:
            cleaned = m2.group(1).strip()

    safe_dict = {
        "__builtins__": None,
        "math": math,
        "sqrt": math.sqrt,
        "pi": math.pi,
        "true": True,
        "false": False,
        "null": None,
    }

    data = None
    try:
        data = json.loads(cleaned)
    except Exception:
        try:
            data = eval(cleaned, safe_dict, {})
        except Exception:
            data = None

    ans = None
    coords = {}
    if isinstance(data, dict):
        for k in ["final_answer", "answer", "result", "dap_an", "ket_qua"]:
            if k in data and data[k]:
                ans = str(data[k]).strip()
                break
        for ck in ["coordinates", "toado", "points", "toa_do"]:
            if ck in data and isinstance(data[ck], dict):
                coords = data[ck]
                break

    # Regex fallbacks if not found
    if not ans:
        ans_m = re.search(r'"(?:final_answer|answer|result|dap_an)"\s*:\s*"([^"]+)"', raw)
        if ans_m:
            ans = ans_m.group(1).strip()
    if not coords:
        coords_m = re.search(r'"(?:coordinates|toado|points)"\s*:\s*(\{[^{}]*?(?:\{[^{}]*?\}[^{}]*?)*\})', raw, re.DOTALL)
        if coords_m:
            try:
                coords = eval(coords_m.group(1), safe_dict, {})
            except Exception:
                pass

    return ans, coords


async def evaluate_single_question(
    q: Dict[str, Any],
    tool_geom_cache: Dict[str, Any],
    llm_svc: Any,
    solver_agent: DeepMathSolverAgent,
    sem: asyncio.Semaphore,
    progress: Dict[str, Any],
) -> Dict[str, Any]:
    qid = q["id"]
    year = q["year"]
    topic = q["topic"]
    subtopic = q.get("subtopic", "")
    difficulty = q.get("difficulty", "medium")
    raw_text = q.get("raw_text", "").strip()
    gt_answer = q.get("answer", "").strip()
    geom_gt = q.get("ground_truth_geometry", {})

    result: Dict[str, Any] = {
        "id": qid,
        "year": year,
        "topic": topic,
        "subtopic": subtopic,
        "difficulty": difficulty,
        "raw_text": raw_text,
        "ground_truth_answer": gt_answer,
        "no_tool": {},
        "tool": {},
    }

    async with sem:
        async def run_notool():
            notool_prompt = f"""Bạn là một chuyên gia toán hình học. Hãy giải bài toán hình học sau ĐỘC LẬP (tuyệt đối không dùng tool, code hay công cụ ngoài).
YÊU CẦU: Giải thật NGẮN GỌN, đi thẳng vào kết quả và toạ độ. Trả về DUY NHẤT một khối JSON hợp lệ theo đúng cấu trúc sau (không kèm văn bản giải thích thừa ngoài JSON):
{{
  "final_answer": "<đáp án toán học cuối cùng: số, phân số, phương trình mp/đt, toạ độ điểm>",
  "coordinates": {{
    "A": [x, y, z],
    "B": [x, y, z]
  }}
}}

Đề bài:
{raw_text}"""
            t0 = time.time()
            try:
                notool_raw = await llm_svc.acomplete(
                    model=MODEL_NAME,
                    messages=[{"role": "user", "content": notool_prompt}],
                    temperature=0.1,
                    max_tokens=3072,
                    timeout=150,
                    reasoning_effort="low",
                    thinking_budget=0,
                )
                t_notool = round((time.time() - t0) * 1000, 1)

                notool_ans, notool_coords = parse_notool_json(notool_raw)
                notool_ans_match = match_answer(notool_ans or "", gt_answer)
                notool_geom_eval = evaluate_coordinates(notool_coords, geom_gt, is_3d=True)

                result["no_tool"] = {
                    "status": "OK",
                    "latency_ms": t_notool,
                    "predicted_answer": notool_ans,
                    "answer_match": notool_ans_match,
                    "predicted_coordinates": notool_coords,
                    "geometry_eval": notool_geom_eval,
                }
            except Exception as e:
                t_notool = round((time.time() - t0) * 1000, 1)
                result["no_tool"] = {
                    "status": "ERROR",
                    "latency_ms": t_notool,
                    "error": str(e),
                    "answer_match": False,
                    "geometry_eval": {"is_valid": False, "csr_pct": 0.0},
                }

        async def run_tool():
            t0 = time.time()
            try:
                # 2a. Solve math answer using DeepMathSolverAgent with SymPy sandbox execution
                solver_res = await solver_agent.solve(problem_text=raw_text)
                t_tool_math = round((time.time() - t0) * 1000, 1)

                tool_ans = solver_res.get("answer") or ""
                tool_ans_match = match_answer(tool_ans, gt_answer)

                # 2b. Retrieve / Validate Tool Geometry from MathSolver Core Benchmark
                cached_geom = tool_geom_cache.get(qid, {})
                tool_geom_stages = cached_geom.get("stages", {})
                engine_stage = tool_geom_stages.get("geometry_engine", {})
                val_stage = tool_geom_stages.get("validation", {})

                tool_coords = engine_stage.get("coordinates", {})
                tool_geom_eval = evaluate_coordinates(tool_coords, geom_gt, is_3d=True)

                result["tool"] = {
                    "status": "OK",
                    "latency_ms": t_tool_math + engine_stage.get("latency_ms", 0.0),
                    "math_latency_ms": t_tool_math,
                    "engine_latency_ms": engine_stage.get("latency_ms", 0.0),
                    "predicted_answer": tool_ans,
                    "answer_match": tool_ans_match,
                    "predicted_coordinates": tool_coords,
                    "geometry_eval": tool_geom_eval,
                }
            except Exception as e:
                t_tool = round((time.time() - t0) * 1000, 1)
                result["tool"] = {
                    "status": "ERROR",
                    "latency_ms": t_tool,
                    "error": str(e),
                    "answer_match": False,
                    "geometry_eval": {"is_valid": False, "csr_pct": 0.0},
                }

        await asyncio.gather(run_notool(), run_tool())

        progress["completed"] += 1
        _log_comparison(progress, qid, year, result)
        return result


def _log_comparison(p: Dict[str, Any], qid: str, year: int, res: Dict[str, Any]):
    nt = res["no_tool"]
    tl = res["tool"]

    nt_ans_icon = "🎯" if nt.get("answer_match") else "❌"
    tl_ans_icon = "🎯" if tl.get("answer_match") else "❌"

    nt_geo_icon = "📐" if nt.get("geometry_eval", {}).get("is_valid") else "⚠️"
    tl_geo_icon = "📐" if tl.get("geometry_eval", {}).get("is_valid") else "⚠️"

    print(
        f"[{p['completed']:3d}/{p['total']}] {year} | {qid[:22]:<22} | "
        f"NoTool: Ans {nt_ans_icon} Geo {nt_geo_icon} (CSR: {nt.get('geometry_eval', {}).get('csr_pct', 0.0):5.1f}%) | "
        f"Tool: Ans {tl_ans_icon} Geo {tl_geo_icon} (CSR: {tl.get('geometry_eval', {}).get('csr_pct', 0.0):5.1f}%)"
    )


def calculate_metrics(results: List[Dict[str, Any]]) -> Dict[str, Any]:
    n = len(results)
    if n == 0:
        return {}

    # 1. Answer Accuracy
    notool_ans_correct = sum(1 for r in results if r.get("no_tool", {}).get("answer_match"))
    tool_ans_correct = sum(1 for r in results if r.get("tool", {}).get("answer_match"))

    # 2. Geometry Strict Validity
    notool_geo_valid = sum(1 for r in results if r.get("no_tool", {}).get("geometry_eval", {}).get("is_valid"))
    tool_geo_valid = sum(1 for r in results if r.get("tool", {}).get("geometry_eval", {}).get("is_valid"))

    # 3. Constraint Satisfaction Rates
    notool_csrs = [r.get("no_tool", {}).get("geometry_eval", {}).get("csr_pct", 0.0) for r in results]
    tool_csrs = [r.get("tool", {}).get("geometry_eval", {}).get("csr_pct", 0.0) for r in results]

    # 4. Oxyz Coordinate MAE
    notool_maes = [
        r["no_tool"]["geometry_eval"]["oxyz_mae"]
        for r in results
        if r.get("topic") == "oxyz" and r.get("no_tool", {}).get("geometry_eval", {}).get("oxyz_mae") is not None
    ]
    tool_maes = [
        r["tool"]["geometry_eval"]["oxyz_mae"]
        for r in results
        if r.get("topic") == "oxyz" and r.get("tool", {}).get("geometry_eval", {}).get("oxyz_mae") is not None
    ]

    # 5. Latency Stats
    notool_lats = [r["no_tool"]["latency_ms"] for r in results if "latency_ms" in r.get("no_tool", {})]
    tool_lats = [r["tool"]["latency_ms"] for r in results if "latency_ms" in r.get("tool", {})]

    def lat_summary(arr):
        if not arr:
            return {}
        a = np.array(arr)
        return {
            "mean_ms": round(float(np.mean(a)), 1),
            "median_ms": round(float(np.median(a)), 1),
            "p95_ms": round(float(np.percentile(a, 95)), 1),
        }

    # Breakdown by topic
    topics = sorted(list(set(r.get("topic", "unknown") for r in results)))
    by_topic = {}
    for t in topics:
        t_items = [r for r in results if r.get("topic") == t]
        t_n = len(t_items)
        by_topic[t] = {
            "total": t_n,
            "notool_math_acc_pct": round(sum(1 for r in t_items if r["no_tool"].get("answer_match")) / t_n * 100, 1),
            "tool_math_acc_pct": round(sum(1 for r in t_items if r["tool"].get("answer_match")) / t_n * 100, 1),
            "notool_geo_valid_pct": round(sum(1 for r in t_items if r["no_tool"].get("geometry_eval", {}).get("is_valid")) / t_n * 100, 1),
            "tool_geo_valid_pct": round(sum(1 for r in t_items if r["tool"].get("geometry_eval", {}).get("is_valid")) / t_n * 100, 1),
            "notool_avg_csr_pct": round(float(np.mean([r["no_tool"].get("geometry_eval", {}).get("csr_pct", 0.0) for r in t_items])), 1),
            "tool_avg_csr_pct": round(float(np.mean([r["tool"].get("geometry_eval", {}).get("csr_pct", 0.0) for r in t_items])), 1),
        }

    return {
        "total_questions": n,
        "summary": {
            "no_tool": {
                "math_accuracy_count": notool_ans_correct,
                "math_accuracy_pct": round(notool_ans_correct / n * 100, 2),
                "geometric_validity_count": notool_geo_valid,
                "geometric_validity_pct": round(notool_geo_valid / n * 100, 2),
                "avg_constraint_satisfaction_pct": round(float(np.mean(notool_csrs)), 2),
                "oxyz_coordinate_mae": round(float(np.mean(notool_maes)), 4) if notool_maes else None,
                "latency": lat_summary(notool_lats),
            },
            "tool": {
                "math_accuracy_count": tool_ans_correct,
                "math_accuracy_pct": round(tool_ans_correct / n * 100, 2),
                "geometric_validity_count": tool_geo_valid,
                "geometric_validity_pct": round(tool_geo_valid / n * 100, 2),
                "avg_constraint_satisfaction_pct": round(float(np.mean(tool_csrs)), 2),
                "oxyz_coordinate_mae": round(float(np.mean(tool_maes)), 4) if tool_maes else None,
                "latency": lat_summary(tool_lats),
            },
        },
        "by_topic": by_topic,
    }


async def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("--limit", type=int, default=None, help="Number of questions to evaluate")
    parser.add_argument("--start-idx", type=int, default=0, help="Starting question index")
    parser.add_argument("--concurrency", type=int, default=1, help="Concurrency limit")
    parser.add_argument("--output-results", type=str, default=None, help="Output results path")
    parser.add_argument("--output-metrics", type=str, default=None, help="Output metrics path")
    args, _ = parser.parse_known_args()

    results_path = os.path.abspath(args.output_results) if args.output_results else OUTPUT_RESULTS_PATH
    metrics_path = os.path.abspath(args.output_metrics) if args.output_metrics else OUTPUT_METRICS_PATH
    sem_limit = args.concurrency

    with open(DATASET_PATH, "r", encoding="utf-8") as f:
        all_questions = [json.loads(l) for l in f]

    if args.limit is not None:
        target_questions = all_questions[args.start_idx : args.start_idx + args.limit]
    else:
        target_questions = all_questions

    print(f"=== COMPARATIVE BENCHMARK: TOOL vs NO-TOOL (N={len(target_questions)}) ===")
    print(f"Model: {MODEL_NAME} Low Thinking (budget=0)")
    print(f"Dataset: {DATASET_PATH}")
    print(f"Concurrency: {sem_limit}, Results Path: {results_path}\n")

    # Load tool geometry cache
    tool_geom_cache: Dict[str, Any] = {}
    if os.path.exists(TOOL_GEOM_RESULTS_PATH):
        with open(TOOL_GEOM_RESULTS_PATH, "r", encoding="utf-8") as f:
            t_items = json.load(f)
            for item in t_items:
                if "id" in item:
                    tool_geom_cache[item["id"]] = item
        print(f"Loaded {len(tool_geom_cache)} tool geometry records from cache.")

    # Check for existing comparative results for resumption
    results_dict: Dict[str, Dict[str, Any]] = {}
    if os.path.exists(results_path):
        try:
            with open(results_path, "r", encoding="utf-8") as f:
                loaded = json.load(f)
                for item in loaded:
                    if "id" in item and item.get("no_tool", {}).get("status") == "OK":
                        results_dict[item["id"]] = item
            print(f"Found {len(results_dict)} existing completed records. Resuming...")
        except Exception:
            pass

    remaining_questions = [q for q in target_questions if q["id"] not in results_dict]
    print(f"Remaining questions to evaluate: {len(remaining_questions)} / {len(target_questions)}")

    llm_svc = get_llm_service()
    solver_agent = DeepMathSolverAgent()
    sem = asyncio.Semaphore(sem_limit)
    progress = {"completed": len(results_dict), "total": len(target_questions)}
    results_lock = asyncio.Lock()

    async def wrapped_eval(q):
        try:
            res = await evaluate_single_question(q, tool_geom_cache, llm_svc, solver_agent, sem, progress)
            async with results_lock:
                results_dict[q["id"]] = res
                with open(results_path, "w", encoding="utf-8") as f_out:
                    ordered = [results_dict.get(orig["id"]) for orig in target_questions if orig["id"] in results_dict]
                    json.dump(ordered, f_out, ensure_ascii=False, indent=2)
        except Exception as err:
            logger.error(f"Error evaluating {q.get('id')}: {err}", exc_info=True)
            print(f"Error evaluating {q.get('id')}: {err}")

    t_start = time.time()
    if remaining_questions:
        tasks = []
        for q in remaining_questions:
            t = asyncio.create_task(wrapped_eval(q))
            tasks.append(t)
            await asyncio.sleep(1.0)
        await asyncio.gather(*tasks)

    total_time = time.time() - t_start
    print(f"\n========================================================")
    print(f"EVALUATION COMPLETE IN {total_time:.1f}s ({total_time/60:.2f} min)")

    ordered_results = [results_dict[orig["id"]] for orig in target_questions if orig["id"] in results_dict]
    metrics = calculate_metrics(ordered_results)

    with open(metrics_path, "w", encoding="utf-8") as f:
        json.dump(metrics, f, ensure_ascii=False, indent=2)
    print(f"Saved comparative metrics to {metrics_path}")

    # Pretty print summary table
    s = metrics["summary"]
    nt_s = s["no_tool"]
    tl_s = s["tool"]

    print("\n================== TOOL vs NO-TOOL COMPARISON ==================")
    print(f"{'Metric':<35} | {'No-Tool (Pure LLM)':<20} | {'Tool (MathSolver)':<20}")
    print("-" * 80)
    print(f"{'Math Answer Accuracy':<35} | {nt_s['math_accuracy_count']}/{metrics['total_questions']} ({nt_s['math_accuracy_pct']}%) | {tl_s['math_accuracy_count']}/{metrics['total_questions']} ({tl_s['math_accuracy_pct']}%)")
    print(f"{'Strict Geometric Validity (0 err)':<35} | {nt_s['geometric_validity_count']}/{metrics['total_questions']} ({nt_s['geometric_validity_pct']}%) | {tl_s['geometric_validity_count']}/{metrics['total_questions']} ({tl_s['geometric_validity_pct']}%)")
    print(f"{'Avg Constraint Satisfaction (CSR)':<35} | {nt_s['avg_constraint_satisfaction_pct']:>18.1f}% | {tl_s['avg_constraint_satisfaction_pct']:>18.1f}%")
    print(f"{'Oxyz Coordinate MAE (lower better)':<35} | {str(nt_s['oxyz_coordinate_mae']):>19} | {str(tl_s['oxyz_coordinate_mae']):>19}")
    print(f"{'Median Latency':<35} | {nt_s['latency'].get('median_ms', 0):>17.1f}ms | {tl_s['latency'].get('median_ms', 0):>17.1f}ms")

    print("\n--- Breakdown By Topic ---")
    for top, stat in metrics["by_topic"].items():
        print(f"  * {top:<16} | Math Acc: NoTool={stat['notool_math_acc_pct']:5.1f}% vs Tool={stat['tool_math_acc_pct']:5.1f}% | CSR: NoTool={stat['notool_avg_csr_pct']:5.1f}% vs Tool={stat['tool_avg_csr_pct']:5.1f}% (N={stat['total']})")


if __name__ == "__main__":
    asyncio.run(main())