{ "cells": [ { "cell_type": "markdown", "id": "47d10a9a", "metadata": {}, "source": [ "# Using feature selection to build an efficient potential" ] }, { "cell_type": "markdown", "id": "85a81755", "metadata": {}, "source": [ "In the previous examples (`MLIP_example` and `zundel_i-PI`) we have seen how to use the tools provided by `librascal` in order to build a working potential. We have already used a \"sparsification\" technique to select a basis set of environments for the fit, as well as to keep the computational cost of both fitting and evaluating the potential manageable.\n", "\n", "Now we will introduce a technique to further optimize the computational cost of the model, and it consists of applying the same sparsification technique along the _features_ (columns) of the feature matrix, in addition to the _samples_ (rows). This technique was introduced in Ref. [1] and has since been used in numerous applications. Feature selection and sparse feature computation are both available directly in librascal, and this example aims to show the user how to apply them in a realistic setting.\n", "\n", "A detail on the selection algorithm itself: `librascal` uses the FPS and CUR algorithms implemented in [`scikit-cosmo`](https://scikit-cosmo.readthedocs.io/en/latest/selection.html). You will therefore need to have `scikit-cosmo` installed (`pip install skcosmo`) to run this notebook. Note also that further selection algorithms, notably PCovCUR and PCovFPS, are implemented in `skcosmo` but not yet included directly in `librascal`.\n", "\n", "[1]: G. Imbalzano, A. Anelli, D. Giofré, S. Klees, J. Behler, and M. Ceriotti, Automatic Selection of Atomic Fingerprints and Reference Configurations for Machine-Learning Potentials, J. Chem. Phys. 148, 241730 (2018). [doi:10.1063/1.5024611](https://aip.scitation.org/doi/full/10.1063/1.5024611)\n", "\n", "----" ] }, { "cell_type": "markdown", "id": "cd8d9cdd", "metadata": {}, "source": [ "Let us begin with the Zundel cation example introduced in the previous notebook:" ] }, { "cell_type": "markdown", "id": "a298ec8b", "metadata": {}, "source": [ "## Initialization" ] }, { "cell_type": "code", "execution_count": 1, "id": "ee33dff6", "metadata": { "init_cell": true }, "outputs": [], "source": [ "%matplotlib inline\n", "from matplotlib import pylab as plt\n", "\n", "import ase\n", "from ase.io import read, write\n", "from ase.build import make_supercell\n", "from ase.visualize import view\n", "import numpy as np\n", "# If installed -- not essential, though\n", "try:\n", " from tqdm.notebook import tqdm\n", "except ImportError:\n", " tqdm = (lambda i, **kwargs: i)\n", "\n", "from time import time\n", "# Print some extra information during the selection process\n", "import logging\n", "logging.basicConfig(level=logging.INFO)\n", "\n", "from rascal.models import Kernel, train_gap_model, compute_KNM\n", "from rascal.representations import SphericalInvariants\n", "from rascal.utils import from_dict, to_dict, CURFilter, FPSFilter, dump_obj, load_obj" ] }, { "cell_type": "code", "execution_count": 2, "id": "c0320d69", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/Max/Work/codes/librascal/examples/i-PI/zundel\n" ] } ], "source": [ "cd i-PI/zundel/" ] }, { "cell_type": "code", "execution_count": 3, "id": "1a8da102", "metadata": {}, "outputs": [], "source": [ "# Load the first N structures of the zundel dataset\n", "N_dataset = 1000\n", "frames = read('zundel_dataset.xyz', index=':{}'.format(N_dataset))\n", "energies = np.loadtxt('zundel_energies.txt')[:N_dataset]\n", "forces = np.concatenate([frame.arrays['forces'] for frame in frames])" ] }, { "cell_type": "code", "execution_count": 4, "id": "2d198335", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(7000, 3)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "forces.shape" ] }, { "cell_type": "code", "execution_count": 5, "id": "0fea5828", "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# Number of structures to train the model with\n", "n_train = 800\n", "\n", "global_species = [1, 8]\n", "\n", "# Select randomly n structures for training the model\n", "ids = list(range(N_dataset))\n", "np.random.seed(10)\n", "np.random.shuffle(ids)\n", "\n", "train_ids = ids[:n_train]\n", "frames_train = [frames[ii] for ii in ids[:n_train]]\n", "e_train = np.array([energies[ii] for ii in ids[:n_train]])\n", "f_train = np.concatenate([frame.arrays['forces'] for frame in frames_train])\n", "# Test on the remaining 200 molecules\n", "test_ids = [int(i) for i in ids[n_train:]] \n", "frames_test = [frames[ii] for ii in test_ids]\n", "e_test = np.array([energies[ii] for ii in test_ids])\n", "f_test = np.concatenate([frame.arrays['forces'] for frame in frames_test])" ] }, { "cell_type": "markdown", "id": "bc1b367a", "metadata": {}, "source": [ "## Set hyperparameters" ] }, { "cell_type": "code", "execution_count": 6, "id": "ee7dcb0b", "metadata": {}, "outputs": [], "source": [ "# Atomic energy baseline - assuming all configurations have the same number of atoms\n", "atom_energy_baseline = np.mean(energies)/len(frames[0])\n", "energy_baseline = {int(species): atom_energy_baseline for species in global_species}" ] }, { "cell_type": "code", "execution_count": 7, "id": "f41ac68a", "metadata": {}, "outputs": [], "source": [ "# define the parameters of the spherical expansion\n", "hypers = dict(soap_type=\"PowerSpectrum\",\n", " interaction_cutoff=3.0, \n", " max_radial=8, \n", " max_angular=6, \n", " gaussian_sigma_constant=0.5,\n", " gaussian_sigma_type=\"Constant\",\n", " cutoff_function_type=\"RadialScaling\",\n", " cutoff_smooth_width=0.5,\n", " cutoff_function_parameters=\n", " dict(\n", " rate=1,\n", " scale=3.5,\n", " exponent=4\n", " ),\n", " radial_basis=\"GTO\",\n", " normalize=True,\n", " optimization=\n", " dict(\n", " Spline=dict(\n", " accuracy=1.0e-05\n", " )\n", " ),\n", " compute_gradients=False\n", " )\n", "\n", "\n", "soap = SphericalInvariants(**hypers)" ] }, { "cell_type": "markdown", "id": "ee0b243f", "metadata": {}, "source": [ "## Compute descriptors" ] }, { "cell_type": "code", "execution_count": 8, "id": "3dd2d30d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Execution: 0.30106329917907715 s\n" ] } ], "source": [ "managers = []\n", "for f in frames_train:\n", " f.wrap(eps=1e-18)\n", "\n", "start = time()\n", "managers = soap.transform(tqdm(frames_train))\n", "print (\"Execution: \", time()-start, \"s\")" ] }, { "cell_type": "markdown", "id": "4a2732fd", "metadata": {}, "source": [ "## Sample selection" ] }, { "cell_type": "markdown", "id": "30c2c3d5", "metadata": {}, "source": [ "Now we select the set of sparse _samples_ (environments). In effect, we are selecting an optimally diverse set of _rows_ of the feature matrix $X$." ] }, { "cell_type": "code", "execution_count": 9, "id": "60d5982b", "metadata": { "code_folding": [], "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:rascal.utils.filter:The number of pseudo points selected by central atom species is: {1: 50, 8: 100}\n", "INFO:rascal.utils.filter:Selecting species: 1\n", "INFO:rascal.utils.filter:Selecting species: 8\n" ] } ], "source": [ "# select the sparse points for the sparse kernel method with FPS on the whole training set\n", "n_sparse_env = {1:50, 8:100}\n", "sample_compressor = FPSFilter(soap, n_sparse_env, act_on='sample per species')\n", "X_sparse = sample_compressor.select_and_filter(managers)" ] }, { "cell_type": "markdown", "id": "8e254859", "metadata": {}, "source": [ "## Model assessment: Sparse samples only" ] }, { "cell_type": "markdown", "id": "51e67152", "metadata": {}, "source": [ "Let's build our potential and get some benchmarks." ] }, { "cell_type": "code", "execution_count": 10, "id": "55bc94ab", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Execution: 4.386488199234009 s\n" ] } ], "source": [ "zeta = 2\n", "\n", "start = time()\n", "hypers['compute_gradients'] = True\n", "soap_grads = SphericalInvariants(**hypers)\n", "kernel = Kernel(soap_grads, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')\n", "\n", "KNM = compute_KNM(tqdm(frames_train, leave=True, desc=\"Computing kernel matrix\"), X_sparse, kernel, soap_grads)\n", "\n", "model = train_gap_model(kernel, frames_train, KNM, X_sparse, e_train, energy_baseline, \n", " grad_train=-f_train, lambdas=[1e-12, 1e-12], jitter=1e-13)\n", "\n", "# save the model to a file in json format for future use\n", "dump_obj('zundel_model.json', model)\n", "np.savetxt('Structure_indices.txt', ids)\n", "print (\"Execution: \", time()-start, \"s\")" ] }, { "cell_type": "code", "execution_count": 11, "id": "4046787a", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# make predictions on the test set\n", "e_pred = []\n", "f_pred = []\n", "for f in tqdm(frames_test):\n", " positions = f.get_positions()\n", " f.wrap(eps=1e-18)\n", " m = soap_grads.transform(f)\n", " e_pred.append(model.predict(m)[0])\n", " f_pred.append(model.predict_forces(m))\n", "\n", "e_pred = np.array(e_pred)\n", "f_pred = np.concatenate(f_pred)" ] }, { "cell_type": "code", "execution_count": 12, "id": "065ae1f5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSE = 37.27416992356848 meV\n", "Test set stdev = 0.274800153468654 eV\n", "Relative RMSE = 13.564100839492538 %\n" ] } ], "source": [ "from rascal.utils import get_score\n", "\n", "score = get_score(e_pred, e_test)\n", "RMSE = score['RMSE']\n", "sigma_test = np.std(e_test)\n", "print(\"RMSE = \", RMSE*1000.0, \"meV\")\n", "print(\"Test set stdev = \", sigma_test, \" eV\")\n", "print(\"Relative RMSE = \", RMSE/sigma_test*100.0, \" %\")" ] }, { "cell_type": "markdown", "id": "fdaa07ef", "metadata": {}, "source": [ "## Now introduce feature selection" ] }, { "cell_type": "markdown", "id": "fecb1e3d", "metadata": {}, "source": [ "This potential uses quite a large number of features to make its predictions:" ] }, { "cell_type": "code", "execution_count": 13, "id": "b10c08f5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1344" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "soap.get_num_coefficients(n_species=2)" ] }, { "cell_type": "markdown", "id": "e761f4b0", "metadata": {}, "source": [ "Since the cost to compute the kernel (and hence evaluate the model) scales asymptotically linearly with the number of features [2], it would be nice if we could get away with using fewer components of the feature vector.\n", "\n", "[2]: F. Musil, M. Veit, A. Goscinski, G. Fraux, M. J. Willatt, M. Stricker, T. Junge, and M. Ceriotti, Efficient Implementation of Atom-Density Representations, J. Chem. Phys. 154, 114109 (2021). [doi:10.1063/5.0044689](https://aip.scitation.org/doi/10.1063/5.0044689)" ] }, { "cell_type": "markdown", "id": "f2381cf3", "metadata": {}, "source": [ "To select features (columns) instead of samples (rows), we use the same selector class as above with a few changes. Let's start with a conservative estimate of half the total number of features:" ] }, { "cell_type": "code", "execution_count": 14, "id": "2fc20224", "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "# select the features with FPS, again the whole training set\n", "n_sparse_feat = 600\n", "feat_compressor = FPSFilter(soap, n_sparse_feat, act_on='feature')\n", "feat_sparse_parameters = feat_compressor.select_and_filter(managers)" ] }, { "cell_type": "markdown", "id": "13c54d51", "metadata": {}, "source": [ "Note that here we use the full feature matrix again, but in real scientific applications this might not be possible or practical due to the size of the full matrix. In this case we might sub-select samples, either randomly or using our existing sparse sample selection*, to obtain a smaller matrix to work with.\n", "\n", "*You might be wondering how we do the sparse sample selection in the first place if the full feature matrix is too big to fit in memory. There are various ways around this, including iterative or hierarchical selection methods, or we could just do feature selection first on a random subset and then do sample selection on the feature-sparsified matrix. The best method to use will in general depend on the application and is out of the scope of this tutorial." ] }, { "cell_type": "code", "execution_count": 15, "id": "55140518", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dict_keys(['coefficient_subselection', 'selected_feature_ids_global', 'selected_feature_ids_global_selection_ordering'])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "feat_sparse_parameters.keys()" ] }, { "cell_type": "markdown", "id": "1d74eb07", "metadata": {}, "source": [ "The result is a dictionary that tells the `SphericalInvariants` class which coefficents to compute. In order to use it, we must update the `hypers` and create a new, sparsified version of the `soap` feature calculator." ] }, { "cell_type": "code", "execution_count": 16, "id": "a19e712a", "metadata": {}, "outputs": [], "source": [ "hypers['coefficient_subselection'] = feat_sparse_parameters['coefficient_subselection']\n", "soap_fsparse = SphericalInvariants(**hypers)" ] }, { "cell_type": "code", "execution_count": 17, "id": "07d21b97", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Execution time: 1.7077140808105469 s\n" ] } ], "source": [ "start = time()\n", "managers_sparsefeat = soap_fsparse.transform(tqdm(frames_train))\n", "print(\"Execution time: {} s\".format(time() - start))" ] }, { "cell_type": "code", "execution_count": 18, "id": "2d2119cc", "metadata": {}, "outputs": [], "source": [ "# Bit of a hack because we can't directly use the old `sample_compressor` on the new features\n", "sample_compressor_sparse = FPSFilter(soap_fsparse, n_sparse_env, act_on='sample per species')\n", "sample_compressor_sparse.selected_sample_ids_by_sp = sample_compressor.selected_sample_ids_by_sp\n", "X_sparse_sparsefeat = sample_compressor_sparse.filter(managers_sparsefeat)" ] }, { "cell_type": "markdown", "id": "6dff3c7e", "metadata": {}, "source": [ "### How do we know how many features to select?" ] }, { "cell_type": "markdown", "id": "7f13fe5d", "metadata": {}, "source": [ "With FPS, one metric we can look at to get a rough estimate of the required number of features is the Hausdorff distance of each selected point (basically, the distance of the selected point to the closest point in the \"blob\" of already-selected points). This tells us how far (in feature space) each selected point is from the selected set, and _very_ roughly, the amount of diversity it adds to the set. So when the decrease in this distance starts to slow down, it indicates that we are getting less of an information return for each point added to the set." ] }, { "cell_type": "code", "execution_count": 19, "id": "2bef5be2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'FPS selection on features')" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5bnA8d+ThCRkBxJISAhhk31HFFdcirjgirZUrVul9l61m/Xq1bZ6tWp7r9ZabatVSq2K+4JUKy4I1gXZQfZFICFAwhYSCCHLc/84JzqmyWSSzOTMTJ7v5zOfzLxzznued8R55n3fc94jqooxxhjTlBivAzDGGBPeLFEYY4zxyxKFMcYYvyxRGGOM8csShTHGGL8sURhjjPHLEoWJaiKyVUTODEG9q0VkYrDrDSYROVFENopIhYhc6HU8JnJZojBt5n4ZV7pfSPWPniJSICLqU7ZVRG7z2e8CEVkuIgdFZI+IvC8iBd61pHEiMlNE7vUtU9WhqvqhRyEF6n+AR1U1RVVfb0tFoUq4JjLEeR2AiRpTVPU93wKfL/0MVa0RkQnA+yKyHNgEPA1cDHwApACTgLp2izj69QZWex0EgIjEqWqN13GY1rEehWk3qvopzhfXMGAU8KWqvq+OclV9RVW3N7aviJwjImtEpFxEdojILT7vnef2TA6IyCciMqKJOmJE5DYR2Swie0XkRRHp6vP+Se7+B0SkUESuFpHpwOXArW6v6E13269+YYtIgog8LCLF7uNhEUlw35soIkUi8jMRKRGRnSJyTVOfkdsTmy0i+0Rkk4hc7/PeXW7MT7ufw2oRGddEPZuBvsCbbtwJIpIuIk+5MewQkXtFJNbdvp+IfOB+LntE5FkRyXDf+zuQ71PXrfXtanBM38/kLhF5WUSeEZGDwNX+Pn8RSXS33et+/otEpEdTn5NpX5YoTLsQx4nAUGAZsBQYJCK/E5HTRCSlmSqeAn6gqqk4ieYDt94xwAzgB0A34HFgdv0XdQM3AxcCpwI9gf3AY249+cDbwB+ALJxEtlxVnwCeBX7rDuFMaaTeO4Dj3X1GAuOBO33ezwbSgVzgOuAxEenSRDtnAUVufFOB+0TkDJ/3zweeBzKA2cCjjVWiqv2A7Tg9vRRVrQL+BtQA/YHROD2477u7CHC/e9zBQC/gLreuKxvU9dsmYm/oAuBlN9Zn8fP5A1fhfEa9cP473gBUBngcE2qqag97tOkBbAUqgAPu43W3vABQt2w/sBa42We/44EXgVLgCDATSGniGNtxkkFag/I/Afc0KFsPnOoT25nu87XAGT7b5QDVOEOwtwOvNXHsmcC9jbS5vt7NwDk+750FbHWfT8T5wovzeb8EOL6R4/QCaoFUn7L7gZnu87uA93zeGwJUNvPfpT7GHkAV0Nnn/WnAvCb2vRBY1lhdPu0q8nO8u4AFDd739/lfC3wCjPD637M9/v1hcxQmWC7UBnMUPjK1kfFpVf0MuAxARI4FXsD5dX57I3VcgvMr/QERWQncps5QVm/gKhG5yWfbeJxfrA31Bl4TEd95kFqcL9FeOF/4rdET2ObzeluD4+9t0P7DOHMyjdWzT1XLG9TlO7y0q0E9iQGO//cGOgE7RaS+LAYoBBCR7sAjwMlAqvve/mbqbE5hIzE09fn/Hee/wfPukNczwB2qWt3GGEwQ2NCTCQuqugh4FWdYqdH3VfUCoDvwOk5PBJwvo1+raobPI0lVZzVSTSFwdoNtE1V1h/tev6bCayb8YpwvwXr5bllLFQNdRSS1QV07WlFXQ4U4PYpMn7anqepQ9/37cdo5QlXTgCtwhqPqNfwMDgFJ9S/cuY6sBts03KfJz19Vq1X1blUdApwAnAd8rw3tNUFkicJ4wp04vt79JYuIDMIZf/+skW3jReRyEUl3f2EexPklCvAX4AYROc6dB0kWkXMbfNnW+zPwaxHp7dabJSIXuO89C5wpIpeJSJyIdBORUe57u3EmhpsyC7jTrS8T+CXOL+IWUdVCnOGX+93J3RE4cxrPtrSuRureCcwFHhSRNHdiuZ+InOpukoo7fCgiucDPG1TR8DPYgNObOVdEOuH09hqbF/LV5OfvzlMNdxPOQZwhqdqmqzLtyRKF8coBnMSwSkQqgH8CrwFNTZReCWx1z6C5AecXL6q6GLgeZ1J3P85pt1c3UcfvcSaA54pIOU5SOs6tZztwDvAzYB+wHGdiGpyJ9CHu2TiNXY9wL7AYWAmswpmov7eR7QIxDWdupxjn8/iVqr7byroa+h7OsNwanM/qZZx5AoC7gTFAGfAPnN6dr/txkuEBEblFVcuA/wCexOnxHMKZhPenyc8fZ8L/ZZwksRaYTyuSrQkNUbUbFxljjGma9SiMMcb4ZYnCGGOMX5YojDHG+GWJwhhjjF9RecFdZmamFhQUeB2GMcZElCVLluxR1YbXw0RnoigoKGDx4sVeh2GMMRFFRLY1Vm5DT8YYY/yyRGGMMcavqEoUIjJFRJ4oKyvzOhRjjIkaUZUoVPVNVZ2enp7udSjGGBM1oipRGGOMCT5LFMYYY/yyRGGMMcYvSxQ+Xl1axLMLGz2N2BhjOqyoShRtPetpzsqdPP95w7s3GmNMxxZViaKtZz0lJ8RRUdXcrYeNMaZjiapE0VYpCXGUH7FEYYwxvixR+EhNjOOQ9SiMMeYbLFH4SEmIo7K6lpraOq9DMcaYsGGJwkdygrOY7qGqWo8jMcaY8GGJwkeqmyjKq6o9jsQYY8JH2CcKEblQRP4iIm+IyKRQHisl0UkUduaTMcZ8zZNEISIzRKRERL5oUD5ZRNaLyCYRuQ1AVV9X1euBq4FvhzKulK+GnixRGGNMPa96FDOByb4FIhILPAacDQwBponIEJ9N7nTfD5n6HoWdImuMMV/zJFGo6gJgX4Pi8cAmVd2iqkeB54ELxPEb4G1VXRrKuOp7FDb0ZIwxXwunOYpcwHf9jCK37CbgTGCqiNzQ1M4iMl1EFovI4tLS0lYF8FWisB6FMcZ8Jc7rAHxII2Wqqo8AjzS3s6o+ISI7gSnx8fFjWxOATWYbY8y/C6ceRRHQy+d1HlDckgravNZTvCUKY4xpKJwSxSJggIj0EZF44DvA7JZU0NbVY2NjhKT4WBt6MsYYH16dHjsL+BQYKCJFInKdqtYANwLvAGuBF1V1dUvqDcY9s1NsBVljjPkGT+YoVHVaE+VvAW+1tl4RmQJM6d+/f2urICUxjnJLFMYY85VwGnpqs2D0KFITbAVZY4zxFVWJoq1zFODevMjmKIwx5itRlSiC0aPIzejMht3lHK2xpcaNMQaiLFEEw+Rh2Rw8UsPnXza8cNwYYzqmqEoUwRh6KshMBmDvoapghWWMMREtqhJFMIaekuJjAbt5kTHG1IuqRBEMSe7V2YeP2oS2McZAlCWKYAw91fcoDh+1HoUxxkCUJYpgDD11io0hPi6GQ9ajMMYYIMoSRbAkx8dy2OYojDEGsETRqKT4OOtRGGOMK6oSRTDmKMCZp6i0OQpjjAGiLFEEY44CICkhjkOWKIwxBoiyRBEszhyFDT0ZYwxYomhUSkIci7ftZ0tphdehGGOM5yxRNOKGif0AeH3ZDo8jMcYY71miaMSY/C6M6pXBJ5v3eh2KMcZ4LqoSRbDOegIYlpvGJht6MsaY6EoUwTrrCSCvSxIHDldTfqQ6CJEZY0zkiqpEEUy9uiQBULS/0uNIjDHGW5YompDXpTMAq4raPoxljDGRzBJFE4b0TGNQdip/mLfR61CMMcZTliia0Ck2hkvG5FG4r5LScrvbnTGm47JE4cfIXhkArCw64HEkxhjjnbBPFCLSV0SeEpGX2/vYg3NSAdiw206TNcZ0XJ4kChGZISIlIvJFg/LJIrJeRDaJyG0AqrpFVa/zIs7UxE5kpSaw2a6nMMZ0YF71KGYCk30LRCQWeAw4GxgCTBORIe0f2jf1y0q2NZ+MMR2aJ4lCVRcA+xoUjwc2uT2Io8DzwAWB1iki00VksYgsLi0tDVqsI/MyWFlUxvJCm6cwxnRM4TRHkQsU+rwuAnJFpJuI/BkYLSK3N7Wzqj6hquNUdVxWVlbQgrrmxD7ExQr/8cwSVDVo9RpjTKQIp0QhjZSpqu5V1RtUtZ+q3u+3giCu9VQvOz2Rn581iOKyI+w+aKfJGmM6nnBKFEVAL5/XeUCxR7F8w6heztpR76ze5XEkxhjT/sIpUSwCBohIHxGJB74DzG5JBcFcFNDX0J7p5KQn8ti8TUGt1xhjIoFXp8fOAj4FBopIkYhcp6o1wI3AO8Ba4EVVXd3CeoM+9ASQ2CmWaePzKSmvoqrG7qVtjOlY4rw4qKpOa6L8LeCtNtT7JvDmuHHjrm9tHU3pnpoAQGl5FXnuyrLGGNMRhNPQU1jrnuYkihJb98kY08FEVaII1dATQPfURAC27T0U9LqNMSacBZwoRCQ5lIEEQ6gmswFy0p1E8ZMXVvCTF5YHvX5jjAlXzSYKETlBRNbgTDAjIiNF5I8hj6wVQtmj6JaSwKzrj2fy0GxeW7aDXWVHgn4MY4wJR4H0KH4HnAXsBVDVFcApoQyqtULZowCY0K8bN57eH4BfvvEFhfsOh+Q4xhgTTgIaelLVwgZFHfYc0SE5aXz3uHw+2riHa2cuoqa2zuuQjDEmpAJJFIUicgKgIhIvIrfgDkN1RDExwn0XDeehy0aysaSC2SvC4uJxY4wJmUASxQ3Af+Is2lcEjHJfh51QzlE0dNbQbIbnpnPPnDVs3WNnQhljoleziUJV96jq5araQ1W7q+oVqrq3PYJrqVDPUfiKiREemTaayupaHl+wOeTHM8YYrwRy1tPfRCTD53UXEZkR2rAiQ5/MZE4ekMW7a0o4VFXjdTjGGBMSgQw9jVDVr+7ao6r7gdGhCymyXDo2jz0VVfzf3PVeh2KMMSERSKKIEZEu9S9EpCserRHVnPaco6g3aWg2Zw/LZvbyYo5Ud9iTwYwxUSyQRPEg8ImI3CMi9wCfAL8NbVit055zFL6+N6GAvYeO8vSnW9v1uMYY0x4Cmcx+GpgK7AZKgItV9e+hDiySTOjXjVOPyeLRDzZRUm5XbBtjokugaz2tA14F3gAqRCQ/dCFFpl+cN4QjNXU8NHeD16EYY0xQBXLW0004vYl3gTnAP9y/xkf/7ilcPDqXN5YXc/BItdfhGGNM0ATSo/gRMFBVh6rqCFUdrqojQh1YJLr8uN5UVtfyypIir0MxxpigCWgJD6D9TiNqAy/OevI1PC+dYwu68Ni8zVTYdRXGmCgRSKLYAnwoIreLyE/rH6EOrDW8OuvJ1x3nDmFPRRWPz7ertY0x0SGQ6yG2u49492H8GNUrg3OGZ/PEgi0s2bafi8fkccmYXETE69CMMaZVmk0Uqnp3ewQSTW6ZNBARYf2ucm55aQX7DlUx/ZR+XodljDGt0myiEJEs4FZgKJBYX66qp4cwrojWNyuFx747BlXlhmeWcN9b6yg+cIS7zh/qdWjGGNNigcxRPItzHUUf4G5gK7AohDFFDRHhd98exVlDezDzk60s2bbf65CMMabFAkkU3VT1KaBaVeer6rXA8SGOK2okxcfx4GWjyExJYOqfP+H0Bz+0hGGMiSiBJIr6q8d2isi5IjIayAthTN8gIsnuUud/EZHL2+u4wZSSEMdLN0zgxtP6c+RoLd//2yK737YxJmIEkijuFZF04GfALcCTwI/bclARmSEiJSLyRYPyySKyXkQ2ichtbvHFwMuqej1wfluO66U+mcn8bNJAnr3+eGrqlO888Rl7Kqq8DssYY5oVSKLYr6plqvqFqp6mqmOBfW087kxgsm+BiMQCjwFnA0OAaSIyBKf3UuhuFvHrePfJTOZ3l41ix4FKFmwo9TocY4xpViCJ4g8BlgVMVRfw78lmPLBJVbeo6lHgeeACnPt01w91NRmviEwXkcUisri0NLy/gCcOzCIhLoaVRRFxwbsxpoNr8vRYEZkAnABkNbgSOw2IDUEsuXzdcwAnQRwHPAI8KiLnAm82tbOqPiEiO4Ep8fHxY0MQX9DExcYwpGcaMz/ZypSROYzt3dXrkIwxpkn+ehTxQApOMkn1eRzEuT9FsDV26bKq6iFVvUZVf6iqz/qrIByW8AjUfRcNp0tSJ+5+cw11dep1OMYY06QmexSqOh+YLyIzVXUbgIjEACmqejAEsRQBvXxe5wHFLalARKYAU/r37x/MuEJicE4av5wyhJ+8sIKL/vQJ15/ch/NG9PQ6LGOM+TeBzFHcLyJpIpIMrAHWi8jPQxDLImCAiPQRkXjgO8DsllQQST0KgAtH5fLAxcPZeaCSm2YtY+l2u77CGBN+AkkUQ9wexIXAW0A+cGVbDiois4BPgYEiUiQi16lqDXAj8A6wFnhRVVe3sF5PlxlvKRHhO+PzmXfLRLolJ/Dj55cz6/Pt1NpQlDEmjIiq/y8lEVkNjAKeAx5V1fkiskJVR7ZHgK0xbtw4Xbx4sddhtMiCDaU88PY61uw8yPF9u3LT6QM4oV83W3XWGNNuRGSJqo5rWB5Ij+JxnPWdkoEFItIbZ0I77ERaj8LXKcdk8Y+bT+J/LhjKppIKLn9yIb+avZrmErkxxoRasz2KRncSiXOHisJSJPYofB0+WsMDb6/j6U+3MW18Ly4d14tReRnExFjvwhgTOk31KPxdR3GFqj7j5252DwUtOvMNSfFx3DVlKBVVNcz6vJBZnxdy8ZhcHrx0pA1FGWPanb/7USS7f1PbI5BgiKTTY5sTEyM8dNkobpk0kMfmbeLZhdtZsKGUzJQELh3Xi+tO6uN1iMaYDqJVQ0/hLtKHnhpSVf7+2TbW7ixnzc6DrCg8wDPXHcdJAzK9Ds0YE0VaM/T0iL8KVfXmYARmmicifG9CAQBHqmuZ/PACrn96Mb/79kgmD8vxNjhjTNTzd9bTEveRCIwBNrqPUYTpKq6RfNZToBI7xfL89An0757Cz19ayb827qGmts7rsIwxUSyQ6yjmAZNUtdp93QmYq6qntUN8rRJtQ0+N2XGgkiufWsiW0kNkpiTwyg8n0LtbcvM7GmNME9pyHUVPvjmhneKWGQ/lZnTmtf84kfsvHk5VTS13vv4FB49UN7+jMca0kL+znuo9ACxzexYApwJ3hSwiE7D0zp2YNj6f8iPV3PfWOkbePZfRvTL4y/fG0S0lwevwjDFRIqCznkQkG+feEAALVXVXSKNqJZ/TY6/fuHGj1+G0G1Xlsy37WPjlXh79YBNJ8bH8+YqxnNDfzooyxgSuqaEnOz02ysxbX8K9c9ZQUl7FnJtOsnkLY0zA2jJHYSLIaQO7M/Oa8cSIcNWMz/lwfYmdFWWMaRN/96C2S38jVK+uScy4ehxHquu4+q+LOO3BD1lZdMDrsIwxEcpfj+JlABF5v51iMUE0tndX5t86kT9MG03FkRoueOxjnvrXlxypDstLYIwxYazJOQoRWQa8Dnwf+F3D91U1bBcF7MhzFI0p3HeYn764nEVb99MlqROXjMnjexMKyO+W5HVoxpgw0po5iu8AR3BOoU1t5BF2OsKV2a3Rq2sSL0yfwDPXHceEft2Y+clWJj08nwfnrrcehjGmWf56FD9S1d+LyC9V9X/aOa42sR6Ff8UHKvnlG6t5b+1uYmOE/5zYj2+Pz6dneqItY25MB9bi02NFZLmqjhKRpao6JuQRBpEliuapKvPWl/D6smJmrygGYHhuOk9dNY7uaYkeR2eM8UJrEsUsYAKQBWz2fQtQVR0RikCDwRJF4FSVJdv2s7zwAA+9uwGAEXnpXDq2F5eMzfM4OmNMe2rxMuOqOs29Ivsd4PxQBme8IyKMK+jKuIKuHN+3Gy8uLuTTzXv52Usr+HjTHm48vT99s1K8DtMY46Hm1noqBVap6rb2CMZ4a1huOsNy06mpreP+t9fx3MLtvLGimFOPyeJ/p46w9aOM6aD8XpmtqrVApojEt1M8JgzExcbwi/OGsODW07j8uHzmrS/h2F+/x+VPfkZpeZXX4Rlj2lkg96N4HOfGRbOBQ/Xl7XUdhYj0Be4A0lV1aiD72BxFcK0pPshbq3byl4+2MCg7laevPY70pE5eh2WMCbK2rPVUDMxxt23RdRQiMkNESkTkiwblk0VkvYhsEpHb/NWhqltU9bpAjmdCY0jPNG45ayC3TBrIiqIyTnvwQ/7zuaWsKLRlQYzpCJq9H4Wq3g0gIqnOS61oQf0zgUeBp+sLRCQWeAz4FlAELBKR2UAscH+D/a9V1ZIWHM+E0PWn9KV/jxRmfryVTzfvZe7qXZx6TBZTx+YxaUg2MTF2DYYx0ajZRCEiw4C/A13d13uA76nq6ub2VdUFIlLQoHg8sElVt7j1PQ9coKr3A+e1KPpvxjkdmA6Qn5/f2mpMM04b2J3TBnanpPwIj36wibmrd/Pe2qVkpyUyeVg2ZwzuTtfkePIykmx4ypgoEcgcxSfAHao6z309EbhPVU8I6ABOopijqsPc11OByar6fff1lcBxqnpjE/t3A36N0wN50k0oftkcRfs5WlPHW6t28taqnczfUEpVjbOkuQicMzyHhy4bSUJcrMdRGmMC0eLrKHwk1ycJAFX9UETacjecxsYnmsxWqroXuCGgir++w10rQzMtFR8Xw4Wjc7lwdC6HqmpYun2/+/cATyzYQtH+Ss4dns3w3AyG5aaRmmi9DGMiTSCJYouI/AJn+AngCuDLNhyzCOjl8zoPZ8K8zVT1TeDNcePGXR+M+kzLJCfEcfKALAAmD8thaM80Hpy7gfveWgdAp1jh8SvHcvqgHl6GaYxpoUCGnroAdwMn4fQGFgB3qer+gA7w70NPccAG4AxgB7AI+G4gcx4BHKtD3jM73O2tqGLljjIeeGsdRfsP8+Mzj2HysGyy0xPpFGs3WTQmXHhyz2x3vaiJQCawG/iVqj4lIucAD+Oc6TRDVX8dzOPaHEV42llWya0vr+SjjXsA6JLUiVvOGsikIdlkpdpV38Z4rTWLAr6J/7mDsFv/yXoUkWFN8UGWFx7ghUXbWVFURqdYYdr4fK49sQ9ZqQkkJwQyImqMCbbWJIpT3acXA9nAM+7racBWVf3vUAQaDNajiAzVtXWsLDrAK0t38OKiQmrqnH+LY/IzuOfCYQztme5xhMZ0LK0eehKRBap6SnNl4cB6FJFr+97DfLSplD3lR3lm4TYOV9Vw9YkFnDawO8Pz0u0UW2PaQVsSxVrgXJ8L5PoAb6nq4JBEGgTWo4hsG3eXc99ba5m3vhSA1MQ4Hrx0JJOGZnscmTHRrS2JYjLwBLDFLSoAfqCq7wQ7yGCxRBEd9h06yqKt+3hw7np27K/k7guGcd6IHBI7We/CmFBo01lPIpIADHJfrlPVsFxr2oaeotPanQe5adYyNpVUMCg7lR+dMYAzh/SwU2uNCbK29CguBf6pquUicifOkuP3qurS0ITadtajiD41tXW8tKSIB95eR1llNf2ykpl+Sl9OOSaLrJQE4ixpGNNmbUkUK1V1hIichLO66/8B/62qx4Um1LazRBG9amrr+Meqnfz2n+vZcaAScOYwJg7szs8nDSS/W5LHERoTudqSKJap6mgRuR/ntqjP1ZeFKtjWsqGnjqOmto6PN++laP9hVhWV8eaKYiqrazl5QBb5XZPISOrElJE9OaZHQLdOMcbQtkQxB2epjTOBsUAl8LmqjgxFoMFgPYqOZ8eBSp79bBv/XL2L/YeOUlZZTafYGK44vjfnjshhTH4Xr0M0Juy1JVEkAZNxehMbRSQHGK6qc0MTattZojAl5Ue4e/Ya3lm9i5o6JTstkZG90rnh1H6MtqRhTKPakigavQuQqm4PUmxBZ4nC1Dtw+CjPLtzO5pIK/rl6F4eP1pKTnsi5w3OYfkpfuqcleh2iMWGjLYliFc6aTwIkAn2A9ao6NBSBtoXNURh/yo9U8/fPtrGysIx31+4mNkY4a2g208b3Ykx+F7s+w3R4QVs9VkTG4Fxw94NgBRds1qMwzdm29xB/+WgLsz4vpLZOyUyJ54R+mdx8xgD6d0/xOjxjPBHUZcZFZKmqjglKZCFgicIEavvew6zaUcbbX+xk/vpSMpI7ccOp/ThraDaZKbb0uelY2jL09FOflzE4F9x1U9Wzghti8FiiMK0xd/UufvT8ciqra0mOj+X+S0bQNzOZvlnJJMXb0ucm+rUlUfzK52UNsBV4RVWPBDXCILJEYVqrpraOjSUV3PDMErbtPQw49wU/oV83Lhqdy5QRPYmJaey278ZEPk/ucOcVSxSmrY5U17K88AAHDh9l0db9vLtmN9v3HeaMQd256oQC+mYlk5vRGRFLGiZ6tKVHkQXcCgzFOesJAFU9PdhBtpWd9WRCpa5OeXzBFn77zjrq/5fJzejMif27cWL/TE7ol2m3czURry2JYi7wAnALcANwFVCqqv8VikCDwXoUJlRKy6vYUlrB+t3lfLp5L59s3ktZZTUAvbslMf2Uvpw5uAc97PoME4HakiiWqOrY+sUB3bL5qnqq3x09ZInCtJfaOmV1cRmfbt7LS0uK2FRSQXxcDMf37caY/AzG9u7CCf0yibV5DRMBmkoUgZzKUe3+3Ski5wLFQF4wgzMmUsXGCCPyMhiRl8G1J/VhdfFBXllSxOdf7uP3729EFQZlp3JCv0xOHZjF4JxUuqdab8NElkB6FOcBHwG9gD8AacDdqjo79OG1jvUoTDgoP1LNswu389rSHWzZU0F1rfP/WkG3JM4f2ZNvj88nN6Ozx1Ea8zU768kYD5VVVrOm+CCri8t4dekO1u46iCpcfUIBlx+XT+9uycTH2c2XjLdanChE5A84azw1SlVvDl54wWWJwoS7rXsOcc+cNby/rgSA5PhYJvTrRm5GZ/K6JNGra2eOLehKN7s63LSj1sxR+H7T3g38qqkNQ0lELgTOBboDj4Xz8ubGBKogM5knrxrHqh1lbCk9xIfrS1i7s5yFW/ZRXlUDOMnjP07rz7Un9qFzvC1YaLwT0NBTa+9oJyIzgPOAElUd5lM+Gfg9EAs8qaoPBFBXF+D/VPW65ra1HoWJVKrKwcoaNu+p4E8fbubdNc4qt3ldOnNCv0yuPbGAAXbXPhMibZqjaO0igCJyCgd7hVgAABOpSURBVFABPF2fKEQkFtgAfAsoAhYB03CSxv0NqrhWVUvc/R4EnlXVpc0d1xKFiRaLtu5j/vpS1u0q5/11uwEY2COVYwu6MmloDzv11gRVW06PbTVVXSAiBQ2KxwObVHWLG9jzwAWqej9O7+MbxFkj4QHgbX9JQkSmA9MB8vMbvdeSMRHn2IKuHFvQFXBu9/r0p1tZv6ucWZ9v5++fbSM5PpaB2akMykljcE4aQ3LSGJOfYUuLmKBqMlGISDlfT2YnicjB+rcAVdW0Vh4zFyj0eV0EHOdn+5tw7tedLiL9VfXPjW2kqk8AT4DTo2hlbMaErdyMztx+9mAADlXVsGBDKQu/3MfanQeZs6KY5xY6N50clpvGpCHZpCTE0SMtkdMHdbc5DtMmTSYKVQ3VQGhjP3X8nV31CPBIQBV/vdZTK0MzJjIkJ8Rx9vAczh6eAzhzGzvLjvDKkiJeXFLIQ+9u+GrblIQ4fjixH9ecWGDLpZtW8eJfTRHOxXv18nCu9jbGtJKI0DOjMzedMYCbzhjAkepaqqrrWLPzIH/8cBP/+8565q7Zzc++dQwnD8i0oSnTIiG/4M6do5jjM5kdhzOZfQawA2cy+7uqujpYx7TJbGO+pqq8vnwHv3h9NRVVNZzQrxunDezOmN5dyOvSme6pCZY4DODRldkiMguYCGQCu4FfqepTInIO8DDOmU4zVPXXQTqeLTNuTBMqj9by7MJt/PXjrew4UPlV+fkje/Lwt0fZDZmMLeFhjPlaSfkRlm8/wKdb9vLXj7dy6jFZnDU0m1MHZtn6Ux2YJ6fHGmPCU/fURCYNzeZbQ3qQmtiJlxYXMn9DKTECN54+gMuPy7d7apivRFWPwoaejGkdVWVTSQX3vbWWeetLSewUw9Un9OGWSccQF2uLFXYUNvRkjAnIllInYby3toTe3ZKYPDSbgdmpnDsih4Q4ux4jmnWIRGE9CmOCo65OeXNlMTM+3sra4oMcra2ja3I8l47L46oJBfS0eYyo1CESRT3rURgTPDW1dby3djevLN3B+2t3IyKcNrA7k4b0YHheOoOyU+302ihhk9nGmFaJi41h8rAcJg/LoXDfYZ75bBuvLtvBe2udRQpz0hO58fT+nD6oOznp1tOIRlHVo7ChJ2PaR01tHdv3HWbBhlJeWlLE6uKDxMUIw/PSGZmXwR3nDqaTTYJHHBt6MsaERHVtHat2lPHykiI2l1Sw8Mt99M1M5s9XjuUYu3dGRLGhJ2NMSHSKjWFMfhfG5HcBYPaKYu6Zs4ZJv1tA38xkenVNok9mMt2S47l4bJ5d0BeBLFEYY4Lq/JE9GZGbzvOLCvlyTwVF+ytZsm0/FVU1PPHRFv5r8iCO79uNPpnJdtOlCGFDT8aYdrF972F++uJyFm/bD0D31ASO7dOVc4fn0C8rhX5ZyXZxn8c6xByFTWYbE97q6pQ1Ow+yZudB5m8o5V8b91BWWQ04N2a689zBX91jw7S/DpEo6lmPwpjIUHm0lk0lFWzYXc5fP/mSL3Yc5KT+mZw0IJNLx+bRLSXB6xA7FEsUxpiwdqiqhsfmbeLtL3bx5Z5DpCTEkd81iZ4ZifTLSmFgdioTB3YnvXMnm9sIEUsUxpiIsWF3OTP+9SV7Kqr4cs8hvtxziDr3qyorNYEHLh7O6YO62xXhQWaJwhgTserqlLlrdlN8oJKnP93K1r2HGZmXzreG9ODMIT0YlJ3mdYhRwRKFMSYq7K2o4qUlRby5opjVxQcBuOL4fG6dPIi0xE4eRxfZOkSisLOejOlYSsur+OOHm/jrx1uJixF+8q1j+OGp/ey2rq3UIRJFPetRGNOxfLGjjN/8cx0fbdxD/+4pnDm4B2N7d2FYbpotVNgCliiMMVGttk6ZvWIHz39eyJJt+6mpU2JjhCuP7801JxbQu1uy1yGGPUsUxpgOY29FFdv2HebZz7Yze8UOqmuVrNQERvfK4H+njiQ9yeYyGmOJwhjTIe0+eIQ3lu9gw+4KXl1aRN+sFC4Y2ZPTB3cnr0sSaYlxdpqtyxKFMabDe2HRdp5fVMiy7Qe+KjtzcA/uuXAo2WmJHT5hRGyiEJHBwI+ATOB9Vf1Tc/tYojDG+LOppIL1u8qZv6GEFxcXAZCdlsikoT04Z3gOx/ft5nGE3vAkUYjIDOA8oERVh/mUTwZ+D8QCT6rqAwHUFQP8RVWva25bSxTGmECtKDzAsu37+WzLPj5YV8LR2jp+O3UEl43r5XVo7c6rRHEKUAE8XZ8oRCQW2AB8CygCFgHTcJLG/Q2quFZVS0TkfOA24FFVfa6541qiMMa0RkVVDec/+i9KDlZx8oBMeqQlMmVkT0b3yugQ12Z4NvQkIgXAHJ9EMQG4S1XPcl/fDqCqDZNEY3X9Q1XPbW47SxTGmNbaXFrBvXPWULi/ksJ9h6mqqSMjqROjemVw3oieXDQ6N2oXJQynW6HmAoU+r4uA45raWEQmAhcDCcBbfrabDkwHyM/PD0acxpgOqF9WCn+9ZjwAZZXVfLBuNx9v2sunm/dyy0sr2FJawa2TB3kcZfvyIlE0loqb7Nao6ofAh81VqqpPiMhOYEp8fPzYVkdnjDGu9M6duGh0HheNzqOmto5rZi7iT/M3s2z7Aa4/pQ+nD+rhdYjtwov7DhYBvrNEeUBxMCpW1TdVdXp6enowqjPGmK/Excbwu2+P4rvj8ykuq+TamYsZe8+73P3maj7aWEq4n0HaFl7MUcThTGafAezAmcz+rqquDsKxbFFAY0zIHa2p4+UlRczfUMK8daUcra1jcE4a/33OICb07Rax9/726qynWcBEnGsgdgO/UtWnROQc4GGcM51mqOqvg3lcm8w2xrSXI9W1vLS4kCc+2kLhvkpSE+O4ZEwe/33OYOLjIithROwFdy1hPQpjjFcOH63h3TW7mbeuhNeXF3P+yJ7cdHp/BvRI9Tq0gHWIRFHPehTGGC89/N4GHnl/I3XqLBFy38XDyExOCPtrMTpEorAehTEmXJQcPMKLiwt56N0N1CnkpCfSv3sKF47K5ZKxeV6H16gOkSjqWY/CGBMulm7fz5Kt+1m8bR8bd1fw5d5DXDQql/86exA90hK9Du8bwumCO2OM6TDG5HdhTH4Xrqcvh6pq+OUbq3l1WRHzN5Ty87MGMmloNl2T470O06+o6lHY0JMxJhKsKirjztdXsaKojKzUBO6aMpQzBncnsVOsp3HZ0JMxxoSRujpl3voSfvH6FxSXHSE2RuiZkUhOemd+eGo/ThvUvd1jskRhjDFhqLZO+WTzHj7/ch/b9h5m6fb9lBys4jdTh3PByNx2PVOqQyQKG3oyxkS6vRVVXP7kQtbtKqdTrDA8N53pp/Rl8rCckB+7QySKetajMMZEsto6ZfaKHazdWc5zC7dTUVXDt4b0YEx+F8YVdGFEXjoJccGfz7BEYYwxEehIdS0/e2kFKwoPULS/EoCCbkncdPoALh6TG9T7fFuiMMaYCKaq7Dp4hM+27OWP8zazsaSCHmkJPHTZKE7snxmUY3SIRGFzFMaYjkBVmb2imD98sImdByq57uS+nD0sm0HZqW3qYXSIRFHPehTGmI5gV9kRfvD3xawoKgOgb1Yyj18xttULEdqV2cYYE2Wy0xN548aT2FNRxT+/2MW7a3aT1yUp6MexRGGMMREuMyWBK47vzRXH9w5J/ZF1Vw1jjDHtzhKFMcYYv6IqUYjIFBF5oqyszOtQjDEmakRVolDVN1V1enp6utehGGNM1IiqRGGMMSb4LFEYY4zxyxKFMcYYvyxRGGOM8Ssql/AQkVJgWyt3zwT2BDEcr0RLO8DaEq6sLeGpLW3prapZDQujMlG0hYgsbmytk0gTLe0Aa0u4sraEp1C0xYaejDHG+GWJwhhjjF+WKP7dE14HECTR0g6wtoQra0t4CnpbbI7CGGOMX9ajMMYY45clCmOMMX5ZonCJyGQRWS8im0TkNq/jaY6IzBCREhH5wqesq4i8KyIb3b9dfN673W3behE5y5uoGycivURknoisFZHVIvIjtzyi2iMiiSLyuYiscNtxt1seUe3wJSKxIrJMROa4ryOyLSKyVURWichyEVnslkVqWzJE5GURWef+PzMh5G1R1Q7/AGKBzUBfIB5YAQzxOq5mYj4FGAN84VP2W+A29/ltwG/c50PcNiUAfdy2xnrdBp+4c4Ax7vNUYIMbc0S1BxAgxX3eCVgIHB9p7WjQpp8CzwFzIvzf2FYgs0FZpLblb8D33efxQEao22I9Csd4YJOqblHVo8DzwAUex+SXqi4A9jUovgDnHxHu3wt9yp9X1SpV/RLYhNPmsKCqO1V1qfu8HFgL5BJh7VFHhfuyk/tQIqwd9UQkDzgXeNKnOCLb0oSIa4uIpOH8SHwKQFWPquoBQtwWSxSOXKDQ53WRWxZpeqjqTnC+fIHubnnEtE9ECoDROL/GI6497lDNcqAEeFdVI7IdroeBW4E6n7JIbYsCc0VkiYhMd8sisS19gVLgr+6Q4JMikkyI22KJwiGNlEXTecMR0T4RSQFeAX6sqgf9bdpIWVi0R1VrVXUUkAeMF5FhfjYP23aIyHlAiaouCXSXRsrCoi2uE1V1DHA28J8icoqfbcO5LXE4Q85/UtXRwCGcoaamBKUtligcRUAvn9d5QLFHsbTFbhHJAXD/lrjlYd8+EemEkySeVdVX3eKIbY87HPAhMJnIbMeJwPkishVnKPZ0EXmGyGwLqlrs/i0BXsMZfonEthQBRW5PFeBlnMQR0rZYonAsAgaISB8RiQe+A8z2OKbWmA1c5T6/CnjDp/w7IpIgIn2AAcDnHsTXKBERnDHXtar6kM9bEdUeEckSkQz3eWfgTGAdEdYOAFW9XVXzVLUA5/+HD1T1CiKwLSKSLCKp9c+BScAXRGBbVHUXUCgiA92iM4A1hLotXs/gh8sDOAfnbJvNwB1exxNAvLOAnUA1zq+G64BuwPvARvdvV5/t73Dbth442+v4G7TlJJzu8Epgufs4J9LaA4wAlrnt+AL4pVseUe1opF0T+fqsp4hrC864/gr3sbr+/+9IbIsb2yhgsfvv7HWgS6jbYkt4GGOM8cuGnowxxvhlicIYY4xfliiMMcb4ZYnCGGOMX5YojDHG+GWJwoQVEbnDXXl1pbvS53HNbD9TRKa24jgFIvLdVuzX6PFEZJAb7zIR6deKen8sIkkt3S/AunuKyMst3OdqEXk0FPGYyGOJwoQNEZkAnIezkuwInAvWCv3v1WoFQIsThR8XAm+o6mhV3dyK/X8MtChRiEhcINuparGqtjiZGlPPEoUJJznAHlWtAlDVPeouvSAiY0Vkvruo2zv1yxX4amobEekvIu+Jc5+Ipe4v/geAk91ewE/cxfz+V0QWub2ZH7j7iog8KiJrROQffL3Ymu9xz8H5ov++iMxzy64Q594Uy0XkcRGJdcv/JCKL5Zv3q7gZ6AnM89m/wqf+qSIy030+U0Qecrf7jXvV8Qw37mUi8m+rHru9py/c51eLyKsi8k9x7l3wW5/trhGRDSIyH2cJj/ryLBF5xT3GIhE50S1/RER+6T4/S0QWiIh9p0Qjr68ytIc96h9ACs5V2RuAPwKnuuWdgE+ALPf1t4EZ7vOZwNRmtlkIXOQ+T8T55T4R92pjt3w6cKf7PAHnytc+wMXAuzj3LOkJHACmNhL7XcAt7vPBwJtAJ/f1H4Hvuc+7un9jcdaCGuG+3orP/RKACp/nU4GZPu2dg3tPAeA+4Ar3eYb72SU3iK0A974lwNXAFiDd/Sy24awFlANsB7Jw7nHwMfCou89zwEnu83ycpVZwP8fVwGk4V/328/rfkD1C8wio62pMe1DVChEZC5yM8+Xzgjh3G1wMDAPedZaFIhZn+RJfAxvbxl3jJ1dVX3OPcQTA3cbXJGCEz/xDOs66OKcAs1S1FigWkQ8CaMoZwFhgkXuczny9SNtl4ixzHYfz5TwEZymGlnjJjac+7vNF5Bb3dSLul7mf/d9X1TIAEVkD9AYygQ9VtdQtfwE4xt3+TGCIz2eWJiKpqlouItcDC4CfaOuG3EwEsERhwor7Bfgh8KGIrMJZ4GwJsFpVJ/jZVRrbRpwbvQRCgJtU9Z0G+59Dy5dlFuBvqnp7g7r6ALcAx6rqfnc4KbGJOnyP2XCbQw2OdYmqrm9BfFU+z2v5+nugqXbGABNUtbKR94YDe3F6WyZK2XiiCRsiMlBEBvgUjcIZGlkPZLmT3YhIJxEZ2mD3RrdR574WRSJyoVue4J5dVI5z29V67wA/FGe5c0TkGHFWGl2As/pmrDvncVoATXkfmCoi3d26uopIbyAN50u+TER64NwboV7DeHaLyGB3zP8iP8d6B7hJ3J/7IjI6gPgasxCYKCLd3M/gUp/35gI31r8QkVHu397Az3BuNHW2NHOGmolclihMOEkB/uZOHK/EGZa5S53b007FmbxdgTOPcYLvjs1scyVws1vnJ0A2znBPjTvB/ROc232uAZa6E7+P4/zSfg1nRc5VwJ+A+c01QlXXAHfi3FFtJc4cR46qrsBZXXY1MANnHqDeE8Db9ZPZODejmQN8wL8Ps/m6B2d+ZqUb9z3NxddEzDtx5lk+Bd4Dlvq8fTMwzp3kXwPc4Camp3DmZYpxVi9+UkSa6iGZCGarxxpjjPHLehTGGGP8skRhjDHGL0sUxhhj/LJEYYwxxi9LFMYYY/yyRGGMMcYvSxTGGGP8+n9Byl4JEDfWAgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fps_dists = feat_compressor.get_fps_distances()\n", "plt.semilogy(fps_dists)\n", "plt.xlabel('Selected feature index')\n", "plt.ylabel('Hausdorff distance')\n", "plt.title('FPS selection on features')" ] }, { "cell_type": "markdown", "id": "da9b6430", "metadata": {}, "source": [ "From this plot, it seems that we might even be able to get away with 200 or 300 features, but we should test this more rigorously. What we really care about is how well these sparsified features reproduce the _kernel_ -- so let's look at the difference in the kernel between the sparse _environments_, with and without _feature_ sparsification." ] }, { "cell_type": "code", "execution_count": 20, "id": "caef3793", "metadata": {}, "outputs": [], "source": [ "K_MM_fullfeat = kernel(X_sparse)\n", "kernel_sparsefeat = Kernel(soap_fsparse, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')\n", "K_MM_sparsefeat = kernel_sparsefeat(X_sparse_sparsefeat)" ] }, { "cell_type": "code", "execution_count": 21, "id": "00de1098", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7171404793446092" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(K_MM_fullfeat - K_MM_sparsefeat, 'fro')" ] }, { "cell_type": "code", "execution_count": 22, "id": "5f5b489a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "104.82383972597799" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(K_MM_fullfeat, 'fro')" ] }, { "cell_type": "markdown", "id": "5ca20ac8", "metadata": {}, "source": [ "So the difference in the norms is less than 1% of the norm of the kernel matrix. Let's see if we can use an even smaller number of features:" ] }, { "cell_type": "code", "execution_count": 23, "id": "aa01d850", "metadata": {}, "outputs": [], "source": [ "def get_feature_sparsified_soap(Nselect, full_soap, managers, sample_compressor):\n", " feat_compressor = FPSFilter(full_soap, Nselect, act_on='feature')\n", " feat_sparse_parameters = feat_compressor.select_and_filter(managers)\n", " soap_hypers = full_soap._get_init_params()\n", " soap_hypers['coefficient_subselection'] = feat_sparse_parameters['coefficient_subselection']\n", " soap_fsparse = SphericalInvariants(**soap_hypers)\n", " managers_sparsefeat = soap_fsparse.transform(managers)\n", " sample_compressor_sparse = FPSFilter(soap_fsparse, n_sparse_env, act_on='sample per species')\n", " sample_compressor_sparse.selected_sample_ids_by_sp = sample_compressor.selected_sample_ids_by_sp\n", " X_sparse_sparsefeat = sample_compressor_sparse.filter(managers_sparsefeat)\n", " return soap_fsparse, managers_sparsefeat, X_sparse_sparsefeat" ] }, { "cell_type": "code", "execution_count": 24, "id": "ab26edf1", "metadata": {}, "outputs": [], "source": [ "soap_smallsparse, managers_sparsefeat, X_sparse_sparsefeat = get_feature_sparsified_soap(\n", " 200, soap, managers, sample_compressor\n", ")" ] }, { "cell_type": "code", "execution_count": 25, "id": "b859802b", "metadata": {}, "outputs": [], "source": [ "kernel_sparsefeat = Kernel(soap_smallsparse, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')\n", "K_MM_sparsefeat = kernel_sparsefeat(X_sparse_sparsefeat)" ] }, { "cell_type": "code", "execution_count": 26, "id": "d2bf5422", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.8514895312542564" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(K_MM_fullfeat - K_MM_sparsefeat, 'fro')" ] }, { "cell_type": "markdown", "id": "4397853b", "metadata": {}, "source": [ "Barely a difference. So let's proceed with 200 features and see how accurate our model is:" ] }, { "cell_type": "markdown", "id": "5ca42c50", "metadata": {}, "source": [ "## Model assessment: Sparse samples and sparse features" ] }, { "cell_type": "code", "execution_count": 27, "id": "0be4e95b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Execution: 1.8869318962097168 s\n" ] } ], "source": [ "zeta = 2\n", "\n", "start = time()\n", "sparse_hypers = soap_smallsparse._get_init_params()\n", "sparse_hypers['compute_gradients'] = True\n", "soap_fsparse_grads = SphericalInvariants(**sparse_hypers)\n", "kernel_fsparse = Kernel(soap_fsparse_grads, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')\n", "\n", "KNM = compute_KNM(tqdm(frames_train, leave=True, desc=\"Computing kernel matrix\"), X_sparse_sparsefeat, kernel_fsparse, soap_fsparse_grads)\n", "\n", "model_fsparse = train_gap_model(\n", " kernel_fsparse,\n", " frames_train,\n", " KNM,\n", " X_sparse_sparsefeat,\n", " e_train,\n", " energy_baseline,\n", " grad_train=-f_train,\n", " lambdas=[1e-12, 1e-12],\n", " jitter=1e-13\n", ")\n", "\n", "# save the model to a file in json format for future use\n", "dump_obj('zundel_model_featsparse.json', model)\n", "np.savetxt('Structure_indices.txt', ids)\n", "print (\"Execution: \", time()-start, \"s\")" ] }, { "cell_type": "code", "execution_count": 28, "id": "22d31fb0", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# make predictions on the test set\n", "e_pred = []\n", "f_pred = []\n", "for f in tqdm(frames_test):\n", " positions = f.get_positions()\n", " f.wrap(eps=1e-18)\n", " m = soap_fsparse_grads.transform(f)\n", " e_pred.append(model_fsparse.predict(m))\n", " f_pred.append(model_fsparse.predict_forces(m))\n", "\n", "e_pred = np.array(e_pred)\n", "f_pred = np.concatenate(f_pred)" ] }, { "cell_type": "markdown", "id": "b320a45b", "metadata": {}, "source": [ "Notice how it is already much faster to train and evaluate a model. Now, how accurate is it?" ] }, { "cell_type": "code", "execution_count": 29, "id": "dff52672", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSE = 26.000018450459176 meV\n", "Test set stdev = 0.274800153468654 eV\n", "Relative RMSE = 9.461427922173616 %\n" ] } ], "source": [ "from rascal.utils import get_score\n", "\n", "score = get_score(e_pred.flat, e_test)\n", "RMSE = score['RMSE']\n", "sigma_test = np.std(e_test)\n", "print(\"RMSE = \", RMSE*1000.0, \"meV\")\n", "print(\"Test set stdev = \", sigma_test, \" eV\")\n", "print(\"Relative RMSE = \", RMSE/sigma_test*100.0, \" %\")" ] }, { "cell_type": "markdown", "id": "3396f60c", "metadata": {}, "source": [ "In fact, it's slightly _more_ accurate than the model built with the full feature set. This might be because the model built with the full features overfits slightly, since it has more than enough features to describe the dataset. If we were to optimize the regularizers for both these fits, it is likely that the feature-sparsified model would have approximately equal or lower accuracy than the full model." ] }, { "cell_type": "markdown", "id": "b03fdbc8", "metadata": {}, "source": [ "## Timings" ] }, { "cell_type": "markdown", "id": "8c6c9c77", "metadata": {}, "source": [ "Now, back to our main motivation for doing feature sparsification. If we use these models to run molecular dynamics (MD), how much of a speedup do we get with the feature-sparsified model?" ] }, { "cell_type": "code", "execution_count": 30, "id": "f3c8abcf", "metadata": {}, "outputs": [], "source": [ "from ase.md import MDLogger\n", "from ase.md.langevin import Langevin\n", "from ase import units\n", "from ase.io.trajectory import Trajectory\n", "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution" ] }, { "cell_type": "code", "execution_count": 31, "id": "7f4c3bb8", "metadata": {}, "outputs": [], "source": [ "from rascal.models.asemd import ASEMLCalculator" ] }, { "cell_type": "markdown", "id": "b9b3ee9d", "metadata": {}, "source": [ "Let's use the ASE calculator (as in `MLIP_example.ipynb`), which we can run directly in the notebook." ] }, { "cell_type": "markdown", "id": "40d996ef", "metadata": {}, "source": [ "Full model:" ] }, { "cell_type": "code", "execution_count": 32, "id": "c29218ba", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elapsed time: 2.3851280212402344 s (4.770256042480469 ms per timestep)\n" ] } ], "source": [ "atoms = frames_test[0].copy()\n", "calc = ASEMLCalculator(model, soap_grads)\n", "T = 200\n", "\n", "start = time()\n", "\n", "MaxwellBoltzmannDistribution(atoms, T * units.kB)\n", "atoms.set_calculator(calc)\n", "dyn = Langevin(atoms, 0.5 * units.fs, units.kB * T, 0.002)\n", "dyn.run(500)\n", "\n", "elapsed = time() - start\n", "print(\"Elapsed time: {} s ({} ms per timestep)\".format(elapsed, elapsed / 500 * 1000))" ] }, { "cell_type": "markdown", "id": "26e48a5b", "metadata": {}, "source": [ "Feature-sparsified model:" ] }, { "cell_type": "code", "execution_count": 33, "id": "43afd338", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elapsed time: 1.2145779132843018 s (2.4291558265686035 ms per timestep)\n" ] } ], "source": [ "atoms = frames_test[0].copy()\n", "calc = ASEMLCalculator(model_fsparse, soap_fsparse_grads)\n", "T = 200\n", "\n", "start = time()\n", "\n", "MaxwellBoltzmannDistribution(atoms, T * units.kB)\n", "atoms.set_calculator(calc)\n", "dyn = Langevin(atoms, 0.5 * units.fs, units.kB * T, 0.002)\n", "dyn.run(500)\n", "\n", "elapsed = time() - start\n", "print(\"Elapsed time: {} s ({} ms per timestep)\".format(elapsed, elapsed / 500 * 1000))" ] }, { "cell_type": "markdown", "id": "1b80dfa8", "metadata": {}, "source": [ "So our feature-sparsified model is about twice as fast as the full model, and with the same accuracy!\n", "\n", "(But since we selected about one-sixth of the features, why isn't it 6x as fast? Well, there are various overheads, especially in the computation of gradients that are needed for MD forces, that keep the computational cost from scaling exactly linearly with the feature vector size. We could probably push the cost down even further by optimizing the feature size more aggressively, but if we want to keep the same accuracy we probably won't be able to get a speedup above about 3x. Note also that the impact of this optimization is much more pronounced for models with much larger feature sizes -- e.g. many atomic species or high `max_radial`/`max_angular` parameters.)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.13" } }, "nbformat": 4, "nbformat_minor": 5 }