Table of Contents
\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Optimizing Hyperparameters\n", + "\n", + "This notebook is intended to _purpose of tutorial_.\n", + "For more information on _nuances of tutorial_, please refer to (among others): \n", + "- [Title (Surname Year)](url) \n", + "- [Title (Surname Year)](url)\n", + "\n", + "Beyond libRascal, the packages used in this tutorial are: [pkg](url), [pkg](url), and [pkg](url)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "# may omit in final version\n", + "%reload_ext autoreload\n", + "%autoreload 2 \n", + "\n", + "import os\n", + "import time\n", + "import sys\n", + "\n", + "import numpy as np\n", + "import ase\n", + "from ase.io import read\n", + "from matplotlib import pyplot as plt\n", + "import json\n", + "\n", + "sys.path.append('./utilities')\n", + "from general_utils import *\n", + "from rascal.representations import SphericalInvariants as SOAP" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": true, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/tutorials/Introduction_to_Atomic_Descriptors.ipynb b/docs/source/tutorials/Introduction_to_Atomic_Descriptors.ipynb new file mode 100644 index 000000000..6c8182007 --- /dev/null +++ b/docs/source/tutorials/Introduction_to_Atomic_Descriptors.ipynb @@ -0,0 +1,128 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Introduction to Atomic Descriptors\n", + "\n", + "This notebook is intended as an introductory to calculating and understanding atomic descriptors, primarily SOAP vectors. For more information on the variable conventions, derivation, and utility of atomic descriptors, please refer to (among others): \n", + "- [On representing chemical environments (Bartók 2013)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.87.184115)\n", + "- [Gaussian approximation potentials: A brief tutorial introduction (Bartók 2015)](https://onlinelibrary.wiley.com/doi/full/10.1002/qua.24927)\n", + "- [Comparing molecules and solids across structural and alchemical space (De 2016)](https://pubs.rsc.org/en/content/articlepdf/2016/cp/c6cp00415f)\n", + "- [Machine Learning of Atomic-Scale Properties Based on Physical Principles (Ceriotti 2018)](https://link.springer.com/content/pdf/10.1007%2F978-3-319-42913-7_68-1.pdf)\n", + "\n", + "Beyond libRascal, the packages used in this tutorial are: [json](https://docs.python.org/2/library/json.html), [numpy](https://numpy.org/), [ipywidgets](https://ipywidgets.readthedocs.io/en/latest/), [matplotlib](https://matplotlib.org/), and [ase](https://wiki.fysik.dtu.dk/ase/index.html)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "# may omit in final version\n", + "%reload_ext autoreload\n", + "%autoreload 2 \n", + "\n", + "import os\n", + "import time\n", + "import sys\n", + "\n", + "import numpy as np\n", + "import ase\n", + "from ase.io import read\n", + "from matplotlib import pyplot as plt\n", + "import json\n", + "\n", + "sys.path.append('./utilities')\n", + "from general_utils import *\n", + "from rascal.representations import SphericalInvariants as SOAP" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Disclaimer: Dirac Notation**\n", + "\n", + "This notebook uses Dirac notation, which may be unfamiliar to many. For a refresher on Dirac notation, check out any of these resources: [Wikipedia](https://en.wikipedia.org/wiki/Bra%E2%80%93ket_notation), [GA Tech](http://vergil.chemistry.gatech.edu/notes/intro_estruc/node5.html).\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Why do we need atomic descriptors?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## What are some simple atomic descriptors?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## What can atomic descriptors be used for?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## What kind of atomic descriptors does LR calculate?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "celltoolbar": "Raw Cell Format", + "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.8" + }, + "toc": { + "base_numbering": 1, + "nav_menu": { + "height": "297px", + "width": "252px" + }, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": "block", + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/tutorials/Kernels.ipynb b/docs/source/tutorials/Kernels.ipynb new file mode 100644 index 000000000..f1107f248 --- /dev/null +++ b/docs/source/tutorials/Kernels.ipynb @@ -0,0 +1,90 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "toc": true + }, + "source": [ + "Table of Contents
\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Similarity Kernels for Environments\n", + "\n", + "This notebook is intended to _purpose of tutorial_.\n", + "For more information on _nuances of tutorial_, please refer to (among others): \n", + "- [Title (Surname Year)](url) \n", + "- [Title (Surname Year)](url)\n", + "\n", + "Beyond libRascal, the packages used in this tutorial are: [pkg](url), [pkg](url), and [pkg](url)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "# may omit in final version\n", + "%reload_ext autoreload\n", + "%autoreload 2 \n", + "\n", + "import os\n", + "import time\n", + "import sys\n", + "\n", + "import numpy as np\n", + "import ase\n", + "from ase.io import read\n", + "from matplotlib import pyplot as plt\n", + "import json\n", + "\n", + "wd=!(pwd)\n", + "sys.path.append('./utilities')\n", + "from general_utils import *\n", + "from rascal.representations import SphericalInvariants as SOAP" + ] + } + ], + "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.8" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": true, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/tutorials/Property_Prediction_Tutorial.ipynb b/docs/source/tutorials/Property_Prediction_Tutorial.ipynb new file mode 100644 index 000000000..4b361e9e5 --- /dev/null +++ b/docs/source/tutorials/Property_Prediction_Tutorial.ipynb @@ -0,0 +1,679 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Machine Learning Materials Properties\n", + "\n", + "This notebook is intended as an introductory how-to on training a model on materials properties based upon SOAP vectors. For more information on the variable conventions, derivation, utility, and calculation of SOAP vectors, please refer to (among others):\n", + "\n", + "- \n",
+ "
- On representing chemical environments (Bartók 2013) \n", + "
- Gaussian approximation potentials: A brief tutorial introduction (Bartók 2015) \n", + "
- Comparing molecules and solids across structural and alchemical space (De 2016) \n", + "
- Machine Learning of Atomic-Scale Properties Based on Physical Principles (Ceriotti 2018) \n", + "
- \n",
+ "
- training percentage: (float) percent of frames to train the model on \n", + "
- $\\zeta$: (int) order the dot kernel \n", + "
- $\\Lambda$: (float) error for the Tikhonov Regularization \n", + "
- jitter: (float) threshold for zero eigenvalues \n", + "
- kernel_type: (string, \"Atom\" or \"Structure\") whether to compare the environmental or average features.\n", + "
| Statistic | Value |
|---|---|
| Mean Average Error | 0.032 |
| Root Mean Squared Error | 0.0448 |
| R2 | 0.9178 |
"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "y = krr.predict(frames[test_idx])\n",
+ "krr.plot(y_known=property_values[test_idx], y=y, property_name=property_name)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Predicting from Another Data Set"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Our weights have been calculated and have shape (80,)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/markdown": [
+ "
"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAEUCAYAAADnbrk+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de5xVdbn/359hgBGESRIURC7GKVIEVBTqaBfsaFQinl91DDupnZPZRfR083qMvCSVHT10O5nd7ER2FUlDU9DqqKAiF1EpCaSQ8BLGzQYd5/n98V1rs2azL2tm7zWz98zzfr3Wa/Zet+8ze/Z65vs83+ciM8NxHKcSGrpbAMdx6h9XJI7jVIwrEsdxKsYVieM4FeOKxHGcinFF4jhOxbgicRynYkoqEkkHd5UgjuPUL+VmJKsk3S3p3yS9qkskchyn7iinSA4BvgQcD/xe0q2STpe0X/aiOY5TLyhtiLykfsAM4HTgrcASMzsjQ9kcx6kTUjtbzewl4HHgCWAH8PqshHIcp74oq0gkHSrp05IeAW6LrplpZkdnLp3jOHVBSdNG0v0EP8lPgJvNbEVXCeY4Tv1QTpG8Cfidea0Bx3FKUNK0MbPfmplJeq2kJZLWAkiaKOmyrhHRcZxaJ62z9VvAxcDLAGa2hrB64ziOk1qRDDCzB/P2tVZbGMdx6pPGlOc9L+k1gAFIejfwl8ykqpADDzzQxowZ091iOE6PY8WKFc+b2dD8/WkVyceAG4Dxkp4GNgI1G4w2ZswYHn744e4Ww3F6HJI2FdqfSpGY2QbgbZIGAg1mtjPv5mea2fcrF9NxnHqkQ2UEzGx3vhKJOL9K8jiOU4dUqx6JqnQfx3HqkGopEg9Yc5xejM9IHMepmGopkvuqdB/HceqQVKs2kvoD/w8Yk7zGzK6Ifn48C+Ecx6kP0saR3ApsB1YAe7ITx3GceiStIhlpZm/PVBLHcarKguWbmL90PXOmj2P21NGZjpXWR3K/pCMzlcRxnKoyf+l6tm5v4StL12c+VlpFcjywQtLvJa2R9KikNVkK5jhOZcyZPo7hzU2cN31c5mOlNW1mZCqF4zhVZ/bU0WVNmmqZP+UaZA2OXu4ssjmOU8dUy/wpNyNZALyLsFpjtA88M+CwikZ3HKfbWLB8Ey/uaaW5qbFi86ekIjGzd0U/x5Y6T9IRZvZYRZI4jtOlzF+6nh0trQxvbqp4Vadaka0/qNJ9HMfpIqrpjE3rbC2H59o4To2Q1oGaxhmbFs/+dZweRlfGj8RUS5E4jlMjdGX8SEy1TJuXqnQfx3EqpJomS1pSKxJJE9k3+/cX0c9pVZfMcZyq0BU5N2nLCHwHmAg8BrRFuw34RSZSOY5TNfJ9JlkolbQzkmlmdnjVRnUcp8uYM30cX1m6nvOmj2unVLpDkTwg6XAze7xqIzuOU1WKmTD5PpNYqVSTtIrkJoIy2UoobCTAzGxiVaVxHKfTpJltZOWITbv8+23gX4G3A6cQ8m9Oqbo0juN0mjnTxzG4qZHde1pZsLxgQ7zMSKtInjOzRWa20cw2xVumkjmO0yFmTx3NgP6N7Ghp7dJgNEhv2qyUtAD4JYmarfHyr+M4tUHSsdqVpFUk+xEUyEmJfb786zg1RncEo0H6JuJnZy2I4zjpiVdoDh7UnzVPb2fmpBFcf/pR3SZPKh+JpJGSbpH0bLT9XNLIrIVzHKcw8QrNqs3baTNYtHpLt8qT1tn6XWARMCLafhntcxynG4gT8yaPbKZBMHPSiG6VR2blKwBIWmVmk8vtqxWmTJliDz/8cHeL4ThVoyt71JRC0gozm5K/P+2M5K+S3i+pT7S9H/hrdUV0nN7LguWbmHbNkqLxH8VqjJS7rqtIq0g+CLwX2Ar8BXg3cFZGMjlOryNWFF9YvK6gYpg2dggNgqljhxS8rqvjRvJJq0hGmtlMMxtqZsPMbBYwKkvBHKc3Efs8DAoqhmUbt9FmsHzjtty+uAr84CpUga+UtD6SR8zs6HL7agX3kTj1yoLlm9oFlMV+EWCf/X/bvYeWVmNwUyNr5p7cJfIV85GUjCOR9AbgjcBQSZ9IHBoM9KmSYJ8ErgWGmtnzBY6fCVwWvb3KzL5fjXEdpxZJBpRNu2ZJbnbywMUn7rM/rrheC5XXy5k2/YD9CQpnUGLbQfCTVISkQwnRsn8qcnwI8FlgKnAc8FlJB1Q6ruPUEoUcphfcvJKt21uAff0isRl06uQRDG9u4sIZ47tU3kKUa5D1G+A3kr6XUZLedcBngFuLHD8ZuMvMtgFIuouQgfyjDGRxnC4hfym3UAWzZyIlAu39ItB9YfClSOtsfVHSlyT9StLSeKtkYEmnAk+b2eoSpx0C/DnxfnO0z3HqlnzFkaz6Hh/r3ygENDWK86aPq5ll3mKkTdr7IfBjQh2Sc4EzgefKXSTpbuDgAocuBS6hfRJgRUg6BzgHYNQoX1ByaovkLGTO9HHMW7yuXd2Q3Xtambd4HdPHD2P5xm2clxd4lvSX1NpsBNKv2qwws2MkrYmrokl6yMyO7dSg0pHAEuDFaNdIYAtwnJltTZz3PuAtZvbh6P03gXvNrKRp46s2Tq0RK4LhzU08cPGJufeDmxrZtaeVtugxjI/nk1zNqefI1pejn3+R9E5JRwFDSl1QCjN7NIpHGWNmYwgmy9FJJRJxJ3CSpAMiJ+tJ0T7HqSuS5ksy/gOgzcLKS6l4kNlTR7dbuak1UyetaXOVpGbgk8BXCMu//5GFQJKmAOea2b+b2TZJVwIPRYeviB2vjlNP5C/r7mhpzSmWzsw0sqoG31lSmTZlbyJdbGbXVEGequCmjVPLVMNM6S5Tp1MBaR3gPUDNKBLHqVUWLN/EvMXrKr5PrS0BV6uJeC0E1zlOl5Dvn+iIv2L+0vXsaGntlgLNWVItRVK5feQ4dUKhALI0Kf6xk7WpUQxuamTq2CE15TCtBJ+ROE4HSa7A5L9PKo/80gDzFq9jR0srBwzsz5q5J7Ns47aaKAFQDcoqkqiQUbkVmp9WSR7HqXnyl2Jh75R83uJ1bN3ewrzF6/YpDQAUVUD1Tllnq5m9EgWGXVfinM9XVSrHqSPimcelt6zNKRRBTtHMW7yO5qZG3jp+GMsSeTO15jCthLSmzX2SvirpBElHx1umkjlOnTBn+jgaRDslYuwNi9/R0sqA/o09ypTJJ+3yb1zk+YrEPgOmV1ccx6k/4lnFFbc9TsvLbfTv25BblcnvfNcdXfC6gqoEpNUaHpDmdAdx/gxAc1MjF84Y32NMl5iKcm0kHSTp25IWR+8Pl/Rv1RbSceqZaYkCRAP6N/Y4JVKKtKbN9wgNsS6N3v+BUFbg2xnI5Dg1S7IcAJCLUr1oxniWrnsWCD6SeCl4/tL1TBs7hGUbt3V7T5osSVtG4CEzO1bSSjM7KtrnDbKcXsWC5Zu4bOFa2ox9lnUF9O/bQMvLbTQ3NbJ67sk5U6dBIcO3uamR/fo31rVCqbSMwG5JryZyTEuaBmyvonyOUxMUq5962MW3c8Vtj9Nm0KAw45iTcJoa0K9PQ66GahzF2tzUyMRDQlvNPa2v9NhVm7SK5BOE3r+vkXQfcBNwXmZSOU43USjcfdHqLbQZtLzcxvDmJq6aNYHZU0fzYF4t1YtmjM8FqiWXfbfu3EObQf/GPgUjYHsCqRSJmT0CvJnQmuLDwBFmtiZLwRynOygU7h7PKCaPbMaABzduY+LcO1m4akvuusFNje2KDsWzkfOmj8t1yXvr+GHtFE1Pmp2kXbVpAuYAVwKfAz4W7XOcHkUy/D1+2J/ZuYcN17yTrTv3sHV7C7eu2sKOltZ2112UaAmRnI3Mnjq6YJe8nhQeD+lNm5uAIwjV0b4avf5BVkI5TpakMSvy22HmMnf7NuyT6j55ZHM752mppL6YQvk69Uza5d8JZnZ44v09kh7PQiDHyZo0ZQrjWcXw5iZmTx2dK4/YkJfnPnlkMws/fnzufX7PGuhZOTXFSDsjeSRaqQFA0lTA11eduiSNWTFn+jiamxrZvaeVC25emZudTDykOddv5vOnTWinRGCvkrps4doe40hNQ1pFcgxwv6SnJD0FPAAcK+lRSe50deqK2KwAipo4s6eOxoAdLa0sjHwiO1paWffMToywAjN76uh2ZlJs/kCIG+kpjtQ0pDVt3l7qoKQDzOyFKsjjOJlRqlVmWtOj5eU2YG+mb/IeseJpbmpkQP/irSV6IqkUSbm+v5IeAbysgFOzJKNSY8WRn5mbfz6QqyNyz7pnQ7p7ohMeUDS7t6f7RPKpVhV5L7Xo1DTzl65vF5UK7Z2g+bOVuCxic1MjxxXIlbng5pVctnAtMyeNaNcZr7cpkBgv/uz0CuZMH0dT3wbajH0iUqG9ibJg+aZcnEjckzeugBbPVOJo10Wrt+xzr95ItWYkjlPTzJ46mssWrgXg1lVbWLZxW7us3KSJkuw702rklIqx14HarzEk6M2cNKLLf5daxKvIO72GmZNG0CDo3yi2bm9h0eot7ZythQLE+vbZ+9Vu6tvAedPHMX/p+lzezfWnH9XVv0ZNkjZE/suSjihxyr7t0x2nm8mPYD1u7BCGDW7i7ROGM7y5KZdDc9Cg/ky7ZgkX3LySadcs4bADB+busV/fPrkgtAMG9Ms5aXtSeHs1SDsjeQK4QdJySedGDcVzeGNvpxYp1shq+cZtPHDxibms3NWbt7N1ewsLV4UZyqrNoUJGPBeZOWlEu0S++UvXM3XsEOZH/hQnffbvjWb2j8AHgDHAGkkLJL01S+EcpxLyM3mTuTPJ430KGOaxCbSjpZV71j27TyJf0iwqRE8rE1CO1D4SSX2A8dH2PLAa+ISkmzOSzXEqIj+Td0dLKwPzaqka0Ni49zGYPLKZ4c1NzJw0gj2tljsnJlY+MyeNYHAUQl9IWfS0MgHlSOsjuQ74PfAO4PNmdoyZfcHMTgHc2+TUDPkh6/HruCbI1KhAcxygtnV7C/36NOT8IM/s3MMDF5/Iso3bMMLMJFkiIFZO159+FAP6NxZtBt7b/Chpa7aeDfzEzHYXONZsZjVVdtFrtvZe4jqpyZqq+a/Pmz6OS25Zm7vm86dNAPYWcp4+fliukPNFJVpKLFi+qddFshar2Zo2jmQ18DqpnTG5HdhUa0rE6d0UClmfOnYIS9c9m/OPzE/MIJoTlc2S/o+4wHMpBdEbygOkJa0i+Tohl2YNwZk9AXgMaJb0ETP7dUbyOU6HyH+4k7VEmhrFpbespY9CGYD+jX24MGG2TBs7hEWrtzDxkGae2bmn15gl1SCts3ULcJSZTTGzYwh+kQ3APwFfrEQASZ+UZJIOLHBssqQHJD0maY2kf6lkLKf3kPSPxLVFWloNI0SrHjCwPxfOGJ9bwl2wfFNuJvLMzj1MHTuEyxau5YKbV3b3r1IXpJ2RvNbMHovfmNnjksab2YY8c6dDSDoUOAn4U5FTXgQ+YGZPShoBrJB0p5n9rdODOr2C2Ez5wuJ17Nc/rK4kiU2cZAmAZFJfnCm8aPUWj15NQdoZyeOSviHpzdH29Whff+DlCsa/DvgMRZL+zOwPZvZk9HoL8CwwtILxnB5AfoxGoZiNeNUk7iXTmviGzZo8ItcBL3a+Jpd15y9dn4t69VyadKRVJGcC64ELom0DcBZBiXQqKE3SqcDTZrY65fnHAf2APxY5fo6khyU9/Nxzz3VGJKdOKBaxGr+Po08PGtSfltb2/6OaGsWyjdtyTtV4xSVe1o2PxZXjfTaSjrKKJApEu9HMvmxmp0XbtWb2opm1mdmuEtfeLWltge1U4BLg8jRCShpOqFp/tpm1FTrHzG6IfDhThg71SUtPplCV9uZEcFisWOJQdwgKZHBTI5efckSuQ16b0S7TN3nvqWOHZBKZ2lMjXtPGkfwfMN3MXqrKoNKRwBKCDwRgJMGhe5yZbc07dzBwLyEQ7mdp7u9xJL2POH5kcFNw+730SluuLKIItvPw5qZcEaKJc+/MFS5aPffkovdLXlNNOat9366i0t6/G4D7JP2npE/EW2eFMbNHzWyYmY0xszHAZuDoAkqkH3ALcFNaJeL0TuKZBIT6IbESgaBEmhrVbjn3ohnjc316S92v2kvAPTXiNe2qzR+jrQEYlJ04IGkKcK6Z/TvwXuBNwKslnRWdcpaZrcpSBqd7KdQbJu15+aZKzEuv2D73KjUXzyrYrKcGsaUybXInSwPM7MXyZ3YvbtrUN/H0v7mpkf36N7arZJZ8COPzGgRXzZqQ843ECJg0spk1T29n5qQR7RynlZgYaRVdT6Qi00bSG6LOeuui95OiJWDHqSq51phRCn+xlP34PBGcpp9d9Bhbt7fQqLC826Aw41jz9HaumjVhn9WXOdPHlczeLUVvy+xNQ1ofyfXAycBfAaIl2zdlJZTTe4nT/V96xXLZt8nCQvnnAQxuauTlV8LMutXg+tOP4qpZE2hQ8UZVs6eOzmXvzlu8rkMrKT3Vz1EJqYs/m9mf86JYX6m+OE5vJ066mzp2SK5/TCHzYU4i+nRg/0YOO3AgqzZvZ/LIULwvvqZY35rkWLv3tHaoUVZP9XNUQtoZyZ8lvREwSX0lfYpQftFxCtLZeIlkvY9ybTX79REi1BhZ+PHjeWreO9v14i3XmjM+Hq/gFFM4PTX2o5qkVSTnAh8DDgGeBiZH7x2nINXyI8Q9Zb6Qtxozb/G6XBJeXDuks7IUqyCf9nonfc3W583sDDM7KIr/eL+Z/TVr4Zz6pdp+hHhtcdZX/48xF92e849AiBspNVuoVBb3iZQnbWTrUOBDhMLPOb+KmX0wM8kqwJd/a4dKlkoXLN/EvMXrePGlVlrbwmrMwlWFO9vVa6RovVFpZOutQDNwN3B7YnOcklRiFsQrM61RkGq+Epk8spnmpsZ2leHBfRrdQdpVmwFmdmGmkjg9kvzSh2mIZzHTopWbgwb1b5eAJ+Dq0yak8mnkn9Obg8myJO2M5DZJ78hUEqdHUs6RWYj8RlbvPfbQdj1hBzc1lrxfqQxed5xmQ1pFcj5BmbRI2iFpp6QdWQrm9A5KFSWKZzHzowpmEIo1F0u0i+83P5oBxbVFLlu4Nnf//KZZbgJVhw7l2tQL7mytH4rlvCRNECBV24cFyzdx6S1rMfYqnDhorZAztt5T+ruDitpRKIS0ngGMNbMro1qrw83swSrL6fQgSvkj8v0gsRkSK46478wVtz3OuitnpDKNkjMXo3B0a1KmzvhvnMJ0pB1FGzAduBLYBXwNODYjuZweQDGnZ3LmcM+6Z1k99+Tc7CAOWY9J1hUpRTLZL9lmItmzJl+mjvpunOKk9ZFMNbOPAS0AZvYCoX6q4xSlWCBX/swhee7UsUPaBZvNmpyu+HK8VHzAwP6snntyOwWRVB4eXJYNaRXJy1HtVoNcgFq6fxVOr6XYik2cwt/c1Mj08cOYds0SHty4jd17WnOxIg1ROYBlG7elcoaWUhDJY51ZRXLKk9a0mU8oeThM0tXAu4HLMpPK6dEkzY071v6Fllbj1lVbcrOT/EJFabJyS2Xk5ps3rkSqTypFYmY/lLQCOJEQDzTLzHLZv5IOiMwdx0lFrCTi+BADmvo20K9PQ7vG3dVyhnZEKTkdpyP1SNYRVUgrwBJCb2DHSUW8YpKMWj1gQL92y7DVrPtRaoXGo10rJ7UiKUPn+3Y6vZL4gb1sYVjmjVtlQjYPdiGlFI/z4p5WdrS0+mylAtI6W8vR86LanKpSKIp0/tL1uX67V82a0M6X0RVh7EnzyldyKqNaisTpZXQ0vDzZ1Du+Ll5NuWrWBIB99uc/2KXG7Ey4ezzOhTPG+0pOhVQlRF7SSjOrmSapHiKfPWnDy/MjWHdHZkT+LCTN/UqdE3fOG9zUyJoCnfO6ip7ub6m0HgmS+kgaIWlUvCUOe6JCLyH+zz9t7JBU5kA8E7l11RZ272ll+vhhueru8wrMToo5Q9OO2d3Out6aXZw21+Y84LPAM+wNRDNgIoCZbctEOqfmyE/xL0ey2vuOllYWrd7CzEkjWL5xGy/s3pNrB7EmLxq1o2NeNGN8TeTN9Nb8nY6UEXidmR1hZkdG28QsBXNqk46m4c+eOpqZk0bkZgptRk4h9GvsA5SfRaRpZlUrEau1IkdXk7odBbC97FlOj6GYkkg+KGmn8cs2bss18m5QaB8Bext5vzUKk0+OlRw/2cyqt5kM9UJaRbIBuFfSxZI+EW9ZCuZ0L2mURLkEuAtuXslhF9/OwYP6M7y5iX6NfXIzEtirlJaue5at21u49Ja9BYji8eOiRJ5sV9ukDUj7U7T1w7N+ewVpbP1ikafxysUz21swYNXm7Xz+tLDEW+qeRnDAxqs8i1ZvybXc7I3mQj3RoeVfSfsDmNmuzCSqAr78273Ey7RJSi3rLli+iS8sXpeLatzR0pqbfaSpjOZ0HRUt/0qaIGkl8BjwmKQVko6otpBO7dCZAK/kMm1DwoOaDH/Pv3c8e7lwxnjWzD25XfvMYo5Lr7Vae6RtkHU/cKmZ3RO9fwvweTN7Y7bidQ6fkVROR+uZLli+qV191POmj2Pe4nUIuDCRzZt/nkHRcYoFd3mt1e6jopqtwMBYiQCY2b2SBlZNOqdmyI9EzZ9JFIvanLd4HW0WlnLj2UQhcySZXxPfu5DfJKlw4sS++H69NVajlkmrSDZI+k/gB9H79xNWcipG0ieBa4GhZvZ8kXMGA48DC83s49UYtzeSJny7VPBXmpoeg6KeM0mFtGzjtnY/YwUV36OUwgFyDtfk+bXsM+npYfKFSLv8+0FgKPCLaBsa7auIqBr9SYQVoVJcCfy20vF6O5Uu6eYfW7B8ExPn3snEuXcyffwwhjc3cVFUdDkZGp/8GSuocg9Y7GeZPLK57pZ9e2OYfNoKaS8AczIY/zrgM4TewgWRdAxwEHAHsI9t5qSnkiXdQsfigsvAPjOYeKwXdu+hpdXo3ygOGNg/tUJYtnEbbQbP7NxTd36Q3mh6lVQkkq43swsk/ZICNUfMbGZnB5Z0KvC0ma0ObXMKntMAfJlgSr2ts2M5gWqaBPntH/IfmnisBcs3dWoJt54fxlo3vbKg3Iwk9olc25mbS7obOLjAoUuBSwhmTSk+CvzKzDYXUzaJsc4BzgEYNWpUyXOdyolnI+VWTjr7UPXGh7GeKalIzGxF9HKymf138pik84HflLm+4CxC0pHAWCCejYwEHpF0nJltTZz6BuAESR8F9gf6SdplZhcVGOsG4AYIy7+l5HLSUchpeMHNK1m4aguNDaEtZrEZQzUcjvltO3ubA7OeSOtsPbPAvrM6O6iZPWpmw8xsjJmNATYDR+cpEczsDDMbFZ3zKeCmQkrEyYZCTsNFq0PfmdY2GNB/7wpNoTKK8bWdDSBL3mPe4nVs3d7CvMXF6o873UlJRSLpfZF/ZKykRYntHiCTGiSSpki6MYt7Ox2j0ArOzEmh811T34bc/kIKJ3ltmlWMQsqm0PjdXbjIKUzJyFZJowkmyDVAciawE1hjZq0FL+xmPLK1aynnUE3jcC0XrdpZp61TXYpFtlalZmut4YqkNLUYMOWKoj6oNGlvmqSHJO2S9JKkVyTtqL6YTldQrYCp/OS7ShLpemtlsZ5CWmfrV4H3AU8C+wH/DnwtK6GcbClXJCitUkgqpN4YzensJXUVeTNbD/Qxs1fM7LvA27MTy8mS2VNH55yghZRFKaWQVDKxQpo6dggv7gmtIOoxgMypnLSK5EVJ/YBVkr4o6T86cK1Tg5RSFqVmLPmJew9cfCLLNm5jR0srA6Pl4GJ4HZGeS1pl8K9AH+DjwG7gUOD/ZSWUkz2llEUhf0Wx3jJxqHyp4LQYN396LmmT9uJ/IX8HPpedOE5XUaqpdjKSNE7/jxtt5yfnJUPl02T0Llq9JVdF3uk5pF21eZeklZK2Sdohaaev2tQ+lfTnvWzhWrZub2HR6i25+quFZjAdqe4eZ/TGVeSdnkNa0+Z6Qpj8q81ssJkNMrPBGcrlVIFipkQxBTNn+jiamxrZ3tKaq2I2c9KIXJ2R2NyJ20xccPPKDi3bekuJnktHGmSttZ4YvdaDKfbgFlMws6eOZr/+wdqNm3xff/pR+yiKuE1EnHeTFo8V6bmkLbX4GeBXkn4D7Il3mtl/ZSKVUxXiB3Z+pDDS1DxNHiv2wM+cNCLXw9dxIL0iuRrYBTThDbLqimJ1VotNLdPUAbn+9KO4/vSjqiilU++kVSQjzGxCppI4HSJtvkz+7CNZnb1UEWfH6QhpfSS/klSumpnThaSNycj3SxRqB1EJHmTmQHpF8hHgDkl/9+Xf2qCjKyD5AWVXzZpQldmIB5k5kMK0UaiFeISZlWsZ4XQhHalpGpdHBFi67lkG9E9r0Zannos0O9Wj7IwkWvK9vQtkcTIif5k2OYPIN006aqr4kq4D6U2bRyQdm6kkTtWJlcLEQ5ppEMyaPKJdk25gn1qobqo4nSHtHHcqcIakTYSkPREmKxMzk8zpNPGKTpwfI2DDNe/MHS80e4hrobqp4nSGtIrk5EylcKpKPKsY3NS4j0M2f9n4ohnj2ykO7yfjdIbU2b+SJgEnRLt+Z2arsxPL6SxxWv/gpkYumjF+H6VQKEDN8x6cSkmb/Xs+8ENgWLT9r6TzshTM6RxxWn+xIkP5y8buE3GqQVpn678BU83scjO7HJgGfCg7sZzOkix/WGj1JX+VxTNynWqQVpEIeCXx/hW8V1GmdDZiNFn+sDORr47TGdIqku8CyyXNlTQXWAZ8OzOpnE6bHMVKIjpOlpR0tkoaa2Ybzey/JN0LHB8dOtvMVmYuXS+ms8uwsQLKL4noOFlSrmXnCjM7RtISM6ubb2Vv67SXX2vVO9Y5WVGs01655d8GSZcAr5X0ifyDXtioNkiaQe7vcLqDcj6S0wmO1UZgUIHNqQG6cuXFywY4hUjVRFzSDDNbXOL4mWb2/apKVgH1aNrUYmPvQky7Zglbt7cwvLnJfTC9kIqaiJdSIhHnd0oqJ0elgWFdNVNIO/vxmUvvolptNz2mpEJKPaBpHsquilBNG3fiEbO9i2opEk/XqJBSD2iah7KYIuqumYFHzPYuUvlIyt5EWmlmNVNWvB59JKVYsHxTp5d03afhVJOKfDfXvmcAABBcSURBVCSSxpbZd18Fgn1Skkk6sMjxUZJ+LekJSY9LGtPZseqVSsLYfWbgdAVp65H8HDg6b9/PgGMAzOzjnRlc0qHASUCperA3AVeb2V2S9gfaOjNWb8XrizhdQbkQ+fHAEUCzpH9OHBpMaJZVKdcRuvjdWmT8w4FGM7sLwMx2VWFMpxPUy/K00z2UM21eB7wLeBVwSmI7mgrLCEg6FXi6TIGk1wJ/k/QLSSslfUlSn0rG7Yl0hUPVV2GcUpQzbd5oZmdLutzMrujozSXdDRxc4NClwCUEs6acfCcARxHMnx8DZ1Eg81jSOcA5AKNGjeqoqHVNsbac1cRruTqlKDcjeUfU12ZWZ25uZm8zswn5G7ABGAuslvQUMJJQqT5f6WwGVpnZBjNrBRayr68mHusGM5tiZlOGDh3aGXFrmlKzjq5wqHrdEqcU5WYkdwAvAPvnddaLq8gP7sygZvYooWRjuFlQJlPM7Pm8Ux8CXiVpqJk9B0wHanZdN0s/Qv6sI38sf8Cd7qTkjMTMPm1mrwJuN7PBiW1QZ5VIOSRNkXRjNP4rwKeAJZIeJSiwb2UxbjXI0o/gtVadWqYqAWm1RncFpFUSOFbLYzlOTLGAtHKFjXZSIvw9q1lJpfS0yFbHqRU6VdjIzAZFF18J/AX4AcG8OAMYnoGcjuPUIWmT9maa2dfNbKeZ7TCzbwCnZimY4zj1Q1pFslvSGZL6SGqQdAahB7DjOE5qRTIbeC/wTLS9J9rndBFeKMipZdJWSHvKzE41swPNbKiZzTKzpzKWzUngy71OLVOtwkZOAao5i/ByAE4tk7aMgNMJqpkD49GrTi3jM5IM8VmE01soV49kn6ZYSbxBVml8FuH0FsqZNnETrNcBxwKLovenAA9mJZTjOPVFucjWzwFI+i1wtJntjN7PBW7PXDrHceqCtD6Sg4CXEu9fivY5juOkViQ3AQ9KmhvNRpYDNdOis9bw4DGnt5E2IO1q4GxCkaMXgLPN7PNZCtYVZPXAe/CY09voyPLvAGCHmf03sLlQr5t6I6sH3pd9nd5GqoA0SZ8FphBWb74L9AX+F/jH7ETLnqwKGvuyr9PbSBvZehqhkvsjAGa2RdKg0pfUPv7AO051SGvavGShlJoBSBqYnUi1iztRHacwaRXJTyR9k1DR/UPA3cCN2YlVm7gT1XEKk3bV5lpCr9+fE/wkl5vZ/CwFq0Xcieo4hUnrbP2CmV0I3FVgX6/BfSqOU5i0ps0/Fdg3o5qCOI5Tv5TL/v0I8FHgNZLWJA4NAu7PUjDHceqHcqbNAmAxcA1wUWL/TjPblplUjuPUFeVadm6ParP+N7DNzDaZ2SagVdLUrhDQcZzaJ62P5BvArsT7XdG+Ho/HjjhOedIqElmit6eZtdFL6r167IjjlCetItkgaY6kvtF2PrAhS8FqBY8dcZzylGwinjtJGgbMB6YTwuSXABeY2bPZitc5vIm442RDsSbiaSNbnzWz081smJkdZGaza1WJ9Cbcf+PUCuXiSD5jZl+U9BWihL0kZjYnM8mcslSzb47jVEI5h+kT0U+3E2qQrOqpOE5HSeUjqTfcR+I42VDMR1LOtPklBUyaGDObWQXBPglcCww1s+cLHP8i8E6CP+cu4HzridrPceqYcqbNtdHPfwYOJpRXBHgf8Eylg0s6FDgJ+FOR428klHOcGO36P+DNwL2Vju04TvUo1yDrNwCSvpw3nfmlpGrYDtcBnwFuLSYC0AT0A0SoFVuxAnMcp7qkDUgbKOmw+E1UQb6icouSTgWeNrPVxc4xsweAe4C/RNudZvZEoXMlnSPpYUkPP/fcc5WI5jhOB0kb5v4fwL2SNhBmBqOBD5e7SNLdBJMon0uBSwhmTanrxwGvB0ZGu+6SdIKZ/S7/XDO7AbgBgrO1nGyO41SPVIrEzO6Q9A/A+GjXOjPbk+K6txXaL+lIYCywWhIERfGIpOPMbGvi1NOAZWa2K7puMfAGYB9F4jhO95HKtJE0APg08PHIFBkl6V2dHdTMHo2iZMeY2RhgM6FJ+da8U/8EvFlSo6S+BEdrQdPGcZzuI62P5LuExuFviN4/DVyVhUCSpkiKK9T/DPgj8CiwGlhtZr/MYlzHcTpP2qS9h81siqSVZnZUtG+1mU3KXMJOIOk5IJmAciCwT4xKjeKyVp96kRNqX9bRZjY0f2daZ+tLkvZjb4Os1wBlfSTdRf4vGivC7pKnI7is1ade5IT6kjVJWkXyWeAO4FBJPyQEiZ2VlVCO49QXZRWJwrLKOkJ06zTC8u/5hcLZHcfpnZRVJGZmkn5lZkcCt3eBTFlwQ3cL0AFc1upTL3JCfcmaI62z9fvAV83soexFchyn3kirSNYB/wA8BewmmDdmZhNLXec4Tu8gbRzJycBhhJqtpwDvin7WNJI+KckkHVjk+BclPSbpCUnzI39Qt5BC1lGSfh3J+rikMV0rYTtZSsoanTNY0mZJX+1K2fJkKCqnpMmSHoj+/msk/Ut3yJiQp9zf/0xJT0bbmV0tXznK1SNpAs4FxhGCwr5tZq1dIVil1FOJgnKyRtwEXG1md0naH2jrEuHySCkrwJXAb7OXqDAp5HwR+ICZPSlpBLBC0p1m9rcuEzIixXd1CGHldAohBGOFpEVm9kLXSVmacjOS7xOEf5TQNPzLmUtUPeISBcVst2SJgv50b4mCkrJKOhxoNLO7AMxsl5m92IXyJSn3uSLpGOAg4NddJVQBSsppZn8wsyej11uAZ4F9Aq26iHKf6cnAXWa2LVIedwFv7yrh0lBu1ebwaLUGSd8GHsxepMpJligoZq2Y2QOS4hIFIjiTuzyPJ42swGuBv0n6BSHZ8W7gIjN7pYvEBNLJKqmB8A/n/UDBpM2sSfmZJs8/jvAP5Y9Zy1Zg7DSyHgL8OfF+c7SvZiinSF6OX5hZaze6EPahK0sUVEqlshL+TicARxGmvz8mBAR+u3pSBqog60eBX5nZ5iy/L1WQM77PcOAHwJlRB8mqUy1ZaxozK7oBrwA7om0n0Jp4vaPUtd21AUcSpqlPRVsr4eE7OO+8TwP/mXh/OfCZGpV1GvCbxPt/Bb5Wo7L+MNr/FCFnZAcwr9bkjM4dDDwCvLvGv6vvA76ZeP9N4H3dIXPR36W7BeiCP9ZTwIEF9v8LwURoJPhHlgCn1KisfQjZz0Oj998FPlaLsuadcxbBZKw5OQmmTNwxstvkSynrEGAjcEC0bQSGdLe8yS3t8m+PoJ5KFCRlteAL+RSwRNKjBJ/Ot7pTviR5n2vNkifne4E3AWdJWhVtk7tRvHbk/f23EVbBHoq2K6J9NUOP7GvjOE7X0qtmJI7jZIMrEsdxKsYVieM4FeOKxHGcinFF4jhOxbgicRynYlyRFEDSXEmfkjQ+ii9YKekYSR9Nce2XotT0L3WFrNGYZ0UZrPH7G6NEPycDJF2S8f1nSrqoive7R9LJefsukPSNqo3hcST7ImkusIsQ9dpoZldF9T9uM7MJZa7dTog6TJVQJ6nRKizNIOle4FNmVo3G7l2GpD5pP6cMZejw5y9pl5ntn5VM1UbSOcAbzOzsxL5lhJSQ6pR66O7Q2lrZCAlUfyDUJfkRIfdmK6EZ2D3AzcDfgVXAl4rcYxEhP2kVIQR/DLAUWEMIxx4Vnfc94H+A5cB/AXMJJRt+R+jH88/AFwlRt3cAfaPrLidENq4l1PYU8G6C0vt9NO5+hJoqU6Jr3hfdZy3whYSsu4CrCVG9y4CDSnw2Q4Gfszey8h+j/XOB70TjbQDmJK55PyFbfBUhN6RPYtwvR+MeD7yDUFx8BTAfuI0wU36SvSkBDcD6+H0B+eLP8+Hob/iuaH8f4EuRzGuAD0f73xJ91ouAP5T4vRdGcj0GnBPtm5f4G/8w2veJ6PNdSxRyH/3t10Wy/YGQg/Q24L7odzuuxLhnEaUWRNfPB+6PPuOSeUGEHLL49/1ctG8IIaenX0K2PxFNJKry/HT3A1wLG3BM9LANICRyrSeEpM8l/KePP/y1Ke61K/H6l4SsUoAPAgsTX47bEg/XXIIC6wtMIhTdmREduwWYFX8hEvf+AVFuEAnFkXwPjIi+MEMJs6uliXtZ4vovApeV+J0WAMdHr0cBTyTkvp9Qz+VA4K/R7/D66HePFeDXCUWE4nHfG71uIqTHj43e/4gw64NQyCd+KE8Cfl5Cvu8RFG4DoSTo5uje58S/VyTjw4QyDG8hlAwdW+ZvOST6uR9BSby6wN84/u4MBPYnKJ2jou9LKyExr4GgkL5DUP6nxt+FIuOeRXtF8tPoHocD60tcdxJ7/8E0EL5jb4qO3QacGr2+CLi2ms+Q+0gCJwC3mNmLZraD8J+qGryB8BBCePCPTxz7qbWf1i82s5cJX8o+hAeD6P2Y6PVbJS2P8m2mA0eUGf9Y4F4ze87C9P2HhPwSCC1Yb4ter0iMUYi3AV+VtIrw2QyOqrQB3G5meyy0J3mWUNDoRMID9lB0zYmEUp0Q/pv/PHo9HthgZhuj9z9KjPkd4APR6w8SEhVL8RMza7NQrGhDdO+TgA9EMiwHXk1QNAAPJsYtxhxJ8Yzt0MS1SY4nfHd2W2h2/wvC9wlgo4U+120EBbPEwpOc/JumYWH0uz1O+HyLcVK0rSRkNY9PyPwj4PTo9em0/6wrJm2DLKf67M57vwfAzNokvRx94SCUVGyMyl5+nTDz+HPkx2mqYPzkGK9Q+rvQAEwzs5bkzqjeSLLjYnwfAd83s4sL3KvFUvhFot/xGUnTgeOAM8pdUuC9gPPM7M48ud/Cvp8/Bc55G8G38GLkh+ro5538bNoS79vo2LOXvE+pIi8CrjGzbxY4ditwnaSjgQFmtqID45fFZySB3wKzJO0naRCFC1vvBAZ18L73s/e/wBkEu7yzxF/i56PZwLtTyPYg8GZJB0rqQ/CX/KYTY/8aOC9+kyJLdgnwbknDovOHSBpd4LzfA4dpbyHr/ALMNwL/y76zt0K8R1KDQjvZw6J73wl8RFLfSI7XShpY5j4xzcALkRIZT6gJE/NyfE/C33SWpAHRvU+jsr9zJdwJfDCeLUo6JP4bRLOlewgzvarORsBnJACY2SOSfkxwAD5LcFbln/NXSfdJWkswQz6d4tbnAd+V9GngOeDsMueXkvFvkr5FsNW35sn4PeB/JP2dYE7F1/wlWka8h/Df6nYzu7UTw88BviZpDeE781tCUfBisj4u6TLg11HpxZeBj9G+sTtm9vdoSf0OSbvZ93NfRDBpypk1EHxBDxJ8XOeaWUuUhj8GeERh+vQcMCvFvSCYludKeoKglJYljt0ArJH0iJmdIel77C1DeqOZrUwoxy7DzH4t6fXAA9FscRfB6f1sdMqPCD630wvfofP48q/TrUja38x2RQ/614Anzey66NgU4DozO6HMPb5HcNL+LHOBnYK4aeN0Nx+KnKGPEcyJbwJEM6mfA4X8LE6N4TOSTiDpSMIqTJI9Zja1O+SpFpIuBd6Tt/unZnZ1d8iTTxbySXo1waeTz4lm9tfO3jfl2GcD5+ftvs/MPlbmupr7/rkicRynYty0cRynYlyROI5TMa5IHMepGFckjuNUjCsSx3Eq5v8DmkBspN/hlSsAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "
| Statistic | Value |
|---|---|
| Mean Average Error | 0.0331 |
| Root Mean Squared Error | 0.0523 |
| R2 | 0.904 |
"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "filename=f'{reference_dir}/inputs/small_molecules-1000.xyz'\n",
+ "new_frames = read(filename,\":400\")\n",
+ "new_property_values = np.array([cc.info[property_name] for cc in new_frames])\n",
+ "y_pred = krr.predict(new_frames)\n",
+ "\n",
+ "krr.plot(y_known=new_property_values, y=y_pred, property_name=property_name)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Training a KRR on Atomic Properties\n",
+ "Now let's consider a property that is atom-centric, like the chemical shift. The workflow of predicting this property does not change much--we can use the same hyperparaters and Kernel class--but we will note the changes."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Load and Read the Input File\n",
+ "We'll now look at the [CSD-500.xyz](https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-06972-x/MediaObjects/41467_2018_6972_MOESM4_ESM.txt) file, which contains the chemical shifts (CS) of 500 molecules. Because this is a property of the atoms, we'll want to use the \"Structure\" kernel given by libRascal.\n",
+ "\n",
+ "Because CSD-500 is a larger file, we don't automatically include it in the rascal repo, but it can be downloaded when needed:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "download_link = \"https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-06972-x/MediaObjects/41467_2018_6972_MOESM4_ESM.txt\"\n",
+ "if(not os.path.exists(f'{reference_dir}/inputs/CSD-500.xyz')):\n",
+ " print(\"Downloading CSD-500\")\n",
+ " import urllib.request\n",
+ " urllib.request.urlretrieve(download_link, f'{reference_dir}/inputs/CSD-500.xyz')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "input_file = f'{reference_dir}/inputs/CSD-500.xyz'\n",
+ "property_name = \"CS\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We then read the trajectory, store the desired property, and create our SOAP calculator.\n",
+ "\n",
+ "**Because the chemical shifts are per-atom properties, they have been stored in ase.arrays, as opposed to ase.info.**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "There are 100 frames and 13233 environments.\n"
+ ]
+ }
+ ],
+ "source": [
+ "frames = np.array(read(input_file,\":100\"))\n",
+ "for frame in frames:\n",
+ " frame.wrap()\n",
+ "number_of_frames = len(frames)\n",
+ "property_values = np.array([cc.arrays[property_name] for cc in frames])\n",
+ "print(f\"There are {number_of_frames} frames and {len(np.concatenate(property_values))} environments.\")\n",
+ "\n",
+ "representation = SOAP(**hyperparameters)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Preparing Data for KRR"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We still need to define the training percentage, $\\zeta$, $\\Lambda$, jitter, and the kernel_type, **however only the kernel_type will differ from the earlier example**. We can then split our dataset into training and testing sets and calculate the SOAP vectors."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "kernel_type = \"Atom\"\n",
+ "train_idx, test_idx = split_dataset(number_of_frames, training_percentage)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "features = representation.transform(frames[train_idx])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Constructing the KRR and Predicting the Properties\n",
+ "\n",
+ "The only change here is the properties array. **Because the properties are stored in lists for each molecule, we will need to pass the concatenated list of the properties to our KRR model.**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "krr = KRR(features, np.concatenate(property_values[train_idx]),\n",
+ " kernel_type=kernel_type, representation=representation,\n",
+ " Lambda=Lambda, jitter=jitter, weights=None,\n",
+ " zeta=zeta, **hyperparameters\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Training the KRR"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "krr.train()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Our weights have been calculated and have shape (10758,)\n"
+ ]
+ }
+ ],
+ "source": [
+ "y = krr.predict(frames[test_idx])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/markdown": [
+ "
"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEGCAYAAAC+UopnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfQElEQVR4nO3df7TcdX3n8eeLRLipZRNSUkgI+YGJUqSE2FtDj+wqCUogSnC3W2P2VKyezbr8aqq7mhSOVq0StdYY2tqmlS7sElNOKxCbIIZE1+NZCQSTIAjqNSRrCIG4mIs/uMGbvPeP73duvvdmZu7cufOd78zc1+Oce2bm8/3O3M9A8s7n5/ujiMDMLA+nFF0BM+tcDjBmlhsHGDPLjQOMmeXGAcbMcjO+6Ark6cwzz4xZs2YVXQ2zjvboo4/+JCKmlLvW0QFm1qxZ7Ny5s+hqmHU0SfsrXXMXycxy4wBjZrkpLMBI6pL0sKQ9kp6Q9NG0fLakHZJ6JP2TpFPT8tPS1z3p9VlF1d3MalNkC+YosDAi5gEXA4slXQJ8CvhcRMwBfgq8N73/vcBP0/LPpfeZWQsrLMBE4ufpy1ekPwEsBP45Lb8DuCZ9vjR9TXp9kSQ1qbpmVodCx2AkjZO0G3ge2Ar8CDgSEf3pLQeAc9Ln5wA/Bkiv9wK/UeYzV0jaKWnn4cOH8/4KZlZFoQEmIo5FxMXAdOD1wPkN+Mz1EdEdEd1TppSdmjezJmmJWaSIOAJ8Hfg9YJKk0vqc6cAz6fNngHMB0usTgf/X5KqadaSVG3dx3urNrNy4q6GfW+Qs0hRJk9LnE4A3A0+SBJrfT2+7Frgvfb4pfU16fXs4mY1ZQ2zac5DjkTw2UpEtmKnA1yU9BjwCbI2IfwU+BLxfUg/JGMsX0/u/CPxGWv5+YFUBdTbrSFfPm8YpSh4bSZ3cCOju7g5vFTDLl6RHI6K73LWWGIMxs8bbsGM/l9y6jQ07Km4VOsmla7Yxa9VmLl2zrSF1cIAxa0O1DMqu297Dod4+btveU/PnHjjSN+hxtBxgzNpQLYOyNy2cw9SJXdy4cM5AWaVWTal80oRkAnf6pK6G1NMBxqwN1TIou3zBTL69ehHLF8wcKKvUqllz/1Mc6u3jeMC+NUu47rI5I+5eleMAY9aG1i6bz95bl7B22fwRva9cqyartPemnu5VOR2dcMrMBlu+YOagFk3JqivP57btPdy4cA4rN+7iud4+usarYiCqlQOMmQ0KPLfc+zgBvHwsygajkXAXycwGaeSiOy+0M7NR8UI7MyuEA4xZGxvJLui8dkxX4wBj1sZGsgs6rx3T1TjAmLWxkQzI5rVjuhoP8prZqHiQ12yMKmLcJcsBxqyDFTHukuUAY9bBihh3yfIYjFnOVm7cxaY9B7l63rQRb06s5tI12zhwpG8gtULp+bdWLWrY76iFx2DMCpRXNyWbHCr7vJEZ6UbLAcYsZ3l1U0otl+mTuk5KEHXgSF/D01/Wo7AukqRzgTuBs0iOjF0fEZ+XNBn4J2AWsA/4g4j4aXpM7OeBq4BfAu+OiO9U+x3uItlYU+o2DbVvzZLcfmerdpH6gQ9ExAXAJcD1ki4gOY5kW0TMBbZx4niSK4G56c8K4AvNr7JZ8cpNPZfKumdN5pQhJ7Y3Kv1lPQrLBxMRzwLPps9/JulJkvOnlwJvSm+7A/gGyVlJS4E708PWHpI0SdLU9HPMOlp2oLg0pnPv7mRMZ+2y+YPGeUr3NHpQuR4tkXBK0ixgPrADOCsTNA6RdKEgCT4/zrztQFo2KMBIWkHSwmHGjBm51dksb9lZooNpvtxSUCnZtOcga5fNPymo1BJYsp+f18xT4QFG0q8D/wKsjIgXk6GWRESEpBENEkXEemA9JGMwjayrWTNlZ4amT+o6aWwlO3Bca1Cp9Pl5KTTASHoFSXC5KyK+nBY/V+r6SJoKPJ+WPwOcm3n79LTMrCOVgkqpBZN1zcWj7/5kPz8vhQWYdFboi8CTEfGXmUulQ+7XpI/3ZcpvkLQRWAD0evzF2l21bkr29cqNuwZ1jxoxttKMBXlFziK9AfhDYKGk3enPVSSB5c2Sfghcnr4G2ALsBXqAvweuK6DOZg1Vazdl7bL5A0eKqOqdraXIWaRvUfm/1UmhNZ09uj7XSpk12Ui6KU/nuJYlL4UP8pqNZc3eN9RsDjBmLa4Z08l58V4ksxbXjOnkvLgFY9YE9bZCshsVi1zyXy8HGLMmqNYKKW0DmDYxWe+SXeKfvb/dukfgLpJZU2RTK2TNXrWZe3cn+4gOHOk7KW9Mpfe1C7dgzJqge9ZkDu45SPesyYPKy+1lOR5w/i1b+PDbXtuWrZYsBxizMho9c5Pd7bxz3wvDrn3p6w9u297D8gUzR/27i+QuklkZjZ65KWW1u+iciTV9dtd4cePCOQ353UVyC8asjEZtBBy6h2j3gd6a3rf4wqlt33oBBxizsurpFpXrVtWb6LuU56XduYtk1iDluj6njat9a+IpStIwFHmOUaO5BWPWANn8uKVu1cqNu3ipv/acZyPJRtcuHGDMGiDbFSqdTTQSeWb9L5K7SGajtGHHfk4dp7rztLRTfpeRcoAxG6V123vo6w/OntjFNRePfOxEHRxhHGDMRummhXPoGi+e7e07Ket/OaV4MmG8OmpAtxyPwdiY0sgVuqXPmjBe9NU4mLtvzRLOW72ZCDh6LNh7a2eOvZS4BWNjSiNX6JY+YyQzRdnxmk5uuZQ4wNiYkl2ZO9pD4UeyyrdrvPjk2y8cNF7TSdPRlRQaYCTdLul5SY9nyiZL2irph+njGWm5JK2T1CPpMUmvK67m1k427NjPJbduY8OO/YO6RbW0YrLvzap0yHwlpc2LNy2cw9SJXR2xz6gWRbdg/geweEjZKmBbRMwFtqWvAa4kOfh+LsnRsF9oUh2tza3b3sOh3j5u294DjCzHSva9l67ZxqxVm0ccXADGC35xtB+Ab69e1BH7jGpRaICJiG8CLwwpXkpy6D3p4zWZ8jsj8RAwKT350ayqoa2Gb61axL41S2oa5M2+d+hRrsPZt2bJwM+xgBf7+rn5nseHfV8nacVZpLMyJzYeAs5Kn58D/Dhz34G0zKc7WlXLF8ysu8WwfMFMHn76BW6593EmjBcv9QfTJ3XRPWsyB9Ip6a4Ks0il1bzTJ3UNJJZq5mHprXAaQdFdpKrSw9ZG9P9E0gpJOyXtPHz4cE41s041dMxlw479Ayktjx6LgZZPdmvAcFPU2RZPM1NftsJpBK3YgnlO0tSIeDbtAj2flj8DnJu5b3paNkhErAfWA3R3dzfzHwxrc9mxlVI2uXXpuA3AtIldnLd6M1fPm8bV86axac9BLjpnIs/97CiHevsq/ktYVAuiGYfbD6cVA8wmkkPv16SP92XKb5C0EVgA9Ga6Umajlv2X/qWX+wdtWJw0YTwHe08k5d5765JB08zVTgYoSivk8y16mvpLwLeB10g6IOm9JIHlzZJ+CFyevgbYAuwFeoC/B64roMrWhipNNQ8tz3ZjjrzUP+jeIy/1D6S9LLdAbu2y+UybmLQYpo2RNS61UDLM0Zm6u7tj586dRVfDCnbJrds41NvH1IldfHv1omHL4eR1LtluzuxVmwmSPUXZA+mzLZ5rLi6+BdMskh6NiO6y1xxgrNNt2LGf27b3cOPCOYNmkzbs2M/N9zxOMPw4SakLdPW8aYM2NJbGOcTg2YhTRMfvMyqpFmBaehbJrBGWL5hZdnHb8gUzB4LCcDMt2WNHSruhlXnf0H+mx8I+o1o4wNiYVhp3mTBenLd6Mys37mLlxl3MWrWZWas2D6TCLI2/TJvYhdLcuU+vWTLw/g5O6TIq7iJZW8t2XUY65pFdiFaaISonm87yvNWbOR7lu0DZI0rcRUq4BWNtLdt1qcXKjbsGWirZhWi1Zv8v3Vfu/rXL5nfcqQCj1YrrYMxqVlrwNtxf6KGzQvfuPjho6X+1MZhL12wbGAAu5X6plAOm004FGC0HGGtr1f5CZ7tP5QJIKbgc7K0+wFvkUvt25y6Sdaxauk8HjlQeeymn1AWqJ7n3WDSiACPpDKmTc6BbJ8muvG3Ufpy1y+aftE3AKqs4iyTpw8DdEfGUpNOArwLzgH5geUQ82Lxq1sezSJZVT6Io6NxD0Rql3lmkdwDfT59fmz5OAd4IfLJx1TPLX73BJbs+xkauWoB5OU40b64ANkbEsYh4Eg8OWwvJprKsdL2e4DJ9UhdHj0XN0+DZKXBLVAswRyVdKGkKcBnwtcy1X8u3Wma1K5dYaTT5c0u+tWpR1R3UQ410Tc5YUK0lshL4Z5Ju0eci4mkASVcBDtHWMsolVhptNrfSTMZI1rXUuiZnLPFWAWtrQ/PO1tNimdg1nsvO/82BZf7D7axuhVy3raSuQV5J708TQA0tf6+klY2soFm9si2VDTv219Vi6e3rH7SHqHvW5KpjKa2Q67ZdVBuD+U/AnWXK/yfwnnyqYzYy2Sx02fy59SiNtQw3ljKaBN7DDUh3mmpjMOMj4ldDCyPiZS+2syJluyglPz/aP6IVueVkdz9XG0sZTbdorLV+qgWYUySdFRHPZQslnVXpDWbNUO4v6dAcuiMxdJYozw2LQwekO308p1qA+QywWdIHgO+kZb+Tlv9F3hUzq2ToX9LRtAY++fYLm3qM69Ag0uktmooBJiLulHQY+BhwIUlWwCeAD0fE/U2qn9lJyv1Ln024XStB4WdEt8LZRXmquiI3DSQtFUwkLQY+D4wD/iEi1gzzFhsDhsvpMlSrJIXqxG5RVlula5A0Dvhr4ErgAuCdki4otlbW6kSyYTHbSvCK2+Zotz1Frwd6ImIvQHrK41Lge4XWynJRbQB06LVKrZfsTujSZ2QTUVm+2i3AnAP8OPP6AMkxsgMkrQBWAMyYMaN5NbOK6p0pqbTHKPu69LxcF6m0E3poQnCntWyeigFG0vurvTEi/rLx1Rm9iFgPrIdkq0DB1TGGnykpBY1JE8bTdep4Lpk9me1PPV/1s0oEXPzRBzjyUv9JAax0AsCmPQcdUApSrQVzevr4GuB3SQ6fB3gb8HCelariGeDczOvpaZm1sOFmSkpB48hL/fBS/8BK2qyLP/oAXaee/Mc1OLEG5sCRPi5ds23g8HlvPizesJsdJX0TWBIRP0tfnw5sjoh/14T6Da3LeOAHwCKSwPIISXa9J8rd782O7SHbgun71TH6+oPxAp0ifnWs/kaoM9E1R7XNjrWMwZwFvJx5/XJa1nQR0S/pBuABkmnq2ysFF2sfQ7s1AMeBvZ+4qubd0ddcPK1sy6eZOn1Vbj1qmaa+E3hY0p9J+jNgB3BHrrWqIiK2RMSrI+JVEfGJouph+Ria4Olbqxaxb82SQdn8961ZMmjaefqkroFk3EVm/e/0Vbn1qCkfjKTXAf82ffnNiGiLhFPuIlkzjdUWzGi7SJCkyHwxIv5R0hRJs0sZ7swsMZaCSq2G7SJJ+gjwIWB1WvQK4H/lWSmzkg079nPJrdvYsGN/0VWxOtTSgnk7MJ90R3VEHExnkswG/NYtW3ipP5gwXjz551c17HPXbe/hUG8fn7r/KdZt7+GmhXMK36BotatlkLd0fEkASHplvlWydjTcofC1KHfsx00L5zB1YhcBHOrt47Y0a121lo1bPa2jlgBzt6S/AyZJ+s/Ag8A/5FstazcTxmvQYz3KpapcvmAm3169iFVXns/UiV3cuHAOcKJlc1uZNJnVrllzDdtFioi/kPRm4EWSVb0fjoitudfM2kq5btGGHftH1K2ptvJ2+YKZgz7jUG8yFfxs78lTwjctnMNt23sGgpEVZ9gAI+lTEfEhYGuZMrOysgvkbtvew998vafiFG52d3M2Ly6MPEjBycHIilNLF+nNZcqubHRFrLNkF5vduHBO1UVo1bL4V+ruaMijtaZqu6n/K3Ad8CpJj2UunQ78n7wrZu0tu8Fx+YKZg1owQ1XrGp19+mkc6u3jrNNPG1T+tPcZtYWKK3klTQTOAG4FVmUu/SwiXmhC3UbNK3mba7juzHArXbNdpVJ6hVLKhVPESd0naw11newYEb0RsY8k/+0LEbE/IvYD/ZIWVHqfjU0bduznlnsfrzp7U66blJ2aLtdVKu1Nmjaxq+ppi9aaahmD+QLw88zrn6dlZgPWbe8ZaGlUmr0pdyJiNqgM3egIDGxiPNjb17Q8uqXTF2ev2jymTmHMQy0reRWZflREHE/zspgNyE4NV5rBKdctKo2/TJvYdVL3qNx9zUgeVWphxZDXNnK1JJz6MvANTrRargMui4hr8q3a6HkMpn200lhLaaxIJEFmrO2OHqnR7qZ+H7AOuIXkv/c20qTaZo3SSuktHUwap6Z8MO3KLZjWMVZzpYwFdbVgJH0wIj4t6TZOdEcHRMRNDayjdThnexubqnWRnkwf3QQYo2ppdZRbu1KurNzJAm7VdL6KASYivpI+FpZ/14pVS6ujNM187+6D3Lv7INMndQ2aUi4FmHIBxK2azldxHYykr0jaVOlnNL9U0n+U9ISk45K6h1xbLalH0vclXZEpX5yW9UhadfKnWqOVW7cyVGntSsmBI31l17PAyfleavl8a2/Vtgq8MX3674GzOZEm853AcxHxJ3X/Uum3SE6m+Dvgv0XEzrT8AuBLJGdQTyPJPfPq9G0/INl4eYDkPKR3RkTVM6k9yNs8tXR3Wmkq2hqnrkHeiPjf6Zs/O+TNX5E0qr+1EfFk+tlDLy0FNkbEUeBpST0kwQZ86H1Lq2UMpZWmoq05alkH80pJ52X+cs8G8kqbeQ7wUOb1gbQMhjn0vkTSCtJ1OjNmzMihilYvHzo/9tQSYP4E+IakvSTpN2YC/2W4N0l6kKRrNdTNEXHfiGo5AhGxHlgPSRcpr99jZsOrJWXmVyXNBc5Pi55KuzDDve/yOupT7XB7H3rf4cpNb1t7q+VcpF8D/jtwQ0TsAWZIemtO9dkELJN0WtoVmws8TDKoO1fSbEmnAsvSe62DVMtsV+7EAWt9taRr+EeSA+9/L339DPDno/mlkt4u6UD6mZslPQCQHmR/N8ng7VeB6yPiWET0A6VD758E7vah9+2rlA5haBqEStPbUD34WOuqZQzmVRHxDknvBIiIX6rM9M9IRMQ9wD0Vrn0COOlQ+4jYAmwZze+11lBpgV21QWDPQLWnWgLMy5ImcOLgtVcBw47BmFVSbtvAcDwD1Z5qCTAfIemunCvpLuANwLvzrJR1Nu87GjuqBpi0K/QUyWreS0imqf84In7ShLqZWZurGmAiIiRtiYjfBjY3qU7WYrzr2epVyyzSdyT9bu41sZblXc9Wr1oCzALgIUk/kvSYpO8OOYjNOpx3PVu9ahnkvWL4W6zTuFtkjVAtZWYXScLvOcB3gS+mC96sw5QLJu4WWSNU6yLdAXSTBJcrgc82pUbWdOWCibtF1gjVukgXpLNHSPoiyZ4g60DlFr65W2SNUC3A/Kr0JCL6R7k7wFqYg4nlpVqAmSfpxfS5gAnpa5Eskfk3udfOzNpatZSZ45pZETPrPD7EfgzzVLTlrZaFdtahPBVteXOAGcM8FW15cxdpDHO3yPLmFoyZ5cYBxsxy4wBjZrlxgDGz3BQSYCR9RtJTaX6ZeyRNylxbLalH0vclXZEpX5yW9UhaVUS9zWxkimrBbAUujIiLgB8AqwEkXUByqNprgcXA30gaJ2kc8Ncku7ovAN6Z3mtmLayQABMRX8vklnmI5ChYgKXAxog4GhFPAz3A69OfnojYGxEvAxvTe82shbXCGMx7gPvT5+cAP85cO5CWVSo/iaQVknZK2nn48OEcqmtmtcptoZ2kB4Gzy1y6OSLuS++5GegH7mrU742I9cB6gO7u7mjU55rZyOUWYCLi8mrXJb0beCuwKCJKgeAZ4NzMbdPTMqqUm1mLKmoWaTHwQeDqiPhl5tImYJmk0yTNBuaSZNJ7BJgrabakU0kGgjc1u95mNjJF7UX6K+A0YGuaKe+hiHhfRDwh6W7geyRdp+sj4hiApBuAB4BxwO0R8UQxVTezWulE76TzdHd3x86dO4uuhllHk/RoRHSXu9YKs0hm1qEcYMwsNw4wZpYbBxgzy40DjJnlxgHGzHLjAGNmuXGAMbPcOMCYWW4cYMwsNz4XqcW0ynGuKzfuYtOeg1w9bxprl80vrB7W3tyCydmla7Yxa9VmLl2zrab7W+U41017DnI8kkezejnA5GykAWO441xXbtzFeas3s3LjrsZUsIKr503jFCWPZvVyFyln0yd1DXR5ajFctyjbssiz67J22Xx3jWzUHGBy1uhxlKvnTRsYGzFrdc4HY2aj4nwwZlYIBxgzy40DjJnlxgFmFJo1ZWzWroo6tuTj6cH3uyV9TdK0tFyS1qUH3D8m6XWZ91wr6Yfpz7VF1HsoL0Yzq66oFsxnIuKiiLgY+Ffgw2n5lSRnIc0FVgBfAJA0GfgIsIDknOqPSDqj6bUewovRzKorZB1MRLyYeflKoDRXvhS4Mz3p8SFJkyRNBd4EbI2IFwAkbQUWA19qXq1PNtrFaK2y78gsL4UttJP0CeBdQC9wWVpc6ZD7SuXlPncFSeuHGTNmNLbSDdYq+47M8pJbF0nSg5IeL/OzFCAibo6Ic0kOvr+hUb83ItZHRHdEdE+ZMqVRH1uzkWxuHG7fkVm7y60FExGX13jrXcAWkjGWZyh/yP0zJN2kbPk3Rl3JHIykVeJukXW6omaR5mZeLgWeSp9vAt6VziZdAvRGxLMkZ1K/RdIZ6eDuW9KyluNWidkJRY3BrJH0GuA4sB94X1q+BbgK6AF+CfwRQES8IOnjwCPpfR8rDfi2GrdKzE4oahbpP1QoD+D6CtduB27Ps17twJnmrJ14JW+b8eI+aycOMG3Gi/usnTgfTAtwt8famfPBtDh3e6xTOcC0AHd7rFO5i2Rmo+IukpkVYsyfKlBtR7N3O5uNzphvwVTbO+TdzmajM+YDTLW9Q95XZDY6HuQ1s1HxIK+ZFcIBxsxy4wBjZrlxgDGz3DjAmFluHGDMLDcOMGaWmzG/VaAcbxEwawy3YMrwFgGzxig0wEj6gKSQdGb6WpLWSeqR9Jik12XuvTY9+P6Hkq7Ns17eImDWGEUeHXsuyflG/zdTfCXJwfdzSQ66/wKwQNJkkoPZuknOsX5U0qaI+GkedXO3yKwximzBfA74ICcOvofkELY7I/EQMEnSVOAKYGtEvJAGla3A4qbX2MxGpKiTHZcCz0TEniGXKh1yX6m83GevkLRT0s7Dhw83sNZmNlK5dZEkPQicXebSzcCfknSPGi4i1gPrIdlNncfvMLPa5BZgIuLycuWSfhuYDeyRBMlB9t+R9HqSQ+7Pzdw+PS17BnjTkPJvNLzSZtZQTe8iRcR3I+I3I2JWRMwi6e68LiIOAZuAd6WzSZcAvRHxLMlB92+RdIakM0haPw80u+5mNjKtttBuC3AV0AP8EvgjgIh4QdLHgUfS+z4WES8UU0Uzq1XhASZtxZSeB3B9hftuB25vUrXMrAE6OmWmpMPA/gqXzwR+0sTqFGWsfE/wdy3KzIiYUu5CRweYaiTtrJRHtJOMle8J/q6tyHuRzCw3DjBmlpuxHGDWF12BJhkr3xP8XVvOmB2DMbP8jeUWjJnlzAHGzHIzZgJMqya3aiRJn5H0VPp97pE0KXNtdfpdvy/pikz54rSsR9KqYmo+ep3yPSDJlSTp65K+J+kJSX+clk+WtDX9c7k13TZT9c9y4SKi439INlA+QLLo7sy07CrgfkDAJcCOtHwysDd9PCN9fkbR36HG7/kWYHz6/FPAp9LnFwB7gNNINpr+CBiX/vwIOA84Nb3ngqK/Rx3fuyO+R+b7TCXZnwdwOvCD9P/hp4FVafmqzP/fsn+WW+FnrLRgxkRyq4j4WkT0py8fItl1Dsl33RgRRyPiaZK9Xq9Pf3oiYm9EvAxsTO9tN53yPQCIiGcj4jvp858BT5LkP1oK3JHedgdwTfq80p/lwnV8gMkzuVWLew/Jv2rQ+d+1U77HSSTNAuYDO4CzIskuAHAIOCt93rLfv/DNjo1QVHKrIlT7rhFxX3rPzUA/cFcz62aNJenXgX8BVkbEi2n+JCDZGCyp5deYdESAiTGU3KrSdy2R9G7grcCiSDvoVP6uVClvJ9W+X1uS9AqS4HJXRHw5LX5O0tSIeDbtAj2flrfu9y96EKiZP8A+TgzyLmHwwNjDaflk4GmSAd4z0ueTi657jd9vMfA9YMqQ8tcyeJB3L8nA6Pj0+WxODI6+tujvUcf37ojvkfk+Au4E1g4p/wyDB3k/nT4v+2e5FX46ogVTp05MbvVXJEFka9pieygi3hcRT0i6myT49APXR8QxAEk3kMywjQNuj4gniql6/SKivxO+R8YbgD8Evitpd1r2p8Aa4G5J7yWZEf2D9FrZP8utwFsFzCw3HT+LZGbFcYAxs9w4wJhZbhxgzCw3DjBmlhsHGCuMpLMlbZT0I0mPStoi6dXpzuDHJX1X0iOSZhddV6vPWF4HYwVSslDnHuCOiFiWls0D3gFMAy6KiOOSpgO/KK6mNhpuwVhRLgN+FRF/WyqIZEPqL4BnI+J4WnYgkl3t1oYcYKwoFwKPlim/G3ibpN2SPitpfpPrZQ3kAGMtJSIOAK8BVgPHgW2SFhVbK6uXx2CsKE8Av1/uQkQcJdm8d7+k50gSK21rYt2sQdyCsaJsB06TtKJUIOkiSW+UNC19fQpwEZXPF7cW582OVpg0kKwFfgfoI0mn8VXg3SS7wgEeBq6LiL4Cqmij5ABjZrlxF8nMcuMAY2a5cYAxs9w4wJhZbhxgzCw3DjBmlhsHGDPLzf8Hw9IMmRTsH24AAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "
| Statistic | Value |
|---|---|
| Mean Average Error | 11.3629 |
| Root Mean Squared Error | 21.0725 |
| R2 | 0.9145 |
"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "krr.plot(y_known=np.concatenate(property_values[test_idx]),\n",
+ " y=y, property_name=property_name)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Raw Cell Format",
+ "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.8"
+ },
+ "toc": {
+ "base_numbering": 1,
+ "nav_menu": {
+ "height": "12px",
+ "width": "252px"
+ },
+ "number_sections": true,
+ "sideBar": true,
+ "skip_h1_title": false,
+ "title_cell": "Table of Contents",
+ "title_sidebar": "Contents",
+ "toc_cell": false,
+ "toc_position": {},
+ "toc_section_display": "block",
+ "toc_window_display": true
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/source/tutorials/TutorialEn.rst b/docs/source/tutorials/TutorialEn.rst
new file mode 100644
index 000000000..b00852ae7
--- /dev/null
+++ b/docs/source/tutorials/TutorialEn.rst
@@ -0,0 +1,342 @@
+.. _TutorialEn:
+
+.. role:: raw-html(raw)
+ :format: html
+
+Energy fitting
+==============
+
+Approach used in Rascal can be used for calculating any of the rotationally
+invariant properties (such as formation energy), as well as rotationally
+covariant properties (such as dipole moment). This example is about fitting the
+formation energies of some small molecules. The dataset is located in the
+`librascal/examples/data/ <../../../../examples/data>`_ folder, the python
+notebook doing the calculation is `librascal/examples/SOAP_example.ipynb
+(outbound)
+`_.
+
+
+Let's take a look at the code and describe what it do step-by-step.
+
+Imports
+*******
+
+Just import of the necessary modules, including Rascal itself.
+
+.. code-block:: python
+
+ from matplotlib import pylab as plt
+ import os, sys
+ sys.path.insert(0,"../build/")
+ import time
+ import rascal
+ import json
+ import ase
+ from ase.io import read, writex
+ from ase.build import make_supercell
+ from ase.visualize import view
+ import numpy as np
+ from rascal.representations import SOAP
+
+Dataset
+*******
+
+Here ASE is used to read the file with information about molecules, then the
+size of the feature matrix shown. The number of columns is the number of atoms
+in the set of molecules, the number of rows is depend on the max_radial and
+max_angular parameters. If the number of columns is big, sparsification needs to
+be done to reduce it.
+
+.. code-block:: python
+
+ frames = read('./data/small_molecules-1000.xyz',':')
+ hypers = dict(soap_type="PowerSpectrum",
+ interaction_cutoff=3.5,
+ max_radial=6,
+ max_angular=6,
+ gaussian_sigma_constant=0.4,
+ gaussian_sigma_type="Constant",
+ cutoff_smooth_width=0.5,
+ )
+ soap = SOAP(**hypers)
+ %time representation = soap.transform(frames)
+ X = representation.get_feature_matrix().T
+ X.shape
+
+
+Functions
+*********
+
+Here we define the functions, which are needed for the further computations.
+Let's describe some of them. compute_representation - creates a feature manager
+class using information about structure compute_kernel - compute the matrix of
+global similarity using global cosine kernel with a power of zeta to itself
+KRR.predict - predicting the energy of the molecules in the set, using trained
+model train_krr_model - train the model basing on the given dataset and using
+global cosine kernel
+
+.. code-block:: python
+
+ # Load the small molecules
+ frames = read('./data/small_molecules-1000.xyz',':600')
+
+ def compute_representation(representation,frames):
+ expansions = soap.transform(frames)
+ return expansions
+
+ def compute_kernel(zeta, rep1, rep2=None):
+ if rep2 is None:
+ kernel = rep1.cosine_kernel_global(zeta)
+ else:
+ kernel = rep1.cosine_kernel_global(rep2,zeta)
+ return kernel
+
+ def extract_energy(frames):
+ prop = [[]]*len(frames)
+ for ii,cc in enumerate(frames):
+ prop[ii] = cc.info['dft_formation_energy_per_atom_in_eV']
+ y = np.array(prop)
+ return y
+
+ def split_dataset(frames, test_fraction, seed=10):
+ N = len(frames)
+ ids = np.arange(N)
+ np.random.seed(seed)
+ np.random.shuffle(ids)
+ Ntrain = int(N*test_fraction)
+ train = ids[:Ntrain]
+ test = ids[Ntrain:]
+ targets = extract_energy(frames)
+ return [frames[ii] for ii in train],targets[train],[frames[ii] for ii in test],targets[test]
+
+ def get_mae(ypred,y):
+ return np.mean(np.abs(ypred-y))
+ def get_rmse(ypred,y):
+ return np.sqrt(np.mean((ypred-y)**2))
+ def get_sup(ypred,y):
+ return np.amax(np.abs((ypred-y)))
+ def get_r2(y_pred,y_true):
+ weight = 1
+ sample_weight = None
+ numerator = (weight * (y_true - y_pred) ** 2).sum(axis=0,dtype=np.float64)
+ denominator = (weight * (y_true - np.average(y_true, axis=0, weights=sample_weight)) ** 2).sum(axis=0,dtype=np.float64)
+ output_scores = 1 - (numerator / denominator)
+ return np.mean(output_scores)
+
+ score_func = dict(
+ MAE=get_mae,
+ RMSE=get_rmse,
+ SUP=get_sup,
+ R2=get_r2,
+ )
+
+ def get_score(ypred,y):
+ scores = {}
+ for k,func in score_func.items():
+ scores[k] = func(ypred,y)
+ return scores
+
+ class KRR(object):
+ def __init__(self,zeta,weights,representation,X):
+ self.weights = weights
+ self.representation = representation
+ self.zeta = zeta
+ self.X = X
+
+ def predict(self,frames):
+ features = compute_representation(self.representation,frames)
+ kernel = compute_kernel(self.zeta , self.X, features)
+ return np.dot(self.weights, kernel)
+
+ def train_krr_model(zeta,Lambda,representation,frames,y,jitter=1e-8):
+ features = compute_representation(representation,frames)
+ kernel = compute_kernel(zeta,features)
+ # adjust the kernel so that it is properly scaled
+ delta = np.std(y) / np.mean(kernel.diagonal())
+ kernel[np.diag_indices_from(kernel)] += Lambda**2 / delta **2 + jitter
+ # train the krr model
+ weights = np.linalg.solve(kernel,y)
+ model = KRR(zeta, weights,representation, features)
+ return model,kernel
+
+
+Full spectrum
+*************
+
+Here the full (with radial and angular parts) energies are computed. Let's
+describe the parameters of the soap descriptor, defined in the ``hypers``
+dictionary.
+
+- **interaction_cutoff**: Maximum pairwise distance for atoms to be considered
+ in expansion
+- **max_radial**: number of radial basis functions
+- **max_angular**: highest angular momentum number in the expansion
+- **gaussian_sigma_constant**: specifies the atomic Gaussian widths, in the
+ case where they're fixed.
+- **gaussian_sigma_type**: how the Gaussian atom sigmas (smearing widths) are
+ allowed to vary -- fixed (``Constant``), by species (``PerSpecies``), or by
+ distance from the central atom (``Radial``)
+- **cutoff_smooth_width**: the distance over which the the interaction is
+ smoothed to zero
+
+.. code-block:: python
+
+ hypers = dict(soap_type="PowerSpectrum",
+ interaction_cutoff=3.5,
+ max_radial=6,
+ max_angular=6,
+ gaussian_sigma_constant=0.4,
+ gaussian_sigma_type="Constant",
+ cutoff_smooth_width=0.5,
+ )
+ soap = SOAP(**hypers)
+
+ frames_train, y_train, frames_test, y_test = split_dataset(frames,0.8)
+
+ zeta = 2
+ Lambda = 5e-3
+ krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)
+
+ y_pred = krr.predict(frames_test)
+ get_score(y_pred, y_test)
+
+ plt.scatter(y_pred, y_test, s=3)
+ plt.axis('scaled')
+ plt.xlabel('DFT energy / (eV/atom)')
+ plt.ylabel('Predicted energy / (eV/atom)')
+
+The result of this block is:
+
+.. image:: ../resources/images/R1.png
+
+The result is quite good. One can try to change the train dataset to see how it
+affects the precision of the result.
+
+Radial spectrum
+***************
+
+Here we compute the energy, supposing the angular component to be zero.
+
+.. code-block:: python
+
+ hypers = dict(soap_type="RadialSpectrum",
+ interaction_cutoff=3.5,
+ max_radial=6,
+ max_angular=0,
+ gaussian_sigma_constant=0.4,
+ gaussian_sigma_type="Constant",
+ cutoff_smooth_width=0.5,
+ )
+ soap = SOAP(**hypers)
+
+ frames_train, y_train, frames_test, y_test = split_dataset(frames,0.8)
+
+ zeta = 2
+ Lambda = 5e-4
+ krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)
+
+ y_pred = krr.predict(frames_test)
+ get_score(y_pred, y_test)
+
+ plt.scatter(y_pred, y_test, s=3)
+ plt.axis('scaled')
+ plt.xlabel('DFT energy / (eV/atom)')
+ plt.ylabel('Predicted energy / (eV/atom)')
+
+Comparison of full and radial spectrum:
+
+.. image:: ../resources/images/Comps.png
+
+It can be seen that the two spectres are quite similar, but the radial spectrum is much more simple to compute (as feature matrix is much smaller and the set of spherical harmonics doesn't have to be computed). It is quite an inteseting fact, but, unfortunately, this feature is probably not generalizable and should be just the feature of this particular dataset.
+
+Map of the dataset
+*******************
+Here we use sklearn to do `kernel principal component analysis (outbound) `_.
+
+.. code-block:: python
+
+ def compute_representation(representation,frames):
+ expansions = soap.transform(frames)
+ return expansions
+
+ def compute_kernel(zeta, rep1, rep2=None):
+ if rep2 is None:
+ kernel = rep1.cosine_kernel_global(zeta)
+ else:
+ kernel = rep1.cosine_kernel_global(rep2,zeta)
+ return kernel
+
+ def link_ngl_wdgt_to_ax_pos(ax, pos, ngl_widget):
+ from matplotlib.widgets import AxesWidget
+ from scipy.spatial import cKDTree
+ r"""
+ Initial idea for this function comes from @arose, the rest is @gph82 and @clonker
+ """
+
+ kdtree = cKDTree(pos)
+ #assert ngl_widget.trajectory_0.n_frames == pos.shape[0]
+ x, y = pos.T
+
+ lineh = ax.axhline(ax.get_ybound()[0], c="black", ls='--')
+ linev = ax.axvline(ax.get_xbound()[0], c="black", ls='--')
+ dot, = ax.plot(pos[0,0],pos[0,1], 'o', c='red', ms=7)
+
+ ngl_widget.isClick = False
+
+ def onclick(event):
+ linev.set_xdata((event.xdata, event.xdata))
+ lineh.set_ydata((event.ydata, event.ydata))
+ data = [event.xdata, event.ydata]
+ _, index = kdtree.query(x=data, k=1)
+ dot.set_xdata((x[index]))
+ dot.set_ydata((y[index]))
+ ngl_widget.isClick = True
+ ngl_widget.frame = index
+
+ def my_observer(change):
+ r"""Here comes the code that you want to execute
+ """
+ ngl_widget.isClick = False
+ _idx = change["new"]
+ try:
+ dot.set_xdata((x[_idx]))
+ dot.set_ydata((y[_idx]))
+ except IndexError as e:
+ dot.set_xdata((x[0]))
+ dot.set_ydata((y[0]))
+ print("caught index error with index %s (new=%s, old=%s)" % (_idx, change["new"], change["old"]))
+
+ # Connect axes to widget
+ axes_widget = AxesWidget(ax)
+ axes_widget.connect_event('button_release_event', onclick)
+
+ # Connect widget to axes
+ ngl_widget.observe(my_observer, "frame", "change")
+
+ # Load the small molecules
+ frames = read('./data/small_molecules-1000.xyz',':600')
+ hypers = dict(soap_type="PowerSpectrum",
+ interaction_cutoff=3.5,
+ max_radial=6,
+ max_angular=6,
+ gaussian_sigma_constant=0.4,
+ gaussian_sigma_type="Constant",
+ cutoff_smooth_width=0.5,
+ )
+ soap = SOAP(**hypers)
+ zeta = 2
+ features = compute_representation(soap, frames)
+ kernel = compute_kernel(zeta,features)
+ from sklearn.decomposition import KernelPCA
+ kpca = KernelPCA(n_components=2,kernel='precomputed')
+ kpca.fit(kernel)
+ X = kpca.transform(kernel)
+ plt.scatter(X[:,0],X[:,1],s=3)
+
+The result of this block is:
+
+.. image:: ../resources/images/PCAs.png
+
+It shows how the structures is located in the abstract 2D map, where similar
+structures are located near to each other, and the very different ones far from
+each other.
diff --git a/docs/source/tutorials/Tutorial_Template.ipynb b/docs/source/tutorials/Tutorial_Template.ipynb
new file mode 100644
index 000000000..816c06a5c
--- /dev/null
+++ b/docs/source/tutorials/Tutorial_Template.ipynb
@@ -0,0 +1,93 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "toc": true
+ },
+ "source": [
+ "; // dim = 1, 2 or 3
+ Eigen::MatrixXd MATRIX(NROWS,NCOLUMNS); //
+ using rascal::VecXi = Eigen::Matrix
+
+Creating a neighborhood manager
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+First, the neighborhood manager is created and initialized using the :cpp:class:`rascal::NeighborhoodManagerCell`.
+
+.. code-block:: cpp
+
+ using Manager_t = rascal::NeighbourhoodManagerCell;
+ Manager_t manager; // initialize the neighborhood manager
+
+While the neighborhood manager can be built for the first time using :cpp:class:`rascal::NeighborhoodManagerCell::build`,
+we recommend updating it using :class:`rascal::NeighborhoodManagerCell::update`. The following inputs should be provided in order:
+
+- the atomic positions matrix column wise as an `Eigen` matrix object
+- an array of the atomic numbers of the species
+- an array of centers ids
+- the cell vectors column wise as an `Eigen` matrix object
+- a boolean array of the periodic boundary conditions over the 3 axes
+- the maximum radial cutoff (in the units of the cell vectors and atomic positions)
+
+The folllowing code is an example of the inputs that need to be customized according to the user's use, as in to provide
+cell vectors, atomic positions and numbers and the ids of the centers.
+
+.. code-block:: cpp
+
+ Eigen::MatrixXd cell(3,3); // the cell vectors matrix
+
+ int number_of_atoms{100}; // the number of atoms in the structure
+ Eigen::MatrixXd positions(3,number_of_atoms); // atomic postions matrix
+ rascal::VecXi numbers(22); // atomic numbers array
+
+ std::vector center_ids{}; // a vector containing the atomic centers ids
+
+ std::array pbc{true, false, false}; // for periodic boundary conditions over the x axis only
+
+ double cutoff{3.5}; // a radial cutoff distance of 3.5 angstrom
+
+ manager.update(positions, numbers, center_ids, cell, pbc, cutoff);
+
+- .. note:: A neighborhood manager can be updated as much as necessary by providing all the upmentioned inputs.
+
+
+Commands of the neighborhood manager:
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In the current implementation of the neighborhood manager, the centers and their neighbors are :cpp:class:`rascal::Neighborhood Cell::ClusterRef` objects of order 1 and 2, respectively.
+In order to access objects of order 2, we need to iterate over the corresponding order 1 object (the corresponding center).
+This implementation has the advantage of not using explicit arrays to deal with center and neighbor properties.
+
+Several methods are implemented for the :cpp:class:`rascal::NeighborhoodCell::ClusterRef` objects (centers and their neighbors), such as retrieving
+the index of the object, its position and its atom type.
+
+- .. note:: The positions of order 2 are given with a certain offset relative to the position of the corresponding center of order 1. If one atom is included in the neighborhood of two different centers, it will have different positions depending on
+ the center being iterated over.
+
+- .. warning:: The relative positions of the neighbors of a center are calculated on the fly and are not stored. If needed they have to be stored manually.
+
+This is an example of code than can be used to print to the screen the positions, types and indices of the centers and their neighbors. It also calculates the relative shift in the positions of every center's neighbors
+ and its norm.
+
+.. code-block:: cpp
+
+ for (auto center : manager) {
+
+ int center_index{center.get_atom_index()}; // get the index of the center atomic
+ Eigen::MatrixXd position{center.get_position()}; // get the position vector of the center
+ auto center_type{center.get_atom_type()}; // get the atome type of the center
+ std::cout << "Neighbors properties : " << endl;
+
+ for (auto neigh : center){
+ int neigh_index{center.get_atom_index()}; // get the index of the neighbor
+ Eigen::MatrixXd position{neigh.get_position()}; // get the position vector of the neighbor
+ auto neigh_type{center.get_atom_type()}; // get the atome type of the neighbor
+ auto relative_shift{position - center.get_position()}; // compute the position offset
+
+ std::cout << "This is the position of atom " << neigh_index << " of a type " << neigh_type << endl;
+ std::cout << "The relative position is : " << neigh_position << endl;
+ std::cout << "The relative shift is : " << relative_shift << endl;
+ std::cout << "The norm of the shif is : " << relative_shift.norm() << " ang" << endl;
+ }
+ }
diff --git a/docs/source/tutorials/utilities/general_utils.py b/docs/source/tutorials/utilities/general_utils.py
new file mode 100644
index 000000000..25f7bd2d5
--- /dev/null
+++ b/docs/source/tutorials/utilities/general_utils.py
@@ -0,0 +1,123 @@
+import os
+import inspect
+import time
+
+from IPython.display import Markdown, display, clear_output, Image, HTML
+import ipywidgets as widgets
+import numpy as np
+import ase
+from ase.io import read
+from matplotlib import pyplot as plt
+from tqdm import tqdm_notebook as tqdm
+import pandas as pd
+import json
+
+# This is for backwards compatibility -- old versions of RASCAL will follow
+# the latter schema, newer versions the former.
+try:
+ from rascal.representations import SphericalInvariants as SOAP
+except:
+ from rascal.representations import SOAP
+
+# These are some imports/helpers for all of the tutorials
+style = json.load(open('./utilities/widget_style.json'))
+hyper_dict = json.load(open("./utilities/hyperparameter_presets.json"))
+hyper_vals = json.load(open("./utilities/hyperparameter_bounds.json"))
+
+
+def mask_body_order(hyperparams):
+ """ This is for backwards compatibility, as we will soon be moving to using
+ "body order" in place of "soap_type"
+
+ Author: R. Cersonsky
+
+ Keywords:
+ hyperparams - full or partial dictionary of hyperparameters
+
+ Returns:
+ dictionary of hyperparameters suitable for a SOAP object
+ """
+
+ soap_types = ["RadialSpectrum", "PowerSpectrum", "BiSpectrum"]
+ return {"soap_type": soap_types[hyperparams["body_order"]-1],
+ **{h: hyperparams[h] for h in hyperparams if h != 'body_order'}}
+
+
+def markdown_table_from_dict(d, headers):
+ """ Function to generate a markdown table from a dictionary
+
+ Author: R. Cersonsky
+
+ Keywords:
+ d - dictionary of type {parameter: value}
+ headers - optional headers for the table
+
+ Returns:
+ markdown string
+ """
+ return '
'
+
+
+def compute_kernel(calculator, features1, features2=None, kernel_type='Structure', **kwargs):
+ """ Function to calculate similarity kernel
+
+ Author: M. Veit or F. Musil, modified by R. Cersonsky
+
+ Keywords:
+ calculator - soap calculator
+ features1 - ASE-Atoms-like object to compute kernel of
+ features2 - ASE-Atoms-like object to compute kernel of.
+ Default is features1
+ kernel_type - "atom" or "structure"
+ kwargs - hyperparameters of soap calculator
+
+ Returns:
+ kernel function
+ """
+ from rascal.models import Kernel
+
+ kernel = Kernel(representation=calculator, name='Cosine', target_type=kernel_type, zeta=2, **kwargs)
+ return kernel(X=features1, Y=features2)
+
+def readme_button():
+ """ Function to generate and display a button to show the package README
+
+ Author: R. Cersonsky
+
+ """
+ def show_readme_for_button(b):
+ clear_output()
+ display(Markdown(str(''.join(open('../README.rst')))))
+
+ button = widgets.Button(description="Show README", style=style)
+ output = widgets.Output()
+ display(button, output)
+
+ button.on_click(show_readme_for_button)
+
+
+def _button_template_(options, description, disp_func, default=None):
+ """ Function to generate and display a ToggleButton instance for a set of
+ options
+
+ Author: R. Cersonsky
+
+ Keywords:
+ options - list of options
+ descriptions - string to describe button
+ disp_func - function object to be called when a button is clicked
+ default - value of the button, needs to be in list of options
+
+ Returns:
+ widgets.ToggleButtons Object
+ """
+ button = widgets.ToggleButtons(
+ options=options,
+ description=description,
+ button_style='',
+ style=style,
+ value=options[0] if default == None else default
+ )
+ display(button)
+ button.observe(disp_func, 'value')
+ return button
diff --git a/docs/source/tutorials/utilities/hyperparameter_bounds.json b/docs/source/tutorials/utilities/hyperparameter_bounds.json
new file mode 100644
index 000000000..3ca1f73fb
--- /dev/null
+++ b/docs/source/tutorials/utilities/hyperparameter_bounds.json
@@ -0,0 +1 @@
+{"body_order": {"options": [1, 2, 3], "name": "Body Order"}, "interaction_cutoff": {"options": [1.75, 5.0], "name": "$r_{cut}$"}, "max_radial": {"options": [0, 10], "name": "$n_{max}$"}, "max_angular": {"options": [0, 10], "name": "$l_{max}$"}, "gaussian_sigma_constant": {"options": [0.01, 1.0], "name": "$\\sigma$"}, "gaussian_sigma_type": {"fixed": "Constant"}, "cutoff_smooth_width": {"fixed": 0.5}}
\ No newline at end of file
diff --git a/docs/source/tutorials/utilities/hyperparameter_presets.json b/docs/source/tutorials/utilities/hyperparameter_presets.json
new file mode 100644
index 000000000..c3200ef85
--- /dev/null
+++ b/docs/source/tutorials/utilities/hyperparameter_presets.json
@@ -0,0 +1 @@
+{"Minimal Power Spectrum": {"body_order": 2, "interaction_cutoff": 3.5, "max_radial": 2, "max_angular": 1, "gaussian_sigma_constant": 0.5, "gaussian_sigma_type": "Constant", "cutoff_smooth_width": 0.0}, "Power Spectrum": {"body_order": 2, "interaction_cutoff": 3.5, "max_radial": 6, "max_angular": 6, "gaussian_sigma_constant": 0.4, "gaussian_sigma_type": "Constant", "cutoff_smooth_width": 0.5}, "Radial Spectrum": {"body_order": 1, "interaction_cutoff": 1.75, "max_radial": 6, "max_angular": 0, "gaussian_sigma_constant": 0.25, "gaussian_sigma_type": "Constant", "cutoff_smooth_width": 0.5}}
\ No newline at end of file
diff --git a/docs/source/tutorials/utilities/learning_utils.py b/docs/source/tutorials/utilities/learning_utils.py
new file mode 100644
index 000000000..4bb4c5985
--- /dev/null
+++ b/docs/source/tutorials/utilities/learning_utils.py
@@ -0,0 +1,644 @@
+import os
+import inspect
+import time
+
+from IPython.display import Markdown, display, clear_output, Image, HTML
+import ipywidgets as widgets
+import numpy as np
+import ase
+from ase.io import read
+from matplotlib import pyplot as plt
+from tqdm import tqdm_notebook as tqdm
+import pandas as pd
+import json
+
+# This is for backwards compatibility -- old versions of RASCAL will follow
+# the latter schema, newer versions the former.
+try:
+ from rascal.representations import SphericalInvariants as SOAP
+except:
+ from rascal.representations import SOAP
+
+from general_utils import (markdown_table_from_dict, \
+ readme_button, _button_template_,
+ compute_kernel, mask_body_order)
+
+# These are some imports/helpers for all of the tutorials
+style = json.load(open('./utilities/widget_style.json'))
+hyper_dict = json.load(open("./utilities/hyperparameter_presets.json"))
+hyper_vals = json.load(open("./utilities/hyperparameter_bounds.json"))
+
+# Properties preset for the learning tutorial
+known_properties = dict(CS="Atom",
+ dft_formation_energy_per_atom_in_eV="Structure")
+ignore = ['Natoms', 'numbers', 'Name',
+ 'positions', 'cutoff', 'nneightol', "NAME"]
+
+def split_dataset(frames, train_fraction, seed=10):
+ """ Function to split a dataset into testing and training portions.
+
+ Author: M. Veit or F. Musil
+
+ Keywords:
+ frames - ASE-Atoms like
+ train_fraction - float in [0.0,1.0], percentage of frames to train
+
+ Returns:
+ indices of training and testing frames
+ """
+
+ N = len(frames)
+ ids = np.arange(N)
+ np.random.seed(seed)
+ np.random.shuffle(ids)
+ Ntrain = int(N*train_fraction)
+ train = ids[:Ntrain]
+ test = ids[Ntrain:]
+ return np.array(train), np.array(test)
+
+
+def get_r2(y_pred, y_true):
+ """ Function to calculate R^2
+ Author: M. Veit or F. Musil
+
+ Keywords:
+ y_pred - list of predicted y values
+ y_true - list of known y values
+
+ Returns:
+ float
+ """
+ numerator = ((y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
+ denominator = ((y_true - np.average(y_true,
+ axis=0)) ** 2).sum(axis=0,
+ dtype=np.float64)
+ output_scores = 1 - numerator / denominator
+ return np.mean(output_scores)
+
+
+def get_score(y_pred, y_true):
+ """ Function to calculate accuracy statistics
+
+ Author: M. Veit or F. Musil
+
+ Keywords:
+ y_pred - list of predicted y values
+ y_true - list of known y values
+
+ Returns:
+ dictionary containing SUP (supremum), MAE (mean squared error),
+ RMSD (root mean squared displacement),
+ and R^2
+ """
+ return {
+ "SUP": [np.amax(np.abs((y_pred-y_true)))],
+ "MAE": [np.mean(np.abs(y_pred-y_true))],
+ "RMSD": [np.sqrt(np.mean((y_pred-y_true)**2))],
+ r"$R^2$": [get_r2(y_pred, y_true)]
+ }
+
+
+class KRR(object):
+ """ Class for Kernel Ridge Regression
+
+ Author: M. Veit or F. Musil
+
+ Keywords:
+ zeta - integer
+ weights - weights of the given by the kernel and the targets given by
+ the training set
+ representation - feature set given by soap vectors
+ kernel_type - "atom" or "Structure"
+ X - training dataset
+
+ Functions:
+ predict - given the KRR, predicts the properties for a set of ASE frames
+ """
+
+ def __init__(self, weights, features, kernel_type, calculator=None, **kwargs):
+ self.weights = weights
+ self.hypers = dict(**kwargs)
+ self.calculator = calculator if calculator!=None else SOAP(**mask_body_order(kwargs))
+ self.X = features
+ self.kernel_type = kernel_type
+
+ def predict(self, frames):
+ features = self.calculator.transform(frames)
+ kernel = compute_kernel(calculator=self.calculator, features1=self.X, \
+ features2=features, kernel_type=self.kernel_type, \
+ **self.hypers)
+ return np.dot(self.weights, kernel)
+
+
+def extract_property(frames, property='energy'):
+ """ Function to read properties from ASE frames object
+
+ Author: M. Veit or F. Musil
+
+ Keywords:
+ frames - ASE Atoms-like object
+ property - property to extract
+
+ Returns:
+ list of the property from the given frames
+
+ Raises error when property is not found in the ASE object
+ """
+ if(property in frames[0].info):
+ return np.array([cc.info[property] for cc in frames])
+ elif(property in frames[0].arrays):
+ return np.array([cc.arrays[property] for cc in frames])
+ else:
+ print(frames[0].info)
+ raise KeyError(
+ "{} is not a property in the given frames".format(property))
+
+
+class learning_tutorial(object):
+ """ Class to wrap the learning tutorial
+
+ Author: R. Cersonsky
+
+ Constructor Keywords:
+ input_file - ASE Atoms-like object containing molecules or
+ crystals with some property to learn
+ training_percentage - float, percentage of frames to include in the
+ training set
+ interactive - boolean, whether to load the tutorial interactive-
+ ly in the jupyter notebook
+ verbose - boolean, whether or not to narrate the actions of
+ the tutorial
+ hyperparameters - dictionary, hyperparameters for SOAP calculation
+ number_of_frames - int, how many frames to include in the dataset
+ property - string, which property to train the dataset on
+
+ """
+
+ def __init__(self, input_file='./data/learning/small_molecules-1000.xyz',
+ training_percentage=0.8, interactive=False, verbose=True,
+ hyperparameters=dict(**hyper_dict['Power Spectrum']),
+ number_of_frames=None, property=None):
+
+ # presets given by M. Veit and F. Musil
+ self.zeta = 2
+ self.Lambda = 5e-3
+
+ self.hyperparameters = hyperparameters
+ self.verbose = verbose
+
+ # The verbosity wrap serves to limit the amount of text output in the
+ # jupyter notebook
+ self.verbosity_wrap = lambda s: (None if not verbose
+ else display(Markdown(s)))
+
+ # Looks in the examples folder for files suitable for the learning
+ # tutorial
+ file_options = ['./data/learning/{}'.format(f) for f in
+ os.listdir('./data/learning/') if f.endswith('xyz')]
+
+ # If the supplied file is not in the examples folder, includes it in the
+ # list of suitable files
+ if(input_file != None):
+ file_options.insert(0, input_file)
+ if(input_file in file_options[1:]):
+ file_options.pop(file_options[1:].index(input_file)+1)
+
+ # Sets the file to the default if not supplied
+ self.input_file = file_options[0] if input_file == None else input_file
+
+ # Reads input file
+ self.frames = np.array(read(self.input_file, ":"))
+ self.train_idx, self.test_idx = None, None
+
+ # This next section is for the jupyter interface. Even if the user has
+ # specified zero interactivity, the sliders serve as containers for the
+ # hyperparameters to change when set
+ self.sliders = {val:
+ widgets.FloatSlider(
+ value=self.hyperparameters.get(val,
+ hyper_vals[val]['options'][0]),
+ min=hyper_vals[val]['options'][0],
+ max=hyper_vals[val]['options'][1],
+ description=hyper_vals[val]['name'],
+ continuous_update=True,
+ step=(hyper_vals[val]['options'][1] -
+ hyper_vals[val]['options'][0])/20.,
+ style=style,
+ )
+ if isinstance(hyper_vals[val]['options'][0], float) else
+ widgets.IntSlider(
+ value=self.hyperparameters.get(val,
+ hyper_vals[val]['options'][0]),
+ min=hyper_vals[val]['options'][0],
+ max=hyper_vals[val]['options'][1],
+ description=hyper_vals[val]['name'],
+ continuous_update=True,
+ step=1,
+ style=style,
+ )
+ if isinstance(hyper_vals[val]['options'][0], int)
+ and hyper_vals[val]['options'][0] != True
+ else widgets.Dropdown(
+ options=hyper_vals[val]['options'],
+ style=style,
+ value=self.hyperparameters.get(val,
+ hyper_vals[val]['options'][0]),
+ description=hyper_vals[val]['name'],
+ )
+ for val in hyper_vals
+ if 'fixed' not in hyper_vals[val]}
+ self.properties = {prop:
+ extract_property(self.frames, prop)
+ for prop in sorted([*list(self.frames[0].info.keys()),
+ *list(self.frames[0].arrays.keys())],
+ key=lambda p: -int(p in known_properties))
+ if prop not in ignore}
+
+ self.sliders['number_of_frames'] = widgets.IntSlider(
+ value=int(len(self.frames)*0.2)
+ if number_of_frames == None
+ else number_of_frames,
+ min=1,
+ max=len(self.frames),
+ description="Number of Frames",
+ step=1,
+ style=style)
+ self.sliders['property_to_ml'] = widgets.Dropdown(
+ value=list(
+ self.properties.keys())[0]
+ if property == None
+ else property,
+ options=list(
+ self.properties.keys()),
+ description="Property to ML",
+ style=style)
+ if(self.sliders['property_to_ml'].value not in known_properties):
+ self.sliders['kernel_type'] = widgets.Dropdown(
+ value='Atom',
+ options=['Atom', 'Structure'],
+ description="Kernel Type",
+ style=style)
+ else:
+ self.sliders['kernel_type'] = widgets.Dropdown(
+ value=known_properties[self.sliders['property_to_ml'].value],
+ options=[
+ known_properties[self.sliders['property_to_ml'].value]],
+ description="Kernel Type",
+ style=style)
+ self.sliders['training_percentage'] = widgets.FloatSlider(
+ value=training_percentage,
+ min=0,
+ max=1,
+ description="Training Percentage",
+ continuous_update=True,
+ step=0.05,
+ style=style,
+ )
+
+ # Show and enable the widgets if interactive==True
+ if(interactive):
+ self.input_button = widgets.Dropdown(options=[*file_options], # , "Other"],
+ style=style,\
+ value=self.input_file,
+ description="Input File: ",
+ )
+ self.input_button.observe(self.get_input, names='value')
+ display(self.input_button)
+
+ self.preset_button = _button_template_(list(hyper_dict.keys()),
+ "SOAP Presets: ", self.preset_func)
+
+ slider_order = ['property_to_ml', 'kernel_type',
+ 'number_of_frames', 'training_percentage',
+ *list(hyper_vals.keys())]
+ for s in slider_order:
+ if(s in self.sliders):
+ self.sliders[s].observe(lambda change:
+ self.change_func(change),
+ names='value')
+ display(self.sliders[s])
+
+ # Set up one KRR for each possible property for caching reasons
+ self.krr = {prop: None for prop in self.properties}
+ self.trained = {prop: False for prop in self.properties}
+
+ # These serve as placeholders for time estimates until any calculation
+ # is run
+ self.est_frames, self.est_times = [], []
+ self.pred_frames, self.pred_times = [], []
+
+ def reset_ML(self, inp_change=False):
+ """ Function to reset the properties, number of frames, and kernel type
+ if any values affecting the machine learning have changed
+
+ Author: R. Cersonsky
+ """
+ if(inp_change):
+ self.properties = {prop:
+ extract_property(self.frames, prop)
+ for prop in sorted([*list(self.frames[0].info.keys()),
+ *list(self.frames[0].arrays.keys())],
+ key=lambda p: -int(p in known_properties))
+ if prop not in ignore}
+
+ self.sliders['number_of_frames'].max = len(self.frames)
+ self.sliders['number_of_frames'].value = int(len(self.frames)*0.2)
+
+ self.sliders['property_to_ml'].options = list(
+ self.properties.keys())
+ self.sliders['property_to_ml'].value = list(
+ self.properties.keys())[0]
+
+ if(self.sliders['property_to_ml'].value in known_properties):
+ self.sliders['kernel_type'].options = [
+ known_properties[self.sliders['property_to_ml'].value]]
+ self.sliders['kernel_type'].value = self.sliders['kernel_type'].options[0]
+ else:
+ self.sliders['kernel_type'].options = ['Atom', 'Structure']
+ self.sliders['kernel_type'].value = 'Atom'
+
+ self.krr = {prop: None for prop in self.properties}
+ self.trained = {prop: False for prop in self.properties}
+
+ def change_func(self, change):
+ """ Function for button events
+
+ Author: R. Cersonsky
+ """
+ change['owner'].value = change['new']
+ for s in self.hyperparameters:
+ if(s in self.sliders):
+ if(self.hyperparameters[s] != self.sliders[s].value):
+ self.est_frames, self.est_times = [], []
+ self.hyperparameters[s] = self.sliders[s].value
+
+ if(change['owner'] == self.sliders['property_to_ml']):
+ if(self.sliders['property_to_ml'].value in known_properties):
+ self.sliders['kernel_type'].options = [
+ known_properties[self.sliders['property_to_ml'].value]]
+ self.sliders['kernel_type'].value = self.sliders['kernel_type'].options[0]
+ else:
+ self.sliders['kernel_type'].options = ['Atom', 'Structure']
+ self.sliders['kernel_type'].value = 'Atom'
+
+ self.reset_ML()
+
+ def preset_func(self, a):
+ """ Function to change hyperparameters based upon presets
+
+ Author: R. Cersonsky
+ """
+ for s in self.sliders:
+ if(s in self.hyperparameters):
+ self.hyperparameters[s] = hyper_dict[self.preset_button.value][s]
+ self.sliders[s].value = self.hyperparameters[s]
+ self.reset_ML()
+
+ def get_input(self, a):
+ """ Function to change input file
+
+ Author: R. Cersonsky
+ """
+ inp = self.input_button.value
+ self.input_file = inp
+ self.frames = np.array(read(self.input_file, ":"))
+ self.reset_ML(True)
+
+ def train_krr_model_func(self, frame_idx, jitter=1e-8, pretend=False):
+ """ Function to trains the model on the given SOAP representation and
+ properties. It is used for both the time estimation and the eventual
+ training. The flag pretend is set to True when it is being used for
+ estimation, false otherwise.
+
+ Author: R. Cersonsky
+
+ Keywords:
+ frame_idx - list of indices of the frames to train on
+ jitter - float, anticipated error to adjust kernel diagonals
+ pretend - boolean, whether to not show the results
+
+ """
+
+ if(pretend == False):
+ self.output_params()
+ verbosity_wrap = lambda s: self.verbosity_wrap(s)
+ else:
+ verbosity_wrap = lambda s: None
+
+ verbosity_wrap("
We will now train a model on {}.".format( + self.sliders['property_to_ml'].value)) + + representation = SOAP(**mask_body_order(self.hyperparameters)) + + if(known_properties[self.sliders['property_to_ml'].value] == 'Atom'): + props = np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[:,0] + else: + props = self.properties[self.sliders['property_to_ml'].value][frame_idx] + + training_properties=props + + verbosity_wrap("First, I am going to separate my dataset:") + verbosity_wrap("
Now we will compute the SOAP representation of \ + our training frames.") + + t=time.time() + features=representation.transform(self.frames[frame_idx]) + + verbosity_wrap('This took {} seconds/frame.'.format( + round((time.time()-t)/len(frame_idx), 8))) + + if(pretend == False): + self.estimate_time(N=max(0, 20-len(self.est_frames)), + p="Estimating time to compute kernel...") + est=int(np.poly1d(np.polyfit(self.est_frames, + self.est_times, + deg=2))(len(frame_idx)))+1 + else: + est=0 + + verbosity_wrap("
Next we find the kernel for our training model.\ +
(This step will take approximately {} minutes and \ + {} seconds.)".format(int(est/60), int(est % 60))) + + time.sleep(0.5) + t=time.time() + kernel = compute_kernel(calculator=representation, \ + features1=features, \ + kernel_type=self.sliders['kernel_type'].value, \ + **self.hyperparameters) + + self.est_frames.append(len(frame_idx)) + self.est_times.append(time.time()-t) + + delta=np.std(training_properties) / np.mean(kernel.diagonal()) + kernel[np.diag_indices_from( + kernel)] += self.Lambda**2 / delta ** 2 + jitter + + verbosity_wrap("
We will adjust the our kernel with the tolerance \ + matrix $\\Lambda = ({})I$.".format(\ + round(self.Lambda**2 / delta ** 2 + jitter, 8))) + verbosity_wrap("
Now we can take this kernel to compute the weights \ + of our KRR.") + + weights=np.linalg.solve(kernel, training_properties) + model=KRR(weights, features, \ + kernel_type=self.sliders['kernel_type'].value, + **self.hyperparameters + ) + if(pretend == False): + self.krr[self.sliders['property_to_ml'].value], k=model, kernel + self.trained[self.sliders['property_to_ml'].value]=True + + def train_krr_model(self, jitter= 1e-8): + """ Wrapper to the train_krr_model_func function for use in the tutorial + + Author: R. Cersonsky + + Keywords: + jitter - float, anticipated error to adjust kernel diagonals + + """ + + training_dict={"Training Set": [int(self.sliders['number_of_frames'].value*(self.sliders['training_percentage'].value)),\ + round(100*(self.sliders['training_percentage'].value))], + "Testing Set": [int(self.sliders['number_of_frames'].value*(1-self.sliders['training_percentage'].value)),\ + round(100*(1-self.sliders['training_percentage'].value))]} + headers=["Partition", "Number of Frames", "Percentage"] + self.verbosity_wrap(markdown_table_from_dict(training_dict, headers)) + + self.train_idx, self.test_idx=split_dataset(self.frames[: self.sliders['number_of_frames'].value], + self.sliders['training_percentage'].value) + self.train_krr_model_func(self.train_idx, jitter= jitter) + + def plot_prediction_func(self, y_known= None, frames = None, + frame_idx= [], pretend = False): + """ Function to predict properties from the given SOAP representation and + KRR model. It is used for both the time estimation and the tutorial + prediction. The flag pretend is set to True when it is being used for + estimation, false otherwise. + + Author: R. Cersonsky + + Keywords: + y_known - list, target properties to train on + frames - frames to predict properties of + frame_idx - list of indices of the frames to train on + pretend - boolean, whether to not show the results + """ + if(len(frame_idx) > 0): + frames=self.frames[frame_idx] + y_known = np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[: , 0] if self.sliders['kernel_type'].value == 'Atom' else self.properties[self.sliders['property_to_ml'].value][frame_idx] + if(pretend == False): + verbosity_wrap=lambda s: self.verbosity_wrap(s) + else: + verbosity_wrap = lambda s: None + if(self.trained[self.sliders['property_to_ml'].value] == False): + verbosity_wrap("Model has not yet been trained, training now...") + self.train_krr_model() + # np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[:,0] if self.sliders['kernel_type'].value=='Atom' else self.properties[self.sliders['property_to_ml'].value][frame_idx] + testing_properties=y_known + + if(pretend == False): + self.estimate_time(x=self.pred_frames, y=self.pred_times, f=self.plot_prediction_func, N=max( + 0, 20-len(self.pred_frames)), ref=len(frames), p="Estimating time to compute prediction...") + est=int(np.poly1d(np.polyfit(self.pred_frames, + self.pred_times, deg=3))(len(frames)))+1 + else: + est=0 + + t=time.time() + verbosity_wrap("Predicting the properties of our data set will take approximately {} minutes and {} seconds.".format( + int(est/60), int(est % 60))) + y_pred=self.krr[self.sliders['property_to_ml'].value].predict(frames) + self.pred_frames.append(len(frame_idx)) + self.pred_times.append(time.time()-t) + verbosity_wrap(markdown_table_from_dict( + get_score(y_pred, testing_properties), headers = ["Statistic", "Value"])) + + if(not pretend): + plt.scatter(y_pred, testing_properties, s=3) + plt.axis('scaled') + plt.xlabel(self.sliders['property_to_ml'].value) + plt.ylabel('Predicted '+self.sliders['property_to_ml'].value) + plt.gca().set_aspect('equal') + plt.show() + + def predict_test_set(self): + """ Wrapper for plot_prediction_func that is used within the tutorial + + Author: R. Cersonsky + """ + self.plot_prediction_func(frame_idx=self.test_idx) + + def predict_new_set(self, filename='./data/small_molecules-1000.xyz', + num_frames=''): + """ Wrapper predicts a new set of data beyond the training set, loaded + from a different file + + Author: R. Cersonsky + + Keywords: + filename - string corresponding to an ASE-like file + num_frames - number of frames to predict + """ + frames=np.array(read(filename, ":{}".format(num_frames))) + properties=extract_property( + frames, self.sliders['property_to_ml'].value) + self.plot_prediction_func(y_known=properties, frames=frames) + + def output_params(self): + """ Function to output the hyperparameters that are currently set + + Author: R. Cersonsky + """ + self.verbosity_wrap('Our input file is {}, of which we are using {} \ + frames.'.format(self.input_file, \ + self.sliders['number_of_frames'].value)) + + hdict={hyper_vals[k]['name']: [self.hyperparameters[k]] + for k in self.hyperparameters + if 'fixed' not in hyper_vals[k]} + headers=["Parameter", "Value"] + + self.verbosity_wrap("
Our hyperparameters are {}".format( + markdown_table_from_dict(hdict, headers))) + + def set(self, value_name, value): + """ Public Setter function for tutorial interface + + Author: R. Cersonsky + + Keywords: + value_name - string, must be in self.sliders + value - value to set + """ + assert value_name in self.sliders + self.sliders[value_name].value=value + + def estimate_time(self, x=None, y=None, f=None, N=20, ref=None, p=None): + """ Helper function for estimating time + + Author: R. Cersonsky + + Keywords: + x - list of numbers of frames computed + y - list, size=size(x), of times to computed frames in x + f - function to estimate the run time of + N - number of trials to run + ref - int, upper limit of trials to run. + Defaults to (number_of_frames in the ultimate computation) / 4 + p - print statement + """ + if(N <= 0): + return + if(x == None or y == None or f == None or ref == None): + x=self.est_frames + y=self.est_times + f=self.train_krr_model_func + ref=int(0.25*self.sliders['number_of_frames'].value * + self.sliders['training_percentage'].value) + + if(p != None and min(N, ref) > 0): + self.verbosity_wrap(p) + for nf in tqdm(np.random.randint(2, ref, size = min(N, ref))): + f(frame_idx=np.random.randint(len(self.frames), size = nf), pretend=True) diff --git a/docs/source/tutorials/utilities/tutorial_utils.py b/docs/source/tutorials/utilities/tutorial_utils.py new file mode 100644 index 000000000..e69de29bb diff --git a/docs/source/tutorials/utilities/vectors_utils.py b/docs/source/tutorials/utilities/vectors_utils.py new file mode 100644 index 000000000..eff2088ea --- /dev/null +++ b/docs/source/tutorials/utilities/vectors_utils.py @@ -0,0 +1,429 @@ +import os +import inspect +import time + +from IPython.display import Markdown, display, clear_output, Image, HTML +import ipywidgets as widgets +import numpy as np +import ase +from ase.io import read +from matplotlib import pyplot as plt +from tqdm import tqdm_notebook as tqdm +import pandas as pd +import json + +# This is for backwards compatibility -- old versions of RASCAL will follow +# the latter schema, newer versions the former. +try: + from rascal.representations import SphericalInvariants as SOAP +except: + from rascal.representations import SOAP + +from general_utils import (markdown_table_from_dict, \ + readme_button, _button_template_, + compute_kernel, mask_body_order) + +# These are some imports/helpers for all of the tutorials +style = json.load(open('./utilities/widget_style.json')) +hyper_dict = json.load(open("./utilities/hyperparameter_presets.json")) +hyper_vals = json.load(open("./utilities/hyperparameter_bounds.json")) + + +def format_positions(positions): + """ Function to format positions for best viewing angle. Does a simple PCA + to find the best direction to view the molecule and rotates position + along these axes + + Author: R. Cersonsky + + Keywords: + positions - list of 3-vectors + + Returns: + positions - list of 3-vectors + """ + from sklearn.decomposition import PCA + pca = PCA(n_components=3) + return pca.fit_transform(positions) + + +class SOAP_tutorial(object): + """ Class to wrap the SOAP vectors tutorial + + Author: R. Cersonsky + + Constructor Keywords: + molecule_file - filename of ASE Atoms-like object + containing molecule + interactive - boolean, whether to load the tutorial interactive- + ly in the jupyter notebook + verbose - boolean, whether or not to narrate the actions of + the tutorial + hyperparameters - dictionary, hyperparameters for SOAP calculation + or string corresponding to one of the SOAP presets + + Functions: + predict - given the KRR, predicts the properties for a set of ASE frames + """ + + def __init__(self, + molecule_file=None, + interactive=False, + verbose=True, + hyperparameters=dict(**hyper_dict['Power Spectrum']), + ): + + # Check if hyperparameters are of type string or dictionary + if(isinstance(hyperparameters, str)): + self.hyperparameters = hyper_dict[hyperparameters] + else: + self.hyperparameters = {h: hyperparameters[h] + if h in hyperparameters + else hyper_dict['Power Spectrum'][h] + for h in hyper_dict['Power Spectrum']} + self.verbose = verbose + self.verbosity_wrap = lambda s: ( + None if not verbose else display(Markdown(s))) + + # Look for and load available molecules + file_options = ['./data/molecules/{}'.format(f) + for f in os.listdir('./data/molecules/') + if f.endswith('xyz')] + + if(molecule_file != None and molecule_file not in file_options): + file_options.append(molecule_file) + + self.frames, self.positions, self.symbols, self.inds, self.labels = {}, {}, {}, {}, {} + self.import_molecules(file_options) + + self.vectors = None + self.interactive = interactive + + # This next section is for the jupyter interface. Even if the user has + # specified zero interactivity, the sliders serve as containers for the + # hyperparameters to change when set + self.sliders = {val: + widgets.FloatSlider( + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + min=hyper_vals[val]['options'][0], + max=hyper_vals[val]['options'][1], + description=hyper_vals[val]['name'], + continuous_update=True, + step=(hyper_vals[val]['options'][1] - + hyper_vals[val]['options'][0])/20., + style=style, + ) + if isinstance(hyper_vals[val]['options'][0], float) else + widgets.IntSlider( + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + min=hyper_vals[val]['options'][0], + max=hyper_vals[val]['options'][1], + description=hyper_vals[val]['name'], + continuous_update=True, + step=1, + style=style, + ) + if isinstance(hyper_vals[val]['options'][0], int) and + hyper_vals[val]['options'][0] != True else + widgets.Dropdown(options=hyper_vals[val]['options'], + style=style, + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + description=hyper_vals[val]['name'], + ) + for val in hyper_vals if 'fixed' not in hyper_vals[val]} + + # Which elements to center smooth gaussians on. Defaults to all elements + # present except for hydrogen. + constituents = list(set([s for k in self.symbols for s in self.symbols[k]])) + self.sliders['center_select']=widgets.SelectMultiple( + options=constituents, + value=[s for s in constituents if s != "H"], + description="Where to Place SOAP Centers", + style=style) + + # Whether to compute the global kernel (average of the environments) or + # atomic kernel. + self.sliders['average'] = widgets.Dropdown( + options=["Environment-Centered", "Average"], + value="Environment-Centered", + description="Type of SOAP Vectors", + style=style) + + # Show and enable the widgets if interactive==True + if(interactive): + self.preset_button = _button_template_(list(hyper_dict.keys()), + "SOAP Presets: ", + self.preset_func) + if(isinstance(hyperparameters, str)): + self.preset_button.value = hyperparameters + + slider_order = ['center_select', + 'average', + *list(hyper_vals.keys())] + + for s in slider_order: + if(s in self.sliders): + display(self.sliders[s]) + self.sliders[s].observe( + lambda change: self.change_func(change), names='value') + + # These buttons allow users to choose up to two molecules to compute + # SOAP vectors and compare. + self.molecule_buttons = [widgets.Dropdown(options=[*self.frames.keys()], # , "Other"], + style=style,\ + value=list( + self.frames.keys())[0], + description="\tMolecule: " + ) for i in range(2)] + + def import_molecules(self, file_list): + """ Function to import ASE-like files and include them in the dictionaries + of available molecules. + + Author: R. Cersonsky + + Keywords: + file_list - list of strings corresponding to ASE-like files + """ + for file in file_list: + for frame in read(file, ":"): + key = str(frame.symbols) + + if(key in self.frames): + print("Overwriting previous entry for {}".format(key)) + + self.frames[key] = frame + self.positions[key] = format_positions(frame.get_positions()) + self.inds[key] = sorted(range(len(self.positions[key])), + key=lambda j: self.positions[key][j][0]) + self.positions[key] = self.positions[key][self.inds[key]] + self.symbols[key] = frame.symbols[self.inds[key]] + super_scripts = ['' if list(self.symbols[key]).count(s) == 1 + else '^{{({})}}'.format(list(self.symbols[key])[:i].count(s)+1) + for i, s in enumerate(list(self.symbols[key]))] + self.labels[key] = np.array([r'${}{}$'.format(self.symbols[key][i], s) + for i, s in enumerate(super_scripts)]) + + def change_func(self, change): + """ Function to import ASE-like files and include them in the dictionaries + of available molecules. + + Author: R. Cersonsky + + Keywords: + file_list - list of strings corresponding to ASE-like files + """ + change['owner'].value = change['new'] + for s in self.hyperparameters: + if(s in self.sliders): + if(self.hyperparameters[s] != self.sliders[s].value): + self.hyperparameters[s] = self.sliders[s].value + self.vectors = None + + def preset_func(self, a): + """ Function to change hyperparameters based upon presets + + Author: R. Cersonsky + """ + for s in self.sliders: + if(s in self.hyperparameters): + self.hyperparameters[s] = hyper_dict[self.preset_button.value][s] + self.sliders[s].value = self.hyperparameters[s] + + def get_soap_vectors(self, key=None, average=False): + """ Get the SOAP vectors for the given molecule, denoted by key + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by \ + import_molecules + average - whether to compute the average SOAP vector for the + molecule + """ + if(key == None): + key = self.input_file.replace('./data/molecules/', '')[:-4] + + representation = SOAP(**mask_body_order(self.hyperparameters)) + all_frames = [self.frames[k] for k in sorted(self.frames.keys())] + features = representation.transform( + all_frames).get_dense_feature_matrix(representation) + + lengths = [len(self.inds[k]) for k in sorted(self.frames.keys())] + + self.vectors = {k: features[sum(lengths[0:i]):sum(lengths[0:i])+lengths[i]] + for i, k in enumerate(sorted(self.frames.keys()))} + vectors = self.vectors[key] + + # if(average): + # return np.mean(vectors, axis=0) + return vectors[self.inds[key]] + + def make_figure(self, key, img_name): + """ Function to generate an image of the given molecule and cache it + under img_name + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by + import_molecules + img_name - filename under which to cache the image + """ + from ase.data import covalent_radii, atomic_numbers + from matplotlib import pyplot, patches + from ase.data.colors import jmol_colors + style = { + "horizontalalignment": "center", + "verticalalignment": "center", + "fontsize": 12 + } + fig, ax = pyplot.subplots(1) + + for i, p in enumerate(self.positions[key]): + an = atomic_numbers[self.symbols[key][i]] + ax.text(*p[:2], + s=self.labels[key][i], + **style, + zorder=p[-1]+0.01 + ) + ax.add_artist(patches.Circle(p[:2], + covalent_radii[an], + facecolor=jmol_colors[an], + edgecolor='k', + alpha=0.95, + zorder=p[-1])) + ax.axis('off') + bounds = [*np.min(self.positions[key], axis=0), * + np.max(self.positions[key], axis=0)] + ax.set_xlim([bounds[0]-1.5, bounds[3]+1.5]) + ax.set_ylim([bounds[1]-1.5, bounds[5]+1.5]) + ax.set_aspect('equal') + fig.savefig(img_name) + fig.clf() + + def show_kernel(self, key1, key2=None, show=True, average=False): + from rascal.models import Kernel + from ase.data import chemical_symbols + """ Function to show pandas dataframe for the comparison kernel + + Author: R. Cersonsky + + Keywords: + key1 - string corresponding to molecule built by + import_molecules + key2 - string corresponding to molecule built by + import_molecules, defaults to key1 + show - boolean, whether to show or return the dataframe + average - boolean, whether to compute the average SOAP vector + for the molecule + """ + key2 = key2 if key2 != None else key1 + average = self.sliders['average'].value == "Average" + + soap=SOAP(**mask_body_order(self.hyperparameters)) + features1=soap.transform(self.frames[key1]) + features2=soap.transform(self.frames[key2]) + + kernel = Kernel(soap, + target_type="Structure" if average else "Atom", + zeta = 2, + **self.hyperparameters) + data = kernel(features2, features1) + + if (average): + + vec1 = np.mean(features1.get_dense_feature_matrix(soap),axis=0) + vec2 = np.mean(features2.get_dense_feature_matrix(soap), axis=0) + if(vec1.shape==vec2.shape): + data = np.dot(vec1, vec2)/(np.linalg.norm(vec1)*np.linalg.norm(vec2)) + df=pd.DataFrame(data=[['%1.3f' %(data)]], + columns =[key2], + index =[key1] + ) + else: + df=pd.DataFrame(data=[['%1.3f' %(data[i][j]) + for i,v in enumerate(features2._frames.numbers) if v!=1] + for j,w in enumerate(features1._frames.numbers) if w!=1], + columns =[chemical_symbols[v] for i,v in enumerate(features2._frames.numbers) if v!=1], + index =[chemical_symbols[v] for i,v in enumerate(features1._frames.numbers) if v!=1] + ) + # + df.name = ("Self-" if key1 == key2 else '') + \ + "Similarity " + ("Kernel" if not average else "") + return df + + def show_molecule(self, key, show=True): + """ Function to either generate or show the image for the given molecule + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by + import_molecules + show - boolean, whether to show or return the dataframe + """ + import os + img_name = './images/cache/{}.png'.format(key) + if(not os.path.exists(img_name)): + self.make_figure(key, img_name) + if(show): + display(Image(img_name)) + return img_name + + def compute(self, a=None): + """ Function to compute the SOAP vectors for the set hyperparameters + + Author: R. Cersonsky + """ + def wrap_compute(a): + self.get_soap_vectors( + key=self.molecule_buttons[0].value, average=self.sliders['average'].value == "Average") + print("SOAP Vectors for {}: \n\t".format( + self.molecule_buttons[0].value), self.vectors[self.molecule_buttons[0].value]) + + compute_button = widgets.Button(description="Compute", + style=style, + value=False, + ) + compute_button.on_click(wrap_compute) + clear_button = widgets.Button(description="Clear Output", + style=style, + value=False, + ) + clear_button.on_click(lambda a: self.compute(a)) + display(self.molecule_buttons[0]) + display(widgets.widgets.HBox((compute_button, clear_button))) + + def compare(self, a=None): + """ Function to compare the SOAP vectors of two molecules + for the set hyperparameters + + Author: R. Cersonsky + """ + clear_output() + + def wrap_compare(a): + df = self.show_kernel( + key1=self.molecule_buttons[0].value, key2=self.molecule_buttons[1].value, show=False) + m1 = self.show_molecule(self.molecule_buttons[0].value, show=False) + m2 = self.show_molecule(self.molecule_buttons[1].value, show=False) + display(HTML("
".format(
+ self.molecule_buttons[0].value, df.name, self.molecule_buttons[1].value, m1, df.to_html(), m2)))
+ compare_button = widgets.Button(description="Compare",
+ style=style,
+ value=False,
+ )
+
+ compare_button.on_click(wrap_compare)
+ clear_button = widgets.Button(description="Clear Output",
+ style=style,
+ value=False,
+ )
+ clear_button.on_click(lambda a: self.compare(a))
+ display(widgets.widgets.HBox(
+ (self.molecule_buttons[0], self.molecule_buttons[1])))
+ display(widgets.widgets.HBox((compare_button, clear_button)))
diff --git a/docs/source/tutorials/utilities/widget_style.json b/docs/source/tutorials/utilities/widget_style.json
new file mode 100644
index 000000000..278465f1c
--- /dev/null
+++ b/docs/source/tutorials/utilities/widget_style.json
@@ -0,0 +1 @@
+{"description_width": "initial"}
\ No newline at end of file
diff --git a/examples/Prediction_example.ipynb b/examples/Prediction_example.ipynb
deleted file mode 100644
index b760399ab..000000000
--- a/examples/Prediction_example.ipynb
+++ /dev/null
@@ -1,330 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Using Rascal and SOAP to Predict Properties\n",
- "\n",
- "This notebook is intended as an introductory how-to on calculating the SOAP vector and train a model for their atomization energies on these vectors. For more information on the variable conventions, derivation, utility, and calculation of SOAP vectors, please refer to (among others): \n",
- "- [On representing chemical environments (Bartók 2013)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.87.184115)\n",
- "- [Gaussian approximation potentials: A brief tutorial introduction (Bartók 2015)](https://onlinelibrary.wiley.com/doi/full/10.1002/qua.24927)\n",
- "- [Comparing molecules and solids across structural and alchemical space (De 2016)](https://pubs.rsc.org/en/content/articlepdf/2016/cp/c6cp00415f)\n",
- "\n",
- "Beyond libRascal, the packages used in this tutorial are: [json](https://docs.python.org/2/library/json.html), [numpy](https://numpy.org/), [matplotlib](https://matplotlib.org/), and [ase](https://wiki.fysik.dtu.dk/ase/index.html)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%matplotlib inline\n",
- "%reload_ext autoreload\n",
- "%autoreload 2\n",
- "from tutorial_utils import *\n",
- "try:\n",
- " from rascal.representations import SphericalInvariants as SOAP\n",
- "except:\n",
- " from rascal.representations import SOAP\n",
- "readme_button()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### First, let's look at how we can use SOAP to represent small molecules.\n",
- "We will play around with the SOAP hyperparameters in later examples, but for now we'll use the default hyperparameters."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP=SOAP_tutorial()\n",
- "mySOAP.output_params()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now let's look at the SOAP representation of our first frame"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "soap = SOAP(**mySOAP.hyperparameters)\n",
- "X = soap.transform(mySOAP.frames[:1])\n",
- "\n",
- "mySOAP.verbosity_wrap(\"Our first frame has {} environments, thus our soap feature matrix has a shape of {}\".format(len(mySOAP.frames[0].positions), X.get_feature_matrix().T.shape, len(X.get_feature_matrix().T[0])))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Now we know how to retrieve SOAP vectors, let's look at the impact of the hyperparameters on training a Kernel Ridge Regression (KRR)\n",
- "This time when we open up the tutorial, you will be able to change the input file, hyperparameters, and property to use for the kernel ridge regression, which are saved to mySOAP as they are changed. We've even included some suggestions for hyperparameters, why not try the Power Spectrum first?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP=SOAP_tutorial(interactive=True)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP.train_krr_model()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP.predict_test_set()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP.predict_new_set(filename='../reference_data/inputs/small_molecules-1000.xyz')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Now that we've explained the workflow, let's strip away the SOAP_tutorial wrapper and run the computation again:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from ase.io import read\n",
- "import numpy as np\n",
- "try:\n",
- " from rascal.representations import SphericalInvariants as SOAP\n",
- "except:\n",
- " from rascal.representations import SOAP\n",
- "\n",
- "def split_dataset(N, training_percentage, seed=20):\n",
- " np.random.seed(seed)\n",
- " ids = list(range(N))\n",
- " np.random.shuffle(ids)\n",
- " return ids[:int(training_percentage*N)], ids[int(training_percentage*N):]\n",
- "\n",
- "class KRR(object):\n",
- " def __init__(self, zeta, weights, representation, X, kernel_type):\n",
- " self.weights = weights\n",
- " self.representation = representation\n",
- " self.zeta = zeta\n",
- " self.X = X\n",
- " self.kernel_type=kernel_type\n",
- "\n",
- " def predict(self,frames):\n",
- " features = self.representation.transform(frames)\n",
- " kernel_function = self.X.cosine_kernel_atomic if self.kernel_type=='atomic' else self.X.cosine_kernel_global\n",
- " kernel = kernel_function(features, zeta) \n",
- " return np.dot(self.weights, kernel)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### These are the parameters you'll want change, such as we did above with the sliders.\n",
- "(Everything else in the workflow is a function of these parameters)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "input_file = '../reference_data/inputs/small_molecules-1000.xyz'\n",
- "hyperparameters = dict(soap_type = 'PowerSpectrum', \\\n",
- " interaction_cutoff = 3.5, \\\n",
- " max_radial = 2, \\\n",
- " max_angular = 1, \\\n",
- " gaussian_sigma_constant = 0.5, \\\n",
- " gaussian_sigma_type = 'Constant', \\\n",
- " cutoff_smooth_width = 0.0\n",
- " )\n",
- "property_to_ml = \"dft_formation_energy_per_atom_in_eV\"\n",
- "\n",
- "training_percentage = 0.8\n",
- "zeta = 2\n",
- "Lambda = 5e-3\n",
- "jitter=1e-8"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Which are then used to compute:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "frames = np.array(read(input_file,\":\"))\n",
- "number_of_frames = len(frames)\n",
- "representation = SOAP(**hyperparameters)\n",
- "property_values = np.array([cc.info[property_to_ml] for cc in frames])\n",
- "train_idx, test_idx = split_dataset(number_of_frames, training_percentage)\n",
- "features = representation.transform(frames[train_idx])"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Then, you can construct the kernel for ML and feed it into a KRR Model:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "kernel_function = features.cosine_kernel_atomic if property_to_ml == \"ENERGY\" else features.cosine_kernel_global\n",
- "kernel = kernel_function(zeta)\n",
- "\n",
- "delta = np.std(property_values[train_idx]) / np.mean(kernel.diagonal())\n",
- "kernel[np.diag_indices_from(kernel)] += Lambda**2 / delta **2 + jitter\n",
- "weights = np.linalg.solve(kernel,property_values[train_idx])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "model = KRR(zeta, weights, representation, features, kernel_type='atomic' if property_to_ml==\"ENERGY\" else \"global\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### This model can be in turn used to predict the data from out testing set:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_pred = model.predict(frames[test_idx])\n",
- "print(dict(\n",
- " mean_average_error= [np.mean(np.abs(y_pred-property_values[test_idx]))],\n",
- " root_mean_squared_error=[np.sqrt(np.mean((y_pred-property_values[test_idx])**2))],\n",
- " R2 = [np.mean(1 - (((property_values[test_idx] - y_pred) ** 2).sum(axis=0,dtype=np.float64) / ((property_values[test_idx] - np.average(property_values[test_idx], axis=0) ** 2).sum(axis=0,dtype=np.float64))))]\n",
- " ))\n",
- "plt.scatter(y_pred, property_values[test_idx], s=3)\n",
- "plt.axis('scaled')\n",
- "plt.xlabel('DFT energy / (eV/atom)')\n",
- "plt.ylabel('Predicted energy / (eV/atom)')\n",
- "plt.gca().set_aspect('equal')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Alternatively we can use the model to predict properties from a new data set"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "filename='../reference_data/inputs/small_molecules-1000.xyz'\n",
- "y_pred = model.predict(np.array(read(filename,\":\")))\n",
- "print(dict(\n",
- " mean_average_error= [np.mean(np.abs(y_pred-property_values[test_idx]))],\n",
- " root_mean_squared_error=[np.sqrt(np.mean((y_pred-property_values[test_idx])**2))],\n",
- " R2 = [np.mean(1 - (((property_values[test_idx] - y_pred) ** 2).sum(axis=0,dtype=np.float64) / ((property_values[test_idx] - np.average(property_values[test_idx], axis=0) ** 2).sum(axis=0,dtype=np.float64))))]\n",
- " ))\n",
- "plt.scatter(y_pred, property_values[test_idx], s=3)\n",
- "plt.axis('scaled')\n",
- "plt.xlabel('DFT energy / (eV/atom)')\n",
- "plt.ylabel('Predicted energy / (eV/atom)')\n",
- "plt.gca().set_aspect('equal')"
- ]
- }
- ],
- "metadata": {
- "celltoolbar": "Raw Cell Format",
- "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.8"
- },
- "toc": {
- "base_numbering": 1,
- "nav_menu": {
- "height": "12px",
- "width": "252px"
- },
- "number_sections": true,
- "sideBar": true,
- "skip_h1_title": false,
- "title_cell": "Table of Contents",
- "title_sidebar": "Contents",
- "toc_cell": false,
- "toc_position": {},
- "toc_section_display": "block",
- "toc_window_display": false
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/SOAP_CSD500.ipynb b/examples/SOAP_CSD500.ipynb
deleted file mode 100644
index 76309c747..000000000
--- a/examples/SOAP_CSD500.ipynb
+++ /dev/null
@@ -1,960 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To install rascal:\n",
- "(NOTE: See the top-level README for the most up-to-date installation instructions.)\n",
- "+ mkdir ../build \n",
- "+ cd build\n",
- "+ cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTS=ON ..\n",
- "+ make -j 4\n",
- "+ make install"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Setup"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%matplotlib notebook\n",
- "from matplotlib import pylab as plt\n",
- "import os, sys"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "tags": [
- "hide_input"
- ]
- },
- "outputs": [],
- "source": [
- "# Ensure the notebook is running in the proper directory\n",
- "# (only needed when compiling from sphinx)\n",
- "if \"docs/source/tutorials\" in os.getcwd():\n",
- " os.chdir(\"../../../examples\")\n",
- "sys.path.insert(0,\"../build/\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import time\n",
- "import rascal\n",
- "import json\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",
- "import sys\n",
- "#import pandas\n",
- "import scipy\n",
- "import json\n",
- "\n",
- "from rascal.representations import SphericalInvariants as SOAP"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import sklearn"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Download data"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Run the cell below to dowload the CSD-500 dataset you'll need for this tutorial.\n",
- "\n",
- "The data is from the paper \"Chemical shifts in molecular solids by machine learning\" by Paruzzo, Hofstetter, Musil, De, Ceriotti and Emsley, available [here](https://www.nature.com/articles/s41467-018-06972-x). (If the cell below doesn't work, try downloading \"Supplementary Data 2\" from the [supplementary info section](https://www.nature.com/articles/s41467-018-06972-x#Sec11) of the paper; be sure to name it `data/CSD-500.xyz` as below)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "download_link = \"https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-06972-x/MediaObjects/41467_2018_6972_MOESM4_ESM.txt\"\n",
- "import urllib.request\n",
- "urllib.request.urlretrieve(download_link, './reference_data/CSD-500.xyz')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "#frames = read('./data/small_molecules-1000.xyz',':600')\n",
- "frames = read('./reference_data/CSD-500.xyz',':3')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "n_atoms = sum(frame.get_number_of_atoms() for frame in frames)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "chem_shifts = []\n",
- "for frame in frames:\n",
- " chem_shifts.append(frame.arrays['CS'])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "chem_shifts[1].shape"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "chem_shifts_atoms = np.concatenate(chem_shifts, axis=0)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "chem_shifts_atoms.shape"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# SOAP: Power spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": []
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=6, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%time representation = soap.transform(frames)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "X = representation.get_features
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "X.shape"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Learning the chemical shifts of a set of crystal structures"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## learning utilities"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": [
- 13,
- 20,
- 31,
- 33,
- 35,
- 37,
- 47,
- 54
- ],
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def compute_representation(representation, frames):\n",
- " result = representation.transform(frames)\n",
- " return result\n",
- "\n",
- "#def compute_representation(frames):\n",
- "# expansions = soap.transform(frames)\n",
- "# return expansions\n",
- "\n",
- "def compute_atomic_kernel(zeta, rep1, rep2=None):\n",
- " if rep2 is not None:\n",
- " kernel = rep1.cosine_kernel_atomic(rep2, zeta)\n",
- " else:\n",
- " kernel = rep1.cosine_kernel_atomic(zeta)\n",
- " return kernel\n",
- "\n",
- "def extract_energy(frames):\n",
- " prop = [[]]*len(frames)\n",
- " for ii,cc in enumerate(frames):\n",
- " #prop[ii] = cc.info['dft_formation_energy_per_atom_in_eV']\n",
- " prop[ii] = cc.info['ENERGY']\n",
- " y = np.array(prop)\n",
- " return y\n",
- "\n",
- "def split_dataset(frames, test_fraction, seed=10):\n",
- " N = len(frames)\n",
- " ids = np.arange(N)\n",
- " np.random.seed(seed)\n",
- " np.random.shuffle(ids)\n",
- " Ntrain = int(N*test_fraction)\n",
- " train = ids[:Ntrain]\n",
- " test = ids[Ntrain:]\n",
- " targets = extract_energy(frames)\n",
- " return [frames[ii] for ii in train],targets[train],[frames[ii] for ii in test],targets[test]\n",
- "\n",
- "def get_mae(ypred,y):\n",
- " return np.mean(np.abs(ypred-y))\n",
- "def get_rmse(ypred,y):\n",
- " return np.sqrt(np.mean((ypred-y)**2))\n",
- "def get_sup(ypred,y):\n",
- " return np.amax(np.abs((ypred-y)))\n",
- "def get_r2(y_pred,y_true):\n",
- " weight = 1\n",
- " sample_weight = None\n",
- " numerator = (weight * (y_true - y_pred) ** 2).sum(axis=0,dtype=np.float64)\n",
- " denominator = (weight * (y_true - np.average(\n",
- " y_true, axis=0, weights=sample_weight)) ** 2).sum(axis=0,dtype=np.float64)\n",
- " output_scores = 1 - (numerator / denominator)\n",
- " return np.mean(output_scores)\n",
- "\n",
- "\n",
- "score_func = dict(\n",
- " MAE=get_mae,\n",
- " RMSE=get_rmse,\n",
- " SUP=get_sup,\n",
- " R2=get_r2,\n",
- ")\n",
- "\n",
- "def get_score(ypred,y):\n",
- " scores = {}\n",
- " for k,func in score_func.items():\n",
- " scores[k] = func(ypred,y)\n",
- " return scores\n",
- "\n",
- "class KRR(object):\n",
- " def __init__(self,zeta,weights,representation,X):\n",
- " self.weights = weights\n",
- " self.representation = representation\n",
- " self.zeta = zeta\n",
- " self.X = X\n",
- " \n",
- " def predict(self,frames):\n",
- " features = compute_representation(self.representation,frames)\n",
- " kernel = compute_atomic_kernel(self.zeta, self.X, features)\n",
- " return np.dot(self.weights, kernel)\n",
- "\n",
- "def train_krr_model(zeta,Lambda,representation,frames,y,jitter=1e-8):\n",
- " features = compute_representation(representation, frames)\n",
- " kernel = compute_atomic_kernel(zeta, features, features) \n",
- " # adjust the kernel so that it is properly scaled\n",
- " delta = np.std(y) / np.mean(kernel.diagonal())\n",
- " kernel[np.diag_indices_from(kernel)] += Lambda**2 / delta **2 + jitter\n",
- " # train the krr model\n",
- " weights = np.linalg.solve(kernel,y)\n",
- " model = KRR(zeta, weights,representation, features)\n",
- " return model,kernel\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "compute_representation(soap, frames)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## With the full power spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=1, \n",
- " max_angular=1, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "frames_train, y_train, frames_test, y_test = split_dataset(frames,0.5)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_train = np.concatenate([frame.arrays['CS'] for frame in frames_train], axis=0)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_train = y_train[:,0]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_test = np.concatenate([frame.arrays['CS'] for frame in frames_test], axis=0)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_test = y_test[:,0]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "n_atoms = sum(frame.get_number_of_atoms() for frame in frames_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "zeta = 2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "features = compute_representation(soap, frames_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "features"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "kernel = compute_atomic_kernel(zeta, features, features)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "kernel.shape"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# WTF?\n",
- "#for i in range(len(kernel)):\n",
- "# kernel[i][0] = np.sqrt(i)/50"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "kernel[191]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# First, try without regularization -- the results will be nonsense\n",
- "weights = np.linalg.solve(kernel,y_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "weights.shape"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "model = KRR(zeta, weights, soap, features)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "y_pred = model.predict(frames_test)\n",
- "get_score(y_pred, y_test)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "Lambda = 10\n",
- "krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "y_pred = krr.predict(frames_test)\n",
- "get_score(y_pred, y_test)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "#sc = plt.scatter(y_pred, y_test, s=3)\n",
- "#ax = plt.axis('scaled')\n",
- "plt.plot(y_test, y_pred, '.')\n",
- "plt.ylabel('DFT shift')\n",
- "plt.xlabel('Predicted shift')\n",
- "plt.savefig('R1.png', dpi=300)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## With just the radial spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"RadialSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=0, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "frames_train, y_train, frames_test, y_testr = split_dataset(frames,0.4)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "Lambda = 5e-3\n",
- "krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "y_predr = krr.predict(frames_test)\n",
- "get_score(y_predr, y_testr)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "#plt.scatter(y_predr, y_testr, s=3)\n",
- "#plt.scatter(y_predr, y_testr, s=3)\n",
- "#ax = plt.axis('scaled')\n",
- "plt.plot(y_pred, y_test, '.b')\n",
- "plt.plot(y_predr, y_testr, '.y')\n",
- "plt.legend(['Full','Radial'])\n",
- "plt.ylabel('DFT energy / (eV/atom)')\n",
- "plt.xlabel('Predicted energy / (eV/atom)')\n",
- "plt.savefig('R1.png', dpi=300)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Make a map of the dataset"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## utils"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def compute_representation(representation,frames):\n",
- " expansions = soap.transform(frames)\n",
- " return expansions\n",
- "\n",
- "def compute_kernel(zeta, rep1, rep2=None):\n",
- " if rep2 is None:\n",
- " kernel = rep1.cosine_kernel_global(zeta)\n",
- " else:\n",
- " kernel = rep1.cosine_kernel_global(rep2,zeta)\n",
- " return kernel"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": [
- 0
- ],
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def link_ngl_wdgt_to_ax_pos(ax, pos, ngl_widget):\n",
- " from matplotlib.widgets import AxesWidget\n",
- " from scipy.spatial import cKDTree\n",
- " r\"\"\"\n",
- " Initial idea for this function comes from @arose, the rest is @gph82 and @clonker\n",
- " \"\"\"\n",
- " \n",
- " kdtree = cKDTree(pos) \n",
- " #assert ngl_widget.trajectory_0.n_frames == pos.shape[0]\n",
- " x, y = pos.T\n",
- " \n",
- " lineh = ax.axhline(ax.get_ybound()[0], c=\"black\", ls='--')\n",
- " linev = ax.axvline(ax.get_xbound()[0], c=\"black\", ls='--')\n",
- " dot, = ax.plot(pos[0,0],pos[0,1], 'o', c='red', ms=7)\n",
- "\n",
- " ngl_widget.isClick = False\n",
- " \n",
- " def onclick(event):\n",
- " linev.set_xdata((event.xdata, event.xdata))\n",
- " lineh.set_ydata((event.ydata, event.ydata))\n",
- " data = [event.xdata, event.ydata]\n",
- " _, index = kdtree.query(x=data, k=1)\n",
- " dot.set_xdata((x[index]))\n",
- " dot.set_ydata((y[index]))\n",
- " ngl_widget.isClick = True\n",
- " ngl_widget.frame = index\n",
- " \n",
- " def my_observer(change):\n",
- " r\"\"\"Here comes the code that you want to execute\n",
- " \"\"\"\n",
- " ngl_widget.isClick = False\n",
- " _idx = change[\"new\"]\n",
- " try:\n",
- " dot.set_xdata((x[_idx]))\n",
- " dot.set_ydata((y[_idx])) \n",
- " except IndexError as e:\n",
- " dot.set_xdata((x[0]))\n",
- " dot.set_ydata((y[0]))\n",
- " print(\"caught index error with index %s (new=%s, old=%s)\" % (_idx, change[\"new\"], change[\"old\"]))\n",
- " \n",
- " # Connect axes to widget\n",
- " axes_widget = AxesWidget(ax)\n",
- " axes_widget.connect_event('button_release_event', onclick)\n",
- " \n",
- " # Connect widget to axes\n",
- " ngl_widget.observe(my_observer, \"frame\", \"change\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## make a map with kernel pca projection"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "# Load the small molecules \n",
- "frames = read('./reference_data/small_molecules-1000.xyz',':600')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=6, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "\n",
- "features = compute_representation(soap, frames)\n",
- "\n",
- "kernel = compute_kernel(zeta,features)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "from sklearn.decomposition import KernelPCA"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "kpca = KernelPCA(n_components=2,kernel='precomputed')\n",
- "kpca.fit(kernel)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "X = kpca.transform(kernel)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "plt.scatter(X[:,0],X[:,1],s=3)\n",
- "#plt.savefig('PCA.png',dpi=300)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## make an interactive map"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# package to visualize the structures in the notebook\n",
- "# https://github.com/arose/nglview#released-version\n",
- "import nglview"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "iwdg = nglview.show_asetraj(frames)\n",
- "# set up the visualization\n",
- "iwdg.add_unitcell()\n",
- "iwdg.add_spacefill()\n",
- "iwdg.remove_ball_and_stick()\n",
- "iwdg.camera = 'orthographic'\n",
- "iwdg.parameters = { \"clipDist\": 0 }\n",
- "iwdg.center()\n",
- "iwdg.update_spacefill(radiusType='covalent',\n",
- " scale=0.6,\n",
- " color_scheme='element')\n",
- "iwdg._remote_call('setSize', target='Widget',\n",
- " args=['%dpx' % (600,), '%dpx' % (400,)])\n",
- "iwdg.player.delay = 200.0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "link_ngl_wdgt_to_ax_pos(plt.gca(), X, iwdg)\n",
- "plt.scatter(X[:,0],X[:,1],s=3)\n",
- "iwdg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "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.7.5rc1"
- },
- "toc": {
- "base_numbering": 1,
- "nav_menu": {
- "height": "12px",
- "width": "252px"
- },
- "number_sections": true,
- "sideBar": true,
- "skip_h1_title": false,
- "title_cell": "Table of Contents",
- "title_sidebar": "Contents",
- "toc_cell": false,
- "toc_position": {},
- "toc_section_display": "block",
- "toc_window_display": false
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/SOAP_example.ipynb b/examples/SOAP_example.ipynb
deleted file mode 100644
index b7859a4e4..000000000
--- a/examples/SOAP_example.ipynb
+++ /dev/null
@@ -1,296 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Using Rascal to Calculate SOAP Vectors\n",
- "\n",
- "This notebook is intended as an introductory how-to on calculating the SOAP vector and train a model for their atomization energies on these vectors. For more information on the variable conventions, derivation, utility, and calculation of SOAP vectors, please refer to (among others): \n",
- "- [On representing chemical environments (Bartók 2013)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.87.184115)\n",
- "- [Gaussian approximation potentials: A brief tutorial introduction (Bartók 2015)](https://onlinelibrary.wiley.com/doi/full/10.1002/qua.24927)\n",
- "- [Comparing molecules and solids across structural and alchemical space (De 2016)](https://pubs.rsc.org/en/content/articlepdf/2016/cp/c6cp00415f)\n",
- "\n",
- "Beyond libRascal, the packages used in this tutorial are: [json](https://docs.python.org/2/library/json.html), [numpy](https://numpy.org/), [matplotlib](https://matplotlib.org/), and [ase](https://wiki.fysik.dtu.dk/ase/index.html)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%matplotlib inline\n",
- "%reload_ext autoreload\n",
- "%autoreload 2\n",
- "import sys\n",
- "sys.path.insert(0, '../build')\n",
- "from tutorial_utils import *\n",
- "try:\n",
- " from rascal.representations import SphericalInvariants as SOAP\n",
- "except:\n",
- " from rascal.representations import SOAP\n",
- "readme_button()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### First, let's look at how we can use SOAP to represent small molecules.\n",
- "We will play around with the SOAP hyperparameters in later examples, but for now we'll use the default hyperparameters."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP=SOAP_tutorial()\n",
- "mySOAP.output_params()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now let's look at the SOAP representation of our first frame"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "soap = SOAP(**mySOAP.hyperparameters)\n",
- "X = soap.transform(mySOAP.frames[:1])\n",
- "\n",
- "mySOAP.verbosity_wrap(\"Our first frame has {} environments, thus our soap feature matrix has a shape of {}\".format(len(mySOAP.frames[0].positions), X.get_features(soap).T.shape, len(X.get_features(soap).T[0])))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Now we know how to retrieve SOAP vectors, let's look at the impact of the hyperparameters on training a Kernel Ridge Regression (KRR)\n",
- "This time when we open up the tutorial, you will be able to change the input file, hyperparameters, and property to use for the kernel ridge regression, which are saved to mySOAP as they are changed. We've even included some suggestions for hyperparameters, why not try the Power Spectrum first?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP=SOAP_tutorial(interactive=True)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP.train_krr_model()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mySOAP.plot_prediction_func()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Now that we've explained the workflow, let's strip away the SOAP_tutorial wrapper and run the computation again:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from ase.io import read\n",
- "import numpy as np\n",
- "try:\n",
- " from rascal.representations import SphericalInvariants as SOAP\n",
- "except:\n",
- " from rascal.representations import SOAP\n",
- "\n",
- "def split_dataset(N, training_percentage, seed=20):\n",
- " np.random.seed(seed)\n",
- " ids = list(range(N))\n",
- " np.random.shuffle(ids)\n",
- " return ids[:int(training_percentage*N)], ids[int(training_percentage*N):]\n",
- "\n",
- "class KRR(object):\n",
- " def __init__(self, zeta, weights, representation, X, kernel_type):\n",
- " self.weights = weights\n",
- " self.representation = representation\n",
- " self.zeta = zeta\n",
- " self.X = X\n",
- " self.kernel_type=kernel_type\n",
- "\n",
- " def predict(self,frames):\n",
- " features = self.representation.transform(frames)\n",
- " kernel_function = self.X.cosine_kernel_atomic if self.kernel_type=='atomic' else self.X.cosine_kernel_global\n",
- " kernel = kernel_function(features, zeta) \n",
- " return np.dot(self.weights, kernel)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### These are the parameters you'll want change, such as we did above with the sliders.\n",
- "(Everything else in the workflow is a function of these parameters)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "input_file = 'reference_data/small_molecules-1000.xyz'\n",
- "hyperparameters = dict(soap_type = 'PowerSpectrum', \\\n",
- " interaction_cutoff = 3.5, \\\n",
- " max_radial = 2, \\\n",
- " max_angular = 1, \\\n",
- " gaussian_sigma_constant = 0.5, \\\n",
- " gaussian_sigma_type = 'Constant', \\\n",
- " cutoff_smooth_width = 0.0\n",
- " )\n",
- "property_to_ml = \"dft_formation_energy_per_atom_in_eV\"\n",
- "\n",
- "training_percentage = 0.8\n",
- "zeta = 2\n",
- "Lambda = 5e-3\n",
- "jitter=1e-8"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Which are then used to compute:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "frames = np.array(read(input_file,\":\"))\n",
- "number_of_frames = len(frames)\n",
- "representation = SOAP(**hyperparameters)\n",
- "property_values = np.array([cc.info[property_to_ml] for cc in frames])\n",
- "train_idx, test_idx = split_dataset(number_of_frames, training_percentage)\n",
- "features = representation.transform(frames[train_idx])"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Then, you can construct the kernel for ML and feed it into a KRR Model:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "kernel_function = features.cosine_kernel_atomic if property_to_ml == \"ENERGY\" else features.cosine_kernel_global\n",
- "kernel = kernel_function(zeta)\n",
- "\n",
- "delta = np.std(property_values[train_idx]) / np.mean(kernel.diagonal())\n",
- "kernel[np.diag_indices_from(kernel)] += Lambda**2 / delta **2 + jitter\n",
- "weights = np.linalg.solve(kernel,property_values[train_idx])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "model = KRR(zeta, weights, representation, features, kernel_type='atomic' if property_to_ml==\"ENERGY\" else \"global\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### This model can be in turn used to predict the data from out testing set:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "y_pred = model.predict(frames[test_idx])\n",
- "print(dict(\n",
- " mean_average_error= [np.mean(np.abs(y_pred-property_values[test_idx]))],\n",
- " root_mean_squared_error=[np.sqrt(np.mean((y_pred-property_values[test_idx])**2))],\n",
- " R2 = [np.mean(1 - (((property_values[test_idx] - y_pred) ** 2).sum(axis=0,dtype=np.float64) / ((property_values[test_idx] - np.average(property_values[test_idx], axis=0) ** 2).sum(axis=0,dtype=np.float64))))]\n",
- " ))\n",
- "plt.scatter(y_pred, property_values[test_idx], s=3)\n",
- "plt.axis('scaled')\n",
- "plt.xlabel('DFT energy / (eV/atom)')\n",
- "plt.ylabel('Predicted energy / (eV/atom)')\n",
- "plt.gca().set_aspect('equal')"
- ]
- }
- ],
- "metadata": {
- "celltoolbar": "Raw Cell Format",
- "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.8"
- },
- "toc": {
- "base_numbering": 1,
- "nav_menu": {
- "height": "12px",
- "width": "252px"
- },
- "number_sections": true,
- "sideBar": true,
- "skip_h1_title": false,
- "title_cell": "Table of Contents",
- "title_sidebar": "Contents",
- "toc_cell": false,
- "toc_position": {},
- "toc_section_display": "block",
- "toc_window_display": false
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/SOAP_small_molecules.ipynb b/examples/SOAP_small_molecules.ipynb
deleted file mode 100644
index e47ebe7b0..000000000
--- a/examples/SOAP_small_molecules.ipynb
+++ /dev/null
@@ -1,694 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To install rascal:\n",
- "(NOTE: See the top-level README for the most up-to-date installation instructions.)\n",
- "+ mkdir ../build \n",
- "+ cd build\n",
- "+ cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTS=ON ..\n",
- "+ make -j 4\n",
- "+ make install"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%matplotlib notebook\n",
- "from matplotlib import pylab as plt\n",
- "\n",
- "import os, sys\n",
- "from ase.io import read\n",
- "sys.path.insert(0,\"../build/\")\n",
- "\n",
- "import sys\n",
- "import time\n",
- "import rascal\n",
- "import json\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",
- "import sys\n",
- "\n",
- "import json\n",
- "\n",
- "from rascal.representations import SOAP"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "frames = read('../reference_data/inputs/reference_data/small_molecules-1000.xyz',':')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# SOAP: Power spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": []
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=4, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "%time representation = soap.transform(frames)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "X = representation.get_feature_matrix().T"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "X.shape"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Learning the formation energies of small molecules"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": []
- },
- "outputs": [],
- "source": [
- "# Load the small molecules \n",
- "frames = read('../reference_data/inputs/small_molecules-1000.xyz',':600')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## learning utilities"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": [
- 13,
- 20,
- 31,
- 33,
- 35,
- 37,
- 47,
- 54
- ],
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def compute_representation(representation,frames):\n",
- " expansions = soap.transform(frames)\n",
- " return expansions\n",
- "\n",
- "def compute_kernel(zeta, rep1, rep2=None):\n",
- " if rep2 is None:\n",
- " kernel = rep1.cosine_kernel_global(zeta)\n",
- " else:\n",
- " kernel = rep1.cosine_kernel_global(rep2,zeta)\n",
- " return kernel\n",
- "\n",
- "def extract_energy(frames):\n",
- " prop = [[]]*len(frames)\n",
- " for ii,cc in enumerate(frames):\n",
- " prop[ii] = cc.info['dft_formation_energy_per_atom_in_eV']\n",
- " y = np.array(prop)\n",
- " return y\n",
- "\n",
- "def split_dataset(frames, test_fraction, seed=10):\n",
- " N = len(frames)\n",
- " ids = np.arange(N)\n",
- " np.random.seed(seed)\n",
- " np.random.shuffle(ids)\n",
- " Ntrain = int(N*test_fraction)\n",
- " train = ids[:Ntrain]\n",
- " test = ids[Ntrain:]\n",
- " targets = extract_energy(frames)\n",
- " return [frames[ii] for ii in train],targets[train],[frames[ii] for ii in test],targets[test]\n",
- "\n",
- "def get_mae(ypred,y):\n",
- " return np.mean(np.abs(ypred-y))\n",
- "def get_rmse(ypred,y):\n",
- " return np.sqrt(np.mean((ypred-y)**2))\n",
- "def get_sup(ypred,y):\n",
- " return np.amax(np.abs((ypred-y)))\n",
- "def get_r2(y_pred,y_true):\n",
- " weight = 1\n",
- " sample_weight = None\n",
- " numerator = (weight * (y_true - y_pred) ** 2).sum(axis=0,dtype=np.float64)\n",
- " denominator = (weight * (y_true - np.average(\n",
- " y_true, axis=0, weights=sample_weight)) ** 2).sum(axis=0,dtype=np.float64)\n",
- " output_scores = 1 - (numerator / denominator)\n",
- " return np.mean(output_scores)\n",
- "\n",
- "\n",
- "score_func = dict(\n",
- " MAE=get_mae,\n",
- " RMSE=get_rmse,\n",
- " SUP=get_sup,\n",
- " R2=get_r2,\n",
- ")\n",
- "\n",
- "def get_score(ypred,y):\n",
- " scores = {}\n",
- " for k,func in score_func.items():\n",
- " scores[k] = func(ypred,y)\n",
- " return scores\n",
- "\n",
- "class KRR(object):\n",
- " def __init__(self,zeta,weights,representation,X):\n",
- " self.weights = weights\n",
- " self.representation = representation\n",
- " self.zeta = zeta\n",
- " self.X = X\n",
- " \n",
- " def predict(self,frames):\n",
- " features = compute_representation(self.representation,frames)\n",
- " kernel = compute_kernel(self.zeta , self.X, features)\n",
- " return np.dot(self.weights, kernel)\n",
- " \n",
- "def train_krr_model(zeta,Lambda,representation,frames,y,jitter=1e-8):\n",
- " features = compute_representation(representation,frames)\n",
- " kernel = compute_kernel(zeta,features) \n",
- " # adjust the kernel so that it is properly scaled\n",
- " delta = np.std(y) / np.mean(kernel.diagonal())\n",
- " kernel[np.diag_indices_from(kernel)] += Lambda**2 / delta **2 + jitter\n",
- " # train the krr model\n",
- " weights = np.linalg.solve(kernel,y)\n",
- " model = KRR(zeta, weights,representation, features)\n",
- " return model,kernel\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## With the full power spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=6, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "frames_train, y_train, frames_test, y_test = split_dataset(frames,0.8)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "Lambda = 5e-3\n",
- "krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "y_pred = krr.predict(frames_test)\n",
- "get_score(y_pred, y_test)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "plt.scatter(y_pred, y_test, s=3)\n",
- "plt.axis('scaled')\n",
- "plt.xlabel('DFT energy / (eV/atom)')\n",
- "plt.ylabel('Predicted energy / (eV/atom)')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## With just the radial spectrum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"RadialSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=0, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "frames_train, y_train, frames_test, y_test = split_dataset(frames,0.8)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "Lambda = 5e-4\n",
- "krr,k = train_krr_model(zeta, Lambda, soap, frames_train, y_train)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "y_pred = krr.predict(frames_test)\n",
- "get_score(y_pred, y_test)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "plt.scatter(y_pred, y_test, s=3)\n",
- "plt.axis('scaled')\n",
- "plt.xlabel('DFT energy / (eV/atom)')\n",
- "plt.ylabel('Predicted energy / (eV/atom)')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Make a map of the dataset"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## utils"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def compute_representation(representation,frames):\n",
- " expansions = soap.transform(frames)\n",
- " return expansions\n",
- "\n",
- "def compute_kernel(zeta, rep1, rep2=None):\n",
- " if rep2 is None:\n",
- " kernel = rep1.cosine_kernel_global(zeta)\n",
- " else:\n",
- " kernel = rep1.cosine_kernel_global(rep2,zeta)\n",
- " return kernel"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "code_folding": [
- 0
- ],
- "hidden": true
- },
- "outputs": [],
- "source": [
- "def link_ngl_wdgt_to_ax_pos(ax, pos, ngl_widget):\n",
- " from matplotlib.widgets import AxesWidget\n",
- " from scipy.spatial import cKDTree\n",
- " r\"\"\"\n",
- " Initial idea for this function comes from @arose, the rest is @gph82 and @clonker\n",
- " \"\"\"\n",
- " \n",
- " kdtree = cKDTree(pos) \n",
- " #assert ngl_widget.trajectory_0.n_frames == pos.shape[0]\n",
- " x, y = pos.T\n",
- " \n",
- " lineh = ax.axhline(ax.get_ybound()[0], c=\"black\", ls='--')\n",
- " linev = ax.axvline(ax.get_xbound()[0], c=\"black\", ls='--')\n",
- " dot, = ax.plot(pos[0,0],pos[0,1], 'o', c='red', ms=7)\n",
- "\n",
- " ngl_widget.isClick = False\n",
- " \n",
- " def onclick(event):\n",
- " linev.set_xdata((event.xdata, event.xdata))\n",
- " lineh.set_ydata((event.ydata, event.ydata))\n",
- " data = [event.xdata, event.ydata]\n",
- " _, index = kdtree.query(x=data, k=1)\n",
- " dot.set_xdata((x[index]))\n",
- " dot.set_ydata((y[index]))\n",
- " ngl_widget.isClick = True\n",
- " ngl_widget.frame = index\n",
- " \n",
- " def my_observer(change):\n",
- " r\"\"\"Here comes the code that you want to execute\n",
- " \"\"\"\n",
- " ngl_widget.isClick = False\n",
- " _idx = change[\"new\"]\n",
- " try:\n",
- " dot.set_xdata((x[_idx]))\n",
- " dot.set_ydata((y[_idx])) \n",
- " except IndexError as e:\n",
- " dot.set_xdata((x[0]))\n",
- " dot.set_ydata((y[0]))\n",
- " print(\"caught index error with index %s (new=%s, old=%s)\" % (_idx, change[\"new\"], change[\"old\"]))\n",
- " \n",
- " # Connect axes to widget\n",
- " axes_widget = AxesWidget(ax)\n",
- " axes_widget.connect_event('button_release_event', onclick)\n",
- " \n",
- " # Connect widget to axes\n",
- " ngl_widget.observe(my_observer, \"frame\", \"change\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "heading_collapsed": true
- },
- "source": [
- "## make a map with kernel pca projection"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "# Load the small molecules \n",
- "frames = read('./reference_data/small_molecules-1000.xyz',':600')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "hypers = dict(soap_type=\"PowerSpectrum\",\n",
- " interaction_cutoff=3.5, \n",
- " max_radial=6, \n",
- " max_angular=6, \n",
- " gaussian_sigma_constant=0.4,\n",
- " gaussian_sigma_type=\"Constant\",\n",
- " cutoff_smooth_width=0.5,\n",
- " )\n",
- "soap = SOAP(**hypers)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "zeta = 2\n",
- "\n",
- "features = compute_representation(soap, frames)\n",
- "\n",
- "kernel = compute_kernel(zeta,features)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "from sklearn.decomposition import KernelPCA"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "kpca = KernelPCA(n_components=2,kernel='precomputed')\n",
- "kpca.fit(kernel)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "X = kpca.transform(kernel)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "hidden": true
- },
- "outputs": [],
- "source": [
- "plt.scatter(X[:,0],X[:,1],s=3)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## make an interactive map"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# package to visualize the structures in the notebook\n",
- "# https://github.com/arose/nglview#released-version\n",
- "import nglview"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "iwdg = nglview.show_asetraj(frames)\n",
- "# set up the visualization\n",
- "iwdg.add_unitcell()\n",
- "iwdg.add_spacefill()\n",
- "iwdg.remove_ball_and_stick()\n",
- "iwdg.camera = 'orthographic'\n",
- "iwdg.parameters = { \"clipDist\": 0 }\n",
- "iwdg.center()\n",
- "iwdg.update_spacefill(radiusType='covalent',\n",
- " scale=0.6,\n",
- " color_scheme='element')\n",
- "iwdg._remote_call('setSize', target='Widget',\n",
- " args=['%dpx' % (600,), '%dpx' % (400,)])\n",
- "iwdg.player.delay = 200.0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "link_ngl_wdgt_to_ax_pos(plt.gca(), X, iwdg)\n",
- "plt.scatter(X[:,0],X[:,1],s=3)\n",
- "iwdg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "celltoolbar": "Initialization Cell",
- "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.8"
- },
- "toc": {
- "base_numbering": 1,
- "nav_menu": {
- "height": "12px",
- "width": "252px"
- },
- "number_sections": true,
- "sideBar": true,
- "skip_h1_title": false,
- "title_cell": "Table of Contents",
- "title_sidebar": "Contents",
- "toc_cell": false,
- "toc_position": {},
- "toc_section_display": "block",
- "toc_window_display": false
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/reference_data/.gitignore b/examples/reference_data/.gitignore
deleted file mode 100644
index e9322b753..000000000
--- a/examples/reference_data/.gitignore
+++ /dev/null
@@ -1,2 +0,0 @@
-# Large file, should be downloaded by user instead
-CSD-500.xyz
diff --git a/shared_hooks/pre-commit b/shared_hooks/pre-commit
index 4fecbbbb6..599161338 100755
--- a/shared_hooks/pre-commit
+++ b/shared_hooks/pre-commit
@@ -9,7 +9,7 @@
# To use this hook, execute:
# git config core.hooksPath shared_hooks
-NBS=$(echo $(git diff --full-index --cached $against --name-only) | sort | grep "ipynb")
+NBS=$(echo $(git diff --full-index --cached $against --name-only) | sort | grep "ipynb" | grep "docs\/source\/tutorials" -v)
if [ ! -z "$NBS" ] ; then
command -v nbstripout --is-installed >/dev/null 2>&1 || { echo >&2 "I require nbstripout but it's not installed. Aborting commit."; exit 1; }
Table of Contents
\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Title of the Tutorial\n", + "\n", + "This notebook is intended to _purpose of tutorial_.\n", + "For more information on _nuances of tutorial_, please refer to (among others): \n", + "- [Title (Surname Year)](url) \n", + "- [Title (Surname Year)](url)\n", + "\n", + "Beyond libRascal, the packages used in this tutorial are: [pkg](url), [pkg](url), and [pkg](url)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import time\n", + "\n", + "import numpy as np\n", + "import ase\n", + "from ase.io import read\n", + "from matplotlib import pyplot as plt\n", + "import json\n", + "from IPython.display import display, Markdown\n", + "\n", + "from rascal.representations import SphericalInvariants as SOAP\n", + "from rascal.models import Kernel\n", + "\n", + "sys.path.append('./utilities')\n", + "from general_utils import *" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": true, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst new file mode 100644 index 000000000..6c6da0757 --- /dev/null +++ b/docs/source/tutorials/index.rst @@ -0,0 +1,38 @@ +.. _tutorials: + +Tutorials +========= + +Here you will find some tutorial examples for getting started with librascal, +both as a new user of the library and for getting started developing a new +representation. + + +Introductory Tutorials +~~~~~~~~~~~~~~~ + +.. toctree:: + :maxdepth: 1 + + Introduction_to_Atomic_Descriptors.ipynb + Hyperparameters.ipynb + Property_Prediction_Tutorial.ipynb + Kernels.ipynb + Benefits_Performance.ipynb + Advanced_SOAP.ipynb + +Converting Your Workflow to Rascal +~~~~~~~~~~~~~~~~~ + +.. toctree:: + :maxdepth: 1 + + Converting_Quip.ipynb + +Rascal Cookbook +~~~~~~~~~~~~~~~~~~~~ + +.. toctree:: + :maxdepth: 1 + + nl-for-user diff --git a/docs/source/tutorials/nl-for-user.rst b/docs/source/tutorials/nl-for-user.rst new file mode 100644 index 000000000..dff3e8e2b --- /dev/null +++ b/docs/source/tutorials/nl-for-user.rst @@ -0,0 +1,109 @@ +.. _nl-for-user: + +How to build a neighbor list? +============================= + +Overview +~~~~~~~~ + + +Here we provide a tutorial on how to use the C++ interface for building neighbor lists. + +Rascal comes with a built-in neighborhood manager. The role of a neighborhood manager is to minimize manual handeling of building the neighborhood lists. +The current implementation avoids the explicit use of embedded lists to describe neighborhoods and makes iterations over them more trivial. +The neighborhood manager interface is universal to all the representations of Rascal. It can also be very easily interfaced with MD/MC simulation packages. +We prefer using `Eigen| {} | '.format(h) for h in headers]))+''.join(['
|---|
| {} | {}{} | '.format(round(v,4)) for v in d[key]])) for key in d])+'
We will now train a model on {}.".format( + self.sliders['property_to_ml'].value)) + + representation = SOAP(**mask_body_order(self.hyperparameters)) + + if(known_properties[self.sliders['property_to_ml'].value] == 'Atom'): + props = np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[:,0] + else: + props = self.properties[self.sliders['property_to_ml'].value][frame_idx] + + training_properties=props + + verbosity_wrap("First, I am going to separate my dataset:") + verbosity_wrap("
Now we will compute the SOAP representation of \ + our training frames.") + + t=time.time() + features=representation.transform(self.frames[frame_idx]) + + verbosity_wrap('This took {} seconds/frame.'.format( + round((time.time()-t)/len(frame_idx), 8))) + + if(pretend == False): + self.estimate_time(N=max(0, 20-len(self.est_frames)), + p="Estimating time to compute kernel...") + est=int(np.poly1d(np.polyfit(self.est_frames, + self.est_times, + deg=2))(len(frame_idx)))+1 + else: + est=0 + + verbosity_wrap("
Next we find the kernel for our training model.\ +
(This step will take approximately {} minutes and \ + {} seconds.)".format(int(est/60), int(est % 60))) + + time.sleep(0.5) + t=time.time() + kernel = compute_kernel(calculator=representation, \ + features1=features, \ + kernel_type=self.sliders['kernel_type'].value, \ + **self.hyperparameters) + + self.est_frames.append(len(frame_idx)) + self.est_times.append(time.time()-t) + + delta=np.std(training_properties) / np.mean(kernel.diagonal()) + kernel[np.diag_indices_from( + kernel)] += self.Lambda**2 / delta ** 2 + jitter + + verbosity_wrap("
We will adjust the our kernel with the tolerance \ + matrix $\\Lambda = ({})I$.".format(\ + round(self.Lambda**2 / delta ** 2 + jitter, 8))) + verbosity_wrap("
Now we can take this kernel to compute the weights \ + of our KRR.") + + weights=np.linalg.solve(kernel, training_properties) + model=KRR(weights, features, \ + kernel_type=self.sliders['kernel_type'].value, + **self.hyperparameters + ) + if(pretend == False): + self.krr[self.sliders['property_to_ml'].value], k=model, kernel + self.trained[self.sliders['property_to_ml'].value]=True + + def train_krr_model(self, jitter= 1e-8): + """ Wrapper to the train_krr_model_func function for use in the tutorial + + Author: R. Cersonsky + + Keywords: + jitter - float, anticipated error to adjust kernel diagonals + + """ + + training_dict={"Training Set": [int(self.sliders['number_of_frames'].value*(self.sliders['training_percentage'].value)),\ + round(100*(self.sliders['training_percentage'].value))], + "Testing Set": [int(self.sliders['number_of_frames'].value*(1-self.sliders['training_percentage'].value)),\ + round(100*(1-self.sliders['training_percentage'].value))]} + headers=["Partition", "Number of Frames", "Percentage"] + self.verbosity_wrap(markdown_table_from_dict(training_dict, headers)) + + self.train_idx, self.test_idx=split_dataset(self.frames[: self.sliders['number_of_frames'].value], + self.sliders['training_percentage'].value) + self.train_krr_model_func(self.train_idx, jitter= jitter) + + def plot_prediction_func(self, y_known= None, frames = None, + frame_idx= [], pretend = False): + """ Function to predict properties from the given SOAP representation and + KRR model. It is used for both the time estimation and the tutorial + prediction. The flag pretend is set to True when it is being used for + estimation, false otherwise. + + Author: R. Cersonsky + + Keywords: + y_known - list, target properties to train on + frames - frames to predict properties of + frame_idx - list of indices of the frames to train on + pretend - boolean, whether to not show the results + """ + if(len(frame_idx) > 0): + frames=self.frames[frame_idx] + y_known = np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[: , 0] if self.sliders['kernel_type'].value == 'Atom' else self.properties[self.sliders['property_to_ml'].value][frame_idx] + if(pretend == False): + verbosity_wrap=lambda s: self.verbosity_wrap(s) + else: + verbosity_wrap = lambda s: None + if(self.trained[self.sliders['property_to_ml'].value] == False): + verbosity_wrap("Model has not yet been trained, training now...") + self.train_krr_model() + # np.concatenate(self.properties[self.sliders['property_to_ml'].value][frame_idx])[:,0] if self.sliders['kernel_type'].value=='Atom' else self.properties[self.sliders['property_to_ml'].value][frame_idx] + testing_properties=y_known + + if(pretend == False): + self.estimate_time(x=self.pred_frames, y=self.pred_times, f=self.plot_prediction_func, N=max( + 0, 20-len(self.pred_frames)), ref=len(frames), p="Estimating time to compute prediction...") + est=int(np.poly1d(np.polyfit(self.pred_frames, + self.pred_times, deg=3))(len(frames)))+1 + else: + est=0 + + t=time.time() + verbosity_wrap("Predicting the properties of our data set will take approximately {} minutes and {} seconds.".format( + int(est/60), int(est % 60))) + y_pred=self.krr[self.sliders['property_to_ml'].value].predict(frames) + self.pred_frames.append(len(frame_idx)) + self.pred_times.append(time.time()-t) + verbosity_wrap(markdown_table_from_dict( + get_score(y_pred, testing_properties), headers = ["Statistic", "Value"])) + + if(not pretend): + plt.scatter(y_pred, testing_properties, s=3) + plt.axis('scaled') + plt.xlabel(self.sliders['property_to_ml'].value) + plt.ylabel('Predicted '+self.sliders['property_to_ml'].value) + plt.gca().set_aspect('equal') + plt.show() + + def predict_test_set(self): + """ Wrapper for plot_prediction_func that is used within the tutorial + + Author: R. Cersonsky + """ + self.plot_prediction_func(frame_idx=self.test_idx) + + def predict_new_set(self, filename='./data/small_molecules-1000.xyz', + num_frames=''): + """ Wrapper predicts a new set of data beyond the training set, loaded + from a different file + + Author: R. Cersonsky + + Keywords: + filename - string corresponding to an ASE-like file + num_frames - number of frames to predict + """ + frames=np.array(read(filename, ":{}".format(num_frames))) + properties=extract_property( + frames, self.sliders['property_to_ml'].value) + self.plot_prediction_func(y_known=properties, frames=frames) + + def output_params(self): + """ Function to output the hyperparameters that are currently set + + Author: R. Cersonsky + """ + self.verbosity_wrap('Our input file is {}, of which we are using {} \ + frames.'.format(self.input_file, \ + self.sliders['number_of_frames'].value)) + + hdict={hyper_vals[k]['name']: [self.hyperparameters[k]] + for k in self.hyperparameters + if 'fixed' not in hyper_vals[k]} + headers=["Parameter", "Value"] + + self.verbosity_wrap("
Our hyperparameters are {}".format( + markdown_table_from_dict(hdict, headers))) + + def set(self, value_name, value): + """ Public Setter function for tutorial interface + + Author: R. Cersonsky + + Keywords: + value_name - string, must be in self.sliders + value - value to set + """ + assert value_name in self.sliders + self.sliders[value_name].value=value + + def estimate_time(self, x=None, y=None, f=None, N=20, ref=None, p=None): + """ Helper function for estimating time + + Author: R. Cersonsky + + Keywords: + x - list of numbers of frames computed + y - list, size=size(x), of times to computed frames in x + f - function to estimate the run time of + N - number of trials to run + ref - int, upper limit of trials to run. + Defaults to (number_of_frames in the ultimate computation) / 4 + p - print statement + """ + if(N <= 0): + return + if(x == None or y == None or f == None or ref == None): + x=self.est_frames + y=self.est_times + f=self.train_krr_model_func + ref=int(0.25*self.sliders['number_of_frames'].value * + self.sliders['training_percentage'].value) + + if(p != None and min(N, ref) > 0): + self.verbosity_wrap(p) + for nf in tqdm(np.random.randint(2, ref, size = min(N, ref))): + f(frame_idx=np.random.randint(len(self.frames), size = nf), pretend=True) diff --git a/docs/source/tutorials/utilities/tutorial_utils.py b/docs/source/tutorials/utilities/tutorial_utils.py new file mode 100644 index 000000000..e69de29bb diff --git a/docs/source/tutorials/utilities/vectors_utils.py b/docs/source/tutorials/utilities/vectors_utils.py new file mode 100644 index 000000000..eff2088ea --- /dev/null +++ b/docs/source/tutorials/utilities/vectors_utils.py @@ -0,0 +1,429 @@ +import os +import inspect +import time + +from IPython.display import Markdown, display, clear_output, Image, HTML +import ipywidgets as widgets +import numpy as np +import ase +from ase.io import read +from matplotlib import pyplot as plt +from tqdm import tqdm_notebook as tqdm +import pandas as pd +import json + +# This is for backwards compatibility -- old versions of RASCAL will follow +# the latter schema, newer versions the former. +try: + from rascal.representations import SphericalInvariants as SOAP +except: + from rascal.representations import SOAP + +from general_utils import (markdown_table_from_dict, \ + readme_button, _button_template_, + compute_kernel, mask_body_order) + +# These are some imports/helpers for all of the tutorials +style = json.load(open('./utilities/widget_style.json')) +hyper_dict = json.load(open("./utilities/hyperparameter_presets.json")) +hyper_vals = json.load(open("./utilities/hyperparameter_bounds.json")) + + +def format_positions(positions): + """ Function to format positions for best viewing angle. Does a simple PCA + to find the best direction to view the molecule and rotates position + along these axes + + Author: R. Cersonsky + + Keywords: + positions - list of 3-vectors + + Returns: + positions - list of 3-vectors + """ + from sklearn.decomposition import PCA + pca = PCA(n_components=3) + return pca.fit_transform(positions) + + +class SOAP_tutorial(object): + """ Class to wrap the SOAP vectors tutorial + + Author: R. Cersonsky + + Constructor Keywords: + molecule_file - filename of ASE Atoms-like object + containing molecule + interactive - boolean, whether to load the tutorial interactive- + ly in the jupyter notebook + verbose - boolean, whether or not to narrate the actions of + the tutorial + hyperparameters - dictionary, hyperparameters for SOAP calculation + or string corresponding to one of the SOAP presets + + Functions: + predict - given the KRR, predicts the properties for a set of ASE frames + """ + + def __init__(self, + molecule_file=None, + interactive=False, + verbose=True, + hyperparameters=dict(**hyper_dict['Power Spectrum']), + ): + + # Check if hyperparameters are of type string or dictionary + if(isinstance(hyperparameters, str)): + self.hyperparameters = hyper_dict[hyperparameters] + else: + self.hyperparameters = {h: hyperparameters[h] + if h in hyperparameters + else hyper_dict['Power Spectrum'][h] + for h in hyper_dict['Power Spectrum']} + self.verbose = verbose + self.verbosity_wrap = lambda s: ( + None if not verbose else display(Markdown(s))) + + # Look for and load available molecules + file_options = ['./data/molecules/{}'.format(f) + for f in os.listdir('./data/molecules/') + if f.endswith('xyz')] + + if(molecule_file != None and molecule_file not in file_options): + file_options.append(molecule_file) + + self.frames, self.positions, self.symbols, self.inds, self.labels = {}, {}, {}, {}, {} + self.import_molecules(file_options) + + self.vectors = None + self.interactive = interactive + + # This next section is for the jupyter interface. Even if the user has + # specified zero interactivity, the sliders serve as containers for the + # hyperparameters to change when set + self.sliders = {val: + widgets.FloatSlider( + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + min=hyper_vals[val]['options'][0], + max=hyper_vals[val]['options'][1], + description=hyper_vals[val]['name'], + continuous_update=True, + step=(hyper_vals[val]['options'][1] - + hyper_vals[val]['options'][0])/20., + style=style, + ) + if isinstance(hyper_vals[val]['options'][0], float) else + widgets.IntSlider( + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + min=hyper_vals[val]['options'][0], + max=hyper_vals[val]['options'][1], + description=hyper_vals[val]['name'], + continuous_update=True, + step=1, + style=style, + ) + if isinstance(hyper_vals[val]['options'][0], int) and + hyper_vals[val]['options'][0] != True else + widgets.Dropdown(options=hyper_vals[val]['options'], + style=style, + value=self.hyperparameters.get(val, + hyper_vals[val]['options'][0]), + description=hyper_vals[val]['name'], + ) + for val in hyper_vals if 'fixed' not in hyper_vals[val]} + + # Which elements to center smooth gaussians on. Defaults to all elements + # present except for hydrogen. + constituents = list(set([s for k in self.symbols for s in self.symbols[k]])) + self.sliders['center_select']=widgets.SelectMultiple( + options=constituents, + value=[s for s in constituents if s != "H"], + description="Where to Place SOAP Centers", + style=style) + + # Whether to compute the global kernel (average of the environments) or + # atomic kernel. + self.sliders['average'] = widgets.Dropdown( + options=["Environment-Centered", "Average"], + value="Environment-Centered", + description="Type of SOAP Vectors", + style=style) + + # Show and enable the widgets if interactive==True + if(interactive): + self.preset_button = _button_template_(list(hyper_dict.keys()), + "SOAP Presets: ", + self.preset_func) + if(isinstance(hyperparameters, str)): + self.preset_button.value = hyperparameters + + slider_order = ['center_select', + 'average', + *list(hyper_vals.keys())] + + for s in slider_order: + if(s in self.sliders): + display(self.sliders[s]) + self.sliders[s].observe( + lambda change: self.change_func(change), names='value') + + # These buttons allow users to choose up to two molecules to compute + # SOAP vectors and compare. + self.molecule_buttons = [widgets.Dropdown(options=[*self.frames.keys()], # , "Other"], + style=style,\ + value=list( + self.frames.keys())[0], + description="\tMolecule: " + ) for i in range(2)] + + def import_molecules(self, file_list): + """ Function to import ASE-like files and include them in the dictionaries + of available molecules. + + Author: R. Cersonsky + + Keywords: + file_list - list of strings corresponding to ASE-like files + """ + for file in file_list: + for frame in read(file, ":"): + key = str(frame.symbols) + + if(key in self.frames): + print("Overwriting previous entry for {}".format(key)) + + self.frames[key] = frame + self.positions[key] = format_positions(frame.get_positions()) + self.inds[key] = sorted(range(len(self.positions[key])), + key=lambda j: self.positions[key][j][0]) + self.positions[key] = self.positions[key][self.inds[key]] + self.symbols[key] = frame.symbols[self.inds[key]] + super_scripts = ['' if list(self.symbols[key]).count(s) == 1 + else '^{{({})}}'.format(list(self.symbols[key])[:i].count(s)+1) + for i, s in enumerate(list(self.symbols[key]))] + self.labels[key] = np.array([r'${}{}$'.format(self.symbols[key][i], s) + for i, s in enumerate(super_scripts)]) + + def change_func(self, change): + """ Function to import ASE-like files and include them in the dictionaries + of available molecules. + + Author: R. Cersonsky + + Keywords: + file_list - list of strings corresponding to ASE-like files + """ + change['owner'].value = change['new'] + for s in self.hyperparameters: + if(s in self.sliders): + if(self.hyperparameters[s] != self.sliders[s].value): + self.hyperparameters[s] = self.sliders[s].value + self.vectors = None + + def preset_func(self, a): + """ Function to change hyperparameters based upon presets + + Author: R. Cersonsky + """ + for s in self.sliders: + if(s in self.hyperparameters): + self.hyperparameters[s] = hyper_dict[self.preset_button.value][s] + self.sliders[s].value = self.hyperparameters[s] + + def get_soap_vectors(self, key=None, average=False): + """ Get the SOAP vectors for the given molecule, denoted by key + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by \ + import_molecules + average - whether to compute the average SOAP vector for the + molecule + """ + if(key == None): + key = self.input_file.replace('./data/molecules/', '')[:-4] + + representation = SOAP(**mask_body_order(self.hyperparameters)) + all_frames = [self.frames[k] for k in sorted(self.frames.keys())] + features = representation.transform( + all_frames).get_dense_feature_matrix(representation) + + lengths = [len(self.inds[k]) for k in sorted(self.frames.keys())] + + self.vectors = {k: features[sum(lengths[0:i]):sum(lengths[0:i])+lengths[i]] + for i, k in enumerate(sorted(self.frames.keys()))} + vectors = self.vectors[key] + + # if(average): + # return np.mean(vectors, axis=0) + return vectors[self.inds[key]] + + def make_figure(self, key, img_name): + """ Function to generate an image of the given molecule and cache it + under img_name + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by + import_molecules + img_name - filename under which to cache the image + """ + from ase.data import covalent_radii, atomic_numbers + from matplotlib import pyplot, patches + from ase.data.colors import jmol_colors + style = { + "horizontalalignment": "center", + "verticalalignment": "center", + "fontsize": 12 + } + fig, ax = pyplot.subplots(1) + + for i, p in enumerate(self.positions[key]): + an = atomic_numbers[self.symbols[key][i]] + ax.text(*p[:2], + s=self.labels[key][i], + **style, + zorder=p[-1]+0.01 + ) + ax.add_artist(patches.Circle(p[:2], + covalent_radii[an], + facecolor=jmol_colors[an], + edgecolor='k', + alpha=0.95, + zorder=p[-1])) + ax.axis('off') + bounds = [*np.min(self.positions[key], axis=0), * + np.max(self.positions[key], axis=0)] + ax.set_xlim([bounds[0]-1.5, bounds[3]+1.5]) + ax.set_ylim([bounds[1]-1.5, bounds[5]+1.5]) + ax.set_aspect('equal') + fig.savefig(img_name) + fig.clf() + + def show_kernel(self, key1, key2=None, show=True, average=False): + from rascal.models import Kernel + from ase.data import chemical_symbols + """ Function to show pandas dataframe for the comparison kernel + + Author: R. Cersonsky + + Keywords: + key1 - string corresponding to molecule built by + import_molecules + key2 - string corresponding to molecule built by + import_molecules, defaults to key1 + show - boolean, whether to show or return the dataframe + average - boolean, whether to compute the average SOAP vector + for the molecule + """ + key2 = key2 if key2 != None else key1 + average = self.sliders['average'].value == "Average" + + soap=SOAP(**mask_body_order(self.hyperparameters)) + features1=soap.transform(self.frames[key1]) + features2=soap.transform(self.frames[key2]) + + kernel = Kernel(soap, + target_type="Structure" if average else "Atom", + zeta = 2, + **self.hyperparameters) + data = kernel(features2, features1) + + if (average): + + vec1 = np.mean(features1.get_dense_feature_matrix(soap),axis=0) + vec2 = np.mean(features2.get_dense_feature_matrix(soap), axis=0) + if(vec1.shape==vec2.shape): + data = np.dot(vec1, vec2)/(np.linalg.norm(vec1)*np.linalg.norm(vec2)) + df=pd.DataFrame(data=[['%1.3f' %(data)]], + columns =[key2], + index =[key1] + ) + else: + df=pd.DataFrame(data=[['%1.3f' %(data[i][j]) + for i,v in enumerate(features2._frames.numbers) if v!=1] + for j,w in enumerate(features1._frames.numbers) if w!=1], + columns =[chemical_symbols[v] for i,v in enumerate(features2._frames.numbers) if v!=1], + index =[chemical_symbols[v] for i,v in enumerate(features1._frames.numbers) if v!=1] + ) + # + df.name = ("Self-" if key1 == key2 else '') + \ + "Similarity " + ("Kernel" if not average else "") + return df + + def show_molecule(self, key, show=True): + """ Function to either generate or show the image for the given molecule + + Author: R. Cersonsky + + Keywords: + key - string corresponding to molecule built by + import_molecules + show - boolean, whether to show or return the dataframe + """ + import os + img_name = './images/cache/{}.png'.format(key) + if(not os.path.exists(img_name)): + self.make_figure(key, img_name) + if(show): + display(Image(img_name)) + return img_name + + def compute(self, a=None): + """ Function to compute the SOAP vectors for the set hyperparameters + + Author: R. Cersonsky + """ + def wrap_compute(a): + self.get_soap_vectors( + key=self.molecule_buttons[0].value, average=self.sliders['average'].value == "Average") + print("SOAP Vectors for {}: \n\t".format( + self.molecule_buttons[0].value), self.vectors[self.molecule_buttons[0].value]) + + compute_button = widgets.Button(description="Compute", + style=style, + value=False, + ) + compute_button.on_click(wrap_compute) + clear_button = widgets.Button(description="Clear Output", + style=style, + value=False, + ) + clear_button.on_click(lambda a: self.compute(a)) + display(self.molecule_buttons[0]) + display(widgets.widgets.HBox((compute_button, clear_button))) + + def compare(self, a=None): + """ Function to compare the SOAP vectors of two molecules + for the set hyperparameters + + Author: R. Cersonsky + """ + clear_output() + + def wrap_compare(a): + df = self.show_kernel( + key1=self.molecule_buttons[0].value, key2=self.molecule_buttons[1].value, show=False) + m1 = self.show_molecule(self.molecule_buttons[0].value, show=False) + m2 = self.show_molecule(self.molecule_buttons[1].value, show=False) + display(HTML("
| {} | {} | {} |
|---|---|---|
| {} |