Wigner_basic.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dad972af",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "#torch.set_default_dtype(torch.float64)\n",
    "import numpy as np\n",
    "import ase.io\n",
    "\n",
    "from pytorch_prototype.code_pytorch import *\n",
    "from pytorch_prototype.utilities import *\n",
    "from pytorch_prototype.miscellaneous import ClebschGordan\n",
    "\n",
    "from rascal.neighbourlist.structure_manager import mask_center_atoms_by_species\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from rascal.models.krr import SparseGPRSolver\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import tqdm\n",
    "import time\n",
    "\n",
    "from rascal.utils import (get_radial_basis_covariance, get_radial_basis_pca, \n",
    "                          get_radial_basis_projections, get_optimal_radial_basis_hypers )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "00841930",
   "metadata": {},
   "outputs": [],
   "source": [
    "METHANE_PATH = 'methane.extxyz'\n",
    "HARTREE_TO_EV = 27.211386245988\n",
    "\n",
    "train_slice = '0:900'\n",
    "test_slice = '900:1000'\n",
    "\n",
    "L_MAX = 6\n",
    "clebsch = ClebschGordan(L_MAX)\n",
    "\n",
    "HYPERS_INITIAL = {\n",
    "    'interaction_cutoff': 6.3,\n",
    "    'max_radial': 20,\n",
    "    'max_angular': L_MAX,\n",
    "    'gaussian_sigma_type': 'Constant',\n",
    "    'gaussian_sigma_constant': 0.2,\n",
    "    'cutoff_smooth_width': 0.3,\n",
    "    'radial_basis': 'DVR'\n",
    "}\n",
    "\n",
    "BATCH_SIZE = 2000\n",
    "DEVICE = 'cuda'\n",
    "\n",
    "\n",
    "\n",
    "structures = process_structures(ase.io.read(METHANE_PATH, index = train_slice))\n",
    "HYPERS = get_optimal_radial_basis_hypers(HYPERS_INITIAL,\n",
    "                                           structures,\n",
    "                                           expanded_max_radial=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "914459e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_wigner_single(first, second):\n",
    "    first_b_size, first_m_size = first.shape[0], first.shape[2]\n",
    "    second_b_size, second_m_size = second.shape[0], second.shape[2]\n",
    "    first = first.transpose(1, 2)\n",
    "    second = second.transpose(1, 2)\n",
    "    first = first.reshape([-1, first.shape[2]])\n",
    "    second = second.reshape([-1, second.shape[2]])\n",
    "    result = torch.matmul(first, second.transpose(0, 1))\n",
    "    result = result.reshape(first_b_size, first_m_size, second_b_size, second_m_size)\n",
    "    return result.transpose(1, 2)\n",
    "\n",
    "def initialize_wigner(first, second):\n",
    "    result = {}\n",
    "    for key in first.keys():\n",
    "        result[str(key)] = initialize_wigner_single(first[key], second[key])\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "51824808",
   "metadata": {},
   "outputs": [],
   "source": [
    "class WignerKernel(torch.nn.Module):\n",
    "    def __init__(self, clebsch, lambda_max, num_iterations):\n",
    "        super(WignerKernel, self).__init__()\n",
    "        main = [WignerCombiningUnrolled(clebsch.precomputed_, lambda_max, algorithm = 'vectorized') \n",
    "                for _ in range(num_iterations)]\n",
    "        self.main = nn.ModuleList(main)\n",
    "        self.last = WignerCombiningUnrolled(clebsch.precomputed_, 0, algorithm = 'vectorized')\n",
    "       \n",
    "            \n",
    "    def forward(self, X):\n",
    "        result = []\n",
    "        wig_now = X\n",
    "        result.append(wig_now['0'][:, 0, 0, None])\n",
    "        for block in self.main:\n",
    "            wig_now = block(wig_now, X)\n",
    "            result.append(wig_now['0'][:, 0, 0, None])\n",
    "        wig_now = self.last(wig_now, X)\n",
    "        result.append(wig_now['0'][:, 0, 0, None])\n",
    "        result = torch.cat(result, dim = -1)\n",
    "        return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "14ff69cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_kernel(model, first, second, batch_size = 1000, device = 'cpu'):\n",
    "    wigner = initialize_wigner(first, second)\n",
    "   \n",
    "    for key in wigner.keys():\n",
    "        initial_shape = [wigner[key].shape[0], wigner[key].shape[1]]\n",
    "        wigner[key] = wigner[key].reshape([-1, wigner[key].shape[2], wigner[key].shape[3]])\n",
    "    ''' for key in wigner.keys():\n",
    "        print(key, wigner[key].shape)'''\n",
    "    \n",
    "    total = initial_shape[0] * initial_shape[1]\n",
    "    result = []\n",
    "    #print(total, batch_size)\n",
    "    #print(initial_shape)\n",
    "    for ind in tqdm.tqdm(range(0, total, batch_size)):\n",
    "        now = {}\n",
    "        for key in wigner.keys():\n",
    "            now[key] = wigner[key][ind : ind + batch_size].to(device)\n",
    "        result_now = model(now).to('cpu')\n",
    "        result.append(result_now)\n",
    "        \n",
    "        \n",
    "    result = torch.cat(result, dim = 0)\n",
    "    return result.reshape(initial_shape + [-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b218587",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 405/405 [03:04<00:00,  2.19it/s]\n",
      "100%|██████████| 45/45 [00:20<00:00,  2.19it/s]\n"
     ]
    }
   ],
   "source": [
    "train_structures = ase.io.read(METHANE_PATH, index = train_slice)\n",
    "test_structures = ase.io.read(METHANE_PATH, index = test_slice)\n",
    "all_species = get_all_species(train_structures + test_structures)\n",
    "\n",
    "for struc in train_structures:\n",
    "    mask_center_atoms_by_species(struc, species_select=[\"C\"])\n",
    "for struc in test_structures:\n",
    "    mask_center_atoms_by_species(struc, species_select=[\"C\"])\n",
    "\n",
    "\n",
    "train_coefs = get_coefs(train_structures, HYPERS, all_species)\n",
    "test_coefs = get_coefs(test_structures, HYPERS, all_species)\n",
    "'''for key in train_coefs.keys():\n",
    "    print(key, train_coefs[key].shape)'''\n",
    "\n",
    "L2_mean = get_L2_mean(train_coefs)\n",
    "#print(L2_mean)\n",
    "for key in train_coefs.keys():\n",
    "    train_coefs[key] /= np.sqrt(L2_mean)\n",
    "    test_coefs[key] /= np.sqrt(L2_mean)\n",
    "\n",
    "'''for key in train_coefs.keys():\n",
    "    train_coefs[key] = train_coefs[key].to(DEVICE)\n",
    "    test_coefs[key] = test_coefs[key].to(DEVICE)\n",
    "'''\n",
    "model = WignerKernel(clebsch, L_MAX, 2)\n",
    "model = model.to(DEVICE)\n",
    "\n",
    "train_train_kernel = compute_kernel(model, train_coefs, train_coefs, batch_size = BATCH_SIZE, device = DEVICE)\n",
    "train_test_kernel = compute_kernel(model, train_coefs, test_coefs, batch_size = BATCH_SIZE, device = DEVICE)\n",
    "\n",
    "train_train_kernel = train_train_kernel.data.cpu().numpy()\n",
    "train_test_kernel = train_test_kernel.data.cpu().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "984879aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.42722014 0.18560453 0.0823369  0.0374895 ]\n",
      "[0.45978022 0.21805915 0.10708366 0.05465367]\n",
      "[0.46098676 0.22094426 0.11126095 0.05959825]\n",
      "[0.43157783 0.19294113 0.09062399 0.04545183]\n",
      "[0.4402474  0.20205282 0.09833136 0.0516623 ]\n",
      "[0.43129173 0.18854591 0.08370247 0.03781608]\n",
      "[0.4258148  0.18502766 0.08266446 0.03834439]\n",
      "[0.44768086 0.20758152 0.10090961 0.05218141]\n",
      "[0.43268204 0.19210142 0.0881346  0.04211764]\n",
      "[0.451247   0.21227995 0.10582463 0.05697034]\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(train_train_kernel[i, i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e60b79d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(900, 900, 4)\n",
      "(900, 100, 4)\n"
     ]
    }
   ],
   "source": [
    "print(train_train_kernel.shape)\n",
    "print(train_test_kernel.shape)\n",
    "train_train_kernel = train_train_kernel[:, :, -1]\n",
    "train_test_kernel = train_test_kernel[:, :, -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0775c7f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_rmse(first, second):\n",
    "    return np.sqrt(np.mean((first - second)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "308cdd47",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_energies = [structure.info['energy'] for structure in train_structures]\n",
    "train_energies = np.array(train_energies) * HARTREE_TO_EV\n",
    "\n",
    "test_energies = [structure.info['energy'] for structure in test_structures]\n",
    "test_energies = np.array(test_energies) * HARTREE_TO_EV\n",
    "\n",
    "mean_e = np.mean(train_energies)\n",
    "train_energies -= mean_e\n",
    "test_energies -= mean_e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a75bab16",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse:  2.384196338104318\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWj0lEQVR4nO3df4wcZ33H8c/X6w1ZJ6gXGpfGmxi7lDoKdWOTK6R11ZKA4kBCcySogdKK/pDcSkUCGh2yC2qCBPKppzZUKqrktgik0JAUzDXFtEfAqVAjJeXMxXFCchTIL9aBHCSXlngh5/O3f9ztZW9vZnd2Z3Znn933S7Lim5ubeTyxP/fcd54f5u4CAIRrQ94NAACkQ5ADQOAIcgAIHEEOAIEjyAEgcBvzuOn555/v27Zty+PWABCsY8eO/dDdNzcezyXIt23bppmZmTxuDQDBMrMnoo5TWgGAwBHkABA4ghwAAkeQA0DgCHIACFwuo1YAhG9qtqLJ6TmdXKhqy0hJ43t3aGx3Oe9m9a1uPi+CHEDbpmYrOnD4hKqLS5KkykJVBw6fkCTCPEK3nxelFQBtm5yeWw2lmurikian53JqUX/r9vMiyAG07eRCta3jw67bz4sgB9C2LSOlto4Pu24/L4IcQNvG9+5QqVhYc6xULGh8746cWtTfuv28eNkJoG21F3SMWkmm28/L8tizc3R01Fk0CwDaY2bH3H208TilFQAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACFxmQW5mBTObNbMvZnVNAEBrWfbI3yfpkQyvBwBIIJNFs8zsQknXSPqYpD/P4poAMCi6vS1eVqsfflzSByW9PO4EM9snaZ8kbd26NaPbAumw7yS6rRfb4qUurZjZtZKecfdjzc5z90PuPuruo5s3b057WyC12j+wykJVrpf+gU3NVvJuGgZIL7bFy6JGvkfSb5vZ45I+K+lKM7stg+sCXcW+k+iFXmyLlzrI3f2Au1/o7tskvVPSUXf/vdQtA7qMfSfRC73YFo9x5Bha7DuJXujFtniZBrm7/6e7X5vlNYFuYd9J9MLY7rIOXr9T5ZGSTFJ5pKSD1+/sy1ErQHDYdxK9Mra73NW/VwQ5hlq3/4EBvUCQA32A8exIgyAHctaLCSMYbIxaAXLGeHakRZADOWM8O9IiyIGcMZ4daRHkQM4Yz460eNkJ5Izx7EiLIAf6AOPZkQalFQAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAEjiAHgMAR5AAQOIIcAAJHkANA4AhyAAgcQQ4AgSPIASBwBDkABI4gB4DAEeQAEDiCHAAClzrIzewiM7vHzL5pZg+b2fuyaBgAIJmNGVzjtKSb3P0bZvZyScfM7G53/2YG1wYAtJC6R+7uT7v7N1Z+/3+SHpFUTntdAEAymdbIzWybpN2S7s/yugCAeJkFuZmdK+nzkt7v7v8b8fl9ZjZjZjPz8/NZ3RYAhl4mQW5mRS2H+Gfc/XDUOe5+yN1H3X108+bNWdwWAKBsRq2YpH+S9Ii7/036JgEA2pFFj3yPpN+XdKWZPbDy660ZXBcAkEDq4Yfu/l+SLIO2AAA6wMxOAAgcQQ4AgSPIASBwBDkABI4gB4DAZbFoFnIwNVvR5PScTi5UtWWkpPG9OzS2eziWuGn8s2/72ZLu++5zWnJXwUzvesNF+ujYzshzh+k5YXiYu/f8pqOjoz4zM9Pz+w6KqdmKDhw+oeri0uqxUrGgg9fvHPiQivqzRznnrILe/rqyPn+sMpTPCYPJzI65+2jjcXrkAZqcnlsXZNXFJU1Ozw1kQNX3qjeYaSlB5+OFF5f0mfueVOOZg/ycMLwI8gCdXKi2dTxkjT3wJCFeE3fmID4nDDeCPEBbRkqqRITRyKai9kwc7et6cLs166ifPtLaMlLK9HpA3hi1EqDxvTtUKhbWHCsWTD/+yWlVFqpySZWFqg4cPqGp2Uo+jYxQ612308a0vefGtSNKxYLG9+5IdU2g3xDkARrbXdbB63eqPFKSSSqPlHTOWRu1eGZtMaFWD+4XzWr7cdL0nkvFgt59+dY1z4kXnRhElFYCNba7vCaQtu8/EnleP9WDO6ntj+/dkWiUikm69cZdDDXEUCLIB0Rc3TyuR5vH+Op22yhptU2T03ORX1szsqm47psbMCworQyIqLp5XD24k1r11GxFeyaOavv+I9ozcbSj2ns7baw3truse/dfqXKTwM9hOgTQNwjyARFVN4+rB7dbq+4k+NO2MUqzwH++uthWW4BBQmllgCQtLbRbq85yAlKa8sfY7rIOHH5Q1cUz6z63ZaTEdHwMLYJ8CDQG3M+UilqI6MHG1aqbBX8vw3NqtqLTZ9bXUIobTFdcvHnNS9HaTw2SCHMMPIJ8wDXOjKwsVFUsmIobbM1wRZN0xcWbI68R95LSTHr/HQ+sfhwXnh+eOqHb738qclGrdkxOz2lxaX2Qn3v2Rt3z6PxQLVsA1KNGPkCiXkhGlUUWl1zFgq2ZLOOSPn+sEln3jnpJKUkRneN1tfYPT53Qbfc9uTq1fsldt933pD48dSK2zXHifjJYOLU4VMsWAI0I8gER90IybsjeqcUzsQtKNWp8SVmw5ntt14fn7fc/FXnO7fc/1fZL1LjSz5aRUtPPAYOOIB8QcS8km0fuenE92NoQwMcmrtGZFmP96sMzbpGrJfe2R880G77Y6dBGYBBQIx8QcQEcF7lm0WOvk/Rg42rmNfXhWWiy7GzcNZp9M5HU9OUqo1YwjAjyAdEqXBu5L/dYGzddSNKDveLizbrtvicjP3feygzLmst/4Tzd+51nE7dLaj3TMy6cmdmJYUVpZUDEvZCMU5uM08nknHsenY88bpJufttr1xx7/EftvWxsNnoGQDR65AMiquzwwk9PR44XN2m17NBuD3ZqthLb83etH7Pd7qiR2uiZ0Ve9gt41kBB7dg6wqP0tTdKvv/oVevxH1bZrya32yzxvU1Gbztq45rqtFruKUx4p6d79V7b9dcAgY8/OPtHLmZBRvfQrLt68ZkPi2pC/mSee1T2PzjdtV7PdeooF0/PVRT13anH1uuOfO64bf/WidRsgJ8H4byA5gryHomZZdnsaeWP5ZM/E0cghf/UbFce1q1m4FjfYulmXi0uuIw8+rYPX71z3zeSeR+eb9tQZ/w0kR5D3UJaLT3Uq6TDFqHbFrdFSbjJi5rlTi7G1+D0TR6On/qv5SocA1mLUSg/1wzTydnq69e2amq3ohRdPrzunuME6Dt2okTYm6d2Xb+VFJ9AGgryH+mEaeTvDFOvb1WzBqrHdZY2UipHXiDsuRa9PfuuNuzpaUAsYZpRWeihq/8leTyNPunVaY7uaLVglSddeekHkJKFrL72gZXvofQPpEOQ91GqKea9GtNTCc/v+I7FT+BsnB7XabzNuklDccQDZySTIzexqSX8rqSDpH919IovrDookAZ3HiJa4cB4pFTU5PacP3PHAantb/TTRD/V/YFilrpGbWUHSJyS9RdIlkt5lZpekve6gSLpUa7srAWYhql5e3GB64cXT69orSTdcVl5dwrZgphsue6ks0g/1f2BYZdEjf72kb7v7dyXJzD4r6TpJ38zg2sFLOuQwTY+205JMVKnn1IunVyf11Lf3lrse1k9Pn1mzQcQd//2Ujjz4tBZOLWpkU3HdrkMsIwv0RhZBXpZUv3vA9yS9IYPrDoSkAd2qBh1laraiW+56eM3Y7nZLMo0vG7fvPxJ5XtT48cUzvhr6z51aVLFgGikV9Xx1kWVkgR7q2ctOM9snaZ8kbd26tVe3zV3SgG53REuzdU+ievxRvXZp/YvXdpfDrbe45DrnZRv1wM1XdfT1ADqTetEsM/s1Sbe4+96Vjw9IkrsfjPuaQV00Ky4soxauci2Pm+501ErcrMj6ezw2cc3qdRvbUCyY5FpXCrnhsvK6tVFM0qazCnrhxdbrpdTu28s1ZYBhEbdoVhZBvlHStyS9SVJF0tcl/a67Pxz3NYMY5FFhWSoWdPD65ckttXHbtRBvPKfdkGs2dFBau3pgq9Bv/LorLt68Zu2VdtS+OcU9C8Ic6FxckKceteLupyW9V9K0pEck3dksxAdVq5ea9+6/UuWRUuINj1tpVjtPOpknysmFqu55dL6jEC8WbHXp2l6PwAGGWSZT9N39S+7+S+7+anf/WBbXDE2Sl5qdjEyZmq1oz8RRbd9/RHsmjq4OW4ybam/2UmjWzm13CGCnNfJzzlqers+YcqC3WGslI0nGUbc71rrZGPTGdUrOWxn+V6uU1Z/bzhDANIW251dGtjCmHOgtgjwjUT3kxhJHknPqxZUobrrzuLbvP6LJ6TmN792hxyau0aazNq55cVl/rrQc9N1WC+p2/5wA0mGtlYy0Wkcl6Tn14koRtUk59WPGm5174PCJyNEopWKh7Z17pPgRL7WgbvfPCSAd9uzsY0lHm5RXesLNzq2NJmkM1yR7akbtxSkR1ECvsWdnD2U1hjpqGF+UkwtV3Xrjrqbnnlyoxi4Z2+oeN7/ttZFfR3AD/YEaecaSLpKVROMLzdqCVY02rBw/eP3O2HPiXjTW7hH3dSOlIoEN9DlKKxmLK4fUT9Cpabfn3mxafv3ko6jJODdcVtY9j87H3qvZhCaCHOgPlFZ6JOkY6k7WH68dv+nO46svPGtqY8dr3ywad62vf9EZdS9eUALhIsgzlnSRrKTL2zYa213WB+54IPJztW8WjbXwPRNHE92LbdeAMFEjz1jSMdRpZj+2O+GGmZbAYCPIMxa1M3xUnTnN7Md2J9ww0xIYbLzszEnal4vtvCiNulfcUroA+hcvO/tMpy8XGwP81ht3tfya+ns1LqVbWajqA3c8oJknntVHx3am/WMByAE98oBEbhCxwXTu2Ru1cCrZ9mpxwyNNSvRNAUB+6JEPgKiRLvX7ZlYWqhr/l+P6yL89HBvscS84feX6BDkQHoK8Q3lsZZZklEljsDeOF2+2JyejWIAwMWqlA1HT8Mc/d1y7PvLldRtAZKmTUSaNO/OM792h6Mn4jGIBQkWQdyCyxLHkWqgupl5fpZlO1/Ou72mP7S7r3ZdvXRfmrBcOhIsg70CSEkS7e1TGbelWb2x3WZuK7f8va+xpf3Rsp269cVfLse4AwkCNvAPN6sz1ktack667MjVbWbcLUKNiwbS4FL3hQz2m4wODgx55B+I2Pm6UtOacdNf5yem5NSHdqDxS0uQ7Lo3saSfp8QMIEz3yDjRO5hnZVNSPf3I6duuzVpKuhdKsh1+7X1RPu5OVFgGEgyDvUGNgphmOmHTFxLjzCmZNa9ydrrQIIAwEeUbS1JyjtnSL6tGP792h8c8dX1NeKRZMk++4tOm9Wf0QGGzUyPvE2XWjUUZKxfgedmOJPMEKC6x+CAw2euQZqZVWKgtVFcy05J5oZcGo9VN+evpM5LmT03PrRq0snvGWJZKkPX4AYaJHnoH6mZ6SVrdhSzIxKOmIFanzEknSNdIBhIkeeQaiwrim1UvFuBCOe/mZ5KVoFMaNA4OLHnkGWvWIoz5fG9cdV+K2lXPqtbszEIDhQI88A61mev5Mqag9E0djd7WPErWsLDvdA4jCxhIZiHphWVPcYJKp6YzMOCbpsYlrMmghgEEQt7EEpZUM1L9MlJYn6EjLLxXPPXtjRyEuMTwQQDKUVjIS9zJx+/4jHV2P2jeApAjyLku6UqIknbepmHjvzTx2KALQn1IFuZlNSnqbpBclfUfSH7r7QgbtGhhR0+qjmEmzf3lVomuyCBaAemlr5HdL+mV3/xVJ35J0IH2TBlCCEnk775zbmUQEYPCl6pG7+5frPrxP0jvSNScsScobUdPqo5TbeLHJIlgA6mU5auWPJP17htfra1EbMEdNx08Sru2+2GQRLAD1Wga5mX3FzB6K+HVd3TkfknRa0meaXGefmc2Y2cz8/Hw2rc9R0vJGXLgWzDpe94QZngDqtSytuPubm33ezP5A0rWS3uRNZhe5+yFJh6TlCUHtNbP/JC1vxK08mGbRKmZ4AqiXdtTK1ZI+KOm33P1UNk0KQ9IFrLoVuiyCBaAm7Tjyv5P0Mkl32/Jsxvvc/U9TtyoA7azxTegC6Ka0o1Z+MauGhIbyBoB+wczOFOhpA+gHBHmPMbUeQNYI8h5iaj2AbiDIe6B+Y+ZGrbaCA4BWCPIua7bpRA1T6wGkwcYSXdZsY+YaptYDSIMg77JWvW2m1gNIiyDvsma97U7WWQGARgR5l8UtcPXxG3fp3v1XEuIAUuNlZ5cxAxRAtxHkPcAMUADdRGkFAAJHkANA4AhyAAgcQQ4AgSPIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAEjiAHgMBlEuRmdpOZuZmdn8X1AADJpQ5yM7tI0lWSnkzfHABAu7Lokd8q6YOSPINrAQDalCrIzew6SRV3P57g3H1mNmNmM/Pz82luCwCos7HVCWb2FUk/H/GpD0n6Cy2XVVpy90OSDknS6OgovXcAyEjLIHf3N0cdN7OdkrZLOm5mknShpG+Y2evd/fuZthIAEKtlkMdx9xOSfq72sZk9LmnU3X+YQbsAAAkxjhwAAtdxj7yRu2/L6loAgOTokQNA4DLrkXfb1GxFk9NzOrlQ1ZaRksb37tDY7nLezQKA3AUR5FOzFR04fELVxSVJUmWhqgOHT0gSYQ5g6AVRWpmcnlsN8Zrq4pImp+dyahEA9I8ggvzkQrWt4wAwTIII8i0jpbaOA8AwCSLIx/fuUKlYWHOsVCxofO+OnFoEAP0jiJedtReajFoBgPWCCHJpOcwJbgBYL4jSCgAgHkEOAIEjyAEgcAQ5AASOIAeAwJl773ddM7N5SU/0/MbRzpfEZhg8hxqewzKew7J+ew6vcvfNjQdzCfJ+YmYz7j6adzvyxnNYxnNYxnNYFspzoLQCAIEjyAEgcAS5dCjvBvQJnsMynsMynsOyIJ7D0NfIASB09MgBIHAEOQAEjiCXZGa3mFnFzB5Y+fXWvNvUS2Z2tZnNmdm3zWx/3u3Ji5k9bmYnVv4OzOTdnl4xs0+a2TNm9lDdsVeY2d1m9j8r/z0vzzZ2W8wzCCYXCPKX3Oruu1Z+fSnvxvSKmRUkfULSWyRdIuldZnZJvq3K1RUrfwf6fuxwhj4l6eqGY/slfdXdXyPpqysfD7JPaf0zkALJBYIcr5f0bXf/rru/KOmzkq7LuU3oIXf/mqRnGw5fJ+nTK7//tKSxXrap12KeQTAI8pe818weXPkRa6B/jGxQlvRU3cffWzk2jFzSl83smJnty7sxOXuluz+98vvvS3plno3JURC5MDRBbmZfMbOHIn5dJ+nvJb1a0i5JT0v66zzbitz8hru/Tstlpj8zs9/Mu0H9wJfHKA/jOOVgciGYrd7Scvc3JznPzP5B0he73Jx+UpF0Ud3HF64cGzruXln57zNm9gUtl52+lm+rcvMDM7vA3Z82swskPZN3g3rN3X9Q+32/58LQ9MibWfmLWvN2SQ/FnTuAvi7pNWa23czOkvROSXfl3KaeM7NzzOzltd9LukrD9feg0V2S3rPy+/dI+tcc25KLkHJhaHrkLfyVme3S8o+Pj0v6k1xb00PuftrM3itpWlJB0ifd/eGcm5WHV0r6gplJy/8u/tnd/yPfJvWGmd0u6Y2Szjez70m6WdKEpDvN7I+1vOT07+TXwu6LeQZvDCUXmKIPAIGjtAIAgSPIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOD+Hxb2/WwziqyJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rgr = SparseGPRSolver(train_train_kernel, regularizer = 1)\n",
    "rgr.fit(train_train_kernel, train_energies)\n",
    "predictions = rgr.predict(train_test_kernel.T)[:, 0]\n",
    "print(\"rmse: \", get_rmse(predictions, test_energies))\n",
    "#print(predictions - test_energies)\n",
    "plt.plot(test_energies, predictions, 'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2476899b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 50/50 [00:55<00:00,  1.11s/it]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbgUlEQVR4nO3de3Sc9Z3f8fdX95st2bJlfEU2BmMHbIwNhBAKIWTrhFtLIIVNy5ZbQk7J6fbknJZttmW3h+3uSS/nbEhCCoSQ7iEEwpLEBid02w0HQgyxfBM2DsGWLVm2jG3dJVu3mW//mMHIQsIzkmae55n5vM6ZY80zj+b5aDSej57n91zM3REREUlFQdABREQkOlQaIiKSMpWGiIikTKUhIiIpU2mIiEjKVBoiIpKyoqADZNqcOXO8vr4+6BgiIpGybdu2E+4+d+z0nC+N+vp6Ghoago4hIhIpZtY83nRtnhIRkZSpNEREJGUqDRERSZlKQ0REUqbSEBGRlKk0REQkZTm/y62ISKYMx+L0DYwwMBJjeMQZjscZiTnDsTjDsTixuBN3iLsTd8dPfw3uzukLUzg4icdPTxr99STzXbtiLsWF07tuoNIQERmj+9QwhzpO0tJxkub2xL+HOk7S0T9E/9AIfQMj9A6OMDQSDzrqx9r18B9RXa7SEBGZViOxOG8d6OClxjb+7973Od47eMbjsytLWDy7gvnVZVSVFVFVOupWVkRZcSHFhQUUFxpFBYl/iwsLKCwwCgsMMygwS97ADCAx3ZLLMDOMDx5LTuPDO6Onp6qypDD9bzoLlYaI5KVY3HnrQDsvN7bxq91Hae8foqKkkOsurOPihdUsmV3BktoKlsyuYEZZcdBxQ0OlISJ556XGI/zFxnc40TdIeXEhn11Zx42r53PtijrKiqf/r/NcotIQkbwxEovzrVfe5fHXmrhkcQ3/5ZZP8JkVdZRnYDNOrlJpiEhe6Ogf4uvPbueNfe3cdeW5/PkNqygp0lEH6VJpiEjO2324m6/+3TaO9w3y325bze3rFwcdKbJUGiKS0/5+Wyv/8WdvU1tZwgsPXMnqRTVBR4o0lYaI5KwnX2/ikZf3cuWyWr7zx2uprSoNOlLkqTREJCe1dp7kv/+fd7l+ZR3f/5frKJrmI6PzlV5FEclJj7y0F8P4y1suUmFMI72SIpJzXvvDcX615ygPXrechTXlQcfJKSoNEckpQyNx/mLjHuprK7jv6qVBx8k5GtMQkZzyg98coOlEPz+8+zJKi3TQ3nSLZGmYWSXwPWAIeNXdnwk4koiEQFv3KR79x/e4fuU8PrOiLug4OWnSm6fMrMzMfmdmu8xsj5n95RSe6ykzO2Zmu8d5bIOZvWtm+8zsoeTkW4EX3P1+4ObJLldEcstfvbyXWNx5+KZVQUfJWVMZ0xgErnP3NcAlwAYz++ToGcyszsxmjJm2fJznehrYMHaimRUC3wU+D6wC7jSzVcAi4FByttgUfgYRyRG/3X+Clxrb+Nq157F4dkXQcXLWpEvDE/qSd4uTt7EXmLoG+LmZlQKY2f3Ao+M812tAxziLuRzY5+5N7j4E/AS4BWglURxT+hlEJDcMx+I8/Is9LJ5dzgPXnBd0nJw2pQ9cMys0s53AMeAf3P2t0Y+7+0+BV4DnzOzLwD3A7WksYiEfrlFAoiwWAi8CXzSzx4BNE2S7ycwe7+7uTmNxIhJFP/rtQd471sd/vvETOrV5hk2pNNw95u6XkPir/3Izu2iceb4FDACPATePWjuZynL73f1ud//aRIPg7r7J3b9SXV091cWJSIjF486Trx/gquW1XL9Sg9+ZNi2bdty9C/g1449LXA1cBPwMeDjNpz4MjD4d5aLkNBERALa3dHK0Z4Db1y3GJnNNVEnLVPaemmtmNcmvy4HPAb8fM89a4HES4xB3A7Vm9kgai9kKnG9mS82sBLgD2DjZzCKSe15qbKOkqIDPai0jK6aypjEf+LWZNZL4cP8Hd39pzDwVwJfcfb+7x4G7gOaxT2RmzwJbgBVm1mpm9wK4+wjwIIlxkb3A8+6+ZwqZRSSHxOPOL3e3ce0Fc3Ud7yyZ9MF97t4IrD3LPG+MuT8MPDHOfHd+zHNsBjZPMqaI5LBtLZ283zPIDavnBx0lb2h3VRGJrJdPb5qaF3SUvKHSEJFIisWdzW+38ZkVc6kqjeQZkSJJpSEikdRwsINjvYPcsHpB0FHyikpDRCLp5bfbKC0q4LMXaq+pbFJpiEjkxOLOL3cf5boL66jUpqmsUmmISORsPdjB8V7tNRUElYaIRM7LjW2UFRdwnTZNZZ1KQ0QiJZY8oO+6C+uoKNGmqWxTaYhIpLx1oJ0TfUPccLH2mgqCSkNEIuXlxjbKiwv5zIVzg46Sl1QaIhIZI7E4r+w5ynUrtWkqKCoNEYmM3x3o4ETfEDderL2mgqLSEJHIeOntNipKCrl2hfaaCopKQ0QiIR53Xkke0Fdeoku6BkWlISKRsO94H+39Q1rLCJhKQ0QiYevBDgDWnzsr4CT5TaUhIpHQcLCTOVWlnFtbEXSUvKbSEJFIaGju4LL6WZhZ0FHymkpDRELvaPcAhzpOsU6bpgKn0hCR0GtoToxnXFY/O+AkotIQkdBrONhJeXEhqxbMDDpK3lNpiEjoNTR3sHZJDcWF+sgKmn4DIhJqfYMjvHOkR7vahoRKQ0RCbUdLJ3GH9RrPCAWVhoiEWsPBTgoM1i6pCTqKoNIQkZBraO7gwnNmMqOsOOgogkpDREJsJBZnR0sXl9VrPCMsVBoiElp723o5ORTTeEaIqDREJLROn6RQaxqhodIQkdBqaO5gYU0586vLg44iSSoNEQkld2frwU6NZ4SMSkNEQulQxymO9w5qPCNkVBoiEkoazwgnlYaIhFJDcwczyoq4oG5G0FFkFJWGiIRSw8FO1p87i4ICXXQpTFQaIhI6nf1DvHesT+MZIVQUdIDJMLNK4HvAEPCquz8TcCQRmUbbmjsBdGbbEJr0moaZLTazX5vZO2a2x8z+7RSe6ykzO2Zmu8d5bIOZvWtm+8zsoeTkW4EX3P1+4ObJLldEwqmhuZPiQmPN4pqgo8gYU9k8NQJ8w91XAZ8E/o2ZrRo9g5nVmdmMMdOWj/NcTwMbxk40s0Lgu8DngVXAncllLAIOJWeLTeFnEJEQajjYwcULqykrLgw6iowx6dJw9zZ33578uhfYCywcM9s1wM/NrBTAzO4HHh3nuV4DOsZZzOXAPndvcvch4CfALUArieKY0s8gIuEzMByjsbVb4xkhNS0fuGZWD6wF3ho93d1/CrwCPGdmXwbuAW5P46kX8uEaBSTKYiHwIvBFM3sM2DRBppvM7PHu7u40FiciQdtzpJuhWJx1Gs8IpSmXhplVAX8P/Km794x93N2/BQwAjwE3u3vfVJfp7v3ufre7f22iQXB33+TuX6murp7q4kQki3YdSvyhd4nGM0JpSqVhZsUkCuMZd39xgnmuBi4CfgY8nOYiDgOLR91flJwmIjmqsbWLeTNLmTezLOgoMo6p7D1lwA+Ave7+PyeYZy3wOIlxiLuBWjN7JI3FbAXON7OlZlYC3AFsnGxmEQm/xsPdrF5UE3QMmcBU1jSuAv4VcJ2Z7UzevjBmngrgS+6+393jwF1A89gnMrNngS3ACjNrNbN7Adx9BHiQxLjIXuB5d98zhcwiEmI9A8M0He9n9UJtVg6rSR/c5+6/AT72+H53f2PM/WHgiXHmu/NjnmMzsHmSMUUkQna3JsYzVms8I7S0u6qIhMauD0pDaxqhpdIQkdB4+3AXS2ZXMKuyJOgoMgGVhoiExq5D3Vy8SGsZYabSEJFQaO8b5HDXKdaoNEJNpSEiodD4wXiGdrcNNZWGiIRCY2s3ZnCRBsFDTaUhIqHQ2NrFeXOrqCqN5GV+8oZKQ0QC5+7sau1mtcYzQk+lISKBa+se4ETfIGs0nhF6Kg0RCdyHg+Ba0wg7lYaIBK6xtYuiAmPl/JlBR5GzUGmISOAaW7tZcc4MXd41AlQaIhIod6extUvHZ0SESkNEAtXcfpKegREdCR4RKg0RCdSu1i5AR4JHhUpDRALV2NpNaVEB58+rCjqKpEClISKBamzt4hMLZlJcqI+jKNBvSUQCE4s7uw/3aNNUhKg0RCQw+471cWo4xprFGgSPCpWGiARGg+DRo9IQkcA0tnYxo7SIpbWVQUeRFKk0RCQwb7d2c9HCagoKLOgokiKVhogEYmgkzt62XlZrPCNSVBoiEojfH+1hKBbX6dAjRqUhIoHYlTwd+sW6vGukqDREJBA7WjqZU1XKolnlQUeRNKg0RCQQO1u6WLukBjMNgkeJSkNEsq6zf4imE/2sXVITdBRJk0pDRLJuZ/KgvkuXzAo2iKRNpSEiWbejuZMC0zXBo0ilISJZt+NQFxeeM5OKkqKgo0iaVBoiklXxuJ8eBJfoUWmISFbtP95H7+AIazWeEUkqDRHJqh0tXQBa04golYaIZNWOQ51UlxfrzLYRpdIQkaza3pwYz9CZbaNJpSEiWdM7MMwfjvWydrHGM6JKpSEiWdPY2o27xjOiTKUhIlmzo6UTgDWLa4INIpOm0hCRrNnR0sXyuiqqy4uDjiKTpNIQkaxwd3Yc6mKt1jIiTaUhIlnR3H6Sjv4hHdQXcSoNEcmKHYcS4xkaBI82lYaIZMWOli4qSwq5YN6MoKPIFKg0RCQrdrR0sWZxDYU6qC/SVBoiknGnhmLsbevRpqkcoNIQkYzbfaSbkbjrSPAcoNIQkYzb3pwYBL9EaxqRp9IQkYzb0dLFktkVzKkqDTqKTJFKQ0Qyyt3Z3tKp8YwcodIQkYxq6x7gWO8gl+qgvpyg0hCRjNKV+nKLSkNEMmp7SyelRQVceM7MoKPINFBpiEhGbW/p5OKF1ZQU6eMmF+i3KCIZMzAcY/fhbtbVazwjV6g0RCRj3j7czXDMWadB8Jyh0hCRjNmWPKhv3bkqjVyh0hCRjNnW3MnSOZXU6qC+nKHSEJGMcHe2N3fq+Iwco9IQkYw42H6S9v4h1msQPKeoNEQkIzSekZtUGiKSEduaO5lRVsTyuVVBR5FppNIQkYzY1tzBpUtmUaAr9eUUlYaITLvuU8P84f0+1mvTVM5RaYjItNvRovGMXKXSEJFpt725k8ICY83imqCjyDRTaYjItGto7mTl/BlUlhYFHUWmmUpDRKbVSCzOzkNdOt9UjlJpiMi0+v3RXk4OxbhU4xk5SaUhItNquwbBc1qkNjiaWSXwPWAIeNXdnwk4koiM0XCwk3NmlrGwpjzoKJIBga9pmNlTZnbMzHaPmb7BzN41s31m9lBy8q3AC+5+P3Bz1sOKyFlta+5k3bmzMNNBfbko8NIAngY2jJ5gZoXAd4HPA6uAO81sFbAIOJScLZbFjCKSgqPdAxzuOqXxjBwWeGm4+2tAx5jJlwP73L3J3YeAnwC3AK0kigM+JruZfcXMGsys4fjx45mILSLj+GA8Q0eC567AS2MCC/lwjQISZbEQeBH4opk9Bmya6Jvd/XF3X+/u6+fOnZvZpCJyWsPBTsqKC1i1YGbQUSRDIjUQ7u79wN1B5xCR8W1r6WT1ohqKC8P696hMVVh/s4eBxaPuL0pOE5GQGhiOsedwt3a1zXFhLY2twPlmttTMSoA7gI0BZxKRj9HY2s1I3DWekeMCLw0zexbYAqwws1Yzu9fdR4AHgVeAvcDz7r4nyJwi8vEamhP7s6zV6UNyWuBjGu5+5wTTNwObsxxHRCZpe3Mny+ZWMruyJOgokkGBr2mISPTF487Wg53aNJUHVBoiMmV7j/bQfWqYK8+rDTqKZJhKQ0SmbMv+dgA+uUylketUGiIyZW82dVBfW8H8ap2kMNepNERkSmJx560D7do0lSdytjTM7CYze7y7uzvoKCI57Z0jPfQOjGjTVJ7I2dJw903u/pXq6uqgo4jktDebEuMZV6o08kLOloaIZMeWpnaWza2kbmZZ0FEkC1QaIjJpI7E4Ww90aNNUHlFpiMik7TnSQ+/giDZN5RGVhohM2pbkeMYVy2YHnESyRaUhIpP2ZlM7y+uqqJuh8Yx8odIQkUkZTo5naNNUflFpiMik7D7cTf9QTIPgeUalISKTovGM/KTSEJFJebOpgwvmVTGnqjToKJJFKg0RSdtwLE7DQY1n5COVhoikrbG1i5Maz8hLKg0RSdubTYnrgV+h0sg7OVsaOsutSOZs2d/OhefM0PXA81DOlobOciuSGUMjcRqadb6pfJWzpSEimbGrtYuB4bguupSnVBoikpYt+9sxgyuW6viMfKTSEJG0vNnUzspzZlJTofGMfKTSEJGU9Q2OsK25U+MZeUylISIpe37rIQZH4tx8yYKgo0hAVBoikpJY3HnqjQNcVj+LSxbXBB1HAqLSEJGUvLLnKK2dp7j308uCjiIBUmmISEqeeL2Jc2sr+NyqeUFHkQCpNETkrLY1d7KjpYt7rlpKYYEFHUcCpNIQkbN68vUmqsuLuX39oqCjSMBUGiLysVraT/LKnqN8+YolVJQUBR1HAqbSEJGP9dQbBygsMP7kU/VBR5EQUGmIyIS6Tw7zfMMhblqzgHkzy4KOIyGQs6WhU6OLTN2zW1s4ORTjPu1mK0k5Wxo6NbrI1AyNxHn6jYNctbyWVQtmBh1HQiJnS0NEpmbz220c7Rngvqu1liEfUmmIyEe4O0+83sTyuiquOX9u0HEkRFQaInKGU0MxHt64hz1Herjv00sp0MF8Mop2uhaR03a0dPKN53fRdKKfu6+q57Z1OphPzqTSEBGGRuJ8+/+9x/de3cf86nJ+fN8VfGr5nKBjSQipNETy3LtHe/l3z+3knbYebl+3iP900ypmlhUHHUtCSqUhkodiceetA+1s3HmEF7cfZmZ5EU/ctV5nsJWzUmmI5Al3Z8+RHn6x8zCbdiV2p60sKeSfr13Iv9+wgtqq0qAjSgSoNCbw4vZWTvQNBh1DpsAYf68fG2eyJScWGFjyvo36usCMwgKS/yZuBWYUFxZQWlRASfL2wdeVJUXMLC9mRmlR4HsfHTjRz6ZdR/jFzsPsP95PcaFxzQV1fPOGlVy/ch7lJYWB5pNoUWlM4Kk3DrD7cE/QMSTizKDqgwIpK6K2qoT51eUsqClnQXUZ82vKWVhTxoKa8mk9g+yRrlO81HiETbvaePtw4lQ6l9fP5p5PL+ULF81nVmXJtC1L8otKYwI//eqniLsHHUMmaaLfnI/6nfrpaR9OcJy4J+bz5GPuiWkxd+JxJxb3018PxeIMjsQZGnUbHInTPzRCz6lhegdG6BkYpudU4t8TfYP85r0TvN87wNi31zkzy1g2t5KlcxK38+ZWcW5tBXNmlDKjtOj02tBYJ4dGaOk4SXP7SQ6c6Ocf9x7jdwc7AFi9qJo/v2ElN6yez/zq8km/niIfUGlMQKvskknDsTjv9wxwpGuAtu5THOo4SdOJfpqOJzYl9QyMnDF/SWEBsyqLmV1ZyuzKYmrKSzjeO8jB9n6O9Z65GfX8uiq+8bkLuGnNAurnVGbzx5I8EMnSMLNlwDeBane/Leg8IukqLixg0awKFs2q+Mhj7k7nyWGajvfR3H6Sjv4h2vuH6OgfpKN/mI7+Qdq6e5hTVco1F8ylfk4lS2ZXUF9byZLaCqrLtbusZE5KpWFmNcCTwEUk1urvcfct6S7MzJ4CbgSOuftFYx7bAPwtUAg86e5/M9HzuHsTcK+ZvZBuBpGwMzNmV5Ywu3I26+tnBx1H5Ayprmn8LfArd7/NzEqAM/48MrM64JS7946attzd9415nqeB7wD/e8z3FwLfBT4HtAJbzWwjiQL56zHPcY+7H0sxt4iITKOzloaZVQP/BPjXAO4+BAyNme0a4AEz+4K7D5rZ/cCtwOdHz+Tur5lZ/TiLuRzYl1yDwMx+Atzi7n9NYs0kbWZ2E3DT8uXLJ/PtIiIyjlTOcrsUOA780Mx2mNmTZnbG6Jq7/xR4BXjOzL4M3APcnkaOhcChUfdbk9PGZWa1ZvZ9YK2Z/dl48+giTCIi0y+V0igCLgUec/e1QD/w0NiZ3P1bwADwGHCzu/dNZ9Axy2p39wfc/bzk2oiIiGRBKqXRCrS6+1vJ+y+QKJEzmNnVJAbKfwY8nGaOw8DiUfcXJaeJiEiInLU03P0ocMjMViQnfRZ4Z/Q8ZrYWeBy4BbgbqDWzR9LIsRU438yWJgfa7wA2pvH9IiKSBaleue/rwDNm1ghcAvzXMY9XAF9y9/3uHgfuAprHPomZPQtsAVaYWauZ3Qvg7iPAgyTGRfYCz7v7nkn8PCIikkHmOX6qjPXr13tDQ0PQMUREIsXMtrn7+o9Mz/XSMLNe4N2gc6RpDnAi6BCToNzZE8XMoNzZNNXM57r73LETI3kakTS9O15bhpmZNUQtMyh3NkUxMyh3NmUqc6pjGiIiIioNERFJXT6UxuNBB5iEKGYG5c6mKGYG5c6mjGTO+YFwERGZPvmwpiEiItNEpSEiIilTaYiISMryrjTMbJmZ/WD0Vf/MrMDM/srMHjWzPwky30QmyL3SzL5vZi+Y2deCzDeeCTJXmtmPzOyJ5Gn0Q8vMrk6+vk+a2W+DzpOO5OvcYGaTuh5NEMZ7v4SRmf2z5Pv3OTP7o6DzpMrMrjWz15Pv6Wsn+zyRKg0ze8rMjpnZ7jHTN5jZu2a2z8w+ctr20dy9yd3vHTP5FhJn1h0mcVbfaZWp3O6+190fAL4EXBWFzCQuzvWCu98P3DydmcfknI78rydf35eAH2Uq65h8U86d9B+A5zOT8qMy+H7JinTyu/vPk+/fB4B/EUTeUfnSed0d6APKmMrnnLtH5kbiCoKXArtHTSsE9gPLgBJgF7AKuJjEf/bRt7pR3/fCqK8fAr46dnrYcyfv3wz8EvjjKGQG/gy4JPn1jyPyXnkemBGV9ziJyybfQeJqmzdGJfdE7/Gw5R/1+P8ALs121im87gXJx+cBz0x2mZE6jYiPf7nY6bhUbCsfXsI2Nh1ZR8tgbtx9I7DRzF4GfjxNkTP9Wi8CdpLBNd3pym9mS4Bud+/NVNbRpiN3ctNDJYkPilNmttkTZ58Ode4gpZPfzPYCfwP80t23ZzfpmdJ83T+4pEUnUDrZZUZq89QEpuNSsS8C/9TMHgVey1jSM005d3Ib5bfN7H8BmzOaNmG6XusvmtljwKaMJR1fWvmT7gV+mLFEqUkrt7t/093/lMQfEU9kujA+xrRfxjnLJsr/deB64DYzeyCIYGcxbm4zuzX5WfF3wHcm++SRWtOYDu7eTmJb5OhpJ0l8OITWBLlfBV4NIk8qJsjcT+JCXZHg7ulehTI03P3poDOkY7z3Sxi5+7eBbwedI13u/iKJP9qmJBfWNKJ6qdgo5o5i5tGiml+5gxHV/BnNnQulEdVLxUYxdxQzjxbV/ModjKjmz2zuIEf+J7GnwLNAGx/uGntvcvoXgD+Q2GPgm0HnzIXcUcycC/mVW/nDnlsnLBQRkZTlwuYpERHJEpWGiIikTKUhIiIpU2mIiEjKVBoiIpIylYaIiKRMpSEiIilTaYiISMpUGiIikrL/D51KcHJBHPMdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha_grid = np.logspace(5, -15, 50)\n",
    "rmse = []\n",
    "for alpha in tqdm.tqdm(alpha_grid):\n",
    "    rgr = SparseGPRSolver(train_train_kernel, regularizer = alpha)\n",
    "    rgr.fit(train_train_kernel, train_energies)\n",
    "    predictions = rgr.predict(train_test_kernel.T)[:, 0]\n",
    "    rmse.append(get_rmse(predictions, test_energies))\n",
    "plt.plot(alpha_grid, rmse)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0a5241f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse:  0.6174583594600753\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAASgUlEQVR4nO3df2zc9X3H8dc75ugubKrDYkWJSxYaIUftaB1qMbpUFWEtpkVdTDuVsnXjj4rwR/lrlaVkQgNNk4jqUfpPhUQ1BBKQDongRiyq2wISEhsIZ6ZNqLCgNEAuKTFLjar1thrz3h93X+d8vvN97+57973PfZ8PyfL5w1ff++jqvvzN+/PL3F0AgHBtSLsDAID2EOQAEDiCHAACR5ADQOAIcgAI3EVpvOnmzZt9x44dabw1AATr+PHj77r7UHV7KkG+Y8cOzc7OpvHWABAsM3uzVjulFQAIHEEOAIEjyAEgcAQ5AASOIAeAwKUyawUAsmR6rqCpmXmdWSxq22Bek+Mjmtg9nNj9CXIA6KDpuYIOHjmh4tKyJKmwWNTBIyckKbEwp7QCAB00NTO/EuKR4tKypmbmE3sPghwAOujMYrGp9lYQ5ADQQdsG8021t4IgB4AOmhwfUT43sKrNJO3dtWbLlJYR5ADQQRO7h/WVTw3LKtpc0hPHC5qeKyTyHgQ5AHTYs68uqPp05CQHPAlyAOiwTg94EuQA0GGdHvAkyAGgw2oNeOZzA5ocH0nk/qzsBIAOi1ZwdmqZPkEOAF0wsXs40f1VKlFaAYDAEeQAEDiCHAACFzvIzexBMztnZicr2u42s4KZvVz++mJnugkAqKeZJ/KHJN1Qo/0+dx8tfx1LplsAgLhiB7m7PyfpfAf7AgBoQRI18jvM7Ofl0sumeheZ2X4zmzWz2YWFhQTeFgAgtR/k90vaKWlU0llJ99a70N0fcPcxdx8bGkpu+0YAyLq2gtzd33H3ZXf/QNL3JV2dTLcAAHG1FeRmtrXix5sknax3LQCgM2Iv0Tezw5KulbTZzE5LukvStWY2qtI+6ack3Z58FwEA64kd5O5+S43mf02wLwCAFrCyEwACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAEjiAHgMDF3o8cACpNzxU0NTOvM4tFbRvMa3J8RBO7h9PuViYR5ACaNj1X0MEjJ1RcWpYkFRaLOnjkhCQR5imgtAKgaVMz8yshHikuLWtqZj6lHmUbQQ6gaWcWi021o7MIcgBN2zaYb6odnUWQA2ja5PiI8rmBVW353IAmx0dS6lG2MdgJoGnRgCazVnoDQQ6gJRO7hwnuHkFpBQACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACFzsIDezB83snJmdrGi71Mx+Ymavlb9v6kw3AQD1NPNE/pCkG6raDkh62t2vkPR0+WcAQBfFDnJ3f07S+armfZIeLr9+WNJEMt0CAMTVbo18i7ufLb/+taQt9S40s/1mNmtmswsLC22+LQAgkthgp7u7JF/nvz/g7mPuPjY0NJTU2wJA5rUb5O+Y2VZJKn8/136XAADNaDfIj0q6tfz6Vkk/bPN+AIAmNTP98LCk/5Q0Ymanzewbkg5J+ryZvSbpc+WfAQBdFPvMTne/pc5/+ouE+gIAaAErOwEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAELvYSfQCdMz1X0NTMvM4sFrVtMK/J8RFN7B5Ou1sIBEEOpGx6rqCDR06ouLQsSSosFnXwyAlJIswRC6UVIGVTM/MrIR4pLi1ramY+pR4hNAQ5kLIzi8Wm2oFqBDmQsm2D+abagWoEOZCyyfER5XMDq9ryuQFNjo+k1COEhsFOIGXRgCazVtAqghzoARO7hwlutIzSCgAEjiAHgMAR5AAQOIIcAAJHkANA4AhyAAgcQQ4AgSPIASBwLAhCprEPOPoBQY7MYh9w9AtKK8gs9gFHvyDIkVnsA45+QZAjs9gHHP2CIEdmsQ84+gWDncgs9gFHvyDIkWnsA45+kEiQm9kpSb+VtCzpfXcfS+K+AIDGknwi3+vu7yZ4PwBADAx2AkDgkgpyl/RjMztuZvtrXWBm+81s1sxmFxYWEnpbAEBSQf4Zd79K0hckfdPMPlt9gbs/4O5j7j42NDSU0NsCABIJcncvlL+fk/SkpKuTuC8AoLG2BzvN7BJJG9z9t+XX10v6p7Z7BiSA3Q2RBUnMWtki6Ukzi+73mLv/KIH7Am1hd0NkRdtB7u5vSPpkAn0BErXe7oYEOfoJ0w/Rt9jdEFlBkKNvsbshsoK9VhC8O6dP6PCLb2vZfaVteDCvvbuG9MTxwqryCrsboh8R5AhO5UyUfG6Dfrf0wZprCotFPXG8oK98aljPvrrArBX0NYIcQameiVIrxCPFpWU9++qCnj9wXbe6B6SCGjmCUmsmynoKDGwiA3giR8+rLKV448tXGSitbwD6GkGOnlZdSmlW5QAo0K8oraCnNVtKqTbMVENkAEGOnhZ38c7G3AYOUkZmEeToaXEX71x80YDu+fKVGh7My1R6Er/ny1cy1RCZQI0cPW1yfCRWjfy94hIHKSOzCHKkqtE2s9HrqZn5dacSsuweWUZpBamJZqQUytMKo21mp+cKq66b2D2s5w9cV3fg0iRq4cg0ghypWW+b2Vomx0fWDGiapL+5ZjslFWQapRWkpt6MlHollMoyC3unABcQ5Oiqypr4BrOaC3asfF2tgGZAE1iLIEfXVK/SrLfq0qW6p/hwBiewFkGOjovCt5kNrGqVXTiDE6iNwU501PRcQX//+MtN70JYazphs4OjQFbwRI7EVJ7UM2CmW/7sMh05flofNLlvVb3phJzBCdTGEzkScef0CT3ywlsrde9ldz3ywlvrHvxQy3rTCTmDE6iNJ3K0pPrpO4ntYocbDF7WWq7PxlgAQY4WRE/fkSRC3KSGR7IxjxyojSBH0w6/+Hbi94xbHmEeObAWNXI0rZ0ncDNpQ9Xpa5RHgPYQ5GhK9YZWzRgezOtX99yo73x1lH3DgQRRWkFs0YKcVkXTBCmPAMniiRyxtXt+JtMEgc4gyBFbnIU3G3OlX6mqMjh1cKCDKK1kXPUmVHt3DenZVxfWTO+bnivU3a2w0qZLPqRfHLiOza2ALjJPYA5ws8bGxnx2drbr74vVqjehqsUk/fnOS/Vfb70Xq6xikn516MbkOglghZkdd/ex6naeyDMsTs3bJT3/y/Ox70kdHOg+auQZ1uyOhHFQBwe6L5EgN7MbzGzezF43swNJ3BOdded069MIAfSWtksrZjYg6XuSPi/ptKSXzOyou/+i3XujOfUGGGsNaD5asVdKkuqd7AOgc5KokV8t6XV3f0OSzOwHkvZJIsi7ZHquoLuPvqLF4tJKW3R6zuyb5/VvL72tpWVfaX+kQyEusTc4kIYkSivDkip3UTpdbkMXRDNPKkM8Ulxa1mMvvrUS4t3AYCfQfV2btWJm+yXtl6Tt27d36237VtxzMJs9nScuk3TRgK36I8GiHyAdSTyRFyRdVvHzR8ptq7j7A+4+5u5jQ0NDCbxtdkVP4Z2YdVLp69ds13dvHlU+N7CqPZ8b0H03j2rqrz7J5ldAD0jiifwlSVeY2eUqBfjXJP11AvdFHe3uedLIpo053fWlj68K5XqrNAluIH1tB7m7v29md0iakTQg6UF3f6XtnqGuTg0obtqY09w/Xr+mnd0Kgd6WSI3c3Y9JOpbEvdDYtsF8R8oqi79bO2AKoPexsjNAk+MjNevWe3Ze2tZ9mXEChIkgD0w0W6W4tKwBK20WGw00Pnrbp1u+r4nl9UCo2DQrINW7FS67y1Ra5DM1My+pFOqtlF1cDFwCoeKJvIdMzxW059AzuvzAv2vPoWfWnI9Za7ZKNIs7Wsm5d9fQmrJLHMOUVYBgEeQ9onJuuOtCMFeGeaPZKsWlZR1+8e01ZZevX7N93f+hWcgDhI0g7xG1nraLS8srJRMp3mBkdILPsvtKQP/zxJX6zs0XTq7ftDGnwXyOhTxAn6BG3iPqPW1Xtu/dNdTUhlfRH4JoHjhhDfQnnsh7RL2n7cr2p352tun7shsh0P8I8h5Rb254Ze261g6HjTA3HOh/lFZ6RFT2SPLkeQYxgWwgyFNU60Sf5w9cV/f6TRtz+k3MZfTDCfwhABAGgrzLKvcRN62dBy7VX5hz4ye2NhzszOcGmIUCZAw18i6q3ke8+syH6umG1Z59dWHd+2/amCPEgQwiyLsozj7i680yaTQDZePFFxHiQAZRWumAeqfZx9kDZb1ZJo22r2WqIZBNBHlC4tS+B8xWVl7WEm2AtefQMzUHKifHR1ZtmlWNqYZANhHkCajelbBe7btRiDca+Ixe3330lTVzyplqCGQXNfIExK19r7fDYNyBz4ndw3r5ruv13Yq9U9gvBcg2nsgTEKc2HdXK1yuNNHNf9k4BEOGJPAGNatPR6TsTu4d11fYPJ3ZfAJB4Im9Z5cyUwY05bZD0QZ1ro9N3pucK+o9fno91/9wGo+YNIBaCvAXVg5uNls1HtfGpmfk1tfB6/vAPmBMOIB6CvAVxBjcj+dyAdvxxXjsPHlt31kq1xZh7qgAAQd6CZhbeXLX9w3o+ZjmlEvVxAHEx2NmCuCE7PJjXC2/8pun7MyccQDMI8hbUOgSiWhTGccopuQHjDE0ALaO00oJaKywvuXhAuYENeq+4tGp/lW89/rO6YW5SIgdIAMg2grwN//f+hQmH//P7ZeVz0n03j64K5Ws+uqlmjXzPzkv16G2fbvge9TbgAoAIQd6iWjNXKk+tj5z679oDo9XttQJb0qppjnEOnwCQPQR5i+rNXKluj3Nd9bz0KLA/dNGGWH8sAGQbg50tqjdzpbo9znX1nu6rdziMsO84gEoEeYtqzVypNW0wznXNBjNzzAFUorTSoqi00WggMs519U7+2bQxp/9d+mDV0zpzzAFUI8jbkNRWsrW2t83nBnTXlz4uqfEfCwDZ1laQm9ndkm6TFB3v/g/ufqzdTvWTegOZ0oWn9UZP7QQ3gPUk8UR+n7v/SwL36UtxpylyUASAVjHY2WFxpykCQKuSeCK/w8z+TtKspG+5e/O7RPWJWot66g1kMvMEQFIaPpGb2U/N7GSNr32S7pe0U9KopLOS7l3nPvvNbNbMZhcWFupdFqyoFl5YLMp1oRa+d9dQrGmKANAq8yYOO1j3RmY7JD3l7n/a6NqxsTGfnZ1N5H17xZ5Dz9R88h4uP5kz8wRAu8zsuLuPVbe3O2tlq7ufLf94k6ST7dwvZOvVwhnIBNBJ7dbIv21moyqdL3xK0u3tdihU1MIBpKWtWSvu/rfufqW7f8Ld/7Li6Txz4i7ZB4CksbIzIXGX7ANA0oIJ8hAOWKAWDiANQQR5nGXuAJBVQazsXG+ZOwBkXRBBzjJ3AKgviCCPexoPAGRREEHO1D4AqC+IwU6m9gFAfUEEucTUPgCoJ4jSCgCgPoIcAAJHkANA4AhyAAgcQQ4AgUvshKCm3tRsQdKbXX/j2jZLejftTvQAPocSPocSPoeSXvsc/sTdh6obUwnyXmJms7WOTsoaPocSPocSPoeSUD4HSisAEDiCHAACR5BLD6TdgR7B51DC51DC51ASxOeQ+Ro5AISOJ3IACBxBDgCBI8glmdndZlYws5fLX19Mu0/dZGY3mNm8mb1uZgfS7k9azOyUmZ0o/w7Mpt2fbjGzB83snJmdrGi71Mx+Ymavlb9vSrOPnVbnMwgmFwjyC+5z99Hy17G0O9MtZjYg6XuSviDpY5JuMbOPpdurVO0t/w70/NzhBD0k6YaqtgOSnnb3KyQ9Xf65nz2ktZ+BFEguEOS4WtLr7v6Gu/9e0g8k7Uu5T+gid39O0vmq5n2SHi6/fljSRDf71G11PoNgEOQX3GFmPy//E6uv/xlZZVjS2xU/ny63ZZFL+rGZHTez/Wl3JmVb3P1s+fWvJW1JszMpCiIXMhPkZvZTMztZ42ufpPsl7ZQ0KumspHvT7CtS8xl3v0qlMtM3zeyzaXeoF3hpjnIW5ykHkwvBHPXWLnf/XJzrzOz7kp7qcHd6SUHSZRU/f6TcljnuXih/P2dmT6pUdnou3V6l5h0z2+ruZ81sq6RzaXeo29z9neh1r+dCZp7I11P+RY3cJOlkvWv70EuSrjCzy83sYklfk3Q05T51nZldYmZ/FL2WdL2y9XtQ7aikW8uvb5X0wxT7koqQciEzT+QNfNvMRlX65+MpSben2psucvf3zewOSTOSBiQ96O6vpNytNGyR9KSZSaX/Xzzm7j9Kt0vdYWaHJV0rabOZnZZ0l6RDkh43s2+otOX0V9PrYefV+QyuDSUXWKIPAIGjtAIAgSPIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOD+H/9qyQVLr4XSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rgr = SparseGPRSolver(train_train_kernel, regularizer = alpha_grid[np.argmin(rmse)])\n",
    "rgr.fit(train_train_kernel, train_energies)\n",
    "predictions = rgr.predict(train_test_kernel.T)[:, 0]\n",
    "print(\"rmse: \", get_rmse(predictions, test_energies))\n",
    "#print(predictions - test_energies)\n",
    "plt.plot(test_energies, predictions, 'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a015644e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('test_saving', [test_energies, predictions])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d26c1316",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}