diff --git a/.gitignore b/.gitignore index 64d49ae..a36be78 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,9 @@ __pycache__/ # C extensions *.so +data/ +data_results/ + # Distribution / packaging .Python build/ diff --git a/README.md b/README.md index 9219cd4..a399ce3 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,15 @@ # forest_python +Tady jsi v branchi youden + +musis mit +data/TEVA +data/TRVA + +a musis mit zscore a zscore_parameters +z projektu +https://github.com/chazzka/depth_estimation + ## Setup (venv) ``` cd forest_python diff --git a/roc_hstree_experiment.py b/roc_hstree_experiment.py new file mode 100644 index 0000000..59fbbf4 --- /dev/null +++ b/roc_hstree_experiment.py @@ -0,0 +1,338 @@ +""" +ROC experiment: HS-Tree on 10 representative ADBench datasets. + +For each dataset: + - Ranges loaded from zscore/ (param_opt per feature) + zscore_parameters/ (mean/std/fallback) + - 15 random seeds, median-seed run selected by ROC AUC + - ROC curve plotted with: + x = fixed threshold operating points (configurable, default shown: 0.3, 0.5, 0.6, 0.7) + o = Youden's index maximising threshold (TPR - FPR) + +Model: 30 trees, max depth 30, batch size min(1024, |train|), c(n) normalisation. +""" + +from __future__ import annotations + +import csv +import os +import random +import re +import sys + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.lines as mlines + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) + +from forest_python import Forest +from forest_python.service import novelty +from forest_python.service.outlier import average_path_length_c + +# --------------------------------------------------------------------------- +# config +# --------------------------------------------------------------------------- +DATASETS = [ + "18_Ionosphere", + "23_mammography", + "26_optdigits", + "27_PageBlocks", + "28_pendigits", + "38_thyroid", + "41_Waveform", + "42_WBC", + "44_Wilt", + "45_wine", +] + +N_TREES = 30 +MAX_DEPTH = 30 +BATCH_SIZE = 1024 +N_SEEDS = 15 +SEED_START = 100 + +# Fixed thresholds to show as cross markers (can be changed freely) +FIXED_THRESHOLDS = [0.5, 0.6, 0.7] + +ROOT = os.path.dirname(os.path.abspath(__file__)) +DATA_DIR = os.path.join(ROOT, "data") +ZSCORE_DIR = os.path.join(ROOT, "zscore") +ZSCORE_PAR_DIR = os.path.join(ROOT, "zscore_parameters") + +# --------------------------------------------------------------------------- +# load optimal ranges +# --------------------------------------------------------------------------- + +def _auc_from_filename(fname: str) -> float: + m = re.search(r"auc([\dp]+)\.csv$", fname) + return float(m.group(1).replace("p", ".")) if m else 0.0 + + +def load_ranges(name: str) -> tuple[list, str]: + """Return (ranges, source_filename) using param_opt from best zscore/ file.""" + files = [f for f in os.listdir(ZSCORE_DIR) if f.startswith(name + "_")] + best_file = max(files, key=_auc_from_filename) + + # param_opt per feature from the best run + param_opt: dict[int, float] = {} + with open(os.path.join(ZSCORE_DIR, best_file)) as f: + for row in csv.DictReader(f): + param_opt[int(row["feature_index"])] = float(row["param_opt"]) + + # mean/std and fallback ranges from zscore_parameters/ + mean_std: dict[int, tuple[float, float]] = {} + fallback: dict[int, tuple[float, float]] = {} + with open(os.path.join(ZSCORE_PAR_DIR, f"{name}.csv")) as f: + for row in csv.DictReader(f): + d = int(row["feature_index"]) + if row["mean"] and row["std"]: + mean_std[d] = (float(row["mean"]), float(row["std"])) + if row["range_l"] != "" and row["range_u"] != "": + fallback[d] = (float(row["range_l"]), float(row["range_u"])) + + n_dims = max(param_opt.keys()) + 1 + ranges = [] + for d in range(n_dims): + k = param_opt.get(d, 0.0) + if k > 0.0 and d in mean_std: + mu, sd = mean_std[d] + ranges.append((mu - k * sd, mu + k * sd)) + elif d in fallback: + ranges.append(fallback[d]) + else: + # constant feature: zero-width range is fine for HST (no split possible) + ranges.append((0.0, 0.0)) + + return ranges, best_file + + +# --------------------------------------------------------------------------- +# data loading +# --------------------------------------------------------------------------- + +def _parse_csv(path: str) -> list: + rows = [] + with open(path) as f: + for line in f: + line = line.strip() + if line: + rows.append(line.split(",")) + return rows + + +def load_dataset(name: str) -> tuple[list, list, list]: + """ + Returns (train, test_features, test_labels). + + TRVA.csv — all normal (label 0): splits TR / VA / TE + TEVA.csv — all anomalies (label 1): splits VA / TE + + Train = TRVA[TR] + Test = TRVA[TE] (label 0) + TEVA[TE] (label 1) + """ + trva = _parse_csv(os.path.join(DATA_DIR, "TRVA", f"{name}_TRVA.csv")) + teva = _parse_csv(os.path.join(DATA_DIR, "TEVA", f"{name}_TEVA.csv")) + + train, test, labels = [], [], [] + + for row in trva: + feats = [float(v) for v in row[:-2]] + if row[-1].strip() == "TR": + train.append(feats) + elif row[-1].strip() == "TE": + test.append(feats) + labels.append(0) + + for row in teva: + feats = [float(v) for v in row[:-2]] + if row[-1].strip() == "TE": + test.append(feats) + labels.append(1) + + return train, test, labels + + +# --------------------------------------------------------------------------- +# model +# --------------------------------------------------------------------------- + +def score_hstree(train: list, test: list, ranges: list, seed: int) -> list[float]: + random.seed(seed) + bs = min(BATCH_SIZE, len(train)) + dp = {"data": train, "ranges": ranges, "batch_size": bs} + forest = Forest.init(N_TREES, dp, novelty.make_split(MAX_DEPTH), novelty.batch) + scores = [] + for point in test: + leaves = Forest.evaluate(forest, point, novelty.decision) + depths = [leaf["depth"] for leaf in leaves] + scores.append(2.0 ** (-novelty.avg(depths) / average_path_length_c(bs))) + return scores + + +# --------------------------------------------------------------------------- +# ROC helpers +# --------------------------------------------------------------------------- + +def roc_curve(scores: list[float], labels: list[int]) -> tuple[list, list]: + n_pos = sum(labels) + n_neg = len(labels) - n_pos + fprs, tprs = [0.0], [0.0] + for t in sorted(set(scores), reverse=True): + tp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 1) + fp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 0) + fprs.append(fp / n_neg if n_neg else 0.0) + tprs.append(tp / n_pos if n_pos else 0.0) + fprs.append(1.0) + tprs.append(1.0) + return fprs, tprs + + +def roc_auc(fprs: list, tprs: list) -> float: + return sum( + (fprs[i] - fprs[i - 1]) * (tprs[i] + tprs[i - 1]) / 2.0 + for i in range(1, len(fprs)) + ) + + +def op_at_threshold(scores: list[float], labels: list[int], + t: float) -> tuple[float, float]: + n_pos = sum(labels) + n_neg = len(labels) - n_pos + tp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 1) + fp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 0) + return (fp / n_neg if n_neg else 0.0), (tp / n_pos if n_pos else 0.0) + + +def youden_op(scores: list[float], labels: list[int]) -> tuple[float, float]: + n_pos = sum(labels) + n_neg = len(labels) - n_pos + best_j, best_fpr, best_tpr = -1.0, 0.0, 0.0 + for t in sorted(set(scores), reverse=True): + tp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 1) + fp = sum(1 for s, l in zip(scores, labels) if s >= t and l == 0) + tpr = tp / n_pos if n_pos else 0.0 + fpr = fp / n_neg if n_neg else 0.0 + if tpr - fpr > best_j: + best_j, best_fpr, best_tpr = tpr - fpr, fpr, tpr + return best_fpr, best_tpr + + +# --------------------------------------------------------------------------- +# experiment +# --------------------------------------------------------------------------- + +def run_dataset(name: str) -> dict: + print(f" {name}") + train, test, labels = load_dataset(name) + ranges, src = load_ranges(name) + print(f" ranges: {src}") + + seed_results = [] + for seed in range(SEED_START, SEED_START + N_SEEDS): + scores = score_hstree(train, test, ranges, seed) + fprs, tprs = roc_curve(scores, labels) + seed_results.append((roc_auc(fprs, tprs), seed, scores, fprs, tprs)) + + seed_results.sort(key=lambda x: x[0]) + auc, seed, scores, fprs, tprs = seed_results[(N_SEEDS - 1) // 2] + print(f" median seed={seed} AUC={auc:.4f}") + + return { + "name": name, + "fprs": fprs, + "tprs": tprs, + "auc": auc, + "fixed_ops": [op_at_threshold(scores, labels, t) for t in FIXED_THRESHOLDS], + "youden_op": youden_op(scores, labels), + } + + +# --------------------------------------------------------------------------- +# plot +# --------------------------------------------------------------------------- + +# marker styles for fixed thresholds +THRESH_MARKERS = ["x", "+", "D", "s"] +THRESH_SIZES = [9, 9, 5, 5 ] + +def plot_all(results: list, out_path: str): + fig, ax = plt.subplots(figsize=(7, 6)) + colours = [c["color"] for c in plt.rcParams["axes.prop_cycle"]] + + # --- random guess --- + ax.plot([0, 1], [0, 1], "b--", lw=1.0, zorder=1) + + # --- ROC curves + markers --- + for i, res in enumerate(results): + col = colours[i % len(colours)] + ax.plot(res["fprs"], res["tprs"], color=col, lw=1.2, zorder=2) + + # fixed threshold markers + for j, (fx, fy) in enumerate(res["fixed_ops"]): + ax.plot(fx, fy, + marker=THRESH_MARKERS[j % len(THRESH_MARKERS)], + color=col, + markersize=THRESH_SIZES[j % len(THRESH_SIZES)], + markeredgewidth=1.6, + linestyle="none", + zorder=3, + clip_on=False) + + # Youden circle + yx, yy = res["youden_op"] + ax.plot(yx, yy, + marker="o", color=col, markersize=6, + markerfacecolor=col, linestyle="none", zorder=3, + clip_on=False) + + # --- legend --- + handles = [ + mlines.Line2D([], [], color="blue", linestyle="--", lw=1.0, + label="Random Guess"), + ] + for i, res in enumerate(results): + col = colours[i % len(colours)] + handles.append(mlines.Line2D([], [], color=col, lw=1.2, + label=res["name"])) + + # threshold marker legend entries + handles.append(mlines.Line2D([], [], linestyle="none", marker="o", + color="gray", markersize=6, + label="Youden optimum")) + for j, t in enumerate(FIXED_THRESHOLDS): + handles.append(mlines.Line2D( + [], [], linestyle="none", + marker=THRESH_MARKERS[j % len(THRESH_MARKERS)], + color="gray", + markersize=THRESH_SIZES[j % len(THRESH_SIZES)], + markeredgewidth=1.6, + label=f"threshold {t}")) + + ax.set_xlabel("False Positive Rate (FPR)") + ax.set_ylabel("True Positive Rate (TPR)") + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.grid(True, linestyle="--", linewidth=0.4, alpha=0.6) + ax.legend(handles=handles, loc="upper left", + bbox_to_anchor=(1.02, 1), borderaxespad=0, + fontsize=7, framealpha=0.9) + + fig.tight_layout() + fig.savefig(out_path, bbox_inches="tight") + print(f"\nSaved: {out_path}") + plt.close(fig) + + +# --------------------------------------------------------------------------- +# main +# --------------------------------------------------------------------------- + +if __name__ == "__main__": + out_path = os.path.join(ROOT, "roc_mean_hstree_all_replicated.pdf") + + results = [] + for name in DATASETS: + results.append(run_dataset(name)) + + plot_all(results, out_path) diff --git a/roc_mean_hstree_all_replicated.pdf b/roc_mean_hstree_all_replicated.pdf new file mode 100644 index 0000000..543c9d6 Binary files /dev/null and b/roc_mean_hstree_all_replicated.pdf differ