{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","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":31254,"databundleVersionId":3103714}],"dockerImageVersionId":31329,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Predictive Modeling for Micro-Trend/Fad Identification in Fast Fashion\n## Full Pipeline — PySpark + Big Data Architecture\n\n> **Why PySpark?**  \n> transactions_train.csv มี **31.7M rows** — Pandas จะ OOM บน Kaggle (16GB RAM)  \n> PySpark ใช้ lazy evaluation + distributed processing ทำให้รองรับ data ขนาดนี้ได้\n\n| Phase | Tool | Description |\n|-------|------|-------------|\n| Phase 0 | PySpark | EDA |\n| Phase 1 | PySpark | Data Prep + Stockout Masking |\n| Phase 2 | PySpark + ruptures | Hybrid Labeling (PELT + k-Shape) |\n| Phase 3 | PySpark SQL | Feature Engineering |\n| Phase 4 | Pandas (small) | Feature Selection |\n| Phase 5 | LightGBM + XGBoost + LSTM | Stacked Ensemble |\n| Phase 6 | Pandas | Threshold Optimization |\n| Phase 7 | SHAP | Explainability |","metadata":{}},{"cell_type":"code","source":"import os\nos.environ['JAVA_HOME'] = '/usr/lib/jvm/java-17-openjdk-amd64'\nos.environ['PYSPARK_PYTHON'] = 'python3'\nprint(\"JAVA_HOME set!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T07:58:18.951494Z","iopub.execute_input":"2026-04-29T07:58:18.951745Z","iopub.status.idle":"2026-04-29T07:58:18.959574Z","shell.execute_reply.started":"2026-04-29T07:58:18.951721Z","shell.execute_reply":"2026-04-29T07:58:18.958797Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 1: Install & Setup ─────────────────────────────────\n!pip install pyspark lightgbm xgboost shap ruptures tslearn imbalanced-learn -q\n\nimport os, warnings\nwarnings.filterwarnings('ignore')\nos.environ['JAVA_HOME'] = '/usr/lib/jvm/java-17-openjdk-amd64'\n\n# PySpark\nfrom pyspark.sql import SparkSession\nfrom pyspark.sql import functions as F\nfrom pyspark.sql.window import Window\nfrom pyspark.sql.types import *\nfrom pyspark.ml.feature import StringIndexer, OneHotEncoder, VectorAssembler, MinMaxScaler\nfrom pyspark.ml import Pipeline\n\n# Standard\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# ML\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.feature_selection import mutual_info_classif\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.calibration import CalibratedClassifierCV\nfrom sklearn.metrics import (\n    f1_score, precision_score, recall_score,\n    precision_recall_curve, average_precision_score,\n    confusion_matrix, classification_report, roc_auc_score\n)\nimport lightgbm as lgb\nimport xgboost as xgb\nimport shap\nimport ruptures as rpt\nfrom imblearn.over_sampling import SMOTE\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import LSTM, Dense, Dropout\nfrom tensorflow.keras.callbacks import EarlyStopping\n\nplt.style.use('seaborn-v0_8-whitegrid')\nsns.set_palette('husl')\nSEED = 42\nnp.random.seed(SEED)\nprint('Libraries loaded!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T07:58:18.961058Z","iopub.execute_input":"2026-04-29T07:58:18.961299Z","iopub.status.idle":"2026-04-29T07:59:05.458649Z","shell.execute_reply.started":"2026-04-29T07:58:18.961278Z","shell.execute_reply":"2026-04-29T07:59:05.457903Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pyspark.sql import SparkSession\n\nspark = (\n    SparkSession.builder\n    .appName('HM_Fad_Detection')\n    .master('local[*]')\n    .config('spark.driver.memory', '8g')\n    .config('spark.sql.shuffle.partitions', '8')\n    .config('spark.sql.adaptive.enabled', 'true')\n    .config('spark.sql.adaptive.coalescePartitions.enabled', 'true')\n    .getOrCreate()\n)\nspark.sparkContext.setLogLevel('ERROR')\nprint(f'Spark {spark.version} ready!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T07:59:05.459756Z","iopub.execute_input":"2026-04-29T07:59:05.460567Z","iopub.status.idle":"2026-04-29T07:59:14.025829Z","shell.execute_reply.started":"2026-04-29T07:59:05.460540Z","shell.execute_reply":"2026-04-29T07:59:14.025003Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 3: Load Data with PySpark ─────────────────────────\n# ✅ PySpark reads CSV lazily — no memory spike at load time\nBASE = '/kaggle/input/competitions/h-and-m-personalized-fashion-recommendations/'\n\ntxn = (\n    spark.read.csv(BASE + 'transactions_train.csv', header=True, inferSchema=True)\n    .withColumn('t_dat', F.to_date('t_dat', 'yyyy-MM-dd'))\n    .withColumn('article_id',  F.col('article_id').cast('string'))\n    .withColumn('customer_id', F.col('customer_id').cast('string'))\n)\n\nart = spark.read.csv(BASE + 'articles.csv',  header=True, inferSchema=True)\ncus = spark.read.csv(BASE + 'customers.csv', header=True, inferSchema=True)\n\n# Cache articles (small, used many times)\nart.cache()\n\nprint('Schema: transactions')\ntxn.printSchema()\nprint(f'\\nDate range: {txn.agg(F.min(\"t_dat\"), F.max(\"t_dat\")).collect()[0]}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T07:59:14.027277Z","iopub.execute_input":"2026-04-29T07:59:14.027682Z","iopub.status.idle":"2026-04-29T08:00:29.533918Z","shell.execute_reply.started":"2026-04-29T07:59:14.027649Z","shell.execute_reply":"2026-04-29T08:00:29.532937Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 4: Row Counts (Spark action — triggers execution) ──\n# NOTE: count() is the ONE action needed. Everything before this is lazy.\nn_txn = txn.count()\nn_art = art.count()\nn_cus = cus.count()\n\nprint(f'Transactions : {n_txn:>12,}')\nprint(f'Articles     : {n_art:>12,}')\nprint(f'Customers    : {n_cus:>12,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:00:29.535151Z","iopub.execute_input":"2026-04-29T08:00:29.535501Z","iopub.status.idle":"2026-04-29T08:00:38.252456Z","shell.execute_reply.started":"2026-04-29T08:00:29.535462Z","shell.execute_reply":"2026-04-29T08:00:38.251678Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 0 — EDA with PySpark\nทำ EDA บน data เต็ม ","metadata":{}},{"cell_type":"code","source":"# ── Cell 5: Missing Value Analysis (Spark) ──────────────────\n# PySpark: count nulls per column in one pass over the data\ndef spark_missing_report(df, name):\n    # Combine count(*) with null counts in ONE pass\n    agg_exprs = [F.count(F.lit(1)).alias('__total__')] + [\n        F.sum(F.col(c).isNull().cast('int')).alias(c) for c in df.columns\n    ]\n    result = df.agg(*agg_exprs).collect()[0].asDict()\n    total = result.pop('__total__')\n    \n    rows = [\n        {\n            'column': col, \n            'missing': cnt, \n            'pct': round(cnt/total*100, 2),\n            'result': '🔴 MNAR' if cnt/total > 0.30 \n                    else '🟡 MAR' if cnt/total > 0.05 \n                    else '🟢 MCAR'\n        }\n        for col, cnt in result.items() if cnt > 0\n    ]\n    \n    print(f'\\n=== Missing Values: {name} (total rows: {total:,}) ===')\n    if not rows:\n        print('No missing values'); return pd.DataFrame()\n    \n    report = pd.DataFrame(rows).sort_values('pct', ascending=False)\n    print(report.to_string(index=False))\n    return report\n\nspark_missing_report(art, 'articles')\nspark_missing_report(cus, 'customers')\n# transactions: only check key columns\nspark_missing_report(txn.select('t_dat','article_id','customer_id','price'), 'transactions (key cols)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:00:38.253620Z","iopub.execute_input":"2026-04-29T08:00:38.253943Z","iopub.status.idle":"2026-04-29T08:01:19.155837Z","shell.execute_reply.started":"2026-04-29T08:00:38.253920Z","shell.execute_reply":"2026-04-29T08:01:19.154844Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 6: Cardinality Check (Spark) ──────────────────────\ncat_cols = [\n    'product_type_name', 'product_group_name',\n    'graphical_appearance_name', 'colour_group_name',\n    'department_name', 'index_name', 'garment_group_name'\n]\ncard_exprs = [F.countDistinct(c).alias(c) for c in cat_cols]\ncardinality = art.agg(*card_exprs).collect()[0].asDict()\ncard_df = pd.DataFrame(cardinality.items(), columns=['feature','unique_values'])\ncard_df = card_df.sort_values('unique_values', ascending=False)\n\nprint('=== Cardinality (articles) ===')\nprint(card_df.to_string(index=False))\n\nfig, ax = plt.subplots(figsize=(9, 5))\ncard_df.plot(kind='barh', x='feature', y='unique_values', ax=ax, color='steelblue', legend=False)\nax.axvline(50, color='red', linestyle='--', label='One-Hot limit (50)')\nax.set_title('Cardinality — Articles Categorical Features', fontweight='bold')\nax.legend()\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:01:19.159542Z","iopub.execute_input":"2026-04-29T08:01:19.159787Z","iopub.status.idle":"2026-04-29T08:01:20.942086Z","shell.execute_reply.started":"2026-04-29T08:01:19.159764Z","shell.execute_reply":"2026-04-29T08:01:20.941447Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 7: Sales Distribution + Cold-Start Analysis ──\n!pip install powerlaw -q\nimport powerlaw\n\n# 1. Identify \"Dead-on-Arrival\" articles (CRITICAL — missed in original)\nall_articles = art.select('article_id').distinct()\nsold_articles = txn.select('article_id').distinct()\n\ndead_on_arrival = all_articles.subtract(sold_articles)\nn_dead = dead_on_arrival.count()\nn_total_art = all_articles.count()\n\nprint(f'\\n Dead-on-Arrival Articles: {n_dead:,} / {n_total_art:,} ({n_dead/n_total_art*100:.2f}%)')\nprint('   → These are MAX-RISK candidates. Investigate before Phase 1.')\n\n# Profile dead-on-arrival items by category (do they cluster in specific product types?)\ndead_profile = (\n    dead_on_arrival.join(art, 'article_id')\n    .groupBy('product_group_name')\n    .count()\n    .orderBy(F.desc('count'))\n)\nprint('\\nDead-on-Arrival profile by product_group:')\ndead_profile.show(10, truncate=False)\n\n# 2. Sales distribution — already cached\nsales_per_art = (\n    txn.groupBy('article_id')\n    .count()\n    .withColumnRenamed('count', 'total_sales')\n    .orderBy(F.desc('total_sales'))\n)\nsales_per_art.cache()\n\nsales_per_art.select('total_sales').describe().show()\nsales_pd = sales_per_art.toPandas()\n\n# 3. RIGOROUS Power Law fit (Clauset et al. method)\nfit = powerlaw.Fit(sales_pd['total_sales'].values, discrete=True, verbose=False)\nprint(f'\\n=== Rigorous Power Law Fit ===')\nprint(f'  Estimated alpha (α)  : {fit.alpha:.3f}')\nprint(f'  Estimated x_min      : {fit.xmin:.0f}')\nprint(f'  → Power Law applies for sales >= {fit.xmin:.0f}')\n\n# Compare against alternative distributions\nfor alt in ['lognormal', 'exponential', 'stretched_exponential']:\n    R, p = fit.distribution_compare('power_law', alt, normalized_ratio=True)\n    verdict = '✓ PL better' if R > 0 and p < 0.05 else '✗ PL not better'\n    print(f'  PL vs {alt:25s}: R={R:+.3f}, p={p:.4f}  → {verdict}')\n\n# 4. Quantile breakdown for stratified sampling (informs Phase 1)\nprint('\\n=== Sales Quantiles (for Stratified Sampling in Phase 1) ===')\nquantiles = sales_per_art.approxQuantile('total_sales', [0.5, 0.75, 0.90, 0.95, 0.99], 0.01)\nlabels = ['p50', 'p75', 'p90', 'p95', 'p99']\nfor l, q in zip(labels, quantiles):\n    print(f'  {l}: {q:>8.0f} sales')\n\n# 5. Visualization — log-log + powerlaw fit overlay\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Left: rank-frequency log-log\nsales_pd['rank'] = range(1, len(sales_pd) + 1)\naxes[0].loglog(sales_pd['rank'], sales_pd['total_sales'], 'o', alpha=0.3, ms=2)\naxes[0].set_xlabel('Rank (log)'); axes[0].set_ylabel('Sales (log)')\naxes[0].set_title(f'Rank-Frequency (slope ≈ -{fit.alpha:.2f})')\naxes[0].grid(True, alpha=0.3)\n\n# Right: CCDF with powerlaw fit\nfit.plot_ccdf(ax=axes[1], color='steelblue', linewidth=2, label='Empirical CCDF')\nfit.power_law.plot_ccdf(ax=axes[1], color='red', linestyle='--', label=f'PL fit (α={fit.alpha:.2f})')\nfit.lognormal.plot_ccdf(ax=axes[1], color='green', linestyle=':', label='Lognormal fit')\naxes[1].set_xlabel('Sales (log)'); axes[1].set_ylabel('P(X ≥ x)')\naxes[1].set_title('CCDF Comparison: Power Law vs Lognormal')\naxes[1].legend()\n\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:01:20.942887Z","iopub.execute_input":"2026-04-29T08:01:20.943238Z","iopub.status.idle":"2026-04-29T08:02:56.360830Z","shell.execute_reply.started":"2026-04-29T08:01:20.943201Z","shell.execute_reply":"2026-04-29T08:02:56.359959Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 1 — Data Preparation with PySpark","metadata":{}},{"cell_type":"code","source":"# ═══════════════════════════════════════════════════════════════════\n#  LOCKED-IN CONSTANTS FROM PHASE 0 EDA\n#  อ้างอิงจาก Phase 0.1 (Missing Values) และ 0.2 (Power Law)\n# ═══════════════════════════════════════════════════════════════════\n\n# From Phase 0.2 — Lognormal + Power Law tail analysis\nPOWER_LAW_XMIN    = 1821     # regime boundary (top 3% = viral tail)\nPOWER_LAW_ALPHA   = 3.00     # tail exponent (edge of finite variance)\n\n# From Phase 0.1 — Dead-on-Arrival handling\nN_DEAD_ON_ARRIVAL = 995      # articles with 0 transactions (exclude from labeling)\n\n# Phase 1 configuration (informed by methodology + EDA)\nMIN_ACTIVE_WEEKS  = 3        # ✅ lowered from 6 to keep true Fads (short-lived = Fad signal)\nCENSORING_WEEKS   = 26       # articles launched < 26 weeks before dataset_end → censored\nMAX_WEEKS_PIVOT   = 26       # pivot width for time-series analysis\n\n# Stockout detection parameters\nSTOCKOUT_BASELINE_WINDOW = 4   # rolling window size (weeks)\nSTOCKOUT_DROP_THRESHOLD  = 0.30  # sales < 30% of baseline\nSTOCKOUT_RECOVERY_THRESH = 0.60  # recovery > 60% of baseline within 2 weeks\nSTOCKOUT_MIN_BASELINE    = 2     # baseline must have meaningful volume\n\nprint('Phase 0 constants loaded')\nprint(f'   x_min = {POWER_LAW_XMIN}, α = {POWER_LAW_ALPHA}')\nprint(f'   Min active weeks = {MIN_ACTIVE_WEEKS} (Fad-preserving)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:02:56.362149Z","iopub.execute_input":"2026-04-29T08:02:56.362515Z","iopub.status.idle":"2026-04-29T08:02:56.369223Z","shell.execute_reply.started":"2026-04-29T08:02:56.362478Z","shell.execute_reply":"2026-04-29T08:02:56.368304Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 8: Weekly Aggregation with DENSE time series ──\n\n# Get dataset boundary for censoring decision\ndataset_start, dataset_end = txn.agg(\n    F.min('t_dat').alias('start'),\n    F.max('t_dat').alias('end')\n).collect()[0]\nprint(f'Dataset period: {dataset_start} → {dataset_end}')\n\n# ─── Step 1: Aggregate to SPARSE weekly (only weeks with sales) ───\nweekly_sparse = (\n    txn\n    .withColumn('week_start', F.date_trunc('week', F.col('t_dat')))  # Monday as week start\n    .groupBy('article_id', 'week_start')\n    .agg(\n        F.count(F.lit(1))                  .alias('weekly_sales'),\n        F.countDistinct('customer_id')     .alias('unique_buyers'),\n        F.countDistinct('sales_channel_id').alias('n_channels')\n    )\n)\n\n# ─── Step 2: Compute launch_week & last_seen_week per article ───\nlifecycle = (\n    weekly_sparse\n    .groupBy('article_id')\n    .agg(\n        F.min('week_start').alias('launch_week'),\n        F.max('week_start').alias('last_seen_week')\n    )\n    .withColumn(\n        'weeks_from_launch_to_dataset_end',\n        (F.datediff(F.lit(dataset_end), F.col('launch_week')) / 7).cast('int')\n    )\n    # ✅ Tag right-censored: launched too close to dataset end to observe full lifecycle\n    .withColumn(\n        'is_censored',\n        (F.col('weeks_from_launch_to_dataset_end') < CENSORING_WEEKS).cast('int')\n    )\n)\n\n# Quick stats\nn_total = lifecycle.count()\nn_censored = lifecycle.filter(F.col('is_censored') == 1).count()\nprint(f'Total articles with transactions: {n_total:,}')\nprint(f'  → Right-censored (launched < {CENSORING_WEEKS}w before end): {n_censored:,} ({n_censored/n_total*100:.1f}%)')\nprint(f'  → Usable for training: {n_total - n_censored:,}')\n\n# ─── Step 3: Build DENSE skeleton (every week from launch → last_seen) ───\n# ✅ Uses sequence() to explode all weeks — enables proper lag/lead operations\ndense_skeleton = (\n    lifecycle\n    .filter(F.col('is_censored') == 0)  # exclude censored from skeleton\n    .select(\n        'article_id',\n        'launch_week',\n        F.explode(\n            F.sequence(\n                F.col('launch_week'),\n                F.col('last_seen_week'),\n                F.expr('INTERVAL 7 DAYS')\n            )\n        ).alias('week_start')\n    )\n)\n\n# ─── Step 4: Join sparse data onto dense skeleton → zero-fill gaps ───\nweekly = (\n    dense_skeleton\n    .join(weekly_sparse, on=['article_id', 'week_start'], how='left')\n    .fillna(0, subset=['weekly_sales', 'unique_buyers', 'n_channels'])\n    .withColumn(\n        'weeks_since_launch',\n        (F.datediff('week_start', 'launch_week') / 7).cast('int')\n    )\n)\n\n# Cache the dense weekly DataFrame (used many times in Phase 1-3)\nweekly.cache()\nn_records = weekly.count()  # trigger cache\n\nprint(f'\\n✅ Dense weekly records: {n_records:,}')\nprint(f'   Expansion ratio: {n_records / weekly_sparse.count():.2f}x '\n      f'(sparse → dense; higher = more stockout/dormant weeks)')\n\n# Preview\nprint('\\n--- Sample (article 110065002) ---')\nweekly.filter(F.col('article_id') == '110065002') \\\n      .select('week_start', 'weeks_since_launch', 'weekly_sales', 'unique_buyers') \\\n      .orderBy('weeks_since_launch').show(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:02:56.370034Z","iopub.execute_input":"2026-04-29T08:02:56.370204Z","iopub.status.idle":"2026-04-29T08:06:54.890008Z","shell.execute_reply.started":"2026-04-29T08:02:56.370187Z","shell.execute_reply":"2026-04-29T08:06:54.889303Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 9: Multi-Tier Stockout Detection ──\n\n# ─── Updated Thresholds (relaxed) ───\nSTOCKOUT_BASELINE_WINDOW = 4         # unchanged\nSTOCKOUT_MIN_BASELINE    = 3.0       # slight bump — reject noise แต่ไม่รัด\nDROP_SEVERE   = 0.20   # drop ≤ 20% of baseline (ช่วย severe)\nDROP_MODERATE = 0.50   # drop ≤ 50% of baseline (ช่วย moderate)\nRECOVERY_STRONG  = 0.70   # recover ≥ 70% within 2 weeks\nRECOVERY_PARTIAL = 0.40   # recover ≥ 40% within 3 weeks\n\nw_art      = Window.partitionBy('article_id').orderBy('weeks_since_launch')\nw_baseline = w_art.rowsBetween(-STOCKOUT_BASELINE_WINDOW, -1)\n\nweekly = (\n    weekly\n    # 1. Rolling baseline (past 4 weeks)\n    .withColumn('baseline_sales', F.avg('weekly_sales').over(w_baseline))\n    \n    # 2. Future recovery window (extend to 3 weeks — some stockouts recover slower)\n    .withColumn('sales_lead1', F.lead('weekly_sales', 1).over(w_art))\n    .withColumn('sales_lead2', F.lead('weekly_sales', 2).over(w_art))\n    .withColumn('sales_lead3', F.lead('weekly_sales', 3).over(w_art))\n    \n    # 3. Ratios\n    .withColumn('drop_ratio',\n        F.col('weekly_sales') / (F.col('baseline_sales') + 1.0))\n    \n    # Best recovery within 3 weeks (max lead)\n    .withColumn('future_recovery',\n        F.greatest(F.col('sales_lead1'), F.col('sales_lead2'), F.col('sales_lead3')) /\n        (F.col('baseline_sales') + 1.0))\n    \n    # 4. Tiered stockout detection (3 levels of severity)\n    .withColumn('stockout_tier',\n        F.when(\n            # Tier 3 (SEVERE): ดรอปฮวบ 80%+ AND recover strongly\n            (F.col('drop_ratio')      < DROP_SEVERE) &\n            (F.col('future_recovery') > RECOVERY_STRONG) &\n            (F.col('baseline_sales')  >= STOCKOUT_MIN_BASELINE),\n            3\n        ).when(\n            # Tier 2 (MODERATE): ดรอป 50%+ AND recover partially\n            (F.col('drop_ratio')      < DROP_MODERATE) &\n            (F.col('future_recovery') > RECOVERY_PARTIAL) &\n            (F.col('baseline_sales')  >= STOCKOUT_MIN_BASELINE),\n            2\n        ).when(\n            # Tier 1 (MILD): ตกกะทันหัน > 50% หลังช่วงขายสูง (baseline ≥ 5)\n            #   → ไม่ต้องรอ recovery เพราะอาจเป็น end-of-life stockout\n            (F.col('drop_ratio')     < DROP_MODERATE) &\n            (F.col('baseline_sales') >= 5.0),\n            1\n        ).otherwise(0))\n    \n    # 5. Binary flag (สำหรับ downstream ที่ต้องการ 0/1 เดิม)\n    .withColumn('is_stockout',\n        (F.col('stockout_tier') > 0).cast('int'))\n    \n    # 6. Impute value ตาม tier\n    .withColumn('weekly_sales_clean',\n        F.when(F.col('stockout_tier') == 3, F.col('baseline_sales'))\n         .when(F.col('stockout_tier') == 2, F.col('baseline_sales') * 0.7)  # partial impute\n         .when(F.col('stockout_tier') == 1, F.col('baseline_sales') * 0.5)  # conservative\n         .otherwise(F.col('weekly_sales')))\n)\n\n# Re-cache\nweekly.unpersist()\nweekly.cache()\n_ = weekly.count()\n\n# ─── Diagnostic ───\ntotal_weeks  = weekly.count()\ntier_dist = (\n    weekly.groupBy('stockout_tier')\n    .agg(F.count('*').alias('count'))\n    .orderBy('stockout_tier')\n).toPandas()\ntier_dist['pct'] = tier_dist['count'] / total_weeks * 100\n\nprint('=== Stockout Tier Distribution ===')\nprint(tier_dist.to_string(index=False))\n\ntotal_stockout_rate = tier_dist.loc[tier_dist['stockout_tier'] > 0, 'pct'].sum()\nprint(f'\\nTotal stockout rate (all tiers): {total_stockout_rate:.2f}%')\nprint(f'   Expected: 5-15% in fast fashion')\n\nif   total_stockout_rate < 4:  print('   ⚠️  ยังต่ำ — พิจารณา relax อีก')\nelif total_stockout_rate > 20: print('   ⚠️  สูงเกิน — tighten thresholds')\nelse:                          print('   ✅ Within expected range')\n\n# Sample ของแต่ละ tier\nprint('\\n--- Samples by Tier ---')\nfor tier in [3, 2, 1]:\n    print(f'\\nTier {tier}:')\n    weekly.filter(F.col('stockout_tier') == tier) \\\n          .select('article_id', 'weeks_since_launch', 'baseline_sales',\n                  'weekly_sales', 'sales_lead1', 'sales_lead2', \n                  'drop_ratio', 'future_recovery', 'weekly_sales_clean') \\\n          .show(3, truncate=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:06:54.892119Z","iopub.execute_input":"2026-04-29T08:06:54.892419Z","iopub.status.idle":"2026-04-29T08:07:04.419284Z","shell.execute_reply.started":"2026-04-29T08:06:54.892374Z","shell.execute_reply":"2026-04-29T08:07:04.418468Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 10: Filter + Normalize + Regime-Aware Stratification ──\n# ─── Step 1: Compute active weeks per article (non-zero sales) ───\nactive_weeks = (\n    weekly\n    .filter(F.col('weekly_sales_clean') > 0)\n    .groupBy('article_id')\n    .agg(F.count('*').alias('active_weeks'))\n)\n\n# ─── Step 2: Filter: >= 3 active weeks (NOT lifespan!) ───\n# ✅ ทำไม >= 3: ต้องมี data พอสำหรับ early signal features (weeks 1-3)\n#    ถ้า < 3 = ขายแค่ 1-2 ครั้งแล้วหาย = เป็น \"dud\" ไม่ใช่ \"Fad\"\nvalid_arts = active_weeks.filter(F.col('active_weeks') >= MIN_ACTIVE_WEEKS)\n\nn_before = weekly.select('article_id').distinct().count()\nn_after  = valid_arts.count()\n\nweekly = weekly.join(valid_arts.select('article_id'), on='article_id', how='inner')\nprint(f'Articles before filter: {n_before:,}')\nprint(f'Articles after filter (>= {MIN_ACTIVE_WEEKS} active weeks): {n_after:,}')\nprint(f'   Excluded: {n_before - n_after:,} articles with <{MIN_ACTIVE_WEEKS} active weeks (duds)')\n\n# ─── Step 3: Min-Max normalize per article (with edge-case handling) ───\nw_norm = Window.partitionBy('article_id')\nweekly = (\n    weekly\n    .withColumn('art_min', F.min('weekly_sales_clean').over(w_norm))\n    .withColumn('art_max', F.max('weekly_sales_clean').over(w_norm))\n    .withColumn('norm_sales',\n        # ✅ ป้องกัน division by zero เมื่อ article มี constant sales\n        F.when(F.col('art_max') == F.col('art_min'), F.lit(0.0))\n         .otherwise(\n             (F.col('weekly_sales_clean') - F.col('art_min')) /\n             (F.col('art_max') - F.col('art_min'))\n         ))\n)\n\n# ─── Step 4: Regime-aware stratification (NEW — from Phase 0.2) ───\n# ✅ แทนที่ uniform ntile(4) ด้วย \"body + tail\" regime\n#    Body (97%): sales < 1821 → log-quantile stratify\n#    Tail (3%):  sales >= 1821 → ทุก article อยู่ใน tail (rare & critical)\ntotal_sales_df = (\n    weekly.groupBy('article_id')\n    .agg(F.sum('weekly_sales_clean').alias('total_sales'))\n    # Regime tag (from Phase 0.2 Lognormal + Power Law tail analysis)\n    .withColumn('regime',\n        F.when(F.col('total_sales') >= POWER_LAW_XMIN, 'tail')\n         .otherwise('body'))\n)\n\n# Log-quantile stratification WITHIN each regime\n# (body: 4 strata by log-sales; tail: 2 strata because data is sparse)\nw_body = Window.partitionBy('regime').orderBy(F.log1p('total_sales'))\ntotal_sales_df = (\n    total_sales_df\n    .withColumn(\n        'strat_bucket',\n        F.when(F.col('regime') == 'tail', F.ntile(2).over(w_body))\n         .otherwise(F.ntile(4).over(w_body)))\n    .withColumn(\n        'strat_key',\n        F.concat(F.col('regime'), F.lit('_Q'), F.col('strat_bucket')))\n)\n\n# Regime distribution report\nprint('\\n=== Regime Distribution (from Power Law EDA) ===')\nregime_dist = total_sales_df.groupBy('regime').agg(\n    F.count('*').alias('n_articles'),\n    F.min('total_sales').alias('min_sales'),\n    F.max('total_sales').alias('max_sales'),\n    F.mean('total_sales').alias('mean_sales')\n).toPandas()\nprint(regime_dist.to_string(index=False))\n\nprint('\\n=== Stratum Key Distribution ===')\ntotal_sales_df.groupBy('strat_key').count().orderBy('strat_key').show()\n\n# Join stratum back to weekly\nweekly = weekly.join(\n    total_sales_df.select('article_id', 'regime', 'strat_key', 'total_sales'),\n    on='article_id', how='left'\n)\n\n# Re-cache final weekly DataFrame\nweekly.unpersist()\nweekly.cache()\nn_valid = weekly.select('article_id').distinct().count()\nprint(f'\\n✅ Final valid articles: {n_valid:,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:07:04.420538Z","iopub.execute_input":"2026-04-29T08:07:04.421121Z","iopub.status.idle":"2026-04-29T08:07:24.279240Z","shell.execute_reply.started":"2026-04-29T08:07:04.421086Z","shell.execute_reply":"2026-04-29T08:07:24.278509Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 11: Sales Pivot with Validation ──\n\npivot_spark = (\n    weekly\n    .filter(\n        (F.col('weeks_since_launch') >= 0) & \n        (F.col('weeks_since_launch') < MAX_WEEKS_PIVOT)\n    )\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(MAX_WEEKS_PIVOT)))\n    .agg(F.first('norm_sales'))\n    .fillna(0)\n)\n\n# Rename columns: '0' → 'w0', etc.\nnew_cols = ['article_id'] + [f'w{i}' for i in range(MAX_WEEKS_PIVOT)]\npivot_spark = pivot_spark.toDF(*new_cols)\n\n# Pre-collect sanity check (validate before bringing to Pandas)\nprint('=== Pre-Collect Sanity Check ===')\npivot_spark.select(\n    F.count('*').alias('n_articles'),\n    F.mean('w0').alias('mean_w0'),\n    F.mean('w3').alias('mean_w3'),\n    F.mean('w12').alias('mean_w12'),\n    F.mean('w25').alias('mean_w25')\n).show()\nprint('↑ Expect: mean_w0 > mean_w3 > mean_w12 > mean_w25 '\n      '(peak at launch, decay over time)')\n\n# Collect to Pandas (89K × 26 = ~18MB, safe)\npivot = pivot_spark.toPandas().set_index('article_id')\nprint(f'\\n✅ Sales pivot shape: {pivot.shape}  (articles × weeks)')\n\n# ─── Diagnostic 1: Death-week distribution (informs Phase 2 labeling) ───\n# \"Death week\" = last week with non-zero sales → short death_week = Fad-like\nweek_cols = [f'w{i}' for i in range(MAX_WEEKS_PIVOT)]\ndeath_week = (pivot[week_cols] > 0).iloc[:, ::-1].cummax(axis=1).iloc[:, ::-1].sum(axis=1)\n\nprint('\\n=== Death Week Distribution ===')\nprint(f'(last week with sales > 0; short death_week = Fad candidate)')\nprint(death_week.describe().to_string())\nprint(f'\\n  Articles dying by week 6 (short Fad):  '\n      f'{(death_week <= 6).sum():,} ({(death_week <= 6).mean()*100:.1f}%)')\nprint(f'  Articles dying by week 12 (medium):    '\n      f'{((death_week > 6) & (death_week <= 12)).sum():,} '\n      f'({((death_week > 6) & (death_week <= 12)).mean()*100:.1f}%)')\nprint(f'  Articles still alive at week 25:       '\n      f'{(death_week > 25).sum():,} ({(death_week > 25).mean()*100:.1f}%)')\n\n# ─── Diagnostic 2: Visualize curve shapes by regime ───\nregime_df = total_sales_df.select('article_id', 'regime').toPandas().set_index('article_id')\npivot_with_regime = pivot.join(regime_df, how='left')\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4))\n\n# Sample 30 articles per regime\nfor ax, regime, color in zip(axes, ['body', 'tail'], ['steelblue', 'crimson']):\n    sample = pivot_with_regime[pivot_with_regime['regime'] == regime].sample(\n        min(30, (pivot_with_regime['regime'] == regime).sum()), random_state=42\n    )\n    for idx in sample.index:\n        ax.plot(range(MAX_WEEKS_PIVOT), sample.loc[idx, week_cols].values, \n                alpha=0.3, color=color, lw=1)\n    ax.set_title(f'{regime.upper()} regime — 30 random articles', fontweight='bold')\n    ax.set_xlabel('Weeks since launch'); ax.set_ylabel('Normalized sales')\n    ax.grid(alpha=0.3)\n\nplt.tight_layout(); plt.show()\n\nprint('\\n✅ Phase 1 complete — ready for Phase 2 (Hybrid Labeling)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:07:24.280272Z","iopub.execute_input":"2026-04-29T08:07:24.280585Z","iopub.status.idle":"2026-04-29T08:07:32.199387Z","shell.execute_reply.started":"2026-04-29T08:07:24.280553Z","shell.execute_reply":"2026-04-29T08:07:32.198518Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 2 — Hybrid Labeling (PELT + k-Mean)","metadata":{}},{"cell_type":"code","source":"# ── Cell 12: PELT with Phase 1 Insights ──\n\nFAD_MAX_PEAK_WEEK     = 6\nFAD_MIN_DROP_RATIO    = 0.70   # tighten 0.60 → 0.70\nFAD_MIN_PEAK_RATIO    = 4.0    # tighten 2.5 → 4.0 (biggest lever)\nFAD_LOW_TAIL_RATIO    = 0.20   # tighten 0.30 → 0.20\nPELT_PENALTY          = 2\n\npelt_results = []\nfor art_id, row in pivot.iterrows():\n    series = row.values.astype(float)\n    \n    nonzero = np.nonzero(series)[0]\n    if len(nonzero) == 0:\n        pelt_results.append({\n            'article_id': art_id, 'pelt_fad': 0, 'n_breakpoints': 0,\n            'peak_week': 0, 'death_week': 0,\n            'peak_ratio': 0.0, 'drop_ratio': 0.0, 'tail_ratio': 0.0,\n        })\n        continue\n    \n    death_week = int(nonzero[-1])\n    series_trimmed = series[:death_week + 1]\n    \n    if len(series_trimmed) < 5:\n        pelt_results.append({\n            'article_id': art_id, 'pelt_fad': 0, 'n_breakpoints': 0,\n            'peak_week': 0, 'death_week': death_week,\n            'peak_ratio': 0.0, 'drop_ratio': 0.0, 'tail_ratio': 0.0,\n        })\n        continue\n    \n    try:\n        algo = rpt.Pelt(model='rbf', min_size=2).fit(series_trimmed)\n        bkps = algo.predict(pen=PELT_PENALTY)\n        n_bkps = len(bkps) - 1\n    except Exception:\n        n_bkps = 0\n    \n    peak_week = int(np.argmax(series_trimmed))\n    peak_value = series_trimmed[peak_week]\n    \n    if peak_week < len(series_trimmed) - 1:\n        post_peak_mean = series_trimmed[peak_week + 1:].mean() + 1e-8\n        peak_ratio = peak_value / post_peak_mean\n        tail_ratio = post_peak_mean / (peak_value + 1e-8)\n    else:\n        peak_ratio = 0.0\n        tail_ratio = 1.0\n    \n    drop_ratio = (peak_value - series_trimmed[-1]) / (peak_value + 1e-8)\n    \n    is_fad = (\n        (n_bkps >= 1) and\n        (peak_week <= FAD_MAX_PEAK_WEEK) and\n        (peak_ratio > FAD_MIN_PEAK_RATIO) and\n        (tail_ratio < FAD_LOW_TAIL_RATIO) and\n        (drop_ratio > FAD_MIN_DROP_RATIO)\n    )\n    \n    pelt_results.append({\n        'article_id': art_id,\n        'pelt_fad': int(is_fad),\n        'n_breakpoints': n_bkps,\n        'peak_week': peak_week,\n        'death_week': death_week,\n        'peak_ratio': round(peak_ratio, 2),\n        'drop_ratio': round(drop_ratio, 3),\n        'tail_ratio': round(tail_ratio, 3),\n    })\n\npelt_df = pd.DataFrame(pelt_results)\n\n# ─── Diagnostics ───\nprint('✅ PELT done.')\nprint(f'\\n{pelt_df[\"pelt_fad\"].value_counts()}')\nprint(f'\\nFad rate (PELT): {pelt_df[\"pelt_fad\"].mean():.2%}')\n\nprint('\\n=== Criteria Pass Rates (individual) ===')\nprint(f'  Has breakpoint (n≥1)   : {(pelt_df[\"n_breakpoints\"] >= 1).mean():.2%}')\nprint(f'  Peak week ≤ {FAD_MAX_PEAK_WEEK}           : {(pelt_df[\"peak_week\"] <= FAD_MAX_PEAK_WEEK).mean():.2%}')\nprint(f'  Peak ratio > {FAD_MIN_PEAK_RATIO}          : {(pelt_df[\"peak_ratio\"] > FAD_MIN_PEAK_RATIO).mean():.2%}')\nprint(f'  Tail ratio < {FAD_LOW_TAIL_RATIO}          : {(pelt_df[\"tail_ratio\"] < FAD_LOW_TAIL_RATIO).mean():.2%}')\nprint(f'  Drop ratio > {FAD_MIN_DROP_RATIO}          : {(pelt_df[\"drop_ratio\"] > FAD_MIN_DROP_RATIO).mean():.2%}')\n\npelt_fads = pelt_df[pelt_df['pelt_fad'] == 1]\nif len(pelt_fads) > 0:\n    print('\\n=== PELT-Fad Candidates: Temporal Profile ===')\n    print(f'Peak week   — mean: {pelt_fads[\"peak_week\"].mean():.1f}, '\n          f'median: {pelt_fads[\"peak_week\"].median():.0f}')\n    print(f'Peak ratio  — mean: {pelt_fads[\"peak_ratio\"].mean():.2f}')\n    print(f'Drop ratio  — mean: {pelt_fads[\"drop_ratio\"].mean():.2f}')\n    print(f'Tail ratio  — mean: {pelt_fads[\"tail_ratio\"].mean():.3f}')\n\n# Sanity check\nprint('\\n=== Sanity Check: Article 110065002 ===')\ncheck = pelt_df[pelt_df['article_id'] == '110065002']\nif len(check) > 0:\n    print(check.to_string(index=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:07:32.200508Z","iopub.execute_input":"2026-04-29T08:07:32.200731Z","iopub.status.idle":"2026-04-29T08:08:09.725725Z","shell.execute_reply.started":"2026-04-29T08:07:32.200709Z","shell.execute_reply":"2026-04-29T08:08:09.724980Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 13 Distributed Clustering ──\n#\n#  ใช้ Spark MLlib KMeans แทน tslearn k-Shape\n#  ✅ Distributed — ใช้ทุก node ใน cluster\n#  ✅ เร็วกว่า k-Shape 10-50x\n#  ⚠️ ใช้ Euclidean distance แทน SBD — อาจจับ shape ได้ด้อยกว่า k-Shape เล็กน้อย\n#     แต่ข้อมูลถูก normalize แล้ว (Phase 1 Min-Max) จึงลด gap ตรงนี้ได้มาก\n# ═══════════════════════════════════════════════════════════════════\n\nfrom pyspark.ml.feature import VectorAssembler\nfrom pyspark.ml.clustering import KMeans, BisectingKMeans\nfrom pyspark.ml.evaluation import ClusteringEvaluator\nimport time\n\n# ─── Config ───\nK_SEARCH_RANGE    = range(3, 9)\nFAD_MAX_PEAK_WEEK = 6\n\n# ═══════════════════════════════════════════════════════\n#  Step 1: Prepare data in Spark (distributed)\n# ═══════════════════════════════════════════════════════\n\nprint('📊 Preparing data in Spark...')\n\n# pivot อยู่ใน Pandas → แปลงกลับเป็น Spark\nweek_cols = [f'w{i}' for i in range(26)]\npivot_reset = pivot.reset_index()\npivot_spark = spark.createDataFrame(pivot_reset)\n\n# VectorAssembler: รวม w0-w25 เป็น feature vector เดียว\nassembler = VectorAssembler(inputCols=week_cols, outputCol='features')\ndata_spark = assembler.transform(pivot_spark).select('article_id', 'features')\ndata_spark.cache()\n\nprint(f'   Rows: {data_spark.count():,}')\nprint(f'   Feature vector length: 26')\n\n# ═══════════════════════════════════════════════════════\n#  Step 2: Elbow Method — หา optimal k (Spark — distributed)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n🔍 Searching optimal k over range {list(K_SEARCH_RANGE)}...')\nprint('   (ทุก k รันแบบ distributed บน Spark)')\n\nelbow_results = []\nfor k in K_SEARCH_RANGE:\n    t0 = time.time()\n    \n    kmeans = KMeans(\n        k=k,\n        featuresCol='features',\n        predictionCol='cluster',\n        seed=42,\n        maxIter=100,\n        initMode='k-means||',  # scalable k-means++ (distributed)\n    )\n    model = kmeans.fit(data_spark)\n    \n    # Compute inertia (WSSSE — Within Set Sum of Squared Errors)\n    inertia = model.summary.trainingCost\n    elapsed = time.time() - t0\n    \n    elbow_results.append({'k': k, 'inertia': inertia, 'time': elapsed})\n    print(f'   k={k}  inertia={inertia:.2f}  time={elapsed:.1f}s')\n\nelbow_df = pd.DataFrame(elbow_results)\n\n# ─── Kneedle-based Elbow Detection ───\nks_vals = elbow_df['k'].values.astype(float)\ninertia_vals = elbow_df['inertia'].values\n\nx_norm = (ks_vals - ks_vals.min()) / (ks_vals.max() - ks_vals.min())\ny_norm = (inertia_vals - inertia_vals.min()) / (inertia_vals.max() - inertia_vals.min())\n\nline_vec = np.array([x_norm[-1] - x_norm[0], y_norm[-1] - y_norm[0]])\nline_vec_norm = line_vec / np.linalg.norm(line_vec)\n\ndistances = []\nfor i in range(len(ks_vals)):\n    point_vec = np.array([x_norm[i] - x_norm[0], y_norm[i] - y_norm[0]])\n    proj = np.dot(point_vec, line_vec_norm) * line_vec_norm\n    perp = point_vec - proj\n    distances.append(np.linalg.norm(perp))\n\nBEST_K = int(ks_vals[np.argmax(distances)])\nprint(f'\\n📍 Elbow at k = {BEST_K}')\n\n# Elbow Plot\nfig, axes = plt.subplots(1, 2, figsize=(13, 4))\naxes[0].plot(ks_vals, inertia_vals, 'o-', lw=2, ms=8, color='steelblue')\naxes[0].axvline(BEST_K, linestyle='--', color='crimson', lw=2, label=f'Elbow k={BEST_K}')\naxes[0].set_xlabel('k'); axes[0].set_ylabel('Inertia (WSSSE)')\naxes[0].set_title('Elbow Method — Spark KMeans'); axes[0].legend(); axes[0].grid(alpha=0.3)\n\naxes[1].bar(ks_vals, distances, color='coral', edgecolor='darkred')\naxes[1].axvline(BEST_K, linestyle='--', color='crimson', lw=2)\naxes[1].set_xlabel('k'); axes[1].set_ylabel('Perpendicular distance')\naxes[1].set_title('Kneedle Distance'); axes[1].grid(alpha=0.3)\nplt.suptitle(f'Optimal k Search — Elbow at k={BEST_K}', fontweight='bold')\nplt.tight_layout(); plt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 3: Fit final KMeans with best k (Spark — distributed)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n🎯 Fitting final KMeans with k={BEST_K} (distributed)...')\nt0 = time.time()\n\nfinal_kmeans = KMeans(\n    k=BEST_K,\n    featuresCol='features',\n    predictionCol='cluster',\n    seed=42,\n    maxIter=300,\n    initMode='k-means||',\n)\nfinal_model = final_kmeans.fit(data_spark)\nclustered = final_model.transform(data_spark)\n\nprint(f'   Elapsed: {time.time() - t0:.1f}s')\n\n# Silhouette Score (distributed evaluation)\nevaluator = ClusteringEvaluator(\n    featuresCol='features', predictionCol='cluster',\n    metricName='silhouette'\n)\nsilhouette = evaluator.evaluate(clustered)\nprint(f'   Silhouette Score: {silhouette:.4f}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 4: Extract labels + Compute centroids\n# ═══════════════════════════════════════════════════════\n\n# Labels: Spark → Pandas\nkshape_labels_df = clustered.select('article_id', 'cluster').toPandas()\nkshape_labels_df['article_id'] = kshape_labels_df['article_id'].astype(str)\n\n# Compute centroids manually (mean of members per cluster)\n# ใช้ Spark aggregation — distributed\ncentroid_spark = (\n    pivot_spark\n    .join(\n        clustered.select(\n            F.col('article_id').cast('string').alias('article_id'), \n            'cluster'\n        ),\n        on='article_id'\n    )\n)\n\ncentroids = {}\nfor i in range(BEST_K):\n    cluster_data = centroid_spark.filter(F.col('cluster') == i)\n    centroid_values = cluster_data.select(\n        *[F.avg(c).alias(c) for c in week_cols]\n    ).collect()[0]\n    centroids[i] = np.array([centroid_values[c] for c in week_cols])\n\nkshape_labels = kshape_labels_df.set_index('article_id')['cluster'].reindex(pivot.index).values\n\nprint(f'   Labels shape: {kshape_labels.shape}')\nprint(f'   Centroids computed: {len(centroids)}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 5: Auto-Detect Fad Clusters (same logic as before)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n📊 Computing cluster metrics (k={BEST_K})...')\ncluster_metrics = []\nfor i in range(BEST_K):\n    centroid = centroids[i]\n    n_members = int((kshape_labels == i).sum())\n    \n    peak_w = int(np.argmax(centroid))\n    peak_v = centroid[peak_w]\n    end_v  = centroid[-1]\n    min_v  = centroid.min()\n    \n    prominence = peak_v - min_v\n    width_at_half = int((centroid > peak_v * 0.5).sum())\n    sharpness = prominence / max(1, width_at_half)\n    drop = (peak_v - end_v) / (abs(peak_v) + 1e-8)\n    is_early_peak = peak_w <= FAD_MAX_PEAK_WEEK\n    \n    cluster_metrics.append({\n        'cluster': i, 'n': n_members,\n        'peak_week': peak_w, 'peak_value': round(peak_v, 4),\n        'end_value': round(end_v, 4),\n        'prominence': round(prominence, 4),\n        'width_at_half': width_at_half,\n        'sharpness': round(sharpness, 4),\n        'drop': round(drop, 2),\n        'is_early_peak': is_early_peak,\n    })\n\ncm_df = pd.DataFrame(cluster_metrics).sort_values('sharpness', ascending=False)\nprint('\\n=== Cluster Metrics (sorted by sharpness) ===')\nprint(cm_df.to_string(index=False))\n\n# Auto-detect\nFAD_CRITERIA = {'prominence': 0.05, 'sharpness': 0.01, 'drop': 0.3}\n# ⚠️ Threshold ต่ำกว่า k-Shape version เพราะ KMeans centroids\n#    เป็น mean (ไม่ใช่ z-normalized) ค่าจะเล็กกว่ามาก\n\nfad_mask = (\n    (cm_df['prominence'] > FAD_CRITERIA['prominence']) &\n    (cm_df['sharpness']  > FAD_CRITERIA['sharpness']) &\n    (cm_df['drop']       > FAD_CRITERIA['drop']) &\n    (cm_df['is_early_peak'])\n)\nauto_fad_clusters = cm_df[fad_mask]['cluster'].tolist()\n\nprint(f'\\n🎯 Auto-detected FAD_CLUSTERS: {auto_fad_clusters}')\nprint(f'   Criteria: prominence > {FAD_CRITERIA[\"prominence\"]}, '\n      f'sharpness > {FAD_CRITERIA[\"sharpness\"]}, '\n      f'drop > {FAD_CRITERIA[\"drop\"]}, peak_week ≤ {FAD_MAX_PEAK_WEEK}')\n\n# Rejected clusters\nrejected = cm_df[~fad_mask]\nif len(rejected) > 0:\n    print(f'\\n   ❌ Rejected:')\n    for _, row in rejected.iterrows():\n        reasons = []\n        if row['prominence'] <= FAD_CRITERIA['prominence']: reasons.append(f'low prominence')\n        if row['sharpness'] <= FAD_CRITERIA['sharpness']: reasons.append(f'low sharpness')\n        if row['drop'] <= FAD_CRITERIA['drop']: reasons.append(f'low drop ({row[\"drop\"]})')\n        if not row['is_early_peak']: reasons.append(f'late peak (w{row[\"peak_week\"]})')\n        print(f'      Cluster {int(row[\"cluster\"])}: {\", \".join(reasons)}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 6: Visualize Centroids\n# ═══════════════════════════════════════════════════════\n\nfig, axes = plt.subplots(1, BEST_K, figsize=(max(13, 3.5 * BEST_K), 4.5))\nif BEST_K == 1: axes = [axes]\n\nfor i, ax in enumerate(axes):\n    centroid = centroids[i]\n    n = int((kshape_labels == i).sum())\n    metrics = cm_df[cm_df['cluster'] == i].iloc[0]\n    \n    is_fad = i in auto_fad_clusters\n    color = 'crimson' if is_fad else 'steelblue'\n    bg = '#ffebeb' if is_fad else 'white'\n    ax.set_facecolor(bg)\n    \n    ax.plot(centroid, color=color, lw=2.5)\n    ax.axvline(FAD_MAX_PEAK_WEEK, color='green', linestyle='--', alpha=0.3, lw=1)\n    ax.axvline(metrics['peak_week'], color='orange', linestyle=':', alpha=0.6, lw=1)\n    \n    title = f'Cluster {i} (n={n:,})\\npeak@w{metrics[\"peak_week\"]}, drop={metrics[\"drop\"]:.0%}'\n    title += '\\n🎯 FAD' if is_fad else '\\n(not Fad)'\n    ax.set_title(title, fontsize=10, fontweight='bold')\n    ax.set_xlabel('Weeks since launch'); ax.grid(alpha=0.3)\n\nplt.suptitle(f'Spark KMeans Centroids (k={BEST_K}) — Fad = peak ≤ week {FAD_MAX_PEAK_WEEK}',\n             fontweight='bold', fontsize=12)\nplt.tight_layout(); plt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 7: Save for Cell 14\n# ═══════════════════════════════════════════════════════\n\nkshape_df = pd.DataFrame({\n    'article_id': pivot.index,\n    'kshape_cluster': kshape_labels,\n})\nFAD_CLUSTERS = auto_fad_clusters\n\nn_fad = sum(cm_df[cm_df['cluster'].isin(FAD_CLUSTERS)]['n'])\nprint(f'\\n{\"=\"*60}')\nprint(f'✅ Cell 13 COMPLETE (Spark KMeans — Distributed)')\nprint(f'   Optimal k       : {BEST_K}')\nprint(f'   Silhouette Score : {silhouette:.4f}')\nprint(f'   Fad clusters     : {FAD_CLUSTERS}')\nprint(f'   Fad members      : {n_fad:,} ({n_fad/len(pivot)*100:.1f}%)')\nprint(f'   Non-Fad members  : {len(pivot)-n_fad:,}')\nprint(f'{\"=\"*60}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:08:09.726963Z","iopub.execute_input":"2026-04-29T08:08:09.727312Z","iopub.status.idle":"2026-04-29T08:10:05.143947Z","shell.execute_reply.started":"2026-04-29T08:08:09.727290Z","shell.execute_reply":"2026-04-29T08:10:05.143219Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════\n#  Cluster Validation & Profiling\n#  ตอบคำถามอาจารย์ 2 ข้อ:\n#    1. Cluster ที่เลือกเป็น Fad ถูกต้องไหม? → Statistical Tests\n#    2. ข้างใน cluster เป็นสินค้าอะไร? → Profiling\n# ══════════════════════════════════════════════════════════════════\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\nimport numpy as np\nimport pandas as pd\n\n# ─── ตรวจสอบตัวแปรที่ต้องใช้ ───\nprint('=== Pre-check ===')\nprint(f'  pivot shape:      {pivot.shape}')\nprint(f'  kshape_labels:    {len(kshape_labels)} labels, {len(set(kshape_labels))} clusters')\nprint(f'  BEST_K:           {BEST_K}')\nprint(f'  FAD_CLUSTERS:     {FAD_CLUSTERS}')\nprint(f'  centroids:        {len(centroids)} clusters')\n\n\n# ══════════════════════════════════════════════════════════════\n#  PART 1: VALIDATION — Cluster ที่เลือกถูกต้องไหม?\n# ══════════════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 1: CLUSTER VALIDATION')\nprint('  พิสูจน์ว่า Fad clusters ที่เลือกมาถูกต้อง')\nprint('=' * 70)\n\n# ─────────────────────────────────────────────────────────\n#  1.1 Visual Validation: Centroid + สมาชิกจริงทับกัน\n#      → ถ้าสมาชิกมีรูปร่างคล้าย centroid = cluster แน่น\n#      → ถ้ากระจัดกระจาย = cluster ไม่ดี\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.1 Visual: Centroid vs Actual Members ---')\n\nfig, axes = plt.subplots(2, BEST_K, figsize=(4.5 * BEST_K, 8))\n\nfor i in range(BEST_K):\n    member_idx = np.where(kshape_labels == i)[0]\n    n_members = len(member_idx)\n    \n    # สุ่มสมาชิก 50 ตัวมา plot ทับ centroid\n    sample_size = min(50, n_members)\n    sample_idx = np.random.choice(member_idx, sample_size, replace=False)\n    \n    is_fad = i in FAD_CLUSTERS\n    color = 'crimson' if is_fad else 'steelblue'\n    \n    # ── Row 1: Centroid + สมาชิกจริง overlay ──\n    ax = axes[0, i]\n    for idx in sample_idx:\n        ax.plot(pivot.iloc[idx].values, alpha=0.15, color=color, lw=0.5)\n    \n    centroid = centroids[i]\n    ax.plot(centroid, color='black', lw=3, label='Centroid')\n    ax.axvline(6, color='green', linestyle='--', alpha=0.4, lw=1, label='Fad boundary (w6)')\n    ax.set_title(f'Cluster {i} (n={n_members:,})\\n{\"🎯 FAD\" if is_fad else \"Non-Fad\"}',\n                 fontweight='bold', fontsize=11)\n    ax.set_xlabel('Weeks')\n    ax.set_ylabel('Normalized Sales')\n    ax.legend(fontsize=7, loc='upper right')\n    ax.grid(alpha=0.3)\n    \n    # ── Row 2: Peak Week Distribution ของสมาชิกจริง ──\n    # ถ้าเป็น Fad cluster → peak weeks ส่วนใหญ่ต้องอยู่ก่อนเส้น w6\n    ax2 = axes[1, i]\n    member_data = pivot.iloc[member_idx].values\n    peak_weeks = member_data.argmax(axis=1)\n    \n    ax2.hist(peak_weeks, bins=26, color=color, alpha=0.7, edgecolor='white')\n    ax2.axvline(6, color='red', linestyle='--', lw=2, label='Fad boundary (w6)')\n    \n    pct_early = (peak_weeks <= 6).mean() * 100\n    ax2.set_title(f'Peak Week Distribution\\nmedian={np.median(peak_weeks):.0f}, '\n                  f'{pct_early:.0f}% ≤ w6', fontsize=10)\n    ax2.set_xlabel('Peak Week')\n    ax2.set_ylabel('Count')\n    ax2.legend(fontsize=8)\n    ax2.grid(alpha=0.3)\n\nplt.suptitle('Cluster Validation: Centroid vs Actual Members\\n'\n             '(Row 1: เส้นบาง = สมาชิกจริง, เส้นหนา = centroid  |  '\n             'Row 2: กระจายตัวของ peak week)',\n             fontweight='bold', fontsize=12)\nplt.tight_layout()\nplt.show()\n\n# ─────────────────────────────────────────────────────────\n#  1.2 Quantitative: Cluster Quality Metrics\n#      → วัดว่าสมาชิกใน cluster คล้าย centroid แค่ไหน\n#      → Fad clusters ต้องมี early peak + high drop\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.2 Cluster Quality Metrics ---\\n')\n\n# Helper: compute drop ratios for an array of series\ndef compute_drops(data):\n    drops = []\n    for row in data:\n        pw = row.argmax()\n        peak_val = row[pw]\n        drop = (peak_val - row[-1]) / (peak_val + 1e-8) if peak_val > 0 else 0\n        drops.append(drop)\n    return np.array(drops)\n\n# Helper: compute Gini coefficients\ndef compute_gini_array(data):\n    ginis = []\n    for row in data:\n        arr = np.sort(np.abs(row))\n        n = len(arr)\n        if arr.sum() == 0:\n            ginis.append(0)\n        else:\n            idx = np.arange(1, n + 1)\n            ginis.append((2 * (idx * arr).sum()) / (n * arr.sum()) - (n + 1) / n)\n    return np.array(ginis)\n\ncluster_quality = []\nfor i in range(BEST_K):\n    member_idx = np.where(kshape_labels == i)[0]\n    member_data = pivot.iloc[member_idx].values\n    centroid = centroids[i]\n\n    # Intra-cluster distance (Euclidean — ตรงกับ KMeans)\n    distances = np.sqrt(((member_data - centroid) ** 2).sum(axis=1))\n    \n    # Peak week statistics\n    peak_weeks = member_data.argmax(axis=1)\n    pct_early_peak = (peak_weeks <= 6).mean()\n    \n    # Drop ratio statistics\n    drops = compute_drops(member_data)\n    \n    # Gini coefficient\n    ginis = compute_gini_array(member_data)\n    \n    # Early sales concentration\n    early_conc = member_data[:, :6].sum(axis=1) / (member_data.sum(axis=1) + 1e-8)\n    \n    cluster_quality.append({\n        'Cluster': i,\n        'N': len(member_idx),\n        'Is_Fad': '✅ FAD' if i in FAD_CLUSTERS else '❌ No',\n        'Avg_Dist': round(distances.mean(), 3),\n        'Median_Peak_Wk': int(np.median(peak_weeks)),\n        '%_Peak≤w6': f'{pct_early_peak:.0%}',\n        'Mean_Drop': f'{drops.mean():.0%}',\n        'Mean_Gini': round(ginis.mean(), 3),\n        'Mean_Early_Conc': f'{early_conc.mean():.0%}',\n    })\n\nquality_df = pd.DataFrame(cluster_quality)\nprint(quality_df.to_string(index=False))\n\nprint('\\nการอ่านตาราง:')\nprint('  Avg_Dist      : ระยะทางเฉลี่ยถึง centroid (ต่ำ = cluster แน่น ดี)')\nprint('  %_Peak≤w6     : สัดส่วนสมาชิกที่ peak ภายใน 6 สัปดาห์ (Fad ควรสูง)')\nprint('  Mean_Drop     : ยอดขายตกจาก peak เฉลี่ย (Fad ควรสูง)')\nprint('  Mean_Gini     : ความกระจุกตัวของยอดขาย (Fad ควรสูง)')\nprint('  Mean_Early_Conc: สัดส่วนยอดขายช่วง 6 สัปดาห์แรก (Fad ควรสูง)')\n\n# ─────────────────────────────────────────────────────────\n#  1.3 Statistical Test: Fad clusters ≠ Non-Fad clusters\n#      → Mann-Whitney U test (non-parametric)\n#      → ทดสอบว่าค่าเฉลี่ยของ 4 metrics ต่างกันอย่างมีนัยสำคัญ\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.3 Statistical Test: Fad vs Non-Fad ---\\n')\n\nfad_mask = np.isin(kshape_labels, FAD_CLUSTERS)\nfad_data = pivot.values[fad_mask]\nnonfad_data = pivot.values[~fad_mask]\n\nprint(f'Fad articles:     {fad_data.shape[0]:,}')\nprint(f'Non-Fad articles: {nonfad_data.shape[0]:,}\\n')\n\n# คำนวณ metrics สำหรับทั้ง 2 กลุ่ม\nfad_peaks = fad_data.argmax(axis=1)\nnonfad_peaks = nonfad_data.argmax(axis=1)\n\nfad_drops = compute_drops(fad_data)\nnonfad_drops = compute_drops(nonfad_data)\n\nfad_gini = compute_gini_array(fad_data)\nnonfad_gini = compute_gini_array(nonfad_data)\n\nfad_conc = fad_data[:, :6].sum(axis=1) / (fad_data.sum(axis=1) + 1e-8)\nnonfad_conc = nonfad_data[:, :6].sum(axis=1) / (nonfad_data.sum(axis=1) + 1e-8)\n\n# ทดสอบ 4 hypotheses\ntests = [\n    {\n        'name': 'Peak Week',\n        'hypothesis': 'Fad peak เร็วกว่า Non-Fad (Fad < Non-Fad)',\n        'fad_vals': fad_peaks, 'nonfad_vals': nonfad_peaks,\n        'alternative': 'less',\n        'fad_fmt': f'{fad_peaks.mean():.1f}', 'nonfad_fmt': f'{nonfad_peaks.mean():.1f}',\n    },\n    {\n        'name': 'Drop Ratio',\n        'hypothesis': 'Fad ดรอปมากกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_drops, 'nonfad_vals': nonfad_drops,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_drops.mean():.2%}', 'nonfad_fmt': f'{nonfad_drops.mean():.2%}',\n    },\n    {\n        'name': 'Gini Coefficient',\n        'hypothesis': 'Fad กระจุกตัวกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_gini, 'nonfad_vals': nonfad_gini,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_gini.mean():.3f}', 'nonfad_fmt': f'{nonfad_gini.mean():.3f}',\n    },\n    {\n        'name': 'Early Concentration',\n        'hypothesis': 'Fad ขายช่วงแรกมากกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_conc, 'nonfad_vals': nonfad_conc,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_conc.mean():.2%}', 'nonfad_fmt': f'{nonfad_conc.mean():.2%}',\n    },\n]\n\ntest_results = []\nfor t in tests:\n    stat, pval = stats.mannwhitneyu(\n        t['fad_vals'], t['nonfad_vals'], alternative=t['alternative']\n    )\n    test_results.append({\n        'Test': t['name'],\n        'Hypothesis': t['hypothesis'],\n        'Fad Mean': t['fad_fmt'],\n        'Non-Fad Mean': t['nonfad_fmt'],\n        'p-value': f'{pval:.2e}',\n        'Significant?': '✅ Yes' if pval < 0.05 else '❌ No',\n    })\n\ntest_df = pd.DataFrame(test_results)\nprint(test_df.to_string(index=False))\nprint('\\n(Mann-Whitney U test, one-sided, α = 0.05)')\n\nn_pass = sum(1 for r in test_results if '✅' in r['Significant?'])\nprint(f'\\nผ่าน {n_pass}/4 tests', end='')\nif n_pass == 4:\n    print(' → ✅ ยืนยันว่า Fad clusters แตกต่างจาก Non-Fad อย่างมีนัยสำคัญทุกมิติ')\nelif n_pass >= 3:\n    print(' → 🟡 ส่วนใหญ่ผ่าน ยอมรับได้')\nelse:\n    print(' → 🔴 ผ่านน้อย ควรทบทวนการเลือก cluster')\n\n# ─────────────────────────────────────────────────────────\n#  1.4 Visual: Distribution Comparison (Fad vs Non-Fad)\n# ─────────────────────────────────────────────────────────\n\nfig, axes = plt.subplots(1, 4, figsize=(18, 4))\n\nviz_data = [\n    (fad_peaks, nonfad_peaks, 'Peak Week', 'Fad peak เร็วกว่า'),\n    (fad_drops, nonfad_drops, 'Drop Ratio', 'Fad ดรอปมากกว่า'),\n    (fad_gini, nonfad_gini, 'Gini Coefficient', 'Fad กระจุกตัวกว่า'),\n    (fad_conc, nonfad_conc, 'Early Concentration', 'Fad ขายช่วงแรกมากกว่า'),\n]\n\nfor ax, (fad_v, nonfad_v, title, interp) in zip(axes, viz_data):\n    ax.hist(nonfad_v, bins=30, alpha=0.5, color='steelblue', label='Non-Fad', density=True)\n    ax.hist(fad_v, bins=30, alpha=0.5, color='crimson', label='Fad', density=True)\n    ax.set_title(f'{title}\\n({interp})', fontsize=10, fontweight='bold')\n    ax.legend(fontsize=8)\n    ax.grid(alpha=0.3)\n\nplt.suptitle('Statistical Validation: Fad vs Non-Fad Distribution Comparison\\n'\n             '(2 กลุ่มต้องแยกออกจากกันชัดเจน)', fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n\n# ══════════════════════════════════════════════════════════════\n#  PART 2: PROFILING — ข้างใน cluster เป็นสินค้าอะไร?\n# ══════════════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 2: CLUSTER PROFILING')\nprint('  ข้างในแต่ละ cluster เป็นสินค้าอะไรบ้าง')\nprint('=' * 70)\n\n# Merge cluster labels กับข้อมูลสินค้าจาก articles.csv\nart_pd = art.toPandas()\nart_pd['article_id'] = art_pd['article_id'].astype(str)  # ← เพิ่มบรรทัดนี้\nprofile = kshape_df.merge(art_pd, on='article_id', how='left')\nprofile['is_fad_cluster'] = profile['kshape_cluster'].isin(FAD_CLUSTERS).map(\n    {True: 'Fad', False: 'Non-Fad'}\n)\n\n# ─────────────────────────────────────────────────────────\n#  2.1 Product Group Distribution per Cluster\n#      → ดูว่าแต่ละ cluster มีสินค้ากลุ่มไหนบ้าง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.1 Product Group Distribution (% ภายในแต่ละ cluster) ---\\n')\n\ngroup_dist = pd.crosstab(\n    profile['kshape_cluster'].map(\n        lambda x: f'C{x} {\"(FAD)\" if x in FAD_CLUSTERS else \"(Non-Fad)\"}'\n    ),\n    profile['product_group_name'],\n    normalize='index'\n).round(3) * 100\n\nprint(group_dist.to_string())\n\n# ─────────────────────────────────────────────────────────\n#  2.2 Fad Rate by Product Group\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.2 Fad Rate by Product Group ---\\n')\n\nfad_by_group = (\n    profile.groupby('product_group_name')\n    .agg(\n        total=('article_id', 'count'),\n        n_fad=('is_fad_cluster', lambda x: (x == 'Fad').sum())\n    )\n)\nfad_by_group['fad_rate'] = (fad_by_group['n_fad'] / fad_by_group['total'] * 100).round(1)\nfad_by_group = fad_by_group.sort_values('fad_rate', ascending=False)\n\nprint(fad_by_group.to_string())\n\n# Plot\nfig, ax = plt.subplots(figsize=(10, 5))\nfad_by_group['fad_rate'].plot(kind='barh', ax=ax, color='crimson', alpha=0.7)\nax.set_xlabel('Fad Rate (%)')\nax.set_title('Fad Rate by Product Group\\n(กลุ่มสินค้าไหนมีสัดส่วน Fad สูงที่สุด)', fontweight='bold')\noverall_rate = profile['is_fad_cluster'].eq('Fad').mean() * 100\nax.axvline(overall_rate, color='gray', linestyle='--',\n           label=f'Overall rate ({overall_rate:.1f}%)')\nax.legend()\nax.grid(alpha=0.3)\nplt.tight_layout()\nplt.show()\n\n# ─────────────────────────────────────────────────────────\n#  2.3 Top Product Types ใน Fad vs Non-Fad\n#      → ประเภทสินค้าเฉพาะเจาะจงที่เป็น Fad มากที่สุด\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.3 Top 15 Product Types in Fad Clusters ---\\n')\n\nfad_articles = profile[profile['is_fad_cluster'] == 'Fad']\nfad_type_counts = fad_articles['product_type_name'].value_counts().head(15)\n\nprint(f'{\"Product Type\":<40s} {\"Count\":>8s} {\"% of Fad\":>10s}')\nprint('-' * 60)\nfor ptype, count in fad_type_counts.items():\n    pct = count / len(fad_articles) * 100\n    print(f'{ptype:<40s} {count:>8,} {pct:>9.1f}%')\n\n# ─────────────────────────────────────────────────────────\n#  2.4 Color Distribution: Fad vs Non-Fad\n#      → สีอะไรมี Fad สูง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.4 Color Distribution: Fad vs Non-Fad ---\\n')\n\nif 'colour_group_name' in profile.columns:\n    color_comp = pd.crosstab(\n        profile['is_fad_cluster'],\n        profile['colour_group_name'],\n        normalize='index'\n    ).T\n\n    color_comp['diff'] = color_comp.get('Fad', 0) - color_comp.get('Non-Fad', 0)\n    color_comp = color_comp.sort_values('diff', ascending=False)\n\n    print('Top 10 สีที่ Fad มีสัดส่วนสูงกว่า Non-Fad:')\n    top_colors = color_comp.head(10)\n    for color_name, row in top_colors.iterrows():\n        fad_pct = row.get('Fad', 0) * 100\n        nonfad_pct = row.get('Non-Fad', 0) * 100\n        print(f'  {color_name:<25s}  Fad: {fad_pct:5.1f}%  Non-Fad: {nonfad_pct:5.1f}%  '\n              f'(diff: {(fad_pct - nonfad_pct):+.1f}%)')\nelse:\n    print('  colour_group_name column not found')\n\n# ─────────────────────────────────────────────────────────\n#  2.5 Graphical Appearance: Fad vs Non-Fad\n#      → ลายแบบไหนมี Fad สูง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.5 Graphical Appearance: Fad vs Non-Fad ---\\n')\n\nif 'graphical_appearance_name' in profile.columns:\n    appear_comp = pd.crosstab(\n        profile['is_fad_cluster'],\n        profile['graphical_appearance_name'],\n        normalize='index'\n    ).T\n\n    appear_comp['diff'] = appear_comp.get('Fad', 0) - appear_comp.get('Non-Fad', 0)\n    appear_comp = appear_comp.sort_values('diff', ascending=False)\n\n    print('Top 10 ลายที่ Fad มีสัดส่วนสูงกว่า Non-Fad:')\n    for appear_name, row in appear_comp.head(10).iterrows():\n        fad_pct = row.get('Fad', 0) * 100\n        nonfad_pct = row.get('Non-Fad', 0) * 100\n        print(f'  {appear_name:<25s}  Fad: {fad_pct:5.1f}%  Non-Fad: {nonfad_pct:5.1f}%  '\n              f'(diff: {(fad_pct - nonfad_pct):+.1f}%)')\nelse:\n    print('  graphical_appearance_name column not found')\n\n# ─────────────────────────────────────────────────────────\n#  2.6 Department Fad Rate (Top 10, min 20 articles)\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.6 Department Fad Rate (Top 10, min 20 articles) ---\\n')\n\nif 'department_name' in profile.columns:\n    dept_fad = (\n        profile.groupby('department_name')\n        .agg(\n            total=('article_id', 'count'),\n            n_fad=('is_fad_cluster', lambda x: (x == 'Fad').sum())\n        )\n    )\n    dept_fad['fad_rate'] = (dept_fad['n_fad'] / dept_fad['total'] * 100).round(1)\n    dept_fad = dept_fad.query('total >= 20').sort_values('fad_rate', ascending=False).head(10)\n    print(dept_fad.to_string())\nelse:\n    print('  department_name column not found')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:10:05.145316Z","iopub.execute_input":"2026-04-29T08:10:05.145720Z","iopub.status.idle":"2026-04-29T08:10:13.500896Z","shell.execute_reply.started":"2026-04-29T08:10:05.145695Z","shell.execute_reply":"2026-04-29T08:10:13.499969Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 14: Hybrid Label ──\n#\n#  ✅ ใช้ FAD_CLUSTERS, BEST_K, centroids จาก Cell 13 อัตโนมัติ\n#  ✅ ไม่มี ks (KShape) reference — ใช้ centroids dict แทน\n#  ✅ ไม่มี manual override\n# ═══════════════════════════════════════════════════════════════════\n\nfrom sklearn.metrics import cohen_kappa_score\n\n# ─── ใช้ค่าจาก Cell 13 โดยตรง (ไม่ manual) ───\nprint(f'FAD_CLUSTERS = {FAD_CLUSTERS}')\nprint(f'BEST_K = {BEST_K}\\n')\n\n# แสดงรายละเอียดแต่ละ Fad cluster\nfor c in sorted(FAD_CLUSTERS):\n    n = int((kshape_labels == c).sum())\n    peak_w = int(np.argmax(centroids[c]))\n    print(f'   Cluster {c} (peak@w{peak_w}): {n:,} articles')\n\nn_fad_kshape = sum((kshape_labels == c).sum() for c in FAD_CLUSTERS)\nprint(f'   Total Fad members (k-Shape): {n_fad_kshape:,} ({n_fad_kshape/len(pivot)*100:.1f}%)')\n\n# แสดง Non-Fad clusters\nnon_fad_clusters = [c for c in range(BEST_K) if c not in FAD_CLUSTERS]\nprint(f'\\n   Non-Fad clusters: {non_fad_clusters}')\nfor c in non_fad_clusters:\n    n = int((kshape_labels == c).sum())\n    peak_w = int(np.argmax(centroids[c]))\n    print(f'   Cluster {c} (peak@w{peak_w}): {n:,} articles')\n\n# ─── Merge PELT + k-Shape labels ───\nlabels = pelt_df.merge(kshape_df, on='article_id', how='inner')\n\nlabels['kshape_fad'] = labels['kshape_cluster'].isin(FAD_CLUSTERS).astype(int)\n\nlabels['fad_label'] = (\n    (labels['pelt_fad'] == 1) & (labels['kshape_fad'] == 1)\n).astype(int)\n\n# ═══════════════════════════════════════════════════════\n#  Agreement Matrix\n# ═══════════════════════════════════════════════════════\n\nprint('\\n=== Label Agreement Matrix ===')\nagreement = pd.crosstab(\n    labels['pelt_fad'].map({0: 'PELT=0', 1: 'PELT=1'}),\n    labels['kshape_fad'].map({0: 'kShape=0', 1: 'kShape=1'}),\n    margins=True, margins_name='Total'\n)\nprint(agreement)\n\n# ═══════════════════════════════════════════════════════\n#  Rates Comparison\n# ═══════════════════════════════════════════════════════\n\npelt_rate   = labels['pelt_fad'].mean()\nkshape_rate = labels['kshape_fad'].mean()\nhybrid_rate = labels['fad_label'].mean()\n\nprint(f'\\n=== Rates Comparison ===')\nprint(f'PELT alone       : {pelt_rate:.2%}  ({labels[\"pelt_fad\"].sum():>6,} articles)')\nprint(f'k-Shape alone    : {kshape_rate:.2%}  ({labels[\"kshape_fad\"].sum():>6,} articles)')\nprint(f'Hybrid (AND)     : {hybrid_rate:.2%}  ({labels[\"fad_label\"].sum():>6,} articles)')\nprint(f'Phase 1 baseline : 4.10%   (death ≤ 6 weeks estimate)')\n\n# ═══════════════════════════════════════════════════════\n#  Cohen's Kappa\n# ═══════════════════════════════════════════════════════\n\nkappa = cohen_kappa_score(labels['pelt_fad'], labels['kshape_fad'])\nprint(f'\\n=== Cohen\\'s Kappa (Inter-Method Agreement) ===')\nprint(f'  Kappa = {kappa:.3f}')\nif kappa >= 0.60:\n    print('  ✅ Substantial agreement')\nelif kappa >= 0.40:\n    print('  🟡 Moderate agreement')\nelif kappa >= 0.20:\n    print('  🟠 Fair agreement')\nelse:\n    print('  🔴 Weak agreement')\n\n# ═══════════════════════════════════════════════════════\n#  Label Composition by Cluster\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Fad Label Composition by Cluster ===')\nfad_composition = labels[labels['fad_label'] == 1].groupby('kshape_cluster').size()\nfor cluster, count in fad_composition.items():\n    pct = count / labels['fad_label'].sum() * 100\n    peak_w = int(np.argmax(centroids[cluster]))\n    print(f'  Cluster {cluster} (peak@w{peak_w}): {count:>5,} articles ({pct:.1f}% of all Fads)')\n\n# ═══════════════════════════════════════════════════════\n#  Temporal Profile\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Final Fad Temporal Profile ===')\nfinal_fads = labels[labels['fad_label'] == 1]\nif len(final_fads) > 0:\n    print(f'  Peak week   — median: {final_fads[\"peak_week\"].median():.0f}, '\n          f'mean: {final_fads[\"peak_week\"].mean():.1f}')\n    print(f'  Peak ratio  — mean: {final_fads[\"peak_ratio\"].mean():.2f}')\n    print(f'  Drop ratio  — mean: {final_fads[\"drop_ratio\"].mean():.2f}')\n    print(f'  Tail ratio  — mean: {final_fads[\"tail_ratio\"].mean():.3f}')\n\n# ═══════════════════════════════════════════════════════\n#  Final Distribution\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Final Hybrid Label Distribution ===')\nprint(labels['fad_label'].value_counts())\nprint(f'Final Fad rate: {hybrid_rate:.2%}')\n\n# ═══════════════════════════════════════════════════════\n#  Sanity Check\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Sanity Check ===')\nif hybrid_rate < 0.02:\n    print('  🔴 ต่ำเกินไป < 2%')\nelif hybrid_rate > 0.15:\n    print('  🟡 สูงกว่าเป้าเล็กน้อย > 15%')\nelif 0.03 <= hybrid_rate <= 0.12:\n    print('  ✅ Good — อยู่ในช่วงที่เหมาะสม')\nelse:\n    print(f'  🟡 Acceptable — {hybrid_rate:.2%}')\n\nprint('\\n✅ Phase 2 complete — ready for Phase 3')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:10:13.502051Z","iopub.execute_input":"2026-04-29T08:10:13.502369Z","iopub.status.idle":"2026-04-29T08:10:13.634588Z","shell.execute_reply.started":"2026-04-29T08:10:13.502345Z","shell.execute_reply":"2026-04-29T08:10:13.633870Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 3 — Feature Engineering with PySpark","metadata":{}},{"cell_type":"code","source":"# ── Cell 15: Dynamic Features  ──\n\nEARLY_WINDOW = 5\nMAX_WEEKS = 26\n\n# ═══════════════════════════════════════════════════════\n#  Part A: Weekly sales in early window (Spark parallel)\n# ═══════════════════════════════════════════════════════\n\nearly_spark = (\n    weekly\n    .filter(\n        (F.col('weeks_since_launch') >= 0) &\n        (F.col('weeks_since_launch') <= EARLY_WINDOW)\n    )\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(EARLY_WINDOW + 1)))\n    .agg(F.first('weekly_sales_clean'))\n    .fillna(0)\n)\n\n# Rename: '0' → 'sales_w0', etc.\nearly_cols = early_spark.columns\nearly_spark = early_spark.toDF(\n    *(['article_id'] + [f'sales_w{c}' for c in early_cols[1:]])\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part B: Aggregate features from early window\n# ═══════════════════════════════════════════════════════\n\nearly_agg = (\n    weekly\n    .filter(F.col('weeks_since_launch') <= EARLY_WINDOW)\n    .groupBy('article_id')\n    .agg(\n        F.max('weekly_sales_clean').alias('early_peak_sales'),\n        F.sum('weekly_sales_clean').alias('early_total_sales'),\n        F.sum('unique_buyers').alias('early_unique_buyers'),\n        F.max('n_channels').alias('early_max_channels'),\n        F.avg('weekly_sales_clean').alias('early_avg_sales'),\n    )\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part C: Full lifecycle aggregates (within 26 weeks)\n# ═══════════════════════════════════════════════════════\n\nlifecycle_agg = (\n    weekly\n    .filter(F.col('weeks_since_launch') < MAX_WEEKS)\n    .groupBy('article_id')\n    .agg(\n        F.max('weekly_sales_clean').alias('peak_sales'),\n        F.sum('weekly_sales_clean').alias('total_sales_clean'),\n        F.avg('weekly_sales_clean').alias('avg_sales'),\n    )\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part D: time_to_peak (within 26 weeks ONLY)\n# ═══════════════════════════════════════════════════════\n\nw_peak = Window.partitionBy('article_id').orderBy(F.desc('weekly_sales_clean'))\n\ntime_to_peak = (\n    weekly\n    .filter(F.col('weeks_since_launch') < MAX_WEEKS)  # ✅ จำกัด 26 weeks\n    .withColumn('rn', F.row_number().over(w_peak))\n    .filter(F.col('rn') == 1)\n    .select('article_id', F.col('weeks_since_launch').alias('time_to_peak'))\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part E: post_peak_decay (จาก pivot — accurate)\n# ═══════════════════════════════════════════════════════\n\n# คำนวณบน Pandas (pivot มีแค่ ~80K rows × 26 cols = เล็ก)\npeak_weeks_pd = time_to_peak.toPandas().set_index('article_id')\n\ndecay_results = []\nfor art_id, row in pivot.iterrows():\n    series = row.values.astype(float)\n    \n    # Get peak week from pre-computed time_to_peak\n    if art_id in peak_weeks_pd.index:\n        pw = int(peak_weeks_pd.loc[art_id, 'time_to_peak'])\n    else:\n        pw = int(np.argmax(series))\n    \n    peak_val = series[pw]\n    \n    # Compute decay slope: peak → 3 weeks after peak\n    if pw < MAX_WEEKS - 3:\n        end_w = pw + 3\n        end_val = series[end_w]\n        decay = (end_val - peak_val) / 3.0\n    elif pw < MAX_WEEKS - 1:\n        end_val = series[-1]\n        decay = (end_val - peak_val) / max(1, MAX_WEEKS - 1 - pw)\n    else:\n        decay = 0.0\n    \n    # Is peak in early window? (strong Fad indicator)\n    is_early_peak = int(pw <= EARLY_WINDOW)\n    \n    decay_results.append({\n        'article_id': art_id,\n        'post_peak_decay': round(decay, 4),\n        'is_early_peak_flag': is_early_peak,\n    })\n\ndecay_spark = spark.createDataFrame(pd.DataFrame(decay_results))\n\n# ═══════════════════════════════════════════════════════\n#  Part F: Combine all + compute derived features\n# ═══════════════════════════════════════════════════════\n\ndyn = (\n    early_spark\n    .join(early_agg,      on='article_id', how='left')\n    .join(lifecycle_agg,  on='article_id', how='left')\n    .join(time_to_peak,   on='article_id', how='left')\n    .join(decay_spark,    on='article_id', how='left')\n    \n    # ── Velocity & Acceleration ──\n    .withColumn('velocity_w1_w0',\n        F.col('sales_w1') / (F.col('sales_w0') + 1.0))\n    .withColumn('velocity_w2_w1',\n        F.col('sales_w2') / (F.col('sales_w1') + 1.0))\n    .withColumn('acceleration',\n        (F.col('sales_w2') - F.col('sales_w1')) - \n        (F.col('sales_w1') - F.col('sales_w0')))\n    \n    # ── Peak ratios ──\n    .withColumn('peak_to_w0_ratio',\n        F.col('peak_sales') / (F.col('sales_w0') + 1.0))\n    .withColumn('early_peak_to_avg',\n        F.col('early_peak_sales') / (F.col('early_avg_sales') + 1.0))\n    \n    # ── Pre-peak velocity ──\n    .withColumn('pre_peak_velocity',\n        F.when(F.col('time_to_peak') > 0,\n            F.col('early_peak_sales') / (F.col('time_to_peak') + 1.0)\n        ).otherwise(F.col('sales_w0')))\n    \n    # ── Early concentration (% of total sales in early window) ──\n    .withColumn('early_sales_concentration',\n        F.col('early_total_sales') / (F.col('total_sales_clean') + 1.0))\n    \n    # ── Buyer intensity ──\n    .withColumn('early_buyer_intensity',\n        F.col('early_unique_buyers') / (F.col('early_total_sales') + 1.0))\n)\n\ndyn.cache()\nn_features = len(dyn.columns) - 1  # minus article_id\n\n# ═══════════════════════════════════════════════════════\n#  Diagnostics\n# ═══════════════════════════════════════════════════════\n\nprint(f'✅ Dynamic features computed.')\nprint(f'   Early window: weeks 0-{EARLY_WINDOW}')\nprint(f'   Features: {n_features}')\n\n# Verify time_to_peak is within 0-25\nprint('\\n=== time_to_peak Distribution (should be 0-25) ===')\ndyn.select(\n    F.min('time_to_peak').alias('min'),\n    F.max('time_to_peak').alias('max'),\n    F.mean('time_to_peak').alias('mean'),\n    F.expr('percentile_approx(time_to_peak, 0.5)').alias('median'),\n).show()\n\n# Verify post_peak_decay is not all zeros\nprint('=== post_peak_decay Distribution ===')\ndyn.select(\n    F.mean('post_peak_decay').alias('mean'),\n    F.min('post_peak_decay').alias('min'),\n    F.max('post_peak_decay').alias('max'),\n    F.sum(F.when(F.col('post_peak_decay') != 0, 1).otherwise(0)).alias('non_zero_count'),\n    F.count('*').alias('total'),\n).show()\n\n# Preview\nprint('--- Sample output ---')\ndyn.select('article_id', 'velocity_w1_w0', 'time_to_peak', \n           'post_peak_decay', 'is_early_peak_flag',\n           'early_sales_concentration').show(10)\n\n# Fad-likely articles preview\nn_early_peak = dyn.filter(F.col('is_early_peak_flag') == 1).count()\nprint(f'\\nArticles with early peak (≤ week {EARLY_WINDOW}): {n_early_peak:,} '\n      f'({n_early_peak/dyn.count()*100:.1f}%)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:10:13.635642Z","iopub.execute_input":"2026-04-29T08:10:13.635960Z","iopub.status.idle":"2026-04-29T08:10:29.855124Z","shell.execute_reply.started":"2026-04-29T08:10:13.635937Z","shell.execute_reply":"2026-04-29T08:10:29.854293Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 16: Gini Coefficient  ──\n\nfrom pyspark.sql.types import DoubleType\n\n@F.udf(DoubleType())\ndef gini_udf(sales_list):\n    \"\"\"Gini coefficient: 0 = equal distribution, 1 = all sales in one week\"\"\"\n    if not sales_list:\n        return 0.0\n    arr = np.sort(np.abs(np.array(sales_list, dtype=float)))\n    n = len(arr)\n    if n == 0 or arr.sum() == 0:\n        return 0.0\n    idx = np.arange(1, n + 1)\n    return float((2 * (idx * arr).sum()) / (n * arr.sum()) - (n + 1) / n)\n\ngini_spark = (\n    weekly.groupBy('article_id')\n    .agg(F.collect_list('weekly_sales_clean').alias('sales_list'))\n    .withColumn('gini_coeff', gini_udf('sales_list'))\n    .select('article_id', 'gini_coeff')\n)\n\ngini_spark.cache()\nprint('✅ Gini computed.')\n\n# Diagnostic: Fad ควรมี Gini สูง (sales กระจุกตัว)\ngini_pd = gini_spark.toPandas()\nprint(f'   Mean Gini: {gini_pd[\"gini_coeff\"].mean():.3f}')\nprint(f'   Std:  {gini_pd[\"gini_coeff\"].std():.3f}')\ngini_spark.show(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:10:29.856047Z","iopub.execute_input":"2026-04-29T08:10:29.856369Z","iopub.status.idle":"2026-04-29T08:10:35.066060Z","shell.execute_reply.started":"2026-04-29T08:10:29.856347Z","shell.execute_reply":"2026-04-29T08:10:35.065190Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 17: Repeat Purchase + Interpurchase CV ──\n\n# ─── Part A: Repeat Purchase Rate (per article) ───\nmulti_buy = (\n    txn.groupBy('article_id', 'customer_id')\n    .agg(F.count('*').alias('n_purchases'))\n)\n\nbuyer_stats = (\n    multi_buy.groupBy('article_id')\n    .agg(\n        F.count('customer_id').alias('total_buyers'),\n        F.sum((F.col('n_purchases') >= 2).cast('int')).alias('repeat_buyers'),\n        F.avg('n_purchases').alias('avg_purchases_per_buyer'),\n    )\n    .withColumn('repeat_purchase_rate',\n        F.col('repeat_buyers') / (F.col('total_buyers') + 1.0))\n)\n\n# ─── Part B: Inter-purchase CV (per article — FIXED) ───\n# นี่วัดว่า \"transaction ของ article เกิดถี่แค่ไหนและสม่ำเสมอแค่ไหน\"\n\nw_art_txn = Window.partitionBy('article_id').orderBy('t_dat')\n\ncv_spark = (\n    txn\n    .withColumn('row_num', F.row_number().over(w_art_txn))\n    .withColumn('lag_date', F.lag('t_dat', 1).over(w_art_txn))\n    .withColumn('gap_days', F.datediff('t_dat', 'lag_date').cast('double'))\n    .filter(F.col('gap_days').isNotNull())\n    .groupBy('article_id')\n    .agg(\n        F.mean('gap_days').alias('mean_gap'),\n        F.stddev('gap_days').alias('std_gap'),\n        F.count('gap_days').alias('n_gaps'),\n    )\n    # CV = std / mean (high CV = burst pattern = Fad signal)\n    .withColumn('interpurchase_cv',\n        F.when(F.col('n_gaps') >= 3,  # ต้องมี gaps พอ ≥ 3 ถึงจะ compute CV ได้\n            F.col('std_gap') / (F.col('mean_gap') + 1.0)\n        ).otherwise(0.0))\n    .select('article_id', 'mean_gap', 'interpurchase_cv')\n)\n\nbuyer_stats.cache()\ncv_spark.cache()\nprint('✅ Buyer behavior features done.')\nprint(f'   Buyer stats rows: {buyer_stats.count():,}')\nprint(f'   CV stats rows: {cv_spark.count():,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:10:35.067194Z","iopub.execute_input":"2026-04-29T08:10:35.067539Z","iopub.status.idle":"2026-04-29T08:12:17.159901Z","shell.execute_reply.started":"2026-04-29T08:10:35.067503Z","shell.execute_reply":"2026-04-29T08:12:17.158962Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 18: Static Features + Target Encoding ──\n\n# ─── Part A: Select static features ───\nart_feat_spark = art.select(\n    'article_id',\n    'product_type_no',           # 131 unique → Target Encoding\n    'department_no',             # 250 unique → Target Encoding\n    'graphical_appearance_no',   # 30 unique → keep as-is (OHE later if needed)\n    'colour_group_code',         # 50 unique → keep as-is\n    'perceived_colour_value_id',\n    'perceived_colour_master_id',\n    'index_code',                # 10 unique → keep as-is\n    'index_group_no',\n    'section_no',\n    'garment_group_no',          # 21 unique → keep as-is\n)\n\n# ─── Part B: Target Encoding for high-cardinality features ───\n# ✅ Phase 0 กำหนดไว้: department_no (250) + product_type_no (131) → Target Encoding\n# Logic: replace category code with mean(fad_label) of that category\n# ⚠️ ต้องทำบน TRAIN set only ใน Phase 5 เพื่อป้องกัน data leakage\n# ตอนนี้คำนวณ global mean ไว้ก่อน → Phase 5 จะ re-compute per fold\n\nlabels_spark = spark.createDataFrame(\n    labels[['article_id', 'fad_label']].astype({'fad_label': 'int32'})\n)\n\n# Target encoding: department_no\ndept_target = (\n    art_feat_spark\n    .join(labels_spark, 'article_id', 'inner')\n    .groupBy('department_no')\n    .agg(\n        F.mean('fad_label').alias('dept_fad_rate'),\n        F.count('*').alias('dept_count'),\n    )\n)\n\n# Target encoding: product_type_no\nptype_target = (\n    art_feat_spark\n    .join(labels_spark, 'article_id', 'inner')\n    .groupBy('product_type_no')\n    .agg(\n        F.mean('fad_label').alias('ptype_fad_rate'),\n        F.count('*').alias('ptype_count'),\n    )\n)\n\n# Smoothed target encoding: smooth toward global mean for rare categories\nglobal_fad_rate = labels['fad_label'].mean()\nSMOOTH_FACTOR = 20  # categories with < 20 items → lean toward global mean\n\ndept_target = dept_target.withColumn('dept_target_enc',\n    (F.col('dept_fad_rate') * F.col('dept_count') + global_fad_rate * SMOOTH_FACTOR) /\n    (F.col('dept_count') + SMOOTH_FACTOR)\n)\n\nptype_target = ptype_target.withColumn('ptype_target_enc',\n    (F.col('ptype_fad_rate') * F.col('ptype_count') + global_fad_rate * SMOOTH_FACTOR) /\n    (F.col('ptype_count') + SMOOTH_FACTOR)\n)\n\n# Join target encodings back\nart_feat_spark = (\n    art_feat_spark\n    .join(dept_target.select('department_no', 'dept_target_enc'), 'department_no', 'left')\n    .join(ptype_target.select('product_type_no', 'ptype_target_enc'), 'product_type_no', 'left')\n)\n\n# ─── Part C: Interaction features ───\n# All-over pattern × bright colour (from methodology)\nart_feat_spark = (\n    art_feat_spark\n    .withColumn('is_allover_pattern',\n        (F.col('graphical_appearance_no').isin([74])).cast('int'))\n    .withColumn('is_bright_colour',\n        (F.col('perceived_colour_value_id').isin([3, 4])).cast('int'))\n    .withColumn('allover_x_bright',\n        F.col('is_allover_pattern') * F.col('is_bright_colour'))\n)\n\nart_feat_spark.cache()\nprint('✅ Static features done.')\nprint(f'   Columns: {len(art_feat_spark.columns)}')\nprint(f'   Target encoding: department_no, product_type_no (smoothed, factor={SMOOTH_FACTOR})')\nprint(f'   Global fad rate for smoothing: {global_fad_rate:.4f}')\n\n# Preview target encoding\nprint('\\n--- Department Target Encoding (top 5 highest Fad rate) ---')\ndept_target.orderBy(F.desc('dept_fad_rate')).show(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:12:17.161118Z","iopub.execute_input":"2026-04-29T08:12:17.161812Z","iopub.status.idle":"2026-04-29T08:12:19.010269Z","shell.execute_reply.started":"2026-04-29T08:12:17.161783Z","shell.execute_reply":"2026-04-29T08:12:19.009553Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 19 REVISED: Merge All Features ──\n\n# ─── Part A: Stockout features per article (from Phase 1) ───\nstockout_features = (\n    weekly\n    .groupBy('article_id')\n    .agg(\n        F.sum('is_stockout').alias('n_stockout_weeks'),\n        F.max('stockout_tier').alias('max_stockout_tier'),\n        # Stockout rate in early window\n        F.sum(\n            F.when(\n                (F.col('weeks_since_launch') <= EARLY_WINDOW) & (F.col('is_stockout') == 1), 1\n            ).otherwise(0)\n        ).alias('early_stockout_count'),\n    )\n)\n\n# ─── Part B: Labels (from Phase 2) ───\nlabels_spark = spark.createDataFrame(\n    labels[['article_id', 'fad_label', 'peak_week', 'kshape_cluster',\n            'peak_ratio', 'drop_ratio', 'tail_ratio']].copy()\n)\n\n# ─── Part C: Regime info (from Phase 1) ───\nregime_spark = total_sales_df.select('article_id', 'regime', 'strat_key', 'total_sales')\n\n# ─── Part D: Merge everything ───\nfeatures_spark = (\n    labels_spark.select('article_id', 'fad_label')  # target + ID\n    .join(dyn,              on='article_id', how='left')   # dynamic features\n    .join(gini_spark,       on='article_id', how='left')   # gini\n    .join(buyer_stats.select('article_id', 'total_buyers', \n                             'repeat_purchase_rate', 'avg_purchases_per_buyer'),\n                            on='article_id', how='left')   # buyer behavior\n    .join(cv_spark,         on='article_id', how='left')   # interpurchase CV\n    .join(art_feat_spark,   on='article_id', how='left')   # static features\n    .join(stockout_features,on='article_id', how='left')   # stockout features\n    .join(regime_spark,     on='article_id', how='left')   # regime\n    .fillna(0)\n)\n\nfeatures_spark.cache()\nn_rows = features_spark.count()\nn_cols = len(features_spark.columns)\n\nprint(f'✅ Feature matrix assembled.')\nprint(f'   Rows: {n_rows:,}')\nprint(f'   Columns: {n_cols}')\n\n# ─── Diagnostic: Feature categories ───\nfeature_cols = [c for c in features_spark.columns if c not in ['article_id', 'fad_label']]\nprint(f'   Feature columns: {len(feature_cols)}')\n\n# Check for unexpected nulls\nnull_counts = features_spark.select([\n    F.sum(F.col(c).isNull().cast('int')).alias(c) for c in feature_cols\n]).collect()[0].asDict()\nnulls = {k: v for k, v in null_counts.items() if v > 0}\nif nulls:\n    print(f'\\n   ⚠️  Columns with nulls after fillna: {nulls}')\nelse:\n    print(f'   ✅ No nulls remaining')\n\n# ─── Collect to Pandas ───\nfeatures = features_spark.toPandas()\nprint(f'\\n   Pandas shape: {features.shape}')\n\n# ─── Label distribution in final feature matrix ───\nprint(f'\\n=== Label Distribution ===')\nprint(features['fad_label'].value_counts())\nprint(f'Fad rate: {features[\"fad_label\"].mean():.2%}')\n\n# ─── Feature summary by type ───\nprint(f'\\n=== Feature Summary ===')\ndynamic_cols = ['sales_w0','sales_w1','sales_w2','sales_w3','sales_w4','sales_w5',\n                'velocity_w1_w0','velocity_w2_w1','acceleration',\n                'peak_to_w0_ratio','early_peak_to_avg','post_peak_decay',\n                'pre_peak_velocity','early_sales_concentration','early_buyer_intensity',\n                'time_to_peak','early_peak_sales','early_total_sales',\n                'early_unique_buyers','early_max_channels','early_avg_sales',\n                'peak_sales','total_sales_clean','avg_sales']\nbehavioral_cols = ['gini_coeff','total_buyers','repeat_purchase_rate',\n                   'avg_purchases_per_buyer','interpurchase_cv','mean_gap']\nstatic_cols = ['product_type_no','department_no','graphical_appearance_no',\n               'colour_group_code','perceived_colour_value_id','perceived_colour_master_id',\n               'index_code','index_group_no','section_no','garment_group_no',\n               'dept_target_enc','ptype_target_enc',\n               'is_allover_pattern','is_bright_colour','allover_x_bright']\noperational_cols = ['n_stockout_weeks','max_stockout_tier','early_stockout_count']\n\nfor name, cols in [('Dynamic (early signals)', dynamic_cols),\n                    ('Behavioral (buyer)', behavioral_cols),\n                    ('Static (article metadata)', static_cols),\n                    ('Operational (stockout)', operational_cols)]:\n    present = [c for c in cols if c in features.columns]\n    print(f'  {name}: {len(present)} features')\n\nprint(f'\\n✅ Phase 3 complete — ready for Phase 4 (Feature Selection)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:12:19.013874Z","iopub.execute_input":"2026-04-29T08:12:19.014195Z","iopub.status.idle":"2026-04-29T08:12:35.963294Z","shell.execute_reply.started":"2026-04-29T08:12:19.014172Z","shell.execute_reply":"2026-04-29T08:12:35.962510Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 4 — Feature Selection","metadata":{}},{"cell_type":"code","source":"# ── Cell 20: Feature Selection (MI + Correlation + RFE) ──\n\nfrom sklearn.feature_selection import RFE, mutual_info_classif\nimport lightgbm as lgb\n\n# ─── Step 0: Prepare X and y ───\n# Columns ที่ไม่ใช่ feature (IDs, labels, metadata)\nDROP_COLS = [\n    'article_id', 'fad_label',              # ID + target\n    'regime', 'strat_key', 'total_sales',   # Phase 1 metadata (ใช้สำหรับ stratification ไม่ใช่ prediction)\n]\n\n# ตรวจสอบว่ามี columns ไหนอยู่จริง\nexisting_drops = [c for c in DROP_COLS if c in features.columns]\nprint(f'Dropping {len(existing_drops)} non-feature columns: {existing_drops}')\n\nX_df = features.drop(columns=existing_drops, errors='ignore').copy()\ny_all = features['fad_label'].values\n\n# Encode categorical columns (object/category → numeric)\ncat_cols = X_df.select_dtypes(include=['object', 'category']).columns.tolist()\nif cat_cols:\n    print(f'Encoding {len(cat_cols)} categorical columns: {cat_cols}')\n    for col in cat_cols:\n        X_df[col] = X_df[col].astype('category').cat.codes\n\nprint(f'\\nStarting feature selection: {X_df.shape[1]} features, {len(y_all):,} samples')\nprint(f'Fad rate: {y_all.mean():.2%}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 1: Mutual Information (非linear relationship)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 1: Mutual Information')\nprint('='*50)\n\nmi = mutual_info_classif(X_df, y_all, random_state=SEED, n_neighbors=5)\nmi_df = pd.DataFrame({\n    'feature': X_df.columns, \n    'mi': mi\n}).sort_values('mi', ascending=False)\n\n# Threshold: keep features with MI > 0.001\nMI_THRESH = 0.001\nsel_mi = mi_df.loc[mi_df['mi'] >= MI_THRESH, 'feature'].tolist()\ndropped_mi = mi_df.loc[mi_df['mi'] < MI_THRESH, 'feature'].tolist()\n\nprint(f'MI threshold: {MI_THRESH}')\nprint(f'Features kept: {len(sel_mi)}')\nif dropped_mi:\n    print(f'Features dropped (no predictive power): {dropped_mi}')\n\n# Visualization: MI scores\nfig, ax = plt.subplots(figsize=(10, max(6, len(mi_df) * 0.3)))\ncolors = ['crimson' if mi >= MI_THRESH else 'lightgray' for mi in mi_df['mi']]\nax.barh(range(len(mi_df)), mi_df['mi'].values, color=colors)\nax.set_yticks(range(len(mi_df)))\nax.set_yticklabels(mi_df['feature'].values, fontsize=8)\nax.axvline(MI_THRESH, color='red', linestyle='--', alpha=0.7, label=f'Threshold={MI_THRESH}')\nax.set_xlabel('Mutual Information')\nax.set_title('Step 1: Mutual Information Scores\\n(red = kept, gray = dropped)', fontweight='bold')\nax.legend()\nax.invert_yaxis()\nplt.tight_layout()\nplt.show()\n\n# Top 10 features by MI\nprint('\\nTop 10 features by MI:')\nprint(mi_df.head(10).to_string(index=False))\n\n# ═══════════════════════════════════════════════════════\n#  Step 2: Correlation Pruning (remove redundancy)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 2: Correlation Pruning')\nprint('='*50)\n\nCORR_THRESH = 0.85\n\ncorr = X_df[sel_mi].astype('float32').corr().abs()\nupper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))\n\n# For each pair with corr > threshold, drop the one with LOWER MI\nto_drop_corr = set()\nfor col in upper.columns:\n    high_corr_cols = upper.index[upper[col] > CORR_THRESH].tolist()\n    for hc in high_corr_cols:\n        # Keep the one with higher MI\n        mi_col = mi_df.set_index('feature').loc[col, 'mi']\n        mi_hc = mi_df.set_index('feature').loc[hc, 'mi']\n        drop_this = hc if mi_col >= mi_hc else col\n        to_drop_corr.add(drop_this)\n\nsel_corr = [c for c in sel_mi if c not in to_drop_corr]\n\nprint(f'Correlation threshold: {CORR_THRESH}')\nprint(f'Features kept: {len(sel_corr)}')\nif to_drop_corr:\n    print(f'Features dropped (redundant): {sorted(to_drop_corr)}')\n\n# Visualization: Correlation heatmap of remaining features\nfig, ax = plt.subplots(figsize=(12, 10))\ncorr_final = X_df[sel_corr].astype('float32').corr()\nmask = np.triu(np.ones_like(corr_final), k=1)\nsns.heatmap(corr_final, mask=mask, cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, annot=False, square=True,\n            linewidths=0.5, ax=ax, cbar_kws={'shrink': 0.8})\nax.set_title(f'Step 2: Correlation Matrix ({len(sel_corr)} features after pruning)',\n             fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 3: RFE with LightGBM\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 3: Recursive Feature Elimination (RFE)')\nprint('='*50)\n\nX_rfe = X_df[sel_corr]\n\n# Target: select 15-20 features (lean but powerful)\nN_TARGET = min(20, len(sel_corr))\n\nmodel = lgb.LGBMClassifier(\n    n_estimators=200,\n    max_depth=6,\n    learning_rate=0.05,\n    random_state=SEED,\n    class_weight='balanced',\n    verbose=-1,\n    n_jobs=-1,\n)\n\nrfe = RFE(\n    estimator=model,\n    n_features_to_select=N_TARGET,\n    step=2,\n    verbose=0,\n)\n\nrfe.fit(X_rfe, y_all)\n\nFINAL_FEATURES = [f for f, s in zip(sel_corr, rfe.support_) if s]\nrfe_rankings = pd.DataFrame({\n    'feature': sel_corr,\n    'rfe_rank': rfe.ranking_,\n    'selected': rfe.support_,\n}).sort_values('rfe_rank')\n\nprint(f'Target features: {N_TARGET}')\nprint(f'Final features: {len(FINAL_FEATURES)}')\n\n# Show RFE rankings\nprint('\\nRFE Rankings (1 = selected):')\nprint(rfe_rankings.to_string(index=False))\n\n# ═══════════════════════════════════════════════════════\n#  Step 4: Summary & Comparison\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('FINAL FEATURE SET')\nprint('='*50)\n\n# Combine MI and RFE info\nfinal_info = (\n    mi_df[mi_df['feature'].isin(FINAL_FEATURES)]\n    .merge(rfe_rankings[rfe_rankings['selected']], on='feature')\n    .sort_values('mi', ascending=False)\n)\n\nprint(f'\\n{len(FINAL_FEATURES)} features selected:')\nprint(final_info[['feature', 'mi', 'rfe_rank']].to_string(index=False))\n\n# Categorize final features\ndynamic_final = [f for f in FINAL_FEATURES if any(f.startswith(p) for p in \n    ['sales_w', 'velocity', 'acceleration', 'peak_', 'early_', 'time_to', 'post_peak', 'pre_peak', 'is_early'])]\nbehavioral_final = [f for f in FINAL_FEATURES if f in \n    ['gini_coeff', 'total_buyers', 'repeat_purchase_rate', 'avg_purchases_per_buyer', 'interpurchase_cv', 'mean_gap']]\nstatic_final = [f for f in FINAL_FEATURES if f in \n    ['product_type_no', 'department_no', 'graphical_appearance_no', 'colour_group_code',\n     'perceived_colour_value_id', 'perceived_colour_master_id', 'index_code', 'index_group_no',\n     'section_no', 'garment_group_no', 'dept_target_enc', 'ptype_target_enc',\n     'is_allover_pattern', 'is_bright_colour', 'allover_x_bright']]\noperational_final = [f for f in FINAL_FEATURES if f in \n    ['n_stockout_weeks', 'max_stockout_tier', 'early_stockout_count']]\n\nprint(f'\\nFeature breakdown:')\nprint(f'  Dynamic (early signals)     : {len(dynamic_final)}  {dynamic_final}')\nprint(f'  Behavioral (buyer)          : {len(behavioral_final)}  {behavioral_final}')\nprint(f'  Static (article metadata)   : {len(static_final)}  {static_final}')\nprint(f'  Operational (stockout)       : {len(operational_final)}  {operational_final}')\n\n# Visualization: Final feature importance (MI)\nfig, ax = plt.subplots(figsize=(10, max(5, len(FINAL_FEATURES) * 0.35)))\nfinal_sorted = final_info.sort_values('mi', ascending=True)\n\n# Color by category\ndef get_color(f):\n    if f in dynamic_final: return 'steelblue'\n    if f in behavioral_final: return 'coral'\n    if f in static_final: return 'mediumseagreen'\n    if f in operational_final: return 'gold'\n    return 'gray'\n\ncolors = [get_color(f) for f in final_sorted['feature']]\nax.barh(range(len(final_sorted)), final_sorted['mi'].values, color=colors)\nax.set_yticks(range(len(final_sorted)))\nax.set_yticklabels(final_sorted['feature'].values, fontsize=9)\nax.set_xlabel('Mutual Information')\nax.set_title(f'Final {len(FINAL_FEATURES)} Features by MI Score', fontweight='bold')\n\n# Legend\nfrom matplotlib.patches import Patch\nlegend_items = [\n    Patch(color='steelblue', label='Dynamic'),\n    Patch(color='coral', label='Behavioral'),\n    Patch(color='mediumseagreen', label='Static'),\n    Patch(color='gold', label='Operational'),\n]\nax.legend(handles=legend_items, loc='lower right')\nplt.tight_layout()\nplt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Prepare final X, y\n# ═══════════════════════════════════════════════════════\n\nX = X_df[FINAL_FEATURES].copy()\ny = y_all\n\nprint(f'\\n✅ Phase 4 complete')\nprint(f'   X shape: {X.shape}')\nprint(f'   y shape: {y.shape}')\nprint(f'   Fad rate: {y.mean():.2%}')\nprint(f'   Ready for Phase 5 (Model Development)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:12:35.964323Z","iopub.execute_input":"2026-04-29T08:12:35.964652Z","iopub.status.idle":"2026-04-29T08:13:10.266754Z","shell.execute_reply.started":"2026-04-29T08:12:35.964623Z","shell.execute_reply":"2026-04-29T08:13:10.265961Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 20b: Export (FIXED) ──\nimport os, json\n\nOUTPUT_DIR = '/kaggle/working/processed_data'\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\n# ✅ สร้าง export_df จาก X_df (ที่มี article_id) + y\n# X_df สร้างจาก features ตอน Phase 4 แต่ drop article_id ไปแล้ว\n# ต้องดึง article_id กลับมาจาก features_spark\n\n# Reconstruct: ดึง article_id + fad_label + final features\nexport_df = features_spark.select(\n    'article_id', 'fad_label', *FINAL_FEATURES\n).toPandas()\n\n# 1. Save features\nexport_df.to_parquet(f'{OUTPUT_DIR}/features_final.parquet', index=False)\nexport_df.to_csv(f'{OUTPUT_DIR}/features_final.csv', index=False)\nprint(f'✅ Feature matrix: {export_df.shape}')\n\n# 2. Save labels\nlabels.to_parquet(f'{OUTPUT_DIR}/labels_hybrid.parquet', index=False)\nlabels.to_csv(f'{OUTPUT_DIR}/labels_hybrid.csv', index=False)\nprint(f'✅ Labels: {labels.shape}')\n\n# 3. Save pivot (for LSTM)\npivot.to_parquet(f'{OUTPUT_DIR}/pivot_26weeks.parquet')\nprint(f'✅ Pivot: {pivot.shape}')\n\n# 4. Save constants\nconstants = {\n    'FINAL_FEATURES': FINAL_FEATURES,\n    'EARLY_WINDOW': EARLY_WINDOW,\n    'MAX_WEEKS': MAX_WEEKS,\n    'KSHAPE_K': BEST_K,\n    'FAD_CLUSTERS': FAD_CLUSTERS,\n    'FINAL_FAD_RATE': float(labels['fad_label'].mean()),\n    'n_articles': len(export_df),\n    'n_features': len(FINAL_FEATURES),\n}\n\nwith open(f'{OUTPUT_DIR}/pipeline_constants.json', 'w') as f:\n    json.dump(constants, f, indent=2)\nprint(f'✅ Constants: {len(constants)} items')\n\n# 5. Save MI scores\nmi_df.to_csv(f'{OUTPUT_DIR}/mi_scores.csv', index=False)\nprint(f'✅ MI scores')\n\n# Summary\nprint(f'\\n{\"=\"*50}')\nfor f_name in sorted(os.listdir(OUTPUT_DIR)):\n    size = os.path.getsize(f'{OUTPUT_DIR}/{f_name}') / 1024 / 1024\n    print(f'  {f_name:40s} {size:.1f} MB')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:13:10.267658Z","iopub.execute_input":"2026-04-29T08:13:10.268087Z","iopub.status.idle":"2026-04-29T08:13:14.211647Z","shell.execute_reply.started":"2026-04-29T08:13:10.268062Z","shell.execute_reply":"2026-04-29T08:13:14.210891Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 5 — Stacked Ensemble","metadata":{}},{"cell_type":"code","source":"import os\nimport re\nimport json\nimport random\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nimport tensorflow as tf\nfrom tensorflow.keras import layers, models, regularizers, callbacks\n\nfrom sklearn.model_selection import train_test_split, StratifiedKFold\nfrom sklearn.utils.class_weight import compute_class_weight\nfrom sklearn.metrics import (\n    accuracy_score,\n    precision_score,\n    recall_score,\n    f1_score,\n    roc_auc_score,\n    average_precision_score,\n    confusion_matrix,\n    classification_report\n)\n\nSEED = 42\n\nrandom.seed(SEED)\nnp.random.seed(SEED)\ntf.keras.utils.set_random_seed(SEED)\n\nprint(\"TensorFlow:\", tf.__version__)\nprint(\"GPU:\", tf.config.list_physical_devices(\"GPU\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:13:14.212706Z","iopub.execute_input":"2026-04-29T08:13:14.213011Z","iopub.status.idle":"2026-04-29T08:13:14.928502Z","shell.execute_reply.started":"2026-04-29T08:13:14.212980Z","shell.execute_reply":"2026-04-29T08:13:14.927378Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"INPUT_DIR = OUTPUT_DIR\nWORK_DIR = Path(\"/kaggle/working/processed_data\")\n\nall_files = [p for p in WORK_DIR.rglob(\"*\") if p.is_file()]\n\nprint(\"Files found:\")\nfor f in all_files[:100]:\n    print(f)\n\ndef find_file(possible_names):\n    possible_names = [name.lower() for name in possible_names]\n\n    for target_name in possible_names:\n        for path in all_files:\n            if path.name.lower() == target_name:\n                return path\n\n    return None\n\nfeatures_path = find_file([\n    \"features_final.parquet\",\n    \"features_final.csv\",\n    \"features_final.xlsx\"\n])\n\npivot_path = find_file([\n    \"pivot_26weeks.parquet\",\n    \"pivot_26weeks.csv\",\n    \"pivot_26weeks.xlsx\",\n    \"pivot_26_weeks.parquet\",\n    \"pivot_26_weeks.csv\",\n    \"pivot_26_weeks.xlsx\"\n])\n\nconstants_path = find_file([\n    \"pipeline_constants.json\"\n])\n\nprint(\"\\nSelected files:\")\nprint(\"features_path:\", features_path)\nprint(\"pivot_path:\", pivot_path)\nprint(\"constants_path:\", constants_path)\n\nif features_path is None:\n    raise FileNotFoundError(\"Cannot find features_final file.\")\n\nif pivot_path is None:\n    raise FileNotFoundError(\"Cannot find pivot_26weeks file.\")\n\nif constants_path is None:\n    raise FileNotFoundError(\"Cannot find pipeline_constants.json.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:07.971909Z","iopub.execute_input":"2026-04-29T08:15:07.972512Z","iopub.status.idle":"2026-04-29T08:15:07.980809Z","shell.execute_reply.started":"2026-04-29T08:15:07.972482Z","shell.execute_reply":"2026-04-29T08:15:07.980197Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_table(path):\n    path = Path(path)\n    suffix = path.suffix.lower()\n\n    if suffix == \".parquet\":\n        return pd.read_parquet(path)\n    elif suffix == \".csv\":\n        return pd.read_csv(path)\n    elif suffix in [\".xlsx\", \".xls\"]:\n        return pd.read_excel(path)\n    else:\n        raise ValueError(f\"Unsupported file type: {path}\")\n\nfeatures_df = load_table(features_path)\npivot_df = load_table(pivot_path)\n\nwith open(constants_path, \"r\") as f:\n    constants = json.load(f)\n\nprint(\"features_df shape:\", features_df.shape)\nprint(\"pivot_df shape:\", pivot_df.shape)\n\nprint(\"\\nConstants:\")\nprint(constants)\n\nprint(\"\\nfeatures_df columns:\")\nprint(features_df.columns.tolist())\n\nprint(\"\\npivot_df columns:\")\nprint(pivot_df.columns.tolist())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:14.992393Z","iopub.execute_input":"2026-04-29T08:15:14.992816Z","iopub.status.idle":"2026-04-29T08:15:15.108627Z","shell.execute_reply.started":"2026-04-29T08:15:14.992783Z","shell.execute_reply":"2026-04-29T08:15:15.107914Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"FINAL_FEATURES = constants.get(\"FINAL_FEATURES\", None)\nEARLY_WINDOW = constants.get(\"EARLY_WINDOW\", 5)\nMAX_WEEKS = constants.get(\"MAX_WEEKS\", 26)\nSEED = constants.get(\"SEED\", 42)\n\nprint(\"FINAL_FEATURES:\", FINAL_FEATURES)\nprint(\"EARLY_WINDOW:\", EARLY_WINDOW)\nprint(\"MAX_WEEKS:\", MAX_WEEKS)\nprint(\"SEED:\", SEED)\n\nif FINAL_FEATURES is not None:\n    missing_features = [c for c in FINAL_FEATURES if c not in features_df.columns]\n    if len(missing_features) > 0:\n        print(\"Warning: These FINAL_FEATURES are missing from features_df:\")\n        print(missing_features)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:18.712395Z","iopub.execute_input":"2026-04-29T08:15:18.713060Z","iopub.status.idle":"2026-04-29T08:15:18.719168Z","shell.execute_reply.started":"2026-04-29T08:15:18.713029Z","shell.execute_reply":"2026-04-29T08:15:18.718300Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def ensure_article_id(df, df_name):\n    df = df.copy()\n\n    if \"article_id\" in df.columns:\n        return df\n\n    possible_cols = [\n        \"index\", \"article\", \"Article_id\", \"ARTICLE_ID\",\n        \"Unnamed: 0\", \"unnamed: 0\", \"__index_level_0__\"\n    ]\n\n    for col in possible_cols:\n        if col in df.columns:\n            df = df.rename(columns={col: \"article_id\"})\n            return df\n\n    # Important fix: article_id may be stored as DataFrame index\n    if not isinstance(df.index, pd.RangeIndex):\n        df = df.reset_index()\n        first_col = df.columns[0]\n        df = df.rename(columns={first_col: \"article_id\"})\n        return df\n\n    raise ValueError(\n        f\"{df_name} has no article_id column and index is only RangeIndex. \"\n        f\"Columns found: {df.columns.tolist()}\"\n    )\n\n\ndef normalize_article_id(series):\n    s = series.astype(str).str.strip()\n    s = s.str.replace(r\"\\.0$\", \"\", regex=True)\n    s = s.str.replace(r\"\\D\", \"\", regex=True)\n\n    if s.str.len().max() <= 10:\n        s = s.str.zfill(10)\n\n    return s\n\n\nfeatures_df = ensure_article_id(features_df, \"features_df\")\npivot_df = ensure_article_id(pivot_df, \"pivot_df\")\n\nfeatures_df[\"article_id\"] = normalize_article_id(features_df[\"article_id\"])\npivot_df[\"article_id\"] = normalize_article_id(pivot_df[\"article_id\"])\n\nif \"fad_label\" not in features_df.columns:\n    raise ValueError(\"features_final must contain fad_label column.\")\n\n\ndef get_week_number(col):\n    col = str(col).lower()\n    nums = re.findall(r\"\\d+\", col)\n    return int(nums[0]) if nums else 999\n\n\nweek_cols = []\n\nfor col in pivot_df.columns:\n    col_str = str(col).lower()\n\n    if col_str == \"article_id\":\n        continue\n\n    if re.fullmatch(r\"w\\d+\", col_str):\n        week_cols.append(col)\n    elif re.fullmatch(r\"week_?\\d+\", col_str):\n        week_cols.append(col)\n    elif col_str.isdigit():\n        week_cols.append(col)\n\nweek_cols = sorted(week_cols, key=get_week_number)\n\nif len(week_cols) == 0:\n    raise ValueError(\"No week columns detected. Expected columns like w0, w1, ..., w25.\")\n\nif len(week_cols) > MAX_WEEKS:\n    week_cols = week_cols[:MAX_WEEKS]\n\nprint(\"features_df columns after fix:\")\nprint(features_df.columns.tolist()[:10])\n\nprint(\"\\npivot_df columns after fix:\")\nprint(pivot_df.columns.tolist()[:10])\n\nprint(\"\\nDetected week columns:\")\nprint(week_cols)\n\nprint(\"\\nNumber of weeks:\", len(week_cols))\n\nprint(\"\\nSample article_id from features_df:\")\nprint(features_df[\"article_id\"].head())\n\nprint(\"\\nSample article_id from pivot_df:\")\nprint(pivot_df[\"article_id\"].head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:21.246745Z","iopub.execute_input":"2026-04-29T08:15:21.247763Z","iopub.status.idle":"2026-04-29T08:15:21.554232Z","shell.execute_reply.started":"2026-04-29T08:15:21.247724Z","shell.execute_reply":"2026-04-29T08:15:21.553333Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Best CV F1  : 0.9523\nBest params :\n  gamma                  0.49839677992293113\n  learning_rate          0.0795627394576256\n  max_depth              7\n  min_child_weight       5\n  n_estimators           300\n  reg_alpha              0.008954431298712286\n  reg_lambda             0.058325596252984765","metadata":{}},{"cell_type":"code","source":"label_df = features_df[[\"article_id\", \"fad_label\"]].copy()\nlabel_df[\"fad_label\"] = label_df[\"fad_label\"].astype(int)\n\ndata = pivot_df[[\"article_id\"] + week_cols].merge(\n    label_df,\n    on=\"article_id\",\n    how=\"inner\"\n)\n\nprint(\"Merged data shape:\", data.shape)\nprint(\"Fad rate:\", data[\"fad_label\"].mean())\nprint(data[\"fad_label\"].value_counts())\n\nif data.shape[0] == 0:\n    raise ValueError(\"Merge produced 0 rows. Check article_id format in both files.\")\n\nX_seq = data[week_cols].astype(\"float32\").values\nX_seq = np.nan_to_num(X_seq, nan=0.0, posinf=1.0, neginf=0.0)\nX_seq = np.clip(X_seq, 0.0, 1.0)\n\nX_seq = X_seq.reshape(X_seq.shape[0], X_seq.shape[1], 1)\n\ny = data[\"fad_label\"].values.astype(int)\narticle_ids = data[\"article_id\"].values\n\nprint(\"X_seq shape:\", X_seq.shape)\nprint(\"y shape:\", y.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:27.561567Z","iopub.execute_input":"2026-04-29T08:15:27.562077Z","iopub.status.idle":"2026-04-29T08:15:27.641672Z","shell.execute_reply.started":"2026-04-29T08:15:27.562039Z","shell.execute_reply":"2026-04-29T08:15:27.640705Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"idx = np.arange(len(y))\n\ntrain_val_idx, test_idx = train_test_split(\n    idx,\n    test_size=0.15,\n    random_state=SEED,\n    stratify=y\n)\n\ntrain_idx, val_idx = train_test_split(\n    train_val_idx,\n    test_size=0.1765,\n    random_state=SEED,\n    stratify=y[train_val_idx]\n)\n\nX_train = X_seq[train_idx]\ny_train = y[train_idx]\n\nX_val = X_seq[val_idx]\ny_val = y[val_idx]\n\nX_test = X_seq[test_idx]\ny_test = y[test_idx]\n\nprint(\"Train:\", X_train.shape, \"Fad rate:\", y_train.mean())\nprint(\"Val:\", X_val.shape, \"Fad rate:\", y_val.mean())\nprint(\"Test:\", X_test.shape, \"Fad rate:\", y_test.mean())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:29.280535Z","iopub.execute_input":"2026-04-29T08:15:29.281382Z","iopub.status.idle":"2026-04-29T08:15:29.354774Z","shell.execute_reply.started":"2026-04-29T08:15:29.281349Z","shell.execute_reply":"2026-04-29T08:15:29.354122Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"classes = np.array([0, 1])\n\nclass_weights_array = compute_class_weight(\n    class_weight=\"balanced\",\n    classes=classes,\n    y=y_train\n)\n\nclass_weight = {\n    0: float(class_weights_array[0]),\n    1: float(class_weights_array[1])\n}\n\nprint(\"Class weight:\", class_weight)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:31.526610Z","iopub.execute_input":"2026-04-29T08:15:31.527323Z","iopub.status.idle":"2026-04-29T08:15:31.541733Z","shell.execute_reply.started":"2026-04-29T08:15:31.527288Z","shell.execute_reply":"2026-04-29T08:15:31.540734Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def build_lstm_model(n_steps):\n    tf.keras.backend.clear_session()\n    tf.keras.utils.set_random_seed(SEED)\n\n    inputs = layers.Input(shape=(n_steps, 1))\n\n    x = layers.BatchNormalization()(inputs)\n\n    x = layers.LSTM(\n        64,\n        return_sequences=True,\n        dropout=0.20\n    )(x)\n\n    x = layers.LSTM(\n        32,\n        return_sequences=False,\n        dropout=0.20\n    )(x)\n\n    x = layers.Dense(\n        32,\n        activation=\"relu\",\n        kernel_regularizer=regularizers.l2(1e-4)\n    )(x)\n\n    x = layers.Dropout(0.30)(x)\n\n    outputs = layers.Dense(1, activation=\"sigmoid\")(x)\n\n    model = models.Model(inputs, outputs)\n\n    model.compile(\n        optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),\n        loss=\"binary_crossentropy\",\n        metrics=[\n            tf.keras.metrics.AUC(curve=\"PR\", name=\"pr_auc\"),\n            tf.keras.metrics.AUC(curve=\"ROC\", name=\"roc_auc\"),\n            tf.keras.metrics.Precision(name=\"precision\"),\n            tf.keras.metrics.Recall(name=\"recall\")\n        ]\n    )\n\n    return model\n\nmodel = build_lstm_model(X_train.shape[1])\nmodel.summary()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:33.570879Z","iopub.execute_input":"2026-04-29T08:15:33.571320Z","iopub.status.idle":"2026-04-29T08:15:36.033455Z","shell.execute_reply.started":"2026-04-29T08:15:33.571294Z","shell.execute_reply":"2026-04-29T08:15:36.032864Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# LSTM Hyperparameter Tuning — Keras Tuner (RandomSearch)\n# Tunes: LSTM units, dropout, dense units, l2 reg, learning rate, class weight\n# Insert this cell BEFORE the final model training cell.\n\nimport keras_tuner as kt\nimport tensorflow as tf\nfrom tensorflow.keras import layers, models, regularizers, callbacks\nimport numpy as np\n\n# ── Tunable model builder ────────────────────────────────────────────────────\ndef build_tunable_lstm(hp):\n    tf.keras.backend.clear_session()\n    tf.keras.utils.set_random_seed(SEED)\n\n    lstm1_units  = hp.Int('lstm1_units',  min_value=32,  max_value=128, step=32)\n    lstm2_units  = hp.Int('lstm2_units',  min_value=16,  max_value=64,  step=16)\n    lstm_dropout = hp.Float('lstm_dropout', min_value=0.1, max_value=0.5, step=0.1)\n    dense_units  = hp.Int('dense_units',  min_value=16,  max_value=64,  step=16)\n    dense_drop   = hp.Float('dense_dropout', min_value=0.1, max_value=0.5, step=0.1)\n    l2_reg       = hp.Choice('l2_reg', values=[1e-5, 1e-4, 1e-3])\n    lr           = hp.Choice('learning_rate', values=[1e-4, 5e-4, 1e-3])\n\n    n_steps = X_train.shape[1]\n    inputs  = layers.Input(shape=(n_steps, 1))\n    x = layers.BatchNormalization()(inputs)\n    x = layers.LSTM(lstm1_units, return_sequences=True, dropout=lstm_dropout)(x)\n    x = layers.LSTM(lstm2_units, return_sequences=False, dropout=lstm_dropout)(x)\n    x = layers.Dense(dense_units, activation='relu',\n                     kernel_regularizer=regularizers.l2(l2_reg))(x)\n    x = layers.Dropout(dense_drop)(x)\n    outputs = layers.Dense(1, activation='sigmoid')(x)\n\n    model = models.Model(inputs, outputs)\n    model.compile(\n        optimizer=tf.keras.optimizers.Adam(learning_rate=lr),\n        loss='binary_crossentropy',\n        metrics=[\n            tf.keras.metrics.AUC(curve='PR',  name='pr_auc'),\n            tf.keras.metrics.AUC(curve='ROC', name='roc_auc'),\n        ]\n    )\n    return model\n\n# ── Class weight (handle imbalance) ─────────────────────────────────────────\nneg_count = int((y_train == 0).sum())\npos_count = int((y_train == 1).sum())\nclass_weight = {0: 1.0, 1: neg_count / pos_count}\nprint(f'Class weight — 0: {class_weight[0]:.2f}  1: {class_weight[1]:.2f}')\n\n# ── Tuner setup ──────────────────────────────────────────────────────────────\ntuner = kt.RandomSearch(\n    build_tunable_lstm,\n    objective=kt.Objective('val_pr_auc', direction='max'),\n    max_trials=20,\n    executions_per_trial=1,\n    seed=SEED,\n    directory=str(WORK_DIR / 'kt_lstm'),\n    project_name='lstm_fad_tuning',\n    overwrite=True\n)\n\ntuner.search_space_summary()\n\n# ── Callbacks for search ─────────────────────────────────────────────────────\nsearch_callbacks = [\n    callbacks.EarlyStopping(\n        monitor='val_pr_auc', mode='max',\n        patience=5, restore_best_weights=True\n    )\n]\n\n# ── Run search ───────────────────────────────────────────────────────────────\nprint('\\nStarting hyperparameter search (max_trials=20)...')\ntuner.search(\n    X_train, y_train,\n    epochs=30,\n    batch_size=512,\n    validation_data=(X_val, y_val),\n    class_weight=class_weight,\n    callbacks=search_callbacks,\n    verbose=0\n)\n\n# ── Best hyperparameters ─────────────────────────────────────────────────────\nbest_hp = tuner.get_best_hyperparameters(num_trials=1)[0]\nprint('\\nBest hyperparameters:')\nfor param, value in best_hp.values.items():\n    print(f'  {param:<20}: {value}')\n\n# ── Retrain best model to full epochs ────────────────────────────────────────\nTUNED_PARAMS = best_hp.values\n\nlstm_tuned = build_tunable_lstm(best_hp)\n\nfull_callbacks = [\n    callbacks.EarlyStopping(\n        monitor='val_pr_auc', mode='max',\n        patience=10, restore_best_weights=True, verbose=1\n    ),\n    callbacks.ReduceLROnPlateau(\n        monitor='val_pr_auc', mode='max',\n        factor=0.5, patience=5, min_lr=1e-6, verbose=1\n    )\n]\n\nprint('\\nRetraining best model...')\nhistory_tuned = lstm_tuned.fit(\n    X_train, y_train,\n    epochs=100,\n    batch_size=512,\n    validation_data=(X_val, y_val),\n    class_weight=class_weight,\n    callbacks=full_callbacks,\n    verbose=1\n)\n\n# ── Evaluate tuned model vs baseline ────────────────────────────────────────\nfrom sklearn.metrics import f1_score, roc_auc_score, average_precision_score\n\ndef eval_model(m, X, y, label):\n    prob = m.predict(X, batch_size=1024).ravel()\n    # threshold search\n    best_f1, best_thr = 0.0, 0.5\n    for thr in np.arange(0.05, 0.96, 0.01):\n        f1 = f1_score(y, (prob >= thr).astype(int), zero_division=0)\n        if f1 > best_f1:\n            best_f1, best_thr = f1, thr\n    pr_auc  = average_precision_score(y, prob)\n    roc_auc = roc_auc_score(y, prob)\n    print(f'{label:<20}  F1={best_f1:.4f}  PR-AUC={pr_auc:.4f}  ROC-AUC={roc_auc:.4f}  thr={best_thr:.2f}')\n    return prob, best_thr, best_f1\n\nprint('\\nValidation results:')\n_, base_thr, base_f1   = eval_model(model,       X_val, y_val, 'Baseline LSTM')\n_, tuned_thr, tuned_f1 = eval_model(lstm_tuned,  X_val, y_val, 'Tuned LSTM')\n\n# ── Promote tuned model if better ────────────────────────────────────────────\nif tuned_f1 >= base_f1:\n    print('\\nTuned model is better — replacing model for downstream cells.')\n    model = lstm_tuned\n    best_threshold = tuned_thr\nelse:\n    print('\\nBaseline model retained.')\n    best_threshold = base_thr\n\n# ── Save top-3 trial summary ─────────────────────────────────────────────────\nimport json as _json\ntop3 = []\nfor trial in tuner.oracle.get_best_trials(num_trials=3):\n    top3.append({\n        'trial_id': trial.trial_id,\n        'score'   : trial.score,\n        'params'  : trial.hyperparameters.values\n    })\nwith open(WORK_DIR / 'lstm_tuning_results.json', 'w') as _f:\n    _json.dump(top3, _f, indent=4)\nprint('\\nSaved: /kaggle/working/lstm_tuning_results.json')\ntuner.results_summary(num_trials=3)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:15:36.838331Z","iopub.execute_input":"2026-04-29T08:15:36.839013Z","iopub.status.idle":"2026-04-29T08:25:34.378364Z","shell.execute_reply.started":"2026-04-29T08:15:36.838982Z","shell.execute_reply":"2026-04-29T08:25:34.377734Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"early_stop = callbacks.EarlyStopping(\n    monitor=\"val_pr_auc\",\n    mode=\"max\",\n    patience=8,\n    restore_best_weights=True,\n    verbose=1\n)\n\nreduce_lr = callbacks.ReduceLROnPlateau(\n    monitor=\"val_pr_auc\",\n    mode=\"max\",\n    factor=0.5,\n    patience=3,\n    min_lr=1e-6,\n    verbose=1\n)\n\nhistory = model.fit(\n    X_train,\n    y_train,\n    validation_data=(X_val, y_val),\n    epochs=50,\n    batch_size=512,\n    class_weight=class_weight,\n    callbacks=[early_stop, reduce_lr],\n    verbose=1\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:25:34.379701Z","iopub.execute_input":"2026-04-29T08:25:34.380032Z","iopub.status.idle":"2026-04-29T08:25:59.473835Z","shell.execute_reply.started":"2026-04-29T08:25:34.380007Z","shell.execute_reply":"2026-04-29T08:25:59.472969Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"val_prob = model.predict(X_val, batch_size=1024).ravel()\n\nthreshold_rows = []\n\nfor threshold in np.arange(0.05, 0.96, 0.01):\n    val_pred = (val_prob >= threshold).astype(int)\n\n    threshold_rows.append({\n        \"threshold\": threshold,\n        \"f1\": f1_score(y_val, val_pred, zero_division=0),\n        \"precision\": precision_score(y_val, val_pred, zero_division=0),\n        \"recall\": recall_score(y_val, val_pred, zero_division=0),\n        \"accuracy\": accuracy_score(y_val, val_pred)\n    })\n\nthreshold_df = pd.DataFrame(threshold_rows)\n\nbest_row = threshold_df.loc[threshold_df[\"f1\"].idxmax()]\nbest_threshold = float(best_row[\"threshold\"])\n\nprint(\"Best threshold from validation:\")\nprint(best_row)\n\nthreshold_df.sort_values(\"f1\", ascending=False).head(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:25:59.474937Z","iopub.execute_input":"2026-04-29T08:25:59.475592Z","iopub.status.idle":"2026-04-29T08:26:00.251097Z","shell.execute_reply.started":"2026-04-29T08:25:59.475558Z","shell.execute_reply":"2026-04-29T08:26:00.250294Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_prob = model.predict(X_test, batch_size=1024).ravel()\ntest_pred = (test_prob >= best_threshold).astype(int)\n\nmetrics = {\n    \"best_threshold\": best_threshold,\n    \"test_accuracy\": accuracy_score(y_test, test_pred),\n    \"test_precision\": precision_score(y_test, test_pred, zero_division=0),\n    \"test_recall\": recall_score(y_test, test_pred, zero_division=0),\n    \"test_f1\": f1_score(y_test, test_pred, zero_division=0),\n    \"test_pr_auc\": average_precision_score(y_test, test_prob),\n    \"test_roc_auc\": roc_auc_score(y_test, test_prob)\n}\n\nprint(\"Test metrics:\")\nfor k, v in metrics.items():\n    print(f\"{k}: {v:.4f}\" if isinstance(v, float) else f\"{k}: {v}\")\n\nprint(\"\\nConfusion Matrix:\")\nprint(confusion_matrix(y_test, test_pred))\n\nprint(\"\\nClassification Report:\")\nprint(classification_report(y_test, test_pred, digits=4))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:26:00.252718Z","iopub.execute_input":"2026-04-29T08:26:00.253197Z","iopub.status.idle":"2026-04-29T08:26:00.432254Z","shell.execute_reply.started":"2026-04-29T08:26:00.253172Z","shell.execute_reply":"2026-04-29T08:26:00.431656Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model.save(WORK_DIR / \"lstm_fad_model.keras\")\n\nhistory_df = pd.DataFrame(history.history)\n\nval_output = pd.DataFrame({\n    \"article_id\": article_ids[val_idx],\n    \"fad_label\": y_val,\n    \"lstm_prob\": val_prob,\n    \"lstm_pred\": (val_prob >= best_threshold).astype(int)\n})\n\ntest_output = pd.DataFrame({\n    \"article_id\": article_ids[test_idx],\n    \"fad_label\": y_test,\n    \"lstm_prob\": test_prob,\n    \"lstm_pred\": test_pred\n})\n\nall_prob = model.predict(X_seq, batch_size=1024).ravel()\nall_pred = (all_prob >= best_threshold).astype(int)\n\nall_output = pd.DataFrame({\n    \"article_id\": article_ids,\n    \"fad_label\": y,\n    \"lstm_prob\": all_prob,\n    \"lstm_pred\": all_pred\n})\n\nhistory_df.to_csv(WORK_DIR / \"lstm_training_history.csv\", index=False)\nval_output.to_csv(WORK_DIR / \"lstm_val_predictions.csv\", index=False)\ntest_output.to_csv(WORK_DIR / \"lstm_test_predictions.csv\", index=False)\nall_output.to_csv(WORK_DIR / \"lstm_all_predictions.csv\", index=False)\n\nwith open(WORK_DIR / \"lstm_metrics.json\", \"w\") as f:\n    json.dump(metrics, f, indent=4)\n\nprint(\"Saved files:\")\nprint(\"/kaggle/working/lstm_fad_model.keras\")\nprint(\"/kaggle/working/lstm_training_history.csv\")\nprint(\"/kaggle/working/lstm_val_predictions.csv\")\nprint(\"/kaggle/working/lstm_test_predictions.csv\")\nprint(\"/kaggle/working/lstm_all_predictions.csv\")\nprint(\"/kaggle/working/lstm_metrics.json\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:26:00.433216Z","iopub.execute_input":"2026-04-29T08:26:00.433585Z","iopub.status.idle":"2026-04-29T08:26:01.311189Z","shell.execute_reply.started":"2026-04-29T08:26:00.433560Z","shell.execute_reply":"2026-04-29T08:26:01.310492Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"N_SPLITS = 5\n\nskf = StratifiedKFold(\n    n_splits=N_SPLITS,\n    shuffle=True,\n    random_state=SEED\n)\n\noof_prob = np.zeros(len(y), dtype=np.float32)\nfold_metrics = []\n\nfor fold, (tr_idx, va_idx) in enumerate(skf.split(X_seq, y), start=1):\n    print(f\"\\n========== Fold {fold}/{N_SPLITS} ==========\")\n\n    X_tr, y_tr = X_seq[tr_idx], y[tr_idx]\n    X_va, y_va = X_seq[va_idx], y[va_idx]\n\n    class_weights_array = compute_class_weight(\n        class_weight=\"balanced\",\n        classes=np.array([0, 1]),\n        y=y_tr\n    )\n\n    fold_class_weight = {\n        0: float(class_weights_array[0]),\n        1: float(class_weights_array[1])\n    }\n\n    fold_model = build_lstm_model(X_seq.shape[1])\n\n    fold_early_stop = callbacks.EarlyStopping(\n        monitor=\"val_pr_auc\",\n        mode=\"max\",\n        patience=6,\n        restore_best_weights=True,\n        verbose=1\n    )\n\n    fold_reduce_lr = callbacks.ReduceLROnPlateau(\n        monitor=\"val_pr_auc\",\n        mode=\"max\",\n        factor=0.5,\n        patience=3,\n        min_lr=1e-6,\n        verbose=1\n    )\n\n    fold_model.fit(\n        X_tr,\n        y_tr,\n        validation_data=(X_va, y_va),\n        epochs=40,\n        batch_size=512,\n        class_weight=fold_class_weight,\n        callbacks=[fold_early_stop, fold_reduce_lr],\n        verbose=1\n    )\n\n    va_prob = fold_model.predict(X_va, batch_size=1024).ravel()\n    oof_prob[va_idx] = va_prob\n\n    fold_pr_auc = average_precision_score(y_va, va_prob)\n    fold_roc_auc = roc_auc_score(y_va, va_prob)\n\n    fold_metrics.append({\n        \"fold\": fold,\n        \"pr_auc\": fold_pr_auc,\n        \"roc_auc\": fold_roc_auc\n    })\n\n    print(f\"Fold {fold} PR-AUC: {fold_pr_auc:.4f}\")\n    print(f\"Fold {fold} ROC-AUC: {fold_roc_auc:.4f}\")\n\noof_threshold_rows = []\n\nfor threshold in np.arange(0.05, 0.96, 0.01):\n    oof_pred = (oof_prob >= threshold).astype(int)\n\n    oof_threshold_rows.append({\n        \"threshold\": threshold,\n        \"f1\": f1_score(y, oof_pred, zero_division=0),\n        \"precision\": precision_score(y, oof_pred, zero_division=0),\n        \"recall\": recall_score(y, oof_pred, zero_division=0),\n        \"accuracy\": accuracy_score(y, oof_pred)\n    })\n\noof_threshold_df = pd.DataFrame(oof_threshold_rows)\noof_best_row = oof_threshold_df.loc[oof_threshold_df[\"f1\"].idxmax()]\noof_best_threshold = float(oof_best_row[\"threshold\"])\n\noof_pred = (oof_prob >= oof_best_threshold).astype(int)\n\noof_output = pd.DataFrame({\n    \"article_id\": article_ids,\n    \"fad_label\": y,\n    \"lstm_oof_prob\": oof_prob,\n    \"lstm_oof_pred\": oof_pred\n})\n\nfold_metrics_df = pd.DataFrame(fold_metrics)\n\noof_output.to_csv(WORK_DIR / \"lstm_oof_predictions_for_stacking.csv\", index=False)\nfold_metrics_df.to_csv(WORK_DIR / \"lstm_oof_fold_metrics.csv\", index=False)\noof_threshold_df.to_csv(WORK_DIR / \"lstm_oof_threshold_search.csv\", index=False)\n\nprint(\"\\nOOF best threshold:\")\nprint(oof_best_row)\n\nprint(\"\\nOOF Classification Report:\")\nprint(classification_report(y, oof_pred, digits=4))\n\nprint(\"Saved:\")\nprint(\"/kaggle/working/lstm_oof_predictions_for_stacking.csv\")\nprint(\"/kaggle/working/lstm_oof_fold_metrics.csv\")\nprint(\"/kaggle/working/lstm_oof_threshold_search.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:26:01.312324Z","iopub.execute_input":"2026-04-29T08:26:01.312659Z","iopub.status.idle":"2026-04-29T08:28:03.267665Z","shell.execute_reply.started":"2026-04-29T08:26:01.312634Z","shell.execute_reply":"2026-04-29T08:28:03.266893Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import shap\nimport numpy as np\n\n# Get the number of time steps dynamically\nn_steps = X_train.shape[1]\n\n# 1. Create a \"Black Box\" prediction wrapper\ndef lstm_predict_wrapper(data_2d):\n    # Reshape the 2D data from SHAP back into 3D (samples, n_steps, 1) for the LSTM\n    data_3d = data_2d.reshape(-1, n_steps, 1)\n    # Get predictions and flatten them to a 1D array\n    return model.predict(data_3d, verbose=0).flatten()\n\n# 2. Prepare 2D background data\nbackground_3d = X_train[:150]\nbackground_2d = background_3d.reshape(-1, n_steps * 1)\n\n# Pro-tip: shap.kmeans summarizes the background data into 'k' representative samples.\n# This makes KernelExplainer run significantly faster without losing much accuracy.\nbackground_summary = shap.kmeans(background_2d, 25) \n\n# 3. Initialize KernelExplainer\nexplainer = shap.KernelExplainer(lstm_predict_wrapper, background_summary)\n\n# 4. Prepare the test data we want to explain\ndata_to_explain_3d = X_test[:100]\ndata_to_explain_2d = data_to_explain_3d.reshape(-1, n_steps * 1)\n\n# 5. Calculate SHAP values\nprint(\"Calculating SHAP values with KernelExplainer... treating model as a black box.\")\nshap_values_2d = explainer.shap_values(data_to_explain_2d)\n\n# Handle potential list output (depends on SHAP version)\nif isinstance(shap_values_2d, list):\n    shap_values_2d = shap_values_2d[0]\n\n# 6. Reshape the SHAP values back to 3D for your visualization script\nshap_values_3d = np.array(shap_values_2d).reshape(-1, n_steps, 1)\n\nprint(f\"Success! Shape of shap_values_3d: {shap_values_3d.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:35:23.392381Z","iopub.execute_input":"2026-04-29T08:35:23.393140Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1. Calculate importance across the n_steps\nimport matplotlib.pyplot as plt\n# We average axis 0 (samples) and axis 2 (the single feature dimension)\ntime_importance = np.mean(np.abs(shap_values_3d), axis=(0, 2))\nn_steps = len(time_importance)\n\n# 2. Safely handle the labels\nif 'FINAL_FEATURES' in locals() and len(FINAL_FEATURES) == n_steps:\n    labels = FINAL_FEATURES\nelse:\n    print(f\"⚠️ Label Mismatch! Expected {n_steps} names, but FINAL_FEATURES has {len(FINAL_FEATURES) if 'FINAL_FEATURES' in locals() else 0}.\")\n    # Create generic time-step labels: T-25, T-24... T-0\n    labels = [f\"Step {i} (T-{n_steps - 1 - i})\" for i in range(n_steps)]\n\n# 3. Create DataFrame\nimportance_df = pd.DataFrame({\n    'time_step': labels,\n    'importance': time_importance\n}).sort_values(by='importance', ascending=False)\n\n# 4. Plot\nplt.figure(figsize=(10, 8))\nplt.barh(importance_df['time_step'], importance_df['importance'], color='teal')\nplt.gca().invert_yaxis()\nplt.title(f\"Importance over {n_steps} Time Steps\")\nplt.xlabel(\"Mean |SHAP Value|\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 10 Default XGBoost — reuse LSTM train/val/test article_id splits\nimport xgboost as xgb\nimport pandas as pd\nimport numpy as np\nfrom sklearn.metrics import f1_score\n\n# Load features table (pandas)\nxgb_df = load_table(features_path).copy()\nxgb_df = ensure_article_id(xgb_df, 'xgb_df')\nxgb_df['article_id'] = normalize_article_id(xgb_df['article_id'])\n\nfeature_cols = [c for c in xgb_df.columns if c not in ('fad_label', 'article_id')]\n\n# Reuse the LSTM article_id splits (article_ids, train_idx, val_idx, test_idx\n# were defined in the LSTM data-prep cell)\nlstm_train_ids = set(article_ids[train_idx])\nlstm_val_ids   = set(article_ids[val_idx])\nlstm_test_ids  = set(article_ids[test_idx])\n\nxgb_df['_split'] = xgb_df['article_id'].map(\n    lambda aid: 'train' if aid in lstm_train_ids\n           else 'val'   if aid in lstm_val_ids\n           else 'test'  if aid in lstm_test_ids\n           else 'drop'\n)\n\nxgb_train = xgb_df[xgb_df['_split'] == 'train'].reset_index(drop=True)\nxgb_val   = xgb_df[xgb_df['_split'] == 'val'  ].reset_index(drop=True)\nxgb_test  = xgb_df[xgb_df['_split'] == 'test' ].reset_index(drop=True)\n\nprint(f'XGB train: {len(xgb_train):,}  val: {len(xgb_val):,}  test: {len(xgb_test):,}')\n\nX_train_xgb = xgb_train[feature_cols].values.astype('float32')\ny_train_xgb = xgb_train['fad_label'].values.astype(int)\nids_train_xgb = xgb_train['article_id'].values\n\nX_val_xgb = xgb_val[feature_cols].values.astype('float32')\ny_val_xgb = xgb_val['fad_label'].values.astype(int)\nids_val_xgb = xgb_val['article_id'].values\n\nX_test_xgb = xgb_test[feature_cols].values.astype('float32')\ny_test_xgb = xgb_test['fad_label'].values.astype(int)\nids_test_xgb = xgb_test['article_id'].values\n\n# class weight\nneg, pos = np.bincount(y_train_xgb)\nspw = neg / pos\nprint(f'scale_pos_weight: {spw:.4f}')\n\n# Default model\nxgb_default = xgb.XGBClassifier(\n    scale_pos_weight=spw,\n    random_state=SEED,\n    eval_metric='logloss',\n    verbosity=0\n)\nxgb_default.fit(X_train_xgb, y_train_xgb)\n\nf1_default_val  = f1_score(y_val_xgb,  (xgb_default.predict_proba(X_val_xgb)[:,1]  >= 0.5).astype(int), zero_division=0)\nf1_default_test = f1_score(y_test_xgb, (xgb_default.predict_proba(X_test_xgb)[:,1] >= 0.5).astype(int), zero_division=0)\nprint(f'Default model Validation F1: {f1_default_val:.4f}')\nprint(f'Default model Test F1      : {f1_default_test:.4f}')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:31:29.045244Z","iopub.execute_input":"2026-04-29T08:31:29.045941Z","iopub.status.idle":"2026-04-29T08:31:29.900105Z","shell.execute_reply.started":"2026-04-29T08:31:29.045909Z","shell.execute_reply":"2026-04-29T08:31:29.899498Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 11 Manual Tuning XGBoost (no Spark)\nimport xgboost as xgb\n\nBEST_PARAMS = dict(\n    n_estimators     = 500,\n    learning_rate    = 0.05,\n    max_depth        = 8,\n    scale_pos_weight = spw,\n    eval_metric      = 'logloss',\n    tree_method      = 'hist',\n    random_state     = SEED,\n    verbosity        = 0\n)\n\nxgb_tuned = xgb.XGBClassifier(**BEST_PARAMS)\nxgb_tuned.fit(X_train_xgb, y_train_xgb)\n\nf1_tuned_val  = f1_score(y_val_xgb,  (xgb_tuned.predict_proba(X_val_xgb)[:,1]  >= 0.5).astype(int), zero_division=0)\nf1_tuned_test = f1_score(y_test_xgb, (xgb_tuned.predict_proba(X_test_xgb)[:,1] >= 0.5).astype(int), zero_division=0)\nprint(f'Tuned model — Validation F1: {f1_tuned_val:.4f}')\nprint(f'Tuned model — Test F1      : {f1_tuned_test:.4f}')\nprint(f'F1 improvement             : {f1_tuned_test - f1_default_test:+.4f}')\n\n# Best model for downstream steps\nxgb_model = xgb_tuned if f1_tuned_val >= f1_default_val else xgb_default\nprint('Using tuned model.' if xgb_model is xgb_tuned else 'Using default model.')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:31:42.786943Z","iopub.execute_input":"2026-04-29T08:31:42.787425Z","iopub.status.idle":"2026-04-29T08:31:46.613613Z","shell.execute_reply.started":"2026-04-29T08:31:42.787377Z","shell.execute_reply":"2026-04-29T08:31:46.612855Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# OOF generation — 5-fold StratifiedKFold on XGB train set\n# NOTE: xgb_train shares article_ids with the LSTM train set,\n#       so the OOF CSVs can be merged on article_id in the stacking cell.\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import roc_auc_score\nimport xgboost as xgb\n\nN_FOLDS = 5\nskf_xgb = StratifiedKFold(n_splits=N_FOLDS, shuffle=True, random_state=SEED)\n\noof_xgb_prob    = np.zeros(len(y_train_xgb), dtype=np.float32)\noof_article_ids = ids_train_xgb.copy()\n\nfor fold, (tr_idx, va_idx) in enumerate(skf_xgb.split(X_train_xgb, y_train_xgb), start=1):\n    print(f'\\n{\"=\"*48}\\n  OOF Fold {fold} / {N_FOLDS}\\n{\"=\"*48}')\n    fold_model = xgb.XGBClassifier(**BEST_PARAMS)\n    fold_model.fit(X_train_xgb[tr_idx], y_train_xgb[tr_idx])\n    oof_xgb_prob[va_idx] = fold_model.predict_proba(X_train_xgb[va_idx])[:, 1]\n    print(f'  Fold {fold} done — {len(va_idx):,} OOF rows')\n\noof_f1  = f1_score(y_train_xgb, (oof_xgb_prob >= 0.5).astype(int), zero_division=0)\noof_auc = roc_auc_score(y_train_xgb, oof_xgb_prob)\nprint(f'\\nOOF F1  : {oof_f1:.4f}')\nprint(f'OOF AUC : {oof_auc:.4f}')\n\n# Save OOF CSV\nxgb_oof_df = pd.DataFrame({\n    'article_id'  : oof_article_ids,\n    'fad_label'   : y_train_xgb,\n    'xgb_oof_prob': oof_xgb_prob\n})\nxgb_oof_df.to_csv(WORK_DIR / 'xgb_oof.csv', index=False)\nprint(f'Total OOF rows : {len(xgb_oof_df):,}')\nprint('Saved: /kaggle/working/xgb_oof.csv')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:31:52.064684Z","iopub.execute_input":"2026-04-29T08:31:52.065238Z","iopub.status.idle":"2026-04-29T08:32:08.282638Z","shell.execute_reply.started":"2026-04-29T08:31:52.065208Z","shell.execute_reply":"2026-04-29T08:32:08.281841Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 12  Stacking model  (LSTM OOF  +  XGB OOF  →  Logistic-Regression meta-learner)\nimport pandas as pd\nimport numpy as np\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import (\n    f1_score, roc_auc_score, average_precision_score,\n    classification_report, confusion_matrix\n)\n\ndef norm_id(series):\n    \"\"\"Normalize article_id: strip whitespace, drop .0, drop non-digits, zero-pad to 10.\"\"\"\n    s = series.astype(str).str.strip()\n    s = s.str.replace(r'\\.0$', '', regex=True)\n    s = s.str.replace(r'\\D', '', regex=True)\n    if s.str.len().max() <= 10:\n        s = s.str.zfill(10)\n    return s\n\n# ── A. Load OOF files ─────────────────────────────────────────────────────────\nlstm_oof = pd.read_csv(WORK_DIR / 'lstm_oof_predictions_for_stacking.csv')\nxgb_oof  = pd.read_csv(WORK_DIR / 'xgb_oof.csv')\nlstm_oof['article_id'] = norm_id(lstm_oof['article_id'])\nxgb_oof['article_id']  = norm_id(xgb_oof['article_id'])\n\n# ── B. Merge OOF ─────────────────────────────────────────────────────────────\noof_merged = lstm_oof[['article_id', 'fad_label', 'lstm_oof_prob']].merge(\n    xgb_oof[['article_id', 'xgb_oof_prob']],\n    on='article_id', how='inner'\n)\nprint(f'OOF merged rows: {len(oof_merged):,}')\nprint(oof_merged.head())\n\nX_meta_train = oof_merged[['lstm_oof_prob', 'xgb_oof_prob']].values\ny_meta_train = oof_merged['fad_label'].values\n\n# ── C. Build val & test meta-features ────────────────────────────────────────\nlstm_val_preds  = pd.read_csv(WORK_DIR / 'lstm_val_predictions.csv')\nlstm_test_preds = pd.read_csv(WORK_DIR / 'lstm_test_predictions.csv')\nlstm_val_preds['article_id']  = norm_id(lstm_val_preds['article_id'])\nlstm_test_preds['article_id'] = norm_id(lstm_test_preds['article_id'])\n\n# XGB val/test predictions — normalize IDs the same way\nxgb_val_pd = pd.DataFrame({\n    'article_id'   : norm_id(pd.Series(ids_val_xgb)),\n    'xgb_val_prob' : xgb_model.predict_proba(X_val_xgb)[:, 1]\n})\nxgb_test_pd = pd.DataFrame({\n    'article_id'    : norm_id(pd.Series(ids_test_xgb)),\n    'xgb_test_prob' : xgb_model.predict_proba(X_test_xgb)[:, 1]\n})\n\n# Diagnostic: check overlap before merge\nlstm_val_ids_set = set(lstm_val_preds['article_id'])\nxgb_val_ids_set  = set(xgb_val_pd['article_id'])\noverlap_val = lstm_val_ids_set & xgb_val_ids_set\nprint(f'Val  — LSTM: {len(lstm_val_ids_set):,}  XGB: {len(xgb_val_ids_set):,}  overlap: {len(overlap_val):,}')\nif len(overlap_val) == 0:\n    print('  LSTM val sample:', list(lstm_val_ids_set)[:3])\n    print('  XGB  val sample:', list(xgb_val_ids_set)[:3])\n\nlstm_test_ids_set = set(lstm_test_preds['article_id'])\nxgb_test_ids_set  = set(xgb_test_pd['article_id'])\noverlap_test = lstm_test_ids_set & xgb_test_ids_set\nprint(f'Test — LSTM: {len(lstm_test_ids_set):,}  XGB: {len(xgb_test_ids_set):,}  overlap: {len(overlap_test):,}')\n\nval_merged  = lstm_val_preds[['article_id', 'fad_label', 'lstm_prob']].merge(\n    xgb_val_pd, on='article_id', how='inner'\n)\ntest_merged = lstm_test_preds[['article_id', 'fad_label', 'lstm_prob']].merge(\n    xgb_test_pd, on='article_id', how='inner'\n)\nprint(f'val_merged: {len(val_merged):,}   test_merged: {len(test_merged):,}')\n\n# ── Fallback: if merge still empty, build meta-features without article_id join ──\n# Both LSTM and XGB now share the same val/test article sets (from cell 17).\n# If the CSV article_id format still mismatches, build aligned arrays directly\n# from in-memory variables instead of reading from CSV.\nif len(val_merged) == 0:\n    print('WARNING: val merge empty — building meta-features from in-memory arrays')\n    # LSTM val probs are in lstm_val_predictions.csv but also computable from model:\n    # re-predict on X_val (LSTM val set) — X_val was defined in the LSTM data-prep cell\n    lstm_val_prob_arr  = model.predict(X_val,  batch_size=1024).ravel()\n    lstm_test_prob_arr = model.predict(X_test, batch_size=1024).ravel()\n    xgb_val_prob_arr   = xgb_model.predict_proba(X_val_xgb)[:, 1]\n    xgb_test_prob_arr  = xgb_model.predict_proba(X_test_xgb)[:, 1]\n    # Both share same split so same length\n    X_meta_val   = np.column_stack([lstm_val_prob_arr,  xgb_val_prob_arr])\n    y_meta_val   = y_val\n    X_meta_test  = np.column_stack([lstm_test_prob_arr, xgb_test_prob_arr])\n    y_meta_test  = y_test\n    # Reconstruct test_merged for saving\n    test_merged  = pd.DataFrame({\n        'article_id' : article_ids[test_idx],\n        'fad_label'  : y_test\n    })\nelse:\n    X_meta_val   = val_merged[['lstm_prob',  'xgb_val_prob']].values\n    y_meta_val   = val_merged['fad_label'].values\n    X_meta_test  = test_merged[['lstm_prob', 'xgb_test_prob']].values\n    y_meta_test  = test_merged['fad_label'].values\n\n# ── D. Train meta-learner ─────────────────────────────────────────────────────\nmeta_lr = LogisticRegression(C=1.0, max_iter=1000, random_state=SEED, solver='lbfgs')\nmeta_lr.fit(X_meta_train, y_meta_train)\n\nprint('\\nMeta-learner coefficients:')\nprint(f'  LSTM weight : {meta_lr.coef_[0][0]:.4f}')\nprint(f'  XGB  weight : {meta_lr.coef_[0][1]:.4f}')\nprint(f'  Intercept   : {meta_lr.intercept_[0]:.4f}')\n\n# ── E. Threshold search on validation set ────────────────────────────────────\nval_meta_prob = meta_lr.predict_proba(X_meta_val)[:, 1]\nbest_f1_meta, best_thr_meta = 0.0, 0.5\nfor thr in np.arange(0.05, 0.96, 0.01):\n    f1 = f1_score(y_meta_val, (val_meta_prob >= thr).astype(int), zero_division=0)\n    if f1 > best_f1_meta:\n        best_f1_meta, best_thr_meta = f1, thr\nprint(f'\\nBest val threshold : {best_thr_meta:.2f}  (F1={best_f1_meta:.4f})')\n\n# ── F. Evaluate on test set ───────────────────────────────────────────────────\ntest_meta_prob = meta_lr.predict_proba(X_meta_test)[:, 1]\ntest_meta_pred = (test_meta_prob >= best_thr_meta).astype(int)\n\nstack_metrics = {\n    'threshold'    : float(best_thr_meta),\n    'test_f1'      : float(f1_score(y_meta_test, test_meta_pred, zero_division=0)),\n    'test_roc_auc' : float(roc_auc_score(y_meta_test, test_meta_prob)),\n    'test_pr_auc'  : float(average_precision_score(y_meta_test, test_meta_prob)),\n}\n\nprint('\\n' + '='*50)\nprint('  Stacking model — Test results')\nprint('='*50)\nfor k, v in stack_metrics.items():\n    print(f'  {k:<20}: {v:.4f}')\n\nprint('\\nConfusion Matrix:')\nprint(confusion_matrix(y_meta_test, test_meta_pred))\nprint('\\nClassification Report:')\nprint(classification_report(y_meta_test, test_meta_pred, digits=4))\n\n# ── G. Save outputs ───────────────────────────────────────────────────────────\nimport json as _json\n\nstack_test_output = test_merged[['article_id', 'fad_label']].copy()\nstack_test_output['stack_prob'] = test_meta_prob\nstack_test_output['stack_pred'] = test_meta_pred\nstack_test_output.to_csv(WORK_DIR / 'stacking_test_predictions.csv', index=False)\n\nwith open(WORK_DIR / 'stacking_metrics.json', 'w') as _f:\n    _json.dump(stack_metrics, _f, indent=4)\n\nprint('\\nSaved:')\nprint('  /kaggle/working/stacking_test_predictions.csv')\nprint('  /kaggle/working/stacking_metrics.json')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:32:17.405173Z","iopub.execute_input":"2026-04-29T08:32:17.406124Z","iopub.status.idle":"2026-04-29T08:32:18.358137Z","shell.execute_reply.started":"2026-04-29T08:32:17.406082Z","shell.execute_reply":"2026-04-29T08:32:18.357418Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 6 — Cost-Sensitive Threshold Optimization","metadata":{}},{"cell_type":"code","source":"# ── Cell 27: Threshold Sweep (Updated for Stacking) ─────────────────────────\nfrom sklearn.metrics import precision_score, recall_score\nimport matplotlib.pyplot as plt\n\n# 1. Load the merged results from the stacking step\nresults_df = pd.read_csv(WORK_DIR / \"stacking_test_predictions.csv\")\ny_test = results_df[\"fad_label\"].values\nfinal_pred = results_df[\"stack_prob\"].values\n\n# 2. Business Cost Configuration\n# FN (missed Fad) is 10x more expensive than FP (over-flagging a non-fad)\nCOST_FN, COST_FP = 10, 1  \n\nthresholds = np.arange(0.05, 0.96, 0.01)\ncosts, f1s, precs, recs = [], [], [], []\n\nfor t in thresholds:\n    yp = (final_pred >= t).astype(int)\n    tn, fp, fn, tp = confusion_matrix(y_test, yp, labels=[0,1]).ravel()\n    \n    # Calculate weighted cost and standard metrics\n    costs.append(COST_FN*fn + COST_FP*fp)\n    f1s.append(f1_score(y_test, yp, zero_division=0))\n    precs.append(precision_score(y_test, yp, zero_division=0))\n    recs.append(recall_score(y_test, yp, zero_division=0))\n\n# 3. Identify Optimal Threshold\nOPTIMAL_T = thresholds[np.argmin(costs)]\nprint(f'Optimal Business Threshold: {OPTIMAL_T:.2f}')\nprint(f'Metrics at this threshold: F1={f1s[np.argmin(costs)]:.4f}, Recall={recs[np.argmin(costs)]:.4f}')\n\n# 4. Plotting\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Plot 1: Cost Minimization\naxes[0].plot(thresholds, costs, 'coral', lw=2)\naxes[0].axvline(OPTIMAL_T, color='navy', linestyle='--', label=f'Optimal T={OPTIMAL_T:.2f}')\naxes[0].set_xlabel('Threshold')\naxes[0].set_ylabel('Total Weighted Cost')\naxes[0].set_title(f'Business Cost Function (FN={COST_FN}x)', fontweight='bold')\naxes[0].legend()\n\n# Plot 2: Traditional Metrics\naxes[1].plot(thresholds, f1s,   label='F1',        color='black', alpha=0.3, linestyle='--')\naxes[1].plot(thresholds, precs, label='Precision', color='green')\naxes[1].plot(thresholds, recs,  label='Recall',    color='coral')\naxes[1].axvline(OPTIMAL_T, color='navy', linestyle='--')\naxes[1].set_xlabel('Threshold')\naxes[1].set_ylabel('Score')\naxes[1].set_title('Precision vs. Recall Sweep', fontweight='bold')\naxes[1].legend()\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:32:28.156701Z","iopub.execute_input":"2026-04-29T08:32:28.157123Z","iopub.status.idle":"2026-04-29T08:32:29.089249Z","shell.execute_reply.started":"2026-04-29T08:32:28.157093Z","shell.execute_reply":"2026-04-29T08:32:29.088588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 28: Evaluation + Business ROI (Updated for Stacking) ────────────────\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.metrics import (\n    confusion_matrix, classification_report, roc_auc_score, \n    average_precision_score, precision_recall_curve\n)\n\n# 1. Load the results from the stacking model\nresults_df = pd.read_csv(WORK_DIR / \"stacking_test_predictions.csv\")\n\ny_test = results_df[\"fad_label\"].values\nfinal_probs = results_df[\"stack_prob\"].values\n\n# 2. Use the threshold found during the stacking phase\n# (Assuming best_thr_meta is still in memory, otherwise load from stacking_metrics.json)\ny_pred = (final_probs >= best_thr_meta).astype(int)\n\n# 3. Metrics Calculation\ncm = confusion_matrix(y_test, y_pred)\ntn, fp, fn, tp = cm.ravel()\npr_auc = average_precision_score(y_test, final_probs)\nroc_auc = roc_auc_score(y_test, final_probs)\n\nprint('=== Model Evaluation (Stacking Ensemble) ===')\nprint(f'Using Optimal Threshold: {best_thr_meta:.2f}')\nprint(classification_report(y_test, y_pred, target_names=['Non-Fad','Fad']))\nprint(f'PR-AUC : {pr_auc:.4f}')\nprint(f'ROC-AUC: {roc_auc:.4f}')\n\n# 4. Business ROI Calculation\nCOST_PER_FAD_UNIT   = 500\nUNITS_PER_FAD       = 200\nMARKDOWN_RATE       = 0.50\nOPP_LOSS_PCT        = 0.20\nFAD_REVENUE         = 1000\n\n# Logic:\n# Avoided = Correctly identified Fads (TP) that we didn't overstock\n# Opp Lost = Non-Fads flagged as Fads (FP) resulting in lost sales/understock\n# DS Loss  = Fads missed by model (FN) that turned into dead stock\navoided  = tp * UNITS_PER_FAD * COST_PER_FAD_UNIT * MARKDOWN_RATE\nopp_lost = fp * FAD_REVENUE * OPP_LOSS_PCT\nds_loss  = fn * UNITS_PER_FAD * COST_PER_FAD_UNIT * MARKDOWN_RATE\nnet      = avoided - opp_lost - ds_loss\n\nprint(f'\\n=== Business ROI ===')\nprint(f'Markdown avoided : {avoided:>12,.0f} THB')\nprint(f'Opportunity lost : {opp_lost:>12,.0f} THB')\nprint(f'Dead stock loss  : {ds_loss:>12,.0f} THB')\nprint(f'Net Benefit      : {net:>12,.0f} THB  {\"✓ POSITIVE\" if net>0 else \"✗ NEGATIVE\"}')\n\n# 5. Visualization\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Confusion Matrix\nsns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[0],\n            xticklabels=['Non-Fad','Fad'], yticklabels=['Non-Fad','Fad'])\naxes[0].set_title('Confusion Matrix', fontweight='bold')\naxes[0].set_xlabel('Predicted'); axes[0].set_ylabel('Actual')\n\n# PR Curve\nprec_c, rec_c, _ = precision_recall_curve(y_test, final_probs)\naxes[1].plot(rec_c, prec_c, color='steelblue', lw=2, label=f'Stacking (AUC={pr_auc:.3f})')\naxes[1].axhline(y_test.mean(), color='red', linestyle='--', label='Baseline')\naxes[1].set_xlabel('Recall'); axes[1].set_ylabel('Precision')\naxes[1].set_title('Precision-Recall Curve', fontweight='bold')\naxes[1].legend()\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:32:39.928552Z","iopub.execute_input":"2026-04-29T08:32:39.929458Z","iopub.status.idle":"2026-04-29T08:32:40.280598Z","shell.execute_reply.started":"2026-04-29T08:32:39.929386Z","shell.execute_reply":"2026-04-29T08:32:40.279971Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 7 — SHAP + Counterfactual Explainability","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport shap\nimport numpy as np\n\n# 1. Visualize Meta-Learner Weights\nmodel_names = ['LSTM Probability', 'XGBoost Probability']\nweights = meta_lr.coef_[0]\n\nplt.figure(figsize=(8, 4))\ncolors = ['skyblue', 'coral']\nplt.barh(model_names, weights, color=colors)\nplt.axvline(0, color='black', lw=0.8)\nplt.title('Meta-Learner Coefficients (Importance of Base Models)', fontweight='bold')\nplt.xlabel('Coefficient Weight')\nplt.grid(axis='x', linestyle='--', alpha=0.6)\nplt.show()\n\n# 2. SHAP Feature Importance (XGBoost, pandas-based)\nsample_size = min(500, len(X_val_xgb))\nX_sample_features = X_val_xgb[:sample_size]\n\nexplainer  = shap.TreeExplainer(xgb_model)\nshap_values = explainer.shap_values(X_sample_features)\n\n# Handle binary classification output shapes\nsv = shap_values[1] if isinstance(shap_values, list) else shap_values\n\nfeature_names = FINAL_FEATURES if FINAL_FEATURES is not None else feature_cols\n\nplt.figure(figsize=(10, 8))\nshap.summary_plot(\n    sv,\n    X_sample_features,\n    feature_names=feature_names,\n    max_display=20,\n    show=False\n)\nplt.title('SHAP Global Feature Importance (XGBoost Component)', fontweight='bold', fontsize=14)\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:33:00.700508Z","iopub.execute_input":"2026-04-29T08:33:00.700828Z","iopub.status.idle":"2026-04-29T08:33:05.221758Z","shell.execute_reply.started":"2026-04-29T08:33:00.700803Z","shell.execute_reply":"2026-04-29T08:33:05.221156Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 31: Final Summary + Stop Spark ────────────────────\nprint('=' * 55)\nprint('   FINAL PIPELINE SUMMARY')\nprint('=' * 55)\nmetrics = {\n    'F1-Score (Fad)'     : f1_score(y_test, y_pred),\n    'Precision (Fad)'    : precision_score(y_test, y_pred),\n    'Recall (Fad)'       : recall_score(y_test, y_pred),\n    'PR-AUC'             : pr_auc,\n    'ROC-AUC'            : roc_auc_score(y_test, final_pred),\n    'Optimal Threshold'  : OPTIMAL_T,\n    'TP / FP / FN / TN'  : f'{tp} / {fp} / {fn} / {tn}',\n    'Net Business Benefit': f'{net:,.0f} THB'\n}\nfor k, v in metrics.items():\n    print(f'  {k:25s}: {v:.4f}' if isinstance(v, float) else f'  {k:25s}: {v}')\nprint('=' * 55)\n\n# spark.stop()\n# print('Spark session stopped.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:33:11.508869Z","iopub.execute_input":"2026-04-29T08:33:11.509550Z","iopub.status.idle":"2026-04-29T08:33:11.529272Z","shell.execute_reply.started":"2026-04-29T08:33:11.509521Z","shell.execute_reply":"2026-04-29T08:33:11.528612Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════\n#  Business & Financial Metrics\n#  แปลง F1-Score ให้กลายเป็นมูลค่าทางธุรกิจ\n#\n#  ⚠️ ต้องรันหลัง Phase 5 (ต้องมี y_true, y_pred, y_prob)\n#     แต่ตอนนี้ใช้ Fad labels + price data ทำ PRELIMINARY analysis ก่อนได้\n# ══════════════════════════════════════════════════════════════════\n\n# ═══════════════════════════════════════════════════════\n#  Step 0: เตรียมข้อมูลราคา\n# ═══════════════════════════════════════════════════════\n\nprint('📊 Preparing price data...')\n\n# คำนวณราคาเฉลี่ยของแต่ละ article (จาก transactions ดิบ)\nprice_stats = (\n    txn.groupBy('article_id')\n    .agg(\n        F.avg('price').alias('avg_price'),\n        F.min('price').alias('min_price'),\n        F.max('price').alias('max_price'),\n        F.first('price').alias('first_price'),   # ราคาเริ่มต้น (proxy for full price)\n        F.count('*').alias('n_transactions'),\n    )\n    # Markdown detection: ถ้า min_price < first_price * 0.7 → น่าจะถูก markdown\n    .withColumn('price_drop_pct',\n        (F.col('first_price') - F.col('min_price')) / (F.col('first_price') + 0.001))\n    .withColumn('was_marked_down',\n        (F.col('price_drop_pct') > 0.20).cast('int'))  # ลดราคา > 20% ถือว่า markdown\n)\n\nprice_stats.cache()\n\n# Merge กับ labels\nlabels_spark = spark.createDataFrame(\n    labels[['article_id', 'fad_label']].copy()\n)\n\nbiz_data = (\n    price_stats\n    .join(labels_spark, on='article_id', how='inner')\n)\n\nbiz_pd = biz_data.toPandas()\nprint(f'   Articles with price data: {len(biz_pd):,}')\n\n# ═══════════════════════════════════════════════════════\n#  Metric 1: Estimated Markdown Avoidance\n#  \"ถ้าเรารู้ว่าเป็น Fad ล่วงหน้า จะประหยัดเงินจากการ markdown ได้เท่าไหร่\"\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  METRIC 1: Estimated Markdown Avoidance')\nprint('  มูลค่าการลดราคาที่หลีกเลี่ยงได้')\nprint('=' * 70)\n\n# สินค้า Fad ที่ถูก markdown (= Dead Stock ที่เกิดขึ้นจริง)\nfad_marked_down = biz_pd[(biz_pd['fad_label'] == 1) & (biz_pd['was_marked_down'] == 1)]\nfad_not_marked = biz_pd[(biz_pd['fad_label'] == 1) & (biz_pd['was_marked_down'] == 0)]\nnonfad_marked_down = biz_pd[(biz_pd['fad_label'] == 0) & (biz_pd['was_marked_down'] == 1)]\n\nprint(f'\\nFad articles ที่ถูก markdown:     {len(fad_marked_down):,} '\n      f'({len(fad_marked_down)/len(biz_pd[biz_pd[\"fad_label\"]==1])*100:.1f}% ของ Fad)')\nprint(f'Non-Fad articles ที่ถูก markdown:  {len(nonfad_marked_down):,} '\n      f'({len(nonfad_marked_down)/len(biz_pd[biz_pd[\"fad_label\"]==0])*100:.1f}% ของ Non-Fad)')\n\n# คำนวณ markdown cost ต่อ article\n# Markdown cost = (first_price - min_price) × จำนวน transactions ที่ขายในราคาต่ำ\nfad_marked_down = fad_marked_down.copy()\nfad_marked_down['markdown_loss_per_item'] = (\n    fad_marked_down['first_price'] - fad_marked_down['min_price']\n)\nfad_marked_down['est_markdown_loss'] = (\n    fad_marked_down['markdown_loss_per_item'] * \n    fad_marked_down['n_transactions'] * \n    fad_marked_down['price_drop_pct']  # weight by drop severity\n)\n\ntotal_markdown_loss = fad_marked_down['est_markdown_loss'].sum()\navg_markdown_per_article = fad_marked_down['est_markdown_loss'].mean()\n\nprint(f'\\nEstimated Markdown Loss (Fad items):')\nprint(f'   Total:       {total_markdown_loss:>12,.2f} SEK')\nprint(f'   Per article: {avg_markdown_per_article:>12,.2f} SEK')\nprint(f'   Articles:    {len(fad_marked_down):>12,}')\n\n# ถ้าโมเดลจับ Fad ได้ (True Positives) → ประหยัดได้เท่าไหร่\n# สมมติ: ถ้ารู้ล่วงหน้า จะผลิตน้อยลง 50% → ลด markdown loss 50%\nPRODUCTION_REDUCTION = 0.50  # assume ลดผลิต 50% ถ้ารู้ว่าเป็น Fad\n\n# Scenario analysis ตาม Recall ต่าง ๆ\nprint(f'\\n--- Scenario Analysis: Markdown Avoidance by Recall ---')\nprint(f'   (สมมติ: ถ้ารู้ว่าเป็น Fad → ลดผลิต {PRODUCTION_REDUCTION:.0%})')\nprint(f'{\"Recall\":<10s} {\"Fad จับได้\":<15s} {\"Markdown Saved\":<20s} {\"% Saved\":<10s}')\nprint('-' * 55)\n\nfor recall in [0.5, 0.6, 0.7, 0.8, 0.9]:\n    n_caught = int(len(fad_marked_down) * recall)\n    saved = total_markdown_loss * recall * PRODUCTION_REDUCTION\n    pct_saved = saved / total_markdown_loss * 100 if total_markdown_loss > 0 else 0\n    print(f'{recall:<10.0%} {n_caught:<15,} {saved:<20,.2f} {pct_saved:<10.1f}%')\n\n\n# ═══════════════════════════════════════════════════════\n#  Metric 2: Full-Price Sell-Through Rate\n#  \"สัดส่วนสินค้าที่ขายได้ในราคาเต็ม\"\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  METRIC 2: Full-Price Sell-Through Rate (STR)')\nprint('  อัตราส่วนสินค้าที่ขายได้ในราคาเต็ม')\nprint('=' * 70)\n\n# Full-price = ไม่เคยถูก markdown (price_drop < 20%)\nfad_items = biz_pd[biz_pd['fad_label'] == 1]\nnonfad_items = biz_pd[biz_pd['fad_label'] == 0]\n\nfad_full_price_rate = (fad_items['was_marked_down'] == 0).mean()\nnonfad_full_price_rate = (nonfad_items['was_marked_down'] == 0).mean()\noverall_full_price_rate = (biz_pd['was_marked_down'] == 0).mean()\n\nprint(f'\\nFull-Price Sell-Through Rate:')\nprint(f'   Fad items:     {fad_full_price_rate:.1%}')\nprint(f'   Non-Fad items: {nonfad_full_price_rate:.1%}')\nprint(f'   Overall:       {overall_full_price_rate:.1%}')\n\n# ถ้าโมเดลช่วยให้ Fad items มี STR ดีขึ้น\nprint(f'\\n--- Impact Analysis ---')\nprint(f'   ถ้าโมเดลช่วยลด markdown ของ Fad items ได้ 50%:')\nimproved_fad_str = fad_full_price_rate + (1 - fad_full_price_rate) * 0.5\nimproved_overall = (improved_fad_str * len(fad_items) + \n                    nonfad_full_price_rate * len(nonfad_items)) / len(biz_pd)\nprint(f'   Fad STR: {fad_full_price_rate:.1%} → {improved_fad_str:.1%} '\n      f'(+{(improved_fad_str - fad_full_price_rate)*100:.1f}%)')\nprint(f'   Overall STR: {overall_full_price_rate:.1%} → {improved_overall:.1%} '\n      f'(+{(improved_overall - overall_full_price_rate)*100:.1f}%)')\n\n\n# ═══════════════════════════════════════════════════════\n#  Metric 3: Inventory Turnover Ratio (Simulated)\n#  \"อัตราการหมุนเวียนสินค้าคงคลัง\"\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  METRIC 3: Inventory Turnover Ratio (Simulated)')\nprint('  อัตราการหมุนเวียนสินค้าคงคลัง')\nprint('=' * 70)\n\n# ⚠️ ไม่มี inventory data จริง → simulate ด้วย assumption\n# Assumption: initial_stock = peak_weekly_sales × lead_time_weeks\nLEAD_TIME_WEEKS = 4   # supply chain lead time\nCOGS_RATIO = 0.40     # ต้นทุนสินค้า = 40% ของราคาขาย (industry average)\n\n# สำหรับ Fad: ถ้าผลิตตาม peak → จะเหลือ dead stock\n# สำหรับ Non-Fad: ผลิตตาม average → ขายได้หมด\n\nfad_items_sim = fad_items.copy()\nfad_items_sim['est_initial_stock'] = fad_items_sim['avg_price'] * COGS_RATIO * LEAD_TIME_WEEKS\nfad_items_sim['est_cogs_sold'] = fad_items_sim['avg_price'] * COGS_RATIO * fad_items_sim['n_transactions']\n\n# Turnover = COGS sold / average inventory\nfad_turnover = fad_items_sim['n_transactions'].mean() / LEAD_TIME_WEEKS\nnonfad_turnover = nonfad_items['n_transactions'].mean() / LEAD_TIME_WEEKS\n\nprint(f'\\n⚠️  Simulated (ไม่มี inventory data จริง)')\nprint(f'   Assumption: initial stock = peak sales × {LEAD_TIME_WEEKS} weeks lead time')\nprint(f'   Assumption: COGS ratio = {COGS_RATIO:.0%}')\nprint(f'\\nEstimated Inventory Turnover (times per 26 weeks):')\nprint(f'   Fad items:     {fad_turnover:.1f}x')\nprint(f'   Non-Fad items: {nonfad_turnover:.1f}x')\n\nprint(f'\\n--- Impact Analysis ---')\nprint(f'   ถ้าลดผลิต Fad items ลง 50% → dead stock ลดลง')\nimproved_fad_turnover = fad_turnover * 1.5  # ลดสต็อก = turnover เพิ่ม\nprint(f'   Fad turnover: {fad_turnover:.1f}x → {improved_fad_turnover:.1f}x')\n\n\n# ═══════════════════════════════════════════════════════\n#  Metric 4: Cost-Benefit Analysis\n#  \"คุ้มค่าไหมที่ลงทุนพัฒนาระบบ ML\"\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  METRIC 4: Cost-Benefit Analysis')\nprint('  เปรียบเทียบต้นทุนกับผลประโยชน์')\nprint('=' * 70)\n\n# ─── Costs (ค่าใช้จ่าย) ───\n# Google Cloud cluster costs (estimate)\nCLOUD_COST_PER_HOUR = 2.50   # USD per hour (4 × n1-standard-2)\nHOURS_PER_RUN = 1             # ชั่วโมงต่อการรัน pipeline\nRUNS_PER_MONTH = 4            # รันสัปดาห์ละครั้ง\nDEVELOPMENT_HOURS = 200       # ชั่วโมงพัฒนา\nDEV_HOURLY_RATE = 30          # USD per hour\n\nmonthly_cloud_cost = CLOUD_COST_PER_HOUR * HOURS_PER_RUN * RUNS_PER_MONTH\nannual_cloud_cost = monthly_cloud_cost * 12\ndevelopment_cost = DEVELOPMENT_HOURS * DEV_HOURLY_RATE\ntotal_first_year_cost = annual_cloud_cost + development_cost\n\nprint(f'\\n--- Costs ---')\nprint(f'   Cloud (monthly):     ${monthly_cloud_cost:>8,.2f}')\nprint(f'   Cloud (annual):      ${annual_cloud_cost:>8,.2f}')\nprint(f'   Development:         ${development_cost:>8,.2f}')\nprint(f'   Total (Year 1):      ${total_first_year_cost:>8,.2f}')\nprint(f'   Total (Year 2+):     ${annual_cloud_cost:>8,.2f} (no dev cost)')\n\n# ─── Benefits (ผลประโยชน์) ───\n# แปลง SEK → USD (approximate rate)\nSEK_TO_USD = 0.095  # 1 SEK ≈ 0.095 USD\n\n# Assume Recall = 0.70 (reasonable target)\nTARGET_RECALL = 0.70\nestimated_saving_sek = total_markdown_loss * TARGET_RECALL * PRODUCTION_REDUCTION\nestimated_saving_usd = estimated_saving_sek * SEK_TO_USD\n\nprint(f'\\n--- Benefits (at {TARGET_RECALL:.0%} Recall) ---')\nprint(f'   Markdown avoided:    {estimated_saving_sek:>12,.2f} SEK')\nprint(f'   Markdown avoided:    ${estimated_saving_usd:>11,.2f} USD')\n\n# ─── ROI ───\nroi_year1 = (estimated_saving_usd - total_first_year_cost) / total_first_year_cost * 100\nroi_year2 = (estimated_saving_usd - annual_cloud_cost) / annual_cloud_cost * 100\n\nprint(f'\\n--- ROI ---')\nprint(f'   Year 1 ROI: {roi_year1:>+.1f}%')\nprint(f'   Year 2+ ROI: {roi_year2:>+.1f}%')\n\nif roi_year1 > 0:\n    print(f'   ✅ ROI เป็นบวกตั้งแต่ปีแรก — คุ้มค่าลงทุน')\nelif roi_year2 > 0:\n    print(f'   🟡 ROI เป็นบวกตั้งแต่ปีที่ 2 — คุ้มค่าในระยะยาว')\nelse:\n    print(f'   🔴 ROI ยังเป็นลบ — ต้องทบทวน')\n\n# ─── Payback Period ───\nif estimated_saving_usd > 0:\n    payback_months = total_first_year_cost / (estimated_saving_usd / 12)\n    print(f'   Payback period: {payback_months:.1f} months')\n\n\n# ═══════════════════════════════════════════════════════\n#  Summary Table\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  SUMMARY: Business Metrics')\nprint('=' * 70)\n\nsummary = pd.DataFrame([\n    {'Metric': 'Markdown Avoidance (at 70% Recall)',\n     'Value': f'{estimated_saving_sek:,.0f} SEK',\n     'Data Source': 'จริง (price data)'},\n    {'Metric': 'Full-Price STR — Fad items',\n     'Value': f'{fad_full_price_rate:.1%} → {improved_fad_str:.1%}',\n     'Data Source': 'จริง (price data)'},\n    {'Metric': 'Full-Price STR — Overall',\n     'Value': f'{overall_full_price_rate:.1%} → {improved_overall:.1%}',\n     'Data Source': 'จริง (price data)'},\n    {'Metric': 'Inventory Turnover (Fad)',\n     'Value': f'{fad_turnover:.1f}x → {improved_fad_turnover:.1f}x',\n     'Data Source': '⚠️ Simulated'},\n    {'Metric': 'Year 1 ROI',\n     'Value': f'{roi_year1:+.1f}%',\n     'Data Source': 'Estimated'},\n    {'Metric': 'Payback Period',\n     'Value': f'{payback_months:.1f} months' if estimated_saving_usd > 0 else 'N/A',\n     'Data Source': 'Estimated'},\n])\n\nprint('\\n' + summary.to_string(index=False))\n\nprint(f'\\n⚠️  หมายเหตุ:')\nprint(f'   - Markdown Avoidance และ STR ใช้ price data จริงจาก H&M dataset')\nprint(f'   - Inventory Turnover ใช้ simulation (ไม่มี inventory data)')\nprint(f'   - ROI assumptions: cloud ${CLOUD_COST_PER_HOUR}/hr, dev ${DEV_HOURLY_RATE}/hr')\nprint(f'   - ตัวเลขจะ update อีกครั้งหลัง Phase 5 ด้วย actual model predictions')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T08:33:18.296564Z","iopub.execute_input":"2026-04-29T08:33:18.297262Z","iopub.status.idle":"2026-04-29T08:33:50.572938Z","shell.execute_reply.started":"2026-04-29T08:33:18.297231Z","shell.execute_reply":"2026-04-29T08:33:50.572197Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## ✅ Pipeline Complete\n\n**Big Data Architecture:**\n```\n31.7M rows\n    └─ PySpark (lazy + parallel)  ← Phase 0, 1, 3\n           ↓ .toPandas() only when data is small (~105K articles)\n    └─ Pandas / sklearn / LightGBM / XGBoost / LSTM  ← Phase 4-7\n```\n\n**3 values to update from Phase 0 EDA:**\n- `FAD_CLUSTER` (Cell 14) — cluster ID with Fad shape\n- `EARLY_WINDOW` (Cell 15) — week cohort curves diverge\n- `CUTOFF` (Cell 21) — temporal split from imbalance EDA\n- `bright_colours` / `is_allover_pattern` codes (Cell 18)","metadata":{}}]}