simulated.py

# Public modules
import json
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
from sklearn.metrics import classification_report, confusion_matrix
from pybocd import BOCDGMM

# Self modules
from utils.simulated import (
    plot_original_data,
    plot_map,
    plot_np,
    plot_comprehensive,
    plot_sigma,
    plot_pi,
)

# ==========================
# === Prepare parameters ===
# ==========================

# 1. Paths
DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
DATA_PATH = os.path.join(DATA_DIR, f"simulated.csv")
PLT_DIR = os.path.join(os.path.dirname(__file__), "plots", f"simulated")
LOG_DIR = os.path.join(os.path.dirname(__file__), "logs", f"simulated")
LOG_PATH = os.path.join(LOG_DIR, f"validating.txt")
if not os.path.exists(DATA_DIR):
    os.makedirs(DATA_DIR, exist_ok=True)
if not os.path.exists(PLT_DIR):
    os.makedirs(PLT_DIR, exist_ok=True)
if not os.path.exists(LOG_DIR):
    os.makedirs(LOG_DIR, exist_ok=True)

# 2. Data source
BOCD_START = 0
BOCD_END = 1000

# 3. Log
LOG_LEVEL = 2

# 4. Hyper parameters
ALPHA_0 = 20
BETA_0 = 19
M_0 = 4.5
KAPPA_0 = 0.1
ALPHA_P_0 = 20
BETA_P_0 = 190
MU_P_0 = -2.217
SIGMA_P_SQ_0 = 0.013

# 5. Parameters for algorithm
JITTER_MU = 0.01
JITTER_SIGMA_SQ = 0.01
JITTER_TAU_SQ = 0.01
JITTER_PI = 0.01
LAMBDA = 100
M = 400
N = 400
INIT_PARTICLE_N = 5
PROB_THRESHOLD = 0.1

if __name__ == "__main__":

    # Log (Lv.0)
    with open(LOG_PATH, "w") as f:
        f.write("----------\n")
        f.write("Start validating\n")

    # =================
    # === Load Data ===
    # =================

    # Raw data
    df = pd.read_csv(DATA_PATH)[BOCD_START:BOCD_END]

    # Colums
    data = df["data"].to_numpy()
    is_level_shift = df["is_level_shift"].to_numpy()
    is_anomaly = df["is_anomaly"].to_numpy()
    is_outlier = is_level_shift | is_anomaly
    baseline = df["baseline"].to_numpy()

    # Inferences
    level_shifts = np.where(is_level_shift)[0]
    anomalies = np.where(is_anomaly)[0]
    outliers = np.where(is_outlier)[0]

    # ============
    # === Init ===
    # ============

    # Create instance
    bocd = BOCDGMM(
        alpha_0=ALPHA_0,
        beta_0=BETA_0,
        m_0=M_0,
        kappa_0=KAPPA_0,
        alpha_p_0=ALPHA_P_0,
        beta_p_0=BETA_P_0,
        mu_p_0=MU_P_0,
        sigma_p_sq_0=SIGMA_P_SQ_0,
        jitter_mu=JITTER_MU,
        jitter_sigma_sq=JITTER_SIGMA_SQ,
        jitter_tau_sq=JITTER_TAU_SQ,
        jitter_pi=JITTER_PI,
        l=LAMBDA,
        m=M,
        n=N,
        init_particle_n=INIT_PARTICLE_N,
    )

    # MAP records
    max_path = np.zeros(1, dtype=int)

    # Parameters records
    mu_records = [np.mean(bocd.model.mu_particles)]
    sigma_sq_records = [np.mean(bocd.model.sigma_sq_particles)]
    tau_sq_records = [np.mean(bocd.model.tau_sq_particles)]
    pi_records = [np.mean(bocd.model.pi_particles)]

    # Prob records
    prob_records = np.array([], dtype=float)

    # =====================
    # === Run Algorithm ===
    # =====================

    map_mask = np.array([0], dtype=int)

    # Gradually insert data x
    for i in range(BOCD_END - BOCD_START):

        # === Before Op 1. Find Proportion ===
        probs_norm_map = bocd.model.pdf_norm(data[i])[map_mask]
        probs_out_map = bocd.model.pdf_out(data[i])[map_mask]
        pi_particles_map = bocd.model.pi_particles[map_mask]
        probs_map = (
            (1 - pi_particles_map)
            * probs_norm_map
            / (
                (1 - pi_particles_map) * probs_norm_map
                + pi_particles_map * probs_out_map
            )
        )
        prob = np.average(probs_map, weights=bocd.weights[map_mask])

        # === Main Operation: Add Data ===

        # 1. Update weights (n -> 2n + 1)
        # 2. Update parameters (n -> 2n + 1)
        # 3. Prune (≤ N)
        bocd.add_data(x=data[i])

        # === Later Op 1. Export States ===

        # Export particles' weights & run length
        weights = bocd.weights
        run_length = bocd.run_length

        # Export particles' parameters
        mu_particles = bocd.model.mu_particles
        sigma_sq_particles = bocd.model.sigma_sq_particles
        tau_sq_particles = bocd.model.tau_sq_particles
        pi_particles = bocd.model.pi_particles

        # === Later Op 2. Find MAP and Parameters ===

        # Find MAP
        combined_weights = np.bincount(run_length, weights=weights)
        r = np.argmax(combined_weights)
        map_mask = run_length == r

        # Find corresponding parameters (weighted average under MAP)
        mu = np.average(mu_particles[map_mask], weights=weights[map_mask])
        sigma_sq = np.average(sigma_sq_particles[map_mask], weights=weights[map_mask])
        tau_sq = np.average(tau_sq_particles[map_mask], weights=weights[map_mask])
        pi = np.average(pi_particles[map_mask], weights=weights[map_mask])

        # === Later Op 3. Add to Records ===

        # MAP records
        max_path = np.append(max_path, r)

        # Parameters records
        mu_records = np.append(mu_records, mu)
        sigma_sq_records = np.append(sigma_sq_records, sigma_sq)
        tau_sq_records = np.append(tau_sq_records, tau_sq)
        pi_records = np.append(pi_records, pi)

        # Prob records
        prob_records = np.append(prob_records, prob)

    # ===============================
    # === Compute Outlier Indices ===
    # ===============================

    outliers_predicted = np.where(prob_records < PROB_THRESHOLD)
    is_outlier_predicted = np.full_like(data, False, dtype=bool)
    is_outlier_predicted[outliers_predicted] = True
    cm = confusion_matrix(is_outlier, is_outlier_predicted, labels=[1, 0])
    cm_df = pd.DataFrame(
        cm,
        index=["Actual Positive (Outlier)", "Actual Negative (Normal)"],
        columns=["Predicted Positive (Outlier)", "Predicted Negative (Normal)"],
    )
    cr = classification_report(
        is_outlier, is_outlier_predicted, target_names=["Normal", "Anomaly"]
    )
    # Log (Lv.0)
    with open(LOG_PATH, "a") as f:
        f.write("----------\n")
        f.write(cm_df.to_string() + "\n")
        f.write("----------\n")
        f.write(cr + "\n")

    # ===================================
    # === Compute Level Shift Indices ===
    # ===================================

    is_level_shift_predicted = np.full_like(data, False, dtype=bool)
    for i in range(1, BOCD_END - BOCD_START):
        if max_path[i] < max_path[i - 1] - 10:
            flag = True
            for j in range(i, min(BOCD_END - BOCD_START, i + 5)):
                if max_path[j] >= max_path[i - 1]:
                    flag = False
                    break
            if flag:
                is_level_shift_predicted[i - max_path[i]] = True
    level_shifts_predicted = np.where(is_level_shift_predicted)

    # ================
    # === Plotting ===
    # ================

    # 1. Plot Original Data + Credible Intervals
    plot_original_data(
        data=data,
        mu_records=mu_records,
        sigma_sq_records=sigma_sq_records,
        tau_sq_records=tau_sq_records,
        outliers=outliers,
        outliers_predicted=outliers_predicted,
        level_shifts=level_shifts,
        level_shifts_predicted=level_shifts_predicted,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        pth=os.path.join(PLT_DIR, "original_data_ci.svg"),
    )

    # 2. Plot MAP
    plot_map(
        max_path=max_path,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        level_shifts=level_shifts,
        level_shifts_predicted=level_shifts_predicted,
        pth=os.path.join(PLT_DIR, "map.svg"),
    )

    # 3. Plot proportion
    plot_np(
        prob_records=prob_records,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        pth=os.path.join(PLT_DIR, "np.svg"),
    )

    # 4. Plot sigma
    plot_sigma(
        sigma_sq_records=sigma_sq_records,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        pth=os.path.join(PLT_DIR, "sigma.svg"),
    )

    # 5. Plot pi
    plot_pi(
        pi_records=pi_records,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        pth=os.path.join(PLT_DIR, "pi.svg"),
    )

    # 6. Plot comprehensive
    plot_comprehensive(
        data=data,
        mu_records=mu_records,
        sigma_sq_records=sigma_sq_records,
        tau_sq_records=tau_sq_records,
        outliers=outliers,
        outliers_predicted=outliers_predicted,
        level_shifts=level_shifts,
        level_shifts_predicted=level_shifts_predicted,
        max_path=max_path,
        prob_records=prob_records,
        bocd_start=BOCD_START,
        bocd_end=BOCD_END,
        pth=os.path.join(PLT_DIR, "comprehension.svg"),
    )