GP_model_performance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "546c8954-eee1-4da9-932a-f5432143afad",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-16 13:30:48.452734: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "2024-06-16 13:30:49.485078: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /share/apps/22.04/openmpi/4.1.5/lib:/share/apps/22.04/slurm/23.02.7/lib\n",
      "2024-06-16 13:30:49.485144: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n",
      "2024-06-16 13:30:49.651620: E tensorflow/stream_executor/cuda/cuda_blas.cc:2981] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered\n",
      "2024-06-16 13:30:53.115710: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer.so.7'; dlerror: libnvinfer.so.7: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /share/apps/22.04/openmpi/4.1.5/lib:/share/apps/22.04/slurm/23.02.7/lib\n",
      "2024-06-16 13:30:53.117201: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer_plugin.so.7'; dlerror: libnvinfer_plugin.so.7: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /share/apps/22.04/openmpi/4.1.5/lib:/share/apps/22.04/slurm/23.02.7/lib\n",
      "2024-06-16 13:30:53.117215: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Cannot dlopen some TensorRT libraries. If you would like to use Nvidia GPU with TensorRT, please make sure the missing libraries mentioned above are installed properly.\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import gpflow\n",
    "import pandas as pd\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import time\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "import statsmodels.api as sm\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "import statsmodels.api as sm\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.svm import SVR\n",
    "from sklearn import preprocessing\n",
    "from xgboost import XGBRegressor\n",
    "import joblib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e051a2ba-77da-4f5c-8513-c1c3da59c992",
   "metadata": {},
   "outputs": [],
   "source": [
    "pfts = [\"Tropical\", \"Temperate_Broadleaf\",\"Temperate_Deciduous\"\n",
    "        ,\"Temperate_Coniferous\",\"Boreal\"]\n",
    "\n",
    "var_list = ['GPP', 'NPP', 'VEGC', 'SOILORGC']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9a493b28-bfe3-4e26-9710-107c00fd874c",
   "metadata": {},
   "outputs": [],
   "source": [
    "pft = pfts[0]\n",
    "all_model_data = pd.read_csv('/home/smmrrr/TEM_Analysis/TEM_Analysis/model_data/merged_outliers_removed/'+pft+'_cru_hist_canesm5_ssp245.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "588f7c6a-e2ba-4da3-b9dc-2966315a8429",
   "metadata": {},
   "outputs": [],
   "source": [
    "### get 15% hold out\n",
    "pred_vars = ['ordinal_stand_age'\n",
    "                         ,'nirr_avg', 'nirr_min', 'nirr_max', 'prec_avg',\n",
    "       'prec_min', 'prec_max', 'tair_avg', 'tair_min', 'tair_max',\n",
    "       'trange_avg', 'trange_min', 'trange_max', 'vpr_avg',\n",
    "       'vpr_min', 'vpr_max', 'wind_avg', 'wind_min', 'wind_max','s1', 's2', 'co2']\n",
    "# X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.15, random_state=42)\n",
    "# Get the indices and split them\n",
    "X = np.array(all_model_data[pred_vars])\n",
    "Y = np.array(all_model_data[var])\n",
    "\n",
    "# Split the data into training and testing sets (80% train, 20% test)\n",
    "# randomness can be controlled with a defined `random_state`\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.99, random_state=42)\n",
    "# Y = np.array(\n",
    "#     model_data.loc[:, 'monthly_mean_ratio']\n",
    "# )\n",
    "# X = np.array(\n",
    "#      # model_data.loc[model_data['year']>=2020,['ordinal_stand_age']]\n",
    "#      model_data.loc[:,['ordinal_stand_age','tmp_ratio']]\n",
    "# )\n",
    "Y = Y.reshape(len(Y),1)\n",
    "Y_train = Y_train.reshape(len(Y_train),1)\n",
    "Y_test = Y_test.reshape(len(Y_test),1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "7f94ea16-05a7-49be-86b7-d091c850efcf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13531265   136679\n"
     ]
    }
   ],
   "source": [
    "print(len(X_test), ' ',len(X_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "232db33c-a170-4005-954c-f97be181e8f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03092781016851379"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "410038/13257906"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cf15921d-a223-4609-8134-fde44bf296f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "var = var_list[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "e4ac0d6e-30c7-443e-a856-49ae804c296f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_train = Y_train.reshape(len(Y_train),1)\n",
    "Y_test = Y_test.reshape(len(Y_test),1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "aaeeff3c-082f-46d3-bdac-c1130c96f494",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.     36.678 112.805 158.25  193.93  226.752 324.7  ]\n"
     ]
    }
   ],
   "source": [
    "percentiles = [0, 10, 25, 50, 75, 90, 100]\n",
    "percentile_values = np.percentile(Y_train, percentiles)\n",
    "print(percentile_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "cb958130-d733-4f69-9722-0eacb4cd98ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.     36.678 112.805 158.25  193.93  226.752 324.7  ]\n",
      "[     4     10     12 ... 136657 136673 136676]\n",
      "[     0      5     15 ... 136641 136662 136663]\n",
      "[     1      6      8 ... 136659 136666 136670]\n",
      "[     7     13     14 ... 136675 136677 136678]\n",
      "[    18     19     53 ... 136651 136671 136672]\n",
      "[     2      3     16 ... 136655 136658 136668]\n",
      "Shape of subset_indices: (240, 1)\n",
      "Shape of inducing_variable: (240, 22)\n"
     ]
    }
   ],
   "source": [
    "rng = np.random.default_rng(12345)\n",
    "bin_size = 40\n",
    "\n",
    "# Calculate the percentile values on Y\n",
    "percentiles = [0, 10, 25, 50, 75, 90, 100]\n",
    "percentile_values = np.percentile(Y_train, percentiles)\n",
    "print(percentile_values)\n",
    "\n",
    "# Initialize a list to store the indices of the subsets\n",
    "subset_indices = []\n",
    "\n",
    "# For each pair of percentile values, select the indices of `bin_size` points from within that range\n",
    "for i in range(len(percentile_values) - 1):\n",
    "    # Get the value range for the current bin\n",
    "    low, high = percentile_values[i], percentile_values[i+1]\n",
    "\n",
    "    # Find the indices where Y falls within the current percentile range\n",
    "    in_range_indices = np.where((Y_train >= low) & (Y_train < high))[0]\n",
    "    print(in_range_indices)\n",
    "\n",
    "    # If there are fewer than bin_size indices, take all of them\n",
    "    if len(in_range_indices) < bin_size:\n",
    "        selected_indices = in_range_indices\n",
    "    else:\n",
    "        # Otherwise, randomly select `bin_size` indices within this range\n",
    "        selected_indices = rng.choice(in_range_indices, size=bin_size, replace=False)\n",
    "\n",
    "    # Extend the subset_indices list with the selected indices\n",
    "    subset_indices.extend(selected_indices)\n",
    "\n",
    "# Convert the list of indices to a numpy array and reshape if necessary\n",
    "subset_indices = np.array(subset_indices).reshape(-1, 1)\n",
    "print(\"Shape of subset_indices:\", subset_indices.shape)\n",
    "\n",
    "# Using the selected indices, we can now reference the corresponding rows from X\n",
    "inducing_variable = X_train[subset_indices.flatten(), :]\n",
    "print(\"Shape of inducing_variable:\", inducing_variable.shape)\n",
    "inducing_variable = X_train[subset_indices.flatten(),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "3871de0d-916f-47d2-9bdf-fb7f963d4e53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.  , 229.16, 184.32, ...,  33.94,  24.03, 301.76],\n",
       "       [  3.  , 265.77, 210.31, ...,  67.62,  19.16, 575.47],\n",
       "       [  5.  , 256.41, 221.03, ...,  49.04,  18.09, 299.84],\n",
       "       ...,\n",
       "       [ 22.  , 196.81, 143.49, ...,  39.26,  20.28, 579.54],\n",
       "       [  6.  , 264.15, 167.78, ...,  44.63,  30.32, 588.38],\n",
       "       [  5.  , 212.68, 173.08, ...,  45.69,  21.81, 406.94]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inducing_variable.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "4adc96b7-7d2f-46c7-8be0-2008888e78bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 40996)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(np.where((Y_train >= low) & (Y_train < high))).shape \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "96d1d84f-473d-4af9-8094-a192b0e64896",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(410038, 1)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "eb269308-416a-465e-85a7-595815c99fb6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([], dtype=float64)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subset_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "960e43a5-3ac1-48db-9586-54a8097bb4b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(410038, 22)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "d0f7bd39-e6d4-48bb-ab70-81147daf374b",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = gpflow.models.SGPR(\n",
    "    (X_train, Y_train),\n",
    "    kernel=gpflow.kernels.RBF(),\n",
    "    # kernel=gpflow.kernels.SquaredExponential(),\n",
    "    inducing_variable=inducing_variable,\n",
    "        # noise_variance=1.0  # Set to some reasonable initial value\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "90667603-a376-4979-9830-21a5a8265ca4",
   "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>SGPR.kernel.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>SGPR.kernel.lengthscales</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>SGPR.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",
       "<tr><td>SGPR.inducing_variable.Z</td><td>Parameter</td><td>Identity        </td><td>       </td><td>True       </td><td>(240, 22)</td><td>float64</td><td>[[1., 229.16, 184.32...</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n"
     ]
    }
   ],
   "source": [
    "print(gpflow.utilities.print_summary(model, \"notebook\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7baa1ed8-f876-4e65-ae43-02ec16a5433f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "Elapsed time: 0.16 minutes\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "\n",
    "start_time = time.time()\n",
    "# Instantiate an Adam optimizer with a set learning rate\n",
    "adam_optimizer = tf.optimizers.Adam(learning_rate=0.01)\n",
    "\n",
    "# Run the optimization loop\n",
    "# Perform a number of optimization steps\n",
    "num_iterations = 1000  # Choose an appropriate number of iterations for your problem\n",
    "for i in range(num_iterations):\n",
    "    adam_optimizer.minimize(model.training_loss, model.trainable_variables)\n",
    "    if (i % 100 ==0):\n",
    "        print(i)\n",
    "        elapsed_time_minutes = (time.time() - start_time) / 60\n",
    "        print(f\"Elapsed time: {elapsed_time_minutes:.2f} minutes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "8e1ee2f5-1791-4131-b0f2-3c50cc2d08b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0, Elapsed time: 1.08 minutes\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_2191542/110709553.py\u001b[0m in \u001b[0;36m<cell line: 23>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     25\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGradientTape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mtape\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m             \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtraining_loss\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_batch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY_batch\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---> 27\u001b[0;31m         \u001b[0mgradients\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtape\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgradient\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\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     28\u001b[0m         \u001b[0madam_optimizer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_gradients\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgradients\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/condaforge/lib/python3.10/site-packages/tensorflow/python/eager/backprop.py\u001b[0m in \u001b[0;36mgradient\u001b[0;34m(self, target, sources, output_gradients, unconnected_gradients)\u001b[0m\n\u001b[1;32m   1111\u001b[0m                           for x in output_gradients]\n\u001b[1;32m   1112\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1113\u001b[0;31m     flat_grad = imperative_grad.imperative_grad(\n\u001b[0m\u001b[1;32m   1114\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_tape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1115\u001b[0m         \u001b[0mflat_targets\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/condaforge/lib/python3.10/site-packages/tensorflow/python/eager/imperative_grad.py\u001b[0m in \u001b[0;36mimperative_grad\u001b[0;34m(tape, target, sources, output_gradients, sources_raw, unconnected_gradients)\u001b[0m\n\u001b[1;32m     65\u001b[0m         \"Unknown value for unconnected_gradients: %r\" % unconnected_gradients)\n\u001b[1;32m     66\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 67\u001b[0;31m   return pywrap_tfe.TFE_Py_TapeGradient(\n\u001b[0m\u001b[1;32m     68\u001b[0m       \u001b[0mtape\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_tape\u001b[0m\u001b[0;34m,\u001b[0m  \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m       \u001b[0mtarget\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/condaforge/lib/python3.10/site-packages/tensorflow/python/eager/backprop.py\u001b[0m in \u001b[0;36m_gradient_function\u001b[0;34m(op_name, attr_tuple, num_inputs, inputs, outputs, out_grads, skip_input_indices, forward_pass_name_scope)\u001b[0m\n\u001b[1;32m    158\u001b[0m       \u001b[0mgradient_name_scope\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mforward_pass_name_scope\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    159\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname_scope\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgradient_name_scope\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--> 160\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mgrad_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmock_op\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mout_grads\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    161\u001b[0m   \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mgrad_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmock_op\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mout_grads\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/condaforge/lib/python3.10/site-packages/tensorflow/python/ops/math_grad.py\u001b[0m in \u001b[0;36m_RealDivGrad\u001b[0;34m(op, grad)\u001b[0m\n\u001b[1;32m   1457\u001b[0m   \u001b[0msx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marray_ops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\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[1;32m   1458\u001b[0m   \u001b[0msy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marray_ops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1459\u001b[0;31m   \u001b[0mrx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mry\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgen_array_ops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbroadcast_gradient_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msy\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   1460\u001b[0m   \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath_ops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconj\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[1;32m   1461\u001b[0m   \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath_ops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/condaforge/lib/python3.10/site-packages/tensorflow/python/ops/gen_array_ops.py\u001b[0m in \u001b[0;36mbroadcast_gradient_args\u001b[0;34m(s0, s1, name)\u001b[0m\n\u001b[1;32m    755\u001b[0m   \u001b[0;32mif\u001b[0m \u001b[0mtld\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_eager\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    756\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[0;32m--> 757\u001b[0;31m       _result = pywrap_tfe.TFE_Py_FastPathExecute(\n\u001b[0m\u001b[1;32m    758\u001b[0m         _ctx, \"BroadcastGradientArgs\", name, s0, s1)\n\u001b[1;32m    759\u001b[0m       \u001b[0m_result\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_BroadcastGradientArgsOutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_make\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_result\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Mini-batch training can be handled by SVGP in GPflow, which allows online training with mini-batches\n",
    "model = gpflow.models.SVGP(\n",
    "    kernel=gpflow.kernels.RBF(),\n",
    "    likelihood=gpflow.likelihoods.Gaussian(),\n",
    "    inducing_variable=inducing_variable,\n",
    "    num_data=X_train.shape[0]  # Provide the total number of data points for scaling\n",
    ")\n",
    "start_time = time.time()\n",
    "\n",
    "# Instantiate an Adam optimizer with a set learning rate\n",
    "initial_learning_rate = 0.1\n",
    "adam_optimizer = tf.optimizers.Adam(learning_rate=initial_learning_rate)\n",
    "\n",
    "# Perform optimization using the Adam optimizer\n",
    "num_iterations = 1000\n",
    "batch_size = 1000\n",
    "\n",
    "# Data pipeline for mini-batch training\n",
    "train_dataset = tf.data.Dataset.from_tensor_slices((X_train, Y_train))\n",
    "train_dataset = train_dataset.shuffle(buffer_size=len(X_train)).batch(batch_size)\n",
    "\n",
    "# Optimization loop\n",
    "for i in range(num_iterations):\n",
    "    for (X_batch, Y_batch) in train_dataset:\n",
    "        with tf.GradientTape() as tape:\n",
    "            loss = model.training_loss((X_batch, Y_batch))\n",
    "        gradients = tape.gradient(loss, model.trainable_variables)\n",
    "        adam_optimizer.apply_gradients(zip(gradients, model.trainable_variables))\n",
    "    \n",
    "    # Log progress every 100 iterations\n",
    "    if i % 100 == 0:\n",
    "        elapsed_time_minutes = (time.time() - start_time) / 60\n",
    "        print(f\"Iteration {i}, Elapsed time: {elapsed_time_minutes:.2f} minutes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc9f6b8-4d7d-432f-8648-ca180f02074f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gpflow.utilities.print_summary(model, \"notebook\"))\n",
    "###save model\n",
    "model.compiled_predict_f = tf.function(\n",
    "    lambda Xnew: model.predict_f(Xnew, full_cov=False),\n",
    "    input_signature=[tf.TensorSpec(shape=[None, X_train.shape[1]], dtype=tf.float64)],\n",
    ")\n",
    "model.compiled_predict_y = tf.function(\n",
    "    lambda Xnew: model.predict_y(Xnew, full_cov=False),\n",
    "    input_signature=[tf.TensorSpec(shape=[None, X_train.shape[1]], dtype=tf.float64)],\n",
    ")\n",
    "\n",
    "save_dir = \"/home/smmrrr/TEM_Analysis/TEM_Analysis/models/\"\n",
    "save_name = \"test_gp\"\n",
    "tf.saved_model.save(model, save_dir+save_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18d69bda-bbd1-40e9-8072-cce34da79775",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_train_sub[var+'_pred'] = xgb_model.predict(X_train_sub)\n",
    "Y_test[var+'_pred'] = xgb_model.predict(X_test)\n",
    "\n",
    "r2_train = r2_score(Y_train_sub[var], Y_train_sub[var+'_pred'])\n",
    "r2_test = r2_score(Y_test[var], Y_test[var+'_pred'])\n",
    "rmse_train = np.sqrt(mean_squared_error(Y_train_sub[var], Y_train_sub[var+'_pred']))\n",
    "rmse_test = np.sqrt(mean_squared_error(Y_test[var], Y_test[var+'_pred']))\n",
    "bias_train = np.mean(Y_train_sub[var+'_pred'] - Y_train_sub[var])\n",
    "bias_test = np.mean(Y_test[var+'_pred'] - Y_test[var])\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "condaforge",
   "language": "python",
   "name": "condaforge"
  },
  "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}