notebooks/Sparse GP Regression for Big Molecular Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "213d6cca",
   "metadata": {},
   "source": [
    "# Sparse GP Regression on Molecules #\n",
    "\n",
    "An example notebook for sparse GP regression to enable scalability to large molecular datasests.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9b75e9c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") # Turn off Graphein warnings\n",
    "\n",
    "import time\n",
    "\n",
    "from botorch import fit_gpytorch_model\n",
    "import gpytorch\n",
    "from mordred import Calculator, descriptors\n",
    "import numpy as np\n",
    "from rdkit import Chem\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import torch\n",
    "\n",
    "from gauche.dataloader import MolPropLoader\n",
    "from gauche.dataloader.data_utils import transform_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "856e6d07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We define our sparse GP model using and inducing point kernel wrapped around the RQ kernel\n",
    "\n",
    "from gpytorch.means import ConstantMean\n",
    "from gpytorch.kernels import ScaleKernel, InducingPointKernel, RQKernel\n",
    "from gpytorch.distributions import MultivariateNormal\n",
    "\n",
    "class SparseGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(SparseGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = ConstantMean()\n",
    "        self.base_covar_module = ScaleKernel(RQKernel())\n",
    "        self.covar_module = InducingPointKernel(self.base_covar_module, inducing_points=train_x[:100, :].clone(), likelihood=likelihood)\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return MultivariateNormal(mean_x, covar_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97f34aff",
   "metadata": {},
   "source": [
    "## Sparse GP Regression on the Lipophilicity Dataset ##\n",
    "\n",
    "We define our experiment parameters. In this case we are working on the large lipophilicity dataset [1] containing 4200 molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bb55fb31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regression experiments parameters, number of random splits and split size\n",
    "\n",
    "n_trials = 20\n",
    "test_set_size = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4869a3a8",
   "metadata": {},
   "source": [
    "Load the Lipophilicity Dataset via the DataLoaderMP class which contains several molecular property prediction benchmark datasets!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a3b1f1de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mordred descriptor computation takes 639.975240945816 seconds\n"
     ]
    }
   ],
   "source": [
    "# Load the Lipophilicity dataset\n",
    "\n",
    "loader = MolPropLoader()\n",
    "loader.load_benchmark(\"Lipophilicity\")\n",
    "\n",
    "# Mordred descriptor computation is expensive\n",
    "calc = Calculator(descriptors, ignore_3D=False)\n",
    "mols = [Chem.MolFromSmiles(smi) for smi in loader.features]\n",
    "t0 = time.time()\n",
    "X_mordred = [calc(mol) for mol in mols]\n",
    "t1 = time.time()\n",
    "print(f'Mordred descriptor computation takes {t1 - t0} seconds')\n",
    "X_mordred = np.array(X_mordred).astype(np.float64)\n",
    "y = loader.labels\n",
    "\n",
    "\"\"\"Collect nan indices\"\"\"\n",
    "\n",
    "nan_dims = []\n",
    "\n",
    "for i in range(len(X_mordred)):\n",
    "    nan_indices = list(np.where(np.isnan(X_mordred[i, :]))[0])\n",
    "    for dim in nan_indices:\n",
    "        if dim not in nan_dims:\n",
    "            nan_dims.append(dim)\n",
    "            \n",
    "X_mordred = np.delete(X_mordred, nan_dims, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74d8d83b",
   "metadata": {},
   "source": [
    "## Model Evaluation ##\n",
    "\n",
    "Here we define a training/evaluation loop assessing performance using the root mean-square error (RMSE), mean average error (MAE), and $R^2$ metrics. The `evaluate_model` function also computes the GP confidence-error curve which will be explained below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0826fcff",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") # Turn off GPyTorch warnings\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "def evaluate_model(X, y):\n",
    "    \"\"\"Helper function for model evaluation.\n",
    "    \n",
    "    Args:\n",
    "        X: n x d NumPy array of inputs representing molecules\n",
    "        y: n x 1 NumPy array of output labels\n",
    "    Returns:\n",
    "        regression metrics and confidence-error curve plot.\n",
    "    \"\"\"\n",
    "\n",
    "    # initialise performance metric lists\n",
    "    r2_list = []\n",
    "    rmse_list = []\n",
    "    mae_list = []\n",
    "    \n",
    "    # We pre-allocate array for plotting confidence-error curves\n",
    "\n",
    "    _, _, _, y_test = train_test_split(X, y, test_size=test_set_size)  # To get test set size\n",
    "    n_test = len(y_test)\n",
    "\n",
    "    mae_confidence_list = np.zeros((n_trials, n_test))\n",
    "    \n",
    "    print('\\nBeginning training loop...')\n",
    "\n",
    "    for i in range(0, n_trials):\n",
    "        \n",
    "        print(f'Starting trial {i}')\n",
    "                \n",
    "        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_set_size, random_state=i)\n",
    "        \n",
    "        scaler = StandardScaler()\n",
    "        X_train = scaler.fit_transform(X_train)\n",
    "        X_test = scaler.transform(X_test)\n",
    "        pca_mordred = PCA(n_components=51)\n",
    "        X_train = pca_mordred.fit_transform(X_train)\n",
    "        X_test = pca_mordred.transform(X_test)\n",
    "\n",
    "        #  We standardise the outputs\n",
    "        _, y_train, _, y_test, y_scaler = transform_data(X_train, y_train, X_test, y_test)\n",
    "\n",
    "        # Convert numpy arrays to PyTorch tensors and flatten the label vectors\n",
    "        X_train = torch.tensor(X_train.astype(np.float64))\n",
    "        X_test = torch.tensor(X_test.astype(np.float64))\n",
    "        y_train = torch.tensor(y_train).flatten()\n",
    "        y_test = torch.tensor(y_test).flatten()\n",
    "\n",
    "        # initialise GP likelihood and model\n",
    "        likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "        model = SparseGPModel(X_train, y_train, likelihood)\n",
    "\n",
    "        # Find optimal model hyperparameters\n",
    "        # \"Loss\" for GPs - the marginal log likelihood\n",
    "        mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "        # Use the BoTorch utility for fitting GPs in order to use the LBFGS-B optimiser (recommended)\n",
    "        fit_gpytorch_model(mll)\n",
    "\n",
    "        # Get into evaluation (predictive posterior) mode\n",
    "        model.eval()\n",
    "        likelihood.eval()\n",
    "\n",
    "        # mean and variance GP prediction\n",
    "        f_pred = model(X_test)\n",
    "\n",
    "        y_pred = f_pred.mean\n",
    "        y_var = f_pred.variance\n",
    "\n",
    "        # Transform back to real data space to compute metrics and detach gradients. Must unsqueeze dimension\n",
    "        # to make compatible with inverse_transform in scikit-learn version > 1\n",
    "        y_pred = y_scaler.inverse_transform(y_pred.detach().unsqueeze(dim=1))\n",
    "        y_test = y_scaler.inverse_transform(y_test.detach().unsqueeze(dim=1))\n",
    "        \n",
    "        # Compute scores for confidence curve plotting.\n",
    "\n",
    "        ranked_confidence_list = np.argsort(y_var.detach(), axis=0).flatten()\n",
    "\n",
    "        for k in range(len(y_test)):\n",
    "\n",
    "            # Construct the MAE error for each level of confidence\n",
    "\n",
    "            conf = ranked_confidence_list[0:k+1]\n",
    "            mae = mean_absolute_error(y_test[conf], y_pred[conf])\n",
    "            mae_confidence_list[i, k] = mae\n",
    "\n",
    "        # Output Standardised RMSE and RMSE on Train Set\n",
    "        y_train = y_train.detach()\n",
    "        y_pred_train = model(X_train).mean.detach()\n",
    "        train_rmse_stan = np.sqrt(mean_squared_error(y_train, y_pred_train))\n",
    "        train_rmse = np.sqrt(mean_squared_error(y_scaler.inverse_transform(y_train.unsqueeze(dim=1)), \n",
    "                                                y_scaler.inverse_transform(y_pred_train.unsqueeze(dim=1))))\n",
    "\n",
    "        # Compute R^2, RMSE and MAE on Test set\n",
    "        score = r2_score(y_test, y_pred)\n",
    "        rmse = np.sqrt(mean_squared_error(y_test, y_pred))\n",
    "        mae = mean_absolute_error(y_test, y_pred)\n",
    "\n",
    "        r2_list.append(score)\n",
    "        rmse_list.append(rmse)\n",
    "        mae_list.append(mae)\n",
    "        \n",
    "    r2_list = np.array(r2_list)\n",
    "    rmse_list = np.array(rmse_list)\n",
    "    mae_list = np.array(mae_list)\n",
    "        \n",
    "    print(\"\\nmean R^2: {:.4f} +- {:.4f}\".format(np.mean(r2_list), np.std(r2_list)/np.sqrt(len(r2_list))))\n",
    "    print(\"mean RMSE: {:.4f} +- {:.4f}\".format(np.mean(rmse_list), np.std(rmse_list)/np.sqrt(len(rmse_list))))\n",
    "    print(\"mean MAE: {:.4f} +- {:.4f}\\n\".format(np.mean(mae_list), np.std(mae_list)/np.sqrt(len(mae_list)))) \n",
    "    \n",
    "    # Plot confidence-error curves\n",
    "\n",
    "    # 1e-14 instead of 0 to for numerical reasons!\n",
    "    confidence_percentiles = np.arange(1e-14, 100, 100/len(y_test))  \n",
    "\n",
    "    # We plot the Mean-absolute error confidence-error curves\n",
    "\n",
    "    mae_mean = np.mean(mae_confidence_list, axis=0)\n",
    "    mae_std = np.std(mae_confidence_list, axis=0)\n",
    "\n",
    "    mae_mean = np.flip(mae_mean)\n",
    "    mae_std = np.flip(mae_std)\n",
    "\n",
    "    # 1 sigma errorbars\n",
    "\n",
    "    lower = mae_mean - mae_std\n",
    "    upper = mae_mean + mae_std\n",
    "\n",
    "    plt.plot(confidence_percentiles, mae_mean, label='mean')\n",
    "    plt.fill_between(confidence_percentiles, lower, upper, alpha=0.2)\n",
    "    plt.xlabel('Confidence Percentile')\n",
    "    plt.ylabel('MAE (nm)')\n",
    "    plt.ylim([0, np.max(upper) + 1])\n",
    "    plt.xlim([0, 100 * ((len(y_test) - 1) / len(y_test))])\n",
    "    plt.yticks(np.arange(0, np.max(upper) + 1, 5.0))\n",
    "    plt.show()\n",
    "    \n",
    "    return rmse_list, mae_list\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba9bcaf5",
   "metadata": {},
   "source": [
    "Check the perfomance achieved by our sparse GP model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fa469e07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Beginning training loop...\n",
      "Starting trial 0\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[14], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m rmse_mordred, mae_mordred \u001b[38;5;241m=\u001b[39m \u001b[43mevaluate_model\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX_mordred\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[13], line 63\u001b[0m, in \u001b[0;36mevaluate_model\u001b[0;34m(X, y)\u001b[0m\n\u001b[1;32m     60\u001b[0m mll \u001b[38;5;241m=\u001b[39m gpytorch\u001b[38;5;241m.\u001b[39mmlls\u001b[38;5;241m.\u001b[39mExactMarginalLogLikelihood(likelihood, model)\n\u001b[1;32m     62\u001b[0m \u001b[38;5;66;03m# Use the BoTorch utility for fitting GPs in order to use the LBFGS-B optimiser (recommended)\u001b[39;00m\n\u001b[0;32m---> 63\u001b[0m \u001b[43mfit_gpytorch_model\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmll\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     65\u001b[0m \u001b[38;5;66;03m# Get into evaluation (predictive posterior) mode\u001b[39;00m\n\u001b[1;32m     66\u001b[0m model\u001b[38;5;241m.\u001b[39meval()\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/fit.py:164\u001b[0m, in \u001b[0;36mfit_gpytorch_model\u001b[0;34m(mll, optimizer, optimizer_kwargs, exclude, max_retries, **kwargs)\u001b[0m\n\u001b[1;32m    158\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m (\n\u001b[1;32m    159\u001b[0m     nullcontext()\n\u001b[1;32m    160\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m exclude \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m    161\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m requires_grad_ctx(mll, assignments\u001b[38;5;241m=\u001b[39m{name: \u001b[38;5;28;01mFalse\u001b[39;00m \u001b[38;5;28;01mfor\u001b[39;00m name \u001b[38;5;129;01min\u001b[39;00m exclude})\n\u001b[1;32m    162\u001b[0m ):\n\u001b[1;32m    163\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 164\u001b[0m         mll \u001b[38;5;241m=\u001b[39m \u001b[43mfit_gpytorch_mll\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    165\u001b[0m \u001b[43m            \u001b[49m\u001b[43mmll\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    166\u001b[0m \u001b[43m            \u001b[49m\u001b[43moptimizer\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptimizer\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    167\u001b[0m \u001b[43m            \u001b[49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    168\u001b[0m \u001b[43m            \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    169\u001b[0m \u001b[43m        \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    170\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m ModelFittingError \u001b[38;5;28;01mas\u001b[39;00m err:\n\u001b[1;32m    171\u001b[0m         warn(\u001b[38;5;28mstr\u001b[39m(err), \u001b[38;5;167;01mRuntimeWarning\u001b[39;00m)\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/fit.py:105\u001b[0m, in \u001b[0;36mfit_gpytorch_mll\u001b[0;34m(mll, closure, optimizer, closure_kwargs, optimizer_kwargs, **kwargs)\u001b[0m\n\u001b[1;32m    102\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m optimizer \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:  \u001b[38;5;66;03m# defer to per-method defaults\u001b[39;00m\n\u001b[1;32m    103\u001b[0m     kwargs[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moptimizer\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m optimizer\n\u001b[0;32m--> 105\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mFitGPyTorchMLL\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    106\u001b[0m \u001b[43m    \u001b[49m\u001b[43mmll\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    107\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;28;43mtype\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mmll\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlikelihood\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    108\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;28;43mtype\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mmll\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmodel\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    109\u001b[0m \u001b[43m    \u001b[49m\u001b[43mclosure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mclosure\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    110\u001b[0m \u001b[43m    \u001b[49m\u001b[43mclosure_kwargs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mclosure_kwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    111\u001b[0m \u001b[43m    \u001b[49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    112\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    113\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/utils/dispatcher.py:93\u001b[0m, in \u001b[0;36mDispatcher.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     91\u001b[0m func \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__getitem__\u001b[39m(types\u001b[38;5;241m=\u001b[39mtypes)\n\u001b[1;32m     92\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m---> 93\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     94\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m MDNotImplementedError:\n\u001b[1;32m     95\u001b[0m     \u001b[38;5;66;03m# Traverses registered methods in order, yields whenever a match is found\u001b[39;00m\n\u001b[1;32m     96\u001b[0m     funcs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdispatch_iter(\u001b[38;5;241m*\u001b[39mtypes)\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/fit.py:252\u001b[0m, in \u001b[0;36m_fit_fallback\u001b[0;34m(mll, _, __, closure, optimizer, closure_kwargs, optimizer_kwargs, max_attempts, warning_handler, caught_exception_types, **ignore)\u001b[0m\n\u001b[1;32m    250\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m catch_warnings(record\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m) \u001b[38;5;28;01mas\u001b[39;00m warning_list, debug(\u001b[38;5;28;01mTrue\u001b[39;00m):\n\u001b[1;32m    251\u001b[0m     simplefilter(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124malways\u001b[39m\u001b[38;5;124m\"\u001b[39m, category\u001b[38;5;241m=\u001b[39mOptimizationWarning)\n\u001b[0;32m--> 252\u001b[0m     \u001b[43moptimizer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmll\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mclosure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mclosure\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    254\u001b[0m \u001b[38;5;66;03m# Resolved warnings and determine whether or not to retry\u001b[39;00m\n\u001b[1;32m    255\u001b[0m done \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/optim/fit.py:92\u001b[0m, in \u001b[0;36mfit_gpytorch_mll_scipy\u001b[0;34m(mll, parameters, bounds, closure, closure_kwargs, method, options, callback, timeout_sec)\u001b[0m\n\u001b[1;32m     89\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m closure_kwargs \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     90\u001b[0m     closure \u001b[38;5;241m=\u001b[39m partial(closure, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mclosure_kwargs)\n\u001b[0;32m---> 92\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mscipy_minimize\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m     93\u001b[0m \u001b[43m    \u001b[49m\u001b[43mclosure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mclosure\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     94\u001b[0m \u001b[43m    \u001b[49m\u001b[43mparameters\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mparameters\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     95\u001b[0m \u001b[43m    \u001b[49m\u001b[43mbounds\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mbounds\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     96\u001b[0m \u001b[43m    \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     97\u001b[0m \u001b[43m    \u001b[49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     98\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcallback\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcallback\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     99\u001b[0m \u001b[43m    \u001b[49m\u001b[43mtimeout_sec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtimeout_sec\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    100\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    101\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m result\u001b[38;5;241m.\u001b[39mstatus \u001b[38;5;241m!=\u001b[39m OptimizationStatus\u001b[38;5;241m.\u001b[39mSUCCESS:\n\u001b[1;32m    102\u001b[0m     warn(\n\u001b[1;32m    103\u001b[0m         \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m`scipy_minimize` terminated with status \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mresult\u001b[38;5;241m.\u001b[39mstatus\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m, displaying\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    104\u001b[0m         \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m original message from `scipy.optimize.minimize`: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mresult\u001b[38;5;241m.\u001b[39mmessage\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    105\u001b[0m         OptimizationWarning,\n\u001b[1;32m    106\u001b[0m     )\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/optim/core.py:109\u001b[0m, in \u001b[0;36mscipy_minimize\u001b[0;34m(closure, parameters, bounds, callback, x0, method, options, timeout_sec)\u001b[0m\n\u001b[1;32m    101\u001b[0m         result \u001b[38;5;241m=\u001b[39m OptimizationResult(\n\u001b[1;32m    102\u001b[0m             step\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mnext\u001b[39m(call_counter),\n\u001b[1;32m    103\u001b[0m             fval\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mfloat\u001b[39m(wrapped_closure(x)[\u001b[38;5;241m0\u001b[39m]),\n\u001b[1;32m    104\u001b[0m             status\u001b[38;5;241m=\u001b[39mOptimizationStatus\u001b[38;5;241m.\u001b[39mRUNNING,\n\u001b[1;32m    105\u001b[0m             runtime\u001b[38;5;241m=\u001b[39mmonotonic() \u001b[38;5;241m-\u001b[39m start_time,\n\u001b[1;32m    106\u001b[0m         )\n\u001b[1;32m    107\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m callback(parameters, result)  \u001b[38;5;66;03m# pyre-ignore [29]\u001b[39;00m\n\u001b[0;32m--> 109\u001b[0m raw \u001b[38;5;241m=\u001b[39m \u001b[43mminimize_with_timeout\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    110\u001b[0m \u001b[43m    \u001b[49m\u001b[43mwrapped_closure\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    111\u001b[0m \u001b[43m    \u001b[49m\u001b[43mwrapped_closure\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mstate\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mif\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01mis\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43;01melse\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mastype\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnp_float64\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcopy\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    112\u001b[0m \u001b[43m    \u001b[49m\u001b[43mjac\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    113\u001b[0m \u001b[43m    \u001b[49m\u001b[43mbounds\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mbounds_np\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    114\u001b[0m \u001b[43m    \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    115\u001b[0m \u001b[43m    \u001b[49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    116\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcallback\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mwrapped_callback\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    117\u001b[0m \u001b[43m    \u001b[49m\u001b[43mtimeout_sec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtimeout_sec\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    118\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    120\u001b[0m \u001b[38;5;66;03m# Post-processing and outcome handling\u001b[39;00m\n\u001b[1;32m    121\u001b[0m wrapped_closure\u001b[38;5;241m.\u001b[39mstate \u001b[38;5;241m=\u001b[39m asarray(raw\u001b[38;5;241m.\u001b[39mx)  \u001b[38;5;66;03m# set parameter state to optimal values\u001b[39;00m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/optim/utils/timeout.py:80\u001b[0m, in \u001b[0;36mminimize_with_timeout\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options, timeout_sec)\u001b[0m\n\u001b[1;32m     77\u001b[0m     wrapped_callback \u001b[38;5;241m=\u001b[39m callback\n\u001b[1;32m     79\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m---> 80\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43moptimize\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mminimize\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m     81\u001b[0m \u001b[43m        \u001b[49m\u001b[43mfun\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfun\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     82\u001b[0m \u001b[43m        \u001b[49m\u001b[43mx0\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mx0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     83\u001b[0m \u001b[43m        \u001b[49m\u001b[43margs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     84\u001b[0m \u001b[43m        \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     85\u001b[0m \u001b[43m        \u001b[49m\u001b[43mjac\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mjac\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     86\u001b[0m \u001b[43m        \u001b[49m\u001b[43mhess\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mhess\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     87\u001b[0m \u001b[43m        \u001b[49m\u001b[43mhessp\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mhessp\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     88\u001b[0m \u001b[43m        \u001b[49m\u001b[43mbounds\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mbounds\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     89\u001b[0m \u001b[43m        \u001b[49m\u001b[43mconstraints\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mconstraints\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     90\u001b[0m \u001b[43m        \u001b[49m\u001b[43mtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtol\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     91\u001b[0m \u001b[43m        \u001b[49m\u001b[43mcallback\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mwrapped_callback\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     92\u001b[0m \u001b[43m        \u001b[49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     93\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     94\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m OptimizationTimeoutError \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[1;32m     95\u001b[0m     msg \u001b[38;5;241m=\u001b[39m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mOptimization timed out after \u001b[39m\u001b[38;5;132;01m{\u001b[39;00me\u001b[38;5;241m.\u001b[39mruntime\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m seconds.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_minimize.py:710\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    707\u001b[0m     res \u001b[38;5;241m=\u001b[39m _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,\n\u001b[1;32m    708\u001b[0m                              \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n\u001b[1;32m    709\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124ml-bfgs-b\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m--> 710\u001b[0m     res \u001b[38;5;241m=\u001b[39m \u001b[43m_minimize_lbfgsb\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfun\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mjac\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbounds\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    711\u001b[0m \u001b[43m                           \u001b[49m\u001b[43mcallback\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcallback\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    712\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtnc\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m    713\u001b[0m     res \u001b[38;5;241m=\u001b[39m _minimize_tnc(fun, x0, args, jac, bounds, callback\u001b[38;5;241m=\u001b[39mcallback,\n\u001b[1;32m    714\u001b[0m                         \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_lbfgsb_py.py:365\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    359\u001b[0m task_str \u001b[38;5;241m=\u001b[39m task\u001b[38;5;241m.\u001b[39mtobytes()\n\u001b[1;32m    360\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m task_str\u001b[38;5;241m.\u001b[39mstartswith(\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFG\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[1;32m    361\u001b[0m     \u001b[38;5;66;03m# The minimization routine wants f and g at the current x.\u001b[39;00m\n\u001b[1;32m    362\u001b[0m     \u001b[38;5;66;03m# Note that interruptions due to maxfun are postponed\u001b[39;00m\n\u001b[1;32m    363\u001b[0m     \u001b[38;5;66;03m# until the completion of the current minimization iteration.\u001b[39;00m\n\u001b[1;32m    364\u001b[0m     \u001b[38;5;66;03m# Overwrite f and g:\u001b[39;00m\n\u001b[0;32m--> 365\u001b[0m     f, g \u001b[38;5;241m=\u001b[39m \u001b[43mfunc_and_grad\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    366\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m task_str\u001b[38;5;241m.\u001b[39mstartswith(\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNEW_X\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[1;32m    367\u001b[0m     \u001b[38;5;66;03m# new iteration\u001b[39;00m\n\u001b[1;32m    368\u001b[0m     n_iterations \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:285\u001b[0m, in \u001b[0;36mScalarFunction.fun_and_grad\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    283\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39marray_equal(x, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx):\n\u001b[1;32m    284\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_x_impl(x)\n\u001b[0;32m--> 285\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_update_fun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    286\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_grad()\n\u001b[1;32m    287\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:251\u001b[0m, in \u001b[0;36mScalarFunction._update_fun\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    249\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_update_fun\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m    250\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf_updated:\n\u001b[0;32m--> 251\u001b[0m         \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_update_fun_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    252\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf_updated \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:155\u001b[0m, in \u001b[0;36mScalarFunction.__init__.<locals>.update_fun\u001b[0;34m()\u001b[0m\n\u001b[1;32m    154\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mupdate_fun\u001b[39m():\n\u001b[0;32m--> 155\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf \u001b[38;5;241m=\u001b[39m \u001b[43mfun_wrapped\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:137\u001b[0m, in \u001b[0;36mScalarFunction.__init__.<locals>.fun_wrapped\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    133\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnfev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m    134\u001b[0m \u001b[38;5;66;03m# Send a copy because the user may overwrite it.\u001b[39;00m\n\u001b[1;32m    135\u001b[0m \u001b[38;5;66;03m# Overwriting results in undefined behaviour because\u001b[39;00m\n\u001b[1;32m    136\u001b[0m \u001b[38;5;66;03m# fun(self.x) will change self.x, with the two no longer linked.\u001b[39;00m\n\u001b[0;32m--> 137\u001b[0m fx \u001b[38;5;241m=\u001b[39m \u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcopy\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    138\u001b[0m \u001b[38;5;66;03m# Make sure the function returns a true scalar\u001b[39;00m\n\u001b[1;32m    139\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39misscalar(fx):\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_optimize.py:77\u001b[0m, in \u001b[0;36mMemoizeJac.__call__\u001b[0;34m(self, x, *args)\u001b[0m\n\u001b[1;32m     75\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__call__\u001b[39m(\u001b[38;5;28mself\u001b[39m, x, \u001b[38;5;241m*\u001b[39margs):\n\u001b[1;32m     76\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\" returns the function value \"\"\"\u001b[39;00m\n\u001b[0;32m---> 77\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_compute_if_needed\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     78\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_value\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/scipy/optimize/_optimize.py:71\u001b[0m, in \u001b[0;36mMemoizeJac._compute_if_needed\u001b[0;34m(self, x, *args)\u001b[0m\n\u001b[1;32m     69\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39mall(x \u001b[38;5;241m==\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx) \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_value \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mjac \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     70\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(x)\u001b[38;5;241m.\u001b[39mcopy()\n\u001b[0;32m---> 71\u001b[0m     fg \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     72\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mjac \u001b[38;5;241m=\u001b[39m fg[\u001b[38;5;241m1\u001b[39m]\n\u001b[1;32m     73\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_value \u001b[38;5;241m=\u001b[39m fg[\u001b[38;5;241m0\u001b[39m]\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/optim/closures/core.py:150\u001b[0m, in \u001b[0;36mNdarrayOptimizationClosure.__call__\u001b[0;34m(self, state, **kwargs)\u001b[0m\n\u001b[1;32m    147\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstate \u001b[38;5;241m=\u001b[39m state\n\u001b[1;32m    149\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 150\u001b[0m     value_tensor, grad_tensors \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mclosure\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    151\u001b[0m     value \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mas_array(value_tensor)\n\u001b[1;32m    152\u001b[0m     grads \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_gradient_ndarray(fill_value\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfill_value)\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/botorch/optim/closures/core.py:66\u001b[0m, in \u001b[0;36mForwardBackwardClosure.__call__\u001b[0;34m(self, **kwargs)\u001b[0m\n\u001b[1;32m     64\u001b[0m values \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mforward(\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m     65\u001b[0m value \u001b[38;5;241m=\u001b[39m values \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mreducer \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mreducer(values)\n\u001b[0;32m---> 66\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbackward\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvalue\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     68\u001b[0m grads \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(param\u001b[38;5;241m.\u001b[39mgrad \u001b[38;5;28;01mfor\u001b[39;00m param \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mparameters\u001b[38;5;241m.\u001b[39mvalues())\n\u001b[1;32m     69\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcallback:\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/torch/_tensor.py:492\u001b[0m, in \u001b[0;36mTensor.backward\u001b[0;34m(self, gradient, retain_graph, create_graph, inputs)\u001b[0m\n\u001b[1;32m    482\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m has_torch_function_unary(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m    483\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m handle_torch_function(\n\u001b[1;32m    484\u001b[0m         Tensor\u001b[38;5;241m.\u001b[39mbackward,\n\u001b[1;32m    485\u001b[0m         (\u001b[38;5;28mself\u001b[39m,),\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    490\u001b[0m         inputs\u001b[38;5;241m=\u001b[39minputs,\n\u001b[1;32m    491\u001b[0m     )\n\u001b[0;32m--> 492\u001b[0m \u001b[43mtorch\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mautograd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbackward\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    493\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgradient\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mretain_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcreate_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minputs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minputs\u001b[49m\n\u001b[1;32m    494\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/gauche/lib/python3.11/site-packages/torch/autograd/__init__.py:251\u001b[0m, in \u001b[0;36mbackward\u001b[0;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)\u001b[0m\n\u001b[1;32m    246\u001b[0m     retain_graph \u001b[38;5;241m=\u001b[39m create_graph\n\u001b[1;32m    248\u001b[0m \u001b[38;5;66;03m# The reason we repeat the same comment below is that\u001b[39;00m\n\u001b[1;32m    249\u001b[0m \u001b[38;5;66;03m# some Python versions print out the first line of a multi-line function\u001b[39;00m\n\u001b[1;32m    250\u001b[0m \u001b[38;5;66;03m# calls in the traceback and some print out the last line\u001b[39;00m\n\u001b[0;32m--> 251\u001b[0m \u001b[43mVariable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_execution_engine\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_backward\u001b[49m\u001b[43m(\u001b[49m\u001b[43m  \u001b[49m\u001b[38;5;66;43;03m# Calls into the C++ engine to run the backward pass\u001b[39;49;00m\n\u001b[1;32m    252\u001b[0m \u001b[43m    \u001b[49m\u001b[43mtensors\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    253\u001b[0m \u001b[43m    \u001b[49m\u001b[43mgrad_tensors_\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    254\u001b[0m \u001b[43m    \u001b[49m\u001b[43mretain_graph\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    255\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcreate_graph\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    256\u001b[0m \u001b[43m    \u001b[49m\u001b[43minputs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    257\u001b[0m \u001b[43m    \u001b[49m\u001b[43mallow_unreachable\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    258\u001b[0m \u001b[43m    \u001b[49m\u001b[43maccumulate_grad\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m    259\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "rmse_mordred, mae_mordred = evaluate_model(X_mordred, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81ff2536",
   "metadata": {},
   "source": [
    "## References \n",
    "\n",
    "[1] Anna Gaulton, Louisa J Bellis, A Patricia Bento, Jon Chambers, Mark Davies, Anne Hersey,\n",
    "Yvonne Light, Shaun McGlinchey, David Michalovich, Bissan Al-Lazikani, et al. [ChEMBL:\n",
    "A large-scale bioactivity database for drug discovery](https://pubmed.ncbi.nlm.nih.gov/21948594/). Nucleic Acids Research, 2012.\n",
    "\n",
    "[2] Bajusz, D., Rácz, A. and Héberger, K., 2015. [Why is Tanimoto index an appropriate choice for fingerprint-based similarity calculations?](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-015-0069-3). Journal of Cheminformatics], 7(1), pp.1-13."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}