diff --git a/notebooks/Gallery.ipynb b/notebooks/Gallery.ipynb index 47a615ce..644e323c 100644 --- a/notebooks/Gallery.ipynb +++ b/notebooks/Gallery.ipynb @@ -35,12 +35,14 @@ "\n", "These tutorials can be run on a 4GB GPU using relatively low volumes of data (3-10GB). They will also work in HPC environments.\n", "\n", - "| Title | Description | Image | Notebooks | Last Tested |\n", + "| Topic | Description | Image | Notebooks | Last Tested |\n", "|-------|--------------|-------|-------------|-------------|\n", "| **Simplified weather model** | Train a reduced-size weather model on a standard GPU with fetchable dataset | ![Image showing FourCastMini prediction outputs](https://pyearthtools.readthedocs.io/en/latest/_images/notebooks_tutorial_FourCastMini_Demo_18_1.png) | [Train and run a simplified global weather model (low hardware and data requirements)](./tutorial/FourCastMini_Demo.ipynb) | 18 Aug 2025 |\n", "| **MLX Demo** | Shows how to integrate PyEarthTools with a non-PyTorch framework (Apple MLX) optimised for M-series chips | ![Image showing weather model outputs from MLX demo](https://pyearthtools.readthedocs.io/en/latest/_images/notebooks_tutorial_MLX-Demo-Custom-Arch_13_1.png) | [MLX Framework Example](./tutorial/MLX-Demo-Custom-Arch.ipynb) | 8 Jun 2025 | \n", "| **Convolutional Neural Net on ERA5** | Shows all steps to train a CNN on ERA5, running on CPU or a standard GPU | ![Image showing weather model outputs](https://pyearthtools.readthedocs.io/en/latest/_images/notebooks_tutorial_CNN-Model-Training_44_1.png) | [End-to-end CNN Training Example](./tutorial/CNN-Model-Training.ipynb) | 25 Aug 2025 |\n", - "| **Radar Visualisation** | Shows how to visualise radar data as a time-series, in 2D and in 3D | ![Image showing a top down view of radar data](https://pyearthtools.readthedocs.io/en/latest/_images/notebooks_RadarVisualisation_10_1.png) | [Radar Visualisation](./RadarVisualisation.ipynb) | 23 Aug 2025 |\n" + "| **Radar Visualisation** | Shows how to visualise radar data as a time-series, in 2D and in 3D | ![Image showing a top down view of radar data](https://pyearthtools.readthedocs.io/en/latest/_images/notebooks_RadarVisualisation_10_1.png) | [Radar Visualisation](./RadarVisualisation.ipynb) | 23 Aug 2025 |\n", + "| **LUCIE Climate Model** | Train a climate model | (no image) | [LUCIE-Training](./tutorial/LUCIE/LUCIE-Training.ipynb) | 13 Nov 2025 |\n", + "| **LUCIE Climate Model** | Make predictions from a climate model | (no image) | [LUCIE-Inference](./tutorial/LUCIE/LUCIE-Inference.ipynb) | 13 Nov 2025 |\n" ] }, { diff --git a/notebooks/tutorial/LUCIE/LUCIE-Inference.ipynb b/notebooks/tutorial/LUCIE/LUCIE-Inference.ipynb new file mode 100644 index 00000000..b8be428a --- /dev/null +++ b/notebooks/tutorial/LUCIE/LUCIE-Inference.ipynb @@ -0,0 +1,147 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "175e2165-f568-48e2-bedb-1245603b1ab5", + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "import lucie\n", + "import lucie.inference\n", + "from pathlib import Path\n", + "import numpy as np\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "bc8460bd-8691-403f-8cbb-dbeb4e39875a", + "metadata": {}, + "outputs": [], + "source": [ + "device = torch.device(\"mps\" if torch.backends.mps.is_available() else \"cpu\")\n", + "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else device)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e5000929-4fec-4b6c-82c1-a5769287aa38", + "metadata": {}, + "outputs": [], + "source": [ + "regridded_path = Path.home() / 'dev/data/lucie' / 'era5_T30_regridded.npz'\n", + "regridded_data = lucie.train.load_data(regridded_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "3facc4a4-ec3f-4fc8-be29-c244ec4268e2", + "metadata": {}, + "outputs": [], + "source": [ + "preprocessed_path = Path.home() / 'dev/data/lucie' / 'era5_T30_preprocessed.npz'\n", + "preprocessed_data = np.load(preprocessed_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "fb5d5b70-681c-4cc4-b973-7b259bb6e81d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 1min 49s, sys: 30.3 s, total: 2min 19s\n", + "Wall time: 1min 54s\n" + ] + } + ], + "source": [ + "%%time\n", + "%%capture\n", + "\n", + "# Note - these timings were obtained on a laptop, not on a high-performance GPU.\n", + "\n", + "predictions = lucie.inference.load_data_and_predict(device, regridded_data, preprocessed_data,model_weights_pth='model.pth')" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "9f9e517e-a6e6-4cff-bc82-fe0cff6c89ec", + "metadata": {}, + "outputs": [], + "source": [ + "da = xr.DataArray(predictions)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e0f847a4-edd0-4dc4-9a2c-278c5df6127e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGxCAYAAADCo9TSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcsxJREFUeJztnQl8HcWd56u736nbkizJxjY2R7CJIRBCwJAhLDA4hCVcuwmBcIWFhDWEY0OABBI2hDE5SciCmckSw04gzDATAnEYwhXMEMxlMMEcDhiDDbYs27Juvau791MtJOv/63Y/yc9+T8/v9/Xnfax6fVVXVbf+qqrfrwzXdV1FCCGEEFIkzGJdiBBCCCGEwQchhBBCig57PgghhBBSVBh8EEIIIaSoMPgghBBCSFFh8EEIIYSQosLggxBCCCFFhcEHIYQQQopKRO3mOI6j1q9fr2pra5VhGKXODiGEkAmM9t3s7e1VU6dOVaa56/4+T6VSKpPJFHyeWCymEomEKjd2++BDBx7Tp08vdTYIIYSUEevWrVPTpk3bZYHHrD1rVHuHXfC52tra1Jo1a8ouANntgw/d46H5zCe/qSJW3Ps50xAT+0RSsgEYOcd3nky9PCZbLSPiyKA8xrVkL4sdl+lcQqYN/yWVE4F9wAk/MijT2SqZJysjt1tpmU7Xy/3dgNaA1zCzcrvhyO2DTZCHNJwQ/pAws353fywbzIOCDiwnCl/A7k5UpqP9cgcbj/fqC/IJ7wjHgiy54dfI1MhrWAF/8OTg3eFCthLdbmg5uVC2DtRnHI7H9uWdA46JDGD9u+HPQQrbJNQltEGNkef9G+2zQ+siB8+WmZPXSGyWjXBgStKfB2jHuSQUJhSVYYevSoF1Y8Fz4wGnMOF5zVXD8w/3ncXt+JjkcH9/FiKDeeofnj18T1mwPZIKfw8GtYl4jxP6bCFYN7EeWTDRfnnjZsbfwAwbbmTUbeTstHr6tVtGfnfsCjKZjBd4rFm+p6qr3fHelZ5eR8065H3vfAw+JhjDQy068IhEht7uThSCjxwEH8ofCeAxbhReurnwh86AX3BuLH/wYeQLPuAl68TgFz+83SL4go3lDz4seMma8C4xYbuFeXDzBB/4xgwoG8yD7xdBnuDDgODDwoAnIPhQsrp9v/AMKCvfSxmuYeE9+a+o3Hh48GHF3NBywuDDd99wPLYv7xxwXxa0MWxT2OYiUFcO3nfAUlJmLvy+I1EIPvBdDfWHbSoC9xmJJvIGHy7cF7Y5DHDytuGAJoZtxoSy8b0jMADG7djuMRiF9uXly85T/1CfmGd8B0Xs/MGHC/UViY4v+FD47oX2EYlA8OEEBB/4jg9ol8UYpq+rNQsKPsqZ3b7ngxBCCJmI2K6j8nSi5T2+XGHwQQghhJQAR7nep5Djy5XK7O8hhBBCSMlgzwchhBBSAhzvX2HHlysMPgghhJASYLuu9ynk+HKlYoIPM2cr0x2a9Vz17la5EWY1O9X+aeHJfqmNjCeklCDTKI8ZmGSGzkTP1sJMdZjx7+UD1AqGI4/JgnwTZbAoQU01GqHSyzRI4Lx9usJnu6MKoGqTEyq9xePdgBnlsb5wKaxfrhkuQe7bA6R5vfLoRFeQugmk0dAkEiBBRSk1DsXWrIcZ+QP+a+ZAKo15QOl01UZ5jmyNlAlYaSdcat3glxWgZByPQWltvNsJlX/WrpcNO9Lr15zayUhoOdiooOqVZZkAaSVK5R1QR8R6/A/bQGs0tA1hWabrrdBnC6XaboCCw6dOslHGhQeEpwOEY4JEZ/5fVJF+mbaj4bJmbB8Dk+WNNrwNWl5d/l2ospPnTHTKY7KNUhpdvaZHnhB+ARtbuuX26ipfHlx4fztV2+Rtrq/gya6gYoIPQgghZCLhVPCEUwYfhBBCSAlwlKvsCg0+qHYhhBBCSFFhzwchhBBSAhwOuxBCCCGkmNhUuxBCCCGkmDgffQo5vlypmGGXbH1cuZEh/Zs9GVa1hMWPUE6mcWFFtRzIUnFxLDsWLr3DVoOy2qBz+q4RDZfuoVww2h++WFrtOv8CTCj3xDzg6p64QFfjm9lQWV2qyd8Eo31O+KJeIHvM1kZC89TwjryvWK+8iQwcHyStjAwaoXLNXJWUGBqw4BZuD1ppDCXELixOF9sq5d5O1AotW8zDYHMkv/zTkl/mktiGZLlkQDKObTBTI7/I1ASsboarFEPZ4Eq5qUaZx0gKpLiwsirOycvWWnnbHLb7vim4kqDvFPKScAksF02uOnxBPbta5snIwbOYBouAKCwCCftn6v2Zxuc5CipVBBenq+6Q16x9b0DuELSQ4ABKaaUUNtsg389mRpZDaqpcbdaCVcmtAKsEXx4G5XvJTm57mdqw0CjZNXDCKSGEEFIC7I/ULoV8xsOiRYvUgQceqOrq6rzPvHnz1H/8x3+MbE+lUmrBggWqqalJ1dTUqNNPP11t3LhRnGPt2rXqxBNPVFVVVaqlpUVdddVVKpcLMKrKA4MPQgghpATYbuGf8TBt2jR18803q+XLl6uXXnpJHXPMMerkk09Wr7/+urf9iiuuUH/4wx/U/fffr5YuXarWr1+vTjvttG35tW0v8MhkMurZZ59Vd999t7rrrrvUd7/73XHfe8UMuxBCCCGVzEknnSTSN910k9cb8txzz3mByZ133qnuvfdeLyjRLF68WM2ZM8fbfvjhh6tHH31UvfHGG+rxxx9Xra2t6qCDDlI33nijuvrqq9UNN9ygYrFtTrH5YM8HIYQQUsIJp04Bnx1F92Lcd999qr+/3xt+0b0h2WxWHXfccSP7zJ49W82YMUMtW7bMS+v/DzjgAC/wGGb+/Pmqp6dnpPdkrLDngxBCCCkBjjKUXcBaMvp4jf7lP5p4PO59gnjttde8YEPP79DzOh544AG1//77qxUrVng9Fw0NDWJ/HWi0t7d7P+v/Rwcew9uHt40H9nwQQgghZcz06dNVfX39yGfhwoXb3Xe//fbzAo3nn39eXXzxxercc8/1hlKKDXs+CCGEkBLguEOfQo7XrFu3zlOvDLO9Xg+N7t3YZ599vJ8POeQQ9eKLL6pf/OIX6ktf+pI3kbSrq0v0fmi1S1tbm/ez/v+FF14Q5xtWwwzvM1YqJvjYdEBcWR9VCGruUevuBMyZwWMwHeTTMRorNf4+KJ83CB6SCb+GBT4A/qXn5YihCb4SmhgsX+6Al4SLngz9dqj3hAPLoyc3+pdYN2FJ9FwSvCeqrdDl7KP94JcBxZBuQIMU/31XbZCFmauBJdf7wCcAlm3HPGIeUgHL2eOS6Fj2gy2yQcS34lLy6LkC5wMvCxP212RhOXv0qxhstkI9OVx4o6SgbqrbA7xkoI1EwacjB3mKDDihz2K0V9ZNFuqu5v3+QB+gMP8ZrBu8JpZ1vCt8f42dkGknAc9nnzwo3Srvy+q3xuf7gY3Qe7Zk2kyHvwvrV4f7+phvvS+3RwJ+xSRkWZu1siDMfnmfuVqZie5Zsj7j3fCsWfL8TkAWsN3mRlmN2LoMnlFFwS5w2GX42GHp7I7gOI5Kp9NeIBKNRtUTTzzhSWw1q1at8qS1ephGo//Xk1Q7Ojo8ma3mscce866th27GQ8UEH4QQQkglc+2116oTTjjBm0Ta29vrKVueeuop9ac//ckbrrngggvUlVdeqRobG72A4tJLL/UCDq100Rx//PFekHH22WerH/3oR948j+uuu87zBgnrbQmCwQchhBBSAuyd1PMxVnSPxTnnnKM2bNjgBRvacEwHHn//93/vbb/llluUaZpez4fuDdFKlttvv33keMuy1JIlS7y5Ijooqa6u9uaMfP/73x933hl8EEIIISXAcQ3vU8jx40H7eISRSCTUbbfd5n22x5577qkefvhhVSgMPgghhJAK6PmYSFBqSwghhJCiwp4PQgghpATYyvQ+O358+cLggxBCCCkBboFzPvTx5UrFBB/pJqXMYTm5CVp1W1agDXp7TbQ3vJJ93iHgwWGDnj5X5ebx5Ai4hhF+zcigTFtpeY3EVhkn5xJmqMeDxolY4R4a4AMR7QXTFPRQgXKwGyJ5NfgW+I+gz0NikzQncKIyUwb4eFiDMlPpSX6JmGuZoXmyEzLfTtwMrSsrI/McAT8FL595FmpAXw68LycKnhumTGdrsL4D2pwb7kVhQbvGxSVivfIENliqZKv9bQy9ImI94BUDZY+eGw7cZ2pyPPSadnKUqcNHpOvM0GcFyyU1Se5vwO42+GOM9pEYxoI2YMeM8DYEvh5GFnw94FHyeZNE/O81K4XtWoXe12Cz3D+5STaA7EF7i3R0i99TxY3LRmFk5TkGp8jCGpgMPj/4DoFnD/OcrfVlwXefo+vLyefJRHYKFRN8EEIIIRMJu4InnDL4IIQQQkqA7ZreZ8ePV2UL1S6EEEIIKSrs+SCEEEJKgKMM5RTQB+DgZKQygsEHIYQQUgLsCp7zwWEXQgghhBSViun5sJOOcpNDki7DkdGimQlfinxoH5mO9apQeRf2hqHszSer86+wrmxQgFogATOz4fKxTB0usW2FSi0DpZ4QWKOcM9ovD0o3ShmdNeiEXiNorhXKUlFqmW8589hmKe9zIyjFAwly1m/VYwzKwo1shMLOynQM5KJubbVI55qkfNDqR82qUna11Gc6cVlf8T5YUr1rQJ7AlUtq2wmQSffJ+0w1+x//nulmaJvK1mHDlkkzA89WNvw50lRvkOmtH4NygGxGB8IlyGYO2j08W7lG/8OGcu7uWeGvRpTB5qrDZe9YjoFy3BontCxRWotlj7LZ2vfk9myNPw+ZhvB8O9HwdM9MeROJTngYA/4wR6n7wJRkqCQ83RD+nsI82bX5yx7bVK5m27PhROwymnDqqnKlYoIPQgghZOLN+TAKOr5cYfBBCCGElACnQHv1cp5wyjkfhBBCCCkq7PkghBBCSoDNOR+EEEIIKfawi8NhF0IIIYSQXQ+HXQghhJASYLuG9ynk+HKlYoIPwzaUAfr/YXLVaD7h3y9TD7tE0CtEhXtmoPzdCdfbB54DPRVAjh4ZCPcWQb28zzfEzu81gt4EFngRpOMyHQE/Ews8VayU31wEfT1ySdD523DNRpmpqq1QEDEU9UufAbM/YH37nj6Zpz7pHWLEoCCS0kzAqZUFF+mR18jVxfP6H0Q3dMsdTJgfbssKMzPg49EkDRB6p0t/i1SzLwt5PVTcKPiZ+HxgYLsFHhxp/xz3XLU8SaJDbrey+ZanhzYH1Z9qUqHeI16+oZ0asE+8R6azkIdoX3ge0VdCY0GzM7rMUP8KN+KOy0uofxrsHiAvwHzlq3/0M4oMuqHP6mAzPCf6ea1PhL5bffcNeUo1QxuDNofvB6xLjZU2tt8uM8XTYdgFql1sql0IIYQQQsZGxfR8EEIIIRMJxzW9z44fX74+Hww+CCGEkBJgc9iFEEIIIaQ4sOeDEEIIKQFOgYqVoLVAywUGH4QQQkhZmoyZqlyZUDm/+eablWEY6vLLLx/5LpVKqQULFqimpiZVU1OjTj/9dLVx48aS5pMQQgjZWfbqdgGfcmXC9Hy8+OKL6h//8R/VgQceKL6/4oor1B//+Ed1//33q/r6enXJJZeo0047Tf3lL38Z1/mdSVmlkh/5HAxKvwMTvCqC6tMGLxC7Vm43QTeOYM8a7h/tN/Jq8B2ZbZWuz+PT4YT7CqD/Afp+aJJb5Gzq1CSZTyeKHhzy+Eyt3B7vlueLmv77TtfL72I9oOOHCd5mVn6R3qNOpO2ErNBIv8yklZL+Gt416pPyGgOysH1zzA2ZZzshK2+wTZ4vW+1vZJEUeKC0oVmE3D44ORLabrHs0Usm0+w3dsF26cTgTuPQqMA7xwAfD/RcCFoBHK+RbobnEY6Jgf1JBNptpi683Qe1cweaAHqF5PU/ge0+3x+/rYuyE/k8UvCc4XXje2+hLQiczztHUt5ItkqmrQF5Ujsh87C1xggtFzPAWwm9grBN5GqcUK8YF9qgkZV5dGtlZRo9Ab/mjIDfD8M/DwYYg5CdzoQIm/r6+tRZZ52lfvWrX6lJkyaNfN/d3a3uvPNO9bOf/Uwdc8wx6pBDDlGLFy9Wzz77rHruuedKmmdCCCGkEBxlFPwpVyZE8KGHVU488UR13HHHie+XL1+ustms+H727NlqxowZatmyZSXIKSGEELJzsDnsUjruu+8+9fLLL3vDLkh7e7uKxWKqoUH2F7e2tnrbgkin095nmJ4e8EUmhBBCSOX2fKxbt05ddtll6p577lGJBAyA7iALFy705oYMf6ZPn75TzksIIYTsCpMxu4BPuVLSnOthlY6ODvXJT35SRSIR77N06VJ16623ej/rHo5MJqO6urrEcVrt0tbWFnjOa6+91psrMvzRAQ4hhBAy0XBco+BPuVJStcuxxx6rXnvtNfHd+eef783ruPrqq71ei2g0qp544glPYqtZtWqVWrt2rZo3b17gOePxuPchhBBCyMSkpMFHbW2tmjt3rviuurra8/QY/v6CCy5QV155pWpsbFR1dXXq0ksv9QKPww8/vES5JoQQQgrHKXDopJxNxiaMz8f2uOWWW5Rpml7Ph55IOn/+fHX77beP+zxmxFZmdMjXwAEvAvTTUKC390C9OphNoA+IkQO/g6w8PjIQ7o/h5Rm9BCCfJsjRc1XhOv9sXb7t/vu2q1AQH+6hgMqvxBa4BpwvqNfQiRihOn/cbqXAq6BWGhrEO2VB2Qm53QgyQIDvjHRABYlMyDx17i99PTLgC2OBD4QmMoh5gO1p8MNoMEKvkauW6WwjmFkElL0NHguG6YYfg14TUWggA1ZoXXrfJcBDB5419K/Igg9MpF/uYKVwe2hVDV0TOkuzNSrca6R3fL4fQfWNfiO5avD+6Q1/3u0kvlPC84TtI+gc6NOiwKclO8kOfQ/iey+QuDyH0StfbE4VXMMxwtNwPtyOXiZDF0Vjlu38POFXtTVVuTLhgo+nnnpKpPVE1Ntuu837EEIIIaT8mXDBByGEEFIJ2MrwPoUcX64w+CCEEEJKgMNhF0IIIYQUE7vA3os8M9EmNOU7W4UQQgghZQmHXQghhJAS4HDYZfenvmFQWR8tFx1vlpLDqCU7rzI5f0zWn5brQPdulbpWAySFKgkdYlmQk0VVaDpIxhqTRq8qU59HepvMs8w3qB6j/UZeCSJKBjGPKBlGGXO6If+S24mtUhqXgWW7LZCcOjEzXHpbJzNhDcL5Ybsm1iPbSG5SIo9MORKax9wUeQ92gPov3RguxcTlzNPNIFkMKEuRZeznxAag/xqpkvftotIWLmH3yoZr1shG6GRQJxvQ2YrSyRgs6x6T9+lE5bOWq5HbbVhC3bXM0OcmqB2j+hrbcQ6etSjIXE1UpMJ2Tc0GeZ8mKKGRNEiMs7lwqS0SJOdHyb8dxzYmG4DVD9LpJllwkTr5grAD6tuxwYagGQp/IBIu/4V260IHvlUFbTBA/mvVyWtakW11YZvoH7DrF5bbUQo5ttSUb84JIYQQMq61zw499FDP4LOlpUWdcsopnmv4aFavXq1OPfVUNXnyZM/Y84tf/KK3pMloOjs71VlnneVt1wu/ajPQvr6+sWeEwQchhBBSGlxlKKeAjz5+POi10xYsWKCee+459dhjj6lsNquOP/541d8/5MSn/9dpwzDUk08+qf7yl79466uddNJJynG29Q7pwOP111/3zrFkyRL19NNPq4suumhceeGcD0IIIaQE2EUednnkkUdE+q677vJ6QPQir0cddZQXbLz33nvqlVde8Xo1NHfffbeaNGmSF4wcd9xx6s033/TO8+KLL6pPfepT3j6//OUv1ec//3n1k5/8RE2dOnVMeeGwCyGEEFLG9PT0iI9eimQs6JXfNXrtNI0+Tvd6jF6cVbuM6yVOnnnmGS+9bNkyb6hlOPDQ6KBE7/P888+POc8MPgghhJAS4LhGwR+NXgG+vr5+5KPnduS9tuOoyy+/XB155JEjC7nqBVv14q56VfmBgQFvGOab3/ymsm1bbdiwwdunvb3d6y0ZTSQS8QIYvW2scNiFEEIIKQF2gavaDh+7bt26kWESzeiei+2h536sXLlypEdDoyeZ3n///eriiy9Wt956q9eb8eUvf1l98pOf9H7emTD4IIQQQsqYuro6EXzk45JLLhmZKDpt2jSxTU841YqXzZs3ez0aeoilra1N7bXXXt52/XNHR4c4JpfLeQoYvW2sVEzwURcfVJGPlu6ui8nxsJ5M/iixvkqud57NwZLrsBy9DVr2LC4bHZX7m1n/rGX06bATeZbYllYkPt8AtHUw0FcgwHsisSX8GpjHCHhTxLrkRaOwvHm8x3/RKHhsVLXLfQZaZX0ZDvh+xM1Q7wKjWm6P9vvzkGqMhi5PHuuGwoOyjaRgqXHwssBy804RkcdkA5ZAlydVoUuHGxFIw/51jVAZnv8IlE1EmkP09slGaIGvQ748BliL+MrOxSXRM2C6YULZQtqtlXWTseTxVr//LzgX34RYX3XoNSLP4Vpy/+Sm/B4bg83yHMktsr7sGLYZI9QXBr1IfAQII1JNefbB+sIK7JIvhGy1vNG6Zr/8Et+d7kdDB8PkLFkOsaZcaDt2sL3A+cyARofv59b63m3Xj6TValUcnFFDJzt6/HhwXVddeuml6oEHHvBWkJ81a9Z2921ubvb+1xNNdbDxhS98wUvPmzdPdXV1eZNUDznkkJF99DDOYYcdNua8VEzwQQghhEwkHGV6n0KOHw96qOXee+9VDz74oOf1MTxHQ88TSSaH/iJavHixmjNnjjcEoyeXXnbZZeqKK65Q++23n7ddb/vc5z6nLrzwQnXHHXd4cl3dk3LGGWeMWemiYfBBCCGElADbNbxPIcePh0WLFnn/H3300eJ7HXCcd9553s/adOzaa6/1hlFmzpypvvOd73jBx2juueceL+A49thjvbkgp59+ujdHZDww+CCEEEIqABfXSwjg5ptv9j5haGWL7kEpBAYfhBBCSAlwijznYyLB4IMQQggpAW6Bq9rq48uV8s05IYQQQsoS9nwQQgghJcBWhvcp5PhypWKCj/p4WkXjQ5NttqSqxLaelPQuSEazfg/8AWnMkIhnQ3XjOdCyYxuxq6WW3ckFNCL4ygVvEPQrsFJGuJcIjg9CFm1Ie9+BH4Ul7U6UBb4eFi4pAPObklvBLyHgmoYtyyZbHQkvS/D1QD8EK+OGei70TfE/BujTEeuVeeqfCoYnQO8MqAu4hA2eHBo3BoWVgIwOgrdMDbTTXHhHZuPkbV4Gmn0nbfbt82G/NCrqS0tPlUjUDvVQSPXHwj086vzPlgv5NmPyGi48W77joV27Wbl/vEU22ky/9HAJOsYcQG8RmbTB9yPrQt1AnqrX+yf6ZaJyn/5W8A4xx+fjkalXoc9JrjrA76IG2hh4ZqCniq9NQt3EauULYDDlf06qq+Q+ff2yjdXXwksmjxcNtkGcB9GQ8J8vasn7yI06p4umPrsQxy1s3gZYHJUVHHYhhBBCSFGpmJ4PQgghZCLhFDjhtJBjSw2DD0IIIaQEOMrwPoUcX64w+CCEEEIqwOF0IlG+fTaEEEIIKUvY80EIIYSUAIdzPnZ/utNxFYkMSbqippRZNVXJpcW7Uv71zltqpUxx66CU6+KK2bm0FS4xHAiXjw4dA8u22yDfrM6GStAiG6TMzU66oUu4Ox9JkcU5ukBCaIcvRR4Bqe1Aq8xzpk4eYKUDJIi18pqZWnkOVMLhUuTJDlluuWp5vt7p4ZJGTbZGXnNgsjyHLdXZyomGS5B96r2g3lKUNaYhnzF5EhP3j8rt8aRsHzGQF3Zl4CaUUrUxWYG9KSmDzGWtUOmtC9staPdugIrRgPuKJzLyGlb4feMS7fFoLvS+O1357GpyGXkOB2XLIDE1emU7ztXL7ZE+uX1gsr/C8VlJNfp2kdeEskOZup1wQ6XbdpW/8C14h2B9xmOyLCdVD4j03nVSrp0BPXBn2l/WMXj/bk3I921Hd61IVydlQSUiMk8oVcX3eSbAQ2C0tFYTMUeVzeifizHnw63MOR8cdiGEEEJIUeGwCyGEEFIC3ALVLvr4coXBByGEEFICnApe1ZbDLoQQQggpKuz5IIQQQkqAQ7ULIYQQQoobfBgcdiGEEEIIKQYVM+ySsSPKsYdu98iWd8W2TZkakW5OSC27JvXRsYG6cKVUe29t6DLfCnw/jKwR7umhAS8B9IFAh4xEvTSXcPaWfgkOLJdtp+T5jW7/UuMRWI0aPTHMbLhXAR6P3gQq5p8wlQP7iVgf3DccMtgkM5Xa1wxdatyEos7W+XX90V55jnQLHATLeGOmzEEoKMizix4dmhzslJD5iiRlHiI+DwZZ2BHwt0hGZGXNqtniy8Lq3maRrolnQtu1AeWQjYOPS0TmYXrTVt818XlD/5F8k+oGc7LdZnIyDza0e58/ileWMp8ZrN8B8O2Jy7oxqmRdpPZGQxz/PVgJ8KvYlAhtI5FeK9yXBy4R2UOWa00MHlal1L5N0qdjUkwe87GqjSKdAkObLdlqkU6Dz0dLXPojaToz8pgtKekFUpuU77E01Cf+1qqOyTbaEIPnwGeyo1RPVpb1mq3bXlz2gPRH2pU4XNuFEEIIIcXEqeBhl4rp+SCEEEImEk4FBx+U2hJCCCGkqLDngxBCCCkBTgX3fDD4IIQQQkqAU8HBB4ddCCGEEFJU2PNBCCGElAD3I7ltIceXKxUTfEyp7lXR6rT38ytbp4ltDTHwxwhoDKgd78okRToZlZr9ZEu3SKdAq57NSc2+Bb4hGrNZNq0Y+DZ0bGiQ1+iJyxNY4I+RkR1dBvZ7BTwDNnhuRHtkOgZpByTyNR+AV0VaprNV/s631CT5XX+rzFj1Rnlfg5Pl8ZlJcN8RN9TvxEn6y97niACeC/H6tLzmJtke1CTpPRCJy7rL9vm9BCLVcEwEvCTAe6Kptl+kp1RJT4VZ1dLDIQ4GJxvTdb48VIEXiAnXxHbcOyAbyKRJ/aHtvDedyNt1PK1aPjtJS+ap35Zl15mWPhGD4OvQn5H7R8F7RJOIymvIHChlw33YW+BZg+qvqpPvlGTc77GRs2U774NyyA3Iss5V26E+IejbU1uVCn1/BHlg2JCH7py8sUlRWb8JMPrpyNaF+mlo6iIyX/vWbxLpzrT0AekYkD5MVVH5nMyq7RTpiGGHeo8Etbm9G7d53mTjGfWGKg4Oh10IIYQQQopDxfR8EEIIIRMJp4J7Phh8EEIIISXAqeDgg2oXQgghhBQV9nwQQgghJcCp4J4PBh+EEEJICXBdw78C+jiPL1cqJvjozcZUJBsPlBMiGRvXfVdqzaBcK742KqWWk6v6QiWKW9PJ0Ig1bsGS7Z5EMB66NHhds7xmz2YpSYsm5H1moboTIO1MR+Vy2d45psn7zIHkMP2elMXVvifzGB2UUj4HijZd7x/5y0nlpHJA1dj5cZm2Yel5F9IqCw8oLqmOy6cHSGUNR+YzC7JGBUusOxl5o9msFVp33jUgHxbks7VGSmk/0fChSDdH5XbHlXk2QVbZZ0PBKqUa4/KYjsFakZ5WK0Wo3THZPtI2yD1hew7yNJQvN3SfuqiUZmZhOz7POagrXHK9LpHKm4fGqoFQSeo7sWaRTsBy9TMbtop0X9YvrW4DafT7iUkincxzX80JmceujJS1ZuE91grvqKAl7zfDcvcoz260pNR2dapF5tmUZR2HctFsytSGXmNGlSy7KUnZ5gZBao1tTsE7Jgp1p6mB9/douXYu599/V+EooyCfj0KOLTWc80EIIYSQolIxPR+EEELIRMLhnA9CCCGEFBO3gud8cNiFEEIIIUWFwy6EEEJICXA47EIIIYSQYuJy2IUQQgghpDhUzLDLpPigisKy5sN0DEp/jIEATX4N+BV0gW9Hc1Lq37th6XD08UDNfsT0a8uTsHR0BvTsNXGZp/2a5NLUH/TVyzw0yjx09Mr7PnCm9I0Ioge8RxItctn21bOk/0FXP5RlCrwnUv77duLSc8G1ZNpIw1SlCPh0gB+KUSe9Blwo+6pav+9DY430ULDhmCaob/SmQF+IRCQXujS9ZkpVD5xDtte2uNxebw2KdGeuOnS587d69xDpxpi8B01TTN53HyyJ3hSXXhGpPL4eDTGZxyAiprzPyTF5jeZoX2i5YB6bE/77CnveNS1JeY0YeE+gZ8qUGbIuBm3pkZO0ZNnPqPL7XfTkZL73rtu2rLsmBedEr5BOeAfVQdnH4H2HdaeZU71epLfGZRuaFJFlaQV4ZoymK1sV6skSVN91Efn8bUxLHxAEPZKwHaNvSBTaS1C+GmLb8pDNyvfuru75cIo44XThwoXqd7/7nXrrrbdUMplURxxxhPrhD3+o9ttvv5F92tvb1VVXXaUee+wx1dvb6237zne+o04//fSRfTo7O9Wll16q/vCHPyjTNL1tv/jFL1RNjf/Z2h6ccEoIIYSUANcLIAr4jPN6S5cuVQsWLFDPPfecF1xks1l1/PHHq/7+bQHcOeeco1atWqUeeugh9dprr6nTTjtNffGLX1SvvPLKyD5nnXWWev31171zLFmyRD399NPqoosuGldeKqbngxBCCKlkHnnkEZG+6667VEtLi1q+fLk66qijvO+effZZtWjRIvXpT3/aS1933XXqlltu8fY5+OCD1Ztvvumd58UXX1Sf+tSnvH1++ctfqs9//vPqJz/5iZo6deqY8sKeD0IIIaQEOB/Zqxfy0fT09IhPOi2H4LZHd/eQdX1j47blQ/RQzL/8y794QyuO46j77rtPpVIpdfTRR3vbly1bphoaGkYCD81xxx3nDb88//zzY753Bh+EEEJICdUubgEfzfTp01V9ff3IR8/tyIcOLC6//HJ15JFHqrlz5458/6//+q/ecExTU5OKx+Pqa1/7mnrggQfUPvvsMzInRPeWjCYSiXgBjN42VjjsQgghhJQAxzWUUcCE0+HJquvWrVN1dXUj3+ugIR967sfKlSvVM888I76//vrrVVdXl3r88cdVc3Oz+v3vf+/N+fjP//xPdcABB6idBYMPQgghpIypq6sTwUc+LrnkkpGJotOmTRv5fvXq1er//J//4wUlH//40PLhn/jEJ7zA47bbblN33HGHamtrUx0dHeJ8uVzOG6bR28pi2EVPajnwwANHCm7evHnqP/7jP0a263EmHZ3p7h8t4dFyno0bN5Yyy4QQQshOwXUL/4zveq4XeOhhlCeffFLNmjVLbB8YGJLa6/kbo7Esyxum0ejf07pnRE9AHUafS28/7LDDyqPnQ0dcN998s9p33329Qrn77rvVySef7El6dNR1xRVXqD/+8Y/q/vvv98axdKFp2c9f/vKXcV+rKTqgYrEhrX2fLfXyWwaktt22/TFZOsCXYTQbHalNj1hSWx5VMp11LJGuivh9ABqSUv+egWPqonK7CRr8+VOlhr850ivSqwamiPSgI30FNNMSW0W6PS2j602gyd9jUpe8hzaZx2nJrtA8e9dISX+SzSnpHdDRK6/Z2y39DiJxqfO3IvIaUxuGJlkNs2+d9CoJyhfeZ11U+ldEYX+sm6QpvQOCtP014HeQMOR91FvSg6PbluVSa8nj7Y8mow1zQK30cenOJfPmIVkr870hLesGaU3KNtYYlR4Mfba/Oxh9GGosOVkuDn4l6F8yq3pzqNcE+j7EwGciiJ6sLJt9a+RfepaS9d0DZZm0MqGeHpoo5MOGNlEVSYf6eqA3UAN4tNRAHvap8v/hNuDEQttY1rVC0+jBMiUmn2874O9b9EzptWXZzKluD63/7hx6icC7FNpPe0CbbY3LdpobVXZpx/8u3l0cThcsWKDuvfde9eCDD6ra2tqRORr696v2/Zg9e7Y3t0PP89DKFf2Hvx52GZbUaubMmaM+97nPqQsvvNDrCdHzQ/Tv5jPOOGPMSpeSBx8nnXSSSN90001eb4jWIOvA5M477/QK6phjjvG2L1682Ltxvf3www8vUa4JIYSQ8mPRokXe/8PKlWH079bzzjtPRaNR9fDDD6trrrnG+/3c19fnBSO6Y0BLaYe55557vIDj2GOPHTEZu/XWW8tzzodt214PhzY70d06uktHR1RawjOMjspmzJjhSX0YfBBCCCln3CL3fOgRhnzokYh///d/D91HK1t0x0AhlDz40A5qOtjQ8zv0vA49FrX//vurFStWqFgs5umJR9Pa2hoq59H65tEaZ615JoQQQnZXtUs5UnKfD+0brwMNbU5y8cUXq3PPPVe98cYbO3w+rW8erXfW+mdCCCGETBxKHnzo3g09pnTIIYd4gYOW9egFarRkJ5PJeLNqR6PVLmFynmuvvdZzbRv+aP0zIYQQUulql4lEyYMPRMt19LCJDkb05JcnnnhiZJte7Gbt2rXeMM320OYqw9Ld8WqfCSGEkGLhegFEIQ6n5VtXJZ3zoXspTjjhBG8SqV66V09geeqpp9Sf/vQnb8jkggsuUFdeeaU3uUUHEXoJXx14cLIpIYQQUr6UNPjQLml6+d4NGzZ4wYY2HNOBx9///d972/VKesMyHt0bMn/+fHX77bcXfN21fZNEuiYudeHVUamP1zTGpf59ZYf0yEAGMtIzI1krtePTa+Rw0l5Vfq+JTRnpLdEN3gPItIQ8Z0tUTrbdmpN+JtMTnaF6e00NeEe0gU1D0pL3dUCd9L9ATFgEelpM5kGzR1zex+uW1I5v6qsR6cnNPaH1OQPKGr1LWmL+ScnoHTA13h3qqVFlymta4PsRN2Q5NYCfgiYGvh7okVBryGuuyzWJ9Mas9DNohmsML0K1vbrVVMN9fOA2yu3godCckD4eaVu+UqZUy3JLuf5XDvqZoPcE+oDUW+FtDH1CouDzUQ0+EUHXSDsyn63wLKGHyiTwM0lAfadifg+dfH4VXfC8z23YIK8ZkfXbnpG9vJNj0suiL+D5npv8QJ4D2hD6E3WBt4zvOYnKZ6vX8V+zM1cT+iyhxw62B/R5wecCn8Va8EsaOmd8u9eIwjttd1K7TCRKGnxoH48wEomEZ+mqP4QQQsjuhPvRp5Djy5WSS20JIYSQSsSt4J6PCTfhlBBCCCG7N+z5IIQQQkqBW7njLgw+CCGEkFLgFjbsoo8vVzjsQgghhJCiUjE9H02xXpX4SO62b52MuTampEQtY0v5mMY0ZP/W1HopIUznZFGmctHQ4/dIdocuRa2pj0hJ4dLN+4p0U7wvVC6I0rrWaHfo/vslpJQvSFqH8kCU+6H0Du8LJYgoSdU0R6REcDLc5wEtMp8zk1tCl/mOgpQT5YCbcn4jur1icgl1B+L0BlPed60p5Xz9PnlguJw0UJ7pRkPzsEdE3kdbpCu07posKQfdYkvptcaCdhqN26FtpgkkpmtTUpr7er+USc+u9rexxoisLxPyjbJky/K3GZnHcMlyVcwvpcf6aI7KNmgpec29QCL+XqY5NE9B0up89d8ajYRKULEu8PlGyWka2pOmzZLH9IMEFZ9nLEtcWwSPzwRIq1HijWWb709klIM3WrL97AHPwdqcbJOa9dlJ278cPKu7ErdAl1KajBFCCCFknMGDQbULIYQQQkgxqJhhF0IIIWRC4RqFTRot4wmnDD4IIYSQEuBW8JwPql0IIYQQUlTY80EIIYSUApcmY4QQQggpZuzhVq7apWJ6PiZF+lUyEglcYhs9OHBJbs3GtPSCOKr5HZH+MN0g0hsG5dLUNRF5zjbwv1iblsujB+n4cfnyFlgyG5ennxzxLxUfpnWvjfqXKk850hvgwKq1It1rJ0OXbZ8Mnh02PCzoK6HZlKsNLasm8GDAcorDktvTo1tCvQtwSe8gHw/06WgHbxBcYt125YhmQuVCfUE0m2y51HiDKeujw64N9YVA/wMs2wzcNy5NHrSEOg7Moh/GxxLtoZ4dWfB5SAX4PqA3BPo+4H2hfwn6W2CbrbUGQ5dTD2qnH2QaQ9sYtmP0hcE8BrVzPCcuNY9+NOj7MjO6WaRXZ1pEekZUvg/WZv1+F3/LtIp0DPKEZYftutdNhD7/WHdB94m+LPXo6+LK9lBlZkLbMb7PPxaVdaOZCe+ErlHvsf4cPAO7GldVJJzzQQghhJCiUjE9H4QQQshEwuWwCyGEEEKKG32oil3VlsMuhBBCCCkqHHYhhBBCSoLx0aeQ48sTBh+EEEJIKXA57EIIIYQQUhQqpufjgPgHqjoxNMVlwI2F6sa35qp9x+9dtUmkP121WqTfjrSJ9L7JjlAvgl5b6uM3G9LDYSgf0ivANJxQnwD0r7Bgf5/3BOjjN4F3habRkr4NneBFgedAj4aEIcu2y8GylRr/IJ0/5hv9D9CvpCrAWyAsz02G9E/RZJUsy/fSzSL9sdhGkd4C94V53mRb4X4aSqkeR3qmmDCbDH09kF5HtqkqI7wc0A/jo4uG+lnMq347NM8IejigD0jQs4DtGvNp5mlj1eCp0wPnrwX/lCAfjinQxrrBtwO9SdDPogHacFDdrcs2javd4jnfyzaHtusVqRkivW9cerJo2rPSn6gOPG/wHYKeOe25+tBynGz15PX5sKHRbYRzoicOvuf8zxr65chy07SBf0nK2far0IZ73qW4ldvzUTHBByGEEDKhcCt3VVuqXQghhBBSVNjzQQghhJQA1x36FHJ8ucLggxBCCCkFbuXO+Rj3sMurr76qfvCDH6jbb79dbd4sFzbq6elRX/3qV3dm/gghhJDde86HW8CnEoKPRx99VH36059W9913n/rhD3+oZs+erf785z+PbB8cHFR33333rsgnIYQQQnYTxjXscsMNN6hvfvOb6qabblKu66of//jH6gtf+IK6//771ec+9zk1kVmba1TJ7NDt/l3yfbEtAXJAlKgFyb1QJrc3SC8zsHS4A3FeW6QrVCapmRGXsrZ3Uy2hEmGU6qFEMQUS431jUnq3yfZLbTFfkyNSOrcFZG24JHd7riE0j0FLjffa8dAl0FHui3nCZdx98mCQPaYMv+QU89UC17DBWRDPibJHbC+YZ80eka2hEkJ/Hp3Qa+By6Ai2B03WlmW3KScl4HvGNodKULGusA2ijDo4X9FQKeZmkIRvhLrYNy6fxV6QA2MbHbqmGVrfWJ94X9jG8NnsdGsCrinP4cBfsXhNlPeuB5nsADw3KO9/ZWCmLw/7QFm9MrCnSH8sId8RA45sM3Fo9yitDSrraTH5XkMcqIt8EmSU3vqebyhnzdvZxHbf1wOOP8+7CsMd+hRyfEUEH6+//rr653/+Z+9nwzDUt771LTVt2jT13/7bf/N6Qw499NBdlU9CCCFk96KC53yMK/iIx+Oqq0v+xX7mmWcq0zTVl770JfXTn/50Z+ePEEIIIbsZ4wo+DjroIG+OxyGHHCK+P+OMM7xhmHPPPXdn548QQgjZPXEr12RsXMHHxRdfrJ5++unAbV/+8pe9AORXv/rVzsobIYQQsvviVu6wy7jULqeeeqq65ZZbtrtdD8GMVr/89re/Vf39/nUzCCGEEFK57FJ79a997Wtq40Y5m5oQQgghalvPRyGfMmWXBh96GIYQQgghpQ8+Fi5c6KlSa2trVUtLizrllFPUqlWrRra/9957npI16KMtNYZZu3atOvHEE1VVVZV3nquuukrlcv4VysOoGHv1aZGtqjo6pAd/Od0Wui8uj63ZI9oZ6vuBngzo84HLeqNmv9HyD0+hDwfq19ErwuexAD4OuJT4a6npIr1XrMOXB8eQE5qqYfnyqCW9BNrt+tBlwDOgycd00DG4xLoFS6ojWLa4rDv6Riioq6FjZBtozONP0QW+Dw1Q1jOj0h+jP8BjI8iPIKz+7Tx+CLg/lktDQJvbBB4a6DWBfhZIc6Q31Bfig2yj7xj02ZkKy9mj1wzWZxzMDtCrogaWia81/W0OQc8UBP0t8HlP5cLrMsgjZVqsM/SdsSE7KfS52Cche5lNeE7Q9yPIUwN9PfBdmIDnH5+lARc9efzlYII/DV6jy67O44eC71orvM3m/G0W2/7oa6Dny+7E0qVL1YIFC7wARAcL3/72t9Xxxx+v3njjDVVdXa2mT5+uNmzYII75p3/6J8/T64QTTvDStm17gUdbW5t69tlnvf3POeccFY1G1T/8wz+MOS8VE3wQQgghlax2eeSRR0T6rrvu8nouli9fro466ihlWZYXVIzmgQceUF/84hdVTU3NiNO5DlYef/xx1dra6qlgb7zxRnX11Vd7RqSxmP8PqyB23xCPEEIImcAYbuGfQujuHuoNa2z090hqdFCyYsUKdcEFF4x8t2zZMnXAAQd4gccw8+fP99Z200akY4U9H4QQQkgZS217enp8hqD6E4bjOOryyy9XRx55pJo7d27gPnfeeaeaM2eOOuKII0a+a29vF4GHZjitt42VXdrzseeee3rjQIQQQgjZNUyfPl3V19ePfPTE0nzouR8rV670lkYJQi8Ue++994pej51JwT0ffX19XgQ1mrq6oUlr+sYIIYQQsutYt27dyO9dTb5ej0suuUQtWbLEMw3V67MF8W//9m9qYGDAm0w6Gj0n5IUXXhDfDVtq4HyRnd7zsWbNGm+2q54dq6OsSZMmeZ+Ghgbvf0IIIYSEYxQ67+Oj8+jAY/Rne8GHtr/QgYeeRPrkk0+qWbNmbTdveshFr1o/efJk8f28efPUa6+9pjo6tqkjH3vsMe+6+++//67t+fjKV77i3cSvf/1rb6xHa4AJIYQQMnFZsGCBN5Ty4IMPel4fw3M0dCdCMpkc2e+dd97xekUefvhh3zm0NFcHGWeffbb60Y9+5J3juuuu886dr8el4ODj1Vdf9WbB7rfffqpc6HaqVNaxAvXxSJUpteyaLXZtqL4dPRXezzSLNF4TvSxqLekLEeyZ0R/qydAU6ZN5BBkW6uUxD12O9B3w9jHlPh1QDkiLJX0e1oKvA2rwgzw78vl0pMBTA8+JvgFBfhajmQp+KZpN4DWAfhU94D2RBa+BTUp6C5gwq2xmgOfCivQUkW6y+kLrp8nsD62rAfASiUE5ojfFUD5lfdRb6dCyxbpCsO6Cnj30p2nPyXbfkZHt/GPJ9tBz4n1iHlen5WS5IH8SbFPY5moj0jukI1sXmqcgj42UK6+xPtsQ6mfhyzM8v43QXhBsox7wd2NUhRtF4X1tgvueHZceER2u/30Rg3P02tt+6QWB79aYKfPYAe/BLLTrtoi/7PFZGv0eyucjVM5S20WLFnn/H3300eL7xYsXq/POO28krTsW9HCMDjQQLcfVQzZ6rTfdC6JHQPSist///vfHlZcdCj60QYkeYyqn4IMQQgip5IXl3DG6jmuzsDDDMC0mCeoV2eXBx//9v/9Xff3rX1cffvihJ9FBRcuBBx5YUKYIIYQQsvuyQ8HHpk2b1OrVq9X5558/8p2e96GjKv2/tl8lhBBCyMTp+Sj74OOrX/2qOvjgg9Vvf/tbTjglhBBCdgCjQJfSQh1Oyy74eP/999VDDz2k9tlnn52fI0IIIYTs1uyQz8cxxxzjKV4IIYQQUuCwi1vAp5J6Pk466SR1xRVXeEYjeoEZnHCqjUkmGhtzdSqZjQQuRW9B31UTyEU1Lw/OFOlOszpUSlcLy3ijVG89LI/dHCD3bIHvUBo5PboldDlrvK9eOz4uad7QOaTsbFiuvL0l1nvNZKiMsg7KBaWYmlqVCpXF4Tlwme9GkByjxBile+05vxwQ9+l0akJlzFjWtabM49sZ6fy3R0B97xHZGrrUfL8j51KZWL8g/62GckGCZLIopU2BXNcGeSjKJrE+sRxRDj50Tvk3UBcsNd+QlJLSqb5yAslqbtK433oor/dJhC2QmEKbQakukk9OqunLyfprjoJ0HnSx2M7xWZwdk7LXVnhuNKvgPYTvGF/9gWwd5eDYZoP4EOoH21wGJMF4X9he8B2F79q3slLCHiSlHi0hTrlFnLPocs7HuNBKF02QrpcTTgkhhJD8GJzzMT5wLRdCCCGEkKItLEcIIYSQie9wWpbBx6233qouuugilUgkvJ/D+MY3vrEz8kYIIYTsvric85GXW265RZ111lle8KF/3h56zgeDD0IIIYQU3POxZs2awJ8JIYQQMn4MTjjNz5VXXjm2wjQM9dOf/pTtkBBCCAnD5bBLXl555RWRfvnll1UulxtZ2fZvf/ubt9TuIYccMiEb21a7Sg3aQ/pvG7zVUKtuG/5JPOtScmn4z9S9LdIfZJpCl6/GdDzPMvBBXhGoof8wJ5fgrjLkORzwBZgd6w71iWgKyMN7sFw1LqmNecSynR7bErrcNWryvXMYZuhS8eg1sC7XGJqHPSx531vANySr/EuNxxTe56BI9zrSt6ENroHeIvvH14v0h1CuQfe1CfwN2sB/pgd8OmqNPP4oUFdBy9tjWaNnBt431l8D3AO2D/SiCTrHjGinSPeD1wjWF3rm4H2hh0eQvw16nuB9Yv1jO8bnG/PQHeDzgfedg3NshOXq0468ZiOUdactvYdMWBq+w5J1G+TjgaB3DDIDlqvPwLOXgnIK8lBpNPtCvWNWZ1vgeDvUcwefvW7wjdFk3fR2zzno5PcqIYUz5mGXP//5zyM//+xnP1O1tbXq7rvvVpMmDT34W7du9Raa+7u/+7udkC1CCCFkN8ctcH0Wt8Ls1fWwysKFC0cCD43++Qc/+AGHXAghhJCx4FauvfoOBR89PT1q06ZNvu/1d7294TbDhBBCCKlsdij4OPXUU70hlt/97nfqgw8+8D7//u//ri644AJ12mmn7fxcEkIIIbsbbuX2fOyQw+kdd9yhvvnNb6ozzzxTZbNDE5YikYgXfPz4xz/e2XkkhBBCdjsMSm3HR1VVlbr99tu9QGP16tXed3vvvbeqrpazrQkhhBBCdsqwyzA62DjwwAO9z44EHnrS6qGHHuopZ1paWtQpp5yiVq1aJfZJpVJqwYIFqqmpSdXU1KjTTz9dbdy4sZBsE0IIIaRSF5ZbunSpF1joAER7hnz7299Wxx9/vHrjjTdGgpkrrrhC/fGPf1T333+/qq+vV5dccok3r+Qvf/nLuK7luKb3CdLc99lSy14foIdvi/WI9OsDe4j0HvGtIj2Qi4V6DbRG0RdC7h/kJYDEQC+P+nbcnnHD/RHWg/dIkCbfMpzQ7SbEs9WGvO8M+GeklN/nAz1P8L62gJ8B+gTUwfEd4JexR0TWZT+cP4guO9zXwwZPlQbwhag1ZTkNOLG8ngrYJtCfJApl6UDZd4G/wQB4WWDdeXlwE6H72JCnyVCW6BOCeVqf83suoIcGttN8zwF6VaDPC97DpgCPFfTpwDaH9ZUGj44BOxbqoTMYUN/ow7EpLdvptIR8pzjgZ/FK354ivV/VBpFe1ruPSO+V9IsE9op1iPT0qLzme9mmUA+Ofqgb9B7amKv3XXMqXKPdlvsMv6eH+SAjfXxqLOkdY0E5Yjvvhfd7oFfIqDaTLuZibS5NxkrCI488ItJ33XWX1wOyfPlyddRRR6nu7m515513qnvvvVcdc8wx3j6LFy9Wc+bMUc8995w6/PDDS5NxQgghpECMCp7zUdCwy85GBxuaxsahSFcHIXpC63HHHTeyz+zZs9WMGTPUsmXLSpZPQgghhJTpsMtoHMdRl19+uTryyCPV3Llzve/a29tVLBZTDQ2yK6+1tdXbFkQ6nfY+oz1JCCGEkAmJqyqSCdPzoed+rFy5Ut13330FnUdPYtVzQ4Y/06dP32l5JIQQQnYabuX6fEyI4ENPIl2yZIm3fsy0adNGvm9ra1OZTEZ1dXWJ/bXaRW8L4tprr/WGb4Y/69at2+X5J4QQQkiZBB+u63qBxwMPPKCefPJJNWvWLLFdr5AbjUbVE088MfKdluKuXbtWzZs3L/Cc8Xhc1dXViQ8hhBAyUSecGgV8ypVIqYdatJLlwQcf9Lw+hudx6OGSZDLp/a9dU6+88kpvEqoOJC699FIv8KDShRBCSFnjUmpbEhYtWuT9f/TRR4vvtZz2vPPO836+5ZZblGmanrmYnkg6f/58z111vHi6/I+0+Qnws6gC3fhbg1MCjpdxWndW+j7UROQ5kMaI9KJIgU/A9OiWvP4FXY685ha7NtQfIwMafEyjl0GVksd718jVhF6jF/LUAB4pXeBNYYOGvtb0lxv6MDRY/aFeI+g9gH4YMyLSk6MLfADQkyHIO6DF6gs9ps2S6U5bplPgXdADnh5BnhgJ8EjZAn4l2D5aLLmo48diG0P9FNCTQ/NhblJoOWAawTaF3iWNUI5BYJvB+kXQx6ML2g96OmAeNVtz1aFeEujrgd4R+E5BBmz/NS0YsG+IDoR6i0RNWfamkve9OtUSet8dmbq8ecD7wDbWAe8c9CrBNozvnKB3Sh2UtWz12iMpHdoGu8HPBj1b6iP+dr45K/MQH+XDk8qF1yXZDXo+9LBLPhKJhLrtttu8DyGEELK7YFSwz8eEkdoSQgghFYXLYRdCCCGEMPioHKktIYQQQioHDrsQQgghJcDgnA9CCCGEFBWXcz4qijjIyYKWN0dmJqQUtsOSsrU1A5NFempCurKuz8r1aarMTF45IS7T3g8SUZScoTQTZW5tsPw5sj4rZZZB0tlOkHuiVLYrj+zNBGkfnk8TB3kfLi1vgiT4Q1uWbRtIa+thcNEGMd8WW8qFNYlR0ruhYySOI+smasn9m0F6u9kOlyRqXuqXJnuz4ptC5dgoOe6HcsKyrzLSoftrmkCui2X/XqZ5XHlCsH0ElUWtORgqz0Y5aJddF3rfuD8u2a5pjsj77gQ5qAPP4rRYp0hvzPqXjh9NleWXb0bgmcfneUsW5L8R2W5nJjaL9Jv9U0W6JyflwFMS8rkIKqu9I/K+TFhd/o1MS6gsGmXPM2OyDWt64XnD99qmXLiFgF86DWULTXDAludHaa2vjcA2smvgsAshhBBSClz2fBBCCCGkiBgVPOeDahdCCCGEFBUOuxBCCCGlwOWwCyGEEEKKiMFhF0IIIYSQ4sA5H4QQQkgph13cAj7jYOHCherQQw9VtbW1qqWlRZ1yyilq1apVvv2WLVumjjnmGFVdXa3q6urUUUcdpQYHt0ngOzs71VlnneVta2hoUBdccIHq68u/YnVFzvnos+MqG7CstabKkjrySWNYgvnk+uUi/WD3ISK9ZkD6IexdJfXubaB/x6WpgzT06MOBXiHob4F6+C7Q16N+Hj09gjT3eE3U6G+E5cxx/z5Yihy3e99F5XcpWM48Y4CPA/iCRMFrYsC3erI0L8gqeT5NkyG9Jj6E+9oDPFPqTOnJ8FpGHt/l1IT6ZWgaI/2hnhkrU9NFep94e+hS8egbge0jG7Dc+Ye5SaG+DfvCNdvBv8axZdna4KnRGPG/oNBTA+vThHJ4LzM5dH9sUz6/E2j3QWWNy7Bvzcn63ZCR910TkX43feCx0Zfze6rsldwc6j+EPh9r+ptEehK0l4aozHNXVj7vHw7KPGvScdkGZkZkvj+0B0PbA3p2NFqyfjeBB0vQOTKOFerDgu0D6c4lQ32bcnB+TVOsb7vvmFS4VU1Zz/lYunSpWrBggReA5HI59e1vf1sdf/zx6o033vACjeHA43Of+5y69tpr1S9/+UsViUTUq6++qkxzW73owGPDhg3qscceU9lsVp1//vnqoosuUvfee++Y81IxwQchhBAykTB8fwqN//jx8Mgjj4j0XXfd5fWALF++3Ovd0FxxxRXqG9/4hrrmmmtG9ttvv/1Gfn7zzTe987z44ovqU5/6lPedDlI+//nPq5/85Cdq6lRpdrc9OOxCCCGElDE9PT3ik077e/eC6O4e6g1tbGz0/u/o6FDPP/+8F5AcccQRqrW1VX32s59VzzzzzMgxumdED7UMBx6a4447zusZ0ceOFQYfhBBCSBnP+Zg+fbqqr68f+ei5HflwHEddfvnl6sgjj1Rz5871vnv33Xe9/2+44QZ14YUXej0cn/zkJ9Wxxx6r3n77bW9be3u7F5yMRg/N6ABGbxsrHHYhhBBCylhqu27dOm/y5zDxuH89G0TP/Vi5cqXo1dABieZrX/uaN49Dc/DBB6snnnhC/frXvx5TUDNWGHwQQgghZUxdXZ0IPvJxySWXqCVLlqinn35aTZs2beT7KVOmeP/vv//+Yv85c+aotWvXej+3tbV5wzOj0ZNXtQJGbxsrHHYhhBBCKkBq67quF3g88MAD6sknn1SzZsmVtGfOnOlNGEX57d/+9je15557ej/PmzdPdXV1eZNUh9Hn0r0mhx122Jjzwp4PQgghpFS4xbuUHmrRctgHH3zQ8/oYnqOh54kkk0llGIa66qqr1Pe+9z31iU98Qh100EHq7rvvVm+99Zb6t3/7t5FeEC3F1XNC7rjjDk9qqwOaM844Y8xKl4oNPtJOJNQHAP0WNM3RXpEecKWWfN+EnGizR0yOub07ODn0fKiPD/Jp2Csuu7qQtdmhGcvDxMDnw4aOLrzvXvDg0JgwIIk+HgnQ7E+OyPt6Py39LFqj0mui1/Ffcz14RyDoX1JrSo+FTZDHJvCWsF0pUGvP+a/3FqQPSbwn0lOt8E7DDHiHZKCssZw0neAlUWuCx0JU3jcyAJ4rM6ObQ8saPVqC2gR6rKDXDN6H31MDPTf8rxz0l1kH7RivEYM2h9fE8+F9d9vSP0djg2gR2zV6aqCXBN5n2pTl1hCVdRlY31Yq1O/ChN9Sawbls9UPXiIDkE5Y8p40U2Lyefw9eInsBX5E6NuC3jLvZ5vz+vhgG/owDf5FVibUxwPPiR4rAxl531nIs6YjI/2LakZ5PaXt3XdAYNGiRd7/Rx99tPh+8eLF6rzzzvN+1pNQU6mUJ7nVQyk6CNF+HnvvvffI/vfcc48XcOiJqFrlcvrpp6tbb711XHmpyOCDEEIIqbS1XVyf4WIw2uNjtM8HopUt4zEUC4LBByGEEFIK3Mpd1Xb37V8ihBBCyISEPR+EEEJIBQy7TCQYfBBCCCGlwK3cYRcGH4QQQkgJMNjzsfsTN3Mqbg7J6db0SznY9KqtoVLOIDnYock18hhYahyX/Z4c6w2V0Xba4ctGj2V5e5S9OSApRSkmygWbA+Sf3bBkNi4V7oAsbm26MTQP+WTP3jEwFSlth8fIKJPUdR1WnygfDqI/J8sqDZLTTHJoDYRhpkPZ1RkyT+tdKeetDljWPWtZofneBGXvk4iDXHuLLaWcvY6syyBwuXmULeNS8pjHzVloo6MkjNuTXvrvQ8pau0AaixJUbKOrU3LdiUnR/tDjgzBN+Xx2ZOpC5Z1Yn1ZEHr8Znl3NByn5TnGiss3URUF6C+08C0vF92Vlm52S7Ak9n+a13j1E+hMt74t0Cp5Py3BCpbetke7Q9hPUhlBai+8hlHu/NyjlwI1RKa2OmLJNRiHPmo60rI9+a1vZZdL+Nkp2Puz5IIQQQkqBy2EXQgghhDD4KAqU2hJCCCGkqHDYhRBCCCkBBiecEkIIIaSouJU754PDLoQQQggpKhx2IYQQQkqA4brep5Djy5WKCT7eH2hSUSMW6D3xSqfUuvfU+pd5n5aQXiAvDs4S6SbwWEBNPi6X/UFG+mFsNOp918RzoKa+z5a6frwv9DPApaXrwKtgRWqGLw+4FDguR49eIZgH9NQI8vVA+mAp8Hz3EaTjH00HeE9Ug69AkP+BCefckJH188/9R4r0QbVrRfrYqlUi3WAOhHpwBNXvXwdlfdRAG9qYlnlaoyaH+iNgG2yNSk8GzfuZ5ryeN6NZ0TsjtNywjfWAV4mmOyt9OvZIynyZ4ImzISXvuxb9MKANdmb9ZY3U4rLsNnjgRPtC/VD6wUOnz/bfJ1KNHijQLt/snSLSCSsr0jloL5jeM7kF8uz3r5gZ3yTSPeAFs8WtDX3ee8FjBbcHgfWDfiX4bsT6joGPT2e2KrQug9452E5LhsthF0IIIYSQolAxPR+EEELIRMKg2oUQQgghRcWt3GEX9nwQQgghJcCo4J4PSm0JIYQQUlTY80EIIYSUApfDLoQQQggpIkYFD7tUTM/H1nRCRSLxQO+JrpTUqv9n116+42e3dIh0U7xfpKfEpTeBBTryuJkN9U/oAb289102EapXR718xJTn7MpI/ftATvo+NMal90QMjtes6W8S6ZqI9CZIgjfBIPgjoPcAemwEafDjoOPHffCcHwzWyDxG5TUiUBd92Xiob4h3DfAeSNmy7GKWzOMfNx4g0h82TBLpY+tel+dz/X4Im3PSU6EjI9NbTelX8Xav9PXAdt0Qk/X7vi39E9ZFZR6DzoHlgGzJVId6MPSCr0dnWrZJTQLKclVPi0inbFn/WVvmyTKd0HvoTcs8VMdkG9b0pOQ+U2t7RHpOXXvo84tt1Pd85/zPdx88K+iB0pmWx6Rt2R5qo/I+co5sxyu6p4v0x+vW+/Kwd3xjqDfM64PTRHrdoGwzM8FLpD4ifYE6c36PFayfLeDD0p+Tz2cc2sdWeK8hH/ZJX5CmpHwONH3Z7fuR5FL+9kF2PhUTfBBCCCETCpfDLoQQQggpMkYZD50UAtUuhBBCCCkqHHYhhBBCSoHrDn0KOb5MYfBBCCGElACjgtUuHHYhhBBCSFGpmJ6PnlRSWdaQhKuzR0q1GuukFCub9csLt6TkMW9ubA29XlOtlOJOq5FS3HpYqn5TWspFg2WtUq77WueUUEnaQIicTNM5UBW6tLWmOiZlq7YjZWyJCEiILTtUFonXQMmqxgUJMUrlMnBOlBz75KIgpY2ANHNtn19yinLdnoyU/yHt3XUivWWwOnS59KCl6juy8hxZV97n2r7G0PvMJ5PNgBx0TZ+UUQeV5ZYBuI9YuIy5Pp4KlW/j+TWboKzwPrAN5WxZn201vSLdC3WF10yB5FxjQLZWfiifrQ9660Pb6Kdb3xfpTWkpi62D512zpqcp9HlN56Cdg5Q2FZP30TMgpboD8A6anJDlpNkKy9F3Qboa3kEmLCbytz75HpwE8u4gKT3WR3O8L/SYQZC5bxmUeUxG5TuoGp7dNV3yudEY+OyMui97MFxevlNxqXYhhBBCSBExnKFPIceXKxXT80EIIYRMKNzK7fngnA9CCCGEFBX2fBBCCCElwKhgtQuDD0IIIaQUuJXr88FhF0IIIYQUFfZ8EEIIISXA4LDL7k9vKqYsc0j/b4OPB/pG1FT7l1Tu6Ja6/VwuvNOoa0Auh22DRr86VhW6rLjm3QG/D8NoBjJS/542ZSzZNwBLx6fldsOUXXZOxq9v74J9TAu0XeCPEIvL+6hNSt+HnkHpRRBEMp4N9T+wwSegFpZIR9+IvnQ81OcD/RS8ayRlfXWlZH02VfWH+j5gm/rDBweI9NFT3vZd01QyX93ZZKjXyFbwO5hcLf0S2vtlmx3MyvYSA08W77tILtRToxvadX3VYKgfRiN4tKRt/987+GxkII1LxWP9rVw7VaSrasKXRI8E3PfAoCzb+lp5Xw2JwdCyR8+UPqirHLQnjQMPD/agW/DsDcLz29crn6VYUj43m3vlc/B4z36+PEyuk22mCjwyplV3iXRjTLb7npzMQ100lbe++2z5PPfnZFl1pmW+u9PyGj150tg+BtNBvi7udp9fO1XEoQyXahdCCCGEkKLAYRdCCCGkBBgcdiGEEEJIUXGpdiGEEEIIKQocdiGEEEJKgFHBwy70+SCEEEJKgbsTPuNg4cKF6tBDD1W1tbWqpaVFnXLKKWrVqlVin6OPPloZhiE+X//618U+a9euVSeeeKKqqqryznPVVVepXM6v2JywwcfTTz+tTjrpJDV16lTvBn//+9+L7a7rqu9+97tqypQpKplMquOOO069/bZfokgIIYSUa8+HUcBnPCxdulQtWLBAPffcc+qxxx5T2WxWHX/88aq/X0qoL7zwQrVhw4aRz49+9KORbbZte4FHJpNRzz77rLr77rvVXXfd5f2uLpthF33Dn/jEJ9RXv/pVddppp/m26xu+9dZbvZubNWuWuv7669X8+fPVG2+8oRKJ/H4Ro0n1JZRpf3SMLfX1W9ulH4IR869TbEbgO6h0B/wQ0PdhcEBq212nTl4TNP3eObPgb5CQkaUNvhxOWqYNuE/lhvsKGFnYX+8Tl/dtWzJPhiVPkoZrDnZKXwgjJ7dbk/yeDFvhmHiD9A5ABlKybLPg4xIH75HuQbl/PCG9DTQbuutCvUe2DEgvgoZq6QPRvqFBnhDK6RlrL981967bItIxU+a7BjwYNvXVyGv21ob6RGTBzwT9EIL8aSbX9IX6NgxkZFk62O5z0dDnIsh/BPOJ9ZkZlPtbUfB16Jd5isSkr4cNz6omB89S30eeQNvL0wC0IWyD6QGZx40xWTdBfkIm/CZBfwrXMULvC99B+E7Bewyq72wcyh78atBTB/2JOgZr83rJbAXPHPTc+HBrQ2i5YHvAconC847tJeiY0e9zZ7CMxzLy8Mgjj4i0Dhp0z8Xy5cvVUUcdNfK97tFoa2sLPMejjz7q/Q5+/PHHVWtrqzrooIPUjTfeqK6++mp1ww03qFhMPgsTsufjhBNOUD/4wQ/Uqaee6tumez1+/vOfq+uuu06dfPLJ6sADD1T/7//9P7V+/XpfDwkhhBBSdjhu4Z8C6O7u9v5vbGwU399zzz2qublZzZ07V1177bVqYGCbWeCyZcvUAQcc4AUew+hOgZ6eHvX666+X/4TTNWvWqPb2dm+oZZj6+np12GGHeTd/xhlnlDR/hBBCyERwOO3p6RFfx+Nx7xOG4zjq8ssvV0ceeaQXZAxz5plnqj333NObDvHXv/7V69HQ80J+97vfedv17+XRgYdmOK23lX3wMXwTQTcZdoPpdNr7DIOVQgghhOxOTJ8+XaS/973veUMgYei5HytXrlTPPPOM+P6iiy4a+Vn3cOg5l8cee6xavXq12nvvvXdanids8LGj6Nm8//t//+9SZ4MQQggJxShQLjs8c2XdunWqrm7bPLV8vR6XXHKJWrJkiSf6mDZtWui+erRB884773jBh54L8sILL4h9Nm7c6P2/vXkiZSW1Hb6J4ZsaRqfDblCPT+lxrOGPrhRCCCFkwjqcugV89IJ+dXXis73gQ8+l1IHHAw88oJ588klPyJGPFStWeP/rHhDNvHnz1GuvvaY6OjpG9tHKGX3d/fffv/yDD10oOsh44oknxBDK888/79389tCFjhVBCCGEVDoLFixQv/nNb9S9997reX3oKQz6Mzg4pNbTQytauaLVL++995566KGH1DnnnOMpYbToQ6OluTrIOPvss9Wrr76q/vSnP3nCEH3ufD0uE2bYpa+vz+vKGT3JVEdZeubtjBkzvMkwWg2z7777jkht9SQYbYxCCCGElDNGkR1OFy1aNGIkNprFixer8847z5PJagmtVppqKww9l+T000/3gothLMvyhmwuvvhiryOgurpanXvuuer73//+uPJS0uDjpZdeUv/lv/yXkfSVV17p/a9vROuPv/Wtb3kFoCfAdHV1qc985jOeTnm8Hh8aN2soN/LRCBl4DZiDoI/vCyiWllSopt7slNpmOwEeCmgTEnVD/S+8fSKg00+BTh8OMTLQkQXHqwAvEXE9fxaUypdP2O6i1wD4W5jV2VCvkqBjfDr9EI2+dzz4PvTD8Sb4Iwz0+ttTBLwCBpxw7fqAioWWNXoZdPZV+a9pOKGeCn1ZeY0cljUUS8ySviCDtiyHTVv83hN19dskdZoP8nguxKKynAbB7wLJgEeDJgWeGXa/zKcB9RWDNpSB/fG5yKbheYf25Z0zmQ31gUBfD3tQviNs8NiJJmW5ZPv85dKTAx+fuLzPDBxjRu1QzxT09UDfHwW+QZr+HvDIaZZ+NdXgw5Nx5Dk/7K0X6UlJ2X7WdYPfjVKqtzsZeh8W1HcG7gM9lxwox3QG/voOep9X+f1HRgjwgZnoapexooddwtDBhjYiy4dWwzz88MOqEEoafOjoK6wwtOupjqbGG1ERQgghZOKy26ldCCGEkHLAcF3vU8jx5QqDD0IIIaQUOP4h+XEfX6Yw+CCEEEJKgFHBPR8TVmpLCCGEkN0T9nwQQgghpcAtrtplIlExwYc5aClTDUm23JgcKDNAVuck/TIsawPIMRulNM+pz4bK/VyQuRlpmbZADqixk+Mc0ENpbR5cK7/cV8F3KMc1BkDOlwA5INwCShQNO0BijDLVgci47kNhniDpWlD/ARJklFLjfZhwDkwnajKhS6ynYbl0zRazWqS7YOnxfD2slumELnePS8n7lhXX1+yQ8ttIVS5UGp2C+8ZzZlORvGXtbpHSSAufx6xMZ6A9GLDdrZFt0KqWdWGBVFOTy7NMu4vtFO8DyjbbCe8LkN5qnHQkXDIM2PAOMeGcDt4DSFCDlPQulEUG8v1uV2vou9Kqk2W7Zb2U3hrwrg3ChXeM3S/LxQXpvE9CjHXjjuG9iPYKo85hDPrl4LsMd5tL6Q4fX6Zw2IUQQgghRaViej4IIYSQSnY4nUgw+CCEEEJKgcthF0IIIYSQosCeD0IIIaQEGM7Qp5DjyxUGH4QQQkgpcDnsQgghhBBSFCqz5wP08k7cCfXg0NhxWBq+W2rynYQTOg3ZwO3oTRBQE2YG8gnnQI2+L994G3BNC/wwHNDTD30JGnrcB68B5zRQMg/Lhvs0+p5vgxnqHWL0QmEl3PDlsoMMDkbhBviIRBtTch8oBxs8FGKxXKj8PprA7f5MZXKysOLQhjJZ8IUAn4hEMhO6VH2QrweCvivoy4L+Fk4KljuH+nfrwd9iq7+sXahuAy0xIN++rma8LchDLi29K+wa8OTxygY8NDrl8+3Cs2dmwIskBs87Pt9mgHcE7gPX8D0H+fwudkD54LsGgO8Ia1DmOZcyw589N7+vR9D7Fj2axP75fD2gjQY8aj6ViDPKC8TnG7QrcWkyRgghhJAiYlTw2i6V2fNBCCGElBqXcz4IIYQQQooCez4IIYSQUuDqCScFHl+mMPgghBBCSoBRwXM+uLAcIYQQQooKez4IIYSQUuAG6PLHe3yZUjHBh9Z6D+u90SfArZHadFfaAgx9B9ry0brwQO06+EAo1MOD9twJsAEwbfAOAG8BMxUJ1eBna8GbACT4eA/oK+Ltk8zvgSLOAb4BxgBs3xpwo3lwImaoJ4rPvgLK3gBvEbc7FurZoMmmquQ54CIueAn0W7Hw/etyoe3FOwZ8HNJOXKQjyWzoi2egOxm6XaEvxBhw0TJlEMrWDPfs8PllBHkoQFniw2HmwsvehXZsYNlC2uyQ5RqUbzynic8v5hiqxlf/Ac8W3rcB+/jvA/wuoNjcKHyBdRfgqeOCPwnWjwMTElwL7isG7wfwffH5AunbhncIlj2+x1zZhHz1j+9SvM8gC3IH8i3Okd8OZ+fhUu1CCCGEEFIUKqbngxBCCJlQOAX2tHBhOUIIIYSMB6OC1S7s+SCEEEJKgcs5H4QQQgghRYE9H4QQQkgpcCu356Nigg8zayjzI5mYARI0c72USQZhJ8Mlaaj+8i/jDPIvbDMBy53jNSL9UpOW2Cz3j/XIdC4h98/UwyVBwhYZ9GVBZepgiXRUa8KNJDfI+4j1yt2TnbLwzZz/4cGys+PmuKS4TkSeINUkm3m0D64XoEBFuZ6VlvmM9kP9Q5sabJLHD06WbSwHqlhNthalz7C8fVUkXLaaR3LoxmF2WoD807cPNmyQWqoE3Dgswe7E4fggxSnIWFEy6rtPfFZgO8ri8VmzayHPmgDpszhntTzG6LdCyy3gcfZhZLGhG+H3jVWBeYiMU4Ic8LvLwVdhnjyZUN/YBl1nDG0ML2niSUJ3zyvFDprQ6ZPjjj5HHjuBnYpbucEHHU4JIYQQUlQqpueDEEIImVA4lNoSQgghpIgYFSy15bALIYQQQooKh10IIYSQUuBW7oRTBh+EEEJIKXDcAOnjOI8vUzjsQgghhJCiUjE9H9EeQ1kf+SbEuuU2K6VCPR2ClpJGL4l0Q/j10UMj1guafNS2K6VqNkhvgRx6JgDomYF5Nt+T23unyeq3/SuNqxrwDokMynNUr8e1xMHnY6ssXHMgE7q/l+9YJLRr0QVfDzcivQbspEzXrbFDjw/668FOyHNYKTt0O5Z1vDu8vdSu9XsdZKtlvrI1kE3Id6oJ8lwFy5uD2YTVZYX6vHjfYTuEtAHeIyoTCfWSMHMqr9cEeoE44POBS7D7lkgftELLwefjEmTCgcu690EbA08U9JJAbxoDvSLQeyLgvtCPArc7sXD/E1/ZuvnLHj2PHPAG8f1Vjo8O1JXvPhP+du7zYcnl8W3Bohw0w8s+Txv09oE8jL4PX/valbgcdiGEEEJIUXELnLdRvsMuFdPzQQghhEwo3Mrt+eCcD0IIIYQUFfZ8EEIIIaXA0T0Xlal2YfBBCCGElALXGfoUcnyZwmEXQgghpAJYuHChOvTQQ1Vtba1qaWlRp5xyilq1alXgvq7rqhNOOEEZhqF+//vfi21r165VJ554oqqqqvLOc9VVV6lcLkBWFAKDD0IIIaSUE07dAj7jYOnSpWrBggXqueeeU4899pjKZrPq+OOPV/39/b59f/7zn3uBB2Lbthd4ZDIZ9eyzz6q7775b3XXXXeq73/3uuPJSMcMuDe84KhId6qJywJMhsUVGbNkaNAbQ2m/ZvRXvlsfgOdG/It0gz5nskP4YTjTA7yIqY8NEpzwmVx0J1cMb4PuB6Xg3+kL4suDzPDFgjDHaI/Nk5ORJrHUb5Qkj0OQsf/xrRKUBhZuMibS5bos8oCoptzfVyfNlZF0ZfWDsEg8wvID6s+sSMp2U92FloCxtWU6RfvAegO1eNrbKfbK1ss2kJpmhHg12Anwh4LZysphUpM+XBZUL8JsR18waoR4aLvhEGKn8f99Yg+BHMgAnzXcKKMp4J/q+4PUCzgG37UAzRa8I9OUxs2aoZ060fzsrmo5OBjRDsT2O3hRyuy0fE7+/ScAvKhPq0xqANgaHRAag/tH2ZRI8B13+XzFWCp6tJL5jwtsH+p3g+aw0HD/gy4Kq2rT991ou66r31O455+ORRx4RaR006J6L5cuXq6OOOmrk+xUrVqif/vSn6qWXXlJTpkwRxzz66KPqjTfeUI8//rhqbW1VBx10kLrxxhvV1VdfrW644QYVi0FD3A7s+SCEEELKmJ6eHvFJpyEC2w7d3UOOm42NjSPfDQwMqDPPPFPddtttqq2tzXfMsmXL1AEHHOAFHsPMnz/fu+7rr78+5jwz+CCEEELKeNhl+vTpqr6+fuSj53bkw3Ecdfnll6sjjzxSzZ07d+T7K664Qh1xxBHq5JNPDjyuvb1dBB6a4bTeNlYqZtiFEEIImVC4BRqFfXTounXrVF3dtuHmeDxgrQxAz/1YuXKleuaZZ0a+e+ihh9STTz6pXnnlFbWrYc8HIYQQUsbU1dWJT77g45JLLlFLlixRf/7zn9W0adNGvteBx+rVq1VDQ4OKRCLeR3P66aero48+2vtZD8Vs3Cjn8g2ng4ZptgeDD0IIIaQC1C6u63qBxwMPPOAFGrNmzRLbr7nmGvXXv/7Vm3A6/NHccsstavHixd7P8+bNU6+99prq6OgYOU4rZ3TQs//++485Lxx2IYQQQkqBo6U9BRiFgQpzLEMt9957r3rwwQc9r4/hORp6nkgymfR6LoJ6L2bMmDESqGhprg4yzj77bPWjH/3IO8d1113nnXsswz2VF3yMUjRZGYgWrfwyyMiAEypbxS4kw5b7R0FqibJIF5Zc947plfo+qw+WrwfsmBm6DLyZlXlKvi81h5m2Wt85UVrrwDUiqz+UB5i4NjksPd/TK8/f0uy7Zm71GnmKA+fIPEyXk50U1JfZI+/LTYB0t0bKZt1YJK/UNlct5WPRnkzo/vYkuX+kV9ZdutEvR0O5dbpOlmUGqsfKhr+Hon3h0stctb+dW/3hy7yjKtAnQYVywGtaoHIeyodMx4Ym4G87hxmejoB0FuWfic3heQ6SLUdS4cu6owQ9W22GHo/S66B9UEqN8v0slFOuBs4H9W3CsvEotQ6qn2gvWgaoUAkxymLjW8Ll3oEy5ky4VBbrBuXBCLaxaIDUNlsjzxHvckMtB3aXheUWLVrk/T88hDKM7tU477zzxnQOy7K8IZuLL77Y6wWprq5W5557rvr+978/rrxUTvBBCCGEVDDuDgQ6Qcfsueee6uGHHy4oLww+CCGEkFLgFrfnYyLB4IMQQggpBU7lrmpLtQshhBBCigp7PgghhJAS4LqO9ynk+HKFwQchhBBSCly3sKGTMp7zwWEXQgghhBSViun5cGKGsmND2u5Yj+yqytTIGCwKnh6aVKMUxFevt0N9PzAi9W0HDb81GNB9BhFxpjHcwCVXbYUuh21CHgb3rBfpqrc3B5xUmiI4jbWheXRbGkQ6PUUubx/f2BfqC+Gd8rOfFGljUBpa2InwZhvrkl4ibrX09XDAs8PIgJGAznerNEWI9MlyMOGYbIOsGzMj69NOWHm9ZKws1HetLJsYeDCgh0K8S6YHm8P9EuJb/GVvZvIsy55vf6ga9McwA6xqIuDDgPeFvg84P8+CPPjARzPAJgK/8vuX5PG36HZCn7Vcwn9R9PWId8lz5Kqw8MCPCF8ZTj4/DV8W8npaYF2gBwuSzwdEE+mX6Wh/eL79ZY0eKnD+QSj7qoB2ntt+OsgHZpfhFjjhtIx7Piom+CCEEEImFI5TmKtZGc/54LALIYQQQooKez4IIYSQUuBy2IUQQgghxYw9HEe5BQy7UGpLCCGEkPFGD6pSJ5xyzgchhBBCikpZBB+33XabmjlzpkokEuqwww5TL7zwQqmzRAghhBSG4xb+KVMm/ITTf/mXf1FXXnmluuOOO7zA4+c//7maP3++WrVqlWppaRnzeaJ9jopEnUANvpWWY25uxK8Lj/WBBj8Z7ttgDYI/xkceI9vT8FtR/zWr1w7Kc8STod4D8a3ZcH8D8BaJ9sn9c63S90OTqZdC/9hWaRZgz5ku0pHulEzDNZx4NNQfI8hDw+yXRg/Wh1vkOdsa5X1MlWknCmWdtvOG4LEt8ppmGrxGauKh92lk5TVyNdJbJLFB+p1ospNk/db3ynOaWWiDNdFQb5H6mLwxF3wijIAXlx2Xx5jgPRKFPOWq5SskC14ziBPwbKGXRLQfPFLgPmI98tmy4L6x3eOz6jP1CDgHXhPLJdof7vNjQHd4POB3RKYhGvo8J7bAfablfcS35kI9VdDnJ6js8V2IbSJXJa8Z7c2F1r8DRW2l/PMZsrWR0HaN+c5hm4Q8RwZkXbjwvEf7xvALenTRwLO7S3F13gqR2pZv8DHhez5+9rOfqQsvvFCdf/75av/99/eCkKqqKvXrX/+61FkjhBBCyO4WfGQyGbV8+XJ13HHHjXxnmqaXXrZsWUnzRgghhBSC67gFf8qVCT3ssnnzZmXbtmptbRXf6/Rbb70VeEw6nfY+w/T09OzyfBJCCCE75lDq7HjB0eF04rBw4UJVX18/8pk+Xc5JIIQQQkhpmdDDLs3NzcqyLLVx40bxvU63tbUFHnPttdeq7u7ukc+6deuKlFtCCCFk7LgVPOwyoYOPWCymDjnkEPXEE0+MfOc4jpeeN29e4DHxeFzV1dWJDyGEEDLhcJ3CP2XKhJ7zodEy23PPPVd96lOfUp/+9Kc9qW1/f7+nfhkL7kdSpFw2tV35oAEyWTdozW3AgHOYuLQ8LEWfQ5lkFmSQcD5vH1vKPXM5mS/MppsLlxz6pJUg0wpy+c2B7MzMSamti7pFyLONUj4bygG2D10DZIxwTsORaceW8l4nJ/V+jgFlDecPkqu5Jsj77FzoffnzjPfphN7T0D7h7c6Ec+TgmrjdMXdAaov3jfWXA6ltTr5Cctk8UtuAZwubnQHPig31Z8Kzla/d+/IUUMx4DrwmlouRG5/UNsjEEp8t3/Pse2fI+7By45TaBpQ9WgT4pLZwTcP3XgOpLVYF1k3AMXmltnnapIK6cPFvarjHQEYVTS6XEr87diU5lS3I4NQ7vkyZ8MHHl770JbVp0yb13e9+V7W3t6uDDjpIPfLII75JqNujt7fX+3/5Izft4pySkrGeZU8I2bno3x163uCu6tVva2tTz7Q/XPC59Hn0+coNwy1GeFdC9DDN+vXrvSh2xowZ3hwQDsUUhlYQ6Ym8LEuW40SBbZJlubPQvyt04DF16lTP2mFXkUqlPDuJQtGBh3b/LjcmfM9HoejGM23atBHJLeeB7DxYlizHiQbbJMtyZ7CrejxGk0gkyjJoqIgJp4QQQgjZ/WDwQQghhJCiUjHBh5bgfu973/P+JyzLiQDbJMtyIsJ2SYrBbj/hlBBCCCETi4rp+SCEEELIxIDBByGEEEKKCoMPQgghhBSVigk+brvtNjVz5kxPV33YYYepF154odRZmvCrAx966KGqtrZWtbS0qFNOOUWtWrXKZ5KzYMEC1dTUpGpqatTpp5/uWwSQSG6++WZlGIa6/PLLWY47wIcffqi+8pWveG0umUyqAw44QL300ksj2/UUNu2GPGXKFG/7cccdp95++202Q8C2bXX99derWbNmeeW09957qxtvvFFYirMsyS7FrQDuu+8+NxaLub/+9a/d119/3b3wwgvdhoYGd+PGjaXO2oRl/vz57uLFi92VK1e6K1ascD//+c+7M2bMcPv6+kb2+frXv+5Onz7dfeKJJ9yXXnrJPfzww90jjjiipPmeyLzwwgvuzJkz3QMPPNC97LLLRr5nOY6Nzs5Od88993TPO+889/nnn3ffffdd909/+pP7zjvvjOxz8803u/X19e7vf/9799VXX3W/8IUvuLNmzXIHBwd3QY2WLzfddJPb1NTkLlmyxF2zZo17//33uzU1Ne4vfvGLkX1YlmRXUhHBx6c//Wl3wYIFI2nbtt2pU6e6CxcuLGm+yomOjg79J5G7dOlSL93V1eVGo1HvpTXMm2++6e2zbNmyEuZ0YtLb2+vuu+++7mOPPeZ+9rOfHQk+WI5j5+qrr3Y/85nPbHe74zhuW1ub++Mf/3jkO12+8Xjc/e1vf1tQ/e1unHjiie5Xv/pV8d1pp53mnnXWWd7PLEuyq9nth120d/7y5cu97tfRlus6vWzZspLmrZzo7u72/m9sbPT+12WazWZFuc6ePdtbP4fl6kcPT5144omivFiO4+Ohhx7yVrf+7//9v3tDgQcffLD61a9+NbJ9zZo13uKTo8tY22TrYVa2SckRRxyhnnjiCfW3v/3NS7/66qvqmWeeUSeccALLkhSF3X5tl82bN3vjm7gKrk6/9dZbJctXuS3Op+coHHnkkWru3Lned/olrxc0amho8JWr3ka2cd9996mXX35Zvfjii75iYTmOnXfffVctWrRIXXnllerb3/62V57f+MY3vHZ47rnnjrS7oGedbVJyzTXXeOtd6T8YLMvy3pE33XSTOuuss0baJcuS7Ep2++CD7Jy/2leuXOn9ZUTGh17597LLLlOPPfZYRS8itbOCYN3z8Q//8A9eWvd86HZ5xx13eMEHGTv/+q//qu655x517733qo9//ONqxYoV3h8YeiVXliUpBrv9sEtzc7MX2aMKQ6fb2tpKlq9y4ZJLLlFLlixRf/7zn73VgYfRZaeHtLq6usT+LFeJHp7q6OhQn/zkJ1UkEvE+S5cuVbfeeqv3s/6rnOU4NrSCZf/99xffzZkzR61du3akTQ63QbbJcK666iqv9+OMM87wFENnn322uuKKKzyVG8uSFIPdPvjQXbKHHHKIN745+i8onZ43b15J8zaR0ZORdeDxwAMPqCeffNKT5I1Gl2k0GhXlqqW4+hcBy3Ubxx57rHrttde8vyyHP/qvd929Pfwzy3Fs6GE/lHvrOQt77rmn97NuozoAGd0m9dDC888/zzYJDAwMeHPfRqP/SNPvRpYlKQpuhUht9Yz3u+66y33jjTfciy66yJPatre3lzprE5aLL77Ykyw+9dRT7oYNG0Y+AwMDQiKq5bdPPvmkJ7WdN2+e9yHhjFa7sBzHJ1WORCKeTPTtt99277nnHreqqsr9zW9+I+Sh+tl+8MEH3b/+9a/uySefTKltAOeee667xx57jEhtf/e737nNzc3ut771LZYlKQoVEXxofvnLX3q/KLXfh5bePvfcc6XO0oRGx6VBH+39MYz2Tvif//N/upMmTfJ+CZx66qlegELGF3ywHMfOH/7wB3fu3LneHxOzZ892/+mf/kls1xLR66+/3m1tbfX2OfbYY91Vq1axSQI9PT1eG9TvxEQi4e61117ud77zHTedTrMsSVHgqraEEEIIKSq7/ZwPQgghhEwsGHwQQgghpKgw+CCEEEJIUWHwQQghhJCiwuCDEEIIIUWFwQchhBBCigqDD0IIIYQUFQYfhBBCCCkqDD4IKXOOPvpob0VSzcyZM9XPf/7zUmeJEEJCYfBByG7Eiy++qC666KJddv5nnnnGW+CtqalJJZNJNXv2bHXLLbfssusRQnZPIqXOACFk5zF58uRdWpzV1dXeascHHnig97MORr72ta95P+/KoIcQsnvBng9Cyoj+/n51zjnnqJqaGjVlyhT105/+VGzHYRfDMNQ//uM/qv/6X/+rqqqqUnPmzFHLli1T77zzjjdco4OGI444Qq1evXpM1z/44IPVl7/8ZfXxj3/cu9ZXvvIVNX/+fPWf//mfO/1eCSG7Lww+CCkjrrrqKrV06VL14IMPqkcffVQ99dRT6uWXXw495sYbb/QClhUrVnjDJGeeeabXW3Httdeql156Sa9s7fVm7AivvPKKevbZZ9VnP/vZHbwjQkglwmEXQsqEvr4+deedd6rf/OY36thjj/W+u/vuu9W0adNCjzv//PPVF7/4Re/nq6++Ws2bN09df/31Xo+F5rLLLvP2GQ/6mps2bVK5XE7dcMMN6n/8j/+xw/dFCKk8GHwQUibooZFMJqMOO+ywke8aGxvVfvvtF3qcnp8xTGtrq/f/AQccIL5LpVKqp6dH1dXVjSkvephFB0PPPfecuuaaa9Q+++zjDccQQshYYPBByG5ONBoVc0C2953jOGM+56xZs0aCmI0bN3q9Hww+CCFjhXM+CCkT9t57by9oeP7550e+27p1q/rb3/5W0nzpoCWdTpc0D4SQ8oI9H4SUCVrhcsEFF3iTTrXPRktLi/rOd76jTLN4f0PcdtttasaMGd7EVc3TTz+tfvKTn6hvfOMbRcsDIaT8YfBBSBnx4x//2JtrcdJJJ6na2lr1v/7X/1Ld3d1F7eXQKpk1a9aoSCTi9cb88Ic/9NQzhBAyVgxX6+wIIYQQQooE53wQQgghpKgw+CCEjKCdS/XckqDPPffcw5IihOwUOOxCCBnh/fffV9lsNrBEtB+InmdCCCGFwuCDEEIIIUWFwy6EEEIIKSoMPgghhBBSVBh8EEIIIaSoMPgghBBCSFFh8EEIIYSQosLggxBCCCFFhcEHIYQQQooKgw9CCCGEqGLy/wFUKCuQ/MePmwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Please note - this image was generated from only a few samples of training and does not represent the final model\n", + "da[5][0].plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25598079-8d9b-4893-a807-cfe1c50d35b8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/tutorial/LUCIE/LUCIE-Training.ipynb b/notebooks/tutorial/LUCIE/LUCIE-Training.ipynb new file mode 100644 index 00000000..8aaf970a --- /dev/null +++ b/notebooks/tutorial/LUCIE/LUCIE-Training.ipynb @@ -0,0 +1,194 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a3575c36-ed8d-4bae-ab90-aefe441949f9", + "metadata": {}, + "source": [ + "# Training the LUCIE model\n", + "\n", + "LUCIE is a climate model developed by Haiwen Guan, Troy Arcomano, Ashesh Chattopadhyay and Romit Maulik (2024). See their preprint at https://doi.org/10.48550/arXiv.2405.16297 and the archive of their training data, code and results here https://doi.org/10.5281/zenodo.15164648.\n", + "\n", + "The code in PyEarthTools was based on their code repository at https://github.com/ISCLPennState/LUCIE, which is made available under the MIT license (see the PyEarthTools NOTICE file for full information on this point)\n", + "\n", + "LUCIE is a model which of interest to climate researchers due to its long-term stability for rollouts for many decades. This model is licensed in a compatible fashion, so we are able to provide a bundled, customised version of LUCIE which can be used within the PyEarthTools framework, integrated with its data pipelines and configurable to work flexibly.\n", + "\n", + "We have only just begun the process of this integration, and so for now the model does not make extensive use of the PyEarthTools classes. This is expected to change fairly quickly, and as this happens, this notebook will be updated. However, in the interests of providing the bundled version to the community as soon as possible for those already seeking to work with the model, we present it in a \"work in progress\" fashion.\n", + "\n", + "You need to manually download the original published dataset from Zenodo, and update the paths in this notebook to point to them. The initial focus will be on reproducing the paper fairly closely using the same data and only slightly modified code (changes to support more devices and updates for compatibility), true enough to the original. Subsequently, we will develop the code further to be adaptable to new data sources.\n", + "\n", + "The intention is to:\n", + " - [done] Supply the source code to train and run the model in PyEarthTools\n", + " - [done] Validate that the model can train without obvious code-level errors\n", + " - Validate inference and reproduce the training results to ensure the trained model is valid\n", + " - Support library updates and other changes\n", + " - Support multiple ML backends beyond CUDA\n", + " - Support connection to multiple data sources through PET data accessors\n", + " - Move the normalisation into a PET pipeline so it can be easily modified and experimented with\n", + "\n", + "If you would like to know more, or get involved with this work, please [let us know on the issue tracker](https://github.com/ACCESS-Community-Hub/PyEarthTools/issues/211)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "e5068eca-cfcc-4dec-bf88-8b1fb870dc3b", + "metadata": {}, + "outputs": [], + "source": [ + "import lucie\n", + "import torch" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f69a338a-ff4e-465f-a664-cd76630baa52", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b81e2c08-bf62-49fc-9090-0595cbfd24ab", + "metadata": {}, + "outputs": [], + "source": [ + "device = torch.device(\"mps\" if torch.backends.mps.is_available() else \"cpu\")\n", + "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else device)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "4180dd8c-ff64-466b-b3bc-9771b2053a57", + "metadata": {}, + "outputs": [], + "source": [ + "regridded_path = Path.home() / 'dev/data/lucie' / 'era5_T30_regridded.npz'" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "7f5ca64a-87c8-4cae-a2e3-3a4788066a73", + "metadata": {}, + "outputs": [], + "source": [ + "regridded_data = lucie.train.load_data(regridded_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b53d4754-4303-4325-801a-afa626aac582", + "metadata": {}, + "outputs": [], + "source": [ + "preprocessed_path = Path.home() / 'dev/data/lucie' / 'era5_T30_preprocessed.npz'\n", + "preprocessed_data = np.load(preprocessed_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "22148013-b8d6-40c7-8c11-9f8545295b85", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Starting Training\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 0%| | 0/2 [00:00 0: - out = torch.sum( - ugrid[..., polar_opt:-polar_opt, :] * quad_weights[polar_opt:-polar_opt] * dlon * radius**2, dim=(-2, -1) - ) - else: - out = torch.sum(ugrid * quad_weights * dlon * radius**2, dim=(-2, -1)) - return out - - -def l2loss_sphere(prd, tar, relative=False, squared=True): - loss = integrate_grid((prd - tar) ** 2, dimensionless=True).sum(dim=-1) - if relative: - loss = loss / integrate_grid(tar**2, dimensionless=True).sum(dim=-1) - - if not squared: - loss = torch.sqrt(loss) - loss = loss.mean() - - return loss - - -def train_model( - model, - train_loader, - val_loader, - optimizer, - scheduler=None, - nepochs=20, - nfuture=0, - num_examples=256, - num_valid=8, - reg_rate=0, -): - - infer_bias = 1e80 - recall_count = 0 - for epoch in tqdm(range(nepochs)): - if epoch < 149: - if scheduler is not None: - scheduler.step() - else: - for param_group in optimizer.param_groups: - param_group["lr"] = 1e-6 - - optimizer.zero_grad() - - model.train() - batch_num = 0 - for inp, tar in train_loader: - batch_num += 1 - loss = 0 - - inp = inp.to(device) - tar = tar.to(device) - prd = model(inp) - - loss_delta = l2loss_sphere(prd[:, :5, :, :], tar[:, :5, :, :], relative=True) - loss_tp = torch.mean((prd[:, 5:, :, :] - tar[:, 5:, :, :]) ** 2) - loss = loss_delta + loss_tp / tar.shape[1] - - lat_index = np.r_[7:15, 32:40] - # lat_index = np.r_[0:48] - # quad_weight_reg = quad_weights.reshape(1,1,48,1)[:,:,lat_index,:] - out_fft = torch.mean(torch.abs(torch.fft.rfft(prd[:, :, lat_index, :], dim=3)), dim=2) - target_fft = torch.mean(torch.abs(torch.fft.rfft(tar[:, :, lat_index, :], dim=3)), dim=2) - loss_reg = 0.05 * torch.mean(torch.abs(out_fft - target_fft)) - - if epoch > 150: - loss = loss + loss_reg - - optimizer.zero_grad() - loss.backward() - optimizer.step() - - if epoch % 10 == 0: - rollout_steps = 2920 - rollout = torch.tensor( - inference( - model, - rollout_steps, - data_inp[0:1].to(device), - data_inp[:1460, -2:].to(device), - 1, - prog_means, - prog_stds, - diag_means, - diag_stds, - diff_stds, - ) - ).to(device) - rollout_clim = torch.mean(rollout[1460:], dim=0) - clim_bias = torch.mean(torch.abs(rollout_clim - true_clim)) - print("2 year rollout bias", clim_bias) - if epoch > 60: - if clim_bias <= infer_bias: - infer_bias = clim_bias - torch.save(model.state_dict(), "regular_training_checkpoint.pth") - recall_count = 0 - else: - state_pth = torch.load("regular_training_checkpoint.pth") - model.load_state_dict(state_pth) - recall_count += 1 - if recall_count > 3: - break - - -data = load_data("era5_T30_regridded.npz")[..., :6] -true_clim = torch.tensor(np.mean(data, axis=0)).to(device).permute(2, 0, 1) - -data = np.load("era5_T30_preprocessed.npz") # standardized data with mean and stds generated from dataset_generator.py -data_inp = torch.tensor(data["data_inp"], dtype=torch.float32) # input data -data_tar = torch.tensor(data["data_tar"], dtype=torch.float32) -raw_means = torch.tensor(data["raw_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) -raw_stds = torch.tensor(data["raw_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) -prog_means = raw_means[:, :5] -prog_stds = raw_stds[:, :5] -diag_means = torch.tensor(data["diag_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) -diag_stds = torch.tensor(data["diag_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) -diff_means = torch.tensor(data["diff_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) -diff_stds = torch.tensor(data["diff_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) - -ntrain = 16000 -nval = 100 - -train_set = TensorDataset(data_inp[:ntrain], data_tar[:ntrain]) -val_set = TensorDataset(data_inp[ntrain : ntrain + nval], data_tar[ntrain : ntrain + nval]) - -train_loader = DataLoader(train_set, batch_size=16, shuffle=True) -val_loader = DataLoader(val_set, batch_size=4, shuffle=False) - - -grid = "legendre-gauss" -nlat = 48 -nlon = 96 -hard_thresholding_fraction = 0.9 -lmax = ceil(nlat / 1) -mmax = lmax -modes_lat = int(nlat * hard_thresholding_fraction) -modes_lon = int(nlon // 2 * hard_thresholding_fraction) -modes_lat = modes_lon = min(modes_lat, modes_lon) -sht = RealSHT(nlat, nlon, lmax=modes_lat, mmax=modes_lon, grid=grid, csphase=False) -radius = 6.37122e6 -cost, quad_weights = legendre_gauss_weights(nlat, -1, 1) -quad_weights = (torch.as_tensor(quad_weights).reshape(-1, 1)).to(device) - -model = SphericalFourierNeuralOperatorNet( - params={}, - spectral_transform="sht", - filter_type="linear", - operator_type="dhconv", - img_shape=(48, 96), - num_layers=8, - in_chans=7, - out_chans=6, - scale_factor=1, - embed_dim=72, - activation_function="silu", - big_skip=True, - pos_embed="latlon", - use_mlp=True, - normalization_layer="instance_norm", - hard_thresholding_fraction=hard_thresholding_fraction, - mlp_ratio=2.0, -).to(device) - -optimizer = torch.optim.Adam(model.parameters(), lr=1e-4, weight_decay=0) -scheduler = CosineAnnealingLR(optimizer, T_max=150, eta_min=1e-5) -train_model(model, train_loader, val_loader, optimizer, scheduler=scheduler, nepochs=500) -torch.save(model.state_dict(), "model.pth") diff --git a/packages/bundled_models/lucie/README.md b/packages/bundled_models/lucie/README.md index 3eebfac5..dca21111 100644 --- a/packages/bundled_models/lucie/README.md +++ b/packages/bundled_models/lucie/README.md @@ -1,13 +1,17 @@ # LUCIE: Lightweight Uncoupled ClImate Emulator -Please note - this is a fork of https://github.com/ISCLPennState/LUCIE which has been adapted included in PyEarthTools for the purposes of maintenance, compatbility and to supply an integrated approach to using the LUCIE model within the PyEarthTools framework. +Please note - this is a adaptation of https://github.com/ISCLPennState/LUCIE which has been modified for inclusion in PyEarthTools for the purposes of maintenance, compatbility and to supply an integrated approach to using the LUCIE model within the PyEarthTools framework. + +This code was copied from the LUCIE repository from commit hash 19a1d6ebe844f49893f92e8b377ebdca8f6aa0e6 (Jul 9th, 2025). --- ## Paper & Data -- [arXiv Preprint: arxiv.org/abs/2405.16297](https://arxiv.org/abs/2405.16297) -- [Zenodo Archive: zenodo.org/records/15164648](https://zenodo.org/records/15164648) +These are the links for the original paper, code and data published by the LUCIE authors. The code was published to Zenodo under a Creative Commons license but the license in their github repository was MIT to allow improved code re-use. + +- [arXiv Preprint: https://doi.org/10.48550/arXiv.2405.16297](https://doi.org/10.48550/arXiv.2405.16297) +- [Zenodo Archive: [https://doi.org/10.5281/zenodo.15164648](https://doi.org/10.5281/zenodo.15164648) --- @@ -22,4 +26,4 @@ This repository prvides the following: 5. The data generator file that precprocesses the regridded ERA5 data. ## Note -Please refer to the zenodo link for the regridded ERA5 data. The link also includes the preprocessed data from the data generator file. +Please refer to the LUCIE zenodo link for the regridded ERA5 data. The link also includes the preprocessed data from the data generator file. diff --git a/packages/bundled_models/lucie/pyproject.toml b/packages/bundled_models/lucie/pyproject.toml new file mode 100644 index 00000000..34fbb687 --- /dev/null +++ b/packages/bundled_models/lucie/pyproject.toml @@ -0,0 +1,55 @@ +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + + +[project] +name = "pyearthtools-bundled-lucie" +version = "0.6.0" +description = "LUCIE Bundled Model" +readme = "README.md" +requires-python = ">=3.11, <3.14" +keywords = ["lucie"] +maintainers = [ + {name = "Tennessee Leeuwenburg", email = "tennessee.leeuwenburg@bom.gov.au"} +] +classifiers = [ + "License :: OSI Approved :: Apache Software License", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", +] +dependencies = [ + 'pyearthtools.training[lightning]>=0.5.1', + 'pyearthtools.zoo>=0.5.1', + 'pyearthtools.data>=0.5.1', + 'pyearthtools.pipeline>=0.5.1', + 'torch_optimizer', + 'timm', +] + + +[project.urls] +homepage = "https://pyearthtools.readthedocs.io/" +documentation = "https://pyearthtools.readthedocs.io/" +repository = "https://github.com/ACCESS-Community-Hub/PyEarthTools" + +[tool.isort] +profile = "black" + +[tool.black] +line-length = 120 + +[tool.mypy] +warn_return_any = true +warn_unused_configs = true + +[[tool.mypy.overrides]] +ignore_missing_imports = true + +[tool.hatch.version] +path = "src/lucie/__init__.py" + +[tool.hatch.build.targets.wheel] +packages = ["src/lucie/"] diff --git a/packages/bundled_models/lucie/src/lucie/__init__.py b/packages/bundled_models/lucie/src/lucie/__init__.py new file mode 100644 index 00000000..95a13f6c --- /dev/null +++ b/packages/bundled_models/lucie/src/lucie/__init__.py @@ -0,0 +1,2 @@ +from lucie import train +from lucie import torch_harmonics_local diff --git a/packages/bundled_models/lucie/dataset_generator.py b/packages/bundled_models/lucie/src/lucie/dataset_generator.py similarity index 100% rename from packages/bundled_models/lucie/dataset_generator.py rename to packages/bundled_models/lucie/src/lucie/dataset_generator.py diff --git a/packages/bundled_models/lucie/LUCIE_inference.py b/packages/bundled_models/lucie/src/lucie/inference.py similarity index 86% rename from packages/bundled_models/lucie/LUCIE_inference.py rename to packages/bundled_models/lucie/src/lucie/inference.py index a36c102f..93dcd756 100644 --- a/packages/bundled_models/lucie/LUCIE_inference.py +++ b/packages/bundled_models/lucie/src/lucie/inference.py @@ -28,13 +28,14 @@ # import torch_harmonics.distributed as thd # from torch_harmonics import * +# from torch._C import float32 import torch.fft from tqdm import tqdm import torch -from torch_harmonics_local import * +from lucie.torch_harmonics_local import * device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu") @@ -42,8 +43,9 @@ torch.cuda.set_device(0) -def inference( - model, steps, initial_frame, forcing, initial_forcing_idx, prog_means, prog_stds, diag_means, diag_stds, diff_stds +def infer(device, + model, steps, initial_frame, forcing, initial_forcing_idx, + prog_means, prog_stds, diag_means, diag_stds, diff_stds ): inf_data = [] model.eval() @@ -76,14 +78,17 @@ def inference( return inf_data -if __name__ == "__main__": +def load_data_and_predict( + device, + regridded_data, + preprocessed_data, # standardised data generated by dataset_generator.py + model_weights_pth='model.pth', + ): - data = load_data("era5_T30_regridded.npz")[..., :6] - true_clim = torch.tensor(np.mean(data, axis=0)).to(device).permute(2, 0, 1) + regridded_data = regridded_data[..., :6] + true_clim = torch.tensor(np.mean(regridded_data, axis=0)).to(device).permute(2, 0, 1) - data = np.load( - "era5_T30_preprocessed.npz" - ) # standardized data with mean and stds generated from dataset_generator.py + data = preprocessed_data # dictionary-like numpy array data_inp = torch.tensor(data["data_inp"], dtype=torch.float32) # input data data_tar = torch.tensor(data["data_tar"], dtype=torch.float32) raw_means = torch.tensor(data["raw_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) @@ -107,7 +112,8 @@ def inference( sht = RealSHT(nlat, nlon, lmax=modes_lat, mmax=modes_lon, grid=grid, csphase=False) radius = 6.37122e6 cost, quad_weights = legendre_gauss_weights(nlat, -1, 1) - quad_weights = (torch.as_tensor(quad_weights).reshape(-1, 1)).to(device) + quad_weights = (torch.as_tensor(quad_weights).reshape(-1, 1)).to(torch.float32).to(device) + # quad_weights = (torch.as_tensor(quad_weights).reshape(-1, 1)).to(device) model = SphericalFourierNeuralOperatorNet( params={}, @@ -129,7 +135,7 @@ def inference( mlp_ratio=2.0, ).to(device) - path = torch.load("regular_8x72_fftreg_baseline.pth") + path = torch.load(model_weights_pth) model.load_state_dict(path) forcing = data_inp[:1460, -2:] # repeating tisr and constant oro @@ -137,7 +143,8 @@ def inference( rollout_step = 14600 initial_frame_idx = 16000 + 100 forcing_initial_idx = (16000 + 100) % 1460 + 1 - rollout = inference( + rollout = infer( + device, model, rollout_step, data_inp[initial_frame_idx].unsqueeze(0).to(device), @@ -149,3 +156,5 @@ def inference( diag_stds, diff_stds, ) + + return rollout diff --git a/packages/bundled_models/lucie/torch_harmonics_local.py b/packages/bundled_models/lucie/src/lucie/torch_harmonics_local.py similarity index 99% rename from packages/bundled_models/lucie/torch_harmonics_local.py rename to packages/bundled_models/lucie/src/lucie/torch_harmonics_local.py index 18024be2..e7b4e08e 100644 --- a/packages/bundled_models/lucie/torch_harmonics_local.py +++ b/packages/bundled_models/lucie/src/lucie/torch_harmonics_local.py @@ -11,7 +11,7 @@ # from torch_harmonics import * import torch.nn.functional as F import torch.fft -from torch.cuda import amp +from torch import amp # was from torch.cuda import amp import math import logging @@ -1158,7 +1158,7 @@ def forward(self, x): # pragma: no cover x = x.float() B, C, H, W = x.shape - with amp.autocast(enabled=False): + with amp.autocast(str(device), enabled=False): x = self.forward_transform(x) if self.scale_residual: x = x.contiguous() @@ -1179,7 +1179,7 @@ def forward(self, x): # pragma: no cover # x = self._contract(x, self.weight, separable=self.separable, operator_type=self.operator_type) # x = x.contiguous() - with amp.autocast(enabled=False): + with amp.autocast(str(device), enabled=False): x = self.inverse_transform(x) if hasattr(self, "bias"): diff --git a/packages/bundled_models/lucie/src/lucie/train.py b/packages/bundled_models/lucie/src/lucie/train.py new file mode 100644 index 00000000..623ae39f --- /dev/null +++ b/packages/bundled_models/lucie/src/lucie/train.py @@ -0,0 +1,282 @@ +# MIT License + +# Copyright (c) 2025 ISCLPennState + +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: + +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. + +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +from math import ceil +import torch + +# import torch_harmonics as th +# import torch_harmonics.distributed as thd + +# from torch_harmonics import * +import torch.fft +from tqdm import tqdm + +import torch + +from torch.utils.data import TensorDataset, DataLoader + +from lucie.torch_harmonics_local import * + +from torch.optim.lr_scheduler import CosineAnnealingLR + +from lucie import inference + + +def integrate_grid(ugrid, nlon, quad_weights, dimensionless=False, polar_opt=0): + + dlon = 2 * torch.pi / nlon + radius = 1 if dimensionless else radius + if polar_opt > 0: + out = torch.sum( + ugrid[..., polar_opt:-polar_opt, :] * quad_weights[polar_opt:-polar_opt] * dlon * radius**2, dim=(-2, -1) + ) + else: + out = torch.sum(ugrid * quad_weights * dlon * radius**2, dim=(-2, -1)) + return out + + +def l2loss_sphere(prd, tar, nlon, quad_weights, relative=False, squared=True): + loss = integrate_grid((prd - tar) ** 2, nlon, quad_weights, dimensionless=True).sum(dim=-1) + if relative: + loss = loss / integrate_grid(tar**2, nlon, quad_weights, dimensionless=True).sum(dim=-1) + + if not squared: + loss = torch.sqrt(loss) + loss = loss.mean() + + return loss + + +def train_model( + device, + model, + train_loader, + val_loader, + optimizer, + data_inp=None, + prog_means=None, + prog_stds=None, + diag_means=None, + diag_stds=None, + diff_stds=None, + nlon=96, + scheduler=None, + nepochs=20, + debug_sample_limit=5, + quad_weights=None, + true_clim=None, + nfuture=0, + num_examples=256, + num_valid=8, + reg_rate=0, +): + """ + Train your own weights for the LUCIE model + """ + + infer_bias = 1e80 + recall_count = 0 + + print("Starting Training") + for epoch in tqdm(range(nepochs)): + + if epoch < 149: + if scheduler is not None: + scheduler.step() + else: + for param_group in optimizer.param_groups: + param_group["lr"] = 1e-6 + + optimizer.zero_grad() + + model.train() + + batch_num = 0 + + zz = 0 + + for inp, tar in train_loader: + batch_num += 1 + loss = 0 + + zz += 1 + if zz > debug_sample_limit: + break + + inp = inp.to(device) + tar = tar.to(device) + prd = model(inp) + + loss_delta = l2loss_sphere(prd[:, :5, :, :], tar[:, :5, :, :], nlon, quad_weights, relative=True) + loss_tp = torch.mean((prd[:, 5:, :, :] - tar[:, 5:, :, :]) ** 2) + loss = loss_delta + loss_tp / tar.shape[1] + + lat_index = np.r_[7:15, 32:40] + # lat_index = np.r_[0:48] + # quad_weight_reg = quad_weights.reshape(1,1,48,1)[:,:,lat_index,:] + out_fft = torch.mean(torch.abs(torch.fft.rfft(prd[:, :, lat_index, :], dim=3)), dim=2) + target_fft = torch.mean(torch.abs(torch.fft.rfft(tar[:, :, lat_index, :], dim=3)), dim=2) + loss_reg = 0.05 * torch.mean(torch.abs(out_fft - target_fft)) + + if epoch > 150: + loss = loss + loss_reg + + optimizer.zero_grad() + loss.backward() + optimizer.step() + + if epoch % 10 == 0: + rollout_steps = 2920 # Per paper + # rollout_steps = 50 # Testing + rollout = torch.tensor( + inference.infer( + device, + model, + rollout_steps, + data_inp[0:1].to(device), + data_inp[:1460, -2:].to(device), + 1, + prog_means, + prog_stds, + diag_means, + diag_stds, + diff_stds, + ) + ).to(device) + rollout_clim = torch.mean(rollout[1460:], dim=0) + clim_bias = torch.mean(torch.abs(rollout_clim - true_clim)) + print("2 year rollout bias", clim_bias) + if epoch > 60: + if clim_bias <= infer_bias: + infer_bias = clim_bias + torch.save(model.state_dict(), "regular_training_checkpoint.pth") + recall_count = 0 + else: + state_pth = torch.load("regular_training_checkpoint.pth") + model.load_state_dict(state_pth) + recall_count += 1 + if recall_count > 3: + break + + +def load_data_and_train( + device, + regridded_data, + preprocessed_data, + *, + debug_sample_limit: int | None = None, + n_epochs: int | None = 500, + ntrain: int | None = 16000, + nval: int | None = 100, +): + """ + + args: + unprocessed_data + reprocessed_data: dictionary or numpy collection containing 'diagn_means', 'diag_stds', 'diff_means' and 'diff_stds' + + """ + + regridded_data = regridded_data[..., :6] + true_clim = torch.tensor(np.mean(regridded_data, axis=0)).to(device).permute(2, 0, 1) + + data = preprocessed_data # dictionary-like numpy array + data_inp = torch.tensor(data["data_inp"], dtype=torch.float32) # input data + data_tar = torch.tensor(data["data_tar"], dtype=torch.float32) + raw_means = torch.tensor(data["raw_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + raw_stds = torch.tensor(data["raw_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + prog_means = raw_means[:, :5] + prog_stds = raw_stds[:, :5] + diag_means = torch.tensor(data["diag_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + diag_stds = torch.tensor(data["diag_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + diff_means = torch.tensor(data["diff_means"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + diff_stds = torch.tensor(data["diff_stds"], dtype=torch.float32).reshape(1, -1, 1, 1).to(device) + + train_set = TensorDataset(data_inp[:ntrain], data_tar[:ntrain]) + val_set = TensorDataset(data_inp[ntrain : ntrain + nval], data_tar[ntrain : ntrain + nval]) + + train_loader = DataLoader(train_set, batch_size=16, shuffle=True) + val_loader = DataLoader(val_set, batch_size=4, shuffle=False) + + grid = "legendre-gauss" + nlat = 48 + nlon = 96 + hard_thresholding_fraction = 0.9 + lmax = ceil(nlat / 1) + mmax = lmax + modes_lat = int(nlat * hard_thresholding_fraction) + modes_lon = int(nlon // 2 * hard_thresholding_fraction) + modes_lat = modes_lon = min(modes_lat, modes_lon) + # sht = RealSHT(nlat, nlon, lmax=modes_lat, mmax=modes_lon, grid=grid, csphase=False) + radius = 6.37122e6 + _cost, quad_weights = legendre_gauss_weights(nlat, -1, 1) + + # mps only supports float32, todo only do this if mps + # That said, most of the data seems to actually be in float32, so it is unclear if the weights + # benefit from being in float64 even on supported devices + # TODO: Experimentally verify the impact of using float32 here by default vs float64 performance + quad_weights = (torch.as_tensor(quad_weights).reshape(-1, 1)).to(torch.float32).to(device) + + model = SphericalFourierNeuralOperatorNet( + params={}, + spectral_transform="sht", + filter_type="linear", + operator_type="dhconv", + img_shape=(48, 96), + num_layers=8, + in_chans=7, + out_chans=6, + scale_factor=1, + embed_dim=72, + activation_function="silu", + big_skip=True, + pos_embed="latlon", + use_mlp=True, + normalization_layer="instance_norm", + hard_thresholding_fraction=hard_thresholding_fraction, + mlp_ratio=2.0, + ).to(device) + + optimizer = torch.optim.Adam(model.parameters(), lr=1e-4, weight_decay=0) + scheduler = CosineAnnealingLR(optimizer, T_max=150, eta_min=1e-5) + train_model( + device, + model, + train_loader, + val_loader, + optimizer, + prog_means=prog_means, + prog_stds=prog_stds, + diag_means=diag_means, + diag_stds=diag_stds, + diff_stds=diff_stds, + true_clim=true_clim, + data_inp=data_inp, + nlon=nlon, + quad_weights=quad_weights, + scheduler=scheduler, + nepochs=n_epochs, + debug_sample_limit=debug_sample_limit, + ) + + return model