{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":31254,"databundleVersionId":3103714,"isSourceIdPinned":false}],"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-11T17:10:46.237469Z","iopub.execute_input":"2026-04-11T17:10:46.237913Z","iopub.status.idle":"2026-04-11T17:10:46.242227Z","shell.execute_reply.started":"2026-04-11T17:10:46.237884Z","shell.execute_reply":"2026-04-11T17:10:46.241526Z"}},"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-11-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-11T17:01:11.432656Z","iopub.execute_input":"2026-04-11T17:01:11.433186Z","iopub.status.idle":"2026-04-11T17:01:15.098804Z","shell.execute_reply.started":"2026-04-11T17:01:11.433139Z","shell.execute_reply":"2026-04-11T17:01:15.097784Z"}},"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-11T17:11:00.711406Z","iopub.execute_input":"2026-04-11T17:11:00.711852Z","iopub.status.idle":"2026-04-11T17:11:13.280327Z","shell.execute_reply.started":"2026-04-11T17:11:00.711822Z","shell.execute_reply":"2026-04-11T17:11:13.279449Z"}},"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-11T17:11:31.733162Z","iopub.execute_input":"2026-04-11T17:11:31.733612Z","iopub.status.idle":"2026-04-11T17:12:43.54307Z","shell.execute_reply.started":"2026-04-11T17:11:31.733578Z","shell.execute_reply":"2026-04-11T17:12:43.542167Z"}},"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-11T17:13:00.041482Z","iopub.execute_input":"2026-04-11T17:13:00.042378Z","iopub.status.idle":"2026-04-11T17:13:09.029865Z","shell.execute_reply.started":"2026-04-11T17:13:00.042345Z","shell.execute_reply":"2026-04-11T17:13:09.028344Z"}},"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    total = df.count()\n    miss_exprs = [\n        F.sum(F.col(c).isNull().cast('int')).alias(c)\n        for c in df.columns\n    ]\n    miss_counts = df.agg(*miss_exprs).collect()[0].asDict()\n    \n    rows = [\n        {'column': col, 'missing': cnt, 'pct': round(cnt/total*100, 2)}\n        for col, cnt in miss_counts.items() if cnt > 0\n    ]\n    \n    print(f'\\n=== Missing Values: {name} ===')\n    if len(rows) == 0:\n        print('  No missing values')\n        return pd.DataFrame()\n    \n    report = pd.DataFrame(rows).sort_values('pct', ascending=False)\n    print(report.to_string(index=False))\n    return report\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-11T17:15:18.583862Z","iopub.execute_input":"2026-04-11T17:15:18.584338Z","iopub.status.idle":"2026-04-11T17:16:01.83676Z","shell.execute_reply.started":"2026-04-11T17:15:18.584307Z","shell.execute_reply":"2026-04-11T17:16:01.835943Z"}},"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-11T17:16:09.865842Z","iopub.execute_input":"2026-04-11T17:16:09.86662Z","iopub.status.idle":"2026-04-11T17:16:11.644585Z","shell.execute_reply.started":"2026-04-11T17:16:09.86659Z","shell.execute_reply":"2026-04-11T17:16:11.643943Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 7: Sales Distribution per Article (Spark) ──────────\n# Aggregate on full 31.7M rows — parallel in Spark\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\n# Describe distribution\nsales_per_art.select('total_sales').describe().show()\n\n# Bring to Pandas for plotting (only summary stats, not full data)\nsales_pd = sales_per_art.toPandas()\nsales_pd['rank'] = range(1, len(sales_pd) + 1)\n\n# Log-log Power Law test\nfrom scipy import stats\nlog_rank  = np.log10(sales_pd['rank'])\nlog_sales = np.log10(sales_pd['total_sales'])\nslope, intercept, r, p, _ = stats.linregress(log_rank, log_sales)\n\nfig, ax = plt.subplots(figsize=(8, 5))\nax.scatter(log_rank, log_sales, alpha=0.2, s=2, color='steelblue')\nx_line = np.linspace(log_rank.min(), log_rank.max(), 100)\nax.plot(x_line, slope*x_line + intercept, 'r-', lw=2,\n        label=f'slope={slope:.2f}  R²={r**2:.3f}')\nax.set_xlabel('log10(Rank)'); ax.set_ylabel('log10(Sales)')\nax.set_title('Power Law Test — Sales per Article (full data)', fontweight='bold')\nax.legend()\nplt.tight_layout(); plt.show()\nprint(f'Power Law slope={slope:.3f}. Straight line with R²>0.9 → Power Law confirmed')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:16:15.526723Z","iopub.execute_input":"2026-04-11T17:16:15.527533Z","iopub.status.idle":"2026-04-11T17:16:38.853176Z","shell.execute_reply.started":"2026-04-11T17:16:15.527503Z","shell.execute_reply":"2026-04-11T17:16:38.852456Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 1 — Data Preparation with PySpark","metadata":{}},{"cell_type":"code","source":"# ── Cell 8: Weekly Aggregation (PySpark — parallel) ─────────\n# weekofyear + year to avoid year-boundary bugs\ntxn_w = (\n    txn\n    .withColumn('year',    F.year('t_dat'))\n    .withColumn('week_of_year', F.weekofyear('t_dat'))\n    # ISO week start (Monday)\n    .withColumn('week_start',\n        F.date_trunc('week', F.col('t_dat')))\n)\n\nweekly = (\n    txn_w.groupBy('article_id', 'week_start')\n    .agg(\n        F.count('customer_id')         .alias('weekly_sales'),\n        F.countDistinct('customer_id') .alias('unique_buyers'),\n        F.countDistinct('sales_channel_id').alias('n_channels')\n    )\n    .orderBy('article_id', 'week_start')\n)\n\n# Add launch_week and weeks_since_launch\nlaunch = weekly.groupBy('article_id').agg(F.min('week_start').alias('launch_week'))\nweekly = (\n    weekly\n    .join(launch, on='article_id', how='left')\n    .withColumn('weeks_since_launch',\n        (F.datediff('week_start', 'launch_week') / 7).cast('int'))\n)\n\n# Cache for repeated use\nweekly.cache()\nprint(f'Weekly records: {weekly.count():,}')\nweekly.show(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:18:38.166873Z","iopub.execute_input":"2026-04-11T17:18:38.167246Z","iopub.status.idle":"2026-04-11T17:20:41.825382Z","shell.execute_reply.started":"2026-04-11T17:18:38.16722Z","shell.execute_reply":"2026-04-11T17:20:41.824655Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 9: Stockout Detection (PySpark Window Functions) ───\nw_art = Window.partitionBy('article_id').orderBy('week_start')\n\nweekly = (\n    weekly\n    .withColumn('sales_lag1',  F.lag('weekly_sales', 1).over(w_art))\n    .withColumn('sales_lead1', F.lead('weekly_sales', 1).over(w_art))\n    .withColumn('pct_drop',\n        (F.col('weekly_sales') - F.col('sales_lag1')) /\n        (F.col('sales_lag1') + 1e-6))\n    .withColumn('pct_recover',\n        (F.col('sales_lead1') - F.col('weekly_sales')) /\n        (F.col('weekly_sales') + 1e-6))\n    .withColumn('is_stockout',\n        ((F.col('pct_drop') < -0.80) &\n         (F.col('pct_recover') > 0.50)).cast('int'))\n    # Interpolate stockout weeks\n    .withColumn('weekly_sales_clean',\n        F.when(F.col('is_stockout') == 1,\n               (F.col('sales_lag1') + F.col('sales_lead1')) / 2)\n         .otherwise(F.col('weekly_sales')))\n)\n\nn_stockout = weekly.filter(F.col('is_stockout') == 1).count()\nprint(f'Stockout weeks detected: {n_stockout:,}')\nweekly.cache()\nprint('Stockout masking done.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:20:52.819311Z","iopub.execute_input":"2026-04-11T17:20:52.820026Z","iopub.status.idle":"2026-04-11T17:20:55.606807Z","shell.execute_reply.started":"2026-04-11T17:20:52.819999Z","shell.execute_reply":"2026-04-11T17:20:55.605599Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 10: Filter & Normalize (PySpark) ──────────────────\n# Keep articles with >= 6 weeks of data\nlifespan = (\n    weekly.groupBy('article_id')\n    .agg(F.max('weeks_since_launch').alias('lifespan_weeks'))\n)\nvalid_arts = lifespan.filter(F.col('lifespan_weeks') >= 6)\nweekly = weekly.join(valid_arts.select('article_id'), on='article_id', how='inner')\n\n# Min-Max normalize per article using Window\nw_norm = Window.partitionBy('article_id')\nweekly = (\n    weekly\n    .withColumn('art_min', F.min('weekly_sales_clean').over(w_norm))\n    .withColumn('art_max', F.max('weekly_sales_clean').over(w_norm))\n    .withColumn('norm_sales',\n        (F.col('weekly_sales_clean') - F.col('art_min')) /\n        (F.col('art_max') - F.col('art_min') + 1e-8))\n)\n\n# Sales quartile for stratified sampling\ntotal_sales_df = (\n    weekly.groupBy('article_id')\n    .agg(F.sum('weekly_sales_clean').alias('total_sales'))\n)\nquartile_expr = F.ntile(4).over(Window.orderBy('total_sales'))\ntotal_sales_df = total_sales_df.withColumn('sales_quartile', quartile_expr)\n\nweekly = weekly.join(total_sales_df.select('article_id','sales_quartile'), on='article_id', how='left')\nweekly.cache()\n\nn_valid = weekly.select('article_id').distinct().count()\nprint(f'Valid articles (>= 6 weeks): {n_valid:,}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:20:59.480529Z","iopub.execute_input":"2026-04-11T17:20:59.481344Z","iopub.status.idle":"2026-04-11T17:21:13.186632Z","shell.execute_reply.started":"2026-04-11T17:20:59.481314Z","shell.execute_reply":"2026-04-11T17:21:13.182711Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 11: Build Sales Pivot (Spark → Pandas) ─────────────\n# Limit to first 26 weeks per article → then pivot to wide format\nMAX_WEEKS = 26\n\npivot_spark = (\n    weekly\n    .filter(F.col('weeks_since_launch') < MAX_WEEKS)\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(MAX_WEEKS)))\n    .agg(F.first('norm_sales'))\n    .fillna(0)\n)\n\n# Rename columns\nold_cols = pivot_spark.columns\nnew_cols = ['article_id'] + [f'w{c}' for c in old_cols[1:]]\npivot_spark = pivot_spark.toDF(*new_cols)\n\n# Collect to Pandas (only ~105K articles × 26 cols = small)\npivot = pivot_spark.toPandas().set_index('article_id')\nprint(f'Sales pivot shape: {pivot.shape}  (articles × weeks)')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:21:31.650248Z","iopub.execute_input":"2026-04-11T17:21:31.650739Z","iopub.status.idle":"2026-04-11T17:21:35.64051Z","shell.execute_reply.started":"2026-04-11T17:21:31.650713Z","shell.execute_reply":"2026-04-11T17:21:35.639654Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 2 — Hybrid Labeling (PELT + k-Shape)","metadata":{}},{"cell_type":"code","source":"# ── Cell 12: PELT Change Point Detection ────────────────────\n# Runs on Pandas pivot (105K × 26 = manageable)\npelt_results = []\n\nfor art_id, row in pivot.iterrows():\n    series = row.values.astype(float)\n    # Trim trailing zeros\n    nonzero = np.nonzero(series)[0]\n    if len(nonzero) == 0:\n        pelt_results.append({'article_id': art_id, 'pelt_fad': 0,\n                              'n_breakpoints': 0, 'peak_week': 0})\n        continue\n    series = series[:nonzero[-1] + 1]\n\n    if len(series) < 5:\n        pelt_results.append({'article_id': art_id, 'pelt_fad': 0,\n                              'n_breakpoints': 0, 'peak_week': 0})\n        continue\n\n    algo = rpt.Pelt(model='rbf').fit(series)\n    bkps = algo.predict(pen=2)\n    n_bkps    = len(bkps) - 1\n    peak_week = int(np.argmax(series))\n    drop      = (series[peak_week] - series[-1]) / (series[peak_week] + 1e-8)\n\n    is_fad = (\n        (n_bkps == 1) and\n        (peak_week <= 0.4 * len(series)) and\n        (drop > 0.5)\n    )\n    pelt_results.append({'article_id': art_id, 'pelt_fad': int(is_fad),\n                          'n_breakpoints': n_bkps, 'peak_week': peak_week})\n\npelt_df = pd.DataFrame(pelt_results)\nprint('PELT done.')\nprint(pelt_df['pelt_fad'].value_counts())\nprint(f'Fad rate (PELT): {pelt_df[\"pelt_fad\"].mean():.2%}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:21:39.410328Z","iopub.execute_input":"2026-04-11T17:21:39.410587Z","iopub.status.idle":"2026-04-11T17:22:10.836018Z","shell.execute_reply.started":"2026-04-11T17:21:39.410567Z","shell.execute_reply":"2026-04-11T17:22:10.835337Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 13: k-Shape Clustering ─────────────────────────────\nfrom tslearn.clustering import KShape\nfrom tslearn.preprocessing import TimeSeriesScalerMeanVariance\n\nX_ts = pivot.values.reshape(pivot.shape[0], pivot.shape[1], 1).astype(np.float64)\nX_ts = TimeSeriesScalerMeanVariance().fit_transform(X_ts)\n\nN_CLUSTERS = 5\nks = KShape(n_clusters=N_CLUSTERS, n_init=3, random_state=SEED)\nkshape_labels = ks.fit_predict(X_ts)\n\nfig, axes = plt.subplots(1, N_CLUSTERS, figsize=(16, 3))\nfor i, ax in enumerate(axes):\n    ax.plot(ks.cluster_centers_[i].ravel(), color='steelblue', lw=2)\n    n = (kshape_labels == i).sum()\n    ax.set_title(f'Cluster {i}\\nn={n:,}', fontsize=9, fontweight='bold')\n    ax.set_xlabel('Weeks')\nplt.suptitle('k-Shape Centroids — identify Fad cluster (sharp rise & fall)', fontweight='bold')\nplt.tight_layout(); plt.show()\n\nkshape_df = pd.DataFrame({'article_id': pivot.index, 'kshape_cluster': kshape_labels})","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-11T17:22:56.855034Z","iopub.execute_input":"2026-04-11T17:22:56.855354Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 14: Hybrid Label ───────────────────────────────────\n# ⚠️  Set FAD_CLUSTER after inspecting Cell 13 centroid plots\nFAD_CLUSTER = 0   # <-- update this!\n\nlabels = (\n    pelt_df\n    .merge(kshape_df, on='article_id', how='inner')\n)\nlabels['kshape_fad'] = (labels['kshape_cluster'] == FAD_CLUSTER).astype(int)\nlabels['fad_label']  = ((labels['pelt_fad'] == 1) & (labels['kshape_fad'] == 1)).astype(int)\n\nprint('=== Hybrid Label ===')\nprint(labels['fad_label'].value_counts())\nprint(f'Final Fad rate: {labels[\"fad_label\"].mean():.2%}')","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 3 — Feature Engineering with PySpark","metadata":{}},{"cell_type":"code","source":"# ── Cell 15: Dynamic Features (PySpark — all 31.7M rows) ────\n# ✅ All aggregations run in parallel on Spark\nEARLY_WINDOW = 3  # set from Cohort EDA\n\n# Sales in early weeks\nearly_spark = (\n    weekly\n    .filter(F.col('weeks_since_launch') < EARLY_WINDOW + 2)\n    .groupBy('article_id')\n    .pivot('weeks_since_launch', list(range(EARLY_WINDOW + 2)))\n    .agg(F.first('weekly_sales_clean'))\n    .fillna(0)\n)\n# Rename\nearly_cols = early_spark.columns\nearly_spark = early_spark.toDF(*(['article_id'] + [f'sales_w{c}' for c in early_cols[1:]]))\n\n# Peak sales\npeak_spark = (\n    weekly.groupBy('article_id')\n    .agg(F.max('weekly_sales_clean').alias('peak_sales'))\n)\n\n# Combine and compute derived features\ndyn = (\n    early_spark.join(peak_spark, on='article_id', how='left')\n    .withColumn('velocity_ratio',\n        F.col('sales_w1') / (F.col('sales_w0') + 1e-6))\n    .withColumn('acceleration',\n        F.col('sales_w2') / (F.col('sales_w1') + 1e-6))\n    .withColumn('peak_to_w0_ratio',\n        F.col('peak_sales') / (F.col('sales_w0') + 1e-6))\n)\n\ndyn.cache()\nprint('Dynamic features computed.')\ndyn.show(3)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 16: Gini Coefficient (PySpark UDF) ─────────────────\nfrom pyspark.sql.types import DoubleType\n\n@F.udf(DoubleType())\ndef gini_udf(sales_list):\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\n# Collect weekly sales as list per article → compute Gini\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)\nprint('Gini computed.')\ngini_spark.show(3)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 17: Repeat Purchase Rate & Interpurchase CV (Spark) -\n# Multi-buy count per (article, customer)\nmulti_buy_spark = (\n    txn.groupBy('article_id', 'customer_id')\n    .agg(F.count('*').alias('n_purchases'))\n)\n\nbuyer_stats = (\n    multi_buy_spark.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    )\n    .withColumn('repeat_purchase_rate',\n        F.col('repeat_buyers') / (F.col('total_buyers') + 1e-6))\n)\n\n# Interpurchase time CV using PySpark Window\nw_txn = Window.partitionBy('article_id').orderBy('t_dat')\ncv_spark = (\n    txn\n    .withColumn('lag_date', F.lag('t_dat', 1).over(w_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    )\n    .withColumn('interpurchase_cv',\n        F.col('std_gap') / (F.col('mean_gap') + 1e-6))\n    .select('article_id', 'interpurchase_cv')\n)\n\nprint('Buyer behavior features done.')","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 18: Static Features + Interaction Terms (Spark) ────\nart_feat_spark = (\n    art.select(\n        'article_id', 'product_type_no', 'product_group_name',\n        'graphical_appearance_no', 'colour_group_code',\n        'perceived_colour_value_id', 'perceived_colour_master_id',\n        'department_no', 'index_code', 'index_group_no',\n        'section_no', 'garment_group_no'\n    )\n    # Domain-informed interaction from Fisher's Exact Test\n    .withColumn('is_allover_pattern',\n        (F.col('graphical_appearance_no') == 74).cast('int'))  # update code from EDA\n    .withColumn('is_bright_colour',\n        F.col('colour_group_code').isin([9, 10, 11, 73]).cast('int'))  # update codes\n    .withColumn('allover_x_bright',\n        F.col('is_allover_pattern') * F.col('is_bright_colour'))\n)\n\n# StringIndexer for product_group_name\nindexer = StringIndexer(inputCol='product_group_name', outputCol='product_group_enc',\n                        handleInvalid='keep')\nart_feat_spark = indexer.fit(art_feat_spark).transform(art_feat_spark)\nart_feat_spark = art_feat_spark.drop('product_group_name')\n\nprint('Static features done.')","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 19: Merge All Features (Spark joins) ───────────────\nfeatures_spark = (\n    spark.createDataFrame(labels[['article_id','fad_label','peak_week']])\n    .join(dyn,           on='article_id', how='left')\n    .join(gini_spark,    on='article_id', how='left')\n    .join(buyer_stats.select('article_id','total_buyers','repeat_purchase_rate'),\n                         on='article_id', how='left')\n    .join(cv_spark,      on='article_id', how='left')\n    .join(art_feat_spark,on='article_id', how='left')\n    .join(total_sales_df.select('article_id','sales_quartile'),\n                         on='article_id', how='left')\n    .fillna(0)\n)\n\nfeatures_spark.cache()\nprint(f'Feature rows: {features_spark.count():,}')\nprint(f'Columns: {len(features_spark.columns)}')\n\n# Bring to Pandas for sklearn/LightGBM (feature matrix is now small ~105K × 30)\nfeatures = features_spark.toPandas()\nprint(f'Pandas shape: {features.shape}')","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## Phase 4 — Feature Selection","metadata":{}},{"cell_type":"code","source":"# ── Cell 20: Mutual Information + Correlation Pruning + RFE ─\nfrom sklearn.feature_selection import RFE\n\nDROP = ['article_id','fad_label','peak_week','sales_quartile']\nfeature_cols = [c for c in features.columns if c not in DROP]\n\nX_all = features[feature_cols].values\ny_all = features['fad_label'].values\n\n# Step 1: MI\nmi = mutual_info_classif(X_all, y_all, random_state=SEED)\nmi_df = pd.DataFrame({'feature': feature_cols, 'mi': mi}).sort_values('mi', ascending=False)\nMI_THRESH = 0.001\nsel_mi = mi_df[mi_df['mi'] >= MI_THRESH]['feature'].tolist()\nprint(f'Step 1 (MI >= {MI_THRESH}): {len(sel_mi)} features')\n\n# Step 2: Correlation pruning\ncorr = features[sel_mi].corr().abs()\nupper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))\nto_drop = [c for c in upper.columns if any(upper[c] > 0.85)]\nsel_corr = [c for c in sel_mi if c not in to_drop]\nprint(f'Step 2 (corr < 0.85): {len(sel_corr)} features')\n\n# Step 3: RFE\nX_rfe = features[sel_corr].values\nrfe = RFE(lgb.LGBMClassifier(n_estimators=100, random_state=SEED, class_weight='balanced', verbose=-1),\n          n_features_to_select=min(20, len(sel_corr)), step=2)\nrfe.fit(X_rfe, y_all)\nFINAL_FEATURES = [f for f, s in zip(sel_corr, rfe.support_) if s]\nprint(f'Step 3 (RFE): {len(FINAL_FEATURES)} final features')\nprint(FINAL_FEATURES)\n\n# Plot MI scores\nfig, ax = plt.subplots(figsize=(10, 5))\nmi_df.head(25).sort_values('mi').plot(\n    kind='barh', x='feature', y='mi', ax=ax, color='steelblue', legend=False)\nax.axvline(MI_THRESH, color='red', linestyle='--', label=f'MI threshold={MI_THRESH}')\nax.set_title('Mutual Information — Top 25 Features', fontweight='bold')\nax.legend()\nplt.tight_layout(); plt.show()\n\nX = features[FINAL_FEATURES].values\ny = features['fad_label'].values","metadata":{},"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":{},"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":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 23: LightGBM (5-fold CV) ──────────────────────────\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)\noof_lgb  = np.zeros(len(X_train_res))\npred_lgb = np.zeros(len(X_test))\n\nlgb_params = dict(n_estimators=500, learning_rate=0.05, num_leaves=31,\n                  subsample=0.8, colsample_bytree=0.8, min_child_samples=20,\n                  class_weight='balanced', random_state=SEED, verbose=-1)\n\nfor fold, (tr, val) in enumerate(cv.split(X_train_res, y_train_res)):\n    m = lgb.LGBMClassifier(**lgb_params)\n    m.fit(X_train_res[tr], y_train_res[tr],\n          eval_set=[(X_train_res[val], y_train_res[val])],\n          callbacks=[lgb.early_stopping(50, verbose=False), lgb.log_evaluation(0)])\n    oof_lgb[val]  = m.predict_proba(X_train_res[val])[:, 1]\n    pred_lgb     += m.predict_proba(X_test)[:, 1] / 5\n    print(f'  Fold {fold+1} F1: {f1_score(y_train_res[val], oof_lgb[val]>0.5):.4f}')\n\nlgb_model = m\nprint(f'LightGBM OOF F1: {f1_score(y_train_res, oof_lgb>0.5):.4f}')","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 24: XGBoost (5-fold CV) ───────────────────────────\nspw = (y_train_res==0).sum() / (y_train_res==1).sum()\noof_xgb  = np.zeros(len(X_train_res))\npred_xgb = np.zeros(len(X_test))\n\nxgb_params = dict(n_estimators=500, learning_rate=0.05, max_depth=6,\n                  subsample=0.8, colsample_bytree=0.8, scale_pos_weight=spw,\n                  random_state=SEED, eval_metric='logloss', use_label_encoder=False)\n\nfor fold, (tr, val) in enumerate(cv.split(X_train_res, y_train_res)):\n    m_x = xgb.XGBClassifier(**xgb_params)\n    m_x.fit(X_train_res[tr], y_train_res[tr],\n            eval_set=[(X_train_res[val], y_train_res[val])],\n            early_stopping_rounds=50, verbose=False)\n    oof_xgb[val]  = m_x.predict_proba(X_train_res[val])[:, 1]\n    pred_xgb     += m_x.predict_proba(X_test)[:, 1] / 5\n    print(f'  Fold {fold+1} F1: {f1_score(y_train_res[val], oof_xgb[val]>0.5):.4f}')\n\nxgb_model = m_x\nprint(f'XGBoost OOF F1: {f1_score(y_train_res, oof_xgb>0.5):.4f}')","metadata":{},"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":{},"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":{},"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":{},"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":{},"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":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Cell 30: Local SHAP + Counterfactual ────────────────────\ndef explain_and_counterfactual(idx):\n    pred_p = final_pred[idx]\n    actual = y_test[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) else ev\n    shap.waterfall_plot(\n        shap.Explanation(values=sv[idx], base_values=base_val,\n                         data=X_test[idx], feature_names=FINAL_FEATURES),\n        max_display=15, show=False)\n    plt.title('SHAP Waterfall (Per-SKU)', fontweight='bold')\n    plt.tight_layout(); 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[idx][FINAL_FEATURES.index(feat)]\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\ntp_idx = np.where((y_test==1) & (y_pred==1))[0]\nfn_idx = np.where((y_test==1) & (y_pred==0))[0]\n\nif len(tp_idx)>0:\n    print('\\n=== True Positive (correctly caught Fad) ===')\n    explain_and_counterfactual(tp_idx[0])\nif len(fn_idx)>0:\n    print('\\n=== False Negative (missed Fad → Dead Stock risk) ===')\n    explain_and_counterfactual(fn_idx[0])","metadata":{},"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":{},"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":{}}]}