{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","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":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"},{"sourceId":292296638,"sourceType":"kernelVersion"},{"sourceId":292448872,"sourceType":"kernelVersion"},{"sourceId":268942,"sourceType":"modelInstanceVersion","modelInstanceId":230141,"modelId":251887},{"sourceId":711299,"sourceType":"modelInstanceVersion","modelInstanceId":540408,"modelId":553604},{"sourceId":719660,"sourceType":"modelInstanceVersion","modelInstanceId":547330,"modelId":560139}],"dockerImageVersionId":31260,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🌾 CSIRO Image2Biomass: Comprehensive EDA & Image Explorer\n\n**Competition:** [CSIRO Image2Biomass Prediction](https://www.kaggle.com/competitions/csiro-biomass)\n\n---\n\n## 📋 **Notebook Overview**\n\nThis notebook provides a **complete exploratory data analysis** for the CSIRO Image2Biomass competition, with a special focus on understanding the relationship between images and their corresponding biomass measurements.\n\n### **What You'll Get:**\n\n1. **📊 Data Quality Assessment**\n   - Missing value analysis\n   - Outlier detection (IQR & Z-score methods)\n   - Physical constraint validation\n\n2. **📈 Statistical Analysis**\n   - Target variable distributions\n   - Correlation patterns\n   - Key statistics and metrics\n\n3. **🗺️ Metadata Exploration**\n   - Geographic distribution (States)\n   - Temporal patterns (Seasons, Years)\n   - Species composition\n   - NDVI and Height relationships\n\n4. **🖼️ Interactive Image Viewer** ⭐\n   - Sort images by any metric\n   - Filter by State, Season, Species\n   - View detailed information\n   - Compare high/low biomass samples\n\n5. **🔍 Insights for Modeling**\n   - Feature importance indicators\n   - Data imbalances\n   - Potential challenges\n\n---\n\n### **Key Features:**\n\n- ✅ Clean, modular, reusable code\n- ✅ Comprehensive visualizations\n- ✅ Interactive image exploration\n- ✅ Physical constraint checking\n- ✅ Easy to adapt and extend\n\n---\n\n### **Competition Context:**\n\n**Goal:** Predict 5 biomass components from pasture images:\n- `Dry_Green_g` (10% weight)\n- `Dry_Dead_g` (10% weight)\n- `Dry_Clover_g` (10% weight)\n- `GDM_g` (20% weight)\n- `Dry_Total_g` (50% weight) ⭐\n\n**Evaluation:** Weighted R² score\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 1️⃣ Setup & Configuration","metadata":{}},{"cell_type":"code","source":"# Create a folder to hold the extracted wheels\n!mkdir -p /kaggle/working/timm_wheels\n\n# Unzip the archive from your input dataset\n!unzip -q /kaggle/input/install-dependency-in-another-notebook/timm_1_0_22_wheels.zip -d /kaggle/working/timm_wheels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:02:37.610237Z","iopub.execute_input":"2026-01-17T23:02:37.610548Z","iopub.status.idle":"2026-01-17T23:03:39.917360Z","shell.execute_reply.started":"2026-01-17T23:02:37.610519Z","shell.execute_reply":"2026-01-17T23:03:39.916559Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Point --find-links to the folder where you just unzipped the files\n!pip install -U --no-deps timm==1.0.22 \\\n    --no-index \\\n    --find-links=/kaggle/working/timm_wheels/kaggle/temp/wheels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:03:39.919363Z","iopub.execute_input":"2026-01-17T23:03:39.919644Z","iopub.status.idle":"2026-01-17T23:03:45.639163Z","shell.execute_reply.started":"2026-01-17T23:03:39.919613Z","shell.execute_reply":"2026-01-17T23:03:45.638158Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Import libraries\nimport os\nimport gc\nimport random\nimport cv2\nimport torch\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom copy import deepcopy\nfrom dataclasses import dataclass\nfrom pathlib import Path\nfrom tqdm.auto import tqdm\nfrom PIL import Image\nfrom scipy import stats\nimport warnings\n\n# Sklearn & Models\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.cross_decomposition import PLSRegression\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.ensemble import GradientBoostingRegressor, HistGradientBoostingRegressor\nfrom lightgbm import LGBMRegressor\nfrom catboost import CatBoostRegressor\nfrom transformers import AutoModel, AutoImageProcessor, AutoTokenizer\n\nwarnings.filterwarnings('ignore')\nos.environ[\"TOKENIZERS_PARALLELISM\"] = \"false\"\n\n# Set plotting style\nplt.style.use('seaborn-v0_8-darkgrid')\nsns.set_palette(\"husl\")\n\n# Display settings\npd.set_option('display.max_columns', 100)\npd.set_option('display.max_rows', 100)\npd.set_option('display.width', 1000)\npd.set_option('display.float_format', '{:.2f}'.format)\n\nprint(\"✅ Libraries imported successfully\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:03:45.640650Z","iopub.execute_input":"2026-01-17T23:03:45.641279Z","iopub.status.idle":"2026-01-17T23:04:36.340204Z","shell.execute_reply.started":"2026-01-17T23:03:45.641245Z","shell.execute_reply":"2026-01-17T23:04:36.339453Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Configuration\nCSV_PATH = \"/kaggle/input/csiro-biomass/train.csv\"\nIMAGE_DIR = \"/kaggle/input/csiro-biomass\"\n\n# Target columns\nTARGET_COLS = ['Dry_Green_g', 'Dry_Dead_g', 'Dry_Clover_g', 'GDM_g', 'Dry_Total_g']\n\n# Competition weights\nWEIGHTS = {\n    'Dry_Green_g': 0.1,\n    'Dry_Dead_g': 0.1,\n    'Dry_Clover_g': 0.1,\n    'GDM_g': 0.2,\n    'Dry_Total_g': 0.5\n}\n\nprint(\"✅ Configuration set\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.341986Z","iopub.execute_input":"2026-01-17T23:04:36.342666Z","iopub.status.idle":"2026-01-17T23:04:36.347500Z","shell.execute_reply.started":"2026-01-17T23:04:36.342642Z","shell.execute_reply":"2026-01-17T23:04:36.346769Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2️⃣ Data Loading & Preprocessing\n\n### What we'll do:\n1. Load the CSV data (long format: 5 rows per image)\n2. Convert to wide format (1 row per image)\n3. Add temporal features (season, year, month)\n4. Add derived features (ratios, physical constraints)","metadata":{}},{"cell_type":"code","source":"def load_and_preprocess_data(csv_path: str) -> tuple[pd.DataFrame, pd.DataFrame]:\n    \"\"\"\n    Load and preprocess biomass data\n    \n    Returns:\n        Tuple of (long_format_df, wide_format_df)\n    \"\"\"\n    # Load data\n    print(\"📂 Loading data...\")\n    df_long = pd.read_csv(csv_path)\n    print(f\"✅ Loaded {len(df_long):,} rows × {len(df_long.columns)} columns\")\n    \n    # Extract image ID\n    df_long['image_id'] = df_long['sample_id'].str.split('__').str[0]\n    \n    # Convert to wide format\n    print(\"\\n🔄 Converting to wide format...\")\n    meta_cols = ['image_path', 'Sampling_Date', 'State', 'Species', \n                 'Pre_GSHH_NDVI', 'Height_Ave_cm']\n    \n    df_wide = df_long.pivot_table(\n        index=['image_id'] + meta_cols,\n        columns='target_name',\n        values='target',\n        aggfunc='first'\n    ).reset_index()\n    \n    n_images = len(df_wide)\n    print(f\"✅ Wide format: {n_images:,} images\")\n    \n    # Add temporal features\n    print(\"\\n📅 Adding temporal features...\")\n    df_wide['date'] = pd.to_datetime(df_wide['Sampling_Date'], format='%Y/%m/%d')\n    df_wide['year'] = df_wide['date'].dt.year\n    df_wide['month'] = df_wide['date'].dt.month\n    df_wide['day'] = df_wide['date'].dt.day\n    df_wide['dayofweek'] = df_wide['date'].dt.dayofweek\n    df_wide['dayofyear'] = df_wide['date'].dt.dayofyear\n    \n    # Add season (Southern Hemisphere)\n    def get_season(month: int) -> str:\n        if month in [9, 10, 11]:\n            return \"Spring\"\n        elif month in [12, 1, 2]:\n            return \"Summer\"\n        elif month in [3, 4, 5]:\n            return \"Autumn\"\n        else:\n            return \"Winter\"\n    \n    df_wide['season'] = df_wide['month'].apply(get_season)\n    \n    # Add derived features\n    print(\"🧮 Adding derived features...\")\n    \n    # Component ratios\n    total = df_wide['Dry_Total_g'] + 1e-8\n    df_wide['green_ratio'] = df_wide['Dry_Green_g'] / total\n    df_wide['dead_ratio'] = df_wide['Dry_Dead_g'] / total\n    df_wide['clover_ratio'] = df_wide['Dry_Clover_g'] / total\n    \n    # Physical constraint check\n    df_wide['sum_components'] = (\n        df_wide['Dry_Green_g'] + \n        df_wide['Dry_Dead_g'] + \n        df_wide['Dry_Clover_g']\n    )\n    df_wide['constraint_diff'] = df_wide['Dry_Total_g'] - df_wide['sum_components']\n    df_wide['constraint_error'] = np.abs(df_wide['constraint_diff'])\n    \n    print(\"✅ Preprocessing complete!\\n\")\n    \n    return df_long, df_wide\n\n# Execute\ndf_long, df = load_and_preprocess_data(CSV_PATH)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.348479Z","iopub.execute_input":"2026-01-17T23:04:36.348838Z","iopub.status.idle":"2026-01-17T23:04:36.513109Z","shell.execute_reply.started":"2026-01-17T23:04:36.348816Z","shell.execute_reply":"2026-01-17T23:04:36.512391Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Quick peek at the data\nprint(\"=\"*80)\nprint(\"📊 DATA OVERVIEW\")\nprint(\"=\"*80)\nprint(f\"\\nShape: {df.shape[0]:,} rows × {df.shape[1]} columns\")\nprint(f\"Memory: {df.memory_usage(deep=True).sum() / 1024**2:.2f} MB\")\nprint(f\"\\nDate range: {df['date'].min().date()} to {df['date'].max().date()}\")\nprint(f\"Number of states: {df['State'].nunique()}\")\nprint(f\"Number of species combinations: {df['Species'].nunique()}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.514016Z","iopub.execute_input":"2026-01-17T23:04:36.514287Z","iopub.status.idle":"2026-01-17T23:04:36.524951Z","shell.execute_reply.started":"2026-01-17T23:04:36.514258Z","shell.execute_reply":"2026-01-17T23:04:36.524240Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Display first few rows\ndf.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.526115Z","iopub.execute_input":"2026-01-17T23:04:36.526442Z","iopub.status.idle":"2026-01-17T23:04:36.565046Z","shell.execute_reply.started":"2026-01-17T23:04:36.526421Z","shell.execute_reply":"2026-01-17T23:04:36.564452Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 3️⃣ Data Quality Assessment\n\nLet's check for:\n- Missing values\n- Data types\n- Basic statistics","metadata":{}},{"cell_type":"code","source":"def analyze_data_quality(df: pd.DataFrame) -> None:\n    \"\"\"Comprehensive data quality analysis\"\"\"\n    \n    print(\"\\n\" + \"=\"*80)\n    print(\"🔍 DATA QUALITY ANALYSIS\")\n    print(\"=\"*80)\n    \n    # Missing values\n    print(\"\\n【Missing Values】\")\n    missing = pd.DataFrame({\n        'Missing_Count': df.isnull().sum(),\n        'Missing_Pct': (df.isnull().sum() / len(df) * 100).round(2)\n    })\n    missing = missing[missing['Missing_Count'] > 0].sort_values('Missing_Count', ascending=False)\n    \n    if len(missing) == 0:\n        print(\"✅ No missing values!\")\n    else:\n        print(missing)\n    \n    # Data types\n    print(\"\\n【Data Types】\")\n    print(df.dtypes.value_counts())\n    \n    # Duplicates\n    n_duplicates = df.duplicated(subset='image_id').sum()\n    print(f\"\\n【Duplicates】\")\n    print(f\"Duplicate image IDs: {n_duplicates}\")\n\nanalyze_data_quality(df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.565836Z","iopub.execute_input":"2026-01-17T23:04:36.566149Z","iopub.status.idle":"2026-01-17T23:04:36.578972Z","shell.execute_reply.started":"2026-01-17T23:04:36.566126Z","shell.execute_reply":"2026-01-17T23:04:36.578088Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4️⃣ Target Variable Analysis\n\n### Key Questions:\n- What are the distributions of each target?\n- Are they normally distributed?\n- What are the correlations between targets?\n- Are there any outliers?","metadata":{}},{"cell_type":"code","source":"def analyze_target_statistics(df: pd.DataFrame, target_cols: list) -> pd.DataFrame:\n    \"\"\"Compute comprehensive target statistics\"\"\"\n    \n    print(\"\\n\" + \"=\"*80)\n    print(\"📈 TARGET VARIABLE STATISTICS\")\n    print(\"=\"*80)\n    \n    stats = df[target_cols].describe().T\n    \n    # Additional metrics\n    stats['cv'] = stats['std'] / stats['mean']  # Coefficient of variation\n    stats['skewness'] = df[target_cols].skew()\n    stats['kurtosis'] = df[target_cols].kurtosis()\n    stats['comp_weight'] = [WEIGHTS[col] for col in target_cols]\n    \n    return stats.round(3)\n\ntarget_stats = analyze_target_statistics(df, TARGET_COLS)\ntarget_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.579830Z","iopub.execute_input":"2026-01-17T23:04:36.580126Z","iopub.status.idle":"2026-01-17T23:04:36.621312Z","shell.execute_reply.started":"2026-01-17T23:04:36.580099Z","shell.execute_reply":"2026-01-17T23:04:36.620581Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Interpretation\nprint(\"\\n💡 Key Observations:\")\nprint(f\"• Dry_Total_g has the highest mean ({target_stats.loc['Dry_Total_g', 'mean']:.1f}g) \" +\n      \"and is the primary target (50% weight)\")\nprint(f\"• Dry_Clover_g has high variability (CV: {target_stats.loc['Dry_Clover_g', 'cv']:.2f})\")\nprint(f\"• All targets show positive skewness (right-tailed distributions)\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.623831Z","iopub.execute_input":"2026-01-17T23:04:36.624236Z","iopub.status.idle":"2026-01-17T23:04:36.629226Z","shell.execute_reply.started":"2026-01-17T23:04:36.624214Z","shell.execute_reply":"2026-01-17T23:04:36.628477Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5️⃣ Target Distribution Visualization\n\nLet's visualize the distributions with:\n- Histograms (with mean/median lines)\n- Box plots (for outlier detection)\n- Q-Q plots (for normality check)","metadata":{}},{"cell_type":"code","source":"def plot_target_distributions(df: pd.DataFrame, target_cols: list) -> None:\n    \"\"\"Plot comprehensive target distributions\"\"\"\n    \n    fig, axes = plt.subplots(3, 5, figsize=(20, 12))\n    \n    for i, target in enumerate(target_cols):\n        data = df[target].dropna()\n        \n        # Histogram\n        axes[0, i].hist(data, bins=50, edgecolor='black', alpha=0.7, color='steelblue')\n        axes[0, i].axvline(data.mean(), color='red', linestyle='--', \n                          linewidth=2, label=f'Mean: {data.mean():.1f}g')\n        axes[0, i].axvline(data.median(), color='green', linestyle='--', \n                          linewidth=2, label=f'Median: {data.median():.1f}g')\n        axes[0, i].set_title(f'{target}\\nDistribution', fontsize=11, fontweight='bold')\n        axes[0, i].set_xlabel('Value (g)')\n        axes[0, i].set_ylabel('Frequency')\n        axes[0, i].legend(fontsize=8)\n        axes[0, i].grid(True, alpha=0.3)\n        \n        # Box plot\n        bp = axes[1, i].boxplot(data, vert=True, patch_artist=True)\n        bp['boxes'][0].set_facecolor('lightblue')\n        axes[1, i].set_title(f'{target}\\nBox Plot', fontsize=11, fontweight='bold')\n        axes[1, i].set_ylabel('Value (g)')\n        axes[1, i].grid(True, alpha=0.3)\n        \n        # Q-Q plot\n        stats.probplot(data, dist=\"norm\", plot=axes[2, i])\n        axes[2, i].set_title(f'{target}\\nQ-Q Plot', fontsize=11, fontweight='bold')\n        axes[2, i].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n\nplot_target_distributions(df, TARGET_COLS)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:36.630211Z","iopub.execute_input":"2026-01-17T23:04:36.630450Z","iopub.status.idle":"2026-01-17T23:04:39.147237Z","shell.execute_reply.started":"2026-01-17T23:04:36.630430Z","shell.execute_reply":"2026-01-17T23:04:39.146270Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"\\n💡 Distribution Insights:\")\nprint(\"• All targets show right-skewed distributions (more low values than high)\")\nprint(\"• Q-Q plots show deviations from normality, especially at the tails\")\nprint(\"• Consider log transformation for modeling\")\nprint(\"• Several outliers present in all targets\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:39.148474Z","iopub.execute_input":"2026-01-17T23:04:39.148836Z","iopub.status.idle":"2026-01-17T23:04:39.153817Z","shell.execute_reply.started":"2026-01-17T23:04:39.148808Z","shell.execute_reply":"2026-01-17T23:04:39.153003Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6️⃣ Correlation Analysis\n\nUnderstanding relationships between targets is crucial for:\n- Multi-task learning strategies\n- Feature engineering\n- Model architecture design","metadata":{}},{"cell_type":"code","source":"def plot_correlation_analysis(df: pd.DataFrame, target_cols: list) -> None:\n    \"\"\"Plot correlation heatmap and key insights\"\"\"\n    \n    corr = df[target_cols].corr()\n    \n    plt.figure(figsize=(10, 8))\n    sns.heatmap(corr, annot=True, fmt='.3f', cmap='coolwarm', \n               center=0, square=True, linewidths=2,\n               cbar_kws={\"shrink\": 0.8}, vmin=-1, vmax=1)\n    plt.title('Target Variable Correlations', fontsize=14, fontweight='bold', pad=20)\n    plt.tight_layout()\n    plt.show()\n    \n    # Print key correlations\n    print(\"\\n💡 Key Correlations:\")\n    \n    # Get upper triangle\n    corr_flat = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))\n    corr_pairs = corr_flat.stack().sort_values(ascending=False)\n    \n    print(\"\\nStrongest positive correlations:\")\n    for (var1, var2), corr_val in corr_pairs.head(3).items():\n        print(f\"  • {var1} ↔ {var2}: {corr_val:.3f}\")\n    \n    print(\"\\nWeakest correlations:\")\n    for (var1, var2), corr_val in corr_pairs.tail(3).items():\n        print(f\"  • {var1} ↔ {var2}: {corr_val:.3f}\")\n\nplot_correlation_analysis(df, TARGET_COLS)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:39.154892Z","iopub.execute_input":"2026-01-17T23:04:39.155184Z","iopub.status.idle":"2026-01-17T23:04:39.456459Z","shell.execute_reply.started":"2026-01-17T23:04:39.155161Z","shell.execute_reply":"2026-01-17T23:04:39.455834Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7️⃣ Metadata Exploration\n\nLet's explore the metadata features:\n- Geographic (State)\n- Temporal (Season, Year)\n- Physical (NDVI, Height)\n- Biological (Species)","metadata":{}},{"cell_type":"code","source":"def plot_metadata_distributions(df: pd.DataFrame) -> None:\n    \"\"\"Plot metadata feature distributions\"\"\"\n    \n    fig, axes = plt.subplots(2, 3, figsize=(18, 10))\n    axes = axes.flatten()\n    \n    # State distribution\n    state_counts = df['State'].value_counts()\n    bars = axes[0].bar(state_counts.index, state_counts.values, \n                      edgecolor='black', color='skyblue', alpha=0.8)\n    axes[0].set_title('State Distribution', fontsize=12, fontweight='bold')\n    axes[0].set_xlabel('State')\n    axes[0].set_ylabel('Sample Count')\n    axes[0].grid(True, alpha=0.3, axis='y')\n    for bar, count in zip(bars, state_counts.values):\n        axes[0].text(bar.get_x() + bar.get_width()/2, count, \n                    str(count), ha='center', va='bottom', fontweight='bold')\n    \n    # Top 10 species\n    species_counts = df['Species'].value_counts().head(10)\n    axes[1].barh(range(len(species_counts)), species_counts.values, \n                edgecolor='black', color='lightgreen', alpha=0.8)\n    axes[1].set_yticks(range(len(species_counts)))\n    axes[1].set_yticklabels([s[:30] for s in species_counts.index], fontsize=9)\n    axes[1].set_title('Top 10 Species Combinations', fontsize=12, fontweight='bold')\n    axes[1].set_xlabel('Sample Count')\n    axes[1].invert_yaxis()\n    axes[1].grid(True, alpha=0.3, axis='x')\n    \n    # NDVI distribution\n    ndvi_data = df['Pre_GSHH_NDVI'].dropna()\n    axes[2].hist(ndvi_data, bins=30, edgecolor='black', alpha=0.7, color='green')\n    axes[2].axvline(ndvi_data.mean(), color='red', linestyle='--', \n                   linewidth=2, label=f'Mean: {ndvi_data.mean():.2f}')\n    axes[2].axvline(ndvi_data.median(), color='orange', linestyle='--', \n                   linewidth=2, label=f'Median: {ndvi_data.median():.2f}')\n    axes[2].set_title('NDVI Distribution', fontsize=12, fontweight='bold')\n    axes[2].set_xlabel('NDVI Value')\n    axes[2].set_ylabel('Frequency')\n    axes[2].legend()\n    axes[2].grid(True, alpha=0.3)\n    \n    # Height distribution\n    height_data = df['Height_Ave_cm'].dropna()\n    axes[3].hist(height_data, bins=30, edgecolor='black', alpha=0.7, color='orange')\n    axes[3].axvline(height_data.mean(), color='red', linestyle='--', \n                   linewidth=2, label=f'Mean: {height_data.mean():.1f}cm')\n    axes[3].axvline(height_data.median(), color='blue', linestyle='--', \n                   linewidth=2, label=f'Median: {height_data.median():.1f}cm')\n    axes[3].set_title('Height Distribution', fontsize=12, fontweight='bold')\n    axes[3].set_xlabel('Height (cm)')\n    axes[3].set_ylabel('Frequency')\n    axes[3].legend()\n    axes[3].grid(True, alpha=0.3)\n    \n    # Season distribution\n    season_counts = df['season'].value_counts()\n    colors = ['#90EE90', '#FFD700', '#FF6347', '#87CEEB']\n    wedges, texts, autotexts = axes[4].pie(\n        season_counts.values, labels=season_counts.index, \n        autopct='%1.1f%%', startangle=90, colors=colors\n    )\n    for autotext in autotexts:\n        autotext.set_color('white')\n        autotext.set_fontweight('bold')\n    axes[4].set_title('Season Distribution', fontsize=12, fontweight='bold')\n    \n    # Yearly sample counts\n    year_counts = df['year'].value_counts().sort_index()\n    axes[5].plot(year_counts.index, year_counts.values, \n                marker='o', linewidth=2, markersize=8, color='purple')\n    axes[5].set_title('Samples by Year', fontsize=12, fontweight='bold')\n    axes[5].set_xlabel('Year')\n    axes[5].set_ylabel('Sample Count')\n    axes[5].grid(True, alpha=0.3)\n    for x, y in zip(year_counts.index, year_counts.values):\n        axes[5].text(x, y, str(y), ha='center', va='bottom', fontweight='bold')\n    \n    plt.tight_layout()\n    plt.show()\n\nplot_metadata_distributions(df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:39.457505Z","iopub.execute_input":"2026-01-17T23:04:39.458236Z","iopub.status.idle":"2026-01-17T23:04:40.384350Z","shell.execute_reply.started":"2026-01-17T23:04:39.458208Z","shell.execute_reply":"2026-01-17T23:04:40.383513Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"\\n💡 Metadata Insights:\")\nprint(f\"• Data collected across {df['State'].nunique()} Australian states\")\nprint(f\"• {df['Species'].nunique()} unique species combinations\")\nprint(f\"• NDVI range: {df['Pre_GSHH_NDVI'].min():.2f} to {df['Pre_GSHH_NDVI'].max():.2f}\")\nprint(f\"• Height range: {df['Height_Ave_cm'].min():.1f} to {df['Height_Ave_cm'].max():.1f} cm\")\nprint(f\"• Data spans {df['year'].nunique()} years ({df['year'].min()}-{df['year'].max()})\")\nprint(f\"• Most samples from: {df['State'].value_counts().index[0]}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:40.385477Z","iopub.execute_input":"2026-01-17T23:04:40.385860Z","iopub.status.idle":"2026-01-17T23:04:40.392588Z","shell.execute_reply.started":"2026-01-17T23:04:40.385830Z","shell.execute_reply":"2026-01-17T23:04:40.391770Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8️⃣ Target vs Metadata Relationships\n\nUnderstanding how targets vary across:\n- States\n- Seasons\n- NDVI values\n- Height measurements\n- Species","metadata":{}},{"cell_type":"code","source":"def plot_target_by_metadata(df: pd.DataFrame, target: str = 'Dry_Total_g') -> None:\n    \"\"\"Plot target relationships with metadata\"\"\"\n    \n    fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n    axes = axes.flatten()\n    \n    # By State\n    df.boxplot(column=target, by='State', ax=axes[0])\n    axes[0].set_title(f'{target} by State', fontsize=12, fontweight='bold')\n    axes[0].set_xlabel('State')\n    axes[0].set_ylabel(f'{target} (g)')\n    plt.sca(axes[0])\n    plt.xticks(rotation=0)\n    \n    # By Season\n    season_order = ['Spring', 'Summer', 'Autumn', 'Winter']\n    season_data = [df[df['season'] == s][target].dropna() for s in season_order]\n    bp = axes[1].boxplot(season_data, labels=season_order, patch_artist=True)\n    colors = ['#90EE90', '#FFD700', '#FF6347', '#87CEEB']\n    for patch, color in zip(bp['boxes'], colors):\n        patch.set_facecolor(color)\n    axes[1].set_title(f'{target} by Season', fontsize=12, fontweight='bold')\n    axes[1].set_xlabel('Season')\n    axes[1].set_ylabel(f'{target} (g)')\n    axes[1].grid(True, alpha=0.3)\n    \n    # vs NDVI\n    corr_ndvi = df[[target, 'Pre_GSHH_NDVI']].corr().iloc[0, 1]\n    axes[2].scatter(df['Pre_GSHH_NDVI'], df[target], alpha=0.5, s=30, edgecolor='black')\n    axes[2].set_title(f'{target} vs NDVI (r={corr_ndvi:.3f})', fontsize=12, fontweight='bold')\n    axes[2].set_xlabel('NDVI')\n    axes[2].set_ylabel(f'{target} (g)')\n    axes[2].grid(True, alpha=0.3)\n    \n    # vs Height\n    corr_height = df[[target, 'Height_Ave_cm']].corr().iloc[0, 1]\n    axes[3].scatter(df['Height_Ave_cm'], df[target], alpha=0.5, s=30, \n                   edgecolor='black', color='orange')\n    axes[3].set_title(f'{target} vs Height (r={corr_height:.3f})', fontsize=12, fontweight='bold')\n    axes[3].set_xlabel('Height (cm)')\n    axes[3].set_ylabel(f'{target} (g)')\n    axes[3].grid(True, alpha=0.3)\n    \n    # By Top 5 Species\n    top_species = df['Species'].value_counts().head(5).index\n    df_top = df[df['Species'].isin(top_species)]\n    species_data = [df_top[df_top['Species'] == s][target].dropna() for s in top_species]\n    bp = axes[4].boxplot(species_data, labels=[s[:20] for s in top_species])\n    axes[4].set_title(f'{target} by Top 5 Species', fontsize=12, fontweight='bold')\n    axes[4].set_xlabel('Species')\n    axes[4].set_ylabel(f'{target} (g)')\n    plt.sca(axes[4])\n    plt.xticks(rotation=45, ha='right')\n    axes[4].grid(True, alpha=0.3)\n    \n    # Monthly trend\n    monthly_mean = df.groupby('month')[target].mean()\n    axes[5].plot(monthly_mean.index, monthly_mean.values, \n                marker='o', linewidth=2, markersize=8, color='green')\n    axes[5].set_title(f'{target} Monthly Average', fontsize=12, fontweight='bold')\n    axes[5].set_xlabel('Month')\n    axes[5].set_ylabel(f'{target} Average (g)')\n    axes[5].set_xticks(range(1, 13))\n    axes[5].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n\n# Plot for main targets\nfor target in ['Dry_Total_g', 'Dry_Green_g', 'GDM_g']:\n    plot_target_by_metadata(df, target)\n    print(f\"\\n💡 {target} shows variation across states, seasons, and metadata\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:40.393665Z","iopub.execute_input":"2026-01-17T23:04:40.394017Z","iopub.status.idle":"2026-01-17T23:04:43.286715Z","shell.execute_reply.started":"2026-01-17T23:04:40.393985Z","shell.execute_reply":"2026-01-17T23:04:43.285924Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 9️⃣ Physical Constraint Validation\n\n### Important Physical Relationship:\n```\nDry_Total = Dry_Green + Dry_Dead + Dry_Clover\n```\n\nLet's check if this constraint holds in the data!","metadata":{}},{"cell_type":"code","source":"def check_physical_constraints(df: pd.DataFrame, threshold: float = 5.0) -> None:\n    \"\"\"Check and visualize physical constraint violations\"\"\"\n    \n    violations = df[df['constraint_error'] > threshold]\n    \n    print(\"\\n\" + \"=\"*80)\n    print(\"🔍 PHYSICAL CONSTRAINT CHECK\")\n    print(\"=\"*80)\n    print(f\"\\nConstraint: Dry_Total ≈ Dry_Green + Dry_Dead + Dry_Clover\")\n    print(f\"Threshold: {threshold}g\")\n    print(f\"Violations: {len(violations):,} / {len(df):,} ({len(violations)/len(df)*100:.1f}%)\")\n    \n    print(f\"\\nConstraint Error Statistics:\")\n    print(df['constraint_error'].describe())\n    \n    if len(violations) > 0:\n        print(f\"\\n⚠️ Top 10 Violations:\")\n        top_viol = violations.nlargest(10, 'constraint_error')[[\n            'image_id', 'Dry_Total_g', 'Dry_Green_g', 'Dry_Dead_g', \n            'Dry_Clover_g', 'sum_components', 'constraint_error'\n        ]]\n        print(top_viol.to_string(index=False))\n    \n    # Visualization\n    fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n    \n    # Error distribution\n    axes[0].hist(df['constraint_error'], bins=50, edgecolor='black', alpha=0.7, color='red')\n    axes[0].axvline(threshold, color='blue', linestyle='--', linewidth=2, \n                   label=f'Threshold: {threshold}g')\n    axes[0].set_title('Physical Constraint Error Distribution', fontsize=12, fontweight='bold')\n    axes[0].set_xlabel('Error (g)')\n    axes[0].set_ylabel('Frequency')\n    axes[0].legend()\n    axes[0].grid(True, alpha=0.3)\n    \n    # Dry_Total vs sum_components\n    axes[1].scatter(df['Dry_Total_g'], df['sum_components'], alpha=0.5, s=20)\n    axes[1].plot([df['Dry_Total_g'].min(), df['Dry_Total_g'].max()],\n                [df['Dry_Total_g'].min(), df['Dry_Total_g'].max()],\n                'r--', linewidth=2, label='Perfect Match')\n    axes[1].set_title('Dry_Total vs Component Sum', fontsize=12, fontweight='bold')\n    axes[1].set_xlabel('Dry_Total (g)')\n    axes[1].set_ylabel('Green + Dead + Clover (g)')\n    axes[1].legend()\n    axes[1].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n\ncheck_physical_constraints(df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:43.287903Z","iopub.execute_input":"2026-01-17T23:04:43.288235Z","iopub.status.idle":"2026-01-17T23:04:43.692886Z","shell.execute_reply.started":"2026-01-17T23:04:43.288211Z","shell.execute_reply":"2026-01-17T23:04:43.692145Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🔟 Outlier Detection\n\nWe'll use two methods:\n1. **IQR Method**: Q1 - 1.5×IQR to Q3 + 1.5×IQR\n2. **Z-Score Method**: |Z| > 3.0","metadata":{}},{"cell_type":"code","source":"def detect_outliers(df: pd.DataFrame, target: str) -> dict:\n    \"\"\"Detect outliers using IQR and Z-score methods\"\"\"\n    \n    print(f\"\\n🔍 Outlier Detection: {target}\")\n    print(\"=\"*60)\n    \n    # IQR method\n    Q1 = df[target].quantile(0.25)\n    Q3 = df[target].quantile(0.75)\n    IQR = Q3 - Q1\n    lower_bound = Q1 - 1.5 * IQR\n    upper_bound = Q3 + 1.5 * IQR\n    \n    iqr_outliers = df[(df[target] < lower_bound) | (df[target] > upper_bound)]\n    \n    print(f\"\\n【IQR Method】\")\n    print(f\"Bounds: [{lower_bound:.2f}g, {upper_bound:.2f}g]\")\n    print(f\"Outliers: {len(iqr_outliers):,} ({len(iqr_outliers)/len(df)*100:.1f}%)\")\n    \n    # Z-score method\n    mean = df[target].mean()\n    std = df[target].std()\n    z_scores = np.abs((df[target] - mean) / std)\n    zscore_outliers = df[z_scores > 3.0]\n    \n    print(f\"\\n【Z-Score Method】(threshold=3.0)\")\n    print(f\"Mean: {mean:.2f}g, Std: {std:.2f}g\")\n    print(f\"Outliers: {len(zscore_outliers):,} ({len(zscore_outliers)/len(df)*100:.1f}%)\")\n    \n    # Visualization\n    fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n    \n    # Box plot\n    bp = axes[0].boxplot(df[target], vert=True, patch_artist=True)\n    bp['boxes'][0].set_facecolor('lightblue')\n    axes[0].axhline(upper_bound, color='red', linestyle='--', linewidth=2, label='Upper')\n    axes[0].axhline(lower_bound, color='red', linestyle='--', linewidth=2, label='Lower')\n    axes[0].set_title(f'{target} - IQR Outlier Detection', fontsize=12, fontweight='bold')\n    axes[0].set_ylabel(f'{target} (g)')\n    axes[0].legend()\n    axes[0].grid(True, alpha=0.3)\n    \n    # Z-score plot\n    axes[1].scatter(range(len(df)), z_scores, alpha=0.5, s=20)\n    axes[1].axhline(3.0, color='red', linestyle='--', linewidth=2, label='Threshold: 3.0')\n    axes[1].set_title(f'{target} - Z-Score Plot', fontsize=12, fontweight='bold')\n    axes[1].set_xlabel('Sample Index')\n    axes[1].set_ylabel('|Z-Score|')\n    axes[1].legend()\n    axes[1].grid(True, alpha=0.3)\n    \n    plt.tight_layout()\n    plt.show()\n    \n    return {\n        'iqr_outliers': iqr_outliers,\n        'zscore_outliers': zscore_outliers\n    }\n\n# Detect outliers for all targets\noutlier_results = {}\nfor target in TARGET_COLS:\n    outlier_results[target] = detect_outliers(df, target)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:43.693976Z","iopub.execute_input":"2026-01-17T23:04:43.694309Z","iopub.status.idle":"2026-01-17T23:04:45.168067Z","shell.execute_reply.started":"2026-01-17T23:04:43.694276Z","shell.execute_reply":"2026-01-17T23:04:45.167399Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1️⃣1️⃣ 🖼️ Interactive Image Viewer\n\n**This is the most powerful feature of this notebook!**\n\nYou can:\n- Sort images by any metric (Dry_Total, NDVI, Height, etc.)\n- Filter by State, Season, Species\n- View detailed information for each image\n- Compare high/low biomass samples visually","metadata":{}},{"cell_type":"code","source":"def show_images_by_condition(\n    df: pd.DataFrame,\n    image_base_dir: str,\n    sort_by: str = 'Dry_Total_g',\n    ascending: bool = False,\n    n_images: int = 9,\n    filters: dict = None,\n    ncols: int = 3\n) -> None:\n    \"\"\"\n    Display images with filtering and sorting\n    \n    Args:\n        df: DataFrame with image info\n        image_base_dir: Base directory for images\n        sort_by: Column to sort by\n        ascending: Sort order\n        n_images: Number of images to display\n        filters: Dictionary of filters (e.g., {'State': 'Vic'})\n        ncols: Number of columns\n    \"\"\"\n    # Apply filters\n    df_filtered = df.copy()\n    if filters:\n        for col, value in filters.items():\n            if col in df_filtered.columns:\n                df_filtered = df_filtered[df_filtered[col] == value]\n                print(f\"✅ Filter: {col} == {value}\")\n    \n    print(f\"📊 Filtered samples: {len(df_filtered):,}\")\n    \n    if len(df_filtered) == 0:\n        print(\"❌ No samples match filters\")\n        return\n    \n    # Sort\n    df_sorted = df_filtered.sort_values(sort_by, ascending=ascending)\n    \n    # Display\n    nrows = (n_images + ncols - 1) // ncols\n    fig, axes = plt.subplots(nrows, ncols, figsize=(ncols*5, nrows*5))\n    axes = axes.flatten() if n_images > 1 else [axes]\n    \n    for idx, (i, row) in enumerate(df_sorted.head(n_images).iterrows()):\n        if idx >= len(axes):\n            break\n            \n        img_path = Path(image_base_dir) / row['image_path']\n        \n        if img_path.exists():\n            img = Image.open(img_path)\n            axes[idx].imshow(img)\n            axes[idx].axis('off')\n            \n            title = f\"ID: {row['image_id']}\\n{sort_by}: {row[sort_by]:.1f}g\\n\" + \\\n                   f\"NDVI: {row['Pre_GSHH_NDVI']:.2f} | Height: {row['Height_Ave_cm']:.1f}cm\\n\" + \\\n                   f\"{row['State']} | {row['season']}\"\n            axes[idx].set_title(title, fontsize=9)\n        else:\n            axes[idx].text(0.5, 0.5, 'Image not found', \n                         ha='center', va='center', color='red')\n            axes[idx].axis('off')\n    \n    for idx in range(n_images, len(axes)):\n        axes[idx].axis('off')\n    \n    plt.tight_layout()\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:45.169035Z","iopub.execute_input":"2026-01-17T23:04:45.169281Z","iopub.status.idle":"2026-01-17T23:04:45.178720Z","shell.execute_reply.started":"2026-01-17T23:04:45.169259Z","shell.execute_reply":"2026-01-17T23:04:45.178000Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def show_single_image_detail(\n    df: pd.DataFrame,\n    image_base_dir: str,\n    image_id: str\n) -> None:\n    \"\"\"Display single image with detailed info\"\"\"\n    \n    matching = df[df['image_id'] == image_id]\n    if len(matching) == 0:\n        print(f\"❌ Image ID '{image_id}' not found\")\n        return\n    \n    row = matching.iloc[0]\n    img_path = Path(image_base_dir) / row['image_path']\n    \n    fig = plt.figure(figsize=(16, 6))\n    gs = fig.add_gridspec(1, 2, width_ratios=[1.5, 1])\n    ax_img = fig.add_subplot(gs[0])\n    ax_info = fig.add_subplot(gs[1])\n    \n    # Image\n    if img_path.exists():\n        img = Image.open(img_path)\n        ax_img.imshow(img)\n        ax_img.axis('off')\n        ax_img.set_title(f\"Image ID: {row['image_id']}\", fontsize=14, fontweight='bold')\n    else:\n        ax_img.text(0.5, 0.5, 'Not found', ha='center', va='center', color='red')\n        ax_img.axis('off')\n    \n    # Info\n    ax_info.axis('off')\n    \n    info_text = f\"\"\"\nIMAGE INFORMATION\n{'='*50}\n\n【Metadata】\nDate: {row['Sampling_Date']}\nState: {row['State']}\nSpecies: {row['Species']}\nNDVI: {row['Pre_GSHH_NDVI']:.3f}\nHeight: {row['Height_Ave_cm']:.2f} cm\nSeason: {row['season']}\n\n【Target Values】\nDry_Total: {row['Dry_Total_g']:.2f} g (100%)\n  ├─ Dry_Green: {row['Dry_Green_g']:.2f} g ({row['green_ratio']*100:.1f}%)\n  ├─ Dry_Dead: {row['Dry_Dead_g']:.2f} g ({row['dead_ratio']*100:.1f}%)\n  └─ Dry_Clover: {row['Dry_Clover_g']:.2f} g ({row['clover_ratio']*100:.1f}%)\nGDM: {row['GDM_g']:.2f} g\n\n【Constraint Check】\nSum(G+D+C): {row['sum_components']:.2f} g\nError: {row['constraint_error']:.2f} g\n    \"\"\"\n    \n    ax_info.text(0.1, 0.5, info_text, fontsize=10, \n                verticalalignment='center', family='monospace',\n                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\n    \n    plt.tight_layout()\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:45.179711Z","iopub.execute_input":"2026-01-17T23:04:45.180068Z","iopub.status.idle":"2026-01-17T23:04:45.204307Z","shell.execute_reply.started":"2026-01-17T23:04:45.180035Z","shell.execute_reply":"2026-01-17T23:04:45.203741Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 📸 Image Viewer Examples","metadata":{}},{"cell_type":"code","source":"# Example 1: Highest Dry_Total images\nprint(\"=\"*80)\nprint(\"【Example 1】Top 9 Images by Dry_Total\")\nprint(\"=\"*80)\nshow_images_by_condition(\n    df, IMAGE_DIR,\n    sort_by='Dry_Total_g',\n    ascending=False,\n    n_images=9,\n    ncols=3\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:45.205063Z","iopub.execute_input":"2026-01-17T23:04:45.205326Z","iopub.status.idle":"2026-01-17T23:04:48.558743Z","shell.execute_reply.started":"2026-01-17T23:04:45.205295Z","shell.execute_reply":"2026-01-17T23:04:48.557555Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example 2: Lowest Dry_Total images\nprint(\"=\"*80)\nprint(\"【Example 2】Bottom 9 Images by Dry_Total\")\nprint(\"=\"*80)\nshow_images_by_condition(\n    df, IMAGE_DIR,\n    sort_by='Dry_Total_g',\n    ascending=True,\n    n_images=9,\n    ncols=3\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:48.560415Z","iopub.execute_input":"2026-01-17T23:04:48.560863Z","iopub.status.idle":"2026-01-17T23:04:51.245179Z","shell.execute_reply.started":"2026-01-17T23:04:48.560820Z","shell.execute_reply":"2026-01-17T23:04:51.243900Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example 3: High NDVI in Spring\nprint(\"=\"*80)\nprint(\"【Example 3】Top 6 Spring Images by NDVI\")\nprint(\"=\"*80)\nshow_images_by_condition(\n    df, IMAGE_DIR,\n    sort_by='Pre_GSHH_NDVI',\n    ascending=False,\n    n_images=6,\n    filters={'season': 'Spring'},\n    ncols=3\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:51.246665Z","iopub.execute_input":"2026-01-17T23:04:51.247106Z","iopub.status.idle":"2026-01-17T23:04:53.004358Z","shell.execute_reply.started":"2026-01-17T23:04:51.247061Z","shell.execute_reply":"2026-01-17T23:04:53.003384Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example 4: Victoria state, high height\nprint(\"=\"*80)\nprint(\"【Example 4】Top 6 Victoria Images by Height\")\nprint(\"=\"*80)\nshow_images_by_condition(\n    df, IMAGE_DIR,\n    sort_by='Height_Ave_cm',\n    ascending=False,\n    n_images=6,\n    filters={'State': 'Vic'},\n    ncols=3\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:53.005510Z","iopub.execute_input":"2026-01-17T23:04:53.005922Z","iopub.status.idle":"2026-01-17T23:04:54.899896Z","shell.execute_reply.started":"2026-01-17T23:04:53.005890Z","shell.execute_reply":"2026-01-17T23:04:54.898603Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example 5: Detailed view of specific image\nprint(\"=\"*80)\nprint(\"【Example 5】Detailed View\")\nprint(\"=\"*80)\nfirst_image_id = df.iloc[0]['image_id']\nshow_single_image_detail(df, IMAGE_DIR, first_image_id)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:54.901558Z","iopub.execute_input":"2026-01-17T23:04:54.902410Z","iopub.status.idle":"2026-01-17T23:04:55.506624Z","shell.execute_reply.started":"2026-01-17T23:04:54.902371Z","shell.execute_reply":"2026-01-17T23:04:55.505836Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example 6: Largest constraint violations\nprint(\"=\"*80)\nprint(\"【Example 6】Top 6 Constraint Violations\")\nprint(\"=\"*80)\nshow_images_by_condition(\n    df, IMAGE_DIR,\n    sort_by='constraint_error',\n    ascending=False,\n    n_images=6,\n    ncols=3\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:55.507656Z","iopub.execute_input":"2026-01-17T23:04:55.507906Z","iopub.status.idle":"2026-01-17T23:04:57.246201Z","shell.execute_reply.started":"2026-01-17T23:04:55.507882Z","shell.execute_reply":"2026-01-17T23:04:57.245180Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1️⃣2️⃣ Key Insights & Recommendations\n\n### 📊 Data Characteristics:\n1. **Target Distributions**\n   - All targets are right-skewed\n   - Consider log transformation: `log(1 + target)`\n   - Outliers present in all targets\n\n2. **Strong Correlations**\n   - High correlation between related targets\n   - Consider multi-task learning\n   - GDM = Dry_Green + Dry_Clover (by definition)\n\n3. **Physical Constraints**\n   - Dry_Total ≈ Green + Dead + Clover\n   - Can be used as regularization\n   - Small violations exist (<5g for most samples)\n\n### 🗺️ Metadata Patterns:\n1. **Geographic Variation**\n   - Different biomass levels across states\n   - Consider state as categorical feature\n\n2. **Temporal Patterns**\n   - Seasonal variation exists\n   - Year effects present\n   - Consider temporal features\n\n3. **Physical Measurements**\n   - NDVI shows positive correlation with biomass\n   - Height correlates with total biomass\n   - Both are valuable features\n\n### 🎯 Modeling Recommendations:\n\n1. **Target Transformation**\n   ```python\n   y_transformed = np.log1p(y)  # log(1 + y)\n   ```\n\n2. **Cross-Validation Strategy**\n   - Stratify by State or Season\n   - GroupKFold by date/location\n   - Watch for temporal leakage\n\n3. **Feature Engineering**\n   - Use NDVI, Height as features\n   - One-hot encode State\n   - Cyclical encoding for month/season\n   - Consider species embeddings\n\n4. **Model Architecture**\n   - Multi-head architecture (one head per target)\n   - Shared backbone with task-specific heads\n   - Add physical constraint loss:\n     ```python\n     constraint_loss = |Dry_Total - (Green + Dead + Clover)|\n     ```\n\n5. **Data Augmentation**\n   - Standard image augmentations\n   - Be careful with spatial augmentations\n   - Consider Mixup/CutMix\n\n### ⚠️ Potential Challenges:\n1. Small dataset (357 images)\n2. Class imbalance in some species\n3. Outliers need handling\n4. Test set may have different distribution\n\n---\n\n**Good luck with your modeling! 🚀**","metadata":{}},{"cell_type":"markdown","source":"# **Modelling**","metadata":{}},{"cell_type":"code","source":"# 1. CONFIGURATION & SEEDING\n# =========================================================================================\n@dataclass\nclass Config:\n    DATA_PATH: Path = Path(\"/kaggle/input/csiro-biomass/\")\n    SPLIT_PATH: Path = Path(\"/kaggle/input/csiro-datasplit/csiro_data_split.csv\")\n    SIGLIP_PATH: str = \"/kaggle/input/google-siglip-so400m-patch14-384/transformers/default/1\"\n    \n    SEED: int = 42\n    DEVICE: str = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n    PATCH_SIZE: int = 520\n    OVERLAP: int = 16\n    \n    # Target definitions\n    TARGET_NAMES = ['Dry_Clover_g', 'Dry_Dead_g', 'Dry_Green_g', 'Dry_Total_g', 'GDM_g']\n    TARGET_MAX = {\n        \"Dry_Clover_g\": 71.7865,\n        \"Dry_Dead_g\": 83.8407,\n        \"Dry_Green_g\": 157.9836,\n        \"Dry_Total_g\": 185.70,\n        \"GDM_g\": 157.9836,\n    }\n\ncfg = Config()\n\ndef seeding(SEED):\n    np.random.seed(SEED)\n    random.seed(SEED)\n    os.environ['PYTHONHASHSEED'] = str(SEED)\n    torch.manual_seed(SEED)\n    if torch.cuda.is_available(): \n        torch.cuda.manual_seed(SEED)\n        torch.backends.cudnn.deterministic = True\n        torch.backends.cudnn.benchmark = False\n\nseeding(cfg.SEED)\n\n# =========================================================================================\n# 2. DATA LOADING & PRE-PROCESSING\n# =========================================================================================\ndef pivot_table(df: pd.DataFrame) -> pd.DataFrame:\n    if 'target' in df.columns.tolist():\n        # Train data\n        df_pt = pd.pivot_table(\n            df, \n            values='target', \n            index=['image_path', 'Sampling_Date', 'State', 'Species', 'Pre_GSHH_NDVI', 'Height_Ave_cm'], \n            columns='target_name', \n            aggfunc='mean'\n        ).reset_index()\n    else:\n        # Test data\n        df['target'] = 0\n        df_pt = pd.pivot_table(\n            df, \n            values='target', \n            index='image_path', \n            columns='target_name', \n            aggfunc='mean'\n        ).reset_index()\n    return df_pt\n\ndef melt_table(df: pd.DataFrame) -> pd.DataFrame:\n    melted = df.melt(\n        id_vars='image_path',\n        value_vars=cfg.TARGET_NAMES,\n        var_name='target_name',\n        value_name='target'\n    )\n    # Create sample_id matching submission format\n    melted['sample_id'] = (\n        melted['image_path']\n        .str.replace(r'^.*/', '', regex=True)\n        .str.replace('.jpg', '', regex=False)\n        + '__' + melted['target_name']\n    )\n    return melted[['sample_id', 'target']]\n\ndef post_process_biomass(df_preds):\n    \"\"\"\n    Derive GDM_g and Dry_Total_g from primary predictions.\n    IMPORTANT: Keeps Dry_Clover_g fixed at 0.0 and does NOT modify Dry_Green_g or Dry_Dead_g.\n    \"\"\"\n    ordered_cols = [\"Dry_Green_g\", \"Dry_Clover_g\", \"Dry_Dead_g\", \"GDM_g\", \"Dry_Total_g\"]\n    \n    # Ensure cols exist\n    for c in ordered_cols:\n        if c not in df_preds.columns:\n            df_preds[c] = 0.0\n    \n    df_out = df_preds.copy()\n    \n    # Keep Dry_Clover_g fixed at 0.0 (should already be 0.0 from model prediction)\n    df_out['Dry_Clover_g'] = 0.0\n    \n    # Keep Dry_Green_g and Dry_Dead_g as predicted (NO MODIFICATION)\n    # Only derive GDM_g and Dry_Total_g\n    df_out['GDM_g'] = df_out['Dry_Green_g'] + df_out['Dry_Clover_g']\n    df_out['Dry_Total_g'] = df_out['GDM_g'] + df_out['Dry_Dead_g']\n    \n    # Clip only derived targets to non-negative\n    df_out['GDM_g'] = df_out['GDM_g'].clip(lower=0.0)\n    df_out['Dry_Total_g'] = df_out['Dry_Total_g'].clip(lower=0.0)\n    \n    return df_out\n\nprint(\"Loading Data...\")\n# Load Train (Metadata Split)\ntrain_df = pd.read_csv(cfg.SPLIT_PATH) \n\n# --- FIX: Remove pre-existing embedding columns to prevent duplication ---\ncols_to_keep = [c for c in train_df.columns if not c.startswith('emb')]\ntrain_df = train_df[cols_to_keep]\n# -----------------------------------------------------------------------\n\n# Ensure train paths match local environment\nif not str(train_df['image_path'].iloc[0]).startswith('/'):\n     train_df['image_path'] = train_df['image_path'].apply(lambda p: str(cfg.DATA_PATH / 'train' / os.path.basename(p)))\n\n# Load Test\ntest_df_raw = pd.read_csv(cfg.DATA_PATH / 'test.csv')\ntest_df = pivot_table(test_df_raw)\ntest_df['image_path'] = test_df['image_path'].apply(lambda p: str(cfg.DATA_PATH / p))\n\n# =========================================================================================\n# 3. FEATURE EXTRACTION: SIGLIP IMAGE EMBEDDINGS\n# =========================================================================================\ndef split_image(image, patch_size=520, overlap=16):\n    h, w, c = image.shape\n    stride = patch_size - overlap\n    patches = []\n    for y in range(0, h, stride):\n        for x in range(0, w, stride):\n            y2 = min(y + patch_size, h)\n            x2 = min(x + patch_size, w)\n            y1 = max(0, y2 - patch_size) # Ensure fixed size\n            x1 = max(0, x2 - patch_size)\n            patch = image[y1:y2, x1:x2, :]\n            patches.append(patch)\n    return patches\n\ndef compute_embeddings(model_path, df):\n    print(f\"Computing Embeddings for {len(df)} images...\")\n    model = AutoModel.from_pretrained(model_path, local_files_only=True).eval().to(cfg.DEVICE)\n    processor = AutoImageProcessor.from_pretrained(model_path)\n    \n    EMBEDDINGS = []\n    \n    for _, row in tqdm(df.iterrows(), total=len(df)):\n        try:\n            img = cv2.imread(row['image_path'])\n            if img is None: raise ValueError(\"Image not found\")\n            img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n            patches = split_image(img, patch_size=cfg.PATCH_SIZE, overlap=cfg.OVERLAP)\n            images = [Image.fromarray(p) for p in patches]\n            \n            # Batch process patches\n            inputs = processor(images=images, return_tensors=\"pt\").to(cfg.DEVICE)\n            with torch.no_grad():\n                features = model.get_image_features(**inputs)\n            \n            # Average pooling of patches\n            avg_embed = features.mean(dim=0).cpu().numpy()\n            EMBEDDINGS.append(avg_embed)\n        except Exception as e:\n            print(f\"Error processing {row['image_path']}: {e}\")\n            # Fallback zero embedding\n            EMBEDDINGS.append(np.zeros(1152))\n        \n    torch.cuda.empty_cache()\n    return np.stack(EMBEDDINGS)\n\n# Compute Features\ntrain_embeddings = compute_embeddings(cfg.SIGLIP_PATH, train_df)\ntest_embeddings = compute_embeddings(cfg.SIGLIP_PATH, test_df)\n\n# Create Feature DataFrames\nemb_cols = [f\"emb{i}\" for i in range(train_embeddings.shape[1])]\ntrain_feat_df = pd.concat([train_df, pd.DataFrame(train_embeddings, columns=emb_cols)], axis=1)\ntest_feat_df = pd.concat([test_df, pd.DataFrame(test_embeddings, columns=emb_cols)], axis=1)\n\n# Double check column counts\nprint(f\"Train Features Shape: {train_feat_df.shape}\")\nprint(f\"Test Features Shape: {test_feat_df.shape}\")\n\n# =========================================================================================\n# 4. FEATURE EXTRACTION: SEMANTIC FEATURES (TEXT PROBING)\n# =========================================================================================\ndef generate_semantic_features(image_embeddings_np, model_path):\n    print(\"Generating Semantic Features...\")\n    model = AutoModel.from_pretrained(model_path).to(cfg.DEVICE)\n    tokenizer = AutoTokenizer.from_pretrained(model_path)\n    \n    # Anchors\n    concept_groups = {\n        \"bare\": [\"bare soil\", \"dirt ground\", \"sparse vegetation\", \"exposed earth\"],\n        \"sparse\": [\"low density pasture\", \"thin grass\", \"short clipped grass\"],\n        \"medium\": [\"average pasture cover\", \"medium height grass\", \"grazed pasture\"],\n        \"dense\": [\"dense tall pasture\", \"thick grassy volume\", \"high biomass\", \"overgrown vegetation\"],\n        \"green\": [\"lush green vibrant pasture\", \"photosynthesizing leaves\", \"fresh growth\"],\n        \"dead\": [\"dry brown dead grass\", \"yellow straw\", \"senesced material\", \"standing hay\"],\n        \"clover\": [\"white clover\", \"trifolium repens\", \"broadleaf legume\", \"clover flowers\"],\n        \"grass\": [\"ryegrass\", \"blade-like leaves\", \"fescue\", \"grassy sward\"]\n    }\n    \n    # Encode Concepts\n    concept_vectors = {}\n    with torch.no_grad():\n        for name, prompts in concept_groups.items():\n            inputs = tokenizer(prompts, padding=\"max_length\", return_tensors=\"pt\").to(cfg.DEVICE)\n            emb = model.get_text_features(**inputs)\n            emb = emb / emb.norm(p=2, dim=-1, keepdim=True)\n            concept_vectors[name] = emb.mean(dim=0, keepdim=True)\n            \n    # Compute Scores\n    img_tensor = torch.tensor(image_embeddings_np, dtype=torch.float32).to(cfg.DEVICE)\n    img_tensor = img_tensor / img_tensor.norm(p=2, dim=-1, keepdim=True)\n    \n    scores = {}\n    for name, vec in concept_vectors.items():\n        scores[name] = torch.matmul(img_tensor, vec.T).cpu().numpy().flatten()\n        \n    df_scores = pd.DataFrame(scores)\n    # Ratios\n    df_scores['ratio_greenness'] = df_scores['green'] / (df_scores['green'] + df_scores['dead'] + 1e-6)\n    df_scores['ratio_clover'] = df_scores['clover'] / (df_scores['clover'] + df_scores['grass'] + 1e-6)\n    df_scores['ratio_cover'] = (df_scores['dense'] + df_scores['medium']) / (df_scores['bare'] + df_scores['sparse'] + 1e-6)\n    \n    torch.cuda.empty_cache()\n    return df_scores.values\n\n# Combine for semantic generation to ensure consistency\nall_emb = np.vstack([train_embeddings, test_embeddings])\nall_semantic = generate_semantic_features(all_emb, cfg.SIGLIP_PATH)\n\nsem_train = all_semantic[:len(train_df)]\nsem_test = all_semantic[len(train_df):]\n\n# =========================================================================================\n# 5. SUPERVISED EMBEDDING ENGINE\n# =========================================================================================\nclass SupervisedEmbeddingEngine:\n    def __init__(self, n_pca=0.80, n_pls=8, n_gmm=6, random_state=42):\n        self.scaler = StandardScaler()\n        self.pca = PCA(n_components=n_pca, random_state=random_state)\n        self.pls = PLSRegression(n_components=n_pls, scale=False)\n        self.gmm = GaussianMixture(n_components=n_gmm, covariance_type='diag', random_state=random_state)\n        self.pls_fitted_ = False\n\n    def fit(self, X, y=None, X_semantic=None):\n        # Scale\n        X_scaled = self.scaler.fit_transform(X)\n        \n        # Unsupervised\n        self.pca.fit(X_scaled)\n        self.gmm.fit(X_scaled)\n        \n        # Supervised\n        if y is not None:\n            self.pls.fit(X_scaled, y)\n            self.pls_fitted_ = True\n        return self\n\n    def transform(self, X, X_semantic=None):\n        X_scaled = self.scaler.transform(X)\n        \n        features = [self.pca.transform(X_scaled)]\n        \n        if self.pls_fitted_:\n            features.append(self.pls.transform(X_scaled))\n            \n        features.append(self.gmm.predict_proba(X_scaled))\n        \n        if X_semantic is not None:\n            # Normalize semantic\n            sem_norm = (X_semantic - np.mean(X_semantic, axis=0)) / (np.std(X_semantic, axis=0) + 1e-6)\n            features.append(sem_norm)\n            \n        return np.hstack(features)\n\n# =========================================================================================\n# 6. TRAINING & INFERENCE (5-FOLD CV)\n# =========================================================================================\ndef cross_validate_predict(model_cls, model_params, train_data, test_data, sem_tr, sem_te, feature_engine):\n    target_max_arr = np.array([cfg.TARGET_MAX[t] for t in cfg.TARGET_NAMES], dtype=float)\n    y_pred_test_accum = np.zeros([len(test_data), len(cfg.TARGET_NAMES)], dtype=float)\n    \n    # Ensure n_splits is integer\n    n_splits = int(train_data['fold'].nunique())\n    \n    # Pre-extract raw columns to avoid indexing overhead\n    # Force float32 to save memory and ensure compatibility\n    X_train_full = train_data[emb_cols].values.astype(np.float32)\n    X_test_raw = test_data[emb_cols].values.astype(np.float32)\n    y_train_full = train_data[cfg.TARGET_NAMES].values.astype(np.float32)\n    \n    for fold in range(n_splits):\n        print(f\"Processing Fold {fold}...\")\n        # Split\n        train_mask = train_data['fold'] != fold\n        \n        X_tr = X_train_full[train_mask]\n        y_tr = y_train_full[train_mask] / target_max_arr # Max Scaling\n        \n        sem_tr_fold = sem_tr[train_mask]\n        \n        # Feature Engineering (Fit on fold train)\n        engine = deepcopy(feature_engine)\n        engine.fit(X_tr, y=y_tr, X_semantic=sem_tr_fold)\n        \n        x_tr_eng = engine.transform(X_tr, X_semantic=sem_tr_fold)\n        x_te_eng = engine.transform(X_test_raw, X_semantic=sem_te)\n        \n        # Train & Predict per target\n        fold_test_pred = np.zeros([len(test_data), len(cfg.TARGET_NAMES)])\n        \n        for k in range(len(cfg.TARGET_NAMES)):\n            target_name = cfg.TARGET_NAMES[k]\n            \n            # Dry_Clover_g: Always predict 0.0 (no training needed)\n            if target_name == 'Dry_Clover_g':\n                fold_test_pred[:, k] = 0.0\n            else:\n                # Train and predict for other targets\n                model = model_cls(**model_params)\n                model.fit(x_tr_eng, y_tr[:, k])\n                pred_raw = model.predict(x_te_eng)\n                fold_test_pred[:, k] = pred_raw * target_max_arr[k] # Inverse Scale\n            \n        y_pred_test_accum += fold_test_pred\n        \n    return y_pred_test_accum / n_splits\n\n# Model Parameters (Optimized)\nparams_cat = {\n    'iterations': 1900, 'learning_rate': 0.045, 'depth': 4, 'l2_leaf_reg': 0.56, \n    'random_strength': 0.045, 'bagging_temperature': 0.98, 'verbose': 0, 'random_state': 42,\n    'allow_writing_files': False\n}\nparams_xgb = { # Using GradientBoostingRegressor as proxy\n    'n_estimators': 1354, 'learning_rate': 0.010, 'max_depth': 3, 'subsample': 0.60, \n    'random_state': 42\n}\nparams_lgbm = {\n    'n_estimators': 807, 'learning_rate': 0.014, 'num_leaves': 48, 'min_child_samples': 19, \n    'subsample': 0.745, 'colsample_bytree': 0.745, 'reg_alpha': 0.21, 'reg_lambda': 3.78,\n    'verbose': -1, 'random_state': 42\n}\nparams_hist = {\n    'max_iter': 300, 'learning_rate': 0.05, 'max_depth': None, 'l2_regularization': 0.44,\n    'random_state': 42\n}\n\nfeat_engine = SupervisedEmbeddingEngine(n_pca=0.80, n_pls=8, n_gmm=6)\n\nprint(\"Training & Inferring Models...\")\n\n# 1. HistGradientBoosting\nprint(\"Model: HistGradientBoosting\")\npred_hist = cross_validate_predict(\n    HistGradientBoostingRegressor, params_hist, \n    train_feat_df, test_feat_df, sem_train, sem_test, feat_engine\n)\n\n# 2. GradientBoosting\nprint(\"Model: GradientBoosting\")\npred_gb = cross_validate_predict(\n    GradientBoostingRegressor, params_xgb, \n    train_feat_df, test_feat_df, sem_train, sem_test, feat_engine\n)\n\n# 3. CatBoost\nprint(\"Model: CatBoost\")\npred_cat = cross_validate_predict(\n    CatBoostRegressor, params_cat, \n    train_feat_df, test_feat_df, sem_train, sem_test, feat_engine\n)\n\n# 4. LightGBM\nprint(\"Model: LightGBM\")\npred_lgbm = cross_validate_predict(\n    LGBMRegressor, params_lgbm, \n    train_feat_df, test_feat_df, sem_train, sem_test, feat_engine\n)\n\n# =========================================================================================\n# 7. ENSEMBLING & SUBMISSION\n# =========================================================================================\nprint(\"Ensembling and Post-processing...\")\n# Simple Average Ensemble\nfinal_pred = (pred_hist + pred_gb + pred_cat + pred_lgbm) / 4.0\n\n# Assign to dataframe\ntest_feat_df[cfg.TARGET_NAMES] = final_pred\n\n# Post-process (Mass Balance)\ntest_processed = post_process_biomass(test_feat_df)\n\n# Create Submission File\nsub_df = melt_table(test_processed)\noutput_path = \"submission_siglip.csv\"\nsub_df.to_csv(output_path, index=False)\n\nprint(f\"✓ Siglip submission generated: {output_path}\")\nprint(sub_df.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:04:57.247577Z","iopub.execute_input":"2026-01-17T23:04:57.247843Z","iopub.status.idle":"2026-01-17T23:19:15.525126Z","shell.execute_reply.started":"2026-01-17T23:04:57.247817Z","shell.execute_reply":"2026-01-17T23:19:15.524434Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile csiro_infer.py\n\nfrom PIL import Image\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import KFold, GroupKFold, StratifiedGroupKFold\nfrom tqdm.auto import tqdm\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.optim import AdamW\nfrom transformers import get_cosine_schedule_with_warmup\n\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision import transforms as T\n\nimport os\nfrom pathlib import Path\nimport timm\nimport warnings \n\nwarnings.filterwarnings('ignore')\ntqdm.pandas()\n\nclass RegressionDataset(Dataset):\n    def __init__(self, data, transform=None):\n        self.data = data\n        self.transform = transform\n    def __len__(self):\n        return self.data.shape[0]\n\n    def __getitem__(self, idx):\n        item = self.data.iloc[idx]\n        image = item.image\n        targets = [item['Dry_Green_g'], item['Dry_Clover_g'], item['Dry_Dead_g']]\n        width, height = image.size\n        mid_point = width // 2\n        left_image = image.crop((0, 0, mid_point, height))\n        right_image = image.crop((mid_point, 0, width, height))\n\n        if self.transform is not None:\n            left_image = self.transform(left_image)\n            right_image = self.transform(right_image)\n\n        return left_image, right_image, targets\n\ndef get_test_dataloaders(data, image_size, batch_size):\n    res = []\n    for trans in [None, T.RandomHorizontalFlip(p=1.0), T.RandomVerticalFlip(p=1.0)]:\n        transform = T.Compose([\n            T.Resize(image_size),\n            T.ToTensor(),\n            T.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])\n        ])\n        if trans:\n            transform = T.Compose([\n                T.Resize(image_size),\n                trans,\n                T.ToTensor(),\n                T.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])\n            ])\n        dataset = RegressionDataset(data, transform=transform)\n        res.append(DataLoader(dataset, batch_size=batch_size, shuffle=False, num_workers=4))\n    return res\n\n\nclass FiLM(nn.Module):\n    def __init__(self, feat_dim):\n        super().__init__()\n        self.mlp = nn.Sequential(\n            nn.Linear(feat_dim, feat_dim // 2), \n            nn.ReLU(inplace=True), \n            nn.Linear(feat_dim // 2, feat_dim * 2)\n        )\n\n    def forward(self, context):\n        gamma_beta = self.mlp(context)\n        return torch.chunk(gamma_beta, 2, dim=1)\n\nclass CSIROModelRegressor(nn.Module):\n    def __init__(self, model_name, pretrained=True, num_classes=3, dropout=0.0, freeze_backbone=False):\n        super().__init__()\n        self.backbone = timm.create_model(model_name, pretrained=pretrained, num_classes=0, global_pool='avg')\n\n        self.film = FiLM(self.backbone.num_features)\n\n        self.dropout = nn.Dropout(dropout)\n\n        def make_head():\n            return nn.Sequential(\n                nn.Linear(self.backbone.num_features * 2, 8),\n                nn.ReLU(inplace=True),\n                nn.Dropout(dropout),\n                nn.Linear(8, 1),\n            )\n\n        self.head_green = make_head()\n        self.head_clover = make_head()\n        self.head_dead = make_head()\n\n        self.softplus = nn.Softplus(beta=1.0)\n\n        if freeze_backbone:\n            for p in self.backbone.parameters():\n                p.requires_grad = False\n\n\n    def forward(self, left_img, right_img):\n        left_feat = self.backbone(left_img)\n        right_feat = self.backbone(right_img)\n\n        context = (left_feat + right_feat) / 2\n        gamma, beta = self.film(context)\n\n        left_feat_modulated = left_feat * (1 + gamma) + beta\n        right_feat_modulated = right_feat * (1 + gamma) + beta\n\n        combined = torch.cat([left_feat_modulated, right_feat_modulated], dim=1)\n\n        green = self.softplus(self.head_green(combined))   \n        clover = self.softplus(self.head_clover(combined))\n        dead = self.softplus(self.head_dead(combined)) \n\n        logits = torch.cat([green, clover, dead], dim=1)\n\n        return logits\n\ndef predict(model, dataloader, device):\n    model.to(device)\n    model.eval()\n\n    all_outputs = []\n    with torch.no_grad():\n        for left_images, right_images, targets in dataloader:\n            left_images = left_images.to(device)\n            right_images = right_images.to(device)\n\n            outputs = model(left_images, right_images)\n            all_outputs.append(outputs.detach().cpu())\n\n    outputs = torch.cat(all_outputs).numpy()\n    return outputs\n\ndef predict_loaders(model, dataloaders, device):\n    all_outputs = []\n    for dataloader in dataloaders:\n        outputs = predict(model, dataloader, device)\n        all_outputs.append(outputs)\n    avg_outputs = np.mean(all_outputs, axis=0)\n    return avg_outputs\n\ndef predict_folds(dataloaders,models_dir):\n    device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')\n\n    all_preds = []\n    for model_file in Path(models_dir).glob('*.pth'):\n        model = CSIROModelRegressor(CFG.MODEL_NAME, pretrained=False, num_classes=3)\n        model.load_state_dict(torch.load(model_file))\n        preds = predict_loaders(model, dataloaders, device)\n        all_preds.append(preds)\n\n    avg_preds = np.mean(all_preds, axis=0)\n    return avg_preds\n\n\nclass CFG:\n    DATA_PATH=\"/kaggle/input/csiro-biomass/\"\n    TEST_DATA_PATH=\"/kaggle/input/csiro-biomass/test.csv\"\n    MODEL_NAME=\"vit_large_patch16_dinov3_qkvb\"\n    MODELS_DIR ='/kaggle/input/modelv3/pytorch/default/1/models_retrained'\n    IMG_SIZE=(512,512)\n\n\ntest_df = pd.read_csv(CFG.TEST_DATA_PATH)\n\ntest_df['target'] = 0.0\ntest_df[['sample_id_prefix', 'sample_id_suffix']] = test_df.sample_id.str.split('__', expand=True)\n\ntest_data_df = test_df.groupby(['sample_id_prefix', 'image_path']).apply(lambda df: df.set_index('target_name').target)\ntest_data_df.reset_index(inplace=True)\ntest_data_df.columns.name = None\n\ntest_data_df['image'] = test_data_df.image_path.progress_apply(lambda path: Image.open(CFG.DATA_PATH + path).convert('RGB'))\n\ntest_loaders = get_test_dataloaders(test_data_df, CFG.IMG_SIZE, 32)\npreds = predict_folds(test_loaders,models_dir=CFG.MODELS_DIR)\n\n\n\ntest_data_df[['Dry_Green_g', 'Dry_Clover_g', 'Dry_Dead_g']] = preds\ntest_data_df['Dry_Green_g'] = test_data_df['Dry_Green_g'].apply(lambda x: 0.0 if x < 0.09 else x)\ntest_data_df['Dry_Clover_g'] = test_data_df['Dry_Clover_g'].apply(lambda x: 0.0 if x < 0.09 else x)\ntest_data_df['Dry_Dead_g'] = test_data_df['Dry_Dead_g'].apply(lambda x: 0.0 if x < 0.09 else x)\ntest_data_df['GDM_g'] = test_data_df.Dry_Green_g + test_data_df.Dry_Clover_g\ntest_data_df['Dry_Total_g'] = test_data_df.GDM_g + test_data_df.Dry_Dead_g\n\ncols = [ 'Dry_Green_g', 'Dry_Dead_g', 'Dry_Clover_g', 'GDM_g' , 'Dry_Total_g']\nsub_df = test_data_df.set_index('sample_id_prefix')[cols].stack()\nsub_df = sub_df.reset_index()\nsub_df.columns = ['sample_id_prefix', 'target_name', 'target']\n\nsub_df['sample_id'] = sub_df.sample_id_prefix + '__' + sub_df.target_name\n\ncols = ['sample_id', 'target']\nsub_df[cols].to_csv('submission_dinov2026.csv', index=False)\n\nprint(sub_df[cols])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:19:15.528301Z","iopub.execute_input":"2026-01-17T23:19:15.528554Z","iopub.status.idle":"2026-01-17T23:19:15.536803Z","shell.execute_reply.started":"2026-01-17T23:19:15.528530Z","shell.execute_reply":"2026-01-17T23:19:15.536141Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!python csiro_infer.py","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:19:15.537719Z","iopub.execute_input":"2026-01-17T23:19:15.538364Z","iopub.status.idle":"2026-01-17T23:21:03.594966Z","shell.execute_reply.started":"2026-01-17T23:19:15.538342Z","shell.execute_reply":"2026-01-17T23:21:03.594209Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# CONFIGURATION\n# =============================================================================\n# Weights for the ensemble (Must sum to 1.0)\n# Adjust based on which model had better local CV or Public LB score.\n# Example: If Siglip was better, give it 0.6.\nW_SIGLIP = 0.35\nW_DINO   = 0.65\n\nFILES = {\n    'siglip': 'submission_siglip.csv',\n    'dino':   'submission_dinov2026.csv'\n}\n\nOUTPUT_FILE = 'submission.csv'\n\n# Target definitions required for Mass Balance\nALL_TARGETS = ['Dry_Green_g', 'Dry_Clover_g', 'Dry_Dead_g', 'GDM_g', 'Dry_Total_g']\n\n# =============================================================================\n# HELPER FUNCTIONS\n# =============================================================================\n\ndef enforce_mass_balance(df_wide, fixed_clover=None):\n    \"\"\"\n    Applies Orthogonal Projection to enforce biological constraints:\n    1. Dry_Green_g + Dry_Clover_g = GDM_g\n    2. GDM_g + Dry_Dead_g = Dry_Total_g\n    \n    If fixed_clover is True, Dry_Clover_g is kept fixed and only other targets are adjusted.\n    \n    This finds the closest set of values to the predictions that satisfy \n    the constraints (minimizing Euclidean distance modification).\n    \"\"\"\n    # 1. Ensure columns are in the specific order for the matrix math\n    # Vector x = [Green, Clover, Dead, GDM, Total]\n    ordered_cols = ['Dry_Green_g', 'Dry_Clover_g', 'Dry_Dead_g', 'GDM_g', 'Dry_Total_g']\n    \n    # Extract values: Shape (5, N_samples)\n    Y = df_wide[ordered_cols].values.T\n    \n    if fixed_clover:\n        # Keep Dry_Clover_g fixed, adjust only other targets\n        # We have: Green + Clover_fixed = GDM, so GDM = Green + Clover_fixed\n        # And: GDM + Dead = Total, so Total = GDM + Dead = Green + Clover_fixed + Dead\n        # Extract fixed Clover values\n        clover_fixed = Y[1, :].copy()  # Dry_Clover_g is index 1\n        # Adjust: GDM = Green + Clover_fixed\n        Y[3, :] = Y[0, :] + clover_fixed  # GDM_g\n        # Adjust: Total = GDM + Dead = Green + Clover_fixed + Dead\n        Y[4, :] = Y[3, :] + Y[2, :]  # Dry_Total_g\n        # Keep Green and Dead as is (or apply minimal adjustment if needed)\n        # Clover stays fixed (already set)\n        Y_reconciled = Y\n    else:\n        # Original method: adjust all targets\n        # 2. Define Constraint Matrix C where Cx = 0\n        # Eq 1: 1*Gr + 1*Cl + 0*De - 1*GDM + 0*Tot = 0\n        # Eq 2: 0*Gr + 0*Cl + 1*De + 1*GDM - 1*Tot = 0\n        C = np.array([\n            [1, 1, 0, -1,  0],\n            [0, 0, 1,  1, -1]\n        ])\n        \n        # 3. Calculate Projection Matrix P = I - C^T * (C * C^T)^-1 * C\n        C_T = C.T\n        try:\n            inv_CCt = np.linalg.inv(C @ C_T)\n            P = np.eye(5) - C_T @ inv_CCt @ C\n        except np.linalg.LinAlgError:\n            # Fallback if singular (unlikely with this specific matrix)\n            print(\"Warning: Singular matrix in projection. Skipping constraint enforcement.\")\n            return df_wide\n\n        # 4. Apply Projection\n        Y_reconciled = P @ Y\n    \n    # 5. Transpose back to (N_samples, 5) and clip negatives\n    Y_reconciled = Y_reconciled.T\n    Y_reconciled = np.maximum(0, Y_reconciled) \n    \n    # 6. Update DataFrame\n    df_out = df_wide.copy()\n    df_out[ordered_cols] = Y_reconciled\n    \n    return df_out\n\ndef robust_ensemble(file_paths, weights):\n    print(f\"--- Starting Ensemble ---\")\n    print(f\"Weights: {weights}\")\n    print(\"NOTE: Using DINO-only for Dry_Clover_g (better detection)\")\n    \n    dfs = []\n    for name, path in file_paths.items():\n        if not os.path.exists(path):\n            raise FileNotFoundError(f\"Missing file: {path}\")\n        \n        # Read and sort by sample_id to ensure alignment\n        df = pd.read_csv(path).sort_values('sample_id').reset_index(drop=True)\n        dfs.append(df)\n        print(f\"Loaded {name}: {len(df)} rows\")\n\n    # 1. Check alignment\n    base_ids = dfs[0]['sample_id']\n    if not all(df['sample_id'].equals(base_ids) for df in dfs[1:]):\n        raise ValueError(\"Sample IDs do not match between submission files!\")\n\n    # 2. Split by target_name to handle Dry_Clover_g separately\n    # Split sample_id into image_id and target_name\n    dfs_split = []\n    for df in dfs:\n        df_split = df.copy()\n        df_split[['image_id', 'target_name']] = df_split['sample_id'].str.rsplit('__', n=1, expand=True)\n        dfs_split.append(df_split)\n    \n    # Get DINO and SigLIP dataframes (identify by order in file_paths)\n    dino_idx = list(file_paths.keys()).index('dino')\n    siglip_idx = list(file_paths.keys()).index('siglip')\n    dino_df = dfs_split[dino_idx]\n    siglip_df = dfs_split[siglip_idx]\n    \n    # 3. Separate targets: Dry_Clover_g uses DINO only, others use ensemble\n    ensemble_results = []\n    \n    for target in ALL_TARGETS:\n        # Filter for this target\n        dino_target = dino_df[dino_df['target_name'] == target].copy()\n        siglip_target = siglip_df[siglip_df['target_name'] == target].copy()\n        \n        if target == 'Dry_Clover_g':\n            # Use DINO only for Clover\n            print(f\"Using DINO-only for {target}\")\n            ensemble_results.append(dino_target[['sample_id', 'target']])\n        else:\n            # Ensemble for other targets\n            # Merge on sample_id to align\n            merged = pd.merge(\n                dino_target[['sample_id', 'target']].rename(columns={'target': 'dino_pred'}),\n                siglip_target[['sample_id', 'target']].rename(columns={'target': 'siglip_pred'}),\n                on='sample_id'\n            )\n            \n            # Weighted average\n            w_vec = np.array([weights['dino'], weights['siglip']])\n            w_vec = w_vec / w_vec.sum()\n            \n            merged['target'] = merged['dino_pred'] * w_vec[0] + merged['siglip_pred'] * w_vec[1]\n            ensemble_results.append(merged[['sample_id', 'target']])\n    \n    # 4. Combine all targets\n    ensemble_df = pd.concat(ensemble_results, ignore_index=True)\n    print(\"Weighted average complete (DINO-only for Clover).\")\n\n    # 5. Prepare for Mass Balance (Convert Long -> Wide)\n    ensemble_df[['image_id', 'target_name']] = ensemble_df['sample_id'].str.rsplit('__', n=1, expand=True)\n    \n    # Pivot\n    wide_df = ensemble_df.pivot(index='image_id', columns='target_name', values='target').reset_index()\n\n    # 6. Apply Robust Constraints (with Dry_Clover_g fixed)\n    print(\"Applying Mass Balance Constraints (Dry_Clover_g fixed to DINO values)...\")\n    wide_balanced = enforce_mass_balance(wide_df, fixed_clover=True)\n\n    # 7. Convert back (Wide -> Long)\n    long_balanced = wide_balanced.melt(\n        id_vars='image_id', \n        value_vars=ALL_TARGETS,\n        var_name='target_name',\n        value_name='target'\n    )\n\n    # Reconstruct sample_id\n    long_balanced['sample_id'] = long_balanced['image_id'] + '__' + long_balanced['target_name']\n\n    # 8. Final Formatting\n    final_submission = long_balanced[['sample_id', 'target']].sort_values('sample_id').reset_index(drop=True)\n\n    return final_submission\n\n# =============================================================================\n# EXECUTION\n# =============================================================================\nif __name__ == \"__main__\":\n    \n    # Define weights\n    ensemble_weights = {\n        'siglip': W_SIGLIP,\n        'dino': W_DINO\n    }\n    \n    try:\n        # Run Ensemble\n        submission = robust_ensemble(FILES, ensemble_weights)\n        \n        # Save\n        submission.to_csv(OUTPUT_FILE, index=False)\n        print(f\"\\nSuccess! Saved to {OUTPUT_FILE}\")\n        print(submission.head())\n        \n        # Sanity Check Stats\n        print(\"\\nStats:\")\n        print(submission['target'].describe())\n        \n    except Exception as e:\n        print(f\"\\nError during ensembling: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-17T23:21:03.596493Z","iopub.execute_input":"2026-01-17T23:21:03.596744Z","iopub.status.idle":"2026-01-17T23:21:03.663736Z","shell.execute_reply.started":"2026-01-17T23:21:03.596715Z","shell.execute_reply":"2026-01-17T23:21:03.663126Z"}},"outputs":[],"execution_count":null}]}