#!/usr/bin/env python3
"""Reproducible analysis of the frozen MCP probe dataset.

Reads mcp-probe-2026-09-14.csv from this script's own directory, verifies the
expected SHA-256 digest and column header before calculating anything, then
prints the figures cited in the whitepaper "MCP Tool Readiness".

Standard library only. No network access. No writes.

Dataset: Cracked Engineering, "The state of public MCP servers, September 2026",
public probe export retrieved 14 September 2026 and licensed CC BY 4.0.

Usage:
    python analyze_probe.py
"""

import csv
import hashlib
import os
import sys

DATA_FILENAME = "mcp-probe-2026-09-14.csv"

EXPECTED_SHA256 = (
    "4e84191c8feda3276fb9493f9f0b34644a90b88c634ca0a99b381c0e0b2475f7"
)

EXPECTED_COLUMNS = [
    "server",
    "title",
    "url",
    "website",
    "tools",
    "list_ok",
    "call_ok",
    "list_ms",
    "call_ms",
    "probed_tool",
    "probed_at",
    "streak",
]

# Bin edges used by the published report, in milliseconds.
LATENCY_BINS = [
    ("under 250 ms", 0, 250),
    ("250-499 ms", 250, 500),
    ("500-999 ms", 500, 1000),
    ("1.0-1.9 seconds", 1000, 2000),
    ("2.0-4.9 seconds", 2000, 5000),
    ("5.0-9.9 seconds", 5000, 10000),
    ("10 seconds or more", 10000, None),
]


class AnalysisError(Exception):
    """Raised when the dataset does not match the frozen research artifact."""


def script_directory():
    return os.path.dirname(os.path.abspath(__file__))


def data_path():
    return os.path.join(script_directory(), DATA_FILENAME)


def sha256_of_file(path):
    digest = hashlib.sha256()
    with open(path, "rb") as handle:
        for block in iter(lambda: handle.read(65536), b""):
            digest.update(block)
    return digest.hexdigest()


def verify_checksum(path):
    if not os.path.isfile(path):
        raise AnalysisError(
            "Dataset not found next to this script.\n"
            "  expected file: %s\n"
            "  Place the frozen CSV in the same directory as analyze_probe.py."
            % path
        )
    actual = sha256_of_file(path)
    if actual != EXPECTED_SHA256:
        raise AnalysisError(
            "SHA-256 mismatch. This is not the frozen research artifact.\n"
            "  expected: %s\n"
            "  actual:   %s\n"
            "  The upstream export is a rolling file. Results printed from a\n"
            "  different copy are not the results cited in the paper."
            % (EXPECTED_SHA256, actual)
        )
    return actual


def verify_columns(fieldnames):
    if fieldnames != EXPECTED_COLUMNS:
        raise AnalysisError(
            "Column header mismatch. The schema this analysis depends on has "
            "changed.\n"
            "  expected: %s\n"
            "  actual:   %s"
            % (",".join(EXPECTED_COLUMNS), ",".join(fieldnames or []))
        )


def read_rows(path):
    with open(path, "r", encoding="utf-8", newline="") as handle:
        reader = csv.DictReader(handle)
        verify_columns(reader.fieldnames)
        return [row for row in reader]


def as_bool(value, column, line):
    text = (value or "").strip().lower()
    if text == "true":
        return True
    if text == "false":
        return False
    raise AnalysisError(
        "Unexpected value in column %s at data line %d: %r. "
        "Expected 'true' or 'false'." % (column, line, value)
    )


def as_int(value, column, line, allow_blank=False):
    text = (value or "").strip()
    if text == "":
        if allow_blank:
            return None
        raise AnalysisError(
            "Empty value in column %s at data line %d. An integer was expected."
            % (column, line)
        )
    try:
        return int(text)
    except ValueError:
        raise AnalysisError(
            "Non-integer value in column %s at data line %d: %r."
            % (column, line, value)
        )


def median(values):
    if not values:
        return None
    ordered = sorted(values)
    count = len(ordered)
    middle = count // 2
    if count % 2 == 1:
        return float(ordered[middle])
    return (ordered[middle - 1] + ordered[middle]) / 2.0


def nearest_rank_percentile(values, percentile):
    """Nearest-rank percentile: index = ceil(p / 100 * N), 1-based."""
    if not values:
        return None
    ordered = sorted(values)
    count = len(ordered)
    rank = -((-percentile * count) // 100)
    rank = max(1, min(count, int(rank)))
    return ordered[rank - 1]


def linear_percentile(values, percentile):
    """N-minus-one linear-interpolation percentile, reported for comparison."""
    if not values:
        return None
    ordered = sorted(values)
    count = len(ordered)
    if count == 1:
        return float(ordered[0])
    position = (percentile / 100.0) * (count - 1)
    lower = int(position)
    upper = min(lower + 1, count - 1)
    weight = position - lower
    return ordered[lower] * (1 - weight) + ordered[upper] * weight


def rate(numerator, denominator):
    if denominator == 0:
        return None
    return 100.0 * numerator / denominator


def fmt_rate(value):
    if value is None:
        return "not defined (empty denominator)"
    return "%.3f%%" % value


def fmt_ms(value):
    if value is None:
        return "not defined"
    if float(value).is_integer():
        return "%d ms" % int(value)
    return "%.1f ms" % value


def bin_label_for(value):
    for label, low, high in LATENCY_BINS:
        if value >= low and (high is None or value < high):
            return label
    return LATENCY_BINS[-1][0]


def analyse(rows):
    facts = {}

    facts["row_count"] = len(rows)
    facts["unique_servers"] = len({row["server"] for row in rows})
    facts["unique_urls"] = len({row["url"] for row in rows})

    timestamps = sorted(row["probed_at"] for row in rows if row["probed_at"])
    facts["first_timestamp"] = timestamps[0] if timestamps else None
    facts["last_timestamp"] = timestamps[-1] if timestamps else None
    facts["timestamp_count"] = len(timestamps)
    facts["probe_dates"] = sorted({stamp[:10] for stamp in timestamps})
    facts["span_seconds"] = span_seconds(
        facts["first_timestamp"], facts["last_timestamp"]
    )

    list_ok = []
    call_ok = []
    list_latencies = []
    call_latencies = []
    tools_values = []
    streaks = []
    blank_probed_tool_rows = 0
    blank_probed_tool_with_list_ok = 0

    for offset, row in enumerate(rows):
        line = offset + 2
        listed = as_bool(row["list_ok"], "list_ok", line)
        called = as_bool(row["call_ok"], "call_ok", line)
        list_ok.append(listed)
        call_ok.append(called)

        list_ms = as_int(row["list_ms"], "list_ms", line, allow_blank=True)
        call_ms = as_int(row["call_ms"], "call_ms", line, allow_blank=True)
        if listed and list_ms is not None:
            list_latencies.append(list_ms)
        if called and call_ms is not None:
            call_latencies.append(call_ms)

        tools_values.append(as_int(row["tools"], "tools", line, allow_blank=True))
        streaks.append(as_int(row["streak"], "streak", line, allow_blank=True))

        if not (row["probed_tool"] or "").strip():
            blank_probed_tool_rows += 1
            if listed:
                blank_probed_tool_with_list_ok += 1

    total = len(rows)
    listed_count = sum(1 for value in list_ok if value)
    called_count = sum(1 for value in call_ok if value)
    listed_then_failed = sum(
        1 for listed, called in zip(list_ok, call_ok) if listed and not called
    )

    facts["list_ok_count"] = listed_count
    facts["list_ok_rate"] = rate(listed_count, total)
    facts["call_ok_count"] = called_count
    facts["call_ok_rate_all"] = rate(called_count, total)
    facts["call_ok_rate_among_listed"] = rate(called_count, listed_count)
    facts["listed_then_call_failed"] = listed_then_failed
    facts["call_failures_total"] = total - called_count
    facts["list_failures_total"] = total - listed_count

    facts["successful_call_sample"] = len(call_latencies)
    facts["successful_list_sample"] = len(list_latencies)
    facts["median_call_ms"] = median(call_latencies)
    facts["median_list_ms"] = median(list_latencies)
    facts["p95_call_nearest_rank_ms"] = nearest_rank_percentile(call_latencies, 95)
    facts["p95_call_linear_ms"] = linear_percentile(call_latencies, 95)
    facts["calls_at_least_1s"] = sum(1 for value in call_latencies if value >= 1000)
    facts["fastest_call_ms"] = min(call_latencies) if call_latencies else None
    facts["slowest_call_ms"] = max(call_latencies) if call_latencies else None

    bins = {label: 0 for label, _, _ in LATENCY_BINS}
    for value in call_latencies:
        bins[bin_label_for(value)] += 1
    facts["call_latency_bins"] = bins

    observed_streaks = [value for value in streaks if value is not None]
    facts["dead_streak_rows"] = sum(1 for value in observed_streaks if value >= 3)
    streak_breakdown = {}
    for value in observed_streaks:
        if value >= 3:
            streak_breakdown[value] = streak_breakdown.get(value, 0) + 1
    facts["dead_streak_breakdown"] = dict(sorted(streak_breakdown.items()))

    observed_tools = [value for value in tools_values if value is not None]
    facts["tools_sum"] = sum(observed_tools)
    facts["tools_max"] = max(observed_tools) if observed_tools else None
    facts["tools_at_cap"] = sum(1 for value in observed_tools if value == 100)
    facts["tools_median"] = median(observed_tools)
    facts["tools_one"] = sum(1 for value in observed_tools if value == 1)
    facts["tools_over_ten"] = sum(1 for value in observed_tools if value > 10)
    facts["tools_51_to_100"] = sum(
        1 for value in observed_tools if 51 <= value <= 100
    )
    facts["empty_lists_with_list_ok"] = sum(
        1
        for listed, value in zip(list_ok, tools_values)
        if listed and value is not None and value == 0
    )
    facts["blank_website_rows"] = sum(
        1 for row in rows if not (row["website"] or "").strip()
    )
    facts["blank_list_ms_rows"] = sum(
        1 for row in rows if not (row["list_ms"] or "").strip()
    )
    facts["blank_call_ms_rows"] = sum(
        1 for row in rows if not (row["call_ms"] or "").strip()
    )

    facts["blank_probed_tool_rows"] = blank_probed_tool_rows
    facts["blank_probed_tool_with_list_ok"] = blank_probed_tool_with_list_ok

    tool_names = {}
    for row in rows:
        name = (row["probed_tool"] or "").strip()
        if name:
            tool_names[name] = tool_names.get(name, 0) + 1
    facts["distinct_probed_tools"] = len(tool_names)
    facts["top_probed_tools"] = sorted(
        tool_names.items(), key=lambda item: (-item[1], item[0])
    )[:8]

    return facts


def span_seconds(first, last):
    """Seconds between two ISO 8601 Zulu timestamps, without dateutil."""
    if not first or not last:
        return None
    return to_epoch(last) - to_epoch(first)


def to_epoch(stamp):
    text = stamp.strip()
    if not text.endswith("Z") or len(text) != 20 or text[10] != "T":
        raise AnalysisError(
            "Unexpected timestamp format in probed_at: %r. "
            "Expected YYYY-MM-DDTHH:MM:SSZ." % stamp
        )
    try:
        year = int(text[0:4])
        month = int(text[5:7])
        day = int(text[8:10])
        hour = int(text[11:13])
        minute = int(text[14:16])
        second = int(text[17:19])
    except ValueError:
        raise AnalysisError("Unparseable timestamp in probed_at: %r." % stamp)
    days = days_from_civil(year, month, day)
    return days * 86400 + hour * 3600 + minute * 60 + second


def days_from_civil(year, month, day):
    """Days since 1970-01-01, Howard Hinnant's civil calendar algorithm."""
    year -= month <= 2
    era = (year if year >= 0 else year - 399) // 400
    year_of_era = year - era * 400
    day_of_year = (153 * (month + (-3 if month > 2 else 9)) + 2) // 5 + day - 1
    day_of_era = year_of_era * 365 + year_of_era // 4 - year_of_era // 100 + day_of_year
    return era * 146097 + day_of_era - 719468


def report(facts, digest, path):
    out = sys.stdout.write

    out("MCP probe dataset analysis\n")
    out("==========================\n\n")

    out("Artifact\n")
    out("--------\n")
    out("  file:              %s\n" % os.path.basename(path))
    out("  directory:         %s\n" % os.path.dirname(path))
    out("  sha256 (verified): %s\n" % digest)
    out("  columns (%d):       %s\n" % (len(EXPECTED_COLUMNS), ",".join(EXPECTED_COLUMNS)))
    out("\n")

    out("Population\n")
    out("----------\n")
    out("  rows:                     %d\n" % facts["row_count"])
    out("  unique server names:      %d\n" % facts["unique_servers"])
    out("  unique endpoint URLs:     %d\n" % facts["unique_urls"])
    out("  probe date(s):            %s\n" % ", ".join(facts["probe_dates"]))
    out("  first probed_at:          %s\n" % facts["first_timestamp"])
    out("  last probed_at:           %s\n" % facts["last_timestamp"])
    out("  wall-clock span:          %d seconds\n" % facts["span_seconds"])
    out("  note: 986 records in 103 seconds is far shorter than sequential\n")
    out("        probing would require, so substantial concurrency is likely.\n")
    out("        The effect of that concurrency on latency is unknown.\n")
    out("\n")

    out("Completion funnel\n")
    out("-----------------\n")
    out("  list_ok = true:                       %d of %d (%s)\n"
        % (facts["list_ok_count"], facts["row_count"], fmt_rate(facts["list_ok_rate"])))
    out("  call_ok = true:                       %d of %d (%s)\n"
        % (facts["call_ok_count"], facts["row_count"], fmt_rate(facts["call_ok_rate_all"])))
    out("  call_ok among list_ok:                %d of %d (%s)\n"
        % (facts["call_ok_count"], facts["list_ok_count"],
           fmt_rate(facts["call_ok_rate_among_listed"])))
    out("  listed, then the chosen call failed:  %d\n" % facts["listed_then_call_failed"])
    out("  listing failures:                     %d\n" % facts["list_failures_total"])
    out("  call failures (all rows):             %d\n" % facts["call_failures_total"])
    out("  note: call_ok = false is a completion outcome for one keyless\n")
    out("        synthetic call. It is not downtime and not an availability SLI.\n")
    out("\n")

    out("Latency of successful operations\n")
    out("--------------------------------\n")
    out("  successful listings sampled:          %d\n" % facts["successful_list_sample"])
    out("  median successful listing:            %s\n" % fmt_ms(facts["median_list_ms"]))
    out("  successful calls sampled:             %d\n" % facts["successful_call_sample"])
    out("  median successful call:               %s\n" % fmt_ms(facts["median_call_ms"]))
    out("  P95 successful call (nearest rank):   %s\n"
        % fmt_ms(facts["p95_call_nearest_rank_ms"]))
    out("  P95 successful call (interpolated):   %s\n"
        % fmt_ms(facts["p95_call_linear_ms"]))
    out("  successful calls at or above 1 s:     %d\n" % facts["calls_at_least_1s"])
    out("  fastest successful call:              %s\n" % fmt_ms(facts["fastest_call_ms"]))
    out("  slowest successful call:              %s\n" % fmt_ms(facts["slowest_call_ms"]))
    out("  distribution of successful calls:\n")
    for label, _, _ in LATENCY_BINS:
        out("      %-20s %d\n" % (label, facts["call_latency_bins"][label]))
    out("  note: one vantage point, one short window, mixed warm and cold paths.\n")
    out("\n")

    out("Dead-streak records\n")
    out("-------------------\n")
    out("  rows with streak >= 3:                %d\n" % facts["dead_streak_rows"])
    for value, count in facts["dead_streak_breakdown"].items():
        out("      streak %-4d %d rows\n" % (value, count))
    out("  note: streaks above 3 show that removal from an active catalogue is\n")
    out("        not removal from probing or from the Official MCP Registry.\n")
    out("\n")

    out("Tool-count projection facts\n")
    out("---------------------------\n")
    out("  sum of the tools column:              %d\n" % facts["tools_sum"])
    out("  maximum value in the tools column:    %d\n" % facts["tools_max"])
    out("  rows at the value 100:                %d\n" % facts["tools_at_cap"])
    out("  median tools per row:                 %s\n" % ("%g" % facts["tools_median"]))
    out("  rows with exactly 1 tool:             %d\n" % facts["tools_one"])
    out("  rows with more than 10 tools:         %d\n" % facts["tools_over_ten"])
    out("  rows with 51 to 100 tools:            %d\n" % facts["tools_51_to_100"])
    out("  successful listings with 0 tools:     %d\n" % facts["empty_lists_with_list_ok"])
    out("  note: this column is a capped public projection. The published report\n")
    out("        cites a larger total and catalogues above 100 tools. Those\n")
    out("        figures cannot be derived from this file and are treated in the\n")
    out("        paper as report assertions, not measured facts.\n")
    out("\n")

    out("Missing values\n")
    out("--------------\n")
    out("  blank website values:                 %d\n" % facts["blank_website_rows"])
    out("  blank list_ms values:                 %d\n" % facts["blank_list_ms_rows"])
    out("  blank call_ms values:                 %d\n" % facts["blank_call_ms_rows"])
    out("\n")

    out("Probe-tool selection\n")
    out("--------------------\n")
    out("  rows with no probed_tool:             %d\n" % facts["blank_probed_tool_rows"])
    out("  of those, listing succeeded:          %d\n"
        % facts["blank_probed_tool_with_list_ok"])
    out("  distinct probed tool names:           %d\n" % facts["distinct_probed_tools"])
    out("  most frequently selected names:\n")
    for name, count in facts["top_probed_tools"]:
        out("      %-26s %d\n" % (name, count))
    out("  note: arguments are not recorded, so a success may show only that a\n")
    out("        self-test, ping or catalogue method worked.\n")
    out("\n")

    out("Not present in this file\n")
    out("------------------------\n")
    for item in [
        "registry snapshot, record version and remote index",
        "MCP protocol version offered or negotiated",
        "call arguments and tool-selection scores",
        "HTTP status codes, error text and failure classes",
        "isError state for calls that returned a result",
        "category labels, descriptions and full tool lists",
        "authentication declarations and tenant context",
        "probe region, observer identity and concurrency",
        "raw responses and repeat observations per endpoint",
    ]:
        out("  - %s\n" % item)
    out("\n")

    out("Done. All figures above are derived only from the verified file.\n")


def main():
    path = data_path()
    try:
        digest = verify_checksum(path)
        rows = read_rows(path)
        facts = analyse(rows)
    except AnalysisError as error:
        sys.stderr.write("analyze_probe.py: %s\n" % error)
        return 1
    report(facts, digest, path)
    return 0


if __name__ == "__main__":
    sys.exit(main())