Random Forest experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This file recreates the random forest experiments the anonymous NeurIPS submission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "import random\n",
    "import PIL\n",
    "import glob\n",
    "import time\n",
    "from collections import Counter\n",
    "from sklearn.ensemble import RandomForestClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def reseed(seed=0):\n",
    "    random.seed(seed)\n",
    "    np.random.seed(seed)\n",
    "    \n",
    "def rand(maxi=10**9):\n",
    "    return np.random.randint(maxi)\n",
    "\n",
    "def flatten(X):\n",
    "    return X.reshape(X.shape[0], -1)\n",
    "\n",
    "def viz(images, num_images=30, cols = 10, seed=None, shuffle=True):    \n",
    "    if seed:\n",
    "        reseed(seed)        \n",
    "    m = min(len(images), num_images)\n",
    "    if shuffle:\n",
    "        pi = np.random.permutation(len(images))[:m]\n",
    "    else:\n",
    "        pi = list(range(m))\n",
    "    n = (num_images + cols-1)//cols    \n",
    "    fig = plt.figure(figsize=(10,(10*n)//10))\n",
    "    for j, i in enumerate(pi):\n",
    "      plt.subplot(n,cols,j+1)\n",
    "      plt.tight_layout()\n",
    "      plt.imshow(images[i], cmap=\"gray\")\n",
    "      plt.xticks([])\n",
    "      plt.yticks([])\n",
    "    \n",
    "def Classifier():\n",
    "    return RandomForestClassifier(n_estimators=100, random_state=rand()) # n_estimators is now default 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data and take most common letters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['balanced', 'byclass', 'bymerge', 'digits', 'letters', 'mnist']"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import emnist \n",
    "emnist.list_datasets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((814255, 28, 28), (814255,))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emnist_map = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'\n",
    "emnist_inv_map = {c: i for i, c in enumerate(emnist_map)}\n",
    "emnist_X_train, emnist_y_train = emnist.extract_training_samples('byclass')\n",
    "emnist_X_test, emnist_y_test = emnist.extract_test_samples('byclass')\n",
    "emnist_X = np.vstack([emnist_X_train, emnist_X_test])\n",
    "emnist_y = np.hstack([emnist_y_train, emnist_y_test])\n",
    "emnist_X.shape, emnist_y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "lower_counts = Counter([emnist_map[i] for y in (emnist_y_train, emnist_y_test) for i in y if emnist_map[i].islower()])\n",
    "upper_counts = Counter([emnist_map[i] for y in (emnist_y_train, emnist_y_test) for i in y if emnist_map[i].isupper()])\n",
    "# plt.scatter([lower_counts[c] for c in sorted(options)], [upper_counts[c.upper()] for c in sorted(options)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(28723, 'e', 5785, 'E'),\n",
       " (21227, 't', 11396, 'T'),\n",
       " (17853, 'l', 5886, 'L'),\n",
       " (16425, 'r', 5882, 'R'),\n",
       " (13316, 'n', 9588, 'N'),\n",
       " (11860, 'd', 5341, 'D'),\n",
       " (11677, 'a', 7469, 'A'),\n",
       " (10217, 'h', 3673, 'H'),\n",
       " (6012, 'b', 4526, 'B'),\n",
       " (4276, 'g', 2964, 'G'),\n",
       " (3499, 'q', 3018, 'Q'),\n",
       " (3378, 'v', 5433, 'V'),\n",
       " (3312, 'u', 14604, 'U'),\n",
       " (3292, 'x', 3203, 'X'),\n",
       " (3286, 'c', 11833, 'C'),\n",
       " (3215, 'o', 29139, 'O'),\n",
       " (3176, 'z', 3165, 'Z'),\n",
       " (3164, 'w', 5501, 'W'),\n",
       " (3152, 'i', 13994, 'I'),\n",
       " (3136, 's', 24272, 'S'),\n",
       " (3109, 'm', 10487, 'M'),\n",
       " (2961, 'f', 10622, 'F'),\n",
       " (2957, 'k', 2850, 'K'),\n",
       " (2816, 'p', 9744, 'P'),\n",
       " (2746, 'y', 5541, 'Y'),\n",
       " (2213, 'j', 4388, 'J')]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted([(l, c, upper_counts[c.upper()], c.upper()) for c, l in lower_counts.most_common()], key = lambda z: (-z[0],-z[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('etlrndah', 'ETLRNDAH', 8, (10217, 3673))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_per_class = 10*1000, 3*1000\n",
    "chars = \"\".join([c for c, l in lower_counts.most_common() if l>=N_per_class[0] and upper_counts[c.upper()]>N_per_class[1]])\n",
    "CHARS = chars.upper()\n",
    "min_counts = min(lower_counts[c] for c in chars), min(upper_counts[c] for c in CHARS)\n",
    "chars, CHARS, len(chars), min_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract(train, test, chars, num_per_class=None):\n",
    "    ids = {emnist_inv_map[c] for c in chars}\n",
    "    ans = [(x, chars.index(emnist_map[y])) for x, y in zip(train, test) if y in ids]\n",
    "    if num_per_class is not None:\n",
    "        random.shuffle(ans)\n",
    "        counts = {}\n",
    "        def helper(z):\n",
    "            x, y = z\n",
    "            counts[y] = counts.setdefault(y, 0) + 1\n",
    "            return counts[y]<=num_per_class\n",
    "        ans = list(filter(helper, ans))\n",
    "    return np.array([x for x, y in ans]), np.array([y for x, y in ans])            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "reseed(0)\n",
    "\n",
    "# P is lowercase letters\n",
    "PX_all, Py_all = extract(emnist_X, emnist_y, chars, num_per_class=N_per_class[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save some P examples for later use in Q\n",
    "n_donate = 3000*len(chars) # donate from P to Q\n",
    "assert n_donate < len(PX_all)\n",
    "\n",
    "PX, PX_donate, Py, Py_donate = train_test_split(PX_all, Py_all, test_size=n_donate, shuffle=True, stratify=Py_all, random_state=rand())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error rates: P_train, P_test\n",
      "[0.0, 0.03642857142857143]\n"
     ]
    }
   ],
   "source": [
    "PX_train, PX_test, Py_train, Py_test = train_test_split(PX, Py, test_size=0.1, shuffle=True, random_state=rand())\n",
    "clf = Classifier()\n",
    "clf.fit(flatten(PX_train), Py_train)\n",
    "print(\"error rates: P_train, P_test\")\n",
    "print([np.mean(clf.predict(flatten(X))!=y) for X, y in [(PX_train, Py_train), (PX_test, Py_test)]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def go(QX, Qy, clf, plot_max_p_rej=0.5, fig_filename=None, fig_title=None, y_ticks=True, legend=False):\n",
    "    QX_train, QX_test, Qy_train, Qy_test = train_test_split(QX, Qy, test_size=0.1, shuffle=True, random_state=rand())\n",
    "\n",
    "    dist = Classifier()\n",
    "    print(f\"Training distinguisher on {len(PX_train):,} P vs {len(QX_train):,} Q\")\n",
    "    dist.fit(flatten(np.vstack([PX_train] + [QX_train])), [0 for _ in PX_train]+[1 for _ in QX_train])\n",
    "    print(dist)\n",
    "    print(\"Distinguisher prob of test: [PX_train, QX_train, PX_test, QX_test]\")\n",
    "    print([np.mean(dist.predict(flatten(X))) for X in [PX_train, QX_train, PX_test, QX_test]])\n",
    "    print(\"Distinguisher mean probabilities: [PX_train, QX_train, PX_test, QX_test]\")\n",
    "    print([np.mean(dist.predict_proba(flatten(X))[:, 1]) for X in [PX_train, QX_train, PX_test, QX_test]])\n",
    "    print()\n",
    "    errors = (clf.predict(flatten(QX_test))!=Qy_test)\n",
    "    # sel = (dist.predict(flatten(QX_test))==0)\n",
    "    # P_rej_rate = np.mean(dist.predict(flatten(PX_test))==1)\n",
    "    # print(1-np.mean(sel), P_rej_rate, np.mean(errors), np.mean(errors[sel]), np.sum(errors[sel])/len(sel))\n",
    "    Pdist_hat = dist.predict_proba(flatten(PX_test))[:, 1]\n",
    "    Qdist_hat = dist.predict_proba(flatten(QX_test))[:, 1]\n",
    "\n",
    "    thresholds = np.linspace(0.0, 1, 100)\n",
    "    p_rejs = [np.mean(Pdist_hat>tau) for tau in thresholds]\n",
    "    thresholds = [t for t, r in zip(thresholds, p_rejs) if r<plot_max_p_rej]\n",
    "\n",
    "    p_rejs = [np.mean(Pdist_hat>tau) for tau in thresholds]\n",
    "    \n",
    "    rejs = [np.mean(Qdist_hat>tau) for tau in thresholds]\n",
    "#     errs = [np.sum(errors[Qdist_hat<=tau])/len(Qdist_hat) for tau in thresholds]\n",
    "    normalized_errs = [np.mean(errors[Qdist_hat<=tau]) for tau in thresholds]\n",
    "#     plt.plot(rejs, errs, label=\"Unnormalized error\")\n",
    "    plt.figure(figsize=(2.8,2))\n",
    "    plt.ylim((-0.03,0.55))\n",
    "    if not y_ticks:\n",
    "        ax = plt.gca()\n",
    "        ax.set_yticklabels([])\n",
    "    if fig_title:\n",
    "        plt.title(fig_title)\n",
    "    plt.plot(*zip(*[(r, e) for r, e in zip(rejs, normalized_errs) if not np.isnan(e)]), label=r'Err on $Q$')\n",
    "    plt.plot(rejs, p_rejs, label=r'Rej on $P$')    \n",
    "#     tots = [np.sum(errors[Qdist_hat<=tau])/len(Qdist_hat) + np.mean(Pdist_hat>tau) for tau in thresholds]\n",
    "#     plt.plot(rejs, tots, label=\"$err_P+rej_Q$\")\n",
    "    plt.xlabel(r'Rejection rate on $Q$')  \n",
    "    # plt.ylabel(\"Error rate\")\n",
    "    if legend:\n",
    "        plt.legend()#loc=2, prop={'size': 6})\n",
    "    if fig_filename:\n",
    "        plt.savefig(fig_filename, bbox_inches = \"tight\")\n",
    "        \n",
    "    return clf, dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Below are errors (1 per class) on held out letters:  ['a', 'd', 'e', 'h', 'l', 'n', 'r', 't']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAA5CAYAAADUdDpEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGt9JREFUeJztnXlQ1OcZxz+/3cVdlWuX1YiooIIHhxK1naICxqhBMRqctEZbZ5KZpJHqjLW1qa3JJLFJOk2cqDOmtqk52kYdnfFChdQjAx5VKVEQRcFY5BKR+3ZhYfvHr783i5oAXrDr+/nHLBL8/R7e43mf5/s8r+JwOJBIJBKJRCJxN3Q9/QASiUQikUgkDwPp5EgkEolEInFLpJMjkUgkEonELZFOjkQikUgkErdEOjkSiUQikUjcEunkSCQSiUQicUukkyORSCQSicQtkU6ORCKRSCQSt0Q6ORKJRCKRSNwSQ3e+WVGUXtce2eFwKD39DM5IG3WOtFHnSBt1DWmnzpE26hxpo85xVRvJSI5EIpFIJBK3pFuRHImkt6DT6fDz8wPAx8eHmpoa6uvrAbDZbD35aI8ERVHE+3t6elJXVyfev7W1tScfrVeh1+tpa2vr6ceQSCQ9RK92cvR6PQ6Hg/b29p5+lF5NUFCQ+O+SkhK33+R0Oh2xsbEkJiYCEBYWRnp6OqmpqQCkpKRQXl6OO18+26dPH6KiogAYMWIEmZmZZGdnA1BVVdWTj9bjKIoawQ4MDGT+/Pns27cPgGvXrvXgU7k/Hh4eAAwePJjGxkYqKysB3HoeSno/vdLJMRqNAPj7+9PQ0PBYndC7y8CBA/n0009paWkBYOnSpW6/mFutVubNm8e8efMAMBgMBAUFMXv2bACmTJnC73//eyoqKnryMR8qAQEBvPHGGwCMHDmSDRs2UFJSAkgn54UXXgBg7dq1jBgxgkWLFgHquJBRnYfHwIEDAVi3bh11dXW89957APz3v/+Vjo6kx5CaHIlEIpFIJG5Jr4rkeHh44OPjQ3x8PACJiYmcO3eOAwcOAGoaQqauVPR6PQALFy5k8uTJFBQUAGpUw50xGAwsXLiQRYsWidNhSkoKXl5ehIWFARAdHY2vr6/bRnJ0Oh3jx48Xacr+/fuLFM3jjsFgYOXKlYCaxnM4HBQXFwMybfIw0ev1eHl5AeDn50dcXJwYk6+99prbzsX7wVlX2K9fP7GGSx4svWJH1Dbs0NBQ5s6dy/LlywF1shQUFIhcr+RbtAXEYrHQp0+fHn6ah4/2vkajkaCgIMxmMzk5OQC8/fbbzJs3D39/f0BdPNwZnU5HUFAQ/fv3B6CtrY36+nrsdnsPP1nPYzQasVqtgDpmKioq2Lp1K4A8ID0kdDod06ZNE+v2j370I0wmEzExMQAMGzZMOjn/Z8SIEQAkJCQQFRVFdHQ0AHa7ncjISMrLy3vy8R4qJpOJ0NBQYmNjAcjMzCQtLe2hz8sed3L69evHzJkzAXjppZcYN24cAwYMACA/P59t27Zx6tQpwP0XqQdRCaI5jO6Gs04rODgYu91OZmYmAKWlpR1O6e5+Yvfw8GDIkCHCmSstLeXYsWOUlZV1+D6dTve9Dp+7OUU6nY7p06cLZxegpqZGOMOSe+f2CLHz2LFarbzyyitMnToV+HYN0n4PsbGxnD9/3mXHm8FgQFEUfH19AbWa05na2lpqamq+t+DDarUybNgwoaObMWMGffv2FfPTbrczdOhQt3NytHV79OjRJCQksGTJEjEudu7cSXZ2NrW1tYCqM/y+TITdbqewsLDbfoB7H3klEolEIpE8tvRoJMdoNBISEsJzzz0HqGFOHx8f4RGfPn2arKws4em5C5r3rv05ZMgQQPXuk5KSaGhoEJUKBoOBxsZGYZP6+nra29tFtKKurg5AnDJCQ0PJzc19dC/zCNDr9TzzzDMAPPfcc8yYMYOysjI++ugjQE3XhIWFiRNWYWGhy54avw9tvHh6ehIUFCRSeC0tLR3SVTqdDqPRiKen5x2nTg273c6NGzdEVZ47REkVRWHkyJHi9Ajqe8mKqntHr9eLUnxPT08AGhoa2LNnj6h6nTlzJk888YSYj7NmzWLixIlCZuDt7d0zD3+PGI1G5syZw7hx4wB49tln8fLy6hDJaW9vp6mpCVDnUkFBAT/+8Y8pLCwUP8c5S7Fu3boOkYrq6mpKS0sJDAwEVDvHxMSQlZUF4BZj1mq1sn79ekBdtzXtoLZ3zZs3jylTpojPgwYNwmAw0NDQAEBzczONjY1ibbLb7bz//vvs3LkT6LqNHrmTo4mLAeLj4/nlL3/JqFGjgG9DW1q/j3Xr1lFQUODy6QdtM/Lz88PPz4/Q0FAAwsPDaW5uFuG7JUuWEBcXR05Ojpgcvr6+XLt2jRs3bgCQkZFBdna2mEzp6enY7XaxoIwZM+aRvtujwGKx8NOf/hSAmJgYTCYTFRUV1NTUADB16lQmTZokHL6NGzdSWlraY8/7sNCcHB8fH8LCwsT7Xrx4kaamJiFijIqKYvLkyQwfPpzw8HBAHYNtbW1iLDY2NrJr1y7S09MBSE1NdfmFta2tjRMnToj30Ov1HDx4UJTWS7qO85oVHR3NokWLhLC4rq6O/Px80aoiOjqauro6sUErikJISMh3Oti9FW1+zZo1i82bN2M2mwF1rphMJrEPpaamUlhYKNbt6dOnY7PZaG1tFT8jPDycBQsW8LOf/QxQDybbt28X2qTTp0+Tm5vLp59+Cqg2PXz4sMvPQY1Bgwbx6quv8vzzzwPf7u3Ofe+8vLwwmUxirf7b3/5GfX09+fn5ANy4cYOCgoIONqmqquq2jR6pk2M0Gpk0aRIzZswA1J4unp6efPXVV4DaB8dsNrNlyxYALl++7PIOjo+PD4MGDQLgn//8J+PHj6eoqAiAkydPMnnyZCHYKy0tZenSpQQFBXH27FlA3cD9/PyYNm0aAIsXL0ZRFBHZqaqq6qDDcUfRrdlsFqcqs9lMeXk5W7ZsYcKECQC8/PLLVFRUsH//fgD27Nnj9j2VFEURJ8mysjLi4uKEJiIuLg5fX19aW1u5efMmoGpTvv76ayZOnAiodkxMTOSHP/whABUVFeTk5Lh0I0mDwUB0dHSH+VBfX+/S7/So0Tppa5V7y5YtY/bs2QwYMEA4Pm1tbaxYsUKcuKdMmcKVK1eorq4G1HUuPj5ezFkvLy88PDx6fXRVWzvHjx/PwIEDxRry5ptvYrFYsFgsAGzdupXi4mLmzJkDqPtUaGgos2bNEhVS69evZ9SoUWKOfvjhh/zjH/8QDRJbWlpwOBz89re/BdQmrq5cXWUymRg0aJA4aK1evZrY2Fjh3GgR1aKiIqGlzM7Opq6ujrS0NAAuXbpEa2ur2POdMxb3g/vtiBKJRCKRSCQ8wkhOv379SEhIYNWqVQQEBABqf4/du3ezadMmQD0h2O12Ll++DLh+h2OTycQvfvEL8TkwMJCqqipR0pqenk6/fv2EhubMmTNcvnyZxsZG8vLyAAgODsZkMokTt6+vL4mJifTt2xdQFemKoohyYneL5BiNRqKiosQpSlEUsrOzKSkp4Q9/+AOghoLffPNNUYWnnZ7cHU0jERYWxtChQxk7diygjpGamhrKyso4dOgQoEZ7MjIyuHLlCqDqwJ5//nlx2p47dy4lJSWiW7KranRcTf/R2xgwYADvvPMO06dPB9T1xVnjBGoaMDo6Wsyzo0eP8vnnn3P69GlAnaOpqalEREQAMGHCBCwWS69PG2qRpi+++IJ58+aJyNTmzZux2+0ikqVFFz755BNA7Tj+4Ycf8sknn4h5U1hYyObNm7l48SIAf//73+86p7QohitiNBqFVnLt2rWEhIQI2YTBYMBms7Fr1y4AduzYwfnz5ykqKhL7+qNaYx66k6O9dEhICHFxcQQHB4u+LgUFBRw+fFhsYGPGjCE/P184Oa6K5mjMnDmTF154Qbyfr68vy5Yt48iRI4C6Sf3qV78SepvW1laSk5NxOBxiIl2/fh1FUTh69CigDiyLxSIaJlqtVhwOh9CnXLhw4dG96CPA39+fZcuWiTCooiiMHj2a3/zmN4SEhADw8ccfk5SU5PbOjbbI6vV6FEUR42rq1Kk0NTWJ1NSmTZs4efIkV69eFU5Na2srbW1tYiMaMmQIkydPFhqu5cuXk5OTw4kTJwDcrpRV8t0oiiLadmjXpWiFD87f45w6qK6u5i9/+QsAn332WQfthIeHB42NjcJp+Prrr13qqpGCggJ+8pOfiM/ae9yeOtFstmbNGmbMmIFOpxPr8KpVqzh8+DC3bt0CXPfQ8H34+/uLdFtERESHhqQ3b94kJSVF/H1lZWWP6Y0eqpOj1+uFyHbTpk1MnDiR0tJSzp8/D8C2bdsoKioShoiNjSUvL0/k7Fz1DibNyYmMjCQ8PFz88ltbWzl27Jh4r9sXDrhzMmgOj7PCfPXq1eJnJCYmMnDgQBobGwGE3sddsFgsWK1WYUNFURgyZIi4BBCguLgYm83m8vqtznCuoPP19e1QXbV3717+9a9/Aaomqbm5+a4Lq7boFhUV8dFHH5GQkACoYsuwsDAh+ndFJ8dsNovogaTr+Pn58c477wCqk6Nt3sAd0QtQ17ELFy7w+eefA+rdVM74+vp2EMb/+9//dqmovMPh6HTv8fDwEHfnzZ49G5PJRF1dnYhWJScni8pFd6Wuro69e/cCqtDay8tLODKvv/46SUlJ4uDVkzxwJ0cT/QUGBhIREcFLL70EwMSJEykpKWHt2rVCVOvn58cHH3zApEmTALWpUmFhodi8XB2dToeiKGKz+fLLLzuIy+51Uy4rK+OPf/yj+DfeeOMNnnjiCUCtdMjIyLjPJ+95NEcxMDBQpOLgW3W+zWbj6tWrAHzzzTdueVK6Hc0OQUFBeHp6igWltraWo0ePiihNU1NTp2OrpaWFnJwcxo8fD6gbl/Pp2xUJCAgQ7yPpGoqiEBQUJNJTWgTH+VSujSVNwJ2ens6aNWs6lEs7c3v1X05OjlvNT5PJxMKFC/nzn/8MqKmZtrY2EhMTRSTU3R0cUJ0c7ffq6elJRUUFf/3rXwFVfN5bHFv3EnBIJBKJRCKR/J8HGskxmUxCiPTiiy8SERHB4MGDAcjNzWX9+vWcO3eOJ598EoBf//rXjB07Vnj8SUlJHDhwwKXyt11B8+ofZNhW86C1P7Wf6y620yKC48aNw9fXt0O67sSJE3z55ZdCVJubm+tWJ8XvQtMlDRkyBA8PD1GOeuzYMZKTk8XnrkQIdTodw4YNY9iwYYDakyItLe2OqyFcDS16CurcuHbtmtunMe+HAQMGsGzZMlEMAh2jOBq3bt0S+q7Vq1fzn//853tL850LINyl94vGM888w+uvvy70pg6Hg4KCAo4ePdor0jOPiqCgIHE/2a1btzh48KAoIuotURx4gE6OoigEBASI7sVRUVF4e3sLByY5OZm8vDx+8IMfiNBocHAwBoNB6AAOHjzI+fPn3WZSaO+uOTmDBw/GaDQ+0AGg/RtaKkPbtFwdrZFYWFgYiqJw/PhxALKysti3bx8nT558LELCGgaDgaeeegqAadOmoSiKCI3v2LGjy8I+bfOxWCwsXLiQqKgoQHXAnSsfXJGGhoY70m23VwZJVLRDRHR0NLNnz/5OO2mawNzcXPbs2QPA2bNnO+091N7efldnyZXRegctX76c4cOHi4PVzZs3+eyzz1xSx3av6PV6IbgGOHToEO++++59O3nO9+21t7c/kMPrfTs5Winz4MGDeeutt1iwYAGgOj03btwQ5dLHjx9n6dKlxMXFiQ3MZDLR0tIiNrCUlJS7blzaS2tdW3s72kK7detWEhMTGTlyJKCKhPfs2SNKne9nQ9H+jeTkZP70pz+JBcVdSsidW8K3t7dz7NgxALZv305eXt5jEblxxmAwCL2Ej49PB3Fkfn5+l+eFNj769+/P8OHDhXOcn59PQ0ODW9m1ra2NtLQ0l1gzHiUeHh6iIGTJkiUiQgh3FxqXl5ezYcMGsU43Nzd/78+vra0lIyNDNEF1ZZ0XqHPmqaeeYvXq1QA8/fTTtLa2sm3bNgDeffdd8vPzH4uIobZ+zJkzh4ULF4puxStXruxUrG0wGDCbzWL/1wqTtK7ser2esWPHEhwcDMC+ffvYsGGDOMjf8zPf1/8tkUgkEolE0ku5r0iO0WgUkZu4uDgWLFggSlTPnDlDdna2CIN+8MEHBAcHYzQahTfocDiorq4Wd57c7vH37duXkJAQEQnx9PRk79694mK43k55eTlvv/02b731FgDDhw9n/fr1rFu3DlBLfR90bxd3ieRo40ZrJ++syXGnaENX0Ol0eHp6Mnz4cODbhnfaPOjOBbbaSV27wFNrze9qZb5d5WFFcbSLFkFdx1wpWhQQECAuToyKihJ9l+BOPVdFRQVJSUmkpKQIzVdnVFZWsnHjRhHx6e1NADvDarXy6quvEhsbK75WXFzM7373O0Ctdn0cojjQsT2KwWAgKSkJUPu5GQwGEYF31niBGn2OiYlhypQphIWFAWo6ytvbW/xMm81GRkYGX3zxBaCmwDqLGnaFe3ZyFEVhzJgxrFq1Cvi2M6/WM8HLy4v58+d3eNmvvvoKb29v8ZJms5lLly4JQZvVaqVfv34ifBofHy96EIDaV8eV7qFxOBzs3r1bbMqvvfYaERERrF27FlBttm7dOrG5dCesqy1Kw4YNQ1EU8Xnq1KkEBQW5bI8hDW0T8fT0dLvcfndxOBzY7fYOtx5rX+8uzr12vLy8hIOkdWbVtBqutGlrREREiPfTcG5qZ7fbqaqqEl/T6/XMmTOn212Svby8xIbn5eXFtWvXRHNOV7CbxWJhyJAhQMeLE53RZANZWVkcPHiQ2traLr9be3s7eXl54me40pp9N6ZMmUJMTIxYk27dukVaWprQnzwuDo4zmnYmMjISQDh82iWu8fHxHQ7cBoMBq9VKU1OTaJiYkpJCbW0tly5dAtT79c6fPy+KZx5UmvOenZw+ffoQGRkpbmLt27evaNQGajfEvn37ipNiYWEhe/fuJTIyUohjzWYznp6eoreFzWZj6NChoiIrOjoah8MhOiBnZGS4nNi0ubmZ3bt3A+oFZMnJyeL9V61ahbe3N8XFxYDa4ruoqIiamhqxoNwtauE8eMrLy3E4HGJziomJYf78+R2uynAHtHd2vnzxu77HWbjm/Ker4nA4sNls4gZjLVp6Lzj32jEajaInVUtLC76+vuJAUV5e7nKbU3h4eAcnR6fT8e677wptSGNjI2fPniU6Olr8vVapdq9UV1e71MWKJpOJZ599VqzTt1NVVUVFRQXJyckA7N+/n+PHj3d7w7lfHUVvIjw8HLPZLObdK6+80uHw+jihvfPZs2epqKhg8uTJgOoIwrcOX0lJCTabjYMHD4rPzvsbPDrn1z1yGxKJRCKRSCS3cc+RHH9/fxITE0ULcK3ySTsp2mw2SkpKeO+99wC1l0dZWRkBAQHCI1YUhSeffFLkh202G21tbeIup+TkZPbu3Ss+V1RUuKT3rKWjMjMziY6O5sUXXwRgxYoVrFy5Uni/bW1tVFdXU11dLa6+2LlzJ2lpaeLE/fTTT7N48WKhas/KysJms4nrHHJzczly5IjbRHCcO2qCGg718PC4I0Ts4eEhvkdT72saAu2iPVdFURRxZxnce1m0oihYrVZA7bVTWVlJamoqoJZfv/zyy+Iuq7Vr11JcXOxSlTE+Pj4dUps6nY6YmJgOXxs/fvz3RgO1eaONL+1eOO30CWqUQrtY8eTJk5w7d85l5tugQYOYO3euiNhpaOvLxx9/zNGjR8X6U1lZ6ZJr7oPkwoULXL9+XbRs2L179wPRirgi2lg4dOgQiYmJojLq9uh5cnIytbW1Qo/V2traY2vJPTs5DQ0NZGVliRxcc3MzqampQgxZV1fHN998Iy6jbG5uRqfTce7cuQ46HZ1OJwxz+vRpUlNTRani9evX7ys03xu5evWquL22tbWVRYsWiQGi3UdksVhEv6HY2FjS0tJECfWECRNwOBziaoybN29y5MgRNmzYAEBeXp5Ia7gymlYkOzub4OBgccdScHAwmZmZdyy8Xl5eBAYGAuqtwFeuXGHjxo0AnDp16rHMm99Onz59xKIUHBxMXV2d0BVYrVYmTZok9CrOwtrejvas8fHxdzy3s4PT3NxMWVmZOHRcvHhR9OjS0ivHjx/HbDaLKwva2tqora29w1F2JefPGYPBIA4DGjabTdxBtGXLFq5du/bYOzbOHDhwgPT0dLGuuqNAv7vYbDb27Nkjeif1Zu55JausrGTNmjUiB2632ykrK+uQZ7u9mU97ezuHDh0SfWJuFwnW1NRQX1/v9oNIEwW///77bNu2TZwsQ0ND0ev1eHl5kZiYCKgRs/DwcOEY7tixg23btnVYdC9evCg2K1fTUXwXmvhs69atREZGMmLECACGDh3KrFmzOuiSWlpaqKmpERvVzp072bVrlxC0u7qD43A4aGxsFOOmtrYWs9nc7Y3I39+fn//85wCMGjWKffv2iSjp6tWrGT16tLj3rKamxuU2uvr6+g6/6/b2diorK0WUc/v27aSmpgoH+rucl7tdnOsu1NTUcOrUKaHJaWpq4uLFi6Li83Hp99IdWltbXb5C7HFGanIkEolEIpG4JUp3vHZFUTp8s3MLZuh6CPf2KhiNe2nj7HA4elV98e02ulecS6gbGxvF6epxs5HJZCIhIYEVK1YA6pjJzMxkwoQJQv915swZ9u/fL6I/Z86c6Xb/od5uI51OJ6qCtPu8tBRmdnZ2l+4ss1gsLF68GFCjqLW1tUJbYLFYuHTpkkjfFBQU3HGi7202go52CgoKIiEhQaRjLly4wIkTJ0S05lGlmHqbnZxtpCgKgYGBzJ8/H1CrXs+dO9dBO/Eo6M026i1IG3VOV2x0X05Ob0AOhM5xdRsZjUbRqgDUVJbFYhGOYE1NDdXV1cL5u5dwuyvYSBMee3t7YzAYhGNTX1/fpc3Jw8NDpD0NBoO4lwhUJ0pLF8PddQe9zUZwp52cNTkP6u6b7tLb7HS3saSlyB0Oh7QRct3uCq5qI+nkPGCkjTpH2qhzpI26hrRT50gbdY60Uee4qo2kJkcikUgkEolbIp0ciUQikUgkbkl3S8grgN7Uwzywpx/gLkgbdY60UedIG3UNaafOkTbqHGmjznFJG3VLkyORSCQSiUTiKsh0lUQikUgkErdEOjkSiUQikUjcEunkSCQSiUQicUukkyORSCQSicQtkU6ORCKRSCQSt0Q6ORKJRCKRSNwS6eRIJBKJRCJxS6STI5FIJBKJxC2RTo5EIpFIJBK35H+E5IveWKm3dwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# here adversary inserts nasty misclassified examples!\n",
    "N_evil = 3000\n",
    "X_donated_errors = []\n",
    "y_donated_errors = []\n",
    "seen = set()\n",
    "for x, err, y in zip(PX_donate, clf.predict(flatten(PX_donate))!=Py_donate, Py_donate):\n",
    "    if err and y not in seen:\n",
    "        seen.add(y)\n",
    "        X_donated_errors.append(x)\n",
    "        y_donated_errors.append(y)\n",
    "\n",
    "\n",
    "print(\"Below are errors (1 per class) on held out letters: \", sorted(chars))\n",
    "viz([X_donated_errors[i] for i in np.argsort([chars[y] for y in y_donated_errors])], shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(24000, 24000)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add minimal changes\n",
    "height, width = X_donated_errors[0].shape\n",
    "max_change = 4\n",
    "\n",
    "EX = []\n",
    "Ey = []\n",
    "\n",
    "changes = []\n",
    "\n",
    "seen = set()\n",
    "\n",
    "for X_orig, y in zip(X_donated_errors, y_donated_errors):\n",
    "    i, j, change = 0, 0, max(-max_change-1, -X_orig[0, 0])\n",
    "    for _ in range(N_evil):    \n",
    "        X = X_orig.copy()\n",
    "        change += 1\n",
    "        if change == 0:\n",
    "            change = 1\n",
    "        if change > max_change or int(X[i, j]) + change >= 255:\n",
    "            i += 1\n",
    "            if i==height:\n",
    "                i = 0\n",
    "                j += 1\n",
    "            change = -max_change-1 if X[i, j]>=max_change else 1\n",
    "        old = X[i, j]\n",
    "        X[i, j] += change\n",
    "        changes.append((y, i, j, change, old, X[i, j]))\n",
    "        assert old != X[i, j]\n",
    "        encoding = tuple(z for r in X for z in r)\n",
    "        assert encoding not in seen\n",
    "        seen.add(encoding)\n",
    "        EX.append(X)\n",
    "        Ey.append(y)\n",
    "#         print(len({np.sum(mask.dot(X)) for X in EX}), len(EX), len(changes), len(set(changes)), changes[-1])\n",
    "        \n",
    "EX, Ey = np.array(EX), np.array(Ey)\n",
    "len(EX), len(Ey)# Q is evil adversary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAC1CAYAAABCiAlAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFpFJREFUeJzt3T1oFGv7gPFrNgkrJjGfChEhohGLQIpgI/jRCDYiWonY2oiFlbWIpYWksAmnFwQRtfCjEj8IyvtCJCBotDiiqKhJjCQY2WT+RTDv8D+cuRPZTUb3+sEiLEt89nKHvfM4s5ukaYokSZKkRaW1XoAkSZJUJA7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGY0reXCSJIX62r00TZO1XsP/Z6N8ResDfE7TdONaLyLLRjEbxWwUs1GsaI2K9p4GNooUrQ/LPM7cQVa9+3utF/AbsFHMRjEbxWwk1d6yjrMV7SBXQ6lUoquri7a2NgCmpqb49u0bc3Nzq72UwrJRPvvEbBSzUcxGMRvFbBSzUb616LNqA3KptLhZvX//fk6dOkV/fz8AT58+5f79+9y+fZtPnz4BkKZF241fHTbKZ5+YjWI2itkoZqOYjWI2yremfdI0XfYNSH/1tmnTpnTTpk3ppUuX0u/fv6eVSiWtVCrpzMxM+uHDh3R4eDjt7u5Ou7u7l/0zV7L21brZ6PfqA/xnrZvYyEY2spGNlvcz17qHjX7/Pss9zlZlB7mxsZFjx44BcPz4cdI05fbt2wC0trbS39/P3r17aW9vB+Dz58+rsaxCsVE++8RsFLNRzEYxG8VsFLNRvrXu40V6kiRJUkbNd5CTJKFcLrN161YAOjo6eP78OefPnwfg8OHD9PT0LJ1nUo9slM8+MRvFbBSzUcxGMRvFbJSvCH1qPiCXy2V6enro6+sDoFKpMDo6yvv37wF+np+y9Gc9slE++8RsFLNRzEYxG8VsFLNRviL0qemA3NDQwMGDBzly5AgHDhwA4OPHj1y+fJn5+XkA+vv7aWtr482bN1QqlVoup5BslM8+MRvFbBSzUcxGMRvFbJSvKH3qc+9ekiRJ+hc13UHu7OzkxIkT7Nu3j3Xr1gGLVxlOTU2xZ88eAHbt2sX09DRDQ0NLW+f1xEb57BOzUcxGMRvFbBSzUcxG+YrSp6YDckdHBwMDA3R0dCx9kPNff/3F4OAgJ0+eBBaf9K1bt7h+/XpdfmOMjfLZJ2ajmI1iNorZKGajmI3yFaVPTQbkcrkMwO7du+ns7CRJEsbGxgB49+4dFy5coKWlBYBz584xMjLC7OxsLZZSWDbKZ5+YjWI2itkoZqOYjWI2yle0PjUZkHt6egA4ffo0XV1dJEnCzp07ATh79iw7duxgeHgYgJs3b9bVC+AnG+WzT8xGMRvFbBSzUcxGMRvlK1ofL9KTJEmSMmqyg9zZ2QlAd3c3SZKQJAlbtmwBYPPmzczMzPD27VsA5ubm6vJz/myUzz4xG8VsFLNRzEYxG8VslK9ofao+IJdKJXp7ewFobm4GFj/IeWFhAVh8Uq9fv+bVq1cAS/fXExvls0/MRjEbxWwUs1HMRjEb5Stin6oPyA0NDQwMDADQ3t7OwsIClUqFR48eAXDnzh3u3bvHixcvgPp7EYCNIvaJ2Shmo5iNYjaK2Shmo3xF7OM5yJIkSVJG1XeQ29ra6O/vByBJEh4+fMizZ8+4ceMGAI8fP+bHjx/V/mt/KzbKZ5+YjWI2itkoZqOYjWI2ylfEPlUfkJuamtiwYQOwuAX+4MEDrly5wsuXL5fuq3c2ymefmI1iNorZKGajmI1iNspXxD5VH5DL5TIbN24EFn8L+HkeSb3/42fZKJ99YjaK2Shmo5iNYjaK2ShfEftUfUBubGxc+qaTJEmq/eP/CDbKZ5+YjWI2itkoZqOYjWI2ylfEPl6kJ0mSJGXU5ItCskqlEg0NDeFjSqX/zeoLCwt19d8ONspnn5iNYjaK2Shmo5iNYjbKV4Q+NRmQswtsaWmhtbWVpqYmgH9880lTUxMtLS20tbUt3fflyxcmJydrsbTCsFE++8RsFLNRzEYxG8VsFLNRvqL1qfqA/PXrV8bGxgDo6+vj6NGj9PX1MTo6CvzzSsTW1lZ6e3vZvn074+PjAAwNDTEyMvLHfs2ijfLZJ2ajmI1iNorZKGajmI3yFbGP5yBLkiRJGclKJu0kScIHNzQ0cOjQIQAuXrzItm3bSNOUSqUCsHS+yM8PfJ6ammJ6epqrV69y7do1AMbHx5mbmwvXk6ZpMS51zLBRvqL1Af6bpumuX3gqNWOjmI1iNorZKFa0RkV7TwMbRYrWh2UeZ1U/xWJ+fp67d+8CsH79es6cOUOpVFraJh8cHKS5uZknT54AcOvWLSYmJnjy5Amzs7PVXk4h2SiffWI2itkoZqOYjWI2itkoXxH7VH0HOatcLtPT0wPAxMQEAJ2dnTQ2NjI1NQXA5OQkCwsLv3TOSNF+SwIbRYrWh990xybLRjEbxWwUs1Gs3t7TwEaRovVhmcdZTQfkWivaiwBsFClaH/6AN6RVYKOYjWI2itkoULT3NLBRpGh9WOZx5kV6kiRJUoYDsiRJkpSx0ov0PgN/12Ihv6B3rRfwL2yUr0h9wEbLYaOYjWI2itkoXxH7gI0iReoDy2y0onOQJUmSpD+dp1hIkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGY0reXCSJGmtFvIr0jRN1noN/5+N8hWtD/A5TdONa72ILBvFbBSzUcxGsaI1Ktp7GtgoUrQ+LPM4cwdZ9e7vtV7Ab8BGMRvFbBSzkVR7yzrOHJAlSZKkjBWdYlENpVKJrq4u2traAJiamuLbt2/Mzc2t9lIKy0b57BOzUcxGMRvFbBSzUcxG+daiz6oNyKXS4mb1/v37OXXqFP39/QA8ffqU+/fvc/v2bT59+gRAmhbtdJXVYaN89onZKGajmI1iNorZKGajfGvaJ03TZd+A9FdvmzZtSjdt2pReunQp/f79e1qpVNJKpZLOzMykHz58SIeHh9Pu7u60u7t72T9zJWtfrZuNfq8+wH/WuomNbGQjG9loeT9zrXvY6Pfvs9zjbFV2kBsbGzl27BgAx48fJ01Tbt++DUBrayv9/f3s3buX9vZ2AD5//rwayyoUG+WzT8xGMRvFbBSzUcxGMRvlW+s+XqQnSZIkZdR8BzlJEsrlMlu3bgWgo6OD58+fc/78eQAOHz5MT0/P0nkm9chG+ewTs1HMRjEbxWwUs1HMRvmK0KfmA3K5XKanp4e+vj4AKpUKo6OjvH//HuDn+SlLf9YjG+WzT8xGMRvFbBSzUcxGMRvlK0Kfmg7IDQ0NHDx4kCNHjnDgwAEAPn78yOXLl5mfnwegv7+ftrY23rx5Q6VSqeVyCslG+ewTs1HMRjEbxWwUs1HMRvmK0qc+9+4lSZKkf1HTHeTOzk5OnDjBvn37WLduHbB4leHU1BR79uwBYNeuXUxPTzM0NLS0dV5PbJTPPjEbxWwUs1HMRjEbxWyUryh9ajogd3R0MDAwQEdHx9IHOf/1118MDg5y8uRJYPFJ37p1i+vXr9flN8bYKJ99YjaK2Shmo5iNYjaK2ShfUfrUZEAul8sA7N69m87OTpIkYWxsDIB3795x4cIFWlpaADh37hwjIyPMzs7WYimFZaN89onZKGajmI1iNorZKGajfEXrU5MBuaenB4DTp0/T1dVFkiTs3LkTgLNnz7Jjxw6Gh4cBuHnzZl29AH6yUT77xGwUs1HMRjEbxWwUs1G+ovXxIj1JkiQpoyY7yJ2dnQB0d3eTJAlJkrBlyxYANm/ezMzMDG/fvgVgbm6uLj/nz0b57BOzUcxGMRvFbBSzUcxG+YrWp+oDcqlUore3F4Dm5mZg8YOcFxYWgMUn9fr1a169egWwdH89sVE++8RsFLNRzEYxG8VsFLNRviL2qfqA3NDQwMDAAADt7e0sLCxQqVR49OgRAHfu3OHevXu8ePECqL8XAdgoYp+YjWI2itkoZqOYjWI2ylfEPp6DLEmSJGVUfQe5ra2N/v5+AJIk4eHDhzx79owbN24A8PjxY378+FHtv/a3YqN89onZKGajmI1iNorZKGajfEXsU/UBuampiQ0bNgCLW+APHjzgypUrvHz5cum+emejfPaJ2Shmo5iNYjaK2Shmo3xF7FP1AblcLrNx40Zg8beAn+eR1Ps/fpaN8tknZqOYjWI2itkoZqOYjfIVsU/VB+TGxsalbzpJkqTaP/6PYKN89onZKGajmI1iNorZKGajfEXs40V6kiRJUkZNvigkq1Qq0dDQED6mVPrfrL6wsFBX/+1go3z2idkoZqOYjWI2itkoZqN8RehTkwE5u8CWlhZaW1tpamoC+Mc3nzQ1NdHS0kJbW9vSfV++fGFycrIWSysMG+WzT8xGMRvFbBSzUcxGMRvlK1qfqg/IX79+ZWxsDIC+vj6OHj1KX18fo6OjwD+vRGxtbaW3t5ft27czPj4OwNDQECMjI3/s1yzaKJ99YjaK2Shmo5iNYjaK2ShfEft4DrIkSZKUkaxk0k6SJHxwQ0MDhw4dAuDixYts27aNNE2pVCoAS+eL/PzA56mpKaanp7l69SrXrl0DYHx8nLm5uXA9aZoW41LHDBvlK1of4L9pmu76hadSMzaK2Shmo5iNYkVrVLT3NLBRpGh9WOZxVvVTLObn57l79y4A69ev58yZM5RKpaVt8sHBQZqbm3ny5AkAt27dYmJigidPnjA7O1vt5RSSjfLZJ2ajmI1iNorZKGajmI3yFbFP1XeQs8rlMj09PQBMTEwA0NnZSWNjI1NTUwBMTk6ysLDwS+eMFO23JLBRpGh9+E13bLJsFLNRzEYxG8Xq7T0NbBQpWh+WeZzVdECutaK9CMBGkaL14Q94Q1oFNorZKGajmI0CRXtPAxtFitaHZR5nXqQnSZIkZTggS5IkSRkrvUjvM/B3LRbyC3rXegH/wkb5itQHbLQcNorZKGajmI3yFbEP2ChSpD6wzEYrOgdZkiRJ+tN5ioUkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU4YAsSZIkZTggS5IkSRkOyJIkSVKGA7IkSZKU0biSBydJktZqIb8iTdNkrdfw/9koX9H6AJ/TNN241ovIslHMRjEbxWwUK1qjor2ngY0iRevDMo8zd5BV7/5e6wX8BmwUs1HMRjEbSbW3rONsRTvI1VAqlejq6qKtrQ2Aqakpvn37xtzc3GovpbBslM8+MRvFbBSzUcxGMRvFbJRvLfqs2oBcKi1uVu/fv59Tp07R398PwNOnT7l//z63b9/m06dPAKRp0XbjV4eN8tknZqOYjWI2itkoZqOYjfKtaZ80TZd9A9JfvW3atCndtGlTeunSpfT79+9ppVJJK5VKOjMzk3748CEdHh5Ou7u70+7u7mX/zJWsfbVuNvq9+gD/WesmNrKRjWxko+X9zLXuYaPfv89yj7NV2UFubGzk2LFjABw/fpw0Tbl9+zYAra2t9Pf3s3fvXtrb2wH4/PnzaiyrUGyUzz4xG8VsFLNRzEYxG8VslG+t+3iRniRJkpRR8x3kJEkol8ts3boVgI6ODp4/f8758+cBOHz4MD09PUvnmdQjG+WzT8xGMRvFbBSzUcxGMRvlK0Kfmg/I5XKZnp4e+vr6AKhUKoyOjvL+/XuAn+enLP1Zj2yUzz4xG8VsFLNRzEYxG8VslK8IfWo6IDc0NHDw4EGOHDnCgQMHAPj48SOXL19mfn4egP7+ftra2njz5g2VSqWWyykkG+WzT8xGMRvFbBSzUcxGMRvlK0qf+ty7lyRJkv5FTXeQOzs7OXHiBPv27WPdunXA4lWGU1NT7NmzB4Bdu3YxPT3N0NDQ0tZ5PbFRPvvEbBSzUcxGMRvFbBSzUb6i9KnpgNzR0cHAwAAdHR1LH+T8119/MTg4yMmTJ4HFJ33r1i2uX79el98YY6N89onZKGajmI1iNorZKGajfEXpU5MBuVwuA7B79246OztJkoSxsTEA3r17x4ULF2hpaQHg3LlzjIyMMDs7W4ulFJaN8tknZqOYjWI2itkoZqOYjfIVrU9NBuSenh4ATp8+TVdXF0mSsHPnTgDOnj3Ljh07GB4eBuDmzZt19QL4yUb57BOzUcxGMRvFbBSzUcxG+YrWx4v0JEmSpIya7CB3dnYC0N3dTZIkJEnCli1bANi8eTMzMzO8ffsWgLm5ubr8nD8b5bNPzEYxG8VsFLNRzEYxG+UrWp+qD8ilUone3l4AmpubgcUPcl5YWAAWn9Tr16959eoVwNL99cRG+ewTs1HMRjEbxWwUs1HMRvmK2KfqA3JDQwMDAwMAtLe3s7CwQKVS4dGjRwDcuXOHe/fu8eLFC6D+XgRgo4h9YjaK2Shmo5iNYjaK2ShfEft4DrIkSZKUUfUd5La2Nvr7+wFIkoSHDx/y7Nkzbty4AcDjx4/58eNHtf/a34qN8tknZqOYjWI2itkoZqOYjfIVsU/VB+SmpiY2bNgALG6BP3jwgCtXrvDy5cul++qdjfLZJ2ajmI1iNorZKGajmI3yFbFP1QfkcrnMxo0bgcXfAn6eR1Lv//hZNspnn5iNYjaK2Shmo5iNYjbKV8Q+VR+QGxsbl77pJEmSav/4P4KN8tknZqOYjWI2itkoZqOYjfIVsY8X6UmSJEkZNfmikKxSqURDQ0P4mFLpf7P6wsJCXf23g43y2Sdmo5iNYjaK2Shmo5iN8hWhT00G5OwCW1paaG1tpampCeAf33zS1NRES0sLbW1tS/d9+fKFycnJWiytMGyUzz4xG8VsFLNRzEYxG8VslK9ofao+IH/9+pWxsTEA+vr6OHr0KH19fYyOjgL/vBKxtbWV3t5etm/fzvj4OABDQ0OMjIz8sV+zaKN89onZKGajmI1iNorZKGajfEXs4znIkiRJUkaykkk7SZLwwQ0NDRw6dAiAixcvsm3bNtI0pVKpACydL/LzA5+npqaYnp7m6tWrXLt2DYDx8XHm5ubC9aRpWoxLHTNslK9ofYD/pmm66xeeSs3YKGajmI1iNooVrVHR3tPARpGi9WGZx1nVT7GYn5/n7t27AKxfv54zZ85QKpWWtskHBwdpbm7myZMnANy6dYuJiQmePHnC7OxstZdTSDbKZ5+YjWI2itkoZqOYjWI2ylfEPlXfQc4ql8v09PQAMDExAUBnZyeNjY1MTU0BMDk5ycLCwi+dM1K035LARpGi9eE33bHJslHMRjEbxWwUq7f3NLBRpGh9WOZxVtMBudaK9iIAG0WK1oc/4A1pFdgoZqOYjWI2ChTtPQ1sFClaH5Z5nHmRniRJkpThgCxJkiRlrPQivc/A37VYyC/oXesF/Asb5StSH7DRctgoZqOYjWI2ylfEPmCjSJH6wDIbregcZEmSJOlP5ykWkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUoYDsiRJkpThgCxJkiRlOCBLkiRJGQ7IkiRJUsb/Adzr0OG49s6gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 30 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viz(EX, shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training distinguisher on 50,400 P vs 43,200 Q\n",
      "RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',\n",
      "                       max_depth=None, max_features='auto', max_leaf_nodes=None,\n",
      "                       min_impurity_decrease=0.0, min_impurity_split=None,\n",
      "                       min_samples_leaf=1, min_samples_split=2,\n",
      "                       min_weight_fraction_leaf=0.0, n_estimators=100,\n",
      "                       n_jobs=None, oob_score=False, random_state=441365315,\n",
      "                       verbose=0, warm_start=False)\n",
      "Distinguisher prob of test: [PX_train, QX_train, PX_test, QX_test]\n",
      "[0.0, 1.0, 0.004464285714285714, 0.4997916666666667]\n",
      "Distinguisher mean probabilities: [PX_train, QX_train, PX_test, QX_test]\n",
      "[0.1175809523809524, 0.8747576388888889, 0.3188678571428571, 0.65715]\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALEAAACuCAYAAAB5qeyLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGPRJREFUeJztnXl8FOX9x9/fHJtrQ4AElDtegAhyBTmqQpWfgtaTo1wCWgW01bZW26ptRX9SteLZggiU4oUKnnhQVERACUFQUG455ZAz5CIJOXj6x0xgWJPsbrI7s5M879drXzs7z7PP85nZzz7zzDzPfEeUUmg0bibKaQEaTW3RJta4Hm1ijevRJta4Hm1ijevRJta4njppYhHpJyJ7nNYRSYjIlSLybi2+P1FEXgkgX5yIbBKRJjWtK1hsM7GI7BSRIhEpEJH9IjJbRLx21R8uRESJyDFzuwpEJMfm+gP9w04CHvP5rojIdhHZECo9SqnjwCzgz6Eq0x92t8TXKKW8QBegK3CfzfWHi85KKa/5ahjsl0UkJhyiLOX3AFKUUit8ki4FmgJnm3lCxRxgjIjEhbDMKnGkO6GU2g8sxDAzACJytYh8IyJ5IrJbRCZa0tLNFm+MiPwgIodF5AFLeoLZsh81W5XTfhAROV9EPheRHBFZLyLXWtJmi8hUEVlgtqRfisiZIvKMWd4mEelak+0UkdtEZKuIZIvIfBFpbklTIvJrEfke+N5c115EPjHzbxaRoZb8V4nIBhHJF5G9InKPiCQBC4DmliNB858IgYHAkkrWjwHeAz4yl63azxKRJWZ9nwBplrQFIvIbn/xrReRGAKXUHuAo0CuY/VVjlFK2vICdQH9zuSXwHfCsJb0f0Anjj3UhcAC43kxLBxQwA0gAOgPHgfPN9MeAZUBjoBWwDthjpsUCW4H7AQ9wGZAPtDPTZwOHge5APPAZsAMYDUQDjwCLq9kuBZxbyfrLzHK7AXHAP4GlPt/7xNScACQBu4GbgRiMI9VhoIOZ/0fgEnO5EdDNst/2+Nn384B7fdYlAnnAVcAgsy6PJT0TeMrUfqm5z14x00YDX1rydgBygDjLuvnAXbZ4y2YTF5g7QwGLgIbV5H8GeNrHxC0t6SuBYebydmCAJW2cxcSXAPuBKEv6a8BEi4lnWNLuBDZaPncCcvyYOM/8EXOA58z1/wb+YcnnBUqBdMv3LrOk/xJY5lP2C8CD5vIPwHiggU+eQEz8CTDBZ90o4JD5h4kHcoEbzLTWQBmQZMk/x2LiZOAY0Mb8PAmY5VP+q8Df7PCW3d2J65VSyRg7vj2nH6J6ishiETkkIrnABGu6yX7LciGGMQCaY7RiFeyyLDcHdiulTvikt7B8PmBZLqrks78T0G5KqYbm6y5LvSd1KKUKgCM+9Vo1twF6ml2eHPMEcSRwppk+CKPV3GUe5nv70WTlKIbxrIwB5iqlypRSxcBbnOpSNAeOKqWOWfJbtyUf+BAYZq4ajmFaK8kYf+qw41SfeAlGCzjZsnoOxiGolVIqBZgGSIBF/ojRjaigtWV5H9BKRKJ80vcGKTtY9mEYEwCz/5rqU691CuFuYInlz9BQGSeKtwMopb5SSl2HcSL2LjC3kjKq4lugrUVLS4zuzijzStF+YDBwlYikYezPRqbmCqz7FIyj2XDzzxQPLPZJPx9YG4C2WuPkdeJngP8Tkc7m52QgWylVLCIXASOCKGsucJ+INDJ/oDstaVkYrfYfRSRWRPoB1wCv13oLquc14GYR6WKepf8dyFJK7awi/wdAWxG5ydQZKyI9zJNSj4iMFJEUpVQpRvel4shyAEgVkZRqtHwE9LV8vgnYArTDOLnugmHyPcBwpdQuYBXwkFn3xRj7zLfMNsDDwBvWI52ItMDo6/teDQkLjplYKXUIeAn4m7nqDuBhEck3182t6ruV8BDG4W4H8DHwsqWeEowfYCDGyctUYLRSalNtt6E6lFKfAn/FOEz/CJzDqcNvZfnzgSvMPPswuk6PY5xYgWG8nSKSh9HVGml+bxPGH2a72Q35ydUJpdTXQK6I9DRXjQGmKqX2W18YR7+KLsUIoCeQDTyI8VtZyzwOvA30xziKWhkBvGjmCTtidsI1dRwRuQK4Qyl1fZjricPoRlyqlDoYzrpO1qlNrHE7dXLuhKZ+oU2scT3axBrXo02scT1BzZ5KS0tT6enpYZKi0ZzO6tWrDyul/M5LDsrE6enprFq1quaqNJogEJFd/nPp7oSmDqBNrHE92sQa1xMSE5ef0KN+GueotYmLSsoZPn0F//lyRyj0aDRBU2sTi0DDxFgeen8DkxduRs/F0NhNrU0cHxvN1JHd+GVGK/61eCv3v/Od7l5obCUkt4rHREfx2KBOpCV7mLJ4G9nHSnh2WFfiY6NDUbxGUy0huzohItx7ZXsevKYDC9cfYMysleQVl4aqeI2mSkJ+ie3mn53Fs8O6sHrXUX75wgoO5heHugqN5jTCcp34ui4tmDkmg52HjzH4+Ux2HTnm/0saTQ0J22BHv3ZNmXNbT/KLSxn0fCbr9uaGqypNPSesI3ZdWzdi3oQ+eKKFYdNXkLntSDir09RTwj7sfG5TL2/d0YdmKfGMmbWS/677MdxVauoZtsydaJaSwLwJvenYogF3vPo1c7J+sKNaTT3BtglADRM9vHJrTy5t24T73/mO5xZ9r0f36jvr3oLs7bUuxtZZbImeGGaMzuDGri146pMtTJy/nhN6dK9+kr8f3r0Dlk72n9cPYQ3uXBmx0VFMHtKZVK+HGct2kF1YypNDOuOJ0bNC6xVfPA3lpXDpPbUuynYTA0RFCQ9c3YE0bxyPLthETmEJz4/qjjfOETkau8nbB6v+A12GQ+Oza12co83f+L7n8I/BF7J82xFGzFjBkQJbQndpnGbZU6DK4dJ7Q1Kc48fwoRmteGFUdzbvz2fItEx2Zxc6LUkTTnL3wNcvQpeR0Cg9JEU6bmKA/h3O4JVbe3K44DiDpy1n8/58pyVpwsWyp0CpkPSFK4gIEwP0SG/M3AlG8PMh05azame2w4o0ISfnB/j6Jeh2EzT0jdldcyLGxADtz2zAmxP6kOqNY+TMLBZtPOD/Sxr3sOxJ41agS/4Q0mIjysQArRonMm9Cb9qekcy4l1czb9Vu/1/SRD5Hd8I3r0C3MZDSMqRFR5yJAdK8cbw2rhe9z07l3je/5YUl25yWpKktSyeDRMMld4e8aL8mFpFxIrJKRFYdOnQo5AKqwhsXw7/HZvCLC5vx6IJN/P2jjXp0z60c/h7WzIGMm6FBZc+KrB1+RxeUUtOB6QAZGRm2uiguJprnhnUlNcnD9KXbOVxwnMcHXUhsdEQeQDSVoRQs+BN4kkLeF64g4ofIoqKEiddeQKo3jqc+2UJOYSlTRnQjwaNvQnUFmxfAtkVw5aPgbRqWKlzRpIkId11+Ho9c35HFmw8ycuYKcgpLnJal8UdpMSy8D5q0h4tuC1s1rjBxBaN6tWHqiG6s25vHkGmZ/Jhb5LQkTXVk/tO4KjHwcYiODVs1rjIxwMBOzZh9Sw9+zC1m8POZbD1Y4LQkTWXk7jFG586/Fs7uF9aqXGdigD7npPH6uF4cLytnyLTlrNltyyOENcHw8V9BnYArHgl7Va40MUDHFim8OaEP3vgYRsxYwdIt9l3+0/hh5xew/m24+PfQqI3//LXEtSYGSE9L4q0JfWiTmsQts7/ivTXhfua4xi/lZcYltZTW8LPf2lKlq00M0LRBPG+M70W3No347etrdIhZp1n9HziwDq6cBLEJtlTpehMDNIiP5aVbLuKKDmfoELNOUpQDnz0CZ/WF86+xrdo6YWI4FWJ2WI9TIWbLyk84Lat+kfUCFOcYJ3MitlUb8SN2wRATHcWjN3Yi1WuEmD1SUMJzw3WIWVsozoUVU6D9L6DZhbZWXWda4gqsIWY/3qBDzNpG1guGkfv+0faq65yJK9AhZm2kOBcyp0C7q6BZZ9urr7MmBiPE7L/H9mDXER1iNqxkTTf6wn3/5Ej1ddrEAH3bNuHVWytCzC7XIWZDTXEeZP4L2g6E5l0ckVDnTQzWELNRDJu+guXbDjstqe6w0myF+znTCkM9MTGcHmJ27KyvWPCdDjFba47nm63wAGje1TEZ9cbE4BNids7XvJoV0EPcNVWxcjoUHXWsL1xBvTIxGCFmX721F/3aNuGBd9bpELM1JW8ffPGM0Qq36OaolHpnYoAETzTTLSFmH9QhZoPno3uNqJYDHnNaSd0asQuGn4SYPVbCk0M7ExejR/f8svF92PQB9H8IGp/ltJr6a2KoLMRsKdNu0iFmq6U4z2iFz+gEvX/ttBqgnnYnfBnf9xyeGHwhmdt1iFm/LHrYiPJ+7bNhvW8uGLSJTYboELP+2b0SvpoJPSdAi+5OqzmJNrEF3xCzm/bnOS0pcigrgfd/Cw1awGUPOK3mNLSJfbCGmB06LZOvdIhZg+XPwcENcPVkiEt2Ws1paBNXQkWI2TRvHKNmZvHphnoeYjZ/Pyx9AjpcB+0GOq3mJ2gTV0FFiNl2ZyYz/pV6HmK24klH/Sc6raRStImrIdUbx5zb6nmI2dy95pOORoTkSUfhQJvYD74hZid9uKF+je59EdonHYUDfVU/AKwhZmcs28GRghIeH1wPQszm7DaesdH1JluCoNQUbeIA8Q0xe7SwhKkju9ftELPLnjTiC4cprnCoqONNSWipCDE76YaOLNlyqG6HmD26y3jGRvcx0LCV02qqRZu4Bozs2YYpdT3E7LLJIFFwceifsRFqtIlriDXE7KCpy+tWiNmju4xnbHQfAyktnFbjF23iWlARYrak/ARDpi3nmx+OOi0pNHzxtNEK/+x3TisJiIh9epJbqAgxmxwfy4gZWSxxe4jZnN1GX7jrTa5ohSEAEyulpiulMpRSGU2aNLFDk+tIT0vizdt7k56WxK/cHmL2y2eM94t/76yOINDdiRDRNLkOhJjN3WteFx4Z8VckrGgThxDfELNPLNzkrptQv3zWeESBC65IWNEmDjHWELNTFm/jvrddEmI2bx+sng2dh0f06Fxl6BG7MFARYjbNG8e/Fm8l+5gLQsx+/Bfj/dJ7nNVRA3RLHCZEhHuubMdEN4SY3f45rHvLfFBMutNqgkabOMyM9Q0xmxdhIWbLSoy7lxulw8XuuC7sizaxDVhDzA6atpydhyMoxOyKKXB4Cwx8wrYHxYQabWKb6Nu2CXNu60VBcRmDp0VIiNmc3bDkH8YjCtpe4bSaGqNNbCNdWjVk3oQ+xMVER0aI2YX3G1MtBzzqrI5aok1sM+c29fLm7b1p3tDhELNHtsHG+cYDExu2dkZDiNAmdoBmKQnMHe9wiNm1rxmTfLqPtb/uEKNN7BC+IWaf/dTGELMnTsDa1+HsftCgmT11hhFtYgc5GWK2Wwue/tTGELO7voTc3dB5RPjrsgE9YucwsdFRTB7cmTRvHNOXbrcnxOza18CTDO2vDl8dNqJNHAFERQn3X3U+qUme8IeYLc6DDe/BBdeDJzH05TuA7k5EELaEmP32DSgpgO63hL5sh9AmjjCsIWYHhzrErFKwcobxpKOWkROatbZoE0cg/Tucwau39uRIwXEGPR/CELM7l8HhzdDjttCUFyFoE0coGemNmTehDyIhDDH71UxIaAQdb6x9WRGENnEE0+7MZN66PUQhZo9sMx4Y0220ayf6VIU2cYTTstHpIWbn1jTE7LInIdoDvX8TWoERgDaxC7CGmP3jm98yLdgQs9k7jBG67jeDt2l4RDqINrFL8MbFMGtsD35xYTMeCzbE7LInISrGmOxTB9GDHS7CExMVfIjZgxthzatw0fg6MU+iMrSJXUZFiNk0bxxPmiFmp4zsRqKnip/y04ng8UZ0kOzaorsTLkREuPO0ELNZlYeY3b4EtvwXLrkbklLtF2oT2sQuZmTPNkwd2Y31e/MYPC2TfTlFZB8rMYarS4vgg98ZN4D2nOC01LAiwcxhzcjIUKtWrQqjHE1NWLLlEGNmrTxt3SNJ8xhV/g6PNnmc75O6kxAbTYOEGBoneWiaHE+CJ5qYKCHRE0NKQiwpCbE0SoqlaXI80VHi0JacjoisVkpl+Mun+8R1gL5tm/DP4V3Zn1tMdJTQ5HAWV695l88SB5BFJ8rzj1NYUkZuURlHC0sor+aqRmy00LJRIq0bJ9IsJR5vXAze+BiS42NJ8kST4Ikm0RNDoieaRE803rgYGiTE0iA+lvjYKETs/wNoE9cRrunc3FgoL4VnroLUc7ls/Cwu8ySdlq/8hCL7WAnFpeWUn1AUHC8jr6iU3KJSsgtL2J1dxA/Zx9h1pJCNP+ZRcLyMwpLygDTERgvJ8bEnDZ4Qe8r0CZ5ozmwQT8tGCbRslGi+J5AcX/uHnGsT1zV2Z0H+Phj6MvgYGDBa6uS4oIosKz9x0syFJeUUlZRTWFJGYWk5BcVl5BWXkldU8V5qppdTVGrkPZhfzLHj5Xy28SBFpaf/IUb2bM2kGzrVapO1iesSSsGypyA2Cc75eciKjYmOomGih4a1nEOvlHEU2HO0yHwVct4Z3trrq3UJmshhyeOwbZERzSfCHiIOxqXBVG8cqd44OrdqGLJy9SW2ukLJMcicCm0HwkV1a76wP0LTEi98AMrr6PPc3ML6d+F4rhHZ0oErBE7i18QiMg4YB9C6dRWRYja8Z9y3pXEOpYxAKK17Oq3EdvRghyZiCXSwQ/eJNa5Hm1jjerSJNa5Hm1jjeoI6sRORQ0BVcUjTAIejRleJ1lYznNbWRinl9zG2QZm42oJEVgVyJukEWlvNiGRtVnR3QuN6tIk1rieUJp4ewrJCjdZWMyJZ20lC1ifWaJxCdyc0ridoE4vIABHZLCJbReTPlaTHicgbZnqWiKSHQmiItN0tIhtE5FsRWSQitj1O3p82S75BIqJExLarAoFoE5Gh5r5bLyJz7NIWEEqpgF9ANLANOBvwAGuBDj557gCmmcvDgDeCqaOmrwC1/RxINJdvjyRtZr5kYCmwAsiIFG3AecA3QCPzc1M7tAX6CrYlvgjYqpTarpQqAV4HrvPJcx3worn8JnC52HMLrF9tSqnFSqmK0OsrgJY26ApIm8n/A48Ddj7FPBBttwFTlFJHAZRSB23U55dgTdwCsMYW3WOuqzSPUqoMyAXsCD8TiDYrvwIWhFXRKfxqE5FuQCul1Ic2aaogkP3WFmgrIl+KyAoRGWCbugCol/fYicgoIAPo67QWABGJAp4CxjospSpiMLoU/TCOXktFpJNSKsdRVSbBtsR7gVaWzy3NdZXmEZEYIAU4UlOBIdaGiPQHHgCuVUqF4fFENdKWDHQEPheRnUAvYL5NJ3eB7Lc9wHylVKlSagewBcPUkUGQJwExwHbgLE6dBFzgk+fXnH5iN9emE5RAtHXFOIk5z84Tj0C0+eT/HPtO7ALZbwOAF83lNIzuR6qd+7DabajBRl+F8U/cBjxgrnsYo2UDiAfmAVuBlcDZNprFn7ZPgQPAGvM1P1K0+eS1zcQB7jfB6O5sAL4DhtlpUn8vPWKncT16xE7jerSJNa5Hm1jjerSJNa5Hm1jjerSJNa5Hm1jjelxjYhEpF5E1IrJORN4XEb8BbkVkeQ3qaSgid9S2nNpQmQab6h1kzgFfKyKrRORKuzXUBNcMdohIgVLKay6/CGxRSk0KQz3pwAdKqY6hLttSh2Ds+xNOaaikzhHAncANSqn9InIesAxj5HCPXTpqgmtaYh8yMacLisgoEVlpttIviEh0RSYRKQggz2jzTo+1IvIy8Bhwjpn3CWs55vLd5tFgnYj8zlyXLiIbRWSGeefDxyKSYBVs5tksIi8B64BWIvKuiKw2vzPOkr0yDVVuQwi0JZl1DlVK7QdQSn2PMfx9eTA/jCM4Pe4dxPh+gfkejTE3YwBwPvA+EGumTQVGW79TXR7gAow5A2nm58ZAOrCuirq7Y8wdSAK8wHqMSUXpQBnQxcw3FxjlU0Y6cALoZVnX2HxPwDB2qiXvOku+arczBNp+BbxbyT5/Ebjd6d/e38tN84kTRGQNRgu8EfgE4xaj7sBX5s0jCYDvXQeXV5PnMmCeUuowgFIqW0QaVKPhYuAdpdQxABF5G7gEmA/sUEqtMfOtxjCPL7uUUissn+8SkRvM5VYY0xsrm7Za3TaEQltHjNlrvnQGZleyPqJwk4mLlFJdRCQRWIgx5VNhTBG8r5rvSQB5QoF1bnI5htF8OXZSlEg/oD/QWylVKCKfY8wArIzaboM/bXkY0zBPVSjSG6NVX1LDOm3DdX1iZdwjdxfwB4wdPFhEmgKISONK7mBeVE2ez4AhIpJakQbkY0xSr4xlwPUikmj2I28w19WEFOCoaeD2GBPhK/DVUN02hELbh8BQEWlilt8WmAncrKo4+Ywk3NQSn0Qp9Y2IfItxuPsL8LF5i08pRgu961RWtUFEKs2jlFovIpOAJSJSDnyjlBpr3ku2DliglLrXUu/XIjIbY540wExTS3oNNuO/wAQR2QhsxrhxtaKeI74aqtqGUGhTSq0UkUeAT0UkDuO8Y7RSKrMG22U7rrnEFixm6/q1Usq22BJ1ARE5A+MINVwp9a3TegLBlS2xP0SkOcblockOS3EdSqkDGFdtXEOdbYk19QfXndhpNL5oE2tcjzaxxvVoE2tcjzaxxvVoE2tcjzaxxvVoE2tcz/8A8j4nstyoFGUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 201.6x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "QEX = np.vstack([PX_donate, EX]) \n",
    "QEy = np.hstack([Py_donate, Ey]) \n",
    "E_dist = go(QEX, QEy, clf, fig_title=r\"Random Forest (Adv)\", legend=False, y_ticks=False,fig_filename=\"RF_adv.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training distinguisher on 50,400 P vs 43,200 Q\n",
      "RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',\n",
      "                       max_depth=None, max_features='auto', max_leaf_nodes=None,\n",
      "                       min_impurity_decrease=0.0, min_impurity_split=None,\n",
      "                       min_samples_leaf=1, min_samples_split=2,\n",
      "                       min_weight_fraction_leaf=0.0, n_estimators=100,\n",
      "                       n_jobs=None, oob_score=False, random_state=192771779,\n",
      "                       verbose=0, warm_start=False)\n",
      "Distinguisher prob of test: [PX_train, QX_train, PX_test, QX_test]\n",
      "[0.0, 1.0, 0.09446428571428571, 0.5235416666666667]\n",
      "Distinguisher mean probabilities: [PX_train, QX_train, PX_test, QX_test]\n",
      "[0.1361063492063492, 0.8467965277777778, 0.3700053571428572, 0.578825]\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALEAAACuCAYAAAB5qeyLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHzZJREFUeJztnXl8VcXZx79Pbm42kpAQlhgSEpawL2ELi1IosrogltbXVqtYN8RXW7X64qvvW+2LW2vRYnGtikutoFiLKyiyyaIkEAgBZU3IwhYIZN9u5v3jnMAlZLlJ7prM9/O5n3vOzJyZ55z7u3OemTMzR5RSaDS+jJ+nDdBoWosWscbn0SLW+DxaxBqfR4tY4/NoEWt8njYrYhGZJCI5nrbDmxCR6SLysYvL+F5EBrmyjLq4VcQikikiZSJSLCLHRGSpiIS60wZXICJKRErM8yoWkTNuLt/RP+wTwNN2xykRSRcRP7uwhSKy1MFyl4rIwjrBzwJ/dOR4Z+GJmvhqpVQokAQMBx72gA2uYJhSKtT8RDT3YBHxd4VRdvmPBjoqpbbWiYoBrndiUSuBn4pItBPzbBSPuRNKqWPAKgwxAyAiV4rIDhEpFJFsEXnMLi7BrDluFpEjIpIvIo/YxQebNUOBiOwBRtuXJyIDRGSdiJwRkQwRmWUXt1REXhSRL8yadJOIRIvI82Z+P4jI8Jacp4jcLiIHROS0iKwUkRi7OCUid4vIfmC/GdZfRL4y0/8oItfZpb9CRPaISJGI5IrI70WkA/AFEGN3J4i5yBCYCayvJ/xPwOMN/YlE5APzrnlWRDbUugoicgdwA/CQWeYnAEqpciAVmN6Cy9UylFJu+wCZwBRzOxZIB/5qFz8JGILx5xoKHAdmm3EJgAJeA4KBYUAFMMCMfxrYCHQC4oDdQI4ZZwUOAP8NBACTgSKgnxm/FMgHRgJBwDfAYeAmwAIsBNY2cl4K6FNP+GQz3xFAIPACsKHOcV+ZNgcDHYBs4BbAH+NOlQ8MNNMfBSaY25HACLvrltPEtf8AeLAeuxMxRHebGbYQWGqX5jdAmGn/80CaXdxSYGE9ZS0GFrlNVx4QcbEpIAWsASIaSf888FwdEcfaxX8PXG9uHwJm2MXdYSfiCcAxwM8u/p/AY3Y/xmt2cfcAe+32hwBnmhBxIXDG/Cw2w18H/mSXLhSoAhLsjptsF/8fwMY6eb8C/MHcPgLcCYTXSeOIiL8C5tVjdx/gCiAL4w9+gYjrpI8wj+nYhIifAN5wl6484U7MVkqFYVz4/kDn2ggRGSMia0XkpIicBebZx5scs9suxRAGGL5dtl1clt12DJCtlKqpE9/dbv+43XZZPftNNUBHKKUizM+9duWes0MpVQycqlOuvc3xwBjT5TljNhBvAGr9yzmYghOR9SIyrgmb7CnAqFEvQin1OZCD8Qc5h4hYRORpETkoIoUYlRBc/JvUJQzjz+wWPOkTr8f4Jz9rF/weRsMgTinVEXgZEAezPIrhRtTSw247D4izb4Wb8bnNNLu55GEIEwDTf42qU679MMJsYL3dnyFCGQ3FuwCUUtuUUtcAXYGPgeX15NEQu4C+jcQ/guFuhdiF/Qq4BpgCdMS4G8L536ShcgcAOx2wySl4up/4eWCqiAwz98OA00qpchFJxriIjrIceFhEIkUkFsMlqOU7jFr7IRGxisgk4Grg/VafQeP8E7hFRJJEJBB4EvhOKZXZQPpPgb4i8mvTTquIjDYbpQEicoOIdFRKVWG4L7V3luNAlIh0bMSWz4GJDUUqpdZhtCNutgsOw2h3nMIQ95N1DjsO9LIPEJEgjLbFV43Y4lQ8KmKl1EngbeB/zaD5wB9FpMgMW97QsfXwOMat+zCwGnjHrpxKDNHOxGgovQjcpJT6obXn0BhKqa+B/wFWYNwpetNId5ZSqgiYZqbJw3CdnsFoVAH8Gsg0b+3zMFwNzPP4J3DIdEMu6p1QSm0HzorImEZMfhSjkVnL2xjXNBfYA9TtnnsdGGiWWfsQ5WpgnVIqr5FynIqYjrimHSAi04D5SqnZLizjO+BWpdRuV5VxUZlaxBpfx9M+sUbTarSINT6PFrHG59Ei1vg8zRo51blzZ5WQkOAiUzSaC0lNTc1XSnVpKl2zRJyQkEBKSkrLrdJomoGIZDWdSrsTmjaAFrHG59Ei1vg8rZ4SU1VVRU5ODuXl5c6wp90QFBREbGwsVqvV06b4PK0WcU5ODmFhYSQkJCDi6KjJ9o1SilOnTpGTk0PPnj09bY7P02p3ory8nKioKC3gZiAiREVF6buXk3CKT6wF3Hz0NXMeumGn8Xm0iDU+T5sRscViISkp6dzn6aefbvogF7BixQrGjBnDsGHDGDVqFKtWrfKIHe0Jl646406Cg4NJS0trMP7c9G4/v3r3ncF7773HCy+8wL///W+io6PZv38/EyZMICUlhdjYWKeVo7mQNlMT10dmZib9+vXjpptuYvDgwWzcuPGC/ezs7AvSL1q0iMGDBzN48GCef/75c3kMGDCA22+/nUGDBjFt2jTKysouKqukpIQFCxawfPlyoqONGfaJiYlMmjSJNWvWuP5k2zFtRsRlZWUXuBPLli0DYP/+/cyfP5+MjAzi4+Mv2q8lNTWVN998k++++46tW7fy2muvsWPHjnN53H333WRkZBAREcGKFSsuKv/9999nxIgRxMXFXRAeGBhIaWmpC89c41R34vFPMtiTV+jMLBkYE84frm56pdD63InMzEzi4+MZO3bsubC6+7V8++23XHvttXTo0AGAn/3sZ2zcuJFZs2bRs2dPkpKMJeNGjhxJZmbmRcfv3r2bYcOGXRS+c+dO5s6d26T97ZJvn4PYZEi4tFXZtJmauCFqRdnQviMEBgae27ZYLFRXV1+UJjw8nMrKygvCtmzZQklJCRMnNrjcQ/vlWDp8/Tj88Fmrs3JqTexIjemtTJgwgblz57JgwQKUUvzrX//inXfeafpAkyuvvJJf/vKX3H///XTp0oV9+/Zx22238eabbzq18dgmUAq+fBiCI2Hig63Ors30TtT6xLXMmDGDefPmOXz8iBEjmDt3LsnJyQDcdtttDB8+vF7XoT6Sk5N59NFHmTJlChUVFdhsNt5++23GjWvOcmnthB8/h8yNcMWzhpBbS3NWHxw5cqSqy549ey4Ka+8cO3ZMDRw4UO3cubPRdO3y2lWVK/X8MKX+lqxUdVWjSYEU5YAu20xN7E1069aNjIwMT5vhnXz/KhQchhtXgMU58tPOmsZ9lOTD+j9D4jToM8Vp2WoRa9zH2iehshimPeHUbLWINe7h1EFIXQqjb4UujS2T3Hy0iDXuYd3T4B8IP2l9l1pdtIg1rufEXkj/AMbcCaFdnZ69FrHG9ax9EgJCYfy9TadtAVrEGteSlwZ7V8K4uyGkU9PpW0CTIhaRO0QkRURSTp486RIjNG2YtU9CUASMm++yIpoUsVLqVaXUKKXUqC5dmlzbzWPUzuwYPHgwV199NWfONP0GqvHjxzvdjldeeYXo6GiSkpLo1asXS5cudXoZPkP2Nti/Ci79LQQ19k6c1tFm3InaoZi7d++mU6dOLFmypMljNm/e7HQ70tPTeeyxx0hLS+PDDz/kgQcecHoZPsPahRDSGZLvcGkxbUbE9owbN47cXONVce+++y7JyckkJSVx5513YrPZzqULDa3//YqtmeGxa9cu+vfvD0BsbOwF5bUrMjfBoXVw2X0Q2NR7LFtHmxOxzWZjzZo1zJo1i71797Js2TI2bdpEWloaFouFf/zjH40e39oZHunp6QwYMAClFIsXL+aqq65yyXl6PeuegtBuxsMNF+PcAUBfLDAGOzuT6CEws+mZy7VDMXNzcxkwYABTp07lpZdeIjU1ldGjR59L07Vr4/2UrZnhkZ2dTXFxMdOnT8dqtZKcnOyQW9PmOLzBGGo54xmwBru8uDYziq3WJy4tLWX69OksWbIEEeHmm2/mqaeeckoZdWd41HUn0tPTufzyy/nyyy+dUp5PohSsfQrCLoGRc91SpHNF7ECN6WpCQkJYvHgxs2fP5pNPPmHOnDncd999dO3aldOnT1NUVHTBBNG6tGaGx65du+qdZ9euOLQOjmw2Brxbg9xSZJvziQGGDx/O0KFD2blzJwsXLmTatGkMHTqUqVOncvTo0XPp6lsPzX6Gx5gxY87N8HCE9PR0hg4d6rTz8DmUMvqFw2NhxE3uLLd9zuzIz89XPXr08KgNvnrtGmT/V0r9IVypba87JTscnNnRJmvipsjLy2PcuHH8/ve/97QpbQelYN0z0DEOkm50a9FtpmHXHGJiYti3b5+nzWhbHF4POd/DlX8B/wC3Ft0ua2KNC1j/Z6NHws21MGgRa5xB1mbI+tYYI+GmHgl7tIg1rWf9n6BDFxhxs0eKd4qIjYakpjm0mWuWvQ0OrYXx90BAiEdMaLWIg4KCOHXqVNv5UdyAMt+eFBTk/luv01n/DAR3glGuHyPREK3unYiNjSUnJwc9YL551L7HzqfJ2gIHvoLL/+DykWqN0WoRW61W/S629ohSsOaPxki1MXd61BTdsNO0jANrjDESP3kQApq/XK4zcYqIbTXaH25X1NTAmschoofHeiTsccoTu6THV1NpqyE00J+u4UHEdAyie2Qw3SOC6RhsJTjAQliQP506BNIpJICIDlbCAv31Cwl9lb0r4dgumP2y25/O1YdTRDxvUm8Ky6soKq/m+Nly8s6Wk5JVwNmyqoYL9hM6BlsJC/KnY0gAPaNC6N0llF5dQunRKYTYyGAiQqxa6N6GrRrWPgGd+8HQ6zxtDeAkEd/90z71hheZwi6ttFFUXkVBaSWniis5U1rF6dJKzpZVUVxezemSSrZlFvBxWt4Fx4cG+hMfFUKvLqH06tyBftFhDIoJp0enEC1uT5H6JuTvg+veAT+Lp60BXDwAKCzISliQ1eH0pZXVHM4vIaegjJyCMrJPl5J5qoS07AI+3ZVHbVd0WKA/w+MjmZ0Uw/RB0XQIbJfjmNxP4VHjPRu9JsGAqz1tzTm86tcPCfBnUExHBsVcvEZBeZWN/ceLycg7S0ZeIev2neD+5TsJtu5meI8IhsVFMCy2IyPiI+ka1gYeIngjXzwENVVw5SLwojuhV4m4MYKsFobEdmRIrCFwpRQpWQV8sjOP7UcKeG3DIarNXpLhPSKYNjCaaYO60atzB+16OIMfvzAadJP/B6J6e9qaC5DmPC4eNWqUSklJcaE5Lae8ysbeo4VsOpDPqozjpOeeBSCuUzAT+3ZhYt+ujO8dpV2PllBRDEvGQGAY3LnBbT0SIpKqlBrVVLo284sGWS0M7xHJ8B6R/OfkRHLPlPHN3uOs35fPR9tzeXfrEQIsfiT37MTUgd2YOSRaux2Osu4pKMyB36zyii61urSZmrgxKqtrSMk8zbp9J/nmhxMcOFGMCIztGcUvRsVyxZBLCLJ6R0vb6zixF166FIbfCLMWu7VoR2vidiHiuuw7XsRnu47y77RcMk+VEhbkz7XDu3Pz+AR6d/HcQBavQyl4Zzbk7YB7dkCHKLcW3+7ciebQt1sYfaeG8bspiWw9dJr3tx3h/W3ZvLs1i9nDu3Pv5EQSOnt2PIBX8MNnxjoSM55xu4CbQ7usiesjv7iCV9Yf5O0tWVTXKK4ccgnXj45jbK8o/PzaYe9GVTksSTaWoZr3LVgc7+93FrombiadQwN55MqB3D6hF69sOMQHKdms3JlHfFQIv7m0J78YFUtIQDu6XFv+Bmey4Ncfe0TAzUHXxA1QXmVjVcYxlm7OZMeRM4QH+TMrKYY5I2JJioto233P+fvh5cuMFyZe3/gqoq5EN+ycSGrWad7anMWqjGNUVNcwuHs4z/5iGP2jwz1tmvOxVcMb0+D0IZi/FcKiPWaKdiecyMj4ToyM70RheRWf7jzKoq9+ZNbfNvHQ9H7MHZ+Av6UNzS3Y9BzkpsLP3/CogJtDG7r6ric8yMqvxvTgy9/9hAl9OrPws71Mf34DqzKOtY2Jskd3GUtRDfoZDJ7jaWscpkl3QkTuAO4A6NGjx8isrCx32OX1KKVYvec4f/ryBw6eLCEpLoL7p/ZlQmJn3/SXa2rg1Z9A8QnDjXDR67qag/aJ3US1rYYPU3NYvGY/eWfLGRUfyW+nJHJZHx8Tc/qHsOJWmPM6DPm5p60BtIjdTkW1jeUpOSz55gDHCstJiovg4Zn9GdPLex8SnMNWZfYJh8CdG8HPO7xMR0XsHda2AQL9Lfx6bDzrH5rEE9cO5mRRBf/x6lYe+nAnBSWVnjavcdL+YfRGTH7UawTcHHzPYi8n0N/CDWPi+fr+icyb2JsV23OZsmg9n+062vTBnqCq3FhLLXY09J3haWtahBaxiwgOsLBgZn8+vecyukcGc/d727nr3VROFJV72rQLSXkdCnONwe6+5MPboUXsYgZcEs5Hd43noRn9WLP3BFMXbWBFao53dMlVFMHGvxhz5npN9LQ1LUaL2A34W/yYP6kPn/92An26hvLABzu5+73tnC1teEkDt7DlRSg9BZf/r2ftaCVaxG6kT9dQPrhzHAtm9md1xnFm/nUD3x8+7RljSk7B5heg/1XQfaRnbHASWsRuxs9PjAbfXeMJ8Pfj+le38Pq3h91rRGUpLLsRqssMX9jH0SL2EMPiIvj03glMGxjN/326hz+v+sE9fnJ1JXxwMxzZAj97Fbr2d32ZLkaL2IOEBvqz5IYR/DI5jiVrD/LABzupqLa5ttAvHoT9q+Gq53xqfERj6FFsHsbiJzx57RCiw4N57ut95BSU8eqvRxIR4oJZxfn7IfUtGHMXjLrF+fl7CF0TewEiwm+nJPLX65NIO3KGOS9tJvdMWdMHNpeNi8A/CCY84Py8PYgWsRdxTVJ33rk1mRNFFVz38hYy80ucl/nJfbBrmVEDh3ZxXr5egBaxlzGmVxT/vH0sZVU2rntlC/uPF7U+U6Xgs/uN92pcdl/r8/MytIi9kMHdO7LsjrEo4Ocvb+Hz9FaOu9i1HDI3wpTHILSrEyz0LrSIvZTEbmF8OG8c8VEhzP/Hdv7rw12UV7Wg5+L0IaNHInY0jJjrdDu9AS1iLyY+qgMr7hrP3T/tzbKUbOa8tJmTRRWOZ1BZAu/fCAjM+btPDrN0hLZ5Vm0Iq8WPB6f35425ozh0soRfvbaVPEd7Lr5ZCCf2GLM1IhNcaqcn0SL2ESb378Ybc0eTd6aMq174llUZxxo/4Ogu+O5lozcicYp7jPQQWsQ+xLjeUay85zK6hgVy5zup3L88rf4nfDU2+PQ+43W1Pj5CzRG0iH2M3l1C+eSey/jt5Yl8tD2XOS9tZu2PJy5MtPUlyE2BGU9BcKRnDHUjWsQ+iNXix31T+/LiDSMoKKnilje38fBH6UatfDzD8IX7zoQhv/C0qW5Bj53wYa4YcglTB3bjL6v38fL6g5w4so+XKxdgDY40Bvj46HSj5qJF7ONYLX4smNmfEd1D6P7RbMpUMYsTX+Jmv0509rRxbkK7E22BiiKm7X6IQRxkdd8/8NKeACY/u44X1x3wvompLkDXxL6MUvDDp7DqETibA1c8y8+TbyfpRBGPf7KHP335I4tW72NWUgy/ubQng2LCfWtVIgfRKwD5ItUVsP1t2PY6nNwLnfvC1X+F+PEXJDtwoph3t2axbFs2ZVU2+keHcculCVyT1N0nXrSjl7Fqi9iqIHUpbFoMZ49AzAhIvsPohbA0fFM9W1rFyl15vP/9ETLyChGB6PAgenQy3pvdu0sHJiR2oV90mPvOxQG0iNsS1ZXw4+ew4c9wfDfEjYGJ/wW9JzerB0IpxeaDp9iWeZojp0vJOlXKwZPFnDGXDkjsGspP+3dlWGwEw+I60j0i2KPuh15k29epsUHudsj4lzGYvTQfInrA9e9Bvyta1H0mIlzapzOX9rmw3+J4YTmrM47x6a6jLN2USaWtBoCuYYGMjI9kZHwkvbp0IDTQSodACxEhAcR0DPIa/1rXxN5CYR6cPmy8M+7AV4aAKwrBzwr9ZsKIm4ya18+1vmxldQ0/HCskLfsM27MKSMkqIKfg4gFHESFWhsdFMLxHJP2jw+geGUxsRAjhwf5OE7d2J7wBpQwhlhca+zXVUH4Wys9A2RkoPg5Hd0Lmt8abimrpOhB6jDMaar0ne3zB6+OF5Rw9W05xeTXFFVXkF1eSnnOWHdkF7D9RjL2EAix+WC2CxU/wt/gRGuhPZIiViJAAIkOsBFkt+PkJ/n6CnwijEzpx5dBL6i1XuxOuxFYNZQXGp6rEWIykMNeoSQsOQ0EmFGRByQlDuI0R0hl6jIWx86FLX4hKhIg4t5yGo3QLD6JbeP3vwS4sryIzv4TcgjJyz5SRX1yJraaG6hpFtU1RVF5FQWkVBaWVHMovprK6BluNwlajqK5RWPykQRE7iveKuKYGbBWGCGpsoGqMrqWqUmOwd3W5sW+rPB+vakDZzH3zu6rMSFddbuRnqzp/nK3SOKb2eFvl+Xxrv6tKjZq0osiMrzDyaYiwGOjU01ikLywaQqIgKBwQED8I6gjBEcbAnJAoCO3m04+Hw4OsDI2NYGhshMdscI6IVz1i/Ni2qvNCUDXG7dRWZfzoVWVGmmo7AdRUGy1vW8WFwquugBoXLrZnCQBLoPGSQT+LIS7xM8Ktwca3fxD4Bxoi69zPmGTpH2QcYw0xxRkBAR2MY8JjICIerPXXWBrX0aSI67x4pv5EGR8bwrVYTQEEgp+ZtcVqCMYabIxv9Q88X/OI5bww/PxNQVnA3xSRJcCIE7/z4dYOEBAC/sHGviXAOLZWiOJ3Ph8/f0NUtYKsFa4P13yai9ENO43Xot/ZoWk3aBFrfB4tYo3Po0Ws8Xma1bATkZNAQ+/F7QzkO8MoF6Btaxmeti1eKdXk6ofNEnGjGYmkONKS9ATatpbhzbbZo90Jjc+jRazxeZwp4ledmJez0ba1DG+27RxO84k1Gk+h3QmNz9NsEYvIDBH5UUQOiMiCeuIDRWSZGf+diCQ4w1An2Xa/iOwRkV0iskZE4r3FNrt0c0REiYjbegUcsU1ErjOvXYaIvOcu2xxCKeXwB7AAB4FeQACwExhYJ8184GVz+3pgWXPKaOnHQdt+CoSY23d5k21mujBgA7AVGOUttgGJwA4g0tzv6g7bHP00tyZOBg4opQ4ppSqB94Fr6qS5BnjL3P4QuFzcM6OwSduUUmuVUqXm7lYg1g12OWSbyf8BzwDuXLbHEdtuB5YopQoAlFJ1luH0LM0VcXcg224/xwyrN41Sqho4C0S11EAn22bPrcAXLrXoPE3aJiIjgDil1GdusqkWR65bX6CviGwSka0iMsNt1jmA905PciEiciMwCpjoaVsARMQPWATM9bApDeGP4VJMwrh7bRCRIUqpMx61yqS5NXEuYD+LMdYMqzeNiPgDHYFTLTXQybYhIlOAR4BZSqlmvMXFpbaFAYOBdSKSCYwFVrqpcefIdcsBViqlqpRSh4F9GKL2DprZCPAHDgE9Od8IGFQnzd1c2LBb7qYGiiO2DcdoxCS6s+HhiG110q/DfQ07R67bDOAtc7szhvsR5c5r2Og5tOCkr8D4Jx4EHjHD/ohRswEEAR8AB4DvgV5uFEtTtn0NHAfSzM9Kb7GtTlq3idjB6yYY7s4eIB243p0ibeqjn9hpfB79xE7j82gRa3weLWKNz6NFrPF5tIg1Po8Wscbn0SLW+Dw+I2IRsYlImojsFpFPRKTJtURFZHMLyokQkfmtzac11GeDm8qdY44B3ykiKSIy3d02tASfedghIsVKqVBz+y1gn1LqCReUkwB8qpQa7Oy87coQjGtf4ykb6inzV8A9wLVKqWMikghsxHhymOMuO1qCz9TEddiCOVxQRG4Uke/NWvoVETn3UgsRKXYgzU3mTI+dIvIO8DTQ20z7Z/t8zO37zbvBbhH5nRmWICJ7ReQ1c+bDahEJtjfYTPOjiLwN7AbiRORjEUk1j7nDLnl9NjR4Dk6wrYNZ5nVKqWMASqn9GI+/L2/OD+MRPP3cuxnP94vNbwvG2IwZwADgE8Bqxr0I3GR/TGNpgEEYYwY6m/udgARgdwNlj8QYO9ABCAUyMAYVJQDVQJKZbjlwY508EoAaYKxdWCfzOxhD2FF2aXfbpWv0PJ1g263Ax/Vc87eAuzz92zf18aXxxMEikoZRA+8FvsKYYjQS2GZOHgkG6s46uLyRNJOBD5RS+QBKqdMiEt6IDZcB/1JKlQCIyEfABGAlcFgplWamS8UQT12ylFJb7fbvFZFrze04jOGN9Q1bbewcnGHbYIzRa3UZBiytJ9yr8CURlymlkkQkBFiFMeRTYQwRfLiR48SBNM7AfmyyDUNodSk5Z5TIJGAKME4pVSoi6zBGANZHa8+hKdsKMYZhni9QZBxGrb6+hWW6DZ/ziZUxR+5e4AGMC/xzEekKICKd6pnBvKaRNN8AvxCRqNo4oAhjkHp9bARmi0iI6Udea4a1hI5AgSng/hgD4Wupa0Nj5+AM2z4DrhORLmb+fYG/A7eoBhqf3oQv1cTnUErtEJFdGLe7R4HV5hSfKowaOut8UrVHROpNo5TKEJEngPUiYgN2KKXmmnPJdgNfKKUetCt3u4gsxRgnDfB305aEFpzGl8A8EdkL/IgxcbW2nFN1bWjoHJxhm1LqexFZCHwtIoEY7Y6blFJbWnBebsdnutiai1m7bldKuW1tibaAiHTDuEP9Uim1y9P2OIJP1sRNISIxGN1Dz3rYFJ9DKXUco9fGZ2izNbGm/eBzDTuNpi5axBqfR4tY4/NoEWt8Hi1ijc+jRazxebSINT6PFrHG5/l/WGbQPiz+KXYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 201.6x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# QU is lower/upper adversary\n",
    "UX, Uy = extract(emnist_X, emnist_y, CHARS, num_per_class=N_per_class[1])    \n",
    "\n",
    "QUX = np.vstack([PX_donate, UX]) \n",
    "QUy = np.hstack([Py_donate, Uy]) \n",
    "U_dist = go(QUX, QUy, clf, fig_title=r\"Random Forest (Nat)\", legend=True, y_ticks=False, fig_filename=\"RF_nat.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error rates: P_train, P_test, U, E\n",
      "[0.0, 0.03642857142857143, 0.5359583333333333, 0.9943333333333333]\n"
     ]
    }
   ],
   "source": [
    "print(\"error rates: P_train, P_test, U, E\")\n",
    "print([np.mean(clf.predict(flatten(X))!=y) for X, y in [(PX_train, Py_train), (PX_test, Py_test), (UX, Uy), (EX, Ey)]])"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}