Ablation-KAN-Nophy-LSTM_GRU.ipynb

{"cells":[{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"_4FezhF-dilH","outputId":"1f49f923-4cdd-46d1-fd17-5cfc71ce1d64","executionInfo":{"status":"ok","timestamp":1718410594837,"user_tz":240,"elapsed":75537,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}}},"outputs":[{"output_type":"stream","name":"stdout","text":["Requirement already satisfied: torchvision in /usr/local/lib/python3.10/dist-packages (0.18.0+cu121)\n","Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from torchvision) (1.25.2)\n","Requirement already satisfied: torch==2.3.0 in /usr/local/lib/python3.10/dist-packages (from torchvision) (2.3.0+cu121)\n","Requirement already satisfied: pillow!=8.3.*,>=5.3.0 in /usr/local/lib/python3.10/dist-packages (from torchvision) (9.4.0)\n","Requirement already satisfied: filelock in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (3.14.0)\n","Requirement already satisfied: typing-extensions>=4.8.0 in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (4.12.2)\n","Requirement already satisfied: sympy in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (1.12.1)\n","Requirement already satisfied: networkx in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (3.3)\n","Requirement already satisfied: jinja2 in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (3.1.4)\n","Requirement already satisfied: fsspec in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (2023.6.0)\n","Collecting nvidia-cuda-nvrtc-cu12==12.1.105 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cuda_nvrtc_cu12-12.1.105-py3-none-manylinux1_x86_64.whl (23.7 MB)\n","Collecting nvidia-cuda-runtime-cu12==12.1.105 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cuda_runtime_cu12-12.1.105-py3-none-manylinux1_x86_64.whl (823 kB)\n","Collecting nvidia-cuda-cupti-cu12==12.1.105 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cuda_cupti_cu12-12.1.105-py3-none-manylinux1_x86_64.whl (14.1 MB)\n","Collecting nvidia-cudnn-cu12==8.9.2.26 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cudnn_cu12-8.9.2.26-py3-none-manylinux1_x86_64.whl (731.7 MB)\n","Collecting nvidia-cublas-cu12==12.1.3.1 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cublas_cu12-12.1.3.1-py3-none-manylinux1_x86_64.whl (410.6 MB)\n","Collecting nvidia-cufft-cu12==11.0.2.54 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cufft_cu12-11.0.2.54-py3-none-manylinux1_x86_64.whl (121.6 MB)\n","Collecting nvidia-curand-cu12==10.3.2.106 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_curand_cu12-10.3.2.106-py3-none-manylinux1_x86_64.whl (56.5 MB)\n","Collecting nvidia-cusolver-cu12==11.4.5.107 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cusolver_cu12-11.4.5.107-py3-none-manylinux1_x86_64.whl (124.2 MB)\n","Collecting nvidia-cusparse-cu12==12.1.0.106 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_cusparse_cu12-12.1.0.106-py3-none-manylinux1_x86_64.whl (196.0 MB)\n","Collecting nvidia-nccl-cu12==2.20.5 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_nccl_cu12-2.20.5-py3-none-manylinux2014_x86_64.whl (176.2 MB)\n","Collecting nvidia-nvtx-cu12==12.1.105 (from torch==2.3.0->torchvision)\n","  Using cached nvidia_nvtx_cu12-12.1.105-py3-none-manylinux1_x86_64.whl (99 kB)\n","Requirement already satisfied: triton==2.3.0 in /usr/local/lib/python3.10/dist-packages (from torch==2.3.0->torchvision) (2.3.0)\n","Collecting nvidia-nvjitlink-cu12 (from nvidia-cusolver-cu12==11.4.5.107->torch==2.3.0->torchvision)\n","  Downloading nvidia_nvjitlink_cu12-12.5.40-py3-none-manylinux2014_x86_64.whl (21.3 MB)\n","\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m21.3/21.3 MB\u001b[0m \u001b[31m61.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n","\u001b[?25hRequirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from jinja2->torch==2.3.0->torchvision) (2.1.5)\n","Requirement already satisfied: mpmath<1.4.0,>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from sympy->torch==2.3.0->torchvision) (1.3.0)\n","Installing collected packages: nvidia-nvtx-cu12, nvidia-nvjitlink-cu12, nvidia-nccl-cu12, nvidia-curand-cu12, nvidia-cufft-cu12, nvidia-cuda-runtime-cu12, nvidia-cuda-nvrtc-cu12, nvidia-cuda-cupti-cu12, nvidia-cublas-cu12, nvidia-cusparse-cu12, nvidia-cudnn-cu12, nvidia-cusolver-cu12\n","Successfully installed nvidia-cublas-cu12-12.1.3.1 nvidia-cuda-cupti-cu12-12.1.105 nvidia-cuda-nvrtc-cu12-12.1.105 nvidia-cuda-runtime-cu12-12.1.105 nvidia-cudnn-cu12-8.9.2.26 nvidia-cufft-cu12-11.0.2.54 nvidia-curand-cu12-10.3.2.106 nvidia-cusolver-cu12-11.4.5.107 nvidia-cusparse-cu12-12.1.0.106 nvidia-nccl-cu12-2.20.5 nvidia-nvjitlink-cu12-12.5.40 nvidia-nvtx-cu12-12.1.105\n","Collecting pyDOE\n","  Downloading pyDOE-0.3.8.zip (22 kB)\n","  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n","Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from pyDOE) (1.25.2)\n","Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from pyDOE) (1.11.4)\n","Building wheels for collected packages: pyDOE\n","  Building wheel for pyDOE (setup.py) ... \u001b[?25l\u001b[?25hdone\n","  Created wheel for pyDOE: filename=pyDOE-0.3.8-py3-none-any.whl size=18168 sha256=aebfe7529f9065ef3c9a1c9f3c6e00fadfd438589a631a0d902de258c2c5c835\n","  Stored in directory: /root/.cache/pip/wheels/ce/b6/d7/c6b64746dba6433c593e471e0ac3acf4f36040456d1d160d17\n","Successfully built pyDOE\n","Installing collected packages: pyDOE\n","Successfully installed pyDOE-0.3.8\n"]}],"source":["!pip install torchvision\n","!pip install pyDOE"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"KqmQpWxud27e"},"outputs":[],"source":["import sys\n","sys.path.insert(0, '.')\n","import argparse\n","import os\n","import torch\n","from collections import OrderedDict\n","import math\n","import numpy as np\n","import matplotlib.pyplot as plt\n","import scipy.io\n","from scipy.interpolate import griddata\n","#from plotting import newfig, savefig\n","from mpl_toolkits.axes_grid1 import make_axes_locatable\n","import matplotlib.gridspec as gridspec\n","import torchvision.transforms as transforms\n","from torchvision.utils import save_image\n","from torch.utils.data import TensorDataset, DataLoader\n","from torchvision import datasets\n","from torch.autograd import Variable\n","import torch.nn as nn\n","import torch.nn.functional as F\n","import torch\n","import seaborn as sns\n","import pylab as py\n","import time\n","from pyDOE import lhs\n","import scipy\n","import warnings\n","sys.path.insert(0, '.')\n","warnings.filterwarnings('ignore')\n","np.random.seed(1234)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"-IP7Q93Ld6zJ","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410608431,"user_tz":240,"elapsed":6423,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"42b2ca3c-6217-4ccd-d0ef-b124cdccfee2"},"outputs":[{"output_type":"stream","name":"stdout","text":["Requirement already satisfied: pandas in /usr/local/lib/python3.10/dist-packages (2.0.3)\n","Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.10/dist-packages (from pandas) (2.8.2)\n","Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2023.4)\n","Requirement already satisfied: tzdata>=2022.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2024.1)\n","Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas) (1.25.2)\n","Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.2->pandas) (1.16.0)\n"]}],"source":["!pip install pandas"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"OAW9qOstd_HL","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410636897,"user_tz":240,"elapsed":28487,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"ee3356ed-1e08-4889-86e6-cbaf70f8c7f5"},"outputs":[{"output_type":"stream","name":"stdout","text":["Mounted at /content/drive\n"]}],"source":["from google.colab import drive\n","drive.mount('/content/drive')"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"xseWEYkAetXL"},"outputs":[],"source":["import pandas as pd"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"mjIeTqaWeGRm","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410641359,"user_tz":240,"elapsed":4472,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"3c422629-fb5a-40a3-86a3-9a03e1e5a27d"},"outputs":[{"output_type":"stream","name":"stdout","text":["Index(['nav_lat', 'nav_lon', 'time_centered', 'sithic'], dtype='object')\n","(2221440, 4)\n"]}],"source":["data1 = pd.read_csv(\"/content/drive/NewReq_Arc1-features.csv\")\n","df_features = pd.DataFrame(data1, columns = data1.columns)\n","#y = data1[\"sithic_value\"].values\n","X = df_features\n","X\n","X['time_centered'] = pd.to_datetime(df_features['time_centered']).astype(int)\n","print(df_features.columns)\n","print(X.shape)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"9PlElEJ0fFXl"},"outputs":[],"source":["df_features = pd.DataFrame(data1, columns = data1.columns)\n","y = data1[\"sithic\"].values"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"V7c7_2v_fPfz"},"outputs":[],"source":["#Set the 75/25% train/test split, set a random state to get train/test for future analysis\n","X = df_features\n","X\n","X['time_centered'] = pd.to_datetime(df_features['time_centered']).astype(int)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"_94Ir8XafRLn","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410641753,"user_tz":240,"elapsed":14,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"fb495db1-34d0-4a36-d75c-e969d91f0e79"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["(2221440, 4)"]},"metadata":{},"execution_count":9}],"source":["X.shape"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"h1ribEkXfT53","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410641753,"user_tz":240,"elapsed":12,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"4cfb21de-9ddb-46d7-9b0d-8555919f7b60"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["array([-58.420696, -53.43733 , 148.11154 , ...,  72.442535,  72.458015,\n","        72.47296 ])"]},"metadata":{},"execution_count":10}],"source":["B=X[\"nav_lon\"].values\n","B"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Zl9KxbZRf1le","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410641753,"user_tz":240,"elapsed":8,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"095931af-12b7-4dfb-f5fb-169fff69f41a"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["nav_lat          float64\n","nav_lon          float64\n","time_centered      int64\n","sithic           float64\n","dtype: object"]},"metadata":{},"execution_count":11}],"source":["x= data1[\"nav_lon\"].values\n","t=data1[\"time_centered\"].values\n","t\n","df_features.dtypes"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"vuwNzZQXfXfF","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410642531,"user_tz":240,"elapsed":784,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"7ffad63e-debc-47f5-b798-d2283c6d7408"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["Index([1975, 1975, 1975, 1975, 1975, 1975, 1975, 1975, 1975, 1975,\n","       ...\n","       2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013],\n","      dtype='int32', length=2165904)"]},"metadata":{},"execution_count":12}],"source":["import pandas as pd\n","####Test differentiation#########\n","\n","# Define the array of timestamps\n","time_strings = t\n","\n","# Convert the strings to datetime objects\n","time_objects = pd.to_datetime(time_strings)\n","\n","# Extract the dates\n","dates = time_objects.year\n","\n","# Create a DataFrame with the year information\n","year_info = pd.DataFrame({'year': dates})\n","\n","\n","# Concatenate X DataFrame and year_info DataFrame horizontally\n","X_with_year = pd.concat([X, year_info], axis=1)\n","\n","# Create a Pandas DataFrame\n","#df_date = pd.DataFrame(data={'dates': dates})\n","\n","# Print the DataFrame\n","\n","X_with_year\n","#nan_indices = year_info['year'].isnull().sum()\n","#nan_indices\n","# Filter X for data corresponding to the year 2007\n","X_2007 = X_with_year[X_with_year['year'] == 2014]\n","data1['time_derivative'] =year_info.diff()\n","\n","data1['time_derivative']\n","# Display the updated DataFrame with extracted features for 2007\n","#print(X_2007)\n","time_objects_excluding_2007 = dates[dates != 2014]\n","time_objects_excluding_2007"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"9faRKyq5gkpa","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410642531,"user_tz":240,"elapsed":7,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"695fc2c1-d103-4c58-e29e-5462392f57d7"},"outputs":[{"output_type":"stream","name":"stdout","text":["Minimum value of 'time_derivative': 0.0\n","Maximum value of 'time_derivative': 1.0\n"]}],"source":["min_derivative = data1['time_derivative'].min()\n","max_derivative = data1['time_derivative'].max()\n","\n","# Print the minimum and maximum values\n","print(\"Minimum value of 'time_derivative':\", min_derivative)\n","print(\"Maximum value of 'time_derivative':\", max_derivative)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"esjQ3QkWhc4x"},"outputs":[],"source":["from sklearn.model_selection import train_test_split\n","from sklearn.linear_model import LinearRegression"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"A_bkq9lrhOzY"},"outputs":[],"source":["# Filter X_with_year to exclude data from the year 2007\n","X_except_2007 = X_with_year[X_with_year['year'] != 2014]\n","y_except_2007 = y[X_with_year['year'] != 2014]\n","\n","# Split the data excluding 2007 into training and testing sets\n","X_train_except_2007, X_test_except_2007, y_train_except_2007, y_test_except_2007 = train_test_split(\n","    X_except_2007, y_except_2007, test_size=0.2, random_state=0\n",")\n","\n","# Filter X_with_year for data corresponding to the year 2007\n","X_2007 = X_with_year[X_with_year['year'] == 2014]\n","y_2007 = y[X_with_year['year'] == 2014]\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"SORcUnZjwMjl","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410643662,"user_tz":240,"elapsed":22,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"349781e4-fcab-48ac-cbee-a7828c0c00d0"},"outputs":[{"output_type":"stream","name":"stdout","text":["[[ 4.9599170e+01 -5.8420696e+01  1.3898736e+18  4.5051612e-02\n","   2.0140000e+03]\n"," [ 5.0527065e+01 -5.3437330e+01  1.3898736e+18  3.5513207e-01\n","   2.0140000e+03]\n"," [ 5.1792973e+01  1.4811154e+02  1.3898736e+18  7.2821360e-02\n","   2.0140000e+03]\n"," ...\n"," [ 6.9997670e+01  7.2442535e+01  1.4187312e+18  1.0973377e+00\n","   2.0140000e+03]\n"," [ 6.9552220e+01  7.2458015e+01  1.4187312e+18  1.2507676e+00\n","   2.0140000e+03]\n"," [ 6.9107480e+01  7.2472960e+01  1.4187312e+18  1.4002302e+00\n","   2.0140000e+03]]\n","float64\n"]}],"source":["print(X_2007.values)\n","print(y_2007.dtype)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"tjj9wjO_y6IH"},"outputs":[],"source":["X_2007=X_2007.values[:,[0,1,4]]\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"GEsUdWMojGmO","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410643662,"user_tz":240,"elapsed":14,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"d1ec06d9-2d34-431a-89e5-2a21a5c06292"},"outputs":[{"output_type":"stream","name":"stdout","text":["(1732723, 5)\n","(433181, 5)\n","(1732723,)\n","(433181,)\n"]}],"source":["print(X_train_except_2007.shape)\n","print(X_test_except_2007.shape)\n","print(y_train_except_2007.shape)\n","print(y_test_except_2007.shape)\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"0uWkebokXZmf","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410643663,"user_tz":240,"elapsed":12,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"f1129806-bb2b-4abc-fd9d-7ef9a49865a6"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["(1732723, 3)"]},"metadata":{},"execution_count":19}],"source":["X_train=X_train_except_2007.values[:,[0,1,4]]\n","X_test=X_test_except_2007.values[:,[0,1,4]]\n","X_train.shape"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"e9KYcqi2oD5v"},"outputs":[],"source":["X_star=X_test\n","u_star = y_test_except_2007\n","X_u_train = X_train\n","u_train=y_train_except_2007"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"DuH7sDUYqAdg"},"outputs":[],"source":["# Choose a random sample size (\n","sample_size = 10000\n","\n","# Select random indices from the entire data range (0 to length-1)\n","random_indices = np.random.choice(X_star.shape[0], size=sample_size, replace=False)\n","\n","# Select random data using the indices\n","X_star_random = X_star[random_indices]\n","u_star_random = u_star[random_indices]\n","\n","# Similarly, for training data\n","random_indices_train = np.random.choice(X_u_train.shape[0], size=sample_size, replace=False)\n","X_u_train_random = X_u_train[random_indices_train]\n","u_train_random = u_train[random_indices_train]\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"7gBaGj5c2sQl"},"outputs":[],"source":["import numpy as np\n","\n","# Assuming X_star might have overlapping indices with previously selected points\n","\n","# Number of elements not yet selected (might be less than 5000)\n","n_remaining = X_star.shape[0] - len(random_indices)\n","\n","# Select a random sample size (up to available elements)\n","sample_size = min(1000, n_remaining)\n","\n","# Randomly choose indices from the remaining set (if any)\n","if n_remaining > 0:\n","  # Ensure no duplicates between random_indices and chosen indices\n","  next_5000_indices = np.random.choice(np.setdiff1d(np.arange(X_star.shape[0]), random_indices), size=sample_size, replace=False)\n","  X_star_next_5000 = X_star[next_5000_indices]\n","else:\n","  # Handle the case where no elements are remaining (all might have been selected previously)\n","  print(\"No elements in X_star haven't been selected yet\")\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"jkQUn8Zp2zLg","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410644188,"user_tz":240,"elapsed":51,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"bea596e9-29d2-4fee-e98d-825b4a35551e"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["array([4.6202254 , 3.4149413 , 0.15174903, ..., 2.9091349 , 3.061613  ,\n","       1.5146224 ])"]},"metadata":{},"execution_count":23}],"source":["X_f_train =X_star_next_5000\n","X_u_train_random\n","u_train_random"]},{"cell_type":"code","source":["print(X_f_train.shape)\n","print(X_u_train_random.shape)\n","print(u_train_random.shape)\n","print(u_star_random.shape)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"y45hq3bXU735","executionInfo":{"status":"ok","timestamp":1718410644188,"user_tz":240,"elapsed":46,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"204723e4-338e-46b8-d2b7-723326f27ff8"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["(1000, 3)\n","(10000, 3)\n","(10000,)\n","(10000,)\n"]}]},{"cell_type":"code","execution_count":null,"metadata":{"id":"bo7YvizHMShe","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1718410644189,"user_tz":240,"elapsed":45,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"75f640d2-1c6d-4b2c-981c-42cd658f2536"},"outputs":[{"output_type":"stream","name":"stdout","text":["10000\n"]}],"source":["print(len(X_u_train_random))"]},{"cell_type":"code","source":["import torch\n","import torch.nn as nn\n","import torch.optim as optim\n","class Spline(nn.Module):\n","    def __init__(self, num_features, spline_order=4):\n","        super(Spline, self).__init__()\n","        self.num_features = num_features\n","        self.spline_order = spline_order\n","        self.coefficients = nn.Parameter(torch.randn(num_features, spline_order))\n","\n","    def forward(self, x):\n","        # For simplicity, using a linear combination of input and spline coefficients\n","        spline_basis = torch.stack([x**i for i in range(self.spline_order)], dim=-1)\n","        out = torch.sum(spline_basis * self.coefficients, dim=-1)\n","        return out\n","\n","class KANBlock(nn.Module):\n","    def __init__(self, in_features, out_features, spline_order=4):\n","        super(KANBlock, self).__init__()\n","        self.linear = nn.Linear(in_features, out_features)\n","        self.spline = Spline(out_features, spline_order)\n","\n","    def forward(self, x):\n","        out = self.linear(x)\n","        out = self.spline(out)\n","        return out\n","\n","class Net(nn.Module):\n","    def __init__(self, input_dim, hidden_dim, output_dim, num_kan_layers=3):\n","        super(Net, self).__init__()\n","\n","        self.rnn = nn.RNN(input_dim, hidden_dim, batch_first=True)\n","                # Subsequent GRU layers\n","        self.lstm_layers = nn.ModuleList([\n","            nn.LSTM(hidden_dim, hidden_dim, batch_first=True),\n","            nn.LSTM(hidden_dim, hidden_dim, batch_first=True),\n","            nn.LSTM(hidden_dim, hidden_dim, batch_first=True),\n","        ])\n","\n","\n","        # KAN blocks\n","        self.kan_layers = nn.ModuleList([KANBlock(hidden_dim, hidden_dim) for _ in range(num_kan_layers)])\n","\n","        # Apply dropout between layers (optional)\n","        self.dropout = nn.Dropout(0.2)  # Adjust dropout probability as needed\n","\n","        # Activation function (choose a suitable activation function)\n","        self.activation = nn.LeakyReLU()  # Example: ReLU activation\n","\n","        # Final linear layer\n","        self.fc = nn.Linear(hidden_dim, output_dim)\n","\n","    def forward(self, x):\n","        # Reshape to 3D if necessary\n","        if len(x.shape) == 2:\n","            x = x.view(-1, 1, x.shape[1])\n","\n","        # Pass input through RNN layer\n","        out, _ = self.rnn(x)\n","\n","        h = None  # Initialize hidden state for GRU layers\n","        # Pass output from RNN through KAN layers with dropout and activation\n","        for layer in self.kan_layers:\n","            for lstm in self.lstm_layers:\n","                out, h = lstm(out, h)\n","            out = self.dropout(self.activation(out))\n","\n","        # Use output from last timestep and apply linear layer\n","        out = self.fc(out[:, -1])\n","        return out"],"metadata":{"id":"QOt8khju66oU"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["import numpy as np\n","import torch\n","import torch.nn as nn\n","import torch.nn.functional as F\n","import scipy.io\n","from scipy.interpolate import griddata\n","from torch.utils.data import TensorDataset, DataLoader\n","import time\n","import torch.nn as nn\n","\n","def set_seed(seed_value):\n","    \"\"\"Set seed for reproducibility.\"\"\"\n","    np.random.seed(seed_value)\n","    torch.manual_seed(seed_value)\n","    torch.cuda.manual_seed_all(seed_value)\n","\n","# Set seed value for reproducibility\n","seed_value = 1234\n","set_seed(seed_value)\n","\n","\n","class BPINN(nn.Module):\n","    def __init__(self, x_u, y_u, x_f, X_star, u_star, net, nepochs, Omega_mse, noise = 0.1):\n","        super(BPINN, self).__init__()\n","        self.device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n","        # Normalize data\n","        self.Xmean, self.Xstd = x_f.mean(0), x_f.std(0)\n","        self.x_f = (x_f - self.Xmean) / self.Xstd\n","        self.x_u = (x_u - self.Xmean) / self.Xstd\n","        self.X_star_norm = (X_star - self.Xmean) / self.Xstd\n","        self.u_star = u_star\n","\n","        #Jacobian of the PDE because of normalization\n","        self.Jacobian_X = 1 / self.Xstd[0]\n","        self.Jacobian_Y = 1 / self.Xstd[1]\n","        self.Jacobian_T = 1 / self.Xstd[2]\n","\n","        self.y_u = y_u + noise * np.std(y_u)*np.random.randn(y_u.shape[0])\n","\n","        self.net = net\n","\n","\n","        self.net_optim = torch.optim.Adam(self.net.parameters(), lr=0.001, betas = (0.9, 0.999))\n","\n","\n","        #self.device = device\n","\n","        # numpy to tensor\n","        self.train_x_u = torch.tensor(self.x_u, requires_grad=True).float()\n","        self.train_y_u = torch.tensor(self.y_u, requires_grad=True).float()\n","        self.train_x_f = torch.tensor(self.x_f, requires_grad=True).float()\n","        self.X_star_norm = torch.tensor(self.X_star_norm, requires_grad=True).float()\n","        self.u_star = torch.tensor(self.u_star, requires_grad=True).float()\n","\n","        self.nepochs = nepochs\n","        self.Omega_mse = Omega_mse\n","\n","        self.batch_size = 100\n","        num_workers = 4\n","        shuffle = True\n","        self.train_loader = DataLoader(\n","            list(zip(self.train_x_u,self.train_y_u)), batch_size=self.batch_size, shuffle=shuffle\n","        )\n","\n","        self.val_loader = DataLoader(\n","            list(zip(self.X_star_norm,self.u_star)), batch_size=self.batch_size, shuffle=shuffle\n","        )\n","\n","\n","\n","    def train(self, patience=10):\n","\n","        MSE_loss = np.zeros(self.nepochs)\n","        VAL_loss = np.zeros(self.nepochs)\n","\n","        epoch_times = []\n","        start_time = time.perf_counter()\n","        min_val_loss = 9999\n","        no_improvement_count = 0\n","        # Define the number of batches to accumulate gradients over\n","        accum_iter = 4  # Adjust this value based on your memory limitations\n","        for epoch in range(self.nepochs):\n","            epoch_loss = 0\n","            epoch_val_loss = 0\n","\n","            for i, (x, y) in enumerate(self.train_loader):\n","                self.net_optim.zero_grad()\n","                y_pred = self.net.forward(x)\n","\n","                mse_loss = torch.nn.functional.mse_loss(y, y_pred)\n","\n","                loss =  mse_loss\n","\n","                loss.backward(retain_graph=True)\n","                # Update weights after 'accum_iter' or last batch\n","                if (i + 1) % accum_iter == 0 or (i + 1) == len(self.train_loader):\n","                  self.net_optim.step()\n","                  self.net_optim.zero_grad()  # Reset gradients for next iteration\n","                self.net_optim.step()\n","\n","\n","                MSE_loss[epoch] += mse_loss.detach().cpu().numpy()\n","\n","\n","            MSE_loss[epoch] /= len(self.train_loader)\n","\n","            # Calculate validation loss\n","            #with torch.no_grad():\n","            self.net.eval()\n","            for x_val, y_val in self.val_loader:\n","                y_pred_val = self.net.forward(x)\n","                val_mse_loss = torch.nn.functional.mse_loss(y_val, y_pred_val)\n","                epoch_val_loss += val_mse_loss.detach().cpu().numpy()\n","            VAL_loss[epoch] = epoch_val_loss / len(self.train_loader)\n","\n","            if min_val_loss > val_mse_loss**0.2:#0.5\n","                min_val_loss = val_mse_loss**0.2\n","                torch.save(burgers.state_dict(), \"/content/drive/MyDrive/No-physics-KDD_GRU-RNN-seq30.pt\")\n","                no_improvement_count = 0\n","            else:\n","                no_improvement_count += 1\n","\n","            if no_improvement_count == patience:\n","                break\n","\n","            if epoch % 10 == 0:\n","                print(\n","                    \"[Epoch %d/%d]  [mse_loss: %f] [val_loss: %f]\"\n","                    % (epoch, self.nepochs,  MSE_loss[epoch], VAL_loss[epoch])\n","                )\n","\n","        total_training_time = time.perf_counter() - start_time\n","        print(\"Total Training Time: {} seconds\".format(total_training_time))\n","        plt.figure(figsize=(6, 6))\n","        plt.plot(MSE_loss, label='MSE loss')\n","        plt.plot(VAL_loss, label='VAL loss')\n","        plt.ylim(0, 9)\n","        plt.xlabel('Epoch')\n","        plt.ylabel('Loss')\n","        plt.legend()\n","        plt.grid(False)\n","        plt.show()"],"metadata":{"id":"a78azAn8t4MR"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"vaK2b4Sk4UIr"},"outputs":[],"source":["num_epochs = 150\n","Omega_mse = 1\n","noise = 0.01"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"iUkW-i7s30NH"},"outputs":[],"source":["net=Net(input_dim=3, hidden_dim=64, output_dim=1,  num_kan_layers=3)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"_mcmx4-p34z8"},"outputs":[],"source":["burgers = Burgers_APINN(X_u_train_random, u_train_random, X_f_train, X_star_random, u_star_random, net, num_epochs, Omega_mse, noise)"]},{"cell_type":"code","source":["burgers.train()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":634},"id":"R0VPwO5zpmi3","executionInfo":{"status":"ok","timestamp":1718415034760,"user_tz":240,"elapsed":369370,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"298d0f0c-6b0b-4baf-fb2a-751a5d16a84d"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["[Epoch 0/150]  [mse_loss: 2.427872] [val_loss: 2.159724]\n","[Epoch 10/150]  [mse_loss: 2.160010] [val_loss: 2.156984]\n","[Epoch 20/150]  [mse_loss: 2.160406] [val_loss: 2.154501]\n","Total Training Time: 369.64177335100067 seconds\n"]},{"output_type":"display_data","data":{"text/plain":["<Figure size 600x600 with 1 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAgsAAAISCAYAAABChoydAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA20UlEQVR4nO3de3wU1d3H8e9kN1kCJBvuSSThVjSEmyiXIla0oIAUUSgqjYJCa9EoIi1V9EFRawP2qVJrn3gpBC8glRaUWhUBlZsiN0FURFQKAQ1YMdlw2yS78/yRZM0mYbIJIZNkP+/Xa15kZ052fido5ss5Z2cM0zRNAQAAnEaE3QUAAID6jbAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAs2RoW8vPzNW3aNHXo0EHR0dG66KKLtGXLFjtLAgAA5dgaFn75y19q1apVeuGFF7Rr1y5dccUVGjp0qA4dOmRnWQAAoAzDrgdJnTx5UjExMXr11Vc1cuTIwP4LL7xQI0aM0O9//3s7ygIAAOU47TpxUVGRfD6fmjRpErQ/OjpaGzZsqPR7vF6vvF5v4LXf79fRo0fVqlUrGYZxVusFAKAxMU1T+fn5SkxMVEREFRMNpo0GDhxoDh482Dx06JBZVFRkvvDCC2ZERIR57rnnVtr+gQceMCWxsbGxsbGx1dKWnZ1d5fXatmkISfryyy81adIkrVu3Tg6HQxdccIHOPfdcbdu2Tbt3767QvvzIQl5enpKTk5Wdna3Y2Ni6LB0AgAbN4/EoKSlJubm5crvdlm1tm4aQpC5dumjt2rU6fvy4PB6PEhISdN1116lz586Vtne5XHK5XBX2x8bGEhYAAKiBUKbx68V9Fpo1a6aEhAR9//33WrlypUaPHm13SQAAoIStIwsrV66UaZo677zz9MUXX2jGjBlKSUnRzTffbGdZAACgDFtHFvLy8pSenq6UlBRNmDBBF198sVauXKnIyEg7ywIAAGXYusDxTHk8HrndbuXl5bFmAQDqgM/nU2Fhod1lIAQOh0NOp/O0axKqcw21dRoCANBwHDt2TAcPHlQD/jdm2GnatKkSEhIUFRV1Ru9DWAAAVMnn8+ngwYNq2rSp2rRpw43w6jnTNFVQUKBvv/1W+/btU9euXau+8ZIFwgIAoEqFhYUyTVNt2rRRdHS03eUgBNHR0YqMjNT+/ftVUFBQ4Y7J1VEvPjoJAGgYGFFoWM5kNCHofWrlXQAAQKNFWAAAAJYICwAAnGUdO3bUvHnz7C6jxggLAIBG66abbpJhGJoyZUqFY+np6TIMQzfddFNg37fffqtbb71VycnJcrlcio+P17Bhw7Rx48ZAm44dO8owjArbnDlz6qJLtuDTEACARi0pKUlLlizR448/Hvgkx6lTp7R48WIlJycHtR07dqwKCgr03HPPqXPnzjp8+LDWrFmj7777LqjdQw89pF/96ldB+2JiYs5uR2zEyAIAoNpM09SJgiJbtureFOqCCy5QUlKSli1bFti3bNkyJScnq0+fPoF9ubm5Wr9+vebOnavLLrtMHTp0UP/+/TVz5kxdddVVQe8ZExOj+Pj4oK1Zs2Yh13TgwAGNHj1azZs3V2xsrK699lodPnw4cHznzp267LLLFBMTo9jYWF144YXaunWrJGn//v0aNWqUWrRooWbNmql79+56/fXXq/UzqS5GFgAA1Xay0KfU+1facu5PHxqmplHVu3xNmjRJWVlZSktLkyQtWLBAN998s959991Am+bNm6t58+Z65ZVX9OMf/1gul6s2yw7w+/2BoLB27VoVFRUpPT1d1113XaCetLQ09enTR5mZmXI4HNqxY0fguUnp6ekqKCjQunXr1KxZM3366adq3rz5Wam1FCMLAIBG74YbbtCGDRu0f/9+7d+/Xxs3btQNN9wQ1MbpdGrhwoV67rnnFBcXp0GDBunee+/VRx99VOH97r777kC4KN3Wr18fUi1r1qzRrl27tHjxYl144YUaMGCAnn/+ea1du1ZbtmyRVDzyMHToUKWkpKhr164aN26cevfuHTg2aNAg9ezZU507d9bPfvYzXXLJJWf4E7LGyAIAoNqiIx369KFhtp27utq0aaORI0dq4cKFMk1TI0eOVOvWrSu0Gzt2rEaOHKn169dr06ZNeuONN/Too4/qb3/7W9BCyBkzZgS9lqRzzjknpFp2796tpKQkJSUlBfalpqYqLi5Ou3fvVr9+/TR9+nT98pe/1AsvvKChQ4dq3Lhx6tKliyRp6tSpuvXWW/XWW29p6NChGjt2rHr16lXtn0l1MLIAAKg2wzDUNMppy1bTu0hOmjQpMHIwadKk07Zr0qSJLr/8cs2aNUvvvfeebrrpJj3wwANBbVq3bq0f/ehHQVtt3gZ79uzZ+uSTTzRy5Ei9/fbbSk1N1fLlyyVJv/zlL/XVV1/pxhtv1K5du9S3b1/95S9/qbVzV4awAAAIC8OHD1dBQYEKCws1bFjooyKpqak6fvx4rdXRrVs3ZWdnKzs7O7Dv008/VW5urlJTUwP7zj33XN1111166623NGbMGGVlZQWOJSUlacqUKVq2bJl+85vf6Nlnn621+irDNAQAICw4HA7t3r078HV53333ncaNG6dJkyapV69eiomJ0datW/Xoo49q9OjRQW3z8/OVk5MTtK9p06aKjY2tso6hQ4eqZ8+eSktL07x581RUVKTbbrtNgwcPVt++fXXy5EnNmDFDP//5z9WpUycdPHhQW7Zs0dixYyVJ06ZN04gRI3Tuuefq+++/1zvvvKNu3brV9McSEsICACBsWF3MmzdvrgEDBujxxx/Xl19+qcLCQiUlJelXv/qV7r333qC2999/v+6///6gfb/+9a/11FNPVVmDYRh69dVXdccdd+iSSy5RRESEhg8fHphKcDgc+u677zRhwgQdPnxYrVu31pgxY/Tggw9KKn5ceHp6ug4ePKjY2FgNHz5cjz/+eHV/FNVimNX9wGo94vF45Ha7lZeXF1KaAwDUzKlTp7Rv3z516tTpjB51jLpl9fdWnWsoaxYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALVs9uzZOv/88+0uo9YQFgAAjdKoUaM0fPjwSo+tX79ehmHoo48+Cuz79a9/LYfDoaVLl1Zo39gu/tVFWAAANEqTJ0/WqlWrdPDgwQrHsrKy1LdvX/Xq1UuSdOLECS1ZskS/+93vtGDBgroutd4jLAAAqs80pYLj9mwhPv/wZz/7mdq0aaOFCxcG7T927JiWLl2qyZMnB/YtXbpUqampuueee7Ru3TplZ2fX5k9Lfr9fDz30kNq3by+Xy6Xzzz9fb775ZuB4QUGBbr/9diUkJKhJkybq0KGDMjIyJEmmaWr27NlKTk6Wy+VSYmKipk6dWqv1VYVHVAMAqq/whPSHRHvOfe/XUlSzKps5nU5NmDBBCxcu1H333SfDMCQVBwOfz6fx48cH2s6fP1833HCD3G63RowYoYULF2rWrFm1VvKf//xn/elPf9LTTz+tPn36aMGCBbrqqqv0ySefqGvXrnriiSe0YsUKvfzyy0pOTlZ2dnYgsPzzn//U448/riVLlqh79+7KycnRzp07a622UDCyAABotCZNmqQvv/xSa9euDezLysrS2LFj5Xa7JUl79+7Vpk2bdN1110mSbrjhBmVlZckMcQQjFP/7v/+ru+++W9dff73OO+88zZ07V+eff77mzZsnSTpw4IC6du2qiy++WB06dNDFF18cCDMHDhxQfHy8hg4dquTkZPXv31+/+tWvaq22UDCyAACovsimxf/Ct+vcIUpJSdFFF12kBQsW6NJLL9UXX3yh9evX66GHHgq0WbBggYYNG6bWrVtLkq688kpNnjxZb7/9toYMGXLG5Xo8Hn399dcaNGhQ0P5BgwYFRghuuukmXX755TrvvPM0fPhw/exnP9MVV1whSRo3bpzmzZunzp07a/jw4bryyis1atQoOZ11dwlnZAEAUH2GUTwVYMdWMp0QqsmTJ+uf//yn8vPzlZWVpS5dumjw4MGSJJ/Pp+eee07//ve/5XQ65XQ61bRpUx09erROFzpecMEF2rdvnx5++GGdPHlS1157rX7+859LkpKSkrRnzx793//9n6Kjo3XbbbfpkksuUWFhYZ3VR1gAADRq1157rSIiIrR48WI9//zzmjRpUmD9wuuvv678/Hx9+OGH2rFjR2B76aWXtGzZMuXm5p7x+WNjY5WYmKiNGzcG7d+4caNSU1OD2l133XV69tln9fe//13//Oc/dfToUUlSdHS0Ro0apSeeeELvvvuu3n//fe3ateuMawsV0xAAgEatefPmuu666zRz5kx5PB7ddNNNgWPz58/XyJEj1bt376DvSU1N1V133aVFixYpPT1dknTy5Ent2LEjqF1MTIy6dOlSZQ0zZszQAw88oC5duuj8889XVlaWduzYoUWLFkmSHnvsMSUkJKhPnz6KiIjQ0qVLFR8fr7i4OC1cuFA+n08DBgxQ06ZN9eKLLyo6OlodOnQ4sx9MNRAWAACN3uTJkzV//nxdeeWVSkws/hTH4cOH9e9//1uLFy+u0D4iIkLXXHON5s+fHwgLn3/+ufr06RPUbsiQIVq9enWV5586dary8vL0m9/8RkeOHFFqaqpWrFihrl27SioOHY8++qj27t0rh8Ohfv366fXXX1dERITi4uI0Z84cTZ8+XT6fTz179tS//vUvtWrV6kx/LCEzzNpc7lnHPB6P3G638vLyFBsba3c5ANBonTp1Svv27VOnTp3UpEkTu8tBiKz+3qpzDWXNAgAAsGRrWPD5fJo1a5Y6deqk6OhodenSRQ8//HCtfrYVAACcGVvXLMydO1eZmZl67rnn1L17d23dulU333yz3G53nd/KEgAAVM7WsPDee+9p9OjRGjlypCSpY8eOeumll7R582Y7ywIAAGXYOg1x0UUXac2aNfr8888lSTt37tSGDRs0YsSIStt7vV55PJ6gDQBQd5gmblhq6+/L1pGFe+65Rx6PRykpKXI4HPL5fHrkkUeUlpZWafuMjAw9+OCDdVwlAMDhcEgqfjpidHS0zdUgVCdOnJAkRUZGntH72BoWXn75ZS1atEiLFy9W9+7dtWPHDk2bNk2JiYmaOHFihfYzZ87U9OnTA689Ho+SkpLqsmQACEult0H+9ttvFRkZqYgIPkxXn5mmqRMnTujIkSOKi4sLhL2asvU+C0lJSbrnnnsCN7yQpN///vd68cUX9dlnn1X5/dxnAQDqTkFBgfbt2ye/3293KQhRXFyc4uPjA7e3Lqs611BbRxZOnDhRIZ06HA7+QwSAeigqKkpdu3ZVQUGB3aUgBJGRkWc8olDK1rAwatQoPfLII0pOTlb37t314Ycf6rHHHtOkSZPsLAsAcBoRERHcwTEM2ToNkZ+fr1mzZmn58uU6cuSIEhMTNX78eN1///2Kioqq8vuZhgAAoGaqcw3l2RAAAIQhng0BAABqDWEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAkq1hoWPHjjIMo8KWnp5uZ1kAAKAMp50n37Jli3w+X+D1xx9/rMsvv1zjxo2zsSoAAFCWrWGhTZs2Qa/nzJmjLl26aPDgwZW293q98nq9gdcej+es1gcAAOrRmoWCggK9+OKLmjRpkgzDqLRNRkaG3G53YEtKSqrjKgEACD+GaZqm3UVI0ssvv6xf/OIXOnDggBITEyttU9nIQlJSkvLy8hQbG1tXpQIA0OB5PB653e6QrqG2TkOUNX/+fI0YMeK0QUGSXC6XXC5XHVYFAADqRVjYv3+/Vq9erWXLltldCgAAKKderFnIyspS27ZtNXLkSLtLAQAA5dgeFvx+v7KysjRx4kQ5nfVioAMAAJRhe1hYvXq1Dhw4oEmTJtldCgAAqITt/5S/4oorVE8+kAEAACph+8gCAACo3wgLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsGR7WDh06JBuuOEGtWrVStHR0erZs6e2bt1qd1kAAKCE086Tf//99xo0aJAuu+wyvfHGG2rTpo327t2rFi1a2FkWAAAow9awMHfuXCUlJSkrKyuwr1OnTjZWBAAAyrN1GmLFihXq27evxo0bp7Zt26pPnz569tlnT9ve6/XK4/EEbQAA4OyyNSx89dVXyszMVNeuXbVy5Urdeuutmjp1qp577rlK22dkZMjtdge2pKSkOq4YAIDwY5imadp18qioKPXt21fvvfdeYN/UqVO1ZcsWvf/++xXae71eeb3ewGuPx6OkpCTl5eUpNja2TmoGAKAx8Hg8crvdIV1DbR1ZSEhIUGpqatC+bt266cCBA5W2d7lcio2NDdoAAMDZZWtYGDRokPbs2RO07/PPP1eHDh1sqggAAJRna1i46667tGnTJv3hD3/QF198ocWLF+uZZ55Renq6nWUBAIAybA0L/fr10/Lly/XSSy+pR48eevjhhzVv3jylpaXZWRYAACjD1gWOZ6o6izMAAMAPGswCRwAAUP8RFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADAEmEBAABYIiwAAABLhAUAAGDJ1rAwe/ZsGYYRtKWkpNhZEgAAKMdpdwHdu3fX6tWrA6+dTttLAgAAZdh+ZXY6nYqPj7e7DAAAcBq2r1nYu3evEhMT1blzZ6WlpenAgQOnbev1euXxeII2AABwdtkaFgYMGKCFCxfqzTffVGZmpvbt26ef/OQnys/Pr7R9RkaG3G53YEtKSqrjigEACD+GaZqm3UWUys3NVYcOHfTYY49p8uTJFY57vV55vd7Aa4/Ho6SkJOXl5Sk2NrYuSwUAoEHzeDxyu90hXUNtX7NQVlxcnM4991x98cUXlR53uVxyuVx1XBUAAOHN9jULZR07dkxffvmlEhIS7C4FAACUsDUs/Pa3v9XatWv1n//8R++9956uueYaORwOjR8/3s6yAABAGbZOQxw8eFDjx4/Xd999pzZt2ujiiy/Wpk2b1KZNGzvLAgAAZdgaFpYsWWLn6QEAQAjq1ZoFAABQ/xAWAACAJcICAACwRFgAAACWCAsAAMASYQEAAFgiLAAAAEuEBQAAYKlGYSE7O1sHDx4MvN68ebOmTZumZ555ptYKAwAA9UONwsIvfvELvfPOO5KknJwcXX755dq8ebPuu+8+PfTQQ7VaIAAAsFeNwsLHH3+s/v37S5Jefvll9ejRQ++9954WLVqkhQsX1mZ9AADAZjUKC4WFhXK5XJKk1atX66qrrpIkpaSk6Jtvvqm96gAAgO1qFBa6d++up556SuvXr9eqVas0fPhwSdLXX3+tVq1a1WqBAADAXjUKC3PnztXTTz+tSy+9VOPHj1fv3r0lSStWrAhMTwAAgMbBME3TrMk3+nw+eTwetWjRIrDvP//5j5o2baq2bdvWWoFWPB6P3G638vLyFBsbWyfnBACgMajONbRGIwsnT56U1+sNBIX9+/dr3rx52rNnT50FBQAAUDdqFBZGjx6t559/XpKUm5urAQMG6E9/+pOuvvpqZWZm1mqBAADAXjUKC9u3b9dPfvITSdI//vEPtWvXTvv379fzzz+vJ554olYLBAAA9qpRWDhx4oRiYmIkSW+99ZbGjBmjiIgI/fjHP9b+/ftrtUAAAGCvGoWFH/3oR3rllVeUnZ2tlStX6oorrpAkHTlyhIWGAAA0MjUKC/fff79++9vfqmPHjurfv78GDhwoqXiUoU+fPrVaIAAAsFeNPzqZk5Ojb775Rr1791ZERHHm2Lx5s2JjY5WSklKrRZ4OH50EAKBmqnMNddb0JPHx8YqPjw88fbJ9+/bckAkAgEaoRtMQfr9fDz30kNxutzp06KAOHTooLi5ODz/8sPx+f23XCAAAbFSjkYX77rtP8+fP15w5czRo0CBJ0oYNGzR79mydOnVKjzzySK0WCQAA7FOjNQuJiYl66qmnAk+bLPXqq6/qtttu06FDh2qtQCusWQAAoGbO+u2ejx49WukixpSUFB09erQmbwkAAOqpGoWF3r1768knn6yw/8knn1SvXr3OuCgAAFB/1GjNwqOPPqqRI0dq9erVgXssvP/++8rOztbrr79eqwUCAAB71WhkYfDgwfr88891zTXXKDc3V7m5uRozZow++eQTvfDCC7VdIwAAsFGNb8pUmZ07d+qCCy6Qz+errbe0xAJHAABq5qwvcAQAAOGDsAAAACwRFgAAgKVqfRpizJgxlsdzc3PPpBYAAFAPVSssuN3uKo9PmDDhjAoCAAD1S7XCQlZW1tmqAwAA1FOsWQAAAJYICwAAwBJhAQAAWCIsAAAAS/UmLMyZM0eGYWjatGl2lwIAAMqoF2Fhy5Ytevrpp3m8NQAA9ZDtYeHYsWNKS0vTs88+qxYtWli29Xq98ng8QRsAADi7bA8L6enpGjlypIYOHVpl24yMDLnd7sCWlJRUBxUCABDebA0LS5Ys0fbt25WRkRFS+5kzZyovLy+wZWdnn+UKAQBAte7gWJuys7N15513atWqVWrSpElI3+NyueRyuc5yZQAAoCzDNE3TjhO/8soruuaaa+RwOAL7fD6fDMNQRESEvF5v0LHKeDweud1u5eXlKTY29myXDABAo1Gda6htIwtDhgzRrl27gvbdfPPNSklJ0d13311lUAAAAHXDtrAQExOjHj16BO1r1qyZWrVqVWE/AACwj+2fhgAAAPWbbSMLlXn33XftLgEAAJTDyAIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlW8NCZmamevXqpdjYWMXGxmrgwIF644037CwJAACUY2tYaN++vebMmaNt27Zp69at+ulPf6rRo0frk08+sbMsAABQhmGapml3EWW1bNlSf/zjHzV58uQq23o8HrndbuXl5Sk2NrYOqgMAoHGozjXUWUc1Vcnn82np0qU6fvy4Bg4cWGkbr9crr9cbeO3xeOqqPAAAwpbtCxx37dql5s2by+VyacqUKVq+fLlSU1MrbZuRkSG32x3YkpKS6rhaAADCj+3TEAUFBTpw4IDy8vL0j3/8Q3/729+0du3aSgNDZSMLSUlJTEMAAFBN1ZmGsD0slDd06FB16dJFTz/9dJVtWbMAAEDNVOcaavs0RHl+vz9o9AAAANjL1gWOM2fO1IgRI5ScnKz8/HwtXrxY7777rlauXGlnWQAAoAxbw8KRI0c0YcIEffPNN3K73erVq5dWrlypyy+/3M6yAABAGbaGhfnz59t5egAAEIJ6t2YBAADUL4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWAAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwkIlThb4VOTz210GAAD1AmGhnAkLNqvH7JXakZ1rdykAANQLhIVyohyGfH5THx3Ms7sUAADqBcJCOb3ax0mSPjqYa2sdAADUF4SFcnq1d0sSIwsAAJQgLJRTOrLw1X+PK+9kob3FAABQDxAWymnZLErtW0RLkj45xOgCAACEhUr0Lhld2MlUBAAAhIXK/LBuIdfeQgAAqAcIC5XoySJHAAACCAuV6HmOW4YhHco9qe+Oee0uBwAAWxEWKhHTJFKdWzeTxOgCAACEhdPoFVjkmGtrHQAA2I2wcBqlixx3MbIAAAhzTrsLqHdOfi8d/lRD8rfK5XxXzfeb2r9qmI627q8TzZNlGEZxO0MyZMgwJEOSYZT9urhB0LGS/aXfE2CaMnxemY4oyYgo3fXDYZV5UeFYRaXniSg5Sdlzlu4vX2PZfUV+U0U+U4U+f8nXfrmjI9W1XUxNf6IAgAaOsFDe/GHSf/coWdIvnCq+Im98Rx0kFZgOnZJLJxWlIjkkSaZ+uPKbplHpBbwyhqRow6tYHVeU4ZPfNJSvaOWrqUzTkGFUfKfSd//hz8r3lz0WqK2kVrOkVenXkuSUT01UIJcK5ZchnxwqUoR8cqhQDn1nxmrTqL/px/36hdg7AEBjQlgor113qfCk1C5VG/Pb6fPD+bpAu9XNv1dRRpGidEKxOlH595a/QldDhGHKrRNy68QZvc/ZkGgc1Xt735IICwAQlggL5Y15VnIU/1gGlWySpMJT0onvioNE4QnJ9JXMCZT8a96UKp8YsBDZVGrillzNpSKvdCpP8np+eJtAaCiTHgwjeN9pX5f9PvOHWiv8qeL+OqMlZ1Txfr9P8hdJ/iJ98uIMdT/2nowiPkIKAOGKsFCe4zQ/ksgmkvucs3vu5m3P7vvXwPGoVpIk01dgcyUAALvwaQhYi4iUJBlFhAUACFe2hoWMjAz169dPMTExatu2ra6++mrt2bPHzpJQjumIKv7TT1gAgHBla1hYu3at0tPTtWnTJq1atUqFhYW64oordPz4cTvLQhmlYcFgGgIAwpataxbefPPNoNcLFy5U27ZttW3bNl1yySU2VYUggbDAAkcACFf1aoFjXl7x3RJbtmxZ6XGv1yuv94eLlsfjqZO6wprTJUkyfIU2FwIAsEu9WeDo9/s1bdo0DRo0SD169Ki0TUZGhtxud2BLSkqq4yrDj1EyshDBmgUACFv1Jiykp6fr448/1pIlS07bZubMmcrLywts2dnZdVhheDKcpWGBkQUACFf1Yhri9ttv12uvvaZ169apffv2p23ncrnkcrnqsDLIycgCAIQ7W8OCaZq64447tHz5cr377rvq1KmTneWgEhElaxYYWQCA8GVrWEhPT9fixYv16quvKiYmRjk5OZIkt9ut6OhoO0tDCcPZRJLkICwAQNiydc1CZmam8vLydOmllyohISGw/f3vf7ezLJQRUToNYRIWACBc2T4NgfrNEVk8DeFkZAEAwla9+TQE6qeI0rBgssARAMIVYQGWSkcWHGaRzZUAAOxCWICliMjiNQtOMQ0BAOGKsABLzqjiT0M4WeAIAGGLsABLzsjisBDJyAIAhC3CAiw5o4rXLESyZgEAwhZhAZZKRxaiRFgAgHBFWIAlZ8mzOFxGoUy/3+ZqAAB2ICzAUqTrh9tuFxRyrwUACEeEBViKivzhKZ8F3pM2VgIAsAthAZbKjiwUFTCyAADhiLAASw5npHymIUkq9J6yuRoAgB0IC6hSYcnzxooKCQsAEI4IC6hSgREpSSosICwAQDgiLKBKRSUjC75Cr82VAADsQFhAlQpVPLJQVEBYAIBwRFhAlYqM0pEFpiEAIBwRFlClIqP4MdU+RhYAICwRFlClopIFjr4iwgIAhCPCAqrkK5mG8LPAEQDCEmEBVSqKKJ6G8DOyAABhibCAKvlLpiEICwAQnggLqJIvojgsmExDAEBYIiygSmZE6cgCD5ICgHBEWECVfCVrFkzCAgCEJcICquQvCQvyMQ0BAOGIsIAqmQ5GFgAgnBEWUKXSNQuGj7AAAOGIsICqOUqnIQgLABCOCAuoUuk0BCMLABCeCAuoWiAssMARAMIRYQFVc7gkSYa/0OZCAAB2ICygSoaTaQgACGeEBVStJCxEMLIAAGGJsIAqGc7iaYgIPyMLABCOCAuoksHIAgCENcICqhThbCJJchAWACAsERZQpYjSkQWTsAAA4YiwgCpFRBavWXASFgAgLBEWUCVHZPE0hJMFjgAQlggLqFLpNIRDRTZXAgCwg61hYd26dRo1apQSExNlGIZeeeUVO8vBaTiiiqchIpmGAICwZGtYOH78uHr37q2//vWvdpaBKjhYswAAYc1p58lHjBihESNGhNze6/XK6/3hYUYej+dslIVyIqOK1yxEirAAAOGoQa1ZyMjIkNvtDmxJSUl2lxQWfhhZYM0CAISjBhUWZs6cqby8vMCWnZ1td0lhoXRkIYoFjgAQlmydhqgul8sll8tldxlhx+kqDgsuo1Cm3y8jokFlTADAGeK3PqrkLBlZkKTCQtYtAEC4ISygSi7XD2GhoOCkjZUAAOxg6zTEsWPH9MUXXwRe79u3Tzt27FDLli2VnJxsY2UoK7LMyEKR1yvF2FgMAKDO2RoWtm7dqssuuyzwevr06ZKkiRMnauHChTZVhfIczkj5TEMOw1RhwSm7ywEA1DFbw8Kll14q0zTtLAEhKpRTDhUSFgAgDLFmASEpMCIlSUWFhAUACDeEBYSkUKVhwVtFSwBAY0NYQEiKSmasfAWEBQAIN4QFhKSoZBrCx5oFAAg7hAWEpMhgGgIAwhVhASEJjCwQFgAg7BAWEBKfUbxmwV9UYHMlAIC6RlhASIoioiRJfj46CQBhh7CAkPhKpiFMRhYAIOwQFhASf0RxWPCzZgEAwg5hASHxl0xDmEWEBQAIN4QFhKR0ZIFpCAAIP4QFhKR0ZEE+RhYAINwQFhAS08HIAgCEK8ICQmKWjCwYPsICAIQbwgJCYjpcxX8SFgAg7BAWEJqSaQhGFgAg/BAWEJLSkQWDBY4AEHYICwiJ4ShZs+AvtLkSAEBdIywgNM7isBDBNAQAhB3CAkJiOBlZAIBwRVhAaErWLET4GVkAgHBDWEBIIiJLwwIjCwAQbggLCEmEszgsOEzCAgCEG8ICQhIRWbxmwcHIAgCEHcICQsLIAgCEL8ICQlK6ZsHJAkcACDuEBYTEUTqyoCKbKwEA1DXCAkISEVUcFiKZhgCAsENYQEickU2K/yQsAEDYISwgJM7SkQWmIQAg7BAWEBJGFgAgfBEWEJLIkpGFKEYWACDsEBYQEkdU8ciCyyiU6ffbXA0AoC4RFhCSSFeTwNdFRUxFAEA4ISwgJFFRP4SFAu9JGysBANQ1wgJCEuWKDnxdVOC1sRIAQF0jLCAkDqdTPtOQJBV6T9lcDQCgLhEWELICRUqSCgsJCwAQTggLCFmh4ZQkFRUQFgAgnNSLsPDXv/5VHTt2VJMmTTRgwABt3rzZ7pJQicLAyAJrFgAgnNgeFv7+979r+vTpeuCBB7R9+3b17t1bw4YN05EjR+wuDeUUqXhkwccCRwAIK4ZpmqadBQwYMED9+vXTk08+KUny+/1KSkrSHXfcoXvuucfyez0ej9xut/Ly8hQbG1sX5Ya1gw+ep/Zmjt7vlK4m8efZXQ4AhJ0el10fuKPumarONdRZK2esoYKCAm3btk0zZ84M7IuIiNDQoUP1/vvvV2jv9Xrl9f7wr9q8vDxJxR3G2XfUG6lY01T3z56UPrO7GgAIP4e7DVZsXMtaea/Sa2coYwa2hoX//ve/8vl8ateuXdD+du3a6bPPKl6NMjIy9OCDD1bYn5SUdNZqBACg3pjTqdbfMj8/X26327KNrWGhumbOnKnp06cHXvv9fh09elStWrWSYRi1cg6Px6OkpCRlZ2c32qkN+tg4NPY+Nvb+SfSxMWjI/TNNU/n5+UpMTKyyra1hoXXr1nI4HDp8+HDQ/sOHDys+Pr5Ce5fLJZcreK4mLi7urNQWGxvb4P7iq4s+Ng6NvY+NvX8SfWwMGmr/qhpRKGXrpyGioqJ04YUXas2aNYF9fr9fa9as0cCBA22sDAAAlLJ9GmL69OmaOHGi+vbtq/79+2vevHk6fvy4br75ZrtLAwAAqgdh4brrrtO3336r+++/Xzk5OTr//PP15ptvVlj0WFdcLpceeOCBCtMdjQl9bBwaex8be/8k+tgYNPb+lbL9PgsAAKB+s/0OjgAAoH4jLAAAAEuEBQAAYImwAAAALBEWymksj8vOyMhQv379FBMTo7Zt2+rqq6/Wnj17gtqcOnVK6enpatWqlZo3b66xY8dWuEFWQzJnzhwZhqFp06YF9jWGPh46dEg33HCDWrVqpejoaPXs2VNbt24NHDdNU/fff78SEhIUHR2toUOHau/evTZWHDqfz6dZs2apU6dOio6OVpcuXfTwww8H3au+ofVv3bp1GjVqlBITE2UYhl555ZWg46H05+jRo0pLS1NsbKzi4uI0efJkHTt2rA57Yc2qj4WFhbr77rvVs2dPNWvWTImJiZowYYK+/vrroPdoyH0sb8qUKTIMQ/PmzQvaX9/7WB2EhTIa0+Oy165dq/T0dG3atEmrVq1SYWGhrrjiCh0/fjzQ5q677tK//vUvLV26VGvXrtXXX3+tMWPG2Fh1zW3ZskVPP/20evXqFbS/offx+++/16BBgxQZGak33nhDn376qf70pz+pRYsWgTaPPvqonnjiCT311FP64IMP1KxZMw0bNkynTp2ysfLQzJ07V5mZmXryySe1e/duzZ07V48++qj+8pe/BNo0tP4dP35cvXv31l//+tdKj4fSn7S0NH3yySdatWqVXnvtNa1bt0633HJLXXWhSlZ9PHHihLZv365Zs2Zp+/btWrZsmfbs2aOrrroqqF1D7mNZy5cv16ZNmyq9ZXJ972O1mAjo37+/mZ6eHnjt8/nMxMREMyMjw8aqaseRI0dMSebatWtN0zTN3NxcMzIy0ly6dGmgze7du01J5vvvv29XmTWSn59vdu3a1Vy1apU5ePBg88477zRNs3H08e677zYvvvji0x73+/1mfHy8+cc//jGwLzc313S5XOZLL71UFyWekZEjR5qTJk0K2jdmzBgzLS3NNM2G3z9J5vLlywOvQ+nPp59+akoyt2zZEmjzxhtvmIZhmIcOHaqz2kNVvo+V2bx5synJ3L9/v2majaePBw8eNM855xzz448/Njt06GA+/vjjgWMNrY9VYWShROnjsocOHRrYZ/W47Iam9HHeLVsWP9p027ZtKiwsDOpvSkqKkpOTG1x/09PTNXLkyKC+SI2jjytWrFDfvn01btw4tW3bVn369NGzzz4bOL5v3z7l5OQE9dHtdmvAgAENoo8XXXSR1qxZo88//1yStHPnTm3YsEEjRoyQ1PD7V14o/Xn//fcVFxenvn37BtoMHTpUERER+uCDD+q85tqQl5cnwzACz/JpDH30+/268cYbNWPGDHXv3r3C8cbQx7Jsv4NjfVHdx2U3JH6/X9OmTdOgQYPUo0cPSVJOTo6ioqIqPIirXbt2ysnJsaHKmlmyZIm2b9+uLVu2VDjWGPr41VdfKTMzU9OnT9e9996rLVu2aOrUqYqKitLEiRMD/ajsv9uG0Md77rlHHo9HKSkpcjgc8vl8euSRR5SWliZJDb5/5YXSn5ycHLVt2zbouNPpVMuWLRtkn0+dOqW7775b48ePDzxoqTH0ce7cuXI6nZo6dWqlxxtDH8siLISB9PR0ffzxx9qwYYPdpdSq7Oxs3XnnnVq1apWaNGlidzlnhd/vV9++ffWHP/xBktSnTx99/PHHeuqppzRx4kSbqztzL7/8shYtWqTFixere/fu2rFjh6ZNm6bExMRG0b9wV1hYqGuvvVamaSozM9PucmrNtm3b9Oc//1nbt2+XYRh2l1MnmIYoUd3HZTcUt99+u1577TW98847at++fWB/fHy8CgoKlJubG9S+IfV327ZtOnLkiC644AI5nU45nU6tXbtWTzzxhJxOp9q1a9fg+5iQkKDU1NSgfd26ddOBAwckKdCPhvrf7YwZM3TPPffo+uuvV8+ePXXjjTfqrrvuUkZGhqSG37/yQulPfHx8hUXVRUVFOnr0aIPqc2lQ2L9/v1atWhX0+OaG3sf169fryJEjSk5ODvzu2b9/v37zm9+oY8eOkhp+H8sjLJRobI/LNk1Tt99+u5YvX663335bnTp1Cjp+4YUXKjIyMqi/e/bs0YEDBxpMf4cMGaJdu3Zpx44dga1v375KS0sLfN3Q+zho0KAKH3n9/PPP1aFDB0lSp06dFB8fH9RHj8ejDz74oEH08cSJE4qICP415HA45Pf7JTX8/pUXSn8GDhyo3Nxcbdu2LdDm7bfflt/v14ABA+q85pooDQp79+7V6tWr1apVq6DjDb2PN954oz766KOg3z2JiYmaMWOGVq5cKanh97ECu1dY1idLliwxXS6XuXDhQvPTTz81b7nlFjMuLs7Mycmxu7Rqu/XWW023222+++675jfffBPYTpw4EWgzZcoUMzk52Xz77bfNrVu3mgMHDjQHDhxoY9VnruynIUyz4fdx8+bNptPpNB955BFz79695qJFi8ymTZuaL774YqDNnDlzzLi4OPPVV181P/roI3P06NFmp06dzJMnT9pYeWgmTpxonnPOOeZrr71m7tu3z1y2bJnZunVr83e/+12gTUPrX35+vvnhhx+aH374oSnJfOyxx8wPP/ww8EmAUPozfPhws0+fPuYHH3xgbtiwwezatas5fvx4u7pUgVUfCwoKzKuuusps3769uWPHjqDfP16vN/AeDbmPlSn/aQjTrP99rA7CQjl/+ctfzOTkZDMqKsrs37+/uWnTJrtLqhFJlW5ZWVmBNidPnjRvu+02s0WLFmbTpk3Na665xvzmm2/sK7oWlA8LjaGP//rXv8wePXqYLpfLTElJMZ955pmg436/35w1a5bZrl070+VymUOGDDH37NljU7XV4/F4zDvvvNNMTk42mzRpYnbu3Nm87777gi4qDa1/77zzTqX/702cONE0zdD6891335njx483mzdvbsbGxpo333yzmZ+fb0NvKmfVx3379p32988777wTeI+G3MfKVBYW6nsfq4NHVAMAAEusWQAAAJYICwAAwBJhAQAAWCIsAAAAS4QFAABgibAAAAAsERYAAIAlwgIAALBEWABQ7xiGoVdeecXuMgCUICwACHLTTTfJMIwK2/Dhw+0uDYBNnHYXAKD+GT58uLKysoL2uVwum6oBYDdGFgBU4HK5FB8fH7S1aNFCUvEUQWZmpkaMGKHo6Gh17txZ//jHP4K+f9euXfrpT3+q6OhotWrVSrfccouOHTsW1GbBggXq3r27XC6XEhISdPvttwcd/+9//6trrrlGTZs2VdeuXbVixYqz22kAp0VYAFBts2bN0tixY7Vz506lpaXp+uuv1+7duyVJx48f17Bhw9SiRQtt2bJFS5cu1erVq4PCQGZmptLT03XLLbdo165dWrFihX70ox8FnePBBx/Utddeq48++khXXnml0tLSdPTo0TrtJ4ASdj/2EkD9MnHiRNPhcJjNmjUL2h555BHTNIsffz5lypSg7xkwYIB56623mqZpms8884zZokUL89ixY4Hj//73v82IiAgzJyfHNE3TTExMNO+7777T1iDJ/J//+Z/A62PHjpmSzDfeeKPW+gkgdKxZAFDBZZddpszMzKB9LVu2DHw9cODAoGMDBw7Ujh07JEm7d+9W79691axZs8DxQYMGye/3a8+ePTIMQ19//bWGDBliWUOvXr0CXzdr1kyxsbE6cuRITbsE4AwQFgBU0KxZswrTArUlOjo6pHaRkZFBrw3DkN/vPxslAagCaxYAVNumTZsqvO7WrZskqVu3btq5c6eOHz8eOL5x40ZFRETovPPOU0xMjDp27Kg1a9bUac0Aao6RBQAVeL1e5eTkBO1zOp1q3bq1JGnp0qXq27evLr74Yi1atEibN2/W/PnzJUlpaWl64IEHNHHiRM2ePVvffvut7rjjDt14441q166dJGn27NmaMmWK2rZtqxEjRig/P18bN27UHXfcUbcdBRASwgKACt58800lJCQE7TvvvPP02WefSSr+pMKSJUt02223KSEhQS+99JJSU1MlSU2bNtXKlSt15513ql+/fmratKnGjh2rxx57LPBeEydO1KlTp/T444/rt7/9rVq3bq2f//zndddBANVimKZp2l0EgIbDMAwtX75cV199td2lAKgjrFkAAACWCAsAAMASaxYAVAszl0D4YWQBAABYIiwAAABLhAUAAGCJsAAAACwRFgAAgCXCAgAAsERYAAAAlggLAADA0v8D1ZVoPZ+TjY0AAAAASUVORK5CYII=\n"},"metadata":{}}]},{"cell_type":"code","source":["def smape(y_true, y_pred):\n","    \"\"\"\n","    Calculate Symmetric Mean Absolute Percentage Error (SMAPE).\n","    \"\"\"\n","    denominator = (np.abs(y_true) + np.abs(y_pred)) / 2.0\n","    diff = np.abs(y_true - y_pred) / denominator\n","    diff[denominator == 0] = 0.0  # Handle division by zero\n","    return np.mean(diff)*100/ len(y_true)\n","\n","def mape(y_true, y_pred):\n","    \"\"\"\n","    Calculate Mean Absolute Percentage Error (MAPE).\n","    \"\"\"\n","    diff = np.abs(y_true - y_pred) / np.abs(y_true)\n","    diff[np.isnan(diff)] = 0.0  # Handle NaN values\n","    diff[np.isinf(diff)] = 0.0  # Handle infinite values\n","    return np.mean(diff) * 100 / len(y_true)\n","\n","def theil_coefficient(y_true, y_pred):\n","    \"\"\"\n","    Calculate Theil's Inequality Coefficient.\n","    \"\"\"\n","    n = len(y_true)\n","\n","    # Calculate three components\n","    mean_true = np.mean(y_true)\n","    mean_pred = np.mean(y_pred)\n","    diff_true = y_true - mean_true\n","    diff_pred = y_pred - mean_pred\n","\n","    # Calculate Theil's Inequality Coefficient\n","    num = np.sqrt(np.mean(diff_true**2)) + np.sqrt(np.mean(diff_pred**2))\n","    denom = np.sqrt(np.mean(y_true**2)) + np.sqrt(np.mean(y_pred**2))\n","\n","    return num / denom"],"metadata":{"id":"c-UxwkxPuvSo"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["batch_size = 64  # You can adjust this value as needed\n","\n","# Convert X_2007 and y_2007 to tensors\n","X_2007_tensor = torch.tensor(X_2007, dtype=torch.float32)\n","y_2007_tensor = torch.tensor(y_2007, dtype=torch.float32)\n","\n","# Calculate the total number of batches\n","num_samples = X_2007_tensor.size(0)\n","num_batches = (num_samples + batch_size - 1) // batch_size\n","\n","# Initialize lists to store predictions and MSE losses\n","predictions = []\n","mse_losses = []\n","smape_values = []\n","mape_values = []\n","theil_coefficients = []\n","\n","# Iterate over batches\n","for i in range(num_batches):\n","    # Calculate start and end indices for the current batch\n","    start_idx = i * batch_size\n","    end_idx = min((i + 1) * batch_size, num_samples)\n","\n","    # Extract the current batch from X_2007_tensor and y_2007_tensor\n","    X_batch = X_2007_tensor[start_idx:end_idx]\n","    y_batch = y_2007_tensor[start_idx:end_idx]\n","\n","    # Pass X_batch through the model to obtain predictions\n","    u_pred_batch = net.forward(X_batch)\n","    predictions.append(u_pred_batch)\n","\n","    # If you have ground truth data for 2007, calculate MSE loss for the current batch\n","    if y_batch is not None:\n","        mse_loss = torch.nn.functional.mse_loss(y_batch, u_pred_batch)\n","        mse_losses.append(mse_loss.item())\n","\n","                # Calculate SMAPE\n","        smape_value = smape(y_batch.detach().numpy(), u_pred_batch.detach().numpy())\n","        smape_values.append(smape_value)\n","        mape_value = mape(y_batch.detach().numpy(), u_pred_batch.detach().numpy())\n","        mape_values.append(mape_value)\n","\n","\n","        # Calculate Theil coefficient\n","        theil_coefficient_value = theil_coefficient(y_batch.detach().numpy(), u_pred_batch.detach().numpy())\n","        theil_coefficients.append(theil_coefficient_value)\n","\n","# Concatenate predictions and calculate overall MSE loss if applicable\n","u_pred_2007 = torch.cat(predictions, dim=0)\n","val_mse_loss = None\n","if y_2007_tensor is not None:\n","    val_mse_loss = sum(mse_losses) / len(mse_losses)\n","    print(\"Validation MSE for 2007:\", val_mse_loss)\n","    print(\"Validation SMAPE for 2007:\", np.mean(smape_values))\n","    print(\"Validation MAPE for 2007:\", np.mean(mape_values))\n","    print(\"Validation Theil coefficient for 2007:\", np.mean(theil_coefficients))\n"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"6ZL87oEbpToC","executionInfo":{"status":"ok","timestamp":1718415039125,"user_tz":240,"elapsed":4380,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"}},"outputId":"8f907d25-1a57-4b53-ac50-450e51a7ac40"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Validation MSE for 2007: 1.9324569610990412\n","Validation SMAPE for 2007: 0.981034422557967\n","Validation MAPE for 2007: 6.8397884123465795\n","Validation Theil coefficient for 2007: 0.1814141\n"]}]},{"cell_type":"markdown","source":["        #without KAN LSTM:\n","                self.net_optim = torch.optim.Adam(self.net.parameters(), lr=0.001, betas = (0.9, 0.999))\n","        Validation MSE for 2007: 2.028034624192984\n","Validation SMAPE for 2007: 1.1119537088753921\n","Validation MAPE for 2007: 4.363998836181778\n","        \n","        "],"metadata":{"id":"cjKVJqTQQ50S"}},{"cell_type":"markdown","source":["With KAN LSTM:\n","        #self.net_optim = torch.optim.Adam(self.net.parameters(), lr=0.001, betas = (0.9, 0.999))\n","        Validation MSE for 2007: 1.9362252502775137\n","Validation SMAPE for 2007: 0.9812685286289313\n","Validation MAPE for 2007: 6.853347201321986\n","Validation Theil coefficient for 2007: 0.1813422\n","     \n","      "],"metadata":{"id":"KH9HC1RtX2aB"}},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":24,"status":"ok","timestamp":1718415039125,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"},"user_tz":240},"id":"pbYRUXjk1bE9","outputId":"ea90a3f1-0b9e-43ad-ddb3-98dda5855bc9"},"outputs":[{"output_type":"stream","name":"stdout","text":["(55536, 3)\n","(55536,)\n"]}],"source":["print(X_2007.shape)\n","print(y_2007.shape)"]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":23,"status":"ok","timestamp":1718415039126,"user":{"displayName":"Akila Sampath","userId":"11204871320919696111"},"user_tz":240},"id":"VUSdMLDi6bxd","outputId":"5ec64efe-ddc8-4d97-abf9-ae48f7b7e8d9"},"outputs":[{"output_type":"stream","name":"stdout","text":["torch.Size([55536, 1])\n"]}],"source":["print(u_pred_2007.shape)"]}],"metadata":{"colab":{"provenance":[],"machine_shape":"hm"},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":0}