wind_power_forecasting_SVGP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "wind power forecasting-SVGP.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "3DJxaCyC_wnx"
      },
      "source": [
        "!pip install gpytorch\n",
        "!pip install tqdm\n",
        "!pip install -U scikit-learn\n",
        "!pip install properscoring"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "e3MI7b9P_2jM"
      },
      "source": [
        "import math\n",
        "import torch\n",
        "import gpytorch\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "from matplotlib import pyplot as plt\n",
        "import tqdm\n",
        "import sklearn\n",
        "import copy\n",
        "import properscoring as ps\n",
        "%matplotlib inline\n",
        "%load_ext autoreload\n",
        "%autoreload 2\n",
        "from torch.distributions.transformed_distribution import TransformedDistribution\n",
        "from torch.distributions.transforms import SigmoidTransform, PowerTransform\n",
        "from torch.distributions.normal import Normal"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Nhn0p0U7XbrD"
      },
      "source": [
        "def fill_ndarray(t1):\n",
        "    for i in range(t1.shape[1]):  \n",
        "        temp_col = t1[:, i]  \n",
        "        nan_num = np.count_nonzero(temp_col != temp_col)\n",
        "        if nan_num != 0:  \n",
        "            temp_not_nan_col = temp_col[temp_col == temp_col]  \n",
        "            temp_col[np.isnan(temp_col)] = temp_not_nan_col.mean()  \n",
        "    return t1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "moYRZggM0Trg"
      },
      "source": [
        "def get_data(df):\n",
        "  data = df['TARGETVAR'].values.reshape((-1,1))\n",
        "  return data\n",
        "\n",
        "def get_nwp(d):\n",
        "  cls = d.columns\n",
        "  data = []\n",
        "  for i in range(4):\n",
        "    data.append(d[cls[i+3]].values.reshape((-1,1)))\n",
        "  data = np.hstack(data)\n",
        "  return data\n",
        "\n",
        "data = pd.read_csv('Task15_W_Zone1.csv',delimiter=',')\n",
        "\n",
        "power = get_data(data); x = get_nwp(data)\n",
        "y = fill_ndarray(power).flatten()\n",
        "y_true = copy.deepcopy(y)\n",
        "\n",
        "v = 1\n",
        "for i in range(len(y)):\n",
        "    if y[i] == 0:\n",
        "        y[i] = 1e-3\n",
        "        y_true[i] = 1e-3\n",
        "        y[i] = np.log(np.power(y[i],v)/(1-np.power(y[i],v)))\n",
        "    elif y[i] == 1:\n",
        "        y[i] = 1-1e-3\n",
        "        y_true[i] = 1-1e-3\n",
        "        y[i] = np.log(np.power(y[i],v)/(1-np.power(y[i],v)))\n",
        "    else:\n",
        "        y[i] = np.log(np.power(y[i],v)/(1-np.power(y[i],v)))\n",
        "y = torch.tensor(y).float()\n",
        "x = torch.tensor(x).float()\n",
        "\n",
        "samples = math.floor(x.shape[0]*0.8)\n",
        "\n",
        "train_y = y[:samples]\n",
        "train_x = x[:samples,:]\n",
        "\n",
        "\n",
        "y_true = y_true[samples:]\n",
        "test_y = y[samples:]\n",
        "test_x = x[samples:,:]\n",
        "\n",
        "if torch.cuda.is_available():\n",
        "    train_x, train_y, test_x, test_y = train_x.cuda(), train_y.cuda(), test_x.cuda(), test_y.cuda()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hRahv5KkhrV7"
      },
      "source": [
        "from torch.utils.data import TensorDataset, DataLoader\n",
        "train_dataset = TensorDataset(train_x, train_y)\n",
        "train_loader = DataLoader(train_dataset, batch_size=1024, shuffle=True)\n",
        "\n",
        "test_dataset = TensorDataset(test_x, test_y)\n",
        "test_loader = DataLoader(test_dataset, batch_size=1024, shuffle=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7xQQtzFSrG4o"
      },
      "source": [
        "from torch import nn\n",
        "class MLPMean(gpytorch.means.Mean):\n",
        "    def __init__(self, input_size):\n",
        "        super(MLPMean, self).__init__()\n",
        "        self.mlp = nn.Sequential(\n",
        "            nn.Linear(input_size, 512),\n",
        "            nn.ReLU(),  \n",
        "            nn.Linear(512, 512),\n",
        "            nn.ReLU(),\n",
        "            nn.Linear(512, 512),\n",
        "            nn.ReLU(), \n",
        "            nn.Linear(512, 1))\n",
        "\n",
        "        count = 0\n",
        "        for n, p in self.mlp.named_parameters():\n",
        "            self.register_parameter(name = 'mlp' + str(count), parameter = p)\n",
        "            count += 1\n",
        "    \n",
        "    def forward(self, x):\n",
        "        m = self.mlp(x)\n",
        "        return m.squeeze()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H0XpKmVaiGEu"
      },
      "source": [
        "from gpytorch.models import ApproximateGP\n",
        "from gpytorch.variational import CholeskyVariationalDistribution\n",
        "from gpytorch.variational import VariationalStrategy\n",
        "\n",
        "class GPModel(ApproximateGP):\n",
        "    def __init__(self, inducing_points):\n",
        "        variational_distribution = CholeskyVariationalDistribution(inducing_points.size(0))\n",
        "        variational_strategy = VariationalStrategy(self, inducing_points, variational_distribution, learn_inducing_locations=True)\n",
        "        super(GPModel, self).__init__(variational_strategy)\n",
        "        self.scale_to_bounds = gpytorch.utils.grid.ScaleToBounds(-1., 1.)\n",
        "        \n",
        "        self.mean_module = MLPMean(input_size=4)\n",
        "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(ard_num_dims=4))\n",
        "        \n",
        "\n",
        "    def forward(self, x):\n",
        "        x = self.scale_to_bounds(x)\n",
        "        mean_x = self.mean_module(x)\n",
        "        covar_x = self.covar_module(x)\n",
        "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
        "\n",
        "inducing_points = train_x[torch.randperm(train_x.size(0))[:1000]]\n",
        "model = GPModel(inducing_points=inducing_points)\n",
        "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
        "\n",
        "if torch.cuda.is_available():\n",
        "    model = model.cuda()\n",
        "    likelihood = likelihood.cuda()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NJp4xXw_ia4j"
      },
      "source": [
        "num_epochs = 200\n",
        "losses = []\n",
        "\n",
        "model.train()\n",
        "likelihood.train()\n",
        "\n",
        "optimizer = torch.optim.Adam([\n",
        "    {'params': model.parameters()},\n",
        "    {'params': likelihood.parameters()},\n",
        "], lr=0.001)\n",
        "scheduler = torch.optim.lr_scheduler.MultiStepLR(optimizer=optimizer, milestones=[50,100], gamma=0.5)\n",
        "\n",
        "mll = gpytorch.mlls.PredictiveLogLikelihood(likelihood, model, num_data=train_y.size(0))\n",
        "\n",
        "\n",
        "epochs_iter = tqdm.notebook.tqdm(range(num_epochs), desc=\"Epoch\")\n",
        "for i in epochs_iter:\n",
        "    minibatch_iter = tqdm.notebook.tqdm(train_loader, desc=\"Minibatch\", leave=False)\n",
        "    for x_batch, y_batch in minibatch_iter:\n",
        "        optimizer.zero_grad()\n",
        "        output = likelihood(model(x_batch))\n",
        "        loss = -mll(output, y_batch)\n",
        "        minibatch_iter.set_postfix(loss=loss.item())\n",
        "        loss.backward()\n",
        "        optimizer.step()\n",
        "    scheduler.step()\n",
        "    losses.append(loss.item())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hAz-gWm4ipUx"
      },
      "source": [
        "model.eval()\n",
        "likelihood.eval()\n",
        "\n",
        "mean = torch.tensor([0.]); std = torch.tensor([0.]); f_std = torch.tensor([0.])\n",
        "with torch.no_grad():\n",
        "    for x_batch, y_batch in test_loader:\n",
        "        preds = likelihood(model(x_batch))\n",
        "        mean = torch.cat([mean, preds.mean.cpu()])\n",
        "        std = torch.cat([std, preds.stddev.cpu()])\n",
        "        f_std = torch.cat([f_std, model(x_batch).stddev.cpu()])\n",
        "mean = mean[1:]; std = std[1:]; f_std = f_std[1:]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w9YtjqCsCJQJ"
      },
      "source": [
        "mean = mean.cpu()\n",
        "std = std.cpu()\n",
        "f_std = f_std.cpu()\n",
        "\n",
        "train_x = train_x.cpu()\n",
        "train_y = train_y.cpu()\n",
        "test_x = test_x.cpu()\n",
        "test_y = test_y.cpu()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-oWVymEQkD_X",
        "outputId": "a45cdfb1-a058-4340-ebc8-b4cd40d6ead2"
      },
      "source": [
        "length = len(y_true)\n",
        "score = 0 \n",
        "for i in range(length):\n",
        "  m = torch.distributions.normal.Normal(mean[i],std[i])\n",
        "  transforms = [SigmoidTransform(),PowerTransform(1)]\n",
        "  logit_normal = TransformedDistribution(m, transforms)\n",
        "  sample = logit_normal.sample([1000]).numpy()\n",
        "  score += ps.crps_ensemble(y_true[i] , sample)\n",
        "\n",
        "score = score/length\n",
        "print(score)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.10304097427849318\n"
          ]
        }
      ]
    }
  ]
}