analysis/covariance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentially Private Covariance\n",
    "\n",
    "SmartNoise offers three different functionalities within its `covariance` function:\n",
    "\n",
    "1. Covariance between two vectors\n",
    "2. Covariance matrix of a matrix\n",
    "3. Cross-covariance matrix of a pair of matrices, where element $(i,j)$ of the returned matrix is the covariance of column $i$ of the left matrix and column $j$ of the right matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load libraries\n",
    "import os\n",
    "import opendp.smartnoise.core as sn\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# establish data information\n",
    "data_path = os.path.join('.', 'data', 'PUMS_california_demographics_1000', 'data.csv')\n",
    "var_names = [\"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"]\n",
    "\n",
    "data = np.genfromtxt(data_path, delimiter=',', names=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functionality\n",
    "\n",
    "Below we show the relationship between the three methods by calculating the same covariance in each. We use a much larger $\\epsilon$ than would ever be used in practice to show that the methods are consistent with one another.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scalar covariance:\n",
      "[[94601.1327106]]\n",
      "\n",
      "covariance matrix:\n",
      "[[3.14917795e+02 9.46406149e+04]\n",
      " [9.46406149e+04 2.65091456e+09]]\n",
      "\n",
      "cross-covariance matrix:\n",
      "[[3.14928928e+02 9.50137258e+04]\n",
      " [9.46392159e+04 2.65184060e+09]]\n"
     ]
    }
   ],
   "source": [
    "with sn.Analysis() as analysis:\n",
    "    wn_data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "    \n",
    "    # get scalar covariance\n",
    "    age_income_cov_scalar = sn.dp_covariance(left = sn.to_float(wn_data['age']),\n",
    "                                        right = sn.to_float(wn_data['income']),\n",
    "                                        privacy_usage = {'epsilon': 5000},\n",
    "                                        left_lower = 0.,\n",
    "                                        left_upper = 100.,\n",
    "                                        left_rows = 1000,\n",
    "                                        right_lower = 0.,\n",
    "                                        right_upper = 500_000.,\n",
    "                                        right_rows = 1000)\n",
    "    \n",
    "    # get full covariance matrix\n",
    "    age_income_cov_matrix = sn.dp_covariance(data = sn.to_float(wn_data['age', 'income']),\n",
    "                                           privacy_usage = {'epsilon': 5000},\n",
    "                                           data_lower = [0., 0.],\n",
    "                                           data_upper = [100., 500_000],\n",
    "                                           data_rows = 1000)\n",
    "\n",
    "    # get cross-covariance matrix\n",
    "    cross_covar = sn.dp_covariance(left = sn.to_float(wn_data['age', 'income']),\n",
    "                                   right = sn.to_float(wn_data['age', 'income']),\n",
    "                                   privacy_usage = {'epsilon': 5000},\n",
    "                                   left_lower = [0., 0.],\n",
    "                                   left_upper = [100., 500_000.],\n",
    "                                   left_rows = 1_000,\n",
    "                                   right_lower = [0., 0.],\n",
    "                                   right_upper = [100., 500_000.],\n",
    "                                   right_rows = 1000)\n",
    "\n",
    "# analysis.release()\n",
    "print('scalar covariance:\\n{0}\\n'.format(age_income_cov_scalar.value))\n",
    "print('covariance matrix:\\n{0}\\n'.format(age_income_cov_matrix.value))    \n",
    "print('cross-covariance matrix:\\n{0}'.format(cross_covar.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DP Covariance in Practice\n",
    "  \n",
    "We now move to an example with a much smaller $\\epsilon$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Non-DP Correlation Matrix:\n",
      "          0         1         2         3         4\n",
      "0  1.000000  0.055088  0.025901  0.103524  0.196072\n",
      "1  0.055088  1.000000 -0.023111 -0.211106 -0.073125\n",
      "2  0.025901 -0.023111  1.000000  0.305238  0.053860\n",
      "3  0.103524 -0.211106  0.305238  1.000000  0.153143\n",
      "4  0.196072 -0.073125  0.053860  0.153143  1.000000\n",
      "\n",
      "\n",
      "DP Correlation Matrix:\n",
      "          0         1         2         3         4\n",
      "0  1.000000  0.575078 -0.369702  1.090716  0.173756\n",
      "1  0.575078  1.000000 -0.068965  0.844628  0.199602\n",
      "2 -0.369702 -0.068965  1.000000 -1.054197  0.167272\n",
      "3  1.090716  0.844628 -1.054197  1.000000 -1.385910\n",
      "4  0.173756  0.199602  0.167272 -1.385910  1.000000\n"
     ]
    }
   ],
   "source": [
    "with sn.Analysis() as analysis:\n",
    "    wn_data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "    # get full covariance matrix\n",
    "    cov = sn.dp_covariance(data = sn.to_float(wn_data['age', 'sex', 'educ', 'income', 'married']),\n",
    "                                          privacy_usage = {'epsilon': 1.},\n",
    "                                          data_lower = [0., 0., 1., 0., 0.],\n",
    "                                          data_upper = [100., 1., 16., 500_000., 1.],\n",
    "                                          data_rows = 1000)\n",
    "analysis.release()\n",
    "\n",
    "# store DP covariance and correlation matrix\n",
    "dp_cov = cov.value\n",
    "dp_corr = dp_cov / np.outer(np.sqrt(np.diag(dp_cov)), np.sqrt(np.diag(dp_cov)))\n",
    "\n",
    "# get non-DP covariance/correlation matrices\n",
    "age = list(data[:]['age'])\n",
    "sex = list(data[:]['sex'])\n",
    "educ = list(data[:]['educ'])\n",
    "income = list(data[:]['income'])\n",
    "married = list(data[:]['married'])\n",
    "non_dp_cov = np.cov([age, sex, educ, income, married])\n",
    "non_dp_corr = non_dp_cov / np.outer(np.sqrt(np.diag(non_dp_cov)), np.sqrt(np.diag(non_dp_cov)))\n",
    "\n",
    "print('Non-DP Correlation Matrix:\\n{0}\\n\\n'.format(pd.DataFrame(non_dp_corr)))\n",
    "print('DP Correlation Matrix:\\n{0}'.format(pd.DataFrame(dp_corr)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0, 0.5, 'age'),\n",
       " Text(0, 1.5, 'sex'),\n",
       " Text(0, 2.5, 'educ'),\n",
       " Text(0, 3.5, 'income'),\n",
       " Text(0, 4.5, 'married')]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAFACAYAAAC1Gdb6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XvcbnOd//HXeztW5DRylgoz6SShaMQv20RyapwTGtoqRkVFkWQopbOamXYlZCoq1S47BpOiybR3ipJyymE75yxn+/3747tuXa593afrfK/7/fRYj31d61rXWp+93etzf9f3KNtERERE1NGMQQcQERER0Ssp6ERERERtpaATERERtZWCTkRERNRWCjoRERFRWynoRERERG2loBMREyJpW0l/knStpCNbfP4OSb+T9FtJl0jaYBBxRkQ0UubRiYjxSFoMuBrYBlgAzAP2sv2HhmOea/uB6vWOwLtsbzuIeCMiRqRGJyImYlPgWtvX234c+DawU+MBI4WcynOAPEVFxMAtPugAIqI7rvnHN7RdsFj/F/99EDCrYdds27Mb3q8B3NzwfgHw6ubzSDoYOAxYEnh9u/FExGDcM++ytvPIiptspG7G0i0p6EQEVaFm9rgHjn+eLwFfkrQ3cDSwX6fnjIjoRJquImIibgHWani/ZrVvNN8Gdu5pRBERE5AanYi6UE+fW+YB60l6AaWAsyew9zMuL61n+5rq7fbANUTE1DJjKFufOpKCTkRdqHcJyvaTkg4BzgMWA06xfaWk44D5tucAh0iaCTwB3EuarSKmHPUwjwxKCjoRddHjJzHbc4G5TfuOaXj97p4GEBHRhvTRiYiIiNpKjU5ETai3fXQiYjqYUb88koJORF3UsBNhRESnUtCJiIiIIp2RI2Jo1TBBRUSf1TCPpKATURc1bFuPiOhUCjoRdVHDJ7GI6LMaPjDV728UERERUUmNTkRN1HFG04jorzrmkdToRNTFjBntbxERfSBpW0l/knStpCNbfH6YpD9IukLShZKe3+k1k+Ei6kJqf4uIgDIfV7vbOCQtBnwJ2A7YANhL0gZNh/0G2Nj2y4HvAp/s+K/U6QkiIiIiJmBT4Frb19t+HPg2sFPjAbZ/avvh6u2lwJqdXjR9dCLqIjMjR0SnOlhKRtIsYFbDrtm2Zze8XwO4ueH9AuDVY5zyAOAnbQdUSUEnIiIiOlYVamaPe+AESNoH2BjYstNzpaATURdZ1DMiOtXbmuFbgLUa3q9Z7XsGSTOBo4AtbT/W6UVT0ImoCaXpKiKG2zxgPUkvoBRw9gT2bjxA0iuBLwPb2r6zGxdNQSeiLjJ6KiI61cM8YvtJSYcA5wGLAafYvlLSccB823OAk4BlgO9Uc/rcZHvHTq6bgk5EXaTpKiKGnO25wNymfcc0vJ7Z7WumoBNRF2m6iohO1fCBKQWdiIiIAOrZ169+RbeIiIiISmp0IuoinZEjIhaRgk5ETSiLc0ZEp2r4wJSCTkRd1DBBRUR0KgWdiLpIjU5EdKqGeSQFnYi6SI1ORHSqhnkkBZ2IuqhhgoqI6FT96qgiIiIiKqnRiYiIiKKGEwamoBNRExleHhGdUpaAiIihlT46ERGLSEEnoi5qWOUcEX1WwwemFHQi6qKGVc4REZ1KQSeiLlKjExGdqmEeySNgRERE1FZqdCLqooZt6xHRZzVsAq/f3ygiIiKikoJOdETSsZLO6OD7V0raqosh9Z2kLST9afBxzGh7i+gGSadKOr6D7z8k6YXdjKnfJL1F0n8POo62zVD725Ca1hlO0g2S7pT0nIZ9B0q6qEfXu0jSo5IelPSApF9LOlLSUg3HHCvpieqGv0/S/0rabIxzri/pO5L+Iul+SVdIOkzSYr34O3SiVRK0/RLbF/XgWhdJsqRXNO3/frV/qwmex5LWHesY2xfb/vsOwu2OGiaoOqvyzyNVPhi519+hhpJndc88XuWDeySdL+kfxjjnppLmVue7R9KvJL2tP3+jyanu0QMb99lexvb1PbjWDdW/49817f9NdY+vM4FzrFMdO2aXD9v/ZfufOos4umlaF3QqiwHv7uP1DrG9LLAacDiwJzBXekYHizNtLwOsDFwCnN30OQCSXgT8H3Az8DLbywG7ARsDy04mqFY373g39BRwNbDvyBtJKwGbAXd16wJD9W8ktb/FoOxQ5YPnAycCRwBfazrmk1U+WBO4Ezi11YmqB6L/AX4GrAusBLwT2G6yQdU0H/wZ2GvkjaSXAc/u5gVq8G+EpLa3YZWCDpwEvE/S8q0+lLS5pHlVbck8SZs3fHaRpH+T9Ivqqey/m58YRmP7r1VNxo6UX77btzjmCeA0YFVK0mr2UeB/bR9m+7bqO3+yvbft+6oYd6yah+6r4n1xQ/w3SDpC0hXAXyUtPsq+1SV9T9Jdkv4s6dDR/l5V7dLt1b/XzyW9pNo/C3gL8IHq6fRHDTHMrF4vJelzkm6tts+N1HZJ2krSAkmHq9TC3TaBJ9X/AvZoqN3aC/g+8HhDvJtK+mX173ObpC9KWrL67OfVYZdXMe/REMcRkm4Hvj6yr/rOi6on6Y2q96tX/25bjRNr52bMaH+LgbJ9v+05wB7AfpJe2uKYh4FvAot8VjkJOM32J2z/xcWvbe8+coCkt0u6tvoZnSNp9YbPLOlgSdcA14yx7x9UapbukfQnSbvTgqQVJP24+vm/t3q9ZvXZCcAWwBere+uLDddbt3q9nKTTq+/fKOloVbVdkvaXdImkT1Xn/rOk8Qp036DhwQfYDzi9KebtVWp5HpB0s6RjGz4eyQf3VTFvVsXxC0mflXQ3cOxIbNX5NlepbV+rev+KKt5Ra+UGroZ5ZHgj65/5wEXA+5o/kLQicA7wBUpB4zPAOSo1AyP2Bt4GPA9YstV5xmL7piqGLVpcfylgf+Bm239p8fWZwHdHO7ek9YFvAe+h1A7NBX408ou8shelkLW87Seb9wELgR8BlwNrAFsD75H0hlEu+xNgPcq/x2WUwga2Z1evP1lVT+/Q4rtHAa8BNgReAWwKHN3w+arAclUcBwBfkrTCaH9/4FbgD8BINfK+NCU24CngvcDfUQqcWwPvqmJ+XXXMK6qYz2yIY0XKU/isxpPZvo7yVH6GpGcDX6f88rlojDgjALD9K2ABrfPBMpSHhd+0+OzZlJ/fsfLB64GPA7tTapRvBL7ddNjOwKuBDVrtU2nmP59S4HoepUb63yVtwKJmUH7+nw+sDTwCfLH6ex4FXEyp4V7G9iEtvn8y5X5/IbAl5f5tfLh5NfAnyr37SeBr0pjVCpcCz5X04urhZ0+guX/hX6vrLE/Jge+UtHP12Ug+WL6K+ZcNcVwPrAKc0Hgy2/8LfBk4TdKzqut92PYfx4gzuiwFneIY4F8lrdy0f3vgGtvfsP2k7W8BfwQaf0l/3fbVth8BzqL8kp6sWym/OEfsLuk+SpPUq4BdRvneSsBtY5x3D+Ac2+dXtUOfAp4FbN5wzBds31zF32rfJsDKto+z/XjVfv4VSpJYhO1TbD9o+zHgWOAVkpYbI8ZGbwGOs32n7bsoNVZvbfj8ierzJ2zPBR4Cxusbczqwb/UEtXxDchqJ99e2L63+/95ASUpbjnPOhcBHbD/W9O82cs6vANdSmhVXoxTgei9NV3XRnA/eV+WDa4FlKA8/zVag5POx8sFbgFNsX1bdnx8ENtMz+6d83PY9TT/XjfveBNxg++vVPfMb4HuUJvNnsH237e/Zftj2g5RCwHj3FgANBZEPVvnkBuDTPDMf3Gj7K7afotR8r0YpbIxlpFZnG+Aq4JammC+y/TvbC21fQXlQHC/mW22fXP17LJIPKHlwOeBX1fW+NM75osumfHtiN9j+vaQfA0dSfvhHrE556ml0I6VGYcTtDa8fpiQiJP0nsE+1/2O2PzZGCGsA/9vw/izb+4x2cIO7KTf3aJ4Rv+2Fkm7mmfHf3OJ7jfueD6xeJdoRi1Gexp6hSk4nUJLeypQCAZQnrvvHiLNlvNXr1Rve391Q6wQN/95jOJuSIO+mJLnmmNen1NRtTGmvXxz49TjnvMv2o+Mc8xVgDjCr+qXSc0qn4rpYA7in4f2nbB892sGVeyn322qUh7FWVqfUsgJg+6GquWUN4IZq90Tywaub8sHitL63ng18FtiWUhADWFbSYlXhZCx/ByzBovmgZe61/XBVmTNePvgGpQnqBSxau4ukV1P6Sr2UUkO/FPCdcc7Z6t/sabafkHQqpWXgMNse53yDVcMHn9To/M1HgLfzzBvpVsqN3Whtmp4CWrH9jqp6c5mxCjlV2+2raFFwmIALgH8e4/NnxF9V667FM+NvddM17rsZ+LPt5Ru2ZW2/scX39gZ2ojSpLQesM3LpMa41aryUf+tbx/nOmKp+DT+hdMpcJBkD/0H5xbCe7ecCH2qId9TTjvVh1cTwOUqn0mOrJtCIcUnahJKDLpnM96qf818yuXzwHEqt8GTzwc+a8sEytt/Z4nuHU2pcX13dWyNNPxPJB3+h1OA254Nxc+9YbN9I6ZT8RspDULNvUh5Q1nIZ3PGfE4h3vHywBuX3y9eBT6thlO10JGnbqm/XtZKObPH5UpLOrD7/v6Yax7akoFOxfS1wJtDY0XYusL6kvVU65e5Babv+cafXk/RsSVsCP6RUac5t4zQfATaXdJKkVavzrivpDJXO1WcB20vaWtISlMTzGM+sPRrPr4AHVTrfPkvSYpJeWiXkZstW57+bUjvSXMC7g9LePppvAUdLWlmlU/cxLNqG3o4PAVtW1d/NlgUeAB6qmreaE/Z4MbfyeWC+7QMpfbz+c5Lfb49mtL/FQEl6rqQ3UfrMnGH7d22c5gPA/pLeP9KPsOr8OtIP51vA2yRtWP2y/Rjwf6PcF6P5MSUnvlXSEtW2iRoGOTRYltIv576qsP+Rps9HvbeqGp+zgBMkLSvp+cBhdCcfHAC83vZfR4n5HtuPStqU8vA24i5KrdmE80H1cHkq5aHnAErT4r+1GXd/9LAJvKr1/xJlJOAGwF4t+ncdANxre11KjeAnOv0rJcM903HA03Pq2L6b0iZ9OOWX9weAN43SMXiivijpQcpN/jlK+/a2theO/bVFVR1fN6PUnFwp6f7qfPOBB23/idJ8djLlCWkHynDWx1ufseU1nqL8G2xIeRL6C/BVSo1Ns9Mp1cu3UDoBX9r0+dcoHRrvk/SDFt8/vor9CuB3lGr2ticfa/g73Gp7tCfk91GS2YOU5qYzmz4/ltKR8D6NMrqkkaSdKFX1IwWmw4CNJL2lndgnJX10pqIfVfngZkpfrs/wzA63E1Z1fH19tV0v6R5gNtVDlO0LgA9TcsRtwIsYpa/dGNd4kNK5f09KDdHtlF9ErWopPkfpE/gXSi44t+nzzwO7VqOQvtDi+/9K6Rx8PaWG65vAKZOJd5S/w3W254/y8buA46r/J8dQClsj33uY0jT/iyofvGYClzuU0mn7w1WT1dsohc1FOpsPjd6OutoUuNb29dXvoW9TWgEa7UTpcwWlc/3W43QyH5eGvbkwIibmpsM+1PbNvPZnPpbSTkTw4AMPtJ1Hln3uc8fMI5J2pTzYH1i9fyulafOQhmN+Xx0zMmXHddUxbVcwpEYnIiIiOiZplqT5Ddus8b/Vexl1FVEX6WsTER1a2EErkct8abPHOOQWyoCYEWuyaAfzkWMWqMw0vRyl60jbkhkjakIz1PY2ofOPP1riMEl/UFlv7cKqA2lExIh5wHqSXqAyce2elFFujeZQZq0G2BX4n06H5KegE1EXgx8t8RtgY9svp3Qi/GSX/4YR0WML3f42nmoetEOA8yhz1p1l+0pJx0nasTrsa8BKkq6lDOZY5KFqstJ0FRET8fRoCYBqyPJOlNF1ANj+acPxl/K3CTMjYopY2OMBStWs9nOb9h3T8PpRWsy03Ym6FHQydCzqauIN5h0sqld1GmzsODi7am8fsQbPnAF2AWWNn9EcQJmocSpJHom6mtajKutS0ImI3nYinEQY2oeypMaE1jWKiOFRxylnUtCJiImYyGgJJM2kTHy3Zb/W+IqIGEsKOhF10dsZjp8eLUEp4OzJM6fHR9IrKau/b2v7zl4GExG9UcMKnRR0IupCHfTRGY/tJyWNjJZYDDhlZLQEZV2vOcBJlNWjv1PN2H6T7R1HPWlERB+koBNRFz1es2oCoyVm9jSAiOi5Xo+6GoQUdCLqYoIT/0VETCcp6ERERASQUVcRMcyy1lVEdKiOBZ1kxoiIiKit1OhE1EX66ERELCIFnYi66PGoq4iov4kszjnVpKATURNKH52IiEWkoBNRF2m6iogO1bEzcgo6EXWRpquI6NBC6lfQSV13RERE1FYKOhEREVFbQ9t0JUmuY2NhRK/0cFHPqSg5JGLy6njLDE1mlLS1pH0lzQIYL0FJmiVpvqT5s2fP7k+QEcNMan+rgcnmkOo7ySMRNTcUNTqS3gCcBBwGnCPpObY/O9Z3bM8GRjJT/YqgEZOkaTzqqp0cAskjEc1qWKEz2IKOysQfSwBvBfYCVgN+C5w1yLgipqRpOI9OckhEdy2sYUln0DU6i9t+TNINwN7Aa4F9bd8iaR/gHttzBxphxFRRkyaoSUoOiYgxDewRUNI/AMdJWhq4CTgCeKftP0l6ZfX+gUHFFzHlzFD72xSUHBLRfbbb3obVQGp0JG0D7Au8CngI+BiwFnCapMuBjYCjbF8yiPgiYrglh0TERPW9oCNpY+CrlDb1+cALgWNtf1jS9yjt7U/Y/m2Gh0ZEs+SQiN5JH53uWBM43/bPJV0MvB44WtJHgJNt3zNyYBJUxCRMn87IySERMWF9y4ySXiRpfeCXwCskbePiQuDPwNrAG6pjp2angYgB0gy1vU0FySERvWe3vw2rvtToSNoBOB64Abgf+B6wi6RVgSuBDYBLgE2Bb+UpLKINNf7dnhwSUW+SVgTOBNah3Oe727636ZgNgf8Angs8BZxg+8zxzt3zgo6k1wDHANtU2+eBh4GfA4cAdwOzKPNf/IukZwGPJlFFTFJNl4BIDononwHeNkcCF9o+UdKR1fsjmo55mDJ9xDWSVgd+Lek82/eNdeJ+1OgsAN4FbAi8B3g1pUS2NmUW08uBLYDPAnvafqQPMUXUT00LOiSHRPTNADsj7wRsVb0+DbiIpoKO7asbXt8q6U5gZWDMgk7PM6PtBbbnAVsCZ9i+DjgD+DvgL8BjlAS2i+0reh1PREwtySERU0Pj2nHVNmsSX1/F9m3V69uBVca51qbAksB14524n6OufgccJGkJ4M3Ae21fCyDpU7YX9jGWiJh6kkMihljT2nGLkHQBsGqLj45qOo8ljVq1JGk14BvAfhO57/tZ0JkLLAXsSOlA9AsooyOSoCK6oMadkSvJIRE91ss+OrZnjvaZpDskrWb7tqogc+coxz0XOIcyIeilE7lu3wo6th+gzFr6X7afHJnIKx0GI7pjqgwTb1dySEStzQH2A06s/vxh8wGSlgS+D5xu+7sTPfEgei8+BZnIK6LrNKP9bWpJDonoEXewdehEYBtJ1wAzq/dI2ljSV6tjdgdeB+wv6bfVtuF4J+77zMhJThE9Uv+mKyA5JKKXBjXqyvbdwNYt9s8HDqxen0EZiDApA1nUMyJ6oOZNVxER7ahNQefey3836BCetsIrXjboECKiDffMu2zQITxtxU02GnQIEbVQm4JOxLQ39fraRMSQqWPLcDJjRERE1FZqdCJqou7DyyOi9wa4BETPpKATURfTZNRVRPRODcs5KehE1EZ9F/WMiGhbCjoRdZEanYjoUDojR0REREwhqdGJqIvU6EREh9IZOSKGltJHJyJiESnoREREBFDPPjop6ETURZquIiIWkYJORF1kwsCI6NDC+lXopKATURtZ6yoiOmTqV9JJZoyIiIjaSo1ORF2k6SoiYhFDVdCRtIzthwYdR8RUpHRGBpJHIjpRx1FXQ9N0JWlL4JDq9bhxSZolab6k+bNnz+55fBFDTzPa32oieSQimg1Tjc66wKYAtheOd7Dt2cBIZvK9l/+uh6FFhKRtgc8DiwFftX1i0+evAz4HvBzY0/Z3+x9lZ3nknnmX9TC0iOGXmZF7oHoCux+YD7yqxedyHevSIrqth310JC0GfAnYBlgAzJM0x/YfGg67CdgfeF/PAhk9vuSRiC6o410ysIJOVa38LOAw4DFgeWCmpLspTWqXAFfZvmFQMUZMKb3to7MpcK3t68ul9G1gJ+Dpgs7IvSpp3JqUbkkeiYjxDLJxfmXbf7W9k+3dgUOBXwHXUQpgW0ENB/RHDKHGvirVNqvpkDWAmxveL6j2DVrySESMaSA1OpIOBraXdA1wv+1jbP9R0sXAY7aPkLSE7ScGEV/ElNTBop5NfVWmhOSRiO6rYwtv32t0JO0J7AbMAlYCXli1/wM8CmwAkOQUMUlS+9v4bgHWani/ZrVvIJJHImKi+lqjI2kZ4HHgcOCNwPOAN9p+StK6wKmUER0RMUnq7YSB84D1JL2AUsDZE9i7lxccTfJIRO9k1FUHJL0LWIrytHUq8CvbM6vP3g6sBxyVJ7CINvVwPhzbT0o6BDiPUog4xfaVko4D5tueI2kT4PvACsAOkj5q+yXdjCN5JKK3BtV0JWlF4ExgHeAGYHfb945y7HMpAyF+YPuQ8c7dl4KOpIMow053sX2LpHWADSStDbwJOAjYO8kpYnjZngvMbdp3TMPreZQmrZ5IHomotSOBC22fKOnI6v0Roxz7b8DPJ3rinhd0JD0L2A74MPCYpHdQqp03BI4HlqYkpyt7HUtErdV4ravkkYj+WDi4lqudKKMkAU4DLqJFQUfSq4BVgHOBjSdy4p4XdGw/ImkucCJlSOpVwPXAN4FjgSdsP9nrOCJi6koeiRh+1bQUjVNTzK5GdE7EKrZvq17fTinMNJ9/BvBpYB9g5kTj6lcfndOB3wDX2b5H0t6UCciU5BTRJfVf1DN5JKLHOumjM940FZIuAFZt8dFRTeexpFaBvAuYa3vBZBYx7ktBx/ajlCnjZ0g6AHgPsJfth/tx/YhpoUaLc7aSPBIxtY0MHGhF0h2SVrN9m6TVgDtbHLYZsEU1KGEZYElJD9k+cqzr9nvCwKWBhZTe1Ff1+doRtdbj4eXDJHkkokcGOGHgHGA/SvP0fsAPmw+w/ZaR15L2BzYer5ADfZ4wsHryOjXJKaIHejth4NBIHomopROBbaqZzmdW75G0saSvdnLivi8BkRWEI3qkgyUgpprkkYjeWDigpeFs3w1s3WL/fODAFvtPpcylNa6BrV4eERERw6WOjxC1Keis8IqXDTqEiJjiVtxko0GHEBFdVpuCzv0LFgw6hKctt+aaLLjz7kGH8bQ1n7fSoEOIfphifW2G0cVXXjPoEJ62xUvW4/a7W86APxCrrrTCoEOIaEttCjoR0970GXUVET1Sx+5vKehE1IRqPo9OREQ7UtCJqIs0XUVEhxamRicihlaariKiQ2m6iojhlaariIhFJDNGREREbaVGJyIiIgBYWL+WqxR0ImojfXQiIhaRgk5ETSijriKiQ+mMHBHDaxot6hkRvZGCTkQMr9ToREQsIgWdiLpIQSciOlTHCQNT1x0RERG1lRqdiLpIH52I6FANK3RSoxMRERH1lRqdiJrI8PKI6NRC6lelM/CCjiS5Gs8maUnbjw86pogpaRpPGJg8EhGjGWjTVVNy2hc4WNKECl+SZkmaL2n+7NmzexpnxJSgGe1vU1jySET32G57G1YDrdFpSE7vAg4E9rT95AS/OxsYyUy+f8GC3gQZMVVM0xqdbuaRi6+8pjdBRkwRw1xgaddAHuVUdSZQsTywBbCX7aslLdHq2IiIRskjETERAynojDyBubiviuP1khaz/QSApC0lLes6Fi8jekFqf5uCkkciYiIG1jgv6S2Svlm9nQusCfxj9dkewAeAZw0ovIipZxr20UkeieiuhW5/G1Z9y3Atqo7nAK+VdJLt04BHgMMl/Rj4IHCk7Tv7FV9EDL/kkYiYrL50Rpa0RENV8vrAg7Zvk/RS4DdVVfNhklYE1gdusH17P2KLqAvVvDNy8khE79WxlbfnNTpVQjpB0mqSngd8EthG0iq2HwReCewv6Wu277F9aZJTRBtq3EcneSSiPwY1vFzSipLOl3RN9ecKoxy3tqT/lnSVpD9IWme8c/e0oFNVMz8PWBI4GLgXOAWYSek0uFqVpE4Gtpa0SkZHRLRpxoz2tyGWPBIxLRwJXGh7PeDC6n0rpwMn2X4xsCkwbtN0z5quGibxukTSi4FtKG3mHweeBPYE1pa0NKWa+TW27+hVPBEx9SSPREwbOwFbVa9PAy4Cjmg8QNIGwOK2zwew/dBETtyzR7mGSbzeC+wDPAhsBhxLKa19BTCwCfCxVDNHdKiGTVfJIxH9tdBue2ucabzaZk3i0qvYvq16fTuwSotj1gfuk3S2pN9IOknSYuOduOs1OiNPYFXV8XLAdsDOtu+VtCWwG2XI5xdsXyzpsyMdDCOiAzXqjJw8EjH1NM00vghJFwCrtvjoqKbzWFKrTj+LUyYGfSVwE3AmsD/wtbHi6mqNTuOaM8DqwKOUv9RWALZ/BiwAdgfeXZXEJjRVe0SMTZrR9jZMkkciBsdufxv/3J5p+6Utth8Cd0haDaD6s1XfmwXAb21fXy3z8gNgo/Gu29UanYZq5kMoT2D/DdwK7CzpEdvnAjcDFwNftv1UN68fEVNf8kjE4Cwc3PDyOcB+wInVnz9sccw8YHlJK9u+C3g9MH+8E/ei6WpnypPWTsD3KBN4/S/wGUl7U6qd3pgOgxFdNsR9bSYreSRi2jkROEvSAcCNlPsfSRsD77B9oO2nJL0PuLBq1v41pZ/emHox6mo54HOUBPWE7e2rYK8Ang0cbfumHlw3YnqrUR8dkkciBsIMpkbH9t3A1i32zwcObHh/PvDyyZy7FwWdGyhzXNxqewsASYcBKwEftf14D64ZEfVyA8kjEdEFvSjo/JrStrZQ0lbA2sBbgP16mZyWW3PNXp26LWs+b6VBhxDTzZB1Ku7QQPLIFi9Zr1enbsuqK7WcHDaiZ+q4BETXCzq2H5J0ErAj8H7gbmB/278eYSIIAAAf00lEQVTv9rUiokGNmq6SRyKiW3oyM3I16c+XJZ1Sve/5/BZ3/egnvb7EhK28w3YcftoPBh3G0z69384AXPOPbxhwJH+z3iXnDTqE2qnbqgeDyCP3zLus15eYsBU32Wjo4gF451e+M+BI/uY/3r7boEOonYX1q9Dp7erlmcAroo+GfM2qdiWPREQnelrQiYiIiKkjfXQiYnjVrOkqIvqvjgWdetZ1R0TXSdpW0p8kXSvpyBafLyXpzOrz/5O0Tv+jjIh4phR0Iuqih6uXV+tJfYmyJMMGwF6SNmg67ADgXtvrAp8FPtHlv2FExKSloBNRFzNmtL+Nb1Pg2moxvceBb1NmLW60E3Ba9fq7wNaq21CwiJpbaLe9DasUdCJqYqHU9iZplqT5DdusptOvQVlIc8SCal/LY6qVhe+nzGQcETEw6YwcUROdzH9hezYwu2vBRMSUNMw1M+1KQSciJuIWYK2G92tW+1ods0DS4pSFOe/uT3gR0Q0ZdRUR09U8YD1JL5C0JLAnMKfpmDnAftXrXYH/cR2zZkRMKanRiaiJXlY5235S0iHAecBiwCm2r5R0HDDf9hzga8A3JF0L3EMpDEVEDFQKOhE10evKE9tzgblN+45peP0okMWHIqawrHUVEUMrjUQREYtKQSeiJuo4WiIi+quO3epS0ImoiTomqIjorzrmkRR0ImqijgkqIqJTAx9e3jxF/ESnjG+cyXX27MxzFjGdJY9ExGgGXqMzMs+GpIMpU8ivJOk4282TkTV/r3EmV9/1o5/0NtCIGFrdyiP3zLust4FGDLk69vUbeI0OgKR3AbsAXwc2Bw4fbEQRU89Ct7/VQfJIRLQy8BqdyuqUmVQPAG4EPiBpaUC2HxloZBFTRProJI9EdKqOaaTvBR1JM2wvbHi/GKWq+RxKcnpzNQvrQcCjkk7JNPIR41vI9LlNkkcieiNNV10wkpwk7SLpdcCKwPHAC4CLbT8uaT/gYOBnSU4RE2O77W2qSR6JiInqW42OJDV0GNwLOAH4BfAwcDqwHfBNSZsDLwR2s31tv+KLmOqmw6/y5JGI3vKAaoYlrQicCawD3ADsbvveFsd9EtieUlFzPvDu8R5k+lKj05ScVgXWBrYC3gVcBvwLpdC1IfAOYEfbV/YjtoiYGpJHImrtSOBC2+sBF1bvn6F6gHkt8HLgpcAmwJbjnbjnNTpNyelwSpCvAubZ/h9JI4sEHkFZEflc4K+9jiuiburYtj4ieSSiPwbYyrsT5cEF4DTgIsr93MjA0sCSgIAlgDvGO3HPa3QaktMOwEzg3cCPgOMlrW77ZuBc4CfA5b2OJyKmnuSRiOHXOAFntc2axNdXsX1b9fp2YJXmA2z/EvgpcFu1nWf7qvFO3Jc+OpL+HpgF3FElpEMkfQH4jqQ9bd8o6Ru2n+xHPBF1VPf+tskjEb3XybxaTRNwLkLSBcCqLT46quk8lrRIJJLWBV4MrFntOl/SFrYvHiuuntTotJh+/VbgB8AakvYFsH0o8Efg1Gpo6FO9iCViulhot70No+SRiHqxPdP2S1tsPwTukLQaQPXnnS1OsQtwqe2HbD9EqcHdbLzrdr1Gp6ktfQ9gOeB+4CzgSeC1khbaPsP2AZJWsZ3kFNGhIS2vtCV5JGIwBlgzPAfYDzix+vOHLY65CXi7pI9T+uhsCXxuvBN3vUanITkdCBxd7T6WUuV8CWUo6A5V8oLWpbaImKQ6zaOTPBIxGAOsGT4R2EbSNZR+eCcCSNpY0lerY74LXAf8jtIX73LbPxrvxF2r0Rl5ApM0A1gG2AF4h+1fSDoX+AqlWvkLwBPAz+FvCS0iInkkYnqyfTewdYv984EDq9dPAQdN9txdqdFprGYGlrb9AHAtsJ6kZ9u+CfggsE01o+mZDb2rI6ILpnofneSRiOiFrhR0GqqZ3wN8qeoU+GfKmPiXVU9nLwIWSlrSDWvURER3TPWmq+SRiMGb6nmklW42XR0I7AnsX1UvfVHSRyhrzSwPrAwcZPvxbl0zIuoleSQiuk3dKoVJ+ihwke2fSlqmGvqFpDWBlYC7bN/alYstaniLkhGdaR5iPaorb7yl7fvgJc9fY8LX6aXkkYiemPD9ffx3z2v7Pjh61zcMRR5p1laNTlNb+oi1KWvN/LQhOe0I/MJ2z2cqvflDx/X6EhO21seO4RdXDc86gq998boA3HvF7wccyd+s8PKXctsXR51Xqu9WO2QyE3gOp2HpazNRw5hH7rtmeO7b5ddbl/v+ePWgw3ja8v+wPgD3XXv9gCP5m+XXfSFbHfvFQYfxtIuOPWTQIXRsquWRiZh0H52m+S3eKOlNkpYAPgQsVo1vR9I+wCeAZbsZcES0NpXa1pNHIqJfJlyjU81SqpEOgJLeBhwO3AP8E3A2cALweUlnA6sDu9q+odtBR8TUlDwSEf02maarxVytISNpe+DNwMuqcxwLbEcZ7jlT0tLAUrbv73K8ETGKKVLlnDwSMcSGefRUuybUdCVpZeBsFYsDG1OmXt7U9hOUKZifAN4haXPbjyY5RfSX3f7WD8kjETEIEyro2L6LMuRzJrCM7Y8CJwNHSnp59fnJwALKBF8R0WfD3kcneSRi+A37A1M7JtwZ2fbDwHOA30tajrL+zDzgI5I2sn0HcLztrDkTES0lj0QMt4W47W1YTWrUle0fUNacmA88F/g4cBVwuKSlyDwUEQMzVZaASB6JiH6a9Dw6ts+V9G7gUmAz20dLWsn2Y90PLyImaip1IkweiRhOUymPTFRbEwbanitpSeBCSRu7rDoaETFhySMR0Q9tr3Vl+weSLmgxs2lEDMDCKXgnJo9EDJcpMk3FpHS0qOfIFO0RMXieol1bkkciope6tnp5RAxWKkUiolN1TCMp6ETURB2rnCMiOpWCTkRERAD1rBkemoKOpKUytDSifTXMT5OWPBLRmTrWDE9qwsBekfQPwLslPX/QsUTE1JQ8EhGtDEVBB1gHeAHwz5LWGHAsEVPSsK911QfrkDwSEU0GWtCRNAPKLKnAdcAmwFskrTqB786SNF/S/NmzZ/c40ojhN1WWgOi25JGI7qnjA9NA++jYXggg6VDgTZQVi7cFFkr6tu0FY3x3NjCSmXzzh47rdbgRQ22YE00vdTOP3HdNFk2PqJuBd0aunrp2AHa1/RdJ2wM7AntI+pbtWwcbYcTUMBVnRu6W5JGI7qjj81Lfm64kqWnXyPo2OwPYPge4Hng7pa194IWxiKmgjlXOo0keieiNOjaB97WgI0kja9pIeo2kjYHnAF8H1pK0XXXotcBlwJm2n+xnjBEx3JJHImIy+vKUM5KYGpLT+ylt6HcADwNzgAeA90maBbwY2MX2nf2ILyKGX/JIRH1J2g04lnLfbmp7/ijHbQt8HlgM+KrtE8c7d9+rcyW9Cnid7a0lfRpY1/YcScsBZwIvBf5g+6Z+xxYxlU3FJqh2JY9E9MYAFwf+PfBm4MujHSBpMeBLwDbAAmCepDm2/zDWiXvadKXi5cBZDbufAq6W9AlgA2D3av/GthfYPjfJKWLyFuK2t2GWPBJRf7avsv2ncQ7bFLjW9vW2Hwe+Dew03rl7WtCpapmvAP5V0j9JWppSalsb2BzYx/Zjkg4CTpC0Yi/jiagzu/2tE5JWlHS+pGuqP1cY5bhzJd0n6ceT+3slj0T0SyedkRvnpaq2WV0Obw3g5ob3C6p9Y+pZ09XIqIiqPf1eSpXUFylPX6dRSmGfknQDsBuwh+17ehVPRN0NsOnqSOBC2ydKOrJ6f0SL404Cng0cNNETJ49E9FcnaaRpXqpFSLoAaDWR51G2f9j+lcfWs4JOQ4fBfwXeYPtNkk4D/o9S/fQHSkfCGcCbbV/dq1gipoMBDu/cCdiqen0acBEtCjq2L5S0VfP+sSSPRNSH7ZkdnuIWYK2G92tW+8bU9YJO09DPN1Em8ToQwPZ+kk4HLgX+0fa/d/v6EdNVJzU6VRVzYzXz7OrpbCJWsX1b9fp2YJW2A/lbPMkjEQMw5IMa5gHrSXoBpYCzJ7D3eF/qah+dpuS0EfA8YAvKQnsA2N4XuBH42ch3uhlDREye7dm2N27YnlHIkXSBpN+32HZqOo+hs97NySMR04+kXSQtADYDzpF0XrV/dUlzAar5sA4BzgOuAs6yfeV45+5qjU5DcjoA2Ad4K/Ax4P2SHrR9WXXcrpJWa/xORAyvsaqcJd0haTXbt1X3dUfz1iSPRAzOoJrAbX8f+H6L/bcCb2x4PxeYO5lzd33UlaTXAfsDe7sspvdV4ALgaEmvbjj09m5fO2I6W+j2tw7NAfarXu8HdNypMHkkIrql44LOSJVxNdfFc4BXAOsAuwBUbfffprStvVvSUtX+PIFFdNFo61hNZOvQicA2kq4BZlbvkbSxpK+OHCTpYuA7wNaSFkh6Q8NnySMRQ2CAeaRnOmq6amxLB5ay/VfgZEmPA5tKerPts23fLulrwJO2H+s06IhY1KASje27ga1b7J9P1YG4er9Fq+8nj0REL6kbyVHSoZQORA8Dp9v+WTV516uAn9r+VscXGdvwFiUjOjPhTrbfv/S3bd8Hu7xmw4F35k0eieiZCd/f7/zKd9q+D/7j7bsNPI+00o2mq4Mpk3h9CFgJOE3SDra/DFwJbCZp2U6vM14Y3diqpNqVc9U1pmGLZxhj6nI8EzaomZG7IXlkSv1M1i6eYYxpcHmkfk1Xky7oSJrR8Hopypoz/0xJUgspE4V9XtL2tj8PHGP7wS7F22vdnq66G4YtpmGLB4YvpmGLZ+gkj/TdsMU0bPHA8MU0bPFMWZPuo2N7IUA1f8YjlKethygTeu1m+65q8rETJV1k+75uBhwRrQ374pyNkkciol8mXKPTOCGXpD0pS6m/HvgMZWbCXwGrSToQuByYWXUqjIgAkkcihl0dm64mVKPTNFPp2pROe6+1fZ2kvSmL+C0BPAq8BdjF9h09irmXJjrlfT8NW0zDFg8MX0wDiWeYEw0kjwzYsMU0bPHA8MU0bPFMWeMWdJqS06GUBLQs8BlJt9j+pqSHgC8BvwFOmkJt6c/gia/t0zfDFtOwxQPDF9Og4hnmgk7yyGANW0zDFg8MX0yDimeAiwP3zLgFnYbktDOwMWU69gOBlwGvkXSJ7TmSlgYun6rJKWKq68IMxz2TPBIxNQxxGmnbhProSFoD+ALwhO2rgWOAByijJP6fpMVtn2X7+t6FGhFjGfa29eSRiBiECRV0bN8CvAfYTtJeth8FPgo8AbwBWLJ3IUbU08gyBtNF8khEd023HNKuCQ8vt322pMeAj0vC9rckfQBYwfbDvQsxon4k/QOwo6Qzbd/YjXMOcx+dEckjEd3RixwCUyOPTNak5tGxfY6khcBsSU/a/g5wV29CGx5Na/EMMo5lbD806Dhgkc6lS9p+fNAxNWv+/zYs/x8r6wAvAP65SlS3dHrCqdKJMHlkoDEMTQ6B5JEOrUOXc0hdtTNh4E8k/QtwXQ/iGQqStgbWAJa2PXsYfqglbUlZB+hESTNGJlwbUCyNyWlfYCVJJ9t+clAxtdIQ48GU/58rSTpukAlh5P+d7XMlvRTYBFhc0um2bx9UXP2WPDKQeIYmh1TxJI+0odc5ZKo8ME1GW2td2T6/rh0GJb0B+CxwK2UF5fcOOKQR6wKbwt9mlR2Uhhv/XZQ+F+cMW3IaUcW4C/B1YHPg8EHG47/NCHwo8Cbgr8C2wD6S1uzs3FNrravkkb4bmhxSxZA80oZe5pDq/EM9qKEdHS/qWReSZlQdu94K7FXt/i1w1uCiKk9hkjYE5gOLlNYbZ5rtQywa+VPS8sAWwF62r5a0xKDiGsfqwK7AjsCNwAckLS3pWYMKSNKqlKUOdrV9IPBpYD1gD0mrt3vehXbbW3THMOaRYcohjddLHmlfr3JIXU266arGFrf9mKQbKFPRvxbY1/YtkvYB7rE9t1/BqCx6+CzgMOAxYHlgpqS7KQXUS4CrbN/Qr5hGnsCqP++rYny9pGttP1HFvSVwmQcwD0pzdbykxShVzedQktObbT+psirwo5JO6UdzQos2/burP3cGvlr1WXkp8HbgcUn/0c6TrWs5A8aUMzR5ZBhzCCSPtBlTX3IIDPd8XO1KjQ5P914/TmWyspsoKye/0/afJL2yev9An8Na2fZfbe9ke3fgUMo6QNdRCqhbMYC5nSS9RdI3q7dzgTWBf6w+2wP4ACW59l1Dle4ukl4HrAgcT+mwd7HtxyXtBxwM/KzfhRxJr5G0MfAcShX4WpK2qw69FrgMOHNYq+9jbEOYR4Yyh0DyyGQkh3Ru2tfoSNoG2Bd4FWX15I8BawGnSboc2Ag4yvYlfYzpYGB7SdcA99s+xvYfJV0MPGb7CElLjDz99DiW5ieJOcDHJJ1k+/2SjgYOl/R+SrJ6q+07ex3XaDFK2gs4AfgF8DBwOrAd8E1JmwMvpKyOfW0/YmqI6/2UdvQ7qrjmUH7pvU9lle4XU9Z2avvfbpjbyOtu2PLIMOWQKp7kkTbj6WcOgXrmkWld0KlKxl+ltKfPp/zwHmv7w5K+R1lg8Anbv21xo/Yqpj2B3YB9gBMpPfwXs/0UZbHDDQD6VMh5OhFKWh940PZtVRXpb6q4DpO0IrA+cIP7PHKoKTmtCqxNeVK9l9J08C/AvwMbUiakW9r23a3P1rMYXwW8zvbWkj4NrOuy3MFywJnAS4E/2L6pk+vUscp5Khi2PDJMOaSKJ3mk8/j6kkPqaloXdChPDufb/nn1pPN64GhJHwFOtn3PyIF9KuQsAzxO6dH/RuB5wBttPyVpXeBUYLFex1HFsj5woKTPAk8BnwTOlnSe7TuqqvgbJS1n+wDg0n7E1RRjY3I6nNIf4lXAPNv/I2mkL8QRwCm2z6WMUOhpTJT1mz5M+WUD5d/vakmfoPyS2bnav7HtC4EFvYwpem5o8sgw5ZAqnuSRNuJhgDmkjjU607KPjqQXVTfgL4FXSNqmqiG8EPgzpTT/hurYvvT6Vxm++HZgFeCnwO62/8ml09vbgVnATS5rBPU6FlES5JKUduh7gVOAmZROg6u5dBI8Gdha0ir9+ndq1JCcdqhiezfwI+B4Savbvhk4F/gJcHm/YrJ9BfCvkv5Jpb/G7yk/U5sD+7h0Vj0IOKF6iu3WtWs3LHSYDVseGaYcUsWTPNJmPIPKIXU17Wp0qh/m44EbgPuB7wG7VNWVV1JKy5dQ5pv4Vp9qcg4C9qe0r94iaR1gA0lrU+ZJOAjYu0/NVSNPN5dIejGwDfBB4OPAk8CewNrVzbc+8Brbd/Q6rjHi/XtKAr+jSkiHSPoC8B1Je9q+UdI33IfOeSNJuvr3uxd4M/BFys/UacBOwKdURuTsBuzR+LTfqRRY+mfY8sgw5ZAqnuSR9uIYaA6BTBg45Ul6DWXF5G0o81q8kTKN9iXAOykLDM4CzgdWl/SsXj9hqMzDsB2lmvIxSe+gVD1vSEmkW1ES1JW9jGNEw9PNeylt/A9SZlM9FrgQ+AplpMYmwMcG0ZbetOtW4AfAGiqzq2L7UOCPwKkqQ0Of6kds1ZOYJf0r8D3b76A87f8f5Wnw48A8ytDQN3f7/2nm0emPYcsjw5ZDIHmkXYPOISWGqTXx6ERMtxqdBcC7KAngPcCrgf+gVAkeRqmW3IIyo+meth/pdUC2H6nagE+s4rsKuB74JiUpPNGv2ojqBhOwHCVx7mz7XpU5LXajDPn8gu2LJX22X0+HzTFWr/eo4ryf8svmSeC1khbaPsP2AZJWcemA2c+43kSZyOtAANv7STqd0vfgH23/e6/jiZ4bqjwyLDkEkke6FNO0zCGSdqP8vL4Y2NT2/BbHrEUZAbcKpaA82/bnxzv3tKrRsb3A9jxgS+AM29cBZwB/B/yFMqnWhpTq3yv6GNrplCfA/Wx/gDJkcFNA/SzkVG9Xp4zMWJXyJIjtn1ES6O7Au6unm77P09CQCA4Ejq52H0v5t7uEMhR0hyp5AfR8eGpTgtqI0idhC8qcGyNx70uZaOxnI9/pRSyp0emPIc0jA80hkDzSrmHKIQP2e0pT3c/HOOZJ4HDbGwCvAQ6WtMF4J55uNTojfgccpDLd+JuB97qaD0HSp9zndWBsPwrMU5k+/gDKU+Jeth/u0/VHbrJDKE9g/02pyt1Z0iMuowxuBi4GvtyPWpJGDU+JM4BlKE8777D9C0nnUqrBnwK+ADxBdaP0ul9E4zWq/2/7UIYYfwx4v6QHbV9WHberpNV6GVf66PTd0OSRQeeQKobkkTYMUw6Bwc2wbvsqgLHKcLZvA26rXj8o6SrKrNV/GOvc06pGp8FcyjDLzYATbP8Cnr4RBrnY3dLAQspoid/388KSdqY8ae1D6fD2FPC/wGeqatPjKUNl+9phsOkpcWnbD1BmAF1P0rNd5o34ILBN9f/uzOpm6GeMr6N0BN3b9gLKnCoXUIYYv7rh0J72Q1jo9rdoyzDmkYHlEEge6SC+ocghnZI0S9L8hm1WD6+1DvBKSv+lMU3LGp3qh/w0Sf/lMvTyGTNQDjCuhyWdOqA4lgM+R0lOT9jeHkDSFcCzgaM9gMmoGp523kMZwnsgZejuVsBVkuYBLwIWSlrS9uO9jqmpH8KzgVdQOqPuAvy7y2Ro36ZMYf9uSW+z/digf76iu4Yxjww4h0DyyIQMcw7ppCnb9mxg9mifS7qA0pzZ7CjbP5zodVTmi/oe8J7qPhzTtCzoNHgK+tPEMVEDjOUGyhwXt9reAkDSYcBKwEf7UYAYTZWU9gT2r6q7v6gyGdvBlIUKVwYO6mchp3q7lO2/AidLehzYVNKbbZ9t+3ZJXwOetP1Yr+OCNF0N0FDlkQHHcQPJI+PFMbQ5BHqbR2zP7PQcVVPx94D/sn32RL4zrQs6w5KYhsSvgR9Snmi2oowgeQulc+PAklNlLeCDLmv1LGP7IdsflbQmJYHeZfvWfgTS8GR4KLCZpIeB021/uWpb3lbSUra/5T6v1ZMf58FIHnmG5JFxDHMOGXZVDdjXgKtsf2ai35uufXSiie2HgJOAa4D3U2YI3X8AfYVa9URbm7LWzEicSNoReMT25f0q5IxQWTDxzcCHKAnyNEk72P4yZbK4zSQt28+YIoZB8siE40sOaaKyWvwCSp+3cySdV+1fXX9bhuO1lM7ar5f022p743jnntY1OvFMVee7L0s6pXo/yPkt3kgpiJ9HSQaflvRx2x+UtA9wFGVkR88X1pM0o+qgiKSlKE0V/0zpPLiQsgbO51Xm3vi8pOVdprbvqwwTj2GQPNIypimRQ2BwgxNsfx/4fov9t1Im5cT2JcCkh9anoBOLGERiosz3MZII3kZZlPAe4J+As4ETKIngbMocHbvavqEf8TXEtRPwCOWJ6yHK8NTdbN9VjS44UdJFtu/rR1yLxjmIq0a0ljzyN1Mlh9RVCjoxDBZzNamZpO0pVbovo/x8Hkt54jrT9kyVtXGWsn1/r4NqejLckzKa5FRga+DfgV8Bq1XJ63LgE1XHwoFIjU5Mc0OXR6ZaDoF69vVLH50YKEkrA2erWBzYmDLj7KbVE+HnKJN3vUPS5rYfHUAhZ23KdOOvtX0kZWr/91KexnalzL9xqge4KCGUib7a/S9iKhvGPDIVcwiUgk6727BKQScGyvZdlCGfM4FlbH8UOBk4UtLLq89Ppkwdf20/YmpKUIcC3wE+Avw/SUvb/iZl6vhlgN8AG7q/S4ZERINhyyPJIcMlBZ0YOJdp6p8D/F7ScpQEMA/4iKSNqqec4/s11LIhQe1MeTJ8K/BjSjX4ayQtbnsOpf3/8kF1GmxWxyexiIkapjwyVXMI1HPNvBR0YijY/gFltd75wHOBj1NWYT68GqXQ17tI0hpUa97Yvho4hrJQ4j9TnsoWt32W7ev7GVdEjG6Y8khyyPBIQSeGhsuif+8GLgWWs300cKgHsHyC7VsoCyNuJ2kvl0UTP0pp538DsGQ/45mIrHUVMTx5ZCrmECijN9vdhlVGXcVQsT1X0pLAhZI2tt3zeXLGiOVsSY8BH5eE7W9J+gCwgvu4KvREpQkqohiWPDLVckhdpaATQ8f2DyRd0O9anFFiOUfSQmC2pCdtfwe4a9BxtTLMbeQR/TYseWQq5RCoZx5JQSeGkqsp2oeB7Z9I+hfgukHHMpYhKBdGDJVhySNTJYfUVQo6ERNg+/xBxzCelHMihtdUyCFQzwemFHQiIiICoJYTiGbUVURERNRWanQiaqKOnQgjIjqVgk5ETdSxbT0i+quO82qloBNREynnREQsKgWdiJpYWMNOhBHRX3WsGU5BJ6Im6pigIqK/6phHMuoqoiYGteqwpBUlnS/pmurPFVocs6GkX0q6UtIVkvbo6KIREROUgk5EdOpI4ELb6wEXVu+bPQzsa/slwLbA5yQt38cYI2KaStNVRHRqJ2Cr6vVpwEXAEY0H2L664fWtku4EVgbu60+IETERdZymIjU6ETVht79JmiVpfsM2axKXXsX2bdXr24FVxjpY0qbAkmTdn4jog9ToRNREJ50Ibc8GZo/2uaQLgFVbfHRU03ksadRAJK0GfAPYz/bCNsONiB6pYYVOCjoRddHLKmfbM0f7TNIdklazfVtVkLlzlOOeC5wDHGX70h6FGhEdSNNVRAwt221vHZoD7Fe93g/4YfMBkpYEvg+cbvu7nV4wImKiUtCJqIlO+uh06ERgG0nXADOr90jaWNJXq2N2B14H7C/pt9W2YcdXjoiuGuADU8+k6SoiOmL7bmDrFvvnAwdWr88AzuhzaBERaJhLYRERERGdSNNVRERE1FYKOhEREVFbKehEREREbaWgExEREbWVgk5ERETUVgo6ERERUVsp6ERERERtpaATERERtZWCTkRERNRWCjoRERFRW/8fSs0LhIC7WU8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x792 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax_1, ax_2) = plt.subplots(1, 2, figsize = (9, 11))\n",
    "\n",
    "# generate a mask for the upper triangular matrix\n",
    "mask = np.triu(np.ones_like(non_dp_corr, dtype = np.bool))\n",
    "\n",
    "# generate color palette\n",
    "cmap = sns.diverging_palette(220, 10, as_cmap = True)\n",
    "\n",
    "# get correlation plots\n",
    "ax_1.title.set_text('Non-DP Correlation Matrix')\n",
    "sns.heatmap(non_dp_corr, mask=mask, cmap=cmap, vmax=.3, center=0,\n",
    "                          square=True, linewidths=.5, cbar_kws={\"shrink\": .5}, ax = ax_1)\n",
    "ax_1.set_xticklabels(labels = ['age', 'sex', 'educ', 'income', 'married'], rotation = 45)\n",
    "ax_1.set_yticklabels(labels = ['age', 'sex', 'educ', 'income', 'married'], rotation = 45)\n",
    "\n",
    "\n",
    "ax_2.title.set_text('DP Correlation Matrix')\n",
    "sns.heatmap(dp_corr, mask=mask, cmap=cmap, vmax=.3, center=0,\n",
    "                          square=True, linewidths=.5, cbar_kws={\"shrink\": .5}, ax = ax_2)\n",
    "ax_2.set_xticklabels(labels = ['age', 'sex', 'educ', 'income', 'married'], rotation = 45)\n",
    "ax_2.set_yticklabels(labels = ['age', 'sex', 'educ', 'income', 'married'], rotation = 45)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the differentially private correlation matrix contains values outside of the feasible range for correlations, $[-1, 1]$. This is not uncommon, especially for analyses with small $\\epsilon$, and is not necessarily indicative of a problem. In this scenario, we will not use these correlations for anything other than visualization, so we will leave our result as is.\n",
    "\n",
    "Sometimes, you may get a result that does cause problems for downstream analysis. For example, say your differentially private covariance matrix is not positive semi-definite. There are a number of ways to deal with problems of this type.\n",
    "\n",
    "1. Relax your original plans: For example, if you want to invert your DP covariance matrix and are unable to do so, you could instead take the pseudoinverse.\n",
    "2. Manual Post-Processing: Choose some way to change the output such that it is consistent with what you need for later analyses. This changed output is still differentially private (we will use this idea again in the next section). For example, map all negative variances to small positive value.\n",
    "3. More releases: You could perform the same release again (perhaps with a larger $\\epsilon$) and combine your results in some way until you have a release that works for your purposes.  Note that additional $\\epsilon$ from will be consumed everytime this happens.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Post-Processing of DP Covariance Matrix: Regression Coefficient\n",
    "\n",
    "Differentially private outputs are \"immune\" to post-processing, meaning functions of differentially private releases are also differentially private (provided that the functions are independent of the underlying data in the dataset). This idea provides us with a relatively easy way to generate complex differentially private releases from simpler ones.\n",
    "\n",
    "Say we wanted to run a linear regression of the form $income = \\alpha + \\beta \\cdot educ$ and want to find an differentially private estimate of the slope, $\\hat{\\beta}_{DP}$. We know that \n",
    "$$ \\beta = \\frac{cov(income, educ)}{var(educ)}, $$ \n",
    "and so \n",
    "$$ \\hat{\\beta}_{DP} = \\frac{\\hat{cov}(income, educ)_{DP}}{ \\hat{var}(educ)_{DP} }. $$\n",
    "\n",
    "We already have differentially private estimates of the necessary covariance and variance, so we can plug them in to find $\\hat{\\beta}_{DP}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "income = alpha + beta * educ\n",
      "DP coefficient: -10900.063151950144\n",
      "Non-DP Coefficient: 4601.803740280991\n"
     ]
    }
   ],
   "source": [
    "'''income = alpha + beta * educ'''\n",
    "# find DP estimate of beta\n",
    "beta_hat_dp = dp_cov[2,3] / dp_cov[2,2]\n",
    "beta_hat = non_dp_cov[2,3] / non_dp_cov[2,2]\n",
    "\n",
    "print('income = alpha + beta * educ')\n",
    "print('DP coefficient: {0}'.format(beta_hat_dp))\n",
    "print('Non-DP Coefficient: {0}'.format(beta_hat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "This result is implausible, as it would suggest that an extra year of education is associated with, on average, a decrease in annual income of nearly $11,000. It's not uncommon for this to be the case for DP releases constructed as post-processing from other releases, especially when they involve taking ratios. \n",
    "\n",
    "If you find yourself in such as situation, it is often worth it to spend some extra privacy budget to estimate your quantity of interest using an algorithm optimized for that specific use case."
   ]
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "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.11"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 2
}