ai-privacy-toolkit/notebooks/membership_inference_diffpriv_nursery.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running membership inference attacks on the Nursery data and defending using differential privacy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we will show how to run black-box membership attacks on both a vanilla model and models trained with differential privacy. The attacks are taken from the ART library's inference module (https://github.com/Trusted-AI/adversarial-robustness-toolbox) and the differentially private model implementation is taken fron the Differential Privacy Library (https://github.com/IBM/differential-privacy-library).\n",
    "\n",
    "This will be demonstrated on the Nursery dataset (original dataset can be found here: https://archive.ics.uci.edu/ml/datasets/nursery). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already preprocessed the dataset such that all categorical features are one-hot encoded, and the data was scaled using sklearn's StandardScaler."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "sys.path.insert(0, os.path.abspath('..'))\n",
    "\n",
    "from art.utils import load_nursery\n",
    "\n",
    "train_size = 200\n",
    "(x_train, y_train), (x_test, y_test), _, _ = load_nursery(test_set=0.5)\n",
    "x_train = x_train[:train_size]\n",
    "y_train = y_train[:train_size]\n",
    "x_test = x_test[:train_size]\n",
    "y_test = y_test[:train_size]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train logistic regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Base model accuracy:  0.91\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from art.estimators.classification.scikitlearn import ScikitlearnLogisticRegression\n",
    "\n",
    "model = LogisticRegression(solver=\"lbfgs\", max_iter=1000)\n",
    "model.fit(x_train, y_train)\n",
    "\n",
    "art_classifier = ScikitlearnLogisticRegression(model)\n",
    "base_model_accuracy = model.score(x_test, y_test)\n",
    "\n",
    "print('Base model accuracy: ', base_model_accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Attack\n",
    "### Black-box attack\n",
    "The black-box attack basically trains an additional classifier (called the attack model) to predict the membership status of a sample. It can use as input to the learning process probabilities/logits or losses, depending on the type of model and provided configuration.\n",
    "#### Train attack model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from art.attacks.inference.membership_inference import MembershipInferenceBlackBox\n",
    "\n",
    "attack_train_ratio = 0.5\n",
    "attack_train_size = int(len(x_train) * attack_train_ratio)\n",
    "attack_test_size = int(len(x_test) * attack_train_ratio)\n",
    "\n",
    "attack = MembershipInferenceBlackBox(art_classifier, attack_model_type='rf') \n",
    "\n",
    "#train attack model\n",
    "attack.fit(x_train[:attack_train_size], y_train[:attack_train_size],\n",
    "           x_test[:attack_test_size], y_test[:attack_test_size])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Infer membership and check accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "attack accuracy on training data:  0.5\n",
      "attack accuracy on test data:  0.5800000000000001\n",
      "overall attack accuracy:  0.54\n"
     ]
    }
   ],
   "source": [
    "# infer attacked feature\n",
    "inferred_train = attack.infer(x_train[attack_train_size:], y_train[attack_train_size:])\n",
    "inferred_test = attack.infer(x_test[attack_test_size:], y_test[attack_test_size:])\n",
    "\n",
    "# check accuracy\n",
    "train_acc = np.sum(inferred_train) / len(inferred_train)\n",
    "test_acc = 1 - (np.sum(inferred_test) / len(inferred_test))\n",
    "acc = (train_acc * len(inferred_train) + test_acc * len(inferred_test)) / (len(inferred_train) + len(inferred_test))\n",
    "print('attack accuracy on training data: ', train_acc)\n",
    "print('attack accuracy on test data: ', test_acc)\n",
    "print('overall attack accuracy: ', acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that for 54% of the data, membership status is inferred correctly (a little better than a random coin flip)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "precision and recall:  (0.5434782608695652, 0.5)\n"
     ]
    }
   ],
   "source": [
    "def calc_precision_recall(predicted, actual, positive_value=1):\n",
    "    score = 0  # both predicted and actual are positive\n",
    "    num_positive_predicted = 0  # predicted positive\n",
    "    num_positive_actual = 0  # actual positive\n",
    "    for i in range(len(predicted)):\n",
    "        if predicted[i] == positive_value:\n",
    "            num_positive_predicted += 1\n",
    "        if actual[i] == positive_value:\n",
    "            num_positive_actual += 1\n",
    "        if predicted[i] == actual[i]:\n",
    "            if predicted[i] == positive_value:\n",
    "                score += 1\n",
    "    \n",
    "    if num_positive_predicted == 0:\n",
    "        precision = 1\n",
    "    else:\n",
    "        precision = score / num_positive_predicted  # the fraction of predicted “Yes” responses that are correct\n",
    "    if num_positive_actual == 0:\n",
    "        recall = 1\n",
    "    else:\n",
    "        recall = score / num_positive_actual  # the fraction of “Yes” responses that are predicted correctly\n",
    "\n",
    "    return precision, recall\n",
    "\n",
    "# rule-based\n",
    "print('precision and recall: ', calc_precision_recall(np.concatenate((inferred_train, inferred_test)), \n",
    "                            np.concatenate((np.ones(len(inferred_train)), np.zeros(len(inferred_test))))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train differentially private model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP model accuracy:  0.56\n"
     ]
    }
   ],
   "source": [
    "import diffprivlib.models as dp\n",
    "\n",
    "dp_model = dp.LogisticRegression(epsilon=5, data_norm=5)\n",
    "dp_model.fit(x_train, y_train)\n",
    "# print('norm: ', np.linalg.norm(x_train) )\n",
    "\n",
    "dp_art_model = ScikitlearnLogisticRegression(dp_model)\n",
    "print('DP model accuracy: ', dp_model.score(x_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Black-box attack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "attack accuracy on training data:  0.4\n",
      "attack accuracy on test data:  0.62\n",
      "overall attack accuracy:  0.51\n",
      "precision and recall:  (0.5128205128205128, 0.4)\n"
     ]
    }
   ],
   "source": [
    "dp_attack = MembershipInferenceBlackBox(dp_art_model, attack_model_type='rf')\n",
    "\n",
    "# train attack model\n",
    "dp_attack.fit(x_train[:attack_train_size].astype(np.float32), y_train[:attack_train_size],\n",
    "              x_test[:attack_test_size].astype(np.float32), y_test[:attack_test_size])\n",
    "\n",
    "# infer \n",
    "dp_inferred_train = dp_attack.infer(x_train.astype(np.float32)[attack_train_size:], y_train[attack_train_size:])\n",
    "dp_inferred_test = dp_attack.infer(x_test.astype(np.float32)[attack_test_size:], y_test[attack_test_size:])\n",
    "\n",
    "# check accuracy\n",
    "dp_train_acc = np.sum(dp_inferred_train) / len(dp_inferred_train)\n",
    "dp_test_acc = 1 - (np.sum(dp_inferred_test) / len(dp_inferred_test))\n",
    "dp_acc = (dp_train_acc * len(dp_inferred_train) + dp_test_acc * len(dp_inferred_test)) / (len(dp_inferred_train) + len(dp_inferred_test))\n",
    "print('attack accuracy on training data: ', dp_train_acc)\n",
    "print('attack accuracy on test data: ', dp_test_acc)\n",
    "print('overall attack accuracy: ', dp_acc)\n",
    "\n",
    "print('precision and recall: ', calc_precision_recall(np.concatenate((dp_inferred_train, dp_inferred_test)), \n",
    "                            np.concatenate((np.ones(len(dp_inferred_train)), np.zeros(len(dp_inferred_test))))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differentially private model reduces the attack accuracy, but also completely ruins the model's accuracy.\n",
    "Let's see if we can find a sweet spot where both acceptable accuracy and privacy can be acheived."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01\n",
      "0.1\n",
      "0.5\n",
      "1.0\n",
      "5.0\n",
      "10.0\n",
      "25.0\n",
      "75.0\n",
      "100.0\n",
      "200.0\n"
     ]
    }
   ],
   "source": [
    "accuracy = []\n",
    "attack_accuracy = []\n",
    "epsilons = [0.01, 0.1, 0.5, 1.0, 5.0, 10.0, 25.0, 75.0, 100.0, 200.0]\n",
    "\n",
    "for eps in epsilons:\n",
    "    print(eps)\n",
    "    dp_clf = dp.LogisticRegression(epsilon=eps, data_norm=5)\n",
    "    dp_clf.fit(x_train, y_train)\n",
    "    accuracy.append(dp_clf.score(x_test, y_test))\n",
    "    dp_art_classifier = ScikitlearnLogisticRegression(dp_clf)\n",
    "    dp_attack = MembershipInferenceBlackBox(dp_art_classifier, attack_model_type='rf')\n",
    "    dp_attack.fit(x_train[:attack_train_size].astype(np.float32), y_train[:attack_train_size].astype(np.float32),\n",
    "                  x_test[:attack_test_size].astype(np.float32), y_test[:attack_test_size].astype(np.float32))\n",
    "    dp_inferred_train = dp_attack.infer(x_train.astype(np.float32)[attack_train_size:], y_train.astype(np.float32)[attack_train_size:])\n",
    "    dp_inferred_test = dp_attack.infer(x_test.astype(np.float32)[attack_train_size:], y_test.astype(np.float32)[attack_train_size:])\n",
    "    dp_train_acc = np.sum(dp_inferred_train) / len(dp_inferred_train)\n",
    "    dp_test_acc = 1 - (np.sum(dp_inferred_test) / len(dp_inferred_test))\n",
    "    dp_acc = (dp_train_acc * len(dp_inferred_train) + dp_test_acc * len(dp_inferred_test)) / (len(dp_inferred_train) + len(dp_inferred_test))\n",
    "    attack_accuracy.append(dp_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3xcdZ3/8dcnSZO2adP7NU3TAm2hrdCWFnBRFq3K/aIiVFERLywKruyqC15QXHdX0VVXV1ysisKitKCIVUHwxypekeZSoC23Utokbeg1adJbLpPP749zkk6mk8xJOtOZJO/n4zGPzLnMmc+cmXw/53y/3/M95u6IiIj0Ji/bAYiISO5TshARkZSULEREJCUlCxERSUnJQkREUlKyEBGRlJQscpiZ3Wlmt8ZNf9jMdpjZfjObYGZnm9lL4fTl2Yw1jG+DmZ0bcV03s5PC5z8ys3/LaHBH3rfbPs2WdH1mM7vazB7r52sjf18DmZl92sy+n+04BjrTdRbZYWZbgClAOxADNgL3ACvdvSPJ+sOAJuAsd386nPc4sMbdv3m84o6L50dAnbt/tp+vd2COu2861m1lQxo+/zG9PtffTwYfnVlk1yXuPhooB74M3Az8oId1pwDDgQ1x88oTpiMzs4L+vG4gM7P8bMcw2KT7dzQUf5cDhrvrkYUHsAV4U8K8M4AOYGE4/SPg34C5wAHAgf3A/wEvh+seCucVAWMIkk09sC18bX64rfcBfwa+AewNlxUB/wnUADuAO4ER4frnAnXAx4Gd4TavDZddB7QBreF7/zLxM4Wf5a9AY/jabwOFcZ/VgZPiP2f4fD1BEu1cbxiwG1iUZB92xvjpcJ0twNVxy38E/A/wcLj/3pTwXs8BF8etXxBuZ0k4/QDwKrAP+AOwIMXnnw78DNgFvAL8Yy/ff1cc4fSHgE3hd7MGmB637C3AC2Ec3wGeAD4Y973+KXxu4fe7M1z3GWBhxO8rP9yPLwPNQCVQliTuWeF39wGC380fwvnvD/dnA/AoUN6H+Pvyu5wI/Irgd7UX+COQFy67meB33xy+3/Jw/m3AvXHxXEpwkNUI/B44JeH/8hPhvtsHrAaGZ7u8yIWHzixyiLs/RVD4vT5h/ovAgnByrLu/0d1PJPhnusTdR7l7C3A3QbXWScBign/SD8Zt6kxgMzAZ+HfgdoJEtCh8TSnwubj1pxIkoFKCwuEOMxvn7iuBHwNfCd/7kiQfJwb8E8E/92uB5cBHIuyGe4B3x01fCNS7+7oe1p8avkcpcA2w0szmxS1/V/hZRwN/SnjtfcA746bPA3a7e1U4/Qgwh2B/VRF8ZpJ9fjPLA34JPB3Gshy4yczOS/WBzeyNwJeAK4FpwFZgVbhsIvBT4FPABIJC8O962NRbgHMIvtOxwFXAnojf1z+H++JCoISg8D/YS9h/D5wCnBe2l30aeBswiaAAv68P8ffld/lxgv+RSQRn258GPPzObwSWeXC2fh5Bwd+Nmc0NY7sp3MbDwC/NrDButSuB84HZwKkECW3IU7LIPduB8X19kZlNAS4AbnL3A+6+k+BobUX8tt39v929HThMcDT7T+6+192bgf9IWL8N+Fd3b3P3hwmOSuML4h65e6W7P+nu7e6+BfguQQGTyr3AhWZWEk6/B/jfFK+51d1b3P0J4NcE/+ydfuHuf3b3Dnc/nPC6nwCXmtnIcPpd4bzOz3CXuzeHifg24DQzG9NDDMuASe7+r+7e6u6bge/RfX/25GrgLnevCt/rU8BrzWwWQeG9wd0fDL+3bxGc7STTRpAUTyZoj3zO3esjvD8EBxWfdfcXPPC0u+/pZf3bwt/ZIeAfgC+F79dO8DtaZGblEePvy++yjSChloe/yz+6uxMcnBQB881smLtvcfeXk8R9FfBrd/+tu7cRnMGMoHsC+5a7b3f3vQQHAIui7MDBTski95QSnF73VTlBlU29mTWaWSNBAT05bp3auOeTgJFAZdz6vwnnd9oT/gN3OgiMihKMmc01s1+Z2atm1kTwDz8x1evcfTtBtcTbzWwsQQL8cS8vaXD3A3HTWwmqgzrV0gN330RQdXJJmDAuJUwWZpZvZl82s5fD+LeEL+vpM5QD0zv3Zbg/P01w9JvK9DDuzrj2A3sIfgvT4z9DWDDW9fB5/o+guu8OYIeZrYxLuqmUEVRBRRW/X8uBb8Z97r0EVWJR4+/L7/KrBNV1j5nZZjO7JdzuJoKzhduAnWa2yszifwedEvd1R/j+pXHrxCezyL/5wU7JIoeY2TKCH21idUkUtUALMNHdx4aPEndfELdOfNe33QTtHQvi1h/j7lH/MVJ1o/sf4HmCHk8lBAWnRdz23QRVUe8A/uru23pZd5yZFcdNzyQ4O4saZ2dV1GXAxrDQgeAs4zKCdo4xBHX1cOQzJG63Fnglbl+OdffR7n5hivcnjLe8cyL8PBMI6t/rgRlxyyx+OpG7f8vdTyeotpwLfLKHeBPVAidGiLXrrRJe+w8Jn32Eu/8lYvyRf5fhmd7H3f0E4BLgn81sebjsJ+7+OoJ96QTVWYkS97URJMrefmOCkkVOMLMSM7uYoJ76Xnd/tq/bCKsbHgO+Fm4vz8xONLOkVT/hEdX3gG+Y2eQwjtIodeyhHcAJvSwfTdDVd7+ZnQx8OOpnAR4ClgAfI2jDSOULZlZoZq8HLiZomI5qFUFd/4eJq4IiiL+F4Ah/JMGZUbzEz/8U0GRmN5vZiPDMZGF4AJDKT4BrzWyRmRWF7/W3sPru18BrzOzysKfQDQTtNEcxs2VmdmbYzfoAQZVOrId4E30f+KKZzbHAqWY2IULsEDRAf8rMFoRxjDGzd4TLIscPqX+XZnaxmZ0UFvJN4eeLmdk8M3tjuP8OEyScWJK3uB+4yMyWh/vp4wTf818iftYhS8kiu35pZs0ER2afAb4OXHsM23svUEhwzUYDQcPitF7Wv5nglP7JsKrl/xGxTYKg19X8sKrgoSTLP0FwdN5M8M+/OuJ2CevBf0bQwPhgitVfJfis2wmqq6539+f78F71BL22/i4hxnsIqiu2EezPJxNe2u3zu3uM4Eh3EUFPqN0EBXBPbRzxMTwO3ErwmesJjvBXhMt2E5xhfYUgcc0HKggKuEQlBPu6IYx9D0Gd/FHxJnnt1wkK0scICuEfENTlp+TuPyc4il8V/o7WE1Qf9jX+Tr39LueE0/sJvrfvuPvvCdorvkyw318lqH79dJJYXyA4a/3vcN1LCDqJtEb5rEOZLsqTnGRmnwPmuvu7e1nnXIIzsR6rZQabsNdVHUEX4d9lO56+GujxD2U6s5CcY2bjCbrqrsx2LLnAzM4zs7FhFUtn20/imU7OGujxSyBjycLM7jKznWa2voflZmbfMrNNZvaMmS3JVCwycJjZhwiq5R5x9z9kO54c8VqCnkqd1SaXh1V1A8VAj1/IYDWUmZ1DUK94j7svTLL8QuCjBP2wzwS+6e5nZiQYERE5Jhk7swiPCnu7XuAygkTi7v4kMNbMemuMFRGRLMnmoF2ldL8Ypy6cd9QVp2Z2HcH4NhQXF59+8sknH5cARUQGi8rKyt3uPin1msllM1kku0AraZ2YB2PbrARYunSpV1RUZDIuEZFBx8y2pl6rZ9nsDVVHcOVkpxl0v/JWRERyRDaTxRrgvWGvqLOAfX0Y9ExERI6jjFVDmdl9BPcbmGhmdcDnCQa6w93vJBga+EKCKzUPcmxXLouISAZlLFm4+ztTLHeCcWJERCTH6QpuERFJSclCRERSUrIQEZGUlCxERCQlJQsREUlJyUJERFJSshARkZSULEREJKVsDiQoIiJp5O60xjpobQ8ebTEPnsc6jnnbShYiIn3UHgsK4LZ2pyUW614wh4Vz59+2hOlg3fgCvYOWpPOdloRtxC9P9h5tsczczA6ULEQkR3V0ePcCMiycW2MxWtv9qIK3pT1h3WTzO5/HOuK2Eeu5oE9SOLfFOuhIc5lcWJBHUX4ehQV5DAv/xj8vys9j+LA8SoYXdFtemH/k77CE6cRtXHb7scWoZCEyhLk7sQ4/6si3exVG8sK5NbEA7pzX0/wkhfNRBXrce6T7KLkgz7oXxgl/h+UHy0sKh1EYPi/Mz0taOA9LKKSTF/RGUUEehfn5DCuwrm0VJRTiBXmGWbLb++QWJQuR46DzKDmxWqIt1llQercCNrH6odv8uEK2s4A9so3uVSItSbbR7Ug51oGnsUw2IyhQkxzZxhesnUfJ3Y6eO58nFLxFcc+7F/RGYX5+t4K+KMn7dcaTl5f7BXIuG3jJYvdL8MOLgufX/vrI/M55mj9o5jdc+XMeWreNXz1Tz+d3f7Jr/r+M/lLX8680fypn5rvD7c2fosMdBz5gt3UVzPcW/GvX+itab+16vqrwi/2ab8C13NZVgH63/XPkmWEGnx17e1cB+uWGW7rmf6f8m10F6Yde/ihmkGfGQ6et7CpUL6r8YNf8p8793671Fz1+NXkGhlF3+QNdBfrkn709mG9GxzW/OnKUnAO/H81PMv8YDLxkIYOa4zQdamdn82Eu+o/HaY11sLC0hKJh+V3rlE8Y2fW86HBuzS9uzcfMyDO4bO70rsJ2xsYRXfO/cMaCriPfk/46ijyCwvzet5zZdTQ+5+ExXYXwX698Y1fhXPyTO7Cw0N547flHdtwPv9H19JfXvi5ufknX0/9asThufnHX048un3Nk/kvDu55e+JppR+b/eVjX05OnHtkm+Ud63+fnqyf+YGaeznPQ40D34B6c6vcd4qcVdayuqKWu4RBjRgzjrYtLuWpZGadMK0m9ARHplZlVuvvS/r5eZxaSNW2xDh5/bier19bwxIu76HA4+6QJfPK8eZy3YCrD484mRCS7lCzkuNu8az+rK2r5WWUdu/e3MqWkiI+cexJXLi1jZlyVjojkDiULOS4OtcZ4+Nl6Vq+t5akte8nPM5afPJkVZ5RxzpxJFKi+WySnKVlIxrg767c1sWptDWvWbae5pZ3ZE4u55YKTeduSUiaPHp56IyKSE5QsJO32HWzjoXXbWLW2lufqmxg+LI8LF07jqmVlnDF7/IC4AElEulOykLTo6HCefGUP96+t5eH1r9LaHnR5/eLlC7n0tOmMGTEs9UZEJGcpWcgx2dF0mJ9W1nF/RS1b9xxk9PACViwr48qlZSwsHZPt8EQkTZQspM/aYx387oVdrF5bw+9e2EWswznrhPHc9KY5XLBwmrq8igxCShYS2ZbdB7q6vO5sbmHS6CL+4ZwTuHJpGbMmFqfegIgMWEoW0qvDbTF+s/5VVq2t4cnNQZfXN8ybxFXLZvKGeeryKjJUKFlIUuu37eP+iloeqt5G0+F2yieM5JPnzeOK02cwpURdXkWGGiUL6bLvUBtrnt7O6rU1rN/WRGFBHhcunMpVy2Zy5uzxGuJZZAhTshji3J2nXtnL6rW1PLy+nsNtHZwyrYQvXLqAyxeVMmakuryKiJLFkLWz+TAPVm3j/rW1bN59gNFFBbx9yQxWLJvJwtISXTgnIt0oWQwh7bEO/vDSLlY9Vcvjz+8k1uGcMWs8N7zhJC58zTRGFKrLq4gkp2QxBNTsOcj9FbU8UFnLjqYWJo4q5IOvn82VS8s4cdKobIcnIgOAksUgdbgtxqMbXuX+ilr+vGkPeQbnzpvMFy4tY/kpkxmmLq8i0gdKFoPMc/VNrF5by8+rt7HvUBtl40fw8TfP5YqlM5g2ZkS2wxORAUrJYhBoPtzGL5+uZ/XaGp6u20dhfh7nLZzKimVlvPaECeryKiLHLKPJwszOB74J5APfd/cvJywfA9wLzAxj+U93/2EmYxos3J3KrQ2sWlvLr5+p51BbjHlTRvP5S+Zz+aJSxhUXZjtEERlEMpYszCwfuAN4M1AHrDWzNe6+MW61G4CN7n6JmU0CXjCzH7t7a6biGuh272/hwao6Vq+t5eVdByguzOfyxdO5atlMTpsxRl1eRSQjMnlmcQawyd03A5jZKuAyID5ZODDaghJuFLAXaM9gTANSrMP5w0u7uH9tLb/duIP2Duf08nF85YoTueg10yguUm2iiGRWJkuZUqA2broOODNhnW8Da4DtwGjgKnfvSNyQmV0HXAcwc+bMjASbi2r3HuSByjoeqKilft9hxhcXcu3Zs7hqWRknTR6d7fBEZAjJZLJIVh/iCdPnAeuANwInAr81sz+6e1O3F7mvBFYCLF26NHEbg0pLe4zfbtzB6rW1/GnTbgDOmTOJz108n+WnTKGwQF1eReT4y2SyqAPK4qZnEJxBxLsW+LK7O7DJzF4BTgaeymBcOenFHc2sXlvLg1V1NBxso3TsCD62fA7vWFpG6Vh1eRWR7MpkslgLzDGz2cA2YAXwroR1aoDlwB/NbAowD9icwZhyyv6Wdn719HZWV9RSXdPIsHzjLfOnctWyMs4+aSL56vIqIjkiY8nC3dvN7EbgUYKus3e5+wYzuz5cfifwReBHZvYsQbXVze6+O1Mx5QJ3p7q2kdVP1fLLZ7ZzsDXGnMmj+OxFp/DWxaVMGFWU7RBFRI6S0W407v4w8HDCvDvjnm8H3pLJGHLF3gOtXV1eX9q5n5GF+Vxy6nSuXFbGkplj1eVVRHKa+lxmUEeH86dNu1ldUctjG16lLeYsnjmWL7/tNVx82nRGqcuriAwQKq0yYGfzYe77Wy33V9SyrfEQ40YO4z1nBV1e501Vl1cRGXiULDLgfXetZWN9E6+fM5FPXXgyb54/haIC3StCRAYuJYs027O/hY31TXzyvHnc8IaTsh2OiEha6AqvNFtX2wjA0vJxWY5ERCR9lCzSrLqmkfw849QZY7MdiohI2ihZpFlVTQOnTBut+1mLyKCiZJFGsQ7n6dpGFpepCkpEBhclizR6aWczB1pjLClXFZSIDC5KFmlUtTVo3NaZhYgMNkoWaVRd08D44kLKJ4zMdigiImmlZJFG1bWNLC7TOE8iMvgoWaTJvoNtbNq5n8Uz1V4hIoOPkkWarKsL2ytmqr1CRAYfJYs0qa5pwAxOK9OZhYgMPkoWaVJV08i8KaM17LiIDEpKFmnQ0eGsq2lQe4WIDFpKFmmwefcBmg63q71CRAYtJYs0qKppAGCJzixEZJBSskiD6ppGSoYXcMLEUdkORUQkI5Qs0qC6poFFM8eRl6eL8URkcFKyOEb7W9p5cUczi9VlVkQGsZTJwswqzOwGM1PrbRLP1DbS4bBEd8YTkUEsypnFCmA6sNbMVpnZeabBj7p0Nm4v0p3xRGQQS5ks3H2Tu38GmAv8BLgLqDGzL5jZ+EwHmOuqaxo5cVIxY0YOy3YoIiIZE6nNwsxOBb4GfBX4GXAF0AT8X+ZCy33uTnVtI0t0fYWIDHIpx6Yws0qgEfgBcIu7t4SL/mZmZ2cyuFy3dc9B9h5o1cV4IjLoRRnI6B3uvjnZAnd/W5rjGVCqa4P2Cg3zISKDXZRqqA+aWVdpaGbjzOzfMhjTgFFd00hxYT5zp4zOdigiIhkVJVlc4O6NnRPu3gBcmLmQBo6qmgZOKxtLvi7GE5FBLkqyyDezos4JMxsBFPWy/pBwqDXGc/XNqoISkSEhSpvFvcDjZvZDwIH3A3dnNKoB4Nlt+4h1uHpCiciQkDJZuPtXzOxZYDlgwBfd/dGMR5bjui7G0zAfIjIERLqtm7s/AjyS4VgGlOqaBsonjGTCqCFfIyciQ0CUsaHOMrO1ZrbfzFrNLGZmTVE2bmbnm9kLZrbJzG7pYZ1zzWydmW0wsyf6+gGywd2pqtHFeCIydEQ5s/g2wfhQDwBLgfcCJ6V6kZnlA3cAbwbqCMaWWuPuG+PWGQt8Bzjf3WvMbHLfP8Lxt63xELuaW9S4LSJDRqThPtx9E5Dv7jF3/yHwhggvOwPY5O6b3b0VWAVclrDOu4AH3b0mfJ+d0UPPnuqaoCfx4jKdWYjI0BAlWRw0s0JgnZl9xcz+CSiO8LpSoDZuui6cF28uMM7Mfm9mlWb23mQbMrPrwqHSK3bt2hXhrTOruqaR4cPyOHmaLsYTkaEhSrJ4T7jejcABoAx4e4TXJbtSzROmC4DTgYuA84BbzWzuUS9yX+nuS9196aRJkyK8dWZV1TRwaulYhuXr3lEiMjT02mYRtjv8u7u/GzgMfKEP264jSCydZgDbk6yz290PAAfM7A/AacCLfXif46qlPcbG7U1ce/asbIciInLc9Hpo7O4xYFJYDdVXa4E5ZjY7fP0KYE3COr8AXm9mBWY2EjgTeK4f73XcbNjeRGusQyPNisiQEqU31Bbgz2a2hqAaCgB3/3pvL3L3djO7EXgUyAfucvcNZnZ9uPxOd3/OzH4DPAN0AN939/X9+yjHR9VWjTQrIkNPlGSxPXzkAX1q0XX3h4GHE+bdmTD9VYKbKg0I1bWNlI4dwZSS4dkORUTkuIky3Edf2ikGvXU1jTqrEJEhJ8qd8n7H0b2YcPc3ZiSiHLaj6TDbGg/x/tfNznYoIiLHVZRqqE/EPR9O0G22PTPh5LbqGrVXiMjQFKUaqjJh1p8HyhhO6VZd00hhfh4LppdkOxQRkeMqSjXU+LjJPIKL6KZmLKIcVlXTwILSEooK8rMdiojIcRWlGqqSoM3CCKqfXgE+kMmgclFbrINn6vZx9Znl2Q5FROS4i1INpdZc4Pn6ZlraO1hSrvYKERl6otzP4oZwKPHO6XFm9pHMhpV7qroat3XltogMPVFGwvuQuzd2Trh7A/ChzIWUm6prGphSUsT0MboYT0SGnijJIs/MukaQDQcX7M9YUQNadW0ji8vGEbcrRESGjCjJ4lHgfjNbbmZvBO4DfpPZsHLL7v0tbN1zUNdXiMiQFaU31M3AdcCHCXpEPQZ8P5NB5Zp14Z3xlpSrvUJEhqYoyWIE8L3OAQDDaqgi4GAmA8slVTUNFOQZC6ePyXYoIiJZEaUa6nGChNFpBPD/MhNObqquaeSUaSWMKNTFeCIyNEVJFsPdfX/nRPh8ZOZCyi2xDufpukaWqL1CRIawKMnigJkt6Zwws9OBQ5kLKbe88GozB1tjur5CRIa0KG0WNwEPmFnn/bOnAVdlLqTcUl2rkWZFRKIM97HWzE4G5hH0hnre3dsyHlmOqK5pZEJxITPHD5maNxGRo0Q5s4AgUcwnuJ/FYjPD3e/JXFi5o6qmgcUzx+piPBEZ0qIMUf554FyCZPEwcAHwJ2DQJ4vGg61s3nWAty+Zke1QRESyKkoD9xXAcuBVd78WOI3gOotBb11tcDGe2itEZKiLkiwOuXsH0G5mJcBO4ITMhpUbqmoayTM4dYaShYgMbVHaLCrCIcq/R3AjpP3AUxmNKkdU1zQwd8poRhVFbdoRERmcovSG6rx3xZ1m9hugxN2fyWxY2dfR4ayrbeSS06ZnOxQRkazr0yGzu2/JUBw55+Vd+2k+3M7iMlVBiYhEabMYkqprOhu3deW2iIiSRQ+qaxsYM2IYJ0wsznYoIiJZ12M1lJmN7+2F7r43/eHkjqqtjSwqG0teni7GExHprc2iEnCCIT4SOYO4+2zz4TZe3NnMBa+Zmu1QRERyQo/Jwt1nH89A0qH5cBsPVW/j6jPLj+mM4Jm6fbjDErVXiIgAEdosLPBuM7s1nJ5pZmdkPrS++2llHbf+YkPXSLH9VbU1eP1p6gklIgJEa+D+DvBa4F3hdDNwR8YiOgYVYSG/YXvTMW2nuraRkyaPYsyIYekIS0RkwIuSLM509xuAwwDu3gAUZjSqfuo8I1i/bV+/t+HuVNc06M54IiJxoiSLNjPLJ2jUxswmAR0Zjaoftjceon7fYeDYziy27DlIw8E2XV8hIhInSrL4FvBzYLKZ/TvB8OT/kdGo+qGzCup1J03kxR3NtLb3L59V1+jOeCIiiVImC3f/MfAvwJeAeuByd38gysbN7Hwze8HMNpnZLb2st8zMYmZ2RdTAE1VtbWBkYT5vP72Utpjz0s7mfm2nuqaRUUUFzJk8ur+hiIgMOlEvytsJ3Be/LNVFeWHV1R3Am4E6YK2ZrXH3jUnWux14tO/hH1GxdS+LysZ2DSe+YVsTC6aP6fN2qmoaOK1sDPm6GE9EpEtvZxaVQEX4dxfwIvBS+LwywrbPADa5+2Z3bwVWAZclWe+jwM8IElK/HGhp57n6Zk4vH8fsCcUUF+azYXvfG7kPtrbz/KvNLC5Te4WISLwek4W7z3b3EwiO+C9x94nuPgG4GHgwwrZLgdq46bpwXhczKwXeCtzZ24bM7DozqzCzil27dh21/OnaRmIdzunl48jLM06ZVtKvRu5n6/YR63CWlKu9QkQkXpQG7mXu/nDnhLs/Avx9hNf1NExIvP8Cbnb3WG8bcveV7r7U3ZdOmjTpqOWVWxswOzJC7ILpJWysbyLWkfh2vasKR5pdpDMLEZFuoiSL3Wb2WTObZWblZvYZYE+E19UBZXHTM4DtCessBVaZ2RaCe31/x8wuj7Dtbiq2NjB38uiui+gWlI7hYGuMLXsO9Gk71TUNzJ5YzPjinLyMREQka6Iki3cCkwi6zz4ETA7npbIWmGNms82sEFgBrIlfIazqmuXus4CfAh9x94f6ED8dHU5VTQNLyo+cDSyYXgL07XoLd6eqplE3OxIRSSLKbVX3Ah8zsxKgw933R9mwu7eb2Y0EbR75wF3uvsHMrg+X99pOEdWm8I52p8clizmTRzMs39iwfR+XRrwtal3DIXbvb9H1FSIiSaRMFmb2GuAeYHw4vRu4xt3Xp3pt2NbxcMK8pEnC3d8XId6jVGwJLqJbGpcsCgvymDd1NBu2RT+zqK7VnfFERHoSpRrqu8A/u3u5u5cDHwdWZjas6Cq3NjChuJDyCSO7zV8wbQwbtu/DPVojd9XWBoYPy+PkqboYT0QkUZRkUezuv+uccPffAzlzr9HKrXs5vXwcZt07Xy0oLaHhYFvXeFGpVNc2cuqMsRTk606zIiKJopSMm83s1rA31Cwz+yzwSqYDi2L3/ha27DnYrb2iU18auQ+3xdi4fZ9udiQi0oMoyeL9BL2hHiToETUJuDaTQUVVGQ4emCxZnDKtBLNow5Vv2L6PtpircVtEpAdRekM1AP94HGLps6qtDRTm57Gw9OgxoEYWFnDCxOJIZxbVNZ2N20oWIn7c+vgAAAymSURBVCLJ9DaQ4JqelgG4+6XpDye1l3cd6blbubWBhaUlDB+Wn3TdBdPHULGl1/EOgSBZzBg3gsmjh6ctThGRwaS3M4vXEoztdB/wN5IP33HcdY7g0dIe45lt+3jf383qcd2FpSWseXo7ew+09npVdlVNA0tnje9xuYjIUNdbm8VU4NPAQuCbBEON73b3J9z9ieMRXG/Wb2uitb2j10bpziHKexuBtn5fcIc9XbktItKz3kadjbn7b9z9GuAsYBPwezP76HGLLonDbTH++NIuKrcG1UvJGrc7RekRtS5sr1jSy3ZERIa6Xhu4zawIuIhgLKhZBLdYjTI8eUa95wdPcd6CKZRPGMmk0UU9rjd2ZCGlY0f0miyqahooLMhj/rSSTIQqIjIo9NbAfTdBFdQjwBeiDO9xPFVubeCcOUcPV55owfQSNvTSfba6ppGF00soLNDFeCIiPemthHwPMBf4GPAXM2sKH81m1vc7C6XZ7v2tnD4rddXRguljeGXPAQ60tB+1rLW9g2e36WI8EZFUejyzcPecP9Turb2i04LpJbjDc/VNR/V4eq6+iZb2Dg0eKCKSQs4nhJ6MLipg7uTUg/51XrCXrN2iuia4AlwX44mI9G7AJovF4f22U5lSUsSE4sKkw35U1zYytWQ408eOyESIIiKDxoBNFmedEO0iOjNj/vSSpGcWVTUNOqsQEYlgQCaLd55RxvvPnh15/QXTx/DSzmZa2mNd83Y1t1C795CShYhIBAMyWSwsHdPjeFDJ1y+hLea8tOPIuFLrwjvjqSeUiEhqAzJZDOvjDYqSDftRVdNAQZ4lHbFWRES6G5DJoqiPF9CVjx/JqKKCbu0W1TUNzJ/e84i1IiJyxIBMFoV9PLPIyzPmTzvSyN0e6+CZOl2MJyIS1YBMFn2thgKYP72EjdubiHU4L+xo5mBrTI3bIiIRDchk0Z9xnBZML+FQW4xXdh84cme8Mp1ZiIhEMSCTRX/OLOIbuatrGpk4qpCy8boYT0QkigGZLPpzZjFnyigK8/PYsL2J6poGFpWNwywnbv4nIpLzBmSy6GtvKAjORuZNHc2fN+1m8+4Daq8QEemDAZks+lMNBeG9LcIeUeoJJSIS3YBMFv29UdGC8AK8PINTZ+hiPBGRqAZkshiW37+2hs57cp88tYTiol7vKCsiInEGZLLo75nFKVNLKMgzlpSrvUJEpC8G5OF1UX7/hugYUZjP3e8/gzmTR6U5IhGRwW1AJothBf3v8nr2SRPTGImIyNAwMKuh+tkbSkRE+mdAlrr5EW6nKiIi6ZPRZGFm55vZC2a2ycxuSbL8ajN7Jnz8xcxOi7jd9AcrIiI9yliyMLN84A7gAmA+8E4zm5+w2ivA37v7qcAXgZWptltcOCCbWUREBrRMnlmcAWxy983u3gqsAi6LX8Hd/+LuDeHkk8CMVBs9YVJx2gMVEZHeZTJZlAK1cdN14byefAB4JNkCM7vOzCrMrGLXrl1pDFFERKLIZLJI1rDgSVc0ewNBsrg52XJ3X+nuS9196aRJk9IYooiIRJHJBoA6oCxuegawPXElMzsV+D5wgbvvyWA8IiLST5k8s1gLzDGz2WZWCKwA1sSvYGYzgQeB97j7ixmMRUREjkHGzizcvd3MbgQeBfKBu9x9g5ldHy6/E/gcMAH4Ttgdtt3dl2YqJhER6R9zT9qMkLOWLl3qFRUV2Q5DRGRAMbPKYzkYH5BXcIuIyPGlZCEiIikpWYiISEpKFiIikpKShYiIpKRkISIiKSlZiIhISkoWIiKSkpKFiIikpGQhIiIpKVmIiEhKShYiIpKSkoWIiKSkZCEiIikpWYiISEpKFiIikpKShYiIpKRkISIiKSlZiIhISkoWIiKSkpKFiIikpGQhIiIpKVmIiEhKShYiIpKSkoWIiKSkZCEiIikpWYiISEpKFiIikpKShYiIpKRkISIiKSlZiIhISkoWIiKSkpKFiIikpGQhIiIpKVmIiEhKGU0WZna+mb1gZpvM7JYky83MvhUuf8bMlmQyHhER6Z+MJQszywfuAC4A5gPvNLP5CatdAMwJH9cB/5OpeEREpP8yeWZxBrDJ3Te7eyuwCrgsYZ3LgHs88CQw1symZTAmERHph4IMbrsUqI2brgPOjLBOKVAfv5KZXUdw5gGw38xeSG+oGTER2J3tICJQnOk1EOIcCDGC4ky3ecfy4kwmC0syz/uxDu6+EliZjqCOFzOrcPel2Y4jFcWZXgMhzoEQIyjOdDOzimN5fSaroeqAsrjpGcD2fqwjIiJZlslksRaYY2azzawQWAGsSVhnDfDesFfUWcA+d69P3JCIiGRXxqqh3L3dzG4EHgXygbvcfYOZXR8uvxN4GLgQ2AQcBK7NVDxZMFCqzRRneg2EOAdCjKA40+2Y4jT3o5oIREREutEV3CIikpKShYiIpKRkcYzMrMzMfmdmz5nZBjP7WDj/NjPbZmbrwseFORDrFjN7NoynIpw33sx+a2YvhX/HZTnGeXH7bJ2ZNZnZTbmwP83sLjPbaWbr4+b1uP/M7FPhUDYvmNl5WY7zq2b2fDiszs/NbGw4f5aZHYrbr3dmOc4ev+cc25+r42LcYmbrwvlZ2Z+9lEPp+326ux7H8ACmAUvC56OBFwmGN7kN+ES240uIdQswMWHeV4Bbwue3ALdnO8642PKBV4HyXNifwDnAEmB9qv0X/gaeBoqA2cDLQH4W43wLUBA+vz0uzlnx6+XA/kz6Pefa/kxY/jXgc9ncn72UQ2n7ferM4hi5e727V4XPm4HnCK5CHyguA+4On98NXJ7FWBItB152963ZDgTA3f8A7E2Y3dP+uwxY5e4t7v4KQY+/M7IVp7s/5u7t4eSTBNc0ZVUP+7MnObU/O5mZAVcC9x2PWHrSSzmUtt+nkkUamdksYDHwt3DWjeFp/13Zrt4JOfCYmVWGQ6gATPHw2pbw7+SsRXe0FXT/J8y1/Qk977+ehrLJBe8HHombnm1m1Wb2hJm9PltBxUn2Pefq/nw9sMPdX4qbl9X9mVAOpe33qWSRJmY2CvgZcJO7NxGMoHsisIhgrKuvZTG8Tme7+xKC0X5vMLNzsh1QTyy4kPNS4IFwVi7uz95EGsrmeDOzzwDtwI/DWfXATHdfDPwz8BMzK8lWfPT8Pefk/gTeSfcDmqzuzyTlUI+rJpnX6/5UskgDMxtG8AX92N0fBHD3He4ec/cO4Hscp1Pm3rj79vDvTuDnBDHtsHCk3/DvzuxF2M0FQJW774Dc3J+hnvZfzg1lY2bXABcDV3tYcR1WQ+wJn1cS1F3PzVaMvXzPubg/C4C3Aas752VzfyYrh0jj71PJ4hiFdZY/AJ5z96/HzY8fav2twPrE1x5PZlZsZqM7nxM0eK4nGHLlmnC1a4BfZCfCo3Q7Ysu1/Rmnp/23BlhhZkVmNpvgni1PZSE+ILgRGXAzcKm7H4ybP8mCe89gZicQxLk5O1H2+j3n1P4MvQl43t3rOmdka3/2VA6Rzt/n8W61H2wP4HUEp2/PAOvCx4XA/wLPhvPXANOyHOcJBL0fngY2AJ8J508AHgdeCv+Oz4F9OhLYA4yJm5f1/UmQvOqBNoIjsw/0tv+AzxAcWb4AXJDlODcR1FF3/kbvDNd9e/h7eBqoAi7Jcpw9fs+5tD/D+T8Crk9YNyv7s5dyKG2/Tw33ISIiKakaSkREUlKyEBGRlJQsREQkJSULERFJSclCRERSUrIQSRMzu9TMbgmf32Zmn8h2TCLpkrHbqooMNe6+hqPvMy8yKOjMQiRkZu82s6fC+xB818zyzWy/mX3NzKrM7HEzmxSu+49mtjEc8G5VOO99ZvbtJNtdZGZP2pF7SYwL5//ezG4P3/PFHBnETyQpJQsRwMxOAa4iGGxxERADrgaKCcaoWgI8AXw+fMktwGJ3PxW4PsXm7wFuDtd9Nm4bENxj4gzgpoT5IjlF1VAigeXA6cDaYJgdRhAMutbBkYHi7gU6B2h7BvixmT0EPNTTRs1sDDDW3Z8IZ93NkZF0idteJcGNc0Ryks4sRAIG3O3ui8LHPHe/Lcl6nePjXATcQZBgKsMRSPujJfwbQwdvksOULEQCjwNXmNlk6Lp3cTnB/8gV4TrvAv5kZnlAmbv/DvgXYCwwKtlG3X0f0BDXHvEeguoskQFFRzIigLtvNLPPEtxJMI9ghNEbgAPAAjOrBPYRtGvkA/eGVUwGfMPdG8Pqq2SuAe40s5EEw1Vfm9lPI5J+GnVWpBdmtt/dk541iAwlqoYSEZGUdGYhIiIp6cxCRERSUrIQEZGUlCxERCQlJQsREUlJyUJERFL6/zZABsCzGKyYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(epsilons, accuracy)\n",
    "plt.plot(epsilons, np.ones_like(epsilons) * base_model_accuracy, dashes=[2,2], label=\"base model\")\n",
    "plt.title(\"Differentially private logistic regression\")\n",
    "plt.xlabel(\"epsilon\")\n",
    "plt.ylabel(\"Model accuracy\")\n",
    "plt.ylim(0, 1)\n",
    "plt.xlim(0.1, 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZgcZbn38e8vk8wkmezJJIQkBBICIUFACCAobmAgyHoOKouKeBTxdQG3w+JhOaKv26seRQRREVAQ9bAryiabCkLCnkAgCYGEQDLZM1lmvd8/qibpmcxMdSbT0z3J73Ndc01X1VNVd1dX1931PE9VKSIwMzPrSK9iB2BmZqXPycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5NFCZN0taSLc4Y/K2mppBpJwyW9U9Ir6fBJxYw1jW+2pPfmWTYk7Zm+vk7SNwsa3Jb1ttimxdJV71nSGZLu7eS8eX9ePZmkiyT9sthx9HTydRbFIWkhMApoABqBOcANwDUR0dRG+T7AWuAdEfFsOu4B4M6I+HF3xZ0Tz3XA4oj4r07OH8CkiJi3vcsqhi54/9s1f6mvz3Y8PrMoruMjYiAwHvgOcD7wq3bKjgL6ArNzxo1vNZw3Sb07M19PJqms2DHsaLp6P9oZ98seIyL8V4Q/YCFwVKtxhwBNwL7p8HXAN4G9gPVAADXA34D5admN6bgKYDBJsnkTeCOdtyxd1ieAfwA/Alam0yqA/we8DiwFrgb6peXfCywGvgIsS5d5VjrtbKAeqEvXfVfr95S+l8eA1em8PwXKc95rAHvmvs/09QskSbS5XB9gOXBAG9uwOcaL0jILgTNypl8HXAXcnW6/o1qt60XguJzyvdPlHJgO/xF4C1gDPAJMzXj/uwK3ANXAq8AXO/j8N8eRDn8amJd+NncCu+ZMmw7MTeP4GfAw8Kmcz/Xv6Wuln++ytOxzwL55fl5l6XacD6wDZgHj2oh79/Sz+w+S/eaRdPwn0+25CrgHGL8N8W/LfjkC+BPJfrUSeBTolU47n2S/X5eu78h0/GXAb3PiOYHkR9Zq4CFgn1bfy6+m224N8Hugb7GPF6Xw5zOLEhIRT5Ac/I5oNf5lYGo6OCQi3h8RE0m+TMdHxICIqAWuJ6nW2hN4O8mX9FM5izoUWACMBL4FfJckER2QzjMGuCSn/C4kCWgMycHhSklDI+Ia4Ebge+m6j2/j7TQCXyL5ch8GHAn8nzw2ww3AR3OGjwXejIhn2im/S7qOMcCZwDWS9s6Zfnr6XgcCf2817++A03KGjwaWR8RT6fBfgEkk2+spkvdMW+9fUi/gLuDZNJYjgfMkHZ31hiW9H/g28GFgNPAacHM6bQTwv8CFwHCSg+Dh7SxqOvBuks90CPARYEWen9eX021xLDCI5OC/oYOw3wPsAxydtpddBPwbUEVyAP/dNsS/LfvlV0i+I1UkZ9sXAZF+5p8HDo7kbP1okgN/C5L2SmM7L13G3cBdkspzin0YOAbYA9iPJKHt9JwsSs8SYNi2ziRpFDADOC8i1kfEMpJfa6fmLjsiroiIBmATya/ZL0XEyohYB/zfVuXrgW9ERH1E3E3yqzT3QNyuiJgVEY9HRENELAR+TnKAyfJb4FhJg9LhjwG/yZjn4oiojYiHgT+TfNmb3RER/4iIpojY1Gq+m4ATJPVPh09PxzW/h2sjYl2aiC8D9pc0uJ0YDgaqIuIbEVEXEQuAX9Bye7bnDODaiHgqXdeFwGGSdic5eM+OiFvTz+0nJGc7baknSYqTSdojX4yIN/NYPyQ/Kv4rIuZG4tmIWNFB+cvS/Wwj8Bng2+n6Gkj2owMkjc8z/m3ZL+tJEur4dL98NCKC5MdJBTBFUp+IWBgR89uI+yPAnyPivoioJzmD6UfLBPaTiFgSEStJfgAckM8G3NE5WZSeMSSn19tqPEmVzZuSVktaTXKAHplTZlHO6yqgPzArp/xf0/HNVqRf4GYbgAH5BCNpL0l/kvSWpLUkX/gRWfNFxBKSaol/lzSEJAHe2MEsqyJifc7wayTVQc0W0Y6ImEdSdXJ8mjBOIE0WksokfUfS/DT+hels7b2H8cCuzdsy3Z4Xkfz6zbJrGndzXDXACpJ9Ydfc95AeGBe3837+RlLddyWwVNI1OUk3yziSKqh85W7X8cCPc973SpIqsXzj35b98vsk1XX3Slog6YJ0ufNIzhYuA5ZJullS7n7QrPW2bkrXPyanTG4yy3uf39E5WZQQSQeT7LStq0vysQioBUZExJD0b1BETM0pk9v1bTlJe8fUnPKDIyLfL0ZWN7qrgJdIejwNIjlwKs9lX09SFfUh4LGIeKODskMlVeYM70ZydpZvnM1VUScCc9KDDiRnGSeStHMMJqmrhy3vofVyFwGv5mzLIRExMCKOzVg/abzjmwfS9zOcpP79TWBszjTlDrcWET+JiINIqi33Ar7WTrytLQIm5hHr5lW1mvczrd57v4j4Z57x571fpmd6X4mICcDxwJclHZlOuyki3kWyLYOkOqu11ttaJImyo33McLIoCZIGSTqOpJ76txHx/LYuI61uuBf4Qbq8XpImSmqz6if9RfUL4EeSRqZxjMmnjj21FJjQwfSBJF19ayRNBj6b73sBbgcOBM4lacPI8t+SyiUdARxH0jCdr5tJ6vo/S04VFEn8tSS/8PuTnBnlav3+nwDWSjpfUr/0zGTf9AdAlpuAsyQdIKkiXde/0uq7PwNvk3RS2lPocyTtNFuRdLCkQ9Nu1utJqnQa24m3tV8Cl0uapMR+kobnETskDdAXSpqaxjFY0ofSaXnHD9n7paTjJO2ZHuTXpu+vUdLekt6fbr9NJAmnsY1V/AH4oKQj0+30FZLP+Z95vtedlpNFcd0laR3JL7OvAz8EztqO5X0cKCe5ZmMVScPi6A7Kn09ySv94WtVyP3m2SZD0upqSVhXc3sb0r5L8Ol9H8uX/fZ7LJa0Hv4WkgfHWjOJvkbzXJSTVVedExEvbsK43SXptHd4qxhtIqiveINmej7eatcX7j4hGkl+6B5D0hFpOcgBur40jN4YHgItJ3vObJL/wT02nLSc5w/oeSeKaAswkOcC1NohkW69KY19BUie/VbxtzPtDkgPpvSQH4V+R1OVniojbSH7F35zuRy+QVB9ua/zNOtovJ6XDNSSf288i4iGS9orvkGz3t0iqXy9qI9a5JGetV6RljyfpJFKXz3vdmfmiPCtJki4B9oqIj3ZQ5r0kZ2LtVsvsaNJeV4tJugg/WOx4tlVPj39n5jMLKzmShpF01b2m2LGUAklHSxqSVrE0t/20PtMpWT09fksULFlIulbSMkkvtDNdkn4iaZ6k5yQdWKhYrOeQ9GmSarm/RMQjxY6nRBxG0lOpudrkpLSqrqfo6fEbBayGkvRuknrFGyJi3zamHwt8gaQf9qHAjyPi0IIEY2Zm26VgZxbpr8KOrhc4kSSRREQ8DgyR1FFjrJmZFUkxb9o1hpYX4yxOx211xamks0nub0NlZeVBkydP7pYAzcx2FLNmzVoeEVXZJdtWzGTR1gVabdaJRXJvm2sApk2bFjNnzixkXGZmOxxJr2WXal8xe0MtJrlystlYWl55a2ZmJaKYyeJO4ONpr6h3AGu24aZnZmbWjQpWDSXpdyTPGxghaTFwKcmN7oiIq0luDXwsyZWaG9i+K5fNzKyACpYsIuK0jOlBcp8YMzMrcb6C28zMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPLVNBkIekYSXMlzZN0QRvTB0u6S9KzkmZLOquQ8ZiZWecULFlIKgOuBGYAU4DTJE1pVexzwJyI2B94L/ADSeWFisnMzDqnkGcWhwDzImJBRNQBNwMntioTwEBJAgYAK4GGAsZkZmadUMhkMQZYlDO8OB2X66fAPsAS4Hng3Ihoar0gSWdLmilpZnV1daHiNTOzdhQyWaiNcdFq+GjgGWBX4ADgp5IGbTVTxDURMS0iplVVVXV9pGZm1qFCJovFwLic4bEkZxC5zgJujcQ84FVgcgFjMjOzTihksngSmCRpj7TR+lTgzlZlXgeOBJA0CtgbWFDAmMzMrBN6F2rBEdEg6fPAPUAZcG1EzJZ0Tjr9auBy4DpJz5NUW50fEcsLFZOZmXVOwZIFQETcDdzdatzVOa+XANMLGYOZmW0/X8FtZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlikzWUi6RdIHJTmxmJntpPJJAFcBpwOvSPqOpMkFjsnMzEpMZrKIiPsj4gzgQGAhcJ+kf0o6S1KfQgdoZmbFl1fVkqThwCeATwFPAz8mSR73FSwyMzMrGb2zCki6FZgM/AY4PiLeTCf9XtLMQgZnZmalITNZAD+NiL+1NSEipnVxPGZmVoLyqYbaR9KQ5gFJQyX9nwLGZGZmJSafZPHpiFjdPBARq4BPFy4kMzMrNfkki16S1DwgqQwoL1xIZmZWavJps7gH+IOkq4EAzgH+WtCozMyspOSTLM4HPgN8FhBwL/DLQgZlZmalJTNZREQTyVXcVxU+HDMzK0X5XGcxCfg2MAXo2zw+IiYUMC4zMysh+TRw/5rkrKIBeB9wA8kFepkkHSNprqR5ki5op8x7JT0jabakh/MN3MzMuk8+yaJfRDwAKCJei4jLgPdnzZT2mroSmEFyVnKapCmtygwBfgacEBFTgQ9tY/xmZtYN8mng3pTenvwVSZ8H3gBG5jHfIcC8iFgAIOlm4ERgTk6Z04FbI+J1gIhYti3Bm5lZ98jnzOI8oD/wReAg4KPAmXnMNwZYlDO8OB2Xay9gqKSHJM2S9PG2FiTpbEkzJc2srq7OY9VmZtaVOjyzSKuSPhwRXwNqgLO2YdlqY1y0sf6DgCOBfsBjkh6PiJdbzBRxDXANwLRp01ovw8zMCqzDZBERjZIOkqSI2NaD9GJgXM7wWGBJG2WWR8R6YL2kR4D9gZcxM7OSkU+bxdPAHZL+CKxvHhkRt2bM9yQwSdIeJO0cp5K0UeS6A/ippN4ktxA5FPhRnrGbmVk3ySdZDANW0LIHVAAdJouIaEgbxO8ByoBrI2K2pHPS6VdHxIuS/go8BzQBv4yIFzrxPszMrIC07bVLxTVt2rSYOdPPXDIz2xaSZm3PM4jyuYL712zdME1EfLKzKzUzs64REazZWM+ydbUsW1vLsnWbtnq9fF3tdq8nn2qoP+W87guczNYN1T1CXUMT59/yHLsPr+Qz75lA3z5lxQ7JzKxNjU3Bipra5MC/blN68K+lel3LhFBdU0tdQ9NW8/cvL2PkwApGDuzLPrsO4sHtjCefGwnekjss6XfA/du53qL4n/tf5ran3wDg1qcXc9kJU3nf3vlcX2hm1jVqGxo3H+SXra2luo0zgWXrallRU0tTG60EQ/r3oWpABSMHVXDIHsMYObCCqoEVjBzUN00OyesBFS0P7z87Y/vizufMorVJwG7bt9ru98SrK7nq4fl8ZNo4jt9/Vy658wXO+vWTHD11FJccP5UxQ/oVO0RL1Tc28frKDcxfVsOC5etZtHIDU3cdzFFTRjJyYN/sBZgVQU1tA8vWbjnYL1u7KT0LaHlmsGZj/Vbz9hIMH1Cx+WC/766DGTmoIk0EfXNeV1DRuzg1IpkN3JLW0bLN4i3gwtZnHN2lMw3c6zbVc8z/PEpZL3H3uUcwoKI3tQ2N/PLRV7nib68gxBeO3JNPvWsC5b3zuajdusKajfUsqK5hfvV65lfXMH9ZDfOra3htxQYacn5SDajoTU1tAxK8fdwQPjBlF6ZPHcXEqgFFjN52BhHBqg31LQ72za9zq4Oq19Wyoa5xq/nLy3qlv/orNlcJjdw83DeZNrCC4QMqKOvV1nXMXWd7G7h3it5QX/nDs9z29GL+eM7hHDR+aItpi1dt4Bt3zeHeOUuZWFXJ5Sfuy+F7jujKkHdqTU3BG6s3JsmgRVJYz/KaLY1uvXuJ3UdUMmFEJRNHDmBi1QAmVlUyoWoAg/r25uWlNdw7+y3unbOU599YA8DEqkqmT92FD0wZxQFjh9CrwF8223E0NDaxvKZu6ySwuU2gluq1m6iuqaW+cetj5ICK3m1W/1Q1J4Q0OQzu14ecp1IXVcGThaSTgb9FxJp0eAjw3oi4vbMr3R7Tdh8cMy99VzJw1p+3TPj1B7e8zhm/4qdH8cqyGsYM6ce4Lz/Ybvm/vbSUy+6cw/dqLmR4ZTnjh1dS/qm/ZC7f4xON1x7LpvpGNtY1csPkn21OCv+96j9p3sVOrbuYwf36MLGqku/XXES/8jL69ilj9YdvY9yw/vQp65XXepec/L/c/+JS7p29lC+8fi4B9CnrxS37/ZzpU0Zx2MThVPzmhJLaPh7fPeObrj2WusYm6hua+Pu7rt+cBP7t2c8k4xub+CSXsWJ9HRFwc/nlm+c9te5ihlWWM3JgBT/a+HX6lPWivKwX9x967eYzgf3uO50+vXtRJpXE+92W8QXvOgtcGhG3NQ9ExGpJlwJFSRbbYunaTSxavp7K8t6MGdpxm8T7J4/i8IkjqL6iH0vWbGT1hnqe+/urnHnYeHqXuWoKklPy+sYmNtY1srG+kavvnM386hoWVK/nBxtWbi53xZJXGDe0PxOrKtmloS99y8vo16eMWR87imGV5ckvrV9vqUIato3VSbsO6cfHD9udjx+2Ow2/GsrqDfWs3FDHHU+/wU3/ep0BFb25tf86hvUvZ0j/Pp1qmLPSEQSNTcHCZTUsW5e0Axy8ZiP1jU3UNTRxyTWPbz4r+EXTlv3wszc+BUBZL/H+inrK04P/B/YelbQDDKxgr1kDKS/rRZ8y8fJZM7ZUQ/960OblfPJde2wJZifuQZnPmcVzEbFfq3HPR8TbChpZO/KthmpqCs789RM8uXAlf/7iEdtUv71w+XouvXM2D79czeRdBvLNk/Zl2u7DtifsHqW2oZHXVmzY0p6QtiXMr15PTW3D5nL9y8s2VxdNrBrAxJEDmFBVye7DK4vSLXlTfSOPzV/BvXPe4r45y1heU0ufMvGOCcOZPmUUR00ZxejB7shQKpqagpUb6lr0AqpOG4aX5TQMV6+rZVP91l1DK3r32vyLP7cXUHM7QHObwLDK8oK3B/QE3VENdS2wmuRBRgF8ARgaEZ/o7Eq3R77J4rp/vMpld83h8pP25WPvGL/N64kI7pn9Ft+4aw5L1mzilIPGcsGMyYwYUNGZsEvSqvV1aRJomRReX7mhRZe90YP7MrEqSQRJchjAxJGV7DKob8nUx7bW1BQ8vWh1kjhmL2XB8uS2ZvuPHcwHpoxi+tRdmDRyQMnG35PVNzZt6QWUc+CvXteym+jymtoWHRmaDezbe0tj8KBWbQFp43DVwL4M6tvbn9826I5kUQlcDByVjroX+FZ6p9hul0+yeGXpOo674u8cPnE4137i4O3aoTbUNfCTB+bxy0cX0L+8jK8dM5nTD9mtx/xSaWhsYvGqjVuSwrL1m1+v2rClC195717sMbySiSNzEkLVAPaoqtyqv3ZPNG9ZDffOeYt7Zy/lmUWrAdh9eP/NDeQH7ja0x3ymxbKhrqHF9QGtrxRuThAr19dtNa8EwyvLGTGg1fUALRqIkzOBfuU7b1VPIe10vaEqRk+K2jdfaXd6XUMTJ135D95au4m/nndEl/XLn7dsHRffPpvHFqxgv7GDufzEfdl/3JAuWXZXWLepngXVWxJB8+uFyzdQ17jlFH54ZfnmM4PcpDBmaL+d5mC5dO2mzQ3k/5y/nPrGYHhlOUftM4rpU0fxzj1H7DRX90cEazc25FwMlts7qOW1ArlVkM1699Lmap/c6wFadxEdPqA86cBgRdMdZxb3AR+KiNXp8FDg5og4urMr3R4VoydFzeK57e543/3rS1z10Hyu+dhBTJ+6S5euOyK489klfPPPL7K8ppbTDtmN/zx6b4b0L+/S9bSnqSl4c+2m5GK13K6o1TUsXbulG2pZLzF+WH8mtEgKlUwYMYChld0Ta0+xblM9D82t5t45S3nopWWsq22gX58y3rNXFdOnjuL9k0d22+fblRqbghXra1teD9BOF9G2bhXRr09ZiwN/Vc6Bf3NV0IAKhvYvd5flHqI7ksXTEfH2rHHdpWL0pFj88gtUDdy67eDZRas56Wf/4MMHjeO7p+zXxtxdY92men503ytc/9hCBvfrwwXHTOaUg8Z22ZdmU30jry5fv1W10YLq9Wys33Lhz8CK3psblZvPEPYcWcluwyp9cWEn1DU08fiC5gbypSxdW0tZL3HI7sOYPnUUH5gyirFD+xc1xtqGxpz2gKQNYMtVwluSwor1dTS20R4wuF+fFr/4m9sDWl8fMKDC7QE7mu5IFrOAkyPi9XR4PHBbRBzY2ZVuj4rRk+KFZ55i0qiBW0374u+e5sG5y3jswiO7pZ59zpK1XHLHC8x8bRUHjR/K5Sfuy5RdB2XPSHKWsrymrs22hDdWbyT3YxkzpF96oVrLBuaqARX+QhdIU1Pw/BtrNrdzvLKsBoCpuw5i+pSknWOf0QO7bPvX1Da03ROo1dnA6g0d3yoitydQi9tFpNN2luo121p3JItjSJ5//XA66t3A2RFxT2dXuj0qRk+KR/7xOIdOGN5ifPW6Wg7/zgOcceh4LjtharfF09QU3PLUYr79l5dYvaGOMw/fnS99YC8G9e0DJD1DXluxYau2hPnLali7aUsdcN8+vZgwYsBWSWGPEZVu8CsBry5fz31p4pj1+ioiYOzQfpsTx8G7D93qepyIYPWG+jbaAtJqoJxG4o5uFbE5AQzaui1gZNo11NcCWZZuaeCWNAJ4ByDgsYhY3tkVbq+K0ZPijvse4Zh9R7cYf+WD8/j+PXO5/8vvYc+R3X/PoNUb6vj+PXO56YnXGTGggv3HDmHB8hpeb3Wfo5EDK7ZqYJ5QVcmug/u57reHqF5XywMvLuW+OUt5dN5y6hqaGNK/D+/ccwR1DU1pF9HkL7dzQbPK8rKtrgdo0TCctgcM6V86t4qwnq+7ksVQkrvNbu5aFBGPdHal26Ni9KS47o4HOO2QLTe+bWwK3v29B9ltWH9+d/Y7ihHWZs8uWs237n6Rlevr2DP32oS0baH5jMN2DOtrG3jk5Wrum7OUf726Mrln0KAtbQBbEsKWLqKVO0BXZOt5uuNJeZ8CzgXGAs+QnGE8Rstncner1v24H355GW+s3shFx+5TpIi22H/cEP7wmcOKHYZ1k8qK3sx422hmvG10dmGzHiyfis5zgYOB1yLifcDbgeqCRpVhVatk8dvHX6dqYAXTp44qUkRmZju2fJLFpojYBCCpIiJeAvYubFgdW7lhS7JYtHIDD85dxqkHj/NFP2ZmBZJP5eni9LbktwP3SVpFkZ/BnVsNddMTryNo0YZhZmZdK59ncJ+cvrxM0oPAYOCvBY2qA2W9xIqaJFnUNjTyhycXceQ+o9jVj0U1MyuYbeqWEREPZ5cqrEF9+7AifcLaX194ixXr6/hoJ+4qa2Zm+etxlfy9e4nl6+uICG58/HXGD+/PEX4MqplZQfW8ZFEm6hqamPXaKp5YuJLTD9nNF7OZmRVYZrKQNKONcecUJpxsvdPE8OMHXqG8dy8+NG1csUIxM9tp5HNmcbGkzRfgSTofOLFwIXWsrFcS8qOvLOe4t41mmG+5bWZWcPk0cJ8A/EnS14BjgMnpuKLoXSaan9xwhhu2zcy6RT5dZ5dLOgG4H5gFnBJFfLxeczXUPqMHceBupfOkOjOzHVm7yULSOiBI7jQbQDkwAThFUkREfg9u6GK9y3oxdddBfO59e/qOnGZm3aTdZBERWz9dqAQI+PMXjyh2GGZmO5V8ekOdLGlwzvAQSScVNiwzMysl+fSGujQi1jQPRMRq4NLChWRmZqUmn2TRVhk/vcXMbCeST7KYKemHkiZKmiDpRyS9ojJJOkbSXEnzJF3QQbmDJTVKOiXfwM3MrPvkkyy+ANQBvwf+CGwCPpc1k6Qy4EpgBjAFOE3SlHbKfRe4J/+wzcysO+VzncV6oN2zgg4cAsyLiAUAkm4mufJ7TqtyXwBuIXkan5mZlaB8nsFdBfwnMBXo2zw+IrKewT0GWJQzvBg4tNWyxwAnkzzPu91kIels4GyA3XbzQ47MzLpbPtVQNwIvAXsA/w0sBJ7MY762rphrfeX3/wDnR0RjRwuKiGsiYlpETKuqqspj1WZm1pXy6dU0PCJ+Jenc9OFHD0vK5yFIi4HcW8KOZevHsU4Dbk6vxB4BHCupISJuz2P5ZmbWTfJJFvXp/zclfZDkgD82j/meBCZJ2gN4AzgVOD23QETs0fxa0nXAn5wozMxKTz7J4pvpFdxfAa4ABgHnZc0UEQ2SPk/Sy6kMuDYiZjc/CyMiru582GZm1p3ySRar0iu41wDvA5D0znwWHhF3A3e3GtdmkoiIT+SzTDMz6375NHBfkec4MzPbQXV0i/LDgMOBKklfzpk0iKRayczMdhIdVUOVAwPSMrm3K18L+LYcZmY7kY6eZ9HcTXZjRHwvd5qkDwGvFDo4MzMrDfm0WZzaxrgLuzoQMzMrXR21WcwAjgXGSPpJzqRBQEOhAzMzs9LRUZvFEmAmcAItb0m+jjyuszAzsx1HR20WzwLPSropIpqv4kbSu4Bvkcdtys3MbMeQzy3K6yUdQHKrjg8DrwK3FjowMzMrHR21WexF0rh9GrCC5OFHioj3dVNsZmZWIjo6s3gJeBQ4PiLmAUj6UrdEZWZmJaWjrrP/DrwFPCjpF5KOpO1nVJiZ2Q6u3WQRETw+QBMAAAlkSURBVLdFxEeAycBDwJeAUZKukjS9m+IzM7MSkHlRXkSsj4gbI+I4kudYPEPnnsltZmY9VD5XcG8WESsj4ud5PH/bzMx2INuULMzMbOfkZGFmZpmcLMzMLJOThZmZZXKyMDOzTE4WZmaWycnCzMwyOVmYmVkmJwszM8vkZGFmZpmcLMzMLJOThZmZZXKyMDOzTE4WZmaWycnCzMwyOVmYmVkmJwszM8vkZGFmZpmcLMzMLFNBk4WkYyTNlTRP0gVtTD9D0nPp3z8l7V/IeMzMrHMKliwklQFXAjOAKcBpkqa0KvYq8J6I2A+4HLimUPGYmVnnFfLM4hBgXkQsiIg64GbgxNwCEfHPiFiVDj4OjC1gPGZm1kmFTBZjgEU5w4vTce35D+AvbU2QdLakmZJmVldXd2GIZmaWj0ImC7UxLtosKL2PJFmc39b0iLgmIqZFxLSqqqouDNHMzPLRu4DLXgyMyxkeCyxpXUjSfsAvgRkRsaKA8ZiZWScV8sziSWCSpD0klQOnAnfmFpC0G3Ar8LGIeLmAsZiZ2XYo2JlFRDRI+jxwD1AGXBsRsyWdk06/GrgEGA78TBJAQ0RMK1RMZmbWOYposxmhZE2bNi1mzpxZ7DDMzHoUSbO258e4r+A2M7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDIVNFlIOkbSXEnzJF3QxnRJ+kk6/TlJBxYyHjMz65yCJQtJZcCVwAxgCnCapCmtis0AJqV/ZwNXFSoeMzPrvEKeWRwCzIuIBRFRB9wMnNiqzInADZF4HBgiaXQBYzIzs07oXcBljwEW5QwvBg7No8wY4M3cQpLOJjnzAKiRNLdrQy2IEcDyYgeRB8fZtXpCnD0hRnCcXW3v7Zm5kMlCbYyLTpQhIq4BrumKoLqLpJkRMa3YcWRxnF2rJ8TZE2IEx9nVJM3cnvkLWQ21GBiXMzwWWNKJMmZmVmSFTBZPApMk7SGpHDgVuLNVmTuBj6e9ot4BrImIN1svyMzMiqtg1VAR0SDp88A9QBlwbUTMlnROOv1q4G7gWGAesAE4q1DxFEFPqTZznF2rJ8TZE2IEx9nVtitORWzVRGBmZtaCr+A2M7NMThZmZpbJyWI7SRon6UFJL0qaLencdPxlkt6Q9Ez6d2wJxLpQ0vNpPDPTccMk3SfplfT/0CLHuHfONntG0lpJ55XC9pR0raRlkl7IGdfu9pN0YXorm7mSji5ynN+X9FJ6W53bJA1Jx+8uaWPOdr26yHG2+zmX2Pb8fU6MCyU9k44vyvbs4DjUdftnRPhvO/6A0cCB6euBwMsktze5DPhqseNrFetCYESrcd8DLkhfXwB8t9hx5sRWBrwFjC+F7Qm8GzgQeCFr+6X7wLNABbAHMB8oK2Kc04He6evv5sS5e265EtiebX7OpbY9W03/AXBJMbdnB8ehLts/fWaxnSLizYh4Kn29DniR5Cr0nuJE4Pr09fXASUWMpbUjgfkR8VqxAwGIiEeAla1Gt7f9TgRujojaiHiVpMffIcWKMyLujYiGdPBxkmuaiqqd7dmektqezSQJ+DDwu+6IpT0dHIe6bP90suhCknYH3g78Kx31+fS0/9piV++kArhX0qz0FioAoyK9tiX9P7Jo0W3tVFp+CUtte0L726+9W9mUgk8Cf8kZ3kPS05IelnREsYLK0dbnXKrb8whgaUS8kjOuqNuz1XGoy/ZPJ4suImkAcAtwXkSsJbmD7kTgAJJ7Xf2giOE1e2dEHEhyt9/PSXp3sQNqj5ILOU8A/piOKsXt2ZG8bmXT3SR9HWgAbkxHvQnsFhFvB74M3CRpULHio/3PuSS3J3AaLX/QFHV7tnEcardoG+M63J5OFl1AUh+SD+jGiLgVICKWRkRjRDQBv6CbTpk7EhFL0v/LgNtIYlqq9E6/6f9lxYuwhRnAUxGxFEpze6ba234ldysbSWcCxwFnRFpxnVZDrEhfzyKpu96rWDF28DmX4vbsDfwb8PvmccXcnm0dh+jC/dPJYjuldZa/Al6MiB/mjM+91frJwAut5+1OkiolDWx+TdLg+QLJLVfOTIudCdxRnAi30uIXW6ltzxztbb87gVMlVUjag+SZLU8UIT4geRAZcD5wQkRsyBlfpeTZM0iaQBLnguJE2eHnXFLbM3UU8FJELG4eUazt2d5xiK7cP7u71X5H+wPeRXL69hzwTPp3LPAb4Pl0/J3A6CLHOYGk98OzwGzg6+n44cADwCvp/2ElsE37AyuAwTnjir49SZLXm0A9yS+z/+ho+wFfJ/llOReYUeQ455HUUTfvo1enZf893R+eBZ4Cji9ynO1+zqW0PdPx1wHntCpblO3ZwXGoy/ZP3+7DzMwyuRrKzMwyOVmYmVkmJwszM8vkZGFmZpmcLMzMLJOThVkXkXSCpAvS15dJ+mqxYzLrKgV7rKrZziYi7mTr58yb7RB8ZmGWkvRRSU+kzyH4uaQySTWSfiDpKUkPSKpKy35R0pz0hnc3p+M+IemnbSz3AEmPa8uzJIam4x+S9N10nS+XyE38zNrkZGEGSNoH+AjJzRYPABqBM4BKkntUHQg8DFyaznIB8PaI2A84J2PxNwDnp2Wfz1kGJM+YOAQ4r9V4s5LiaiizxJHAQcCTyW126Edy07Umttwo7rdA8w3angNulHQ7cHt7C5U0GBgSEQ+no65ny510yVneLJIH55iVJJ9ZmCUEXB8RB6R/e0fEZW2Ua74/zgeBK0kSzKz0DqSdUZv+b8Q/3qyEOVmYJR4ATpE0EjY/u3g8yXfklLTM6cDfJfUCxkXEg8B/AkOAAW0tNCLWAKty2iM+RlKdZdaj+JeMGRARcyT9F8mTBHuR3GH0c8B6YKqkWcAaknaNMuC3aRWTgB9FxOq0+qotZwJXS+pPcrvqswr7bsy6nu86a9YBSTUR0eZZg9nOxNVQZmaWyWcWZmaWyWcWZmaWycnCzMwyOVmYmVkmJwszM8vkZGFmZpn+P+uTLnXaZD2mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(epsilons, attack_accuracy)\n",
    "plt.plot(epsilons, np.ones_like(epsilons) * acc, dashes=[2,2], label=\"base model\")\n",
    "plt.title(\"Differentially private logistic regression\")\n",
    "plt.xlabel(\"epsilon\")\n",
    "plt.ylabel(\"Attack accuracy\")\n",
    "plt.ylim(0, 1)\n",
    "plt.xlim(0.1, 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like epsilon=25 is a good choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP model accuracy with eps=100:  0.795\n",
      "DP model attack accuracy with eps=100:  0.52\n"
     ]
    }
   ],
   "source": [
    "dp_model = dp.LogisticRegression(epsilon=25, data_norm=5)\n",
    "dp_model.fit(x_train, y_train)\n",
    "\n",
    "dp_art_classifier = ScikitlearnLogisticRegression(dp_model)\n",
    "print('DP model accuracy with eps=100: ', dp_model.score(x_test, y_test))\n",
    "\n",
    "dp_attack = MembershipInferenceBlackBox(dp_art_classifier, attack_model_type='rf')\n",
    "dp_attack.fit(x_train[:attack_train_size].astype(np.float32), y_train[:attack_train_size].astype(np.float32),\n",
    "              x_test[:attack_test_size].astype(np.float32), y_test[:attack_test_size].astype(np.float32))\n",
    "dp_inferred_train = dp_attack.infer(x_train.astype(np.float32)[attack_train_size:], y_train.astype(np.float32)[attack_train_size:])\n",
    "dp_inferred_test = dp_attack.infer(x_test.astype(np.float32)[attack_test_size:], y_test.astype(np.float32)[attack_test_size:])\n",
    "dp_train_acc = np.sum(dp_inferred_train) / len(dp_inferred_train)\n",
    "dp_test_acc = 1 - (np.sum(dp_inferred_test) / len(dp_inferred_test))\n",
    "dp_acc = (dp_train_acc * len(dp_inferred_train) + dp_test_acc * len(dp_inferred_test)) / (len(dp_inferred_train) + len(dp_inferred_test))\n",
    "    \n",
    "print('DP model attack accuracy with eps=100: ', dp_acc)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}