{"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":"gpu","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"},{"sourceId":4533,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":3325,"modelId":986},{"sourceId":4537,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":3329,"modelId":986}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Topic: https://www.kaggle.com/competitions/csiro-biomass/discussion/654028","metadata":{}},{"cell_type":"code","source":"import os\nfrom tqdm import tqdm\n\nimport pandas as pd\nimport numpy as np\nimport numpy.linalg as la\nimport torch\nimport matplotlib.pyplot as plt\nfrom PIL import Image\n\nfrom transformers import AutoImageProcessor, AutoModel\n\nfrom sklearn.metrics.pairwise import cosine_similarity\n%matplotlib inline","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:42.658018Z","iopub.execute_input":"2025-12-07T07:47:42.658275Z","iopub.status.idle":"2025-12-07T07:47:43.507692Z","shell.execute_reply.started":"2025-12-07T07:47:42.658248Z","shell.execute_reply":"2025-12-07T07:47:43.506233Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ROOT_DIR = '/kaggle/input/csiro-biomass/'\ndevice = 'cuda'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.508054Z","iopub.status.idle":"2025-12-07T07:47:43.508279Z","shell.execute_reply.started":"2025-12-07T07:47:43.508172Z","shell.execute_reply":"2025-12-07T07:47:43.508182Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_train = pd.read_csv('/kaggle/input/csiro-biomass/train.csv')\ndf_train = df_train.pivot(\n    index=[\"image_path\", \"Sampling_Date\", \"State\", \"Species\", \"Pre_GSHH_NDVI\", \"Height_Ave_cm\"],\n    columns=\"target_name\",\n    values=\"target\",\n).reset_index()\ndf_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.509409Z","iopub.status.idle":"2025-12-07T07:47:43.509645Z","shell.execute_reply.started":"2025-12-07T07:47:43.509529Z","shell.execute_reply":"2025-12-07T07:47:43.509542Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Define Backbone Embedding Model\nEach image will be cropped into m_crop by n_crop square regions","metadata":{}},{"cell_type":"code","source":"processor = AutoImageProcessor.from_pretrained('/kaggle/input/dinov2/pytorch/small/1', use_fast=True)\nmodel = AutoModel.from_pretrained('/kaggle/input/dinov2/pytorch/small/1').to(device)\nm_crop = 4\nn_crop = 2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.510707Z","iopub.status.idle":"2025-12-07T07:47:43.511012Z","shell.execute_reply.started":"2025-12-07T07:47:43.510845Z","shell.execute_reply":"2025-12-07T07:47:43.510857Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Helper Functions","metadata":{}},{"cell_type":"code","source":"def embed_img(images, model, processor):\n    with torch.no_grad():\n        inputs = processor(images=images, return_tensors=\"pt\", device=device)\n        outputs = model(**inputs)\n        last_hidden_states = outputs.last_hidden_state\n        output = last_hidden_states[:, 0, :]\n    return output\n\ndef crop_to_mxn(img: Image.Image, m: int = 4, n: int = 2):\n    \"\"\"\n    Split a PIL image into an m x n grid of *square* crops.\n\n    - m: number of tiles horizontally (x direction)\n    - n: number of tiles vertically (y direction)\n    - Tiles are all squares of the same side length.\n    - Grid is aligned with the top-left corner of the image.\n    - Side length is chosen so that the full m x n grid fits\n      inside the image (no overflow), maximizing the used area.\n    - Any leftover strip on the right/bottom is ignored.\n\n    Returns a list of length m*n in row-major order:\n    [ (row=0, col=0), (0,1), ..., (0,m-1),\n      (1,0), (1,1), ..., ...\n      (n-1,0), ..., (n-1,m-1) ]\n    \"\"\"\n    w, h = img.size\n\n    # Max square side that allows an m x n grid to fit\n    side = min(w // m, h // n)\n    if side <= 0:\n        raise ValueError(\n            f\"Image ({w}x{h}) is too small to fit a {m}x{n} grid of square tiles.\"\n        )\n\n    crops: List[Image.Image] = []\n    for row in range(n):        # vertical (y)\n        for col in range(m):    # horizontal (x)\n            left   = col * side\n            upper  = row * side\n            right  = left + side\n            lower  = upper + side\n            crops.append(img.crop((left, upper, right, lower)))\n\n    return crops\n\n\ndef make_grid_embeddings(images, model, processor, m: int = 4, n: int = 2):\n    \"\"\"\n    images: list of length n_s, each a PIL Image\n    m: number of splits per row\n    n: number of splits per col\n\n    returns:\n        patch_embeds: np.ndarray of shape (n_s, m*n, d)\n                      where dim 1 indexes the patches in row-major order.\n    \"\"\"\n    n_s = len(images)\n    num_patches_per_img = m * n\n\n    # 1) build list of all crops\n    all_crops = []\n    for img in images:\n        crops = crop_to_mxn(img, m=m, n=n)   # list of m*n PIL images\n        all_crops.extend(crops)           # total length = n_s * m * n\n\n    # 2) embed all crops at once\n    crop_feats = embed_img(all_crops, model, processor)\n\n    # 3) reshape back to (n_s, m*n, d), convert to numpy if desired\n    crop_feats = crop_feats.detach().cpu().numpy()\n    d = crop_feats.shape[1]\n    patch_embeds = crop_feats.reshape(n_s, num_patches_per_img, d)\n\n    return patch_embeds  # shape: (n_s,m*n, d)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.512034Z","iopub.status.idle":"2025-12-07T07:47:43.512350Z","shell.execute_reply.started":"2025-12-07T07:47:43.512191Z","shell.execute_reply":"2025-12-07T07:47:43.512206Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Embed Image Crops","metadata":{}},{"cell_type":"code","source":"embedded_patches = []\npaths = df_train['image_path'].to_list()\nb = 16\nfor i in tqdm(range(0, len(paths), b)):\n    batch_paths = paths[i:i+b]\n    imgs = [Image.open(os.path.join(ROOT_DIR, p)) for p in batch_paths]\n    batch_embedded_patches = make_grid_embeddings(imgs, model, processor, m_crop, n_crop)\n    embedded_patches.extend(list(batch_embedded_patches))\n\nembedded_patches = np.stack(embedded_patches, axis=0)\nembedded_patches.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.513815Z","iopub.status.idle":"2025-12-07T07:47:43.514117Z","shell.execute_reply.started":"2025-12-07T07:47:43.513987Z","shell.execute_reply":"2025-12-07T07:47:43.514002Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Task Setup See Description\n\n### Graph Laplacian\n\nLet \\(W \\in \\mathbb{R}^{N \\times N}\\) be the patch similarity matrix and define:\n\n$$\nD_{pp} = \\sum_q W_{pq}, \\quad\nL = D - W.\n$$\n\nThen the smoothness term can be written as:\n\n$$\n\\frac{1}{2} \\sum_{p,q} W_{pq} (x_p - x_q)^2\n= x^\\top L x.\n$$\n\nOptionally define a strictly positive-definite matrix:\n\n$$\nQ = L + \\varepsilon I, \\quad \\varepsilon > 0.\n$$\n\n### Constraint matrix\n\nDefine \\(A \\in \\mathbb{R}^{n_s \\times N}\\) such that row \\(i\\) selects the patches of image \\(i\\):\n\n$$\nA_{i,p} =\n\\begin{cases}\n1 & \\text{if patch } p \\text{ belongs to image } i, \\\\\n0 & \\text{otherwise}.\n\\end{cases}\n$$\n\nThen the constraints\n$$\n\\sum_{j=1}^{n_p} x_{i,j} = y_i\n$$\nbecome\n$$\nA x = y.\n$$\n\n### Final matrix form\n\n$$\n\\min_{x \\in \\mathbb{R}^N} \\quad \\frac{1}{2} x^\\top Q x\n\\quad \\text{s.t.} \\quad A x = y.\n$$","metadata":{}},{"cell_type":"code","source":"n_s, n_p, d = embedded_patches.shape   # n_p = m*n\nN = n_s * n_p\n\n# A: (n_s, N)\nA = np.zeros((n_s, N), dtype=float)\nfor i in range(n_s):\n    start = i * n_p\n    end   = (i + 1) * n_p\n    A[i, start:end] = 1.0\n\npatch_feats = embedded_patches.reshape(N, d)  # (N, d)\n\n\n\n# 1) cosine similarity matrix\nC = cosine_similarity(patch_feats)  # (N, N), values in [-1, 1]\n\n# 2) make non-negative weights; choice 1: clip negatives to 0\nW = np.abs(C)\n\nnp.fill_diagonal(W, 0.0)  # no self-loops (optional)\n\n# 3) Laplacian\nD = np.diag(W.sum(axis=1))\nL = D - W\n\neps = 1e-3\nQ = L + eps * np.eye(N)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.514853Z","iopub.status.idle":"2025-12-07T07:47:43.515142Z","shell.execute_reply.started":"2025-12-07T07:47:43.514995Z","shell.execute_reply":"2025-12-07T07:47:43.515009Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Solve the Optimization Problem in Closed Form","metadata":{}},{"cell_type":"code","source":"all_x = {}\nfor i, label in enumerate(['Dry_Clover_g', 'Dry_Green_g', 'Dry_Dead_g']):\n    y = df_train[label].to_numpy()\n    \n    Q_inv = la.inv(Q)\n    Q_inv_A_T = Q_inv @ A.T\n    AQA_inv = la.inv(A@Q_inv_A_T)\n    x = Q_inv_A_T @ (AQA_inv @ y)\n    \n    x = x.reshape(n_s, -1)\n    all_x[label] = x","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.516556Z","iopub.status.idle":"2025-12-07T07:47:43.516890Z","shell.execute_reply.started":"2025-12-07T07:47:43.516722Z","shell.execute_reply":"2025-12-07T07:47:43.516737Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Visualize","metadata":{}},{"cell_type":"code","source":"def show_image_and_crops(idx, df, label, all_x, root_dir, m: int = 4, n: int = 2, figsize=(10, 6)):\n    \"\"\"\n    idx      : row index into df (and x)\n    df       : dataframe with column 'image_path'\n    x        : np array of per-patch values, shape (n_samples, m*n)\n               we'll use x[idx] -> length m*n\n    root_dir : base directory for images\n    n        : split image into m x n crops\n    figsize  : matplotlib figure size\n    \"\"\"\n    # --- load image ---\n    img_path = os.path.join(root_dir, df.loc[idx, 'image_path'])\n    img = Image.open(img_path)\n\n    # --- crops & patch values for this image ---\n    crops = crop_to_mxn(img, m=m, n=n)        # length m*n\n    x = all_x[label]\n    x_idx = x[idx]                         # shape (m*n,)\n\n    # --- figure with a big panel on the left + grid on the right ---\n    fig = plt.figure(figsize=figsize)\n\n    # grid: n rows, n+1 columns; left column is the big image\n    gs = fig.add_gridspec(nrows=n, ncols=m+1,\n                          width_ratios=[m] + [1]*m)\n\n    # left: original image spanning all rows of first column\n    ax_big = fig.add_subplot(gs[:, 0])\n    ax_big.imshow(img)\n    ax_big.set_title(f\"{df.loc[idx, label]:.2f}\")\n    ax_big.axis('off')\n\n    # right: m x n small crops\n    axes = []\n    for r in range(n):\n        for c in range(1, m+1):\n            axes.append(fig.add_subplot(gs[r, c]))\n\n    for i, ax in enumerate(axes):\n        if i < len(crops):\n            ax.imshow(crops[i])\n            # assume scalar per patch; format to 2 decimals\n            ax.set_title(f\"{x_idx[i]:.2f}\")\n            ax.axis('off')\n        else:\n            ax.axis('off')\n\n    plt.tight_layout()\n    plt.show()\n    return fig\n\n\nfig = show_image_and_crops(32, df_train, 'Dry_Clover_g', all_x, ROOT_DIR, m=m_crop, n=n_crop)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.517888Z","iopub.status.idle":"2025-12-07T07:47:43.518179Z","shell.execute_reply.started":"2025-12-07T07:47:43.518027Z","shell.execute_reply":"2025-12-07T07:47:43.518040Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Collect the labels into a Dataset","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\noutput = {\n    'image_path':[],\n    'crop_path': [],\n    'target_name': [],\n    'target': []\n}\n\nout_dir = 'crops'\nos.makedirs(out_dir, exist_ok=True)\nfor i, row in tqdm(df_train.iterrows(), total=len(df_train)):\n    img_path = row['image_path']\n    img = Image.open(os.path.join(ROOT_DIR, img_path))\n    crops = crop_to_mxn(img, m_crop, n_crop)\n    for j, crop in enumerate(crops):\n        crop_name = f\"{Path(row['image_path']).stem}_{j}.jpg\"\n        crop_path = os.path.join(out_dir, crop_name)\n        crop.save(crop_path)\n\n        for target_name in ['Dry_Clover_g', 'Dry_Green_g', 'Dry_Dead_g']:\n            x = all_x[target_name][i,j]\n            output['image_path'].append(img_path)\n            output['crop_path'].append(crop_path)\n            output['target_name'].append(target_name)\n            output['target'].append(x)\n\n        x_Dry_Total_g = sum([x_target[i,j] for x_target in all_x.values()])\n        output['image_path'].append(img_path)\n        output['crop_path'].append(crop_path)\n        output['target_name'].append('Dry_Total_g')\n        output['target'].append(x_Dry_Total_g)\n\n        x_GDM_g = all_x['Dry_Clover_g'][i,j] + all_x['Dry_Green_g'][i,j]\n        output['image_path'].append(img_path)\n        output['crop_path'].append(crop_path)\n        output['target_name'].append('GDM_g')\n        output['target'].append(x_GDM_g)\n\noutput = pd.DataFrame(output)\noutput.to_csv('crop.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.518920Z","iopub.status.idle":"2025-12-07T07:47:43.519195Z","shell.execute_reply.started":"2025-12-07T07:47:43.519052Z","shell.execute_reply":"2025-12-07T07:47:43.519065Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"p = df_train.loc[10, 'image_path']\nx = [x_t[10] for x_t in all_x.values()]\nprint(x)\noutput.loc[output['image_path']==p]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-07T07:47:43.520142Z","iopub.status.idle":"2025-12-07T07:47:43.520365Z","shell.execute_reply.started":"2025-12-07T07:47:43.520255Z","shell.execute_reply":"2025-12-07T07:47:43.520266Z"}},"outputs":[],"execution_count":null}]}