{"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-27T20:12:50.748669Z","iopub.execute_input":"2026-04-27T20:12:50.749247Z","iopub.status.idle":"2026-04-27T20:12:50.753213Z","shell.execute_reply.started":"2026-04-27T20:12:50.749216Z","shell.execute_reply":"2026-04-27T20:12:50.752520Z"}},"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-27T20:12:52.116888Z","iopub.execute_input":"2026-04-27T20:12:52.117472Z","iopub.status.idle":"2026-04-27T20:13:55.373650Z","shell.execute_reply.started":"2026-04-27T20:12:52.117441Z","shell.execute_reply":"2026-04-27T20:13:55.372993Z"}},"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-27T20:14:58.307051Z","iopub.execute_input":"2026-04-27T20:14:58.308266Z","iopub.status.idle":"2026-04-27T20:15:09.594858Z","shell.execute_reply.started":"2026-04-27T20:14:58.308230Z","shell.execute_reply":"2026-04-27T20:15:09.594014Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 3: Load Data with PySpark ─────────────────────────\n# ✅ PySpark reads CSV lazily — no memory spike at load time\nBASE = '/kaggle/input/competitions/h-and-m-personalized-fashion-recommendations/'\n\ntxn = (\n    spark.read.csv(BASE + 'transactions_train.csv', header=True, inferSchema=True)\n    .withColumn('t_dat', F.to_date('t_dat', 'yyyy-MM-dd'))\n    .withColumn('article_id',  F.col('article_id').cast('string'))\n    .withColumn('customer_id', F.col('customer_id').cast('string'))\n)\n\nart = spark.read.csv(BASE + 'articles.csv',  header=True, inferSchema=True)\ncus = spark.read.csv(BASE + 'customers.csv', header=True, inferSchema=True)\n\n# Cache articles (small, used many times)\nart.cache()\n\nprint('Schema: transactions')\ntxn.printSchema()\nprint(f'\\nDate range: {txn.agg(F.min(\"t_dat\"), F.max(\"t_dat\")).collect()[0]}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:15:52.094122Z","iopub.execute_input":"2026-04-27T20:15:52.094937Z","iopub.status.idle":"2026-04-27T20:17:09.375388Z","shell.execute_reply.started":"2026-04-27T20:15:52.094904Z","shell.execute_reply":"2026-04-27T20:17:09.374556Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 4: Row Counts (Spark action — triggers execution) ──\n# NOTE: count() is the ONE action needed. Everything before this is lazy.\nn_txn = txn.count()\nn_art = art.count()\nn_cus = cus.count()\n\nprint(f'Transactions : {n_txn:>12,}')\nprint(f'Articles     : {n_art:>12,}')\nprint(f'Customers    : {n_cus:>12,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:17:44.026238Z","iopub.execute_input":"2026-04-27T20:17:44.026777Z","iopub.status.idle":"2026-04-27T20:17:52.639580Z","shell.execute_reply.started":"2026-04-27T20:17:44.026751Z","shell.execute_reply":"2026-04-27T20:17:52.638501Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 0 — EDA with PySpark\nทำ EDA บน data เต็ม ","metadata":{}},{"cell_type":"code","source":"# ── Cell 5: Missing Value Analysis (Spark) ──────────────────\n# PySpark: count nulls per column in one pass over the data\ndef spark_missing_report(df, name):\n    # Combine count(*) with null counts in ONE pass\n    agg_exprs = [F.count(F.lit(1)).alias('__total__')] + [\n        F.sum(F.col(c).isNull().cast('int')).alias(c) for c in df.columns\n    ]\n    result = df.agg(*agg_exprs).collect()[0].asDict()\n    total = result.pop('__total__')\n    \n    rows = [\n        {\n            'column': col, \n            'missing': cnt, \n            'pct': round(cnt/total*100, 2),\n            'result': '🔴 MNAR' if cnt/total > 0.30 \n                    else '🟡 MAR' if cnt/total > 0.05 \n                    else '🟢 MCAR'\n        }\n        for col, cnt in result.items() if cnt > 0\n    ]\n    \n    print(f'\\n=== Missing Values: {name} (total rows: {total:,}) ===')\n    if not rows:\n        print('No missing values'); return pd.DataFrame()\n    \n    report = pd.DataFrame(rows).sort_values('pct', ascending=False)\n    print(report.to_string(index=False))\n    return report\n\nspark_missing_report(art, 'articles')\nspark_missing_report(cus, 'customers')\n# transactions: only check key columns\nspark_missing_report(txn.select('t_dat','article_id','customer_id','price'), 'transactions (key cols)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:17:58.270616Z","iopub.execute_input":"2026-04-27T20:17:58.270889Z","iopub.status.idle":"2026-04-27T20:18:36.245432Z","shell.execute_reply.started":"2026-04-27T20:17:58.270867Z","shell.execute_reply":"2026-04-27T20:18:36.243821Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 6: Cardinality Check (Spark) ──────────────────────\ncat_cols = [\n    'product_type_name', 'product_group_name',\n    'graphical_appearance_name', 'colour_group_name',\n    'department_name', 'index_name', 'garment_group_name'\n]\ncard_exprs = [F.countDistinct(c).alias(c) for c in cat_cols]\ncardinality = art.agg(*card_exprs).collect()[0].asDict()\ncard_df = pd.DataFrame(cardinality.items(), columns=['feature','unique_values'])\ncard_df = card_df.sort_values('unique_values', ascending=False)\n\nprint('=== Cardinality (articles) ===')\nprint(card_df.to_string(index=False))\n\nfig, ax = plt.subplots(figsize=(9, 5))\ncard_df.plot(kind='barh', x='feature', y='unique_values', ax=ax, color='steelblue', legend=False)\nax.axvline(50, color='red', linestyle='--', label='One-Hot limit (50)')\nax.set_title('Cardinality — Articles Categorical Features', fontweight='bold')\nax.legend()\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:18:48.036465Z","iopub.execute_input":"2026-04-27T20:18:48.037124Z","iopub.status.idle":"2026-04-27T20:18:49.713137Z","shell.execute_reply.started":"2026-04-27T20:18:48.037087Z","shell.execute_reply":"2026-04-27T20:18:49.712482Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 7: Sales Distribution + Cold-Start Analysis ──\n!pip install powerlaw -q\nimport powerlaw\n\n# 1. Identify \"Dead-on-Arrival\" articles (CRITICAL — missed in original)\nall_articles = art.select('article_id').distinct()\nsold_articles = txn.select('article_id').distinct()\n\ndead_on_arrival = all_articles.subtract(sold_articles)\nn_dead = dead_on_arrival.count()\nn_total_art = all_articles.count()\n\nprint(f'\\n Dead-on-Arrival Articles: {n_dead:,} / {n_total_art:,} ({n_dead/n_total_art*100:.2f}%)')\nprint('   → These are MAX-RISK candidates. Investigate before Phase 1.')\n\n# Profile dead-on-arrival items by category (do they cluster in specific product types?)\ndead_profile = (\n    dead_on_arrival.join(art, 'article_id')\n    .groupBy('product_group_name')\n    .count()\n    .orderBy(F.desc('count'))\n)\nprint('\\nDead-on-Arrival profile by product_group:')\ndead_profile.show(10, truncate=False)\n\n# 2. Sales distribution — already cached\nsales_per_art = (\n    txn.groupBy('article_id')\n    .count()\n    .withColumnRenamed('count', 'total_sales')\n    .orderBy(F.desc('total_sales'))\n)\nsales_per_art.cache()\n\nsales_per_art.select('total_sales').describe().show()\nsales_pd = sales_per_art.toPandas()\n\n# 3. RIGOROUS Power Law fit (Clauset et al. method)\nfit = powerlaw.Fit(sales_pd['total_sales'].values, discrete=True, verbose=False)\nprint(f'\\n=== Rigorous Power Law Fit ===')\nprint(f'  Estimated alpha (α)  : {fit.alpha:.3f}')\nprint(f'  Estimated x_min      : {fit.xmin:.0f}')\nprint(f'  → Power Law applies for sales >= {fit.xmin:.0f}')\n\n# Compare against alternative distributions\nfor alt in ['lognormal', 'exponential', 'stretched_exponential']:\n    R, p = fit.distribution_compare('power_law', alt, normalized_ratio=True)\n    verdict = '✓ PL better' if R > 0 and p < 0.05 else '✗ PL not better'\n    print(f'  PL vs {alt:25s}: R={R:+.3f}, p={p:.4f}  → {verdict}')\n\n# 4. Quantile breakdown for stratified sampling (informs Phase 1)\nprint('\\n=== Sales Quantiles (for Stratified Sampling in Phase 1) ===')\nquantiles = sales_per_art.approxQuantile('total_sales', [0.5, 0.75, 0.90, 0.95, 0.99], 0.01)\nlabels = ['p50', 'p75', 'p90', 'p95', 'p99']\nfor l, q in zip(labels, quantiles):\n    print(f'  {l}: {q:>8.0f} sales')\n\n# 5. Visualization — log-log + powerlaw fit overlay\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Left: rank-frequency log-log\nsales_pd['rank'] = range(1, len(sales_pd) + 1)\naxes[0].loglog(sales_pd['rank'], sales_pd['total_sales'], 'o', alpha=0.3, ms=2)\naxes[0].set_xlabel('Rank (log)'); axes[0].set_ylabel('Sales (log)')\naxes[0].set_title(f'Rank-Frequency (slope ≈ -{fit.alpha:.2f})')\naxes[0].grid(True, alpha=0.3)\n\n# Right: CCDF with powerlaw fit\nfit.plot_ccdf(ax=axes[1], color='steelblue', linewidth=2, label='Empirical CCDF')\nfit.power_law.plot_ccdf(ax=axes[1], color='red', linestyle='--', label=f'PL fit (α={fit.alpha:.2f})')\nfit.lognormal.plot_ccdf(ax=axes[1], color='green', linestyle=':', label='Lognormal fit')\naxes[1].set_xlabel('Sales (log)'); axes[1].set_ylabel('P(X ≥ x)')\naxes[1].set_title('CCDF Comparison: Power Law vs Lognormal')\naxes[1].legend()\n\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:19:18.824405Z","iopub.execute_input":"2026-04-27T20:19:18.825063Z","iopub.status.idle":"2026-04-27T20:20:45.210589Z","shell.execute_reply.started":"2026-04-27T20:19:18.825033Z","shell.execute_reply":"2026-04-27T20:20:45.209856Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 1 — Data Preparation with PySpark","metadata":{}},{"cell_type":"code","source":"# ═══════════════════════════════════════════════════════════════════\n#  LOCKED-IN CONSTANTS FROM PHASE 0 EDA\n#  อ้างอิงจาก Phase 0.1 (Missing Values) และ 0.2 (Power Law)\n# ═══════════════════════════════════════════════════════════════════\n\n# From Phase 0.2 — Lognormal + Power Law tail analysis\nPOWER_LAW_XMIN    = 1821     # regime boundary (top 3% = viral tail)\nPOWER_LAW_ALPHA   = 3.00     # tail exponent (edge of finite variance)\n\n# From Phase 0.1 — Dead-on-Arrival handling\nN_DEAD_ON_ARRIVAL = 995      # articles with 0 transactions (exclude from labeling)\n\n# Phase 1 configuration (informed by methodology + EDA)\nMIN_ACTIVE_WEEKS  = 3        # ✅ lowered from 6 to keep true Fads (short-lived = Fad signal)\nCENSORING_WEEKS   = 26       # articles launched < 26 weeks before dataset_end → censored\nMAX_WEEKS_PIVOT   = 26       # pivot width for time-series analysis\n\n# Stockout detection parameters\nSTOCKOUT_BASELINE_WINDOW = 4   # rolling window size (weeks)\nSTOCKOUT_DROP_THRESHOLD  = 0.30  # sales < 30% of baseline\nSTOCKOUT_RECOVERY_THRESH = 0.60  # recovery > 60% of baseline within 2 weeks\nSTOCKOUT_MIN_BASELINE    = 2     # baseline must have meaningful volume\n\nprint('Phase 0 constants loaded')\nprint(f'   x_min = {POWER_LAW_XMIN}, α = {POWER_LAW_ALPHA}')\nprint(f'   Min active weeks = {MIN_ACTIVE_WEEKS} (Fad-preserving)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:21:16.039589Z","iopub.execute_input":"2026-04-27T20:21:16.040221Z","iopub.status.idle":"2026-04-27T20:21:16.046063Z","shell.execute_reply.started":"2026-04-27T20:21:16.040187Z","shell.execute_reply":"2026-04-27T20:21:16.045400Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 8: Weekly Aggregation with DENSE time series ──\n\n# Get dataset boundary for censoring decision\ndataset_start, dataset_end = txn.agg(\n    F.min('t_dat').alias('start'),\n    F.max('t_dat').alias('end')\n).collect()[0]\nprint(f'Dataset period: {dataset_start} → {dataset_end}')\n\n# ─── Step 1: Aggregate to SPARSE weekly (only weeks with sales) ───\nweekly_sparse = (\n    txn\n    .withColumn('week_start', F.date_trunc('week', F.col('t_dat')))  # Monday as week start\n    .groupBy('article_id', 'week_start')\n    .agg(\n        F.count(F.lit(1))                  .alias('weekly_sales'),\n        F.countDistinct('customer_id')     .alias('unique_buyers'),\n        F.countDistinct('sales_channel_id').alias('n_channels')\n    )\n)\n\n# ─── Step 2: Compute launch_week & last_seen_week per article ───\nlifecycle = (\n    weekly_sparse\n    .groupBy('article_id')\n    .agg(\n        F.min('week_start').alias('launch_week'),\n        F.max('week_start').alias('last_seen_week')\n    )\n    .withColumn(\n        'weeks_from_launch_to_dataset_end',\n        (F.datediff(F.lit(dataset_end), F.col('launch_week')) / 7).cast('int')\n    )\n    # ✅ Tag right-censored: launched too close to dataset end to observe full lifecycle\n    .withColumn(\n        'is_censored',\n        (F.col('weeks_from_launch_to_dataset_end') < CENSORING_WEEKS).cast('int')\n    )\n)\n\n# Quick stats\nn_total = lifecycle.count()\nn_censored = lifecycle.filter(F.col('is_censored') == 1).count()\nprint(f'Total articles with transactions: {n_total:,}')\nprint(f'  → Right-censored (launched < {CENSORING_WEEKS}w before end): {n_censored:,} ({n_censored/n_total*100:.1f}%)')\nprint(f'  → Usable for training: {n_total - n_censored:,}')\n\n# ─── Step 3: Build DENSE skeleton (every week from launch → last_seen) ───\n# ✅ Uses sequence() to explode all weeks — enables proper lag/lead operations\ndense_skeleton = (\n    lifecycle\n    .filter(F.col('is_censored') == 0)  # exclude censored from skeleton\n    .select(\n        'article_id',\n        'launch_week',\n        F.explode(\n            F.sequence(\n                F.col('launch_week'),\n                F.col('last_seen_week'),\n                F.expr('INTERVAL 7 DAYS')\n            )\n        ).alias('week_start')\n    )\n)\n\n# ─── Step 4: Join sparse data onto dense skeleton → zero-fill gaps ───\nweekly = (\n    dense_skeleton\n    .join(weekly_sparse, on=['article_id', 'week_start'], how='left')\n    .fillna(0, subset=['weekly_sales', 'unique_buyers', 'n_channels'])\n    .withColumn(\n        'weeks_since_launch',\n        (F.datediff('week_start', 'launch_week') / 7).cast('int')\n    )\n)\n\n# Cache the dense weekly DataFrame (used many times in Phase 1-3)\nweekly.cache()\nn_records = weekly.count()  # trigger cache\n\nprint(f'\\n✅ Dense weekly records: {n_records:,}')\nprint(f'   Expansion ratio: {n_records / weekly_sparse.count():.2f}x '\n      f'(sparse → dense; higher = more stockout/dormant weeks)')\n\n# Preview\nprint('\\n--- Sample (article 110065002) ---')\nweekly.filter(F.col('article_id') == '110065002') \\\n      .select('week_start', 'weeks_since_launch', 'weekly_sales', 'unique_buyers') \\\n      .orderBy('weeks_since_launch').show(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:21:20.267810Z","iopub.execute_input":"2026-04-27T20:21:20.268142Z","iopub.status.idle":"2026-04-27T20:25:04.150653Z","shell.execute_reply.started":"2026-04-27T20:21:20.268113Z","shell.execute_reply":"2026-04-27T20:25:04.149918Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 9: Multi-Tier Stockout Detection ──\n\n# ─── Updated Thresholds (relaxed) ───\nSTOCKOUT_BASELINE_WINDOW = 4         # unchanged\nSTOCKOUT_MIN_BASELINE    = 3.0       # slight bump — reject noise แต่ไม่รัด\nDROP_SEVERE   = 0.20   # drop ≤ 20% of baseline (ช่วย severe)\nDROP_MODERATE = 0.50   # drop ≤ 50% of baseline (ช่วย moderate)\nRECOVERY_STRONG  = 0.70   # recover ≥ 70% within 2 weeks\nRECOVERY_PARTIAL = 0.40   # recover ≥ 40% within 3 weeks\n\nw_art      = Window.partitionBy('article_id').orderBy('weeks_since_launch')\nw_baseline = w_art.rowsBetween(-STOCKOUT_BASELINE_WINDOW, -1)\n\nweekly = (\n    weekly\n    # 1. Rolling baseline (past 4 weeks)\n    .withColumn('baseline_sales', F.avg('weekly_sales').over(w_baseline))\n    \n    # 2. Future recovery window (extend to 3 weeks — some stockouts recover slower)\n    .withColumn('sales_lead1', F.lead('weekly_sales', 1).over(w_art))\n    .withColumn('sales_lead2', F.lead('weekly_sales', 2).over(w_art))\n    .withColumn('sales_lead3', F.lead('weekly_sales', 3).over(w_art))\n    \n    # 3. Ratios\n    .withColumn('drop_ratio',\n        F.col('weekly_sales') / (F.col('baseline_sales') + 1.0))\n    \n    # Best recovery within 3 weeks (max lead)\n    .withColumn('future_recovery',\n        F.greatest(F.col('sales_lead1'), F.col('sales_lead2'), F.col('sales_lead3')) /\n        (F.col('baseline_sales') + 1.0))\n    \n    # 4. Tiered stockout detection (3 levels of severity)\n    .withColumn('stockout_tier',\n        F.when(\n            # Tier 3 (SEVERE): ดรอปฮวบ 80%+ AND recover strongly\n            (F.col('drop_ratio')      < DROP_SEVERE) &\n            (F.col('future_recovery') > RECOVERY_STRONG) &\n            (F.col('baseline_sales')  >= STOCKOUT_MIN_BASELINE),\n            3\n        ).when(\n            # Tier 2 (MODERATE): ดรอป 50%+ AND recover partially\n            (F.col('drop_ratio')      < DROP_MODERATE) &\n            (F.col('future_recovery') > RECOVERY_PARTIAL) &\n            (F.col('baseline_sales')  >= STOCKOUT_MIN_BASELINE),\n            2\n        ).when(\n            # Tier 1 (MILD): ตกกะทันหัน > 50% หลังช่วงขายสูง (baseline ≥ 5)\n            #   → ไม่ต้องรอ recovery เพราะอาจเป็น end-of-life stockout\n            (F.col('drop_ratio')     < DROP_MODERATE) &\n            (F.col('baseline_sales') >= 5.0),\n            1\n        ).otherwise(0))\n    \n    # 5. Binary flag (สำหรับ downstream ที่ต้องการ 0/1 เดิม)\n    .withColumn('is_stockout',\n        (F.col('stockout_tier') > 0).cast('int'))\n    \n    # 6. Impute value ตาม tier\n    .withColumn('weekly_sales_clean',\n        F.when(F.col('stockout_tier') == 3, F.col('baseline_sales'))\n         .when(F.col('stockout_tier') == 2, F.col('baseline_sales') * 0.7)  # partial impute\n         .when(F.col('stockout_tier') == 1, F.col('baseline_sales') * 0.5)  # conservative\n         .otherwise(F.col('weekly_sales')))\n)\n\n# Re-cache\nweekly.unpersist()\nweekly.cache()\n_ = weekly.count()\n\n# ─── Diagnostic ───\ntotal_weeks  = weekly.count()\ntier_dist = (\n    weekly.groupBy('stockout_tier')\n    .agg(F.count('*').alias('count'))\n    .orderBy('stockout_tier')\n).toPandas()\ntier_dist['pct'] = tier_dist['count'] / total_weeks * 100\n\nprint('=== Stockout Tier Distribution ===')\nprint(tier_dist.to_string(index=False))\n\ntotal_stockout_rate = tier_dist.loc[tier_dist['stockout_tier'] > 0, 'pct'].sum()\nprint(f'\\nTotal stockout rate (all tiers): {total_stockout_rate:.2f}%')\nprint(f'   Expected: 5-15% in fast fashion')\n\nif   total_stockout_rate < 4:  print('   ⚠️  ยังต่ำ — พิจารณา relax อีก')\nelif total_stockout_rate > 20: print('   ⚠️  สูงเกิน — tighten thresholds')\nelse:                          print('   ✅ Within expected range')\n\n# Sample ของแต่ละ tier\nprint('\\n--- Samples by Tier ---')\nfor tier in [3, 2, 1]:\n    print(f'\\nTier {tier}:')\n    weekly.filter(F.col('stockout_tier') == tier) \\\n          .select('article_id', 'weeks_since_launch', 'baseline_sales',\n                  'weekly_sales', 'sales_lead1', 'sales_lead2', \n                  'drop_ratio', 'future_recovery', 'weekly_sales_clean') \\\n          .show(3, truncate=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:25:30.358091Z","iopub.execute_input":"2026-04-27T20:25:30.358749Z","iopub.status.idle":"2026-04-27T20:25:40.623055Z","shell.execute_reply.started":"2026-04-27T20:25:30.358723Z","shell.execute_reply":"2026-04-27T20:25:40.622163Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 10: Filter + Normalize + Regime-Aware Stratification ──\n# ─── Step 1: Compute active weeks per article (non-zero sales) ───\nactive_weeks = (\n    weekly\n    .filter(F.col('weekly_sales_clean') > 0)\n    .groupBy('article_id')\n    .agg(F.count('*').alias('active_weeks'))\n)\n\n# ─── Step 2: Filter: >= 3 active weeks (NOT lifespan!) ───\n# ✅ ทำไม >= 3: ต้องมี data พอสำหรับ early signal features (weeks 1-3)\n#    ถ้า < 3 = ขายแค่ 1-2 ครั้งแล้วหาย = เป็น \"dud\" ไม่ใช่ \"Fad\"\nvalid_arts = active_weeks.filter(F.col('active_weeks') >= MIN_ACTIVE_WEEKS)\n\nn_before = weekly.select('article_id').distinct().count()\nn_after  = valid_arts.count()\n\nweekly = weekly.join(valid_arts.select('article_id'), on='article_id', how='inner')\nprint(f'Articles before filter: {n_before:,}')\nprint(f'Articles after filter (>= {MIN_ACTIVE_WEEKS} active weeks): {n_after:,}')\nprint(f'   Excluded: {n_before - n_after:,} articles with <{MIN_ACTIVE_WEEKS} active weeks (duds)')\n\n# ─── Step 3: Min-Max normalize per article (with edge-case handling) ───\nw_norm = Window.partitionBy('article_id')\nweekly = (\n    weekly\n    .withColumn('art_min', F.min('weekly_sales_clean').over(w_norm))\n    .withColumn('art_max', F.max('weekly_sales_clean').over(w_norm))\n    .withColumn('norm_sales',\n        # ✅ ป้องกัน division by zero เมื่อ article มี constant sales\n        F.when(F.col('art_max') == F.col('art_min'), F.lit(0.0))\n         .otherwise(\n             (F.col('weekly_sales_clean') - F.col('art_min')) /\n             (F.col('art_max') - F.col('art_min'))\n         ))\n)\n\n# ─── Step 4: Regime-aware stratification (NEW — from Phase 0.2) ───\n# ✅ แทนที่ uniform ntile(4) ด้วย \"body + tail\" regime\n#    Body (97%): sales < 1821 → log-quantile stratify\n#    Tail (3%):  sales >= 1821 → ทุก article อยู่ใน tail (rare & critical)\ntotal_sales_df = (\n    weekly.groupBy('article_id')\n    .agg(F.sum('weekly_sales_clean').alias('total_sales'))\n    # Regime tag (from Phase 0.2 Lognormal + Power Law tail analysis)\n    .withColumn('regime',\n        F.when(F.col('total_sales') >= POWER_LAW_XMIN, 'tail')\n         .otherwise('body'))\n)\n\n# Log-quantile stratification WITHIN each regime\n# (body: 4 strata by log-sales; tail: 2 strata because data is sparse)\nw_body = Window.partitionBy('regime').orderBy(F.log1p('total_sales'))\ntotal_sales_df = (\n    total_sales_df\n    .withColumn(\n        'strat_bucket',\n        F.when(F.col('regime') == 'tail', F.ntile(2).over(w_body))\n         .otherwise(F.ntile(4).over(w_body)))\n    .withColumn(\n        'strat_key',\n        F.concat(F.col('regime'), F.lit('_Q'), F.col('strat_bucket')))\n)\n\n# Regime distribution report\nprint('\\n=== Regime Distribution (from Power Law EDA) ===')\nregime_dist = total_sales_df.groupBy('regime').agg(\n    F.count('*').alias('n_articles'),\n    F.min('total_sales').alias('min_sales'),\n    F.max('total_sales').alias('max_sales'),\n    F.mean('total_sales').alias('mean_sales')\n).toPandas()\nprint(regime_dist.to_string(index=False))\n\nprint('\\n=== Stratum Key Distribution ===')\ntotal_sales_df.groupBy('strat_key').count().orderBy('strat_key').show()\n\n# Join stratum back to weekly\nweekly = weekly.join(\n    total_sales_df.select('article_id', 'regime', 'strat_key', 'total_sales'),\n    on='article_id', how='left'\n)\n\n# Re-cache final weekly DataFrame\nweekly.unpersist()\nweekly.cache()\nn_valid = weekly.select('article_id').distinct().count()\nprint(f'\\n✅ Final valid articles: {n_valid:,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:25:45.225036Z","iopub.execute_input":"2026-04-27T20:25:45.225346Z","iopub.status.idle":"2026-04-27T20:26:05.162251Z","shell.execute_reply.started":"2026-04-27T20:25:45.225320Z","shell.execute_reply":"2026-04-27T20:26:05.160797Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 11: Sales Pivot with Validation ──\n\npivot_spark = (\n    weekly\n    .filter(\n        (F.col('weeks_since_launch') >= 0) & \n        (F.col('weeks_since_launch') < MAX_WEEKS_PIVOT)\n    )\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(MAX_WEEKS_PIVOT)))\n    .agg(F.first('norm_sales'))\n    .fillna(0)\n)\n\n# Rename columns: '0' → 'w0', etc.\nnew_cols = ['article_id'] + [f'w{i}' for i in range(MAX_WEEKS_PIVOT)]\npivot_spark = pivot_spark.toDF(*new_cols)\n\n# Pre-collect sanity check (validate before bringing to Pandas)\nprint('=== Pre-Collect Sanity Check ===')\npivot_spark.select(\n    F.count('*').alias('n_articles'),\n    F.mean('w0').alias('mean_w0'),\n    F.mean('w3').alias('mean_w3'),\n    F.mean('w12').alias('mean_w12'),\n    F.mean('w25').alias('mean_w25')\n).show()\nprint('↑ Expect: mean_w0 > mean_w3 > mean_w12 > mean_w25 '\n      '(peak at launch, decay over time)')\n\n# Collect to Pandas (89K × 26 = ~18MB, safe)\npivot = pivot_spark.toPandas().set_index('article_id')\nprint(f'\\n✅ Sales pivot shape: {pivot.shape}  (articles × weeks)')\n\n# ─── Diagnostic 1: Death-week distribution (informs Phase 2 labeling) ───\n# \"Death week\" = last week with non-zero sales → short death_week = Fad-like\nweek_cols = [f'w{i}' for i in range(MAX_WEEKS_PIVOT)]\ndeath_week = (pivot[week_cols] > 0).iloc[:, ::-1].cummax(axis=1).iloc[:, ::-1].sum(axis=1)\n\nprint('\\n=== Death Week Distribution ===')\nprint(f'(last week with sales > 0; short death_week = Fad candidate)')\nprint(death_week.describe().to_string())\nprint(f'\\n  Articles dying by week 6 (short Fad):  '\n      f'{(death_week <= 6).sum():,} ({(death_week <= 6).mean()*100:.1f}%)')\nprint(f'  Articles dying by week 12 (medium):    '\n      f'{((death_week > 6) & (death_week <= 12)).sum():,} '\n      f'({((death_week > 6) & (death_week <= 12)).mean()*100:.1f}%)')\nprint(f'  Articles still alive at week 25:       '\n      f'{(death_week > 25).sum():,} ({(death_week > 25).mean()*100:.1f}%)')\n\n# ─── Diagnostic 2: Visualize curve shapes by regime ───\nregime_df = total_sales_df.select('article_id', 'regime').toPandas().set_index('article_id')\npivot_with_regime = pivot.join(regime_df, how='left')\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4))\n\n# Sample 30 articles per regime\nfor ax, regime, color in zip(axes, ['body', 'tail'], ['steelblue', 'crimson']):\n    sample = pivot_with_regime[pivot_with_regime['regime'] == regime].sample(\n        min(30, (pivot_with_regime['regime'] == regime).sum()), random_state=42\n    )\n    for idx in sample.index:\n        ax.plot(range(MAX_WEEKS_PIVOT), sample.loc[idx, week_cols].values, \n                alpha=0.3, color=color, lw=1)\n    ax.set_title(f'{regime.upper()} regime — 30 random articles', fontweight='bold')\n    ax.set_xlabel('Weeks since launch'); ax.set_ylabel('Normalized sales')\n    ax.grid(alpha=0.3)\n\nplt.tight_layout(); plt.show()\n\nprint('\\n✅ Phase 1 complete — ready for Phase 2 (Hybrid Labeling)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:26:12.124307Z","iopub.execute_input":"2026-04-27T20:26:12.124609Z","iopub.status.idle":"2026-04-27T20:26:20.300474Z","shell.execute_reply.started":"2026-04-27T20:26:12.124583Z","shell.execute_reply":"2026-04-27T20:26:20.299699Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 2 — Hybrid Labeling (PELT + k-Mean)","metadata":{}},{"cell_type":"code","source":"# ── Cell 12: PELT with Phase 1 Insights ──\n\nFAD_MAX_PEAK_WEEK     = 6\nFAD_MIN_DROP_RATIO    = 0.70   # tighten 0.60 → 0.70\nFAD_MIN_PEAK_RATIO    = 4.0    # tighten 2.5 → 4.0 (biggest lever)\nFAD_LOW_TAIL_RATIO    = 0.20   # tighten 0.30 → 0.20\nPELT_PENALTY          = 2\n\npelt_results = []\nfor art_id, row in pivot.iterrows():\n    series = row.values.astype(float)\n    \n    nonzero = np.nonzero(series)[0]\n    if len(nonzero) == 0:\n        pelt_results.append({\n            'article_id': art_id, 'pelt_fad': 0, 'n_breakpoints': 0,\n            'peak_week': 0, 'death_week': 0,\n            'peak_ratio': 0.0, 'drop_ratio': 0.0, 'tail_ratio': 0.0,\n        })\n        continue\n    \n    death_week = int(nonzero[-1])\n    series_trimmed = series[:death_week + 1]\n    \n    if len(series_trimmed) < 5:\n        pelt_results.append({\n            'article_id': art_id, 'pelt_fad': 0, 'n_breakpoints': 0,\n            'peak_week': 0, 'death_week': death_week,\n            'peak_ratio': 0.0, 'drop_ratio': 0.0, 'tail_ratio': 0.0,\n        })\n        continue\n    \n    try:\n        algo = rpt.Pelt(model='rbf', min_size=2).fit(series_trimmed)\n        bkps = algo.predict(pen=PELT_PENALTY)\n        n_bkps = len(bkps) - 1\n    except Exception:\n        n_bkps = 0\n    \n    peak_week = int(np.argmax(series_trimmed))\n    peak_value = series_trimmed[peak_week]\n    \n    if peak_week < len(series_trimmed) - 1:\n        post_peak_mean = series_trimmed[peak_week + 1:].mean() + 1e-8\n        peak_ratio = peak_value / post_peak_mean\n        tail_ratio = post_peak_mean / (peak_value + 1e-8)\n    else:\n        peak_ratio = 0.0\n        tail_ratio = 1.0\n    \n    drop_ratio = (peak_value - series_trimmed[-1]) / (peak_value + 1e-8)\n    \n    is_fad = (\n        (n_bkps >= 1) and\n        (peak_week <= FAD_MAX_PEAK_WEEK) and\n        (peak_ratio > FAD_MIN_PEAK_RATIO) and\n        (tail_ratio < FAD_LOW_TAIL_RATIO) and\n        (drop_ratio > FAD_MIN_DROP_RATIO)\n    )\n    \n    pelt_results.append({\n        'article_id': art_id,\n        'pelt_fad': int(is_fad),\n        'n_breakpoints': n_bkps,\n        'peak_week': peak_week,\n        'death_week': death_week,\n        'peak_ratio': round(peak_ratio, 2),\n        'drop_ratio': round(drop_ratio, 3),\n        'tail_ratio': round(tail_ratio, 3),\n    })\n\npelt_df = pd.DataFrame(pelt_results)\n\n# ─── Diagnostics ───\nprint('✅ PELT done.')\nprint(f'\\n{pelt_df[\"pelt_fad\"].value_counts()}')\nprint(f'\\nFad rate (PELT): {pelt_df[\"pelt_fad\"].mean():.2%}')\n\nprint('\\n=== Criteria Pass Rates (individual) ===')\nprint(f'  Has breakpoint (n≥1)   : {(pelt_df[\"n_breakpoints\"] >= 1).mean():.2%}')\nprint(f'  Peak week ≤ {FAD_MAX_PEAK_WEEK}           : {(pelt_df[\"peak_week\"] <= FAD_MAX_PEAK_WEEK).mean():.2%}')\nprint(f'  Peak ratio > {FAD_MIN_PEAK_RATIO}          : {(pelt_df[\"peak_ratio\"] > FAD_MIN_PEAK_RATIO).mean():.2%}')\nprint(f'  Tail ratio < {FAD_LOW_TAIL_RATIO}          : {(pelt_df[\"tail_ratio\"] < FAD_LOW_TAIL_RATIO).mean():.2%}')\nprint(f'  Drop ratio > {FAD_MIN_DROP_RATIO}          : {(pelt_df[\"drop_ratio\"] > FAD_MIN_DROP_RATIO).mean():.2%}')\n\npelt_fads = pelt_df[pelt_df['pelt_fad'] == 1]\nif len(pelt_fads) > 0:\n    print('\\n=== PELT-Fad Candidates: Temporal Profile ===')\n    print(f'Peak week   — mean: {pelt_fads[\"peak_week\"].mean():.1f}, '\n          f'median: {pelt_fads[\"peak_week\"].median():.0f}')\n    print(f'Peak ratio  — mean: {pelt_fads[\"peak_ratio\"].mean():.2f}')\n    print(f'Drop ratio  — mean: {pelt_fads[\"drop_ratio\"].mean():.2f}')\n    print(f'Tail ratio  — mean: {pelt_fads[\"tail_ratio\"].mean():.3f}')\n\n# Sanity check\nprint('\\n=== Sanity Check: Article 110065002 ===')\ncheck = pelt_df[pelt_df['article_id'] == '110065002']\nif len(check) > 0:\n    print(check.to_string(index=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:26:45.366873Z","iopub.execute_input":"2026-04-27T20:26:45.367513Z","iopub.status.idle":"2026-04-27T20:27:23.523690Z","shell.execute_reply.started":"2026-04-27T20:26:45.367485Z","shell.execute_reply":"2026-04-27T20:27:23.522814Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 13 Distributed Clustering ──\n#\n#  ใช้ Spark MLlib KMeans แทน tslearn k-Shape\n#  ✅ Distributed — ใช้ทุก node ใน cluster\n#  ✅ เร็วกว่า k-Shape 10-50x\n#  ⚠️ ใช้ Euclidean distance แทน SBD — อาจจับ shape ได้ด้อยกว่า k-Shape เล็กน้อย\n#     แต่ข้อมูลถูก normalize แล้ว (Phase 1 Min-Max) จึงลด gap ตรงนี้ได้มาก\n# ═══════════════════════════════════════════════════════════════════\n\nfrom pyspark.ml.feature import VectorAssembler\nfrom pyspark.ml.clustering import KMeans, BisectingKMeans\nfrom pyspark.ml.evaluation import ClusteringEvaluator\nimport time\n\n# ─── Config ───\nK_SEARCH_RANGE    = range(3, 9)\nFAD_MAX_PEAK_WEEK = 6\n\n# ═══════════════════════════════════════════════════════\n#  Step 1: Prepare data in Spark (distributed)\n# ═══════════════════════════════════════════════════════\n\nprint('📊 Preparing data in Spark...')\n\n# pivot อยู่ใน Pandas → แปลงกลับเป็น Spark\nweek_cols = [f'w{i}' for i in range(26)]\npivot_reset = pivot.reset_index()\npivot_spark = spark.createDataFrame(pivot_reset)\n\n# VectorAssembler: รวม w0-w25 เป็น feature vector เดียว\nassembler = VectorAssembler(inputCols=week_cols, outputCol='features')\ndata_spark = assembler.transform(pivot_spark).select('article_id', 'features')\ndata_spark.cache()\n\nprint(f'   Rows: {data_spark.count():,}')\nprint(f'   Feature vector length: 26')\n\n# ═══════════════════════════════════════════════════════\n#  Step 2: Elbow Method — หา optimal k (Spark — distributed)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n🔍 Searching optimal k over range {list(K_SEARCH_RANGE)}...')\nprint('   (ทุก k รันแบบ distributed บน Spark)')\n\nelbow_results = []\nfor k in K_SEARCH_RANGE:\n    t0 = time.time()\n    \n    kmeans = KMeans(\n        k=k,\n        featuresCol='features',\n        predictionCol='cluster',\n        seed=42,\n        maxIter=100,\n        initMode='k-means||',  # scalable k-means++ (distributed)\n    )\n    model = kmeans.fit(data_spark)\n    \n    # Compute inertia (WSSSE — Within Set Sum of Squared Errors)\n    inertia = model.summary.trainingCost\n    elapsed = time.time() - t0\n    \n    elbow_results.append({'k': k, 'inertia': inertia, 'time': elapsed})\n    print(f'   k={k}  inertia={inertia:.2f}  time={elapsed:.1f}s')\n\nelbow_df = pd.DataFrame(elbow_results)\n\n# ─── Kneedle-based Elbow Detection ───\nks_vals = elbow_df['k'].values.astype(float)\ninertia_vals = elbow_df['inertia'].values\n\nx_norm = (ks_vals - ks_vals.min()) / (ks_vals.max() - ks_vals.min())\ny_norm = (inertia_vals - inertia_vals.min()) / (inertia_vals.max() - inertia_vals.min())\n\nline_vec = np.array([x_norm[-1] - x_norm[0], y_norm[-1] - y_norm[0]])\nline_vec_norm = line_vec / np.linalg.norm(line_vec)\n\ndistances = []\nfor i in range(len(ks_vals)):\n    point_vec = np.array([x_norm[i] - x_norm[0], y_norm[i] - y_norm[0]])\n    proj = np.dot(point_vec, line_vec_norm) * line_vec_norm\n    perp = point_vec - proj\n    distances.append(np.linalg.norm(perp))\n\nBEST_K = int(ks_vals[np.argmax(distances)])\nprint(f'\\n📍 Elbow at k = {BEST_K}')\n\n# Elbow Plot\nfig, axes = plt.subplots(1, 2, figsize=(13, 4))\naxes[0].plot(ks_vals, inertia_vals, 'o-', lw=2, ms=8, color='steelblue')\naxes[0].axvline(BEST_K, linestyle='--', color='crimson', lw=2, label=f'Elbow k={BEST_K}')\naxes[0].set_xlabel('k'); axes[0].set_ylabel('Inertia (WSSSE)')\naxes[0].set_title('Elbow Method — Spark KMeans'); axes[0].legend(); axes[0].grid(alpha=0.3)\n\naxes[1].bar(ks_vals, distances, color='coral', edgecolor='darkred')\naxes[1].axvline(BEST_K, linestyle='--', color='crimson', lw=2)\naxes[1].set_xlabel('k'); axes[1].set_ylabel('Perpendicular distance')\naxes[1].set_title('Kneedle Distance'); axes[1].grid(alpha=0.3)\nplt.suptitle(f'Optimal k Search — Elbow at k={BEST_K}', fontweight='bold')\nplt.tight_layout(); plt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 3: Fit final KMeans with best k (Spark — distributed)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n🎯 Fitting final KMeans with k={BEST_K} (distributed)...')\nt0 = time.time()\n\nfinal_kmeans = KMeans(\n    k=BEST_K,\n    featuresCol='features',\n    predictionCol='cluster',\n    seed=42,\n    maxIter=300,\n    initMode='k-means||',\n)\nfinal_model = final_kmeans.fit(data_spark)\nclustered = final_model.transform(data_spark)\n\nprint(f'   Elapsed: {time.time() - t0:.1f}s')\n\n# Silhouette Score (distributed evaluation)\nevaluator = ClusteringEvaluator(\n    featuresCol='features', predictionCol='cluster',\n    metricName='silhouette'\n)\nsilhouette = evaluator.evaluate(clustered)\nprint(f'   Silhouette Score: {silhouette:.4f}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 4: Extract labels + Compute centroids\n# ═══════════════════════════════════════════════════════\n\n# Labels: Spark → Pandas\nkshape_labels_df = clustered.select('article_id', 'cluster').toPandas()\nkshape_labels_df['article_id'] = kshape_labels_df['article_id'].astype(str)\n\n# Compute centroids manually (mean of members per cluster)\n# ใช้ Spark aggregation — distributed\ncentroid_spark = (\n    pivot_spark\n    .join(\n        clustered.select(\n            F.col('article_id').cast('string').alias('article_id'), \n            'cluster'\n        ),\n        on='article_id'\n    )\n)\n\ncentroids = {}\nfor i in range(BEST_K):\n    cluster_data = centroid_spark.filter(F.col('cluster') == i)\n    centroid_values = cluster_data.select(\n        *[F.avg(c).alias(c) for c in week_cols]\n    ).collect()[0]\n    centroids[i] = np.array([centroid_values[c] for c in week_cols])\n\nkshape_labels = kshape_labels_df.set_index('article_id')['cluster'].reindex(pivot.index).values\n\nprint(f'   Labels shape: {kshape_labels.shape}')\nprint(f'   Centroids computed: {len(centroids)}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 5: Auto-Detect Fad Clusters (same logic as before)\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n📊 Computing cluster metrics (k={BEST_K})...')\ncluster_metrics = []\nfor i in range(BEST_K):\n    centroid = centroids[i]\n    n_members = int((kshape_labels == i).sum())\n    \n    peak_w = int(np.argmax(centroid))\n    peak_v = centroid[peak_w]\n    end_v  = centroid[-1]\n    min_v  = centroid.min()\n    \n    prominence = peak_v - min_v\n    width_at_half = int((centroid > peak_v * 0.5).sum())\n    sharpness = prominence / max(1, width_at_half)\n    drop = (peak_v - end_v) / (abs(peak_v) + 1e-8)\n    is_early_peak = peak_w <= FAD_MAX_PEAK_WEEK\n    \n    cluster_metrics.append({\n        'cluster': i, 'n': n_members,\n        'peak_week': peak_w, 'peak_value': round(peak_v, 4),\n        'end_value': round(end_v, 4),\n        'prominence': round(prominence, 4),\n        'width_at_half': width_at_half,\n        'sharpness': round(sharpness, 4),\n        'drop': round(drop, 2),\n        'is_early_peak': is_early_peak,\n    })\n\ncm_df = pd.DataFrame(cluster_metrics).sort_values('sharpness', ascending=False)\nprint('\\n=== Cluster Metrics (sorted by sharpness) ===')\nprint(cm_df.to_string(index=False))\n\n# Auto-detect\nFAD_CRITERIA = {'prominence': 0.05, 'sharpness': 0.01, 'drop': 0.3}\n# ⚠️ Threshold ต่ำกว่า k-Shape version เพราะ KMeans centroids\n#    เป็น mean (ไม่ใช่ z-normalized) ค่าจะเล็กกว่ามาก\n\nfad_mask = (\n    (cm_df['prominence'] > FAD_CRITERIA['prominence']) &\n    (cm_df['sharpness']  > FAD_CRITERIA['sharpness']) &\n    (cm_df['drop']       > FAD_CRITERIA['drop']) &\n    (cm_df['is_early_peak'])\n)\nauto_fad_clusters = cm_df[fad_mask]['cluster'].tolist()\n\nprint(f'\\n🎯 Auto-detected FAD_CLUSTERS: {auto_fad_clusters}')\nprint(f'   Criteria: prominence > {FAD_CRITERIA[\"prominence\"]}, '\n      f'sharpness > {FAD_CRITERIA[\"sharpness\"]}, '\n      f'drop > {FAD_CRITERIA[\"drop\"]}, peak_week ≤ {FAD_MAX_PEAK_WEEK}')\n\n# Rejected clusters\nrejected = cm_df[~fad_mask]\nif len(rejected) > 0:\n    print(f'\\n   ❌ Rejected:')\n    for _, row in rejected.iterrows():\n        reasons = []\n        if row['prominence'] <= FAD_CRITERIA['prominence']: reasons.append(f'low prominence')\n        if row['sharpness'] <= FAD_CRITERIA['sharpness']: reasons.append(f'low sharpness')\n        if row['drop'] <= FAD_CRITERIA['drop']: reasons.append(f'low drop ({row[\"drop\"]})')\n        if not row['is_early_peak']: reasons.append(f'late peak (w{row[\"peak_week\"]})')\n        print(f'      Cluster {int(row[\"cluster\"])}: {\", \".join(reasons)}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 6: Visualize Centroids\n# ═══════════════════════════════════════════════════════\n\nfig, axes = plt.subplots(1, BEST_K, figsize=(max(13, 3.5 * BEST_K), 4.5))\nif BEST_K == 1: axes = [axes]\n\nfor i, ax in enumerate(axes):\n    centroid = centroids[i]\n    n = int((kshape_labels == i).sum())\n    metrics = cm_df[cm_df['cluster'] == i].iloc[0]\n    \n    is_fad = i in auto_fad_clusters\n    color = 'crimson' if is_fad else 'steelblue'\n    bg = '#ffebeb' if is_fad else 'white'\n    ax.set_facecolor(bg)\n    \n    ax.plot(centroid, color=color, lw=2.5)\n    ax.axvline(FAD_MAX_PEAK_WEEK, color='green', linestyle='--', alpha=0.3, lw=1)\n    ax.axvline(metrics['peak_week'], color='orange', linestyle=':', alpha=0.6, lw=1)\n    \n    title = f'Cluster {i} (n={n:,})\\npeak@w{metrics[\"peak_week\"]}, drop={metrics[\"drop\"]:.0%}'\n    title += '\\n🎯 FAD' if is_fad else '\\n(not Fad)'\n    ax.set_title(title, fontsize=10, fontweight='bold')\n    ax.set_xlabel('Weeks since launch'); ax.grid(alpha=0.3)\n\nplt.suptitle(f'Spark KMeans Centroids (k={BEST_K}) — Fad = peak ≤ week {FAD_MAX_PEAK_WEEK}',\n             fontweight='bold', fontsize=12)\nplt.tight_layout(); plt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 7: Save for Cell 14\n# ═══════════════════════════════════════════════════════\n\nkshape_df = pd.DataFrame({\n    'article_id': pivot.index,\n    'kshape_cluster': kshape_labels,\n})\nFAD_CLUSTERS = auto_fad_clusters\n\nn_fad = sum(cm_df[cm_df['cluster'].isin(FAD_CLUSTERS)]['n'])\nprint(f'\\n{\"=\"*60}')\nprint(f'✅ Cell 13 COMPLETE (Spark KMeans — Distributed)')\nprint(f'   Optimal k       : {BEST_K}')\nprint(f'   Silhouette Score : {silhouette:.4f}')\nprint(f'   Fad clusters     : {FAD_CLUSTERS}')\nprint(f'   Fad members      : {n_fad:,} ({n_fad/len(pivot)*100:.1f}%)')\nprint(f'   Non-Fad members  : {len(pivot)-n_fad:,}')\nprint(f'{\"=\"*60}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:28:08.971527Z","iopub.execute_input":"2026-04-27T20:28:08.972370Z","iopub.status.idle":"2026-04-27T20:30:07.797914Z","shell.execute_reply.started":"2026-04-27T20:28:08.972329Z","shell.execute_reply":"2026-04-27T20:30:07.797331Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════\n#  Cluster Validation & Profiling\n#  ตอบคำถามอาจารย์ 2 ข้อ:\n#    1. Cluster ที่เลือกเป็น Fad ถูกต้องไหม? → Statistical Tests\n#    2. ข้างใน cluster เป็นสินค้าอะไร? → Profiling\n# ══════════════════════════════════════════════════════════════════\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\nimport numpy as np\nimport pandas as pd\n\n# ─── ตรวจสอบตัวแปรที่ต้องใช้ ───\nprint('=== Pre-check ===')\nprint(f'  pivot shape:      {pivot.shape}')\nprint(f'  kshape_labels:    {len(kshape_labels)} labels, {len(set(kshape_labels))} clusters')\nprint(f'  BEST_K:           {BEST_K}')\nprint(f'  FAD_CLUSTERS:     {FAD_CLUSTERS}')\nprint(f'  centroids:        {len(centroids)} clusters')\n\n\n# ══════════════════════════════════════════════════════════════\n#  PART 1: VALIDATION — Cluster ที่เลือกถูกต้องไหม?\n# ══════════════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 1: CLUSTER VALIDATION')\nprint('  พิสูจน์ว่า Fad clusters ที่เลือกมาถูกต้อง')\nprint('=' * 70)\n\n# ─────────────────────────────────────────────────────────\n#  1.1 Visual Validation: Centroid + สมาชิกจริงทับกัน\n#      → ถ้าสมาชิกมีรูปร่างคล้าย centroid = cluster แน่น\n#      → ถ้ากระจัดกระจาย = cluster ไม่ดี\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.1 Visual: Centroid vs Actual Members ---')\n\nfig, axes = plt.subplots(2, BEST_K, figsize=(4.5 * BEST_K, 8))\n\nfor i in range(BEST_K):\n    member_idx = np.where(kshape_labels == i)[0]\n    n_members = len(member_idx)\n    \n    # สุ่มสมาชิก 50 ตัวมา plot ทับ centroid\n    sample_size = min(50, n_members)\n    sample_idx = np.random.choice(member_idx, sample_size, replace=False)\n    \n    is_fad = i in FAD_CLUSTERS\n    color = 'crimson' if is_fad else 'steelblue'\n    \n    # ── Row 1: Centroid + สมาชิกจริง overlay ──\n    ax = axes[0, i]\n    for idx in sample_idx:\n        ax.plot(pivot.iloc[idx].values, alpha=0.15, color=color, lw=0.5)\n    \n    centroid = centroids[i]\n    ax.plot(centroid, color='black', lw=3, label='Centroid')\n    ax.axvline(6, color='green', linestyle='--', alpha=0.4, lw=1, label='Fad boundary (w6)')\n    ax.set_title(f'Cluster {i} (n={n_members:,})\\n{\"🎯 FAD\" if is_fad else \"Non-Fad\"}',\n                 fontweight='bold', fontsize=11)\n    ax.set_xlabel('Weeks')\n    ax.set_ylabel('Normalized Sales')\n    ax.legend(fontsize=7, loc='upper right')\n    ax.grid(alpha=0.3)\n    \n    # ── Row 2: Peak Week Distribution ของสมาชิกจริง ──\n    # ถ้าเป็น Fad cluster → peak weeks ส่วนใหญ่ต้องอยู่ก่อนเส้น w6\n    ax2 = axes[1, i]\n    member_data = pivot.iloc[member_idx].values\n    peak_weeks = member_data.argmax(axis=1)\n    \n    ax2.hist(peak_weeks, bins=26, color=color, alpha=0.7, edgecolor='white')\n    ax2.axvline(6, color='red', linestyle='--', lw=2, label='Fad boundary (w6)')\n    \n    pct_early = (peak_weeks <= 6).mean() * 100\n    ax2.set_title(f'Peak Week Distribution\\nmedian={np.median(peak_weeks):.0f}, '\n                  f'{pct_early:.0f}% ≤ w6', fontsize=10)\n    ax2.set_xlabel('Peak Week')\n    ax2.set_ylabel('Count')\n    ax2.legend(fontsize=8)\n    ax2.grid(alpha=0.3)\n\nplt.suptitle('Cluster Validation: Centroid vs Actual Members\\n'\n             '(Row 1: เส้นบาง = สมาชิกจริง, เส้นหนา = centroid  |  '\n             'Row 2: กระจายตัวของ peak week)',\n             fontweight='bold', fontsize=12)\nplt.tight_layout()\nplt.show()\n\n# ─────────────────────────────────────────────────────────\n#  1.2 Quantitative: Cluster Quality Metrics\n#      → วัดว่าสมาชิกใน cluster คล้าย centroid แค่ไหน\n#      → Fad clusters ต้องมี early peak + high drop\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.2 Cluster Quality Metrics ---\\n')\n\n# Helper: compute drop ratios for an array of series\ndef compute_drops(data):\n    drops = []\n    for row in data:\n        pw = row.argmax()\n        peak_val = row[pw]\n        drop = (peak_val - row[-1]) / (peak_val + 1e-8) if peak_val > 0 else 0\n        drops.append(drop)\n    return np.array(drops)\n\n# Helper: compute Gini coefficients\ndef compute_gini_array(data):\n    ginis = []\n    for row in data:\n        arr = np.sort(np.abs(row))\n        n = len(arr)\n        if arr.sum() == 0:\n            ginis.append(0)\n        else:\n            idx = np.arange(1, n + 1)\n            ginis.append((2 * (idx * arr).sum()) / (n * arr.sum()) - (n + 1) / n)\n    return np.array(ginis)\n\ncluster_quality = []\nfor i in range(BEST_K):\n    member_idx = np.where(kshape_labels == i)[0]\n    member_data = pivot.iloc[member_idx].values\n    centroid = centroids[i]\n\n    # Intra-cluster distance (Euclidean — ตรงกับ KMeans)\n    distances = np.sqrt(((member_data - centroid) ** 2).sum(axis=1))\n    \n    # Peak week statistics\n    peak_weeks = member_data.argmax(axis=1)\n    pct_early_peak = (peak_weeks <= 6).mean()\n    \n    # Drop ratio statistics\n    drops = compute_drops(member_data)\n    \n    # Gini coefficient\n    ginis = compute_gini_array(member_data)\n    \n    # Early sales concentration\n    early_conc = member_data[:, :6].sum(axis=1) / (member_data.sum(axis=1) + 1e-8)\n    \n    cluster_quality.append({\n        'Cluster': i,\n        'N': len(member_idx),\n        'Is_Fad': '✅ FAD' if i in FAD_CLUSTERS else '❌ No',\n        'Avg_Dist': round(distances.mean(), 3),\n        'Median_Peak_Wk': int(np.median(peak_weeks)),\n        '%_Peak≤w6': f'{pct_early_peak:.0%}',\n        'Mean_Drop': f'{drops.mean():.0%}',\n        'Mean_Gini': round(ginis.mean(), 3),\n        'Mean_Early_Conc': f'{early_conc.mean():.0%}',\n    })\n\nquality_df = pd.DataFrame(cluster_quality)\nprint(quality_df.to_string(index=False))\n\nprint('\\nการอ่านตาราง:')\nprint('  Avg_Dist      : ระยะทางเฉลี่ยถึง centroid (ต่ำ = cluster แน่น ดี)')\nprint('  %_Peak≤w6     : สัดส่วนสมาชิกที่ peak ภายใน 6 สัปดาห์ (Fad ควรสูง)')\nprint('  Mean_Drop     : ยอดขายตกจาก peak เฉลี่ย (Fad ควรสูง)')\nprint('  Mean_Gini     : ความกระจุกตัวของยอดขาย (Fad ควรสูง)')\nprint('  Mean_Early_Conc: สัดส่วนยอดขายช่วง 6 สัปดาห์แรก (Fad ควรสูง)')\n\n# ─────────────────────────────────────────────────────────\n#  1.3 Statistical Test: Fad clusters ≠ Non-Fad clusters\n#      → Mann-Whitney U test (non-parametric)\n#      → ทดสอบว่าค่าเฉลี่ยของ 4 metrics ต่างกันอย่างมีนัยสำคัญ\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 1.3 Statistical Test: Fad vs Non-Fad ---\\n')\n\nfad_mask = np.isin(kshape_labels, FAD_CLUSTERS)\nfad_data = pivot.values[fad_mask]\nnonfad_data = pivot.values[~fad_mask]\n\nprint(f'Fad articles:     {fad_data.shape[0]:,}')\nprint(f'Non-Fad articles: {nonfad_data.shape[0]:,}\\n')\n\n# คำนวณ metrics สำหรับทั้ง 2 กลุ่ม\nfad_peaks = fad_data.argmax(axis=1)\nnonfad_peaks = nonfad_data.argmax(axis=1)\n\nfad_drops = compute_drops(fad_data)\nnonfad_drops = compute_drops(nonfad_data)\n\nfad_gini = compute_gini_array(fad_data)\nnonfad_gini = compute_gini_array(nonfad_data)\n\nfad_conc = fad_data[:, :6].sum(axis=1) / (fad_data.sum(axis=1) + 1e-8)\nnonfad_conc = nonfad_data[:, :6].sum(axis=1) / (nonfad_data.sum(axis=1) + 1e-8)\n\n# ทดสอบ 4 hypotheses\ntests = [\n    {\n        'name': 'Peak Week',\n        'hypothesis': 'Fad peak เร็วกว่า Non-Fad (Fad < Non-Fad)',\n        'fad_vals': fad_peaks, 'nonfad_vals': nonfad_peaks,\n        'alternative': 'less',\n        'fad_fmt': f'{fad_peaks.mean():.1f}', 'nonfad_fmt': f'{nonfad_peaks.mean():.1f}',\n    },\n    {\n        'name': 'Drop Ratio',\n        'hypothesis': 'Fad ดรอปมากกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_drops, 'nonfad_vals': nonfad_drops,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_drops.mean():.2%}', 'nonfad_fmt': f'{nonfad_drops.mean():.2%}',\n    },\n    {\n        'name': 'Gini Coefficient',\n        'hypothesis': 'Fad กระจุกตัวกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_gini, 'nonfad_vals': nonfad_gini,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_gini.mean():.3f}', 'nonfad_fmt': f'{nonfad_gini.mean():.3f}',\n    },\n    {\n        'name': 'Early Concentration',\n        'hypothesis': 'Fad ขายช่วงแรกมากกว่า Non-Fad (Fad > Non-Fad)',\n        'fad_vals': fad_conc, 'nonfad_vals': nonfad_conc,\n        'alternative': 'greater',\n        'fad_fmt': f'{fad_conc.mean():.2%}', 'nonfad_fmt': f'{nonfad_conc.mean():.2%}',\n    },\n]\n\ntest_results = []\nfor t in tests:\n    stat, pval = stats.mannwhitneyu(\n        t['fad_vals'], t['nonfad_vals'], alternative=t['alternative']\n    )\n    test_results.append({\n        'Test': t['name'],\n        'Hypothesis': t['hypothesis'],\n        'Fad Mean': t['fad_fmt'],\n        'Non-Fad Mean': t['nonfad_fmt'],\n        'p-value': f'{pval:.2e}',\n        'Significant?': '✅ Yes' if pval < 0.05 else '❌ No',\n    })\n\ntest_df = pd.DataFrame(test_results)\nprint(test_df.to_string(index=False))\nprint('\\n(Mann-Whitney U test, one-sided, α = 0.05)')\n\nn_pass = sum(1 for r in test_results if '✅' in r['Significant?'])\nprint(f'\\nผ่าน {n_pass}/4 tests', end='')\nif n_pass == 4:\n    print(' → ✅ ยืนยันว่า Fad clusters แตกต่างจาก Non-Fad อย่างมีนัยสำคัญทุกมิติ')\nelif n_pass >= 3:\n    print(' → 🟡 ส่วนใหญ่ผ่าน ยอมรับได้')\nelse:\n    print(' → 🔴 ผ่านน้อย ควรทบทวนการเลือก cluster')\n\n# ─────────────────────────────────────────────────────────\n#  1.4 Visual: Distribution Comparison (Fad vs Non-Fad)\n# ─────────────────────────────────────────────────────────\n\nfig, axes = plt.subplots(1, 4, figsize=(18, 4))\n\nviz_data = [\n    (fad_peaks, nonfad_peaks, 'Peak Week', 'Fad peak เร็วกว่า'),\n    (fad_drops, nonfad_drops, 'Drop Ratio', 'Fad ดรอปมากกว่า'),\n    (fad_gini, nonfad_gini, 'Gini Coefficient', 'Fad กระจุกตัวกว่า'),\n    (fad_conc, nonfad_conc, 'Early Concentration', 'Fad ขายช่วงแรกมากกว่า'),\n]\n\nfor ax, (fad_v, nonfad_v, title, interp) in zip(axes, viz_data):\n    ax.hist(nonfad_v, bins=30, alpha=0.5, color='steelblue', label='Non-Fad', density=True)\n    ax.hist(fad_v, bins=30, alpha=0.5, color='crimson', label='Fad', density=True)\n    ax.set_title(f'{title}\\n({interp})', fontsize=10, fontweight='bold')\n    ax.legend(fontsize=8)\n    ax.grid(alpha=0.3)\n\nplt.suptitle('Statistical Validation: Fad vs Non-Fad Distribution Comparison\\n'\n             '(2 กลุ่มต้องแยกออกจากกันชัดเจน)', fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n\n# ══════════════════════════════════════════════════════════════\n#  PART 2: PROFILING — ข้างใน cluster เป็นสินค้าอะไร?\n# ══════════════════════════════════════════════════════════════\n\nprint('\\n' + '=' * 70)\nprint('  PART 2: CLUSTER PROFILING')\nprint('  ข้างในแต่ละ cluster เป็นสินค้าอะไรบ้าง')\nprint('=' * 70)\n\n# Merge cluster labels กับข้อมูลสินค้าจาก articles.csv\nart_pd = art.toPandas()\nart_pd['article_id'] = art_pd['article_id'].astype(str)  # ← เพิ่มบรรทัดนี้\nprofile = kshape_df.merge(art_pd, on='article_id', how='left')\nprofile['is_fad_cluster'] = profile['kshape_cluster'].isin(FAD_CLUSTERS).map(\n    {True: 'Fad', False: 'Non-Fad'}\n)\n\n# ─────────────────────────────────────────────────────────\n#  2.1 Product Group Distribution per Cluster\n#      → ดูว่าแต่ละ cluster มีสินค้ากลุ่มไหนบ้าง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.1 Product Group Distribution (% ภายในแต่ละ cluster) ---\\n')\n\ngroup_dist = pd.crosstab(\n    profile['kshape_cluster'].map(\n        lambda x: f'C{x} {\"(FAD)\" if x in FAD_CLUSTERS else \"(Non-Fad)\"}'\n    ),\n    profile['product_group_name'],\n    normalize='index'\n).round(3) * 100\n\nprint(group_dist.to_string())\n\n# ─────────────────────────────────────────────────────────\n#  2.2 Fad Rate by Product Group\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.2 Fad Rate by Product Group ---\\n')\n\nfad_by_group = (\n    profile.groupby('product_group_name')\n    .agg(\n        total=('article_id', 'count'),\n        n_fad=('is_fad_cluster', lambda x: (x == 'Fad').sum())\n    )\n)\nfad_by_group['fad_rate'] = (fad_by_group['n_fad'] / fad_by_group['total'] * 100).round(1)\nfad_by_group = fad_by_group.sort_values('fad_rate', ascending=False)\n\nprint(fad_by_group.to_string())\n\n# Plot\nfig, ax = plt.subplots(figsize=(10, 5))\nfad_by_group['fad_rate'].plot(kind='barh', ax=ax, color='crimson', alpha=0.7)\nax.set_xlabel('Fad Rate (%)')\nax.set_title('Fad Rate by Product Group\\n(กลุ่มสินค้าไหนมีสัดส่วน Fad สูงที่สุด)', fontweight='bold')\noverall_rate = profile['is_fad_cluster'].eq('Fad').mean() * 100\nax.axvline(overall_rate, color='gray', linestyle='--',\n           label=f'Overall rate ({overall_rate:.1f}%)')\nax.legend()\nax.grid(alpha=0.3)\nplt.tight_layout()\nplt.show()\n\n# ─────────────────────────────────────────────────────────\n#  2.3 Top Product Types ใน Fad vs Non-Fad\n#      → ประเภทสินค้าเฉพาะเจาะจงที่เป็น Fad มากที่สุด\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.3 Top 15 Product Types in Fad Clusters ---\\n')\n\nfad_articles = profile[profile['is_fad_cluster'] == 'Fad']\nfad_type_counts = fad_articles['product_type_name'].value_counts().head(15)\n\nprint(f'{\"Product Type\":<40s} {\"Count\":>8s} {\"% of Fad\":>10s}')\nprint('-' * 60)\nfor ptype, count in fad_type_counts.items():\n    pct = count / len(fad_articles) * 100\n    print(f'{ptype:<40s} {count:>8,} {pct:>9.1f}%')\n\n# ─────────────────────────────────────────────────────────\n#  2.4 Color Distribution: Fad vs Non-Fad\n#      → สีอะไรมี Fad สูง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.4 Color Distribution: Fad vs Non-Fad ---\\n')\n\nif 'colour_group_name' in profile.columns:\n    color_comp = pd.crosstab(\n        profile['is_fad_cluster'],\n        profile['colour_group_name'],\n        normalize='index'\n    ).T\n\n    color_comp['diff'] = color_comp.get('Fad', 0) - color_comp.get('Non-Fad', 0)\n    color_comp = color_comp.sort_values('diff', ascending=False)\n\n    print('Top 10 สีที่ Fad มีสัดส่วนสูงกว่า Non-Fad:')\n    top_colors = color_comp.head(10)\n    for color_name, row in top_colors.iterrows():\n        fad_pct = row.get('Fad', 0) * 100\n        nonfad_pct = row.get('Non-Fad', 0) * 100\n        print(f'  {color_name:<25s}  Fad: {fad_pct:5.1f}%  Non-Fad: {nonfad_pct:5.1f}%  '\n              f'(diff: {(fad_pct - nonfad_pct):+.1f}%)')\nelse:\n    print('  colour_group_name column not found')\n\n# ─────────────────────────────────────────────────────────\n#  2.5 Graphical Appearance: Fad vs Non-Fad\n#      → ลายแบบไหนมี Fad สูง\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.5 Graphical Appearance: Fad vs Non-Fad ---\\n')\n\nif 'graphical_appearance_name' in profile.columns:\n    appear_comp = pd.crosstab(\n        profile['is_fad_cluster'],\n        profile['graphical_appearance_name'],\n        normalize='index'\n    ).T\n\n    appear_comp['diff'] = appear_comp.get('Fad', 0) - appear_comp.get('Non-Fad', 0)\n    appear_comp = appear_comp.sort_values('diff', ascending=False)\n\n    print('Top 10 ลายที่ Fad มีสัดส่วนสูงกว่า Non-Fad:')\n    for appear_name, row in appear_comp.head(10).iterrows():\n        fad_pct = row.get('Fad', 0) * 100\n        nonfad_pct = row.get('Non-Fad', 0) * 100\n        print(f'  {appear_name:<25s}  Fad: {fad_pct:5.1f}%  Non-Fad: {nonfad_pct:5.1f}%  '\n              f'(diff: {(fad_pct - nonfad_pct):+.1f}%)')\nelse:\n    print('  graphical_appearance_name column not found')\n\n# ─────────────────────────────────────────────────────────\n#  2.6 Department Fad Rate (Top 10, min 20 articles)\n# ─────────────────────────────────────────────────────────\n\nprint('\\n--- 2.6 Department Fad Rate (Top 10, min 20 articles) ---\\n')\n\nif 'department_name' in profile.columns:\n    dept_fad = (\n        profile.groupby('department_name')\n        .agg(\n            total=('article_id', 'count'),\n            n_fad=('is_fad_cluster', lambda x: (x == 'Fad').sum())\n        )\n    )\n    dept_fad['fad_rate'] = (dept_fad['n_fad'] / dept_fad['total'] * 100).round(1)\n    dept_fad = dept_fad.query('total >= 20').sort_values('fad_rate', ascending=False).head(10)\n    print(dept_fad.to_string())\nelse:\n    print('  department_name column not found')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:32:51.760211Z","iopub.execute_input":"2026-04-27T20:32:51.760805Z","iopub.status.idle":"2026-04-27T20:33:00.097894Z","shell.execute_reply.started":"2026-04-27T20:32:51.760772Z","shell.execute_reply":"2026-04-27T20:33:00.097203Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 14: Hybrid Label ──\n#\n#  ✅ ใช้ FAD_CLUSTERS, BEST_K, centroids จาก Cell 13 อัตโนมัติ\n#  ✅ ไม่มี ks (KShape) reference — ใช้ centroids dict แทน\n#  ✅ ไม่มี manual override\n# ═══════════════════════════════════════════════════════════════════\n\nfrom sklearn.metrics import cohen_kappa_score\n\n# ─── ใช้ค่าจาก Cell 13 โดยตรง (ไม่ manual) ───\nprint(f'FAD_CLUSTERS = {FAD_CLUSTERS}')\nprint(f'BEST_K = {BEST_K}\\n')\n\n# แสดงรายละเอียดแต่ละ Fad cluster\nfor c in sorted(FAD_CLUSTERS):\n    n = int((kshape_labels == c).sum())\n    peak_w = int(np.argmax(centroids[c]))\n    print(f'   Cluster {c} (peak@w{peak_w}): {n:,} articles')\n\nn_fad_kshape = sum((kshape_labels == c).sum() for c in FAD_CLUSTERS)\nprint(f'   Total Fad members (k-Shape): {n_fad_kshape:,} ({n_fad_kshape/len(pivot)*100:.1f}%)')\n\n# แสดง Non-Fad clusters\nnon_fad_clusters = [c for c in range(BEST_K) if c not in FAD_CLUSTERS]\nprint(f'\\n   Non-Fad clusters: {non_fad_clusters}')\nfor c in non_fad_clusters:\n    n = int((kshape_labels == c).sum())\n    peak_w = int(np.argmax(centroids[c]))\n    print(f'   Cluster {c} (peak@w{peak_w}): {n:,} articles')\n\n# ─── Merge PELT + k-Shape labels ───\nlabels = pelt_df.merge(kshape_df, on='article_id', how='inner')\n\nlabels['kshape_fad'] = labels['kshape_cluster'].isin(FAD_CLUSTERS).astype(int)\n\nlabels['fad_label'] = (\n    (labels['pelt_fad'] == 1) & (labels['kshape_fad'] == 1)\n).astype(int)\n\n# ═══════════════════════════════════════════════════════\n#  Agreement Matrix\n# ═══════════════════════════════════════════════════════\n\nprint('\\n=== Label Agreement Matrix ===')\nagreement = pd.crosstab(\n    labels['pelt_fad'].map({0: 'PELT=0', 1: 'PELT=1'}),\n    labels['kshape_fad'].map({0: 'kShape=0', 1: 'kShape=1'}),\n    margins=True, margins_name='Total'\n)\nprint(agreement)\n\n# ═══════════════════════════════════════════════════════\n#  Rates Comparison\n# ═══════════════════════════════════════════════════════\n\npelt_rate   = labels['pelt_fad'].mean()\nkshape_rate = labels['kshape_fad'].mean()\nhybrid_rate = labels['fad_label'].mean()\n\nprint(f'\\n=== Rates Comparison ===')\nprint(f'PELT alone       : {pelt_rate:.2%}  ({labels[\"pelt_fad\"].sum():>6,} articles)')\nprint(f'k-Shape alone    : {kshape_rate:.2%}  ({labels[\"kshape_fad\"].sum():>6,} articles)')\nprint(f'Hybrid (AND)     : {hybrid_rate:.2%}  ({labels[\"fad_label\"].sum():>6,} articles)')\nprint(f'Phase 1 baseline : 4.10%   (death ≤ 6 weeks estimate)')\n\n# ═══════════════════════════════════════════════════════\n#  Cohen's Kappa\n# ═══════════════════════════════════════════════════════\n\nkappa = cohen_kappa_score(labels['pelt_fad'], labels['kshape_fad'])\nprint(f'\\n=== Cohen\\'s Kappa (Inter-Method Agreement) ===')\nprint(f'  Kappa = {kappa:.3f}')\nif kappa >= 0.60:\n    print('  ✅ Substantial agreement')\nelif kappa >= 0.40:\n    print('  🟡 Moderate agreement')\nelif kappa >= 0.20:\n    print('  🟠 Fair agreement')\nelse:\n    print('  🔴 Weak agreement')\n\n# ═══════════════════════════════════════════════════════\n#  Label Composition by Cluster\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Fad Label Composition by Cluster ===')\nfad_composition = labels[labels['fad_label'] == 1].groupby('kshape_cluster').size()\nfor cluster, count in fad_composition.items():\n    pct = count / labels['fad_label'].sum() * 100\n    peak_w = int(np.argmax(centroids[cluster]))\n    print(f'  Cluster {cluster} (peak@w{peak_w}): {count:>5,} articles ({pct:.1f}% of all Fads)')\n\n# ═══════════════════════════════════════════════════════\n#  Temporal Profile\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Final Fad Temporal Profile ===')\nfinal_fads = labels[labels['fad_label'] == 1]\nif len(final_fads) > 0:\n    print(f'  Peak week   — median: {final_fads[\"peak_week\"].median():.0f}, '\n          f'mean: {final_fads[\"peak_week\"].mean():.1f}')\n    print(f'  Peak ratio  — mean: {final_fads[\"peak_ratio\"].mean():.2f}')\n    print(f'  Drop ratio  — mean: {final_fads[\"drop_ratio\"].mean():.2f}')\n    print(f'  Tail ratio  — mean: {final_fads[\"tail_ratio\"].mean():.3f}')\n\n# ═══════════════════════════════════════════════════════\n#  Final Distribution\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Final Hybrid Label Distribution ===')\nprint(labels['fad_label'].value_counts())\nprint(f'Final Fad rate: {hybrid_rate:.2%}')\n\n# ═══════════════════════════════════════════════════════\n#  Sanity Check\n# ═══════════════════════════════════════════════════════\n\nprint(f'\\n=== Sanity Check ===')\nif hybrid_rate < 0.02:\n    print('  🔴 ต่ำเกินไป < 2%')\nelif hybrid_rate > 0.15:\n    print('  🟡 สูงกว่าเป้าเล็กน้อย > 15%')\nelif 0.03 <= hybrid_rate <= 0.12:\n    print('  ✅ Good — อยู่ในช่วงที่เหมาะสม')\nelse:\n    print(f'  🟡 Acceptable — {hybrid_rate:.2%}')\n\nprint('\\n✅ Phase 2 complete — ready for Phase 3')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:37:37.540921Z","iopub.execute_input":"2026-04-27T20:37:37.541676Z","iopub.status.idle":"2026-04-27T20:37:37.679039Z","shell.execute_reply.started":"2026-04-27T20:37:37.541643Z","shell.execute_reply":"2026-04-27T20:37:37.678402Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 3 — Feature Engineering with PySpark","metadata":{}},{"cell_type":"code","source":"# ── Cell 15: Dynamic Features  ──\n\nEARLY_WINDOW = 5\nMAX_WEEKS = 26\n\n# ═══════════════════════════════════════════════════════\n#  Part A: Weekly sales in early window (Spark parallel)\n# ═══════════════════════════════════════════════════════\n\nearly_spark = (\n    weekly\n    .filter(\n        (F.col('weeks_since_launch') >= 0) &\n        (F.col('weeks_since_launch') <= EARLY_WINDOW)\n    )\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(EARLY_WINDOW + 1)))\n    .agg(F.first('weekly_sales_clean'))\n    .fillna(0)\n)\n\n# Rename: '0' → 'sales_w0', etc.\nearly_cols = early_spark.columns\nearly_spark = early_spark.toDF(\n    *(['article_id'] + [f'sales_w{c}' for c in early_cols[1:]])\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part B: Aggregate features from early window\n# ═══════════════════════════════════════════════════════\n\nearly_agg = (\n    weekly\n    .filter(F.col('weeks_since_launch') <= EARLY_WINDOW)\n    .groupBy('article_id')\n    .agg(\n        F.max('weekly_sales_clean').alias('early_peak_sales'),\n        F.sum('weekly_sales_clean').alias('early_total_sales'),\n        F.sum('unique_buyers').alias('early_unique_buyers'),\n        F.max('n_channels').alias('early_max_channels'),\n        F.avg('weekly_sales_clean').alias('early_avg_sales'),\n    )\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part C: Full lifecycle aggregates (within 26 weeks)\n# ═══════════════════════════════════════════════════════\n\nlifecycle_agg = (\n    weekly\n    .filter(F.col('weeks_since_launch') < MAX_WEEKS)\n    .groupBy('article_id')\n    .agg(\n        F.max('weekly_sales_clean').alias('peak_sales'),\n        F.sum('weekly_sales_clean').alias('total_sales_clean'),\n        F.avg('weekly_sales_clean').alias('avg_sales'),\n    )\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part D: time_to_peak (within 26 weeks ONLY)\n# ═══════════════════════════════════════════════════════\n\nw_peak = Window.partitionBy('article_id').orderBy(F.desc('weekly_sales_clean'))\n\ntime_to_peak = (\n    weekly\n    .filter(F.col('weeks_since_launch') < MAX_WEEKS)  # ✅ จำกัด 26 weeks\n    .withColumn('rn', F.row_number().over(w_peak))\n    .filter(F.col('rn') == 1)\n    .select('article_id', F.col('weeks_since_launch').alias('time_to_peak'))\n)\n\n# ═══════════════════════════════════════════════════════\n#  Part E: post_peak_decay (จาก pivot — accurate)\n# ═══════════════════════════════════════════════════════\n\n# คำนวณบน Pandas (pivot มีแค่ ~80K rows × 26 cols = เล็ก)\npeak_weeks_pd = time_to_peak.toPandas().set_index('article_id')\n\ndecay_results = []\nfor art_id, row in pivot.iterrows():\n    series = row.values.astype(float)\n    \n    # Get peak week from pre-computed time_to_peak\n    if art_id in peak_weeks_pd.index:\n        pw = int(peak_weeks_pd.loc[art_id, 'time_to_peak'])\n    else:\n        pw = int(np.argmax(series))\n    \n    peak_val = series[pw]\n    \n    # Compute decay slope: peak → 3 weeks after peak\n    if pw < MAX_WEEKS - 3:\n        end_w = pw + 3\n        end_val = series[end_w]\n        decay = (end_val - peak_val) / 3.0\n    elif pw < MAX_WEEKS - 1:\n        end_val = series[-1]\n        decay = (end_val - peak_val) / max(1, MAX_WEEKS - 1 - pw)\n    else:\n        decay = 0.0\n    \n    # Is peak in early window? (strong Fad indicator)\n    is_early_peak = int(pw <= EARLY_WINDOW)\n    \n    decay_results.append({\n        'article_id': art_id,\n        'post_peak_decay': round(decay, 4),\n        'is_early_peak_flag': is_early_peak,\n    })\n\ndecay_spark = spark.createDataFrame(pd.DataFrame(decay_results))\n\n# ═══════════════════════════════════════════════════════\n#  Part F: Combine all + compute derived features\n# ═══════════════════════════════════════════════════════\n\ndyn = (\n    early_spark\n    .join(early_agg,      on='article_id', how='left')\n    .join(lifecycle_agg,  on='article_id', how='left')\n    .join(time_to_peak,   on='article_id', how='left')\n    .join(decay_spark,    on='article_id', how='left')\n    \n    # ── Velocity & Acceleration ──\n    .withColumn('velocity_w1_w0',\n        F.col('sales_w1') / (F.col('sales_w0') + 1.0))\n    .withColumn('velocity_w2_w1',\n        F.col('sales_w2') / (F.col('sales_w1') + 1.0))\n    .withColumn('acceleration',\n        (F.col('sales_w2') - F.col('sales_w1')) - \n        (F.col('sales_w1') - F.col('sales_w0')))\n    \n    # ── Peak ratios ──\n    .withColumn('peak_to_w0_ratio',\n        F.col('peak_sales') / (F.col('sales_w0') + 1.0))\n    .withColumn('early_peak_to_avg',\n        F.col('early_peak_sales') / (F.col('early_avg_sales') + 1.0))\n    \n    # ── Pre-peak velocity ──\n    .withColumn('pre_peak_velocity',\n        F.when(F.col('time_to_peak') > 0,\n            F.col('early_peak_sales') / (F.col('time_to_peak') + 1.0)\n        ).otherwise(F.col('sales_w0')))\n    \n    # ── Early concentration (% of total sales in early window) ──\n    .withColumn('early_sales_concentration',\n        F.col('early_total_sales') / (F.col('total_sales_clean') + 1.0))\n    \n    # ── Buyer intensity ──\n    .withColumn('early_buyer_intensity',\n        F.col('early_unique_buyers') / (F.col('early_total_sales') + 1.0))\n)\n\ndyn.cache()\nn_features = len(dyn.columns) - 1  # minus article_id\n\n# ═══════════════════════════════════════════════════════\n#  Diagnostics\n# ═══════════════════════════════════════════════════════\n\nprint(f'✅ Dynamic features computed.')\nprint(f'   Early window: weeks 0-{EARLY_WINDOW}')\nprint(f'   Features: {n_features}')\n\n# Verify time_to_peak is within 0-25\nprint('\\n=== time_to_peak Distribution (should be 0-25) ===')\ndyn.select(\n    F.min('time_to_peak').alias('min'),\n    F.max('time_to_peak').alias('max'),\n    F.mean('time_to_peak').alias('mean'),\n    F.expr('percentile_approx(time_to_peak, 0.5)').alias('median'),\n).show()\n\n# Verify post_peak_decay is not all zeros\nprint('=== post_peak_decay Distribution ===')\ndyn.select(\n    F.mean('post_peak_decay').alias('mean'),\n    F.min('post_peak_decay').alias('min'),\n    F.max('post_peak_decay').alias('max'),\n    F.sum(F.when(F.col('post_peak_decay') != 0, 1).otherwise(0)).alias('non_zero_count'),\n    F.count('*').alias('total'),\n).show()\n\n# Preview\nprint('--- Sample output ---')\ndyn.select('article_id', 'velocity_w1_w0', 'time_to_peak', \n           'post_peak_decay', 'is_early_peak_flag',\n           'early_sales_concentration').show(10)\n\n# Fad-likely articles preview\nn_early_peak = dyn.filter(F.col('is_early_peak_flag') == 1).count()\nprint(f'\\nArticles with early peak (≤ week {EARLY_WINDOW}): {n_early_peak:,} '\n      f'({n_early_peak/dyn.count()*100:.1f}%)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:41:03.762394Z","iopub.execute_input":"2026-04-27T20:41:03.762731Z","iopub.status.idle":"2026-04-27T20:41:19.436256Z","shell.execute_reply.started":"2026-04-27T20:41:03.762705Z","shell.execute_reply":"2026-04-27T20:41:19.435009Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 16: Gini Coefficient  ──\n\nfrom pyspark.sql.types import DoubleType\n\n@F.udf(DoubleType())\ndef gini_udf(sales_list):\n    \"\"\"Gini coefficient: 0 = equal distribution, 1 = all sales in one week\"\"\"\n    if not sales_list:\n        return 0.0\n    arr = np.sort(np.abs(np.array(sales_list, dtype=float)))\n    n = len(arr)\n    if n == 0 or arr.sum() == 0:\n        return 0.0\n    idx = np.arange(1, n + 1)\n    return float((2 * (idx * arr).sum()) / (n * arr.sum()) - (n + 1) / n)\n\ngini_spark = (\n    weekly.groupBy('article_id')\n    .agg(F.collect_list('weekly_sales_clean').alias('sales_list'))\n    .withColumn('gini_coeff', gini_udf('sales_list'))\n    .select('article_id', 'gini_coeff')\n)\n\ngini_spark.cache()\nprint('✅ Gini computed.')\n\n# Diagnostic: Fad ควรมี Gini สูง (sales กระจุกตัว)\ngini_pd = gini_spark.toPandas()\nprint(f'   Mean Gini: {gini_pd[\"gini_coeff\"].mean():.3f}')\nprint(f'   Std:  {gini_pd[\"gini_coeff\"].std():.3f}')\ngini_spark.show(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:46:39.554038Z","iopub.execute_input":"2026-04-27T20:46:39.554860Z","iopub.status.idle":"2026-04-27T20:46:44.919960Z","shell.execute_reply.started":"2026-04-27T20:46:39.554829Z","shell.execute_reply":"2026-04-27T20:46:44.919152Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 17: Repeat Purchase + Interpurchase CV ──\n\n# ─── Part A: Repeat Purchase Rate (per article) ───\nmulti_buy = (\n    txn.groupBy('article_id', 'customer_id')\n    .agg(F.count('*').alias('n_purchases'))\n)\n\nbuyer_stats = (\n    multi_buy.groupBy('article_id')\n    .agg(\n        F.count('customer_id').alias('total_buyers'),\n        F.sum((F.col('n_purchases') >= 2).cast('int')).alias('repeat_buyers'),\n        F.avg('n_purchases').alias('avg_purchases_per_buyer'),\n    )\n    .withColumn('repeat_purchase_rate',\n        F.col('repeat_buyers') / (F.col('total_buyers') + 1.0))\n)\n\n# ─── Part B: Inter-purchase CV (per article — FIXED) ───\n# นี่วัดว่า \"transaction ของ article เกิดถี่แค่ไหนและสม่ำเสมอแค่ไหน\"\n\nw_art_txn = Window.partitionBy('article_id').orderBy('t_dat')\n\ncv_spark = (\n    txn\n    .withColumn('row_num', F.row_number().over(w_art_txn))\n    .withColumn('lag_date', F.lag('t_dat', 1).over(w_art_txn))\n    .withColumn('gap_days', F.datediff('t_dat', 'lag_date').cast('double'))\n    .filter(F.col('gap_days').isNotNull())\n    .groupBy('article_id')\n    .agg(\n        F.mean('gap_days').alias('mean_gap'),\n        F.stddev('gap_days').alias('std_gap'),\n        F.count('gap_days').alias('n_gaps'),\n    )\n    # CV = std / mean (high CV = burst pattern = Fad signal)\n    .withColumn('interpurchase_cv',\n        F.when(F.col('n_gaps') >= 3,  # ต้องมี gaps พอ ≥ 3 ถึงจะ compute CV ได้\n            F.col('std_gap') / (F.col('mean_gap') + 1.0)\n        ).otherwise(0.0))\n    .select('article_id', 'mean_gap', 'interpurchase_cv')\n)\n\nbuyer_stats.cache()\ncv_spark.cache()\nprint('✅ Buyer behavior features done.')\nprint(f'   Buyer stats rows: {buyer_stats.count():,}')\nprint(f'   CV stats rows: {cv_spark.count():,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:46:51.866655Z","iopub.execute_input":"2026-04-27T20:46:51.867277Z","iopub.status.idle":"2026-04-27T20:48:37.451472Z","shell.execute_reply.started":"2026-04-27T20:46:51.867246Z","shell.execute_reply":"2026-04-27T20:48:37.450602Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 18: Static Features + Target Encoding ──\n\n# ─── Part A: Select static features ───\nart_feat_spark = art.select(\n    'article_id',\n    'product_type_no',           # 131 unique → Target Encoding\n    'department_no',             # 250 unique → Target Encoding\n    'graphical_appearance_no',   # 30 unique → keep as-is (OHE later if needed)\n    'colour_group_code',         # 50 unique → keep as-is\n    'perceived_colour_value_id',\n    'perceived_colour_master_id',\n    'index_code',                # 10 unique → keep as-is\n    'index_group_no',\n    'section_no',\n    'garment_group_no',          # 21 unique → keep as-is\n)\n\n# ─── Part B: Target Encoding for high-cardinality features ───\n# ✅ Phase 0 กำหนดไว้: department_no (250) + product_type_no (131) → Target Encoding\n# Logic: replace category code with mean(fad_label) of that category\n# ⚠️ ต้องทำบน TRAIN set only ใน Phase 5 เพื่อป้องกัน data leakage\n# ตอนนี้คำนวณ global mean ไว้ก่อน → Phase 5 จะ re-compute per fold\n\nlabels_spark = spark.createDataFrame(\n    labels[['article_id', 'fad_label']].astype({'fad_label': 'int32'})\n)\n\n# Target encoding: department_no\ndept_target = (\n    art_feat_spark\n    .join(labels_spark, 'article_id', 'inner')\n    .groupBy('department_no')\n    .agg(\n        F.mean('fad_label').alias('dept_fad_rate'),\n        F.count('*').alias('dept_count'),\n    )\n)\n\n# Target encoding: product_type_no\nptype_target = (\n    art_feat_spark\n    .join(labels_spark, 'article_id', 'inner')\n    .groupBy('product_type_no')\n    .agg(\n        F.mean('fad_label').alias('ptype_fad_rate'),\n        F.count('*').alias('ptype_count'),\n    )\n)\n\n# Smoothed target encoding: smooth toward global mean for rare categories\nglobal_fad_rate = labels['fad_label'].mean()\nSMOOTH_FACTOR = 20  # categories with < 20 items → lean toward global mean\n\ndept_target = dept_target.withColumn('dept_target_enc',\n    (F.col('dept_fad_rate') * F.col('dept_count') + global_fad_rate * SMOOTH_FACTOR) /\n    (F.col('dept_count') + SMOOTH_FACTOR)\n)\n\nptype_target = ptype_target.withColumn('ptype_target_enc',\n    (F.col('ptype_fad_rate') * F.col('ptype_count') + global_fad_rate * SMOOTH_FACTOR) /\n    (F.col('ptype_count') + SMOOTH_FACTOR)\n)\n\n# Join target encodings back\nart_feat_spark = (\n    art_feat_spark\n    .join(dept_target.select('department_no', 'dept_target_enc'), 'department_no', 'left')\n    .join(ptype_target.select('product_type_no', 'ptype_target_enc'), 'product_type_no', 'left')\n)\n\n# ─── Part C: Interaction features ───\n# All-over pattern × bright colour (from methodology)\nart_feat_spark = (\n    art_feat_spark\n    .withColumn('is_allover_pattern',\n        (F.col('graphical_appearance_no').isin([74])).cast('int'))\n    .withColumn('is_bright_colour',\n        (F.col('perceived_colour_value_id').isin([3, 4])).cast('int'))\n    .withColumn('allover_x_bright',\n        F.col('is_allover_pattern') * F.col('is_bright_colour'))\n)\n\nart_feat_spark.cache()\nprint('✅ Static features done.')\nprint(f'   Columns: {len(art_feat_spark.columns)}')\nprint(f'   Target encoding: department_no, product_type_no (smoothed, factor={SMOOTH_FACTOR})')\nprint(f'   Global fad rate for smoothing: {global_fad_rate:.4f}')\n\n# Preview target encoding\nprint('\\n--- Department Target Encoding (top 5 highest Fad rate) ---')\ndept_target.orderBy(F.desc('dept_fad_rate')).show(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:50:19.999196Z","iopub.execute_input":"2026-04-27T20:50:19.999557Z","iopub.status.idle":"2026-04-27T20:50:21.801005Z","shell.execute_reply.started":"2026-04-27T20:50:19.999527Z","shell.execute_reply":"2026-04-27T20:50:21.800231Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 19 REVISED: Merge All Features ──\n\n# ─── Part A: Stockout features per article (from Phase 1) ───\nstockout_features = (\n    weekly\n    .groupBy('article_id')\n    .agg(\n        F.sum('is_stockout').alias('n_stockout_weeks'),\n        F.max('stockout_tier').alias('max_stockout_tier'),\n        # Stockout rate in early window\n        F.sum(\n            F.when(\n                (F.col('weeks_since_launch') <= EARLY_WINDOW) & (F.col('is_stockout') == 1), 1\n            ).otherwise(0)\n        ).alias('early_stockout_count'),\n    )\n)\n\n# ─── Part B: Labels (from Phase 2) ───\nlabels_spark = spark.createDataFrame(\n    labels[['article_id', 'fad_label', 'peak_week', 'kshape_cluster',\n            'peak_ratio', 'drop_ratio', 'tail_ratio']].copy()\n)\n\n# ─── Part C: Regime info (from Phase 1) ───\nregime_spark = total_sales_df.select('article_id', 'regime', 'strat_key', 'total_sales')\n\n# ─── Part D: Merge everything ───\nfeatures_spark = (\n    labels_spark.select('article_id', 'fad_label')  # target + ID\n    .join(dyn,              on='article_id', how='left')   # dynamic features\n    .join(gini_spark,       on='article_id', how='left')   # gini\n    .join(buyer_stats.select('article_id', 'total_buyers', \n                             'repeat_purchase_rate', 'avg_purchases_per_buyer'),\n                            on='article_id', how='left')   # buyer behavior\n    .join(cv_spark,         on='article_id', how='left')   # interpurchase CV\n    .join(art_feat_spark,   on='article_id', how='left')   # static features\n    .join(stockout_features,on='article_id', how='left')   # stockout features\n    .join(regime_spark,     on='article_id', how='left')   # regime\n    .fillna(0)\n)\n\nfeatures_spark.cache()\nn_rows = features_spark.count()\nn_cols = len(features_spark.columns)\n\nprint(f'✅ Feature matrix assembled.')\nprint(f'   Rows: {n_rows:,}')\nprint(f'   Columns: {n_cols}')\n\n# ─── Diagnostic: Feature categories ───\nfeature_cols = [c for c in features_spark.columns if c not in ['article_id', 'fad_label']]\nprint(f'   Feature columns: {len(feature_cols)}')\n\n# Check for unexpected nulls\nnull_counts = features_spark.select([\n    F.sum(F.col(c).isNull().cast('int')).alias(c) for c in feature_cols\n]).collect()[0].asDict()\nnulls = {k: v for k, v in null_counts.items() if v > 0}\nif nulls:\n    print(f'\\n   ⚠️  Columns with nulls after fillna: {nulls}')\nelse:\n    print(f'   ✅ No nulls remaining')\n\n# ─── Collect to Pandas ───\nfeatures = features_spark.toPandas()\nprint(f'\\n   Pandas shape: {features.shape}')\n\n# ─── Label distribution in final feature matrix ───\nprint(f'\\n=== Label Distribution ===')\nprint(features['fad_label'].value_counts())\nprint(f'Fad rate: {features[\"fad_label\"].mean():.2%}')\n\n# ─── Feature summary by type ───\nprint(f'\\n=== Feature Summary ===')\ndynamic_cols = ['sales_w0','sales_w1','sales_w2','sales_w3','sales_w4','sales_w5',\n                'velocity_w1_w0','velocity_w2_w1','acceleration',\n                'peak_to_w0_ratio','early_peak_to_avg','post_peak_decay',\n                'pre_peak_velocity','early_sales_concentration','early_buyer_intensity',\n                'time_to_peak','early_peak_sales','early_total_sales',\n                'early_unique_buyers','early_max_channels','early_avg_sales',\n                'peak_sales','total_sales_clean','avg_sales']\nbehavioral_cols = ['gini_coeff','total_buyers','repeat_purchase_rate',\n                   'avg_purchases_per_buyer','interpurchase_cv','mean_gap']\nstatic_cols = ['product_type_no','department_no','graphical_appearance_no',\n               'colour_group_code','perceived_colour_value_id','perceived_colour_master_id',\n               'index_code','index_group_no','section_no','garment_group_no',\n               'dept_target_enc','ptype_target_enc',\n               'is_allover_pattern','is_bright_colour','allover_x_bright']\noperational_cols = ['n_stockout_weeks','max_stockout_tier','early_stockout_count']\n\nfor name, cols in [('Dynamic (early signals)', dynamic_cols),\n                    ('Behavioral (buyer)', behavioral_cols),\n                    ('Static (article metadata)', static_cols),\n                    ('Operational (stockout)', operational_cols)]:\n    present = [c for c in cols if c in features.columns]\n    print(f'  {name}: {len(present)} features')\n\nprint(f'\\n✅ Phase 3 complete — ready for Phase 4 (Feature Selection)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:50:26.526651Z","iopub.execute_input":"2026-04-27T20:50:26.526949Z","iopub.status.idle":"2026-04-27T20:50:43.941545Z","shell.execute_reply.started":"2026-04-27T20:50:26.526923Z","shell.execute_reply":"2026-04-27T20:50:43.940578Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 4 — Feature Selection","metadata":{}},{"cell_type":"code","source":"# ── Cell 20: Feature Selection (MI + Correlation + RFE) ──\n\nfrom sklearn.feature_selection import RFE, mutual_info_classif\nimport lightgbm as lgb\n\n# ─── Step 0: Prepare X and y ───\n# Columns ที่ไม่ใช่ feature (IDs, labels, metadata)\nDROP_COLS = [\n    'article_id', 'fad_label',              # ID + target\n    'regime', 'strat_key', 'total_sales',   # Phase 1 metadata (ใช้สำหรับ stratification ไม่ใช่ prediction)\n]\n\n# ตรวจสอบว่ามี columns ไหนอยู่จริง\nexisting_drops = [c for c in DROP_COLS if c in features.columns]\nprint(f'Dropping {len(existing_drops)} non-feature columns: {existing_drops}')\n\nX_df = features.drop(columns=existing_drops, errors='ignore').copy()\ny_all = features['fad_label'].values\n\n# Encode categorical columns (object/category → numeric)\ncat_cols = X_df.select_dtypes(include=['object', 'category']).columns.tolist()\nif cat_cols:\n    print(f'Encoding {len(cat_cols)} categorical columns: {cat_cols}')\n    for col in cat_cols:\n        X_df[col] = X_df[col].astype('category').cat.codes\n\nprint(f'\\nStarting feature selection: {X_df.shape[1]} features, {len(y_all):,} samples')\nprint(f'Fad rate: {y_all.mean():.2%}')\n\n# ═══════════════════════════════════════════════════════\n#  Step 1: Mutual Information (非linear relationship)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 1: Mutual Information')\nprint('='*50)\n\nmi = mutual_info_classif(X_df, y_all, random_state=SEED, n_neighbors=5)\nmi_df = pd.DataFrame({\n    'feature': X_df.columns, \n    'mi': mi\n}).sort_values('mi', ascending=False)\n\n# Threshold: keep features with MI > 0.001\nMI_THRESH = 0.001\nsel_mi = mi_df.loc[mi_df['mi'] >= MI_THRESH, 'feature'].tolist()\ndropped_mi = mi_df.loc[mi_df['mi'] < MI_THRESH, 'feature'].tolist()\n\nprint(f'MI threshold: {MI_THRESH}')\nprint(f'Features kept: {len(sel_mi)}')\nif dropped_mi:\n    print(f'Features dropped (no predictive power): {dropped_mi}')\n\n# Visualization: MI scores\nfig, ax = plt.subplots(figsize=(10, max(6, len(mi_df) * 0.3)))\ncolors = ['crimson' if mi >= MI_THRESH else 'lightgray' for mi in mi_df['mi']]\nax.barh(range(len(mi_df)), mi_df['mi'].values, color=colors)\nax.set_yticks(range(len(mi_df)))\nax.set_yticklabels(mi_df['feature'].values, fontsize=8)\nax.axvline(MI_THRESH, color='red', linestyle='--', alpha=0.7, label=f'Threshold={MI_THRESH}')\nax.set_xlabel('Mutual Information')\nax.set_title('Step 1: Mutual Information Scores\\n(red = kept, gray = dropped)', fontweight='bold')\nax.legend()\nax.invert_yaxis()\nplt.tight_layout()\nplt.show()\n\n# Top 10 features by MI\nprint('\\nTop 10 features by MI:')\nprint(mi_df.head(10).to_string(index=False))\n\n# ═══════════════════════════════════════════════════════\n#  Step 2: Correlation Pruning (remove redundancy)\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 2: Correlation Pruning')\nprint('='*50)\n\nCORR_THRESH = 0.85\n\ncorr = X_df[sel_mi].astype('float32').corr().abs()\nupper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))\n\n# For each pair with corr > threshold, drop the one with LOWER MI\nto_drop_corr = set()\nfor col in upper.columns:\n    high_corr_cols = upper.index[upper[col] > CORR_THRESH].tolist()\n    for hc in high_corr_cols:\n        # Keep the one with higher MI\n        mi_col = mi_df.set_index('feature').loc[col, 'mi']\n        mi_hc = mi_df.set_index('feature').loc[hc, 'mi']\n        drop_this = hc if mi_col >= mi_hc else col\n        to_drop_corr.add(drop_this)\n\nsel_corr = [c for c in sel_mi if c not in to_drop_corr]\n\nprint(f'Correlation threshold: {CORR_THRESH}')\nprint(f'Features kept: {len(sel_corr)}')\nif to_drop_corr:\n    print(f'Features dropped (redundant): {sorted(to_drop_corr)}')\n\n# Visualization: Correlation heatmap of remaining features\nfig, ax = plt.subplots(figsize=(12, 10))\ncorr_final = X_df[sel_corr].astype('float32').corr()\nmask = np.triu(np.ones_like(corr_final), k=1)\nsns.heatmap(corr_final, mask=mask, cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, annot=False, square=True,\n            linewidths=0.5, ax=ax, cbar_kws={'shrink': 0.8})\nax.set_title(f'Step 2: Correlation Matrix ({len(sel_corr)} features after pruning)',\n             fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Step 3: RFE with LightGBM\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('Step 3: Recursive Feature Elimination (RFE)')\nprint('='*50)\n\nX_rfe = X_df[sel_corr]\n\n# Target: select 15-20 features (lean but powerful)\nN_TARGET = min(20, len(sel_corr))\n\nmodel = lgb.LGBMClassifier(\n    n_estimators=200,\n    max_depth=6,\n    learning_rate=0.05,\n    random_state=SEED,\n    class_weight='balanced',\n    verbose=-1,\n    n_jobs=-1,\n)\n\nrfe = RFE(\n    estimator=model,\n    n_features_to_select=N_TARGET,\n    step=2,\n    verbose=0,\n)\n\nrfe.fit(X_rfe, y_all)\n\nFINAL_FEATURES = [f for f, s in zip(sel_corr, rfe.support_) if s]\nrfe_rankings = pd.DataFrame({\n    'feature': sel_corr,\n    'rfe_rank': rfe.ranking_,\n    'selected': rfe.support_,\n}).sort_values('rfe_rank')\n\nprint(f'Target features: {N_TARGET}')\nprint(f'Final features: {len(FINAL_FEATURES)}')\n\n# Show RFE rankings\nprint('\\nRFE Rankings (1 = selected):')\nprint(rfe_rankings.to_string(index=False))\n\n# ═══════════════════════════════════════════════════════\n#  Step 4: Summary & Comparison\n# ═══════════════════════════════════════════════════════\n\nprint('\\n' + '='*50)\nprint('FINAL FEATURE SET')\nprint('='*50)\n\n# Combine MI and RFE info\nfinal_info = (\n    mi_df[mi_df['feature'].isin(FINAL_FEATURES)]\n    .merge(rfe_rankings[rfe_rankings['selected']], on='feature')\n    .sort_values('mi', ascending=False)\n)\n\nprint(f'\\n{len(FINAL_FEATURES)} features selected:')\nprint(final_info[['feature', 'mi', 'rfe_rank']].to_string(index=False))\n\n# Categorize final features\ndynamic_final = [f for f in FINAL_FEATURES if any(f.startswith(p) for p in \n    ['sales_w', 'velocity', 'acceleration', 'peak_', 'early_', 'time_to', 'post_peak', 'pre_peak', 'is_early'])]\nbehavioral_final = [f for f in FINAL_FEATURES if f in \n    ['gini_coeff', 'total_buyers', 'repeat_purchase_rate', 'avg_purchases_per_buyer', 'interpurchase_cv', 'mean_gap']]\nstatic_final = [f for f in FINAL_FEATURES if f in \n    ['product_type_no', 'department_no', 'graphical_appearance_no', 'colour_group_code',\n     'perceived_colour_value_id', 'perceived_colour_master_id', 'index_code', 'index_group_no',\n     'section_no', 'garment_group_no', 'dept_target_enc', 'ptype_target_enc',\n     'is_allover_pattern', 'is_bright_colour', 'allover_x_bright']]\noperational_final = [f for f in FINAL_FEATURES if f in \n    ['n_stockout_weeks', 'max_stockout_tier', 'early_stockout_count']]\n\nprint(f'\\nFeature breakdown:')\nprint(f'  Dynamic (early signals)     : {len(dynamic_final)}  {dynamic_final}')\nprint(f'  Behavioral (buyer)          : {len(behavioral_final)}  {behavioral_final}')\nprint(f'  Static (article metadata)   : {len(static_final)}  {static_final}')\nprint(f'  Operational (stockout)       : {len(operational_final)}  {operational_final}')\n\n# Visualization: Final feature importance (MI)\nfig, ax = plt.subplots(figsize=(10, max(5, len(FINAL_FEATURES) * 0.35)))\nfinal_sorted = final_info.sort_values('mi', ascending=True)\n\n# Color by category\ndef get_color(f):\n    if f in dynamic_final: return 'steelblue'\n    if f in behavioral_final: return 'coral'\n    if f in static_final: return 'mediumseagreen'\n    if f in operational_final: return 'gold'\n    return 'gray'\n\ncolors = [get_color(f) for f in final_sorted['feature']]\nax.barh(range(len(final_sorted)), final_sorted['mi'].values, color=colors)\nax.set_yticks(range(len(final_sorted)))\nax.set_yticklabels(final_sorted['feature'].values, fontsize=9)\nax.set_xlabel('Mutual Information')\nax.set_title(f'Final {len(FINAL_FEATURES)} Features by MI Score', fontweight='bold')\n\n# Legend\nfrom matplotlib.patches import Patch\nlegend_items = [\n    Patch(color='steelblue', label='Dynamic'),\n    Patch(color='coral', label='Behavioral'),\n    Patch(color='mediumseagreen', label='Static'),\n    Patch(color='gold', label='Operational'),\n]\nax.legend(handles=legend_items, loc='lower right')\nplt.tight_layout()\nplt.show()\n\n# ═══════════════════════════════════════════════════════\n#  Prepare final X, y\n# ═══════════════════════════════════════════════════════\n\nX = X_df[FINAL_FEATURES].copy()\ny = y_all\n\nprint(f'\\n✅ Phase 4 complete')\nprint(f'   X shape: {X.shape}')\nprint(f'   y shape: {y.shape}')\nprint(f'   Fad rate: {y.mean():.2%}')\nprint(f'   Ready for Phase 5 (Model Development)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T20:50:54.529321Z","iopub.execute_input":"2026-04-27T20:50:54.529646Z","iopub.status.idle":"2026-04-27T20:51:27.243643Z","shell.execute_reply.started":"2026-04-27T20:50:54.529619Z","shell.execute_reply":"2026-04-27T20:51:27.243039Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 20b: Export (FIXED) ──\nimport os, json\n\nOUTPUT_DIR = '/kaggle/working/processed_data'\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\n# ✅ สร้าง export_df จาก X_df (ที่มี article_id) + y\n# X_df สร้างจาก features ตอน Phase 4 แต่ drop article_id ไปแล้ว\n# ต้องดึง article_id กลับมาจาก features_spark\n\n# Reconstruct: ดึง article_id + fad_label + final features\nexport_df = features_spark.select(\n    'article_id', 'fad_label', *FINAL_FEATURES\n).toPandas()\n\n# 1. Save features\nexport_df.to_parquet(f'{OUTPUT_DIR}/features_final.parquet', index=False)\nexport_df.to_csv(f'{OUTPUT_DIR}/features_final.csv', index=False)\nprint(f'✅ Feature matrix: {export_df.shape}')\n\n# 2. Save labels\nlabels.to_parquet(f'{OUTPUT_DIR}/labels_hybrid.parquet', index=False)\nlabels.to_csv(f'{OUTPUT_DIR}/labels_hybrid.csv', index=False)\nprint(f'✅ Labels: {labels.shape}')\n\n# 3. Save pivot (for LSTM)\npivot.to_parquet(f'{OUTPUT_DIR}/pivot_26weeks.parquet')\nprint(f'✅ Pivot: {pivot.shape}')\n\n# 4. Save constants\nconstants = {\n    'FINAL_FEATURES': FINAL_FEATURES,\n    'EARLY_WINDOW': EARLY_WINDOW,\n    'MAX_WEEKS': MAX_WEEKS,\n    'KSHAPE_K': BEST_K,\n    'FAD_CLUSTERS': FAD_CLUSTERS,\n    'FINAL_FAD_RATE': float(labels['fad_label'].mean()),\n    'n_articles': len(export_df),\n    'n_features': len(FINAL_FEATURES),\n}\n\nwith open(f'{OUTPUT_DIR}/pipeline_constants.json', 'w') as f:\n    json.dump(constants, f, indent=2)\nprint(f'✅ Constants: {len(constants)} items')\n\n# 5. Save MI scores\nmi_df.to_csv(f'{OUTPUT_DIR}/mi_scores.csv', index=False)\nprint(f'✅ MI scores')\n\n# Summary\nprint(f'\\n{\"=\"*50}')\nfor f_name in sorted(os.listdir(OUTPUT_DIR)):\n    size = os.path.getsize(f'{OUTPUT_DIR}/{f_name}') / 1024 / 1024\n    print(f'  {f_name:40s} {size:.1f} MB')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-16T18:36:09.754599Z","iopub.execute_input":"2026-04-16T18:36:09.755385Z","iopub.status.idle":"2026-04-16T18:36:17.326292Z","shell.execute_reply.started":"2026-04-16T18:36:09.755352Z","shell.execute_reply":"2026-04-16T18:36:17.325438Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 5 — Stacked Ensemble","metadata":{}},{"cell_type":"code","source":"# ── Cell 21: Temporal Train/Test Split ──────────────────────\nlaunch_map = (\n    weekly.groupBy('article_id')\n    .agg(F.min('week_start').alias('launch_week'))\n    .toPandas()\n)\nfeatures = features.merge(launch_map, on='article_id', how='left')\n\nCUTOFF = pd.Timestamp('2020-06-01')  # adjust from EDA temporal imbalance\ntrain_mask = features['launch_week'] < CUTOFF\ntest_mask  = features['launch_week'] >= CUTOFF\n\nX_train, y_train = X[train_mask], y[train_mask]\nX_test,  y_test  = X[test_mask],  y[test_mask]\n\nprint(f'Train: {X_train.shape[0]:,} | Fad rate: {y_train.mean():.2%}')\nprint(f'Test : {X_test.shape[0]:,}  | Fad rate: {y_test.mean():.2%}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-16T18:10:55.710446Z","iopub.execute_input":"2026-04-16T18:10:55.711693Z","iopub.status.idle":"2026-04-16T18:10:59.642718Z","shell.execute_reply.started":"2026-04-16T18:10:55.711659Z","shell.execute_reply":"2026-04-16T18:10:59.641956Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 22: SMOTE ─────────────────────────────────────────\nsmote = SMOTE(random_state=SEED, k_neighbors=5)\nX_train_res, y_train_res = smote.fit_resample(X_train, y_train)\nprint(f'Before SMOTE: {np.bincount(y_train)}')\nprint(f'After  SMOTE: {np.bincount(y_train_res)}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T16:46:30.958717Z","iopub.execute_input":"2026-04-15T16:46:30.959664Z","iopub.status.idle":"2026-04-15T16:46:31.278860Z","shell.execute_reply.started":"2026-04-15T16:46:30.959631Z","shell.execute_reply":"2026-04-15T16:46:31.278192Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 23: LightGBM (5-fold CV) ──────────────────────────\nimport numpy as np\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import f1_score\nimport lightgbm as lgb\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)\n\n# แปลงเป็น numpy array ก่อน loop — เร็วกว่า DataFrame indexing มาก\n# สำหรับ big data ลด overhead ของ pandas ได้ชัดเจน\nif hasattr(X_train_res, 'values'):\n    X_tr_arr = X_train_res.values  # DataFrame → numpy\nelse:\n    X_tr_arr = X_train_res         # ถ้าเป็น numpy อยู่แล้ว\n\nif hasattr(X_test, 'values'):\n    X_te_arr = X_test.values\nelse:\n    X_te_arr = X_test\n\nif hasattr(y_train_res, 'values'):\n    y_tr_arr = y_train_res.values\nelse:\n    y_tr_arr = y_train_res\n\noof_lgb  = np.zeros(len(X_tr_arr))\npred_lgb = np.zeros(len(X_te_arr))\n\nlgb_params = dict(\n    n_estimators=500,\n    learning_rate=0.05,\n    num_leaves=31,\n    subsample=0.8,\n    colsample_bytree=0.8,\n    min_child_samples=20,\n    class_weight='balanced',\n    random_state=SEED,\n    verbose=-1,\n    # big data optimizations\n    n_jobs=-1,           # ใช้ทุก CPU core\n    device='cpu',        # เปลี่ยนเป็น 'gpu' ถ้ามี GPU\n)\n\nfor fold, (tr, val) in enumerate(cv.split(X_tr_arr, y_tr_arr)):\n    m = lgb.LGBMClassifier(**lgb_params)\n    \n    # ใช้ numpy indexing — ไม่มี KeyError + เร็วกว่า\n    m.fit(\n        X_tr_arr[tr], y_tr_arr[tr],\n        eval_set=[(X_tr_arr[val], y_tr_arr[val])],\n        callbacks=[\n            lgb.early_stopping(50, verbose=False),\n            lgb.log_evaluation(0)\n        ]\n    )\n    \n    oof_lgb[val]  = m.predict_proba(X_tr_arr[val])[:, 1]\n    pred_lgb     += m.predict_proba(X_te_arr)[:, 1] / 5\n    print(f'  Fold {fold+1} F1: {f1_score(y_tr_arr[val], oof_lgb[val] > 0.5):.4f}')\n\nlgb_model = m\nprint(f'LightGBM OOF F1: {f1_score(y_tr_arr, oof_lgb > 0.5):.4f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T16:46:36.828362Z","iopub.execute_input":"2026-04-15T16:46:36.829039Z","iopub.status.idle":"2026-04-15T16:46:59.144787Z","shell.execute_reply.started":"2026-04-15T16:46:36.829007Z","shell.execute_reply":"2026-04-15T16:46:59.144264Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 24: XGBoost (5-fold CV) ───────────────────────────\n\n# แปลงเป็น numpy array ก่อน — แก้ KeyError + เร็วกว่าสำหรับ big data\nif hasattr(X_train_res, 'values'):\n    X_tr_arr = X_train_res.values\nelse:\n    X_tr_arr = X_train_res\n\nif hasattr(X_test, 'values'):\n    X_te_arr = X_test.values\nelse:\n    X_te_arr = X_test\n\nif hasattr(y_train_res, 'values'):\n    y_tr_arr = y_train_res.values\nelse:\n    y_tr_arr = y_train_res\n\n# คำนวณ scale_pos_weight จาก numpy array\nspw = (y_tr_arr == 0).sum() / (y_tr_arr == 1).sum()\n\noof_xgb  = np.zeros(len(X_tr_arr))\npred_xgb = np.zeros(len(X_te_arr))\n\nxgb_params = dict(\n    n_estimators=500,\n    learning_rate=0.05,\n    max_depth=6,\n    subsample=0.8,\n    colsample_bytree=0.8,\n    scale_pos_weight=spw,\n    random_state=SEED,\n    eval_metric='logloss',\n    early_stopping_rounds=50,  # ✅ ย้ายมาไว้ตรงนี้แทน\n    tree_method='hist',\n    n_jobs=-1,\n    device='cpu',\n)\n\nfor fold, (tr, val) in enumerate(cv.split(X_tr_arr, y_tr_arr)):\n    m_x = xgb.XGBClassifier(**xgb_params)\n\n    m_x.fit(\n        X_tr_arr[tr], y_tr_arr[tr],\n        eval_set=[(X_tr_arr[val], y_tr_arr[val])],\n        verbose=False  # ลบ early_stopping_rounds ออกจากตรงนี้\n    )\n\n    oof_xgb[val]  = m_x.predict_proba(X_tr_arr[val])[:, 1]\n    pred_xgb     += m_x.predict_proba(X_te_arr)[:, 1] / 5\n    print(f'  Fold {fold+1} F1: {f1_score(y_tr_arr[val], oof_xgb[val] > 0.5):.4f}')\n\nxgb_model = m_x\nprint(f'XGBoost OOF F1: {f1_score(y_tr_arr, oof_xgb > 0.5):.4f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T16:58:02.554159Z","iopub.execute_input":"2026-04-15T16:58:02.554846Z","iopub.status.idle":"2026-04-15T16:58:24.780918Z","shell.execute_reply.started":"2026-04-15T16:58:02.554814Z","shell.execute_reply":"2026-04-15T16:58:24.780316Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 25: LSTM (Time-Series Sequences) ──────────────────\ntrain_arts = features[train_mask]['article_id'].values\ntest_arts  = features[test_mask]['article_id'].values\n\nX_lstm_train = pivot.reindex(train_arts).fillna(0).values.reshape(-1, MAX_WEEKS, 1)\nX_lstm_test  = pivot.reindex(test_arts).fillna(0).values.reshape(-1, MAX_WEEKS, 1)\n\n# ACF-informed weight: reduce LSTM contribution if Fad series have short memory\nLSTM_WEIGHT = 0.2\n\ncw = {0: 1.0, 1: (y_train==0).sum()/(y_train.sum()+1e-6)}\nlstm_model = Sequential([\n    LSTM(64, return_sequences=True, input_shape=(MAX_WEEKS, 1)),\n    Dropout(0.3),\n    LSTM(32),\n    Dropout(0.3),\n    Dense(16, activation='relu'),\n    Dense(1, activation='sigmoid')\n])\nlstm_model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['AUC'])\nes = EarlyStopping(monitor='val_loss', patience=10, restore_best_weights=True)\nlstm_model.fit(X_lstm_train, y_train, validation_split=0.2,\n               epochs=50, batch_size=256, class_weight=cw,\n               callbacks=[es], verbose=1)\n\npred_lstm = lstm_model.predict(X_lstm_test).ravel()\nprint(f'LSTM done. Predictions shape: {pred_lstm.shape}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T16:58:29.435675Z","iopub.execute_input":"2026-04-15T16:58:29.436438Z","iopub.status.idle":"2026-04-15T17:00:05.893852Z","shell.execute_reply.started":"2026-04-15T16:58:29.436407Z","shell.execute_reply":"2026-04-15T17:00:05.893260Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 26: Meta-Learner + Final Blend ────────────────────\nmeta_train = np.column_stack([oof_lgb, oof_xgb])\nmeta_test  = np.column_stack([pred_lgb, pred_xgb])\n\nmeta_lr = LogisticRegression(C=1.0, random_state=SEED)\nmeta_lr.fit(meta_train, y_train_res)\n\nmeta_pred = meta_lr.predict_proba(meta_test)[:, 1]\nfinal_pred = (1 - LSTM_WEIGHT) * meta_pred + LSTM_WEIGHT * pred_lstm\n\n# Calibration\ncal_lgb = CalibratedClassifierCV(lgb_model, method='sigmoid', cv='prefit')\ncal_lgb.fit(X_train_res, y_train_res)\n\nprint('Meta-learner weights:', meta_lr.coef_)\nprint('Ensemble ready.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:01:00.229095Z","iopub.execute_input":"2026-04-15T17:01:00.229515Z","iopub.status.idle":"2026-04-15T17:01:02.779748Z","shell.execute_reply.started":"2026-04-15T17:01:00.229476Z","shell.execute_reply":"2026-04-15T17:01:02.779160Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 6 — Cost-Sensitive Threshold Optimization","metadata":{}},{"cell_type":"code","source":"# ── Cell 27: Threshold Sweep ────────────────────────────────\nCOST_FN, COST_FP = 10, 1  # FN (miss Fad → Dead Stock) costs 10x more than FP\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    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\nOPTIMAL_T = thresholds[np.argmin(costs)]\nprint(f'Optimal threshold: {OPTIMAL_T:.2f}')\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\naxes[0].plot(thresholds, costs, 'coral', lw=2)\naxes[0].axvline(OPTIMAL_T, color='navy', linestyle='--', label=f'Optimal={OPTIMAL_T:.2f}')\naxes[0].set_title(f'Business Cost (FN={COST_FN}x, FP={COST_FP}x)', fontweight='bold')\naxes[0].legend()\n\naxes[1].plot(thresholds, f1s,   label='F1',        color='steelblue')\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_title('Metrics vs Threshold', fontweight='bold')\naxes[1].legend()\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:01:09.049496Z","iopub.execute_input":"2026-04-15T17:01:09.050233Z","iopub.status.idle":"2026-04-15T17:01:09.868070Z","shell.execute_reply.started":"2026-04-15T17:01:09.050203Z","shell.execute_reply":"2026-04-15T17:01:09.867301Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 28: Evaluation + Business ROI ─────────────────────\ny_pred = (final_pred >= OPTIMAL_T).astype(int)\ncm = confusion_matrix(y_test, y_pred)\ntn, fp, fn, tp = cm.ravel()\npr_auc = average_precision_score(y_test, final_pred)\n\nprint('=== Model Evaluation ===')\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_score(y_test, final_pred):.4f}')\n\n# Business ROI\nCOST_PER_FAD_UNIT   = 500\nUNITS_PER_FAD       = 200\nMARKDOWN_RATE       = 0.50\nOPP_LOSS_PCT        = 0.20\nFAD_REVENUE         = 1000\n\navoided  = tp * UNITS_PER_FAD * COST_PER_FAD_UNIT * MARKDOWN_RATE\nopp_lost = fp * FAD_REVENUE * OPP_LOSS_PCT\nds_loss  = fn * UNITS_PER_FAD * COST_PER_FAD_UNIT * MARKDOWN_RATE\nnet      = avoided - opp_lost - ds_loss\n\nprint(f'\\n=== Business ROI ===')\nprint(f'Markdown avoided : {avoided:>12,.0f} THB')\nprint(f'Opportunity lost : {opp_lost:>12,.0f} THB')\nprint(f'Dead stock loss  : {ds_loss:>12,.0f} THB')\nprint(f'Net Benefit      : {net:>12,.0f} THB  {\"✓ POSITIVE\" if net>0 else \"✗ NEGATIVE\"}')\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 5))\nsns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[0],\n            xticklabels=['Non-Fad','Fad'], yticklabels=['Non-Fad','Fad'])\naxes[0].set_title('Confusion Matrix', fontweight='bold')\n\nprec_c, rec_c, _ = precision_recall_curve(y_test, final_pred)\naxes[1].plot(rec_c, prec_c, color='steelblue', lw=2, label=f'PR-AUC={pr_auc:.3f}')\naxes[1].axhline(y_test.mean(), color='red', linestyle='--', label='Random')\naxes[1].set_xlabel('Recall'); axes[1].set_ylabel('Precision')\naxes[1].set_title('Precision-Recall Curve', fontweight='bold')\naxes[1].legend()\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:01:14.408621Z","iopub.execute_input":"2026-04-15T17:01:14.408894Z","iopub.status.idle":"2026-04-15T17:01:14.712981Z","shell.execute_reply.started":"2026-04-15T17:01:14.408872Z","shell.execute_reply":"2026-04-15T17:01:14.712186Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 7 — SHAP + Counterfactual Explainability","metadata":{}},{"cell_type":"code","source":"# ── Cell 29: Global SHAP ────────────────────────────────────\nexplainer   = shap.TreeExplainer(lgb_model)\nshap_values = explainer.shap_values(X_test)\nsv = shap_values[1] if isinstance(shap_values, list) else shap_values\n\nplt.figure(figsize=(10, 6))\nshap.summary_plot(sv, X_test, feature_names=FINAL_FEATURES, max_display=20, show=False)\nplt.title('SHAP Summary — Global Feature Importance', fontweight='bold')\nplt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:01:33.049952Z","iopub.execute_input":"2026-04-15T17:01:33.050485Z","iopub.status.idle":"2026-04-15T17:01:43.346351Z","shell.execute_reply.started":"2026-04-15T17:01:33.050454Z","shell.execute_reply":"2026-04-15T17:01:43.345615Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 30: Local SHAP + Counterfactual ────────────────────\n\n# แปลงเป็น numpy ก่อน ป้องกัน KeyError ทุกจุด\nX_test_arr = X_test.values if hasattr(X_test, 'values') else X_test\ny_test_arr = y_test.values if hasattr(y_test, 'values') else y_test\ny_pred_arr = y_pred.values if hasattr(y_pred, 'values') else y_pred\n\ndef explain_and_counterfactual(idx):\n    pred_p = final_pred[idx]\n    actual = y_test_arr[idx]\n    print(f\"Predicted Fad prob : {pred_p:.3f}  ({'FAD' if pred_p >= OPTIMAL_T else 'Non-Fad'})\")\n    print(f\"Actual label       : {'FAD' if actual == 1 else 'Non-Fad'}\")\n\n    # SHAP waterfall\n    plt.figure(figsize=(10, 4))\n    ev = explainer.expected_value\n    base_val = ev[1] if isinstance(ev, (list, np.ndarray)) else ev\n\n    shap.waterfall_plot(\n        shap.Explanation(\n            values=sv[idx],\n            base_values=base_val,\n            data=X_test_arr[idx],        # ✅ numpy indexing\n            feature_names=FINAL_FEATURES\n        ),\n        max_display=15,\n        show=False\n    )\n    plt.title('SHAP Waterfall (Per-SKU)', fontweight='bold')\n    plt.tight_layout()\n    plt.show()\n\n    # Counterfactual\n    shap_s = pd.Series(sv[idx], index=FINAL_FEATURES)\n    top5 = shap_s.nlargest(5)\n    print('\\nCounterfactual — features that drive FAD prediction:')\n    for feat, val in top5.items():\n        feat_val = X_test_arr[idx][FINAL_FEATURES.index(feat)]  # ✅ numpy indexing\n        print(f'  {feat:35s}: current={feat_val:.3f}  SHAP=+{val:.4f}')\n        print(f'  {\"\":35s}  → reduce to {feat_val * 0.5:.3f} to lower Fad probability')\n\n# ใช้ y_test_arr และ y_pred_arr แทน\ntp_idx = np.where((y_test_arr == 1) & (y_pred_arr == 1))[0]\nfn_idx = np.where((y_test_arr == 1) & (y_pred_arr == 0))[0]\n\nif len(tp_idx) > 0:\n    print('\\n=== True Positive (correctly caught Fad) ===')\n    explain_and_counterfactual(tp_idx[0])\n\nif len(fn_idx) > 0:\n    print('\\n=== False Negative (missed Fad → Dead Stock risk) ===')\n    explain_and_counterfactual(fn_idx[0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:02:02.933945Z","iopub.execute_input":"2026-04-15T17:02:02.934318Z","iopub.status.idle":"2026-04-15T17:02:03.932792Z","shell.execute_reply.started":"2026-04-15T17:02:02.934288Z","shell.execute_reply":"2026-04-15T17:02:03.932181Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 31: Final Summary + Stop Spark ────────────────────\nprint('=' * 55)\nprint('   FINAL PIPELINE SUMMARY')\nprint('=' * 55)\nmetrics = {\n    'F1-Score (Fad)'     : f1_score(y_test, y_pred),\n    'Precision (Fad)'    : precision_score(y_test, y_pred),\n    'Recall (Fad)'       : recall_score(y_test, y_pred),\n    'PR-AUC'             : pr_auc,\n    'ROC-AUC'            : roc_auc_score(y_test, final_pred),\n    'Optimal Threshold'  : OPTIMAL_T,\n    'TP / FP / FN / TN'  : f'{tp} / {fp} / {fn} / {tn}',\n    'Net Business Benefit': f'{net:,.0f} THB'\n}\nfor k, v in metrics.items():\n    print(f'  {k:25s}: {v:.4f}' if isinstance(v, float) else f'  {k:25s}: {v}')\nprint('=' * 55)\n\nspark.stop()\nprint('Spark session stopped.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-15T17:02:14.837686Z","iopub.execute_input":"2026-04-15T17:02:14.838508Z","iopub.status.idle":"2026-04-15T17:02:16.230667Z","shell.execute_reply.started":"2026-04-15T17:02:14.838476Z","shell.execute_reply":"2026-04-15T17:02:16.229921Z"}},"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":{}}]}