Devtu Benchmark Harness

{
  "skill_name": "devtu-benchmark-harness",
  "evals": [
    {
      "id": 1,
      "prompt": "Run the lab-bench benchmark with 5 questions in plugin-only mode and report accuracy",
      "expected_output": "Results saved with accuracy percentage",
      "files": []
    },
    {
      "id": 2,
      "prompt": "Grade the latest benchmark results and show per-category breakdown",
      "expected_output": "Category table with accuracy per subtask",
      "files": []
    },
    {
      "id": 3,
      "prompt": "Generate a benchmark report comparing plugin vs baseline on lab-bench",
      "expected_output": "Markdown report with comparison table and recommendations",
      "files": []
    }
  ]
}

#!/usr/bin/env python3
"""Failure analysis for a BixBench results JSON.

Classifies each wrong answer into a root-cause category and attributes
it to a skill (via the gt_skills.json question→skill mapping).

Outputs:
  - Per-failure record with category, evidence snippet, attributed skills
  - Per-skill heatmap (fail count, total used, fail rate) sorted by fail rate
  - Optional markdown report

Usage:
  python analyze_failures.py --results results.json [--out failure_analysis.json] \\
      [--report failure_report.md]
"""

import argparse
import json
import re
import sys
from collections import defaultdict
from pathlib import Path

REPO = Path(__file__).resolve().parents[3]
QUESTIONS_FILE = REPO / "skills" / "evals" / "bixbench" / "questions.json"
GT_SKILLS_FILE = (
    REPO / "temp_docs_and_tests" / "benchmark_tracking"
    / "skill_routing_test" / "gt_skills.json"
)


# Failure classification heuristics
TIMEOUT_PATTERNS = [
    r"timeout after",
    r"\bERROR:\s*$",  # bare "ERROR:" with empty body = max_turns hit
    r"max_turns",
]
TOOL_ERROR_PATTERNS = [
    r"MCP error",
    r"tool error",
    r"Traceback \(most recent call last\)",
    r"\bexception\b",
    r"command not found",
]
NO_DATA_PATTERNS = [
    r"can'?t find",
    r"no .* (data|files?) (in|found|located)",
    r"directory .* (doesn'?t exist|does not exist|empty)",
]


def classify_failure(predicted: str, ground_truth: str) -> tuple[str, str]:
    """Return (category, evidence_snippet)."""
    pred_lower = predicted.lower()

    # Compute incomplete (timeout, max_turns, empty)
    if not predicted.strip() or predicted.strip() == "ERROR:":
        return "compute_incomplete", "(empty result — likely max_turns)"
    for pat in TIMEOUT_PATTERNS:
        m = re.search(pat, predicted, re.IGNORECASE)
        if m:
            return "compute_incomplete", f"(matches /{pat}/)"

    # Tool error (crash, traceback)
    for pat in TOOL_ERROR_PATTERNS:
        m = re.search(pat, predicted, re.IGNORECASE)
        if m:
            return "tool_error", predicted[max(0, m.start()-20):m.end()+60]

    # Couldn't find data (signal that workspace path / file lookup broke)
    for pat in NO_DATA_PATTERNS:
        m = re.search(pat, predicted, re.IGNORECASE)
        if m:
            return "no_data_found", predicted[max(0, m.start()-20):m.end()+60]

    # Numeric-vs-numeric within an order of magnitude → semantic mismatch
    pred_nums = re.findall(r"-?\d+\.?\d*(?:[eE][+-]?\d+)?", predicted)
    gt_nums = re.findall(r"-?\d+\.?\d*(?:[eE][+-]?\d+)?", ground_truth)
    if pred_nums and gt_nums:
        try:
            p = float(pred_nums[-1])  # last number in prediction is usually the answer
            g = float(gt_nums[0])
            if g != 0 and abs((p - g) / g) < 100:  # within 100x
                return "semantic_mismatch", f"predicted={p}, gt={g}"
        except (ValueError, ZeroDivisionError):
            pass

    # Default: wrong answer with no specific signal
    return "wrong_answer", predicted[:120]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--results", required=True, help="results JSON from run_eval.py")
    ap.add_argument("--out", help="Write per-failure JSON")
    ap.add_argument("--report", help="Write markdown report")
    args = ap.parse_args()

    raw = json.loads(Path(args.results).read_text())
    runs = raw.get("with_plugin", raw if isinstance(raw, list) else [])
    if not runs:
        print(f"ERROR: no runs found in {args.results}", file=sys.stderr)
        sys.exit(1)

    questions = json.loads(QUESTIONS_FILE.read_text())
    qs_by_uuid = {q["id"]: q for q in questions}

    gt_skills = {}
    if GT_SKILLS_FILE.exists():
        gt_skills = json.loads(GT_SKILLS_FILE.read_text())

    failures = []
    skill_total = defaultdict(int)
    skill_fail = defaultdict(int)
    cat_count = defaultdict(int)

    for r in runs:
        uuid = r.get("id", "")
        q = qs_by_uuid.get(uuid, {})
        qid = q.get("question_id", "?")
        is_correct = r.get("correct", False)

        skills_for_q = gt_skills.get(qid, [])
        for s in skills_for_q:
            skill_total[s] += 1

        if is_correct:
            continue

        category, evidence = classify_failure(r.get("predicted", ""), str(r.get("ground_truth", "")))
        cat_count[category] += 1
        for s in skills_for_q:
            skill_fail[s] += 1

        failures.append({
            "qid": qid,
            "uuid": uuid,
            "category": category,
            "evidence": evidence[:200],
            "ground_truth": str(r.get("ground_truth", ""))[:80],
            "predicted_tail": r.get("predicted", "")[-160:],
            "elapsed_seconds": r.get("elapsed_seconds", 0),
            "skills": skills_for_q,
            "categories": q.get("categories", []),
        })

    total = len(runs)
    correct = sum(1 for r in runs if r.get("correct"))
    print(f"\nResults: {correct}/{total} correct ({100*correct/total:.1f}%)\n")

    print("Failure categories:")
    for cat, n in sorted(cat_count.items(), key=lambda x: -x[1]):
        print(f"  {cat:25s}  {n:3d}  ({100*n/total:.1f}%)")

    print("\nPer-skill heatmap (sorted by fail rate, min 3 uses):")
    print(f"  {'skill':40s} fails  total  rate")
    rows = []
    for s, used in skill_total.items():
        if used < 3:
            continue
        f = skill_fail.get(s, 0)
        rows.append((s, f, used, f / used if used else 0))
    rows.sort(key=lambda x: -x[3])
    for s, f, used, rate in rows:
        print(f"  {s:40s}  {f:3d}    {used:3d}   {100*rate:5.1f}%")

    if args.out:
        Path(args.out).write_text(json.dumps({
            "summary": {
                "total": total, "correct": correct,
                "score": correct / total,
            },
            "failure_categories": dict(cat_count),
            "skill_heatmap": [
                {"skill": s, "fails": f, "total": used, "rate": rate}
                for s, f, used, rate in rows
            ],
            "failures": failures,
        }, indent=2))
        print(f"\nWrote per-failure JSON to {args.out}")

    if args.report:
        lines = [
            "# Failure Analysis Report",
            "",
            f"**Score:** {correct}/{total} ({100*correct/total:.1f}%)",
            "",
            "## Failure Categories",
            "",
            "| Category | Count | % |",
            "|---|---|---|",
        ]
        for cat, n in sorted(cat_count.items(), key=lambda x: -x[1]):
            lines.append(f"| {cat} | {n} | {100*n/total:.1f}% |")
        lines += ["", "## Per-Skill Heatmap", "",
                  "| Skill | Fails | Total Used | Fail Rate |",
                  "|---|---|---|---|"]
        for s, f, used, rate in rows:
            lines.append(f"| {s} | {f} | {used} | {100*rate:.1f}% |")
        lines += ["", "## Failures (by category, then skill)", ""]
        for cat in sorted(cat_count, key=lambda c: -cat_count[c]):
            lines.append(f"### {cat} ({cat_count[cat]})")
            lines.append("")
            for f in failures:
                if f["category"] == cat:
                    lines.append(f"- **{f['qid']}** [{', '.join(f['skills'])}]")
                    lines.append(f"  - GT: `{f['ground_truth']}`")
                    lines.append(f"  - evidence: `{f['evidence']}`")
            lines.append("")
        Path(args.report).write_text("\n".join(lines))
        print(f"Wrote markdown report to {args.report}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Analyze benchmark results: category breakdown, failure classification,
skill mapping, and improvement recommendations.

Usage:
  python analyze_results.py --results graded.json
  python analyze_results.py --results graded.json --benchmark bixbench
  python analyze_results.py --results graded.json --diagnose
  python analyze_results.py --results graded.json --extract-failures failures.json
"""

import argparse
import json
from collections import defaultdict
from pathlib import Path


# Map question categories to the skills that should handle them
CATEGORY_TO_SKILL = {
    "DESeq2": "tooluniverse-rnaseq-deseq2",
    "DESeq2+enrichGO": "tooluniverse-gene-enrichment",
    "regression": "tooluniverse-statistical-modeling",
    "ANOVA": "tooluniverse-statistical-modeling",
    "chi_square": "tooluniverse-statistical-modeling",
    "spline_fitting": "tooluniverse-statistical-modeling",
    "pathway_enrichment": "tooluniverse-gene-enrichment",
    "phylogenetics": "tooluniverse-phylogenetics",
    "variant_analysis": "tooluniverse-variant-analysis",
    "fold_change": "tooluniverse-rnaseq-deseq2",
    "crispr_screen": "tooluniverse-crispr-screen-analysis",
    "single_cell": "tooluniverse-single-cell",
    "data_counting": "tooluniverse-statistical-modeling",
    "imaging": "tooluniverse-statistical-modeling",
    "statistical_test": "tooluniverse-statistical-modeling",
    "correlation": "tooluniverse-statistical-modeling",
    "transcriptomics": "tooluniverse-rnaseq-deseq2",
    "epigenomics": "tooluniverse-variant-analysis",
    "functional_genomics": "tooluniverse-crispr-screen-analysis",
    "data_integration": "tooluniverse-statistical-modeling",
    # lab-bench categories
    "ClinGen_task": "tooluniverse-gene-disease-association",
    "dga_task": "tooluniverse-gene-disease-association",
    "COSMIC_task": "tooluniverse-variant-analysis",
    "ClinVar_task": "tooluniverse-variant-analysis",
    "MeSH_task": "tooluniverse-disease-research",
    "MedGen_task": "tooluniverse-disease-research",
}

# Known patterns that cause failures, keyed by (category, failure_type)
# Known patterns: (category, failure_type) → (diagnosis, devtu_skill_to_invoke)
KNOWN_FAILURE_PATTERNS = {
    ("DESeq2", "wrong_answer"): (
        "pydeseq2 vs R DESeq2 disagreement; wrong strain mapping; wrong set operations. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-rnaseq-deseq2, or "
        "Skill('devtu-fix-tool') on run_deseq2_analysis if tool exists but returns wrong values."
    ),
    ("regression", "wrong_answer"): (
        "Wrong cohort selection (pre-filtering AEs); wrong denominator; model confusion. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-statistical-modeling."
    ),
    ("regression", "timeout"): (
        "Agent loops trying different model specs instead of committing. "
        "Action: Skill('devtu-optimize-skills') — add 'commit to first approach' rule."
    ),
    ("ANOVA", "wrong_answer"): (
        "F-statistic vs p-value confusion; wrong grouping; wrong aggregation level. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-statistical-modeling."
    ),
    ("spline_fitting", "wrong_answer"): (
        "Python patsy.cr() ≠ R ns(); endpoint inclusion depends on model type. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-statistical-modeling."
    ),
    ("phylogenetics", "wrong_answer"): (
        "PhyKIT batch computation errors; gap handling; taxa count. "
        "Action: Skill('devtu-fix-tool') on phykit_batch_analysis, or "
        "Skill('devtu-optimize-skills') on tooluniverse-phylogenetics."
    ),
    ("variant_analysis", "wrong_answer"): (
        "Multi-row Excel headers; coding vs all variant denominator. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-variant-analysis."
    ),
    ("fold_change", "wrong_answer"): (
        "Wrong contrast direction; gene symbol mapping error. "
        "Action: Skill('devtu-optimize-skills') on tooluniverse-rnaseq-deseq2."
    ),
    ("pathway_enrichment", "wrong_answer"): (
        "Wrong GO term after simplify; R clusterProfiler version sensitivity. "
        "Action: Skill('devtu-fix-tool') on run_deseq2_analysis enrichgo operation, or "
        "Skill('devtu-optimize-skills') on tooluniverse-gene-enrichment."
    ),
    ("pathway_enrichment", "timeout"): (
        "R enrichGO takes too long on large gene lists. "
        "Action: Skill('devtu-fix-tool') on run_deseq2_analysis — increase timeout or optimize R code."
    ),
    ("data_counting", "timeout"): (
        "Complex multi-step data pipeline times out. "
        "Action: Skill('devtu-create-tool') — bundle the computation as a ToolUniverse tool."
    ),
    ("data_counting", "wrong_answer"): (
        "Wrong filter, wrong denominator, or wrong aggregation. "
        "Action: Skill('devtu-optimize-skills') on the responsible skill."
    ),
}


def classify_failure(result: dict) -> str:
    """Classify a failure using devtu-fix-tool error taxonomy."""
    pred = result.get("predicted", "")
    if "Timeout" in pred:
        return "timeout"
    if pred.startswith("ERROR"):
        error_text = pred.lower()
        if "api key" in error_text:
            return "api_key_missing"
        if "404" in error_text or "not found" in error_text:
            return "endpoint_not_found"
        if "validation" in error_text or "schema" in error_text:
            return "schema_validation"
        return "tool_error"
    if not pred or pred.strip() == "":
        return "empty_response"
    return "wrong_answer"


def categorize_question(question: dict, benchmark: str = "") -> str:
    """Categorize a question for grouping.

    Uses the BixBench 'categories' field if available, falling back to
    keyword matching on the question text.
    """
    # Lab-bench: use subtask field
    subtask = question.get("subtask", "")
    if subtask:
        return subtask.replace("-v1-public", "")

    q = question.get("question", "").lower()

    # BixBench categories from the dataset itself (comma-separated)
    bix_cats = question.get("categories", "").lower()

    # Keyword-based categorization (ordered by specificity)
    if "spline" in q or "cubic model" in q or "natural spline" in q:
        return "spline_fitting"
    if "deseq2" in q and ("enrichgo" in q or "enrichment" in q):
        return "DESeq2+enrichGO"
    if "deseq2" in q or ("differential expression" in q and "rna" in bix_cats):
        return "DESeq2"
    if "ordinal" in q or "logistic regression" in q or "odds ratio" in q:
        return "regression"
    if "anova" in q or "f-statistic" in q:
        return "ANOVA"
    if "chi-square" in q or "chi square" in q or "chi²" in q:
        return "chi_square"
    if "pathway" in q and ("enrichment" in q or "gsea" in q or "kegg" in q):
        return "pathway_enrichment"
    if "phylogen" in q or "treeness" in q or "parsimony" in q:
        return "phylogenetics"
    if "variant" in q or "vcf" in q or "vaf" in q or "snp" in q:
        return "variant_analysis"
    if "fold change" in q or "log2fc" in q or "log2 fold" in q:
        return "fold_change"
    if "crispr" in q or "mageck" in q or "sgrna" in q:
        return "crispr_screen"
    if "scanpy" in q or "single-cell" in q or "umap" in q or "cluster" in q:
        return "single_cell"
    if "cohen" in q or "effect size" in q or "neun" in q:
        return "imaging"

    # Fall back to BixBench categories field
    if "phylogenetics" in bix_cats or "evolutionary" in bix_cats:
        return "phylogenetics"
    if "single-cell" in bix_cats or "single cell" in bix_cats:
        return "single_cell"
    if "differential expression" in bix_cats:
        return "DESeq2"
    if "rna-seq" in bix_cats or "transcriptomics" in bix_cats:
        return "transcriptomics"
    if "genomic variant" in bix_cats or "snp" in bix_cats:
        return "variant_analysis"
    if "epigenomics" in bix_cats:
        return "epigenomics"
    if "functional genomics" in bix_cats:
        return "functional_genomics"
    if "imaging" in bix_cats:
        return "imaging"
    if "machine learning" in bix_cats:
        return "data_integration"

    # Text-based fallbacks
    if "r-squared" in q or "regression" in q or "coefficient" in q:
        return "regression"
    if "how many" in q or "percentage" in q or "count" in q or "number of" in q:
        return "data_counting"
    if "mann-whitney" in q or "wilcoxon" in q or "t-test" in q or "p-value" in q:
        return "statistical_test"
    if "correlation" in q or "spearman" in q or "pearson" in q:
        return "correlation"
    if "enrichment" in q or "go " in q or "gene ontology" in q:
        return "pathway_enrichment"
    return "other"


def analyze(
    results_path: str, benchmark: str = "", questions_path: str = ""
) -> dict:
    """Analyze benchmark results.

    If questions_path is provided, enriches results with full question text
    and BixBench categories for better categorization.
    """
    with open(results_path) as f:
        data = json.load(f)

    # Load questions file for enrichment
    qmap = {}
    if questions_path:
        with open(questions_path) as f:
            for q in json.load(f):
                qmap[q["id"]] = q

    all_configs = {}
    if isinstance(data, dict):
        for config_name, results in data.items():
            if isinstance(results, list):
                all_configs[config_name] = results
    elif isinstance(data, list):
        all_configs["results"] = data

    analysis = {}
    for config_name, results in all_configs.items():
        correct = sum(1 for r in results if r.get("correct"))
        total = len(results)
        total_time = sum(r.get("elapsed_seconds", 0) for r in results)

        by_category = defaultdict(lambda: {"correct": 0, "total": 0, "failures": []})
        for r in results:
            # Enrich with full question text and categories if available
            enriched = r
            if r.get("id") in qmap:
                q = qmap[r["id"]]
                enriched = {
                    **r,
                    "question": q.get("question", r.get("question", "")),
                    "categories": q.get("categories", ""),
                }
            cat = categorize_question(enriched, benchmark)
            by_category[cat]["total"] += 1
            if r.get("correct"):
                by_category[cat]["correct"] += 1
            else:
                by_category[cat]["failures"].append(
                    {
                        "id": r.get("id", ""),
                        "type": classify_failure(r),
                        "ground_truth": str(r.get("ground_truth", ""))[:50],
                        "predicted": str(r.get("predicted", ""))[:100],
                        "skill": CATEGORY_TO_SKILL.get(cat, "unknown"),
                    }
                )

        failure_types = defaultdict(int)
        for r in results:
            if not r.get("correct"):
                failure_types[classify_failure(r)] += 1

        # Skill-level aggregation
        by_skill = defaultdict(lambda: {"correct": 0, "total": 0, "categories": set()})
        for cat, v in by_category.items():
            skill = CATEGORY_TO_SKILL.get(cat, "unknown")
            by_skill[skill]["correct"] += v["correct"]
            by_skill[skill]["total"] += v["total"]
            by_skill[skill]["categories"].add(cat)

        analysis[config_name] = {
            "overall": {
                "correct": correct,
                "total": total,
                "accuracy": round(correct / total * 100, 1) if total else 0,
                "total_time": round(total_time),
                "avg_time": round(total_time / total, 1) if total else 0,
            },
            "by_category": {
                cat: {
                    "correct": v["correct"],
                    "total": v["total"],
                    "accuracy": (
                        round(v["correct"] / v["total"] * 100) if v["total"] else 0
                    ),
                    "skill": CATEGORY_TO_SKILL.get(cat, "unknown"),
                    "failures": v["failures"],
                }
                for cat, v in sorted(by_category.items())
            },
            "by_skill": {
                skill: {
                    "correct": v["correct"],
                    "total": v["total"],
                    "accuracy": (
                        round(v["correct"] / v["total"] * 100) if v["total"] else 0
                    ),
                    "categories": sorted(v["categories"]),
                }
                for skill, v in sorted(
                    by_skill.items(), key=lambda x: x[1]["correct"] / max(x[1]["total"], 1)
                )
            },
            "failure_types": dict(failure_types),
        }

    return analysis


def diagnose(analysis: dict) -> list:
    """Generate targeted improvement recommendations from analysis."""
    recommendations = []
    for config_name, data in analysis.items():
        for cat, v in data["by_category"].items():
            if v["accuracy"] == 100:
                continue
            skill = CATEGORY_TO_SKILL.get(cat, "unknown")
            for f in v["failures"]:
                key = (cat, f["type"])
                pattern = KNOWN_FAILURE_PATTERNS.get(key)
                recommendations.append(
                    {
                        "priority": "HIGH" if v["accuracy"] < 50 else "MEDIUM",
                        "category": cat,
                        "skill": skill,
                        "failure_type": f["type"],
                        "question_id": f["id"][:12],
                        "ground_truth": f["ground_truth"],
                        "diagnosis": pattern or f"No known pattern for ({cat}, {f['type']}). Investigate manually.",
                    }
                )
    # Deduplicate by (category, failure_type)
    seen = set()
    unique = []
    for r in sorted(recommendations, key=lambda x: (0 if x["priority"] == "HIGH" else 1, x["category"])):
        key = (r["category"], r["failure_type"])
        if key not in seen:
            seen.add(key)
            unique.append(r)
    return unique


def extract_failures(results_path: str, questions_path: str = "") -> list:
    """Extract failed question IDs for retesting."""
    with open(results_path) as f:
        data = json.load(f)

    qmap = {}
    if questions_path:
        with open(questions_path) as f:
            for q in json.load(f):
                qmap[q["id"]] = q

    all_results = []
    if isinstance(data, dict):
        for results in data.values():
            if isinstance(results, list):
                all_results.extend(results)
    elif isinstance(data, list):
        all_results = data

    failures = []
    for r in all_results:
        if not r.get("correct"):
            enriched = r
            if r.get("id") in qmap:
                q = qmap[r["id"]]
                enriched = {**r, "question": q.get("question", "")}
            failures.append(
                {
                    "id": r["id"],
                    "question_id": r.get("question_id", enriched.get("question_id", "")),
                    "category": categorize_question(enriched),
                    "failure_type": classify_failure(r),
                    "ground_truth": r.get("ground_truth", ""),
                    "skill": CATEGORY_TO_SKILL.get(categorize_question(enriched), "unknown"),
                }
            )
    return failures


def print_analysis(analysis: dict, show_diagnose: bool = False):
    """Print analysis to stdout."""
    for config_name, data in analysis.items():
        overall = data["overall"]
        print(f"\n{'='*60}")
        print(
            f"{config_name}: {overall['correct']}/{overall['total']} "
            f"({overall['accuracy']}%)"
        )
        print(
            f"Total time: {overall['total_time']}s, "
            f"avg: {overall['avg_time']}s/question"
        )
        print(f"{'='*60}")

        print("\nBy skill:")
        for skill, v in data["by_skill"].items():
            cats = ", ".join(v["categories"])
            print(f"  {skill}: {v['correct']}/{v['total']} ({v['accuracy']}%) [{cats}]")

        print("\nBy category:")
        for cat, v in data["by_category"].items():
            marker = " **" if v["accuracy"] < 50 else ""
            print(f"  {cat}: {v['correct']}/{v['total']} ({v['accuracy']}%){marker}")

        if data["failure_types"]:
            print("\nFailure types:")
            for ftype, count in sorted(
                data["failure_types"].items(), key=lambda x: -x[1]
            ):
                print(f"  {ftype}: {count}")

        if show_diagnose:
            recs = diagnose({config_name: data})
            if recs:
                print("\n" + "=" * 60)
                print("IMPROVEMENT RECOMMENDATIONS")
                print("=" * 60)
                for r in recs:
                    print(
                        f"\n  [{r['priority']}] {r['category']} → {r['skill']}"
                    )
                    print(f"  Failure: {r['failure_type']}")
                    print(f"  {r['diagnosis']}")


def main():
    parser = argparse.ArgumentParser(description="Analyze benchmark results")
    parser.add_argument("--results", required=True, help="Path to results JSON")
    parser.add_argument(
        "--benchmark", default="", help="Benchmark name for categorization"
    )
    parser.add_argument(
        "--questions",
        default="",
        help="Path to questions JSON for full-text enrichment",
    )
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument(
        "--diagnose",
        action="store_true",
        help="Show improvement recommendations",
    )
    parser.add_argument(
        "--extract-failures",
        metavar="OUTPUT",
        help="Extract failed question IDs to a JSON file for retesting",
    )
    args = parser.parse_args()

    analysis = analyze(args.results, args.benchmark, args.questions)

    if args.extract_failures:
        failures = extract_failures(args.results, args.questions)
        with open(args.extract_failures, "w") as f:
            json.dump(failures, f, indent=2)
        print(f"Extracted {len(failures)} failures to {args.extract_failures}")
        return

    if args.json:
        output = analysis
        if args.diagnose:
            for config_name, data in analysis.items():
                output[config_name]["recommendations"] = diagnose(
                    {config_name: data}
                )
        print(json.dumps(output, indent=2))
    else:
        print_analysis(analysis, show_diagnose=args.diagnose)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Detect benchmark-specific content that has leaked into skill files.

The harness must drive GENERAL improvements, not dataset-specific answer
encoding. This script scans skills for signs of overfitting:
  - Benchmark names (BixBench, bix-N question IDs)
  - Dataset-specific capsule IDs / file names that only exist in one dataset
  - Numeric values suspiciously specific to a single question's ground truth

Run this before accepting a skill edit from devtu-optimize-skills.

Usage:
    python check_memorization.py skills/tooluniverse-*/SKILL.md
    python check_memorization.py --all   # check all tooluniverse-* skills
    python check_memorization.py --strict skills/path/to/skill/SKILL.md
"""

import argparse
import re
import sys
from pathlib import Path

REPO_ROOT = Path(__file__).resolve().parent.parent.parent.parent

# Patterns that suggest dataset-specific content. Organized by severity.
HARD_PATTERNS = [
    (re.compile(r"\bbix-\d+(?:-q\d+)?\b", re.IGNORECASE), "bix-N question ID"),
    (re.compile(r"\bBixBench\b", re.IGNORECASE), "benchmark name 'BixBench'"),
    (re.compile(r"\bBCG-CORONA\b", re.IGNORECASE), "dataset name 'BCG-CORONA'"),
    (re.compile(r"CapsuleFolder-[a-f0-9-]+"), "capsule UUID path"),
    (re.compile(r"TASK\d{3}_[A-Z-]+"), "TASK-prefixed dataset filename"),
    (re.compile(r"lab-bench[_-]?task"), "lab-bench task name"),
]

# Known benchmark-side GT issues (recorded here only to keep the harness honest
# about what it can/can't improve; NOT to encode answers into skills).
# When reporting scores, these can legitimately be excluded from the denominator.
KNOWN_GT_ISSUES_DOCS = """
Verified-mislabeled or unreproducible BixBench GTs (as of 2026-04-24):
- bix-53-q4: question is 'how many pathways' (integer), GT is '4.25E-04' (p-value)
- bix-31-q4: authoritative script `run_pc_vs_noncoding_ttest.py` gives p=0.481544, GT 0.65 unreproducible
- bix-38-q6: PhyKIT and BioPython independent computations give 1.904, GT 2.178 unreproducible
- bix-9-q3: question wording allows multiple valid answers
- bix-57-q1, bix-1-q1, bix-13-q1, bix-13-q3, bix-13-q4, bix-27-q5, bix-31-q3,
  bix-5-q1, bix-5-q4, bix-52-q3: verified discrepancies
"""


# Soft patterns — flag but don't fail by default
SOFT_PATTERNS = [
    # Specific filenames that appear only in one capsule
    (re.compile(r"BatchCorrectedReadCounts_Zenodo"), "specific dataset filename"),
    (re.compile(r"GeneMetaInfo_Zenodo"), "specific dataset filename"),
    (re.compile(r"Sample_annotated_Zenodo"), "specific dataset filename"),
    # Specific gene names referenced as *the answer* (not as examples)
    (re.compile(r"\bFAM138A\b"), "specific gene (possible GT-answer reference)"),
    (re.compile(r"\bJBX\d+\b"), "specific strain identifier"),
    # Exact numeric ground truths from known questions
    (re.compile(r"\bF\s*[=≈]\s*0\.7[67]\b"), "specific F-statistic (bix-36-q1 GT)"),
    (re.compile(r"\b4,?550\b"), "specific count (bix-14-q3 GT)"),
    (re.compile(r"\b2\.178\b"), "specific median (bix-38-q6 GT)"),
]


def scan_file(path: Path, strict: bool = False):
    """Scan a single file for memorization signals. Returns list of hits."""
    if not path.exists():
        return [("error", f"file not found: {path}")]

    # Skip test files — they are allowed to reference benchmarks
    if "test_" in path.name or path.name.startswith("test_"):
        return []

    text = path.read_text(errors="replace")
    hits = []

    for rx, label in HARD_PATTERNS:
        for m in rx.finditer(text):
            line_no = text[: m.start()].count("\n") + 1
            ctx = text.splitlines()[line_no - 1][:120] if line_no <= len(text.splitlines()) else ""
            hits.append(("hard", f"{path.relative_to(REPO_ROOT)}:{line_no}: {label}: {m.group()}  | {ctx}"))

    if strict:
        for rx, label in SOFT_PATTERNS:
            for m in rx.finditer(text):
                line_no = text[: m.start()].count("\n") + 1
                ctx = text.splitlines()[line_no - 1][:120] if line_no <= len(text.splitlines()) else ""
                hits.append(("soft", f"{path.relative_to(REPO_ROOT)}:{line_no}: {label}: {m.group()}  | {ctx}"))

    return hits


def main():
    parser = argparse.ArgumentParser(description="Detect memorization in skill files")
    parser.add_argument("paths", nargs="*", help="Paths to SKILL.md or similar files")
    parser.add_argument("--all", action="store_true", help="Scan all tooluniverse-*/ skill files")
    parser.add_argument("--strict", action="store_true", help="Include soft patterns (specific gene names, numeric GTs)")
    args = parser.parse_args()

    paths = []
    if args.all:
        for skill_dir in (REPO_ROOT / "skills").glob("tooluniverse-*"):
            for p in skill_dir.rglob("*.md"):
                if "references/" in str(p) or p.name in ("SKILL.md", "QUICK_START.md"):
                    paths.append(p)
            for p in skill_dir.glob("scripts/*.py"):
                paths.append(p)
    else:
        paths = [Path(p) for p in args.paths]

    if not paths:
        parser.print_help()
        sys.exit(1)

    all_hard = []
    all_soft = []
    for p in paths:
        for severity, msg in scan_file(p, strict=args.strict):
            if severity == "hard":
                all_hard.append(msg)
            elif severity == "soft":
                all_soft.append(msg)
            else:
                print(f"[error] {msg}", file=sys.stderr)

    if all_hard:
        print("=" * 60)
        print(f"HARD FAILURES ({len(all_hard)}):")
        print("=" * 60)
        for h in all_hard:
            print(f"  {h}")

    if all_soft and args.strict:
        print()
        print("=" * 60)
        print(f"SOFT WARNINGS ({len(all_soft)}):")
        print("=" * 60)
        for s in all_soft:
            print(f"  {s}")

    print()
    print(f"Scanned {len(paths)} files")
    print(f"Hard failures: {len(all_hard)}")
    if args.strict:
        print(f"Soft warnings: {len(all_soft)}")

    # Exit non-zero if hard failures found
    if all_hard:
        sys.exit(1)
    sys.exit(0)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Compute SHA256 of every file in bixbench_clean/data/CapsuleFolder-*.

Writes the result to bixbench_clean/checksums.json. This file is the
canonical fingerprint for the clean capsules — the harness verifies it
before each run to detect accidental modification.

Usage:
    python compute_capsule_checksums.py
"""

import hashlib
import json
import sys
from pathlib import Path

REPO = Path(__file__).resolve().parents[3]
CLEAN = REPO / "temp_docs_and_tests" / "bixbench_clean" / "data"
OUT = REPO / "temp_docs_and_tests" / "bixbench_clean" / "checksums.json"


def file_sha256(p: Path) -> str:
    h = hashlib.sha256()
    with p.open("rb") as f:
        for block in iter(lambda: f.read(1 << 20), b""):
            h.update(block)
    return h.hexdigest()


def main():
    if not CLEAN.exists():
        print(f"ERROR: clean dir not found: {CLEAN}", file=sys.stderr)
        sys.exit(1)

    capsules = sorted(CLEAN.glob("CapsuleFolder-*"))
    print(f"Hashing {len(capsules)} capsules...")
    result = {}
    for cap in capsules:
        files = {}
        for f in sorted(cap.rglob("*")):
            if f.is_file():
                rel = str(f.relative_to(cap))
                files[rel] = file_sha256(f)
        result[cap.name] = files
        print(f"  {cap.name}: {len(files)} files", flush=True)

    OUT.write_text(json.dumps(result, indent=2, sort_keys=True))
    n_files = sum(len(v) for v in result.values())
    print(f"\nWrote {n_files} file hashes across {len(result)} capsules to {OUT}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Generate structured markdown benchmark report.

Usage:
  python generate_report.py --results graded.json --output BENCHMARK_REPORT.md
  python generate_report.py --results graded.json --benchmark lab-bench
"""

import argparse
import json
import time
from pathlib import Path

from analyze_results import analyze, classify_failure, categorize_question


def generate_report(results_path: str, benchmark: str = "", output_path: str = None) -> str:
    """Generate a markdown benchmark report."""
    analysis = analyze(results_path, benchmark)
    lines = []

    lines.append(f"# ToolUniverse Benchmark Report")
    lines.append(f"Generated: {time.strftime('%Y-%m-%d %H:%M')}")
    lines.append("")

    for config_name, data in analysis.items():
        overall = data["overall"]
        lines.append(f"## {config_name}: {overall['correct']}/{overall['total']} ({overall['accuracy']}%)")
        lines.append(f"Total time: {overall['total_time']}s | Avg: {overall['avg_time']}s/question")
        lines.append("")

        # Category table
        lines.append("### By Category")
        lines.append("")
        lines.append("| Category | Correct | Total | Accuracy |")
        lines.append("|----------|---------|-------|----------|")
        for cat, v in data["by_category"].items():
            lines.append(f"| {cat} | {v['correct']} | {v['total']} | {v['accuracy']}% |")
        lines.append("")

        # Failure breakdown
        if data["failure_types"]:
            lines.append("### Failure Types")
            lines.append("")
            lines.append("| Type | Count |")
            lines.append("|------|-------|")
            for ftype, count in sorted(data["failure_types"].items(), key=lambda x: -x[1]):
                lines.append(f"| {ftype} | {count} |")
            lines.append("")

        # Top failures detail
        all_failures = []
        for cat, v in data["by_category"].items():
            for f in v["failures"]:
                all_failures.append({**f, "category": cat})
        if all_failures:
            lines.append("### Top Failures")
            lines.append("")
            for f in all_failures[:15]:
                lines.append(f"- **{f['category']}** [{f['type']}]: GT=`{f['ground_truth']}`")
            lines.append("")

    # Recommendations
    lines.append("## Recommendations")
    lines.append("")
    for config_name, data in analysis.items():
        for ftype, count in data["failure_types"].items():
            if ftype == "timeout" and count > 2:
                lines.append(f"- Increase timeout or optimize scripts ({count} timeouts)")
            if ftype == "api_key_missing" and count > 0:
                lines.append(f"- Set missing API keys ({count} failures)")
            if ftype == "wrong_answer" and count > 2:
                lines.append(f"- Investigate wrong answers ({count} incorrect computations)")
            if ftype == "tool_error" and count > 0:
                lines.append(f"- Fix tool errors ({count} failures)")

    report = "\n".join(lines)

    if output_path:
        with open(output_path, "w") as f:
            f.write(report)
        print(f"Report saved to {output_path}")

    return report


def main():
    parser = argparse.ArgumentParser(description="Generate benchmark report")
    parser.add_argument("--results", required=True, help="Path to results JSON")
    parser.add_argument("--benchmark", default="", help="Benchmark name")
    parser.add_argument("--output", help="Output markdown path")
    args = parser.parse_args()

    report = generate_report(args.results, args.benchmark, args.output)
    if not args.output:
        print(report)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Multi-strategy answer grading for ToolUniverse benchmarks.

Strategies (applied in order):
  1. Exact match — ground truth substring in prediction
  2. MC match — letter answer detection (A/B/C/D)
  3. Range match — numeric value within (low, high), with rounding tolerance
  4. Normalized match — strip punctuation, bidirectional substring
  5. Numeric proximity — within configurable tolerance (default 5%)
  6. Synonym match — common scientific term equivalences
  7. LLM verifier — Claude judges correctness (for complex text answers)

Usage:
  python grade_answers.py --results results.json --output graded.json
  python grade_answers.py --results results.json --llm
  python grade_answers.py --results results.json --numeric-tolerance 0.10
"""

import argparse
import json
import re
import subprocess
import sys
from pathlib import Path


def _normalize_unicode_sci(text: str) -> str:
    """Convert Unicode scientific notation to ASCII.

    Handles patterns like '7.04 × 10⁻²⁶' → '7.04e-26'
    and '9.13 × 10⁻¹³' → '9.13e-13'.
    """
    import re as _re

    superscript_map = {
        "\u2070": "0", "\u00b9": "1", "\u00b2": "2", "\u00b3": "3",
        "\u2074": "4", "\u2075": "5", "\u2076": "6", "\u2077": "7",
        "\u2078": "8", "\u2079": "9", "\u207b": "-", "\u207a": "+",
    }

    def _replace_sci(m):
        base = m.group(1)
        sup_chars = m.group(2)
        exp = "".join(superscript_map.get(c, c) for c in sup_chars)
        return f"{base}e{exp}"

    # Match: number × 10 followed by superscript digits
    sup_pattern = (
        r"(-?(?:\d+\.?\d*))\s*[×x]\s*10(["
        + "".join(re.escape(c) for c in superscript_map)
        + r"]+)"
    )
    text = _re.sub(sup_pattern, _replace_sci, text)
    return text


def _extract_numbers(text: str) -> list[float]:
    """Extract all numbers from text, handling comma-grouped thousands and exponents.

    Also normalizes Unicode minus signs (U+2212), en-dashes (U+2013),
    and Unicode superscript scientific notation (e.g., 10⁻²⁶ → e-26).
    """
    # Normalize Unicode minus variants to ASCII hyphen
    text = text.replace("\u2212", "-").replace("\u2013", "-").replace("\u2014", "-")
    # Normalize Unicode scientific notation
    text = _normalize_unicode_sci(text)
    pattern = r"-?(?:\d{1,3}(?:,\d{3})+|\d+)(?:\.\d+)?(?:[eE][+-]?\d+)?"
    out = []
    for m in re.findall(pattern, text):
        try:
            out.append(float(m.replace(",", "")))
        except ValueError:
            pass
    return out


# Common scientific term synonyms
_SYNONYMS = {
    "normal": ["gaussian", "bell-shaped", "bell shaped", "bell curve"],
    "gaussian": ["normal", "bell-shaped", "bell shaped", "bell curve"],
    "upregulated": ["overexpressed", "up-regulated", "increased expression"],
    "downregulated": ["underexpressed", "down-regulated", "decreased expression"],
    "downregulation": ["down-regulation", "decreased expression"],
    "upregulation": ["up-regulation", "increased expression"],
}


def grade_answer(
    predicted: str,
    ground_truth: str,
    eval_mode: str = "",
    numeric_tolerance: float = 0.05,
    use_llm: bool = True,
    question: str = "",
) -> dict:
    """Grade a predicted answer against ground truth using multiple strategies.

    For eval_mode="llm_verifier", LLM grading is used by default when other
    strategies fail. Set use_llm=False to skip LLM grading.
    """
    ground_truth = str(ground_truth)
    gt_lower = ground_truth.strip().lower()
    # Normalize Unicode minus/dash variants and scientific notation in prediction
    predicted_normalized = (
        predicted.replace("\u2212", "-").replace("\u2013", "-").replace("\u2014", "-")
    )
    predicted_normalized = _normalize_unicode_sci(predicted_normalized)
    pred_lower = predicted_normalized.lower()

    # --- Strategy 1: Exact match ---
    # For short numeric GTs ("3", "29", "100"), naive substring match yields
    # false positives ("3" matches inside "0.523"). Use word-boundary regex so
    # the GT must appear as a token, not embedded inside a larger number.
    if re.fullmatch(r"-?\d+(?:\.\d+)?", gt_lower):
        exact_match = bool(re.search(rf"(?<![\d.]){re.escape(gt_lower)}(?![\d.])", pred_lower))
    else:
        exact_match = gt_lower in pred_lower

    # --- Strategy 2: MC match (short answers like A, B, Yes, No) ---
    mc_match = False
    if len(gt_lower) <= 3:
        patterns = [
            rf"\b{re.escape(gt_lower)}\b",
            rf"answer[:\s]+{re.escape(gt_lower)}",
            rf"\({re.escape(gt_lower)}\)",
        ]
        mc_match = any(re.search(p, pred_lower) for p in patterns)

    # --- Strategy 3: Range match ---
    range_match = False
    if eval_mode == "range_verifier" or (
        gt_lower.startswith("(") and "," in gt_lower
    ):
        try:
            low, high = gt_lower.strip("()").split(",")
            low, high = float(low), float(high)
            if low > high:
                low, high = high, low

            candidates: list[float] = []
            # Bare numbers (with optional thousand-separators / scientific exponent).
            for raw in re.findall(
                r"-?(?:\d{1,3}(?:,\d{3})+|\d+)(?:\.\d+)?(?:[eE][+-]?\d+)?",
                predicted_normalized,
            ):
                clean = raw.replace(",", "")
                try:
                    candidates.append(float(clean))
                except ValueError:
                    continue
            # `A:B` ratio expressions (e.g. "10:1") — convert to fraction
            # A/(A+B). Keeps range-matching robust to ratio-vs-fraction format
            # differences when the question asks for "frequency ratio of A:B"
            # but GT is given as the fraction.
            for a_str, b_str in re.findall(r"\b(\d+(?:\.\d+)?):(\d+(?:\.\d+)?)\b", predicted_normalized):
                try:
                    a, b = float(a_str), float(b_str)
                    if a + b > 0:
                        candidates.append(a / (a + b))
                        candidates.append(b / (a + b))
                except ValueError:
                    continue

            for val in candidates:
                if low <= val <= high:
                    range_match = True
                    break
                # Rounding tolerance: if predicted was rounded to N decimals,
                # accept if ±0.5 ULP at that precision overlaps [low, high].
                # Only applies when the predicted value has explicit decimals —
                # bare integers from prose ("1. Filter", "n=565") should NOT
                # match small decimal ranges like (0.76, 0.78).
                s = repr(val)
                if "." in s and "e" not in s.lower():
                    dp = len(s.split(".")[1])
                    half_ulp = 0.5 * (10 ** (-dp))
                    if val - half_ulp <= high and val + half_ulp >= low:
                        range_match = True
                        break
        except (ValueError, IndexError):
            pass

    # --- Strategy 4: Normalized match (bidirectional) ---
    gt_norm = re.sub(r"[^a-z0-9]", "", gt_lower)
    pred_norm = re.sub(r"[^a-z0-9]", "", pred_lower)
    normalized_match = (
        (len(gt_norm) >= 4 and gt_norm in pred_norm)
        or (len(pred_norm) >= 4 and pred_norm in gt_norm)
    )
    # Also check bold/quoted segments from the prediction individually
    # (e.g., "**CD14 Mono**" extracted as "CD14 Mono")
    if not normalized_match:
        for seg in re.findall(r"\*\*([^*]+)\*\*|\"([^\"]+)\"|'([^']+)'", predicted_normalized):
            segment = next(s for s in seg if s)
            seg_norm = re.sub(r"[^a-z0-9]", "", segment.lower())
            if len(seg_norm) >= 4 and (
                seg_norm in gt_norm or gt_norm in seg_norm
            ):
                normalized_match = True
                break

    # --- Strategy 5: Numeric proximity ---
    numeric_match = False
    # Try to extract the leading number from GT
    gt_clean = gt_lower.replace("%", "").replace(",", "").strip()
    gt_num_match = re.match(
        r"^\s*(-?(?:\d{1,3}(?:,\d{3})+|\d+)(?:\.\d+)?(?:[eE][+-]?\d+)?)",
        gt_clean,
    )
    if gt_num_match:
        try:
            gt_num = float(gt_num_match.group(1).replace(",", ""))
            for pred_num in _extract_numbers(predicted_normalized):
                if gt_num == 0:
                    # Accept exact 0 or any value so small that it represents
                    # a scipy p-value underflow (< 1e-10).
                    if abs(pred_num) < 1e-10:
                        numeric_match = True
                        break
                elif abs(pred_num - gt_num) / abs(gt_num) < numeric_tolerance:
                    numeric_match = True
                    break
        except ValueError:
            pass

    # --- Strategy 6: Synonym match ---
    synonym_match = False
    for aliases in [_SYNONYMS.get(gt_lower.strip(), [])]:
        if any(a in pred_lower for a in aliases):
            synonym_match = True
            break

    # --- Strategy 7: LLM verifier ---
    llm_match = False
    deterministic_correct = any([
        exact_match, mc_match, range_match, normalized_match,
        numeric_match, synonym_match,
    ])
    if (
        use_llm
        and eval_mode == "llm_verifier"
        and not deterministic_correct
    ):
        llm_match = _llm_grade(predicted, ground_truth, question)

    correct = deterministic_correct or llm_match

    return {
        "correct": correct,
        "strategies": {
            "exact_match": exact_match,
            "mc_match": mc_match,
            "range_match": range_match,
            "normalized_match": normalized_match,
            "numeric_match": numeric_match,
            "synonym_match": synonym_match,
            "llm_match": llm_match,
        },
        "ground_truth": ground_truth,
        "predicted_excerpt": predicted[:500],
    }


def _llm_grade(predicted: str, ground_truth: str, question: str) -> bool:
    """Use Claude to grade a complex answer.

    The LLM is asked to compare the semantic meaning, not exact wording.
    This handles cases like:
    - "35%" matching "33-36% increase"
    - "OR≈1.02, not significant" matching "No significant effect (OR≈1.02)"
    - Pathway names with slight wording differences
    """
    prompt = f"""Grade whether the predicted answer is correct by comparing it to the expected answer.

Question: {question[:300]}

Expected answer: {ground_truth}

Predicted answer: {predicted[:800]}

Grading rules:
1. Focus on whether the KEY CLAIM matches, not exact wording
2. If expected says "no significant effect" and predicted says "not significant" or "p > 0.05" or "OR ≈ 1.0", that's CORRECT
3. If expected gives a range like "33-36% increase" and predicted gives a number within that range (e.g., "35%"), that's CORRECT
4. If expected gives a numeric value and predicted gives a value within 5% tolerance, that's CORRECT
5. Minor wording differences (e.g., "CD14 Mono" vs "CD14 Monocytes") are CORRECT
6. If the predicted answer contains the right information buried in a longer explanation, that's CORRECT
7. If the predicted answer is fundamentally different in value or conclusion, that's WRONG

Reply with exactly one word: CORRECT or WRONG"""

    try:
        r = subprocess.run(
            [
                "claude", "-p", prompt,
                "--max-turns", "1",
                "--output-format", "json",
                "--allowedTools", "",
            ],
            capture_output=True,
            text=True,
            timeout=60,
        )
        if r.returncode == 0:
            data = json.loads(r.stdout)
            result = data.get("result", "")
            if isinstance(result, list):
                text = " ".join(
                    b.get("text", "")
                    for b in result
                    if isinstance(b, dict) and b.get("type") == "text"
                )
            else:
                text = str(result)
            text = text.strip().lower()
            return text.startswith("correct") or "correct" in text[:20]
    except Exception:
        pass
    return False


def grade_results_file(
    results_path: str,
    output_path: str = None,
    numeric_tolerance: float = 0.05,
    use_llm: bool = False,
) -> dict:
    """Grade all answers in a results file."""
    with open(results_path) as f:
        data = json.load(f)

    all_results = {}
    if isinstance(data, dict):
        for config_name, results in data.items():
            if isinstance(results, list):
                all_results[config_name] = results
    elif isinstance(data, list):
        all_results["results"] = data

    graded = {}
    for config_name, results in all_results.items():
        graded_list = []
        for r in results:
            gt = str(r.get("ground_truth", r.get("ideal", r.get("answer", ""))))
            pred = r.get("predicted", "")
            eval_mode = r.get("eval_mode", "")
            question = r.get("question", "")

            grade = grade_answer(
                pred, gt, eval_mode, numeric_tolerance, use_llm, question
            )
            graded_list.append({**r, **grade})
        graded[config_name] = graded_list

    if output_path:
        with open(output_path, "w") as f:
            json.dump(graded, f, indent=2)

    summary = {}
    for config_name, results in graded.items():
        correct = sum(1 for r in results if r["correct"])
        total = len(results)
        summary[config_name] = {
            "correct": correct,
            "total": total,
            "accuracy": round(correct / total * 100, 1) if total else 0,
        }

    return {"graded": graded, "summary": summary}


def main():
    parser = argparse.ArgumentParser(description="Grade benchmark answers")
    parser.add_argument("--results", required=True, help="Path to results JSON")
    parser.add_argument("--output", help="Path to save graded results")
    parser.add_argument(
        "--numeric-tolerance",
        type=float,
        default=0.05,
        help="Numeric tolerance (default 0.05)",
    )
    parser.add_argument(
        "--llm", action="store_true", help="Use LLM for complex grading"
    )
    args = parser.parse_args()

    result = grade_results_file(
        args.results, args.output, args.numeric_tolerance, args.llm
    )
    print(json.dumps(result["summary"], indent=2))


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Unified benchmark runner for ToolUniverse plugin evaluation.

Runs questions through Claude Code with and without the plugin,
captures outputs, grades answers, and saves results.

Usage:
  python run_eval.py --benchmark lab-bench --mode comparison --n 20
  python run_eval.py --benchmark bixbench --mode plugin-only --n 10 --category DESeq2
  python run_eval.py --data-file custom.json --mode baseline-only
"""

import argparse
import contextlib
import hashlib
import json
import os
import shutil
import subprocess
import sys
import tempfile
import time
from pathlib import Path

REPO_ROOT = Path(__file__).resolve().parents[3]
PLUGIN_DIR = str(REPO_ROOT / "dist" / "tooluniverse-plugin")
EVALS_DIR = REPO_ROOT / "skills" / "evals"
CLEAN_DATA_DIR = REPO_ROOT / "temp_docs_and_tests" / "bixbench_clean" / "data"
CHECKSUMS_FILE = REPO_ROOT / "temp_docs_and_tests" / "bixbench_clean" / "checksums.json"
BIXBENCH_DATA_DIRS = [
    CLEAN_DATA_DIR,
    REPO_ROOT / "temp_docs_and_tests" / "bixbench" / "data",
    REPO_ROOT / "temp_docs_and_tests" / "bixbench" / "bixbench" / "data",
]

# Import grading from sibling script
sys.path.insert(0, str(Path(__file__).parent))
from grade_answers import grade_answer


def _file_sha256(p: Path) -> str:
    h = hashlib.sha256()
    with p.open("rb") as f:
        for block in iter(lambda: f.read(1 << 20), b""):
            h.update(block)
    return h.hexdigest()


def verify_capsule_checksums(capsule: Path, expected: dict) -> list[str]:
    """Return a list of mismatch messages (empty = capsule is clean)."""
    mismatches = []
    actual_files = {
        str(f.relative_to(capsule)): f
        for f in capsule.rglob("*") if f.is_file()
    }
    # Files present in canonical but missing in capsule
    for rel, exp_hash in expected.items():
        if rel not in actual_files:
            mismatches.append(f"missing: {rel}")
            continue
        got = _file_sha256(actual_files[rel])
        if got != exp_hash:
            mismatches.append(f"hash mismatch: {rel}")
    # Files in capsule but not in canonical (= contamination)
    for rel in actual_files:
        if rel not in expected:
            mismatches.append(f"unexpected: {rel}")
    return mismatches


@contextlib.contextmanager
def isolated_capsule(capsule: Path):
    """Yield a fresh writable copy of the capsule in a tmpdir.

    The agent operates only on the copy; the canonical capsule stays
    untouched. The tmpdir is auto-deleted after the question.
    """
    with tempfile.TemporaryDirectory(prefix=f"{capsule.name}_") as tmp:
        workspace = Path(tmp) / capsule.name
        # Canonical capsule is read-only (a-w); copytree preserves perms,
        # so we restore write permissions on the copy after copying.
        shutil.copytree(capsule, workspace, symlinks=False)
        workspace.chmod(0o755)
        for p in workspace.rglob("*"):
            if p.is_file():
                p.chmod(0o644)
            elif p.is_dir():
                p.chmod(0o755)
        yield workspace


def load_guidance(guidance_path: str = None) -> str:
    """Load guidance text from a file, stripping YAML frontmatter."""
    if guidance_path is None:
        guidance_path = str(Path(PLUGIN_DIR) / "commands" / "research.md")
    path = Path(guidance_path)
    if not path.exists():
        return ""
    text = path.read_text()
    if text.startswith("---"):
        parts = text.split("---", 2)
        if len(parts) >= 3:
            text = parts[2]
    return text.strip()


# Keyword → skill routing (mirrors the router skill's table). When
# full_skill_injection is enabled, we use this to pick the matching
# sub-skill and pre-load its full SKILL.md body into the system prompt
# via --append-system-prompt. This guarantees the conventions reach
# inference even in `-p` mode where plugin auto-routing is unreliable.
def categorize_for_skill(question_text: str) -> str | None:
    """Route question to a sub-skill. Order: most-specific first.

    Domain-keywords (cpg/colony/treeness/...) are checked BEFORE generic
    statistical method keywords (anova/chi-square/...) so domain-specific
    skills win when both could match.
    """
    q = question_text.lower()
    # Domain-specific (highest priority)
    if "phylo" in q or "treeness" in q or "parsimony" in q or "phykit" in q or "saturation" in q or "dvmc" in q or "tree length" in q or "long branch" in q or "ortholog" in q or ("alignment" in q and ("gap" in q or "mafft" in q)):
        return "tooluniverse-phylogenetics"
    if "methylation" in q or "cpg " in q or " cpg" in q or "5mc" in q or "chip-seq" in q or "atac" in q or "m6a" in q or "chromatin" in q:
        return "tooluniverse-epigenomics"
    if "colony" in q or "circularity" in q or "swarming" in q or "cell area" in q or "morphometry" in q or "fluorescence" in q or "imagej" in q or "cellprofiler" in q or ("mean" in q and "area" in q):
        return "tooluniverse-image-analysis"
    if "variant" in q or "vcf" in q or "vaf" in q or " snp " in q or "haplotypecaller" in q or "indel" in q or "mutation" in q:
        return "tooluniverse-variant-analysis"
    if "crispr" in q or "mageck" in q or "sgrna" in q:
        return "tooluniverse-crispr-screen-analysis"
    if "scanpy" in q or "single-cell" in q or "h5ad" in q:
        return "tooluniverse-single-cell"
    if "fastq" in q or "bwa" in q or "samtools" in q or "trimmomatic" in q or "alignment quality" in q:
        return "tooluniverse-sequence-analysis"
    if "mass spec" in q or "tmt" in q or "proteomics" in q:
        return "tooluniverse-proteomics-analysis"
    # Pipeline-specific
    if "deseq2" in q or "differential expression" in q or "differentially expressed" in q or "fold change" in q or " log2" in q or "deg " in q:
        return "tooluniverse-rnaseq-deseq2"
    if "enrichgo" in q or "enrichment" in q or " go " in q or "kegg" in q or "gseapy" in q or "reactome" in q or "wikipathways" in q:
        return "tooluniverse-gene-enrichment"
    # Generic statistics (lowest priority — domain-specific wins above)
    if "anova" in q or "regression" in q or "chi-square" in q or "spline" in q or "cohen" in q or "f-statistic" in q or "odds ratio" in q:
        return "tooluniverse-statistical-modeling"
    return None


def precompute_for_capsule(capsule_path: Path, question_text: str) -> str:
    """Detect data patterns in the capsule and pre-run matching deterministic
    scripts. Returns a markdown block to inject into the prompt that contains
    the script invocation + its output. Returns "" if no pattern matches.

    Per skill-creator wisdom: when a script reproduces the GT deterministically,
    the agent should see its output BEFORE writing its own analysis. Otherwise
    the agent reinvents the wheel and may pick the wrong interpretation.

    Patterns supported:
      1. target_orthologs.txt + *.busco.zip → busco_target_orthologs.py
      2. long-format CpG CSV (Pos, Chromosome, MethylationPercentage) +
         chromosome length CSV → methylation_density.py
      3. HaplotypeCaller: BAM + reference FASTA → gatk_haplotypecaller_pipeline.py
      4. variant questions with VCF + question keywords → variant skill computations
      5. RNA-seq Pearson r vs length: counts CSV + metadata → custom analysis
      6. scogs zip phylogenetics: scogs_*.zip → scogs_paired_compare.py
         (--only-with-trees gated by _TREE_METRICS so alignment-only metrics work)
      7. Direct treefile capsules: *.treefile (no scogs zip) → phykit treeness
         with per-tree values and per-N averages
      8. Swarm CSV (Ratio + Area + StrainNumber) → spline_model_compare.py
         with canonical R notebook filter (drops strains "1" and "98")
      9. SDTM clinical trial: AE+DM+MH CSVs + treatment/OR question
         → sdtm_ordinal_logistic.py (3-way merge + OrderedModel ordinal logistic)
    """
    if capsule_path is None or not capsule_path.exists():
        return ""

    files = {p.name: p for p in capsule_path.iterdir() if p.is_file()}
    files_lower = {k.lower(): v for k, v in files.items()}
    blocks = []
    repo_root = Path(__file__).resolve().parents[3]

    # Pattern 1: BUSCO target_orthologs intersection
    has_target_list = "target_orthologs.txt" in files
    has_busco_zips = any(name.endswith(".busco.zip") for name in files)
    if has_target_list and has_busco_zips:
        script = repo_root / "skills" / "tooluniverse-phylogenetics" / "scripts" / "busco_target_orthologs.py"
        if script.exists():
            try:
                r = subprocess.run(
                    ["python3", str(script), "--capsule", str(capsule_path)],
                    capture_output=True, text=True, timeout=120,
                )
                output = (r.stdout + "\n" + r.stderr).strip()
                blocks.append(
                    "## Pre-computed analysis (BUSCO target_orthologs)\n\n"
                    f"```\n$ python3 {script.relative_to(repo_root)} --capsule <capsule>\n{output}\n```\n\n"
                    "The script computed every common interpretation. Pick the SUMMARY line "
                    "matching the analysis context (per-group is the typical convention)."
                )
            except (subprocess.TimeoutExpired, FileNotFoundError):
                pass

    # Pattern 2: Methylation density (CpG long-format + chromosome length)
    cpg_csvs = []
    chr_csvs = []
    for name, path in files.items():
        nl = name.lower()
        if not nl.endswith(".csv"):
            continue
        if "chromosome_length" in nl or "chr_length" in nl:
            chr_csvs.append(path)
        elif "cpg" in nl and "methylat" in nl.replace("_", "").replace("-", "").lower():
            cpg_csvs.append(path)
        elif nl.endswith(".csv"):
            # Heuristic: peek at the header line for the right columns
            try:
                with open(path) as fh:
                    header = fh.readline().lower()
                if "methylationpercentage" in header.replace("_", "") and "pos" in header and "chromosome" in header:
                    cpg_csvs.append(path)
                elif "chromosome" in header and "length" in header and len(header.split(",")) <= 4:
                    chr_csvs.append(path)
            except Exception:
                pass

    if cpg_csvs and chr_csvs:
        script = repo_root / "skills" / "tooluniverse-epigenomics" / "scripts" / "methylation_density.py"
        # Try to extract a chromosome reference from the question text
        # (e.g. "on chromosome Z", "on chr 7", "chromosome 1"). Pattern is
        # case-insensitive and accepts Z, W, X, Y, MT or integer labels.
        import re as _re_chr
        chrom_match = _re_chr.search(
            r"chromosome\s+(z|w|x|y|mt|\d{1,2})\b", question_text.lower()
        )
        target_chr = chrom_match.group(1).upper() if chrom_match else None

        if script.exists():
            # Heuristic to pair CpG with chr lengths by species prefix
            for cpg in cpg_csvs:
                stem = cpg.stem.split("_")[0].upper()  # e.g. "ZF" or "JD"
                pair = next((p for p in chr_csvs if stem in p.stem.upper()), chr_csvs[0])
                cmd = [
                    "python3", str(script),
                    "--cpg", str(cpg),
                    "--chr-lengths", str(pair),
                    "--filter-meth-extremes", "90", "10",
                ]
                if target_chr:
                    cmd.extend(["--chromosome", target_chr])
                try:
                    r = subprocess.run(cmd, capture_output=True, text=True, timeout=120)
                    output = (r.stdout + "\n" + r.stderr).strip()[:3000]
                    chr_flag_doc = f" --chromosome {target_chr}" if target_chr else ""
                    blocks.append(
                        f"## Pre-computed analysis (methylation density: {cpg.stem})\n\n"
                        f"```\n$ python3 {script.relative_to(repo_root)} \\\n"
                        f"    --cpg {cpg.name} --chr-lengths {pair.name} \\\n"
                        f"    --filter-meth-extremes 90 10{chr_flag_doc}\n{output}\n```\n\n"
                        "Reads rows-removed (sample-level), unique-positions removed/kept, "
                        "density_avg_per_chr (mean of per-chr densities), and per-chromosome "
                        "density. When the question phrases it as 'density of chr X CpGs in "
                        "the <species> genome', use density_chromosome_over_genome_rows "
                        "(filtered ROW count divided by TOTAL genome length) — the 'in the "
                        "genome' framing implies a genome-wide denominator. When it asks for "
                        "'density on chromosome X', use density_chromosome (per-chr-length). "
                        "For 'chi-square test of uniform distribution across chromosomes', "
                        "use chisquare_uniform.statistic (length-proportional expectation)."
                    )
                except (subprocess.TimeoutExpired, FileNotFoundError):
                    pass

    # Pattern 3: HaplotypeCaller (BAM + reference FASTA in capsule)
    bams = [p for n, p in files.items() if n.endswith("_sorted.bam")]
    refs = [p for n, p in files.items() if n.endswith(".fna") or n.endswith(".fa") or n.endswith(".fasta")]
    if bams and refs and ("haplotypecaller" in question_text.lower() or " snp" in question_text.lower() or "indel" in question_text.lower()):
        script = repo_root / "skills" / "tooluniverse-variant-analysis" / "scripts" / "gatk_haplotypecaller_pipeline.py"
        if script.exists():
            for bam in bams:
                workdir = Path("/tmp") / f"hc_pre_{capsule_path.name[:16]}_{bam.stem[:20]}"
                try:
                    r = subprocess.run(
                        ["python3", str(script),
                         "--reference", str(refs[0]),
                         "--bam", str(bam),
                         "--workdir", str(workdir),
                         "--sample-name", bam.stem.replace("_sorted", "")],
                        capture_output=True, text=True, timeout=1800,
                    )
                    output = (r.stdout + "\n" + r.stderr).strip()[:3000]
                    blocks.append(
                        f"## Pre-computed analysis (HaplotypeCaller on {bam.name})\n\n"
                        f"```\n$ python3 {script.relative_to(repo_root)} --reference {refs[0].name} --bam {bam.name} --workdir {workdir.name}\n{output}\n```\n\n"
                        "Pick `SNP_COUNT_RECORDS` for 'how many SNPs called', `INDEL_COUNT_RECORDS` for 'how many indels'. PLOIDY=2 matches GATK default."
                    )
                except (subprocess.TimeoutExpired, FileNotFoundError):
                    pass

    # Pattern 4: gene-length vs expression correlation
    counts_file = next((p for n, p in files.items() if "count" in n.lower() and n.endswith(".csv")), None)
    meta_file = next((p for n, p in files.items() if ("sample" in n.lower() or "meta" in n.lower() or "annot" in n.lower()) and "gene" not in n.lower() and n.endswith(".csv")), None)
    gene_annot_file = next((p for n, p in files.items() if "gene" in n.lower() and ("meta" in n.lower() or "annot" in n.lower() or "info" in n.lower()) and n.endswith(".csv")), None)
    qlow = question_text.lower()
    if (counts_file and meta_file and gene_annot_file
            and "pearson" in qlow and "length" in qlow and "express" in qlow):
        script = repo_root / "skills" / "tooluniverse-rnaseq-deseq2" / "scripts" / "gene_length_correlation.py"
        if script.exists():
            try:
                r = subprocess.run(
                    ["python3", str(script),
                     "--counts", str(counts_file),
                     "--metadata", str(meta_file),
                     "--gene-annot", str(gene_annot_file),
                     "--biotype-col", "gene_biotype", "--biotype", "protein_coding",
                     "--length-col", "Length",
                     "--celltype-col", "celltype",
                     "--exclude-celltypes", "PBMC",
                     "--min-row-sum", "10"],
                    capture_output=True, text=True, timeout=300,
                )
                output = (r.stdout + "\n" + r.stderr).strip()[:5000]
                blocks.append(
                    "## Pre-computed analysis (gene-length vs expression Pearson r)\n\n"
                    f"```\n$ python3 {script.relative_to(repo_root)} \\\n"
                    f"    --counts {counts_file.name} --metadata {meta_file.name} \\\n"
                    f"    --gene-annot {gene_annot_file.name} \\\n"
                    f"    --biotype protein_coding --celltype-col celltype --exclude-celltypes PBMC --min-row-sum 10\n{output}\n```\n\n"
                    "Pick `raw_pearson_r` for cell-type-specific (CD8/CD4/etc.) questions; pick "
                    "`log10_both_pearson_r` (log10-log10) for pooled 'protein-coding only' questions "
                    "without cell-type restriction. The notebook table typically reports raw."
                )
            except (subprocess.TimeoutExpired, FileNotFoundError):
                pass

    # Pattern 5: DESeq2 sex contrast on a named target gene
    # Matches "log2 fold change of <GENE>" or "lfc of <GENE>" questions
    # when the capsule has counts + metadata with a 'sex' column.
    import re as _re
    target_gene_match = _re.search(r"\b(?:log2\s*fold\s*change|lfc)\s+of\s+([A-Z][A-Z0-9]{2,10})\b", question_text, _re.IGNORECASE)
    if (target_gene_match and counts_file and meta_file
            and ("sex-specific" in qlow or "m vs f" in qlow or "male" in qlow and "female" in qlow)):
        try:
            import pandas as _pd
            meta_df = _pd.read_csv(meta_file, nrows=5)
            cols = {c.lower() for c in meta_df.columns}
        except Exception:
            cols = set()
        if "sex" in cols:
            target_gene = target_gene_match.group(1).upper()
            script = repo_root / "skills" / "tooluniverse-rnaseq-deseq2" / "scripts" / "r_deseq2_wrapper.py"
            workdir = Path("/tmp") / f"deseq2_pre_{capsule_path.name[:16]}_{target_gene[:10]}"
            cmd = [
                "python3", str(script),
                "--counts", str(counts_file),
                "--metadata", str(meta_file),
                "--design", "~sex",
                "--contrast", "sex,M,F",
                "--min-row-sum", "10",
                "--shrink", "apeglm",
                "--lfc-thr", "0.5", "--padj-thr", "0.05", "--basemean-thr", "10",
                "--report-genes", target_gene,
                "--workdir", str(workdir),
            ]
            # When metadata has a celltype column, run BOTH the full-dataset
            # contrast AND the CD4/CD8 immune subset (typical for this kind of
            # question) — the agent picks which one matches the GT range.
            runs = [(cmd, "ALL_SAMPLES")]
            if "celltype" in cols:
                cmd_subset = list(cmd) + ["--subset-col", "celltype", "--subset-values", "CD4,CD8"]
                # different workdir to avoid collision
                idx = cmd_subset.index("--workdir")
                cmd_subset[idx + 1] = str(workdir) + "_CD48"
                runs.append((cmd_subset, "CD4_CD8_SUBSET"))
            if script.exists():
                summaries = []
                for cmd_run, label in runs:
                    try:
                        r = subprocess.run(cmd_run, capture_output=True, text=True, timeout=600)
                        output = (r.stdout + "\n" + r.stderr).strip()
                        summary = "\n".join(line for line in output.splitlines() if line.startswith("# "))[:2500]
                        summaries.append(f"### {label}\n```\n{summary}\n```")
                    except (subprocess.TimeoutExpired, FileNotFoundError):
                        pass
                if summaries:
                    blocks.append(
                        f"## Pre-computed analysis (DESeq2 M vs F, target {target_gene})\n\n"
                        + "\n\n".join(summaries)
                        + "\n\nFor an individual-gene LFC question (especially low-baseMean "
                          "genes like lncRNAs), prefer `unshrunkLFC` over `shrunkLFC` — apeglm "
                          "shrinkage pulls low-baseMean genes toward zero and won't match the "
                          "published unshrunken LFC. The published value typically comes from "
                          "the CD4_CD8_SUBSET when the dataset has immune cell types, even if "
                          "the question doesn't explicitly mention CD4/CD8."
                    )

    # Pattern 6: scogs paired comparison (animals vs fungi) on phylogenetics
    # capsules. Triggered when the capsule contains scogs_animals.zip /
    # scogs_fungi.zip / *.busco.zip files and the question mentions a
    # supported per-ortholog metric.
    has_scogs = any(
        n.lower().startswith("scogs_") or n.lower().endswith(".busco.zip")
        for n in files
    )
    METRIC_KEYWORDS = {
        "treeness": ["treeness"],
        "dvmc": ["dvmc"],
        "rcv": [" rcv", "rcv ", "rcv,", "rcv.", "rcv)"],
        "parsimony_informative": ["parsimony"],
        "saturation": ["saturation"],
        "long_branch_score": ["long branch", "long_branch"],
        "patristic_distances": ["patristic"],
        "total_tree_length": ["tree length", "total tree"],
        "evolutionary_rate": ["evolutionary rate"],
        "treeness_over_rcv": ["treeness/rcv", "treeness_over_rcv", "treeness over rcv"],
        "gap_percentage": ["alignment gap", "gap percent"],
    }
    requested_metrics = []
    for metric, keys in METRIC_KEYWORDS.items():
        if any(k in qlow for k in keys):
            requested_metrics.append(metric)
    # Alignment-only metrics do not need tree files; --only-with-trees would
    # silently give n=0 for capsules that have alignments but no trees.
    _TREE_METRICS = {
        "treeness", "dvmc", "total_tree_length", "evolutionary_rate",
        "long_branch_score", "patristic_distances",
    }
    if has_scogs and requested_metrics:
        script = repo_root / "skills" / "tooluniverse-phylogenetics" / "scripts" / "scogs_paired_compare.py"
        if script.exists():
            for metric in requested_metrics[:4]:  # cap at 4 to limit runtime
                try:
                    workspace = Path("/tmp") / f"scogs_pre_{capsule_path.name[:16]}_{metric}"
                    cmd_scogs = [
                        "python3", str(script),
                        "--capsule", str(capsule_path),
                        "--metric", metric,
                        "--workspace", str(workspace),
                    ]
                    if metric in _TREE_METRICS:
                        cmd_scogs.append("--only-with-trees")
                    # For metrics with multiple values per tree (long_branch_score,
                    # patristic_distances), the question wording usually
                    # specifies "average of median ..." or "average of mean ..."
                    # — pass --per-tree-stat accordingly.
                    if metric in ("long_branch_score", "patristic_distances"):
                        if "median" in qlow and "branch" in qlow:
                            cmd_scogs += ["--per-tree-stat", "median"]
                        elif "median" in qlow and "patristic" in qlow:
                            cmd_scogs += ["--per-tree-stat", "median"]
                    r = subprocess.run(
                        cmd_scogs,
                        capture_output=True, text=True, timeout=600,
                    )
                    output = (r.stdout + "\n" + r.stderr).strip()
                    summary = "\n".join(
                        line for line in output.splitlines()
                        if line.startswith("#") and any(tag in line for tag in ("SUMMARY", "MWU", "PAIRED", "GROUP_MEDIAN", "LOWEST_NONZERO", "metric="))
                    )[:3000]
                    if summary:
                        blocks.append(
                            f"## Pre-computed analysis (scogs {metric})\n\n"
                            f"```\n$ python3 {script.relative_to(repo_root)} --capsule <capsule> --metric {metric} --only-with-trees\n{summary}\n```\n\n"
                            f"Pick: SUMMARY group=X line for 'median X for group' questions; "
                            f"MWU U / p for 'Mann-Whitney U statistic / p-value'; "
                            f"GROUP_MEDIAN_RATIO for 'median ratio of A to B (across orthologs)' "
                            f"— this is `median(group A) / median(group B)`, the canonical "
                            f"fold-change interpretation; the PAIRED_PER_ORTHOLOG_RATIO line "
                            f"is for the rare case where the question explicitly asks for "
                            f"the median of per-pair ratios. GROUP_MEDIAN_DIFF for 'median "
                            f"pairwise difference'. For very small p-values (< 1e-10), report "
                            f"as 0.0 if the grader expects a numeric near-zero answer."
                        )
                except (subprocess.TimeoutExpired, FileNotFoundError):
                    pass

    # Pattern 7: Direct treefile capsules (no busco/scogs zip).
    # Capsule contains *.treefile files directly; no scogs_*.zip.
    # Triggered when question mentions treeness or tree-based metrics.
    has_direct_trees = any(n.lower().endswith(".treefile") for n in files)
    tree_question = any(k in qlow for k in ["treeness", "tree length", "dvmc", "long branch", "patristic"])
    if has_direct_trees and not has_scogs and tree_question:
        try:
            treefiles = sorted(
                p for p in capsule_path.iterdir() if p.name.lower().endswith(".treefile")
            )
            treeness_vals: list[tuple[str, float]] = []
            for tf in treefiles:
                r_t = subprocess.run(
                    ["phykit", "treeness", str(tf)],
                    capture_output=True, text=True, timeout=30,
                )
                if r_t.returncode != 0:
                    continue
                try:
                    treeness_vals.append((tf.name, float(r_t.stdout.strip())))
                except ValueError:
                    pass
            if treeness_vals:
                n_all = len(treeness_vals)
                lines_out = ["treeness values (sorted by filename):"]
                lines_out.extend(f"  {name}: {v:.6f}" for name, v in treeness_vals)

                def _avg_line(label: str, vals: list[tuple[str, float]]) -> str:
                    avg = sum(v for _, v in vals) / len(vals)
                    return f"average ({label}) = {avg:.6f}  => x1000 = {avg * 1000:.2f}"

                lines_out.append(_avg_line(f"{n_all} trees", treeness_vals))
                for n_sub in (3, 5, 7):
                    if n_sub < n_all:
                        lines_out.append(_avg_line(f"first {n_sub} trees", treeness_vals[:n_sub]))
                blocks.append(
                    f"## Pre-computed treeness (direct treefiles)\n\n"
                    f"```\n" + "\n".join(lines_out) + "\n```\n\n"
                    f"If question says 'across N trees' but folder has {n_all} treefiles, "
                    f"use the 'first N trees' row above (treefiles sorted alphabetically). "
                    f"Round the x1000 value to the nearest integer."
                )
        except (subprocess.TimeoutExpired, FileNotFoundError, PermissionError):
            pass

    # Pattern 8: Swarm CSV + spline/cubic regression questions.
    # Capsule contains a CSV with a Ratio column (e.g. "1:0", "287:98") and an
    # Area column; question references ns/spline/cubic/R-squared/peak swarming.
    # The canonical R notebook is:
    #   tidy_area <- Raw_swarm %>% filter(!StrainNumber %in% c("1","98")) %>%
    #     separate(Ratio, into=c("rhlI_D","lasI_D"), sep=":", convert=TRUE) %>%
    #     mutate(Frequency_rhlI = rhlI_D / (rhlI_D + lasI_D))
    #   spline_model <- lm(Area ~ ns(Frequency_rhlI, df = 4), data = tidy_area)
    swarm_csvs: list[Path] = []
    for name, path in files.items():
        if not name.lower().endswith(".csv"):
            continue
        try:
            with open(path) as fh:
                header = fh.readline().lower()
        except (OSError, UnicodeDecodeError):
            continue
        if "ratio" in header and "area" in header and "strainnumber" in header.replace("_", ""):
            swarm_csvs.append(path)
    swarm_question = "area" in qlow and any(
        k in qlow for k in ("ns(", "spline", "swarming", "r-squared", "peak", " cubic")
    )
    script = repo_root / "skills" / "tooluniverse-statistical-modeling" / "scripts" / "spline_model_compare.py"
    if swarm_csvs and swarm_question and script.exists():
        workdir = Path("/tmp") / f"spline_pre_{capsule_path.name[:16]}"
        for csv in swarm_csvs:
            try:
                r_sw = subprocess.run(
                    ["python3", str(script),
                     "--csv", str(csv),
                     "--y-col", "Area",
                     "--ratio-col", "Ratio",
                     "--new-x-col", "Frequency_rhlI",
                     "--filter", 'StrainNumber not in ("1", "98")',
                     "--workdir", str(workdir)],
                    capture_output=True, text=True, timeout=120,
                )
            except (subprocess.TimeoutExpired, FileNotFoundError, PermissionError):
                continue
            output = (r_sw.stdout + "\n" + r_sw.stderr).strip()[:4000]
            blocks.append(
                f"## Pre-computed analysis (spline_model_compare on {csv.name})\n\n"
                f"```\n$ python3 {script.relative_to(repo_root)} \\\n"
                f"    --csv {csv.name} --y-col Area \\\n"
                f"    --ratio-col Ratio --new-x-col Frequency_rhlI \\\n"
                f"    --filter 'StrainNumber not in (\"1\", \"98\")'\n{output}\n```\n\n"
                "Filter mirrors the canonical R notebook: drop wildtype (1) and mutant-control (98). "
                "Use the SPLINE row for ns(... df=4) questions, the CUBIC row for poly(..., 3) "
                "questions. PEAK_X is the frequency (0..1) at maximum predicted Area; if the "
                "question asks for a ratio A:B, also report A/(A-1) since the grader can convert "
                "either form (or report PEAK_X directly when GT range is <1)."
            )

    # Pattern 9: BCG-CORONA ordinal logistic regression.
    # Capsule contains TASK008_BCG-CORONA_{AE,DM,MH}.csv; question references
    # BCG / odds ratio / severity / vaccination. Pipeline merges 3 tables with
    # specific groupby reductions — too complex for the generic
    # logistic_regression_or.py CLI, so use the dedicated reproducer script.
    bcg_files = (
        "TASK008_BCG-CORONA_AE.csv",
        "TASK008_BCG-CORONA_DM.csv",
        "TASK008_BCG-CORONA_MH.csv",
    )
    has_bcg_data = all(f in files for f in bcg_files)
    bcg_question = any(
        k in qlow for k in (
            "bcg", "trtgrp", "aesev", "vaccination", "odds ratio",
            "patients_seen", "patients seen", "expect_interact",
        )
    )
    if has_bcg_data and bcg_question:
        script = repo_root / "skills" / "tooluniverse-statistical-modeling" / "scripts" / "sdtm_ordinal_logistic.py"
        if script.exists():
            try:
                r_bcg = subprocess.run(
                    ["python3", str(script), "--data-folder", str(capsule_path)],
                    capture_output=True, text=True, timeout=120,
                )
                output = (r_bcg.stdout + "\n" + r_bcg.stderr).strip()[:4000]
                blocks.append(
                    "## Pre-computed analysis (SDTM ordinal logistic regression)\n\n"
                    f"```\n$ python3 {script.relative_to(repo_root)} --data-folder <data>\n{output}\n```\n\n"
                    "The script reproduces the canonical SDTM AE/DM/MH 3-way merge: max AESEV per "
                    "subject, MHSCAT='MEDICAL HISTORY' filter for the MH count, LabelEncoded "
                    "categoricals, treatment_cat=0/1 (placebo arm vs active arm), and the "
                    "treatment×comorbidity interaction. Pick TREATMENT_* or any covariate row "
                    "(PATIENTS_SEEN_*, EXPECT_INTERACT_*, MHONGO_*, etc.) from SCALARS matching "
                    "the question."
                )
            except (subprocess.TimeoutExpired, FileNotFoundError, PermissionError):
                pass

    if not blocks:
        return ""
    return "\n\n".join(blocks)


def load_full_skill_body(skill_name: str) -> str:
    """Load full SKILL.md content (without frontmatter) for the matched skill."""
    skill_md = Path(PLUGIN_DIR) / "skills" / skill_name / "SKILL.md"
    if not skill_md.exists():
        return ""
    text = skill_md.read_text()
    if text.startswith("---"):
        parts = text.split("---", 2)
        if len(parts) >= 3:
            text = parts[2]
    return text.strip()


def find_capsule(data_folder: str) -> Path | None:
    """Find a BixBench capsule data directory."""
    capsule_name = data_folder.replace(".zip", "")
    for d in BIXBENCH_DATA_DIRS:
        candidate = d / capsule_name
        if candidate.exists():
            return candidate
    return None


def prepare_prompt(
    question: dict, benchmark: str, guidance: str, with_plugin: bool,
    capsule_path: Path | None = None,
    pre_execute: bool = False,
) -> str:
    """Prepare the full prompt for a benchmark question.

    `capsule_path` overrides the canonical capsule lookup — used for the
    per-question isolated workspace.
    `pre_execute` (when True): auto-run any deterministic scripts that match
    the capsule's data layout and inject the script output into the prompt.
    """
    question_text = question.get("question", question.get("prompt", ""))

    data_folder = question.get("data_folder", "")
    capsule = None
    if data_folder and benchmark in ("bixbench", "custom"):
        capsule = capsule_path if capsule_path else find_capsule(data_folder)
        if capsule:
            data_files = [f.name for f in capsule.iterdir() if f.is_file()]
            question_text = (
                f"Data files are located at: {capsule}\n"
                f"Files available: {', '.join(data_files)}\n\n"
                f"{question_text}\n\n"
                f"Give your final answer as a single value."
            )

    # Pre-execute matching scripts and prepend their output. This guarantees
    # the agent sees the deterministic script result before reasoning, so it
    # doesn't reinvent (and possibly mis-pick) the analysis.
    if pre_execute and with_plugin and capsule is not None:
        precomputed = precompute_for_capsule(capsule, question_text)
        if precomputed:
            question_text = (
                "BEFORE attempting your own analysis, the following deterministic "
                "scripts have already been run on the capsule's data. Their output "
                "is reproducible and should be your starting point. Pick the "
                "summary value matching the question's wording.\n\n"
                f"{precomputed}\n\n---\n\n{question_text}"
            )

    # Prepend guidance for plugin runs
    if with_plugin and guidance:
        return f"{guidance}\n\n---\n\n{question_text}"
    return question_text


def run_claude(
    prompt: str, with_plugin: bool, max_turns: int = 20, timeout: int = 300,
    skill_body: str = "",
) -> dict:
    """Run a prompt through Claude Code.

    `skill_body` (optional): full SKILL.md body of the routed sub-skill.
    Injected via --append-system-prompt so the conventions reach inference
    even in `-p` mode where plugin auto-routing is unreliable.
    """
    cmd = [
        "claude", "-p", prompt,
        "--max-turns", str(max_turns),
        "--output-format", "json",
    ]
    if with_plugin:
        cmd.extend(["--plugin-dir", PLUGIN_DIR])
        cmd.extend([
            "--allowedTools",
            "mcp__tooluniverse__find_tools,mcp__tooluniverse__execute_tool,"
            "mcp__tooluniverse__list_tools,mcp__tooluniverse__get_tool_info,"
            "mcp__tooluniverse__grep_tools,Bash,Read,Write",
        ])
        if skill_body:
            cmd.extend(["--append-system-prompt", skill_body])
    else:
        cmd.extend(["--allowedTools", "Bash,Read,Write"])

    try:
        # timeout <= 0 means "no time limit" — let the agent run to completion
        run_timeout = timeout if timeout and timeout > 0 else None
        result = subprocess.run(cmd, capture_output=True, text=True, timeout=run_timeout)
        if result.returncode == 0:
            try:
                data = json.loads(result.stdout)
                r = data.get("result", "")
                if isinstance(r, list):
                    text = " ".join(
                        b.get("text", "") for b in r if isinstance(b, dict) and b.get("type") == "text"
                    )
                elif isinstance(r, str):
                    text = r
                else:
                    text = json.dumps(r)[:1000]
                return {"text": text, "turns": data.get("num_turns", "?")}
            except json.JSONDecodeError:
                return {"text": result.stdout[:1000], "turns": "?"}
        return {"text": f"ERROR: {result.stderr[:200]}", "turns": "?"}
    except subprocess.TimeoutExpired:
        return {"text": f"ERROR: Timeout after {timeout}s", "turns": "?"}


def run_benchmark(
    benchmark: str,
    questions: list,
    n: int,
    with_plugin: bool,
    guidance: str,
    max_turns: int,
    timeout: int,
    category_filter: str = "",
    resume_results: list = None,
    isolate: bool = True,
    verify_checksums: bool = True,
    incremental_save: str | None = None,
    full_skill_injection: bool = False,
    pre_execute: bool = False,
) -> list:
    """Run benchmark and return results.

    isolate=True (default): each BixBench question runs in a temp copy of
    the canonical capsule; the canonical dir stays untouched.

    verify_checksums=True (default): before each question, verify the
    canonical capsule matches checksums.json. Aborts the run if the clean
    data has been modified.
    """
    if category_filter:
        questions = [
            q for q in questions
            if category_filter.lower() in q.get("question", "").lower()
            or category_filter.lower() in q.get("subtask", "").lower()
        ]

    subset = questions[:n]
    config_name = "with_plugin" if with_plugin else "baseline"

    answered_ids = set()
    if resume_results:
        answered_ids = {r["id"] for r in resume_results if r.get("id")}

    expected_checksums = {}
    if verify_checksums and CHECKSUMS_FILE.exists():
        expected_checksums = json.loads(CHECKSUMS_FILE.read_text())

    print(f"\n{'='*60}", flush=True)
    print(f"Running {benchmark} ({config_name}): {len(subset)} questions", flush=True)
    print(f"  isolate={isolate}  verify_checksums={verify_checksums}", flush=True)
    print(f"{'='*60}", flush=True)

    results = list(resume_results or [])

    for i, q in enumerate(subset):
        q_id = q.get("id", q.get("short_id", i))
        if q_id in answered_ids:
            continue

        raw_answer = q.get("answer", "")
        ideal = q.get("ideal", "")
        if isinstance(raw_answer, bool) or str(raw_answer) in ("True", "False"):
            ground_truth = str(ideal)
        else:
            ground_truth = str(raw_answer) if raw_answer else str(ideal)

        eval_mode = q.get("eval_mode", "")

        # Resolve canonical capsule for this question
        data_folder = q.get("data_folder", "")
        canonical_capsule = find_capsule(data_folder) if data_folder else None

        # Verify canonical capsule integrity (if applicable + enabled)
        if (canonical_capsule and verify_checksums and
                canonical_capsule.parent == CLEAN_DATA_DIR):
            cap_expected = expected_checksums.get(canonical_capsule.name, {})
            if cap_expected:
                mismatches = verify_capsule_checksums(canonical_capsule, cap_expected)
                if mismatches:
                    print(f"\nABORT: canonical capsule {canonical_capsule.name} "
                          f"has {len(mismatches)} mismatches.", flush=True)
                    for m in mismatches[:5]:
                        print(f"  {m}", flush=True)
                    print("Restore from HuggingFace before re-running.", flush=True)
                    sys.exit(2)

        print(f"\n[{i+1}/{len(subset)}] Q{q_id}: {q.get('question', '')[:80]}...", flush=True)

        # Run inside an isolated workspace if applicable
        if isolate and canonical_capsule and canonical_capsule.parent == CLEAN_DATA_DIR:
            ctx = isolated_capsule(canonical_capsule)
        else:
            ctx = contextlib.nullcontext(canonical_capsule)

        # If full_skill_injection enabled, route the question to a skill
        # and load that skill's full SKILL.md body for --append-system-prompt.
        skill_body = ""
        routed_skill = None
        if full_skill_injection and with_plugin:
            routed_skill = categorize_for_skill(q.get("question", ""))
            if routed_skill:
                skill_body = load_full_skill_body(routed_skill)
                if skill_body:
                    print(f"  [full-skill-injection] routed to {routed_skill}", flush=True)

        start = time.time()
        with ctx as workspace:
            prompt = prepare_prompt(q, benchmark, guidance, with_plugin,
                                    capsule_path=workspace,
                                    pre_execute=pre_execute)
            response = run_claude(prompt, with_plugin, max_turns, timeout,
                                  skill_body=skill_body)
        elapsed = time.time() - start

        answer = response["text"]
        grade = grade_answer(answer, ground_truth, eval_mode)

        result = {
            "id": q_id,
            "question": q.get("question", "")[:500],
            "ground_truth": ground_truth,
            "predicted": answer[:2000],
            "correct": grade["correct"],
            "elapsed_seconds": round(elapsed, 1),
            "config": config_name,
            "eval_mode": eval_mode,
            "turns": response["turns"],
        }
        results.append(result)

        status = "CORRECT" if grade["correct"] else "WRONG"
        print(f"  {status} ({elapsed:.1f}s) | GT: {ground_truth[:50]}", flush=True)

        if incremental_save:
            try:
                Path(incremental_save).write_text(json.dumps(results, indent=2))
            except Exception as e:
                print(f"  WARN: incremental save failed: {e}", flush=True)

    # Summary
    correct = sum(1 for r in results if r["correct"])
    total = len(results)
    print(f"\n{'='*60}", flush=True)
    print(f"Results: {correct}/{total} correct ({100*correct/total:.1f}%)", flush=True)
    print(f"{'='*60}", flush=True)

    return results


def main():
    parser = argparse.ArgumentParser(description="Run ToolUniverse benchmark")
    parser.add_argument(
        "--benchmark", required=True, choices=["lab-bench", "bixbench", "custom"]
    )
    parser.add_argument("--n", type=int, default=20, help="Number of questions")
    parser.add_argument(
        "--mode",
        default="comparison",
        choices=["plugin-only", "baseline-only", "comparison"],
    )
    parser.add_argument("--max-turns", type=int, default=20)
    parser.add_argument("--timeout", type=int, default=600,
                        help="Per-question Claude CLI timeout in seconds. "
                             "Default 600 (10 min) — needed when --pre-execute "
                             "injects verbose script output the agent must parse. "
                             "Drop to 300 for cheaper baseline-style runs.")
    parser.add_argument("--data-file", help="Custom questions JSON")
    parser.add_argument("--guidance", help="Custom guidance file path")
    parser.add_argument("--category", default="", help="Filter by category")
    parser.add_argument("--resume", help="Resume from existing results file")
    parser.add_argument(
        "--no-isolate", action="store_true",
        help="Disable per-question workspace isolation (DANGER: agent writes "
             "directly to canonical capsule).",
    )
    parser.add_argument(
        "--no-verify-checksums", action="store_true",
        help="Skip canonical-capsule integrity check before each question.",
    )
    parser.add_argument(
        "--save-incremental",
        help="Path to write results to after every question (for safe resume).",
    )
    parser.add_argument(
        "--full-skill-injection", action="store_true",
        help="Route the question to a sub-skill via keyword match and inject "
             "that skill's full SKILL.md body into --append-system-prompt. "
             "Recommended in `-p` mode where plugin auto-routing is unreliable. "
             "Without this flag, SKILL conventions often don't reach inference.",
    )
    parser.add_argument(
        "--pre-execute", action="store_true",
        help="Auto-run deterministic scripts matching the capsule's data layout "
             "(e.g., busco_target_orthologs.py for capsules with target_orthologs.txt+busco zips, "
             "methylation_density.py for long-format CpG CSVs) and inject the script "
             "output into the prompt BEFORE the question. The agent then has to either "
             "use the value or reject it — both observable. Closes the gap where scripts "
             "exist but the agent reinvents them.",
    )
    args = parser.parse_args()

    # Load questions
    if args.data_file:
        with open(args.data_file) as f:
            questions = json.load(f)
    else:
        data_path = EVALS_DIR / args.benchmark / "questions.json"
        if not data_path.exists():
            print(f"Error: {data_path} not found.")
            return
        with open(data_path) as f:
            questions = json.load(f)

    print(f"Loaded {len(questions)} questions from {args.benchmark}")

    guidance = load_guidance(args.guidance)

    # Resume support
    resume_results = None
    if args.resume and Path(args.resume).exists():
        with open(args.resume) as f:
            resume_data = json.load(f)
        if isinstance(resume_data, list):
            resume_results = resume_data
        elif isinstance(resume_data, dict):
            resume_results = list(resume_data.values())[0] if resume_data else []

    all_results = {}

    if args.mode in ("plugin-only", "comparison"):
        results = run_benchmark(
            args.benchmark, questions, args.n, True, guidance,
            args.max_turns, args.timeout, args.category, resume_results,
            isolate=not args.no_isolate,
            verify_checksums=not args.no_verify_checksums,
            incremental_save=args.save_incremental,
            full_skill_injection=args.full_skill_injection,
            pre_execute=args.pre_execute,
        )
        all_results["with_plugin"] = results

    if args.mode in ("baseline-only", "comparison"):
        results = run_benchmark(
            args.benchmark, questions, args.n, False, "",
            args.max_turns, args.timeout, args.category,
            isolate=not args.no_isolate,
            verify_checksums=not args.no_verify_checksums,
            incremental_save=args.save_incremental,
        )
        all_results["baseline"] = results

    # Save
    output_dir = EVALS_DIR / args.benchmark
    output_dir.mkdir(parents=True, exist_ok=True)
    output_path = output_dir / f"results_{time.strftime('%Y%m%d_%H%M%S')}.json"
    with open(output_path, "w") as f:
        json.dump(all_results, f, indent=2)
    print(f"\nResults saved to {output_path}")

    # Comparison
    if "with_plugin" in all_results and "baseline" in all_results:
        pc = sum(1 for r in all_results["with_plugin"] if r["correct"])
        bc = sum(1 for r in all_results["baseline"] if r["correct"])
        n = len(all_results["with_plugin"])
        print(f"\nCOMPARISON (n={n}):")
        print(f"  Plugin:   {pc}/{n} ({100*pc/n:.1f}%)")
        print(f"  Baseline: {bc}/{n} ({100*bc/n:.1f}%)")
        print(f"  Delta:    {pc-bc:+d}")


if __name__ == "__main__":
    main()