{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🌿 CSIRO Img2Bio: XGBoost + PCA + RGB Vegetation Indices\n\nThis notebook replaces DINOv2 features with **RGB-based Vegetation Indices** for biomass prediction.\nThe goal is to compare domain-specific handcrafted features against deep learning embeddings.\n\n### Approach\n- Extract **vegetation indices** (NGRDI, ExG, VARI, etc.) from RGB images\n- Add **color space features** (HSV, LAB)\n- Add **texture features** (GLCM, edge density)\n- Apply **PCA** for dimensionality reduction (optional, for fair comparison)\n- Train **XGBoost** regression models for each target\n\n### Key References\n- NGRDI achieves R² ~0.62 for grassland biomass (Hunt et al., 2005)\n- ExG, VARI validated for pasture monitoring (Bareth et al., 2018)","metadata":{}},{"cell_type":"code","source":"import os\nimport glob\nimport random\nfrom typing import Dict, List, Optional, Tuple, Any\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nfrom tqdm.auto import tqdm\n\nimport torch\nimport warnings\n\nwarnings.simplefilter(action='ignore', category=FutureWarning)\nprint(f\"PyTorch: {torch.__version__}\")\nprint(f\"Device: {torch.cuda.get_device_name(0) if torch.cuda.is_available() else 'CPU'}\")","metadata":{"trusted":true,"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 📊 Explore Data","metadata":{}},{"cell_type":"code","source":"PATH_DATA = '/kaggle/input/csiro-biomass'\nPATH_TRAIN_CSV = os.path.join(PATH_DATA, 'train.csv')\nPATH_TRAIN_IMG = os.path.join(PATH_DATA, 'train')\nPATH_TEST_IMG = os.path.join(PATH_DATA, 'test')\n\ndf = pd.read_csv(PATH_TRAIN_CSV)\nprint(f\"Dataset size: {df.shape}\")\ndisplay(df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"TARGET_COLS = [c for c in df.columns if c not in ['image_id', 'Image']]\nprint(f\"Target columns: {TARGET_COLS}\")\nprint(f\"Number of targets: {len(TARGET_COLS)}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Plot target distribution","metadata":{}},{"cell_type":"code","source":"# Exclude non-numeric or identifier columns from histogram plotting\ncols_to_plot = [col for col in TARGET_COLS if col not in ['sample_id', 'image_path', 'State', 'target_name']]\n\nfor col in cols_to_plot:\n    plt.figure(figsize=(8, 3)) # Create a new figure for each histogram\n    plt.hist(df[col].dropna(), bins=50, edgecolor='black', alpha=0.7)\n    plt.xlabel(col, fontsize=12)\n    plt.ylabel('Count', fontsize=12)\n    plt.title(f'{col} Distribution', fontsize=14, fontweight='bold')\n    plt.grid(alpha=0.3)\n    plt.xticks(rotation=45, ha=\"right\") # Rotate x-axis labels\n    plt.tight_layout() # Adjust layout to prevent overlap\n    plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cols_to_plot = ['State', 'target_name']\nn_rows, n_cols = 1, len(cols_to_plot)\nfig, axes = plt.subplots(n_rows, n_cols, figsize=(4 * n_cols, 4 * n_rows))\n\n# Ensure axes is an array even for a single subplot\naxes = axes.flatten()\n\nfor ax, col in zip(axes, cols_to_plot):\n    counts = df[col].value_counts()\n    ax.pie(counts, labels=counts.index, autopct='%1.1f%%', startangle=140)\n    ax.set_title(f'Distribution of {col}', fontsize=14, fontweight='bold')\n    ax.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\n\ngrouped_train_data = df.groupby('target_name')\n\nplt.figure(figsize=(12, 8))\nfor target_name, group_data in grouped_train_data:\n    sns.histplot(data=group_data, x='target', kde=True, label=target_name)\n\nplt.title('Distribution of Target for Each Target Name Class')\nplt.xlabel('Target')\nplt.ylabel('Frequency')\nplt.legend(title='Target Name')\nplt.grid(True) # Added grid here\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Show sample images","metadata":{}},{"cell_type":"code","source":"def show_images(df_sample, n=12, path_img=PATH_DATA):\n    \"\"\"Displays a linear sampling of images sorted by target value.\"\"\"\n\n    # Sort the DataFrame by the 'target' column\n    df_sorted = df_sample.sort_values(by='target').reset_index(drop=True)\n\n    # Perform linear sampling\n    indices_to_show = np.linspace(0, len(df_sorted) - 1, n, dtype=int)\n    df_to_show = df_sorted.iloc[indices_to_show]\n\n    # Determine the number of rows and columns for subplots\n    n_cols = 3  # You can adjust this number\n    n_rows = (n + n_cols - 1) // n_cols\n\n    fig, axes = plt.subplots(n_rows, n_cols, figsize=(16, 4 * n_rows))\n    axes = axes.flatten()\n\n    # Remove unused subplots if any\n    for i in range(n, len(axes)):\n        fig.delaxes(axes[i])\n\n    for i, (idx, row) in enumerate(df_to_show.iterrows()):\n        # Use image_path directly (includes train/ID....jpg)\n        img_path = os.path.join(path_img, row['image_path'])\n\n        if os.path.exists(img_path):\n            img = Image.open(img_path).convert('RGB')\n            axes[i].imshow(img)\n            # Include the target value in the title\n            title = f\"ID: {row['sample_id']}\\nTarget: {row['target']:.2f}\"\n            axes[i].set_title(title, fontsize=10)\n        axes[i].axis('off')\n\n    plt.tight_layout()\n    plt.show()\n\n# Example usage: Show 12 images linearly sampled based on target value\nshow_images(df, n=12)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🌱 RGB Vegetation Indices\n\nThese indices are specifically designed for vegetation monitoring and biomass estimation.\nThey encode domain knowledge about chlorophyll absorption and reflection patterns.\n\n| Index | Formula | Use Case |\n|-------|---------|----------|\n| NGRDI | (G-R)/(G+R) | Green biomass estimation |\n| ExG | 2g-r-b | Vegetation/soil separation |\n| VARI | (G-R)/(G+R-B) | Atmospheric resistant |\n| ExGR | ExG - 1.4r + g | Built-in threshold |\n| GLI | (2G-R-B)/(2G+R+B) | Green leaf index |\n| RGRI | R/G | Dead matter detection |\n| MGRVI | (G²-R²)/(G²+R²) | Modified green-red |","metadata":{}},{"cell_type":"code","source":"def compute_vegetation_indices(img: np.ndarray) -> Dict[str, float]:\n    \"\"\"Compute RGB-based vegetation indices optimized for biomass prediction.\n    \n    Args:\n        img: HxWx3 uint8 RGB array\n        \n    Returns:\n        Dictionary of vegetation index statistics\n    \"\"\"\n    img_f = img.astype(np.float32) / 255.0\n    R, G, B = img_f[:, :, 0], img_f[:, :, 1], img_f[:, :, 2]\n    eps = 1e-7\n    \n    # Normalized chromatic coordinates\n    total = R + G + B + eps\n    r, g, b = R / total, G / total, B / total\n    \n    features = {}\n    \n    # =========================================================================\n    # PRIMARY VEGETATION INDICES\n    # =========================================================================\n    \n    # 1. NGRDI - Normalized Green-Red Difference Index\n    NGRDI = (G - R) / (G + R + eps)\n    features[\"NGRDI_mean\"] = np.mean(NGRDI)\n    features[\"NGRDI_std\"] = np.std(NGRDI)\n    features[\"NGRDI_median\"] = np.median(NGRDI)\n    features[\"NGRDI_p10\"] = np.percentile(NGRDI, 10)\n    features[\"NGRDI_p25\"] = np.percentile(NGRDI, 25)\n    features[\"NGRDI_p75\"] = np.percentile(NGRDI, 75)\n    features[\"NGRDI_p90\"] = np.percentile(NGRDI, 90)\n    features[\"NGRDI_iqr\"] = features[\"NGRDI_p75\"] - features[\"NGRDI_p25\"]\n    \n    # 2. ExG - Excess Green Index\n    ExG = 2 * g - r - b\n    features[\"ExG_mean\"] = np.mean(ExG)\n    features[\"ExG_std\"] = np.std(ExG)\n    features[\"ExG_median\"] = np.median(ExG)\n    features[\"ExG_p10\"] = np.percentile(ExG, 10)\n    features[\"ExG_p90\"] = np.percentile(ExG, 90)\n    features[\"ExG_positive_frac\"] = np.mean(ExG > 0)\n    features[\"ExG_sum\"] = np.sum(np.maximum(ExG, 0))\n    \n    # 3. VARI - Visible Atmospherically Resistant Index\n    VARI = (G - R) / (G + R - B + eps)\n    VARI = np.clip(VARI, -2, 2)\n    features[\"VARI_mean\"] = np.mean(VARI)\n    features[\"VARI_std\"] = np.std(VARI)\n    features[\"VARI_median\"] = np.median(VARI)\n    \n    # 4. ExGR - Excess Green minus Excess Red\n    ExR = 1.4 * r - g\n    ExGR = ExG - ExR\n    features[\"ExGR_mean\"] = np.mean(ExGR)\n    features[\"ExGR_std\"] = np.std(ExGR)\n    features[\"ExGR_positive_frac\"] = np.mean(ExGR > 0)\n    \n    # =========================================================================\n    # SECONDARY INDICES\n    # =========================================================================\n    \n    # 5. RGRI - Red-Green Ratio Index\n    RGRI = R / (G + eps)\n    features[\"RGRI_mean\"] = np.mean(RGRI)\n    features[\"RGRI_std\"] = np.std(RGRI)\n    features[\"RGRI_median\"] = np.median(RGRI)\n    features[\"RGRI_p90\"] = np.percentile(RGRI, 90)\n    \n    # 6. GLI - Green Leaf Index\n    GLI = (2 * G - R - B) / (2 * G + R + B + eps)\n    features[\"GLI_mean\"] = np.mean(GLI)\n    features[\"GLI_std\"] = np.std(GLI)\n    features[\"GLI_median\"] = np.median(GLI)\n    \n    # 7. MGRVI - Modified Green-Red Vegetation Index\n    MGRVI = (G**2 - R**2) / (G**2 + R**2 + eps)\n    features[\"MGRVI_mean\"] = np.mean(MGRVI)\n    features[\"MGRVI_std\"] = np.std(MGRVI)\n    \n    # 8. GBVI - Green-Blue Vegetation Index\n    GBVI = (G - B) / (G + B + eps)\n    features[\"GBVI_mean\"] = np.mean(GBVI)\n    features[\"GBVI_std\"] = np.std(GBVI)\n    \n    # 9. CIVE - Color Index of Vegetation Extraction\n    CIVE = 0.441 * r - 0.811 * g + 0.385 * b + 18.787\n    features[\"CIVE_mean\"] = np.mean(CIVE)\n    features[\"CIVE_std\"] = np.std(CIVE)\n    \n    # 10. VEG - Vegetative Index\n    VEG = g / ((r ** 0.667) * (b ** 0.333) + eps)\n    VEG = np.clip(VEG, 0, 10)\n    features[\"VEG_mean\"] = np.mean(VEG)\n    features[\"VEG_std\"] = np.std(VEG)\n    \n    return features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_color_features(img: np.ndarray) -> Dict[str, float]:\n    \"\"\"Compute color space features from HSV and LAB.\n    \n    Args:\n        img: HxWx3 uint8 RGB array\n        \n    Returns:\n        Dictionary of color features\n    \"\"\"\n    features = {}\n    \n    # =========================================================================\n    # HSV Color Space\n    # =========================================================================\n    hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)\n    H, S, V = hsv[:, :, 0], hsv[:, :, 1], hsv[:, :, 2]\n    \n    features[\"HSV_H_mean\"] = np.mean(H)\n    features[\"HSV_H_std\"] = np.std(H)\n    features[\"HSV_H_median\"] = np.median(H)\n    features[\"HSV_S_mean\"] = np.mean(S)\n    features[\"HSV_S_std\"] = np.std(S)\n    features[\"HSV_S_median\"] = np.median(S)\n    features[\"HSV_V_mean\"] = np.mean(V)\n    features[\"HSV_V_std\"] = np.std(V)\n    \n    # Hue-based vegetation masks\n    green_mask = (H > 35) & (H < 85) & (S > 50)\n    features[\"green_hue_frac\"] = np.mean(green_mask)\n    \n    brown_mask = (H > 10) & (H < 35) & (S > 30)\n    features[\"brown_hue_frac\"] = np.mean(brown_mask)\n    \n    red_mask = ((H < 10) | (H > 170)) & (S > 30)\n    features[\"red_hue_frac\"] = np.mean(red_mask)\n    \n    # =========================================================================\n    # LAB Color Space\n    # =========================================================================\n    lab = cv2.cvtColor(img, cv2.COLOR_RGB2LAB)\n    L, A, B_lab = lab[:, :, 0], lab[:, :, 1], lab[:, :, 2]\n    \n    features[\"LAB_L_mean\"] = np.mean(L)\n    features[\"LAB_L_std\"] = np.std(L)\n    features[\"LAB_A_mean\"] = np.mean(A)\n    features[\"LAB_A_std\"] = np.std(A)\n    features[\"LAB_A_median\"] = np.median(A)\n    features[\"LAB_A_p10\"] = np.percentile(A, 10)\n    features[\"LAB_A_p90\"] = np.percentile(A, 90)\n    features[\"LAB_B_mean\"] = np.mean(B_lab)\n    features[\"LAB_B_std\"] = np.std(B_lab)\n    \n    # =========================================================================\n    # Simple RGB Statistics\n    # =========================================================================\n    R, G, B = img[:, :, 0], img[:, :, 1], img[:, :, 2]\n    \n    features[\"R_mean\"] = np.mean(R)\n    features[\"G_mean\"] = np.mean(G)\n    features[\"B_mean\"] = np.mean(B)\n    features[\"R_std\"] = np.std(R)\n    features[\"G_std\"] = np.std(G)\n    features[\"B_std\"] = np.std(B)\n    \n    features[\"green_dominance\"] = np.mean((G > R) & (G > B))\n    features[\"green_ratio\"] = np.mean(G) / (np.mean(R) + np.mean(B) + 1e-7)\n    features[\"brightness_mean\"] = np.mean(img)\n    features[\"brightness_std\"] = np.std(img)\n    features[\"contrast\"] = np.std(V)\n    \n    return features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_texture_features(img: np.ndarray) -> Dict[str, float]:\n    \"\"\"Compute texture features using edge detection and local statistics.\n    \n    Args:\n        img: HxWx3 uint8 RGB array\n        \n    Returns:\n        Dictionary of texture features\n    \"\"\"\n    features = {}\n    \n    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)\n    \n    # Edge features\n    edges = cv2.Canny(gray, 50, 150)\n    features[\"edge_density\"] = np.mean(edges > 0)\n    features[\"edge_mean\"] = np.mean(edges)\n    \n    # Gradient features\n    sobel_x = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)\n    sobel_y = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)\n    gradient_magnitude = np.sqrt(sobel_x**2 + sobel_y**2)\n    features[\"gradient_mean\"] = np.mean(gradient_magnitude)\n    features[\"gradient_std\"] = np.std(gradient_magnitude)\n    \n    # Laplacian\n    laplacian = cv2.Laplacian(gray, cv2.CV_64F)\n    features[\"laplacian_mean\"] = np.mean(np.abs(laplacian))\n    features[\"laplacian_std\"] = np.std(laplacian)\n    \n    # Local variance\n    kernel_size = 15\n    local_mean = cv2.blur(gray.astype(np.float32), (kernel_size, kernel_size))\n    local_sq_mean = cv2.blur(gray.astype(np.float32)**2, (kernel_size, kernel_size))\n    local_var = local_sq_mean - local_mean**2\n    features[\"local_variance_mean\"] = np.mean(local_var)\n    features[\"local_variance_std\"] = np.std(local_var)\n    \n    # GLCM-like features (simplified)\n    gray_quant = (gray // 16).astype(np.uint8)\n    h, w = gray_quant.shape\n    left = gray_quant[:, :-1].flatten()\n    right = gray_quant[:, 1:].flatten()\n    \n    contrast = np.mean((left.astype(np.float32) - right.astype(np.float32))**2)\n    features[\"glcm_contrast\"] = contrast\n    features[\"glcm_homogeneity\"] = np.mean(1.0 / (1.0 + np.abs(left - right)))\n    \n    diff_hist, _ = np.histogram(np.abs(left - right), bins=16, range=(0, 16))\n    diff_hist = diff_hist / (diff_hist.sum() + 1e-7)\n    features[\"glcm_energy\"] = np.sum(diff_hist**2)\n    features[\"glcm_entropy\"] = -np.sum(diff_hist * np.log(diff_hist + 1e-10))\n    \n    return features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_spatial_features(img: np.ndarray) -> Dict[str, float]:\n    \"\"\"Compute spatial distribution features.\n    \n    Args:\n        img: HxWx3 uint8 RGB array\n        \n    Returns:\n        Dictionary of spatial features\n    \"\"\"\n    features = {}\n    \n    img_f = img.astype(np.float32) / 255.0\n    R, G, B = img_f[:, :, 0], img_f[:, :, 1], img_f[:, :, 2]\n    \n    green_mask = ((G > R) & (G > B)).astype(np.float32)\n    h, w = green_mask.shape\n    \n    # Quadrant features\n    top_half = green_mask[:h//2, :]\n    bottom_half = green_mask[h//2:, :]\n    left_half = green_mask[:, :w//2]\n    right_half = green_mask[:, w//2:]\n    \n    features[\"green_top_frac\"] = np.mean(top_half)\n    features[\"green_bottom_frac\"] = np.mean(bottom_half)\n    features[\"green_left_frac\"] = np.mean(left_half)\n    features[\"green_right_frac\"] = np.mean(right_half)\n    \n    features[\"green_vertical_uniformity\"] = 1 - abs(features[\"green_top_frac\"] - features[\"green_bottom_frac\"])\n    features[\"green_horizontal_uniformity\"] = 1 - abs(features[\"green_left_frac\"] - features[\"green_right_frac\"])\n    \n    # Center vs edges\n    margin = h // 4\n    center = green_mask[margin:h-margin, margin:w-margin]\n    features[\"green_center_frac\"] = np.mean(center)\n    features[\"green_center_ratio\"] = features[\"green_center_frac\"] / (np.mean(green_mask) + 1e-7)\n    \n    # Row/column profiles\n    row_profile = np.mean(G, axis=1)\n    features[\"row_profile_mean\"] = np.mean(row_profile)\n    features[\"row_profile_std\"] = np.std(row_profile)\n    features[\"row_profile_gradient\"] = np.mean(np.diff(row_profile))\n    \n    col_profile = np.mean(G, axis=0)\n    features[\"col_profile_std\"] = np.std(col_profile)\n    \n    return features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def extract_all_features(img_path: str, resize: Optional[Tuple[int, int]] = None) -> Dict[str, float]:\n    \"\"\"Extract all features from a single image.\n    \n    Args:\n        img_path: Path to image file\n        resize: Optional (width, height) to resize image\n        \n    Returns:\n        Dictionary of all features\n    \"\"\"\n    img = cv2.imread(str(img_path))\n    if img is None:\n        raise ValueError(f\"Could not load image: {img_path}\")\n    \n    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n    \n    if resize is not None:\n        img = cv2.resize(img, resize, interpolation=cv2.INTER_AREA)\n    \n    features = {}\n    features.update(compute_vegetation_indices(img))\n    features.update(compute_color_features(img))\n    features.update(compute_texture_features(img))\n    features.update(compute_spatial_features(img))\n    \n    return features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Extract test image features\nimport multiprocessing\nfrom functools import partial\n\ndef _extract_features_for_row_parallel(row_dict, path_data):\n    \"\"\"Helper function to extract features for a single row,\n    designed for multiprocessing. Expects row_dict as a dictionary.\n    \"\"\"\n    image_path = os.path.join(path_data, row_dict['image_path']) # Access as dict key\n    features = extract_all_features(image_path)\n    return row_dict['sample_id'], features # Access as dict key\n\n# Create a partial function to fix path_data argument\n_partial_extract_features = partial(_extract_features_for_row_parallel, path_data=PATH_DATA)\n\n# Number of CPU cores to use. Use os.cpu_count() for better practice.\n# 'os' is already imported in cell t3uqRSVUaCn8\nnum_cores = os.cpu_count()\nprint(f\"Using {num_cores} cores for feature extraction.\")\n\nwith multiprocessing.Pool(processes=num_cores) as pool:\n    # Convert DataFrame rows to a list of dictionaries, which are easily picklable\n    # tqdm is used to show a progress bar\n    results = list(tqdm(pool.imap(_partial_extract_features, df.to_dict('records')),\n                        total=len(df), desc=\"Extracting test features\"))\n\nimage_features = {}\n# Collect results from multiprocessing\nfor sample_id, features in results:\n    image_features[sample_id] = features\n\nprint(\"Image feature extraction complete.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Extract features for all images\n# image_features = {}\n# for index, row in tqdm(df.iterrows(), total=df.shape[0]):\n#     image_path = os.path.join(PATH_DATA, row['image_path'])\n#     # Use the pipeline function\n#     features = extract_all_features(image_path)\n#     # Use the sample_id as the key for the features\n#     image_features[row['sample_id']] = features\n\n# print(\"Image feature extraction complete.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Convert the image features dictionary to a DataFrame\nimage_features_df = pd.DataFrame.from_dict(image_features, orient='index')\nimage_features_df.index.name = 'sample_id'\nimage_features_df.columns = [f'img_feature_{i}' for i in range(image_features_df.shape[1])]\n\n# One-hot encode the 'target_name' column\ndf_one_hot = pd.get_dummies(df['target_name'], prefix='target_name').astype(int) # Convert to int\n\n# Merge the image features and one-hot encoded features with the original DataFrame\ndf_combined = df.merge(image_features_df, on='sample_id', how='left')\ndf_combined = df_combined.merge(df_one_hot, left_index=True, right_index=True)\n\nprint(\"Image features and one-hot encoded target names combined with the original DataFrame.\")\ndisplay(df_combined.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Flatten all image features into a single Series\nimage_features_cols = [col for col in df_combined.columns if col.startswith(\"img_feature_\")]\nall_features_flat = image_features_df[image_features_cols].values.flatten()\nall_features_series = pd.Series(all_features_flat)\n\n# Plot a single histogram for all features\nplt.figure(figsize=(10, 6))\nsns.histplot(all_features_series, kde=True, bins=100) # Adjust bins as needed\nplt.title('Distribution of All Extracted Image Feature Values')\nplt.xlabel('Feature Value')\nplt.ylabel('Frequency')\nplt.grid(True, alpha=0.3)\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🏋️ Train Regressor","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\n\n# Separate image features and other features\nimage_feature_cols = [col for col in df_combined.columns if col.startswith('img_feature_')]\nother_features_cols = [col for col in df_combined.columns if col.startswith('target_name_')]\n\nX_image_features = df_combined[image_feature_cols]\nX_other_features = df_combined[other_features_cols]\n\n# Standardize the image features before applying PCA\nscaler = StandardScaler()\nX_image_scaled = scaler.fit_transform(X_image_features)\n\n# Apply PCA to reduce dimensionality\n# You can adjust the number of components (n_components)\npca = PCA(n_components=60) # Example: Reduce components by 20x\nX_image_pca = pca.fit_transform(X_image_scaled)\n\n# Convert the PCA reduced features to a DataFrame\nX_image_pca_df = pd.DataFrame(\n    X_image_pca, index=df_combined.index,\n    columns=[f'pca_img_feature_{i}' for i in range(pca.n_components)])\n\n# Combine the PCA reduced image features with the other features\nX_combined_pca = pd.concat([X_other_features, X_image_pca_df], axis=1)\n\nprint(f\"Original number of image features: {X_image_features.shape[1]}\")\nprint(f\"Reduced number of image features after PCA: {X_image_pca_df.shape[1]}\")\nprint(f\"Total number of features for regression after PCA: {X_combined_pca.shape[1]}\")\n\n# Display the first few rows of the combined feature DataFrame after PCA\ndisplay(X_combined_pca.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import xgboost as xgb\nfrom sklearn.model_selection import train_test_split, RandomizedSearchCV\nfrom sklearn.metrics import mean_squared_error\nfrom scipy.stats import randint, uniform, loguniform\nimport numpy as np\n\n# Define features (X) and target (y) using the PCA-reduced image features\n# Exclude identifier columns and the original target column\n# Use X_combined_pca which contains the one-hot encoded target_name and PCA image features\nfeatures = X_combined_pca.columns.tolist() # Use all columns from X_combined_pca\nprint(f\"training features: {features}\")\ntarget = 'target'\n\nX = X_combined_pca\ny = df_combined[target] # Target remains the same\n\n# Split data into training and testing sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\n# Initialize and train the XGBoost Regressor\n# The features include binary (from one-hot encoding) and float types.\n# XGBoost can handle this mix of feature types effectively.\n\n# Define a distribution for hyperparameters to sample from\n# Using distributions allows for linear or log-uniform (experimental/exponential) sampling\nparam_dist = {\n    'n_estimators': randint(100, 1000),\n    'learning_rate': loguniform(0.001, 0.2),        # Log-uniform sampling\n    'max_depth': randint(3, 11),\n    'colsample_bytree': uniform(0.5, 0.5),          # Uniform sampling [0.5, 1.0]\n    'reg_alpha': loguniform(1e-4, 10.0),            # Log-uniform sampling\n    'reg_lambda': loguniform(1.0, 10.0)             # Log-uniform sampling\n}\n\nxgb_model = xgb.XGBRegressor(objective='reg:squarederror', random_state=42)\n\n# Use RandomizedSearchCV to find the best parameters\n# n_iter controls the number of parameter settings that are sampled\nrandom_search = RandomizedSearchCV(\n    estimator=xgb_model,\n    param_distributions=param_dist,\n    n_iter=800,\n    scoring='neg_mean_squared_error',\n    cv=3,\n    verbose=2,\n    random_state=42,\n    n_jobs=-1,\n)\n\nrandom_search.fit(X_train, y_train)\n\n# Get the best parameters and the best model\nbest_params = random_search.best_params_\nxgb_model = random_search.best_estimator_\n\nprint(\"Randomized search complete.\")\nprint(f\"Best parameters found: {best_params}\")\nprint(\"XGBoost model training complete with best parameters.\")","metadata":{"trusted":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Get feature importances from the best model found by GridSearchCV\n# If you didn't use GridSearchCV, use xgb_model.feature_importances_\nfeature_importances = xgb_model.feature_importances_\n\n# Create a pandas Series for easier handling and sorting\nfeature_importances_series = pd.Series(feature_importances, index=X_train.columns)\n\n# Sort features by importance\nsorted_feature_importances = feature_importances_series.sort_values(ascending=False)\n\n# Plot the top N feature importances (e.g., top 20)\nn_top_features = 20\nplt.figure(figsize=(10, 8))\nsns.barplot(x=sorted_feature_importances.head(n_top_features).values,\n            y=sorted_feature_importances.head(n_top_features).index,\n            palette='viridis')\nplt.title(f'Top {n_top_features} Feature Importances (XGBoost)', fontsize=14)\nplt.xlabel('Importance', fontsize=12)\nplt.ylabel('Feature', fontsize=12)\nplt.tight_layout()\nplt.show()\n\n# You can also display the sorted feature importances\nprint(\"Sorted Feature Importances:\")\ndisplay(sorted_feature_importances)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Make predictions on the test data\ny_pred = xgb_model.predict(X_test)\n\nprint(\"Predictions on test data complete.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.metrics import mean_absolute_error, r2_score\n\n# Calculate Root Mean Squared Error (RMSE)\nrmse = np.sqrt(mean_squared_error(y_test, y_pred))\nprint(f\"Root Mean Squared Error (RMSE) on the test set: {rmse}\")\n\n# Calculate Mean Absolute Error (MAE)\nmae = mean_absolute_error(y_test, y_pred)\nprint(f\"Mean Absolute Error (MAE) on the test set: {mae}\")\n\n# Calculate R-squared\nr2 = r2_score(y_test, y_pred)\nprint(f\"R-squared on the test set: {r2}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🔮 Test predictions","metadata":{}},{"cell_type":"code","source":"# 1. Load test data\nPATH_TEST_CSV = os.path.join(PATH_DATA, 'test.csv')\ndf_test = pd.read_csv(PATH_TEST_CSV)\nprint(f\"Test Dataset size: {df_test.shape}\")\ndisplay(df_test.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 2. Extract test image features\ntest_image_features = {}\nfor index, row in tqdm(df_test.iterrows(), total=df_test.shape[0]):\n    image_path = os.path.join(PATH_DATA, row['image_path']) # Access as dict key\n    features = extract_all_features(image_path)\n    # Use the sample_id as the key for the features\n    test_image_features[row['sample_id']] = features\n\nprint(\"Test Image feature extraction complete.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 3. Combine test features\ntest_image_features_df = pd.DataFrame.from_dict(test_image_features, orient='index')\ntest_image_features_df.index.name = 'sample_id'\ntest_image_features_df.columns = [f'img_feature_{i}' for i in range(test_image_features_df.shape[1])]\n\n# One-hot encode the 'target_name' column for the test set\ntest_df_one_hot = pd.get_dummies(df_test['target_name'], prefix='target_name').astype(int)\n\n# Merge the test image features and one-hot encoded features with the original test DataFrame\ndf_test_combined = df_test.merge(test_image_features_df, on='sample_id', how='left')\ndf_test_combined = df_test_combined.merge(test_df_one_hot, left_index=True, right_index=True)\n\nprint(\"Test image features and one-hot encoded target names combined.\")\ndisplay(df_test_combined.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 4. Scale and PCA transform test image features\n# Separate image features and other features for the test set\ntest_image_feature_cols = [col for col in df_test_combined.columns if col.startswith('img_feature_')]\ntest_other_features_cols = [col for col in df_test_combined.columns if col.startswith('target_name_')]\n\nX_test_image_features = df_test_combined[test_image_feature_cols]\nX_test_other_features = df_test_combined[test_other_features_cols]\n\n# Apply the *trained* scaler and PCA to the test image features\nX_test_image_scaled = scaler.transform(X_test_image_features) # Use the fitted scaler\nX_test_image_pca = pca.transform(X_test_image_scaled)       # Use the fitted pca\n\n# Convert the PCA reduced features to a DataFrame\nX_test_image_pca_df = pd.DataFrame(\n    X_test_image_pca, index=df_test_combined.index,\n    columns=[f'pca_img_feature_{i}' for i in range(pca.n_components)])\n\nprint(f\"Scaled and PCA-transformed test image features shape: {X_test_image_pca_df.shape}\")\n\n# Display the first few rows of the PCA reduced test image features\ndisplay(X_test_image_pca_df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 5. Prepare test data for inference\n# Combine the PCA reduced image features with the other features (one-hot encoded target_name)\nX_test_final = pd.concat([X_test_other_features, X_test_image_pca_df], axis=1)\n\nprint(f\"Final test feature DataFrame shape: {X_test_final.shape}\")\n\n# Display the first few rows of the final test feature DataFrame\ndisplay(X_test_final.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 6. Make predictions\ntest_predictions = xgb_model.predict(X_test_final)\n\nprint(\"Predictions on test data complete.\")\ndisplay(test_predictions)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 📝 Prepare submission","metadata":{}},{"cell_type":"code","source":"# 7. Generate submission file\nsubmission_df = pd.DataFrame({'sample_id': df_test['sample_id'], 'target': test_predictions})\n\n# Ensure predictions are non-negative if the target represents a physical quantity\nsubmission_df['target'] = submission_df['target'].clip(lower=0)\n\nsubmission_path = 'submission.csv'\nsubmission_df.to_csv(submission_path, index=False)\n\nprint(f\"Submission file generated at: {submission_path}\")\ndisplay(submission_df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!head submission.csv","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}