{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Aggregation: Sum\n", "\n", "This notebook is a deep-dive on transformations for computing the sum with bounded stability.\n", "\n", "---\n", "Any constructors that have not completed the proof-writing and vetting process may still be accessed if you opt-in to \"contrib\".\n", "Please contact us if you are interested in proof-writing. Thank you!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from opendp.mod import enable_features\n", "enable_features(\"contrib\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Known/Unknown Dataset Size\n", "Constructing a sum transformation is easy!\n", "If the dataset size is not public information, then use `make_bounded_sum`:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from opendp.transformations import *\n", "bs_trans = make_bounded_sum(bounds=(-10, 10))\n", "\n", "# compute the sum\n", "assert bs_trans([1, 2, 4]) == 7 # 1 + 2 + 4\n", "\n", "# compute the sensitivity\n", "assert bs_trans.map(d_in=1) == 10 # d_in * max(|L|, U)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "In the integer case, since the sensitivity is scaled by $max(|L|, U)$, the sensitivity won't increase if you were to widen $L$ to $min(L, -U)$, or widen $U$ to $max(-L, U)$.\n", "This doesn't hold for floating-point datasets with unknown size, for reasons we'll cover in the next section.\n", "\n", "If the dataset size is public information, then use `make_sized_bounded_sum`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "sbs_trans = make_sized_bounded_sum(size=3, bounds=(-10, 10))\n", "\n", "# compute the sum\n", "assert sbs_trans([1, 2, 4]) == 7 # 1 + 2 + 4\n", "\n", "# compute the sensitivity\n", "assert sbs_trans.map(d_in=2) == 20 # (d_in // 2) * (U - L)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Since the sensitivity is scaled by $U - L$, the sensitivity won't increase if you shift both bounds by the same fixed constant. Therefore, the sensitivity remains small in situations where $L$ and $U$ share the same sign and are large in magnitude.\n", "\n", "_All_ sum transformations expect a `d_in` in terms of the `SymmetricDistance`.\n", "However, when the dataset size is known, we are operating in the bounded-DP regime where adjacent datasets have the same size.\n", "This means it takes both an addition and a removal to change any single record, to reach any adjacent dataset.\n", "This results in a stair-step relationship between $d_{in}$ and $d_{out}$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEVCAYAAAD5IL7WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAonUlEQVR4nO3deZgU1Rnv8e+PHcWAwkgExNEooGyDDhCjURQUDSpo3NG45RqNiktiJJrE5WqCS66JZiHGBYy7JEYFY0CUiBsKqICiooiKIuACatTI8t4/qmZshp6ZBqZnKP19nmeerq7tvFXd89bpU3WqFBGYmVn2NGroAMzMbP04gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTeEZIekHSgDpe5xhJl9YyzwJJg+qy3PVVSLzFWLahSfqlpD82dBw1yUKMX0VNGjoAK0xEdG/oGKzBdAceaeggapGFGL9yXAM32/h1B2Y1dBC1yEKMXzlO4BsZSedJelvSx5JeljQwHb9A0iBJR0j6JOfvf5KmpPN0kPR3SUslvS5pRJV195E0M133nUCLAsPqK+lFSR9KuklS5XKSdpQ0RdKytJnnoJxpIWn7nPdrNGOk2/RTSbMkLZd0Z5V1VxtvMba1hn1fyHacm27HfyXdIKm9pH+l63pI0uaF7GhJjST9XNISSe9IOhLYHphTyPJ5tmdclXG/l3RNTdtbnzHaBooI/20kf0BX4C2gQ/q+FPhWOrwAGFRl/m8Ac4EfkRyMZwC/ApoB2wHzgcHpvM2AN4CzgabAocAK4NJaYlpA8o+5NbAF8HjFMul6XgXOT9e/N/Ax0DWdHsD2Oesak1teuu6ngQ7puucCp9QWbzG2tZZ9X8h2PAW0BzoCS4CZQB+SA8fDwIXpvH8C/lRDHBel+/ibQOt0eP56fp+2AT4FNkvfNwYWAd+uaXsLWG+dxei/DftzDXzjsgpoDuwkqWlELIiI1/LNKKkRcBswJSL+AvQFSiLikoj4IiLmA38FjkwX+TZJMvtdRKyIiHHAMwXG9YeIeCsiPgAuA47KWWcrYFRa5sPA+JzphbgmIt5J130/UFZAvMXY1oL3fTWujYjFEfE2MBWYFhHPRsTnwD0kyZyI+HFE/DjfCiSVAD8FfhAR70bEcmACMDtnntMl7VBIQBHxBsmB5OB01N7ApxHx1PpubyExFkrSXpI6r+ty9iUn8I1IRLwKnEVSw1ki6Q5JHaqZ/TJgM6Ci6WAboEPalLFM0jKSmnH7dHoH4O2IyL1/8BsFhvZWlWUqYuoAvBURq6tM71jgegHezRn+lOSAULHu6uKt821dx32fz+Kc4c/yvG9F7QYCc6sk0vbktC1HxB8iYt46xHUbXx5Qj07fb8j21hrjOjiR5NeNrScn8I1MRNwWEbuTJKkALq86T9rmeBRwaESsSEe/BbweEW1y/jaLiO+l0xcBHSUpZ1WF1n62rrLMO+nwO8DW6a+B3Olvp8OfApvkTPtmgeVBzfEWZVtr2Pcbsh3roh1J8wsAkpoCw8hJjjnnO56WdG163uGMGtZ5NzBAUieSmvhtFRMK+a6ta4ySTpP0lKRpkvpXzCPp4Zxl/pOeKzkA+JukYwso1/JwAt+ISOoqaW9JzYHPSWpuq6vM0we4FhgWEUtzJj0NfJyemGopqbGkHpL6ptOfBFYCI9J/qEOAfgWGdpqkTpK2AC4A7kzHTyNJbj9L1zkAOBC4I53+HHB0Gst+wJ4FlldbvHW+rbXs+w3ZjnXxMrC7pC6SWgN/JjnwzE5jbEdSW24HtCX51bEHsH91K0y/I1OAm0gOenPTddX6XVvXGCX1A3YHdgUOB85Nl9mR5PwGktqTHADGAzMiYkBE/K2QnWNrcwLfuDQHRgHvkTQtbAn8vMo8Q4HNgcf05ZUo/4qIVSQ1mjLg9XQd15OcZCIivgAOAY4HPgCOAP5RYFy3ARNJThS+RnIisWKdB5IkkPdITtD9ICJeSpc7M52+DBgO/LPA8mqMt0jbWtO+X+/tqErSaEmj802LiEkkB7/pJG32S0mSa0WTSS+SZN4LuCMiPiZpvlhQS7G3AYPIqX1Tw/amV8+cvx4xDiM5QRskB4PP0sV6As+nw2UktfXtc7bL1pPWbCY0s42VpLNIknUpybmHv0s6CvhGeiK7QUm6EngwIiZLuhx4JiLGSfolyUndiZL+SnLSU8A2EfG7Bgw581wDN8uOniS1154kzToAvXOGG9pfgEskPQFEevUPwKPABZL+L8mvh1kkTTE/lPS7Bon0K8I18K+59DKuF6uZvFNEvFmf8RTT12lb7evBCdzMLKPchGJmllFO4GZmGVWvt5Nt165dlJaW1meRZmaZN2PGjPcioqTq+HpN4KWlpUyfPr0+izQzyzxJeW8F4SYUM7OMcgI3M8soJ3Azs4xq8GdirlixgoULF/L55583dChmRdGiRQs6depE06ZNGzoU+4pp8AS+cOFCNttsM0pLS1nz7p9m2RcRvP/++yxcuJBtt922ocOxr5iCmlAknZ3ed3iOpNsltZC0bXrP31eVPMuw2foE8Pnnn9O2bVsnb/tKkkTbtm39C9OKotYELqkjyVNfyiOiB8lz9Y4kufn71RGxPfAhcNL6BuHkbV9l/n5bsRR6ErMJ0FJSE5Inkywieb5exd3GxpLcC9iq+M53vgPAggULuO22L2/HPH36dEaMGFHdYgCMHj2am2++GYAxY8bwzjvv1Dh/Poceeijz58+vcZ7S0lLee++9dV53rilTpnDAAQfUOM+yZcv405/+tEHljBkzhtNPPx1Yc/9UF9MTTzyxQeWti/Hjx/OrX/2q3sozq7UNPCLelnQV8CbJDdonkjwRfFlErExnW0g1z0GUdDJwMkDnzrU/wat05ISCAi/UglFD6nR966oigVQk8KOPPhqA8vJyysvLa1z2lFNOqRweM2YMPXr0oEOHwh/T+MILL7Bq1Sq222679Yi87lUk8B//OO8zfddZ7v7JZ8qUKbRq1aryIFpsQ4YM4Ze//CUjR45kk002qX0BW0td//8XomqO2BhiKFQhTSibkzwFZluSh8VuCuxXaAERcV1ElEdEeUnJWj1BG9x///tfhgwZQu/evenRowd33pk8LWzGjBnsueee7LLLLgwePJhFixYBMGDAAM477zz69etHly5dmDp1KpAky379+lFWVkavXr2YNy952EirVsmzbEeOHMnUqVMpKyvj6quvrqyxrl69mtLSUpYtW1YZ0w477MDixYu56KKLuOqqqxg3bhzTp09n+PDhlJWVMWHCBIYNG1Y5/6RJkzj44IOp6tZbb2Xo0KGV70899VTKy8vp3r07F1544RrzXnHFFfTs2ZN+/frx6quvAnD33XfTo0cPevfuzR577AEk5yxOOOEEevbsSZ8+fXjkkUfWKrci7go9evRgwYIFjBw5ktdee42ysjLOPTd52taVV15J37596dWr11oxVbjpppvo0qUL/fr14/HHH89bzjXXXMNOO+1Er169OPLII1mwYAGjR4/m6quvpqysjKlTp3L//ffTv39/+vTpw6BBg1i8eHHlek488UQGDBjAdtttxzXXXFNZxs0330yvXr3o3bs3xx6bPLpx6dKlfP/736dv37707du3MiZJDBgwgPHjx+fdDrO6VshVKINInqW3FEDSP4DdgDaSmqS18E58+SDbTHnwwQfp0KEDEyYkR93ly5ezYsUKzjjjDO69915KSkq48847ueCCC7jxxhsBWLlyJU8//TQPPPAAF198MQ899BCjR4/mzDPPZPjw4XzxxResWrVqjXJGjRrFVVddVfnPPWXKFAAaNWrE0KFDueeeezjhhBOYNm0a22yzDe3bt69c9tBDD+UPf/gDV111FeXl5UQEP/nJT1i6dCklJSXcdNNNnHjiiWtt2+OPP85RRx1V+f6yyy5jiy22YNWqVQwcOJBZs2bRq1cvAFq3bs3s2bO5+eabOeussxg/fjyXXHIJ//73v+nYsWPlAeaPf/wjkpg9ezYvvfQS++67L6+88kpB+3rUqFHMmTOH5557DoCJEycyb948nn76aSKCgw46iEcffbTyYAGwaNEiLrzwQmbMmEHr1q3Za6+96NOnT951v/766zRv3pxly5bRpk0bTjnlFFq1asVPf/pTAD788EOeeuopJHH99ddzxRVX8Nvf/haAl156iUceeYSPP/6Yrl27cuqpp/LKK69w6aWX8sQTT9CuXTs++OADAM4880zOPvtsdt99d958800GDx7M3LlzgeSX1dSpUzn88MML2idmG6KQNvA3gW9L2iR9yvdAkpviPwIcms5zHHBvcUIsrp49ezJp0iTOO+88pk6dSuvWrXn55ZeZM2cO++yzD2VlZVx66aUsXLiwcplDDjkEgF122YUFCxYAsOuuu/LrX/+ayy+/nDfeeIOWLVsWHMMRRxxRWfO/4447OOKII2qcXxLHHnsst9xyC8uWLePJJ59k//3Xfq7tokWLyP3Vc9ddd7HzzjvTp08fXnjhBV588ctnG1Qk+qOOOoonn3wSgN12243jjz+ev/71r5UHpMcee4xjjjkGgG7durHNNtsUnMCrmjhxIhMnTqRPnz7svPPOvPTSS5W/XCpMmzaNAQMGUFJSQrNmzardN7169WL48OHccsstNGmSv16ycOFCBg8eTM+ePbnyyit54YUXKqcNGTKE5s2b065dO7bccksWL17Mww8/zGGHHUa7du0A2GKLLQB46KGHOP300ykrK+Oggw7io48+4pNPPgFgyy23XK9zFWbro5A28GmSxgEzSZ70/SxwHclz7e6QdGk67oZiBlosXbp0YebMmTzwwAP84he/YODAgRx88MF07969MpFV1bx5cwAaN27MypXJaYCjjz6a/v37M2HCBL73ve/xl7/8hb333rugGHbddVdeffVVli5dyj//+U9+8Ytf1LrMCSecwIEHHkiLFi047LDD8iatli1bVl6+9vrrr3PVVVfxzDPPsPnmm3P88cevcWlb7pUSFcOjR49m2rRpTJgwgV122YUZM2YUtD1NmjRh9eovH3Be3SV0EcHPf/5zfvSjHxW03ppMmDCBRx99lPvvv5/LLruM2bNnrzXPGWecwTnnnMNBBx3ElClTuOiiiyqnVXymsObnms/q1at56qmnaNGixVrTPv/883U6eJttiIKuQomICyOiW0T0iIhjI+J/ETE/IvpFxPYRcVhE/K/YwRbDO++8wyabbMIxxxzDueeey8yZM+natStLly6tTOArVqxYo7aWz/z589luu+0YMWIEQ4cOZdasWWtM32yzzfj444/zLiuJgw8+mHPOOYcdd9yRtm3brjVP1eU7dOhAhw4duPTSSznhhBPyrnfHHXesbM/+6KOP2HTTTWndujWLFy/mX//61xrzVvwCuPPOO9l1110BeO211+jfvz+XXHIJJSUlvPXWW3z3u9/l1ltvBeCVV17hzTffpGvXrmusq7S0lJkzZwIwc+ZMXn/99bzbMHjwYG688cbK2uvbb7/NkiVL1lhX//79+c9//sP777/PihUruPvuu9faztWrV/PWW2+x1157cfnll7N8+XI++eSTtcpbvnw5HTsm59rHjh2bd5/l2nvvvbn77rt5//33ASqbUPbdd1+uvfbayvkqmoQq9kmPHj1qXbdZXWjwnpgNbfbs2Zx77rk0atSIpk2b8uc//5lmzZoxbtw4RowYwfLly1m5ciVnnXUW3bt3r3Y9d911F3/7299o2rQp3/zmNzn//PPXmN6rVy8aN25M7969Of7449dqxz3iiCPo27cvY8aMybv+448/nlNOOYWWLVvy5JNP0rJlS4YPH87SpUvZcccd8y4zZMgQpkyZwqBBg+jduzd9+vShW7dubL311uy2225rzPvhhx/Sq1cvmjdvzu233w7Aueeey7x584gIBg4cSO/evenWrRunnnoqPXv2pEmTJowZM2aN2ivA97//fW6++Wa6d+9O//796dKlCwBt27Zlt912o0ePHuy///5ceeWVzJ07t/KA0apVK2655Ra23HLLynVttdVWXHTRRey66660adOGsrKytbZz1apVHHPMMSxfvpyIYMSIEbRp04YDDzyQQw89lHvvvZdrr72Wiy66iMMOO4zNN9+cvffeu/LAUp3u3btzwQUXsOeee9K4cWP69OnDmDFjuOaaazjttNPo1asXK1euZI899mD06NEAPPLII/zmN7+pcb1mdaVen4lZXl4eVe8HPnfu3GoTkNXs9NNPp0+fPpx0Uv4+VJ999hl77bUXjz/+OI0bN67n6L5+Fi9ezNFHH83kyZPXmubveWE2hkv4NoYYqpI0IyLWuu7YdyPMqF122YVZs2ZVnlDMp2XLllx88cW8/XYmLxDKnDfffLPyqhaz+vC1b0LJqkJPKA4ePLjIkViFvn37NnQI9jXjGriZWUZtFAm8Ptvhzeqbv99WLA2ewFu0aMH777/vL7l9JVXcDzzfNeNmG6rB28A7derEwoULWbp0aUOHYlYUFU/kMatrDZ7AmzZt6ieVmJmthwZvQjEzs/XjBG5mllFO4GZmGeUEbmaWUU7gZmYZ5QRuZpZRTuBmZhlVyEONu0p6LufvI0lnSdpC0iRJ89LXzesjYDMzS9SawCPi5Ygoi4gyYBfgU+AeYCQwOSJ2ACan783MrJ6saxPKQOC1iHgDGApUPJdqLDCsDuMyM7NarGsCPxK4PR1uHxGL0uF3gfZ1FpWZmdWq4AQuqRlwELDWU2UjuZVg3tsJSjpZ0nRJ033DKjOzurMuNfD9gZkRsTh9v1jSVgDp65J8C0XEdRFRHhHlJSUlGxatmZlVWpcEfhRfNp8A3Acclw4fB9xbV0GZmVntCkrgkjYF9gH+kTN6FLCPpHnAoPS9mZnVk4LuBx4R/wXaVhn3PslVKWZm1gDcE9PMLKOcwM3MMsoJ3Mwso5zAzcwyygnczCyjnMDNzDLKCdzMLKOcwM3MMsoJ3Mwso5zAzcwyygnczCyjnMDNzDLKCdzMLKOcwM3MMsoJ3Mwso5zAzcwyygnczCyjCn2kWhtJ4yS9JGmupF0lbSFpkqR56evmxQ7WzMy+VGgN/PfAgxHRDegNzAVGApMjYgdgcvrezMzqSa0JXFJrYA/gBoCI+CIilgFDgbHpbGOBYcUJ0czM8imkBr4tsBS4SdKzkq5Pn1LfPiIWpfO8C7QvVpBmZra2Qp5K3wTYGTgjIqZJ+j1VmksiIiRFvoUlnQycDNC5c+cNDNfsq6105IR6L3PBqCEbTfm2bgqpgS8EFkbEtPT9OJKEvljSVgDp65J8C0fEdRFRHhHlJSUldRGzmZlRQAKPiHeBtyR1TUcNBF4E7gOOS8cdB9xblAjNzCyvQppQAM4AbpXUDJgPnECS/O+SdBLwBnB4cUI0M7N8CkrgEfEcUJ5n0sA6jcbMzArmnphmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZ5QRuZpZRTuBmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZ5QRuZpZRTuBmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZVdATeSQtAD4GVgErI6Jc0hbAnUApsAA4PCI+LE6YZmZW1brUwPeKiLKIqHi02khgckTsAExO35uZWT3ZkCaUocDYdHgsMGyDozEzs4IVmsADmChphqST03HtI2JROvwu0L7OozMzs2oV1AYO7B4Rb0vaEpgk6aXciRERkiLfgmnCPxmgc+fOGxSsmZl9qaAaeES8nb4uAe4B+gGLJW0FkL4uqWbZ6yKiPCLKS0pK6iZqMzOrPYFL2lTSZhXDwL7AHOA+4Lh0tuOAe4sVpJmZra2QJpT2wD2SKua/LSIelPQMcJekk4A3gMOLF6aZmVVVawKPiPlA7zzj3wcGFiMoMzOrnXtimplllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZVTBCVxSY0nPShqfvt9W0jRJr0q6U1Kz4oVpZmZVrUsN/Exgbs77y4GrI2J74EPgpLoMzMzMalZQApfUCRgCXJ++F7A3MC6dZSwwrAjxmZlZNQqtgf8O+BmwOn3fFlgWESvT9wuBjnUbmpmZ1aTWp9JLOgBYEhEzJA1Y1wIknQycDNC5c+d1Xdy+ZkpHTqj3MheMGrLRlG+2Lgqpge8GHCRpAXAHSdPJ74E2kioOAJ2At/MtHBHXRUR5RJSXlJTUQchmZgYFJPCI+HlEdIqIUuBI4OGIGA48AhyaznYccG/RojQzs7VsyHXg5wHnSHqVpE38hroJyczMClFrG3iuiJgCTEmH5wP96j4kMzMrhHtimplllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZVStCVxSC0lPS3pe0guSLk7HbytpmqRXJd0pqVnxwzUzswqF1MD/B+wdEb2BMmA/Sd8GLgeujojtgQ+Bk4oWpZmZraWQp9JHRHySvm2a/gWwNzAuHT8WGFaMAM3MLL+C2sAlNZb0HLAEmAS8BiyLiJXpLAuBjkWJ0MzM8ioogUfEqogoAzqRPIm+W6EFSDpZ0nRJ05cuXbp+UZqZ2VrW6SqUiFgGPALsCrSR1CSd1Al4u5plrouI8ogoLykp2ZBYzcwsRyFXoZRIapMOtwT2AeaSJPJD09mOA+4tUoxmZpZHk9pnYStgrKTGJAn/rogYL+lF4A5JlwLPAjcUMU4zM6ui1gQeEbOAPnnGzydpDzczswbgnphmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZ5QRuZpZRTuBmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZ5QRuZpZRTuBmZhnlBG5mllFO4GZmGeUEbmaWUU7gZmYZVcgzMbeW9IikFyW9IOnMdPwWkiZJmpe+bl78cM3MrEIhNfCVwE8iYifg28BpknYCRgKTI2IHYHL63szM6kmtCTwiFkXEzHT4Y5In0ncEhgJj09nGAsOKFKOZmeWxTm3gkkpJHnA8DWgfEYvSSe8C7es2NDMzq0mtT6WvIKkV8HfgrIj4SFLltIgISVHNcicDJwN07tx5w6L9iisdOaHey1wwashGU76ZrZuCauCSmpIk71sj4h/p6MWStkqnbwUsybdsRFwXEeURUV5SUlIXMZuZGYVdhSLgBmBuRPy/nEn3Acelw8cB99Z9eGZmVp1CmlB2A44FZkt6Lh13PjAKuEvSScAbwOFFidDMzPKqNYFHxGOAqpk8sG7DMTOzQrknpplZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRjmBm5lllBO4mVlGOYGbmWWUE7iZWUYV8kzMGyUtkTQnZ9wWkiZJmpe+bl7cMM3MrKpCauBjgP2qjBsJTI6IHYDJ6XszM6tHtSbwiHgU+KDK6KHA2HR4LDCsbsMyM7ParG8bePuIWJQOvwu0r6N4zMysQBt8EjMiAojqpks6WdJ0SdOXLl26ocWZmVlqfRP4YklbAaSvS6qbMSKui4jyiCgvKSlZz+LMzKyq9U3g9wHHpcPHAffWTThmZlaoQi4jvB14EugqaaGkk4BRwD6S5gGD0vdmZlaPmtQ2Q0QcVc2kgXUci5mZrQP3xDQzyygncDOzjHICNzPLKCdwM7OMcgI3M8soJ3Azs4xyAjczyygncDOzjHICNzPLKCdwM7OMcgI3M8soJ3Azs4xyAjczyygncDOzjHICNzPLKCdwM7OMcgI3M8uoDUrgkvaT9LKkVyWNrKugzMysduudwCU1Bv4I7A/sBBwlaae6CszMzGq2ITXwfsCrETE/Ir4A7gCG1k1YZmZWG0XE+i0oHQrsFxE/TN8fC/SPiNOrzHcycHL6tivw8vqHW612wHtFWG+WYnD5X+/yN4YYXH7xyt8mIkqqjqz1qfQbKiKuA64rZhmSpkdEeTHL2NhjcPlf7/I3hhhcfv2XvyFNKG8DW+e875SOMzOzerAhCfwZYAdJ20pqBhwJ3Fc3YZmZWW3WuwklIlZKOh34N9AYuDEiXqizyNZNUZtoCtTQMbj8r3f50PAxuPx6tt4nMc3MrGG5J6aZWUY5gZuZZZQTuJlZRhX9OvBikNSNpNdnx3TU28B9ETG34aKqP+n2dwSmRcQnOeP3i4gH6ymGfkBExDPpLRT2A16KiAfqo/w88dwcET9ooLJ3J+mZPCciJtZDef2BuRHxkaSWwEhgZ+BF4NcRsbweYhgB3BMRbxW7rGrKr7jy7Z2IeEjS0cB3gLnAdRGxoh5i2A44hORy6lXAK8BtEfFRscuujCFrJzElnQccRdJ1f2E6uhPJh3lHRIxqqNgAJJ0QETcVcf0jgNNIvqhlwJkRcW86bWZE7FyssnNiuJDkHjhNgElAf+ARYB/g3xFxWZHLr3q5qoC9gIcBIuKgIpf/dET0S4f/D8nncQ+wL3B/sb+Dkl4AeqdXgl0HfAqMAwam4w8pZvlpDMuB/wKvAbcDd0fE0mKXm1P+rSTfv02AZUAr4B8k+0ARcVyRyx8BHAA8CnwPeDaN42DgxxExpZjlV4qITP2RHOWa5hnfDJi3EcT3ZpHXPxtolQ6XAtNJkjjAs/W0jbNJLh3dBPgI+EY6viUwqx7KnwncAgwA9kxfF6XDe9ZD+c/mDD8DlKTDmwKz66H8ubn7osq05+rpO/AsSRPsvsANwFLgQeA4YLN6KH9W+toEWAw0Tt+rnr6Ds3PK3ASYkg53rq//w4jIZBPKaqAD8EaV8Vul04pO0qzqJgHti1x8o0ibTSJigaQBwDhJ26Tl14eVEbEK+FTSa5H+ZIyIzyTVx2dQDpwJXACcGxHPSfosIv5TD2UDNJK0OUkCU6Q1z4j4r6SV9VD+nJxfes9LKo+I6ZK6AEVvOkhFRKwGJgITJTUl+VV2FHAVsNZ9O+pYo7QZZVOSBNoa+ABoDjQtctkVmpA0nTQn+QVARLyZ7ot6CyBrzgImS5oHVLS/dQa2B06vbqE61h4YDHxYZbyAJ4pc9mJJZRHxHEBEfCLpAOBGoGeRy67whaRNIuJTYJeKkZJaUw8H0TRxXC3p7vR1MfX7XW4NzCD5vEPSVhGxSFIr6ucg+kPg95J+QXLzpCclvUXy//DDeigfqmxnJG3O9wH3SdqkHsq/AXiJ5JfgBcDdkuYD3yZpXi2264FnJE0DvgtcDiCphORAUi8y1wYOIKkRyUmj3JOYz6S1wvoo/wbgpoh4LM+02yLi6CKW3YmkBvxunmm7RcTjxSo7p5zmEfG/POPbAVtFxOxix1Cl3CHAbhFxfn2WmyeOTYD2EfF6PZX3DWBbkoPXwohYXB/lpmV3iYhX6qu8amLoABAR70hqAwwiacJ8up7K7w7sSHLy+qX6KHOtGLKYwM3MzNeBm5lllhO4mVlGOYFbrSQNkPSdGqYfVPFQa0lj0qc1rcv6z6/yvqgngiUNk/SrYpaxviS1kfTjWuZZ7/0j6ZP0tYOkcRsSR12S1EzSo5KyeGFFg3ECt0IMIOnlthZJTSLivtiwzitrJPCIqPZgUUd+BvypyGWsrzZA3sRZkdzqYv9ExDsRUdOBtto4iiGS5+pOBo6orzK/CpzAG4ikUkkvpTXWVyTdKmmQpMclzUu7qiOpn6QnJT0r6QlJXdPxZ0u6MR3uKWlO1cu3JDWWdFU6bZakM9LxA9P1zZZ0o6Tm6fgFki6WNDOd1k1SKXAKcLak5yR9N415dHoJ1RWSjpf0h5yiB0manm7XAem615hH0vi0Zj8KaJmu+9Z0WkUtUZKuTOOfLemIdPwASVMkjUv34a2SlE4bJenFdHuvyrPfuwD/i4j30veHpet/XtKj6bhHJZXlLPOYpN6SLpI0VtJUSW9IOkTSFWlsDyq9/jfdj79Jt2m6pJ0l/VvSa5JOyVnvuZKeSWO9OB09CvhWuuyV6bZOVdL79MXc/ZMOn5eW/3y6L6tu77bp92e2pEurfP/mpMPdJT2dljlL0g554mglaXLOd2NoznrmSvqrpBckTVTSvR9J20t6KI1tpqRv1bDdAP8EhlfdBqtBffUY8t9aPblKgZUk1243Irmu+EaS62uHAv9M5/sG0CQdHgT8PR1uRNKN92CS3pi75SnjVJIu1hXLbwG0ILleuEs67mbgrHR4AXBGOvxj4Pp0+CLgpznrHQOM58ueaMcDf8iZ9mAa3w4ktztokTtPOt94YEA6/EmVuD9JX79P0lW/Mcm192+SdNgaACwnuYVCI+BJYHegLclDsyuurmqTZ5+cAPw25/1soGPu/CS9CX+XDncBpufsh8dIOor0JunCvn867R5gWM5+PDUdvhqYBWxG0rllcTp+X5IHACjdhvHAHiTfizk58Q0g6bK+bZ79sz9Jv4NNKj7fPNt7H/CDdPi0nGUrywGuBYanw81IetRWjaMJX/a4bQe8msZeSvI9Lkun3QUckw5PAw5Oh1uQdLjJu93pPI2BpQ39v5mlP9fAG9brETE7ko4pLwCTI/kmzyb5x4Ck08jdaW3paqA7VHZmOR74G/CfyH/99yDgLxGxMl3mA6BrWm7FNbxjSRJHhX+krzNyYsjn7qj+uvu7ImJ1RMwD5gPdalhPTXYHbo+IVZFc4/wfoG867emIWJjuh+fSWJcDnwM3SDqEJMFWtRVJt+8KjwNjlNzTpHHFtgEHpDXqE0kOShX+FUmnlYrbCVTcPCz3M4MvHy84m+SmYx9H0mPzf0quWd43/XuW5NYA3UgOePk8HfmvLR9E0h/hU6j8fKvajeReJZB8V/J5EjhfyX2GtomIz/LMI+DXSnohP0TSB6Oi1/HrkXYsI/3eSNqM5MB4Txrb52mc1W53+n36Il3WCuATBg0rtzPM6pz3q/nys/m/wCMRcXDanDElZ5kdgE9Ibi1Q1zGtoubvx39rmFa1c0GQ1NJyKwwt1j20NeTuu1UkvzJWKml6GggcStIzd+8qy31GclBMAos4Rcnd/YYAMyTtEhHvS5pE8kvocHJ6m1aUGxGrJa1ID7iw5meWG1/u55o7n4DfRMRfcoNLP+OqatrXhaixs0dE3KakOWwI8ICkH5EceHMNJ/kFsUtErJC0gC8/w6qfRcsaisu73TmakxyErQCugW/8WpP0NIWkxg1Udlu/hqT23Fb5r/yYBPxI6ckvSVuQNDGUSto+nedYkpptTT4maQIo1GGSGqVtntulZS4AytLxW5P0pK2wQvnvHzEVOEJJW34JybZW28tOSVf21pHc0vZskmaOquaS3HahYplvRcS0iPgVSc1863TS9ST795mIqHrLhLrwb+DENGYkdZS0Jeu2rycBJyg995F+vlU9TnKnTqimfVnJbVHnR8Q1wL1ArzxxtAaWpMl7L2CbmgKLiI+BhZKGpWU0T+OsbruR1BZ4L+rhVrBfFU7gG78rgN9IepY1a3hXA39Mm0JOAkZV/CPkuJ6k3XiWpOeBoyPic5J24LslzSapEY6uJYb7gYPTE1rfLSDmN0kS7b+AU9IyHwdeJzkRdw3Jz+cK16Ux3lplPfeQtB8/T3Kr2J9FnlsI5NgMGJ/+zH8MOCfPPI8CfSRV3MvjyvSk3ByS9uTnASJiBsmdFotya+BI7ht+G8l9TGaTnKvYLCLeBx5XcmL1ylrW8SBJU810Sc8BP80z25nAaWkZHfNMh+RXxpx0HT2Am/PEcStQnq7nByT3IanNscCI9PN4Avhmddudzr8XMKGA9VrKXenta0fS70nu2/1QDfN0IGmu6pa2s1uRSfoHMDIa+B4rWeIauH0d/Zrkioi8JP2A5AqKC5y864eSW8P+08l73bgGbmaWUa6Bm5lllBO4mVlGOYGbmWWUE7iZWUY5gZuZZZQTuJlZRv1/gvrjT4puKvkAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import pandas as pd\n", "pd.DataFrame(\n", " [(d_in, sbs_trans.map(d_in)) for d_in in range(10)], \n", " columns=[\"max contributions (symmetric distance)\", \"sensitivity (absolute distance)\"]\n", ").plot.bar(0, title=\"sized_bounded_sum: $d_{in}$ vs. $d_{out}$\", width=0.9);" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Intuitively, if we say the symmetric distance between adjacent datasets is at most one, and all adjacent datasets differ by an even number of additions and removals, then the sensitivity is zero." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Floating-Point\n", "\n", "Floating-point addition is not closed, that is, \n", "adding two floating point numbers doesn't necessarily result in another floating-point number.\n", "To resolve this, the IEEE-754 floating-point standard requires a rounding to the nearest floating-point number.\n", "Unfortunately, this influences the sensitivity of the summation.\n", "\n", "In the OpenDP Library, stability maps account for the increased sensitivity due to floating-point rounding in intermediate operations by adding an additional constant term that scales with the dataset size." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20.00000000004426\n" ] } ], "source": [ "sbs_trans = make_sized_bounded_sum(size=1000, bounds=(-10., 10.))\n", "\n", "# The sensitivity is now slightly larger than 20 because of the floating-point constant term\n", "print(sbs_trans.map(d_in=2))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Unfortunately, the worst-case sensitivity analysis hits a snag when the dataset size is unknown, as the rounding error becomes unbounded.\n", "\n", "To keep the sensitivity finite, a dataset truncation operation is applied first:\n", "The dataset size is reduced to no greater than $2^{20}$ elements (a little over 1 million records), if necessary, via a simple random sample.\n", "\n", "The dataset truncation also causes a regression in the sensitivity, as it it is now scaled by $max(|L|, U, U - L)$. \n", "This accounts for the case where an adjacent dataset with one additional row needs to drop a random row to preserve the dataset size.\n", "In practice, the worst-case penalty on the sensitivity is when $U = -L$. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20.00000009313226\n" ] } ], "source": [ "# show the worst-case degradation\n", "bs_trans = make_bounded_sum(bounds=(-10., 10.))\n", "\n", "# the sensitivity is now scaled by max(|L|, U, U - L) = max(|-10.|, 10., 10. - -10.) = 20.\n", "print(bs_trans.map(d_in=1))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "If the bounds share the same sign, then the sensitivity remains unchanged, save for the constant term to account for float rounding." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10.00000009313226\n" ] } ], "source": [ "# if the bounds share sign, the sensitivity is unaffected\n", "bs_trans = make_bounded_sum(bounds=(-10., 0.))\n", "\n", "# the sensitivity is now scaled by max(|L|, U, U - L) = max(|-10.|, 0., 0. - -10.) = 10.\n", "print(bs_trans.map(d_in=1))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Specialized Floating-Point Constructors\n", "\n", "In the previous section an arbitrary limit of ($2^{20}$) on dataset size was baked into the constructor,\n", "to help simplify the library interface.\n", "This limit be manipulated by calling the appropriate constructor:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10.00000000000295" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "make_bounded_float_checked_sum(size_limit=100, bounds=(-10., 0.)).map(d_in=1)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The size of the relaxation term varies according to the dataset size, magnitude of the bounds, summation algorithm and floating-point bit depth.\n", "The following visualization shows the effect of dataset size and choice of algorithm.\n", "\n", "To isolate the relaxation term for this visualization, the sensitivity is calculated for the case when datasets differ by zero additions or removals.\n", "A dataset with the same rows but a different row ordering will result in a different answer for the same sum query." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtZklEQVR4nO3dd3gVVfrA8e+bQi+BhB4giHRCS2iCvSJgA1dQEARl1cWyKrq7uoqKu/DT1VUUFQUBC7BiQ8GGylJEIPQOEQKE3hICIf39/XEv2YiBXCCTyb33/TwPT+7MnJl55ybc955zZs4RVcUYY0zwCnE7AGOMMe6yRGCMMUHOEoExxgQ5SwTGGBPkLBEYY0yQC3M7gLMVFRWlMTExbodhjDF+ZdmyZQdVtUZh2/wuEcTExJCQkOB2GMYY41dEZPvptlnTkDHGBDlLBMYYE+QsERhjTJDzuz6CwmRnZ5OcnExGRobboQSccuXKER0dTXh4uNuhGGMcEhCJIDk5mcqVKxMTE4OIuB1OwFBVDh06RHJyMo0aNXI7HGOMQwKiaSgjI4PIyEhLAsVMRIiMjLSaljEBLiASAWBJwCH2vhoT+AImERhjTKDKy1NemLWenYfTHTm+JYIANXfuXH7++ef85bfeeospU6accZ+RI0fy0ksvOR2aMeYsjf0xkXfmb2P+loOOHD8gOovN782dO5dKlSpx0UUXAXDvvfe6HJEx5lz8d/MB/v3DZm5pX4/+neo7cg6rERST48eP07NnT9q2bUvr1q2ZPn06y5Yt49JLLyUuLo5rr72WPXv2ALBs2TLatm1L27ZtGTFiBK1btwZg0qRJDB8+PP+YvXr1Yu7cuQB89913dO3alQ4dOnDrrbdy7NgxwDPkxjPPPEOHDh2IjY1l48aNJCUl8dZbb/HKK6/Qrl075s+f/5tv+++88w4dO3akbdu29OnTh/R0Z6qbxpjzk3wknYemraBZrcq8cHOsY312AVcjePbLdazffbRYj9mybhWe6d3qjGW++eYb6taty6xZswBITU2lR48efPHFF9SoUYPp06fz5JNPMnHiRO666y5ef/11LrnkEkaMGFHk+Q8ePMioUaOYM2cOFStWZMyYMbz88ss8/fTTAERFRbF8+XLGjRvHSy+9xLvvvsu9995LpUqVeOyxxwD44Ycf8o93yy23cM899wDw1FNPMWHCBB544IFzem+MMc7IzMnl/g+Xk5urvDkgjvJlQh07V8AlArfExsby6KOP8sQTT9CrVy+qVavG2rVrufrqqwHIzc2lTp06pKSkkJKSwiWXXALAwIED+frrr8947F9++YX169fTrVs3ALKysujatWv+9ltuuQWAuLg4Pv300yJjXbt2LU899RQpKSkcO3aMa6+99pyu2RjjnJEz17M6OZW3B8bRKKqio+cKuERQ1Dd3pzRt2pTly5cze/ZsnnrqKa644gpatWrFokWLflMuJSXltMcICwsjLy8vf/nk/fuqytVXX83UqVML3a9s2bIAhIaGkpOTU2SsgwcP5vPPP6dt27ZMmjQpv/nJGFM6TFuyg6lLdnDfZY25tlVtx89nfQTFZPfu3VSoUIEBAwYwYsQIFi9ezIEDB/ITQXZ2NuvWrSMiIoKIiAgWLFgAwIcffph/jJiYGFauXEleXh47d+5kyZIlAHTp0oWFCxeSmJgIePojNm/efMZ4KleuTFpaWqHb0tLSqFOnDtnZ2b85vzHGfSt3pvD0F+u4uEkUj13TrETOGXA1AresWbOGESNGEBISQnh4OG+++SZhYWE8+OCDpKamkpOTw8MPP0yrVq147733GDJkCCLCNddck3+Mbt260ahRI1q2bEmLFi3o0KEDADVq1GDSpEn079+fzMxMAEaNGkXTpk1PG0/v3r3p27cvX3zxBWPHjv3Ntueff57OnTtTo0YNOnfufNqEYYwpWQePZXLfB8uoUbksr/VrT2hIyTzQKapaIicqLvHx8XrqxDQbNmygRYsWLkV0fpKSkujVqxdr1651O5TT8uf31xh/kZObx4AJi1mxI4VP7ruI1vWqFuvxRWSZqsYXts1qBMYYUwqM+WYjv2w9zMt/aFvsSaAo1kfgspiYmFJdGzDGOO/LVbt5Z/42BnVtyC0dokv8/JYIjDHGRRv3HuXxGauJb1iNJ3u2dCUGSwTGGOOS1BPZ/PH9ZVQuF8a4OzpQJsydj2TrIzDGGBfk5Sl/nr6S3SknmDasCzWrlHMtFqsRGGOMC177cQs/btzP071aEtewuquxWCJwydNPP82cOXPOaV9fhpQ2xpReP2zYx7/nbKFvXDQDujR0Oxxnm4ZE5DrgVSAUeFdVR5+mXB9gBtBRVRMKKxNonnvuuULX5+bmEhp65sGlbEhpY/xX4v40Hpq2ktb1qjDqptalYhZAx2oEIhIKvAH0AFoC/UXkd13iIlIZeAhY7FQsJSEpKYnmzZtzxx130KJFC/r27Ut6ejrPPfccHTt2pHXr1gwbNoyTD/ANHjyYGTNmAJ5bSJ944gk6dOjAhx9+SFxcHACrVq1CRNixYwcAjRs3Jj09/TdDSr/22mu0bNmSNm3a0K9fP8AzBMWQIUPo1KkT7du354svvijpt8MYU4jU9GzunpxAufBQxg+Mp1y4cyOKng0nawSdgERV3QogItOAG4H1p5R7HhgDFD0esy++/gvsXVMsh8pXOxZ6FFqZ+Y1NmzYxYcIEunXrxpAhQxg3bhzDhw/PHy564MCBfPXVV/Tu3ft3+0ZGRrJ8+XIAxowZw9GjR5k/fz7x8fHMnz+f7t27U7NmTSpUqPCb/UaPHs22bdsoW7Zs/oB2L7zwAldccQUTJ04kJSWFTp06cdVVV1GxorMjGBpjTi8nN4/hU5ezK+UEU+/pQt2I8m6HlM/JPoJ6wM4Cy8nedflEpANQX1VnORhHialfv37+UNEDBgxgwYIF/PTTT3Tu3JnY2Fh+/PFH1q1bV+i+t912W/7riy66iIULFzJv3jz+9re/MW/ePObPn8/FF1/8u/3atGnDHXfcwQcffEBYmCevf/fdd4wePZp27dpx2WWXkZGRkV+rMMa4Y8w3G5m/5SCjbmpNfIy7ncOncu32UREJAV4GBvtQdhgwDKBBgwZnLuzDN3ennNrWJyLcf//9JCQkUL9+fUaOHJk/tPSpCn5bv+SSS5g/fz7bt2/nxhtvZMyYMYgIPXv2/N1+s2bNYt68eXz55Ze88MILrFmzBlXlk08+oVmzkhm50BhzZp8sS85/cvi2jkV8hrnAyRrBLqDgBJvR3nUnVQZaA3NFJAnoAswUkd8NiqSq41U1XlXja9So4WDI52fHjh35w05/9NFHdO/eHfDMIHbs2LH8PoGiXHzxxXzwwQc0adKEkJAQqlevzuzZs/OPd9LJ4aovv/xyxowZQ2pqav5EM2PHjs3vj1ixYkUxXqUx5mys3JnCXz9bQ9cLInmqlztPDhfFyRrBUqCJiDTCkwD6Abef3KiqqUDUyWURmQs85s93DTVr1ow33niDIUOG0LJlS+677z6OHDlC69atqV27Nh07dvTpODExMahq/ixm3bt3Jzk5mWrVqv2mXG5uLgMGDCA1NRVV5cEHHyQiIoK///3vPPzww7Rp04a8vDwaNWrEV199VezXa4w5s31HMxg2JYFaVcoy7o4OhIeWzjv2HR2GWkSuB/6N5/bRiar6gog8BySo6sxTys7Fh0RQWoeh9ofhpM9VaXh/jfE3Gdm53Db+F7bsS+PT+y+iee0qrsbj2jDUqjobmH3KuqdPU/YyJ2MxxpiSoqo8+dlaVu1M4a0Bca4ngaKUznqKH7LhpI0xJ01YsI1Plifz56uacl1r5+ccPl8Bkwj8baY1f2HvqzFnZ97mA/xj9gZ6tK7NA1dc6HY4PgmIRFCuXDkOHTpkH1rFTFU5dOgQ5cq5NyqiMf4kcf8x/vTRcprWqsxLt7YlpITmHD5fATEMdXR0NMnJyRw4cMDtUAJOuXLliI4u+RmTjPE3R45nMXTyUsqGhfDuoHgqlvWfj1f/ifQMwsPDadSokdthGGOCVFZOHn/8YBl7UjOYNqwL0dUqFL1TKRIQTUPGGOMWzx1Ca1iy7TAv9m1DhwbVit6plLFEYIwx52H8vK18vCyZB69swo3t6hW9QylkicAYY87Rt+v2MvqbjfRsU4eHr2zidjjnzBKBMcacg7W7Unl42kraREfwLz+6Q6gwlgiMMeYs7T+awT1TEoioEM47A+NKzQQz5yog7hoyxpiSciIrl3umJJB6IpuP7+1KzSr+/5yNJQJjjPFRXp7y6McrWb0rlfED42lVt6rbIRULaxoyxhgfvfz9Zmav2ctfezTn6pa13A6n2FgiMMYYH0xfuoPXf0qkX8f63HPxBW6HU6wsERhjTBHmbT7A3z5by8VNonj+pta/m5bW31kiMMaYM1i/+yj3f7icJjUrlepZxs5H4F2RMcYUk72pGQyZtJRKZcN4766OVC4X7nZIjrBEYIwxhUjLyOauSUtJy8hm4uCO1Kla3u2QHGO3jxpjzCmyc/P400cr2LwvjYmDO9KybumeavJ8WY3AGGMKUFWe/mIt8zYfYNRNrbm0aQ23Q3KcJQJjjClg3NxfmbpkJ/df1pj+nRq4HU6JsERgjDFeX6zcxYvfbuKGtnV57JpmbodTYiwRGGMMsHjrIUZ8vJpOMdV58dY2fj2a6NmyRGCMCXob9x7l7ikJRFcvz/g74ygb5t+jiZ4tSwTGmKC2K+UEgyYuoUKZUKYM6UREhTJuh1Ti7PZRY0zQOnI8izsnLCY9K5eP7+3qd5POFxerERhjgtKJrFyGTl7KzsMneOfOeJrXDuxnBc7EagTGmKCTk5vHA1NXsGJnCm/c3oEuF0S6HZKrrEZgjAkqqsrfv1jLnA37GNm7FdfH1nE7JNdZIjDGBJVX5mxh6pKd/Onyxgy6KMbtcEoFSwTGmKDxwS/bee2HLfwhPjqoHhgriiUCY0xQ+GbtXp7+Yi1XNK/JP26ODbjJZc6HJQJjTMBbvPUQD05bQZvoCF6/vT1hATi5zPmwd8MYE9DW7krl7skJ1K9WnomDO1KhjN0seSpLBMaYgJW4/xh3TlxClfLhfHB3Z6pXDL6nhn1hicAYE5B2pZzgzgmLCRF4f2ingJ5h7Hw5mghE5DoR2SQiiSLyl0K23ysia0RkpYgsEJGWTsZjjAkOB49lMvDdxaRl5jB5SCcuqFHJ7ZBKNccSgYiEAm8APYCWQP9CPug/UtVYVW0H/B/wslPxGGOCw9GMbAZNXMLu1BNMHNyRVnWruh1SqedkjaATkKiqW1U1C5gG3FiwgKoeLbBYEVAH4zHGBLiM7FzunpzApr1pvDkgjo4x1d0OyS842X1eD9hZYDkZ6HxqIRH5E/AIUAa4orADicgwYBhAgwbBMXWcMebsZOfmcf+Hy1madJjX+rXn8mY13Q7Jb7jeWayqb6hqY+AJ4KnTlBmvqvGqGl+jRuBPJG2MOTt5ecpjH6/ix437GXVTa3q3ret2SH7FyUSwC6hfYDnau+50pgE3ORiPMSYAqSojv1zHFyt38/h1zbijc0O3Q/I7TiaCpUATEWkkImWAfsDMggVEpEmBxZ7AFgfjMcYEGFVl9DcbmbJoO3+85ALuu7Sx2yH5Jcf6CFQ1R0SGA98CocBEVV0nIs8BCao6ExguIlcB2cARYJBT8RhjAs+rP2zh7f9uZWCXhvylR3MbP+gcOfqstarOBmafsu7pAq8fcvL8xpjA9dZ/f+Xfc7Zwa1w0z97QypLAeXC9s9gYY87Wewu3MfrrjdzQti6j+7QhJMSSwPmwRGCM8StTl+zg2S/Xc22rWvzrD20JtSRw3iwRGGP8xmcrkvnbZ2u4vFkNxvbvQLgNJ10s7F00xviFWav38Oh/VtH1gkjeHBBHmTD7+Cou9k4aY0q979fv46FpK4hrWI13B8VTLjzU7ZACiiUCY0ypNm/zAf704XJa1a1iE8s4xKdEICKfikhPEbHEYYwpMfM2H+CeKQk0rlmJyUM6UblcuNshBSRfP9jHAbcDW0RktIg0czAmY4zJTwIX1KjEh3d3JqKCzS7mFJ8SgarOUdU7gA5AEjBHRH4WkbtExFK0MaZYzd/y2yRgU0w6y+emHhGJBAYDdwMrgFfxJIbvHYnMGBOU5m85wN2TLQmUJJ96XUTkM6AZ8D7QW1X3eDdNF5EEp4IzxgSXk0mgUVRFSwIlyNfu93e84wblE5GyqpqpqvEOxGWMCTILthzMTwIf3dPFkkAJ8rVpaFQh6xYVZyDGmOC1YMtBhk5eaknAJWesEYhIbTxTTpYXkfbAyUE9qgAVHI7NGBMELAm4r6imoWvxdBBHAy8XWJ8G/M2hmIwxQeK/mw8wbIo1B7ntjIlAVScDk0Wkj6p+UkIxGWOCwPfr9/GnD5dzYc1KfGAdw64qqmlogKp+AMSIyCOnblfVlwvZzRhjzmjW6j08NG0FretVZfJdnahawR5HclNRTUMVvT8rOR2IMSY4fLo8mcc+XkVcw2pMHNzRho0oBYpqGnrb+3Kcqh4ogXiMMQHso8U7ePLzNVzUOJJ37oy3AeRKCV9vH10oIt+JyFARqeZoRMaYgPTewm387bM1XNa0BhMG2SiipYmvYw01BZ4CWgHLROQrERngaGTGmIDx5txf86eXfHugzSdQ2vg81pCqLlHVR4BOwGFgsmNRGWMCgqryyvebGfONZ6L512/vYDOLlUK+zkdQRUQGicjXwM/AHjwJwRhjCqWqjP56I6/+sIVb46J55bZ2NsdwKeVrI90q4HPgOVW1oSWMMWeUm6c8+dkapi3dycAuDXn2hlaEhEjROxpX+JoILlBVdTQSY0xAyMzJ5c/TVzJ7zV4euOJCHrm6KSKWBEqzoh4o+7eqPgzMFJHfJQJVvcGpwIwx/ud4Zg73frCM+VsO8lTPFtx98QVuh2R8UFSN4H3vz5ecDsQY499S0rMY/N5SVien8GLfNtwaX9/tkIyPinqgbJn3ZTtVfbXgNhF5CPivU4EZY/zHvqMZ3DlhCdsOHufNAXFc26q22yGZs+BrF/6gQtYNLsY4jDF+avuh4/R962eSj6Tz3l0dLQn4oaL6CPoDtwONRGRmgU2V8TxLYIwJYhv2HOXOiUvIyc3jo3u60LZ+hNshmXNQVB/ByWcGooB/FVifBqx2KihjTOm3NOkwQyctpUKZMD76Y1ea1KrsdkjmHBXVR7Ad2A50LZlwjDH+4Ju1e3hw2kqiI8ozeUgn6le3CQv9WVFNQwtUtbuIpAEFbx8VQFW1iqPRGWNKnSmLknhm5jra1Y9gwqCONqFMACiqRtDd+9PqfMYEOVXlxW83MW7ur1zVoiZj+3egfBkbPC4Q+DrWUGMRKet9fZmIPCgiEY5GZowpNbJz83j041WMm/sr/TvV560BcZYEAoivt49+AuSKyIXAeKA+8JFjURljSo1jmTkMnZzAp8t38cjVTfnHzbGE2eBxAcXX32aequYANwNjVXUEUKeonUTkOhHZJCKJIvKXQrY/IiLrRWS1iPwgIg3PLnxjjJMOpGXSf/wvLEw8yJg+sTx4ZRMbNygA+ZoIsr3PFAwCvvKuO+NEoyISCrwB9ABaAv1FpOUpxVYA8araBpgB/J+vgRtjnLX1wDH6vPkzifuP8c6dcdzWsYHbIRmH+JoI7sJzC+kLqrpNRBrxv3GITqcTkKiqW1U1C5gG3FiwgKr+pKrp3sVfgGjfQzfGOGXx1kPc8ubPHMvMYeqwLlzRvJbbIRkH+TQMtaquBx4ssLwNGFPEbvWAnQWWk4HOZyg/FPi6sA0iMgwYBtCggX0rMcZJn61I5vEZq6lfvQLvDe5Iw8iKbodkHOZTIhCRbsBIoKF3n5PPERTLGLPe+Y/jgUsL266q4/F0UhMfH2/zIhjjAFXl33O28OoPW+h6QSRvDYijaoUztgCbAOHrxDQTgD8Dy4BcH/fZhefuopOivet+Q0SuAp4ELlXVTB+PbYwpRpk5uTwxYzWfr9xN37ho/nFzrM0tHER8TQSpqlpos80ZLAWaePsTdgH98Axgl09E2gNvA9ep6v6zPL4xphgcOZ7FH99fxpKkw4y4thn3X9bY7gwKMr4mgp9E5EXgUyD/W7uqLj/dDqqaIyLDgW+BUGCiqq4TkeeABFWdCbwIVAI+9v7h7bBZz4wpOdsOHueu95awOzWDsf3b07ttXbdDMi7wNRGc7OSNL7BOgSvOtJOqzgZmn7Lu6QKvr/Lx/MaYYrZ46yH++MEyQkSYek9n4hpWdzsk4xJf7xq63OlAjDElZ+qSHfz987U0iLQ7g4zvYw3VEpEJIvK1d7mliAx1NjRjTHHLyc3jmS/W8tdP19Dtwig+u7+bJQHj8wNlk/C09Z9sQNwMPOxAPMYYh6SkZzHovSVMXrSdu7s3YuLgjlQtb7eHGt/7CKJU9T8i8lfI7wj29TZSY4zLEvencffkBHanZPBi3zbcGl+/6J1M0PA1ERwXkUi8k9OISBcg1bGojDHF5qeN+3lg6grKhYcydZh1Cpvf8zURPALMBBqLyEKgBtDXsaiMMedNVRk/byujv9lIq7pVGD8wnroR5d0Oy5RCviaCxnhGEa0P9MFzO6mv+xpjSlhGdi5/+3QNn67YRc82dXipb1ubSMaclq+dxX9X1aNANeByYBzwpmNRGWPO2c7D6dwy7mc+W7mLR69uyuv921sSMGfk67f6kx3DPYF3VHWWiIxyKCZjzDmau2k/D01biaoycXBHLm9W0+2QjB/wNRHsEpG3gauBMd75i21EKmNKibw85Y2fEnl5zmaa1arM2wPj7PkA4zNfE8EfgOuAl1Q1RUTqACOcC8sY46ujGdk8Mn0Vczbs46Z2dfnnLW2sKcicFV+HmEjHM+DcyeU9wB6ngjLG+GbT3jTu/WAZOw+nM7J3SwZdFGMjh5qzZnf+GOOnvlq9m8dnrKZi2TCmDutCxxh7PsCcG0sExviZzJxc/jl7I5N+TiK+YTXG3dGBmlXKuR2W8WOWCIzxIzsOpTN86nJWJ6dyd/dGPH5dc5tJzJw3SwTG+Ilv1u5lxIxVCDB+YBzXtKrtdkgmQFgiMKaUy8rJ459fb+C9hUm0ja7K67d3oH71Cm6HZQKIJQJjSrGdh9MZPnUFq3amcFe3GP7ao4U1BZliZ4nAmFLqu3V7eezjVSjw1oAOXNe6jtshmQBlicCYUiYjO5d/zN7AlEXbia1XlTdu70CDSGsKMs6xRGBMKbJpbxoPTl3Bpn1pDO3eiMeva0bZMHtK2DjLEoExpYCq8v4v2xk1awNVyoUzeUgnLm1aw+2wTJCwRGCMyw4dy+TxGav5YeN+LmtWg5dubUtUpbJuh2WCiCUCY1w0f8sBHvnPKlLTs3mmd0sG21hBxgWWCIxxQUZ2Li9/v5nx87ZyYc1KTL6rEy3rVnE7LBOkLBEYU8LW7U7lkemr2LQvjTs6N+Cpni1t2GjjKksExpSQnNw83vrvr/x7zhaqVyzDe4M7cnlzm0HMuM8SgTEl4NcDx3j0P6tYuTOF3m3r8vyNrYioUMbtsIwBLBEY46i8PGXKoiRGf7ORcuGhjO3fnt5t67odljG/YYnAGIfsSjnB4zNWsTDxEJc3q8GYPm1s3gBTKlkiMKaY5eUpU5fu4J+zN6Kq/POWWPp1rG+3hZpSyxKBMcUo6eBx/vLpan7ZephuF0Yy+pY2NmS0KfUsERhTDHLzlAkLtvKv7zZTJiyEMX1i+UO81QKMf7BEYMx52rQ3jcdnrGJVcipXtajFCze3ppb1BRg/YonAmHOUlZPHuLmJvPFTIlXKhTO2f3t6taljtQDjdxyd6khErhORTSKSKCJ/KWT7JSKyXERyRKSvk7EYU5yWJh2m19j5/HvOFq6PrcP3j1xK77Z1LQkYv+RYjUBEQoE3gKuBZGCpiMxU1fUFiu0ABgOPORWHMcXp8PEsRn+9gf8kJFMvojwTBsVzZYtabodlzHlxsmmoE5CoqlsBRGQacCOQnwhUNcm7Lc/BOIw5b6rKjGXJ/GP2BtIycvjjpRfw0JVNqFDGWleN/3Pyr7gesLPAcjLQ+VwOJCLDgGEADRo0OP/IjDkLW/al8eTna1my7TDxDasx6ubWNK9tI4WawOEXX2dUdTwwHiA+Pl5dDscEiRNZuYz9cQvj522lUrkwxvSJ5da4+oSEWD+ACSxOJoJdQP0Cy9HedcaUaqrK12v38sKsDexKOUHfuGj+2qM5kTZrmAlQTiaCpUATEWmEJwH0A2538HzGnLeNe4/y7Mz1LNp6iOa1KzNtWBe6XBDpdljGOMqxRKCqOSIyHPgWCAUmquo6EXkOSFDVmSLSEfgMqAb0FpFnVbWVUzEZczop6Vm88v1m3v9lO1XKh/P8Ta3p37E+YaGO3mFtTKngaB+Bqs4GZp+y7ukCr5fiaTIyxhW5ecrUJTv413ebSD2RzR2dG/LI1U2pVtHmCjDBwy86i41xws+JBxk1awPr9xylc6PqjLyhFS3q2N1AJvhYIjBBZ8u+NP759UZ+3LifehHlef329vSMtaEhTPCyRGCCxv6jGbwyZzPTl+6kYtkw/tKjOYMviqFcuE0cb4KbJQIT8I5n5vDO/K2Mn7eVrJw87uwaw4NXNqG69QMYA1giMAEsOzePGcuSefn7zRxIy+T62No8fm1zYqIquh2aMaWKJQITcPLylC9X7+aV7zeTdCidDg0ieGtAB+IaVnc7NGNKJUsEJmCoKt+t38fL321m0740mteuzDt3xnNVi5rWEWzMGVgiMH5PVZm/5SD/+m4Tq5JTuSCqIq/1b0+v2Do2LpAxPrBEYPza4q2H+Nf3m1my7TD1Isrzf33acEuHevZEsDFnwRKB8TuqyoLEg4z9IZElSYepUbksz97Qin6d6lM2zG4FNeZsWSIwfkNV+XHjfsb+mMjKnSnUrlKOkb1b0q9TA3sWwJjzYInAlHp5ecq36/Yy9sdE1u85SnS18vzj5lj6xNWzGoAxxcASgSm1MnNymblyN+PnbWXL/mM0iqrIi33bcFP7eoRbH4AxxcYSgSl1Uk9k89HiHby3cBv70zJpXrsyr/ZrR682dQm1u4CMKXaWCEypkXwknYkLkpi+dAfHs3LpfmEUL93aloubRNlzAMY4yBKBcd3q5BTenb+NWWv2IEDvtnW5++JGtKpb1e3QjAkKlgiMKzJzcvl6zV4mL0pixY4UKpUNY0i3GO7q1oi6EeXdDs+YoGKJwJSoPakn+GjxDqYu2cHBY1k0iqrIM71b0icumirlwt0Oz5igZInAOE5VWbztMFMWJfHtun3kqXJl85rc2TWG7hdG2TAQxrjMEoFxzMFjmXy2fBfTlu7g1wPHqVo+nKHdGzGwS0PqV6/gdnjGGC9LBKZY5eYp87ccYPrSnczZsI/sXKVDgwj+r08beretS/ky9gCYMaWNJQJTLHYeTmfGsmQ+TtjJ7tQMqlUIZ1DXGG7rWJ8mtSq7HZ4x5gwsEZhzlpKexaw1e/h8xS6WJh1BBLpfGMWTPVtyVcuaNvyDMX7CEoE5KxnZufy0cT+frdjFT5v2k52rXFizEiOubcaN7eoSXc3a/o3xN5YITJGyc/P4+ddDzF69h9lr95CWkUONymUZ1DWGm9rXo1XdKvbkrzF+zBKBKVRmTi4Lthxk9pq9zNmwj9QT2VQqG8Y1rWpxc/t6XNQ4ysb9MSZAWCIw+U5k5TJvywG+XrOHHzbsJy0zhyrlwri6ZW16tK5N9yZRNu6/MQHIEkGQSz6Szk8b9/Pjxv38/OshMnPyiKgQzvWxdegRW5uLGkdRJsyGfDYmkFkiCDK5ecrKnUf4YYPnw3/j3jQAGkZW4I7ODbmyRU06Napu4/0bE0QsEQQ4VWXbweMsTDzIgsSDLPr1EEczcggNETrGVOPJ61twRYuaXBBV0Tp8jQlSlggC0P60DBb9eogFWw6yMPEgu1MzAKgXUZ4erevQvUkUlzStQdXyNsibMcYSgd9TVbYePE5C0mGWJh0hIekwSYfSAahaPpyLGkdy/+VRdL8wioaRFexbvzHmdywR+JljmTms25XKquQUEpKOkLD9CIePZwFQvWIZ4hpWo3+nBnRtHEmrulXtFk9jTJEsEZRi6Vk5rN99lNXJqazZlcrq5BS2HjyOqmd7TGQFrmhek44x1YiPqW7t/MaYc2KJoBTIyc1j++F0Nu9NY9O+NDbvS2PT3jS2HTxOnvdDv1aVssTWi+DGdvWIrVeV2OiqRFUq627gxpiAYImghKgqh49nkXQone2Hjuf/3LLvGIkHjpGVkweACDSsXoGmtSrTM7YObaIjiI2uSq0q5Vy+AmNMoHI0EYjIdcCrQCjwrqqOPmV7WWAKEAccAm5T1SQnY3JKTm4e+9My2ZOawd7UDPaknvD+zGD74eNsP5hOWmZOfvkQgboR5bmwZiUubhJF01qVaVa7Mo1rVLIx+40xJcqxRCAiocAbwNVAMrBURGaq6voCxYYCR1T1QhHpB4wBbnMqpqKoKhnZeRzPyiE9M9fzMyuH45m5pJ7I5kh6FoePZ3HkeBaH07M5cjyLI+lZHDyWyYG0zPxmnJPKhYdQt2p5GkRWIL5hdRpGViAmsiINIysQXa2CPbFrjCkVnKwRdAISVXUrgIhMA24ECiaCG4GR3tczgNdFRFT1lI/U87d50p9I37kSVFFAFe/P/y3nnea0ZYAa3n8AoSFCeEgIYaFCeKgQFhZCmRohlAnz/Csb6vkZGiII3s7bQ95/xhhzrmrHQo/RRZc7S04mgnrAzgLLyUDn05VR1RwRSQUigYMFC4nIMGAYQIMGDc4pmLLhIVAmFPEcz/sTBPH+hJAQITRECBHPz1DxrhMhNPR/H/75H+7GGBMA/KKzWFXHA+MB4uPjz6m20PCOscUakzHGBAonG6l3AfULLEd71xVaRkTCgKpYA4oxxpQoJxPBUqCJiDQSkTJAP2DmKWVmAoO8r/sCPzrRP2CMMeb0HGsa8rb5Dwe+xXP76ERVXScizwEJqjoTmAC8LyKJwGE8ycIYY0wJcrSPQFVnA7NPWfd0gdcZwK1OxmCMMebM7EZ2Y4wJcpYIjDEmyFkiMMaYIGeJwBhjgpz4292aInIA2H6Ou0dxylPLQcCuOTjYNQeH87nmhqpao7ANfpcIzoeIJKhqvNtxlCS75uBg1xwcnLpmaxoyxpggZ4nAGGOCXLAlgvFuB+ACu+bgYNccHBy55qDqIzDGGPN7wVYjMMYYcwpLBMYYE+QCMhGIyHUisklEEkXkL4VsLysi073bF4tIjAthFisfrvkREVkvIqtF5AcRaehGnMWpqGsuUK6PiKiI+P2thr5cs4j8wfu7XiciH5V0jMXNh7/tBiLyk4is8P59X+9GnMVFRCaKyH4RWXua7SIir3nfj9Ui0uG8T6qqAfUPz5DXvwIX4JlueBXQ8pQy9wNveV/3A6a7HXcJXPPlQAXv6/uC4Zq95SoD84BfgHi34y6B33MTYAVQzbtc0+24S+CaxwP3eV+3BJLcjvs8r/kSoAOw9jTbrwe+xjPDbhdg8fmeMxBrBJ2ARFXdqqpZwDTgxlPK3AhM9r6eAVwpIv48EXGR16yqP6lqunfxFzwzxvkzX37PAM8DY4CMkgzOIb5c8z3AG6p6BEBV95dwjMXNl2tWoIr3dVVgdwnGV+xUdR6e+VlO50Zginr8AkSISJ3zOWcgJoJ6wM4Cy8nedYWWUdUcIBWILJHonOHLNRc0FM83Cn9W5DV7q8z1VXVWSQbmIF9+z02BpiKyUER+EZHrSiw6Z/hyzSOBASKSjGf+kwdKJjTXnO3/9yL5xeT1pviIyAAgHrjU7VicJCIhwMvAYJdDKWlheJqHLsNT65snIrGqmuJmUA7rD0xS1X+JSFc8sx62VtU8twPzF4FYI9gF1C+wHO1dV2gZEQnDU508VCLROcOXa0ZErgKeBG5Q1cwSis0pRV1zZaA1MFdEkvC0pc708w5jX37PycBMVc1W1W3AZjyJwV/5cs1Dgf8AqOoioByewdkClU//389GICaCpUATEWkkImXwdAbPPKXMTGCQ93Vf4Ef19sL4qSKvWUTaA2/jSQL+3m4MRVyzqqaqapSqxqhqDJ5+kRtUNcGdcIuFL3/bn+OpDSAiUXiairaWYIzFzZdr3gFcCSAiLfAkggMlGmXJmgnc6b17qAuQqqp7zueAAdc0pKo5IjIc+BbPHQcTVXWdiDwHJKjqTGACnupjIp5OmX7uRXz+fLzmF4FKwMfefvEdqnqDa0GfJx+vOaD4eM3fAteIyHogFxihqn5b2/Xxmh8F3hGRP+PpOB7sz1/sRGQqnmQe5e33eAYIB1DVt/D0g1wPJALpwF3nfU4/fr+MMcYUg0BsGjLGGHMWLBEYY0yQs0RgjDFBzhKBMcYEOUsExhhTihU1CN0pZV8RkZXef5tFJMWXc1giMAFPREaKyGNFlLlJRFoW83ljROT2s9xntohEFGccxu9NAnwaKkRV/6yq7VS1HTAW+NSX/SwRGONxE56RK4tTDHBWiUBVrw/w4SDMWSpsEDoRaSwi34jIMhGZLyLNC9m1PzDVl3NYIjABSUSe9FaNFwDNCqy/R0SWisgqEflERCqIyEXADcCL3ip148LKefe/VUTWetfP864LFZEXveVXi8gfvacbDVzsPeafT4mvjojM825bKyIXe9cniUiUiNxboIq/TUR+8m6/RkQWichyEflYRCo5/maa0mg88ICqxgGPAeMKbhTPfCONgB99OZg9UGYCjojE4alOd8bz9PxyPPNPvCQikSeftBWRUcA+VR0rIpOAr1R1hnfb6cqtAa5T1V0iEqGqKSIyDM+4/6NEpCywELgVaAg8pqq9ConxUaCcqr4gIqF45opI846LFK+qB73lwvH8Z/4/YBGeqn4PVT0uIk8AZVX1ueJ/F01pIp7Js75S1dbe5H8A2FSgSFlVbVGg/BNAtKr6NBJrwA0xYQxwMfDZyfkXRKTgcBOtvR/sEXiG3Pj2NMc4XbmFwCQR+Q//a3+9BmgjIn29y1XxDPSWdYYYlwITvR/0n6vqytOUexXPWFhfikgvPM1XC73DhJTBkxxMcAkBUrz9AKfTD/jT2RzQmGAyCRiuqrHAs3gGKPO5nKreCzyFZ/THZSISiWemqAdOdtKpaiNV/e5MQXjbfS/BM2rkJBG589QyIjIYT63i2ZOrgO8LnKelqg71/dJNIFDVo8A2EbkV8qeubHtyu7e/oBpn8SXBEoEJRPOAm0SkvIhUBnoX2FYZ2OP9Jn5HgfVp3m1nLCcijVV1sao+jad6Xh9PbeE+b1lEpKmIVCzkmBQ4TkM8zU3vAO/imZqw4PaTbb8DCoyr/wvQTUQu9JapKCJNfX5XjF/yDkK3CGgmIskiMhTP3+RQEVkFrOO3s7b1A6adzcB71jRkAo6qLheR6Xjmt92PpxnmpL8Di/F8iC/mfx/U0/CMYPkgnqHJT1fuRRFpgufb+Q/ec6zGc4fQcvG02RzAcxfSaiDX+591kqq+UiCOy4ARIpINHANOrREMB6oDP3mbgRJU9W5vLWGqty8CPLWTzWf3Dhl/oqr9T7Op0FtKVXXk2Z7DOouNMSbIWdOQMcYEOUsExhgT5CwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhjTJD7f5VtPh1CaCl8AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "bounds = (0., 10.)\n", "sizes = list(range(1, 10_000_000, 10_000))\n", "\n", "pd.DataFrame({\n", " \"dataset size\": sizes,\n", " \"sequential\": [make_sized_bounded_float_checked_sum(size, bounds, S=\"Sequential\").map(0) for size in sizes],\n", " \"pairwise\": [make_sized_bounded_float_checked_sum(size, bounds, S=\"Pairwise\").map(0) for size in sizes],\n", "}).plot(0, ylabel=\"sensitivity\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The increase to sensitivity when using a sequential summation algorithm grows on the order of $O(n^2)$, while the pairwise algorithm grows on the order of $O(n \\log_2(n))$.\n", "\n", "OpenDP defaults to the pairwise algorithm, but the ability to configure the algorithm can be useful to calculate the sensitivity in situations where you don't have control over how the summation is computed.\n", "For example, floating-point aggregations in SQLite and MySQL both exhibit increases in sensitivity akin to the sequential algorithm.\n", "\n", "Beware, these relaxation terms grow far more quickly when the data type is adjusted to single-precision floats (`f32`)!\n", "\n", "\n", "\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Specialized Integer Constructors\n", "\n", "Just as in the case for floating-point types, there are specialized constructors for summation over integral types.\n", "\n", "The integral transformations from `make_bounded_sum` and `make_sized_bounded_sum` each use properties of the bounds, data types and input metric to determine which strategy to use to compute the sum in a way that protects the sensitivity from the effects of numerical overflow.\n", "\n", "The following strategies are ordered by computational efficiency:\n", "\n", "* ``checked`` can be used when the dataset size multiplied by the bounds doesn't overflow.\n", "* ``monotonic`` can be used when the bounds share the same sign.\n", "* ``ordered`` can be used when the input metric is ``InsertDeleteDistance``.\n", "* ``split`` separately sums positive and negative numbers, and then adds these sums together.\n", "\n", "``monotonic``, ``ordered`` and ``split`` are implemented with saturation arithmetic. \n", "``checked``, ``monotonic`` and ``split`` protect against underestimating sensitivity by preserving associativity.\n", "\n", "These each have their own uses.\n", "For example, if the dataset is considered ordered-- that is, the dataset distance metric is sensitive to changes in row ordering (`InsertDeleteDistance`), then neighboring datasets share the same row ordering, and it becomes safe to use arithmetic that saturates at the minimum and maximum representable values of the data type." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# because of the way this constructor is called...\n", "make_bounded_sum(bounds=(1, 20), MI=\"InsertDeleteDistance\")\n", "# ...it internally uses this constructor to build the transformation:\n", "make_bounded_int_ordered_sum(bounds=(1, 20))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "If you are trying to model a computation performed outside of OpenDP, you may not have access to saturation arithmetic.\n", "If this is the case, you can use `make_sized_bounded_int_checked_sum` to perform an overflow check at the moment the constructor is called." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Configure the type argument \"T\" to determine the necessary bit depth to avoid overflow\n", "make_sized_bounded_int_checked_sum(1234, (-2, 4), T=\"i32\")" ] } ], "metadata": { "kernelspec": { "display_name": "psi", "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.13" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "3220da548452ac41acb293d0d6efded0f046fab635503eb911c05f743e930f34" } } }, "nbformat": 4, "nbformat_minor": 2 }