{"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:15:38.023795Z","iopub.execute_input":"2026-03-03T23:15:38.024481Z","iopub.status.idle":"2026-03-03T23:15:38.2899Z","shell.execute_reply.started":"2026-03-03T23:15:38.02445Z","shell.execute_reply":"2026-03-03T23:15:38.289128Z"}},"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\nfrom sklearn.ensemble import IsolationForest\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 = 12\nDIFF_STEPS = 50\nDIFF_T = 15\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    le_b, le_s, le_u = LabelEncoder(), LabelEncoder(), LabelEncoder()\n    df[\"b_idx\"] = le_b.fit_transform(df[\"building_id\"])\n    df[\"s_idx\"] = le_s.fit_transform(df[\"site_id\"])\n    df[\"u_idx\"] = le_u.fit_transform(df[\"primary_use\"])\n\n    return df, feat_cols, le_b, le_s, le_u\n\n# =========================================================\n# 2. SEQUENCE WINDOWING\n# =========================================================\ndef create_sequences(df, feat_cols):\n    X, b, s, u, 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        b_idx = group[\"b_idx\"].values\n        s_idx = group[\"s_idx\"].values\n        u_idx = group[\"u_idx\"].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            b.append(b_idx[i+SEQ_LEN-1])\n            s.append(s_idx[i+SEQ_LEN-1])\n            u.append(u_idx[i+SEQ_LEN-1])\n            bid.append(building_id)\n\n    return (\n        torch.tensor(np.array(X), dtype=torch.float32),\n        torch.tensor(b),\n        torch.tensor(s),\n        torch.tensor(u),\n        np.array(y),\n        np.array(bid)\n    )\n\n# =========================================================\n# 3. SIMPLE PIGD MODEL (NO TRANSFORMER)\n# =========================================================\nclass PIGD(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 + ISOLATION FOREST\n# =========================================================\ndef train_kfold(X, y):\n\n    kf = KFold(n_splits=3, shuffle=True, random_state=42)\n    ensemble_scores = np.zeros(len(y))\n\n    alphas = torch.linspace(1e-4, 0.02, DIFF_STEPS).to(device)\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        print(f\"\\n🔥 Fold {fold+1}\")\n\n        model = PIGD(X.shape[2]).to(device)\n        opt = optim.AdamW(model.parameters(), lr=1e-3)\n\n        # Train only on normal\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        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                alpha = alphas[t].unsqueeze(1)\n\n                noisy = torch.sqrt(1 - alpha) * xb + torch.sqrt(alpha) * noise\n                pred = model(noisy)\n\n                loss = nn.MSELoss()(pred, xb)\n\n                opt.zero_grad()\n                loss.backward()\n                torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n                opt.step()\n\n        # Compute reconstruction error for all\n        model.eval()\n        errors = []\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\n                t = torch.full((xb.size(0),), DIFF_T, device=device)\n                noise = torch.randn_like(xb)\n                alpha = alphas[t].unsqueeze(1)\n\n                noisy = torch.sqrt(1 - alpha) * xb + torch.sqrt(alpha) * noise\n                pred = model(noisy)\n\n                err = (pred - xb) ** 2\n                err = (\n                    err[:,0] +\n                    0.75 * err[:,2] +\n                    0.5 * err[:,1]\n                )\n\n                errors.append(err.cpu())\n\n        errors = torch.cat(errors).numpy()\n\n        # Isolation Forest on reconstruction errors\n        iso = IsolationForest(\n            n_estimators=200,\n            contamination=np.mean(y),\n            random_state=42\n        )\n\n        iso.fit(errors.reshape(-1,1))\n        scores = -iso.decision_function(errors.reshape(-1,1))\n\n        # Savgol smoothing per building\n        df_tmp = pd.DataFrame({\"score\": scores})\n        smooth = savgol_filter(df_tmp[\"score\"], 11, 2)\n\n        ensemble_scores += smooth / 3\n\n    return ensemble_scores\n\n# =========================================================\n# 5. EVALUATION\n# =========================================================\ndef evaluate(scores, y):\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# RUN\n# =========================================================\npath = \"/kaggle/input/competitions/energy-anomaly-detection/train_features.csv\"\n\ndf, feats, le_b, le_s, le_u = prepare_data(path)\nX, b, s, u, y, bid = create_sequences(df, feats)\n\nscores = train_kfold(X, y)\nevaluate(scores, y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-03T23:16:32.403539Z","iopub.execute_input":"2026-03-03T23:16:32.40436Z","iopub.status.idle":"2026-03-03T23:27:19.573559Z","shell.execute_reply.started":"2026-03-03T23:16:32.404327Z","shell.execute_reply":"2026-03-03T23:27:19.572924Z"}},"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\nfrom sklearn.preprocessing import LabelEncoder, StandardScaler\nfrom sklearn.ensemble import IsolationForest\nfrom sklearn.model_selection import KFold\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)\n\nfrom scipy.signal import savgol_filter\n\n# =========================================================\n# 1. DATA PREPARATION (PHYSICS FEATURES)\n# =========================================================\nSEQ_LEN = 6   # shorter window works better for IF\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    df[\"rolling_std\"] = (\n        df.groupby(\"building_id\")[\"meter_log\"]\n        .transform(lambda x: x.rolling(6, min_periods=1).std())\n        .fillna(0)\n    )\n\n    feat_cols = [\n        \"shock_scaled\",\n        \"thermal_p\",\n        \"kinetic_s\",\n        \"air_temperature\",\n        \"shock_v\",\n        \"rolling_std\"\n    ]\n\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            window = feats[i:i+SEQ_LEN]\n            last = window[-1]\n            X.append(last)\n            y.append(labels[i+SEQ_LEN-1])\n            bid.append(building_id)\n\n    return np.array(X), np.array(y), np.array(bid)\n\n# =========================================================\n# 3. ISOLATION FOREST K-FOLD\n# =========================================================\ndef train_isolation_forest(X, y, bid):\n\n    kf = KFold(n_splits=3, shuffle=True, random_state=42)\n    ensemble_scores = np.zeros(len(y))\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        print(f\"\\n🔥 Fold {fold+1}\")\n\n        # Train ONLY on normal samples\n        normal_mask = y[train_idx] == 0\n        X_train = X[train_idx][normal_mask]\n\n        scaler = StandardScaler()\n        X_train_scaled = scaler.fit_transform(X_train)\n\n        iso = IsolationForest(\n            n_estimators=400,\n            max_samples='auto',\n            contamination=np.mean(y),\n            random_state=42,\n            n_jobs=-1\n        )\n\n        iso.fit(X_train_scaled)\n\n        # Score all data\n        X_scaled = scaler.transform(X)\n        scores = -iso.decision_function(X_scaled)\n\n        # Savgol smoothing per building\n        df_tmp = pd.DataFrame({\"score\": 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# 4. 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    sample = bid[0]\n    mask = (bid == sample)\n\n    plt.figure(figsize=(12,5))\n    plt.plot(scores[mask], label=\"Predicted\")\n    plt.scatter(np.where(y[mask]==1), scores[mask][y[mask]==1],\n                color=\"red\", label=\"Actual Anomaly\")\n    plt.legend()\n    plt.title(f\"Building {sample}\")\n    plt.show()\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_isolation_forest(X, y, bid)\nevaluate(scores, y, bid)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-03T23:39:06.45258Z","iopub.execute_input":"2026-03-03T23:39:06.453704Z","iopub.status.idle":"2026-03-03T23:41:27.833184Z","shell.execute_reply.started":"2026-03-03T23:39:06.45367Z","shell.execute_reply":"2026-03-03T23:41:27.832365Z"}},"outputs":[],"execution_count":null}]}