util.py

import math
import numpy as np
import krippendorff
from typing import List, Iterable, Union
from collections import Counter
from statistics import median, mean
from scipy.stats import spearmanr, pearsonr, kendalltau


def review_pmf_exact(score_pmfs: List[np.ndarray],
                     weight_pmfs: List[np.ndarray],
                     cutpoints=(1.5, 2.5, 3.5, 4.5)) -> np.ndarray:
    """
    Return a 5-class probability vector for a single review.

    score_pmfs  : list of k arrays length-5   (Pr[S_i = 1..5])
    weight_pmfs : list of k arrays length-3   (Pr[W_i = 0,1,2])
    """
    state = Counter({(0, 0): 1.0})                  # (score_sum , weight_sum) → prob

    for p_s, p_w in zip(score_pmfs, weight_pmfs):
        nxt = Counter()
        for (s_sum, w_sum), p_accum in state.items():
            for s, ps in enumerate(p_s, start=1):
                if ps == 0:
                    continue
                for w, pw in enumerate(p_w):        # w ∈ {0,1,2}
                    if pw == 0:
                        continue
                    nxt[(s_sum + w * s, w_sum + w)] += p_accum * ps * pw
        state = nxt

    bins = np.array([1, *cutpoints, 5.01])
    pmf = np.zeros(5)
    for (s_sum, w_sum), p in state.items():
        r = (s_sum / w_sum) if w_sum else 3.0       # fallback if all weights = 0
        y = np.digitize(r, bins)                    # 1-based index
        pmf[y - 1] += p

    return pmf / pmf.sum()                          # numerical hygiene


def brier_multi(pmf: np.ndarray, true_label: int) -> float:
    """
    Multi-class Brier score for one example.

    pmf         : length-5 probability vector
    true_label  : integer 1..5
    """
    y_onehot = np.zeros(5)
    y_onehot[true_label - 1] = 1
    return np.sum((pmf - y_onehot) ** 2)


def ece_topclass(pmfs: List[np.ndarray],
                 labels: List[int],
                 num_bins: int = 10) -> float:
    """
    Classic ECE on the *top-class* confidence.

    pmfs   : list of N length-5 probability vectors
    labels : list of N integer gold labels (1..5)
    """
    confidences = np.array([p.max() for p in pmfs])
    preds       = np.array([p.argmax() + 1 for p in pmfs])
    labels      = np.array(labels)

    bin_edges = np.linspace(0, 1, num_bins + 1)
    ece = 0.0
    for lo, hi in zip(bin_edges[:-1], bin_edges[1:]):
        mask = (confidences >= lo) & (confidences < hi)
        if mask.any():
            acc  = (preds[mask] == labels[mask]).mean()
            conf = confidences[mask].mean()
            ece += mask.mean() * abs(acc - conf)
    return ece


def print_grounding_acc(ground_gt_list, ground_list):
    assert len(ground_gt_list) == len(ground_list)
    cnt = 0
    for tmp_gt, tmp in zip(ground_gt_list, ground_list):
        if tmp_gt == tmp:
            cnt += 1
    acc_ = cnt / len(ground_gt_list)
    print(f'\n**Grounding Acc**\n{acc_}')


def print_statistic(score_name, gt_list, score_list, N=3):
    assert len(gt_list) == len(score_list)
    print(f'\n**{score_name} Alignment**')
    # print(f'spearman: {round(spearmanr(score_list, gt_list).statistic, 3)}')
    # print(f'kendall-tau: {round(kendalltau(score_list, gt_list).statistic, 3)}')
    try:
        print(f'pearson: {round(pearsonr(score_list, gt_list).statistic, 3)}')
    except:
        print('pearson: nan')


def calc_krippendorff(list_A, list_B, list_C, N=3):
    # ------------------------------------------------------------
    # Compute Krippendorff’s α for interval-scale data
    #    (use 'ordinal' (rank-only) or 'nominal' if that fits your rubric better)
    # ------------------------------------------------------------
    ratings = np.array([list_A, list_B, list_C])
    try:
        alpha_interval = krippendorff.alpha(ratings, level_of_measurement='ordinal')       # interval
        alpha_interval = round(alpha_interval, N)
    except:
        alpha_interval = np.nan
    
    return alpha_interval

def calc_krippendorff_two(list_A, list_B):
    # ------------------------------------------------------------
    # Compute Krippendorff’s α for interval-scale data
    #    (use 'ordinal' (rank-only) or 'nominal' if that fits your rubric better)
    # ------------------------------------------------------------
    ratings = np.array([list_A, list_B])
    try:
        alpha_interval = krippendorff.alpha(ratings, level_of_measurement='ordinal')       # interval
    except:
        alpha_interval = np.nan
    
    return alpha_interval

def calc_avg_spearman(list_A, list_B, list_C):
    try:
        ab = spearmanr(list_A, list_B).statistic
        bc = spearmanr(list_B, list_C).statistic
        ca = spearmanr(list_C, list_A).statistic
        avg = (ab+bc+ca)/3
    except:
        avg = np.nan
    return avg

def weighted_average(scores, importance):
    """
    Compute the weighted average of scores given corresponding importances.
    
    Parameters
    ----------
    scores : list[int]
        List of integer scores from 1 to 5.
    importance : list[int]
        List of integer weights from 0 to 2, same length as `scores`.
    
    Returns
    -------
    float
        Weighted average score. If all importances are zero, returns 0.0.
    """
    if len(scores) != len(importance):
        raise ValueError("scores and importance must have the same length")
    
    if np.nan in scores or np.nan in importance:
        return np.nan

    total_weight = sum(importance)
    if total_weight == 0:
        return mean(scores)     # changed 0609; original: 0.0

    return sum(s * w for s, w in zip(scores, importance)) / total_weight

def nanmean(values: Iterable[Union[int, float, np.floating, np.integer]]) -> float:
    """
    Return the arithmetic mean of `values`, skipping any np.nan entries.

    Parameters
    ----------
    values : iterable of int | float
        A 1-D container (list, tuple, NumPy array, etc.).  Elements may be
        regular numbers or `np.nan`.

    Returns
    -------
    float
        The average of the non-NaN elements.  If *all* elements are NaN,
        the function returns `float('nan')`.

    Examples
    --------
    >>> nanmean([1, 2, 3])
    2.0
    >>> nanmean([1.0, np.nan, 3.0])
    2.0
    >>> nanmean([np.nan, np.nan])
    nan
    """
    total = 0.0
    count = 0

    for v in values:
        # math.isnan handles both Python floats and NumPy scalars
        if v is not None and not (isinstance(v, float) and math.isnan(v)):
            total += float(v)
            count += 1

    return total / count if count else float("nan")


# --------------------------------------------- #
# Paper Parsing
# --------------------------------------------- #

def insert_tables_and_figures(paper: dict, section: str) -> str:
    '''
    Insert tables and figures referenced in the [section]
    '''
    line = ''
    for img in paper['table_figure_positions'][section]:
        if 'Figure' in img:
            figure_caption = paper['figures'][img]

            line += f'{img}: {figure_caption}'
            line += '\n\n'
        
        elif 'Table' in img:
            table_caption = paper['tables'][img]['caption']
            table_content = paper['tables'][img]['content']

            line += f'{img}: {table_caption}'
            line += '\n'
            line += table_content
            line += '\n\n'
    return line


def paper_generation(paper, use_fig=False, use_appendix=False):
    line = f"Title:\n{paper['title']}"
    line += '\n\n'
    line += f"Abstract:\n{paper['ABSTRACT']}"

    for section in paper['table_figure_positions']:
        if section == 'APPENDIX' and not use_appendix:
            continue

        line += '\n\n'
        if len(paper['table_figure_positions'][section]) > 0:
            line += insert_tables_and_figures(paper, section)
        line += f"{section}:\n{paper[f'[{section}]']}"
    return line


def paper_generation_ARR(paper, use_fig=False, use_appendix=False):
    paper = paper['metadata']
    line = f"## Title\n\n{paper['title']}\n\n### Abstract\n\n{paper['abstractText']}"

    for section in paper['sections']:
        if 'appendix' in section['heading'].lower() and not use_appendix:
            continue

        line += '\n\n'
        line += f'### {section['heading']}\n\n{section['text']}'
    
    return line


from pathlib import Path
import json
from typing import List, Tuple, Any

def load_json_dir(
    directory: str,
    *,
    recursive: bool = False
) -> Tuple[List[Any], List[str]]:
    """
    Load every  annotation_{annotator}.json  file in *directory*.

    Returns
    -------
    records : list
        Parsed JSON objects.
    annotators : list[str]
        Annotator names (derived from the file name) in the same order.
    """
    directory = Path(directory).expanduser().resolve()
    pattern = "**/annotation_*.json" if recursive else "annotation_*.json"

    records: List[Any] = []
    annotators: List[str] = []

    for path in sorted(directory.glob(pattern)):
        # ── extract the bit after "annotation_" and before ".json"
        annotator = path.stem.split("annotation_", 1)[1]

        try:
            with path.open(encoding="utf-8") as f:
                records.append(json.load(f))
            annotators.append(annotator)
        except json.JSONDecodeError as e:
            print(f"[WARN] {path} is not valid JSON: {e}")

    return records, annotators


# --------------------------------------------- #
# Alignment metrics
# --------------------------------------------- #
from sklearn.metrics import mean_absolute_error, cohen_kappa_score

def ccc(y_true, y_pred):
    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    mean_true = np.mean(y_true)
    mean_pred = np.mean(y_pred)
    var_true = np.var(y_true)
    var_pred = np.var(y_pred)
    cov = np.mean((y_true - mean_true) * (y_pred - mean_pred))
    return (2 * cov) / (var_true + var_pred + (mean_true - mean_pred) ** 2)

def tol_within_1(y_true, y_pred):
    return np.mean(np.abs(np.asarray(y_true) - np.asarray(y_pred)) <= 1.0)

# ---------- Utilities ----------
def _ensure_int_scale(y_true, y_pred):
    y_true = np.rint(y_true).astype(int)  # np.asarray(y_true).astype(int)
    y_pred = np.rint(y_pred).astype(int)  # np.asarray(y_pred).astype(int)
    y_pred = np.clip(y_pred, a_min=1, a_max=5)
    if y_true.shape != y_pred.shape:
        raise ValueError("y_true and y_pred must have the same shape.")
    return y_true, y_pred

def _quad_weights(k):
    # Quadratic agreement weights in [0,1]; w_ii = 1 (full credit), farther mistakes penalized more
    I, J = np.ogrid[:k, :k]
    return 1.0 - ((I - J) ** 2) / ((k - 1) ** 2)

def _confusion(y_true, y_pred, k, offset=1):
    # offset=1 if categories are 1..k; set 0 for 0..k-1
    cm = np.zeros((k, k), dtype=float)
    for t, p in zip(y_true, y_pred):
        cm[t - offset, p - offset] += 1.0
    return cm

def gwet_ac2_ordinal(y_true, y_pred, k=5):
    """
    Gwet's AC2 for two raters with ordinal (quadratic) weights.
    Reference idea: AC2 = 1 - D_o / D_e, where D_o is weighted disagreement observed,
    and D_e is weighted disagreement expected under chance using the *average marginal*
    category probabilities across raters.
    """
    y_true, y_pred = _ensure_int_scale(y_true, y_pred)
    N = float(len(y_true))
    if N == 0:
        raise ValueError("Empty inputs.")

    # Confusion and weights
    try:
        cm = _confusion(y_true, y_pred, k=k, offset=1)
    except:
        import pdb; pdb.set_trace()
    W = _quad_weights(k)  # agreement weights in [0,1]
    V = 1.0 - W          # convert to *disagreement* weights

    # Observed disagreement probability
    p_ij = cm / N
    Do = np.sum(V * p_ij)

    # Average marginals across the two raters (π_i)
    row_marg = cm.sum(axis=1) / N   # human distribution over categories
    col_marg = cm.sum(axis=0) / N   # model distribution over categories
    pi = 0.5 * (row_marg + col_marg)

    # Expected disagreement under chance using averaged marginals
    # (This avoids the kappa prevalence/bias paradox and is standard for AC2.)
    Pi = pi[:, None] * pi[None, :]
    De = np.sum(V * Pi)

    # Guard against degenerate cases
    if np.isclose(De, 0.0):
        return 1.0 if np.isclose(Do, 0.0) else 0.0

    ac2 = 1.0 - (Do / De)
    return float(ac2)

def icc_a1(y_true, y_pred):
    """
    ICC(A,1) per McGraw & Wong (1996): two-way random effects, absolute agreement, single measure.
    Falls back to a direct ANOVA-based computation if pingouin is unavailable.
    """
    y_true, y_pred = _ensure_int_scale(y_true, y_pred)
    x = np.column_stack([y_true, y_pred]).astype(float)  # shape (n, k=2)
    n, k = x.shape

    # Try pingouin (nice CI & extras if installed)
    try:
        import pandas as pd
        import pingouin as pg
        df = pd.DataFrame({
            "targets": np.repeat(np.arange(n), k),
            "raters":  np.tile(["human", "model"], n),
            "scores":  x.flatten(order="C")
        })
        icc_tbl = pg.intraclass_corr(data=df, targets="targets", raters="raters", ratings="scores")
        # ICC2 == two-way random, absolute agreement, single rater  (Pingouin naming)
        val = float(icc_tbl.loc[icc_tbl["Type"] == "ICC2", "ICC"].values[0])
        return val
    except Exception:
        pass  # fall back to pure Python

    # Pure-Python ANOVA components for ICC(A,1)
    grand_mean = x.mean()
    row_means = x.mean(axis=1, keepdims=True)  # targets
    col_means = x.mean(axis=0, keepdims=True)  # raters

    ss_rows = k * np.sum((row_means - grand_mean) ** 2)
    ss_cols = n * np.sum((col_means - grand_mean) ** 2)
    ss_total = np.sum((x - grand_mean) ** 2)
    ss_error = ss_total - ss_rows - ss_cols

    ms_rows = ss_rows / (n - 1) if n > 1 else 0.0
    ms_cols = ss_cols / (k - 1) if k > 1 else 0.0
    ms_error = ss_error / ((n - 1) * (k - 1)) if (n > 1 and k > 1) else 0.0

    # McGraw & Wong (1996), absolute-agreement, two-way random, single measurement:
    # ICC(A,1) = (MSR - MSE) / (MSR + (k-1)MSE + k*(MSC - MSE)/n)
    numerator = ms_rows - ms_error
    denominator = ms_rows + (k - 1) * ms_error + (k * (ms_cols - ms_error) / n if n > 0 else 0.0)
    if np.isclose(denominator, 0.0):
        return 0.0
    return float(numerator / denominator)

def compute_balanced_metrics(score_name, human_scores, pred, N=3):
    if len(human_scores) == 0 or len(pred) == 0:
        print('Empty')
        return

    # mae = mean_absolute_error(human_scores, pred)
    # ccc_val = ccc(human_scores, pred)
    kappa = cohen_kappa_score(np.rint(human_scores).astype(int),    # np.asarray(human_scores, dtype=int),
                              np.rint(pred).astype(int),
                              weights="quadratic")
    # tol = tol_within_1(human_scores, pred)

    ac2 = gwet_ac2_ordinal(human_scores, pred)
    # icc = icc_a1(human_scores, pred)

    # alpha = krippendorff.alpha(np.array([human_scores, pred]), level_of_measurement='ordinal')       # interval

    out = {
        # 'mae': mae,
        # 'ccc': ccc_val,
        'ac2': round(ac2, N),
        'kappa': round(kappa, N),
        # 'tol': tol,
        # 'icc_a1': round(icc, N),
        # 'alpha': round(alpha,N)
    }
    print(score_name)
    print(out)

def map_1to5_to_binary(gt, pred):
    out_gt, out_pred = [], []
    for data_gt, data_pred in zip(gt, pred):
        if data_gt == 3 or data_pred == 3:
            continue

        if data_gt > 3:
            out_gt.append(0)
        else:
            out_gt.append(1)

        if data_pred > 3:
            out_pred.append(0)
        else:
            out_pred.append(1)
    assert len(out_gt) == len(out_pred)
    return out_gt, out_pred

def agg_1to5_to_binary(xs):
    for x in xs:
        if x < 3:
            return 1
    return 0

def binary_classification_metrics(score_name, gt, pred):
    """
    Compute accuracy, precision, recall, and F1 for binary labels.

    Args:
        gt: Ground-truth labels (each 0 or 1).
        pred: Predicted labels (each 0 or 1).

    Returns:
        dict with keys: accuracy, precision, recall, f1, tp, tn, fp, fn

    Raises:
        ValueError: if lengths mismatch or labels are not in {0,1}.
    """
    print(f'\n**{score_name} Alignment**')

    if len(gt) != len(pred):
        raise ValueError(f"Length mismatch: gt={len(gt)} vs pred={len(pred)}")

    tp = tn = fp = fn = 0
    for g, p in zip(gt, pred):
        # Normalize and validate
        try:
            g_ = int(g)
            p_ = int(p)
        except Exception:
            raise ValueError("Labels must be 0 or 1 (ints or bools).")
        if g_ not in (0, 1) or p_ not in (0, 1):
            raise ValueError("Labels must be 0 or 1 (ints or bools).")

        if g_ == 1 and p_ == 1: tp += 1
        elif g_ == 0 and p_ == 0: tn += 1
        elif g_ == 0 and p_ == 1: fp += 1
        else: fn += 1  # g_ == 1 and p_ == 0

    n = tp + tn + fp + fn
    accuracy = (tp + tn) / n if n > 0 else 0.0
    precision = tp / (tp + fp) if (tp + fp) > 0 else 0.0
    recall = tp / (tp + fn) if (tp + fn) > 0 else 0.0
    f1 = (2 * precision * recall / (precision + recall)) if (precision + recall) > 0 else 0.0

    print({
        # "accuracy": accuracy,
        "precision": round(precision, 3),
        "recall": round(recall, 3),
        "f1": round(f1, 3),
        # "tp": float(tp),
        # "tn": float(tn),
        # "fp": float(fp),
        # "fn": float(fn),
    })


# pdf parsing

from pathlib import Path
import re, os
import json
import difflib
from typing import List, Dict

try:
    import PyPDF2
except Exception as e:
    raise RuntimeError(f"Failed to import PyPDF2: {e}")

def extract_pdf_text(path: Path) -> str:
    reader = PyPDF2.PdfReader(str(path))
    texts = []
    for page in reader.pages:
        try:
            texts.append(page.extract_text() or "")
        except Exception:
            texts.append("")
    return "\n".join(texts)

# Remove page headers like "Under review as a conference paper at ICLR 20XX"
header_pattern = re.compile(r"^Under review as a conference paper at ICLR\s*\d{4}\s*$", re.IGNORECASE | re.MULTILINE)

def strip_headers(text: str) -> str:
    text = header_pattern.sub("", text)
    # Remove page numbers that are a single number on a line
    text = re.sub(r"(?m)^\s*\d+\s*$", "", text)
    return text

# Truncate at the start of Appendix / APPENDIX
def truncate_at_appendix(text: str) -> str:
    m = re.search(r"(?mi)^(appendix|appendices)\b.*$", text)
    if m:
        return text[:m.start()].strip()
    return text.strip()

# Normalize for comparison (ignore line breaks and simple hyphenation artifacts)
def normalize_for_compare(s: str) -> str:
    s2 = re.sub(r"[ \t]+", " ", s)
    s2 = s2.replace("\u00ad", "")  # soft hyphen
    s2 = re.sub(r"-\s*\n\s*", "", s2)  # join line-broken hyphen words
    s2 = s2.replace("\n", " ")
    s2 = re.sub(r"\s+", " ", s2).strip()
    return s2

# Sentence/Unit segmentation (v1 heuristics)
def sentence_tokenize(text: str) -> List[str]:
    parts = re.split(r"\n\s*\n", text)
    units: List[str] = []
    for p in parts:
        p = p.strip()
        if not p:
            continue
        # Keep captions intact: expecting "Table 1:" etc.
        if re.match(r"^(Table|Figure|Algorithm)\s*\d+:", p):
            units.append(p)
            continue
        # Keep section headings that look like all-caps or title-case numerals line (stricter v1 rule)
        if re.match(r"^\d+(\.\d+)*\s+[A-Z][A-Z0-9 \-/]+$", p):
            units.append(p)
            continue
        # Conservative sentence split
        splits = re.split(r"(?<=[.?!])\s+(?=[A-Z0-9(])", p)
        for s in splits:
            s = s.strip()
            if s:
                units.append(s)
    return units

def build_section_map(units: List[str]) -> Dict[int, str]:
    section_map: Dict[int, str] = {}
    current_section = "Unknown section"
    for idx, u in enumerate(units):
        if re.match(r"^\d+(\.\d+)*\s+.+", u):
            current_section = u
        section_map[idx] = current_section
    return section_map

def similarity(a: str, b: str) -> float:
    return difflib.SequenceMatcher(
        a=normalize_for_compare(a).lower(),
        b=normalize_for_compare(b).lower(),
        autojunk=False
    ).ratio()

def is_trivial_change(before_text: str, after_text: str) -> bool:
    # High similarity => trivial (formatting/grammar/paraphrase)
    if similarity(before_text, after_text) >= 0.97:
        return True
    # Ignore renumbering of captions when caption text is effectively the same
    caption_num_change = re.compile(r"^(Table|Figure|Algorithm)\s*\d+:")
    if caption_num_change.match(before_text) and caption_num_change.match(after_text):
        before_caption = re.sub(r"^(Table|Figure|Algorithm)\s*\d+:\s*", "", before_text).strip()
        after_caption  = re.sub(r"^(Table|Figure|Algorithm)\s*\d+:\s*", "", after_text).strip()
        if similarity(before_caption, after_caption) >= 0.98:
            return True
    # Ignore pure equation/footnote numbering differences like "(1)" -> "(2)"
    b_clean = re.sub(r"\(\s*\d+\s*\)", "()", normalize_for_compare(before_text))
    a_clean = re.sub(r"\(\s*\d+\s*\)", "()", normalize_for_compare(after_text))
    if b_clean == a_clean:
        return True
    return False

def make_summary(b: str, a: str) -> str:
    if b and a:
        # if "train" in a.lower() and "valid" in a.lower() and "test" in a.lower():
        #     return "Updated dataset split statistics/descriptions"
        # if "EASY" in a or "HARD" in a or "biased" in a.lower():
        #     return "Revised description of EASY/HARD or bias analysis"
        if "Table" in a or "Figure" in a:
            return "Changed table/figure caption or content"
        # if "RACE" in a or "DREAM" in a:
        #     return "Revised related work or transfer learning details"
        return "Edited prose in the section"
    elif b and not a:
        return "Removed content"
    elif a and not b:
        return "Added new content"
    return "Change"

def find_changes(before_units: List[str], after_units: List[str]) -> list:
    before_norm = [normalize_for_compare(u).lower() for u in before_units]
    after_norm  = [normalize_for_compare(u).lower() for u in after_units]
    sm = difflib.SequenceMatcher(a=before_norm, b=after_norm, autojunk=False)

    changes = []
    for tag, i1, i2, j1, j2 in sm.get_opcodes():
        if tag == "equal":
            continue
        before_block = "\n".join(before_units[i1:i2]).strip()
        after_block  = "\n".join(after_units[j1:j2]).strip()
        if before_block and after_block and is_trivial_change(before_block, after_block):
            continue
        # Determine location using AFTER doc
        after_section_map = build_section_map(after_units)
        loc_idx = j1 if (j1 < len(after_units) and after_units[j1].strip()) else max(j1 - 1, 0)
        location = after_section_map.get(loc_idx, "Unknown section")
        changes.append({
            "summary": make_summary(before_block, after_block),
            "location": location,
            "before": before_block,
            "after": after_block
        })

    # Post-filter: drop tiny punctuation-only diffs
    filtered = []
    for ch in changes:
        b = normalize_for_compare(ch["before"])
        a = normalize_for_compare(ch["after"])
        if b == "" and a == "":
            continue
        if len(b) + len(a) < 40 and similarity(b, a) > 0.9:
            continue
        filtered.append(ch)
    return filtered

def extract_main_body_text_pageaware(pdf_path: Path) -> str:
    import PyPDF2, re
    reader = PyPDF2.PdfReader(str(pdf_path))
    pages_out = []
    cutoff_seen = False
    REF_RE = re.compile(r'(?mi)^\s*references\s*(?:$|\n)')

    for page in reader.pages:
        if cutoff_seen:
            break

        txt = page.extract_text() or ""

        # If a "References" heading appears on this page, keep only the text BEFORE it
        m = REF_RE.search(txt)
        if m:
            pages_out.append(txt[:m.start()].rstrip())
            cutoff_seen = True
        else:
            pages_out.append(txt)

    return "\n".join(pages_out).strip()

# alignment with external signals

from math import sqrt
from collections import Counter
import random

def agreement_metrics(a, b, bootstrap_ci=False, R=2000, seed=0):
    assert len(a) == len(b) and len(a) > 0
    # Confusion counts
    TP = sum(x==1 and y==1 for x,y in zip(a,b))
    TN = sum(x==0 and y==0 for x,y in zip(a,b))
    FP = sum(x==1 and y==0 for x,y in zip(a,b))
    FN = sum(x==0 and y==1 for x,y in zip(a,b))
    n  = TP+TN+FP+FN

    # Raw agreement
    p0 = (TP+TN)/n

    # Marginals for kappa
    pA1, pA0 = (TP+FP)/n, (TN+FN)/n
    pB1, pB0 = (TP+FN)/n, (TN+FP)/n
    pe = pA1*pB1 + pA0*pB0
    # kappa = (p0 - pe) / (1 - pe) if (1 - pe) != 0 else None  # None if undefined
    kappa = cohen_kappa_score(b,    # np.asarray(human_scores, dtype=int),
                              a,
                              weights="quadratic")

    # MCC / phi
    denom = (TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)
    mcc = (TP*TN - FP*FN) / sqrt(denom) if denom > 0 else None

    # Jaccard & F1 (Dice) on positives
    pos_denom = TP + FP + FN
    jaccard = TP / pos_denom if pos_denom > 0 else (1.0 if (TP==FP==FN==0) else None)
    f1 = (2*TP) / (2*TP + FP + FN) if (2*TP + FP + FN) > 0 else (1.0 if (TP==FP==FN==0) else None)

    # PABAK and Hamming distance rate
    pabak = 2*p0 - 1
    hamming_rate = 1 - p0

    # McNemar components (b=FP, c=FN). Report chi2 with continuity correction.
    b, c = FP, FN
    if (b + c) > 0:
        chi2_cc = (abs(b - c) - 1)**2 / (b + c)
    else:
        chi2_cc = None  # no discordant pairs

    out = {
        "n": n, "TP": TP, "TN": TN, "FP": FP, "FN": FN,
        "raw_agreement": p0,
        "kappa": kappa,
        "MCC": mcc,
        "Jaccard": jaccard,
        "F1_on_positives": f1,
        "PABAK": pabak,
        "hamming_rate": hamming_rate,
        "mcnemar_b": b, "mcnemar_c": c, "mcnemar_chi2_cc": chi2_cc,
    }

    # Optional: simple bootstrap CIs for κ and MCC
    if bootstrap_ci:
        rng = random.Random(seed)
        idx = list(range(n))
        def boot_stat(fn):
            vals = []
            for _ in range(R):
                samp = [rng.choice(idx) for _ in range(n)]
                aa = [a[i] for i in samp]; bb = [b[i] for i in samp]
                vals.append(fn(aa, bb))
            vals = [v for v in vals if v is not None]
            if not vals: return None
            vals.sort()
            lo = vals[int(0.025*len(vals))]
            hi = vals[int(0.975*len(vals))-1]
            return (lo, hi)
        def _k(aa, bb):
            TP = sum(x==1 and y==1 for x,y in zip(aa,bb))
            TN = sum(x==0 and y==0 for x,y in zip(aa,bb))
            FP = sum(x==1 and y==0 for x,y in zip(aa,bb))
            FN = sum(x==0 and y==1 for x,y in zip(aa,bb))
            n = TP+TN+FP+FN
            p0 = (TP+TN)/n
            pA1, pA0 = (TP+FP)/n, (TN+FN)/n
            pB1, pB0 = (TP+FN)/n, (TN+FP)/n
            pe = pA1*pB1 + pA0*pB0
            return (p0 - pe)/(1 - pe) if (1 - pe)!=0 else None
        def _m(aa, bb):
            TP = sum(x==1 and y==1 for x,y in zip(aa,bb))
            TN = sum(x==0 and y==0 for x,y in zip(aa,bb))
            FP = sum(x==1 and y==0 for x,y in zip(aa,bb))
            FN = sum(x==0 and y==1 for x,y in zip(aa,bb))
            denom = (TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)
            return (TP*TN-FP*FN)/sqrt(denom) if denom>0 else None
        out["kappa_ci95"] = boot_stat(_k)
        out["MCC_ci95"] = boot_stat(_m)

    return out