{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "16537827-7386-4163-b95f-2997fc020a2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from shutil import copyfile\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a517af1c-4204-45c9-aae4-865a2cb259e9",
   "metadata": {},
   "source": [
    "Data manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "62628e60-28c6-4a9a-8a81-22e5bfd74722",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0fa93674-d4e7-4b36-ab3a-ebb11df12ed3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler, RobustScaler\n",
    "from sklearn.exceptions import NotFittedError"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acb33a41-06b9-4a1d-9ea7-6a2d87b1f4fb",
   "metadata": {},
   "source": [
    "Plotting / Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a42ae056-7511-4e17-b4ba-981fbcfaf922",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f8b23d5-e253-408b-907f-6f9990a98a96",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.rcParams[\"figure.figsize\"] = (15, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90fdac33-eed4-4ab4-b2b1-de0f1f27727b",
   "metadata": {},
   "source": [
    "Gaussian Process Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e629471f-350e-4af3-83df-377794a20a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gpflow\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b33dbca7-4419-4201-8d49-9fe3a5e16a33",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gpflow.utilities import print_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fe8d81a8-f8ec-41d0-8a6c-77ee73619def",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpflow.config.set_default_summary_fmt(\"notebook\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aba0df5-b0e3-4738-bb61-1dad869d1ea3",
   "metadata": {},
   "source": [
    "## Load previously exported data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a8bf0b3f-1236-41c5-ba72-7e274a75d22f",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_exps = ['Exp1', 'Exp3', 'Exp5', 'Exp6']\n",
    "test_exps = ['Exp2', 'Exp4', 'Exp7']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3c8ad21b-8566-4d14-a13c-99e2dc3efc74",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_cols = ['time_h', 'time_m']\n",
    "w_cols = ['SolRad', 'OutsideTemp']\n",
    "u_cols = ['SimulatedHeat']\n",
    "y_cols = ['SimulatedTemp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ad00a16f-7cb6-44a6-90e7-ecc07003ffad",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_lags = 4\n",
    "w_lags = 1\n",
    "u_lags = 3\n",
    "y_lags = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e3225bc5-70b1-4b27-9893-8d0aae750bc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dict_cols = pickle.load(open(Path(\"dict_cols.pkl\"), 'rb'))\n",
    "#dict_cols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "815913ad-73b3-407c-8e05-e3acfa65d73c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_cols = {\n",
    "    't': (t_lags, t_cols),\n",
    "    'w': (w_lags, w_cols),\n",
    "    'u': (u_lags, u_cols),\n",
    "    'y': (y_lags, y_cols)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b1574455-ece7-4a45-92a9-70410f73266e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs_train = pickle.load(open(Path(\"dfs_train.pkl\"), 'rb'))\n",
    "dfs_test = pickle.load(open(Path(\"dfs_test.pkl\"), 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0ebe637c-ad84-4393-9199-6ba088462132",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = pickle.load(open(Path(\"scaler.pkl\"), 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "136ba207-b4a7-427d-bdb2-0321110be9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_scaled_df(df, dict_cols, scaler):\n",
    "    \n",
    "    t_list = dict_cols['t'][1]\n",
    "    w_list = dict_cols['w'][1]\n",
    "    u_list = dict_cols['u'][1]\n",
    "    y_list = dict_cols['y'][1]\n",
    "    \n",
    "    df_local = df[t_list + w_list + u_list + y_list]\n",
    "    df_scaled = df_local.to_numpy()\n",
    "    \n",
    "    try:\n",
    "        df_scaled = scaler.transform(df_scaled)\n",
    "    except NotFittedError:\n",
    "        df_scaled = scaler.fit_transform(df_scaled)\n",
    "        \n",
    "    df_scaled = pd.DataFrame(df_scaled, index = df_local.index, columns = df_local.columns)\n",
    "    \n",
    "    return df_scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d5130037-26f7-4f47-8382-81b12a131e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = pd.concat(dfs_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "a1d9c9f8-6689-40b4-8ae7-75395b343de2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35185.23586737608"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(df_train.to_numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b862da6c-9bca-4647-aeec-ee861c11293c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_sc = get_scaled_df(df_train, dict_cols, scaler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c0c18759-db54-4c46-8875-535b28d9ede6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.478002157732377"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(df_train_sc.to_numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5b39e5c3-94ca-4765-aa6b-1426b226ca27",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs_train_sc = []\n",
    "dfs_test_sc = []\n",
    "for df in dfs_train:\n",
    "    df_sc = get_scaled_df(df, dict_cols, scaler)\n",
    "    dfs_train_sc.append(df_sc)\n",
    "    \n",
    "for df in dfs_test:\n",
    "    df_sc = get_scaled_df(df, dict_cols, scaler)\n",
    "    dfs_test_sc.append(df_sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e8c97d70-2c97-4db8-9e2e-636086bd6271",
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_to_gpr(df, dict_cols):\n",
    "    \n",
    "    t_list = dict_cols['t'][1]\n",
    "    w_list = dict_cols['w'][1]\n",
    "    u_list = dict_cols['u'][1]\n",
    "    y_list = dict_cols['y'][1]\n",
    "    \n",
    "    df_gpr = df[t_list + w_list + u_list + y_list].copy()\n",
    "    \n",
    "    for lags, names in dict_cols.values():\n",
    "        for name in names:\n",
    "            col_idx = df_gpr.columns.get_loc(name)\n",
    "            for lag in range(1, lags + 1):\n",
    "                df_gpr.insert(col_idx + lag, f\"{name}_{lag}\", df_gpr.loc[:, name].shift(lag))\n",
    "\n",
    "    df_gpr.dropna(inplace = True)\n",
    "    \n",
    "    return df_gpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b351b32a-7018-4c90-8aa0-7c0be5461c2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dfs_gpr_train = pickle.load(open(Path(\"dfs_gpr_train.pkl\"), 'rb'))\n",
    "#dfs_gpr_test = pickle.load(open(Path(\"dfs_gpr_test.pkl\"), 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b0453fdf-58f3-49e6-bdea-ebe5940422e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time_h</th>\n",
       "      <th>time_h_1</th>\n",
       "      <th>time_h_2</th>\n",
       "      <th>time_h_3</th>\n",
       "      <th>time_h_4</th>\n",
       "      <th>time_m</th>\n",
       "      <th>time_m_1</th>\n",
       "      <th>time_m_2</th>\n",
       "      <th>time_m_3</th>\n",
       "      <th>time_m_4</th>\n",
       "      <th>...</th>\n",
       "      <th>OutsideTemp</th>\n",
       "      <th>OutsideTemp_1</th>\n",
       "      <th>SimulatedHeat</th>\n",
       "      <th>SimulatedHeat_1</th>\n",
       "      <th>SimulatedHeat_2</th>\n",
       "      <th>SimulatedHeat_3</th>\n",
       "      <th>SimulatedTemp</th>\n",
       "      <th>SimulatedTemp_1</th>\n",
       "      <th>SimulatedTemp_2</th>\n",
       "      <th>SimulatedTemp_3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timestamp</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2017-06-01 20:20:00+02:00</th>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>-0.454545</td>\n",
       "      <td>-0.636364</td>\n",
       "      <td>-0.818182</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.179560</td>\n",
       "      <td>-0.132679</td>\n",
       "      <td>-0.094006</td>\n",
       "      <td>-0.076890</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-06-01 20:25:00+02:00</th>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>-0.090909</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>-0.454545</td>\n",
       "      <td>-0.636364</td>\n",
       "      <td>-0.818182</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.208254</td>\n",
       "      <td>-0.179560</td>\n",
       "      <td>-0.132679</td>\n",
       "      <td>-0.094006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-06-01 20:30:00+02:00</th>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.090909</td>\n",
       "      <td>-0.090909</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>-0.454545</td>\n",
       "      <td>-0.636364</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.580115</td>\n",
       "      <td>-0.222268</td>\n",
       "      <td>-0.208254</td>\n",
       "      <td>-0.179560</td>\n",
       "      <td>-0.132679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-06-01 20:35:00+02:00</th>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.272727</td>\n",
       "      <td>0.090909</td>\n",
       "      <td>-0.090909</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>-0.454545</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.234855</td>\n",
       "      <td>-0.222268</td>\n",
       "      <td>-0.208254</td>\n",
       "      <td>-0.179560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-06-01 20:40:00+02:00</th>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.73913</td>\n",
       "      <td>0.454545</td>\n",
       "      <td>0.272727</td>\n",
       "      <td>0.090909</td>\n",
       "      <td>-0.090909</td>\n",
       "      <td>-0.272727</td>\n",
       "      <td>...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.904139</td>\n",
       "      <td>-0.247166</td>\n",
       "      <td>-0.234855</td>\n",
       "      <td>-0.222268</td>\n",
       "      <td>-0.208254</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 22 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                            time_h  time_h_1  time_h_2  time_h_3  time_h_4  \\\n",
       "timestamp                                                                    \n",
       "2017-06-01 20:20:00+02:00  0.73913   0.73913   0.73913   0.73913   0.73913   \n",
       "2017-06-01 20:25:00+02:00  0.73913   0.73913   0.73913   0.73913   0.73913   \n",
       "2017-06-01 20:30:00+02:00  0.73913   0.73913   0.73913   0.73913   0.73913   \n",
       "2017-06-01 20:35:00+02:00  0.73913   0.73913   0.73913   0.73913   0.73913   \n",
       "2017-06-01 20:40:00+02:00  0.73913   0.73913   0.73913   0.73913   0.73913   \n",
       "\n",
       "                             time_m  time_m_1  time_m_2  time_m_3  time_m_4  \\\n",
       "timestamp                                                                     \n",
       "2017-06-01 20:20:00+02:00 -0.272727 -0.454545 -0.636364 -0.818182 -1.000000   \n",
       "2017-06-01 20:25:00+02:00 -0.090909 -0.272727 -0.454545 -0.636364 -0.818182   \n",
       "2017-06-01 20:30:00+02:00  0.090909 -0.090909 -0.272727 -0.454545 -0.636364   \n",
       "2017-06-01 20:35:00+02:00  0.272727  0.090909 -0.090909 -0.272727 -0.454545   \n",
       "2017-06-01 20:40:00+02:00  0.454545  0.272727  0.090909 -0.090909 -0.272727   \n",
       "\n",
       "                           ...  OutsideTemp  OutsideTemp_1  SimulatedHeat  \\\n",
       "timestamp                  ...                                              \n",
       "2017-06-01 20:20:00+02:00  ...     0.058824       0.058824      -0.904139   \n",
       "2017-06-01 20:25:00+02:00  ...     0.058824       0.058824      -0.904139   \n",
       "2017-06-01 20:30:00+02:00  ...     0.058824       0.058824      -0.904139   \n",
       "2017-06-01 20:35:00+02:00  ...     0.058824       0.058824      -0.904139   \n",
       "2017-06-01 20:40:00+02:00  ...     0.058824       0.058824      -0.904139   \n",
       "\n",
       "                           SimulatedHeat_1  SimulatedHeat_2  SimulatedHeat_3  \\\n",
       "timestamp                                                                      \n",
       "2017-06-01 20:20:00+02:00        -0.580115        -0.580115        -0.580115   \n",
       "2017-06-01 20:25:00+02:00        -0.904139        -0.580115        -0.580115   \n",
       "2017-06-01 20:30:00+02:00        -0.904139        -0.904139        -0.580115   \n",
       "2017-06-01 20:35:00+02:00        -0.904139        -0.904139        -0.904139   \n",
       "2017-06-01 20:40:00+02:00        -0.904139        -0.904139        -0.904139   \n",
       "\n",
       "                           SimulatedTemp  SimulatedTemp_1  SimulatedTemp_2  \\\n",
       "timestamp                                                                    \n",
       "2017-06-01 20:20:00+02:00      -0.179560        -0.132679        -0.094006   \n",
       "2017-06-01 20:25:00+02:00      -0.208254        -0.179560        -0.132679   \n",
       "2017-06-01 20:30:00+02:00      -0.222268        -0.208254        -0.179560   \n",
       "2017-06-01 20:35:00+02:00      -0.234855        -0.222268        -0.208254   \n",
       "2017-06-01 20:40:00+02:00      -0.247166        -0.234855        -0.222268   \n",
       "\n",
       "                           SimulatedTemp_3  \n",
       "timestamp                                   \n",
       "2017-06-01 20:20:00+02:00        -0.076890  \n",
       "2017-06-01 20:25:00+02:00        -0.094006  \n",
       "2017-06-01 20:30:00+02:00        -0.132679  \n",
       "2017-06-01 20:35:00+02:00        -0.179560  \n",
       "2017-06-01 20:40:00+02:00        -0.208254  \n",
       "\n",
       "[5 rows x 22 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfs_gpr_train = []\n",
    "for df_sc in dfs_train_sc:\n",
    "    dfs_gpr_train.append(data_to_gpr(df_sc, dict_cols))\n",
    "df_gpr_train = pd.concat(dfs_gpr_train)\n",
    "df_gpr_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8b866f02-084a-4ec1-b269-818949461530",
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_gpr_train = pd.concat(dfs_gpr_train)\n",
    "\n",
    "df_input_train = df_gpr_train.drop(columns = dict_cols['u'][1] + dict_cols['y'][1])\n",
    "df_output_train = df_gpr_train[dict_cols['y'][1]]\n",
    "\n",
    "np_input_train = df_input_train.to_numpy()\n",
    "np_output_train = df_output_train.to_numpy().reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "288c99da-56c2-4216-819e-b3722f2a502c",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Define Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "dadb3f43-af78-4cb6-98c6-408559f62479",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_dims = np_input_train.shape[1]\n",
    "rational_dims = np.arange(0, (dict_cols['t'][0] + 1) * len(dict_cols['t'][1]), 1)\n",
    "nb_rational_dims = len(rational_dims)\n",
    "squared_dims = np.arange(nb_rational_dims, nb_dims, 1)\n",
    "nb_squared_dims = len(squared_dims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "9bca461d-b4f9-4f9f-a6d6-65f7601fa020",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rational dims: 10\n",
      "squared dims: 10\n"
     ]
    }
   ],
   "source": [
    "print(f\"rational dims: {nb_rational_dims}\")\n",
    "print(f\"squared dims: {nb_squared_dims}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "dc8ef81f-271c-4e7a-8eb7-8981557cf7bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "squared_l = [1e-4] * nb_squared_dims\n",
    "rational_l = [1e-7] * nb_rational_dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "da302d5c-cb00-47e5-9cad-6e4de21eceb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "squared_l = np.linspace(0.01, 1, nb_squared_dims)\n",
    "rational_l = np.linspace(0.01, 1, nb_rational_dims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "61caaaba-7527-4437-a7bf-62b5594c147c",
   "metadata": {},
   "outputs": [],
   "source": [
    "k0 = gpflow.kernels.SquaredExponential(lengthscales = squared_l, active_dims = squared_dims, variance = 2)\n",
    "k1 = gpflow.kernels.Constant()\n",
    "k2 = gpflow.kernels.RationalQuadratic(lengthscales = rational_l, active_dims = rational_dims, variance = 2)\n",
    "k3 = gpflow.kernels.Periodic(k2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "87e3c3a2-4e67-4181-bd2d-22d684b7ce03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>name                                      </th><th>class    </th><th>transform  </th><th>prior  </th><th>trainable  </th><th>shape  </th><th>dtype  </th><th>value               </th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>Product.kernels[0].kernels[0].variance    </td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>2.0                 </td></tr>\n",
       "<tr><td>Product.kernels[0].kernels[0].lengthscales</td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>(10,)  </td><td>float64</td><td>[0.01, 0.12, 0.23...</td></tr>\n",
       "<tr><td>Product.kernels[0].kernels[1].variance    </td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>1.0                 </td></tr>\n",
       "<tr><td>Product.kernels[1].variance               </td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>2.0                 </td></tr>\n",
       "<tr><td>Product.kernels[1].lengthscales           </td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>(10,)  </td><td>float64</td><td>[0.01, 0.12, 0.23...</td></tr>\n",
       "<tr><td>Product.kernels[1].alpha                  </td><td>Parameter</td><td>Softplus   </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>1.0                 </td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k = (k0 + k1) * k2\n",
    "print_summary(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4af25a43-15c9-4543-af73-3c313b5fc7af",
   "metadata": {},
   "source": [
    "## Compile Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "99da48d2-f04e-4ef8-a248-bc9b82fdbabd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>name                                         </th><th>class    </th><th>transform       </th><th>prior  </th><th>trainable  </th><th>shape  </th><th>dtype  </th><th>value               </th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>GPR.kernel.kernels[0].kernels[0].variance    </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>2.0                 </td></tr>\n",
       "<tr><td>GPR.kernel.kernels[0].kernels[0].lengthscales</td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>(10,)  </td><td>float64</td><td>[0.01, 0.12, 0.23...</td></tr>\n",
       "<tr><td>GPR.kernel.kernels[0].kernels[1].variance    </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>1.0                 </td></tr>\n",
       "<tr><td>GPR.kernel.kernels[1].variance               </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>2.0                 </td></tr>\n",
       "<tr><td>GPR.kernel.kernels[1].lengthscales           </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>(10,)  </td><td>float64</td><td>[0.01, 0.12, 0.23...</td></tr>\n",
       "<tr><td>GPR.kernel.kernels[1].alpha                  </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>1.0                 </td></tr>\n",
       "<tr><td>GPR.likelihood.variance                      </td><td>Parameter</td><td>Softplus + Shift</td><td>       </td><td>True       </td><td>()     </td><td>float64</td><td>1.0                 </td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = gpflow.models.GPR(\n",
    "    data = (np_input_train, np_output_train), \n",
    "    kernel = k, \n",
    "    mean_function = None\n",
    "    )\n",
    "print_summary(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08f41235-12df-4e9c-bf63-e7a4390cf21a",
   "metadata": {},
   "source": [
    "## Train Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "9e5e2138-b342-4d44-8987-b7758b0daa6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = gpflow.optimizers.Scipy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "fcfb78e0-2a3b-4a16-a621-6698abdaf3ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "c449b100-acff-414d-a648-c4233879c253",
   "metadata": {},
   "outputs": [
    {
     "ename": "InvalidArgumentError",
     "evalue": " Input matrix is not invertible.\n\t [[node gradient_tape/triangular_solve/MatrixTriangularSolve (defined at /usr/lib/python3.9/site-packages/gpflow/optimizers/scipy.py:173) ]] [Op:__inference__tf_eval_1118]\n\nErrors may have originated from an input operation.\nInput Source operations connected to node gradient_tape/triangular_solve/MatrixTriangularSolve:\n Cholesky (defined at /usr/lib/python3.9/site-packages/gpflow/models/gpr.py:87)\n\nFunction call stack:\n_tf_eval\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mInvalidArgumentError\u001b[0m                      Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-36-5f570574e163>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mstart_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdatetime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mopt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtraining_loss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrainable_variables\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"Finished fitting in {datetime.now() - start_time}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mprint_summary\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/gpflow/optimizers/scipy.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(self, closure, variables, method, step_callback, compile, **scipy_kwargs)\u001b[0m\n\u001b[1;32m     88\u001b[0m             \u001b[0mscipy_kwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcallback\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m         return scipy.optimize.minimize(\n\u001b[0m\u001b[1;32m     91\u001b[0m             \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_params\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjac\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mscipy_kwargs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     92\u001b[0m         )\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    617\u001b[0m                                   **options)\n\u001b[1;32m    618\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'l-bfgs-b'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 619\u001b[0;31m         return _minimize_lbfgsb(fun, x0, args, jac, bounds,\n\u001b[0m\u001b[1;32m    620\u001b[0m                                 callback=callback, **options)\n\u001b[1;32m    621\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'tnc'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36m_minimize_lbfgsb\u001b[0;34m(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)\u001b[0m\n\u001b[1;32m    358\u001b[0m             \u001b[0;31m# until the completion of the current minimization iteration.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    359\u001b[0m             \u001b[0;31m# Overwrite f and g:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 360\u001b[0;31m             \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    361\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mtask_str\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstartswith\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mb'NEW_X'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    362\u001b[0m             \u001b[0;31m# new iteration\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py\u001b[0m in \u001b[0;36mfun_and_grad\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    258\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    259\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_x_impl\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 260\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    261\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    262\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py\u001b[0m in \u001b[0;36m_update_fun\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    224\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_update_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    225\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_updated\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 226\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_fun_impl\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    227\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_updated\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    228\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py\u001b[0m in \u001b[0;36mupdate_fun\u001b[0;34m()\u001b[0m\n\u001b[1;32m    131\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    132\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0mupdate_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 133\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    134\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    135\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_fun_impl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mupdate_fun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py\u001b[0m in \u001b[0;36mfun_wrapped\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    128\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnfev\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 130\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    131\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    132\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0mupdate_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x, *args)\u001b[0m\n\u001b[1;32m     72\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     73\u001b[0m         \u001b[0;34m\"\"\" returns the the function value \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 74\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compute_if_needed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     75\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_value\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     76\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m_compute_if_needed\u001b[0;34m(self, x, *args)\u001b[0m\n\u001b[1;32m     66\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_value\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjac\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 68\u001b[0;31m             \u001b[0mfg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjac\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_value\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/gpflow/optimizers/scipy.py\u001b[0m in \u001b[0;36m_eval\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    111\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    112\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0m_eval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0mTuple\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 113\u001b[0;31m             \u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_tf_eval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconvert_to_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    114\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mloss\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    115\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m    826\u001b[0m     \u001b[0mtracing_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexperimental_get_tracing_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    827\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_name\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mtm\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 828\u001b[0;31m       \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    829\u001b[0m       \u001b[0mcompiler\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"xla\"\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_experimental_compile\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m\"nonXla\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    830\u001b[0m       \u001b[0mnew_tracing_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexperimental_get_tracing_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m_call\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m    860\u001b[0m       \u001b[0;31m# In this case we have not created variables on the first call. So we can\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    861\u001b[0m       \u001b[0;31m# run the first trace but we should fail if variables are created.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 862\u001b[0;31m       \u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stateful_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    863\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_created_variables\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    864\u001b[0m         raise ValueError(\"Creating variables on a non-first call to a function\"\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2940\u001b[0m       (graph_function,\n\u001b[1;32m   2941\u001b[0m        filtered_flat_args) = self._maybe_define_function(args, kwargs)\n\u001b[0;32m-> 2942\u001b[0;31m     return graph_function._call_flat(\n\u001b[0m\u001b[1;32m   2943\u001b[0m         filtered_flat_args, captured_inputs=graph_function.captured_inputs)  # pylint: disable=protected-access\n\u001b[1;32m   2944\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_call_flat\u001b[0;34m(self, args, captured_inputs, cancellation_manager)\u001b[0m\n\u001b[1;32m   1916\u001b[0m         and executing_eagerly):\n\u001b[1;32m   1917\u001b[0m       \u001b[0;31m# No tape is watching; skip to running the function.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1918\u001b[0;31m       return self._build_call_outputs(self._inference_function.call(\n\u001b[0m\u001b[1;32m   1919\u001b[0m           ctx, args, cancellation_manager=cancellation_manager))\n\u001b[1;32m   1920\u001b[0m     forward_backward = self._select_forward_and_backward_functions(\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36mcall\u001b[0;34m(self, ctx, args, cancellation_manager)\u001b[0m\n\u001b[1;32m    553\u001b[0m       \u001b[0;32mwith\u001b[0m \u001b[0m_InterpolateFunctionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    554\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcancellation_manager\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 555\u001b[0;31m           outputs = execute.execute(\n\u001b[0m\u001b[1;32m    556\u001b[0m               \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    557\u001b[0m               \u001b[0mnum_outputs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_num_outputs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.9/site-packages/tensorflow/python/eager/execute.py\u001b[0m in \u001b[0;36mquick_execute\u001b[0;34m(op_name, num_outputs, inputs, attrs, ctx, name)\u001b[0m\n\u001b[1;32m     57\u001b[0m   \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m     \u001b[0mctx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mensure_initialized\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 59\u001b[0;31m     tensors = pywrap_tfe.TFE_Py_Execute(ctx._handle, device_name, op_name,\n\u001b[0m\u001b[1;32m     60\u001b[0m                                         inputs, attrs, num_outputs)\n\u001b[1;32m     61\u001b[0m   \u001b[0;32mexcept\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_NotOkStatusException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mInvalidArgumentError\u001b[0m:  Input matrix is not invertible.\n\t [[node gradient_tape/triangular_solve/MatrixTriangularSolve (defined at /usr/lib/python3.9/site-packages/gpflow/optimizers/scipy.py:173) ]] [Op:__inference__tf_eval_1118]\n\nErrors may have originated from an input operation.\nInput Source operations connected to node gradient_tape/triangular_solve/MatrixTriangularSolve:\n Cholesky (defined at /usr/lib/python3.9/site-packages/gpflow/models/gpr.py:87)\n\nFunction call stack:\n_tf_eval\n"
     ]
    }
   ],
   "source": [
    "start_time = datetime.now()\n",
    "opt.minimize(m.training_loss, m.trainable_variables)\n",
    "print(f\"Finished fitting in {datetime.now() - start_time}\")\n",
    "print_summary(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd49280-bb3f-4903-a339-b7225a56ae16",
   "metadata": {},
   "source": [
    "## Evaluate performance on training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bf833f2-83d7-4af4-8d4c-d3a663c18b53",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_plts = len(dfs_gpr_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "063ef841-2708-421c-9793-f878ac8a6e1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (20, 20))\n",
    "\n",
    "for idx, df_iter in enumerate(dfs_gpr_train):\n",
    "    plt.subplot(nb_plts, 1, idx + 1)\n",
    "    df_input_iter = df_iter.drop(columns = dict_cols['y'][1] + dict_cols['u'][1])\n",
    "    df_output_iter = df_iter[dict_cols['y'][1]]\n",
    "    np_input_iter = df_input_iter.to_numpy()\n",
    "    np_output_iter = df_output_iter.to_numpy().reshape(-1, 1)\n",
    "    \n",
    "    mean, var = m.predict_f(np_input_iter)\n",
    "    \n",
    "    plt.plot(df_iter.index, np_output_iter[:, :], label = 'Measured data')\n",
    "    plt.plot(df_iter.index, mean[:, :], label = 'Gaussian Process Prediction')\n",
    "    plt.fill_between(\n",
    "        df_iter.index, \n",
    "        mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),\n",
    "        mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),\n",
    "        alpha = 0.2\n",
    "    )\n",
    "    plt.title(f\"Model Performance on training data: {train_exps[idx]}\")\n",
    "    plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d7d8ca2-1c2d-42dc-a1d2-1f4af11c9d19",
   "metadata": {},
   "source": [
    "## Evaluate performance on test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc73ac1e-024b-4c62-9f2a-4a61b776ecb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs_gpr_test = []\n",
    "for df_sc in dfs_test_sc:\n",
    "    dfs_gpr_test.append(data_to_gpr(df_sc, dict_cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "095c0b27-3faa-4225-a91c-3d18f2a033f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (20, 20))\n",
    "\n",
    "for idx, df_iter in enumerate(dfs_gpr_test):\n",
    "    plt.subplot(nb_plts, 1, idx + 1)\n",
    "    df_input_iter = df_iter.drop(columns = dict_cols['y'][1] + dict_cols['u'][1])\n",
    "    df_output_iter = df_iter[dict_cols['y'][1]]\n",
    "    np_input_iter = df_input_iter.to_numpy()\n",
    "    np_output_iter = df_output_iter.to_numpy().reshape(-1, 1)\n",
    "    \n",
    "    mean, var = m.predict_f(np_input_iter)\n",
    "    \n",
    "    plt.plot(df_iter.index, np_output_iter[:, :], label = 'Measured data')\n",
    "    plt.plot(df_iter.index, mean[:, :], label = 'Gaussian Process Prediction')\n",
    "    plt.fill_between(\n",
    "        df_iter.index, \n",
    "        mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),\n",
    "        mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),\n",
    "        alpha = 0.2\n",
    "    )\n",
    "    plt.title(f\"Model Performance on test data: {test_exps[idx]}\")\n",
    "    plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}