ex1_genotype_lattice_signals.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import numpy as np\n",
    "import lsp \n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HIV Genotype Lattice Signals\n",
    "\n",
    "The signals s1 and s2 are taken from [1, Appendix Table 4 and Table 7]. s1 maps HIV RT genotypes to their RT inhibitor zidovudine resistance, and, s2 maps HIV PR genotypes to their PR inhibitor indinavir resistance.\n",
    "\n",
    "[1] Beerenwinkel, N., Knupfer, P. and Tresch, A.; Learning Monotonic Genotype-Phenotype Maps; https://www.research-collection.ethz.ch/handle/20.500.11850/30573"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = {}\n",
    "s1[(0,0,0,0,0,0,0)] = -0.06\n",
    "s1[(0,1,0,0,0,0,0)] = 0.13\n",
    "s1[(0,1,1,0,0,0,0)] = 0.13\n",
    "s1[(1,0,0,0,0,0,0)] = 0.15\n",
    "s1[(1,0,0,0,0,1,0)] = 0.80\n",
    "s1[(0,1,0,1,0,0,0)] = 1\n",
    "s1[(0,1,0,1,0,0,1)] = 1\n",
    "s1[(0,1,1,1,0,0,0)] = 1\n",
    "s1[(1,0,0,0,1,1,0)] = 1.23\n",
    "s1[(1,1,0,0,0,0,0)] = 1.3\n",
    "s1[(1,1,0,0,0,1,0)] = 1.3\n",
    "s1[(1,1,1,0,0,0,0)] = 1.3\n",
    "s1[(1,1,1,0,0,1,0)] = 1.3\n",
    "s1[(1,1,0,0,1,1,0)] = 1.67\n",
    "s1[(1,1,0,1,0,0,0)] = 1.67\n",
    "s1[(1,1,0,1,0,1,0)] = 1.67\n",
    "s1[(1,1,0,1,1,1,0)] = 1.67\n",
    "s1[(0,1,1,1,0,0,1)] = 1.83\n",
    "s1[(1,1,0,1,0,0,1)] = 2.06\n",
    "s1[(1,1,0,1,0,1,1)] = 2.06\n",
    "s1[(1,1,0,1,1,1,1)] = 2.06\n",
    "s1[(1,1,1,0,1,1,0)] = 2.24\n",
    "s1[(1,1,1,1,0,0,0)] = 2.24\n",
    "s1[(1,1,1,1,0,0,1)] = 2.24\n",
    "s1[(1,1,1,1,0,1,0)] = 2.24\n",
    "s1[(1,1,1,1,0,1,1)] = 2.24\n",
    "s1[(1,1,1,1,1,1,0)] = 2.24\n",
    "s1[(1,1,1,1,1,1,1)] = 2.24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = {}\n",
    "s2[(0,0,0,0,0,0)] = -0.05\n",
    "s2[(0,0,0,0,0,1)] = 0.59\n",
    "s2[(0,0,1,0,0,0)] = 1.07\n",
    "s2[(0,0,1,1,0,0)] = 1.07 \n",
    "s2[(1,0,0,0,0,0)] = 1.08\n",
    "s2[(1,0,1,0,0,0)] = 1.16\n",
    "s2[(1,0,1,1,0,0)] = 1.16\n",
    "s2[(1,0,0,0,0,1)] = 1.20 \n",
    "s2[(0,0,0,0,1,1)] = 1.23\n",
    "s2[(1,0,0,0,1,1)] = 1.35\n",
    "s2[(0,0,1,0,0,1)] = 1.40\n",
    "s2[(0,0,1,0,1,1)] = 1.40\n",
    "s2[(0,0,1,1,0,1)] = 1.40\n",
    "s2[(0,0,1,1,1,1)] = 1.40\n",
    "s2[(1,0,1,0,0,1)] = 1.53\n",
    "s2[(1,0,1,0,1,1)] = 1.53\n",
    "s2[(1,0,1,1,0,1)] = 1.53\n",
    "s2[(1,0,1,1,1,1)] = 1.53\n",
    "s2[(0,1,0,0,0,0)] = 1.99 \n",
    "s2[(0,1,0,0,0,1)] = 1.99\n",
    "s2[(0,1,0,0,1,1)] = 1.99\n",
    "s2[(0,1,1,0,0,0)] = 1.99\n",
    "s2[(0,1,1,0,0,1)] = 1.99\n",
    "s2[(0,1,1,0,1,1)] = 1.99\n",
    "s2[(0,1,1,1,0,0)] = 1.99\n",
    "s2[(0,1,1,1,0,1)] = 1.99\n",
    "s2[(0,1,1,1,1,1)] = 1.99\n",
    "s2[(1,1,0,0,0,0)] = 2.38\n",
    "s2[(1,1,0,0,0,1)] = 2.38\n",
    "s2[(1,1,0,0,1,1)] = 2.38\n",
    "s2[(1,1,1,0,0,0)] = 2.38\n",
    "s2[(1,1,1,1,0,0)] = 2.38\n",
    "s2[(1,1,1,0,0,1)] = 2.61\n",
    "s2[(1,1,1,0,1,1)] = 2.61\n",
    "s2[(1,1,1,1,0,1)] = 2.61\n",
    "s2[(1,1,1,1,1,1)] = 2.61"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def meet(a, b):\n",
    "    c = np.asarray(a)*np.asarray(b)\n",
    "    c = c.astype(np.int32)\n",
    "    c = tuple(c.tolist())\n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "closing under meet...\n",
      "computing partial order...\n",
      "inverting partial order matrix...\n",
      "closing under meet...\n",
      "computing partial order...\n",
      "inverting partial order matrix...\n"
     ]
    }
   ],
   "source": [
    "sp1 = lsp.LSP(list(s1.keys()), meet)\n",
    "sp2 = lsp.LSP(list(s2.keys()), meet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rt = np.asarray([s1[tuple(l)] for l in sp1.lattice])\n",
    "pr = np.asarray([s2[tuple(l)] for l in sp2.lattice])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 28 artists>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABq8AAAI/CAYAAADp3kNcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdX4il933f8c+3Gtm9aKnTaFqrkrbjUl3UpaFxBsUhN4bEIGmL1DYJSBfxH2oWikVa6M24ARt8Nb1pwdjYbBphuRTbxW2aDbtFxPmDW4qN1kZRLAmRrdiiZUW9sYrS4DRG7bcXmpbpanZnJD075zszrxcMe57z/Pw7X+nsrDy893lOdXcAAAAAAABggj+z6gEAAAAAAADg/xKvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhjbVUvfMcdd/TGxsaqXh4AAFjYt7/97T/s7vVVz8Hx5GdIAAA4Xm72M+TK4tXGxkYuXry4qpcHAAAWVlX/ddUzcHz5GRIAAI6Xm/0M6baBAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAADAoqrq8ar6XlV99wbnq6o+U1WXquqZqnrfYc8IAADMJV4BAACwtC8muf8m5x9Icu/O15kknz+EmQAAgCNCvAIAAGBR3f2NJK/cZMnDSb7Ur/tmkndV1Z2HMx0AADCdeAUAAMBhuyvJS7uOr+w8BwAAIF4BAABw6GqP5/oNi6rOVNXFqrp47dq1QxgLAACYQLwCAADgsF1Jcs+u47uTXL1+UXef7e7N7t5cX18/tOEAAIDVEq8AAAA4bOeSfKhe9/4kr3b3y6seCgAAmGFt1QMAAABwvFTVl5N8IMkdVXUlyaeS3J4k3f2FJBeSPJjkUpIfJPnoaiYFAAAmEq8AAABYVHc/us/5TvLxQxoHAAA4Ytw2EAAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGCMtVUPAJxsG1vnF93v8vbpRfc7SbwXJ89xeM+Pwz8Db85hvOfH5TUATgJ/ngIAHE+uvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxlhb9QAAAADA4dvYOr/ofpe3Ty+6HwAAJ5crrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxlhb9QAAAAAAwHI2ts4vut/l7dOL7gcA+3HlFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY6ytegCAo25j6/yi+13ePr3ofrBKS35/+N6Yz5+HAAAAwBJceQUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjLFvvKqqe6rqd6rq+ap6tqr+0R5rqqo+U1WXquqZqnrfrRkXAAAAAACA42ztAGteS/JPuvs7VfXnk3y7qn6zu5/bteaBJPfufP1kks/v/AoAAAAAAAAHtu+VV939cnd/Z+fx/0jyfJK7rlv2cJIv9eu+meRdVXXn4tMCAAAAAABwrL2pz7yqqo0kP57kW9eduivJS7uOr+SNgQsAAAAAAABu6sDxqqr+XJJ/m+Qfd/cfXX96j/9J77HHmaq6WFUXr1279uYmBQAAAAAA4Ng7ULyqqtvzerj619397/ZYciXJPbuO705y9fpF3X22uze7e3N9ff2tzAsAAAAAAMAxtm+8qqpK8qtJnu/uf36DZeeSfKhe9/4kr3b3ywvOCQAAAAAAwAmwdoA1P53kF5P8flU9vfPcP01yKkm6+wtJLiR5MMmlJD9I8tHlRwUAAAAAAOC42zdedfd/yt6fabV7TSf5+FJDAQAAAAAAcDId6DOvAAAAAAAA4DCIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIyxtuoBAAAAAACut7F1ftH9Lm+fXnQ/AG4dV14BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMMbaqgcAYH8bW+cX3e/y9ulF9wMAADgMfjYC4M3w342jy5VXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAi6qq+6vqhaq6VFVbe5z/SFVdq6qnd74+too5AQCAmdZWPQAAAADHR1XdluRzST6Y5EqSp6rqXHc/d93Sr3b3Y4c+IAAAMJ4rrwAAAFjSfUkudfeL3f3DJF9J8vCKZwIAAI4Q8QoAAIAl3ZXkpV3HV3aeu97PVdUzVfW1qrrncEYDAACOAvEKAACAJdUez/V1x7+RZKO7fyzJ15M8sedGVWeq6mJVXbx27drCYwIAAFOJVwAAACzpSpLdV1LdneTq7gXd/f3u/tOdw19J8hN7bdTdZ7t7s7s319fXb8mwAADAPOIVAAAAS3oqyb1V9Z6qekeSR5Kc272gqu7cdfhQkucPcT4AAGC4tVUPAAAAwPHR3a9V1WNJnkxyW5LHu/vZqvp0kovdfS7JL1XVQ0leS/JKko+sbGAAAGAc8QoAAIBFdfeFJBeue+6Tux5/IsknDnsuAADgaHDbQAAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYIy1VQ8AAAAAAADHwcbW+UX3u7x9etH94Khw5RUAAAAAAABjuPKKE8fffgCOAn9WnSyH8X77PQUAAAAcFa68AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAx1lY9AAAAAACcJBtb5xfb6/L26cX2AoApXHkFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGPvGq6p6vKq+V1XfvcH5D1TVq1X19M7XJ5cfEwAAAAAAgJNg7QBrvpjks0m+dJM1/7G7/84iEwEAAAAAAHBi7XvlVXd/I8krhzALAAAAAAAAJ9xSn3n1U1X1e1X1H6rqby60JwAAAAAAACfMQW4buJ/vJPmr3f3HVfVgkn+f5N69FlbVmSRnkuTUqVMLvDQAAAAAAADHydu+8qq7/6i7/3jn8YUkt1fVHTdYe7a7N7t7c319/e2+NAAAAAAAAMfM245XVfXuqqqdx/ft7Pn9t7svAAAAAAAAJ8++tw2sqi8n+UCSO6rqSpJPJbk9Sbr7C0l+Psk/rKrXkvxJkke6u2/ZxAAAAAAAABxb+8ar7n50n/OfTfLZxSYCAAAAAADgxHrbtw0EAAAAAACApYhXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIyxtuoBAAAAgONpY+v8Yntd3j692F4AAMzmyisAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDJ95BQAAAMMs+VlRic+LAgDgaHHlFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADDG2qoHgN02ts4vut/l7dOL7nfSeD84ig7j963vDQAAAAC4dVx5BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAsKiqur+qXqiqS1W1tcf5d1bVV3fOf6uqNg5/SgAAYCrxCgAAgMVU1W1JPpfkgSTvTfJoVb33umX/IMl/7+6/nuRfJPlnhzslAAAwmXgFAADAku5Lcqm7X+zuHyb5SpKHr1vzcJIndh5/LcnPVFUd4owAAMBga6seAABgso2t84vud3n79KL7AQx0V5KXdh1fSfKTN1rT3a9V1atJfjTJHx7KhAAAwGjiFQAAAEva6wqqfgtrUlVnkpxJklOnTr39yRZyGH+x4TD+ssNxeI3D+GfwF1lmOS7//nxvvNGq/iz0Pc7Sjsvv2+PyvXFc3o+TyG0DAQAAWNKVJPfsOr47ydUbramqtSR/Ickr12/U3We7e7O7N9fX12/RuAAAwDTiFQAAAEt6Ksm9VfWeqnpHkkeSnLtuzbkkH955/PNJfru733DlFQAAcDK5bSAAAACL2fkMq8eSPJnktiSPd/ezVfXpJBe7+1ySX03yr6rqUl6/4uqR1U0MAABMI14BAACwqO6+kOTCdc99ctfj/5nkFw57LgAA4Ghw20AAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHWVj0AAAAAHCWXt0+vegQA4ATz/0U4CVx5BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABj7BuvqurxqvpeVX33Buerqj5TVZeq6pmqet/yYwIAAAAAAHASHOTKqy8muf8m5x9Icu/O15kkn3/7YwEAAAAAAHAS7RuvuvsbSV65yZKHk3ypX/fNJO+qqjuXGhAAAAAAAICTY4nPvLoryUu7jq/sPAcAAAAAAABvytoCe9Qez/WeC6vO5PVbC+bUqVMLvDS7bWydX2yvy9unF9uLW2PJ9zs53u+5f1cAAAAAAEfHEldeXUlyz67ju5Nc3Wthd5/t7s3u3lxfX1/gpQEAAAAAADhOlohX55J8qF73/iSvdvfLC+wLAAAAAADACbPvbQOr6stJPpDkjqq6kuRTSW5Pku7+QpILSR5McinJD5J89FYNCwAAAAAAwPG2b7zq7kf3Od9JPr7YRAAAAAAAAJxYS9w2EAAAAAAAABYhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwxtqqBwAAAAAAAOa4vH161SNwwrnyCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYY23VAwBwcmxsnV90v8vbpxfdDwAAAABYPVdeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAMAiquovVtVvVtUf7Pz6IzdY97+q6umdr3OHPScAADCbeAUAAMBStpL8Vnffm+S3do738ifd/bd3vh46vPEAAICjQLwCAABgKQ8neWLn8RNJ/u4KZwEAAI4o8QoAAICl/OXufjlJdn79SzdY92er6mJVfbOqBC4AAOD/s7bqAQAAADg6qurrSd69x6lffhPbnOruq1X115L8dlX9fnf/lz1e60ySM0ly6tSptzQvvF2Xt0+vegQAgBNHvAIAAODAuvtnb3Suqv5bVd3Z3S9X1Z1JvneDPa7u/PpiVf1ukh9P8oZ41d1nk5xNks3NzV5gfAAA4Ahw20AAAACWci7Jh3cefzjJr1+/oKp+pKreufP4jiQ/neS5Q5sQAAAYT7wCAABgKdtJPlhVf5DkgzvHqarNqvqXO2v+RpKLVfV7SX4nyXZ3i1cAAMD/47aBAAAALKK7v5/kZ/Z4/mKSj+08/s9J/tYhjwYAABwhrrwCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGWFv1AAAAAAAAAEfR5e3Tqx7hWHLlFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjrK16AAAAAACAVbi8fXrVIwCwB1deAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIxxoHhVVfdX1QtVdamqtvY4/5GqulZVT+98fWz5UQEAAAAAADju1vZbUFW3Jflckg8muZLkqao6193PXbf0q9392C2YEQAAAAAAgBPiIFde3ZfkUne/2N0/TPKVJA/f2rEAAAAAAAA4iQ4Sr+5K8tKu4ys7z13v56rqmar6WlXds8h0AAAAAAAAnCgHiVe1x3N93fFvJNno7h9L8vUkT+y5UdWZqrpYVRevXbv25iYFAAAAAADg2DtIvLqSZPeVVHcnubp7QXd/v7v/dOfwV5L8xF4bdffZ7t7s7s319fW3Mi8AAAAAAADH2EHi1VNJ7q2q91TVO5I8kuTc7gVVdeeuw4eSPL/ciAAAAAAAAJwUa/st6O7XquqxJE8muS3J4939bFV9OsnF7j6X5Jeq6qEkryV5JclHbuHMAAAAAAAAHFP7xqsk6e4LSS5c99wndz3+RJJPLDsaAAAAAAAAJ81BbhsIAAAAAAAAh0K8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAYBFV9QtV9WxV/e+q2rzJuvur6oWqulRVW4c5IwAAMJ94BQAAwFK+m+TvJ/nGjRZU1W1JPpfkgSTvTfJoVb33cMYDAACOgrVVDwAAAMDx0N3PJ0lV3WzZfUkudfeLO2u/kuThJM/d8gEBAIAjwZVXAAAAHKa7kry06/jKznMAAABJXHkFAADAm1BVX0/y7j1O/XJ3//pBttjjub7Ba51JciZJTp06deAZAQCAo028AgAA4MC6+2ff5hZXktyz6/juJFdv8Fpnk5xNks3NzT0DFwAAcPy4bSAAAACH6akk91bVe6rqHUkeSXJuxTMBAACDiFcAAAAsoqr+XlVdSfJTSc5X1ZM7z/+VqrqQJN39WpLHkjyZ5Pkk/6a7n13VzAAAwDxuGwgAAMAiuvvXkvzaHs9fTfLgruMLSS4c4mgAAMARIl5xYBtb5xfd7/L26UX3AwAAAAAAjj63DQQAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADG8JlXh8TnRZ0s3m8AAAAAAHhrXHkFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMcaB4VVX3V9ULVXWpqrb2OP/Oqvrqzvlv1f9p725CPavLOIB/H2ZyYxGWU4g6aeVGWvRycVOEEL2oiynI0EUYFLbIqF3SoiQILHrbRGAoWFQSvQ4kVFBQq3AUSU2sSWZyUjQzKhch5tPi/oXLZWbupPfec+75fT6be86Zw/8+wzO/e8/z/845/6qLtrtQAAAAAAAAlm/L8Kqq9iX5epIrklya5NqqunTTaR9O8o/ufn2Sryb5wnYXCgAAAAAAwPKdyZ1XlyU52t0Pd/czSe5IcmjTOYeS3L7a/kGSd1RVbV+ZAAAAAAAAjOBMwqvzkzyyYf/E6thJz+nuZ5P8M8krt6NAAAAAAAAAxlHdffoTqq5O8u7u/shq/4NJLuvuj28454HVOSdW+39enfP3Ta91fZLrk+TgwYNvOX78+Hb+XV6wi2782ba+3rGbr9rW1wMAgL2gqu7u7rWp62CZ1tbW+siRI1OXAcCK99MAeLFON0OeyZ1XJ5JcuGH/giSPnuqcqtqf5OVJntr8Qt19S3evdffagQMHzqR2AAAAAAAABnIm4dVdSS6pqour6qwk1yQ5vOmcw0muW22/P8mveqtbugAAAAAAAGCT/Vud0N3PVtUNSX6eZF+S27r7gar6XJIj3X04ya1Jvl1VR7N+x9U1O1k0AAAAAAAAy7RleJUk3X1nkjs3HfvMhu3/JLl6e0sDAAAAAABgNGfy2EAAAAAAAADYFcIrAAAAAAAAZkN4BQAAAAAAwGwIrwAAAAAAAJgN4RUAAAAAAACzIbwCAAAAAABgNoRXAAAAAAAAzIbwCgAAAAAAgNnYP3UBc3Ds5qumLgEAAAAA9gzvpwGwk9x5BQAAAAAAwGwIrwAAAAAAAJgN4RUAAAAAAACzIbwCAAAAAABgNv+4W04AAAU5SURBVIRXAAAAAAAAzIbwCgAAAAAAgNkQXgEAAAAAADAbwisAAAAAAABmQ3gFAAAAAADAbAivAAAAAAAAmA3hFQAAAAAAALMhvAIAAAAAAGA2hFcAAAAAAADMhvAKAAAAAACA2RBeAQAAAAAAMBvCKwAAAAAAAGZDeAUAAAAAAMBsCK8AAAAAAACYDeEVAAAAAAAAsyG8AgAAAAAAYDaEVwAAAAAAAMyG8AoAAAAAAIDZEF4BAAAAAAAwG8IrAAAAAAAAZqO6e5pvXPW3JMcn+eYv3LlJnpy6CHaNfo9Hz8ei32PR7/Ho+TRe090Hpi6CZTJDsgfo93j0fCz6PRb9Ho+eT+OUM+Rk4dVeVFVHuntt6jrYHfo9Hj0fi36PRb/Ho+fAHPhZNBb9Ho+ej0W/x6Lf49Hz+fHYQAAAAAAAAGZDeAUAAAAAAMBsCK/+P7dMXQC7Sr/Ho+dj0e+x6Pd49ByYAz+LxqLf49Hzsej3WPR7PHo+Mz7zCgAAAAAAgNlw5xUAAAAAAACzIbw6A1X1nqp6qKqOVtWNU9fDzquqY1V1X1XdW1VHpq6H7VVVt1XVE1V1/4Zjr6iqX1bVn1Zfz5myRrbXKXp+U1X9dbXO762qK6eske1TVRdW1a+r6sGqeqCqPrE6bp0v0Gn6bY0DkzFDjscMuWxmyPGYIcdihhyLGXLv8NjALVTVviR/TPLOJCeS3JXk2u7+w6SFsaOq6liSte5+cupa2H5V9fYkTyf5Vne/YXXsi0me6u6bV28wnNPdn5qyTrbPKXp+U5Knu/tLU9bG9quq85Kc1933VNXLktyd5L1JPhTrfHFO0+8PxBoHJmCGHJMZctnMkOMxQ47FDDkWM+Te4c6rrV2W5Gh3P9zdzyS5I8mhiWsCXoTu/k2SpzYdPpTk9tX27Vn/pcVCnKLnLFR3P9bd96y2/53kwSTnxzpfpNP0G2AqZkhYGDPkeMyQYzFDjsUMuXcIr7Z2fpJHNuyfiH/MI+gkv6iqu6vq+qmLYVe8ursfS9Z/iSV51cT1sDtuqKrfrx4J4fb/Baqqi5K8KcnvYp0v3qZ+J9Y4MA0z5JjMkONxbTkm15cLZ4Ycixly3oRXW6uTHPOsxeV7a3e/OckVST62ul0cWJZvJHldkjcmeSzJl6cth+1WVS9N8sMkn+zuf01dDzvrJP22xoGpmCHHZIaE5XN9uXBmyLGYIedPeLW1E0ku3LB/QZJHJ6qFXdLdj66+PpHkx1l/9AfL9vjqmbfPP/v2iYnrYYd19+Pd/d/ufi7JN2OdL0pVvSTrF6Hf6e4frQ5b5wt1sn5b48CEzJADMkMOybXlYFxfLpsZcixmyL1BeLW1u5JcUlUXV9VZSa5JcnjimthBVXX26sP6UlVnJ3lXkvunrYpdcDjJdavt65L8dMJa2AXPX4CuvC/W+WJUVSW5NcmD3f2VDX9knS/QqfptjQMTMkMOxgw5LNeWg3F9uVxmyLGYIfeO6vb0gq1U1ZVJvpZkX5LbuvvzE5fEDqqq12b9f8olyf4k39XzZamq7yW5PMm5SR5P8tkkP0ny/SQHk/wlydXd7cNZF+IUPb8867eCd5JjST76/LOs2duq6m1JfpvkviTPrQ5/OuvPsLbOF+Y0/b421jgwETPkWMyQy2eGHI8ZcixmyLGYIfcO4RUAAAAAAACz4bGBAAAAAAAAzIbwCgAAAAAAgNkQXgEAAAAAADAbwisAAAAAAABmQ3gFAAAAAADAbAivAAAAAAAAmA3hFQAAAAAAALMhvAIAAAAAAGA2/gesGN8RN+0jtQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 2160x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(30, 10))\n",
    "plt.subplot(1, 2, 1).bar(np.arange(len(rt)), rt)\n",
    "plt.subplot(1, 2, 2).bar(np.arange(len(rt)), sp1.ft(rt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11 out of 28 Fourier coefficients are zero.\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(np.abs(sp1.ft(rt)) <= 1e-8), \"out of\", len(sp1.lattice), \"Fourier coefficients are zero.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 36 artists>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABq8AAAI/CAYAAADp3kNcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdX4il933f8c83HiW9cKjTaloLSZtxqS6a0FC7gxzjG0EasLXGuogD8oX/EbM0WMQGX3Ttgg2+2ty4YGSsbrCwXYKdkgR3w24Izj9sU2y8EopsS4RuzbZaJOqtVGQLuzVqv73YKR2PZjWj3WfnfHfm9YJhz3Oe33nme3aMzPDe33OquwMAAAAAAAAT/MyqBwAAAAAAAID/R7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYIy1VX3jW2+9tTc2Nlb17QEAgIU98sgj/72711c9B4eT3yEBAOBwebnfIVcWrzY2NnL+/PlVfXsAAGBhVfVfVj0Dh5ffIQEA4HB5ud8h3TYQAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYY23VAwAAsFobJ89e1+svnjq+0CQvdVhn2znX5PcJ3FyW/G8TAACsip1XAAAAAAAAjCFeAQAAsJiqurOq/qqqnqyq71bVB3dZU1X1qaq6UFWPV9UbVjErAAAwk9sGAgAAsKQXk3y4ux+tqp9P8khVfaW7n9i25q1J7tr6emOSz2z9CQAAYOcVAAAAy+nuZ7r70a3HP0zyZJLbdyy7L8kX+opvJHlNVd12wKMCAABDiVcAAADcEFW1keT1Sb6549TtSZ7adnwpLw1cAADAESVeAQAAsLiqenWSP0ryoe7+wc7Tu7ykd7nGiao6X1XnL1++fCPGBAAABhKvAAAAWFRV3ZIr4er3u/uPd1lyKcmd247vSPL0zkXdfbq7N7t7c319/cYMCwAAjCNeAQAAsJiqqiSfTfJkd3/yKsvOJHl3XfGrSZ7v7mcObEgAAGC0tVUPAAAAwKHy5iTvSvLtqnps67mPJjmWJN39UJJzSe5NciHJj5K8bwVzAgAAQ4lXAAAALKa7v57dP9Nq+5pO8oGDmQgAALjZuG0gAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGPsGa+q6s6q+quqerKqvltVH9xlzT1V9XxVPbb19bEbMy4AAAAAAACH2do+1ryY5MPd/WhV/XySR6rqK939xI51X+vuty0/IgAAAAAAAEfFnjuvuvuZ7n506/EPkzyZ5PYbPRgAAAAAAABHzyv6zKuq2kjy+iTf3OX0m6rqb6rqT6vqlxeYDQAAAAAAgCNmP7cNTJJU1auT/FGSD3X3D3acfjTJL3b3C1V1b5IvJ7lrl2ucSHIiSY4dO3bNQwMAAAAAAHA47WvnVVXdkivh6ve7+493nu/uH3T3C1uPzyW5papu3WXd6e7e7O7N9fX16xwdAAAAAACAw2bPeFVVleSzSZ7s7k9eZc1rt9alqu7euu6zSw4KAAAAAADA4bef2wa+Ocm7kny7qh7beu6jSY4lSXc/lOQdSX67ql5M8uMk93d334B5AQAAAAAAOMT2jFfd/fUktceaB5M8uNRQAAAAAAAAHE37+swrAAAAAAAAOAjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjrK16AACAqTZOnr2u1188dXyx6y15rd2ut6Ql3+eSJv+dAQAAAP+fnVcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAspqoerqrvV9V3rnL+nqp6vqoe2/r62EHPCAAAzLa26gEAAAA4VD6X5MEkX3iZNV/r7rcdzDgAAMDNxs4rAAAAFtPdX03y3KrnAAAAbl7iFQAAAAftTVX1N1X1p1X1y6seBgAAmMVtAwEAADhIjyb5xe5+oaruTfLlJHfttrCqTiQ5kSTHjh07uAkBAICVsvMKAACAA9PdP+juF7Yen0tyS1XdepW1p7t7s7s319fXD3ROAABgdcQrAAAADkxVvbaqauvx3bnye+mzq50KAACYxG0DAQAAWExVfTHJPUlurapLST6e5JYk6e6HkrwjyW9X1YtJfpzk/u7uFY0LAAAMJF4BAACwmO5+5x7nH0zy4AGNAwAA3ITcNhAAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYIw941VV3VlVf1VVT1bVd6vqg7usqar6VFVdqKrHq+oNN2ZcAAAAAAAADrO1fax5McmHu/vRqvr5JI9U1Ve6+4lta96a5K6trzcm+czWnwAAAAAAALBve+686u5nuvvRrcc/TPJkktt3LLsvyRf6im8keU1V3bb4tAAAAAAAABxqr+gzr6pqI8nrk3xzx6nbkzy17fhSXhq4AAAAAAAA4GXtO15V1auT/FGSD3X3D3ae3uUlvcs1TlTV+ao6f/ny5Vc2KQAAAAAAAIfevuJVVd2SK+Hq97v7j3dZcinJnduO70jy9M5F3X26uze7e3N9ff1a5gUAAAAAAOAQW9trQVVVks8mebK7P3mVZWeSPFBVX0ryxiTPd/czy415NG2cPHtdr7946vhi19t5rSUt/T6n8vPcHz/P67/Wko7KzxMAAAAAmGPPeJXkzUneleTbVfXY1nMfTXIsSbr7oSTnktyb5EKSHyV53/KjAgAAAAAAcNjtGa+6++vZ/TOttq/pJB9YaigAAAAAAACOpn195hUAAAAAAAAcBPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYY23VAwAAAAAA7LRx8uw1v/biqeMLTgLAQbPzCgAAAAAAgDHsvAIAAIAVuZ5dBYmdBQAAHE52XgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOsrXoAAAAAAADgp22cPHvNr7146viCk8DBs/MKAAAAAACAMcQrAAAAAAAAxhCvAAAAWFRVPVxV36+q71zlfFXVp6rqQlU9XlVvOOgZAQCAucQrAAAAlva5JG95mfNvTXLX1teJJJ85gJkAAICbhHgFAADAorr7q0mee5kl9yX5Ql/xjSSvqarbDmY6AABgOvEKAACAg3Z7kqe2HV/aeg4AAEC8AgAA4MDVLs/1SxZVnaiq81V1/vLlywcwFgAAMIF4BQAAwEG7lOTObcd3JHl656LuPt3dm929ub6+fmDDAQAAqyVeAQAAcNDOJHl3XfGrSZ7v7mdWPRQAADDD2qoHAAAA4HCpqi8muSfJrVV1KcnHk9ySJN39UJJzSe5NciHJj5K8bzWTAgAAE4lXAAAALKq737nH+U7ygQMaBwAAuMm4bSAAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOsrXoAAODqNk6evebXXjx1fLFr7Xa9JR2V9wkAAADA3vbceVVVD1fV96vqO1c5f09VPV9Vj219fWz5MQEAAAAAADgK9rPz6nNJHkzyhZdZ87XuftsiEwEAAAAAAHBk7bnzqru/muS5A5gFAAAAAACAI27PeLVPb6qqv6mqP62qX17omgAAAAAAABwx+7lt4F4eTfKL3f1CVd2b5MtJ7tptYVWdSHIiSY4dO7bAtwYAAAAAAOAwue6dV939g+5+YevxuSS3VNWtV1l7urs3u3tzfX39er81AAAAAAAAh8x1x6uqem1V1dbju7eu+ez1XhcAAAAAAICjZ8/bBlbVF5Pck+TWqrqU5ONJbkmS7n4oyTuS/HZVvZjkx0nu7+6+YRMDAAAAAABwaO0Zr7r7nXucfzDJg4tNBAAAAAAAwJF13bcNBAAAAAAAgKWIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGGurHgAAAABgVTZOnr3m1148dXzBSVjCUfl5HpX3CcDRZecVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOsrXoAmGzj5Nnrev3FU8cXmuTGup73ebO8x8T73I+b6X0CAAAAAIeTnVcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjLG26gEAAAAAALh+GyfPXvNrL546vuAkANfHzisAAAAAAADGEK8AAAAAAAAYw20DAYBr4nYUAAAAANwIdl4BAAAAAAAwhngFAAAAAADAGOIVAAAAi6qqt1TV31bVhao6ucv591bV5ap6bOvr/auYEwAAmMlnXgEAALCYqnpVkk8n+fUkl5J8q6rOdPcTO5b+QXc/cOADAgAcQdfzudWJz67m4Nl5BQAAwJLuTnKhu7/X3T9J8qUk9614JgAA4CYiXgEAALCk25M8te340tZzO/1GVT1eVX9YVXcezGgAAMDNQLwCAABgSbXLc73j+E+SbHT3ryT58ySf3/VCVSeq6nxVnb98+fLCYwIAAFP5zCsAAACWdCnJ9p1UdyR5evuC7n522+HvJfnd3S7U3aeTnE6Szc3NnQEMADiirufzm3x20zx+nuzGzisAAACW9K0kd1XV66rqZ5Pcn+TM9gVVddu2w7cnefIA5wMAAIaz8woAAIDFdPeLVfVAkj9L8qokD3f3d6vqE0nOd/eZJL9TVW9P8mKS55K8d2UDAwAA44hXAAAALKq7zyU5t+O5j217/JEkHznouQAAgJuD2wYCAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAYa6seAODlbJw8e82vvXjq+IKTAAAAAABwEOy8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGCMtVUPAACHzcbJs9f82ounji84CQAAAADcfOy8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAx1lY9AABci42TZ6/5tRdPHV/sWrtdDwAAAAC4duIVAAAAwCG35D/+AgC40dw2EAAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHWVj0AAEfHxsmz1/zai6eOLzgJAAAAADCVnVcAAAAAAACMYefVwq5nV0Fy8+wsOCrvEwAAAAAAOFh2XgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjLFnvKqqh6vq+1X1naucr6r6VFVdqKrHq+oNy48JAAAAAADAUbCfnVefS/KWlzn/1iR3bX2dSPKZ6x8LAAAAAACAo2jPeNXdX03y3MssuS/JF/qKbyR5TVXdttSAAAAAAAAAHB1LfObV7Ume2nZ8aes5AAAAAAAAeEWWiFe1y3O968KqE1V1vqrOX758eYFvDQAAAAAAwGGyRLy6lOTObcd3JHl6t4Xdfbq7N7t7c319fYFvDQAAAAAAwGGyRLw6k+TddcWvJnm+u59Z4LoAAAAAAAAcMWt7LaiqLya5J8mtVXUpyceT3JIk3f1QknNJ7k1yIcmPkrzvRg0LAAAAAADA4bZnvOrud+5xvpN8YLGJAAAAAAAAOLKWuG0gAAAAAAAALEK8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAFlVVb6mqv62qC1V1cpfzP1dVf7B1/ptVtXHwUwIAAFOJVwAAACymql6V5NNJ3prkl5K8s6p+acey30ryP7r7Hyf5N0l+92CnBAAAJhOvAAAAWNLdSS509/e6+ydJvpTkvh1r7kvy+a3Hf5jk16qqDnBGAABgMPEKAACAJd2e5Kltx5e2ntt1TXe/mOT5JH//QKYDAADGW1v1AAAAABwqu+2g6mtYk6o6keREkhw7duz6J1vIxsmz1/zai6eOv+zx9Vryeku+z8mm/J0lPz3Lktfa7fiw8vM8XKb8PG/03/fk/19Zkve5txv5Ppf+O5v6Ppf+7y2rY+cVAAAAS7qU5M5tx3ckefpqa6pqLcnfTfLczgt19+nu3uzuzfX19Rs0LgAAMI14BQAAwJK+leSuqnpdVf1skvuTnNmx5kyS92w9fkeSv+zul+y8AgAAjia3DQQAAGAx3f1iVT2Q5M+SvCrJw9393ar6RJLz3X0myWeT/LuqupArO67uX93EAADANOIVAAAAi+ruc0nO7XjuY9se/88kv3nQcwEAADcHtw0EAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxvCZVwAAAPAKXDx1fNUjAADAoWbnFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADDG2qoHAAAAADgMLp46vuoRAAAOBTuvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDF85hUAAAAAALBvPueRG83OKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAMn3nFobNx8uw1v9a9WgEAAAAAYLXsvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhjX/Gqqt5SVX9bVReq6uQu599bVZer6rGtr/cvPyoAAAAAAACH3dpeC6rqVUk+neTXk1xK8q2qOtPdT+xY+gfd/cANmBEAAAAAAIAjYj87r+5OcqG7v9fdP0nypST33dixAAAAAAAAOIr2E69uT/LUtuNLW8/t9BtV9XhV/WFV3bnIdAAAAAAAABwp+4lXtctzveP4T5JsdPevJPnzJJ/f9UJVJ6rqfFWdv3z58iubFAAAAAAAgENvP/HqUpLtO6nuSPL09gXd/Wx3/6+tw99L8s93u1B3n+7uze7eXF9fv5Z5AQAAGKqq/l5VfaWq/tPWn79wlXX/u6oe2/o6c9BzAgAAs+0nXn0ryV1V9bqq+tkk9yf5qV8uquq2bYdvT/LkciMCAABwkziZ5C+6+64kf7F1vJsfd/c/2/p6+8GNBwAA3AzW9lrQ3S9W1QNJ/izJq5I83N3frapPJDnf3WeS/E5VvT3Ji0meS/LeGzgzAAAAM92X5J6tx59P8tdJ/tWqhgEAAG5Oe8arJOnuc0nO7XjuY9sefyTJR5YdDQAAgJvMP+zuZ5Kku5+pqn9wlXV/p6rO58o/gDzV3V8+sAkBgJW4eOr4qkfgCPC/s8NjX/EKAAAAkqSq/jzJa3c59a9fwWWOdffTVfWPkvxlVX27u//zLt/rRJITSXLs2LFrmhcAALj5iFcAAADsW3f/i6udq6r/VlW3be26ui3J969yjae3/vxeVf11ktcneUm86u7TSU4nyebmZi8wPgAAcBP4mVUPAAAAwKFxJsl7th6/J8l/2Lmgqn6hqn5u6/GtSd6c5IkDmxAAABjPzisAAACWcirJv6+q30ryX5P8ZpJU1WaSf9nd70/yT5L826r6P7nyDypPdbd4BcCR5TN6AF5KvAIAAGAR3f1skl/b5fnzSd6/9fg/JvmnBzwaAABwE3HbQAAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGWFv1AAAAAAD8tIunjq96hCNvyZ+BnycAvDJ2XgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjiFcAAAAAAACMIV4BAAAAAAAwhngFAAAAAADAGOIVAAAAAAAAY4hXAAAAAAAAjCFeAQAAAAAAMIZ4BQAAAAAAwBjiFQAAAAAAAGOIVwAAAAAAAIwhXgEAAAAAADCGeAUAAAAAAMAY4hUAAAAAAABjrK16AAAAAAAAOAwunjq+6hHgULDzCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7wCAAAAAABgDPEKAAAAAACAMcQrAAAAAAAAxhCvAAAAAAAAGEO8AgAAAAAAYAzxCgAAAAAAgDHEKwAAAAAAAMYQrwAAAAAAABhDvAIAAAAAAGAM8QoAAAAAAIAxxCsAAAAAAADGEK8AAAAAAAAYQ7yC/9ve/YZYel91AP8edpsqUWxrV61J1kQNaCkSdQyKIkVbTdMXa6XRDSgRlLVgpOIbo4LWQiGK/96USqTRKOo2tE272IVYaaX6Ju6mpjabNbqNG7NuSKyxal5oiTm+mGdhWGcms3Xm/p659/OBZe7z3N/ee/Z3+M08Z888vwsAAAAAAMyG5hUAAAAAAACzoXkFAAAAAADAbGheAQAAAAAAMBuaVwAAAAAAAMyG5hUAAAAAAACzoXkFAAAAAADAbGheAQAAAAAAMBuaVwAAAAAAAMyG5hUAAAAAAACzoXkFAAAAAADAbGheAQAAAAAAMBuaVwAAAAAAAMyG5hUAAAAAAACzoXkFAAAAAADAbGheAQAAAAAAMBuaVxfiz/EAAAonSURBVAAAAOyKqrqtqs5U1YtVtbbNuFuq6vGqOldVdy0yRgAAYP40rwAAANgtjyb5wSSf2GpAVR1I8u4kb0ry2iS3V9VrFxMeAACwHxwcHQAAAADLobvPJklVbTfs5iTnuvuJaezxJEeSPLbnAQIAAPuCO68AAABYpGuSPLXh+MJ0DgAAIIk7rwAAALgCVfXnSb5qk6d+sbs/vJOX2ORcb/Fex5IcS5LDhw/vOEYAAGB/07wCAABgx7r7Df/Pl7iQ5LoNx9cmubjFe92T5J4kWVtb27TBBQAALB/bBgIAALBIp5LcWFU3VNVVSY4mOTE4JgAAYEY0rwAAANgVVfWWqrqQ5DuSfKSqHpzOf3VVnUyS7n4hyZ1JHkxyNsn93X1mVMwAAMD82DYQAACAXdHdDyR5YJPzF5PcuuH4ZJKTCwwNAADYR9x5BQAAAAAAwGzsqHlVVbdU1eNVda6q7trk+ZdX1fum5x+qqut3O1AAAAAAAACW30tuG1hVB5K8O8kbk1xIcqqqTnT3YxuG/XiSf+vur6+qo0l+NckP70XAAAAAAABX4vzdbx4dAgBXYCd3Xt2c5Fx3P9Hdn09yPMmRy8YcSXLf9Pj9Sb63qmr3wgQAAAAAAGAV7KR5dU2SpzYcX5jObTqmu19I8u9Jvnw3AgQAAAAAAGB1VHdvP6DqtiTf390/MR3/aJKbu/unN4w5M425MB1/Zhrzr5e91rEkx5Lk8OHD3/rkk0/u5r/lC3b9XR/5gv+uW44BAGBdVT3c3Wuj42A5ra2t9enTp0eHsVLUygAA7KXtasid3Hl1Icl1G46vTXJxqzFVdTDJlyV57vIX6u57unutu9cOHTq0k9gBAAAAAABYITtpXp1KcmNV3VBVVyU5muTEZWNOJLljevzWJB/rl7qlCwAAAAAAAC5z8KUGdPcLVXVnkgeTHEhyb3efqap3Jjnd3SeSvDfJH1bVuazfcXV0L4MGAAAAAABgOb1k8ypJuvtkkpOXnfulDY//K8ltuxsaAAAAAAAAq2Yn2wYCAAAAAADAQmheAQAAAAAAMBuaVwAAAAAAAMyG5hUAAAAAAACzcXB0AHNw/u43jw4BAAAAZkWtDADAKO68AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2ajuHvPGVf+S5Mkhb37lXp3ks6ODWGHmfzw5GMv8jycH48nBWOZ/vP2Sg6/p7kOjg2A5qSG5AuZ/PDkYy/yPJwfjycFY5n+8/ZKDLWvIYc2r/aSqTnf32ug4VpX5H08OxjL/48nBeHIwlvkfTw5gf7FmxzL/48nBWOZ/PDkYTw7GMv/jLUMObBsIAAAAAADAbGheAQAAAAAAMBuaVztzz+gAVpz5H08OxjL/48nBeHIwlvkfTw5gf7FmxzL/48nBWOZ/PDkYTw7GMv/j7fsc+MwrAAAAAAAAZsOdVwAAAAAAAMyG5tU2quqWqnq8qs5V1V2j41lFVXW+qj5dVY9U1enR8ayCqrq3qp6tqkc3nHtVVX20qv5h+vrKkTEusy3m/x1V9c/TOnikqm4dGeMyq6rrqurjVXW2qs5U1dun89bAgmyTA+tgQarqi6rqr6vqU1MOfmU6f0NVPTStg/dV1VWjY11G28z/71fVP25YAzeNjhX4v9SQ46khF08NOZYaciw15HhqyPHUkGMtcw1p28AtVNWBJH+f5I1JLiQ5leT27n5saGArpqrOJ1nr7s+OjmVVVNV3J3k+yR909+umc7+W5Lnuvnsqwl/Z3T83Ms5ltcX8vyPJ89396yNjWwVV9Zokr+nuT1bVlyZ5OMkPJPmxWAMLsU0OfijWwUJUVSW5urufr6qXJfmrJG9P8rNJPtjdx6vqd5J8qrvfMzLWZbTN/L8tyZ929/uHBghsSQ05D2rIxVNDjqWGHEsNOZ4acjw15FjLXEO682prNyc5191PdPfnkxxPcmRwTLDnuvsTSZ677PSRJPdNj+/L+kUAe2CL+WdBuvvp7v7k9Pg/k5xNck2sgYXZJgcsSK97fjp82fSnk3xPkksXvdbBHtlm/oH5U0OyktSQY6khx1JDjqeGHE8NOdYy15CaV1u7JslTG44vxDe+ETrJn1XVw1V1bHQwK+wru/vpZP2iIMlXDI5nFd1ZVX87bQlhu4EFqKrrk3xzkodiDQxxWQ4S62BhqupAVT2S5NkkH03ymSSf6+4XpiGui/bQ5fPf3ZfWwLumNfBbVfXygSECm1NDzoMach5cP4/n2nnB1JDjqSHHUUOOtaw1pObV1mqTc0vRsdxnvrO7vyXJm5L81HQ7PKya9yT5uiQ3JXk6yW+MDWf5VdWXJPlAkp/p7v8YHc8q2iQH1sECdff/dPdNSa7N+p0E37jZsMVGtToun/+qel2Sn0/yDUm+Lcmrkth2BuZHDTkPakhw7bxwasjx1JBjqSHHWtYaUvNqaxeSXLfh+NokFwfFsrK6++L09dkkD2T9mx+L98y0h/ClvYSfHRzPSunuZ6YfQi8m+d1YB3tq2h/4A0n+qLs/OJ22BhZosxxYB2N09+eS/EWSb0/yiqo6OD3lumgBNsz/LdN2KN3d/53k92INwBypIWdADTkbrp8Hcu28WGrI8dSQ86GGHGvZakjNq62dSnJjVd1QVVclOZrkxOCYVkpVXT190GKq6uok35fk0bFRrawTSe6YHt+R5MMDY1k5ly54J2+JdbBnpg+5fG+Ss939mxuesgYWZKscWAeLU1WHquoV0+MvTvKGrO8b//Ekb52GWQd7ZIv5/7sN//lRWd8r3hqA+VFDDqaGnBXXzwO5dl4cNeR4asjx1JBjLXMNWd3u1ttKVd2a5LeTHEhyb3e/a3BIK6WqvjbrvymXJAeT/LEc7L2q+pMkr0/y6iTPJPnlJB9Kcn+Sw0n+Kclt3e0DYffAFvP/+qzf5t5Jzif5yUt7Z7O7quq7kvxlkk8neXE6/QtZ3y/bGliAbXJwe6yDhaiqb8r6h+keyPovOt3f3e+cfi4fz/p2A3+T5Eem3+BiF20z/x9Lcijr25I9kuRtGz6UF5gJNeRYasgx1JBjqSHHUkOOp4YcTw051jLXkJpXAAAAAAAAzIZtAwEAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDY0rwAAAAAAAJgNzSsAAAAAAABmQ/MKAAAAAACA2dC8AgAAAAAAYDb+F3aIZgjEer3nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 2160x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(30, 10))\n",
    "plt.subplot(1, 2, 1).bar(np.arange(len(pr)), pr)\n",
    "plt.subplot(1, 2, 2).bar(np.arange(len(pr)), sp2.ft(pr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12 out of 36 Fourier coefficients are zero.\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(np.abs(sp2.ft(pr)) < 1e-8), \"out of\", len(sp2.lattice), \"Fourier coefficients are zero.\")"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}