SCA_on_NTRU.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of key:  (509,)\n",
      "shape of trace dataset:  (100, 507, 950)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABZYklEQVR4nO2dd3hUxfrHP7PpvYeEdHpLgYQuIlJERWzoFbkK9nYVvdd67fV61WtD7Ijlp2LHggURRIr0Ir0npPfek53fH7O7yaZAIBsSNvN5nn32nJk558zZTd6d884731dIKdFoNBqNfWLo7A5oNBqNpuPQRl6j0WjsGG3kNRqNxo7RRl6j0WjsGG3kNRqNxo5x7OwONCYwMFBGR0d3djc0Go3mtGLz5s15Usqgluq6lJGPjo5m06ZNnd0NjUajOa0QQqS0VmcTd40QwlcI8aUQYq8QYo8QYrQQwl8I8asQ4oDp3c8W19JoNBpN27GVT/4V4Gcp5QAgHtgD3A/8JqXsC/xm2tdoNBrNKaTdRl4I4QOcCSwAkFLWSCmLgAuBD0zNPgAuau+1NBqNRnNi2MInHwPkAguFEPHAZmAu0ENKmWlqkwX0aOlgIcSNwI0AkZGRNuiOpjtRW1tLWloaVVVVnd0VjabDcXV1JTw8HCcnpzYfYwsj7wgMA26XUq4XQrxCE9eMlFIKIVoUyZFSvg28DZCUlKSFdDQnRFpaGl5eXkRHRyOE6OzuaDQdhpSS/Px80tLSiImJafNxtvDJpwFpUsr1pv0vUUY/WwgRCmB6z7HBtTQaK6qqqggICNAGXmP3CCEICAg44afWdht5KWUWkCqE6G8qmgjsBr4DZpvKZgPftvdaGk1LaAOv6S6czN+6reLkbwc+FkI4A4eBa1A/IJ8LIa4DUoDLbXStE6e6FDa+C4MvAb+oTuuGRqPRnGpsEkIppdwmpUySUsZJKS+SUhZKKfOllBOllH2llJOklAW2uFabqSqG7N1Qng87voRlj8GKZ05pFzTdg2uvvZbg4GCGDBliVV5QUMDkyZPp27cvkydPprCwEIDff/+dtWvXWtrNmTOHL7/88rjXcXBwICEhwfJ69tlnATjrrLOIjIykcW6Iiy66CE9PT1vcXqu0td8tsXfvXhISEhg6dCiHDh2yqnvwwQeJiIhoU/+3bdvGjz/+eFJ96Cjef/99MjIyLPvXX389u3fv7rT+2K92zcLz4Y3R8PooSDVNF2Tv7Nw+aeySOXPm8PPPPzcrf/bZZ5k4cSIHDhxg4sSJFqPc1Mi3FTc3N7Zt22Z53X9/Q3yDr68va9asAaCoqIjMzMzWTtMlWLx4MTNmzGDr1q307t3bqu6CCy5gw4YNbTrP6WDk3333XQYNGtRp/bFPI1+WC9k7IHgQlOfAXtMfQWEK6ExYGhtz5pln4u/v36z822+/ZfZsNS01e/ZsFi9eTHJyMm+++SYvvfQSCQkJrFq1CoA//viDMWPG0KtXr5MaHV9xxRUsWrQIgK+//ppLLrmkxXbJyckMGDCAWbNmMXDgQGbMmEFFRQUAmzdvZvz48SQmJnLOOedYfijeeecdhg8fTnx8PJdeeqmlfWMefvhh5syZQ319vVX5tm3bGDVqFHFxcVx88cUUFhby448/8vLLL/PGG28wYcKEZucaNWoUoaGhzcq/+OILhgwZQnx8PGeeeSY1NTU88sgjfPbZZyQkJPDZZ59Zta+vr+eee+5h+PDhxMXF8dZbbwHw0ksvce211wKwY8cOhgwZQkVFBYcOHWLq1KkkJiYybtw49u7dC0B2djYXX3wx8fHxxMfHs3btWpKTk62e3F544QUee+wxvvzySzZt2sSsWbNISEigsrKSs846yyLX8umnnxIbG8uQIUO47777LMd7enry4IMPEh8fz6hRo8jOzm7x+zsZupR2zUmTswe2fQIT/g1ObpC3X5WPvAm+nwvVxWq/plS5cdx8O62rmo7j8e93sTujxKbnHNTTm0cvGHxSx2ZnZ1uMVUhICNnZ2URHR3PzzTfj6enJ3XffDcCCBQvIzMxk9erV7N27l+nTpzNjxoxm56usrCQhIcGy/8ADD/C3v/0NgIkTJ3LDDTdQX1/PokWLePvtt3nyySdb7Ne+fftYsGABY8eO5dprr+X1119n7ty53H777Xz77bcEBQXx2Wef8eCDD/Lee+9xySWXcMMNNwDw0EMPsWDBAm6//XbL+e655x5KS0tZuHBhs4nBq6++mnnz5jF+/HgeeeQRHn/8cV5++eVmn0FbeOKJJ/jll18ICwujqKgIZ2dnnnjiCTZt2sRrr73WrP2CBQvw8fFh48aNVFdXM3bsWKZMmcLcuXM566yz+Oabb3j66ad56623cHd358Ybb+TNN9+kb9++rF+/nltvvZXly5dzxx13MH78eL755hvq6+spKyuzuN6aMmPGDF577TVeeOEFkpKSrOoyMjK477772Lx5M35+fkyZMoXFixdz0UUXUV5ezqhRo3j66ae59957eeedd3jooYfa/NkcC/sw8oUpsPZViBoL/adCcaoqjxoLrj7KsIfEQdZfkLsXIkd1bn813Q4hxDEjIy666CIMBgODBg1qdRRndte0hIODA2eccQaLFi2isrKSY6m5RkREMHbsWAD+/ve/8+qrrzJ16lR27tzJ5MmTATUKNv9A7dy5k4ceeoiioiLKyso455xzLOd68sknGTlyJG+//Xaz6xQXF1NUVMT48eMB9TRz2WWXtdqv4zF27FjmzJnD5Zdf3uqTSmOWLl3KX3/9ZXkyKi4u5sCBA8TExPD+++8TFxfHTTfdxNixYykrK2Pt2rVW/auurgZg+fLlfPjhh4D6nH18fFo18sdi48aNnHXWWQQFKbHIWbNm8ccff3DRRRfh7OzMtGnTAEhMTOTXX3894fO3hn0Y+ciR6v3n+5SRLzgCCPAJh6gzYN8SOOMu+PIaSFmrjbydcrIj7o6iR48eZGZmEhoaSmZmJsHBwa22dXFxsWzLk3QpXnHFFVx88cU89thjx2zX9MdGCIGUksGDB/Pnn382az9nzhwWL15MfHw877//Pr///rulbvjw4WzevJmCgoIWXVa25M0332T9+vUsWbKExMRENm/efMz2UkrmzZtn9aNk5sCBA3h6elp850ajEV9f31Z/RJvi6OiI0Wi07Ld3xbWTk5Ple3FwcKCurq5d52uMffjk3fxUeGRhMnwwHZJXQY/BynUz/VWY9SUMvhgC+8Nvj0NNc5+iRmNrpk+fzgcfKPmmDz74gAsvvBAALy8vSktLbX69cePG8cADDzBz5sxjtjt69KjFmH/yySecccYZ9O/fn9zcXEt5bW0tu3btAqC0tJTQ0FBqa2v5+OOPrc41depU7r//fs4///xm9+Tj44Ofn59l3uGjjz6yjOpPhkOHDjFy5EieeOIJgoKCSE1NPeZnec455/DGG29QW1sLwP79+ykvL6e4uJg77riDP/74g/z8fL788ku8vb2JiYnhiy++ANQPxPbt2wHlCnvjjTcA9YRTXFxMjx49yMnJIT8/n+rqan744QfLdVvr04gRI1i5ciV5eXnU19fz6aeftuvzaCv2YeQBJj0KvSfCkZWQsgYiR6tyj0DoOxmEgIEXqLLUdZ3XT43dMXPmTEaPHs2+ffsIDw9nwYIFANx///38+uuv9O3bl2XLllmiYS644AK++eYbq4nXtmD2yZtfjaNrQI3I7777bgIDA495nv79+zN//nwGDhxIYWEht9xyC87Oznz55Zfcd999xMfHk5CQYIkAMrtkxo4dy4ABA5qd77LLLuOGG25g+vTpVFZWWtV98MEH3HPPPcTFxbFt2zYeeeSR497nvffeS3h4OBUVFYSHh1ueTO655x7LpOWYMWOIj49nwoQJ7N69u8WJ1+uvv55BgwYxbNgwhgwZwk033URdXR133XUXt912G/369WPBggXcf//95OTk8PHHH7NgwQLi4+MZPHgw336r1m++8sorrFixgtjYWBITE9m9ezdOTk488sgjjBgxgsmTJ1t9LnPmzOHmm2+2TLyaCQ0N5dlnn2XChAnEx8eTmJho+eHvSMTJPhp2BElJSbJdSUOqimFeElQWUD/zc7KDxhDk5UJ+WQ0hPq5qUdSzkTDubjj7Qdt1XNNp7Nmzh4EDB3Z2N04bkpOTmTZtGjt36nDi05WW/uaFEJullEkttbcPn7wZVx+4ez9IySvLDvDqguXEBHpwJK+cH24/gyFhPmoCNuXEY5Q1Go3mdMR+3DVmhACDgRX7cgE4klcOwOqDeao+aiykb4K6Gvjpfvjq+s7qqUZzyomOjtaj+G6GXRn5balFvP77QTKLK9mXVYoQMHVwCAAbjphUFUJioa5KhVKufwN2fAFFqbDrGzi0vBN7r9FoNLbHLtw1VbX1HC2o4NHvdrE9tYi1B/OpqTfy7tVJTBrUg399vp1VB9TIHp9w9b71/xpOsP9n+NG0KOPRIvU0oNFoNHaAXRj5X3ZlMXfRNsu+2TWTFK1yh0f4u5FTWk11XT0uvqbsUxveajjButcbtvMPQWCfju6yRqPRnBLswl0T5uvWrKx/Dy983Z2t6i98bQ3rCzzAxVs1GnIphMZDweGGA4uSO7q7Go1Gc8qwDyPv12DkrxkbzXmxIdwxsa+lbHy/IKbH9+RgThm/7cuFKU+pxVETHwWfCOuTFZkkEZJXw5pXVdilRtMKqampTJgwgUGDBjF48GBeeeUVS52WGm6dY0kNT5061RKrfvPNNzcTPWuJpoJhnc3ixYut5IUfeeQRli1b1il9sQsjH+zlatk+e0Awr89K5Py4BhW7YG9XXp05lAh/d9KLKiFxNlz2vkog4t9LNfIOB0c3NfkqJXx+Nfz6MPz1GRpNazg6OvK///2P3bt3s27dOubPn2/559ZSw61zLKnhzz//nO3bt7Nz505yc3Mtq1AbY8tl/x1BUyP/xBNPMGnSpE7pi10YeQeDoFeQBwBR/h6ttgv3c+NofhNJgygl1ETECPXa851KMFKRr8p1TL3mGISGhjJs2DBALWcfOHAg6enpgJYahpOTGvb2Vu7Uuro6ampqLJouZ511FnfeeSdJSUm88sorbN682SL/O3/+/FY/m+eff94iN/zoo48C8M033zBx4kSklGRmZtKvXz+ysrLIzc3l0ksvZfjw4QwfPtzyw1lWVsY111xDbGwscXFxfPXVVwBWT0tffvklc+bMYe3atXz33Xfcc889JCQkcOjQIaunnt9++42hQ4cSGxvLtddeaxFCi46O5tFHH2XYsGHExsZapI7bjZTSJi/AAdgK/GDajwHWAweBzwDn450jMTFRniypBeVy5b6cY7Z57uc9Muq+H+Q5L62UO9KKVGF9vZR7fpDfrNoq/zF/sZSPeje8Xhsh5QsDpDQapcw7KOW7k6U88OtJ91Fje3bv3t2w8+N9Ur53nm1fP97X5r4cOXJERkREyOLiYimllD4+PpY6o9Fo2X/00Ufl888/b6mbPXu2nDFjhqyvr5e7du2SvXv3bvH8BoNBxsfHW16LFi2SUko5fvx4uW7dOhkbGyvr6urk5MmT5ZEjR6SHh0eLfQTk6tWrpZRSXnPNNfL555+XNTU1cvTo0TInR/0PLVq0SF5zzTVSSinz8vIsxz/44IPy1VdftfT7iy++kHfffbe86aabpNFobHa92NhY+fvvv0sppXz44Yfl3LlzW/wMmjJlyhTp6+srZ86cKevq6iz3ecstt1ide+XKlVJKKe+++245ePDgZuf55Zdf5A033CCNRqOsr6+X559/vuWYWbNmyXnz5snzzz9ffvLJJ1JKKWfOnClXrVolpZQyJSVFDhgwQEop5b333mvpu5RSFhQUSCml1Wf8xRdfyNmzZ1t9NmbM+5WVlTI8PFzu27dPSinlVVddJV966SUppZRRUVGWz3b+/Pnyuuuua/GzsfqbNwFskq3YVVuO5OcCexrt/xd4SUrZBygErrPhtZoR7ufOmf2CjtlmREwAAHuzSvlmazqpBRW8/sdh5mf2584f0vn+qCNZYVNUY4MjJF0LpRlK+Gz3YpVhan1zSdVjUpQK69+CrB2w4R2o79qPmZqTo6ysjEsvvZSXX37ZMhJtjC2lhs0vs5Y8tE9qePXq1ezbt88iNZyQkMBTTz1FWloaoKSGx40bR2xsLB9//LFFuAyUrk1xcTFvvvlms/trSWr4jz/+aLVfjfnll1/IzMykurqa5csb1q+Y77moqIiioiLOPPNMAK666qoWz7N06VKWLl3K0KFDGTZsGHv37uXAgQMAzJs3j//85z+4uLhYRN2WLVvGP/7xDxISEpg+fTolJSWUlZWxbNkybrvtNst5/fz82nQfTdm3bx8xMTH069cPaP6ZmJ/AEhMTSU5OPqlrNMUmIZRCiHDgfOBp4J9CfdtnA1eamnwAPAa8YYvrnSyJUX4EeDiTX17Droxinl6yh593ZVm1OVIbQAiAd0+IHqcKj/4JGdvUduZ25bNvayz9T/cpqWMzvpHQr5H0aXUZuHTsBFm34dxnO+WytbW1XHrppcyaNcvKTaKlhtuHq6srF154Id9++61F597Do3V3bEtIKXnggQe46aabmtWlpaVhMBjIzs7GaDRiMBgwGo2sW7cOV1fXFs7WnMafZXvlhqHh78CWcsO2Gsm/DNwLmAWWA4AiKaW5l2lAmI2uddJ4ujiy4cFJTI/vybrDBfy8K4sJ/YPY99RU9j91LqN7BbCl3PQL7R0OQQPA1VepWpoTkZRlQdHRtl80c5v1vjnfLKhkJ8/1gv2/tOe2NJ2IlJLrrruOgQMH8s9//tOqTksNn7jUcFlZmWUuoK6ujiVLlrSofOnr64uvry+rV68GaNYvM+eccw7vvfceZWVlAKSnp5OTk0NdXR3XXnstn376KQMHDuTFF18EYMqUKcybN89yvFlffvLkyVZ+f3OkVI8ePdizZw9Go5FvvvnGUt/ad9y/f3+Sk5M5ePBgmz+T9tJuIy+EmAbkSCmPreDf+vE3CiE2CSE25ebmtrc7x8XBIKxCLm+d0AcXRwecHQ0Mj/Hn9YIkKqe/DRfNB4MBosZAyp/K7WKOxClMbvsFa5tMUjU+9shKqK/WRv40Zs2aNXz00UcsX77cEtpoTiytpYZPXGq4vLyc6dOnExcXR0JCAsHBwdx8880ttl24cCG33XYbCQkJrT79TJkyhSuvvJLRo0cTGxvLjBkzKC0t5ZlnnmHcuHGcccYZvPjii7z77rvs2bOHV199lU2bNhEXF8egQYN48803AZX2sLCw0JJjdsWKFYCKoJo2bRpjxoyxykt7xRVX8PzzzzcLEXV1dWXhwoVcdtllxMbGYjAYWr0/W9FuqWEhxH+Aq4A6wBXwBr4BzgFCpJR1QojRwGNSyuYpWhrRbqnhNvLN1jTu+kwlBDjyn/Msj1xrDuYx6931vH/NcM7qb3q0XvOqCqUESPg7bPs/uHA+DP07GI2weSF4hcCA85tfyFgPTwQAjT7j8BFwvSm11ze3wPZPVMLxW5s/JmuOj5YaPjG01PDpz4lKDbd7JC+lfEBKGS6ljAauAJZLKWcBKwBzNuLZwLftvZatGNs7EHdnB64eHWXlUxsa6YujQbAxuaChsTnEElTaQAdnSN8MFQUq3HLJP2HRlcq33pTKQkCCRzAYnCAsEYrTGupTVHgWObvV+TQajcbGdGSc/H2oSdiDKB/9gg681gkR7O3KX49O4fHp1jlB3Z0dGRLm06BYCRAa17Ad2A/CkmDTe/BcDHwxu6EubUPzC5Wb3E9T/wMPZkGfyVCaqWSOi9OhKKUhW9VRna1K0/FoqeHuh02NvJTydynlNNP2YSnlCCllHynlZVLKalteq704OhhaDGkbEePPxuRC3l2l9GyW7s3nHt8X2TLqVbVY6qL5cNYDDQdc+TkIQ/NFU0fXweumhOHeYeDgaFLAlFCSBr89oepG366eDn68GxaeD6Uth89pWqe9LkeN5nThZP7W7WLFqy25LFFJET+1ZA+lVbV8vP4oX2SF8GLaAH7elc2GYl8Y9y/ln7/4bRUOGRKnJmcbs/kDy2atVxgfrUshy8U0cbvqRTiwVEXuhCfBhH+rp4SU1Ur2WNNmXF1dyc/P14ZeY/dIKcnPz29zeKcZu5AatiV9e3jxwmXx3P3Fdt5ZdYQtKSpUavXBPIuE8ZH/nIe4qNEy6qixsPFdqKsGR1O8c65pXZhXT35MgYcX72R5Xz8WeofD1o9U3WXvg8EBzrgLxt4JL/SF3d/CsKu1pn0bCQ8PJy0tjVMRmaXRdDaurq6Eh4ef0DHayLfAjMRwXvltP6/+plbGDe7pza6MEkt9elEl4X7uDQdEjYZ18yFjq5qcBeV2ib8Sps/joSd/A2BnVgXcv12FTQoDODWSSBZC/VjsXgyrX4Jx1jHXmpZxcnIiJiams7uh0XRZtLumFd65OonnZsTx6syhXD/O2og88f1u3l9zpKEgcrR6T14N3/4DktdAeQ54hfDDrhxKq9SasLyyamoxsK/AyO1f7qWqtp4Xl+7j8rf+JCW/HM55BhDw2+Na4lij0dgEbeRbYUCIN5cnRTA9vieDQn2s6pbuzuax73c3+IE9AiGwP2x+X7li3j8PjHXKyG9Xq/f+NbkfUkJWcRX3f/0X32/PYO2hPF5dfpANRwp49qe94BMGE02LRZY+pEI1NRqNph1oI98G+gZ7EubrRkKEr1V5Sn4FVbX1VNTUqZWxZukDE9KzBxuSC7h0WDhDI5VcQlphJbX1Sv3hg7UplrY/7cyirLoORt8GXqHqB2PxraqytgrKWvE5Swnl+Ta5T41GY39oI98GDAbBynvO4qtbxhAT2CCQtDmlkHNfWcX4539HNl40ZSK9zouC8hpGxPhZpBTSiyoprqwFYOV+Zbj/d1k8AE9+v1tN3N65AyY8CLl7oTwP3p0EL/RRbqCmbFoAz/eCfToqR6PRNEcb+Tbi6GDAwSB4b85w3vx7IkLAppRCjuSVk1taTUnM+XDpAohvEIf696/KiI+ICSDUR4U9pRZUkFXcoFbXK9CDi4eGEeTlwtZUFcmDgxPEmESL9nwP2TvU9kFT+rDcffDeuSoUM6tJnUaj0TRCR9ecIDGBHsQEetDDy5XFW9Mt5Z9szsTPfQSXD/HHsP1TALYWe3DVqCiiA9wRQhDk5cLyvTnU1ktuOrMXB3LKOD82FINBMGdMNM//so/C8hr8PJyh51BwdIXVLzZc/OifUFMBn89WIZqZ2yHa9AShM1hpNJoW0Eb+JIkMcLeSP/jvzypVV8TVgxnp4svyit4MCg/kyYsakgv39HVje2oRAFeOjCQqoMH1M9Tk79+ZUcy4vkHg6Ax9JsHeH8AzBAZfBJveo3bLRziZY/Bry9WiKlD6N43j9DUajQbtrjlpkqLURGqgp7VRffyXFAYWv8pNtXfx30vjrOrCTX55g1AGvzER/iruPqOokVTr5R/B/alw107odRbU15D14/PUO3nB7Vsa2vWIBaS1+JlGo9GgjfxJMzxGZcFxcTTw1S1jeOWKBAb39GZ/dhmurq68O3sE0YHWWWziwlQopqeLI04O1h99iI8rBgHphY2MvMHAvDXZXLFgEykesQBEGHI57J3UoG0PEHeZem8S3aPRaDTaXXOSjO4VwAXxPTlncA8So/xIjPLD0WBg0cajjO8XxMSBPZodMy2+JxuTCxjVK6BZnZODgR7erqQ1Gsnnllbzv1/3A/D9/kD6+M7ELX8XG53P524hYNpLkLYZYi+HXx9RPvpeZ3XYPWs0mo5h4ZojxIb5kBRt2xSKoI38SePq5MC8mUOtys6PC+X8uNBWjoAwXzfenT38mPWNR/KfbWxIM7ghuZCljn/nr9piQgtcuRtUovGka1WDgD5q8nXs3JO6H41G0zkYjZKnluzh5vG9OsTIa3dNFyLcz43UgoZ0gelFlXi7OjJrZCRbUgotPwDZJVWWBVUWosYoeWNjk3KNRtOlKamqpd4o8ffomKAJbeS7ELHhvmQUVzHr3XXc+OEmDuWUExPoQe8gT8qq68gvryHS3x2jSR7BisgxUFXUoH6p0WhOCwrKawDw93DqkPNrI9+FGGF6VFtzMJ+lu7PZkFxAn2Avq8TjI0wTvukm3/321CJ+3Z2tRvKg4+U1mtOMBiOvR/J2z8BQL4K9rL/okTH+DAzxtuxPMk3oml03F85fww0fbqLcrafKQKWNvEZzWmEx8u7OHXL+dk+8CiEigA+BHoAE3pZSviKE8Ac+A6KBZOByKWVhe69nzzg6GFh939nUGyUTXvidrJIqRsT4Exngzp4npiKRGEzJRNKLKskva8iomFtWg0fUGDiySomW6aQjGs1pQWGFych7doyRt8VIvg74l5RyEDAKuE0IMQi4H/hNStkX+M20rzkOzo4G3Jwd+Oi6Ebx9VaIl1t7N2QF3Z0dcnRwI9HQho6iS7JIGI794WzofpIdCWRb89XlndV+j0Zwg+R08km+3kZdSZkopt5i2S4E9QBhwIWBOdPoBcFF7r9Wd6NvDiymDQ1qsC/NzY9HGVP440CA//PKyA8zLGKh2dn19Krqo0WhsQGF5DW5ODrg5O3TI+W3qkxdCRANDgfVADyllpqkqC+XOaemYG4UQm4QQm3SezrZhXjn77E97rcrz8OFg+KVKyEyHUmo0XZ6y6jpyS6sZ7F7UYf+zNjPyQghP4CvgTillSeM6qVIoyZaOk1K+LaVMklImBQUF2ao7ds0TFw7mrkn9Wqw74hkPVcVKsEyj0XRZDuaUMuTRX/huWxof1twFP3eMR9smRl4I4YQy8B9LKc2+gmwhRKipPhTIscW1NCCEYFy/QMu+n3tDfO1uR5PqpY6y0Wi6NCv35wEwUKTgLisgvPXV8O2h3UZeCCGABcAeKWUj8XO+A2abtmcD37b3WpoG4sJ8uCihJ09eNISARkqYacYAFUp5VBt5jaYrczCnFIC/BZuEBaNGd8h1bKFdMxa4CtghhNhmKvs38CzwuRDiOiAFuNwG19KYcHQw8PIVSjvn43UNuWLLa+vVwqgDS3UopUbThUkrrCQu3Ierg9JBRIJPeIdcp91GXkq5GmjNkkxs7/k1x6exjk15dT30HwM7voAN78DIGzuxZxqNpjXSiyrp38MLUjdArwkddh294tUOKK+uB5SKZUVNHcRdoSpyG0XfVBbCpzNVzliNRtOpSCnJKKok2gsoy4aglgMpbIGWGrYD5s8aypK/sjhaUE56URU4u0NIXEMSkYytsOVD2PcjZO+CgRd0boc1mm5OQXkNVbVG+roWqQKfiA67ljbydkBilD+JUf7MXbRVjeQB/GMgY5va/nQmlJqWLJTpICeNprPZlVGCA/X0Nx5SBX4xHXYt7a6xIzxcHCmrMhn58BFQlAIrn1MGfuIjcNYDUFcJ1WXWB66dB4/5QG1l85NqNBqb8vnGVK5+bwOfOj/N4PX3gKMbhMYd/8CTRBt5OyLAw5nCihrqjRICTT6+FU+r9xE3gW+U2i7LbjiorhqWPgTAfxd8bJEwpuAw/N8M+OIaqG2iXa/RaE6In3Zk8u6qw1BfS9TKudzgsYpEp2SVrvPKReDYMTLDoN01dkWwlwtGCfnl1QT7NvLxJV0HLp7gZVKWSN0AR1aq7UYG3Dl1LfuWVhE2MhGSV8PBX00nHgS+EcqX72ydnFyj0RyfWz7eAsBlkWWMLPuNkfymKvpN7fC8zNrI2xFBJi36nJJqgoOiwNkTaivgbDVSx9Nk5BffbHWcRJAuA7jL6SvY/RXsBmLOBJ9IKD4KK55SDatLYcQNp+huNBr7470lq7ircUEH+uLNaHeNHRHk5QpAblm1irC5ez/cnwrupuTAno1ULeOugPFKK6PAsw+/1VsnJefIH9BvCoz7V0NZ8uqO7L5GY5co6S5FXvpB68qIER1+fT2StyPMWaVyzTrzTVwrhdIDP/NO77Oh/1TwCGTeNi+8ardD9a/kSh+CRDEA8w4Fk+Y/mqemhuC09QOlh6NX0Wo0J0RxZS0Ao3sFcLkzyGSBMOs1mgdgHYgeydsRZnfNW38csmSNWrY7m9/3qbDJjSlFDY2jRoOrD4y4gRXFPcgMP4/yvhcyvfopPqybzIr6eH6vi+WznaVsDLwEEq6E8hy1qKox+Ycgd9+puD2N5rRCSsmKfTlklah5rytHRhLvWYLwjYTJT8BFb56Sfmgjb0e4OjkQFeDOodxy5i0/yOHcMq7/cBNzFm4ks7iS9KJKvqofR7V0pN5bTcwajZLMoioCA4Nwv/IDKtxCeKTuGu5yeoiFt05GCNiQXNCwWMO8wMrMvGEwv+MfOTWa041vt2Vwx8KVzF+mBkFBns6QuV2tYRk7FxJmnpJ+aCNvZ/xy55mM6uXPusP5rD9SYCl/6JudpBdWcnftTQyqXsj+bKWAl11aRU29kTA/N4QQxIWrhCTPXByLt6sTA0O82XCkoEE8qaiRka9rSD/YbISv0XQnairUWpMN71iKvDa9xg7X67n9wDWAZPDGByD/AESMOqVd0z55O8PVyYExvQN5adl+lu3OJtDTGS9XJ37bm0NcuA8SA/XAvOUHKK+uZ2CoNwCxYQ3Gfe2hPKYMUpE4I2L8WbTxKLVeiTiB9Ui+PK9huyQD3PzQaLolaRvV+6+PsDVkBk8v2cP9mcvAAP1EKtMM6/BMWaYi3EbedEq7pkfydsiIGH+khN/25jAixp9HLhgEwF9pxZaR+o87sli5P5c3Vx5CCBgY6gVAhL87fxseiaOD+tMYEuZDVa2RtGp3tTKvOK3hQhWNjLyWS9B0Z8xJeoz1fLIuhR0p2cSJQ3xdfwYArznPQ1QWwNkPn5LJ1sZoI2+HJET44uOmskWN7xfEhP7B9ApUkTbmEXtjgr1ccHFsOYlwmK8bgBI+842AoqMNlRX5DdvlOj+vphuTska911dz9Mg+Rhr24CzqKek9ncX1Yxra9Rp/yrumjbwd4urkwJr7z2bDgxP52/BIQI3IAZKi/SyG/uwBwQD0NBnylgj3Mxv5CuWXt3LXNDLyjaUSNJruRF0NpG2CoAEA9CjezofO/wXAGD6SR2vnqHYBfcE38pR3T/vk7RRPF0c8XRq+3n+fN5CRvfw5d0goceG+bE4u5JzBIXyzNY2k6NYfH0N8XDEISC+sVBE2mX81VJY0ct0cWApH18GUp1T0gEZj7/z+rPLFh49Qwn9xl8NvT/A3h99V/aTHCPQKpph0nvR7hoevuaxTuqmNfDchxMeVWSOVQFnvIE96B3kCMGfssQ2yk4OBHt6uyl0TMQC2fABb/w/iZ6pIG1dfJZ9w5A91QM8EOPOeDrwTjaYLUF2mFF5lPRxcBsCBkGlEi2cZ67AL6eqLGDOXkWU1nNEnkDPGDQfP4E7paoe7a4QQU4UQ+4QQB4UQ93f09TS2J8LfnYO5ZdB7AhIB394GB35VrhvfCPAMamhs1rDXaOyZrB3KwAf0AWCHMZrLPk1hj0Gpv4qoMWBQA6T/u34kE/p3joGHDjbyQggHYD5wLjAImCmEGNSR19TYnsQoP3alFzN/lxMjquZTIx3J37VcjeR9IsGj0R+w9s1rugPmAIS/fczaC1Yyo+YxiipqWVXdS5VHju68vjWho0fyI4CDUsrDUsoaYBFwYQdfU2NjRkT7U2eUPP/LPpx9Q9kme1Ow53cVTtl0JF+qjbymG1BsMvK+keyr8qEaZwDerTuPwyMeh8TZndg5azrayIcBjdfBp5nKLAghbhRCbBJCbMrN1WF4XZHE6IZFTjMSwzngGkvf2n1QU8rmiiAKjY2E0EozYfe3UJrVCT3VaE4N1WnbqXDrSZ2DqwpKMFGINz0n3650oboInR5CKaV8W0qZJKVMCgoKOv4BmlOOt6sT8RG+eDg7cMtZvRk07hIA6qXgn5v8eG+/GsUQMRKMtfD51fCznn7R2Al11WA0WhWVHdnIsrIoft6VRUZxJSHerrg4GhgZ44+rU8trTjqLjo6uSQcapyEPN5VpTjM+u3EUtfVGXJ0c6DN8CklL3qAewZC+vZh3oAfX/utW/LLWQup6dUCx/po1dkDqBnjvHIgaC3N+UGX1dfjW5pAiR5GfXEh6YSV9e3jy2pXDcHXq9HFzMzq6RxuBvkKIGCGEM3AF8F0HX1PTAbg6OeDlqlbRerk68cxVZ3PPxWO5/ey+gODltQUNmafAejWsLaivhW//AUv+1WxUpdF0GPt/BmmE5FXwwQWw5SMozcQBIxkygOT8ctKLKgnzdcPHzanVleOdSYeO5KWUdUKIfwC/AA7Ae1LKXR15Tc2pYcpglWWqps6Io0HwwZ8pXD88kYiY8ZC7V03KGo1gsNE4In0zbP1IbQ+/HoIH2ua8Gs2xyNvfsH3kDyjNpt6vFw5AugwkObecvLIai/xHV6TDny2klD9KKftJKXtLKZ/u6OtpTi3OjgaeumgIAHd/fwRmfwfj7ob6amsBs/aSvqVh26wTotF0NGU5ED5cifMBlGZRkZsMQAaBHC2oACDMrxsbeY3987fhEUyP78nW1CKq6+pVWCVYa8+3l6KjamWtV6iST9BoTgWlWSrZ9r/2woSHoLqY2owdAPiF9rI069YjeY39I4RgWlwoNXVGzn5hJR/sqlMVxUePfeCJUJyqtHOC+kPBYdudV6NpSk05fH0jLH9KLe7zDAY3XwhUq1sd0taRL72Ije5pOeRYIn+djTbyGpsw3CRyll5Uyf82mOKGd3wJ2TaaginJAO9QZegba9prNLYmfTP89Rn88TzUVYGXmn8yp8D0ydtKugxkfP+GkO8QH9fO6Gmb0EZeYxP8PJwtWaZK8MAonGDvD/DWeJCy/Rcoz1XRO76RanRVW9X+c2o0LdHUzWiOGjOnwERNuo7pHYCvuxMxgR44OXRdU6pVKDU24+tbxlBeU8fYZ5eT6daHsIo9anFU0VHwizr5E0upJsA8ghr+0X6+XyVFvmE5CGGbG9BooHmyerORdw+wFPWI6IOTg4FV907A0MX//rruz4/mtMPN2YFATxeGRvrylNMdkHiNqvj+DvjsKsjdf+wTtEJyRhbUV5Mjvfk5TcXqs3khZGyBwmTbdF7TvTEa1RqMJXerQYmDS0NdSKx6d3CyFA0aOBhQa0Y8XLr2WFkbeY3NGRHtz0/ZPnwV+k+MLj5w+HfY8x1s+/i4x67Yl0NmcaVV2bKNKprh93TBU2vKrA8w59bUaNpD4RHY+C5sfAeO/gmhcZB0rRqotJCT1TXk9FmnoY28xuaMiFGPtf/6cge7XIc1VGRsaeUIxc70Yq5ZuJF/frZdFRiNYKwnM0NF6fxwsI4s2eQfLnm1zfqt6cY0nswvOAx+0TDtJbjgZQCMRkm9UfKnVGtCiBhxyrt4smgjr7E5QyN9Ldv/FnfA3O0w5FLrJOAtsP5IAQAHcsqUH/71UfCfcIKK1Ug+R/pQ12gaSUaNhe2fQP4h29+EpntRlGK9H95gxOvqjYx7bgW9//0jc6rv5v/O/B1cvU9t/9qBNvIam+Ph4sibf08kIcKXXdmVlLqFmUIf0xt0Z1Y+B4tmweJboUatGqxO3crLTq8RWJ8NeQcgbx/UVjC+agUAeVLJt77f+2VurrmThYaL1bmWPnzK71FjZ6RtBJdGhjtqjGUzq6SK9CLlQqzGmX4x7Qgi6AS69oyB5rRl6pAQvFwdmfXuel5bcZAHAsJVpE1ZNji5wopnlK+zIh8GXgD9z6V/5rdMdFhLTq0f+39YRT/TuQaKFIwYGNQnmjP7h3DB0Ek8tmsZP++RzPFwxWBWvtRoTobiNNjyIfQ7V02yZu+CYJXAbm9WCakFysA7OxgY3TuA+IiuoxXfFrSR13QYZrfNWysPM/uSIHqCCk+ryAckXPw2fHK5WnzS9xyESevmRsclkAJpMpDtxl6c77ABo5s/H17fMLr6+tYxXPL6Wn72n8V5uQugthKcuu6qQ00X5rs71Hvvs2HkjQBIKamqqWfqy6sszZbedSbRgR4tnaFLo901mg7D3dmRr29VhnlLiacqLE5VLhphgOix4N0TNryDfLoHZ9c1TKJWhJ/BL+O+YoNRRTE4GGuszj0s0o8R0f5sL/FSBcfx92s0LVJfq7SQ+p8PI26wFF/x9joGPvKzZT/Aw5moAPfO6GG70SN5TYcyNMKXIC8X/sh2ZRpQ/dVtuMhKGDgdnNyo9gjFpTgV83KS7Ihz6dE3CfcB5zPTty8/Vk2FTR8gZHMN+cRoP35P7ckDTsD3d8J5z0PIkFN4d5pTyqb31LqISY/bbgFc5l9QWw6xMyznzCurZv2RAiYOCCbQ04UePq4Mj/ZDdPFFT62hjbymQxFCMCLGn1UphVQLV2XgASY/AcA3dWNIMGYywKBWGdYOvgxGXQqAOzDjvHOg6tKGhVWNCPN1Y299T6oGnIvrkV/hp/vgmiVIKfnjQB4DQ70I9uq6miKaNlKcBsY6+OEutZ94DfjH2ObcR9U6i2y/YdQUVJBaUMGerFIAbp3Qh8Qov2MdfVqgjbymwxkR7c+SvzLZ4xxOguEgr4Q8w1z/GIora7n/6AhgBLc5LOYep8/xHXCm9cEGB5jxXovnjfR3R2Jgy+jXGFM9E1JWQ8FhdlUGMPu9DYyI8efzm0Z3/A1qOpaXBlvvZ2yxnZFPWQv+vRg5bxfQIKbn6mQgNuz0mmBtDe2T13Q4ZoXKK2v+zRnVL7Mwuw9Go2RLSiEAd03qx+v10xkl38fTt+3J3OMjfBECNqcUKlcNwPd3kpxfDsAGU9y95jSmprx5ma3mX4xGOPon9RFjmlUlRPji7Ggf5lGP5DUdTv8QNTlagSszRvXnwz9TOJhbxn9/3oujQXDjmb0I9nahXw/PEzqvj5sTId6uHMkvh9B4CBoAmdvIiG4wDNV19V0y76amjaRtVO+RY2Do3+GXf9suGU3uXqgspKTHcAAi/N24LDECAZw9MNg21+gCtOunSgjxvBBirxDiLyHEN0II30Z1DwghDgoh9gkhzml3TzWnLQ4GwR0T+zKhfxCzx0QD8MNfmezNKiU23Ac3ZwdmjogkMaq5RsjxCPN1I6Ookp0ZJdSMngtVxdRlqsduf0oo2/oVVOgR/WlLyloViXXlZzB0lso61lQl8mTJ2a3ePFQE1yPTBnPHxL7cPrEvg3vah6sG2u+u+RUYIqWMA/YDDwAIIQYBVwCDganA60IIPZzqxvxzcj8WXjOCXoEe9PB24fUVBwF49ILBxzny2IT5ubE9tZhp81bz351qxaJHntK++bfTJwQsuUEtvNKcnqSsVQuUzDICPpG2Sxpj+rHIFMpFGOjpbJvzdjHaZeSllEullKZcb6wDzKr6FwKLpJTVUsojwEHg9FH00XQYKtomgDqjxM3JgcE926cBEubrRmVtPQDv7zZSJw2UZB2mh7cLowxqpKYTf5+m1NUod01kI5+5T7hy19giEU1hCrj5kV2tvNaBni7HOeD0xJY++WuBz0zbYSijbybNVNYMIcSNwI0AkZGRNuyOpqty16S+9A7yYHBPn3Zn1GmcW7MeB7LwZ4phE4megvCaPOqFIw4Fh5VROE3jnO2a3d8pGerx91tyqJKyVslSD7pYpd+LahQh5RsBNaVQVQRu7QhvTFmrchLEjCevTC20C/LqpkZeCLEMCGmh6kEp5bemNg8CdcDxBcObIKV8G3gbICkpyQY/z5quTq8gT+6c1O/4DdvAqF4BDO7pzfBofw7klJKfG0V89WYoUI/0e4LOZUjO91CeB55tj9zRnCJWvwgZW1WC9jPvUWW/PwtHVqocwdBkJK/yrJL5F3iHNfwwnAhSwrf/UNvDrib3SDWeLo64OtmnR/m4Rl5KOelY9UKIOcA0YKKUlmeodCCiUbNwU5lGY1P6BHuy5I5xDQVrL4almwEoxZ31zqMZwveQf1Ab+a5IaZZ6bxwxU1fV8B7Qx/p78zWZlQ+nq/dHCsFwgk+DR/+EgkMQeznEziBvxxa79cdDO901QoipwL3AeCllRaOq74BPhBAvAj2BvsCG9lxLo2kTo26FvlNAOPDoDykcyC7hOlD/2FF6YVSXoq6mwcg3jpgpbjQeDI23PsYnwno/dw/0aMPkfVUxvDYCBpwHB39TZeeoCfm8smq7ddVA+33yrwEuwK8mXYd1UsqbpZS7hBCfA7tRbpzbpJT17byWRnN8DA7q0R8Y3MeBr/fuprZnP5xS1sK4f3Zy5zRWlKQBpod/c8RMfS2UZsCo29QIftBF1sd4NHkaS1nbNiOf8ieUZSn9G4BeEyxPCHllNfQNPrE1GqcT7TLyUspWHWJSyqeBp9tzfo2mPcSFq1jnXP9h9Dy4CKpKTquMPnaP2bCHxELeQeUrL80EaVQ/1Imzmx/TdPI8d2/brpXXJIn81P80VJVVM7pXwAl0/PTCPtbtajQtEGaKvDnopVY0sv6tTuyNphl5B9R79Dioq1R5BsxlfsfIvjTYlBEssJ9KQtMWTLkKAJUcJFgtgKquq6eootau3TXayGvslh7erjgaBGudxqpVk03zeGo6l6PrwCMYosaq/ed7w/9dor6rnsNaP+6Sd+D+o+AVonRsXomHp0LgmXA4vLLlYyrywasn3LnTSvAus0hN8jYOxbU3tJHX2C0OBsHAUG+2pRVBz6G2WympaT+l2bDjczUZ7ttkMrXxCteWcHACVx/wDIHM7Upjvq5Sxc+btW6aUp4P7gHqWs4q+UdRRQ03faQisXr62q8ktTbyGrtmeLQ/W48WUe8dbjvNE037Wflf9T7oQuV2iTRFPvn3hsQ5Vk1zSqrIKq5qfg6vHg3bnqbtpt9xfS0cWg4Fh8HTWnRs7aF89mWXEujpwuBQ+9GqaYpWodTYNUnRfry35ggFjsEEFf+iV752BWqrVJSLeyAMUQliuPbnVpuPeEaFPCY/e751hdmwe4fBP3fDW+ObK1Tu+xE+v1ptNwmhzSlRPxy/3DkOH3enk7uX0wBt5DV2jTkvZ64hmKC6KijPbTai05xiMrcBEiY/3moTo1Ey9ZU/6B3UENpYXl2Hh0sjk+UVqt7DEtW7bwTkNomiKclo2PYJt6rKKa3G0SDwc7ffhVCg3TUaO8ccYZMmA1WBrbTINSePOelHxMhWm+zLLmV/dhk/7cyylGUWV1o36jsFzn4Ypjyp9n0ilbumsXhZRX7DdvAgq8NzS6sJ9HTBYLDvJztt5DV2jY+bEx7ODhyqMWnVa79852Ne5erZo9Umh3ObZ4TKKa22LnD1hjPvBr9ote8TDrUVKn9Azl71KjeFTvaaADHjm53PnkMnzWh3jcauEUIQ5ufGvkrTP7OOsOl8yrLB0Q1cvFptUlRZ06wst6mRb4o5SqcoGd45W20PmAaB/eHqxS2eL9THfqNqzOiRvMbuCfN140CJAzi4qKXtms6lLFtFxhxjAryoorZZ2dxF23jwmx2tn9esa2M28ADJq8AjsMXm3WUkr428xu7p6etGenGVMixlOZ3dHU1p1jFdNQAllc2NPMDH64+RxNukWWRFVTG4N08rWVVbT15ZNaE+9rsIyow28hq7J8zPjaKKWqpcgqgpyjj+AZqOpSznuEa+qdRA4+3iVn4AcHKD4TeohO7X/dpQ7t5cl8Ycdx/mZ/9GXvvkNXaPOcLm9wwDMeIILYz3NKeS0izoNf6YTYora/Fzd8IglMGvrm0Qsc0oqsTHrZW49vNfUO9SgrMn1JSpePwmpBepSB17XulqRo/kNXZPuGm0liN9CRZFlFfXsWjDUYY+sZScH/8DH13cyT3sRlQVQ3Vxs5j1phRV1uDj5sSKu89i00OT6NejYZI2vbDyGEeaEEItkoIWr2U+R7ive9v7fpqijbzG7jGLT+VKH/xEGbs/vIuP16VQWFFL8IZn1bL3stxO7mU3wbxOoWnyj6bNKmrxcXPG3dkRL1cn5s8axtyJfQEoKG8eedMi5mibRjo41XX1PP79Lj7flIoQEKKjazSa059gL/WP/KdxEJnSn+HpH+JRuIvB4khDo7x9ndS7boY5hNU38pjNSiprrVwyPbxduW5cjDpFaz75pkz9r1K47DXBUrT2YD4L1ySzKaWQSH93nB3t3wTa/x1quj0OBkGkvzvbDQNZMfIdABbJ+1ji8mBDIx11c2owL0Y7hrtGSklBRQ1+TfRkvFwccTCIFmPoWySwD1zzo1V0zYbkAsv22D4th1baGzaZeBVC/At4AQiSUuYJlQvwFeA8oAKYI6XcYotraTQnw09zx1EvJW6OBnLXexMkSqwblJ+cu+ZQbhmXvfknsWE+fHDtCEt5bmk1l7yxhlevGMrQSL/2dN2+KEoBB2elI2+iuLKW819dRX5ZDS5OBl69YihVtUZCm2i8CyHwdnVs+0i+BTYlFzAkzJtnL4mjjx2n/GtMu0fyQogIYArQOID1XFTy7r7AjcAb7b2ORtMePFwc8XZ1wsnRgfRxz7El9G+kB4/nnzU3I4XDCY/k92aVcN+Xf/H0kj0UlNewcn8u+WUNKzKfXrKb1IJKZrz5J++uOszG5AK+3tLNV9tu+Qh2fqP04g3K9Czbnc3cRVtJK6zkvNhQiitrufq9DUBDVFRjfN2dW1wo1VaO5FUwpKcPQ8J8cHVyOOnznE7YYiT/EnAv8G2jsguBD6WUElgnhPAVQoRKKTNtcD2Npl0kTJoJk2ayO6OEr19dxZOOX+BRfmJG/r3VR/h8UxqeVDBcHGWT7Mee7es4I64/siSTVQeUZoq7swNP/7jHopk1La5nt/ADN6OqBL77h0r2MfRfANTVG7n9060IAYlRfjxzyRB2phezL7sUgIQI32an8XZzOqmRfE2dkYM5ZeSVVbf442HPtOuvTQhxIZAupdzepCoMaKwElWYqa+kcNwohNgkhNuXm6ggHzamjf4gXTg6C1BpP5AmO5LelFgHwkOP/8YXLE7zr/BJn/Dod/tcf8c5ZuJRn8N9LY3n7qiQrUcSdGcU2vIPTiFX/U+/TX4OkawFYvC2Dytp6bpvQh69uGYOLowPPXhoLwIzE8BYjX3xP0si/vGw/5726CoAIf/sPm2zMcUfyQohlQEgLVQ8C/0a5ak4aKeXbwNsASUlJ8jjNNRqb4WAQzJ3Yl5wV3gRmpdLWaTgpJWmFlfx9VCTTsiohCyYaNlu1GWHYy/DoSy3hm4PFET50fhbnjz3gjjWt6qnYLanrwTscBl5gKVp7SD3t3HhmL0vZ0Eg/fr3rzFYNsY+bE8n5zRUqj4f5R1ldw/eEjz+dOe5IXko5SUo5pOkLOAzEANuFEMlAOLBFCBECpAONA2HDTWUaTZdiztgYcvHBoTKvzccUV9ZSUVNPdIAHnvXWE7jL/S4DYHb/OmJ8DLj+9jAfTqrjncBFBIhSvKqz4EgryabtEaMRfnsC0jfDkIutRMk2HClg6uAQnByszVDfHl6t+st93Z1OyidvDsd89IJBRHazkfxJu2uklDuklMFSymgpZTTKJTNMSpkFfAdcLRSjgGLtj9d0RTxdHHHw6oFHbYF1soljsCdT+Yx7B3qokEAHFyrdw7it5g6uzbyYAodAhnqXIQ4th3XzOXPr3fQs3UGqQzhVwhVS/uzIW+paHFiqXDXugTCoYWVxWXUdaYWVxEWcWG7VAA8Xiitrqa6rP37jRmQUVTKubyDXjI1BdLP0jx01A/QjaqR/EHgHuLWDrqPRtBvpHoQztVBdcvzGqBGoEJDYQyhtlImPkHL1OpYYRwFQ5RkB+QfV6BXANKn7Reg97BT9IH1Th9xHl2Tls+r9xhUQnmgptsgK+J3YqNqsNZNZ1EJi72OQXlTZ7SZczdjMyJtG9HmmbSmlvE1K2VtKGSul7EZ/1ZrTDUdvpYhoLG3b5OvG5AIGhHjjXW3SpveNsDIghojhkLYBtn9qdZxfn5HsqwnCWHgMuVx7orYSsnZC7GXgZT2tl1ZYAbQcJnkszD8KV76zrs3HKFnhGm3kNZruiou/Wn1ZnH2k1TY5pVU89t0ulu/NZvXBPEbG+DdavRmBl6sTT140hHvO6U/wWbeo8tJGHkpXX4b1DiVdBmKozKe+qqyjbqfrkL4ZjLUw5NJmVdtSizAIFeF0IiRG+TEgxIuM4irLD8Xyvdm88Ms+6o3N3W1v/3GI+SsOAt1DVrgltJHXdHs8epg0UbJaN/Jfb0nn/bXJXPu+eiidNtAXtn6sKk1iW1eNiuK2CX0wBPZSmikAfc+BwH4w5UkG9/Qmz0W1Pbiz7SPR05Ztn6j3Jgm7S6pq+WT9UQb39MHT5cSW6jg7Gnjx8gQAFm1IJbWggps/2sJrKw7yxSbr/L05pVU88+Ne5i1XRr6nHslrNN2TgJ4xGKWgOi+l1TbbG4XgJUb5kZS6EPYtATe/lsMho8ep94SZ8I+NMOxqHB0M3HvjHAAKd9t5hE1VMWz7GPximmVmeuqH3eSX15y0dkz/EC9cnQy8tuIg455bQU29EYD7v95hlQd245FCq+O0u0aj6aaEBfiQjR/GotRW26SaXAMAY3sK+ON5tTz/plUt5yod9y9VN/BCq+LAkAhSRU9GHX4VaiqaH9fVOPArvDZc+ddPhFQlTcDEh5tVrT2UT/8eXtw5qe9JdcnBIKz05QGuHKlULR//fhc3fbSJEU8v49MN1nMf3UFWuCV0ZihNt8fL1YlkhyCciq2N/OcbU9mbVco/p/QjvbCSyYN6kORdwtXFL6sGw2Y3aJY3xdEZQuNarNruO4mIwg/V5Gyvs2x3Ix3BT/dBwSFI3wLRY9t+XMoaMDhCv3OtijOKKkkrrOSRaYPapR3z9EWxfLk5FQeDAV93J245qze/7s5m2Z5sqmrVyD6ntJoz+gSy+qBaA9E0Hr+7oI28RgOUe8bQv+RPFSsvBDV1Ru796i8A4sJ9KKyoJSHCl5vqvodtv0D4cEicc1LX2hN1FecVfIRIWYvo6kbeI1AZ+dw9J2bks3ZC0EBwtg6R3GiS+h0R0zy59okQG+5DbLh1jP1NZ/biqSV7rMqGR/uTFO1HRc2JxdXbE93zp02jaUJBQCJ+FGPM3Q/AjvQiS51Zgzzczw1S1ioDf/0ycGglz+hxCAgMYreMou7wmnb3u8OpM/m4T1RvvzgN/KKaFW89WoS7swMDQ71bOKh9hDeKnpkWFwrA2D4B3DmpH/8+b6DNr3e6oI28RgPUhKsIkO1rfgRgU3LDpN26w/kAhHsCGVshaky7rhXm58ZGY38cMjZBXRsTYJxqfv8v/HCXSroNUJbd9mONRhVe2kJikLTCCiL93XEw2H7VaZgpX6u3qyPPzYjjxzvGkRTdvicGe0AbeY0GOGPESHKlD712zQOjkcO55QR6ujA00pfDueUIIYnb/G8V991eI+/rxnrjQAz1VZDZVMC1g8j8C5Y+BNVtjM///RnY9B6UmYx8YQosfRh+ebD5OeprYdnj8Ofraj9vn1oJHBLb7LRphR238nRQT2/unNSXd2cPx93ZkUE9bf+0cDqiffIaDRDo5coWrwSGla1k9+pvWHXAlzA/Ny4dGoZv/nbODK7Cac9iEAaIHN2ua0UFuLNZ9lc7f30GtRXQc6hVwmmb8+sjcHgFRI6BAecdu21RCytyD69QL1BrABqfI2UtrH5RbSfOhuxdarvnsGanySiqbLc/vjUcDII7J/XrkHOfzuiRvEZjInm4CverXfYUGcVV+Lk78feggyysf4BrMh9XjW5d325j7OXqRFBIBEcdo2HjO/DhdFj64HGPaxfm9IbFrYeJWkjbaL3fdHK46TkytjaqS2v4kWiSrLu0qpaSqrpuG6/eWWgjr9GYmH5GIgXR5xNvOMxwsZd/FT0Nn10FBieY/T3cshaCbDNSHBHjz6zq+6i7+geIORO2fKhcIh1FrSkmvzAFPpgOuxZb12//DJ7rDc/1guVPW9f1HKrefSLBwaW5kW88KVuUCsufUhmgXKxzqKaZRMm668rTzkK7azQaE44OBvwveBLmLeF9rzfwKM5VromEK5UhtiEDQrx4v9aHLL9EwhNmwZE/1Ii4hYiUdiMlFJvSORxYCvkHlKb9gHxwMJmAQ8uVb93FU4VMuvrA5R+pOYjKItXGWKcmU5suGqvIB+EAsl49mch66H12s26YVw13RGSNpnX0SF6jaUxAb4gcjUeNyb0x4UEY+nebX8YslpVeWAn9pqrCtrhSTobyXKg3hULmH2goT22kn1OcBj0GwZg7lPb7gGnQazz0maTcNQF9YPRtavGXuZ+5+6E0GyryoMdgNV+x/2dVN+nxZt34eks6gZ7O9A7y6Jj71LSINvIaTVOmvdywHZbYarP2EGGSzE3OL1ejZlcfyD/UIdeyjLw9laQyTh6AUBOmoEb6BYeUD33UzXDvIbjo9YbjPQLh9s0w5h9KjK04DSoLYf5weP989SPi2QM8glX7uL81eyLZm1XChuQCYsN8ul3Sjs5GG3mNpinBA+Dm1XD7FnDqGL2TqAB3/NydVDy+EBCWpPKgdgTFponQwZeYLj5GjbxT1sIP/4THfZUscvjw45/LN1LFzB82CazlH4Ds3cqNM+wqVebf2/rylbVMfVkl0b7/3O67KKmz0D55jaYlWojxtiVCCJKi/S2raYkarSYsj66HyJHHPvhEKU5T76NvU6t04/4Gm9+HLR80hEVC2+L/zQucdnzRUGasVW6c0f8AJzcYebPVIa+b9NwvHRZ+wvrxmvajR/IaTScxItqflPwKckqqlA8cYPmTtr/Q/l+Un903AqY8CSFDlEGvN622dfJQmu9BbRhle/dU7weWQtCAhvKAPuDootQ3nRt87rX1RhasVjr9z1wyxFZ3pDkB2m3khRC3CyH2CiF2CSGea1T+gBDioBBinxDinPZeR6OxN5Ki/QDYcrQIggdC3ymQvEq5P2xF9m51TnMYpBlzUhOAO7bCdUvB0AZz4G7SgK+vgaD+DU88rSwQ+2T9UeqMkicuHIyL48mrTmpOnna5a4QQE4ALgXgpZbUQIthUPgi4AhgM9ASWCSH6SSm7rxScRtOEmEA14jWnsWPio2qEfOQPFeliC478od4nPWZd7tUDbl3XsN1WGidI8YmA8/6noms8g5s1ra6r59HvduHkIJg5IrJZvebU0N6R/C3As1LKagAppXlVxIXAIilltZTyCHAQGNHOa2k0doWPmxPuzg6kF5kScoQMUYbz6FrbXeToWnXOkBZcJcED1etEcGskSRASC55BrZ5je2oxANeOjem2Wu5dgfZ+8v2AcUKI9UKIlUII8/R8GNA46DfNVNYMIcSNQohNQohNubm57eyORnP6IIQgzNdNxcqbiRqjol5k86TUJ0RFAeTug6wdENZcQ+akcXRu2D6Ohs/uDGXkrzsjxnbX15wwxzXyQohlQoidLbwuRLl7/IFRwD3A5+IEg2CllG9LKZOklElBQUEndRMazelKmJ9bw0ge1ARoeW7LImEnwseXwfwRUHAYfG28itY84ep7bBdMRnEVLo4GgrxcbHt9zQlxXJ+8lHJSa3VCiFuAr6WUEtgghDACgUA60DgvWripTKPRNCLM142tR4uQUqpFQgGmGPOilJOXOKirhvRNDfs+raQoBIoqarj8rT/xdXdm0Q2jMLRF5/26X5UWznHGc+kmWWG9+Klzaa+7ZjEwAUAI0Q9wBvKA74ArhBAuQogYoC+woZ3X0mjsjoGh3hRX1nK0wDT5ajbIx0gqbkX+IaXlXtmQ5IT0LdZtWvLHm3ht+UH2Z5ex4UgBN360mS82Nb/ukbxy/vvzXgrLTSGXrt7gFXLcrqUVVWoxsi5Ae438e0AvIcROYBEwWyp2AZ8Du4Gfgdt0ZI1G0xyztvr6I6ZFUT4RSukxp41hlOteV1ruO75sKDNP3A6+BHoMaVHXHUBKyZIdmZb9ZXuyuefLv6g05UMtqaqluKKW+SsO8sbvh/hm64k9jKd3YIIQTdtpl5GXUtZIKf8upRwipRwmpVzeqO5pKWVvKWV/KeVP7e+qRmN/9AnyxM/diY1mI+/orPRyUtoYYVNdqt7TGrlnUv5UfvPLFsIta1qVZkgrrCSzuIonLhyMt2uD53anacJ0/HMrGP7MMvZmlQANSbjbQll1HXll1VZ5VzWdg45r0mg6EYOhibwBqAibzO1tS9Vnzr1alAIb3oF5iXDw1zZJFGww/bCMiPFnyR3jePfqJABSCyqorqunsKKWmjoj+7LUD8nG5AJkG6J+rn1/I+OfU3IJQyP9jn8Pmg5FG3mNppMZGaPkDbJLqlRB1GilyZ7WhmmsUrORT4UVz0C+0okhsnUjn19WzX9+2sOzP+/Fx82JfsFeRPi7c0ZftdDp++0ZjHrmN0v72npJTKAHeWU1HM4rP2Z30gorWL43h749PLnpzF4dlupP03a0kddoOpn4CF8AdmcqtwjhI5Q2e8qfyngfa/RsTrRdkgaVZr9+ZPOUfY34aksab608TG5pNcOj/SwRNa5ODgR6OrNiXy6FFbVWx1yepCaELW6lVjA/HTwybTAPnDcQZ0dtYjob/Q1oNJ1MeOMEIqCiV0LiYP1b8PIQ2Lyw5QPrqlVUTeMQyau+gbt2qJWorbDNlKEJaDbSDvRsiGm/eGgYvUwJPi5PCsffw5mtR4s4FttTi/B0cdRqk10ILTWs0XQywV6uOBqE9aKoqDGQuU1t/3AXJF7TPC7dnFs1LLEhW1N4y+ohmcWVXLVgAzeO68XRggrG9glg7sR+xEf4WLVzd1YiYg+dP5C/j4qiuLKWgvIaAjxdiA5wJ9Wss9MKuWXV9PB2waEt8faaU4I28hpNJ+NgEIT6upLR2MiHJ1k3yj8IgX2ty0pMIY19p8DuxWq7SfLslPxyPl6vVs8ezCnjPz/tASAu3LdFf/m/zxvIzzuz+PuoKFydHHB1cqCHt4rO6enrxo704mPeS0F5Df4ezsdsozm1aCOv0XQBevo00bDpfTYE9IWxd8B3t6uQSrORlxJKMhoWTIUlKkPfQvLsp5fsYenubMu+2dceE9ByntWkaH+SolueLO0V6MEPf2WSX1ZNgGfLUgWF5bVEB7of73Y1pxDtk9dougDNNGzc/OD2TTD0KvAIso6bX/MyvDQIUlarfd8ImPUFjLrF6pypBRVWBr4xo3oFnHAfzcf89+e9rbbJL6/WI/kuhjbyGk0XINzXjeySKmrrjdYVQii1x78WQYHKsMRfnze8u/lbZWJqzJ+H8gEl9QvQJ9iTRTeO4qtbxhAb7tPiMcdidO8AYgI9Wp18raqtJ6+shp4+egFUV0IbeY2mCxDm54ZRQlZxVfPKoaYE2b8+olw1Zp2a2go1im+FDckF+Hs4c885/bl2bAzPXBzLqF4BJEad3AIlIQQzEsM5kFNGgVnHphHmOQWtV9O10EZeo+kCmA2jlcvGTL8pasS+70cozVQvM/69Wz3nrowS4sJ9cHN24JELBtlkYZL5HGaJg3qjVDlqG/U9TEsZdCm0kddougBmIS+rydfGXPgaGOtg93dqP2Kkeo8e23J7IL2wggg/206CxoX74OxgsLhsPtlwlHHPrSC/rNrSdy1K1rXQ0TUaTRfAPJJPa83Im7Mw/fWZep/ylFoMZTb2TSitqqWkqs7mo2oXRwfC/Nws8fJ/7M+lus7I1FdWUV1bj0FAiE/LgmiazkEbeY2mC+Dq5ECvQA92pBe13MDdH4IGQoZJK96/N3i0HiHzV5qKZ++Ilac9fV1JL6zEaJRsMrltckurGRLmzUUJYTqfaxdDfxsaTRchIcKXneklrTeIMo3mHVyU0T8Gu0xywUNNuji2JCbQg4M5ZezPKbXSuJk7sR/Xj+tl8+tp2oc28hpNFyHcz42c0hbCKM1Ej1PvvpHHTb2XUVSFl4sjvu62j1kfHu1PWXUdH/6ZAsCwSF8Akk4yakfTsWh3jUbTRTCHUWYWVREZ0MKE6aAL4arF4B9z3HOlFVZ0WCjjcNOK2M82phLs5cInN4wirbASP70IqkuijbxG00UYEOINwMK1RzhncEjzVakGB+g94bjnyS+rZtmeHC5M6NkR3aSnrxv/mtyPI3nljO8fhKuTA32CPY9/oKZTaJeRF0IkAG8CrkAdcKuUcoNQ6dlfAc4DKoA5UsotrZ5Io9EwuKc37s4OLFyTzMI1ySQ/e/4JHS+l5GBOGV9uSQNg8qAeHdFNAG6f2Pf4jTRdgvb65J8DHpdSJgCPmPYBzgX6ml43Am+08zoajd3j6GCwWrBkyRTVRv44kMfkl/7grZWH6R3kwbS4jhnJa04v2mvkJeBt2vYBMkzbFwIfSsU6wFcIEdrOa2k0ds9zM+K4a1I/QClItoXskirOe2UVdy7aaimLD/ftiO5pTkPa65O/E/hFCPEC6gfDnFgyDEht1C7NVJZJE4QQN6JG+0RGRrazOxrN6U2wlyv/OLsPH/6ZzHfbM5g9Jvq4WjP/W7qP3ZklnDskhD7BnuSV1TB7TNQp6rGmq3NcIy+EWAaEtFD1IDARuEtK+ZUQ4nJgATDpRDogpXwbeBsgKSnp+KngNRo7x8EgePriWG7+v8088cNuvr2tdemCeqPkxx1ZODsYeH3WMMRxQis13Y/jGnkpZatGWwjxITDXtPsF8K5pOx1oLI8XbirTaDRtYOqQEG45qzfv/HGYmjpjqwmx92SWUFZdxytXJGgDr2mR9vrkM4Dxpu2zgQOm7e+Aq4ViFFAspWzmqtFoNK0TE+hBnVG2LD9swqwGObyVbE4aTXt98jcArwghHIEqTL514EdU+ORBVAjlNe28jkbT7TCrOaYVVbS8OArYcKSAcD83reGuaZV2GXkp5WogsYVyCdzWnnNrNN2dUJOaY2ZRyyP59KJKNiYXcGbfoFPZLc1phl7xqtF0UYK9lZHPLatuVpdVXMX451ZQZ5Qnla9V033QAmUaTRfF08URd2cHckubG/n1R/KpM0r+c0ksFw8L64TeaU4XtJHXaLowwV4u5LRg5A/llGEQMCMxXOu3a46J/uvQaLowwV6ulhyqjUkrqqSHt6s28Jrjov9CNJouTJCXS4s++fTCSp1LVdMmtJHXaLowQV4u5JY0N/IZxZU2z9+qsU+0kddoujChPq6UVtfxV1oRK/blAErKILOoSsfGa9qEDqHUaLowZkM+/bU1APw0dxzZJVXUGSUxAR6d2TXNaYI28hpNFya8iUvm3FdWWbaHx2gpA83x0UZeo+nCDAz1tmw/PG0QIaYFUgGezsQE6pG85vhoI6/RdGFcnRx46PyBVNXWc90Zx0/grdE0RRt5jaaLc/24Xp3dBc1pjI6u0Wg0GjtGG3mNRqOxY7SR12g0GjtGG3mNRqOxY7SR12g0GjtGG3mNRqOxY7SR12g0GjtGG3mNRqOxY4TKud01EELkAikneXggkGfD7pyOdPfPQN+/vv/uev9RUsoWM7p3KSPfHoQQm6SUSZ3dj86ku38G+v71/Xfn+28N7a7RaDQaO0YbeY1Go7Fj7MnIv93ZHegCdPfPQN9/96a733+L2I1PXqPRaDTNsaeRvEaj0WiaoI28RqPR2DF2YeSFEFOFEPuEEAeFEPd3dn86AiFEhBBihRBitxBilxBirqncXwjxqxDigOndz1QuhBCvmj6Tv4QQwzr3DmyDEMJBCLFVCPGDaT9GCLHedJ+fCSGcTeUupv2DpvroTu24DRBC+AohvhRC7BVC7BFCjO5O378Q4i7T3/5OIcSnQgjX7vT9nyynvZEXQjgA84FzgUHATCHEoM7tVYdQB/xLSjkIGAXcZrrP+4HfpJR9gd9M+6A+j76m143AG6e+yx3CXGBPo/3/Ai9JKfsAhcB1pvLrgEJT+Uumdqc7rwA/SykHAPGoz6FbfP9CiDDgDiBJSjkEcACuoHt9/yeHlPK0fgGjgV8a7T8APNDZ/ToF9/0tMBnYB4SaykKBfabtt4CZjdpb2p2uLyAcZcjOBn4ABGqFo2PTvwXgF2C0advR1E509j204959gCNN76G7fP9AGJAK+Ju+zx+Ac7rL99+e12k/kqfhyzeTZiqzW0yPnkOB9UAPKWWmqSoL6GHatsfP5WXgXsBo2g8AiqSUdab9xvdouX9TfbGp/elKDJALLDS5q94VQnjQTb5/KWU68AJwFMhEfZ+b6T7f/0ljD0a+WyGE8AS+Au6UUpY0rpNq2GKXMbFCiGlAjpRyc2f3pZNwBIYBb0gphwLlNLhmALv//v2AC1E/dj0BD2Bqp3bqNMEejHw6ENFoP9xUZncIIZxQBv5jKeXXpuJsIUSoqT4UyDGV29vnMhaYLoRIBhahXDavAL5CCEdTm8b3aLl/U70PkH8qO2xj0oA0KeV60/6XKKPfXb7/ScARKWWulLIW+Br1N9Fdvv+Txh6M/Eagr2mW3Rk1GfNdJ/fJ5gghBLAA2COlfLFR1XfAbNP2bJSv3lx+tSnKYhRQ3Oix/rRDSvmAlDJcShmN+o6XSylnASuAGaZmTe/f/LnMMLU/bUe5UsosIFUI0d9UNBHYTTf5/lFumlFCCHfT/4L5/rvF998uOntSwBYv4DxgP3AIeLCz+9NB93gG6lH8L2Cb6XUeys/4G3AAWAb4m9oLVNTRIWAHKiqh0+/DRp/FWcAPpu1ewAbgIPAF4GIqdzXtHzTV9+rsftvgvhOATaa/gcWAX3f6/oHHgb3ATuAjwKU7ff8n+9KyBhqNRmPH2IO7RqPRaDStoI28RqPR2DHayGs0Go0do428RqPR2DHayGs0Go0do428RqPR2DHayGs0Go0d8/8r/Y8jTCBeMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# example of how to load and use the traces\n",
    "\n",
    "import numpy as np, tqdm, matplotlib.pyplot as plt\n",
    "from h5py import File\n",
    "\n",
    "with File(\"./data/data_1.h5\",\"r\") as f_data:\n",
    "    key = np.array(f_data['key'], dtype='int16') # secret key to be found\n",
    "    traces = np.array(f_data[\"traces\"], dtype='int8') # EM traces\n",
    "\n",
    "print(\"shape of key: \",key.shape)\n",
    "print(\"shape of trace dataset: \",traces.shape)\n",
    "\n",
    "plt.plot(traces[0,9], label=\"10th EM peak of 1st execution\")\n",
    "plt.plot(traces[2,19], label=\"20th EM peak of 3rd execution\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Attack script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np, tqdm, matplotlib.pyplot as plt\n",
    "from h5py import File\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "def cyclic_dot_product(a,b):\n",
    "    assert(len(a)==len(b))\n",
    "    return [np.dot(a,np.roll(b,i)) for i in range(len(b))]\n",
    "\n",
    "def attack(key_numbers, n_at):\n",
    "\n",
    "    for key_number in key_numbers:\n",
    "        # time window for locating the -1's, from first leakage assessment\n",
    "        targeted_window_1 = np.arange(195,250)\n",
    "        window_size_1 = len(targeted_window_1)\n",
    "\n",
    "        # time window for locating the 0's, from second leakage assessment\n",
    "        targeted_window_2 = np.arange(545,580)\n",
    "        window_size_2 = len(targeted_window_2)\n",
    "        \n",
    "        # get the traces\n",
    "        with File(\"./data/data_{}.h5\".format(key_number),\"r\") as f_data:\n",
    "            true_key = np.array(f_data['key'], dtype='int16')\n",
    "            traces = np.array(f_data[\"traces\"], dtype='int8')\n",
    "\n",
    "        # do you want to randomly permute the traces and start over?\n",
    "        # then change the next condition\n",
    "        if False:\n",
    "            traces = traces[np.random.permutation(traces.shape[0])]\n",
    "\n",
    "        # true key used for verification purpose only\n",
    "        true_key.flags.writeable = False\n",
    "        \n",
    "        for n_traces in n_at:\n",
    "\n",
    "            ########################################################\n",
    "            ##### First phase of the attack: locating the -1's #####\n",
    "            ########################################################\n",
    "\n",
    "            # getting ready for clustering\n",
    "            # gathering all the peaks to perform the clustering\n",
    "            X = traces[:n_traces,:,targeted_window_1].reshape((N*n_traces,window_size_1))\n",
    "\n",
    "            y = KMeans(n_clusters=2,random_state=42).fit(X).labels_\n",
    "            \n",
    "            # for a 2-clustering, set label 1 to smallest cluster, that should be the -1's\n",
    "            bcount = np.bincount(y)\n",
    "            if bcount[0] < bcount[1]:\n",
    "                    y = 1-y\n",
    "            \n",
    "            y = y.reshape(n_traces,N)\n",
    "\n",
    "            # choosing a reference vector of labels (minimizing CSR)\n",
    "            ref_trace_index = 0\n",
    "            best_ratio = 1\n",
    "            for i in range(n_traces):\n",
    "                bc = np.bincount(y[i])\n",
    "                ratio = np.abs(np.min(bc)/np.max(bc)-0.5)\n",
    "                if ratio < best_ratio:\n",
    "                    best_ratio = ratio\n",
    "                    ref_trace_index = i\n",
    "\n",
    "            # add two zero's for the missing first and last peaks\n",
    "            ref_labels = np.concatenate((np.zeros(1,dtype='int'), y[ref_trace_index], np.zeros(1,dtype='int')))\n",
    "\n",
    "            # accumulator to apply majority vote\n",
    "            acc_maj_vote = np.zeros((NTRU_N,2), dtype='int')\n",
    "\n",
    "            kept_indices = [] # traces that we do not discard\n",
    "            rotation_indices = [] # used to store the guessed rotation indices\n",
    "\n",
    "            for i in range(n_traces):\n",
    "                # for each vector of labels, find the best rotation that matches ref_labels\n",
    "\n",
    "                # add two zero's for the missing first and last peaks\n",
    "                tmp_y = np.concatenate((np.zeros(1,dtype='int'), y[i], np.zeros(1,dtype='int')))\n",
    "\n",
    "                # perform the rotations and dot products\n",
    "                cdp = cyclic_dot_product(ref_labels, tmp_y)\n",
    "\n",
    "                # get the two maximal values\n",
    "                maxv = np.sort(cdp)[-2:]\n",
    "\n",
    "                # discard if the two maximal values are too close (no peak)\n",
    "                # (threshold that can be changed)\n",
    "                if maxv[1]-maxv[0] >= 10:\n",
    "                    # keep this trace\n",
    "                    kept_indices.append(i)\n",
    "\n",
    "                    # associated assumed rotation index\n",
    "                    rotation_indices.append(np.argmax(cdp))\n",
    "\n",
    "                    # rotate the vector of labels\n",
    "                    shifted_tmp_y = np.roll(tmp_y, rotation_indices[-1])\n",
    "\n",
    "                    # update majority vote with shifted labels\n",
    "                    for j in range(NTRU_N):\n",
    "                        acc_maj_vote[j, shifted_tmp_y[j]] += 1\n",
    "            \n",
    "            # majority vote to discriminate -1 vs 0/1\n",
    "            guessed_key_m1 = np.argmax(acc_maj_vote, axis=1)\n",
    "\n",
    "            # indices where the -1 should be\n",
    "            ind_m1 = np.where(guessed_key_m1 == 1)[0]\n",
    "\n",
    "            #######################################################\n",
    "            ##### second part of the attack: locating the 0's #####\n",
    "            #######################################################\n",
    "\n",
    "            # indices that are assumed to point to 0 and 1 coefficients: [0,508] - ind_m1\n",
    "            ind_01 = np.setdiff1d(np.arange(NTRU_N), ind_m1, assume_unique=True)\n",
    "\n",
    "            # getting the L1 distance of consecutive peaks on the area of interest\n",
    "            l1_dist_traces = {}\n",
    "\n",
    "            for c in ind_01:\n",
    "                l1_dist = []\n",
    "                for i, ki in enumerate(kept_indices):\n",
    "                    current_place = (c-rotation_indices[i])%NTRU_N\n",
    "                    if 507 >= current_place >= 2: # 3 positions for which we don't have diff traces: (0,1,508) / 0 and 508 by choice (difformed peaks), 1 because not possible        \n",
    "                        current_l1_dist = np.sum(np.abs(traces[ki, current_place-1, targeted_window_2]-traces[ki, current_place-2, targeted_window_2]))\n",
    "                        l1_dist.append(current_l1_dist)\n",
    "\n",
    "                # average over the distances for c-th coeff\n",
    "                # hopefully, one average trace per coefficient\n",
    "                if len(l1_dist) > 0:\n",
    "                    l1_dist_traces[c] = np.mean(l1_dist)\n",
    "\n",
    "            # we're gonna classify according to the average of each trace compared to the global average\n",
    "            # the classification threshold is the following one:\n",
    "            mean_l1_dist = np.mean(list(l1_dist_traces.values()))\n",
    "            \n",
    "            # calculating the labels\n",
    "            labels = []\n",
    "            for c in ind_01:\n",
    "                \n",
    "                # if traces are available for c\n",
    "                if c in l1_dist_traces.keys():\n",
    "                    m = l1_dist_traces[c]\n",
    "                    if m < mean_l1_dist:\n",
    "                        labels.append(0)\n",
    "                    else:\n",
    "                        labels.append(1)\n",
    "                        \n",
    "                # otherwise choose randomly\n",
    "                else:\n",
    "                    labels.append(np.random.choice([0,1]))\n",
    "\n",
    "            labels = np.array(labels)\n",
    "\n",
    "            # indices where the 0's and the 1's should be located\n",
    "            ind_0 = ind_01[np.where(labels==0)[0]]\n",
    "            ind_1 = ind_01[np.where(labels==1)[0]]\n",
    "\n",
    "            # building our key candidate\n",
    "            final_key = np.zeros(NTRU_N, dtype='int')\n",
    "            final_key[ind_m1] = 2**11-1\n",
    "            final_key[ind_0] = 0\n",
    "            final_key[ind_1] = 1\n",
    "\n",
    "            #######################\n",
    "            ##### Final check #####\n",
    "            #######################\n",
    "\n",
    "            # we might get the right key up to a final global rotation\n",
    "            # that can be brute forced\n",
    "            global_rot_ind = np.argmax(cyclic_dot_product(true_key, final_key))\n",
    "            nb_matching_coeff = np.sum(true_key == np.roll(final_key, global_rot_ind))\n",
    "\n",
    "            # (brute force correction of wrong coefficients is not implemented)\n",
    "            \n",
    "            if nb_matching_coeff == NTRU_N:\n",
    "                print(\"Key {} found with {} traces.\".format(key_number, n_traces))\n",
    "                break # done!\n",
    "    \n",
    "    print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Key 1 found with 40 traces.\n",
      "Key 2 found with 23 traces.\n",
      "Key 3 found with 57 traces.\n",
      "Key 4 found with 60 traces.\n",
      "Key 5 found with 38 traces.\n",
      "Key 6 found with 28 traces.\n",
      "Key 7 found with 50 traces.\n",
      "Key 8 found with 30 traces.\n",
      "Key 9 found with 49 traces.\n",
      "Key 10 found with 52 traces.\n",
      "Key 11 found with 68 traces.\n",
      "Key 12 found with 28 traces.\n",
      "Key 13 found with 24 traces.\n",
      "Key 14 found with 34 traces.\n",
      "Key 15 found with 37 traces.\n",
      "Key 16 found with 49 traces.\n",
      "Key 17 found with 64 traces.\n",
      "Key 18 found with 41 traces.\n",
      "Key 19 found with 54 traces.\n",
      "Key 20 found with 49 traces.\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "# next parameters should not be modified\n",
    "NTRU_N = 509\n",
    "NTRU_LOGQ = 11\n",
    "NTRU_Q = (1 << NTRU_LOGQ)\n",
    "\n",
    "N = 507 # number of coeff/peaks extracted (do not change)\n",
    "dim_traces = 950\n",
    "\n",
    "# you can tweak the next two parameters\n",
    "Sk = np.arange(1,21) # set of key indices to test, number of keys up to 200, for full set, write: np.arange(1,201)\n",
    "n_at = np.arange(20,70,1) # number of attack traces to try, up to 100 per key (np.arange(2,100))\n",
    "\n",
    "# start attack\n",
    "attack(Sk, n_at)"
   ]
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}