{"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-29T14:22:09.484564Z","iopub.execute_input":"2026-04-29T14:22:09.484951Z","iopub.status.idle":"2026-04-29T14:22:09.489719Z","shell.execute_reply.started":"2026-04-29T14:22:09.484924Z","shell.execute_reply":"2026-04-29T14:22:09.488852Z"}},"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-29T14:22:09.491295Z","iopub.execute_input":"2026-04-29T14:22:09.491673Z","iopub.status.idle":"2026-04-29T14:22:12.963421Z","shell.execute_reply.started":"2026-04-29T14:22:09.491645Z","shell.execute_reply":"2026-04-29T14:22:12.962622Z"}},"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-29T14:22:12.964775Z","iopub.execute_input":"2026-04-29T14:22:12.965147Z","iopub.status.idle":"2026-04-29T14:22:12.976078Z","shell.execute_reply.started":"2026-04-29T14:22:12.965106Z","shell.execute_reply":"2026-04-29T14:22:12.975497Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 3: Load Data with PySpark ─────────────────────────\n\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\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-29T14:22:12.977531Z","iopub.execute_input":"2026-04-29T14:22:12.978119Z","iopub.status.idle":"2026-04-29T14:23:13.627250Z","shell.execute_reply.started":"2026-04-29T14:22:12.978095Z","shell.execute_reply":"2026-04-29T14:23:13.626262Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 4: Row Counts ──\n\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-29T14:23:13.628415Z","iopub.execute_input":"2026-04-29T14:23:13.628712Z","iopub.status.idle":"2026-04-29T14:23:19.194853Z","shell.execute_reply.started":"2026-04-29T14:23:13.628683Z","shell.execute_reply":"2026-04-29T14:23:19.193957Z"}},"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 ──────────────────\n\ndef spark_missing_report(df, name):\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\nspark_missing_report(txn.select('t_dat','article_id','customer_id','price'), 'transactions (key cols)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:23:19.195942Z","iopub.execute_input":"2026-04-29T14:23:19.196490Z","iopub.status.idle":"2026-04-29T14:23:54.785705Z","shell.execute_reply.started":"2026-04-29T14:23:19.196458Z","shell.execute_reply":"2026-04-29T14:23:54.784930Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 6: Cardinality Check ──────────────────────\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-29T14:23:54.786753Z","iopub.execute_input":"2026-04-29T14:23:54.787143Z","iopub.status.idle":"2026-04-29T14:23:55.499636Z","shell.execute_reply.started":"2026-04-29T14:23:54.787111Z","shell.execute_reply":"2026-04-29T14:23:55.499005Z"}},"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\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\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\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\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\nprint('\\n=== Sales Quantiles ===')\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-29T14:23:55.500443Z","iopub.execute_input":"2026-04-29T14:23:55.500679Z","iopub.status.idle":"2026-04-29T14:24:56.194946Z","shell.execute_reply.started":"2026-04-29T14:23:55.500658Z","shell.execute_reply":"2026-04-29T14:24:56.194200Z"}},"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# ═══════════════════════════════════════════════════════════════════\n\n# From Phase 0.1 — Dead-on-Arrival handling\nN_DEAD_ON_ARRIVAL = 995 \n\n# From Phase 0.2 — Lognormal + Power Law tail analysis\nPOWER_LAW_XMIN    = 1821     \nPOWER_LAW_ALPHA   = 3.00       \n\n# Phase 1 configuration \nMIN_ACTIVE_WEEKS  = 3       \nCENSORING_WEEKS   = 26       \nMAX_WEEKS_PIVOT   = 26     \n\n# Stockout detection parameters\nSTOCKOUT_BASELINE_WINDOW = 4   \nSTOCKOUT_DROP_THRESHOLD  = 0.30  \nSTOCKOUT_RECOVERY_THRESH = 0.60 \nSTOCKOUT_MIN_BASELINE    = 2   \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-29T14:24:56.197306Z","iopub.execute_input":"2026-04-29T14:24:56.197551Z","iopub.status.idle":"2026-04-29T14:24:56.202843Z","shell.execute_reply.started":"2026-04-29T14:24:56.197530Z","shell.execute_reply":"2026-04-29T14:24:56.202113Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 8: Weekly Aggregation with DENSE time series ──\n\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    .withColumn(\n        'is_censored',\n        (F.col('weeks_from_launch_to_dataset_end') < CENSORING_WEEKS).cast('int')\n    )\n)\n\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) ───\ndense_skeleton = (\n    lifecycle\n    .filter(F.col('is_censored') == 0)  \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\nweekly.cache()\nn_records = weekly.count()\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-29T14:24:56.204037Z","iopub.execute_input":"2026-04-29T14:24:56.204373Z","iopub.status.idle":"2026-04-29T14:26:37.390100Z","shell.execute_reply.started":"2026-04-29T14:24:56.204323Z","shell.execute_reply":"2026-04-29T14:26:37.389191Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 9: Multi-Tier Stockout Detection ──\n\n# ─── Updated Thresholds ───\nSTOCKOUT_BASELINE_WINDOW = 4       \nSTOCKOUT_MIN_BASELINE    = 3.0     \nDROP_SEVERE   = 0.20 \nDROP_MODERATE = 0.50  \nRECOVERY_STRONG  = 0.70  \nRECOVERY_PARTIAL = 0.40 \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\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    .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)\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)\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)\n            (F.col('drop_ratio')     < DROP_MODERATE) &\n            (F.col('baseline_sales') >= 5.0),\n            1\n        ).otherwise(0))\n    \n    # 5. Binary flag\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)  \n         .when(F.col('stockout_tier') == 1, F.col('baseline_sales') * 0.5)  \n         .otherwise(F.col('weekly_sales')))\n)\n\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-29T14:26:37.391416Z","iopub.execute_input":"2026-04-29T14:26:37.391727Z","iopub.status.idle":"2026-04-29T14:26:46.740469Z","shell.execute_reply.started":"2026-04-29T14:26:37.391696Z","shell.execute_reply":"2026-04-29T14:26:46.739577Z"}},"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 ───\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 ───\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        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 ───\ntotal_sales_df = (\n    weekly.groupBy('article_id')\n    .agg(F.sum('weekly_sales_clean').alias('total_sales'))\n    .withColumn('regime',\n        F.when(F.col('total_sales') >= POWER_LAW_XMIN, 'tail')\n         .otherwise('body'))\n)\n\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-29T14:26:46.741505Z","iopub.execute_input":"2026-04-29T14:26:46.741775Z","iopub.status.idle":"2026-04-29T14:27:05.754863Z","shell.execute_reply.started":"2026-04-29T14:26:46.741741Z","shell.execute_reply":"2026-04-29T14:27:05.754022Z"}},"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\nnew_cols = ['article_id'] + [f'w{i}' for i in range(MAX_WEEKS_PIVOT)]\npivot_spark = pivot_spark.toDF(*new_cols)\n\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\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 ───\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')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:27:05.756542Z","iopub.execute_input":"2026-04-29T14:27:05.757035Z","iopub.status.idle":"2026-04-29T14:27:12.669465Z","shell.execute_reply.started":"2026-04-29T14:27:05.757002Z","shell.execute_reply":"2026-04-29T14:27:12.668689Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 2 — Hybrid Labeling (PELT + k-Mean)","metadata":{}},{"cell_type":"code","source":"# ── Cell 12: PELT ──\n\nFAD_MAX_PEAK_WEEK     = 6\nFAD_MIN_DROP_RATIO    = 0.70  \nFAD_MIN_PEAK_RATIO    = 4.0    \nFAD_LOW_TAIL_RATIO    = 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-29T14:27:12.670711Z","iopub.execute_input":"2026-04-29T14:27:12.671035Z","iopub.status.idle":"2026-04-29T14:27:47.586617Z","shell.execute_reply.started":"2026-04-29T14:27:12.671013Z","shell.execute_reply":"2026-04-29T14:27:47.585739Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 13 Distributed Clustering ──\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\n# ═══════════════════════════════════════════════════════\n\nprint('Preparing data in Spark...')\n\nweek_cols = [f'w{i}' for i in range(26)]\npivot_reset = pivot.reset_index()\npivot_spark = spark.createDataFrame(pivot_reset)\n\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\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n🔍 Searching optimal k over range {list(K_SEARCH_RANGE)}...')\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||', \n    )\n    model = kmeans.fit(data_spark)\n    \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\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\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\nkshape_labels_df = clustered.select('article_id', 'cluster').toPandas()\nkshape_labels_df['article_id'] = kshape_labels_df['article_id'].astype(str)\n\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\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\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-29T14:27:47.587658Z","iopub.execute_input":"2026-04-29T14:27:47.587994Z","iopub.status.idle":"2026-04-29T14:29:32.370508Z","shell.execute_reply.started":"2026-04-29T14:27:47.587972Z","shell.execute_reply":"2026-04-29T14:29:32.369802Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════\n#  Cluster Validation & 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\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('=' * 70)\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    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    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    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\nprint('\\n--- 1.2 Cluster Quality Metrics ---\\n')\n\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\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\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\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\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\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\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\n# ══════════════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 2: CLUSTER PROFILING')\nprint('=' * 70)\n\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\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\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\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\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\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\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-29T14:29:32.371545Z","iopub.execute_input":"2026-04-29T14:29:32.372363Z","iopub.status.idle":"2026-04-29T14:29:39.958093Z","shell.execute_reply.started":"2026-04-29T14:29:32.372318Z","shell.execute_reply":"2026-04-29T14:29:39.957208Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 14: Hybrid Label ──\n\nfrom sklearn.metrics import cohen_kappa_score\n\nprint(f'FAD_CLUSTERS = {FAD_CLUSTERS}')\nprint(f'BEST_K = {BEST_K}\\n')\n\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\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\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\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\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\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\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\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\nprint(f'\\n=== Final Hybrid Label Distribution ===')\nprint(labels['fad_label'].value_counts())\nprint(f'Final Fad rate: {hybrid_rate:.2%}')\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')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:29:39.959242Z","iopub.execute_input":"2026-04-29T14:29:39.959557Z","iopub.status.idle":"2026-04-29T14:29:40.077548Z","shell.execute_reply.started":"2026-04-29T14:29:39.959535Z","shell.execute_reply":"2026-04-29T14:29:40.076877Z"}},"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\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\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\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\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\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) \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\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    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    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_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\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\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\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-29T14:29:40.078591Z","iopub.execute_input":"2026-04-29T14:29:40.079311Z","iopub.status.idle":"2026-04-29T14:29:54.209355Z","shell.execute_reply.started":"2026-04-29T14:29:40.079287Z","shell.execute_reply":"2026-04-29T14:29:54.208396Z"}},"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\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-29T14:29:54.210527Z","iopub.execute_input":"2026-04-29T14:29:54.210840Z","iopub.status.idle":"2026-04-29T14:29:57.974094Z","shell.execute_reply.started":"2026-04-29T14:29:54.210808Z","shell.execute_reply":"2026-04-29T14:29:57.973113Z"}},"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) ───\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    .withColumn('interpurchase_cv',\n        F.when(F.col('n_gaps') >= 3, \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-29T14:29:57.975294Z","iopub.execute_input":"2026-04-29T14:29:57.975600Z","iopub.status.idle":"2026-04-29T14:31:34.725501Z","shell.execute_reply.started":"2026-04-29T14:29:57.975568Z","shell.execute_reply":"2026-04-29T14:31:34.724705Z"}},"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',           \n    'department_no',             \n    'graphical_appearance_no',   \n    'colour_group_code',         \n    'perceived_colour_value_id',\n    'perceived_colour_master_id',\n    'index_code',                \n    'index_group_no',\n    'section_no',\n    'garment_group_no',          \n)\n\n# ─── Part B: Target Encoding for high-cardinality features ───\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  \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\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-29T14:31:34.726636Z","iopub.execute_input":"2026-04-29T14:31:34.726997Z","iopub.status.idle":"2026-04-29T14:31:36.510603Z","shell.execute_reply.started":"2026-04-29T14:31:34.726973Z","shell.execute_reply":"2026-04-29T14:31:36.509874Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 19: Merge All Features ──\n\n# ─── Part A: Stockout features per article ───\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        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 ───\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 ───\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')  \n    .join(dyn,              on='article_id', how='left')  \n    .join(gini_spark,       on='article_id', how='left')  \n    .join(buyer_stats.select('article_id', 'total_buyers', \n                             'repeat_purchase_rate', 'avg_purchases_per_buyer'),\n                            on='article_id', how='left')   \n    .join(cv_spark,         on='article_id', how='left')  \n    .join(art_feat_spark,   on='article_id', how='left')   \n    .join(stockout_features,on='article_id', how='left')   \n    .join(regime_spark,     on='article_id', how='left')  \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')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:31:36.511538Z","iopub.execute_input":"2026-04-29T14:31:36.511940Z","iopub.status.idle":"2026-04-29T14:31:51.974114Z","shell.execute_reply.started":"2026-04-29T14:31:36.511914Z","shell.execute_reply":"2026-04-29T14:31:51.973190Z"}},"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 ───\nDROP_COLS = [\n    'article_id', 'fad_label',             \n    'regime', 'strat_key', 'total_sales',   \n]\n\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\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\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\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%}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:31:51.975353Z","iopub.execute_input":"2026-04-29T14:31:51.975781Z","iopub.status.idle":"2026-04-29T14:32:24.079702Z","shell.execute_reply.started":"2026-04-29T14:31:51.975747Z","shell.execute_reply":"2026-04-29T14:32:24.079124Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 20.5: Export ──\nimport os, json\n\nOUTPUT_DIR = '/kaggle/working/processed_data'\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\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-29T14:32:24.083761Z","iopub.execute_input":"2026-04-29T14:32:24.084194Z","iopub.status.idle":"2026-04-29T14:32:28.377998Z","shell.execute_reply.started":"2026-04-29T14:32:24.084169Z","shell.execute_reply":"2026-04-29T14:32:28.377236Z"}},"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-29T14:32:28.378884Z","iopub.execute_input":"2026-04-29T14:32:28.379231Z","iopub.status.idle":"2026-04-29T14:32:29.345233Z","shell.execute_reply.started":"2026-04-29T14:32:28.379192Z","shell.execute_reply":"2026-04-29T14:32:29.344263Z"}},"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-29T14:32:29.346346Z","iopub.execute_input":"2026-04-29T14:32:29.346701Z","iopub.status.idle":"2026-04-29T14:32:29.365720Z","shell.execute_reply.started":"2026-04-29T14:32:29.346679Z","shell.execute_reply":"2026-04-29T14:32:29.365184Z"}},"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-29T14:32:29.366643Z","iopub.execute_input":"2026-04-29T14:32:29.366987Z","iopub.status.idle":"2026-04-29T14:32:29.502827Z","shell.execute_reply.started":"2026-04-29T14:32:29.366955Z","shell.execute_reply":"2026-04-29T14:32:29.502090Z"}},"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-29T14:32:29.503727Z","iopub.execute_input":"2026-04-29T14:32:29.504031Z","iopub.status.idle":"2026-04-29T14:32:29.509983Z","shell.execute_reply.started":"2026-04-29T14:32:29.504010Z","shell.execute_reply":"2026-04-29T14:32:29.509155Z"}},"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-29T14:32:29.510979Z","iopub.execute_input":"2026-04-29T14:32:29.511219Z","iopub.status.idle":"2026-04-29T14:32:29.779594Z","shell.execute_reply.started":"2026-04-29T14:32:29.511199Z","shell.execute_reply":"2026-04-29T14:32:29.778827Z"}},"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-29T14:32:29.780583Z","iopub.execute_input":"2026-04-29T14:32:29.780932Z","iopub.status.idle":"2026-04-29T14:32:29.851551Z","shell.execute_reply.started":"2026-04-29T14:32:29.780902Z","shell.execute_reply":"2026-04-29T14:32:29.850907Z"}},"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-29T14:32:29.852544Z","iopub.execute_input":"2026-04-29T14:32:29.852981Z","iopub.status.idle":"2026-04-29T14:32:29.916115Z","shell.execute_reply.started":"2026-04-29T14:32:29.852957Z","shell.execute_reply":"2026-04-29T14:32:29.915447Z"}},"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-29T14:32:29.917053Z","iopub.execute_input":"2026-04-29T14:32:29.917381Z","iopub.status.idle":"2026-04-29T14:32:29.930218Z","shell.execute_reply.started":"2026-04-29T14:32:29.917358Z","shell.execute_reply":"2026-04-29T14:32:29.929463Z"}},"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-29T14:32:29.931286Z","iopub.execute_input":"2026-04-29T14:32:29.931523Z","iopub.status.idle":"2026-04-29T14:32:32.825884Z","shell.execute_reply.started":"2026-04-29T14:32:29.931502Z","shell.execute_reply":"2026-04-29T14:32:32.825259Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# LSTM Hyperparameter Tuning — Keras Tuner (RandomSearch)\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-29T14:32:32.826846Z","iopub.execute_input":"2026-04-29T14:32:32.827148Z","iopub.status.idle":"2026-04-29T14:42:13.555582Z","shell.execute_reply.started":"2026-04-29T14:32:32.827128Z","shell.execute_reply":"2026-04-29T14:42:13.554714Z"}},"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-29T14:42:13.556637Z","iopub.execute_input":"2026-04-29T14:42:13.557026Z","iopub.status.idle":"2026-04-29T14:42:37.241894Z","shell.execute_reply.started":"2026-04-29T14:42:13.557000Z","shell.execute_reply":"2026-04-29T14:42:37.240995Z"}},"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-29T14:42:37.242950Z","iopub.execute_input":"2026-04-29T14:42:37.243605Z","iopub.status.idle":"2026-04-29T14:42:37.967667Z","shell.execute_reply.started":"2026-04-29T14:42:37.243579Z","shell.execute_reply":"2026-04-29T14:42:37.966914Z"}},"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-29T14:42:37.968555Z","iopub.execute_input":"2026-04-29T14:42:37.968917Z","iopub.status.idle":"2026-04-29T14:42:38.143144Z","shell.execute_reply.started":"2026-04-29T14:42:37.968853Z","shell.execute_reply":"2026-04-29T14:42:38.142253Z"}},"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-29T14:42:38.144206Z","iopub.execute_input":"2026-04-29T14:42:38.144474Z","iopub.status.idle":"2026-04-29T14:42:39.018742Z","shell.execute_reply.started":"2026-04-29T14:42:38.144451Z","shell.execute_reply":"2026-04-29T14:42:39.017935Z"}},"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-29T14:42:39.019798Z","iopub.execute_input":"2026-04-29T14:42:39.020219Z","iopub.status.idle":"2026-04-29T14:44:38.936184Z","shell.execute_reply.started":"2026-04-29T14:42:39.020184Z","shell.execute_reply":"2026-04-29T14:44:38.935347Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import shap\nimport numpy as np\n\nn_steps = X_train.shape[1]\n\n# 1. Create a \"Black Box\" prediction wrapper\ndef lstm_predict_wrapper(data_2d):\n    data_3d = data_2d.reshape(-1, n_steps, 1)\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\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\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-29T14:44:38.937281Z","iopub.execute_input":"2026-04-29T14:44:38.937591Z","iopub.status.idle":"2026-04-29T14:52:22.981668Z","shell.execute_reply.started":"2026-04-29T14:44:38.937568Z","shell.execute_reply":"2026-04-29T14:52:22.980814Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1. Calculate importance across the n_steps\nimport matplotlib.pyplot as plt\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,"execution":{"iopub.status.busy":"2026-04-29T14:52:22.982802Z","iopub.execute_input":"2026-04-29T14:52:22.983239Z","iopub.status.idle":"2026-04-29T14:52:23.249325Z","shell.execute_reply.started":"2026-04-29T14:52:22.983198Z","shell.execute_reply":"2026-04-29T14:52:23.248496Z"}},"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\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-29T14:52:23.250410Z","iopub.execute_input":"2026-04-29T14:52:23.250735Z","iopub.status.idle":"2026-04-29T14:52:24.138928Z","shell.execute_reply.started":"2026-04-29T14:52:23.250710Z","shell.execute_reply":"2026-04-29T14:52:24.138431Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 11 Manual Tuning XGBoost\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-29T14:52:24.139540Z","iopub.execute_input":"2026-04-29T14:52:24.139744Z","iopub.status.idle":"2026-04-29T14:52:27.817516Z","shell.execute_reply.started":"2026-04-29T14:52:24.139722Z","shell.execute_reply":"2026-04-29T14:52:27.816951Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# OOF generation — 5-fold StratifiedKFold on XGB train set\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-29T14:52:27.818659Z","iopub.execute_input":"2026-04-29T14:52:27.818908Z","iopub.status.idle":"2026-04-29T14:52:43.726880Z","shell.execute_reply.started":"2026-04-29T14:52:27.818885Z","shell.execute_reply":"2026-04-29T14:52:43.726157Z"}},"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\nif len(val_merged) == 0:\n    print('WARNING: val merge empty — building meta-features from in-memory arrays')\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\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\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-29T14:52:43.727966Z","iopub.execute_input":"2026-04-29T14:52:43.728288Z","iopub.status.idle":"2026-04-29T14:52:44.624743Z","shell.execute_reply.started":"2026-04-29T14:52:43.728265Z","shell.execute_reply":"2026-04-29T14:52:44.624088Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 6 — Cost-Sensitive Threshold Optimization","metadata":{}},{"cell_type":"code","source":"# ── Cell 27: Threshold Sweep ─────────────────────────\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-29T14:52:44.625864Z","iopub.execute_input":"2026-04-29T14:52:44.626170Z","iopub.status.idle":"2026-04-29T14:52:45.539804Z","shell.execute_reply.started":"2026-04-29T14:52:44.626147Z","shell.execute_reply":"2026-04-29T14:52:45.539178Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 28: Evaluation + Business Metrics (Real Price Data) ──\n\nimport pandas as pd\nimport numpy as np\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    f1_score, precision_score, recall_score\n)\n\n# ═══════════════════════════════════════════════════════\n#  PART 1: Model Evaluation\n# ═══════════════════════════════════════════════════════\n\nresults_df = pd.read_csv(WORK_DIR / \"stacking_test_predictions.csv\")\ny_test = results_df[\"fad_label\"].values\nfinal_probs = results_df[\"stack_prob\"].values\n\n# ใช้ threshold จาก Cell 27\ny_pred = (final_probs >= OPTIMAL_T).astype(int)\n\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('=' * 70)\nprint('  PART 1: Model Evaluation (Stacking Ensemble)')\nprint('=' * 70)\nprint(f'\\nOptimal Threshold: {OPTIMAL_T:.2f}')\nprint(f'Confusion Matrix: TP={tp:,}, FP={fp:,}, FN={fn:,}, TN={tn:,}')\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\nactual_recall = tp / (tp + fn) if (tp + fn) > 0 else 0\nactual_precision = tp / (tp + fp) if (tp + fp) > 0 else 0\nactual_f1 = f1_score(y_test, y_pred)\n\n# ═══════════════════════════════════════════════════════\n#  PART 2: Business Metrics (Real Price Data)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 2: Business Metrics (จาก Price Data จริง)')\nprint('=' * 70)\n\n# ─── 2.0 เตรียมข้อมูลราคา ───\nprint('\\n--- Preparing price data ---')\n\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'),\n        F.count('*').alias('n_transactions'),\n    )\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'))\n    .withColumn('markdown_loss_per_item',\n        F.col('first_price') - F.col('min_price'))\n)\n\nprice_pd = price_stats.toPandas()\nprice_pd['article_id'] = price_pd['article_id'].astype(str)\n\n# Merge predictions กับ price data\ntest_articles = results_df[['article_id', 'fad_label']].copy()\ntest_articles['article_id'] = test_articles['article_id'].astype(str)\ntest_articles['y_pred'] = y_pred\ntest_articles['y_prob'] = final_probs\n\n# Classify: TP, FP, FN, TN\ntest_articles['outcome'] = 'TN'\ntest_articles.loc[(test_articles['fad_label'] == 1) & (test_articles['y_pred'] == 1), 'outcome'] = 'TP'\ntest_articles.loc[(test_articles['fad_label'] == 0) & (test_articles['y_pred'] == 1), 'outcome'] = 'FP'\ntest_articles.loc[(test_articles['fad_label'] == 1) & (test_articles['y_pred'] == 0), 'outcome'] = 'FN'\n\nbiz = test_articles.merge(price_pd, on='article_id', how='inner')\nprint(f'   Test articles with price data: {len(biz):,}')\nprint(f'   TP={tp:,}, FP={fp:,}, FN={fn:,}, TN={tn:,}')\n\n# ─── 2.1 Metric 1: Markdown Avoidance ───\nprint('\\n' + '-' * 50)\nprint('  Metric 1: Estimated Markdown Avoidance')\nprint('-' * 50)\n\n# TP: Fad ที่จับได้ถูก → ถ้าลดผลิตจะ avoid markdown\ntp_data = biz[biz['outcome'] == 'TP']\ntp_marked_down = tp_data[tp_data['was_marked_down'] == 1]\n\n# FN: Fad ที่พลาด → ยังเกิด markdown loss\nfn_data = biz[biz['outcome'] == 'FN']\nfn_marked_down = fn_data[fn_data['was_marked_down'] == 1]\n\nPRODUCTION_REDUCTION = 0.50\n\n# Markdown ที่ avoid ได้ (จาก TP)\ntp_markdown_loss = (\n    tp_marked_down['markdown_loss_per_item'] *\n    tp_marked_down['n_transactions'] *\n    tp_marked_down['price_drop_pct']\n).sum()\navoided_markdown = tp_markdown_loss * PRODUCTION_REDUCTION\n\n# Markdown ที่ยังเกิด (จาก FN — missed Fads)\nfn_markdown_loss = (\n    fn_marked_down['markdown_loss_per_item'] *\n    fn_marked_down['n_transactions'] *\n    fn_marked_down['price_drop_pct']\n).sum()\n\nprint(f'\\nTP (Fad จับได้): {len(tp_data):,} articles')\nprint(f'   TP ที่ถูก markdown:  {len(tp_marked_down):,} ({len(tp_marked_down)/max(1,len(tp_data))*100:.1f}%)')\nprint(f'   Markdown loss (TP):  {tp_markdown_loss:>12,.2f} SEK')\nprint(f'   Avoided (ลดผลิต {PRODUCTION_REDUCTION:.0%}): {avoided_markdown:>12,.2f} SEK')\n\nprint(f'\\nFN (Fad ที่พลาด): {len(fn_data):,} articles')\nprint(f'   FN ที่ถูก markdown:  {len(fn_marked_down):,} ({len(fn_marked_down)/max(1,len(fn_data))*100:.1f}%)')\nprint(f'   Remaining loss (FN): {fn_markdown_loss:>12,.2f} SEK')\n\nprint(f'\\nNet Markdown Improvement: {avoided_markdown:,.2f} SEK')\n\n# ─── 2.2 Metric 2: Full-Price STR ───\nprint('\\n' + '-' * 50)\nprint('  Metric 2: Full-Price Sell-Through Rate')\nprint('-' * 50)\n\n# Current state (ไม่มีโมเดล)\nfad_test = biz[biz['fad_label'] == 1]\nnonfad_test = biz[biz['fad_label'] == 0]\n\ncurrent_fad_str = (fad_test['was_marked_down'] == 0).mean()\ncurrent_nonfad_str = (nonfad_test['was_marked_down'] == 0).mean()\ncurrent_overall_str = (biz['was_marked_down'] == 0).mean()\n\n# Improved state (โมเดลช่วย)\n# TP: ลดผลิต → markdown ลดลง 50%\n# FN: ไม่รู้ → markdown เท่าเดิม\ntp_current_markdown_rate = tp_data['was_marked_down'].mean() if len(tp_data) > 0 else 0\ntp_improved_markdown_rate = tp_current_markdown_rate * (1 - PRODUCTION_REDUCTION)\n\nn_fad_test = len(fad_test)\nn_tp = len(tp_data)\nn_fn = len(fn_data)\n\n# Improved Fad STR\nif n_fad_test > 0:\n    fn_str = (fn_data['was_marked_down'] == 0).mean() if len(fn_data) > 0 else current_fad_str\n    tp_improved_str = 1 - tp_improved_markdown_rate\n    improved_fad_str = (tp_improved_str * n_tp + fn_str * n_fn) / n_fad_test\nelse:\n    improved_fad_str = current_fad_str\n\nimproved_overall_str = (\n    improved_fad_str * n_fad_test + current_nonfad_str * len(nonfad_test)\n) / len(biz)\n\nprint(f'\\n{\"Category\":<25s} {\"Current\":<12s} {\"With Model\":<12s} {\"Change\":<10s}')\nprint('-' * 60)\nprint(f'{\"Fad items\":<25s} {current_fad_str:<12.1%} {improved_fad_str:<12.1%} '\n      f'+{(improved_fad_str - current_fad_str)*100:.1f}pp')\nprint(f'{\"Non-Fad items\":<25s} {current_nonfad_str:<12.1%} {current_nonfad_str:<12.1%} '\n      f'  0.0pp')\nprint(f'{\"Overall\":<25s} {current_overall_str:<12.1%} {improved_overall_str:<12.1%} '\n      f'+{(improved_overall_str - current_overall_str)*100:.1f}pp')\n\n# ─── 2.3 Metric 3: Inventory Turnover (Simulated) ───\nprint('\\n' + '-' * 50)\nprint('  Metric 3: Inventory Turnover (Simulated)')\nprint('-' * 50)\n\nLEAD_TIME_WEEKS = 4\n\nfad_avg_txn = fad_test['n_transactions'].mean()\nnonfad_avg_txn = nonfad_test['n_transactions'].mean()\n\ncurrent_fad_turnover = fad_avg_txn / LEAD_TIME_WEEKS\ncurrent_nonfad_turnover = nonfad_avg_txn / LEAD_TIME_WEEKS\n\n# TP: ลดสต็อก 50% → turnover เพิ่ม 2x\n# FN: เท่าเดิม\nimproved_fad_turnover = (\n    (current_fad_turnover / (1 - PRODUCTION_REDUCTION)) * (n_tp / max(1, n_fad_test)) +\n    current_fad_turnover * (n_fn / max(1, n_fad_test))\n)\n\nprint(f'\\n⚠️  Simulated (ไม่มี inventory data จริง)')\nprint(f'\\n{\"Category\":<25s} {\"Current\":<12s} {\"With Model\":<12s}')\nprint('-' * 50)\nprint(f'{\"Fad items\":<25s} {current_fad_turnover:<12.1f}x {improved_fad_turnover:<12.1f}x')\nprint(f'{\"Non-Fad items\":<25s} {current_nonfad_turnover:<12.1f}x {current_nonfad_turnover:<12.1f}x')\n\n# ─── 2.4 Metric 4: Cost-Benefit Analysis ───\nprint('\\n' + '-' * 50)\nprint('  Metric 4: Cost-Benefit Analysis')\nprint('-' * 50)\n\n# Costs\nCLOUD_COST_PER_HOUR = 2.50\nHOURS_PER_RUN = 1\nRUNS_PER_MONTH = 4\nDEVELOPMENT_HOURS = 200\nDEV_HOURLY_RATE = 30\nSEK_TO_USD = 0.095\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_year1_cost = annual_cloud_cost + development_cost\n\n# Benefits (จาก actual model performance)\nbenefit_sek = avoided_markdown\nbenefit_usd = benefit_sek * SEK_TO_USD\n\n# FP cost: opportunity loss จากการลดผลิต Non-Fad ที่ถูก flag ผิด\nfp_data = biz[biz['outcome'] == 'FP']\nfp_opp_loss_sek = (fp_data['avg_price'] * fp_data['n_transactions'] * 0.10).sum()\nfp_opp_loss_usd = fp_opp_loss_sek * SEK_TO_USD\n\nnet_benefit_usd = benefit_usd - fp_opp_loss_usd\n\nroi_year1 = (net_benefit_usd - total_year1_cost) / total_year1_cost * 100\nroi_year2 = (net_benefit_usd - annual_cloud_cost) / annual_cloud_cost * 100\n\nprint(f'\\n--- Costs ---')\nprint(f'   Cloud (annual):      ${annual_cloud_cost:>10,.2f}')\nprint(f'   Development:         ${development_cost:>10,.2f}')\nprint(f'   Total Year 1:        ${total_year1_cost:>10,.2f}')\n\nprint(f'\\n--- Benefits (Actual Model at Recall={actual_recall:.0%}) ---')\nprint(f'   Markdown avoided:    {benefit_sek:>12,.2f} SEK (${benefit_usd:>10,.2f})')\nprint(f'   FP opportunity loss: {fp_opp_loss_sek:>12,.2f} SEK (${fp_opp_loss_usd:>10,.2f})')\nprint(f'   Net benefit:         ${net_benefit_usd:>10,.2f}')\n\nprint(f'\\n--- ROI ---')\nprint(f'   Year 1:  {roi_year1:>+.1f}%')\nprint(f'   Year 2+: {roi_year2:>+.1f}%')\n\nif net_benefit_usd > 0 and net_benefit_usd > total_year1_cost:\n    print(f'   ✅ ROI เป็นบวกตั้งแต่ปีแรก')\nelif net_benefit_usd > 0:\n    payback_months = total_year1_cost / (net_benefit_usd / 12)\n    print(f'   🟡 Payback period: {payback_months:.1f} months')\nelse:\n    print(f'   ⚠️ Net benefit ติดลบ — ดู caveat เรื่อง encoded price')\n\n# ═══════════════════════════════════════════════════════\n#  PART 3: Visualization\n# ═══════════════════════════════════════════════════════\n\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\n\n# 3.1 Confusion Matrix\nsns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[0, 0],\n            xticklabels=['Non-Fad', 'Fad'], yticklabels=['Non-Fad', 'Fad'])\naxes[0, 0].set_title(f'Confusion Matrix (T={OPTIMAL_T:.2f})', fontweight='bold')\naxes[0, 0].set_xlabel('Predicted')\naxes[0, 0].set_ylabel('Actual')\n\n# 3.2 PR Curve\nprec_c, rec_c, _ = precision_recall_curve(y_test, final_probs)\naxes[0, 1].plot(rec_c, prec_c, color='steelblue', lw=2,\n                label=f'Stacking (PR-AUC={pr_auc:.3f})')\naxes[0, 1].axhline(y_test.mean(), color='red', linestyle='--', label='Baseline')\naxes[0, 1].set_xlabel('Recall')\naxes[0, 1].set_ylabel('Precision')\naxes[0, 1].set_title('Precision-Recall Curve', fontweight='bold')\naxes[0, 1].legend()\naxes[0, 1].grid(alpha=0.3)\n\n# 3.3 Business Impact: Markdown Avoidance by Outcome\noutcome_counts = biz.groupby('outcome').agg(\n    n_articles=('article_id', 'count'),\n    n_marked_down=('was_marked_down', 'sum'),\n    total_markdown=('markdown_loss_per_item', lambda x: (x * biz.loc[x.index, 'n_transactions'] * biz.loc[x.index, 'price_drop_pct']).sum())\n).reindex(['TP', 'FP', 'FN', 'TN'])\n\ncolors_map = {'TP': 'green', 'FP': 'orange', 'FN': 'red', 'TN': 'steelblue'}\nbars = axes[1, 0].bar(\n    outcome_counts.index,\n    outcome_counts['n_articles'],\n    color=[colors_map.get(x, 'gray') for x in outcome_counts.index],\n    alpha=0.7, edgecolor='gray'\n)\naxes[1, 0].set_title('Articles by Prediction Outcome', fontweight='bold')\naxes[1, 0].set_ylabel('Count')\nfor bar, val in zip(bars, outcome_counts['n_articles']):\n    axes[1, 0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 50,\n                    f'{val:,}', ha='center', fontsize=10)\naxes[1, 0].grid(alpha=0.3, axis='y')\n\n# 3.4 STR Comparison\nstr_data = pd.DataFrame({\n    'Category': ['Fad\\n(Current)', 'Fad\\n(With Model)', 'Non-Fad', 'Overall\\n(Current)', 'Overall\\n(With Model)'],\n    'STR': [current_fad_str, improved_fad_str, current_nonfad_str, current_overall_str, improved_overall_str],\n    'Color': ['#ff7675', '#00b894', '#74b9ff', '#b2bec3', '#00b894']\n})\nbars = axes[1, 1].bar(str_data['Category'], str_data['STR'] * 100, color=str_data['Color'], alpha=0.8, edgecolor='gray')\naxes[1, 1].set_ylabel('Full-Price STR (%)')\naxes[1, 1].set_title('Full-Price Sell-Through Rate', fontweight='bold')\nfor bar, val in zip(bars, str_data['STR']):\n    axes[1, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,\n                    f'{val:.1%}', ha='center', fontsize=10)\naxes[1, 1].grid(alpha=0.3, axis='y')\n\nplt.suptitle(f'Phase 6: Model Evaluation & Business Impact\\n'\n             f'(Recall={actual_recall:.1%}, Precision={actual_precision:.1%}, F1={actual_f1:.4f})',\n             fontweight='bold', fontsize=13)\nplt.tight_layout()\nplt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Summary\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  SUMMARY')\nprint('=' * 70)\n\nsummary = pd.DataFrame([\n    {'Metric': 'Model — F1 Score', 'Value': f'{actual_f1:.4f}'},\n    {'Metric': 'Model — Recall', 'Value': f'{actual_recall:.1%}'},\n    {'Metric': 'Model — Precision', 'Value': f'{actual_precision:.1%}'},\n    {'Metric': 'Model — PR-AUC', 'Value': f'{pr_auc:.4f}'},\n    {'Metric': 'Threshold (Cost-Optimized)', 'Value': f'{OPTIMAL_T:.2f}'},\n    {'Metric': '', 'Value': ''},\n    {'Metric': 'Markdown Avoided (TP)', 'Value': f'{avoided_markdown:,.2f} SEK'},\n    {'Metric': 'Markdown Remaining (FN)', 'Value': f'{fn_markdown_loss:,.2f} SEK'},\n    {'Metric': 'Full-Price STR — Fad', 'Value': f'{current_fad_str:.1%} → {improved_fad_str:.1%}'},\n    {'Metric': 'Full-Price STR — Overall', 'Value': f'{current_overall_str:.1%} → {improved_overall_str:.1%}'},\n    {'Metric': '', 'Value': ''},\n    {'Metric': 'Year 1 ROI', 'Value': f'{roi_year1:+.1f}%'},\n    {'Metric': 'Year 2+ ROI', 'Value': f'{roi_year2:+.1f}%'},\n])\n\nprint('\\n' + summary.to_string(index=False))\n\nprint(f'\\n⚠️  Caveat: H&M price data ถูก encode/normalize ไม่ใช่ราคาจริง')\nprint(f'   ตัวเลข absolute (SEK) ใช้เปรียบเทียบสัดส่วนได้ แต่ไม่ใช่มูลค่าจริง')\nprint(f'   Pattern และ % improvement เชื่อถือได้')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T14:52:45.541072Z","iopub.execute_input":"2026-04-29T14:52:45.541377Z","iopub.status.idle":"2026-04-29T14:53:13.883771Z","shell.execute_reply.started":"2026-04-29T14:52:45.541352Z","shell.execute_reply":"2026-04-29T14:53:13.882945Z"}},"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-29T14:53:13.884677Z","iopub.execute_input":"2026-04-29T14:53:13.885010Z","iopub.status.idle":"2026-04-29T14:53:18.302604Z","shell.execute_reply.started":"2026-04-29T14:53:13.884987Z","shell.execute_reply":"2026-04-29T14:53:18.301890Z"}},"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'{avoided_markdown:,.0f} SEK avoided (test set)'\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\n# ═══════════════════════════════════════════════════════\n#  Full-Scale Projection\n#  (Test set → Full catalog)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  FULL-SCALE PROJECTION')\nprint('  (Scale จาก test set → full catalog)')\nprint('=' * 70)\n\n# Scale factor\nn_test = len(results_df)\nn_full = 80899\nscale_factor = n_full / n_test\n\nprojected_saving_sek = avoided_markdown * scale_factor\nprojected_saving_usd = projected_saving_sek * SEK_TO_USD\nannual_saving_usd = projected_saving_usd * 2  # 2 cycles per year\n\nprint(f'\\n   Test set:     {n_test:,} articles → {avoided_markdown:,.0f} SEK saved')\nprint(f'   Scale factor: {scale_factor:.2f}x')\nprint(f'   Full catalog: {n_full:,} articles → {projected_saving_sek:,.0f} SEK saved')\nprint(f'   Per year (2 cycles): ${annual_saving_usd:,.0f} USD')\n\n# Projected ROI\nproj_roi_year1 = (annual_saving_usd - total_year1_cost) / total_year1_cost * 100\nproj_payback = total_year1_cost / (annual_saving_usd / 12)\n\nprint(f'\\n   Year 1 ROI:   {proj_roi_year1:+.1f}%')\nprint(f'   Payback:      {proj_payback:.0f} months')\nprint(f'   Annual savings vs cloud cost: {annual_saving_usd/annual_cloud_cost:.0f}x return')\n\n# spark.stop()\n# print('Spark session stopped.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-29T15:32:27.361474Z","iopub.execute_input":"2026-04-29T15:32:27.362401Z","iopub.status.idle":"2026-04-29T15:32:27.386918Z","shell.execute_reply.started":"2026-04-29T15:32:27.362368Z","shell.execute_reply":"2026-04-29T15:32:27.386212Z"}},"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":{}}]}