{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":35297,"databundleVersionId":3589144}],"dockerImageVersionId":31287,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-03-03T23:57:08.925707Z","iopub.execute_input":"2026-03-03T23:57:08.926079Z","iopub.status.idle":"2026-03-03T23:57:10.963790Z","shell.execute_reply.started":"2026-03-03T23:57:08.926032Z","shell.execute_reply":"2026-03-03T23:57:10.963104Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =========================================================\n# 0. IMPORTS\n# =========================================================\nimport os\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nimport torch\nfrom torch import nn, optim\nfrom torch.utils.data import DataLoader, TensorDataset\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import (\n    roc_auc_score,\n    average_precision_score,\n    precision_score,\n    recall_score,\n    roc_curve,\n    precision_recall_curve\n)\nfrom sklearn.model_selection import KFold\n\nfrom scipy.signal import savgol_filter\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(\"Device:\", device)\ntorch.cuda.empty_cache()\n\n# =========================================================\n# 1. DATA PREPARATION\n# =========================================================\nSEQ_LEN = 8\nDIFF_STEPS = 100\nMC_SAMPLES = 5   # Monte Carlo passes (PR booster)\n\ndef prepare_data(path):\n    df = pd.read_csv(path)\n\n    df[\"timestamp\"] = pd.to_datetime(df[\"timestamp\"], format=\"mixed\", errors=\"coerce\")\n    df = df.dropna(subset=[\"timestamp\"])\n    df = df.sort_values([\"building_id\", \"timestamp\"])\n\n    df[\"meter_log\"] = np.log1p(df[\"meter_reading\"])\n    df[\"hour\"] = df[\"timestamp\"].dt.hour\n    df[\"dow\"] = df[\"timestamp\"].dt.dayofweek\n\n    norm = (\n        df[df[\"anomaly\"] == 0]\n        .groupby([\"building_id\", \"dow\", \"hour\"])[\"meter_log\"]\n        .agg([\"median\", \"std\"])\n        .reset_index()\n    )\n    norm.columns = [\"building_id\", \"dow\", \"hour\", \"h_med\", \"h_std\"]\n    df = df.merge(norm, on=[\"building_id\", \"dow\", \"hour\"], how=\"left\")\n\n    df[\"h_med\"] = df[\"h_med\"].fillna(\n        df.groupby(\"building_id\")[\"meter_log\"].transform(\"median\")\n    )\n\n    df[\"temp_inertia\"] = (\n        df.groupby(\"building_id\")[\"air_temperature\"]\n        .transform(lambda x: x.rolling(6, min_periods=1).mean())\n    )\n\n    df[\"thermal_p\"] = df[\"air_temperature\"] - df[\"temp_inertia\"]\n\n    df[\"shock\"] = df[\"meter_log\"] - df[\"h_med\"]\n    df[\"shock_v\"] = df.groupby(\"building_id\")[\"shock\"].diff().fillna(0)\n    df[\"kinetic_s\"] = df.groupby(\"building_id\")[\"shock_v\"].diff().fillna(0)\n\n    df[\"shock_scaled\"] = df.groupby(\"building_id\")[\"shock\"].transform(\n        lambda x: (x - x.median()) /\n        (x.quantile(0.75) - x.quantile(0.25) + 1e-6)\n    )\n\n    feat_cols = [\"shock_scaled\", \"thermal_p\", \"kinetic_s\", \"air_temperature\"]\n    df[feat_cols] = df[feat_cols].replace([np.inf, -np.inf], 0).fillna(0)\n\n    return df, feat_cols\n\n# =========================================================\n# 2. SEQUENCE WINDOWING\n# =========================================================\ndef create_sequences(df, feat_cols):\n\n    X, y, bid = [], [], []\n\n    for building_id, group in df.groupby(\"building_id\"):\n        group = group.sort_values(\"timestamp\")\n\n        feats = group[feat_cols].values\n        labels = group[\"anomaly\"].values\n\n        for i in range(len(group) - SEQ_LEN):\n            X.append(feats[i:i+SEQ_LEN])\n            y.append(labels[i+SEQ_LEN-1])\n            bid.append(building_id)\n\n    return (\n        torch.tensor(np.array(X), dtype=torch.float32),\n        np.array(y),\n        np.array(bid)\n    )\n\n# =========================================================\n# 3. VANILLA DIFFUSION MODEL (Predict Noise)\n# =========================================================\nclass VanillaDiffusion(nn.Module):\n    def __init__(self, inp):\n        super().__init__()\n        self.net = nn.Sequential(\n            nn.Linear(inp, 128),\n            nn.SiLU(),\n            nn.Linear(128, 128),\n            nn.SiLU(),\n            nn.Linear(128, inp)\n        )\n\n    def forward(self, x):\n        return self.net(x)\n\n# =========================================================\n# 4. TRAIN + MONTE CARLO SCORING\n# =========================================================\ndef train_kfold(X, y, bid):\n\n    kf = KFold(n_splits=3, shuffle=True, random_state=42)\n    ensemble_scores = np.zeros(len(y))\n\n    betas = torch.linspace(1e-4, 0.02, DIFF_STEPS).to(device)\n    alphas = 1 - betas\n    alpha_bar = torch.cumprod(alphas, dim=0)\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        print(f\"\\n🔥 Fold {fold+1}\")\n\n        model = VanillaDiffusion(X.shape[2]).to(device)\n        opt = optim.AdamW(model.parameters(), lr=1e-3)\n\n        normal_mask = y[train_idx] == 0\n        X_train = X[train_idx][normal_mask][:, -1, :]\n\n        train_loader = DataLoader(\n            TensorDataset(X_train),\n            batch_size=512,\n            shuffle=True\n        )\n\n        # Train to predict noise ε\n        for epoch in range(20):\n            model.train()\n            for (xb,) in train_loader:\n                xb = xb.to(device)\n\n                t = torch.randint(0, DIFF_STEPS, (xb.size(0),), device=device)\n                noise = torch.randn_like(xb)\n\n                a_bar = alpha_bar[t].unsqueeze(1)\n                x_t = torch.sqrt(a_bar) * xb + torch.sqrt(1 - a_bar) * noise\n\n                pred_noise = model(x_t)\n                loss = nn.MSELoss()(pred_noise, noise)\n\n                opt.zero_grad()\n                loss.backward()\n                torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n                opt.step()\n\n        # Inference with Monte Carlo\n        model.eval()\n        scores = []\n\n        with torch.no_grad():\n            for i in range(0, len(X), 512):\n                xb = X[i:i+512][:, -1, :].to(device)\n                total_err = torch.zeros(xb.size(0), device=device)\n\n                for _ in range(MC_SAMPLES):\n                    t = torch.full((xb.size(0),), DIFF_STEPS//2, device=device)\n                    noise = torch.randn_like(xb)\n\n                    a_bar = alpha_bar[t].unsqueeze(1)\n                    x_t = torch.sqrt(a_bar) * xb + torch.sqrt(1 - a_bar) * noise\n\n                    pred_noise = model(x_t)\n\n                    x_recon = (x_t - torch.sqrt(1 - a_bar) * pred_noise) / torch.sqrt(a_bar)\n\n                    err = (x_recon - xb) ** 2\n                    err = err[:,0] + 0.75*err[:,2] + 0.5*err[:,1]\n\n                    total_err += err\n\n                total_err /= MC_SAMPLES\n                scores.append(total_err.cpu())\n\n        fold_scores = torch.cat(scores).numpy()\n\n        # Savgol smoothing per building\n        df_tmp = pd.DataFrame({\"score\": fold_scores, \"b\": bid})\n        smooth = df_tmp.groupby(\"b\")[\"score\"].transform(\n            lambda x: savgol_filter(x, 11, 2) if len(x) > 11 else x\n        )\n\n        ensemble_scores += smooth / 3\n\n    return ensemble_scores\n\n# =========================================================\n# 5. EVALUATION\n# =========================================================\ndef evaluate(scores, y, bid):\n\n    scores = np.nan_to_num(scores)\n\n    auc = roc_auc_score(y, scores)\n    pr = average_precision_score(y, scores)\n\n    thresh = np.percentile(scores, 95)\n    preds = (scores > thresh).astype(int)\n\n    precision = precision_score(y, preds)\n    recall = recall_score(y, preds)\n    fpr = np.sum((preds==1)&(y==0)) / np.sum(y==0)\n\n    print(\"\\n🏆 FINAL RESULTS\")\n    print(\"ROC-AUC :\", auc)\n    print(\"PR-AUC  :\", pr)\n    print(\"Precision:\", precision)\n    print(\"Recall:\", recall)\n    print(\"TPR:\", recall)\n    print(\"FPR:\", fpr)\n\n    fpr_c, tpr_c, _ = roc_curve(y, scores)\n    plt.plot(fpr_c, tpr_c)\n    plt.title(\"ROC Curve\")\n    plt.show()\n\n    p_c, r_c, _ = precision_recall_curve(y, scores)\n    plt.plot(r_c, p_c)\n    plt.title(\"PR Curve\")\n    plt.show()\n\n    # =========================================================\n    # Predicted vs Actual Anomaly Points (Sample Building)\n    # =========================================================\n    \n    sample_building = bid[0]   # pick one building\n    mask = bid == sample_building\n    \n    plt.figure(figsize=(12,5))\n    \n    # predicted anomaly score\n    plt.plot(scores[mask], label=\"Predicted Anomaly Score\")\n    \n    # actual anomaly points\n    plt.scatter(\n        np.where(y[mask] == 1),\n        scores[mask][y[mask] == 1],\n        color=\"red\",\n        label=\"Actual Anomaly\"\n    )\n    \n    plt.title(f\"Predicted vs Actual Anomalies (Building {sample_building})\")\n    plt.xlabel(\"Time Step\")\n    plt.ylabel(\"Anomaly Score\")\n    plt.legend()\n    plt.show()\n\n\n# =========================================================\n# RUN\n# =========================================================\npath = \"/kaggle/input/competitions/energy-anomaly-detection/train_features.csv\"\n\ndf, feats = prepare_data(path)\nX, y, bid = create_sequences(df, feats)\n\nscores = train_kfold(X, y, bid)\nevaluate(scores, y, bid)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-04T11:58:16.792119Z","iopub.execute_input":"2026-03-04T11:58:16.792405Z","iopub.status.idle":"2026-03-04T12:07:47.714124Z","shell.execute_reply.started":"2026-03-04T11:58:16.792380Z","shell.execute_reply":"2026-03-04T12:07:47.713419Z"}},"outputs":[],"execution_count":null}]}