MQRNN end-to-end demo in PyTorch

MQRNN.ipynb

{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nfrom tqdm import trange, tqdm\n\nfrom io import BytesIO\nfrom urllib.request import urlopen\nfrom zipfile import ZipFile\n\nfrom pandas import read_csv\nfrom scipy import stats","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:29.893322Z","iopub.execute_input":"2023-01-10T06:54:29.893743Z","iopub.status.idle":"2023-01-10T06:54:30.618207Z","shell.execute_reply.started":"2023-01-10T06:54:29.893660Z","shell.execute_reply":"2023-01-10T06:54:30.617115Z"},"trusted":true},"execution_count":1,"outputs":[]},{"cell_type":"code","source":"name = 'LD2011_2014.txt'\nsave_name = 'elect'\nnum_covariates = 3\ntrain_start = '2011-01-01 00:00:00'\ntrain_end = '2014-08-31 23:00:00'\ntest_start = '2014-08-24 00:00:00' #need additional 7 days as given info\ntest_end = '2014-09-07 23:00:00'\n\nsave_path = os.path.join('data', save_name)\nif not os.path.exists(save_path):\n    os.makedirs(save_path)\ncsv_path = os.path.join(save_path, name)\nif not os.path.exists(csv_path):\n    zipurl = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00321/LD2011_2014.txt.zip'\n    with urlopen(zipurl) as zipresp:\n        with ZipFile(BytesIO(zipresp.read())) as zfile:\n            zfile.extractall(save_path)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:30.773982Z","iopub.execute_input":"2023-01-10T06:54:30.774640Z","iopub.status.idle":"2023-01-10T06:54:38.642682Z","shell.execute_reply.started":"2023-01-10T06:54:30.774594Z","shell.execute_reply":"2023-01-10T06:54:38.641667Z"},"trusted":true},"execution_count":3,"outputs":[]},{"cell_type":"code","source":"# initial fixing\ntarget_df = pd.read_csv(csv_path, sep=\";\", index_col=0, parse_dates=True, decimal=',')\ntarget_df = target_df.resample('12H',label = 'left',closed = 'right').sum()[train_start:test_end]\ntarget_df.fillna(0, inplace=True)\n\ntrain_target_df = target_df[train_start:train_end]\ntest_target_df = target_df[test_start:test_end]\nhorizon_size = test_target_df.shape[0]","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:38.645423Z","iopub.execute_input":"2023-01-10T06:54:38.645782Z","iopub.status.idle":"2023-01-10T06:54:46.858712Z","shell.execute_reply.started":"2023-01-10T06:54:38.645746Z","shell.execute_reply":"2023-01-10T06:54:46.857768Z"},"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\nscaler.fit(train_target_df)\ntrain_target_df = pd.DataFrame(scaler.transform(train_target_df), index=train_target_df.index, columns=train_target_df.columns)\ntest_target_df = pd.DataFrame(scaler.transform(test_target_df), index=test_target_df.index, columns=test_target_df.columns)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:46.871640Z","iopub.execute_input":"2023-01-10T06:54:46.872534Z","iopub.status.idle":"2023-01-10T06:54:46.998146Z","shell.execute_reply.started":"2023-01-10T06:54:46.872496Z","shell.execute_reply":"2023-01-10T06:54:46.997141Z"},"trusted":true},"execution_count":6,"outputs":[]},{"cell_type":"code","source":"covariate_df = pd.DataFrame(index=target_df.index,\n                                data={'hour':target_df.index.hour,\n                                      'dayofweek':target_df.index.dayofweek,\n                                      'month': target_df.index.month\n                                })\nfor col in covariate_df.columns:\n        covariate_df[col] = (covariate_df[col] - np.mean(covariate_df[col]))/np.std(covariate_df[col])\n\ntrain_covariate_df = covariate_df[train_start:train_end]\ntest_covariate_df = covariate_df[test_start:test_end]","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:46.999682Z","iopub.execute_input":"2023-01-10T06:54:47.000147Z","iopub.status.idle":"2023-01-10T06:54:47.015969Z","shell.execute_reply.started":"2023-01-10T06:54:47.000108Z","shell.execute_reply":"2023-01-10T06:54:47.014937Z"},"trusted":true},"execution_count":7,"outputs":[]},{"cell_type":"code","source":"import torch\nfrom torch.utils.data import Dataset\n\ndesired_quantiles = [0.25, 0.5, 0.75, 0.95]\n\nclass MQRNN_dataset(Dataset):\n    \n    def __init__(self,\n                series_df:pd.DataFrame,\n                covariate_df:pd.DataFrame, \n                horizon_size:int=horizon_size,\n                quantile_size:int=len(desired_quantiles)):\n        \n        self.series_df = series_df\n        self.covariate_df = covariate_df\n        self.horizon_size = horizon_size\n        self.quantile_size = quantile_size\n        full_covariate = []\n        covariate_size = self.covariate_df.shape[1]\n        for i in range(1, self.covariate_df.shape[0] - horizon_size+1):\n            cur_covariate = []\n            cur_covariate.append(self.covariate_df.iloc[i:i+horizon_size,:].to_numpy())\n            full_covariate.append(cur_covariate)\n        full_covariate = np.array(full_covariate)\n        full_covariate = full_covariate.reshape(-1, horizon_size * covariate_size)\n        self.next_covariate = full_covariate\n    \n    def __len__(self):\n        return self.series_df.shape[1]\n    \n    def __getitem__(self,idx):\n        cur_series = np.array(self.series_df.iloc[: -self.horizon_size, idx])\n        cur_covariate = np.array(self.covariate_df.iloc[:-self.horizon_size, :]) # covariate used in generating hidden states\n\n        covariate_size = self.covariate_df.shape[1]\n\n        real_vals_list = []\n        for i in range(1, self.horizon_size+1):\n            real_vals_list.append(np.array(self.series_df.iloc[i: self.series_df.shape[0]-self.horizon_size+i, idx]))\n        real_vals_array = np.array(real_vals_list) #[horizon_size, seq_len]\n        real_vals_array = real_vals_array.T #[seq_len, horizon_size]\n        cur_series_tensor = torch.tensor(cur_series)\n        \n        cur_series_tensor = torch.unsqueeze(cur_series_tensor,dim=1) # [seq_len, 1]\n        cur_covariate_tensor = torch.tensor(cur_covariate) #[seq_len, covariate_size]\n        cur_series_covariate_tensor = torch.cat([cur_series_tensor, cur_covariate_tensor],dim=1)\n        next_covariate_tensor = torch.tensor(self.next_covariate) #[seq_len, horizon_size * covariate_size]\n\n        cur_real_vals_tensor = torch.tensor(real_vals_array)\n        return cur_series_covariate_tensor, next_covariate_tensor, cur_real_vals_tensor","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:47.026796Z","iopub.execute_input":"2023-01-10T06:54:47.027505Z","iopub.status.idle":"2023-01-10T06:54:48.660729Z","shell.execute_reply.started":"2023-01-10T06:54:47.027469Z","shell.execute_reply":"2023-01-10T06:54:48.659715Z"},"trusted":true},"execution_count":9,"outputs":[]},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\n\nclass MQRNN(nn.Module):\n    def __init__(self, \n                horizon_size:int=horizon_size, \n                hidden_size:int=100, \n                quantiles:list=desired_quantiles,\n                dropout:float=0.3,\n                layer_size:int=3,\n                context_size:int=50, \n                covariate_size:int=num_covariates,\n                bidirectional=False,\n                device=torch.device('cuda')):\n        super(MQRNN, self).__init__()\n        self.quantiles = desired_quantiles\n        self.quantile_size = len(quantiles)\n        self.bidirectional = bidirectional\n        self.hidden_size = hidden_size\n        self.horizon_size = horizon_size\n        self.device = device\n        self.covariate_size = covariate_size\n\n        self.encoder = nn.LSTM(input_size= covariate_size+1, \n                            hidden_size=hidden_size, \n                            num_layers=layer_size, \n                            dropout=dropout,\n                            bidirectional=bidirectional)\n        self.global_decoder = nn.Sequential(nn.Linear(in_features= hidden_size + covariate_size*horizon_size, out_features= horizon_size*hidden_size*3),\n                                           nn.ReLU(),\n                                           nn.Linear(in_features= horizon_size*hidden_size*3, out_features= horizon_size*hidden_size*2),\n                                           nn.ReLU(),\n                                           nn.Linear(in_features= horizon_size*hidden_size*2, out_features= (horizon_size+1)*context_size),\n                                           nn.ReLU())\n        self.local_decoder = nn.Sequential(nn.Linear(in_features= horizon_size*context_size + horizon_size* covariate_size + context_size, out_features= horizon_size* context_size),\n                                           nn.ReLU(),\n                                           nn.Linear(in_features= horizon_size* context_size, out_features= horizon_size*self.quantile_size),\n                                           nn.ReLU())\n        self.encoder.double().to(self.device)\n        self.global_decoder.double().to(self.device)\n        self.local_decoder.double().to(self.device)\n        \n    def forward(self, cur_series_covariate_tensor, next_covariate_tensor):\n        seq_len, batch_size = cur_series_covariate_tensor.shape[0], cur_series_covariate_tensor.shape[1]\n        direction_size = 2 if self.bidirectional else 1\n        outputs,_ = self.encoder(cur_series_covariate_tensor)\n        outputs_reshape = outputs.view(seq_len,batch_size,direction_size,self.hidden_size)[:,:,-1,:]\n        encoder_outputs = outputs_reshape.view(seq_len,batch_size,self.hidden_size)\n        if not self.training:\n            encoder_outputs = torch.unsqueeze(encoder_outputs[-1],dim=0) #[1,1,hidden_size]\n\n        global_decoder_output = self.global_decoder(torch.cat([encoder_outputs, next_covariate_tensor], dim=2))\n        local_decoder_output = self.local_decoder(torch.cat([global_decoder_output, next_covariate_tensor], dim=2))\n        \n        seq_len = local_decoder_output.shape[0]\n        batch_size = local_decoder_output.shape[1]\n        return local_decoder_output.view(seq_len, batch_size, self.horizon_size, self.quantile_size)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:48.666205Z","iopub.execute_input":"2023-01-10T06:54:48.666687Z","iopub.status.idle":"2023-01-10T06:54:48.684984Z","shell.execute_reply.started":"2023-01-10T06:54:48.666656Z","shell.execute_reply":"2023-01-10T06:54:48.683862Z"},"trusted":true},"execution_count":10,"outputs":[]},{"cell_type":"code","source":"from torch.utils.data import DataLoader\nfrom tqdm import tqdm, trange\n\ndef train(model, train_dataset, desired_quantiles=desired_quantiles, batch_size=1, num_epochs=1, lr=1e-3, device=torch.device(\"cuda\")):\n    optimizer = torch.optim.Adam(model.parameters(),lr=lr)\n    data_iter = DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True)\n    \n    model.train()\n    for epoch in range(num_epochs):\n        epoch_loss_sum = 0.0\n        total_sample = 0\n        pbar = tqdm(data_iter)\n        for (cur_series_tensor, cur_covariate_tensor, cur_real_vals_tensor) in pbar:\n            batch_size, seq_len, horizon_size  = cur_series_tensor.shape[0], cur_series_tensor.shape[1], cur_covariate_tensor.shape[-1]\n            total_sample += batch_size * seq_len * horizon_size\n            optimizer.zero_grad()\n            cur_series_covariate_tensor = cur_series_tensor.double().permute(1,0,2).to(device) #[seq_len, batch_size, 1+covariate_size]\n            next_covariate_tensor = cur_covariate_tensor.double().permute(1,0,2).to(device) #[seq_len, batch_size, covariate_size * horizon_size]\n            cur_real_vals_tensor = cur_real_vals_tensor.double().permute(1,0,2).to(device) #[seq_len, batch_size, horizon_size]\n            \n            model_output = model(cur_series_covariate_tensor, next_covariate_tensor)\n            \n            losses = []\n            for i, p in enumerate(desired_quantiles):\n                errors = cur_real_vals_tensor - model_output[:,:,:,i]\n                losses.append(torch.max((p-1)*errors, p*errors))\n            total_loss = torch.mean(torch.sum(torch.cat(losses, dim=1), dim=1)).to(device)\n            \n#             total_loss = torch.tensor([0.0]).to(device)\n#             for i in range(len(desired_quantiles)):\n#                 p = desired_quantiles[i]\n#                 errors = cur_real_vals_tensor - model_output[:,:,:,i]\n#                 cur_loss = torch.max((p-1)*errors, p*errors)\n#                 total_loss += torch.sum(cur_loss)\n\n            pbar.set_description(f\"Loss:{total_loss.item()}\")\n            total_loss.backward()\n            optimizer.step()\n            epoch_loss_sum += total_loss.item()\n        epoch_loss_mean = epoch_loss_sum/ total_sample\n        if (epoch+1)%5 == 0:\n            print(f\"epoch_num {epoch+1}, current loss is: {epoch_loss_mean}\")\n\ndef evaluate(model, device=torch.device('cuda'), covariate_size=3):\n    full_covariate_tensor = torch.tensor(train_covariate_df.to_numpy()).to(device)\n    next_covariate_tensor = torch.tensor(test_covariate_df.to_numpy().reshape(-1, horizon_size * covariate_size)).unsqueeze(dim=0).to(device)#[1,1,horizon_size * covariate_size]\n    results = []\n\n    model.eval()\n    with torch.no_grad():\n        for colname in tqdm(train_target_df.columns):\n            input_target_tensor = torch.tensor(train_target_df[[colname]].to_numpy()).to(device)\n            input_target_covariate_tensor = torch.unsqueeze(torch.cat([input_target_tensor, full_covariate_tensor], dim=1), dim= 0).to(device) #[1, seq_len, 1+covariate_size]\n            input_target_covariate_tensor = input_target_covariate_tensor.permute(1,0,2).to(device) #[seq_len, 1, 1+covariate_size]\n            \n            model_output = model(input_target_covariate_tensor, next_covariate_tensor)\n            \n            output_array = model_output.squeeze(0).cpu()\n            results.append(output_array)\n    predictions = torch.cat(results)\n    diff = torch.sum(torch.mul(predictions[:,:,1].T - torch.tensor(test_target_df.to_numpy()), predictions[:,:,1].T - torch.tensor(test_target_df.to_numpy()))).item()\n    test_rmse = np.sqrt(diff/test_target_df.shape[1])\n    return test_rmse","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:48.686660Z","iopub.execute_input":"2023-01-10T06:54:48.687057Z","iopub.status.idle":"2023-01-10T06:54:48.706036Z","shell.execute_reply.started":"2023-01-10T06:54:48.687019Z","shell.execute_reply":"2023-01-10T06:54:48.705030Z"},"trusted":true},"execution_count":11,"outputs":[]},{"cell_type":"code","source":"model = MQRNN()\ntrain_dataset = MQRNN_dataset(train_target_df, train_covariate_df)\ntrain(model, train_dataset, batch_size=16, num_epochs=15, lr=0.0001)\ntest_rmse = evaluate(model)\nprint(test_rmse)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T06:54:48.707318Z","iopub.execute_input":"2023-01-10T06:54:48.707915Z","iopub.status.idle":"2023-01-10T07:24:45.061124Z","shell.execute_reply.started":"2023-01-10T06:54:48.707874Z","shell.execute_reply":"2023-01-10T07:24:45.059971Z"},"trusted":true},"execution_count":12,"outputs":[{"name":"stderr","text":"Loss:0.46324248885113073: 100%|██████████| 24/24 [01:57<00:00,  4.91s/it]\nLoss:0.4486111273566888: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.9046695750178001: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.5590738157245124: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.2502994214335062: 100%|██████████| 24/24 [01:56<00:00,  4.83s/it]\n","output_type":"stream"},{"name":"stdout","text":"epoch_num 5, current loss is: 9.108657407976985e-07\n","output_type":"stream"},{"name":"stderr","text":"Loss:0.17912541354587783: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.22242166081433362: 100%|██████████| 24/24 [01:56<00:00,  4.85s/it]\nLoss:0.23131378937569877: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.22128374589107713: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.17112013396108156: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\n","output_type":"stream"},{"name":"stdout","text":"epoch_num 10, current loss is: 4.746060925913057e-07\n","output_type":"stream"},{"name":"stderr","text":"Loss:0.23409598323551858: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.18404477580326017: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.22592413670993697: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.17513563014407954: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\nLoss:0.18691363430222796: 100%|██████████| 24/24 [01:56<00:00,  4.84s/it]\n","output_type":"stream"},{"name":"stdout","text":"epoch_num 15, current loss is: 4.6011582813914996e-07\n","output_type":"stream"},{"name":"stderr","text":"100%|██████████| 370/370 [00:46<00:00,  7.97it/s]","output_type":"stream"},{"name":"stdout","text":"0.5596733488977182\n","output_type":"stream"},{"name":"stderr","text":"\n","output_type":"stream"}]}]}