standard_security_experiments.py

from __future__ import annotations

import argparse
import csv
import json
import math
import os
import tempfile
from dataclasses import asdict, dataclass
from pathlib import Path
from typing import Any, Callable, Iterable

os.environ.setdefault(
    "MPLCONFIGDIR",
    str(Path(tempfile.gettempdir()) / "roi-matplotlib-cache"),
)

import cv2
import matplotlib
import numpy as np
from scipy.special import erfc, gammaincc, ndtr
from scipy.stats import chi2

matplotlib.use("Agg")
import matplotlib.pyplot as plt

from roi_sensitivity_adapter import (
    algorithm_class,
    default_center_roi,
    full_image_roi,
    roi_pixel_count,
)
from sensitivity import (
    adjacent_correlation,
    channel_names,
    image_entropy,
    npcr_uaci,
    split_channels,
)


DEFAULT_KEY = "roi-test-key"
THEORETICAL_NPCR = 100.0 * 255.0 / 256.0
THEORETICAL_UACI = 100.0 * 257.0 / 768.0


@dataclass
class NistResult:
    name: str
    status: str
    p_value: float | None
    details: dict[str, Any]


def json_value(value: Any) -> Any:
    if isinstance(value, dict):
        return {str(key): json_value(item) for key, item in value.items()}
    if isinstance(value, (list, tuple)):
        return [json_value(item) for item in value]
    if isinstance(value, np.ndarray):
        return value.tolist()
    if isinstance(value, np.generic):
        return value.item()
    if isinstance(value, float) and not math.isfinite(value):
        return None
    return value


def load_image(path: str | None) -> tuple[np.ndarray, str]:
    if path is None:
        default_path = Path(__file__).with_name("natural.jpg")
        if default_path.exists():
            path = str(default_path)

    if path is not None:
        image = cv2.imread(path, cv2.IMREAD_UNCHANGED)
        if image is None:
            raise ValueError(f"Cannot read image: {path}")
        source = str(Path(path).resolve())
    else:
        height, width = 512, 512
        y, x = np.mgrid[0:height, 0:width]
        image = np.stack(
            (
                (x * 255 // max(1, width - 1)),
                (y * 255 // max(1, height - 1)),
                ((x + y) * 255 // max(1, width + height - 2)),
            ),
            axis=2,
        ).astype(np.uint8)
        source = "generated 512x512 gradient"

    if image.dtype != np.uint8:
        image = np.clip(image, 0, 255).astype(np.uint8)
    if image.ndim == 3 and image.shape[2] == 4:
        image = cv2.cvtColor(image, cv2.COLOR_BGRA2BGR)
    if image.ndim not in (2, 3):
        raise ValueError("Only 8-bit grayscale or color images are supported")
    return np.ascontiguousarray(image), source


def save_image(path: Path, image: np.ndarray) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    if not cv2.imwrite(str(path), image):
        raise OSError(f"Failed to save image: {path}")


def parse_box(text: str) -> tuple[int, int, int, int]:
    values = [int(round(float(value.strip()))) for value in text.split(",")]
    if len(values) != 4:
        raise argparse.ArgumentTypeError("ROI box must be x1,y1,x2,y2")
    return tuple(values)  # type: ignore[return-value]


def build_rois(
    mode: str,
    image: np.ndarray,
    ratio: float,
    boxes: Iterable[tuple[int, int, int, int]] | None,
) -> list[tuple[int, int, int, int]] | None:
    if mode == "center":
        return default_center_roi(image, ratio)
    if mode == "full":
        return full_image_roi(image)
    if mode == "none":
        return []
    if mode == "custom":
        return list(boxes or [])
    if mode == "yolo":
        return None
    raise ValueError(f"Unsupported ROI mode: {mode}")


def flip_one_plaintext_bit(
    image: np.ndarray,
    rng: np.random.Generator,
) -> tuple[np.ndarray, dict[str, int]]:
    changed = image.copy()
    flat = changed.reshape(-1)
    byte_index = int(rng.integers(0, flat.size))
    bit_index = int(rng.integers(0, 8))
    before = int(flat[byte_index])
    flat[byte_index] ^= np.uint8(1 << bit_index)
    return changed, {
        "byte_index": byte_index,
        "bit_index": bit_index,
        "before": before,
        "after": int(flat[byte_index]),
    }


def one_bit_key_variant(key: str) -> tuple[str, dict[str, Any]]:
    encoded = bytearray(key.encode("utf-8"))
    for byte_index in range(len(encoded)):
        for bit_index in range(8):
            candidate_bytes = encoded.copy()
            candidate_bytes[byte_index] ^= 1 << bit_index
            try:
                candidate = bytes(candidate_bytes).decode("utf-8")
            except UnicodeDecodeError:
                continue
            if candidate == key or candidate.encode("utf-8") != bytes(candidate_bytes):
                continue
            if not all(character.isprintable() for character in candidate):
                continue
            distance = sum(
                (left ^ right).bit_count()
                for left, right in zip(encoded, candidate_bytes)
            )
            if distance == 1:
                return candidate, {
                    "byte_index": byte_index,
                    "bit_index": bit_index,
                    "xor_mask_hex": f"{1 << bit_index:02x}",
                    "base_utf8_length": len(encoded),
                    "variant_utf8_length": len(candidate_bytes),
                    "hamming_distance_bits": distance,
                }
    raise ValueError(
        "Could not create a printable UTF-8 key differing by exactly one bit. "
        "Use a printable ASCII --key."
    )


def metric_summary(values: list[float]) -> dict[str, float]:
    array = np.asarray(values, dtype=np.float64)
    return {
        "mean": float(np.mean(array)),
        "std": float(np.std(array)),
        "min": float(np.min(array)),
        "max": float(np.max(array)),
    }


def differential_experiment(
    encryptor: algorithm_class,
    original: np.ndarray,
    key: str,
    iterations: int,
    seed: int,
) -> tuple[dict[str, Any], list[dict[str, Any]], np.ndarray]:
    base_cipher = encryptor.encrypt(original, key=key)
    rng = np.random.default_rng(seed)
    rows: list[dict[str, Any]] = []
    npcr_values: list[float] = []
    uaci_values: list[float] = []
    first_changed_cipher = base_cipher

    for trial in range(max(1, iterations)):
        changed_plain, change = flip_one_plaintext_bit(original, rng)
        changed_cipher = encryptor.encrypt(changed_plain, key=key)
        if trial == 0:
            first_changed_cipher = changed_cipher
        npcr, uaci = npcr_uaci(base_cipher, changed_cipher)
        npcr_values.append(npcr)
        uaci_values.append(uaci)
        rows.append(
            {
                "trial": trial + 1,
                **change,
                "npcr_percent": npcr,
                "uaci_percent": uaci,
            }
        )

    return (
        {
            "iterations": max(1, iterations),
            "npcr_percent": metric_summary(npcr_values),
            "uaci_percent": metric_summary(uaci_values),
            "theoretical_npcr_percent": THEORETICAL_NPCR,
            "theoretical_uaci_percent": THEORETICAL_UACI,
        },
        rows,
        first_changed_cipher,
    )


def chi_square_uniformity(image: np.ndarray) -> dict[str, dict[str, float]]:
    result: dict[str, dict[str, float]] = {}
    for name, channel in zip(channel_names(image), split_channels(image)):
        counts = np.bincount(channel.reshape(-1), minlength=256).astype(np.float64)
        expected = channel.size / 256.0
        statistic = float(np.sum((counts - expected) ** 2 / expected))
        result[name] = {
            "statistic": statistic,
            "p_value": float(chi2.sf(statistic, df=255)),
        }
    statistics = [item["statistic"] for item in result.values()]
    p_values = [item["p_value"] for item in result.values()]
    result["mean"] = {
        "statistic": float(np.mean(statistics)),
        "p_value": float(np.mean(p_values)),
    }
    return result


def plot_histograms(original: np.ndarray, cipher: np.ndarray, path: Path) -> None:
    names = channel_names(original)
    colors = ["black"] if original.ndim == 2 else ["blue", "green", "red"]
    figure, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
    for axis, image, title in zip(
        axes,
        (original, cipher),
        ("Original image histogram", "Cipher image histogram"),
    ):
        for name, channel, color in zip(names, split_channels(image), colors):
            counts = np.bincount(channel.reshape(-1), minlength=256)
            axis.plot(np.arange(256), counts, color=color, label=name, linewidth=1.0)
        axis.set_title(title)
        axis.set_ylabel("Frequency")
        axis.grid(alpha=0.2)
        axis.legend()
    axes[-1].set_xlabel("Pixel value")
    figure.tight_layout()
    path.parent.mkdir(parents=True, exist_ok=True)
    figure.savefig(path, dpi=180)
    plt.close(figure)


def adjacent_pairs(channel: np.ndarray, direction: str) -> tuple[np.ndarray, np.ndarray]:
    if direction == "horizontal":
        return channel[:, :-1].reshape(-1), channel[:, 1:].reshape(-1)
    if direction == "vertical":
        return channel[:-1, :].reshape(-1), channel[1:, :].reshape(-1)
    if direction == "diagonal":
        return channel[:-1, :-1].reshape(-1), channel[1:, 1:].reshape(-1)
    raise ValueError(f"Unsupported direction: {direction}")


def sample_xy(
    x: np.ndarray,
    y: np.ndarray,
    count: int,
    rng: np.random.Generator,
) -> tuple[np.ndarray, np.ndarray]:
    if count > 0 and x.size > count:
        indexes = rng.choice(x.size, size=count, replace=False)
        return x[indexes], y[indexes]
    return x, y


def plot_correlation_scatters(
    original: np.ndarray,
    cipher: np.ndarray,
    out_dir: Path,
    sample_count: int,
    seed: int,
) -> list[str]:
    directions = ("horizontal", "vertical", "diagonal")
    output_paths: list[str] = []
    rng = np.random.default_rng(seed)
    for channel_name, original_channel, cipher_channel in zip(
        channel_names(original),
        split_channels(original),
        split_channels(cipher),
    ):
        figure, axes = plt.subplots(2, 3, figsize=(12, 8), sharex=True, sharey=True)
        for column, direction in enumerate(directions):
            for row, (channel, title_prefix) in enumerate(
                ((original_channel, "Original"), (cipher_channel, "Cipher"))
            ):
                x, y = adjacent_pairs(channel, direction)
                x, y = sample_xy(x, y, sample_count, rng)
                axes[row, column].scatter(x, y, s=3, alpha=0.25, rasterized=True)
                axes[row, column].set_title(f"{title_prefix} {direction}")
                axes[row, column].set_xlim(0, 255)
                axes[row, column].set_ylim(0, 255)
                axes[row, column].grid(alpha=0.15)
        figure.suptitle(f"Adjacent pixel correlation: {channel_name}")
        figure.supxlabel("Pixel x")
        figure.supylabel("Adjacent pixel y")
        figure.tight_layout()
        output_path = out_dir / f"correlation_scatter_{channel_name}.png"
        figure.savefig(output_path, dpi=180)
        plt.close(figure)
        output_paths.append(str(output_path.resolve()))
    return output_paths


def ensure_bits(bits: np.ndarray) -> np.ndarray:
    array = np.asarray(bits, dtype=np.uint8).reshape(-1)
    if array.size and not np.all((array == 0) | (array == 1)):
        raise ValueError("NIST input must contain only zero and one")
    return array


def nist_result(
    name: str,
    p_values: float | Iterable[float],
    alpha: float,
    **details: Any,
) -> NistResult:
    values = np.atleast_1d(np.asarray(p_values, dtype=np.float64))
    finite = values[np.isfinite(values)]
    if finite.size == 0:
        return NistResult(name, "SKIP", None, details)
    minimum = float(np.min(finite))
    details["p_values"] = [float(value) for value in finite]
    return NistResult(
        name=name,
        status="PASS" if np.all(finite >= alpha) else "FAIL",
        p_value=minimum,
        details=details,
    )


def nist_skip(name: str, reason: str, **details: Any) -> NistResult:
    return NistResult(name, "SKIP", None, {"reason": reason, **details})


def frequency_monobit(bits: np.ndarray, alpha: float) -> NistResult:
    n = bits.size
    if n < 100:
        return nist_skip("Frequency (Monobit)", "requires at least 100 bits", n=n)
    total = int(np.sum(bits.astype(np.int64) * 2 - 1))
    p_value = float(erfc(abs(total) / math.sqrt(2.0 * n)))
    return nist_result("Frequency (Monobit)", p_value, alpha, n=n, sum=total)


def block_frequency(bits: np.ndarray, alpha: float, block_size: int = 128) -> NistResult:
    n = bits.size
    blocks = n // block_size
    if blocks < 20:
        return nist_skip(
            "Frequency within a Block",
            "requires at least 20 complete blocks",
            n=n,
            block_size=block_size,
        )
    used = bits[: blocks * block_size].reshape(blocks, block_size)
    proportions = np.mean(used, axis=1)
    statistic = float(4.0 * block_size * np.sum((proportions - 0.5) ** 2))
    p_value = float(gammaincc(blocks / 2.0, statistic / 2.0))
    return nist_result(
        "Frequency within a Block",
        p_value,
        alpha,
        n=n,
        block_size=block_size,
        blocks=blocks,
        chi_square=statistic,
    )


def runs_test(bits: np.ndarray, alpha: float) -> NistResult:
    n = bits.size
    if n < 100:
        return nist_skip("Runs", "requires at least 100 bits", n=n)
    proportion = float(np.mean(bits))
    prerequisite = 2.0 / math.sqrt(n)
    if abs(proportion - 0.5) >= prerequisite:
        return NistResult(
            "Runs",
            "FAIL",
            0.0,
            {
                "reason": "frequency prerequisite failed",
                "n": n,
                "proportion_ones": proportion,
                "allowed_deviation": prerequisite,
            },
        )
    runs = 1 + int(np.count_nonzero(bits[1:] != bits[:-1]))
    numerator = abs(runs - 2.0 * n * proportion * (1.0 - proportion))
    denominator = 2.0 * math.sqrt(2.0 * n) * proportion * (1.0 - proportion)
    p_value = float(erfc(numerator / denominator))
    return nist_result(
        "Runs",
        p_value,
        alpha,
        n=n,
        runs=runs,
        proportion_ones=proportion,
    )


def longest_run_test(bits: np.ndarray, alpha: float) -> NistResult:
    n = bits.size
    if n < 128:
        return nist_skip("Longest Run of Ones", "requires at least 128 bits", n=n)
    if n < 6272:
        block_size = 8
        probabilities = np.array([0.2148, 0.3672, 0.2305, 0.1875])

        def category(value: int) -> int:
            return min(max(value - 1, 0), 3)

    elif n < 750000:
        block_size = 128
        probabilities = np.array([0.1174, 0.2430, 0.2493, 0.1752, 0.1027, 0.1124])

        def category(value: int) -> int:
            if value <= 4:
                return 0
            if value >= 9:
                return 5
            return value - 4

    else:
        block_size = 10000
        probabilities = np.array(
            [0.0882, 0.2092, 0.2483, 0.1933, 0.1208, 0.0675, 0.0727]
        )

        def category(value: int) -> int:
            if value <= 10:
                return 0
            if value >= 16:
                return 6
            return value - 10

    blocks = n // block_size
    if blocks == 0:
        return nist_skip("Longest Run of Ones", "no complete blocks", n=n)
    observed = np.zeros(probabilities.size, dtype=np.int64)
    for block in bits[: blocks * block_size].reshape(blocks, block_size):
        padded = np.concatenate(([0], block, [0]))
        changes = np.flatnonzero(padded[1:] != padded[:-1])
        longest = 0
        for start, end in zip(changes[::2], changes[1::2]):
            longest = max(longest, int(end - start))
        observed[category(longest)] += 1
    expected = blocks * probabilities
    statistic = float(np.sum((observed - expected) ** 2 / expected))
    p_value = float(gammaincc((probabilities.size - 1) / 2.0, statistic / 2.0))
    return nist_result(
        "Longest Run of Ones",
        p_value,
        alpha,
        n=n,
        block_size=block_size,
        blocks=blocks,
        observed=observed.tolist(),
        chi_square=statistic,
    )


def gf2_rank_32(matrix_bits: np.ndarray) -> int:
    packed = np.packbits(matrix_bits.reshape(32, 32), axis=1, bitorder="big")
    rows = [int.from_bytes(row.tobytes(), "big") for row in packed]
    rank = 0
    for column in range(31, -1, -1):
        pivot = next(
            (index for index in range(rank, 32) if (rows[index] >> column) & 1),
            None,
        )
        if pivot is None:
            continue
        rows[rank], rows[pivot] = rows[pivot], rows[rank]
        for index in range(32):
            if index != rank and ((rows[index] >> column) & 1):
                rows[index] ^= rows[rank]
        rank += 1
        if rank == 32:
            break
    return rank


def rank_probability(rank: int, rows: int = 32, columns: int = 32) -> float:
    product = 1.0
    for index in range(rank):
        numerator = (1.0 - 2.0 ** (index - rows)) * (
            1.0 - 2.0 ** (index - columns)
        )
        denominator = 1.0 - 2.0 ** (index - rank)
        product *= numerator / denominator
    return 2.0 ** (rank * (rows + columns - rank) - rows * columns) * product


def binary_matrix_rank_test(bits: np.ndarray, alpha: float) -> NistResult:
    matrix_bits = 32 * 32
    blocks = bits.size // matrix_bits
    if blocks < 38:
        return nist_skip(
            "Binary Matrix Rank",
            "requires at least 38 complete 32x32 matrices",
            n=bits.size,
            matrices=blocks,
        )
    frequencies = np.zeros(3, dtype=np.int64)
    used = bits[: blocks * matrix_bits].reshape(blocks, matrix_bits)
    for block in used:
        rank = gf2_rank_32(block)
        if rank == 32:
            frequencies[0] += 1
        elif rank == 31:
            frequencies[1] += 1
        else:
            frequencies[2] += 1
    p32 = rank_probability(32)
    p31 = rank_probability(31)
    probabilities = np.array([p32, p31, 1.0 - p32 - p31])
    expected = blocks * probabilities
    statistic = float(np.sum((frequencies - expected) ** 2 / expected))
    p_value = float(math.exp(-statistic / 2.0))
    return nist_result(
        "Binary Matrix Rank",
        p_value,
        alpha,
        n=bits.size,
        matrices=blocks,
        observed=frequencies.tolist(),
        probabilities=probabilities.tolist(),
        chi_square=statistic,
    )


def discrete_fourier_transform_test(bits: np.ndarray, alpha: float) -> NistResult:
    n = bits.size
    if n < 1000:
        return nist_skip("Discrete Fourier Transform", "requires at least 1000 bits", n=n)
    sequence = bits.astype(np.float64) * 2.0 - 1.0
    magnitudes = np.abs(np.fft.fft(sequence))[: n // 2]
    threshold = math.sqrt(math.log(20.0) * n)
    expected_below = 0.95 * n / 2.0
    observed_below = int(np.count_nonzero(magnitudes < threshold))
    deviation = (observed_below - expected_below) / math.sqrt(
        n * 0.95 * 0.05 / 4.0
    )
    p_value = float(erfc(abs(deviation) / math.sqrt(2.0)))
    return nist_result(
        "Discrete Fourier Transform",
        p_value,
        alpha,
        n=n,
        threshold=threshold,
        expected_below=expected_below,
        observed_below=observed_below,
        normalized_deviation=deviation,
    )


def non_overlapping_template_test(
    bits: np.ndarray,
    alpha: float,
    template: np.ndarray | None = None,
) -> NistResult:
    if template is None:
        template = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1], dtype=np.uint8)
    template = ensure_bits(template)
    m = template.size
    block_count = 8
    block_size = bits.size // block_count
    if block_size < 1000 or block_size < m:
        return nist_skip(
            "Non-overlapping Template Matching",
            "requires eight sufficiently large blocks",
            n=bits.size,
            block_size=block_size,
        )
    counts = np.zeros(block_count, dtype=np.int64)
    for block_index in range(block_count):
        block = bits[block_index * block_size : (block_index + 1) * block_size]
        windows = np.lib.stride_tricks.sliding_window_view(block, m)
        matches = np.flatnonzero(np.all(windows == template, axis=1))
        next_allowed = 0
        count = 0
        for match in matches:
            if int(match) >= next_allowed:
                count += 1
                next_allowed = int(match) + m
        counts[block_index] = count
    mean = (block_size - m + 1) / (2.0**m)
    variance = block_size * (
        1.0 / (2.0**m) - (2.0 * m - 1.0) / (2.0 ** (2 * m))
    )
    statistic = float(np.sum((counts - mean) ** 2 / variance))
    p_value = float(gammaincc(block_count / 2.0, statistic / 2.0))
    return nist_result(
        "Non-overlapping Template Matching",
        p_value,
        alpha,
        n=bits.size,
        template="".join(str(int(value)) for value in template),
        block_size=block_size,
        observed=counts.tolist(),
        mean=mean,
        variance=variance,
        chi_square=statistic,
        note="One representative aperiodic template is tested; official STS uses 148 templates.",
    )


def overlapping_template_test(bits: np.ndarray, alpha: float) -> NistResult:
    template_size = 9
    block_size = 1032
    blocks = bits.size // block_size
    if blocks < 20:
        return nist_skip(
            "Overlapping Template Matching",
            "requires at least 20 complete 1032-bit blocks",
            n=bits.size,
            blocks=blocks,
        )
    used = bits[: blocks * block_size].reshape(blocks, block_size)
    windows = np.lib.stride_tricks.sliding_window_view(
        used,
        template_size,
        axis=1,
    )
    match_counts = np.count_nonzero(np.all(windows == 1, axis=2), axis=1)
    observed = np.bincount(np.minimum(match_counts, 5), minlength=6)
    probabilities = np.array(
        [0.364091, 0.185659, 0.139381, 0.100571, 0.070432, 0.139865]
    )
    expected = blocks * probabilities
    statistic = float(np.sum((observed - expected) ** 2 / expected))
    p_value = float(gammaincc(5.0 / 2.0, statistic / 2.0))
    return nist_result(
        "Overlapping Template Matching",
        p_value,
        alpha,
        n=bits.size,
        block_size=block_size,
        blocks=blocks,
        observed=observed.tolist(),
        chi_square=statistic,
    )


UNIVERSAL_PARAMETERS = (
    (1059061760, 16, 655360, 15.167379, 3.421),
    (496435200, 15, 327680, 14.167488, 3.419),
    (231669760, 14, 163840, 13.167693, 3.416),
    (107560960, 13, 81920, 12.168070, 3.410),
    (49643520, 12, 40960, 11.168765, 3.401),
    (22753280, 11, 20480, 10.170032, 3.384),
    (10342400, 10, 10240, 9.172324, 3.356),
    (4654080, 9, 5120, 8.176424, 3.311),
    (2068480, 8, 2560, 7.183666, 3.238),
    (904960, 7, 1280, 6.196251, 3.125),
    (387840, 6, 640, 5.217705, 2.954),
)


def universal_statistical_test(bits: np.ndarray, alpha: float) -> NistResult:
    n = bits.size
    parameters = next(
        (values for values in UNIVERSAL_PARAMETERS if n >= values[0]),
        None,
    )
    if parameters is None:
        return nist_skip(
            "Maurer's Universal Statistical",
            "requires at least 387840 bits",
            n=n,
        )
    _, block_length, initialization_blocks, expected, variance = parameters
    total_blocks = n // block_length
    test_blocks = total_blocks - initialization_blocks
    if test_blocks <= 0:
        return nist_skip(
            "Maurer's Universal Statistical",
            "not enough test blocks after initialization",
            n=n,
        )
    used = bits[: total_blocks * block_length].reshape(total_blocks, block_length)
    weights = 1 << np.arange(block_length - 1, -1, -1, dtype=np.int64)
    values = used @ weights
    table = np.zeros(1 << block_length, dtype=np.int64)
    for index in range(initialization_blocks):
        table[int(values[index])] = index + 1
    total = 0.0
    for index in range(initialization_blocks, total_blocks):
        value = int(values[index])
        position = index + 1
        total += math.log2(position - int(table[value]))
        table[value] = position
    statistic = total / test_blocks
    correction = (
        0.7
        - 0.8 / block_length
        + (4.0 + 32.0 / block_length)
        * (test_blocks ** (-3.0 / block_length))
        / 15.0
    )
    sigma = correction * math.sqrt(variance / test_blocks)
    p_value = float(erfc(abs(statistic - expected) / (math.sqrt(2.0) * sigma)))
    return nist_result(
        "Maurer's Universal Statistical",
        p_value,
        alpha,
        n=n,
        block_length=block_length,
        initialization_blocks=initialization_blocks,
        test_blocks=test_blocks,
        statistic=statistic,
        expected=expected,
        sigma=sigma,
    )


def berlekamp_massey_length(block: np.ndarray) -> int:
    connection = 1
    previous = 1
    complexity = 0
    last_update = -1
    reversed_prefix = 0
    for index, value in enumerate(block):
        reversed_prefix = (reversed_prefix << 1) | int(value)
        discrepancy = (connection & reversed_prefix).bit_count() & 1
        if discrepancy:
            old_connection = connection
            connection ^= previous << (index - last_update)
            if 2 * complexity <= index:
                complexity = index + 1 - complexity
                previous = old_connection
                last_update = index
    return complexity


def linear_complexity_test(
    bits: np.ndarray,
    alpha: float,
    block_size: int = 500,
) -> NistResult:
    blocks = bits.size // block_size
    if blocks < 200:
        return nist_skip(
            "Linear Complexity",
            "requires at least 200 complete 500-bit blocks",
            n=bits.size,
            blocks=blocks,
        )
    probabilities = np.array(
        [0.01047, 0.03125, 0.12500, 0.50000, 0.25000, 0.06250, 0.020833]
    )
    observed = np.zeros(7, dtype=np.int64)
    mean = (
        block_size / 2.0
        + (9.0 + (-1.0) ** (block_size + 1)) / 36.0
        - (block_size / 3.0 + 2.0 / 9.0) / (2.0**block_size)
    )
    used = bits[: blocks * block_size].reshape(blocks, block_size)
    for block in used:
        complexity = berlekamp_massey_length(block)
        transformed = ((-1.0) ** block_size) * (complexity - mean) + 2.0 / 9.0
        if transformed <= -2.5:
            category = 0
        elif transformed <= -1.5:
            category = 1
        elif transformed <= -0.5:
            category = 2
        elif transformed <= 0.5:
            category = 3
        elif transformed <= 1.5:
            category = 4
        elif transformed <= 2.5:
            category = 5
        else:
            category = 6
        observed[category] += 1
    expected = blocks * probabilities
    statistic = float(np.sum((observed - expected) ** 2 / expected))
    p_value = float(gammaincc(3.0, statistic / 2.0))
    return nist_result(
        "Linear Complexity",
        p_value,
        alpha,
        n=bits.size,
        block_size=block_size,
        blocks=blocks,
        observed=observed.tolist(),
        chi_square=statistic,
    )


def cyclic_pattern_counts(bits: np.ndarray, pattern_size: int) -> np.ndarray:
    extended = np.concatenate((bits, bits[: pattern_size - 1]))
    windows = np.lib.stride_tricks.sliding_window_view(extended, pattern_size)
    weights = 1 << np.arange(pattern_size - 1, -1, -1, dtype=np.int64)
    values = windows @ weights
    return np.bincount(values, minlength=1 << pattern_size)


def psi_square(bits: np.ndarray, pattern_size: int) -> float:
    if pattern_size <= 0:
        return 0.0
    counts = cyclic_pattern_counts(bits, pattern_size).astype(np.float64)
    return float((2.0**pattern_size / bits.size) * np.sum(counts**2) - bits.size)


def serial_test(bits: np.ndarray, alpha: float, pattern_size: int = 16) -> NistResult:
    n = bits.size
    if n < 1000:
        return nist_skip("Serial", "requires at least 1000 bits", n=n)
    maximum = max(3, int(math.floor(math.log2(n))) - 2)
    pattern_size = min(pattern_size, maximum)
    psi_m = psi_square(bits, pattern_size)
    psi_m1 = psi_square(bits, pattern_size - 1)
    psi_m2 = psi_square(bits, pattern_size - 2)
    delta1 = psi_m - psi_m1
    delta2 = psi_m - 2.0 * psi_m1 + psi_m2
    p1 = float(gammaincc(2.0 ** (pattern_size - 2), delta1 / 2.0))
    p2 = float(gammaincc(2.0 ** (pattern_size - 3), delta2 / 2.0))
    return nist_result(
        "Serial",
        [p1, p2],
        alpha,
        n=n,
        pattern_size=pattern_size,
        delta1=delta1,
        delta2=delta2,
    )


def approximate_entropy_test(
    bits: np.ndarray,
    alpha: float,
    pattern_size: int = 10,
) -> NistResult:
    n = bits.size
    if n < 1000:
        return nist_skip("Approximate Entropy", "requires at least 1000 bits", n=n)
    maximum = max(2, int(math.floor(math.log2(n))) - 5)
    pattern_size = min(pattern_size, maximum)

    def phi(size: int) -> float:
        counts = cyclic_pattern_counts(bits, size).astype(np.float64)
        probabilities = counts[counts > 0] / n
        return float(np.sum(probabilities * np.log(probabilities)))

    approximate_entropy = phi(pattern_size) - phi(pattern_size + 1)
    statistic = float(2.0 * n * (math.log(2.0) - approximate_entropy))
    p_value = float(gammaincc(2.0 ** (pattern_size - 1), statistic / 2.0))
    return nist_result(
        "Approximate Entropy",
        p_value,
        alpha,
        n=n,
        pattern_size=pattern_size,
        approximate_entropy=approximate_entropy,
        chi_square=statistic,
    )


def cumulative_sums_p_value(sequence: np.ndarray) -> tuple[float, int]:
    n = sequence.size
    walk = np.cumsum(sequence.astype(np.int64) * 2 - 1)
    maximum = int(np.max(np.abs(walk)))
    if maximum == 0:
        return 1.0, maximum
    root_n = math.sqrt(n)
    first_start = math.floor((-n / maximum + 1.0) / 4.0)
    first_end = math.floor((n / maximum - 1.0) / 4.0)
    first_sum = sum(
        ndtr((4 * index + 1) * maximum / root_n)
        - ndtr((4 * index - 1) * maximum / root_n)
        for index in range(first_start, first_end + 1)
    )
    second_start = math.floor((-n / maximum - 3.0) / 4.0)
    second_end = math.floor((n / maximum - 1.0) / 4.0)
    second_sum = sum(
        ndtr((4 * index + 3) * maximum / root_n)
        - ndtr((4 * index + 1) * maximum / root_n)
        for index in range(second_start, second_end + 1)
    )
    return float(np.clip(1.0 - first_sum + second_sum, 0.0, 1.0)), maximum


def cumulative_sums_test(bits: np.ndarray, alpha: float) -> NistResult:
    if bits.size < 100:
        return nist_skip("Cumulative Sums", "requires at least 100 bits", n=bits.size)
    forward_p, forward_maximum = cumulative_sums_p_value(bits)
    reverse_p, reverse_maximum = cumulative_sums_p_value(bits[::-1])
    return nist_result(
        "Cumulative Sums",
        [forward_p, reverse_p],
        alpha,
        n=bits.size,
        forward_maximum=forward_maximum,
        reverse_maximum=reverse_maximum,
    )


def excursion_cycle_data(bits: np.ndarray) -> tuple[np.ndarray, np.ndarray, int]:
    walk = np.concatenate(
        (
            np.array([0], dtype=np.int64),
            np.cumsum(bits.astype(np.int64) * 2 - 1),
            np.array([0], dtype=np.int64),
        )
    )
    zero_indexes = np.flatnonzero(walk == 0)
    cycles = zero_indexes.size - 1
    positions = np.arange(walk.size)
    cycle_ids = np.searchsorted(zero_indexes, positions, side="right") - 1
    valid = (cycle_ids >= 0) & (cycle_ids < cycles)
    return walk[valid], cycle_ids[valid], cycles


def random_excursions_test(bits: np.ndarray, alpha: float) -> NistResult:
    walk, cycle_ids, cycles = excursion_cycle_data(bits)
    if cycles < 500:
        return nist_skip(
            "Random Excursions",
            "NIST requires at least 500 zero-return cycles",
            n=bits.size,
            cycles=cycles,
        )
    p_values: list[float] = []
    state_details: dict[str, Any] = {}
    for state in (-4, -3, -2, -1, 1, 2, 3, 4):
        visits = np.bincount(
            cycle_ids[walk == state],
            minlength=cycles,
        )
        observed = np.bincount(np.minimum(visits, 5), minlength=6)
        absolute = abs(state)
        base = 1.0 - 1.0 / (2.0 * absolute)
        probabilities = np.empty(6, dtype=np.float64)
        probabilities[0] = base
        for count in range(1, 5):
            probabilities[count] = (base ** (count - 1)) / (4.0 * absolute**2)
        probabilities[5] = (base**4) / (2.0 * absolute)
        expected = cycles * probabilities
        statistic = float(np.sum((observed - expected) ** 2 / expected))
        p_value = float(gammaincc(5.0 / 2.0, statistic / 2.0))
        p_values.append(p_value)
        state_details[str(state)] = {
            "p_value": p_value,
            "observed": observed.tolist(),
            "chi_square": statistic,
        }
    return nist_result(
        "Random Excursions",
        p_values,
        alpha,
        n=bits.size,
        cycles=cycles,
        states=state_details,
    )


def random_excursions_variant_test(bits: np.ndarray, alpha: float) -> NistResult:
    walk, _, cycles = excursion_cycle_data(bits)
    if cycles < 500:
        return nist_skip(
            "Random Excursions Variant",
            "NIST requires at least 500 zero-return cycles",
            n=bits.size,
            cycles=cycles,
        )
    p_values: list[float] = []
    state_details: dict[str, Any] = {}
    for state in tuple(range(-9, 0)) + tuple(range(1, 10)):
        visits = int(np.count_nonzero(walk == state))
        denominator = math.sqrt(2.0 * cycles * (4.0 * abs(state) - 2.0))
        p_value = float(erfc(abs(visits - cycles) / denominator))
        p_values.append(p_value)
        state_details[str(state)] = {
            "p_value": p_value,
            "visits": visits,
        }
    return nist_result(
        "Random Excursions Variant",
        p_values,
        alpha,
        n=bits.size,
        cycles=cycles,
        states=state_details,
    )


def run_nist_suite(bits: np.ndarray, alpha: float) -> list[NistResult]:
    bits = ensure_bits(bits)
    tests: list[Callable[[np.ndarray, float], NistResult]] = [
        frequency_monobit,
        block_frequency,
        runs_test,
        longest_run_test,
        binary_matrix_rank_test,
        discrete_fourier_transform_test,
        non_overlapping_template_test,
        overlapping_template_test,
        universal_statistical_test,
        linear_complexity_test,
        serial_test,
        approximate_entropy_test,
        cumulative_sums_test,
        random_excursions_test,
        random_excursions_variant_test,
    ]
    results: list[NistResult] = []
    for test in tests:
        try:
            results.append(test(bits, alpha))
        except Exception as exc:
            results.append(
                nist_skip(
                    test.__name__,
                    f"implementation error: {type(exc).__name__}: {exc}",
                )
            )
    return results


def write_csv(path: Path, rows: list[dict[str, Any]]) -> None:
    if not rows:
        return
    path.parent.mkdir(parents=True, exist_ok=True)
    fieldnames: list[str] = []
    for row in rows:
        for key in row:
            if key not in fieldnames:
                fieldnames.append(key)
    with path.open("w", newline="", encoding="utf-8") as file:
        writer = csv.DictWriter(file, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(rows)


def format_channel_metrics(
    title: str,
    values: dict[str, Any],
    value_name: str | None = None,
) -> list[str]:
    lines = [title]
    for channel, item in values.items():
        if isinstance(item, dict):
            if value_name is None:
                rendered = ", ".join(
                    f"{name}={float(value):.6f}" for name, value in item.items()
                )
            else:
                rendered = f"{value_name}={float(item[value_name]):.6f}"
        else:
            rendered = f"{float(item):.6f}"
        lines.append(f"  {channel}: {rendered}")
    return lines


def main() -> None:
    parser = argparse.ArgumentParser(
        description=(
            "Complete standard security experiments for the roi2.py encryption "
            "pipeline without modifying the algorithm source."
        )
    )
    parser.add_argument("--image", help="Input image. Defaults to ROI/natural.jpg.")
    parser.add_argument("--key", default=DEFAULT_KEY, help="Base password.")
    parser.add_argument(
        "--algorithm-source",
        default=None,
        help="Algorithm source file/module. Auto-selection resolves to roi2.py.",
    )
    parser.add_argument(
        "--roi-mode",
        choices=["center", "full", "none", "custom", "yolo"],
        default="center",
    )
    parser.add_argument("--roi-ratio", type=float, default=0.25)
    parser.add_argument("--roi-box", action="append", type=parse_box)
    parser.add_argument("--model", default=None, help="YOLO model path.")
    parser.add_argument("--confidence", type=float, default=0.25)
    parser.add_argument("--iterations", type=int, default=20)
    parser.add_argument("--correlation-samples", type=int, default=10000)
    parser.add_argument("--scatter-samples", type=int, default=3000)
    parser.add_argument("--alpha", type=float, default=0.01)
    parser.add_argument(
        "--nist-max-bits",
        type=int,
        default=1_000_000,
        help="Maximum contiguous cipher bits used by the built-in NIST screening.",
    )
    parser.add_argument("--skip-nist", action="store_true")
    parser.add_argument("--skip-plots", action="store_true")
    parser.add_argument(
        "--out-dir",
        default="standard_security_outputs",
    )
    args = parser.parse_args()

    out_dir = Path(args.out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    original, image_source = load_image(args.image)
    fixed_rois = build_rois(
        args.roi_mode,
        original,
        args.roi_ratio,
        args.roi_box,
    )
    use_yolo = args.roi_mode == "yolo"

    production_encryptor = algorithm_class(
        algorithm_source=args.algorithm_source,
        use_image_hash=True,
        use_yolo=use_yolo,
        model_path=args.model,
        confidence=args.confidence,
        roi_boxes=fixed_rois,
    )
    fixed_material_encryptor = algorithm_class(
        algorithm_source=args.algorithm_source,
        use_image_hash=False,
        use_yolo=use_yolo,
        model_path=args.model,
        confidence=args.confidence,
        roi_boxes=fixed_rois,
    )

    cipher, stage_timings = production_encryptor.encrypt_timed(
        original,
        key=args.key,
    )
    algorithm_path = str(Path(production_encryptor.alg.__file__).resolve())
    save_image(out_dir / "original.png", original)
    save_image(out_dir / "cipher.png", cipher)

    original_entropy = image_entropy(original)
    cipher_entropy = image_entropy(cipher)
    original_chi_square = chi_square_uniformity(original)
    cipher_chi_square = chi_square_uniformity(cipher)
    correlation_rng = np.random.default_rng(20260614)
    original_correlation = adjacent_correlation(
        original,
        args.correlation_samples,
        correlation_rng,
    )
    cipher_correlation = adjacent_correlation(
        cipher,
        args.correlation_samples,
        correlation_rng,
    )

    key_variant, key_change = one_bit_key_variant(args.key)
    key_cipher = production_encryptor.encrypt(original, key=key_variant)
    key_npcr, key_uaci = npcr_uaci(cipher, key_cipher)
    save_image(out_dir / "cipher_key_one_bit_changed.png", key_cipher)

    production_differential, production_rows, production_changed_cipher = (
        differential_experiment(
            production_encryptor,
            original,
            args.key,
            args.iterations,
            seed=20260615,
        )
    )
    fixed_differential, fixed_rows, fixed_changed_cipher = differential_experiment(
        fixed_material_encryptor,
        original,
        args.key,
        args.iterations,
        seed=20260615,
    )
    save_image(
        out_dir / "cipher_plaintext_one_bit_changed_production.png",
        production_changed_cipher,
    )
    save_image(
        out_dir / "cipher_plaintext_one_bit_changed_fixed_material.png",
        fixed_changed_cipher,
    )
    for row in production_rows:
        row["mode"] = "production_image_bound_kdf"
    for row in fixed_rows:
        row["mode"] = "fixed_derivation_material"
    write_csv(
        out_dir / "differential_trials.csv",
        production_rows + fixed_rows,
    )

    plot_paths: list[str] = []
    if not args.skip_plots:
        histogram_path = out_dir / "histograms.png"
        plot_histograms(original, cipher, histogram_path)
        plot_paths.append(str(histogram_path.resolve()))
        plot_paths.extend(
            plot_correlation_scatters(
                original,
                cipher,
                out_dir,
                args.scatter_samples,
                seed=20260616,
            )
        )

    nist_results: list[NistResult] = []
    nist_bits_used = 0
    if not args.skip_nist:
        cipher_bits = np.unpackbits(
            np.ascontiguousarray(cipher).reshape(-1),
            bitorder="big",
        )
        if args.nist_max_bits > 0:
            cipher_bits = cipher_bits[: args.nist_max_bits]
        nist_bits_used = int(cipher_bits.size)
        nist_results = run_nist_suite(cipher_bits, args.alpha)
        write_csv(
            out_dir / "nist_results.csv",
            [
                {
                    "test": result.name,
                    "status": result.status,
                    "minimum_p_value": result.p_value,
                    "details_json": json.dumps(
                        json_value(result.details),
                        ensure_ascii=False,
                        sort_keys=True,
                    ),
                }
                for result in nist_results
            ],
        )

    if fixed_rois is None:
        roi_pixels: int | None = None
        roi_ratio: float | None = None
        roi_description: Any = "dynamic YOLO boxes"
    else:
        roi_pixels = roi_pixel_count(original.shape, fixed_rois)
        roi_ratio = roi_pixels / (original.shape[0] * original.shape[1])
        roi_description = fixed_rois

    report = {
        "image": {
            "source": image_source,
            "shape": list(original.shape),
            "dtype": str(original.dtype),
        },
        "algorithm": {
            "source": algorithm_path,
            "pipeline": (
                "AES-CTR key streams + recursive ROI Ring-Feistel-like mixing + "
                "global row/column permutation + global XOR diffusion"
            ),
            "roi_mode": args.roi_mode,
            "roi_boxes": roi_description,
            "roi_pixels": roi_pixels,
            "roi_ratio": roi_ratio,
            "stage_timings_seconds": stage_timings,
        },
        "histogram": {
            "plot": str((out_dir / "histograms.png").resolve())
            if not args.skip_plots
            else None,
        },
        "information_entropy": {
            "original": original_entropy,
            "cipher": cipher_entropy,
        },
        "chi_square_uniformity": {
            "degrees_of_freedom": 255,
            "alpha_reference": 0.05,
            "critical_value_reference": float(chi2.ppf(0.95, df=255)),
            "original": original_chi_square,
            "cipher": cipher_chi_square,
        },
        "adjacent_pixel_correlation": {
            "sample_count_per_direction_channel": args.correlation_samples,
            "original": original_correlation,
            "cipher": cipher_correlation,
            "plots": plot_paths[1:] if plot_paths else [],
        },
        "key_sensitivity_one_bit": {
            **key_change,
            "npcr_percent": key_npcr,
            "uaci_percent": key_uaci,
        },
        "plaintext_sensitivity_one_bit": {
            "production_image_bound_kdf": production_differential,
            "fixed_derivation_material": fixed_differential,
            "interpretation": (
                "The production mode changes image-bound KDF material when one "
                "plaintext bit changes. The fixed-material mode isolates diffusion "
                "inside the encryption pipeline."
            ),
        },
        "nist_sp800_22_single_sequence_screening": {
            "enabled": not args.skip_nist,
            "alpha": args.alpha,
            "cipher_bits_used": nist_bits_used,
            "results": [asdict(result) for result in nist_results],
            "limitation": (
                "This built-in suite is a single-sequence screening. Formal NIST "
                "reporting also requires many independent sequences, pass "
                "proportions, p-value uniformity, and validation with official STS."
            ),
        },
        "theoretical_random_cipher_reference": {
            "npcr_percent": THEORETICAL_NPCR,
            "uaci_percent": THEORETICAL_UACI,
        },
    }
    report_path = out_dir / "security_report.json"
    report_path.write_text(
        json.dumps(json_value(report), ensure_ascii=False, indent=2),
        encoding="utf-8",
    )

    report_lines = [
        "ROI standard security experiments",
        f"Image: {image_source}",
        f"Shape: {original.shape}",
        f"Algorithm source: {algorithm_path}",
        f"ROI mode: {args.roi_mode}",
        f"ROI boxes: {roi_description}",
        f"ROI ratio: {roi_ratio if roi_ratio is not None else 'dynamic'}",
        "",
        *format_channel_metrics("Original entropy:", original_entropy),
        *format_channel_metrics("Cipher entropy:", cipher_entropy),
        "",
        *format_channel_metrics("Original chi-square:", original_chi_square),
        *format_channel_metrics("Cipher chi-square:", cipher_chi_square),
        "",
        "Key sensitivity (exactly one UTF-8 bit changed):",
        f"  changed byte index: {key_change['byte_index']}",
        f"  changed bit index: {key_change['bit_index']}",
        f"  XOR mask: 0x{key_change['xor_mask_hex']}",
        f"  Hamming distance: {key_change['hamming_distance_bits']} bit",
        f"  NPCR: {key_npcr:.6f}%",
        f"  UACI: {key_uaci:.6f}%",
        "",
        "Plaintext sensitivity (one plaintext bit changed):",
        (
            "  production image-bound KDF NPCR/UACI: "
            f"{production_differential['npcr_percent']['mean']:.6f}% / "
            f"{production_differential['uaci_percent']['mean']:.6f}%"
        ),
        (
            "  fixed derivation material NPCR/UACI: "
            f"{fixed_differential['npcr_percent']['mean']:.6f}% / "
            f"{fixed_differential['uaci_percent']['mean']:.6f}%"
        ),
        (
            "  random-cipher reference NPCR/UACI: "
            f"{THEORETICAL_NPCR:.6f}% / {THEORETICAL_UACI:.6f}%"
        ),
        "",
        "NIST SP 800-22 single-sequence screening:",
    ]
    if args.skip_nist:
        report_lines.append("  skipped")
    else:
        report_lines.extend(
            f"  {result.name}: {result.status}, minimum p-value="
            f"{result.p_value if result.p_value is not None else 'N/A'}"
            for result in nist_results
        )
        report_lines.append(
            "  Formal paper results must be cross-checked with official NIST STS "
            "using multiple independent sequences."
        )
    report_lines.extend(
        (
            "",
            f"JSON report: {report_path.resolve()}",
            f"Differential trials: {(out_dir / 'differential_trials.csv').resolve()}",
        )
    )
    if not args.skip_nist:
        report_lines.append(
            f"NIST CSV: {(out_dir / 'nist_results.csv').resolve()}"
        )
    text_path = out_dir / "security_report.txt"
    text_path.write_text("\n".join(report_lines) + "\n", encoding="utf-8")
    print("\n".join(report_lines))


if __name__ == "__main__":
    main()