{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "17da4353",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "✅ Import thành công\n"
     ]
    }
   ],
   "source": [
    "# Import libraries for Element84 Earth Search\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# Element84 Earth Search STAC\n",
    "import pystac_client\n",
    "from odc.stac import load\n",
    "\n",
    "# Machine Learning\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "from sklearn.metrics import classification_report, confusion_matrix, accuracy_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "print(\"✅ Import thành công\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9c063be3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "✅ Element84 Earth Search kết nối thành công\n",
      "   API: earth-search.aws.element84.com\n"
     ]
    }
   ],
   "source": [
    "# Kết nối Element84 Earth Search (hosted trên AWS)\n",
    "def connect_earth_search():\n",
    "    \"\"\"Kết nối đến Element84 Earth Search STAC API\"\"\"\n",
    "    catalog = pystac_client.Client.open(\n",
    "        \"https://earth-search.aws.element84.com/v1\"\n",
    "    )\n",
    "    return catalog\n",
    "\n",
    "catalog = connect_earth_search()\n",
    "print(\"✅ Element84 Earth Search kết nối thành công\")\n",
    "print(f\"   API: earth-search.aws.element84.com\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "83784d01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "✅ Tìm thấy 45 scenes Sentinel-2\n",
      "✅ Sentinel-2 raw: FrozenMappingWarningOnValuesAccess({'y': 8874, 'x': 9902, 'time': 28})\n",
      "   Variables: ['red', 'green', 'blue', 'nir', 'swir16', 'swir22', 'scl']\n"
     ]
    }
   ],
   "source": [
    "# Cấu hình vùng và thời gian\n",
    "date_range = (\"2022-09-01\", \"2023-10-01\")\n",
    "longtitude_range = (105.5, 106.4)\n",
    "latitude_range = (9.2, 10.0)\n",
    "\n",
    "# Tạo bounding box\n",
    "bbox = (longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1])\n",
    "\n",
    "# Query Sentinel-2 từ Element84\n",
    "search = catalog.search(\n",
    "    collections=[\"sentinel-2-l2a\"],\n",
    "    bbox=bbox,\n",
    "    datetime=f\"{date_range[0]}/{date_range[1]}\",\n",
    "    query={\"eo:cloud_cover\": {\"lt\": 30}}\n",
    ")\n",
    "\n",
    "items = search.item_collection()\n",
    "print(f\"✅ Tìm thấy {len(items)} scenes Sentinel-2\")\n",
    "\n",
    "# Load data với odc-stac\n",
    "data_sen2 = load(\n",
    "    items,\n",
    "    bands=[\"red\", \"green\", \"blue\", \"nir\", \"swir16\", \"swir22\", \"scl\"],\n",
    "    bbox=bbox,\n",
    "    resolution=10,\n",
    "    chunks={\"time\": 1, \"x\": 2048, \"y\": 2048},\n",
    "    groupby=\"solar_day\"\n",
    ")\n",
    "\n",
    "print(f\"✅ Sentinel-2 raw: {data_sen2.dims}\")\n",
    "print(f\"   Variables: {list(data_sen2.data_vars)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c3faed92",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "dimension time on 0th function argument to apply_ufunc with dask='parallelized' consists of multiple chunks, but is also a core dimension. To fix, either rechunk into a single array chunk along this dimension, i.e., ``.chunk(dict(time=-1))``, or pass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` but beware that this may significantly increase memory usage.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[4], line 32\u001b[0m\n\u001b[1;32m     30\u001b[0m data_clean \u001b[38;5;241m=\u001b[39m mask_clean(data_sen2)\n\u001b[1;32m     31\u001b[0m data_ndvi \u001b[38;5;241m=\u001b[39m calculate_indices(data_clean, index\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNDVI\u001b[39m\u001b[38;5;124m\"\u001b[39m, satellite_mission\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124ms2\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m---> 32\u001b[0m data_fill \u001b[38;5;241m=\u001b[39m \u001b[43mfill_nan\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata_ndvi\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     33\u001b[0m data_sen2_monthly \u001b[38;5;241m=\u001b[39m data_fill\u001b[38;5;241m.\u001b[39mresample(time\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m1MS\u001b[39m\u001b[38;5;124m\"\u001b[39m)\u001b[38;5;241m.\u001b[39mmean()\u001b[38;5;241m.\u001b[39mcompute()\n\u001b[1;32m     35\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m✅ S2 monthly shape: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdata_sen2_monthly\u001b[38;5;241m.\u001b[39mdims\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
      "Cell \u001b[0;32mIn[4], line 27\u001b[0m, in \u001b[0;36mfill_nan\u001b[0;34m(ds)\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfill_nan\u001b[39m(ds):\n\u001b[1;32m     26\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Fill NaN bằng interpolation theo thời gian\"\"\"\u001b[39;00m\n\u001b[0;32m---> 27\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mds\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minterpolate_na\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mtime\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mlinear\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfill_value\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mextrapolate\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/dataset.py:6773\u001b[0m, in \u001b[0;36mDataset.interpolate_na\u001b[0;34m(self, dim, method, limit, use_coordinate, max_gap, **kwargs)\u001b[0m\n\u001b[1;32m   6655\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Fill in NaNs by interpolating according to different methods.\u001b[39;00m\n\u001b[1;32m   6656\u001b[0m \n\u001b[1;32m   6657\u001b[0m \u001b[38;5;124;03mParameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   6769\u001b[0m \u001b[38;5;124;03m    D        (x) float64 40B 5.0 3.0 1.0 -1.0 4.0\u001b[39;00m\n\u001b[1;32m   6770\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   6771\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mxarray\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcore\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mmissing\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m _apply_over_vars_with_dim, interp_na\n\u001b[0;32m-> 6773\u001b[0m new \u001b[38;5;241m=\u001b[39m \u001b[43m_apply_over_vars_with_dim\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   6774\u001b[0m \u001b[43m    \u001b[49m\u001b[43minterp_na\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6775\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6776\u001b[0m \u001b[43m    \u001b[49m\u001b[43mdim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6777\u001b[0m \u001b[43m    \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6778\u001b[0m \u001b[43m    \u001b[49m\u001b[43mlimit\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlimit\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6779\u001b[0m \u001b[43m    \u001b[49m\u001b[43muse_coordinate\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43muse_coordinate\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6780\u001b[0m \u001b[43m    \u001b[49m\u001b[43mmax_gap\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmax_gap\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6781\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   6782\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   6783\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m new\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/missing.py:222\u001b[0m, in \u001b[0;36m_apply_over_vars_with_dim\u001b[0;34m(func, self, dim, **kwargs)\u001b[0m\n\u001b[1;32m    220\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m name, var \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata_vars\u001b[38;5;241m.\u001b[39mitems():\n\u001b[1;32m    221\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m dim \u001b[38;5;129;01min\u001b[39;00m var\u001b[38;5;241m.\u001b[39mdims:\n\u001b[0;32m--> 222\u001b[0m         ds[name] \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvar\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    223\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    224\u001b[0m         ds[name] \u001b[38;5;241m=\u001b[39m var\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/missing.py:367\u001b[0m, in \u001b[0;36minterp_na\u001b[0;34m(self, dim, use_coordinate, method, limit, max_gap, keep_attrs, **kwargs)\u001b[0m\n\u001b[1;32m    365\u001b[0m     warnings\u001b[38;5;241m.\u001b[39mfilterwarnings(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mignore\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moverflow\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;167;01mRuntimeWarning\u001b[39;00m)\n\u001b[1;32m    366\u001b[0m     warnings\u001b[38;5;241m.\u001b[39mfilterwarnings(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mignore\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minvalid value\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;167;01mRuntimeWarning\u001b[39;00m)\n\u001b[0;32m--> 367\u001b[0m     arr \u001b[38;5;241m=\u001b[39m \u001b[43mapply_ufunc\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    368\u001b[0m \u001b[43m        \u001b[49m\u001b[43minterpolator\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    369\u001b[0m \u001b[43m        \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m    370\u001b[0m \u001b[43m        \u001b[49m\u001b[43mindex\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    371\u001b[0m \u001b[43m        \u001b[49m\u001b[43minput_core_dims\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[43m[\u001b[49m\u001b[43mdim\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43mdim\u001b[49m\u001b[43m]\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    372\u001b[0m \u001b[43m        \u001b[49m\u001b[43moutput_core_dims\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[43m[\u001b[49m\u001b[43mdim\u001b[49m\u001b[43m]\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    373\u001b[0m \u001b[43m        \u001b[49m\u001b[43moutput_dtypes\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    374\u001b[0m \u001b[43m        \u001b[49m\u001b[43mdask\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mparallelized\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m    375\u001b[0m \u001b[43m        \u001b[49m\u001b[43mvectorize\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    376\u001b[0m \u001b[43m        \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    377\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mtranspose(\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdims)\n\u001b[1;32m    379\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m limit \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    380\u001b[0m     arr \u001b[38;5;241m=\u001b[39m arr\u001b[38;5;241m.\u001b[39mwhere(valids)\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/computation.py:1265\u001b[0m, in \u001b[0;36mapply_ufunc\u001b[0;34m(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, on_missing_core_dim, *args)\u001b[0m\n\u001b[1;32m   1263\u001b[0m \u001b[38;5;66;03m# feed DataArray apply_variable_ufunc through apply_dataarray_vfunc\u001b[39;00m\n\u001b[1;32m   1264\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28many\u001b[39m(\u001b[38;5;28misinstance\u001b[39m(a, DataArray) \u001b[38;5;28;01mfor\u001b[39;00m a \u001b[38;5;129;01min\u001b[39;00m args):\n\u001b[0;32m-> 1265\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mapply_dataarray_vfunc\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   1266\u001b[0m \u001b[43m        \u001b[49m\u001b[43mvariables_vfunc\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1267\u001b[0m \u001b[43m        \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1268\u001b[0m \u001b[43m        \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1269\u001b[0m \u001b[43m        \u001b[49m\u001b[43mjoin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mjoin\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1270\u001b[0m \u001b[43m        \u001b[49m\u001b[43mexclude_dims\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mexclude_dims\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1271\u001b[0m \u001b[43m        \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1272\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1273\u001b[0m \u001b[38;5;66;03m# feed Variables directly through apply_variable_ufunc\u001b[39;00m\n\u001b[1;32m   1274\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28many\u001b[39m(\u001b[38;5;28misinstance\u001b[39m(a, Variable) \u001b[38;5;28;01mfor\u001b[39;00m a \u001b[38;5;129;01min\u001b[39;00m args):\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/computation.py:307\u001b[0m, in \u001b[0;36mapply_dataarray_vfunc\u001b[0;34m(func, signature, join, exclude_dims, keep_attrs, *args)\u001b[0m\n\u001b[1;32m    302\u001b[0m result_coords, result_indexes \u001b[38;5;241m=\u001b[39m build_output_coords_and_indexes(\n\u001b[1;32m    303\u001b[0m     args, signature, exclude_dims, combine_attrs\u001b[38;5;241m=\u001b[39mkeep_attrs\n\u001b[1;32m    304\u001b[0m )\n\u001b[1;32m    306\u001b[0m data_vars \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mgetattr\u001b[39m(a, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mvariable\u001b[39m\u001b[38;5;124m\"\u001b[39m, a) \u001b[38;5;28;01mfor\u001b[39;00m a \u001b[38;5;129;01min\u001b[39;00m args]\n\u001b[0;32m--> 307\u001b[0m result_var \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mdata_vars\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    309\u001b[0m out: \u001b[38;5;28mtuple\u001b[39m[DataArray, \u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m] \u001b[38;5;241m|\u001b[39m DataArray\n\u001b[1;32m    310\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m signature\u001b[38;5;241m.\u001b[39mnum_outputs \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n",
      "File \u001b[0;32m/env/lib/python3.12/site-packages/xarray/core/computation.py:764\u001b[0m, in \u001b[0;36mapply_variable_ufunc\u001b[0;34m(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)\u001b[0m\n\u001b[1;32m    762\u001b[0m             \u001b[38;5;28;01mfor\u001b[39;00m axis, dim \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(core_dims, start\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;28mlen\u001b[39m(core_dims)):\n\u001b[1;32m    763\u001b[0m                 \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(data\u001b[38;5;241m.\u001b[39mchunks[axis]) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[0;32m--> 764\u001b[0m                     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    765\u001b[0m                         \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdimension \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdim\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m on \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mn\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124mth function argument to \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    766\u001b[0m                         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mapply_ufunc with dask=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mparallelized\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m consists of \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    767\u001b[0m                         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmultiple chunks, but is also a core dimension. To \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    768\u001b[0m                         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfix, either rechunk into a single array chunk along \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    769\u001b[0m                         \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mthis dimension, i.e., ``.chunk(dict(\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdim\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m=-1))``, or \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    770\u001b[0m                         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    771\u001b[0m                         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbut beware that this may significantly increase memory usage.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    772\u001b[0m                     )\n\u001b[1;32m    773\u001b[0m     dask_gufunc_kwargs[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mallow_rechunk\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m    775\u001b[0m output_sizes \u001b[38;5;241m=\u001b[39m dask_gufunc_kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moutput_sizes\u001b[39m\u001b[38;5;124m\"\u001b[39m, {})\n",
      "\u001b[0;31mValueError\u001b[0m: dimension time on 0th function argument to apply_ufunc with dask='parallelized' consists of multiple chunks, but is also a core dimension. To fix, either rechunk into a single array chunk along this dimension, i.e., ``.chunk(dict(time=-1))``, or pass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` but beware that this may significantly increase memory usage."
     ]
    }
   ],
   "source": [
    "# Tiền xử lý Sentinel-2: cloud mask + NDVI + resampling\n",
    "def mask_clean(ds):\n",
    "    \"\"\"Cloud masking sử dụng SCL band (Scene Classification Layer)\"\"\"\n",
    "    if \"scl\" not in ds:\n",
    "        print(\"⚠️  Không có SCL band, bỏ qua cloud masking\")\n",
    "        return ds\n",
    "    \n",
    "    cloud_mask = (\n",
    "        (ds[\"scl\"] == 3) |   # cloud shadow\n",
    "        (ds[\"scl\"] == 8) |   # cloud medium probability\n",
    "        (ds[\"scl\"] == 9) |   # cloud high probability\n",
    "        (ds[\"scl\"] == 10)    # thin cirrus\n",
    "    )\n",
    "    \n",
    "    ds_masked = ds.where(~cloud_mask)\n",
    "    return ds_masked.drop_vars(\"scl\", errors=\"ignore\")\n",
    "\n",
    "def calculate_indices(ds, index=\"NDVI\", satellite_mission=\"s2\"):\n",
    "    \"\"\"Tính chỉ số NDVI cho Sentinel-2\"\"\"\n",
    "    if index == \"NDVI\":\n",
    "        ndvi = (ds[\"nir\"] - ds[\"red\"]) / (ds[\"nir\"] + ds[\"red\"] + 1e-8)\n",
    "        ds[\"NDVI\"] = ndvi\n",
    "    return ds\n",
    "\n",
    "def fill_nan(ds):\n",
    "    \"\"\"Fill NaN bằng interpolation theo thời gian\"\"\"\n",
    "    return ds.interpolate_na(dim=\"time\", method=\"linear\", fill_value=\"extrapolate\")\n",
    "\n",
    "# Áp dụng tiền xử lý\n",
    "data_clean = mask_clean(data_sen2)\n",
    "data_ndvi = calculate_indices(data_clean, index=\"NDVI\", satellite_mission=\"s2\")\n",
    "data_fill = fill_nan(data_ndvi)\n",
    "data_sen2_monthly = data_fill.resample(time=\"1MS\").mean().compute()\n",
    "\n",
    "print(f\"✅ S2 monthly shape: {data_sen2_monthly.dims}\")\n",
    "print(f\"   Variables: {list(data_sen2_monthly.data_vars)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "569bfebb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tải Sentinel-1 (SAR VV/VH) từ Element84\n",
    "search_s1 = catalog.search(\n",
    "    collections=[\"sentinel-1-grd\"],\n",
    "    bbox=bbox,\n",
    "    datetime=f\"{date_range[0]}/{date_range[1]}\",\n",
    "    query={\n",
    "        \"sat:orbit_state\": {\"eq\": \"descending\"},\n",
    "        \"sar:product_type\": {\"eq\": \"GRD\"}\n",
    "    }\n",
    ")\n",
    "\n",
    "items_s1 = search_s1.item_collection()\n",
    "print(f\"✅ Tìm thấy {len(items_s1)} scenes Sentinel-1\")\n",
    "\n",
    "if len(items_s1) > 0:\n",
    "    data_sen1 = load(\n",
    "        items_s1,\n",
    "        bands=[\"vv\", \"vh\"],\n",
    "        bbox=bbox,\n",
    "        resolution=10,\n",
    "        chunks={\"time\": 1, \"x\": 2048, \"y\": 2048},\n",
    "        groupby=\"solar_day\"\n",
    "    )\n",
    "    data_sen1_monthly = data_sen1.resample(time=\"1MS\").mean().compute()\n",
    "else:\n",
    "    # Tạo dummy data nếu không có S1\n",
    "    print(\"⚠️  Không có Sentinel-1, tạo dummy data\")\n",
    "    data_sen1_monthly = xr.Dataset({\n",
    "        \"vv\": xr.DataArray(\n",
    "            np.zeros_like(data_sen2_monthly[\"red\"].values),\n",
    "            coords=data_sen2_monthly[\"red\"].coords,\n",
    "            dims=data_sen2_monthly[\"red\"].dims\n",
    "        ),\n",
    "        \"vh\": xr.DataArray(\n",
    "            np.zeros_like(data_sen2_monthly[\"red\"].values),\n",
    "            coords=data_sen2_monthly[\"red\"].coords,\n",
    "            dims=data_sen2_monthly[\"red\"].dims\n",
    "        )\n",
    "    })\n",
    "\n",
    "print(f\"✅ S1 monthly shape: {data_sen1_monthly.dims}\")\n",
    "print(f\"   Variables: {list(data_sen1_monthly.data_vars)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc56c2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Ánh xạ nhãn lớp đất\n",
    "label_mapping = {\n",
    "    \"Lua tom\": \"0\", \"Lua\": \"1\", \"CHN\": \"2\", \"CLN\": \"3\",\n",
    "    \"TS\": \"4\", \"Song\": \"5\", \"Dat xay dung\": \"6\", \"Rung\": \"7\",\n",
    "}\n",
    "\n",
    "def load_train_data(csv_path=\"/media/x79/2A7D-FAA0/remote-sensing/train_data.csv\", label_mapping=None):\n",
    "    \"\"\"Load training points từ CSV\"\"\"\n",
    "    df = pd.read_csv(csv_path)\n",
    "    if label_mapping:\n",
    "        df[\"label\"] = df[\"class_name\"].map(label_mapping).astype(str)\n",
    "    return df\n",
    "\n",
    "def get_data_sen1_and_sen2(train_df, sen2_data, sen1_data):\n",
    "    \"\"\"Extract features từ S2 và S1 tại các điểm training\"\"\"\n",
    "    X_list = []\n",
    "    y_list = []\n",
    "    \n",
    "    for idx, row in train_df.iterrows():\n",
    "        try:\n",
    "            lon, lat = float(row[\"longitude\"]), float(row[\"latitude\"])\n",
    "            label = int(row[\"label\"])\n",
    "            \n",
    "            # Extract S2 features\n",
    "            s2_point = sen2_data.sel(x=lon, y=lat, method=\"nearest\")\n",
    "            s2_features = []\n",
    "            \n",
    "            for var in [\"red\", \"green\", \"blue\", \"nir\", \"swir16\", \"swir22\", \"NDVI\"]:\n",
    "                if var in s2_point:\n",
    "                    vals = s2_point[var].values\n",
    "                    if vals.size > 0:\n",
    "                        s2_features.extend([\n",
    "                            np.nanmean(vals), np.nanstd(vals),\n",
    "                            np.nanmin(vals), np.nanmax(vals)\n",
    "                        ])\n",
    "                    else:\n",
    "                        s2_features.extend([0, 0, 0, 0])\n",
    "            \n",
    "            # Extract S1 features\n",
    "            s1_point = sen1_data.sel(x=lon, y=lat, method=\"nearest\")\n",
    "            s1_features = []\n",
    "            \n",
    "            for var in [\"vv\", \"vh\"]:\n",
    "                if var in s1_point:\n",
    "                    vals = s1_point[var].values\n",
    "                    if vals.size > 0:\n",
    "                        s1_features.extend([np.nanmean(vals), np.nanstd(vals)])\n",
    "                    else:\n",
    "                        s1_features.extend([0, 0])\n",
    "            \n",
    "            features = s2_features + s1_features\n",
    "            \n",
    "            if not np.isnan(features).any() and not np.isinf(features).any():\n",
    "                X_list.append(features)\n",
    "                y_list.append(label)\n",
    "        except Exception as e:\n",
    "            continue\n",
    "    \n",
    "    return X_list, y_list\n",
    "\n",
    "def split_train_data(X, y, test_size=0.2, val_size=0.1, random_state=42):\n",
    "    \"\"\"Split data thành train/val/test\"\"\"\n",
    "    X = np.array(X, dtype=np.float32)\n",
    "    y = np.array(y, dtype=np.int64)\n",
    "    \n",
    "    X_train, X_temp, y_train, y_temp = train_test_split(\n",
    "        X, y, test_size=test_size + val_size, random_state=random_state, stratify=y\n",
    "    )\n",
    "    \n",
    "    val_ratio = val_size / (test_size + val_size)\n",
    "    X_val, X_test, y_val, y_test = train_test_split(\n",
    "        X_temp, y_temp, test_size=(1 - val_ratio), random_state=random_state, stratify=y_temp\n",
    "    )\n",
    "    \n",
    "    return X_train, X_val, X_test, y_train, y_val, y_test\n",
    "\n",
    "train_data = load_train_data(label_mapping=label_mapping)\n",
    "X, y = get_data_sen1_and_sen2(train_data, data_sen2_monthly, data_sen1_monthly)\n",
    "X_train, X_val, X_test, y_train, y_val, y_test = split_train_data(X, y, test_size=0.2, val_size=0.1)\n",
    "\n",
    "X_train_np = np.array(X_train, dtype=np.float32)\n",
    "X_val_np   = np.array(X_val,   dtype=np.float32)\n",
    "X_test_np  = np.array(X_test,  dtype=np.float32)\n",
    "y_train_np = np.array(y_train, dtype=np.int64)\n",
    "y_val_np   = np.array(y_val,   dtype=np.int64)\n",
    "y_test_np  = np.array(y_test,  dtype=np.int64)\n",
    "\n",
    "n_features = X_train_np.shape[1]\n",
    "n_classes  = len(np.unique(y_train_np))\n",
    "\n",
    "print(f\"✅ Train: {X_train_np.shape}  Val: {X_val_np.shape}  Test: {X_test_np.shape}\")\n",
    "print(f\"   n_features={n_features}  n_classes={n_classes}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5492528b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "\n",
    "# ── MobileNetV3 + LR-ASPP classifier ──────────────────────────────────────────\n",
    "class MobileNetLRASPPClassifier(nn.Module):\n",
    "    \"\"\"\n",
    "    MobileNetV3-inspired backbone with LR-ASPP (Lite Reduced ASPP) head\n",
    "    for land-use classification on flat feature vectors.\n",
    "    \"\"\"\n",
    "    def __init__(self, n_features, n_classes):\n",
    "        super().__init__()\n",
    "        # Feature extraction backbone\n",
    "        self.feature_extractor = nn.Sequential(\n",
    "            nn.Linear(n_features, 128), nn.BatchNorm1d(128), nn.ReLU(inplace=True), nn.Dropout(0.2),\n",
    "            nn.Linear(128, 256),        nn.BatchNorm1d(256), nn.ReLU(inplace=True), nn.Dropout(0.3),\n",
    "            nn.Linear(256, 512),        nn.BatchNorm1d(512), nn.ReLU(inplace=True), nn.Dropout(0.3),\n",
    "        )\n",
    "        # LR-ASPP Branch 1: global pooling → 128\n",
    "        self.global_pool = nn.AdaptiveAvgPool1d(1)\n",
    "        self.global_conv = nn.Sequential(nn.Linear(512, 128), nn.ReLU(inplace=True))\n",
    "        # LR-ASPP Branch 2: direct 1×1 → 128\n",
    "        self.branch_conv = nn.Sequential(nn.Linear(512, 128), nn.BatchNorm1d(128), nn.ReLU(inplace=True))\n",
    "        # Fusion → n_classes\n",
    "        self.classifier = nn.Sequential(\n",
    "            nn.Linear(256, 128), nn.BatchNorm1d(128), nn.ReLU(inplace=True), nn.Dropout(0.4),\n",
    "            nn.Linear(128, n_classes),\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        feat        = self.feature_extractor(x)\n",
    "        global_feat = self.global_pool(feat.unsqueeze(-1)).squeeze(-1)\n",
    "        global_feat = self.global_conv(global_feat)\n",
    "        branch_feat = self.branch_conv(feat)\n",
    "        fused       = torch.cat([global_feat, branch_feat], dim=1)\n",
    "        return self.classifier(fused)\n",
    "\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "print(f\"✅ Device: {device}\")\n",
    "model = MobileNetLRASPPClassifier(n_features, n_classes).to(device)\n",
    "print(model)\n",
    "total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)\n",
    "print(f\"   Trainable params: {total_params:,}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da40f5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# ── Hyper-parameters ───────────────────────────────────────────────────────────\n",
    "LEARNING_RATE     = 1e-3\n",
    "BATCH_SIZE        = 64\n",
    "N_EPOCHS          = 60\n",
    "PATIENCE          = 10\n",
    "CONV_THRESHOLD    = 1e-4   # cải thiện val_loss < THRESHOLD → đánh dấu hội tụ\n",
    "CONV_WINDOW       = 3      # cần CONV_WINDOW bước liên tiếp thỏa mãn\n",
    "\n",
    "# ── Tensors & DataLoaders ──────────────────────────────────────────────────────\n",
    "X_tr_t = torch.FloatTensor(X_train_np)\n",
    "y_tr_t = torch.LongTensor(y_train_np)\n",
    "X_va_t = torch.FloatTensor(X_val_np)\n",
    "y_va_t = torch.LongTensor(y_val_np)\n",
    "\n",
    "train_loader = DataLoader(TensorDataset(X_tr_t, y_tr_t), batch_size=BATCH_SIZE, shuffle=True)\n",
    "val_loader   = DataLoader(TensorDataset(X_va_t, y_va_t), batch_size=BATCH_SIZE, shuffle=False)\n",
    "\n",
    "# ── Class-weighted loss ────────────────────────────────────────────────────────\n",
    "class_counts  = np.bincount(y_train_np)\n",
    "class_weights = 1.0 / (class_counts + 1e-6)\n",
    "class_weights = class_weights / class_weights.sum() * n_classes\n",
    "criterion = nn.CrossEntropyLoss(weight=torch.FloatTensor(class_weights).to(device))\n",
    "\n",
    "optimizer = optim.Adam(model.parameters(), lr=LEARNING_RATE, weight_decay=1e-4)\n",
    "scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode=\"min\", factor=0.5, patience=5)\n",
    "\n",
    "# ── Training loop ──────────────────────────────────────────────────────────────\n",
    "train_losses, val_losses, val_accs, lr_history = [], [], [], []\n",
    "best_val_loss    = float(\"inf\")\n",
    "best_epoch       = 1\n",
    "best_state_dict  = None\n",
    "patience_counter = 0\n",
    "convergence_epoch = None   # ← điểm hội tụ sẽ được ghi lại ở đây\n",
    "early_stop_epoch  = None\n",
    "\n",
    "print(\"🚀 Training MobileNetV3 + LR-ASPP...\")\n",
    "print(f\"   conv_threshold={CONV_THRESHOLD}  conv_window={CONV_WINDOW}  patience={PATIENCE}\")\n",
    "for epoch in range(1, N_EPOCHS + 1):\n",
    "    # --- Train ---\n",
    "    model.train()\n",
    "    epoch_loss = 0.0\n",
    "    for bx, by in train_loader:\n",
    "        bx, by = bx.to(device), by.to(device)\n",
    "        optimizer.zero_grad()\n",
    "        loss = criterion(model(bx), by)\n",
    "        loss.backward()\n",
    "        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)\n",
    "        optimizer.step()\n",
    "        epoch_loss += loss.item()\n",
    "    avg_train_loss = epoch_loss / len(train_loader)\n",
    "\n",
    "    # --- Validate ---\n",
    "    model.eval()\n",
    "    val_loss = 0.0; correct = 0; total = 0\n",
    "    with torch.no_grad():\n",
    "        for bx, by in val_loader:\n",
    "            bx, by = bx.to(device), by.to(device)\n",
    "            out = model(bx)\n",
    "            val_loss += criterion(out, by).item()\n",
    "            pred = out.argmax(1)\n",
    "            correct += (pred == by).sum().item()\n",
    "            total   += by.size(0)\n",
    "    avg_val_loss = val_loss / len(val_loader)\n",
    "    val_acc = correct / total\n",
    "\n",
    "    scheduler.step(avg_val_loss)\n",
    "    lr = optimizer.param_groups[0][\"lr\"]\n",
    "    lr_history.append(lr)\n",
    "\n",
    "    train_losses.append(avg_train_loss)\n",
    "    val_losses.append(avg_val_loss)\n",
    "    val_accs.append(val_acc)\n",
    "\n",
    "    # ── Phát hiện điểm hội tụ (sliding window trên val_loss) ──────────────\n",
    "    if convergence_epoch is None and epoch >= CONV_WINDOW + 1:\n",
    "        window = val_losses[-(CONV_WINDOW + 1):]\n",
    "        improvements = [abs(window[i] - window[i - 1]) for i in range(1, len(window))]\n",
    "        if all(imp < CONV_THRESHOLD for imp in improvements):\n",
    "            convergence_epoch = epoch - CONV_WINDOW + 1\n",
    "            print(f\"  📍 Hội tụ phát hiện tại epoch {convergence_epoch}  \"\n",
    "                  f\"(val_loss={val_losses[convergence_epoch-1]:.4f})\")\n",
    "\n",
    "    if epoch % 5 == 0 or epoch == 1:\n",
    "        print(f\"  Epoch {epoch:3d}/{N_EPOCHS}  train={avg_train_loss:.4f}  \"\n",
    "              f\"val={avg_val_loss:.4f}  acc={val_acc:.4f}  lr={lr:.2e}\")\n",
    "\n",
    "    # ── Early stopping ─────────────────────────────────────────────────────\n",
    "    if avg_val_loss < best_val_loss:\n",
    "        best_val_loss   = avg_val_loss\n",
    "        best_epoch      = epoch\n",
    "        best_state_dict = {k: v.clone() for k, v in model.state_dict().items()}\n",
    "        patience_counter = 0\n",
    "    else:\n",
    "        patience_counter += 1\n",
    "        if patience_counter >= PATIENCE:\n",
    "            early_stop_epoch = epoch\n",
    "            print(f\"⏹  Early stopping tại epoch {epoch}  (patience={PATIENCE})\")\n",
    "            break\n",
    "\n",
    "# Restore best weights\n",
    "model.load_state_dict(best_state_dict)\n",
    "total_epochs = len(train_losses)\n",
    "print(f\"\\n✅ Training hoàn tất!  Best epoch={best_epoch}  Best val_loss={best_val_loss:.4f}\")\n",
    "if convergence_epoch:\n",
    "    print(f\"   Điểm hội tụ  : epoch {convergence_epoch}\")\n",
    "\n",
    "# ═══════════════════════════════════════════════════════════════════════════════\n",
    "# PHÂN TÍCH ĐIỂM HỘI TỤ — MobileNetV3 + LR-ASPP\n",
    "# EMA-smoothed curves + annotated convergence / best / early-stop markers\n",
    "# ═══════════════════════════════════════════════════════════════════════════════\n",
    "epochs_axis = list(range(1, total_epochs + 1))\n",
    "\n",
    "# Exponential Moving Average smoothing\n",
    "def ema(values, alpha=0.2):\n",
    "    s = [values[0]]\n",
    "    for v in values[1:]:\n",
    "        s.append(alpha * v + (1 - alpha) * s[-1])\n",
    "    return s\n",
    "\n",
    "val_losses_ema  = ema(val_losses,  alpha=0.3)\n",
    "train_losses_ema = ema(train_losses, alpha=0.3)\n",
    "val_accs_ema    = ema(val_accs,    alpha=0.3)\n",
    "\n",
    "# ── Tính ΔVal-loss per epoch ──────────────────────────────────────────────────\n",
    "delta_val = [abs(val_losses_ema[i] - val_losses_ema[i-1])\n",
    "             for i in range(1, len(val_losses_ema))]\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(16, 10))\n",
    "\n",
    "# --- Top-left: Loss curves ---\n",
    "ax = axes[0, 0]\n",
    "ax.plot(epochs_axis, train_losses, \"b-\", alpha=0.25, linewidth=0.8)\n",
    "ax.plot(epochs_axis, train_losses_ema, \"b-\", linewidth=2,   label=\"Train loss (EMA)\")\n",
    "ax.plot(epochs_axis, val_losses,   \"g-\", alpha=0.25, linewidth=0.8)\n",
    "ax.plot(epochs_axis, val_losses_ema,   \"g-\",  linewidth=2,   label=\"Val loss (EMA)\")\n",
    "ax.axvline(x=best_epoch,       color=\"red\",    linestyle=\"--\", linewidth=1.5,\n",
    "           label=f\"Best  epoch={best_epoch}\")\n",
    "if convergence_epoch:\n",
    "    ax.axvline(x=convergence_epoch, color=\"orange\", linestyle=\":\", linewidth=1.5,\n",
    "               label=f\"HỘI TỤ epoch={convergence_epoch}\")\n",
    "if early_stop_epoch:\n",
    "    ax.axvline(x=early_stop_epoch,  color=\"gray\",   linestyle=\"-.\", linewidth=1.5,\n",
    "               label=f\"Early stop epoch={early_stop_epoch}\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "ax.set_ylabel(\"Loss\")\n",
    "ax.set_title(\"Loss Curves (raw + EMA)\")\n",
    "ax.legend(fontsize=8)\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "# --- Top-right: Val accuracy ---\n",
    "ax = axes[0, 1]\n",
    "ax.plot(epochs_axis, val_accs,     \"g-\", alpha=0.3, linewidth=0.8)\n",
    "ax.plot(epochs_axis, val_accs_ema, \"g-\", linewidth=2, label=\"Val acc (EMA)\")\n",
    "ax.axvline(x=best_epoch, color=\"red\", linestyle=\"--\", linewidth=1.5,\n",
    "           label=f\"Best  epoch={best_epoch}  ({val_accs[best_epoch-1]*100:.2f}%)\")\n",
    "if convergence_epoch:\n",
    "    ax.axvline(x=convergence_epoch, color=\"orange\", linestyle=\":\", linewidth=1.5,\n",
    "               label=f\"HỘI TỤ epoch={convergence_epoch}  ({val_accs[convergence_epoch-1]*100:.2f}%)\")\n",
    "if early_stop_epoch:\n",
    "    ax.axvline(x=early_stop_epoch, color=\"gray\", linestyle=\"-.\", linewidth=1.5,\n",
    "               label=f\"Early stop\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "ax.set_ylabel(\"Accuracy\")\n",
    "ax.set_title(\"Val Accuracy Curve\")\n",
    "ax.legend(fontsize=8)\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "# --- Bottom-left: ΔVal-loss (marginal improvement) ---\n",
    "ax = axes[1, 0]\n",
    "ax.bar(epochs_axis[1:], delta_val,\n",
    "       color=[\"green\" if d > CONV_THRESHOLD else \"salmon\" for d in delta_val],\n",
    "       alpha=0.75)\n",
    "ax.axhline(y=CONV_THRESHOLD, color=\"red\", linestyle=\"--\",\n",
    "           label=f\"Threshold = {CONV_THRESHOLD:.0e}\")\n",
    "if convergence_epoch:\n",
    "    ax.axvline(x=convergence_epoch, color=\"orange\", linestyle=\":\", linewidth=1.5,\n",
    "               label=f\"HỘI TỤ epoch={convergence_epoch}\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "ax.set_ylabel(\"|ΔVal Loss|\")\n",
    "ax.set_title(\"Marginal Val-Loss Improvement per Epoch\")\n",
    "ax.legend(fontsize=8)\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "# --- Bottom-right: Learning rate schedule ---\n",
    "ax = axes[1, 1]\n",
    "ax.semilogy(epochs_axis, lr_history, \"purple\", linewidth=2)\n",
    "if convergence_epoch:\n",
    "    ax.axvline(x=convergence_epoch, color=\"orange\", linestyle=\":\", linewidth=1.5,\n",
    "               label=f\"HỘI TỤ epoch={convergence_epoch}\")\n",
    "ax.axvline(x=best_epoch, color=\"red\", linestyle=\"--\", linewidth=1.5,\n",
    "           label=f\"Best epoch={best_epoch}\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "ax.set_ylabel(\"Learning Rate (log)\")\n",
    "ax.set_title(\"Learning Rate Schedule (ReduceLROnPlateau)\")\n",
    "ax.legend(fontsize=8)\n",
    "ax.grid(True, alpha=0.3)\n",
    "\n",
    "plt.suptitle(\"MobileNetV3 + LR-ASPP — Convergence Analysis\", fontsize=13, fontweight=\"bold\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# ── Tổng kết ──────────────────────────────────────────────────────────────────\n",
    "print(f\"\\n{'═'*60}\")\n",
    "print(f\"  Tổng số epoch       : {total_epochs}\")\n",
    "print(f\"  Best epoch          : {best_epoch}  (val_loss={best_val_loss:.4f})\")\n",
    "print(f\"  Best val accuracy   : {val_accs[best_epoch-1]*100:.4f}%\")\n",
    "if convergence_epoch:\n",
    "    print(f\"  Điểm HỘI TỤ        : epoch {convergence_epoch}  \"\n",
    "          f\"(val_acc={val_accs[convergence_epoch-1]*100:.2f}%)\")\n",
    "    wasted = total_epochs - convergence_epoch\n",
    "    print(f\"  Epochs sau hội tụ   : {wasted}  \"\n",
    "          f\"(có thể giảm N_EPOCHS không ảnh hưởng nhiều đến kết quả)\")\n",
    "if early_stop_epoch:\n",
    "    print(f\"  Early stop tại epoch: {early_stop_epoch}\")\n",
    "print(f\"{'═'*60}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc690b92",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import classification_report, confusion_matrix, accuracy_score\n",
    "import seaborn as sns\n",
    "\n",
    "# ── Đánh giá trên tập test ─────────────────────────────────────────────────────\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    X_te_t = torch.FloatTensor(X_test_np).to(device)\n",
    "    logits  = model(X_te_t)\n",
    "    y_pred  = logits.argmax(1).cpu().numpy()\n",
    "\n",
    "acc = accuracy_score(y_test_np, y_pred)\n",
    "print(f\"Test Accuracy : {acc:.4f} ({acc*100:.2f}%)\\n\")\n",
    "print(classification_report(y_test_np, y_pred, digits=4))\n",
    "\n",
    "# ── Confusion matrix ────────────────────────────────────────────────────────────\n",
    "class_names = list(label_mapping.keys())\n",
    "cm = confusion_matrix(y_test_np, y_pred)\n",
    "plt.figure(figsize=(9, 7))\n",
    "sns.heatmap(cm, annot=True, fmt=\"d\", cmap=\"Purples\",\n",
    "            xticklabels=class_names, yticklabels=class_names)\n",
    "plt.xlabel(\"Predicted\")\n",
    "plt.ylabel(\"Actual\")\n",
    "plt.title(\"Confusion Matrix — MobileNetV3 + LR-ASPP\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}