{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"}],"dockerImageVersionId":31154,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"\n# Pasture Biomass Prediction — **Stacked Ensemble (Offline, Kaggle No-Internet)**\n\n\n- Base models: `RandomForestRegressor`, `ExtraTreesRegressor`, `Ridge` (scaled), `GradientBoostingRegressor` (\\*MultiOutput wrapper\n","metadata":{}},{"cell_type":"code","source":"\n# =====================\n# 1) Packages and settings\n# =====================\nimport os, glob, time, math, random\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor\nfrom sklearn.linear_model import Ridge\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.ensemble import GradientBoostingRegressor\n\nimport matplotlib.pyplot as plt\n\nSEED = 42\nrandom.seed(SEED)\nnp.random.seed(SEED)\n\nIMG_MAX_SIDE = 512\nTARGETS = [\"Dry_Green_g\",\"Dry_Dead_g\",\"Dry_Clover_g\",\"GDM_g\",\"Dry_Total_g\"]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:32:07.518599Z","iopub.execute_input":"2025-11-05T02:32:07.518891Z","iopub.status.idle":"2025-11-05T02:32:07.530056Z","shell.execute_reply.started":"2025-11-05T02:32:07.518870Z","shell.execute_reply":"2025-11-05T02:32:07.529506Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 2) Data path\n# =====================\n\ndef guess_base_dir():\n    if os.path.exists('/kaggle/input/csiro-biomass/train.csv') and os.path.exists('/kaggle/input/csiro-biomass/test.csv'):\n        return '.'\n    cands = glob.glob('/kaggle/input/**/train.csv', recursive=True)\n    if cands:\n        return os.path.dirname(cands[0])\n    return '.'\n\nBASE_DIR = '/kaggle/input/csiro-biomass'\nprint('BASE_DIR =', BASE_DIR)\n\nTRAIN_CSV = os.path.join(BASE_DIR, 'train.csv')\nTEST_CSV  = os.path.join(BASE_DIR, 'test.csv')\n\nassert os.path.exists(TRAIN_CSV), f\"train.csv олдсонгүй: {TRAIN_CSV}\"\nassert os.path.exists(TEST_CSV),  f\"test.csv олдсонгүй: {TEST_CSV}\"\n\ntrain_df = pd.read_csv(TRAIN_CSV)\nprint('train_df:', train_df.shape)\nprint(train_df.head(3))\n\ntest_df = pd.read_csv(TEST_CSV)\nprint('test_df:', test_df.shape)\nprint(test_df.head(3))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:34:46.767729Z","iopub.execute_input":"2025-11-05T02:34:46.768036Z","iopub.status.idle":"2025-11-05T02:34:46.791576Z","shell.execute_reply.started":"2025-11-05T02:34:46.768014Z","shell.execute_reply":"2025-11-05T02:34:46.790758Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 3) Img feature\n# =====================\n\ndef resolve_image_path(image_path: str):\n    p = os.path.join(BASE_DIR, image_path)\n    if os.path.exists(p):\n        return p\n    # fallback variations\n    alt1 = os.path.join(BASE_DIR, 'train/')\n    if os.path.exists(alt1):\n        return alt1\n    alt2 = os.path.join(BASE_DIR, 'test/')\n    if os.path.exists(alt2):\n        return alt2\n    if os.path.exists(image_path):\n        return image_path\n    return p\n\n\ndef load_image_rgb(path: str, max_side=IMG_MAX_SIDE):\n    with Image.open(path) as im:\n        im = im.convert('RGB')\n        w,h = im.size\n        if max(w,h) > max_side:\n            scale = max_side / float(max(w,h))\n            im = im.resize((int(round(w*scale)), int(round(h*scale))), resample=Image.BILINEAR)\n        arr = np.asarray(im).astype(np.float32) / 255.0\n    return arr\n\n\ndef rgb_to_hsv_np(rgb):\n    r,g,b = rgb[...,0], rgb[...,1], rgb[...,2]\n    cmax = np.max(rgb, axis=-1)\n    cmin = np.min(rgb, axis=-1)\n    delta = cmax - cmin + 1e-12\n    h = np.zeros_like(cmax)\n    m = delta > 0\n    idx = (cmax==r) & m; h[idx] = ((g[idx]-b[idx]) / delta[idx]) % 6\n    idx = (cmax==g) & m; h[idx] = ((b[idx]-r[idx]) / delta[idx]) + 2\n    idx = (cmax==b) & m; h[idx] = ((r[idx]-g[idx]) / delta[idx]) + 4\n    h = h / 6.0\n    s = np.zeros_like(cmax); s[m] = delta[m] / (cmax[m] + 1e-12)\n    v = cmax\n    return np.stack([h,s,v], axis=-1)\n\n\ndef basic_stats(x, prefix):\n    feats = {}\n    x = x.reshape(-1)\n    x = x[np.isfinite(x)]\n    if x.size==0:\n        x = np.array([0.0])\n    feats[f'{prefix}_mean'] = float(np.mean(x))\n    feats[f'{prefix}_std']  = float(np.std(x))\n    feats[f'{prefix}_min']  = float(np.min(x))\n    for p in [1,5,25,50,75,95,99]:\n        feats[f'{prefix}_p{p:02d}'] = float(np.percentile(x, p))\n    feats[f'{prefix}_max']  = float(np.max(x))\n    return feats\n\n\ndef channel_hist(x, bins=16, prefix='r'):\n    x = x.reshape(-1)\n    x = x[np.isfinite(x)]\n    if x.size==0:\n        hist = np.zeros(bins, dtype=np.float32)\n    else:\n        hist, _ = np.histogram(x, bins=bins, range=(0.0,1.0), density=True)\n    return {f'{prefix}_hist_{i}': float(v) for i,v in enumerate(hist)}\n\n\ndef laplacian_var(gray):\n    c = gray\n    up    = np.roll(c, -1, axis=0)\n    down  = np.roll(c,  1, axis=0)\n    left  = np.roll(c,  1, axis=1)\n    right = np.roll(c, -1, axis=1)\n    lap = (-4.0*c + up + down + left + right)\n    return float(np.var(lap))\n\n\ndef gradient_stats(gray, prefix='grad'):\n    gx = np.gradient(gray, axis=1)\n    gy = np.gradient(gray, axis=0)\n    mag = np.sqrt(gx*gx + gy*gy)\n    feats = basic_stats(mag, prefix)\n    feats[f'{prefix}_mean_gx'] = float(np.mean(np.abs(gx)))\n    feats[f'{prefix}_mean_gy'] = float(np.mean(np.abs(gy)))\n    return feats\n\n\ndef vegetation_indices(rgb):\n    r = rgb[...,0]; g = rgb[...,1]; b = rgb[...,2]\n    eps = 1e-6\n    exg  = 2*g - r - b\n    vari = (g - r) / (g + r - b + eps)\n    ndi  = (g - r) / (g + r + eps)\n    cive = 0.441*r - 0.811*g + 0.385*b + 18.787\n    tgi  = -0.5*(190*(r-g) - 120*(r-b))\n    feats = {}\n    feats.update(basic_stats(exg,  'exg'))\n    feats.update(basic_stats(vari, 'vari'))\n    feats.update(basic_stats(ndi,  'ndi'))\n    feats.update(basic_stats(cive, 'cive'))\n    feats.update(basic_stats(tgi,  'tgi'))\n    feats['green_dom_frac'] = float(np.mean((g>r) & (g>b)))\n    feats['exg_pos_frac']   = float(np.mean(exg>0))\n    return feats\n\n\ndef extract_features_from_image(rgb):\n    H,W,_ = rgb.shape\n    feats = {'img_h':float(H),'img_w':float(W),'img_ratio':float(W/max(H,1))}\n    r,g,b = rgb[...,0], rgb[...,1], rgb[...,2]\n    feats.update(basic_stats(r,'r'))\n    feats.update(basic_stats(g,'g'))\n    feats.update(basic_stats(b,'b'))\n    feats.update(channel_hist(r,16,'r'))\n    feats.update(channel_hist(g,16,'g'))\n    feats.update(channel_hist(b,16,'b'))\n    hsv = rgb_to_hsv_np(rgb)\n    h,s,v = hsv[...,0], hsv[...,1], hsv[...,2]\n    feats.update(basic_stats(h,'h'))\n    feats.update(basic_stats(s,'s'))\n    feats.update(basic_stats(v,'v'))\n    gray = 0.299*r + 0.587*g + 0.114*b\n    feats.update(gradient_stats(gray,'grad'))\n    feats['lap_var'] = laplacian_var(gray)\n    feats.update(vegetation_indices(rgb))\n    return feats\n\n_FEATURE_CACHE = {}\n\ndef compute_image_features(image_path: str):\n    abspath = resolve_image_path(image_path)\n    if abspath in _FEATURE_CACHE:\n        return _FEATURE_CACHE[abspath]\n    if not os.path.exists(abspath):\n        raise FileNotFoundError(f\"Зураг олдсонгүй: {abspath}\")\n    rgb = load_image_rgb(abspath)\n    feats = extract_features_from_image(rgb)\n    _FEATURE_CACHE[abspath] = feats\n    return feats\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:35:14.578026Z","iopub.execute_input":"2025-11-05T02:35:14.578316Z","iopub.status.idle":"2025-11-05T02:35:14.599409Z","shell.execute_reply.started":"2025-11-05T02:35:14.578295Z","shell.execute_reply":"2025-11-05T02:35:14.598788Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 4) Train matrices\n# =====================\n\nuniq_images = train_df[['image_path']].drop_duplicates().reset_index(drop=True)\nrows = []\nstart = time.time()\nfor i, r in uniq_images.iterrows():\n    if (i+1)%100==0:\n        print(f\"Features: {i+1}/{len(uniq_images)}\")\n    feats = compute_image_features(r['image_path'])\n    feats['image_path'] = r['image_path']\n    rows.append(feats)\nprint(f\"Feature extraction finished in {time.time()-start:.1f}s for {len(uniq_images)} images\")\n\nfeatures_df = pd.DataFrame(rows)\n\n# pivot targets to wide\npiv = train_df.pivot_table(index='image_path', columns='target_name', values='target', aggfunc='mean').reset_index()\nfor t in TARGETS:\n    if t not in piv.columns:\n        piv[t] = np.nan\n\nfull_df = features_df.merge(piv, on='image_path', how='inner')\nprint('full_df:', full_df.shape)\n\nY = full_df[TARGETS].values\nY_log = np.log1p(np.clip(Y, a_min=0, a_max=None))\nfeature_cols = [c for c in full_df.columns if c not in (['image_path'] + TARGETS)]\nX = full_df[feature_cols].values.astype(np.float32)\nprint('X:', X.shape, 'Y:', Y.shape)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:35:19.355400Z","iopub.execute_input":"2025-11-05T02:35:19.356182Z","iopub.status.idle":"2025-11-05T02:36:55.418500Z","shell.execute_reply.started":"2025-11-05T02:35:19.356155Z","shell.execute_reply":"2025-11-05T02:36:55.417828Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 5) 5-Fold OOF — Base models & Meta blender (Stacking)\n# =====================\n\nn_splits = 5\nkf = KFold(n_splits=n_splits, shuffle=True, random_state=SEED)\n\nn = X.shape[0]\nM = 4  # number of base models\n# OOF predictions for each base model, log-space\nOOF_LOG = np.zeros((n, len(TARGETS), M), dtype=np.float32)\n\nfold_idx = 0\nfor tr_idx, va_idx in kf.split(X):\n    fold_idx += 1\n    print(f\"\\n===== Fold {fold_idx}/{n_splits} =====\")\n    X_tr, X_va = X[tr_idx], X[va_idx]\n    Y_tr, Y_va = Y_log[tr_idx], Y_log[va_idx]\n\n    # Base-1: RandomForest (multioutput)\n    rf = RandomForestRegressor(\n        n_estimators=500,\n        max_depth=None,\n        min_samples_leaf=2,\n        n_jobs=-1,\n        random_state=SEED+fold_idx,\n    )\n    rf.fit(X_tr, Y_tr)\n    pred_rf = rf.predict(X_va)\n\n    # Base-2: ExtraTrees (multioutput)\n    et = ExtraTreesRegressor(\n        n_estimators=500,\n        max_depth=None,\n        min_samples_leaf=2,\n        n_jobs=-1,\n        random_state=SEED+100+fold_idx,\n    )\n    et.fit(X_tr, Y_tr)\n    pred_et = et.predict(X_va)\n\n    # Base-3: Ridge (scaled)\n    rg = Pipeline([\n        ('scaler', StandardScaler(with_mean=True, with_std=True)),\n        ('ridge', Ridge(alpha=1.0, random_state=SEED+fold_idx))\n    ])\n    rg.fit(X_tr, Y_tr)\n    pred_rg = rg.predict(X_va)\n\n    # Base-4: GradientBoosting (wrapped for multioutput)\n    gbr_base = GradientBoostingRegressor(\n        n_estimators=300,\n        learning_rate=0.05,\n        max_depth=3,\n        subsample=1.0,\n        random_state=SEED+fold_idx,\n    )\n    gbr = MultiOutputRegressor(gbr_base)\n    gbr.fit(X_tr, Y_tr)\n    pred_gb = gbr.predict(X_va)\n\n    # collect OOF base preds (log-space)\n    OOF_LOG[va_idx, :, 0] = pred_rf\n    OOF_LOG[va_idx, :, 1] = pred_et\n    OOF_LOG[va_idx, :, 2] = pred_rg\n    OOF_LOG[va_idx, :, 3] = pred_gb\n\n    # quick fold metric using simple average\n    pred_log_avg = (pred_rf + pred_et + pred_rg + pred_gb) / 4.0\n    pred_avg = np.expm1(pred_log_avg)\n    fold_rmse = [math.sqrt(mean_squared_error(np.expm1(Y_va[:,j]), pred_avg[:,j])) for j in range(len(TARGETS))]\n    print('Fold simple-avg mean RMSE:', np.mean(fold_rmse))\n\n# Train meta blender (one Ridge per target) on OOF base predictions\nMETA = []  # list of (scaler, ridge) pipelines per target\nOOF_final = np.zeros((n, len(TARGETS)), dtype=np.float32)\n\nfor j, t in enumerate(TARGETS):\n    Z = OOF_LOG[:, j, :]  # (n_samples, M)\n    meta = Pipeline([\n        ('scaler', StandardScaler(with_mean=True, with_std=True)),\n        ('ridge', Ridge(alpha=1.0, random_state=SEED))\n    ])\n    meta.fit(Z, Y_log[:, j])\n    META.append(meta)\n    pred_log = meta.predict(Z)\n    OOF_final[:, j] = np.expm1(pred_log)\n\n# OOF metrics for meta-stacked predictions\noof_rmse = []\nfor j, t in enumerate(TARGETS):\n    rmse = math.sqrt(mean_squared_error(np.expm1(Y_log[:, j]), OOF_final[:, j]))\n    oof_rmse.append(rmse)\n    print(f\"OOF {t}: RMSE = {rmse:.4f}\")\nprint(f\"OOF mean RMSE (Stacked) = {np.mean(oof_rmse):.4f}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:37:29.796627Z","iopub.execute_input":"2025-11-05T02:37:29.797095Z","iopub.status.idle":"2025-11-05T02:39:23.184003Z","shell.execute_reply.started":"2025-11-05T02:37:29.797070Z","shell.execute_reply":"2025-11-05T02:39:23.183057Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 6) All data base models full-fit + Test inference (Stacked)\n# =====================\n\n# 6.1 Full-fit base models on entire training set\nrf_final = RandomForestRegressor(\n    n_estimators=700,\n    max_depth=None,\n    min_samples_leaf=2,\n    n_jobs=-1,\n    random_state=SEED,\n)\nrf_final.fit(X, Y_log)\n\net_final = ExtraTreesRegressor(\n    n_estimators=700,\n    max_depth=None,\n    min_samples_leaf=2,\n    n_jobs=-1,\n    random_state=SEED+111,\n)\net_final.fit(X, Y_log)\n\nrg_final = Pipeline([\n    ('scaler', StandardScaler(with_mean=True, with_std=True)),\n    ('ridge', Ridge(alpha=1.0, random_state=SEED))\n])\nrg_final.fit(X, Y_log)\n\ngbr_final = MultiOutputRegressor(GradientBoostingRegressor(\n    n_estimators=400,\n    learning_rate=0.05,\n    max_depth=3,\n    subsample=1.0,\n    random_state=SEED,\n))\ngbr_final.fit(X, Y_log)\n\nprint('Full-fit base models trained.')\n\n# 6.2 Build test features\nuniq_test = test_df[['image_path']].drop_duplicates().reset_index(drop=True)\nrows = []\nstart = time.time()\nfor i, r in uniq_test.iterrows():\n    if (i+1)%100==0:\n        print(f\"Test features: {i+1}/{len(uniq_test)}\")\n    feats = compute_image_features(r['image_path'])\n    feats['image_path'] = r['image_path']\n    rows.append(feats)\nprint(f\"Test feature extraction in {time.time()-start:.1f}s for {len(uniq_test)} images\")\n\nfeat_test = pd.DataFrame(rows)\nfor c in [col for col in feature_cols if col not in feat_test.columns]:\n    feat_test[c] = 0.0\nX_test = feat_test[feature_cols].values.astype(np.float32)\n\n# 6.3 Base predictions (log-space)\nlog_rf = rf_final.predict(X_test)\nlog_et = et_final.predict(X_test)\nlog_rg = rg_final.predict(X_test)\nlog_gb = gbr_final.predict(X_test)\n\n# 6.4 Stack through meta blender per target\nZ_test_all = np.stack([log_rf, log_et, log_rg, log_gb], axis=2)  # (n_test, n_targets, M)\nfinal_log = np.zeros_like(log_rf)\n\nfor j in range(len(TARGETS)):\n    Zj = Z_test_all[:, j, :]  # (n_test, M)\n    final_log[:, j] = META[j].predict(Zj)\n\npred_test = np.expm1(final_log)\n# clip to >= 0 for safety\npred_test = np.clip(pred_test, a_min=0.0, a_max=None)\n\npred_df = pd.DataFrame(pred_test, columns=TARGETS)\npred_df['image_path'] = feat_test['image_path']\n\n# 6.5 Map to submission rows\nsub = test_df.merge(pred_df, on='image_path', how='left')\nname_to_col = {t:t for t in TARGETS}\nvalues = []\nfor _, r in sub.iterrows():\n    values.append(r[name_to_col.get(r['target_name'], TARGETS[-1])])\n\nsubmission = pd.DataFrame({'sample_id': sub['sample_id'], 'target': values})\n\nout_path = 'submission.csv'\nsubmission.to_csv(out_path, index=False)\nprint('Saved:', out_path, 'shape =', submission.shape)\nprint(submission.head(10))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:40:11.719623Z","iopub.execute_input":"2025-11-05T02:40:11.720324Z","iopub.status.idle":"2025-11-05T02:40:50.351796Z","shell.execute_reply.started":"2025-11-05T02:40:11.720298Z","shell.execute_reply":"2025-11-05T02:40:50.350904Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# =====================\n# 7) Fast visualition(Optional)\n# =====================\ntry:\n    samp = test_df.sample(n=min(6, len(test_df)), random_state=SEED)\n    fig, axes = plt.subplots(2,3, figsize=(12,8)); axes = axes.ravel()\n    for i, p in enumerate(samp['image_path'].tolist()[:6]):\n        im = load_image_rgb(resolve_image_path(p), max_side=512)\n        axes[i].imshow(im); axes[i].set_title(os.path.basename(p)); axes[i].axis('off')\n    plt.tight_layout(); plt.show()\nexcept Exception as e:\n    print('Viz skipped:', e)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T02:40:50.353144Z","iopub.execute_input":"2025-11-05T02:40:50.353518Z","iopub.status.idle":"2025-11-05T02:40:51.272399Z","shell.execute_reply.started":"2025-11-05T02:40:50.353500Z","shell.execute_reply":"2025-11-05T02:40:51.271648Z"}},"outputs":[],"execution_count":null}]}