{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Aggregation: Quantile\n", "\n", "This notebook explains the b-ary tree technique for releasing quantiles.\n", "Examples in this notebook will a dataset sampled from $Exponential(\\lambda=20)$." ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHpCAYAAABN+X+UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZAUlEQVR4nO3de1zUdb4/8NdcmBmuw01AFASVwguKqSBqmUc2TNtyszKPpXlcq10ziy5e1rTf7mnRWssubq57um6a5a5Z65otkWmt5A1vmPcLIDBcBGa4DszM9/fHMF8YAeX+ncvr+XjMw/zOd2beY8XLz10mCIIAIiIickhyqQsgIiKitjGoiYiIHBiDmoiIyIExqImIiBwYg5qIiMiBMaiJiIgcGIOaiIjIgTGoO0kQBBgMBnAZOhER9SQGdSdVVlZCq9WisrJS6lKIiMiFMaiJiIgcGIOaiIjIgTGoiYiIHBiDmoiIyIExqImIiBwYg5qIiMiBMaiJiIgcGIOaiIjIgTGoiYiIHBiDmoiIyIExqImIiBwYg5qIiMiBMaiJiIgcGIOaiIjIgTGoiYiIHBiDmoiIyIExqImIiBwYg5qIiMiBMaiJiIgcGIOauuxYXgUWfnwY+86VSF0KEZHLYVBTlxQZ6vDrjw4h/ecizP/wEDYfyJG6JCIil8Kgpk5rMFuwaHMWSqvq4aNWwmwR8LsvsrHnTLHUpRERuQwGNXXaV8cKcDinHL4aJf65eCJmjYkAAGw+kCtxZUREroNBTZ32feOY9LykKEQHe2PhHdHW62eLca3KKGVpREQug0FNnWKxCPjPhVIAwB239AEADA7xxcj+WpgsAr48ViBleURELoNBTZ3yc6EBZdX18FYpMCrSX7w+c3R/AMA/sq5KVBkRkWthUFOn/NjYmh43MAgeiqb/jH45IhweChlOFRhwobhKqvKIiFwGg5o65Yfz1vHpiTHBdtcDvFVIiA4EAGReutbrdRERuRoGNXVYXYMZh66UAwBuvy6oASAhKggAcOhyWa/WRUTkihjU1GEn8/WoN1kQ4qvGoD4+LZ4fGx0AADh4uQyCIPR2eURELoVBTR12vsg69hzb1w8ymazF86MiAuChkEFnqMPV8treLo+IyKUwqKnDbJPEBrfSmgYAT5UCw/tpAVhb1URE1HkMauqwCyWNQR3SelADQEKUdULZoSsMaiKirmBQU4ddLL55UI9tDOqDDGoioi5hUFOHVBtNyK+wjju3J6gvlVSjoqa+V2ojInJFDGrqkEsl1QCAQG8VAr1Vbd6n9fJAZKAXAOsuZkRE1DkMauqQiyU3nkjW3NC+fgCAnwsY1EREncWgpg6xzfgedINub5th4QxqIqKuYlBTh1xox0Qym6G2oGbXNxFRpzGoqUPaszTLxhbUF4qrUNdg7tG6iIhcFYOa2s1ktuBKqXUyWXuCOsxPgwAvD5gsAk/SIiLqJAY1tZvOUAeTRYCHQoa+fpqb3i+TycRW9akCfU+XR0TkkhwiqDds2ICoqChoNBokJibi4MGDN7x/27ZtiI2NhUajQVxcHHbt2mX3/Msvv4zY2Fh4e3sjICAAycnJOHDggN09ZWVlmDNnDvz8/ODv748FCxagqoqtvhsp1NcBAMK0GsjlLff4bg1nfhMRdY3kQf3ZZ58hNTUVq1evRlZWFkaOHImUlBQUFxe3ev/+/fsxe/ZsLFiwAEePHsWMGTMwY8YMZGdni/fccssteOedd3Dy5En8+OOPiIqKwl133YWSkhLxnjlz5uDUqVNIT0/Hzp07sW/fPjz++OM9/n2dWUHjRid9tZ7tfg0nlBERdY1MkPgcwsTERIwdOxbvvPMOAMBisSAiIgKLFy/GsmXLWtw/a9YsVFdXY+fOneK1cePGIT4+Hhs3bmz1MwwGA7RaLb799ltMmTIFp0+fxtChQ3Ho0CGMGTMGALB7925MmzYNV69eRXh4eIv3MBqNMBqNdu8ZEREBvV4PPz+/Lv0ZOIuNey9izddnMCM+HOsfHtWu15zRGTB1/Q/wVStx4uW7Wj1ti4iI2iZpi7q+vh5HjhxBcnKyeE0ulyM5ORmZmZmtviYzM9PufgBISUlp8/76+nps2rQJWq0WI0eOFN/D399fDGkASE5Ohlwub9FFbpOWlgatVis+IiIiOvRdXUGhrUXt3/4WdXSwNxRyGSqNJhRXGm/+AiIisiNpUJeWlsJsNiM0NNTuemhoKHQ6Xauv0el07bp/586d8PHxgUajwRtvvIH09HQEBweL7xESEmJ3v1KpRGBgYJufu3z5cuj1evGRl5fXoe/qCgoax6jDtTefSGajViowIMi6lajtHGsiImo/yceoe8rkyZNx7Ngx7N+/H1OnTsVDDz3U5rh3e6jVavj5+dk93I2uMag7MkYNADGNS7nOF1d2e01ERK5O0qAODg6GQqFAUVGR3fWioiKEhYW1+pqwsLB23e/t7Y3Bgwdj3LhxeO+996BUKvHee++J73F9aJtMJpSVlbX5uQQU6m1d3+1vUQNATIgvAOA811ITEXWYpEGtUqkwevRoZGRkiNcsFgsyMjKQlJTU6muSkpLs7geA9PT0Nu9v/r62yWBJSUmoqKjAkSNHxOe/++47WCwWJCYmdvbruDSjyYzSKutxleEdbVGHWlvUF9j1TUTUYUqpC0hNTcW8efMwZswYJCQkYP369aiursb8+fMBAHPnzkW/fv2QlpYGAFiyZAkmTZqEdevWYfr06di6dSsOHz6MTZs2AQCqq6vxyiuv4N5770Xfvn1RWlqKDRs2ID8/Hw8++CAAYMiQIZg6dSoWLlyIjRs3oqGhAU899RQefvjhVmd8U1O3t8ZDDn8vjw691raL2bniSgiCwJnfREQdIHlQz5o1CyUlJVi1ahV0Oh3i4+Oxe/duccJYbm4u5PKmhv/48eOxZcsWrFy5EitWrEBMTAx27NiB4cOHAwAUCgXOnDmDjz76CKWlpQgKCsLYsWPxww8/YNiwYeL7bN68GU899RSmTJkCuVyOmTNn4q233urdL+9ECipsE8k8Oxy0g/r4QCYDKmoacK26HsE+6p4okYjIJUm+jtpZ2dZmu8s66u1ZV5H6+XFMGByEzb8e1+HXT3ptD3Ku1eDTheOQNCioByokInJNLjvrm7pXYSdnfNvYZn5f4MxvIqIOYVBTu9i2D+3IGurmBnPmNxFRpzCoqV3EFnUHdiVrbrDYomZQExF1BIOa2qXpQI7OtagH9vEGAFxuPM+aiIjah0FN7VLSuE93WGeDOtga1IX6OtTUm7qtLiIiV8egppsymS0oq7FudtLZpVX+XioEeqsAAFdKa7qtNiIiV8egppsqq6mHIAByGRDgper0+0QHs/ubiKijGNR0U9catw4N9FZBIe/8rmJNQc0JZURE7cWgppsqrbKOTwd5d21HMVtQX2KLmoio3RjUdFO2oA727Xy3N9A0oYxd30RE7cegppuydX13dY/uaC7RIiLqMAY13VRJN3V9RwVZg7qipgHl1fVdrouIyB0wqOmmSisbW9Rd7PrWeCjQr3Fns0ucUEZE1C4Marqpa9WNY9RdbFEDzSaUlbD7m4ioPRjUdFPdNZkM4FpqIqKOYlDTTYld312cTAYwqImIOopBTTckCILY9R3UHUHNmd9ERB3CoKYbMtSa0GAWAABB3l3v+m6+ltpiEbr8fkREro5BTTdkW5rlq1FC46Ho8vv1D/CCh0IGo8mCQkNdl9+PiMjVMajphq7ZJpJ1Q7c3ACjkMgxoXE99mTO/iYhuikFNN1Qq7krW9W5vGx7OQUTUfgxquqHuOpCjuYE8nIOIqN0Y1HRD17pxDbUNl2gREbUfg5puqKSbDuRojruTERG1H4OabsjWou6ONdQ2trXUV8trYDSZu+19iYhcEYOabqi8xtqiDvTqvq7vPj5q+KiVsAhAXllNt70vEZErYlDTDZXXNAAAArw8uu09ZTKZ2P19kd3fREQ3xKCmG6pobFH7d2OLGmgap77CCWVERDfEoKY2CYKACluL2rv7WtQAEGUL6msMaiKiG2FQU5sqjSaYGvfjDujmFvVAzvwmImoXBjW1qaLa2prWeMi7ZZ/v5tiiJiJqHwY1tck247u7W9MAEN2433eRwYhqo6nb35+IyFUwqKlN5T00kQwAtF4eCGw8NpM7lBERtY1BTW2q6IGlWc1Fs/ubiOimGNTUpp7s+gaa7fnNCWVERG1iUFObbC1q/x5uUV9mi5qIqE0MampTRW+1qDlGTUTUJgY1tam8h1vUUUHcnYyI6GYY1NSmnpz1DQBRwV6Nn9Mgtt6JiMgeg5ra1NOzvr1USoT5aQCw+5uIqC0MampTT7eoAY5TExHdDIOa2tTTLWqgaStRBjURUesY1NSqepMFVY1be/bUrG+g6XAOBjURUesY1NSqilprt7dMBvh5skVNRCQVBjW1ytbtrfX0gEIu67HPEbcRLa2GIAg99jlERM6KQU2tKq/u2c1ObCIDvSCXAdX1ZpRUGnv0s4iInBGDmlrV05ud2KiUcvQPsK6nZvc3EVFLDGpqlb62d1rUAMepiYhuhEFNreqtFjXQbOY3D+cgImqBQU2tEjc78eyFFnVQY9c3j7skImqBQU2tMtT2Xos6uo8PAHZ9ExG1hkFNrdI3BrWfRtnjnxXdeIpWTlkNzBYu0SIiao5BTa0y1Fp3JdP2Qou6X4AnPBQy1JssKKio7fHPIyJyJgxqapWtRa3twV3JbBRyGSIDrePUVzihjIjIDoOaWtWbQQ0A0cEcpyYiag2DmlplqLONUfdWUHPTEyKi1jhEUG/YsAFRUVHQaDRITEzEwYMHb3j/tm3bEBsbC41Gg7i4OOzatUt8rqGhAUuXLkVcXBy8vb0RHh6OuXPnoqCgwO49oqKiIJPJ7B5r1qzpke/nbCwWQZz1zRY1EZG0JA/qzz77DKmpqVi9ejWysrIwcuRIpKSkoLi4uNX79+/fj9mzZ2PBggU4evQoZsyYgRkzZiA7OxsAUFNTg6ysLLz00kvIysrC9u3bcfbsWdx7770t3uv3v/89CgsLxcfixYt79Ls6i6p6E2yTr3vy5Kzmohpb1FcY1EREdmSCxEcWJSYmYuzYsXjnnXcAABaLBREREVi8eDGWLVvW4v5Zs2ahuroaO3fuFK+NGzcO8fHx2LhxY6ufcejQISQkJCAnJweRkZEArC3qZ555Bs8880y76jQajTAamw6NMBgMiIiIgF6vh5+fX3u/rlO4Wl6DiWv3QKWU49z/3t0rn6nT12FcWgYUchlO/34qVErJ/w5JROQQJP1pWF9fjyNHjiA5OVm8JpfLkZycjMzMzFZfk5mZaXc/AKSkpLR5PwDo9XrIZDL4+/vbXV+zZg2CgoIwatQovPbaazCZTG2+R1paGrRarfiIiIhoxzd0Tr09kQwAQv3U8PRQwGwRkFde02ufS0Tk6CQN6tLSUpjNZoSGhtpdDw0NhU6na/U1Op2uQ/fX1dVh6dKlmD17tl3L9+mnn8bWrVuxZ88ePPHEE/jjH/+IF198sc1aly9fDr1eLz7y8vLa+zWdjhRBLZPJxMM52P1NRNSk57edklBDQwMeeughCIKAd9991+651NRU8Z9HjBgBlUqFJ554AmlpaVCr1S3eS61Wt3rdFfX2RDKbgcHeOF1o4IQyIqJmJG1RBwcHQ6FQoKioyO56UVERwsLCWn1NWFhYu+63hXROTg7S09NvOo6cmJgIk8mEK1eudPyLuBjbrmS9sX1oc1FcokVE1IKkQa1SqTB69GhkZGSI1ywWCzIyMpCUlNTqa5KSkuzuB4D09HS7+20hff78eXz77bcICgq6aS3Hjh2DXC5HSEhIJ7+N65Ci6xvgEi0iotZI3vWdmpqKefPmYcyYMUhISMD69etRXV2N+fPnAwDmzp2Lfv36IS0tDQCwZMkSTJo0CevWrcP06dOxdetWHD58GJs2bQJgDekHHngAWVlZ2LlzJ8xmszh+HRgYCJVKhczMTBw4cACTJ0+Gr68vMjMz8eyzz+KRRx5BQECANH8QDkS6oGaLmojoepIH9axZs1BSUoJVq1ZBp9MhPj4eu3fvFieM5ebmQi5vaviPHz8eW7ZswcqVK7FixQrExMRgx44dGD58OAAgPz8fX331FQAgPj7e7rP27NmDO++8E2q1Glu3bsXLL78Mo9GI6OhoPPvss3bj1u5M3JWs18eorS3qQn0dqo0meKsl/8+TiEhykq+jdlYGgwFardYl11Ev2XoUXx4rwMrpQ/Dr2wf26meP/kM6rlXXY+fiiRjeT9urn01E5Ii4qwS1IJ5F3cstagAY1Mfaqr5YUtXrn01E5IgY1NSCVGPUADAoxLqW+mIxg5qICGBQUyts66h76+Ss5ppa1JxQRkQEMKipFfrGddSStKjZ9U1EZIdBTXYEodkRl17SBfWl0mqYLZznSETEoCY7RpMF9WYLgN7fmQwA+gV4QqWUo95kQX55ba9/PhGRo2FQkx3bRDKFXAYfCdYxK+QyDGw8nIPd30REDGq6jrg0S6OETCaTpAaOUxMRNWFQkx2pTs5qblAftqiJiGwY1GRHys1ObAaFNLaoi7lEi4iIQU12pNzsxIZd30RETRjUZKeyznYWtXRBHd04mexadT3Kq+slq4OIyBEwqMmObYzaV4KlWTbeaiXCtRoAwKVStqqJyL0xqMlOpdHaopYyqAGOUxMR2TCoyU5lnXT7fDfHcWoiIisGNdkx1DpIi5pLtIiIADCo6TqGOtsYtaO0qNn1TUTujUFNdsRZ3xIuzwKaxqhzy2pgNJklrYWISEoMarLT1KKWtus7xFcNH7USZouA3Gs1ktZCRCQlBjXZsbWopQ5qmUzGcWoiIjCo6TqOMusbaBqnvlDMoCYi98WgJlG9yYK6BttZ1A4Q1CEMaiIiBjWJbK1pAPCRuOsbAGIag/pcEYOaiNwXg5pEtvFpH7USCrk0Z1E3d0uoLwDrGLXZIkhcDRGRNBjUJHKUGd82EYFeUCvlMJosyCvjzG8ick8MahI5yoxvG4VcJk4oO1dUKXE1RETSYFCTyJFmfNvcEmoN6vOcUEZEbopBTSJH2ee7uZjGcerzbFETkZtiUJPIUfb5bo4zv4nI3TGoSdS0z7fjtKg585uI3B2DmkSO2KLmzG8icncMahI52qxvgDO/iYgY1CRyxFnfAGd+E5F7Y1CTyBFnfQOc+U1E7o1BTaJKo2O2qDnzm4jcGYOaRI446xvgzG8icm8MahIZah1v1jfAmd9E5N4Y1AQAEATBIWd9A5z5TUTujUFNAIC6BgtMjd3KjjZGDXDmNxG5LwY1AWhamqWQy+ClUkhcTUuc+U1E7opBTQCadiXzUSshk8kkrqYlzvwmInfFoCYAgMFBx6dtOPObiNwVg5oAAFViUDve+DRgP/M7lzO/iciNMKgJgGPu892cQi7D4Mbu77M6jlMTkftgUBOApslkvmrHDGoAuDXM2v3NoCYid8KgJgBAldGxW9QAEGsL6iKDxJUQEfUeBjUBaJpM5uPAQX1rmB8A4Axb1ETkRhjUBMDxJ5MBwJDGFvWV0mrUNZglroaIqHcwqAlA0xi1jwOPUffxVSPAywMWATjP9dRE5CYY1ASgaYzaz4G7vmUymTih7IyO49RE5B4Y1ASg+fIsx+36BoDYxnFqzvwmInfBoCYAztH1DTRbosU9v4nITTCoCQBQ6QTLswA06/pmUBORe2BQE4Cmrm9HXp4FALc27vldUmnEtSqjxNUQEfU8BjUBaFqe5YhnUTfnrVYiMtALAMepicg9OERQb9iwAVFRUdBoNEhMTMTBgwdveP+2bdsQGxsLjUaDuLg47Nq1S3yuoaEBS5cuRVxcHLy9vREeHo65c+eioKDA7j3KysowZ84c+Pn5wd/fHwsWLEBVlXsu+WkwW1DbuC7Z0ceoAXZ/E5F7kTyoP/vsM6SmpmL16tXIysrCyJEjkZKSguLi4lbv379/P2bPno0FCxbg6NGjmDFjBmbMmIHs7GwAQE1NDbKysvDSSy8hKysL27dvx9mzZ3Hvvffavc+cOXNw6tQppKenY+fOndi3bx8ef/zxHv++jqi6cXwacPyub6DZVqIMaiJyAzJBECQ93DcxMRFjx47FO++8AwCwWCyIiIjA4sWLsWzZshb3z5o1C9XV1di5c6d4bdy4cYiPj8fGjRtb/YxDhw4hISEBOTk5iIyMxOnTpzF06FAcOnQIY8aMAQDs3r0b06ZNw9WrVxEeHn7Tug0GA7RaLfR6Pfz8/Drz1R1GXlkNbn91DzQecpz5w91Sl3NTO08U4KktRzEywh9fLpogdTlERD1K0hZ1fX09jhw5guTkZPGaXC5HcnIyMjMzW31NZmam3f0AkJKS0ub9AKDX6yGTyeDv7y++h7+/vxjSAJCcnAy5XI4DBw60+h5GoxEGg8Hu4SoMtpOzHHx82sa2lvp8USUsFkn/nklE1OMkDerS0lKYzWaEhobaXQ8NDYVOp2v1NTqdrkP319XVYenSpZg9e7bY8tXpdAgJCbG7T6lUIjAwsM33SUtLg1arFR8RERHt+o7OQNzn2wnGpwEgKsgLKqUcNfVm5JXXSF0OEVGPknyMuic1NDTgoYcegiAIePfdd7v0XsuXL4derxcfeXl53VSl9Jp2JXOOoFYq5IgJ8QEAnC7kODURuTZJgzo4OBgKhQJFRUV214uKihAWFtbqa8LCwtp1vy2kc3JykJ6ebjeOHBYW1mKymslkQllZWZufq1ar4efnZ/dwFU1nUTtH1zfQbIcyTigjIhcnaVCrVCqMHj0aGRkZ4jWLxYKMjAwkJSW1+pqkpCS7+wEgPT3d7n5bSJ8/fx7ffvstgoKCWrxHRUUFjhw5Il777rvvYLFYkJiY2B1fzak4y/ahzYkzv4tcZ64AEVFrJP/JnJqainnz5mHMmDFISEjA+vXrUV1djfnz5wMA5s6di379+iEtLQ0AsGTJEkyaNAnr1q3D9OnTsXXrVhw+fBibNm0CYA3pBx54AFlZWdi5cyfMZrM47hwYGAiVSoUhQ4Zg6tSpWLhwITZu3IiGhgY89dRTePjhh9s149vVOMv2oc3ZJpSdYdc3Ebk4yX8yz5o1CyUlJVi1ahV0Oh3i4+Oxe/duccJYbm4u5PKmhv/48eOxZcsWrFy5EitWrEBMTAx27NiB4cOHAwDy8/Px1VdfAQDi4+PtPmvPnj248847AQCbN2/GU089hSlTpkAul2PmzJl46623ev4LOyBn2T60udi+1hb15WvVqKk3wUvlPLUTEXWE5OuonZUrraNeueMkPvkpF09PiUHqL26Rupx2G/O/36K0yogvfjseoyIDpC6HiKhHuPSsb2ofZ1ueZTM03PoXpJ8LOU5NRK6LQU1OtzzLZkhj9/dpBjURuTAGNYmTyZxpjBoAhvZtbFEXMKiJyHUxqKlZi9p51lEDTUF9RsetRInIdTGoCVVG217fztWijg72hrpxK9GcMm4lSkSuiUFNTS1qJ5tMplTIxY1P2P1NRK6KQe3mBEFomvXtZF3fADCksfubE8qIyFUxqN1cXYMFpsbxXWebTAZwiRYRuT4GtZuz7fMtkwHeKoXE1XQcZ34TkatjULs5cWmWWgmZTCZxNR0X2xjUOkMdyqrrJa6GiKj7MajdnG0imZ8Tjk8D1r9gDAjyAsBxaiJyTQxqN2ebSOZMR1xeb0gYJ5QRketiULs52xi1s62hbk6cUMZxaiJyQQxqN+eMZ1FfT5xQxhY1EbmgTgX1pUuXursOkkjTWdTOOUYNAEMaW9QXiqtgNJklroaIqHt1KqgHDx6MyZMn45NPPkFdXV1310S9yBW6vsO1Gmg9PWCyCDhfVCV1OURE3apTQZ2VlYURI0YgNTUVYWFheOKJJ3Dw4MHuro16gbOeRd2cTCZj9zcRuaxOBXV8fDzefPNNFBQU4P3330dhYSEmTpyI4cOH4/XXX0dJSUl310k9xFnPor7eEG58QkQuqkuTyZRKJe6//35s27YNa9euxYULF/D8888jIiICc+fORWFhYXfVST2kyuj8y7MAYHg/a1Bn5+slroSIqHt1KagPHz6M3/72t+jbty9ef/11PP/887h48SLS09NRUFCA++67r7vqpB5iEMeonXcyGQDE9dMCsHZ9m3k2NRG5kE41o15//XV88MEHOHv2LKZNm4aPP/4Y06ZNg1xuzf3o6Gh8+OGHiIqK6s5aqQeILWon7/oe2McHnh4K1NSbcbm0CoNDfKUuiYioW3Tqp/O7776L//mf/8Fjjz2Gvn37tnpPSEgI3nvvvS4VRz3PVcaoFXIZhob74UhOObLzDQxqInIZner6Tk9Px9KlS1uEtCAIyM3NBQCoVCrMmzev6xVSj7Itz3LWvb6bG964nvokx6mJyIV0KqgHDRqE0tLSFtfLysoQHR3d5aKo97jCXt82wxvHqTmhjIhcSaeCWhBan6xTVVUFjUbTpYKo95gtAqrrrTt5OXvXN9AU1KcKDLBwQhkRuYgO/XROTU0FYN1gYtWqVfDy8hKfM5vNOHDgAOLj47u1QOo5tolkgPNPJgOAmBAfqJVyVBlNyCmrQXSwt9QlERF1WYd+Oh89ehSAtUV98uRJqFQq8TmVSoWRI0fi+eef794KqcfYxqdVSjnUSoXE1XSdUiFHbF8/HM+rwMl8PYOaiFxCh4J6z549AID58+fjzTffhJ+fX48URb3D1qJ25u1DrxfXzxrUp/L1uHdkuNTlEBF1Wad+Qn/wwQfdXQdJwFWWZjU3PNw6Ts2Z30TkKtr9E/r+++/Hhx9+CD8/P9x///03vHf79u1dLox6nq3r2xXGp22az/wWBAEymUziioiIuqbdP6G1Wq34Q0+r1fZYQdR7xBa12vnXUNvcEuoLlUIOQ50JeWW1iAzyuvmLiIgcWLuDunl3N7u+XYMtqF2pRa1SynFrmC9O5uuRXaBnUBOR0+vUOura2lrU1NSIv8/JycH69evx73//u9sKo54nTiZzoaAGmk7S4jg1EbmCTgX1fffdh48//hgAUFFRgYSEBKxbtw733Xcf3n333W4tkHqOK20f2hx3KCMiV9KpoM7KysLtt98OAPj73/+OsLAw5OTk4OOPP8Zbb73VrQVSz3Gl7UObs838tk0oIyJyZp0K6pqaGvj6Wk8n+ve//437778fcrkc48aNQ05OTrcWSD3HFZdnAcCtYb5QymUor2lAfkWt1OUQEXVJp4J68ODB2LFjB/Ly8vDNN9/grrvuAgAUFxdzExQnYnDByWQAoPFQ4NYw618kT1xl9zcRObdOBfWqVavw/PPPIyoqComJiUhKSgJgbV2PGjWqWwuknlNltI5R+7rYGDUAjIzwBwAcz6uQtA4ioq7qVFA/8MADyM3NxeHDh7F7927x+pQpU/DGG290W3HUs5rWUbtWixoA4huD+iiDmoicXKd/QoeFhSEsLMzuWkJCQpcLot7jqsuzgKagPnlVD5PZAqWiU38nJSKSXKd+QldXV2PNmjXIyMhAcXExLBaL3fOXLl3qluKoZzVNJnO9ru9BfXzgrVKgut6MCyVViA3j3Akick6dCupf//rX2Lt3Lx599FH07duX+yk7IUEQYKi1jVG7XotaIZdhRH9/ZF66hmO5FQxqInJanfoJ/fXXX+Nf//oXJkyY0N31UC+pa7DAZLGuMXbFoAasE8oyL13D8asVeDghUupyiIg6pVMDdwEBAQgMDOzuWqgX2XYlk8kAb5VrBrU4oSy3QtI6iIi6olNB/Yc//AGrVq2y2++bnIuh2a5kcrlrDl3YgvpcUSVq6k3SFkNE1EmdakqtW7cOFy9eRGhoKKKiouDhYT8ZKSsrq1uKo57jqvt8Nxem1SDMTwOdoQ4nr+qRODBI6pKIiDqsU0E9Y8aMbi6Depurbh96vZERWuhO1eH41QoGNRE5pU79lF69enV310G9zBbUrtyiBoD4iAB8c6oIx/O4lSgROadO7wJRUVGB//u//8Py5ctRVlYGwNrlnZ+f323FUc+xdX27Q4saAI5xhzIiclKd+il94sQJJCcnQ6vV4sqVK1i4cCECAwOxfft25ObmimdVk+MyuElQj+jvD5kMyK+oRXFlHUJ8NVKXRETUIZ1qUaempuKxxx7D+fPnodE0/eCbNm0a9u3b123FUc9x5V3JmvNRKxET4gMA7P4mIqfUqaA+dOgQnnjiiRbX+/XrB51O1+WiqOe5y2QyoGmZFk/SIiJn1KmgVqvVMBgMLa6fO3cOffr06XJR1POaur5du0UNNB15yXFqInJGnQrqe++9F7///e/R0GDb3UqG3NxcLF26FDNnzuzWAqlniLO+PV2/RT0qIgCANajNjdumEhE5i04F9bp161BVVYU+ffqgtrYWkyZNwuDBg+Hr64tXXnmlu2ukHlDpRi3qW8N84aNWospowlldpdTlEBF1SKeaU1qtFunp6fjPf/6D48ePo6qqCrfddhuSk5O7uz7qIe40Rq2QyzAq0h8/nC/F4ZwyDA3nSVpE5Dw63KK2WCx4//33cc899+CJJ57Au+++ix9//BEFBQUQhI53K27YsAFRUVHQaDRITEzEwYMHb3j/tm3bEBsbC41Gg7i4OOzatcvu+e3bt+Ouu+5CUFAQZDIZjh071uI97rzzTshkMrvHk08+2eHanVnThieuH9QAMDbKeojMoSvlEldCRNQxHQpqQRBw77334te//jXy8/MRFxeHYcOGIScnB4899hh+9atfdejDP/vsM6SmpmL16tXIysrCyJEjkZKSguLi4lbv379/P2bPno0FCxbg6NGjmDFjBmbMmIHs7GzxnurqakycOBFr16694WcvXLgQhYWF4uPVV1/tUO3Ozp0mkwHAmCjrOPWRK2USV0JE1EFCB7z//vuCr6+v8N1337V4LiMjQ/D19RU++uijdr9fQkKCsGjRIvH3ZrNZCA8PF9LS0lq9/6GHHhKmT59udy0xMVF44oknWtx7+fJlAYBw9OjRFs9NmjRJWLJkSbvrFARBqKurE/R6vfjIy8sTAAh6vb5D7+MILBaLMHD5v4QBS3cKBRU1UpfTK6qNDeJ3vlruHt+ZiFxDh1rUn376KVasWIHJkye3eO6//uu/sGzZMmzevLld71VfX48jR47YjWvL5XIkJycjMzOz1ddkZma2GAdPSUlp8/4b2bx5M4KDgzF8+HAsX778pkd2pqWlQavVio+IiIgOf6ajqG0wi7Of3aVF7aVSYnjj2PRhtqqJyIl0KKhPnDiBqVOntvn83XffjePHj7frvUpLS2E2mxEaGmp3PTQ0tM1NU3Q6XYfub8t///d/45NPPsGePXuwfPly/O1vf8Mjjzxyw9csX74cer1efOTl5XXoMx2JbXxaLgO8VQqJq+k9Y8RxagY1ETmPDs0kKisraxGUzYWGhqK83PEn6zz++OPiP8fFxaFv376YMmUKLl68iEGDBrX6GrVaDbVa3Vsl9qjmS7NkMpnE1fSeMQMC8N6Pl3GYE8qIyIl0qEVtNpuhVLad7QqFAiaTqV3vFRwcDIVCgaKiIrvrRUVFCAsLa/U1YWFhHbq/vRITEwEAFy5c6NL7OAuDGy3Nam5044Sys0WV0Nc2SFwNEVH7dOgntSAIeOyxx9psWRqNxna/l0qlwujRo5GRkYEZM2YAsC79ysjIwFNPPdXqa5KSkpCRkYFnnnlGvJaeno6kpKR2f25rbEu4+vbt26X3cRbuciDH9UJ8NYgK8sKVazXIyi3H5FtDpC6JiOimOhTU8+bNu+k9c+fObff7paamYt68eRgzZgwSEhKwfv16VFdXY/78+eJ79evXD2lpaQCAJUuWYNKkSVi3bh2mT5+OrVu34vDhw9i0aZP4nmVlZcjNzUVBQQEA4OzZswCsrfGwsDBcvHgRW7ZswbRp0xAUFIQTJ07g2WefxR133IERI0a0u3ZnZqh1jyMuWzN6QCCuXKvBkSsMaiJyDh36Sf3BBx9064fPmjULJSUlWLVqFXQ6HeLj47F7925xHDw3NxdyeVPv/Pjx47FlyxasXLkSK1asQExMDHbs2IHhw4eL93z11Vdi0APAww8/DABYvXo1Xn75ZahUKnz77bfiXwoiIiIwc+ZMrFy5slu/myNzt81OmhsbFYB/ZF3lhDIichoyQejEdmIEg8EArVYLvV4PPz/n2pLyL3svIu3rM/jVqH54Y1a81OX0qgvFVUh+fS/USjlOvpwClbJT290TEfUa/pRyQ+7coh7UxxsBXh4wmizILtBLXQ4R0U0xqN2QO52cdT2ZTIbRA6zrqbnxCRE5Awa1G3Knk7NakxhtDeqfLjGoicjxMajdkMFNl2fZJA0KAgAcuHQNDWaLxNUQEd0Yg9oNNXV9u2eLemhfP2g9PVBdb8bJfI5TE5FjY1C7IXfdmcxGLpchaaC1VZ158ZrE1RAR3RiD2g2582Qym/GDrUG9/2KpxJUQEd0Yg9oN2SaTaT3ds0UNAOMbx6kPXylHXYNZ4mqIiNrGoHYzgiCgyujek8kAYFAfH4T4qmE0WZCVy9O0iMhxMajdTE29GWaLdTM6dx2jBqzrqW2tao5TE5EjY1C7GVu3t0Iug6eHQuJqpDV+UDAAYD+DmogcGIPazTRfmiWTySSuRlq29dTH8yrE4QAiIkfDoHYz7r40q7mIQC9EBHrCZBF4mhYROSwGtZsx2FrUavedSNbc+IHW7m+OUxORo2JQuxnx5Cw3XprVHNdTE5GjY1C7GW52Ys+2Q9mpAgPKq+slroaIqCUGtZtx95Ozrhfip0FMiA8EAfgPW9VE5IAY1G7G1qL2Y4taNOmWPgCAvWdLJK6EiKglBrWbYYu6pTtvDQEA7D1XAkEQJK6GiMgeg9rNMKhbGhsdAE8PBYorjThdWCl1OUREdhjUbsZQy8lk11MrFeJ2ot+fK5a4GiIiewxqNyMuz2JQ25l0K8epicgxMajdjKHZFqLU5M5brOPUR3LKxQl3RESOgEHtZjhG3brIIC9EB3vDZBHwnwtcpkVEjoNB7Wa44UnbxGVa59j9TUSOg0HtRgRBEE+J8mOLuoU7G8epvz/LZVpE5DgY1G6kut4MS2P+sEXd0riBQVAr5SjU1+F8cZXU5RARAWBQuxVbt7eHQgaNB//VX0/jocC4xr2/vz/LZVpE5Bj409qNGGptE8k8IJPJJK7GMdnGqb/nMi0ichAMajdSyaVZNzU51rpM69CVMnEpGxGRlBjUboRLs24uOtgbA/t4o8EsYB9nfxORA2BQuxFxsxM1J5LdyC+GhAIAvv25SOJKiIgY1G6FLer2SR5qDeo9Z0tgMlskroaI3B2D2o00BTVb1DdyW2QAArw8oK9twOGccqnLISI3x6B2I5xM1j4KuUycVMbubyKSGoPajegbj7jUerJFfTPiOPXpIu5SRkSSYlC7EQZ1+91+Sx+oFHJcuVaDiyXVUpdDRG6MQe1GDLazqBnUN+WjVmLcIOsuZd+eZvc3EUmHQe1G2KLumF8MsY5TZzCoiUhCDGo3YmBQd8iUxnHqIznluFZllLgaInJXDGo3whZ1x4T7e2JYuB8sAvDdGR7SQUTSYFC7CUEQGNSdkNzYqk7nMi0ikgiD2k3U1JthbjyM2s+T66jb665h1qDed74EtfVmiashInfEoHYTtta0h0IGTw+FxNU4j6F9/RAR6Im6Bgv28pAOIpIAg9pNNO/25lnU7SeTyZAyNAwA8M0pncTVEJE7YlC7CVtQcw11x6UMtwZ1xukiNPCQDiLqZQxqN8GJZJ13W2QAgn1UMNSZ8NOla1KXQ0RuhkHtJmxrqP14claHKeQy/KKx+3t3Nru/iah3MajdBFvUXZMyrGmZlsXCQzqIqPcwqN0EdyXrmvGDguGrVqK40oijeRVSl0NEboRB7SbYou4alVKO/2rc+/vfnP1NRL2IQe0mmmZ9c7OTzkoZ1jhOfUrHM6qJqNcwqN2E7YhLtqg7b9ItfaBSypFzrQZniyqlLoeI3ASD2k2w67vrvNVK3BETDAD4Jpt7fxNR72BQuwlueNI9bN3f3KWMiHoLg9pNsEXdPZKHhEIhl+HnQgPyymqkLoeI3IDkQb1hwwZERUVBo9EgMTERBw8evOH927ZtQ2xsLDQaDeLi4rBr1y6757dv34677roLQUFBkMlkOHbsWIv3qKurw6JFixAUFAQfHx/MnDkTRUWu3ZWp54Yn3SLAW4WEqEAAbFUTUe+QNKg/++wzpKamYvXq1cjKysLIkSORkpKC4uLiVu/fv38/Zs+ejQULFuDo0aOYMWMGZsyYgezsbPGe6upqTJw4EWvXrm3zc5999ln885//xLZt27B3714UFBTg/vvv7/bv5yjqGsyoN1n3qNZ6Mai7aupwdn8TUe+RCRKuM0lMTMTYsWPxzjvvAAAsFgsiIiKwePFiLFu2rMX9s2bNQnV1NXbu3CleGzduHOLj47Fx40a7e69cuYLo6GgcPXoU8fHx4nW9Xo8+ffpgy5YteOCBBwAAZ86cwZAhQ5CZmYlx48a1WqvRaITRaBR/bzAYEBERAb1eDz8/v07/GfSGYkMdEv6YAbkMuPDKNMjlPD2rKwr1tUhK+w4yGXBwRTL6+KqlLomIXJhkLer6+nocOXIEycnJTcXI5UhOTkZmZmarr8nMzLS7HwBSUlLavL81R44cQUNDg937xMbGIjIy8obvk5aWBq1WKz4iIiLa/ZlSaz6RjCHddX21nhjZXwtBsG4pSkTUkyQL6tLSUpjNZoSGhtpdDw0NhU7XepeiTqfr0P1tvYdKpYK/v3+H3mf58uXQ6/XiIy8vr92fKTVOJOt+d3H2NxH1EsknkzkLtVoNPz8/u4ez4ESy7mcbp95/sRSGugaJqyEiVyZZUAcHB0OhULSYbV1UVISwsLBWXxMWFtah+9t6j/r6elRUVHTpfZyJLUjYou4+g/r4YHCIDxrMAvacaX3yIxFRd5AsqFUqFUaPHo2MjAzxmsViQUZGBpKSklp9TVJSkt39AJCent7m/a0ZPXo0PDw87N7n7NmzyM3N7dD7OJPy6sag5ozvbmU7+pLd30TUkyQ9oSE1NRXz5s3DmDFjkJCQgPXr16O6uhrz588HAMydOxf9+vVDWloaAGDJkiWYNGkS1q1bh+nTp2Pr1q04fPgwNm3aJL5nWVkZcnNzUVBQAMAawoC1JR0WFgatVosFCxYgNTUVgYGB8PPzw+LFi5GUlNTmjG9nV1FTDwAIYFB3q6nD+mLDnov4/mwJauvN8FQppC6JiFyQpEE9a9YslJSUYNWqVdDpdIiPj8fu3bvFCWO5ubmQy5sa/ePHj8eWLVuwcuVKrFixAjExMdixYweGDx8u3vPVV1+JQQ8ADz/8MABg9erVePnllwEAb7zxBuRyOWbOnAmj0YiUlBT8+c9/7oVvLI2KxjHqAC+VxJW4luH9/NDP3xP5FbXYe65EHLcmIupOkq6jdmYGgwFardYp1lEv/vQo/nm8AC/dMxQLJkZLXY5LeeVfP+OvP1zGjPhwrH94lNTlEJEL4qxvN2Dr+vbnZLJuN3V4XwDAt6eLYTSZJa6GiFwRg9oNVNQ0dn17M6i726gIf4T5aVBlNOHH86VSl0NELohB7QbKbS1qjlF3O7lcJo5Nf53N2d9E1P0Y1G7A1qJm13fPuLsxqP99SicefkJE1F0Y1C6uwWxBldEEgLO+e8qYqEAE+6hhqDMh89I1qcshIhfDoHZxtta0TGY9lIO6n0IuEzc/2Z1dKHE1RORqGNQuzjbj20/jAQVPzuox0+Kss7+/OVUEk5nd30TUfRjULq7cNuObu5L1qMToQAR4eaCsuh4HL5dJXQ4RuRAGtYur4IzvXqFUyHHXUM7+JqLux6B2ceKMb7aoe9zdcdag3n1KB7OFG/4RUfdgULu4cvFADraoe9r4QcHw1ShRUmnEkZxyqcshIhfBoHZxtgM52KLueSqlHL8Yap39/TVnfxNRN2FQu7imfb7Zou4Ndzfu/b07WwcLu7+JqBswqF1ceTX3+e5Nt8cEw1ulQKG+DseuVkhdDhG5AAa1i6uo5azv3qTxUGDKEGv3979OsPubiLqOQe3iuM937/vlyHAAwM4TBZz9TURdxqB2cZz13fvuuCUYfholigxGHLrCzU+IqGsY1C6O66h7n1qpEI++/OfxAomrISJnx6B2YbX1Zhgbj10M8GaLujfZur+/ztahgXt/E1EXMKhdmK3bWymXwVulkLga95I0MAhB3iqUVddj/0UefUlEncegdmHlzfb5lsl4clZvUirk4olaXx1j9zcRdR6D2oXZ1lAHcg21JO6Nt3Z///uUDnUNZomrISJnxaB2YdeqjQCAIG+1xJW4p9GRAeir1aDSaMLecyVSl0NETopB7cKuVVm7vgN9OJFMCnK5DPeMsHZ/c/Y3EXUWg9qF2VrUwZzxLRnb7O9vTxehymiSuBoickYMahdWVt3YombXt2Ti+mkxsI836hos2HWSW4oSUccxqF0Yu76lJ5PJMPO2/gCAvx+5KnE1ROSMGNQu7Fpji5pd39K6/7Z+kMmAg5fLkHOtWupyiMjJMKhdWFPXN4NaSn21npg4OBgA8I+sfImrISJnw6B2YaVVjcuzfDhGLbUHRlu7v7dnXYWFJ2oRUQcwqF1UvcmCyjrrLOMgtqgllzIsDL5qJa6W1+LAZZ6oRUTtx6B2UbbtQxVyGbQ8i1pyGg8F7mlcqsVJZUTUEQxqF2Xr9g7wUkEu5z7fjsDW/f11diGquaaaiNqJQe2ibBPJ2O3tOG6L9MfAYG/U1Ju5ppqI2o1B7aJsa6iDuIbaYchkMsxsbFV/fjhP4mqIyFkwqF3UNS7NckgPjO4PhVyGQ1fKcUZnkLocInICDGoXVWbb55tLsxxKqJ8Gdw0NBQBs/ilX4mqIyBkwqF2UuH0oW9QO55FxAwAAXxzN56QyIropBrWLYte340oaGIToYG9UGU344ih3KiOiG2NQuyjbrO9gTiZzOHK5TGxVf7j/CgSBO5URUdsY1C7qWuM6ah5x6ZgeGtMfPmolLhRXYd/5UqnLISIHxqB2Ubauby7Pcky+Gg88OMa6VOu9Hy9LXA0ROTIGtQsymszc59sJzB8fDZkM2HeuBOeLKqUuh4gcFIPaBZVUWru9VQo59/l2YJFBXuJSrb/suyRxNUTkqBjULqi4Maj7+Kohk3Gfb0f2mzsHAwB2HM3H1fIaiashIkfEoHZBxYamoCbHFh/hj4mDg2GyCPgrW9VE1AoGtQsqqawDAIQwqJ3Cb+8cBADYeigPxY3/7oiIbBjULsjW9R3qp5G4EmqPpEFBuC3SH0aTBX/ec1HqcojIwTCoXZCt65staucgk8nw3F23AgC2HMhFfkWtxBURkSNhULsgW/dpiB+D2llMGByMpIFBqDdb8HbGeanLISIHwqB2Qbau7xBfdn07k+dTrK3qbUeucl01EYkY1C6oiLO+ndLoAQFIGRYKs0XA//7rtNTlEJGDYFC7GJPZgmuNZ1Gz69v5LL97CDwUMuw9V4I9Z4ulLoeIHACD2sVcq66HIAByGRDEAzmcTlSwN+ZPiAYA/GHnzzCazBJXRERSY1C7mOabnSjk3JXMGT31X4MR7KPGpZJq/GUvN0EhcncOEdQbNmxAVFQUNBoNEhMTcfDgwRvev23bNsTGxkKj0SAuLg67du2ye14QBKxatQp9+/aFp6cnkpOTcf68/UzaqKgoyGQyu8eaNWu6/bv1NnHGNyeSOS0/jQdW/XIoAOCdPRdwubRa4oqISEqSB/Vnn32G1NRUrF69GllZWRg5ciRSUlJQXNz6+Nz+/fsxe/ZsLFiwAEePHsWMGTMwY8YMZGdni/e8+uqreOutt7Bx40YcOHAA3t7eSElJQV2d/a5Pv//971FYWCg+Fi9e3KPftTc0zfhmt7cz++WIvrg9Jhj1JguW/eMELBZB6pKISCKSB/Xrr7+OhQsXYv78+Rg6dCg2btwILy8vvP/++63e/+abb2Lq1Kl44YUXMGTIEPzhD3/AbbfdhnfeeQeAtTW9fv16rFy5Evfddx9GjBiBjz/+GAUFBdixY4fde/n6+iIsLEx8eHt79/TX7XHiZiecSObUZDIZXpkRBy+VAgcul+HD/VekLomIJCJpUNfX1+PIkSNITk4Wr8nlciQnJyMzM7PV12RmZtrdDwApKSni/ZcvX4ZOp7O7R6vVIjExscV7rlmzBkFBQRg1ahRee+01mEymNms1Go0wGAx2D0dU1Nj13Ydd304vMsgLK6YNAQCs3X0GF0uqJK6IiKQgaVCXlpbCbDYjNDTU7npoaCh0Ol2rr9HpdDe83/brzd7z6aefxtatW7Fnzx488cQT+OMf/4gXX3yxzVrT0tKg1WrFR0RERPu/aC/i9qGuZU5iJG6PCYbRZMGizVmoa+AscCJ3I3nXt1RSU1Nx5513YsSIEXjyySexbt06vP322zAaja3ev3z5cuj1evGRl5fXyxW3D0/Oci0ymQx/enAkgn1UOKOrxOovT0ldEhH1MkmDOjg4GAqFAkVFRXbXi4qKEBYW1uprwsLCbni/7deOvCcAJCYmwmQy4cqVK60+r1ar4efnZ/dwRLZdyXhylusI9dPgzYdHQSYDPjuch22HHfMviUTUMyQNapVKhdGjRyMjI0O8ZrFYkJGRgaSkpFZfk5SUZHc/AKSnp4v3R0dHIywszO4eg8GAAwcOtPmeAHDs2DHI5XKEhIR05StJqsFsEceow/09Ja6GutOEwcFITb4FAPDSl9k4o3PMORJE1P2UUheQmpqKefPmYcyYMUhISMD69etRXV2N+fPnAwDmzp2Lfv36IS0tDQCwZMkSTJo0CevWrcP06dOxdetWHD58GJs2bQJg7Sp85pln8L//+7+IiYlBdHQ0XnrpJYSHh2PGjBkArBPSDhw4gMmTJ8PX1xeZmZl49tln8cgjjyAgIECSP4fuoNPXQRAAlVKOYB+V1OVQN1s0eTAO55Rj77kS/PaTLHyxaAK0nh5Sl0VEPUzyoJ41axZKSkqwatUq6HQ6xMfHY/fu3eJksNzcXMjlTQ3/8ePHY8uWLVi5ciVWrFiBmJgY7NixA8OHDxfvefHFF1FdXY3HH38cFRUVmDhxInbv3g2NxtodrFarsXXrVrz88sswGo2Ijo7Gs88+i9TU1N798t3Mdo5xP39PyGTclczVyOUyvDErHve89QMulVZj0eYsfDB/LDwUbjvVhMgtyARB4E4KnWAwGKDVaqHX6x1mvHp71lWkfn4cEwYHYfOvx0ldDvWQUwV6PLgxEzX1Zjw8NgJp98fxL2ZELox/FXchBY0t6nAtx6dd2bBwLd6ePQpyGbD1UB7+so/7gRO5Mga1CxG7vgMY1K5uypBQvHSPdT/wNV+fwa6ThRJXREQ9hUHtQvIrOOPbncyfEI15SQMAAM9sPYb9F0slroiIegKD2oUUNJtMRu7hpXuG4q6hoag3W7Dwo8M4cbVC6pKIqJsxqF2EIAgMajekVMjx1uxRGD8oCNX1Zsx7/yAuFFdKXRYRdSMGtYuoqGlATb11H+gwLXclcycaDwU2zR2DEf21KK9pwKPvHRTnKxCR82NQuwjbD+ZgHzU0HgqJq6He5qNW4sP5CRjUxxuF+jo8+n8HUFrV+r71RORcGNQuooAzvt1eoLcKn/w6Ef38PXGptBrz3j8IQ12D1GURURcxqF1E065k7PZ2Z321nvjbggQEeatwqsCAue8dRCXDmsipMahdBDc7IZuBfXzwtwWJ8PfywLG8Csx7/yCqjCapyyKiTmJQuwhbi5prqAkAhob74ZMFidB6eiArtwKPMayJnBaD2kXkXKsBAEQGeklcCTmK4f20+GRBIvw0ShzOKcf8Dw6immFN5HQY1C5AEARcLq0GAET38Za4GnIkcf21+OTXifDVKHHoSjnmf3gINfUMayJnwqB2ASWVRtTUmyGXAREBbFGTvRH9/fHJgkT4qpU4eLkMj71/iLPBiZwIg9oFXGpsTUcEekGl5L9SamlkhD8+XpBgDesrZZj1l59QXFkndVlE1A78qe4CrjQGdVQQu72pbaMiA7D1iXEI9lHjdKEBD27MRG7j3AYiclwMahcgjk8HM6jpxoaFa/GP3yQhItATOddqMHPjfmTn66Uui4hugEHtAhjU1BEDgrzxjyfHIzbMFyWVRjywcT++PJYvdVlE1AYGtQtgUFNHhfhp8PmTSZh8ax/UNViwZOsx/HHXaZgtgtSlEdF1GNROzmwRkFNmHWdkUFNH+Gk88H/zxmLR5EEAgE37LuGxDw6i2MBJZkSOhEHt5AoqalFvskClkHNXMuowhVyGF1JiseG/b4OnhwI/nC/FL97Yhy+P5UMQ2LomcgQMaid35Zq12zsyyAsKuUziashZTR/RF189NQFx/bTQ1zZgydZj+M0nWTwqk8gBMKidHMenqbvEhPpi+2/HI/UXt0Apl2H3KR0mv/Y9/vz9BdQ1mKUuj8htMaid3IXiKgDAQAY1dQMPhRxPT4nBl09NwPB+fqg0mvDq7rP4rz99j38cucrJZkQSYFA7udOFBgDAkL5+EldCrmRYuBZfLZqI1x8aiXCtBgX6Ojy37TimrPsef/spB7X1bGET9RaZwBkjnWIwGKDVaqHX6+HnJ01IWiwCRvy/f6PKaMI3z9yBW8N8JamDXFtdgxkf7r+CjXsvoqLGukd4oLcKj4wbgEcSIxHip5G4QiLXxqDuJEcI6pxr1Zj02vdQKeU49f9S4KFgBwn1nJp6Ez4/lIf/+/EyrpZbzz9XymWYOjwM88ZHYcyAAMhknNBI1N2UUhdAnWfr9r4l1IchTT3OS6XEYxOi8ci4Adh9SoeP9l/BoSvl2HmiEDtPFGJIXz/MTRqA++LD4aXijxai7sL/m5zYzwXWoB7K8WnqRUqFHPeMCMc9I8JxqkCPv2XmYMexfJwuNGD59pNI23UaD46JwKPjBiCKkxyJuozNMCf2cyGDmqQ1LFyLNTNH4KflU/C7aUMQGegFQ50J7/14GXf+6Xs89sFBHLxcJnWZRE6NLWondrqwEgBnfJP0/L1UWHjHQCyYGI2950rwUeYVfH+2RHyMGxiIp6fEIGlgEMexiTqIQe2kKmrqkV9hndAzJJxBTY5BLpdhcmwIJseG4EppNTb9cAnbDufhp0tl+OnSASREBWLZtFjcFhkgdalEToNd307K1u3dP8ATfhoPiashaikq2Bt//FUc9r4wGfOSBkCllOPglTLc/+f9+O3mI+KuekR0YwxqJ3UsrwIAMDxcK20hRDcR7u+J/3ffcOx7YTJmjYmAXAbsOqnDL17fi9VfZnM/caKbYFA7KdsEnbHRgRJXQtQ+YVoN1j4wAruW3I7Jt/aBySLgo8wc3Pna93jnu/Pc7YyoDQxqJ2S2CDh8pRwAkMigJicTG+aHD+YnYMuvEzG8nx+qjCb86d/ncOef9uDzQ3ncT5zoOgxqJ3S60IAqowm+aiVnfJPTGj84GF8tmog3H45H/wBPFBmMePEfJzD9rR/w/dlinodN1IhB7YQONHZ7j4kK4BnU5NTkchnui++HjOcm4XfThsBPo8QZXSUe++AQHnnvALLz9VKXSCQ5BrUTOnj5GgAgITpI4kqIuodaqcDCOwZi34uTsfD2aKgUcvznwjXc8/aPWPzpUZwvqpS6RCLJMKidjCAI4kSyxIEcnybX4u+lwu+mD0XGc5NwX3w4AOCfxwtw1/p9WLQ5S9zfnsid8PSsTpLq9KwzOgOmrv8Bnh4KHF99F1RK/l2LXNepAj3e+e4Cvs7WidfuuKUP5k+IwqSYPpBz6IfcAHcmczLfZBcBAMYPCmJIk8sbFq7Fu4+MxhmdAe98dwH/OlmIfedKsO9cCQb28cb88VG4/7b+8FbzRxm5LraoO0mqFvXU9ftwRleJ1x4YgQfHRPTa5xI5gtxrNfgo8wo+P5SHSqMJAOCrUeKhMRF4eGwEYkJ9Ja6QqPsxqDtJiqC+UlqNO//0PZRyGQ6vTIa/l6pXPpfI0VQZTfj74Tx8uP8KrlyrEa+PivTHw2MjMH1EOHzYyiYXwaDuJCmC+t3vL2Lt7jO4PSYYf1uQ2CufSeTILBYBe8+V4NODufjuTDFMjZuleKkUuGdEX8waG4HbIgN4Yhc5Nf6V04l8nV0IAJg6PEziSogcQ/PTuoor67A9Kx+fH8rDpdJqfH74Kj4/fBWD+njj/tv64774cPQP8JK6ZKIOY4u6k3q7RX1WV4mU9fsglwEHViSjj6+6xz+TyBkJgoDDOeXYejAPu04WorahaQ/xhKhAzBjVD9Pj+kLrxVPnyDkwqDupt4P6xb8fx+eHr2JaXBj+PGd0j38ekSuorGvA1yd12HEsH5mXrsH2006lkOP2mGCkDA9D8pBQBHpzvgc5LgZ1J/VmUJdUGjFhzXeoN1vwj9+Mx+gBAT36eUSuqFBfi6+OFeCLo/k4o2va6UwuAxKiA3F7TB+MHxSEuH5aKBVc+kiOg0HdSb0Z1K+nn8NbGecxKtIfX/x2Qo9+FpE7OFdUid3ZOuzO1uHn63Y781ErMXpAAIaG+yE2zBdD+vohMtALGg+FRNWSu2NQd1JvBXWxoQ5TXt+LyjoT3vnvUbhnRHiPfRaRO8orq8F3Z4qx/2IpMi9eg6HO1Op9Wk8PhPqp0cdXDS+VEhoPBTw95FArFVDIZfBQyKCQyxt/lUEpl0GpkFt/lcugUMjhIbc+56lSwFfjAV+NEn4apfjPnh4KzlCnFhjUndRbQf3UlizsPFGIEf21+OK3E3haFlEPMlsE/FxgwLGrFThTaMDpQgPO6ipRXW+++Yu7gVIug2+z4G7+z34aD7tQ99V4QOvpgT6+aoT4quHv5cGQd1FcnuXA9pwtxs4ThZDLgD/+Ko4hTdTDFHIZ4vprEddfK14TBAGGWhOKK+tQZDCitMqI2gYzauvNqG0ww2iywGyxwGQRYDILMFsENJgtjb8KMFssaLAIMJsF6z0WC2rqzaisM6GyrkH81SIAJouA8poGlNc0dLh2lUKOPr7WFn+4vwbhWk+E+1sf/QOsvwYwzJ0Sg9pBXSiuxDNbjwEA5k+IxvB+2hu/gIh6hEwmg9bLA1ovjx7bolQQBLvwNtQ1wFBnahHmlc2uGepMqKipR0mlEeU1Dag3W5BfUYv8ilocy2v9czQecoT7e6Jf48MW5Lbfh2k1PEPAAbHru5N6suv7ankNHtqYiQJ9HeIj/PHpwnHwVHEiCxG1zmgyo6TSiOJKI4oNRhTqa1FQUYuCijoxvEsqjTd9H5kM6OOjtoZ3gCfC/DQI9lEjyEeFPo2/BvmoEeSt4uS6XsSg7qSeCuqM00V4bttxVNQ0YGAfb/z9yfFc40lEXWY0maHTW4O7oKIO+eWNYa6vbbxWi7oGS7vfz1etRJCPSgzyIB81gn3U8Pf0gI9aCW+1Ej4aJXzUCnirldAoFdB4KKDxkEPjoYBKIecxpe3kEEG9YcMGvPbaa9DpdBg5ciTefvttJCQktHn/tm3b8NJLL+HKlSuIiYnB2rVrMW3aNPF5QRCwevVq/PWvf0VFRQUmTJiAd999FzExMeI9ZWVlWLx4Mf75z39CLpdj5syZePPNN+Hj49OumrszqAVBQFZuBTbsuYDvzhQDAEb2tx7vF+7v2aX3JiJqD0EQUFZdL7bCCypqUWSoQ2lVPUqrjLhWbcS1xn9uMHdPbKiUcmiU8sYAt4a4WtkU5rZ/bn5N6+mBQG8VgrxV1l99VAj0tv4FwVWDX/Kg/uyzzzB37lxs3LgRiYmJWL9+PbZt24azZ88iJCSkxf379+/HHXfcgbS0NNxzzz3YsmUL1q5di6ysLAwfPhwAsHbtWqSlpeGjjz5CdHQ0XnrpJZw8eRI///wzNBoNAODuu+9GYWEh/vKXv6ChoQHz58/H2LFjsWXLlnbV3R1BXVJpxMa9F7E7W4f8iloA1s0XHhsfjWV3x3KsiIgcjiAIMNSZcK3KiNKq+sZfjWKgG+pMqKprQLXRjCqjCVVGE6qNJtQ1mFFnsk6y6wlyGRDgZQ3vpgC3hniglwcCG7vsbSHv5+nhNK16yYM6MTERY8eOxTvvvAMAsFgsiIiIwOLFi7Fs2bIW98+aNQvV1dXYuXOneG3cuHGIj4/Hxo0bIQgCwsPD8dxzz+H5558HAOj1eoSGhuLDDz/Eww8/jNOnT2Po0KE4dOgQxowZAwDYvXs3pk2bhqtXryI8vOVaZaPRCKOxaYxHr9cjMjISeXl5nQ5qfW0DJr26ByaLAI2HHHcNDcPCO6IRHdy+Vj0RkbMxmS2oM1lgbDCL4W1snD1f32BBndkMY4OlxXPGBgvqTGYYahtQVlOP8uoGlNfUo6zaiMq6zi+f81BY18CrFHKolHJ4KOVQKeRQyGWQy2SQyWSQyQA5ZJDLrJMLo4K9kHb/iG758/D19b35THxBQkajUVAoFMIXX3xhd33u3LnCvffe2+prIiIihDfeeMPu2qpVq4QRI0YIgiAIFy9eFAAIR48etbvnjjvuEJ5++mlBEAThvffeE/z9/e2eb2hoEBQKhbB9+/ZWP3f16tUCAD744IMPPvjotoder79pVkq6PKu0tBRmsxmhoaF210NDQ3HmzJlWX6PT6Vq9X6fTic/brt3onuu71ZVKJQIDA8V7rrd8+XKkpqaKv7dYLCgrK0NQUFCLvw0ZDAZERER0qbUtBdbdu1h373HGmgHW3dukqNvX9+ZL/riOup3UajXUavujJf39/W/4Gj8/P6f6j9SGdfcu1t17nLFmgHX3NkerW9LZSsHBwVAoFCgqKrK7XlRUhLCwsFZfExYWdsP7bb/e7J7i4mK7500mE8rKytr8XCIiIilIGtQqlQqjR49GRkaGeM1isSAjIwNJSUmtviYpKcnufgBIT08X74+OjkZYWJjdPQaDAQcOHBDvSUpKQkVFBY4cOSLe891338FisSAxMbHbvh8REVGX3XQUu4dt3bpVUKvVwocffij8/PPPwuOPPy74+/sLOp1OEARBePTRR4Vly5aJ9//nP/8RlEql8Kc//Uk4ffq0sHr1asHDw0M4efKkeM+aNWsEf39/4csvvxROnDgh3HfffUJ0dLRQW1sr3jN16lRh1KhRwoEDB4Qff/xRiImJEWbPnt0t36murk5YvXq1UFdX1y3v11tYd+9i3b3HGWsWBNbd2xy1bsmDWhAE4e233xYiIyMFlUolJCQkCD/99JP43KRJk4R58+bZ3f/5558Lt9xyi6BSqYRhw4YJ//rXv+yet1gswksvvSSEhoYKarVamDJlinD27Fm7e65duybMnj1b8PHxEfz8/IT58+cLlZWVPfYdiYiIOkPyddRERETUNm59RURE5MAY1ERERA6MQU1EROTAGNREREQOjEHdja5cuYIFCxYgOjoanp6eGDRoEFavXo36+nq7+06cOIHbb78dGo0GERERePXVVyWq2N6GDRsQFRUFjUaDxMREHDx4UOqSRGlpaRg7dix8fX0REhKCGTNm4OzZs3b31NXVYdGiRQgKCoKPjw9mzpzZYuMbqa1ZswYymQzPPPOMeM1R687Pz8cjjzyCoKAgeHp6Ii4uDocPHxafFwQBq1atQt++feHp6Ynk5GScP39esnrNZjNeeuklu////vCHP6D5fFlHqHnfvn345S9/ifDwcMhkMuzYscPu+fbUWFZWhjlz5sDPzw/+/v5YsGABqqqqJKu7oaEBS5cuRVxcHLy9vREeHo65c+eioKDAoeu+3pNPPgmZTIb169dLXrcdKaecu5qvv/5aeOyxx4RvvvlGuHjxovDll18KISEhwnPPPSfeo9frhdDQUGHOnDlCdna28Omnnwqenp7CX/7yFwkrt65nV6lUwvvvvy+cOnVKWLhwoeDv7y8UFRVJWpdNSkqK8MEHHwjZ2dnCsWPHhGnTpgmRkZFCVVWVeM+TTz4pRERECBkZGcLhw4eFcePGCePHj5ewansHDx4UoqKihBEjRghLliwRrzti3WVlZcKAAQOExx57TDhw4IBw6dIl4ZtvvhEuXLgg3rNmzRpBq9UKO3bsEI4fPy7ce++9LfYr6E2vvPKKEBQUJOzcuVO4fPmysG3bNsHHx0d48803HarmXbt2Cb/73e+E7du3CwBaHErUnhqnTp0qjBw5Uvjpp5+EH374QRg8eHC37QPRmborKiqE5ORk4bPPPhPOnDkjZGZmCgkJCcLo0aPt3sPR6m5u+/btwsiRI4Xw8PAWBz9JUXdzDOoe9uqrrwrR0dHi7//85z8LAQEBgtFoFK8tXbpUuPXWW6UoT5SQkCAsWrRI/L3ZbBbCw8OFtLQ0CatqW3FxsQBA2Lt3ryAI1h8UHh4ewrZt28R7Tp8+LQAQMjMzpSpTVFlZKcTExAjp6enCpEmTxKB21LqXLl0qTJw4sc3nLRaLEBYWJrz22mvitYqKCkGtVguffvppb5TYwvTp04X/+Z//sbt2//33C3PmzBEEwTFrvj442lPjzz//LAAQDh06JN7z9ddfCzKZTMjPz5ek7tYcPHhQACDk5OQIguDYdV+9elXo16+fkJ2dLQwYMMAuqB2hbnZ99zC9Xo/AwEDx95mZmbjjjjugUqnEaykpKTh79izKy8ulKBH19fU4cuQIkpOTxWtyuRzJycnIzMyUpKab0ev1ACD+2R45cgQNDQ123yE2NhaRkZEO8R0WLVqE6dOn29UHOG7dX331FcaMGYMHH3wQISEhGDVqFP7617+Kz1++fBk6nc6ubq1Wi8TERMnqHj9+PDIyMnDu3DkAwPHjx/Hjjz/i7rvvdtiar9eeGjMzM+Hv748xY8aI9yQnJ0Mul+PAgQO9XnNb9Ho9ZDKZeHiRo9ZtsVjw6KOP4oUXXsCwYcNaPO8IdfP0rB504cIFvP322/jTn/4kXtPpdIiOjra7z3Ykp06nQ0BAQK/WCHTuuFEpWSwWPPPMM5gwYQKGDx8OwPpnp1KpWpxo1vx4U6ls3boVWVlZOHToUIvnHLXuS5cu4d1330VqaipWrFiBQ4cO4emnn4ZKpcK8efPadZxsb1u2bBkMBgNiY2OhUChgNpvxyiuvYM6cOQDadwSu1HrqmN7eVldXh6VLl2L27NniKVSOWvfatWuhVCrx9NNPt/q8I9TNoG6HZcuWYe3atTe85/Tp04iNjRV/n5+fj6lTp+LBBx/EwoULe7pEt7Jo0SJkZ2fjxx9/lLqUm8rLy8OSJUuQnp4OjUYjdTntZrFYMGbMGPzxj38EAIwaNQrZ2dnYuHEj5s2bJ3F1rfv888+xefNmbNmyBcOGDcOxY8fwzDPPIDw83GFrdkUNDQ146KGHIAgC3n33XanLuaEjR47gzTffRFZWFmQymdTltIld3+3w3HPP4fTp0zd8DBw4ULy/oKAAkydPxvjx47Fp0ya792rrmE7bc1LozHGjUnnqqaewc+dO7NmzB/379xevh4WFob6+HhUVFXb3S/0djhw5guLiYtx2221QKpVQKpXYu3cv3nrrLSiVSoSGhjpk3X379sXQoUPtrg0ZMgS5ubkA2necbG974YUXsGzZMjz88MOIi4vDo48+imeffRZpaWkAHLPm6zn7Mb22kM7JyUF6errdmc6OWPcPP/yA4uJiREZGiv9/5uTk4LnnnkNUVBQAx6ibQd0Offr0QWxs7A0ftjHn/Px83HnnnRg9ejQ++OADyOX2f8RJSUnYt28fGhoaxGvp6em49dZbJen2Bjp33GhvEwQBTz31FL744gt89913LYYPRo8eDQ8PD7vvcPbsWeTm5kr6HaZMmYKTJ0/i2LFj4mPMmDGYM2eO+M+OWPeECRNaLH87d+4cBgwYAKB9x8n2tpqamhb/vykUClgsFgCOWfP1nPmYXltInz9/Ht9++y2CgoLsnnfEuh999FGcOHHC7v/P8PBwvPDCC/jmm28cp+5embLmJq5evSoMHjxYmDJlinD16lWhsLBQfNhUVFQIoaGhwqOPPipkZ2cLW7duFby8vBxiedaNjhuV2m9+8xtBq9UK33//vd2fa01NjXjPk08+KURGRgrfffedcPjwYSEpKUlISkqSsOrWNZ/1LQiOWffBgwcFpVIpvPLKK8L58+eFzZs3C15eXsInn3wi3tOe42R707x584R+/fqJy7O2b98uBAcHCy+++KJD1VxZWSkcPXpUOHr0qABAeP3114WjR4+Ks6OlPqa3M3XX19cL9957r9C/f3/h2LFjdv+PNl/h4mh1t+b6Wd9S1d0cg7obffDBBwKAVh/NHT9+XJg4caKgVquFfv36CWvWrJGoYns3Om5Uam39uX7wwQfiPbW1tcJvf/tbISAgQPDy8hJ+9atf2f0lyVFcH9SOWvc///lPYfjw4YJarRZiY2OFTZs22T3fnuNke5PBYBCWLFkiREZGChqNRhg4cKDwu9/9zi4oHKHmPXv2tPrfsu04X0c9pvdGdV++fLnN/0f37NnjsHW3prWglvpYZB5zSURE5MA4Rk1EROTAGNREREQOjEFNRETkwBjUREREDoxBTURE5MAY1ERERA6MQU1EROTAGNREREQOjEFNRETkwBjUREREDoxBTURE5MD+P49GPetHSW8QAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# privacy settings for all examples:\n", "max_contributions = 1\n", "epsilon = 1.\n", "\n", "import numpy as np\n", "data = np.random.exponential(20., size=1000)\n", "bounds = 0., 100. # a best guess!\n", "\n", "import seaborn as sns\n", "sns.displot(data, kind=\"kde\");" ] }, { "cell_type": "code", "execution_count": 130, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[5.342567934782284, 13.46163987459231, 27.466767771051295]" ] }, "execution_count": 130, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# true quantiles\n", "true_quantiles = np.quantile(data, [0.25, 0.5, 0.75])\n", "true_quantiles.tolist()\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "----\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": 131, "metadata": {}, "outputs": [], "source": [ "import opendp.prelude as dp\n", "dp.enable_features(\"contrib\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Quantile via Histogram\n", "\n", "One approach for releasing quantiles is to estimate the cumulative distribution via a histogram query.\n", "\n", "The basic procedure for estimating an $\\alpha$-quantile is as follows:\n", "\n", "1. bin the data, and count the number of elements in each bin privately\n", "2. divide the counts by the total to get the probability density of landing in each bin\n", "3. sum the densities while scanning from the left until the sum is at least $\\alpha$\n", "4. interpolate the bin edges of the terminal bin" ] }, { "cell_type": "code", "execution_count": 132, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[5.33, 13.411111111111111, 29.875]" ] }, "execution_count": 132, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "quart_alphas = [0.25, 0.5, 0.75]\n", "input_space = dp.vector_domain(dp.atom_domain(T=float)), dp.symmetric_distance()\n", "\n", "def make_hist_quantiles(alphas, d_in, d_out, num_bins=500):\n", "\n", " edges = np.linspace(*bounds, num=num_bins + 1)\n", " bin_names = [str(i) for i in range(num_bins)]\n", "\n", " def make_from_scale(scale):\n", " return (\n", " input_space >>\n", " dp.t.then_find_bin(edges=edges) >> # bin the data\n", " dp.t.then_index(bin_names, \"0\") >> # can be omitted. Just set TIA=\"usize\", categories=list(range(num_bins)) on next line:\n", " dp.t.then_count_by_categories(categories=bin_names, null_category=False) >>\n", " dp.m.then_base_discrete_laplace(scale) >>\n", " # we're really only interested in the function on this transformation- the domain and metric don't matter\n", " dp.t.make_cast_default(dp.vector_domain(dp.atom_domain(T=int)), dp.symmetric_distance(), TOA=float) >>\n", " dp.t.make_quantiles_from_counts(edges, alphas=alphas)\n", " )\n", "\n", " return dp.binary_search_chain(make_from_scale, d_in, d_out)\n", "\n", "hist_quartiles_meas = make_hist_quantiles(quart_alphas, max_contributions, epsilon)\n", "hist_quartiles_meas(data)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "A drawback of using this algorithm is that it can be difficult to choose the number of bins.\n", "\n", "If the number of bins is chosen to be very small, then the postprocessor will need to sum fewer instances of noise before reaching the bin of interest, resulting in a better bin selection.\n", "However, the bin will be wider, so there will be greater error when interpolating the final answer.\n", "\n", "If the number of bins is chosen to be very large, then the same holds in the other direction.\n", "\n", "Estimating quantiles via the next algorithm can help make choosing the number of bins less sensitive.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Quantile via B-ary Tree\n", "\n", "A slightly more complicated algorithm that tends to provide better utility is to privatize a B-ary tree instead of a histogram.\n", "In this algorithm, the raw counts form the leaf nodes, and a complete tree is constructed by recursively summing groups of size `b`.\n", "This results in a structure where each parent node is the sum of its `b` children.\n", "Noise is added to each node in the tree, and then a postprocessor makes all nodes of the tree consistent with each other, and returns the leaf nodes.\n", "\n", "In the histogram approach, the postprocessor would be influenced by a number of noise sources approximately $O(n)$ in the number of scanned bins.\n", "After this modification, the postprocessor is influenced by a number of noise sources approximately $O(log_b(n))$ in the number of scanned bins, and with noise sources of similarly greater magnitude.\n", "\n", "This modification introduces a new hyperparameter, the branching factor.\n", "`choose_branching_factor` provides a heuristic for the ideal branching factor, based on information (or a best guess) of the dataset size." ] }, { "cell_type": "code", "execution_count": 133, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25" ] }, "execution_count": 133, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = dp.t.choose_branching_factor(size_guess=1_500)\n", "b" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "\n", "We now make the following adjustments to the histogram algorithm:\n", "\n", "* insert a stable (Lipschitz) transformation to construct a b-ary tree before the noise mechanism\n", "* replace the cast postprocessor with a consistency postprocessor\n" ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[5.0371403927139795, 13.146207218782838, 26.037387422664267]" ] }, "execution_count": 134, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def make_tree_quantiles(alphas, b, d_in, d_out, num_bins=500):\n", "\n", " edges = np.linspace(*bounds, num=num_bins + 1)\n", " bin_names = [str(i) for i in range(num_bins)]\n", "\n", " def make_from_scale(scale):\n", " return (\n", " input_space >>\n", " dp.t.then_find_bin(edges=edges) >> # bin the data\n", " dp.t.then_index(bin_names, \"0\") >> # can be omitted. Just set TIA=\"usize\", categories=list(range(num_bins)) on next line:\n", " dp.t.then_count_by_categories(categories=bin_names, null_category=False) >>\n", " dp.t.then_b_ary_tree(leaf_count=len(bin_names), branching_factor=b) >>\n", " dp.m.then_base_discrete_laplace(scale) >> \n", " dp.t.make_consistent_b_ary_tree(branching_factor=b) >> # postprocessing\n", " dp.t.make_quantiles_from_counts(edges, alphas=alphas) # postprocessing\n", " )\n", "\n", " return dp.binary_search_chain(make_from_scale, d_in, d_out)\n", "\n", "tree_quartiles_meas = make_tree_quantiles(quart_alphas, b, max_contributions, epsilon)\n", "tree_quartiles_meas(data)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As mentioned earlier, using the tree-based approach can help make the algorithm less sensitive to the number of bins:" ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABoZklEQVR4nO3dd3yNd//H8dfJ3oOQBEEQsYlYKW211aLjRpeqGi2dtPy6tffdeffW3epSLTVaoxMtpUXRIkVSscUKMSJm9j7n+v1xSKVWIuPKeD8fj/NwxjU+53KS8873+l7fr8UwDAMRERERkziYXYCIiIjUbAojIiIiYiqFERERETGVwoiIiIiYSmFERERETKUwIiIiIqZSGBERERFTOZldQHHYbDYOHz6Mt7c3FovF7HJERESkGAzDID09nXr16uHgcOH2jyoRRg4fPkxISIjZZYiIiMhlOHDgAA0aNLjg61UijHh7ewP2N+Pj42NyNSIiIlIcaWlphISEFH6PX0iVCCNnTs34+PgojIiIiFQxl+pioQ6sIiIiYiqFERERETGVwoiIiIiYqkr0GSkOm81GXl6e2WVIDeHi4nLRy9RERKT4qkUYycvLIyEhAZvNZnYpUkM4ODgQGhqKi4uL2aWIiFR5VT6MGIZBUlISjo6OhISE6K9VKXdnBuFLSkqiYcOGGohPRKSUqnwYKSgoICsri3r16uHh4WF2OVJD1KlTh8OHD1NQUICzs7PZ5YiIVGlVvhnBarUCqLlcKtSZz9uZz5+IiFy+Kh9GzlBTuVQkfd5ERMpOtQkjIiIiUjUpjIiIiIipFEZM0rNnT8aOHXvRZSwWC/PmzauQekRERMyiMFKJJSUl0bdv32Itq+AiIiKXIz0nnzV7jptag8JIJRYUFISrq6vZZZwjPz/f7BIuymq1nncAvMsdoVcj+4pIdXXgZBa3T4xm+NT1/JV4yrQ6ql0YMQyDrLwCU26GYZSoVpvNxtNPP02tWrUICgripZdeKvL62a0deXl5jB49muDgYNzc3GjUqBHjx48HoHHjxgAMGDAAi8VS+Bhg4sSJNG3aFBcXF8LDw/nyyy+L7GPHjh306NEDNzc3WrVqxdKlS4vsd9++fVgsFr7++muuvvpq3NzcmDlzJidOnGDQoEHUr18fDw8P2rZty+zZs4tsu2fPnjz66KOMHTsWf39/AgMD+fzzz8nMzOTee+/F29ubZs2asWjRoosep9zcXJ588knq16+Pp6cnXbt2ZcWKFYWvT5s2DT8/P3788UdatWqFq6sriYmJNG7cmFdffZWhQ4fi4+PDAw88AMD3339P69atcXV1pXHjxrzzzjtF9neh9UREqpPY/acY8Mlq4pPT8XV3xsnBvKsEq/ygZ/+UnW+l1Qu/mLLvba/0xsOl+Id0+vTpPP7446xdu5bo6GiGDx9O9+7duf76689Z9oMPPuDHH3/km2++oWHDhhw4cIADBw4AsH79eurWrcvUqVPp06cPjo6OAMydO5cxY8bw/vvv06tXLxYsWMC9995LgwYNuOaaa7BarfTv35+GDRuydu1a0tPTeeKJJ85b67PPPss777xDREQEbm5u5OTkEBkZyTPPPIOPjw8LFy5kyJAhNG3alC5duhR5j08//TTr1q3j66+/5uGHH2bu3LkMGDCA5557jvfee48hQ4aQmJh4wUHrRo8ezbZt25gzZw716tVj7ty59OnTh82bNxMWFgZAVlYWb7zxBpMnT6Z27drUrVsXgLfffpsXXniBF198EYDY2FjuvPNOXnrpJQYOHMiaNWt45JFHqF27NsOHDy/c5z/XExGpTubHHeKp7zaRV2CjVbAPk4d1op6fu2n1WIyS/jlvgrS0NHx9fUlNTcXHx6fIazk5OSQkJBAaGoqbmxtZeQVVIoz07NkTq9XKH3/8Ufhcly5duPbaa3n99dcBe8vI3Llz6d+/P4899hhbt24tbLn4p7OXPaN79+60bt2azz77rPC5O++8k8zMTBYuXMjixYu55ZZbOHDgAEFBQQAsXbqU66+/vnBb+/btIzQ0lPfff58xY8Zc9D3dfPPNtGjRgrfffvu879FqteLr68utt97KjBkzADhy5AjBwcFER0fTrVu3c7aZmJhIkyZNSExMpF69eoXP9+rViy5duvC///2PadOmce+99xIXF0f79u0Ll2ncuDERERHMnTu38LnBgwdz7Ngxfv3118Lnnn76aRYuXMjWrVsvuN4//fNzJyJSFRiGwXtLd/HBsl0A9GoZyIS7OuDpWj5tExf7/j5btWsZcXd2ZNsrvU3bd0m0a9euyOPg4GCOHj163mWHDx/O9ddfT3h4OH369OHmm2/mhhtuuOj2t2/ffs4phu7duzNhwgQA4uPjCQkJKQwiQJFWjbN16tSpyGOr1cr//vc/vvnmGw4dOkReXh65ubnntG6c/R4dHR2pXbs2bdu2LXwuMDAQ4ILve/PmzVitVpo3b17k+dzcXGrXrl342MXF5Zzjeb66t2/fTr9+/Yo81717d95//32sVmthq9I/1xMRqepy8q08+e1GFmxKAuDBq5rwdJ8WOJp4euaMahdGLBZLiU6VmOmfc5pYLJYLzjzcsWNHEhISWLRoEUuXLuXOO++kV69efPfddxVRKp6enkUev/XWW0yYMIH333+ftm3b4unpydixY8/p7Hm+93j2c2daeS70vjMyMnB0dCQ2NrYwKJzh5eVVeN/d3f28LUb/rLu4Lnc9EZHK6Gh6DvfPiGXjgRScHCy8NqANAzs3NLusQlXjW1sA8PHxYeDAgQwcOJDbb7+dPn36cPLkSWrVqoWzs/M586S0bNmS1atXM2zYsMLnVq9eTatWrQAIDw/nwIEDJCcnF7ZQrF+/vli1rF69mn79+nHPPfcA9jCxc+fOwm2XlYiICKxWK0ePHuXKK68s9fbOHJOzrV69mubNm58TdkREqoPtSWmMmLaew6k5+Hk4M3FwJFFNa196xQqkMFJFvPvuuwQHBxMREYGDgwPffvstQUFB+Pn5AfZ+DsuWLaN79+64urri7+/PU089xZ133klERAS9evXip59+4ocffmDp0qUAXH/99TRt2pRhw4bx5ptvkp6ezr///W/g0nOvhIWF8d1337FmzRr8/f159913SU5OLvMw0rx5cwYPHszQoUMLO9AeO3aMZcuW0a5dO2666aYSbe+JJ56gc+fOvPrqqwwcOJDo6Gg++ugjPvnkkzKtW0SkMli6LZnH5mwgK89KkwBPpgzvTGhA5Wv5rXaX9lZX3t7evPnmm3Tq1InOnTuzb98+fv75Zxwc7P+F77zzDkuWLCEkJISIiAgA+vfvz4QJE3j77bdp3bo1kyZNYurUqfTs2ROw9+GYN28eGRkZdO7cmZEjR/L8888DXLJT5r///W86duxI79696dmzJ0FBQUU6z5alqVOnMnToUJ544gnCw8Pp378/69evp2HDkjcxduzYkW+++YY5c+bQpk0bXnjhBV555ZUiV9KIiFR1hmEw+Y+93P9lDFl5Vq5oWpu5j3SvlEEEquHVNFI6q1evpkePHuzevZumTZuaXU6lpc+diFRW+VYbL8zfwux19uEfBnVpyCv9WuPsWPHtDzX2ahopmblz5+Ll5UVYWBi7d+9mzJgxdO/eXUFERKQKSs3K5+GZsazZcwKLBZ6/sSUjeoRe8tS72RRGarj09HSeeeYZEhMTCQgIoFevXueMSCoiIpVfwvFMRkxbz97jmXi6OPLBoAiuaxlodlnFojBSww0dOpShQ4eaXYaIiJRC9J4TPPRVLKnZ+dT3c2fysE60DL7waZHKRmFERESkCvt6fSLPz91Cgc2gQ4gfnw2NpK531erLpjAiIiJSBVltBm8s3sFnv+8F4OZ2wbx9R3vcSjgaeGWgMCIiIlLFZOYWMGZOHEu3JwMw5rowxvYKq/QdVS+kRNf5TJw4kXbt2uHj44OPjw9RUVEXnf592rRpWCyWIjddBikiInL5Dqdkc/un0SzdnoyLkwMT7urA/13fvMoGEShhy0iDBg14/fXXCQsLwzAMpk+fTr9+/diwYQOtW7c+7zo+Pj7Ex8cXPq7KB0tERMRMGw+kMHJGDMfScwnwcmHSkE5ENvI3u6xSK1EYueWWW4o8fu2115g4cSJ//vnnBcOIxWIpMiusiIiIlNzCTUk8/k0cuQU2wgO9mTK8Ew38PS69YhVw2cOxWa1W5syZQ2ZmJlFRURdcLiMjg0aNGhESEkK/fv3YunXrJbedm5tLWlpakVt107NnT8aOHWt2GSIiUskZhsGHy3YxatZf5BbYuCa8Dt89HFVtgghcRhjZvHkzXl5euLq68tBDDzF37twLTo4WHh7OF198wfz58/nqq6+w2WxcccUVHDx48KL7GD9+PL6+voW3kJCQkpZZLRiGQUFBgdllFMrLyzvnOavVis1mK/G2Lnc9EZGaJCffyv99Hcc7S3YCcF/3UCYP64y3m7PJlZWtEoeR8PBw4uLiWLt2LQ8//DDDhg1j27Zt5102KiqKoUOH0qFDB66++mp++OEH6tSpw6RJky66j3HjxpGamlp4O3DgQEnLrNSGDx/OypUrmTBhQmHH3n379rFixQosFguLFi0iMjISV1dXVq1ahc1mY/z48YSGhuLu7k779u357rvvimxzy5Yt9O3bFy8vLwIDAxkyZAjHjx+/aB2rVq3iyiuvxN3dnZCQEB577DEyMzMLX2/cuDGvvvoqQ4cOxcfHhwceeIBp06bh5+fHjz/+SKtWrXB1dSUxMZFTp04xdOhQ/P398fDwoG/fvuzatatwWxdaT0REzu94Ri6DJ69lXtxhHB0s/Ld/G164pRWODtWv72WJw4iLiwvNmjUjMjKS8ePH0759eyZMmFCsdZ2dnYmIiGD37t0XXc7V1bXwip0zt2IzDMjLNOdWzDkHJ0yYQFRUFPfffz9JSUkkJSUVaf159tlnef3119m+fTvt2rVj/PjxzJgxg08//ZStW7fyf//3f9xzzz2sXLkSgJSUFK699loiIiKIiYlh8eLFJCcnc+edd16whj179tCnTx9uu+02Nm3axNdff82qVasYPXp0keXefvtt2rdvz4YNG/jPf/4DQFZWFm+88QaTJ09m69at1K1bl+HDhxMTE8OPP/5IdHQ0hmFw4403kp+fX7it860nIiLn2pmcTv+PVxO7/xTebk5Mv7cL93RrZHZZ5abU44zYbDZyc3OLtazVamXz5s3ceOONpd3theVnwf/qld/2L+a5w+By6emZfX19cXFxwcPD47yde1955RWuv/56wN5/5n//+x9Lly4t7JvTpEkTVq1axaRJk7j66qv56KOPiIiI4H//+1/hNr744gtCQkLYuXMnzZs3P2cf48ePZ/DgwYX9VsLCwvjggw+4+uqrmThxYuEl2Ndeey1PPPFE4Xp//PEH+fn5fPLJJ7Rv3x6AXbt28eOPP7J69WquuOIKAGbOnElISAjz5s3jjjvuADhnPREROdeK+KOMnrWBjNwCGtX2YMqwzjSr62V2WeWqRGFk3Lhx9O3bl4YNG5Kens6sWbNYsWIFv/zyC2Cf56R+/fqMHz8esH+pduvWjWbNmpGSksJbb73F/v37GTlyZNm/k2qkU6dOhfd3795NVlZWYTg5Iy8vj4iICAA2btzI8uXL8fI698O6Z8+e84aRjRs3smnTJmbOnFn4nGEY2Gw2EhISaNmy5Tm1nOHi4kK7du0KH2/fvh0nJye6du1a+Fzt2rUJDw9n+/btF1xPRESKmr5mHy//tBWbAV1CazHpnkj8PV3MLqvclSiMHD16lKFDh5KUlISvry/t2rXjl19+KfyiTExMxMHh7zM/p06d4v777+fIkSP4+/sTGRnJmjVrLtjhtUw4e9hbKMzgXDY9mz09/25dycjIAGDhwoXUr1+/yHKurq6Fy9xyyy288cYb52wrODj4vPvIyMjgwQcf5LHHHjvntYYNG563ljPc3d0va7yYy11PRKS6K7DaePmnbXz5534A7ohswGsD2uLidNkXvVYpJQojU6ZMuejrK1asKPL4vffe47333itxUaVisRTrVInZXFxcsFqtl1zu7M6eV1999XmX6dixI99//z2NGzfGyal4/6UdO3Zk27ZtNGvWrER1n0/Lli0pKChg7dq1hadpTpw4QXx8fPkGTxGRaiA1O5/Rs/7ij13HsVjgmT4tePCqJjXqj7eaEbkqocaNG7N27Vr27dvH8ePHL3iZq7e3N08++ST/93//x/Tp09mzZw9//fUXH374IdOnTwdg1KhRnDx5kkGDBrF+/Xr27NnDL7/8wr333nvBwPPMM8+wZs0aRo8eTVxcHLt27WL+/PnndGAtjrCwMPr168f999/PqlWr2LhxI/fccw/169enX79+Jd6eiEhNkXgii9smruGPXcdxd3Zk4uBIHrq6aY0KIqAwYponn3wSR0dHWrVqRZ06dS56meurr77Kf/7zH8aPH0/Lli3p06cPCxcuJDQ0FIB69eqxevVqrFYrN9xwA23btmXs2LH4+fkVOW12tnbt2rFy5Up27tzJlVdeSUREBC+88AL16l1e59+pU6cSGRnJzTffTFRUFIZh8PPPP+PsXL2uhRcRKSvrEk7S7+NV7D6aQZCPG98+FEWfNjVzxHKLYRTzelQTpaWl4evrS2pq6jmX+ebk5JCQkEBoaKgm4ZMKo8+diJTG97EHefaHTeRbDdrW92XysE4E+lS/3yUX+/4+W6kv7RUREZHisdkM3v41nk9W7AGgb5sg3r2zA+4ujiZXZi6FERERkQqQnWcf2n3x1iMAjLqmKU9cH45DNRxRtaQURkRERMpZcloOI6fHsPlQKs6OFl6/tR23RTYwu6xKQ2FERESkHG05lMrI6TEcScuhlqcLk4ZE0rlxLbPLqlSqTRipAv1wpRrR501EiuOXrUcYOyeO7Hwrzep68cWwzjSsXTYDZFYnVT6MODraO/3k5eXh7u5ucjVSU+Tl5QF/f/5ERM5mGAafrtzLm7/swDDgyrAAPh7cER83DXdwPlU+jDg5OeHh4cGxY8dwdna+4LgaImXFZrNx7NgxPDw8ij3irYjUHHkFNp6bu5nvYg8CMDSqES/c3AonR30/XUiV/01qsVgIDg4mISGB/fv3m12O1BAODg40bNiwxo2SKCIXdzIzj4e+imVdwkkcLPDiLa0ZdkVjs8uq9Kp8GAH7PC9hYWGFTeci5c3FxUWtcCJSxO6jGYyYvp79J7LwcnXio7sj6Ble1+yyqoRqEUbA/peqRsIUEREzrNp1nIdnxpKeU0ADf3e+GN6Z5oHeZpdVZVSbMCIiImKGr/7cz4s/bsVqM4hs5M9nQyKp7eVqdllVisKIiIjIZbDaDP67cBtTV+8DYEBEfcbf2hY3Z11lV1IKIyIiIiWUnpPPY7M3sDz+GABP3tCcUdc0U6f2y6QwIiIiUgIHTmYxcnoM8cnpuDo58O6dHbipXbDZZVVpCiMiIiLFFLv/FA9+GcPxjDzqeLsyeWgn2of4mV1WlacwIiIiUgzz4w7x1HebyCuw0SrYh8nDOlHPTyN/lwWFERERkYswDIP3lu7ig2W7AOjVMpAJd3XA01VfoWVFR1JEROQCcvKtPPntRhZsSgLgwaua8HSfFjg6qKNqWVIYEREROY+j6Tk8MCOWuAMpODlYeG1AGwZ2bmh2WdWSwoiIiMg/bE9KY8S09RxOzcHPw5mJgyOJalrb7LKqLYURERGRsyzbnsxjszeQmWelSYAnU4Z3JjTA0+yyqjWFEREREewdVaesSuC1n7djGHBF09pMHByJr4ez2aVVewojIiJS4+Vbbbwwfwuz1x0AYFCXhrzSrzXOjpqduyIojIiISI2WmpXPwzNjWbPnBBYLPH9jS0b0CNXQ7hVIYURERGqshOOZjJi2nr3HM/F0cWTCXRH0ahVodlk1jsKIiIjUSNF7TvDQV7GkZudT38+dycM60TLYx+yyaiSFERERqXG+Xp/I83O3UGAz6BDix2dDI6nr7WZ2WTWWwoiIiNQYVpvBG4t38NnvewG4uV0wb9/RHjdnR5Mrq9kURkREpEbIzC1gzJw4lm5PBmDMdWGM7RWmjqqVgMKIiIhUe4dTshkxPYbtSWm4ODnw1u3t6NehvtllyWkKIyIiUq1tPJDCyBkxHEvPJcDLhUlDOhHZyN/ssuQsCiMiIlJtLdyUxOPfxJFbYCM80JspwzvRwN/D7LLkHxRGRESk2jEMg49+2807S3YCcE14HT4YFIG3m4Z2r4wURkREpFrJLbDy7PebmbvhEAD3dQ/l+Zta4uigjqqVlcKIiIhUG8czcnnwy1hi95/C0cHCy/9qzT3dGpldllyCwoiIiFQLO5PTuW/aeg6eysbbzYmJgyPpERZgdllSDAojIiJS5a2IP8roWRvIyC2gUW0PpgzrTLO6XmaXJcWkMCIiIlXa9DX7ePmnrdgM6BJai0n3ROLv6WJ2WVICCiMiIlIlFVhtvPzTNr78cz8Ad0Q24LUBbXFxcjC5MikphREREalyUrPzGT3rL/7YdRyLBZ7p04IHr2qiod2rKIURERGpUhJPZHHf9PXsPpqBu7Mj7w3sQJ82QWaXJaWgMCIiIlXG+n0neWBGDKey8gnycWPysE60qe9rdllSSgojIiJSJXwfe5BxP2wmz2qjbX1fJg/rRKCPm9llSRlQGBERkUrNZjN4+9d4PlmxB4C+bYJ4984OuLs4mlyZlJUSdTmeOHEi7dq1w8fHBx8fH6Kioli0aNFF1/n2229p0aIFbm5utG3blp9//rlUBYuISM2RnWflkZl/FQaRUdc05eO7OyqIVDMlCiMNGjTg9ddfJzY2lpiYGK699lr69evH1q1bz7v8mjVrGDRoECNGjGDDhg3079+f/v37s2XLljIpXkREqq/ktBzunBTN4q1HcHa08M4d7XmqdwscNMdMtWMxDMMozQZq1arFW2+9xYgRI855beDAgWRmZrJgwYLC57p160aHDh349NNPi72PtLQ0fH19SU1NxcfHpzTliohIFbDlUCojp8dwJC2HWp4uTBoSSefGtcwuS0qouN/flz0yjNVqZc6cOWRmZhIVFXXeZaKjo+nVq1eR53r37k10dPRFt52bm0taWlqRm4iI1Ay/bD3CHZ9GcyQth2Z1vZj3SHcFkWquxB1YN2/eTFRUFDk5OXh5eTF37lxatWp13mWPHDlCYGBgkecCAwM5cuTIRfcxfvx4Xn755ZKWJiIiVZhhGHy6ci9v/rIDw4ArwwL4eHBHfNyczS5NylmJW0bCw8OJi4tj7dq1PPzwwwwbNoxt27aVaVHjxo0jNTW18HbgwIEy3b6IiFQueQU2nv5uE28stgeRoVGNmDq8s4JIDVHilhEXFxeaNWsGQGRkJOvXr2fChAlMmjTpnGWDgoJITk4u8lxycjJBQRcfKc/V1RVXV9eSliYiIlXQycw8HvoqlnUJJ3GwwIu3tGbYFY3NLksqUKlnE7LZbOTm5p73taioKJYtW1bkuSVLllywj4mIiNQsu49mMOCT1axLOImXqxNfDO+sIFIDlahlZNy4cfTt25eGDRuSnp7OrFmzWLFiBb/88gsAQ4cOpX79+owfPx6AMWPGcPXVV/POO+9w0003MWfOHGJiYvjss8/K/p2IiEiVsmrXcR6eGUt6TgEN/N35Ynhnmgd6m12WmKBEYeTo0aMMHTqUpKQkfH19adeuHb/88gvXX389AImJiTg4/N3YcsUVVzBr1iz+/e9/89xzzxEWFsa8efNo06ZN2b4LERGpUr76cz8v/rgVq80gspE/k4ZEEuCl0/M1VanHGakIGmdERKR6sNoM/rtwG1NX7wNgQER9xt/aFjdnjahaHRX3+1tz04iISIVIz8nnsdkbWB5/DIAnb2jOqGuaYbFoRNWaTmFERETK3YGTWYycHkN8cjquTg68e2cHbmoXbHZZUkkojIiISLmK3X+KB7+M4XhGHnW8XZk8tBPtQ/zMLksqEYUREREpN/PjDvHUd5vIK7DRKtiHycM6Uc/P3eyypJJRGBERkTJnGAbvLd3FB8t2AdCrZSAT7uqAp6u+duRc+lSIiEiZysm38uS3G1mwKQmAB69qwtN9WuDooI6qcn4KIyIiUmaOpufwwIxY4g6k4ORg4bUBbRjYuaHZZUklpzAiIiJlYntSGiOmredwag5+Hs5MHBxJVNPaZpclVYDCiIiIlNqy7ck8NnsDmXlWmgR4MmV4Z0IDPM0uS6oIhREREblshmEwZVUCr/28HcOAK5rWZuLgSHw9nM0uTaoQhREREbks+VYbL8zfyux1iQAM6tKQV/q1xtmx1BPCSw2jMCIiIiWWmpXPwzNjWbPnBBYLPH9jS0b0CNXQ7nJZFEZERKREEo5nMmLaevYez8TTxZEJd0XQq1Wg2WVJFaYwIiIixRa95wQPfRVLanY+9XzdmDK8My2DNZu6lI7CiIiIFMvX6xN5fu4WCmwG7UP8+HxoJHW93cwuS6oBhREREbkoq83gzcU7mPT7XgBubhfM23e0x83Z0eTKpLpQGBERkQvKzC1g7NdxLNmWDMCY68IY2ytMHVWlTCmMiIjIeR1OyWbE9Bi2J6Xh4uTAW7e3o1+H+maXJdWQwoiIiJxj44EURs6I4Vh6LgFeLkwa0onIRv5mlyXVlMKIiIgUsXBTEo9/E0dugY3wQG+mDO9EA38Ps8uSakxhREREAPvQ7h/9tpt3luwE4JrwOnwwKAJvNw3tLuVLYURERMgtsPLs95uZu+EQAPd1D+X5m1ri6KCOqlL+FEZERGq44xm5PPhlLLH7T+HoYOHlf7Xmnm6NzC5LahCFERGRGmxncjr3TVvPwVPZeLs5MXFwJD3CAswuS2oYhRERkRpqRfxRHp21gfTcAhrV9mDKsM40q+tldllSAymMiIjUQNPX7OPln7ZiM6BLaC0m3ROJv6eL2WVJDaUwIiJSgxRYbbyyYBszovcDcEdkA14b0BYXJweTK5OaTGFERKSGSMvJZ9TMv/hj13EsFnimTwsevKqJhnYX0ymMiIjUAIknsrhv+np2H83A3dmR9wZ2oE+bILPLEgEURkREqr31+07ywIwYTmXlE+jjypRhnWlT39fsskQKKYyIiFRj38ceZNwPm8mz2mhb35fPh3YiyNfN7LJEilAYERGphmw2g3eWxPPx8j0A9G0TxLt3dsDdxdHkykTOpTAiIlLNZOdZefybOBZtOQLAqGua8sT14ThoaHeppBRGRESqkeS0HEZOj2HzoVScHS28fms7botsYHZZIhelMCIiUk1sOZTKyOkxHEnLoZanC5OGRNK5cS2zyxK5JIUREZFq4JetRxg7J47sfCvN6nrxxbDONKztYXZZIsWiMCIiUoUZhsGnK/fy5i87MAy4MiyAjwd3xMfN2ezSRIpNYUREpIrKK7Dx/NzNfBt7EIAh3Rrx4i2tcHLU0O5StSiMiIhUQScz83joq1jWJZzEwQIv3tKaYVc0NrsskcuiMCIiUsXsPprBiOnr2X8iCy9XJz66O4Ke4XXNLkvksimMiIhUIat2HefhmbGk5xTQwN+dL4Z3pnmgt9lliZSKwoiISBXx1Z/7efHHrVhtBpGN/Jk0JJIAL1ezyxIpNYUREZFKzmozeG3hdr5YnQDAgIj6jL+1LW7OGtpdqgeFERGRSiw9J58xc+L4bcdRAJ68oTmjrmmGxaKh3aX6UBgREamkDpzMYuT0GOKT03F1cuDdOztwU7tgs8sSKXMKIyIilVDs/lM8+GUMxzPyqOPtyuShnWgf4md2WSLlQmFERKSSmR93iKe+20RegY2WwT5MGdaJen7uZpclUm4URkREKgnDMHh/6S4mLNsFQK+WgUy4qwOervpVLdVbicYMHj9+PJ07d8bb25u6devSv39/4uPjL7rOtGnTsFgsRW5ubm6lKlpEpLrJybfy2Jy4wiDy4FVNmDQkUkFEaoQSfcpXrlzJqFGj6Ny5MwUFBTz33HPccMMNbNu2DU9Pzwuu5+PjUyS0qBe4iMjfjqbn8MCMWOIOpODkYOG1AW0Y2Lmh2WWJVJgShZHFixcXeTxt2jTq1q1LbGwsV1111QXXs1gsBAUFFXs/ubm55ObmFj5OS0srSZkiIlXG9qQ0Rkxbz+HUHPw8nJk4OJKoprXNLkukQpVqasfU1FQAatWqddHlMjIyaNSoESEhIfTr14+tW7dedPnx48fj6+tbeAsJCSlNmSIildKy7cncPnENh1NzaBLgydxHuiuISI1kMQzDuJwVbTYb//rXv0hJSWHVqlUXXC46Oppdu3bRrl07UlNTefvtt/n999/ZunUrDRo0OO8652sZCQkJITU1FR8fn8spV0Sk0jAMgymrEnjt5+0YBlzRtDYTB0fi6+FsdmkiZSotLQ1fX99Lfn9fdhh5+OGHWbRoEatWrbpgqDif/Px8WrZsyaBBg3j11VeLtU5x34yISGWXb7XxwvytzF6XCMCgLiG80q8Nzo6laqgWqZSK+/19Wd20R48ezYIFC/j9999LFEQAnJ2diYiIYPfu3ZezaxGRKis1K5+HZ8ayZs8JLBZ4/saWjOgRqk79UuOVKIobhsHo0aOZO3cuv/32G6GhoSXeodVqZfPmzQQHa0hjEak5Eo5nMuCT1azZcwJPF0c+H9KJkVc2URARoYQtI6NGjWLWrFnMnz8fb29vjhw5AoCvry/u7vbRAYcOHUr9+vUZP348AK+88grdunWjWbNmpKSk8NZbb7F//35GjhxZxm9FRKRyit5zgoe+iiU1O596vm5MGd6ZlsE65SxyRonCyMSJEwHo2bNnkeenTp3K8OHDAUhMTMTB4e8Gl1OnTnH//fdz5MgR/P39iYyMZM2aNbRq1ap0lYuIVAFfr0/k+blbKLAZtA/x4/OhkdT11sCPIme77A6sFUkdWEWkqrHaDN5cvINJv+8F4OZ2wbx9R3vcnB1Nrkyk4pRrB1YREbmwzNwCxn4dx5JtyQCMuS6Msb3C1D9E5AIURkREytDhlGxGTI9he1IaLk4OvHV7O/p1qG92WSKVmsKIiEgZ2XgghZEzYjiWnkuAlwuThnQispG/2WWJVHoKIyIiZWDhpiQe/yaO3AIb4YHeTB7WiZBaHmaXJVIlKIyIiJSCYRh8vHw3b/+6E4BrwuvwwaAIvN00tLtIcSmMiIhcptwCK89+v5m5Gw4BcF/3UJ6/qSWODuqoKlISCiMiIpfheEYuD34ZS+z+Uzg6WHj5X625p1sjs8sSqZIURkRESmhncjr3TVvPwVPZeLs5MXFwJD3CAswuS6TKUhgRESmBFfFHeXTWBtJzC2hU24MpwzrTrK6X2WWJVGkKIyIixTR9zT5e/mkrNgO6hNZi0j2R+Hu6mF2WSJWnMCIicgkFVhuvLNjGjOj9ANwe2YD/DWiLi1OJJj4XkQtQGBERuYi0nHxGzfyLP3Ydx2KBZ/q04MGrmmhod5EypDAiInIBiSeyuG/6enYfzcDd2ZH3BnagT5sgs8sSqXYURkREzmP9vpM8+GUsJzPzCPRxZcqwzrSp72t2WSLVksKIiMg/fB97kHE/bCbPaqNtfV8+H9qJIF83s8sSqbYURkRETrPZDN5ZEs/Hy/cA0Kd1EO8N7IC7i6PJlYlUbwojIiJAdp6Vx7+JY9GWIwA80rMpT94QjoOGdhcpdwojIlLjJaflMHJ6DJsPpeLsaOH1W9txW2QDs8sSqTEURkSkRttyKJWR02M4kpaDv4czk4Z0oktoLbPLEqlRFEZEpMb6ZesRxs6JIzvfSrO6XkwZ1olGtT3NLkukxlEYEZEaxzAMJv2+lzcW78Aw4MqwAD4e3BEfN2ezSxOpkRRGRKRGySuw8fzczXwbexCAId0a8eItrXBy1NDuImZRGBGRGuNkZh4PfRXLuoSTOFjgxVtaM+yKxmaXJVLjKYyISI2w+2gGI6avZ/+JLLxcnfjo7gh6htc1uywRQWFERGqAVbuO8/DMWNJzCmjg786UYZ0JD/I2uywROU1hRESqtZlr9/PC/K1YbQaRjfyZNCSSAC9Xs8sSkbMojIhItWS1Gby2cDtfrE4AYEBEfcbf2hY3Zw3tLlLZKIyISLWTnpPPmDlx/LbjKABPXN+c0dc2w2LR0O4ilZHCiIhUKwdPZTFiWgzxyem4Ojnw7p0duKldsNllichFKIyISLXxV+IpHpgRw/GMPOp4uzJ5aCfah/iZXZaIXILCiIhUC/PjDvHUd5vIK7DRMtiHKcM6Uc/P3eyyRKQYFEZEpEozDIP3l+5iwrJdAPRqGciEuzrg6apfbyJVhX5aRaTKysm38tR3m/hp42EAHriqCc/0aYGjgzqqilQlCiMiUiUdTc/hgRmxxB1IwcnBwmsD2jCwc0OzyxKRy6AwIiJVzvakNEZOj+FQSja+7s58ek8kUU1rm12WiFwmhRERqVKWbU/msdkbyMyz0iTAkynDOxMa4Gl2WSJSCgojIlIlGIbBlFUJvPbzdgwDrmham4mDI/H1cDa7NBEpJYUREan08q02Xpi/ldnrEgEY1CWEV/q1wdnRweTKRKQsKIyISKWWmpXPwzNjWbPnBBYLPH9jS0b0CNXQ7iLViMKIiFRaCcczGTFtPXuPZ+Lh4sgHd0XQq1Wg2WWJSBlTGBGRSil6zwke+iqW1Ox86vm6MXlYZ1rV8zG7LBEpBwojIlLpfLP+AM/N3UyBzaB9iB+fD42krreb2WWJSDlRGBGRSsNqM3hz8Q4m/b4XgJvbBfP2He1xc3Y0uTIRKU8KIyJSKWTmFjD26ziWbEsG4LHrwhh7XRgOGtpdpNpTGBER0yWlZjNiWgzbktJwcXLgrdvb0a9DfbPLEpEKojAiIqbaeCCF+2fEcDQ9lwAvFyYN6URkI3+zyxKRCqQwIiKmWbgpice/iSO3wEZ4oDeTh3UipJaH2WWJSAUr0fCF48ePp3Pnznh7e1O3bl369+9PfHz8Jdf79ttvadGiBW5ubrRt25aff/75sgsWkarPMAw++m0Xo2b9RW6BjZ7hdfju4SgFEZEaqkRhZOXKlYwaNYo///yTJUuWkJ+fzw033EBmZuYF11mzZg2DBg1ixIgRbNiwgf79+9O/f3+2bNlS6uJFpOrJLbDy+DcbefvXnQDc270xk4d2wttNc8yI1FQWwzCMy1352LFj1K1bl5UrV3LVVVedd5mBAweSmZnJggULCp/r1q0bHTp04NNPPz3vOrm5ueTm5hY+TktLIyQkhNTUVHx8NOiRSFV1IiOXB7+MJWb/KRwdLLz8r9bc062R2WWJSDlJS0vD19f3kt/fpZplKjU1FYBatWpdcJno6Gh69epV5LnevXsTHR19wXXGjx+Pr69v4S0kJKQ0ZYpIJbAzOZ3+n6wmZv8pvN2cmH5vFwUREQFKEUZsNhtjx46le/futGnT5oLLHTlyhMDAonNJBAYGcuTIkQuuM27cOFJTUwtvBw4cuNwyRaQSWBF/lNs+WcOBk9k0qu3B3Ee60yMswOyyRKSSuOyraUaNGsWWLVtYtWpVWdYDgKurK66urmW+XRGpeNPX7OPln7ZiM6BLaC0m3ROJv6eL2WWJSCVyWWFk9OjRLFiwgN9//50GDRpcdNmgoCCSk5OLPJecnExQUNDl7FpEqogCq41XFmxjRvR+AG6PbMD/BrTFxalUZ4dFpBoq0W8FwzAYPXo0c+fO5bfffiM0NPSS60RFRbFs2bIizy1ZsoSoqKiSVSoiVUZaTj73TlvPjOj9WCzwbN8WvHV7OwURETmvErWMjBo1ilmzZjF//ny8vb0L+334+vri7u4OwNChQ6lfvz7jx48HYMyYMVx99dW888473HTTTcyZM4eYmBg+++yzMn4rIlIZJJ7I4r7p69l9NAN3Z0feG9iBPm3UEioiF1aiP1MmTpxIamoqPXv2JDg4uPD29ddfFy6TmJhIUlJS4eMrrriCWbNm8dlnn9G+fXu+++475s2bd9FOryJSNa3fd5L+n6xm99EMAn1c+fahKAUREbmkUo0zUlGKe52yiJjn+9iDjPthM3lWG23r+/L50E4E+bqZXZaImKi439+am0ZESsVmM3hnSTwfL98DQJ/WQbw3sAPuLo4mVyYiVYXCiIhctuw8K49/E8eiLfb+Y4/0bMqTN4Tj4GAxuTIRqUoURkTksiSn5TByegybD6Xi7Gjh9VvbcVvkxS/1FxE5H4URESmxLYdSGTk9hiNpOfh7ODNpSCe6hF54WggRkYtRGBGREvll6xHGzokjO99Ks7peTBnWiUa1Pc0uS0SqsJodRlIPQeZRqNUE3HzNrkakUjMMg0m/7+WNxTswDLgyLICP7u6Ir7uz2aWJSBVXs4dDnD0QPusJB9abXYlIpZZXYOPp7zbx+iJ7EBnSrRFTh3dWEBGRMlGzW0Y8Ts8amnXC3DpEKrFTmXk8+FUs6xJO4mCBF29pzbArGptdlohUIzU8jNS2/6swInJeu49mMGL6evafyMLL1YmP7o6gZ3hds8sSkWqmRoeRLGc/PIC89KNoQnORolbtOs7DM2NJzymggb87U4Z1JjzI2+yyRKQaqtF9Rr7fng3A8eTDJlciUrnMXLufYVPXkZ5TQGQjf+aN6q4gIiLlpka3jOS7+kMOkHXS7FJEKgWrzeC1hdv5YnUCAAMi6jP+1ra4OWtodxEpPzU6jFjda0EqOGQrjIik5+QzZk4cv+04CsAT1zdn9LXNsFg0tLuIlK8aHUYsnvaraZxzFUakZjt4KosR02KIT07H1cmBd+/swE3tgs0uS0RqiBodRhw87VfTuOWfMrkSEfP8lXiKB2bEcDwjjzrernw+tBMdQvzMLktEapAaHUZcvOsA4F6QBjYbONTo/rxSA82PO8RT320ir8BGy2AfpgzrRD0/d7PLEpEapkaHEVdfexhxwAY5KeChib6kZjAMg/eX7mLCsl0A9GoZyIS7OuDpWqN/JYiISWr0bx5fTw/SDA98LFn2gc8URqQGyMm38tR3m/hpo/2S9geuasIzfVrg6KCOqiJijhodRvw9XThpeP8dRggzuySRcnU0PYcHZsQSdyAFJwcLrw1ow8DODc0uS0RquBodRvzcnTmFN41J1pDwUu1tT0pj5PQYDqVk4+vuzKf3RBLVtLbZZYmI1Oww4uvhzH7DPqqkLfN4zR6OVqq133Yk8+isDWTmWWkS4MmU4Z0JDfA0uywREaCGhxE/dxdOnQ4juWnH0DUEUt0YhsGUVQn87+ft2AyIalKbifd0xM9DszGJSOVRo8OIi5MD6Y6+AOSlHlUYkWol32rjhflbmb0uEYBBXUJ4pV8bnB3VBigilUuNDiMAOc5+UADWTPUZkeojNSufh2fGsmbPCSwWeP7GlozoEaqh3UWkUqrxYSTf1R8K7H1GRKqDhOOZjJi2nr3HM/FwceSDuyLo1SrQ7LJERC6oxoeRArdakKnJ8qR6iN5zgoe+iiU1O596vm5MHtaZVvV8zC5LROSianwYMdztlzY65SiMSNX2zfoDPDd3MwU2g/Yhfnw+NJK63m5mlyUickk1Pow4eNnDiGueJsuTqslqM3hz8Q4m/b4XgJvbBfP2He1xc3Y0uTIRkeKp8WHE2SsAAFdrJhTkgZMueZSqIzO3gLFfx7FkWzIAj10XxtjrwnDQ0O4iUoXU+DDi5l2bAsMBJ4sNsk+Cd5DZJYkUS1JqNiOmxbAtKQ0XJwfeur0d/TrUN7ssEZESq/FhxM/TlVN4UYc0yDyuMCJVwsYDKdw/I4aj6bkEeLkwaUgnIhv5m12WiMhlURjxsI/CWseSpvlppEr4eXMS//d1HLkFNsIDvZk8rBMhtTzMLktE5LIpjHjYJ8sDFEakUjMMg4+X7+btX3cC0DO8Dh8OisDbzdnkykRESqfGhxF/D2d2GqfHYVAYkUoqt8DKs99vZu6GQwDc270xz9/YEicN7S4i1UCNDyO+Z02Wp5l7pTI6kZHLg1/GErP/FI4OFl7+V2vu6dbI7LJERMqMwoi7MydPn6bJSzuGhoiSymRncjojpq/nwMlsvN2cmDg4kh5hAWaXJSJSpmp8GHFxciDz9My9BRman0Yqj5U7jzF65l+k5xbQqLYHU4Z1plldL7PLEhEpczU+jADkuvhBPtg0c69UEtPX7OPln7ZiM6BL41p8OiSSWp4akE9EqieFEaDAtRbkgyVLLSNirgKrjVcWbGNG9H4Abo9swGsD2uDqpKHdRaT6UhgBDI/akAFOOZqfptIxDPvNofp3LU7LyWfUzL/4Y9dxLBZ4pk8LHryqCRaLhnYXkepNYQSweNony3POO2X/4tMv/8oh+xR8ehV41YFhC8Cl+g7slXgiixHT17PraAbuzo68N7ADfdpoNGARqRmq/5+bxeDoab86wcmWB3mZJlcjhWKnQWoiHIqFRU+ZXU25Wb/vJP0/Wc2uoxkE+rjy7UNRCiIiUqMojABe3r5kG6c7B2rgs8qhIA/WTvr78YavYOMc8+opJ9/HHmTw52s5mZlH2/q+zB/Vgzb1fc0uS0SkQimMYB8S/mThkPDqxFopbJsH6UngFQhXPml/bsH/wbF4U8sqKzabwVu/7OCJbzeSZ7XRp3UQ3zwYRZCvRroRkZpHYYS/J8sDIOukucWIvd/Omg/t97vcD9c8B6FXQX4WfDsc8rJMLa+0svOsjJr1Fx8v3wPAIz2b8sngjri76IoZEamZFEYAP3dnThqaLK/S2LcKjmwCJ3foNAIcHOHWyeBZF45ug8XPmF3hZUtOy2HgZ9Es2nIEZ0cL79zRnqf7tMDBQZ2mRaTmUhjhn6dpFEZMF/2x/d8Od4NHLft970C47XPAAn/NgI1fm1be5dpyKJV+H61m08FU/D2cmTmyG7dFNjC7LBER05U4jPz+++/ccsst1KtXD4vFwrx58y66/IoVK7BYLOfcjhw5crk1l7kip2ky1WfEVMd3w85FAOxuMpSBk6LpO+EPXl2wjeX5rcjvcfqqmgX/B8d2mlhoyfyy9Qh3fBrNkbQcmtX1Yt6o7nQJrWV2WSIilUKJxxnJzMykffv23Hfffdx6663FXi8+Ph4fH5/Cx3Xr1i3prsuNn8ffp2lsWSfUXGSmPz8B4GDdq7lxVhJ5BTYAtielMWVVAm6OHfjOsz1t8jaSPWsILg8tx9G18o4/YhgGk37fyxuLd2AYcGVYAB/d3RFfd2ezSxMRqTRKHEb69u1L3759S7yjunXr4ufnV+L1KoKvuzOnTp+mKUg/jmYAMUnWSYy4WViAJw9eSZ7NxjXhdegfUZ/oPSf4Y9dxDqVkc2/ag/zs+ix1Tu3gu/H3sKzZ83RvFsCVYQE0qu1p9rsolFdg4/m5m/k29iAAQ7o14sVbWuHkqLgrInK2ChuBtUOHDuTm5tKmTRteeuklunfvfsFlc3Nzyc3NLXyclpZWrrU5OzqQ7eQHgFUz95rmyLJPCCrIZrOtMTG04rkbWzCyRxMcHCz061AfwzDYfyKLP3Yf56uN/2ZM0jPczjJ+39aCf2+xf55CarnTo1kdejQLoHuz2vh5mBMtT2Xm8eBXsaxLOImDBV68pTXDrmhsSi0iIpVduYeR4OBgPv30Uzp16kRubi6TJ0+mZ8+erF27lo4dO553nfHjx/Pyyy+Xd2lF5LnWgjxNlmcGwzD4ctVO+sR8Dhb4wbU/34y8go4N/YssZ7FYaBzgSeMAT+j2ILZlx+CPt3jbfSqWWhH8fNiTAyezmb0ukdnrErFYoG19X3o0C6BHswAiG/tXyIRzu49mMGL6evafyMLL1YmP7o6gZ3jlOS0pIlLZWAzDMC57ZYuFuXPn0r9//xKtd/XVV9OwYUO+/PLL875+vpaRkJAQUlNTi/Q7KUuj3vuSj1NHk+fqj8u4feWyDzlXalY+T3+/Ea/t3/KOy6eccgzAYewmfL2LcbrFZoUZ/WDfHxDYlswhi1h3MJs/dh1n1e5j7EzOKLK4m7MDXUJr06NZbXo0q0OLIO8yv6R29e7jPPxVLGk5BTTwd2fKsM6EB3mX6T5ERKqKtLQ0fH19L/n9bcpEeV26dGHVqlUXfN3V1RVXV9cKrAgsngGQCs65qfYvOQcNQFXeNiSe4tHZGzh4KotFrj8D4NdzNJbiBBE4Pf7I5/BpD0jejOfy/3DNLe9zTQt7K0RyWg6rdx9n1a7j/LH7OMfSc/l95zF+33kM2EGAlwtXNA2gR5i9v0mwr3up3s/Mtft5Yf5WrDaDyEb+TBoSSYBXxX6ORUSqIlPCSFxcHMHBwWbs+oIcTs/ca8EGOal/j28hZc4wDKasSuD1RTsosBkM8N1Ny9xEcPbE0ml4yTbmEwy3fgZf3QaxU6FxD2h7OwCBPm7c2rEBt3ZsgGEY7EzOYNXu46zadYy1CSc5npHHjxsP8+PGwwA0reNpP6UTVoduTWrh7Va8K16sNoPXFm7ni9UJAPTvUI/Xb2uHm7MCrYhIcZQ4jGRkZLB79+7CxwkJCcTFxVGrVi0aNmzIuHHjOHToEDNmzADg/fffJzQ0lNatW5OTk8PkyZP57bff+PXXX8vuXZQBH093Ug0PfC1Z9rFGFEbKxanMPJ78diPLdhwF4Ma2Qbxp/QL2AhH3gLv/xTdwPs2ugyufgD/ehp/GQL0IqN20yCIWi4XwIG/Cg7wZ0SOUvAIbGxJPsWr3cf7YdZxNB1PYcyyTPccymR69H0cHCxEhfoVX6bQP8cP5PFfBZOQW8NjsDfx2+v08cX1zRl/bDItFI6qKiBRXicNITEwM11xzTeHjxx9/HIBhw4Yxbdo0kpKSSExMLHw9Ly+PJ554gkOHDuHh4UG7du1YunRpkW1UBv4eLpw0vO1hRKOwlouYfSd5bPYGDqfm4OLkwH9ubsU9TbKwfLIUsEC3hy5/4z3HQWI07F8N3wyDkUvB+cKTzrk4OdC1SW26NqnNEzeEk5qVT/TeE6zafYxVu46z70QWMftPEbP/FBOW7cLL1YluTU73NwmrQ9M6nhxKyWbk9Bh2HEnH1cmBd+/swE3tKleLn4hIVVCqDqwVpbgdYEpj8h976bjkDjo67IaBM6HlzeWyn5rIZjP49Pc9vPPrTqw2g9AATz66O4LW9Xzhx8fgr+nQ4ma4a2bpdpSWZO8/knXcPqfNze9e9qYOnMxi9W57X5M1u49zKiu/yOvBvm7kFtg4mZlHHW9XPh/aiQ4hfqWrX0SkmqnUHVgrI7/TLSOAWkbK0PGMXB7/ZuPpTqPQr0M9XhvQFi9XJ/vpsI1z7AtGjS79znyC4dZJ9v4jMVPs/UfaFH+U4LOF1PLgri4NuatLQ2w2g62H0+z9TXYfY/2+UySl5gDQMtiHKcM6Uc+vdJ1fRURqMoWR0/w9nDlhnE5tGmukTPy59wSPzd7A0fRcXJ0ceKVfa+7sFPJ3f4r1k8GaC/U6QsNuZbPTZr1O9x95x97qEtz+nP4jJeXgYKFtA1/aNvDl4Z5Nyc6zErP/JPtPZDEgoj6ervoxEhEpDf0WPc3Pw5ldhTP3njS3mCrOajP4ePlu3l+6E5sBzep68fHdHYuOt5GfA+s+t9+/YjSUZYfPns/B/mhIXAPfDocRSy7af6Sk3F0cuTKsDleGldkmRURqNE2ScZqv+1kz9+o0zWU7mp7D0C/W8u4SexC5rWMDfhzd/dyBvzZ/Y2+B8g2Blv3KtghHJ7htMnjUhiOb4Nd/l+32RUSkTCmMnObv4czJ0y0jtkydprkcq3Yd58YJq1i9+wTuzo68fUd73rmzPR4u/2iAMwyI/th+v+uD9vBQ1nzrw4DP7PfXfw5b55X9PkREpEwojJzm6+5c2DJi02R5JVJgtfHOr/EM+WItxzNyCQ/05qdHu3N7ZIPzr7BnGRzbAS5e0HFo+RUW1gt6/J/9/o+Pwsm95bcvERG5bAojpzk5OpDjbB9wy9BpmmI7kprD3ZPX8uFvuzEMGNQlhPmju9Os7kXmY1nzkf3fjkPBzbd8C7zm3xDSDXLT7P1HCnIvuYqIiFQshZGzZLsHAuCUfkidWIthRfxRbvzgD9YlnMTTxZEJd3Vg/K2XGAY9eSvsXQ4WB+haikHOisvRCW6fAu61IGmj+o+IiFRCCiNnyfOsx3ZbQyyGFXYsNLucSivfauP1RTsYPnU9JzPzaBXsw4LHrqRfh/qXXjn6E/u/Lf8F/o3Kt9AzfBvAgEn2++s+g23zK2a/IiJSLAojZ/HzcGahtav9wbZ5ptZSWR1Kyeauz/7k05V7ABjSrRE/PHIFoQHFmGk3Pdl+FQ2UzSBnJdH8Bug+xn5//mg4mVCx+xcRkQtSGDmLn4cLi2xd7A/2rtCpmn9Ysi2ZGyf8Qez+U3i7OvHJ4I682r9N8WenXT8ZrHnQoAuEdC7fYs/n2v9ASFd7/5Hv7lX/ERGRSkJh5Cx+7s7sMepz3KMp2Aog/mezS6oU8gpsvLpgG/fPiCE1O592DXxZ+NiV3Ni2BJPC5WXZwwhA1KjyKfRSHJ3h9i/sMwMf3gBLXjCnDhGRysQwIDcdrAWmlaARWM/i5+EMwEafnlyXtcfetyDiHpOrMteBk1mMnr2BjQdSALiveyjP9m2Bi1MJc+ymOZB9EvwaQctbyr7Q4vJtAP0/hdkDYe2n9vlrzKxHRKSsWAsgJxWyT/19y0kp+jj7FGSnnLuMrQAeWgVBbU0pXWHkLH4eLgBEu13FdUyBPcvt/2nufqbWZZbFW5J46rtNpOcU4OPmxNt3tOeG1kEl35DN9nfH1W4Pg0MxT+uUl/A+cMWjsOZDmDfK/sPn39jcmkREwN5KkZ99gSDxz8dnL5NiPwVdGtmnSl//ZVIYOYufu71lJN4aDHVawrHt9lM1He42ubKKlZNvZfzP25kevR+AiIZ+fDgoggb+Hpe3wd1L4MQucPWpPC1N170IiX/CwfXw7b1w3y/g5GJ2VSJSXdhskJt6gSCRcvEWC2sp+7O5+tj/iHb3t9/czrpfePvHc25+4Gze7OMKI2fx97SHkZSsfGjbH1Zstw8jXoPCyL7jmYya9RdbD9sT9oNXN+HJG8JxdixF96Lo04OcRQ4D14sMhlaRHJ3h9qnwaQ84/BcsfRH6jDe7KhGpbApyz39a45ItFqmAcfn7dXC6QJDwOzdEFHnsWz5TbJSzqldxOfJ1t/9lnJKdB636w4rxsOe3GnOq5qeNhxn3w2Yycgvw93Dm3Ts7cE2LuqXbaNImSPgdLI7Q5cGyKbSs+IXAgE9h9l3w5yfQqDu0vNnsqkSkrJ3poFmcIPHPPhf5WaXbt7PnP0KE30VaLM567OJVtrOZV3IKI2ep52efZv7gqWx22OrRok4L+xwq8YugwyCTqys/OflWXv5pG7PXJQLQubE/HwyKINi3DJrszkyI13qA/cu/sgnvax/zJPojmP/I6f4jFTQYm4iUjLWgaB+JC/afOE/QMKyl2LGl+CGiyDJ+4ORa2nddIyiMnCXY150b2wbx8+YjvLk4ni9a9YOVO+wDoFXTMJJbYOX+GTH8ses4FguM6tmMsb3CcCrNaZkz0g7Dlu/s9826nLc4er1k7z9yKMY+/si9i9V/RKS8nN1Bs9hXfKTYlyltB01H14v0nfC7cP8KVx9w0EgY5Ulh5B+e6t2CX7cm89uOo2xs25P2vGE/VZOTWv6TulWwAquNMbPj+GPXcTxcHJl4TyRXN69TdjtY97n9crGGV0D9jmW33bLm6Ax3nO4/cigWlr0MvV8zuyqRys1mtf9eLBIkUs4NEucLGta80u3b1RfcfS/QGnGRFgsTO2jKxSmM/ENogCeDujTkyz/380K0jXkBzbEc3wnxi6H9QLPLKzM2m8GzP2xm8dYjuDg68PnQTnRvFlB2O8jLhJgv7Pcrc6vIGX4Nof9EmHO3/ZRNo+7Q4kazqxIpfwW5Fx9/4kKtFWXVQbNEpz2qbgdNuTj9j57HY9eF8f1fB9l4MJXdHXoRdnyn/VRNNQkjhmHw34Xb+S72IA4W+GBQRNkGEYC4WfZfZP6h9n4ZVUGLm6DbKPjzY5j3MDz0hz2kiFR2Z3fQvGSQSCn6epl10DwTIvyK11pRwzpoysUpjJxHHW9X7r+yCROW7eJ/+5szFWD3MshJAzcfs8srtQ+W7eaL1faJ4t68vT192lzGQGYXY7Par04Be6uI2YOclUSvl+DAn/bTNd/eC/cuUv8RqTjW/HOv5ihui0VpOmhaHOwtDhcNERe4rFQ/H1IGFEYu4P6rmjBz7X6Wn6pDau3G+Gbug52Lod2dZpdWKlNXJ/De0p0AvHBzK26PbFD2O9m5GE7utf+iqmpjtDi52OevmXSVvUOr+o9ISRmGvbWhRFd8nH4uL710+3Z0BY9aF2iN8Ltw0FAHTTGZwsgFeLk6Mea6MP4zfyvfZkUykn32AdCqcBj5LvYgL/+0DYCxvcK4r0do+ezozOW8ne4FF8/y2Ud58m8M/T6Brwfb+4807lF1TjVJ2TnTQbPwao7zXOFxoaBRJh00/S7ed+J8rRXqoClVlMLIRdzVpSFTViXw3cnOjHT9HnYvrbKnan7ZeoRnvt8EwL3dGzPmurDy2dGhv2D/anvntC4PlM8+KkLLm6Hrw7B2Isx9yD6BVGUcJ0UuLT+nGANd/fP1U/af9bLqoFmS0x7qoCk1kD7xF+Hs6MBTvVswalYm+4xgGluTYOcv0O4Os0srkdW7j/PorA1YbQa3RzbgPze1wlJeHcfOtIq0uQ186pXPPirK9a/AgbX24eK/uw/u/dl+GbBUPMOwjzFxySCRcu7rBdml2/c5HTSLMSS3u7+9VVAdNEWKRWHkEm5sG0T7EH8WJHVhtNN8+1U1VSiMbEg8xf0zYsiz2ujdOpDXb22Lg0M5/YJMPQhb59rvV4XLeS/FyeX0+CNXwcF1sOwVuOFVs6uq2qz5517NUawWi5Qy6KDpV/LTHuqgKVIhFEYuwWKxMK5vC175vCujneZj27UEh9z0yjPh20XsOJLG8Knrycqz0qNZAB8MiiibkVUvZO0k+xdG4yshuH357aci+TeGfh/BN0NgzQf2/iPNe5tdlbmKdNAsRv+Js1ssyqqD5jkhwu/iQUMdNEUqNYWRYujWpDZBzTuTkBBIKMn2UzVtbze7rIvafyKTIVPWkZqdT0RDPyYNicTVqRwvsc1Nh9jp9vtRo8tvP2Zo9S/7JH/rJsHcB+39R3zL4SqkivbPDprFHpr7FNjyS7fvwg6aF+k7cb7+FeqgKVItKYwU09N9W/DzR10Z5fAjKTHf4leJw0hyWg73TFnLsfRcWgR5M3V4Zzxdy/m/esNMyE2F2mEQdkP57ssMN7xq7z+SFGfvPzJ8YeXpP5KfU7wRM/+5TE5q6fb7zw6axR1N09VHHTRFpAj9RiimFkE+/Nz8X7D3Rzz2/4aRm46lEp6qOZWZxz2T13LgZDaNanswY0QX/DzK+Zx3kUHOHqmezeFOrnDHNPv4IwfWwm//hetfLrvtF+mgebH+EynnLlPaDpouXmeFCL/iDcmtDpoiUoYURkrgrltuYv/7gTSyJLNl5be0ueE+s0sqIiO3gOFT17HraAaBPq58NaIrdb3dyn/HOxZAyn5wrwXt7ir//ZmlVij860P4dhisft8+f03zf7QCnd1Bs7hzfJxppSiTDprFPe3hpw6aIlJpKIyUQD1/D9bU602jpBmcWvct1l734lheV6aUUE6+lfunx7DxYCr+Hs58NaIrIbU8KmbnZy7n7TwCXCpon2Zp3R/23Q/rP4cfRkJQu6ItFnkZpdu+k9sFQoTfxYOGOmiKSBWmMFJCba8fBjNm0Ck/hvnrdnJrt3CzSyLfamP0rA1E7z2Bl6sT0+/rQlhgBZ1COrDeftrC0QU6318x+zTbDf+1X+qbtBH2/XH+Zdx8L9134nyvq4OmiNRACiMl5B0aSZpbfXxyDhG79BtujHwON2fzJoKz2Qye/m4TS7cn4+LkwOdDO9GugV/FFRD9kf3ftneAd2DF7ddMzm4w+Hv76akz/S2KjE3hW7UmBxQRMZnCSElZLHhE3AbRHxCV+wfT1uzjoaubmlKKYRi89NNW5m44hKODhU/u7khU09oVV8Cp/bD9R/v96jDIWUl41bHPvSMiIqWmk8yXwanNAACudYjji+VbSckq5aRYl+ndJTuZEb0fiwXeuaM9vVpVcMvE2klg2KDJNRDYumL3LSIi1YbCyOWoF4Hh1xAPSy6ReTF8smJPhZcw+Y+9fPjbbgBe+Vdr+kfUr9gCclLhrxn2+9VtkDMREalQCiOXw2LB0qo/ADc5rmXamn0cSinlWA8l8PX6RP67cDsAT/UOZ0hU4wrbd6G/vrQP7V2nBTS7ruL3LyIi1YbCyOVq3R+AXk4bcCjI5p1f4ytktz9vTmLcD5sBeOCqJjzS04T+KtYCWPup/X63RzTwlYiIlIrCyOWq1xH8GuJm5NLTYSNzNxxi2+G0ct3lyp3HGDNnAzYD7uocwri+LbCYEQS2z4fUA+ARAO0GVvz+RUSkWlEYuVwWC7TqB8DI2hsxDHjzlx3ltrvY/Sd56MtY8q0GN7UN5rUBbc0JIoYBa05fztvlfvtlriIiIqWgMFIap/uNdMxZi6dDPivij7Fmz/Ey3822w2kMn7qe7HwrVzevw3sDO5g38uuBtXD4L/tU7p1GmFODiIhUKwojpVE/EnxDcMjP4t/hhwB4fdEObDajzHax91gGQ79YS3pOAZ0a+fPpPZG4OJn433ZmkLP2A+1jbYiIiJSSwkhpnHWqZoBrDJ4ujmw6mMrCzUllsvnDKdkMmbKO4xl5tAr2Ycrwzri7mDiy58m9sH2B/X63GjbImYiIlBuFkdI6farGbe+vPNLDPtbHW7/Ek1dgK9VmT2Tkcs+UtRxKyaZJgCczRnTB1925tNWWzp+fAgY0ux7qtjC3FhERqTZKHEZ+//13brnlFurVq4fFYmHevHmXXGfFihV07NgRV1dXmjVrxrRp0y6j1EqqQSfwaQB5GYwI3kuAlyuJJ7OYvS7xsjeZlpPPsKnr2Hssk3q+bnw5sisBXq5lWPRlyD4FG76y369pQ7+LiEi5KnEYyczMpH379nz88cfFWj4hIYGbbrqJa665hri4OMaOHcvIkSP55ZdfSlxspXTWqRq3nT8xtlcYAB8s20V6Tn6JN5edZ2XktBi2HEqjtqcLX47sSn2/SjCTa+x0yM+Euq2hSU+zqxERkWqkxGGkb9++/Pe//2XAgAHFWv7TTz8lNDSUd955h5YtWzJ69Ghuv/123nvvvRIXW2mdHgCN+MUMjKhDkwBPTmTm8fnve0u0mbwCG4/MjGXdvpN4uzox/b4uNK3jVfb1lpQ13z4PDdhbRTTImYiIlKFy7zMSHR1Nr169ijzXu3dvoqOjL7hObm4uaWlpRW6VWv1O4FMf8tJxTljO033CAfj8jwSOpuUUaxNWm8ET325kefwx3JwdmDK8M23q+5Zn1cW3dS6kHwavQGh7u9nViIhINVPuYeTIkSMEBhadTTYwMJC0tDSys88/n8v48ePx9fUtvIWEhJR3maXj4AAt/2W/v20+vVsHEdHQj+x8KxOW7brk6oZh8J/5W/hp42GcHCxMvCeSLqG1yrnoYjKMvy/n7XI/OJncd0VERKqdSnk1zbhx40hNTS28HThwwOySLq3wVM0iLNY8xvVtCcCc9QfYcyzjoqu+sTieWWsTsVjgvYEduCa8bjkXWwL7V0PSRnByh8j7zK5GRESqoXIPI0FBQSQnJxd5Ljk5GR8fH9zdz98x09XVFR8fnyK3Sq9BF/CuB7lpsOc3uoTWolfLulhtBm8tvvAkehNX7OHTlXsA+N+AttzSvl5FVVw80ac7KncYBJ61za1FRESqpXIPI1FRUSxbtqzIc0uWLCEqKqq8d12xHByg1elTNVvnAfB0nxY4WGDx1iPE7j91zioz1+7njcX2+WzG9W3BoC4NK6ra4jm+G+IX2e93e8TcWkREpNoqcRjJyMggLi6OuLg4wH7pblxcHImJ9nE1xo0bx9ChQwuXf+ihh9i7dy9PP/00O3bs4JNPPuGbb77h//7v/8rmHVQmpwdAI/5nKMileaA3t0c2AOCNRTswjL+Hif9x42H+PW8LAI/0bMqDVzet6Gov7c9PAAOa94GAMLOrERGRaqrEYSQmJoaIiAgiIiIAePzxx4mIiOCFF14AICkpqTCYAISGhrJw4UKWLFlC+/bteeedd5g8eTK9e/cuo7dQiYR0Be/g06dqlgPwf9c3x9XJgXX7TrJs+1EAlu84yuNfx2EYcE+3hjzVO9zMqs8v6yTEzbLfjxptbi0iIlKtOZV0hZ49exb5C/+fzje6as+ePdmwYUNJd1X1nLmqZt0k2DYPwvsQ7OvOfT1CmbhiD28s3oGHqyMPfRVLgc2gX4d6vPKvNlgq47gdMV9AQTYEtYPGPcyuRkREqrFKeTVNlXbmqpod9lM1AA9d3RQ/D2d2Hc1gyJR15BbYuK5FXd6+oz0ODpUwiBTkwrrP7PejRmuQMxERKVcKI2UtpKt9cLDcVNi7EgBfd2dGX9MMsA9u1jW0Fh8P7oizYyU9/Ft+gIxk+ymn1sUbaVdERORyVdJvwyrMwfGsAdDmFT49JKoRUU1q06NZAJOHdcLN2dGc+i6lyCBnD4CTi7n1iIhItVfiPiNSDK37w/rPYccCKHgfnFxwdXJk9gPdzK7s0hJWQvIWcPaATveaXY2IiNQAahkpDw2jwLMu5KTav9yrkjODnEXcA+7+5tYiIiI1gsJIeXBwPGcAtCrhWDzs+hWwQNeHzK5GRERqCIWR8nJmALQdC8Cab2opxfbnJ/Z/W9wEtSvhIGwiIlItKYyUl0ZXnD5Vk1J4VU2llnkcNs6x348aZW4tIiJSoyiMlBcHR2h5i/3+trnm1lIc66dAQQ7U62jv8yIiIlJBFEbKU6t+9n93LKzcp2ryc+xX/4C9VUSDnImISAXSpb3lqVF38AiArOOQ8Ds0u87siv6WdRKSNsLhDfbaMo+BT4O/A5SIiEgFURgpT45O9lM1sVPtA6CZFUayT8HhOEiKs/97eAOk7D93uStGg6NzBRcnIiI1ncJIeWvd3x5Gti+Am94t/y/77FN/t3icCSCn9p1/Wf9QqNcBgjtAg07Q8IryrU1EROQ8FEbKW6Me4FEbsk7Avj+g6bVlt+3C4BF3utVjw0WCR2N76KgXcTqAtNegZiIiUikojJS3wlM10+wDoF1uGMlO+bvF48zpllMJ51/Wr9FZoaODPXh41Lq8/YqIiJQzhZGK0Kq/PYzsOHOq5hKH/UzwOLuPx0WDR4e/Wz0UPEREpIpRGKkIja/8+1TN/lXQpOffr+Wk/n2q5Uyrx8m959+OX8N/nGrpoOAhIiJVnsJIRXB0ghY3w1/TYe1nkLTp71aPk3vOv45vQ3vgOLvVQ8FDRESqIYWRitKqnz2MxC+0387m2xDqtT/rVEsH8KxtRpUiIiIVTmGkooReBU2ugRN7ILjdWadaIhQ8RESkRlMYqSiOzjB0ntlViIiIVDqam0ZERERMpTAiIiIiplIYEREREVMpjIiIiIipFEZERETEVAojIiIiYiqFERERETGVwoiIiIiYSmFERERETKUwIiIiIqZSGBERERFTKYyIiIiIqRRGRERExFQKIyIiImIqJ7MLKA7DMABIS0szuRIREREprjPf22e+xy+kSoSR9PR0AEJCQkyuREREREoqPT0dX1/fC75uMS4VVyoBm83G4cOH8fb2xmKxXHTZtLQ0QkJCOHDgAD4+PhVUYdWmY1ZyOmYlp2NWcjpmJadjVnLlecwMwyA9PZ169erh4HDhniFVomXEwcGBBg0alGgdHx8ffRBLSMes5HTMSk7HrOR0zEpOx6zkyuuYXaxF5Ax1YBURERFTKYyIiIiIqapdGHF1deXFF1/E1dXV7FKqDB2zktMxKzkds5LTMSs5HbOSqwzHrEp0YBUREZHqq9q1jIiIiEjVojAiIiIiplIYEREREVMpjIiIiIipql0Y+fjjj2ncuDFubm507dqVdevWmV2SKV566SUsFkuRW4sWLQpfz8nJYdSoUdSuXRsvLy9uu+02kpOTi2wjMTGRm266CQ8PD+rWrctTTz1FQUFBRb+VcvP7779zyy23UK9ePSwWC/PmzSvyumEYvPDCCwQHB+Pu7k6vXr3YtWtXkWVOnjzJ4MGD8fHxwc/PjxEjRpCRkVFkmU2bNnHllVfi5uZGSEgIb775Znm/tXJzqWM2fPjwcz53ffr0KbJMTTtm48ePp3Pnznh7e1O3bl369+9PfHx8kWXK6udxxYoVdOzYEVdXV5o1a8a0adPK++2Vi+Ics549e57zWXvooYeKLFOTjtnEiRNp165d4cBlUVFRLFq0qPD1Sv8ZM6qROXPmGC4uLsYXX3xhbN261bj//vsNPz8/Izk52ezSKtyLL75otG7d2khKSiq8HTt2rPD1hx56yAgJCTGWLVtmxMTEGN26dTOuuOKKwtcLCgqMNm3aGL169TI2bNhg/Pzzz0ZAQIAxbtw4M95Oufj555+N559/3vjhhx8MwJg7d26R119//XXD19fXmDdvnrFx40bjX//6lxEaGmpkZ2cXLtOnTx+jffv2xp9//mn88ccfRrNmzYxBgwYVvp6ammoEBgYagwcPNrZs2WLMnj3bcHd3NyZNmlRRb7NMXeqYDRs2zOjTp0+Rz93JkyeLLFPTjlnv3r2NqVOnGlu2bDHi4uKMG2+80WjYsKGRkZFRuExZ/Dzu3bvX8PDwMB5//HFj27Ztxocffmg4OjoaixcvrtD3WxaKc8yuvvpq4/777y/yWUtNTS18vaYdsx9//NFYuHChsXPnTiM+Pt547rnnDGdnZ2PLli2GYVT+z1i1CiNdunQxRo0aVfjYarUa9erVM8aPH29iVeZ48cUXjfbt25/3tZSUFMPZ2dn49ttvC5/bvn27ARjR0dGGYdi/dBwcHIwjR44ULjNx4kTDx8fHyM3NLdfazfDPL1abzWYEBQUZb731VuFzKSkphqurqzF79mzDMAxj27ZtBmCsX7++cJlFixYZFovFOHTokGEYhvHJJ58Y/v7+RY7ZM888Y4SHh5fzOyp/Fwoj/fr1u+A6Nf2YGYZhHD161ACMlStXGoZRdj+PTz/9tNG6desi+xo4cKDRu3fv8n5L5e6fx8ww7GFkzJgxF1ynph8zwzAMf39/Y/LkyVXiM1ZtTtPk5eURGxtLr169Cp9zcHCgV69eREdHm1iZeXbt2kW9evVo0qQJgwcPJjExEYDY2Fjy8/OLHKsWLVrQsGHDwmMVHR1N27ZtCQwMLFymd+/epKWlsXXr1op9IyZISEjgyJEjRY6Rr68vXbt2LXKM/Pz86NSpU+EyvXr1wsHBgbVr1xYuc9VVV+Hi4lK4TO/evYmPj+fUqVMV9G4q1ooVK6hbty7h4eE8/PDDnDhxovA1HTNITU0FoFatWkDZ/TxGR0cX2caZZarD779/HrMzZs6cSUBAAG3atGHcuHFkZWUVvlaTj5nVamXOnDlkZmYSFRVVJT5jVWKivOI4fvw4Vqu1yIEECAwMZMeOHSZVZZ6uXbsybdo0wsPDSUpK4uWXX+bKK69ky5YtHDlyBBcXF/z8/IqsExgYyJEjRwA4cuTIeY/lmdequzPv8XzH4OxjVLdu3SKvOzk5UatWrSLLhIaGnrONM6/5+/uXS/1m6dOnD7feeiuhoaHs2bOH5557jr59+xIdHY2jo2ONP2Y2m42xY8fSvXt32rRpA1BmP48XWiYtLY3s7Gzc3d3L4y2Vu/MdM4C7776bRo0aUa9ePTZt2sQzzzxDfHw8P/zwA1Azj9nmzZuJiooiJycHLy8v5s6dS6tWrYiLi6v0n7FqE0akqL59+xbeb9euHV27dqVRo0Z88803Ve4HTKqOu+66q/B+27ZtadeuHU2bNmXFihVcd911JlZWOYwaNYotW7awatUqs0upMi50zB544IHC+23btiU4OJjrrruOPXv20LRp04ous1IIDw8nLi6O1NRUvvvuO4YNG8bKlSvNLqtYqs1pmoCAABwdHc/pHZycnExQUJBJVVUefn5+NG/enN27dxMUFEReXh4pKSlFljn7WAUFBZ33WJ55rbo78x4v9nkKCgri6NGjRV4vKCjg5MmTOo6nNWnShICAAHbv3g3U7GM2evRoFixYwPLly2nQoEHh82X183ihZXx8fKrsHyAXOmbn07VrV4Ain7WadsxcXFxo1qwZkZGRjB8/nvbt2zNhwoQq8RmrNmHExcWFyMhIli1bVviczWZj2bJlREVFmVhZ5ZCRkcGePXsIDg4mMjISZ2fnIscqPj6exMTEwmMVFRXF5s2bi3xxLFmyBB8fH1q1alXh9Ve00NBQgoKCihyjtLQ01q5dW+QYpaSkEBsbW7jMb7/9hs1mK/zFGBUVxe+//05+fn7hMkuWLCE8PLxKn24oroMHD3LixAmCg4OBmnnMDMNg9OjRzJ07l99+++2cU1Bl9fMYFRVVZBtnlqmKv/8udczOJy4uDqDIZ60mHbPzsdls5ObmVo3PWKm7wFYic+bMMVxdXY1p06YZ27ZtMx544AHDz8+vSO/gmuKJJ54wVqxYYSQkJBirV682evXqZQQEBBhHjx41DMN+mVfDhg2N3377zYiJiTGioqKMqKiowvXPXOZ1ww03GHFxccbixYuNOnXqVKtLe9PT040NGzYYGzZsMADj3XffNTZs2GDs37/fMAz7pb1+fn7G/PnzjU2bNhn9+vU776W9ERERxtq1a41Vq1YZYWFhRS5TTUlJMQIDA40hQ4YYW7ZsMebMmWN4eHhU2ctUL3bM0tPTjSeffNKIjo42EhISjKVLlxodO3Y0wsLCjJycnMJt1LRj9vDDDxu+vr7GihUrilyGmpWVVbhMWfw8nrns8qmnnjK2b99ufPzxx1X2MtVLHbPdu3cbr7zyihETE2MkJCQY8+fPN5o0aWJcddVVhduoacfs2WefNVauXGkkJCQYmzZtMp599lnDYrEYv/76q2EYlf8zVq3CiGEYxocffmg0bNjQcHFxMbp06WL8+eefZpdkioEDBxrBwcGGi4uLUb9+fWPgwIHG7t27C1/Pzs42HnnkEcPf39/w8PAwBgwYYCQlJRXZxr59+4y+ffsa7u7uRkBAgPHEE08Y+fn5Ff1Wys3y5csN4JzbsGHDDMOwX977n//8xwgMDDRcXV2N6667zoiPjy+yjRMnThiDBg0yvLy8DB8fH+Pee+810tPTiyyzceNGo0ePHoarq6tRv3594/XXX6+ot1jmLnbMsrKyjBtuuMGoU6eO4ezsbDRq1Mi4//77z/ljoKYds/MdL8CYOnVq4TJl9fO4fPlyo0OHDoaLi4vRpEmTIvuoSi51zBITE42rrrrKqFWrluHq6mo0a9bMeOqpp4qMM2IYNeuY3XfffUajRo0MFxcXo06dOsZ1111XGEQMo/J/xiyGYRilb18RERERuTzVps+IiIiIVE0KIyIiImIqhRERERExlcKIiIiImEphREREREylMCIiIiKmUhgRERERUymMiIiIiKkURkSkRHr27MnYsWPNLqOQYRg88MAD1KpVC4vFUjhHydmmTZt2zvTp//TSSy/RoUOHcqlRRC5OYUREqrTFixczbdo0FixYQFJSEm3atLms7Tz55JPnTAImIhXDyewCRESsVisWiwUHh5L/fXRmNuorrriiVDV4eXnh5eVVqm2IyOVRy4hIFdSzZ08ee+wxnn76aWrVqkVQUBAvvfRS4ev79u0755RFSkoKFouFFStWALBixQosFgu//PILERERuLu7c+2113L06FEWLVpEy5Yt8fHx4e677yYrK6vI/gsKChg9ejS+vr4EBATwn//8h7OnucrNzeXJJ5+kfv36eHp60rVr18L9wt+nTX788UdatWqFq6sriYmJ532vK1eupEuXLri6uhIcHMyzzz5LQUEBAMOHD+fRRx8lMTERi8VC48aNL3rc5s2bR1hYGG5ubvTu3ZsDBw4UvvbP0zTDhw+nf//+vP322wQHB1O7dm1GjRpFfn5+4TKffPJJ4fYCAwO5/fbbL7p/ETk/hRGRKmr69Ol4enqydu1a3nzzTV555RWWLFlS4u289NJLfPTRR6xZs4YDBw5w55138v777zNr1iwWLlzIr7/+yocffnjOvp2cnFi3bh0TJkzg3XffZfLkyYWvjx49mujoaObMmcOmTZu444476NOnD7t27SpcJisrizfeeIPJkyezdetW6tate05thw4d4sYbb6Rz585s3LiRiRMnMmXKFP773/8CMGHCBF555RUaNGhAUlIS69evv+D7zMrK4rXXXmPGjBmsXr2alJQU7rrrrosem+XLl7Nnzx6WL1/O9OnTmTZtGtOmTQMgJiaGxx57jFdeeYX4+HgWL17MVVdddcnjLSLnUSZz/4pIhbr66quNHj16FHmuc+fOxjPPPGMYhmEkJCQYgLFhw4bC10+dOmUAxvLlyw3DsE8FDhhLly4tXGb8+PEGYOzZs6fwuQcffNDo3bt3kX23bNnSsNlshc8988wzRsuWLQ3DMIz9+/cbjo6OxqFDh4rUd9111xnjxo0zDMMwpk6dagBGXFzcRd/nc889Z4SHhxfZ18cff2x4eXkZVqvVMAzDeO+994xGjRpddDtn9vfnn38WPrd9+3YDMNauXWsYhmG8+OKLRvv27QtfHzZsmNGoUSOjoKCg8Lk77rjDGDhwoGEYhvH9998bPj4+Rlpa2kX3LSKXppYRkSqqXbt2RR4HBwdz9OjRUm0nMDAQDw8PmjRpUuS5f263W7duWCyWwsdRUVHs2rULq9XK5s2bsVqtNG/evLAfhpeXFytXrmTPnj2F67i4uJzzHv5p+/btREVFFdlX9+7dycjI4ODBgyV6n05OTnTu3LnwcYsWLfDz82P79u0XXKd169Y4OjoWPj77GF9//fU0atSIJk2aMGTIEGbOnHnO6SwRKR51YBWpopydnYs8tlgs2Gw2gMKOoMZZ/TjO7utwoe1YLJaLbrc4MjIycHR0JDY2tsgXOVCkg6i7u3uRkFEZXexYeHt789dff7FixQp+/fVXXnjhBV566SXWr19/ycuIRaQotYyIVEN16tQBICkpqfC5842/cbnWrl1b5PGff/5JWFgYjo6OREREYLVaOXr0KM2aNStyCwoKKtF+WrZsSXR0dJFQtXr1ary9vWnQoEGJtlVQUEBMTEzh4/j4eFJSUmjZsmWJtnM2JycnevXqxZtvvsmmTZvYt28fv/3222VvT6SmUhgRqYbc3d3p1q0br7/+Otu3b2flypX8+9//LrPtJyYm8vjjjxMfH8/s2bP58MMPGTNmDADNmzdn8ODBDB06lB9++IGEhATWrVvH+PHjWbhwYYn288gjj3DgwAEeffRRduzYwfz583nxxRd5/PHHS3wZsLOzM48++ihr164lNjaW4cOH061bN7p06VKi7ZyxYMECPvjgA+Li4ti/fz8zZszAZrMRHh5+WdsTqcl0mkakmvriiy8YMWIEkZGRhIeH8+abb3LDDTeUybaHDh1KdnY2Xbp0wdHRkTFjxvDAAw8Uvj516lT++9//8sQTT3Do0CECAgLo1q0bN998c4n2U79+fX7++Weeeuop2rdvT61atRgxYsRlBSsPDw+eeeYZ7r77bg4dOsSVV17JlClTSrydM/z8/Pjhhx946aWXyMnJISwsjNmzZ9O6devL3qZITWUxzm7/FBEREalgOk0jIiIiplIYEREREVMpjIiIiIipFEZERETEVAojIiIiYiqFERERETGVwoiIiIiYSmFERERETKUwIiIiIqZSGBERERFTKYyIiIiIqf4frBiZX+ucV+gAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def average_error(num_bins, num_trials):\n", " hist_quantiles_meas = make_hist_quantiles(quart_alphas, max_contributions, epsilon, num_bins)\n", " tree_quantiles_meas = make_tree_quantiles(quart_alphas, b, max_contributions, epsilon, num_bins)\n", "\n", " def sample_error(meas):\n", " return np.linalg.norm(true_quantiles - meas(data))\n", " hist_err = np.mean([sample_error(hist_quantiles_meas) for _ in range(num_trials)])\n", " tree_err = np.mean([sample_error(tree_quantiles_meas) for _ in range(num_trials)])\n", "\n", " return num_bins, hist_err, tree_err\n", "\n", "import pandas as pd\n", "pd.DataFrame(\n", " [average_error(nb, num_trials=25) for nb in [70, 100, 250, 500, 750, 1_000, 3_000]],\n", " columns=[\"number of bins\", \"histogram error\", \"tree error\"]\n", ").plot(0); # type: ignore" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Privately Estimating the Distribution\n", "\n", "Minor note: instead of postprocessing the noisy counts into quantiles, they can be left as counts, which can be used to visualize the distribution.\n", "\n" ] }, { "cell_type": "code", "execution_count": 136, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvUklEQVR4nO3df1RVdb7/8Rc/BNQETZOjhmJlokmg8kO071iJUXEnyaaLjlfJ7LcayOSEjspKx0ErjVSujHXV5S3Tsate07KQ/DEliYJmlmm3H+KoB/Q6chQTHc7+/tHqNOeCxqFzOMB+PtY6K8/nvPfmvT9rEa/12fvs7WMYhiEAAAAT8fV2AwAAAI2NAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEzH39sNNEV2u10nT55Uu3bt5OPj4+12AABAPRiGofPnz6tr167y9b32Gg8BqA4nT55UWFiYt9sAAAANcPz4cd14443XrCEA1aFdu3aSfpjA4OBgL3cDAADqw2azKSwszPF3/FoIQHX48bRXcHAwAQgAgGamPpevcBE0AAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHX9vN2BG4Vlbfrbmu3nJjdAJAADmxAoQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHa8HoLy8PIWHhysoKEjx8fEqLi6+au3nn3+uhx56SOHh4fLx8VFubu4v3icAADAfrwagtWvXKjMzU9nZ2SotLVVUVJSSkpJUUVFRZ/3Fixd10003ad68ebJYLG7ZJwAAMB+vBqCFCxfq8ccf1/jx49W3b1/l5+erTZs2Wr58eZ31sbGxeumllzRq1CgFBga6ZZ+SVF1dLZvN5vQCAAAtl9cC0OXLl1VSUqLExMSfmvH1VWJiooqKihp1nzk5OQoJCXG8wsLCGvTzAQBA8+C1AHTmzBnV1NQoNDTUaTw0NFRWq7VR9zlt2jRVVlY6XsePH2/QzwcAAM0DD0OVFBgYeNVTagAAoOXx2gpQp06d5Ofnp/Lycqfx8vLyq17g7I19AgCAlsdrASggIEADBw5UYWGhY8xut6uwsFAJCQlNZp8AAKDl8eopsMzMTKWlpSkmJkZxcXHKzc1VVVWVxo8fL0kaN26cunXrppycHEk/XOT8xRdfOP594sQJHThwQNddd51uueWWeu0TAADAqwEoNTVVp0+f1qxZs2S1WhUdHa2tW7c6LmIuKyuTr+9Pi1QnT55U//79He9ffvllvfzyyxo6dKh27NhRr30CAAD4GIZheLuJpsZmsykkJESVlZUKDg52+/7Ds7b8bM1385Ld/nMBAGjJXPn77fVHYQAAADQ2AhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdrwegvLw8hYeHKygoSPHx8SouLr5m/bp16xQREaGgoCBFRkbq3Xffdfr8woULmjRpkm688Ua1bt1affv2VX5+vicPAQAANDNeDUBr165VZmamsrOzVVpaqqioKCUlJamioqLO+t27d2v06NGaMGGC9u/fr5SUFKWkpOjQoUOOmszMTG3dulVvvPGGDh8+rIyMDE2aNEmbNm1qrMMCAABNnI9hGIa3fnh8fLxiY2O1ZMkSSZLdbldYWJgmT56srKysWvWpqamqqqrS5s2bHWODBg1SdHS0Y5WnX79+Sk1N1cyZMx01AwcO1H333ac//vGPdfZRXV2t6upqx3ubzaawsDBVVlYqODjYLcf6z8KztvxszXfzkt3+cwEAaMlsNptCQkLq9ffbaytAly9fVklJiRITE39qxtdXiYmJKioqqnOboqIip3pJSkpKcqofPHiwNm3apBMnTsgwDG3fvl1Hjx7VPffcc9VecnJyFBIS4niFhYX9wqMDAABNmdcC0JkzZ1RTU6PQ0FCn8dDQUFmt1jq3sVqtP1u/ePFi9e3bVzfeeKMCAgJ07733Ki8vT7/61a+u2su0adNUWVnpeB0/fvwXHBkAAGjq/L3dgLstXrxYn3zyiTZt2qQePXpo165dmjhxorp27Vpr9ehHgYGBCgwMbOROAQCAt3gtAHXq1El+fn4qLy93Gi8vL5fFYqlzG4vFcs3677//XtOnT9eGDRuUnPzDNTS33367Dhw4oJdffvmqAQgAAJiL106BBQQEaODAgSosLHSM2e12FRYWKiEhoc5tEhISnOolqaCgwFF/5coVXblyRb6+zofl5+cnu93u5iMAAADNlVdPgWVmZiotLU0xMTGKi4tTbm6uqqqqNH78eEnSuHHj1K1bN+Xk5EiS0tPTNXToUC1YsEDJyclas2aN9u3bp2XLlkmSgoODNXToUE2dOlWtW7dWjx49tHPnTq1atUoLFy702nECAICmxasBKDU1VadPn9asWbNktVoVHR2trVu3Oi50Lisrc1rNGTx4sFavXq0ZM2Zo+vTp6tWrlzZu3Kh+/fo5atasWaNp06ZpzJgxOnv2rHr06KG5c+fqqaeeavTjAwAATZNX7wPUVLlyH4GG4D5AAAC4X7O4DxAAAIC3EIAAAIDpEIAAAIDptLgbIbYUXCcEAIDnsAIEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMx+UA9M0333iiDwAAgEbjcgC65ZZbdNddd+mNN97QpUuXPNETAACAR7kcgEpLS3X77bcrMzNTFotFTz75pIqLiz3RGwAAgEe4HICio6P16quv6uTJk1q+fLlOnTqlO+64Q/369dPChQt1+vRpT/QJAADgNg2+CNrf318jR47UunXrNH/+fP3P//yPnnvuOYWFhWncuHE6deqUO/sEAABwmwYHoH379umZZ55Rly5dtHDhQj333HP6+uuvVVBQoJMnT2rEiBHu7BMAAMBt/F3dYOHChVqxYoWOHDmi+++/X6tWrdL9998vX98fslTPnj21cuVKhYeHu7tXAAAAt3A5AC1dulSPPvqoHnnkEXXp0qXOms6dO+s//uM/fnFzAAAAnuByACooKFD37t0dKz4/MgxDx48fV/fu3RUQEKC0tDS3NQkAAOBOLl8DdPPNN+vMmTO1xs+ePauePXu6pSkAAABPcjkAGYZR5/iFCxcUFBT0ixsCAADwtHqfAsvMzJQk+fj4aNasWWrTpo3js5qaGu3Zs0fR0dFubxAAAMDd6h2A9u/fL+mHFaDPPvtMAQEBjs8CAgIUFRWl5557zv0dAgAAuFm9A9D27dslSePHj9err76q4OBgjzUFAADgSS5/C2zFihWe6AMAAKDR1CsAjRw5UitXrlRwcLBGjhx5zdr169e7pTEAAABPqVcACgkJkY+Pj+PfAAAAzVm9AtA/n/biFBgAAGjuXL4P0Pfff6+LFy863h87dky5ubn64IMP3NoYAACAp7gcgEaMGKFVq1ZJks6dO6e4uDgtWLBAI0aM0NKlS93eIAAAgLu5HIBKS0v1//7f/5Mkvf3227JYLDp27JhWrVqlRYsWub1BAAAAd3M5AF28eFHt2rWTJH3wwQcaOXKkfH19NWjQIB07dszlBvLy8hQeHq6goCDFx8eruLj4mvXr1q1TRESEgoKCFBkZqXfffbdWzeHDh/XAAw8oJCREbdu2VWxsrMrKylzuDQAAtEwuB6BbbrlFGzdu1PHjx/X+++/rnnvukSRVVFS4fHPEtWvXKjMzU9nZ2SotLVVUVJSSkpJUUVFRZ/3u3bs1evRoTZgwQfv371dKSopSUlJ06NAhR83XX3+tO+64QxEREdqxY4cOHjyomTNn8pwyAADg4GNc7emmV/H222/rt7/9rWpqajRs2DDHxc85OTnatWuX3nvvvXrvKz4+XrGxsVqyZIkkyW63KywsTJMnT1ZWVlat+tTUVFVVVWnz5s2OsUGDBik6Olr5+fmSpFGjRqlVq1b6z//8z3r3UV1drerqasd7m82msLAwVVZWeuSO1+FZW9yyn+/mJbtlPwAAtAQ2m00hISH1+vvt8grQb37zG5WVlWnfvn3aunWrY3zYsGF65ZVX6r2fy5cvq6SkRImJiT814+urxMREFRUV1blNUVGRU70kJSUlOertdru2bNmiW2+9VUlJSercubPi4+O1cePGa/aSk5OjkJAQxyssLKzexwEAAJoflwOQJFksFvXv31++vj9tHhcXp4iIiHrv48yZM6qpqVFoaKjTeGhoqKxWa53bWK3Wa9ZXVFTowoULmjdvnu6991598MEHevDBBzVy5Ejt3Lnzqr1MmzZNlZWVjtfx48frfRwAAKD5cflZYFVVVZo3b54KCwtVUVEhu93u9Pk333zjtuZc9WMvI0aM0JQpUyRJ0dHR2r17t/Lz8zV06NA6twsMDFRgYGCj9QkAALzL5QD02GOPaefOnRo7dqy6dOnieESGqzp16iQ/Pz+Vl5c7jZeXl8tisdS5jcViuWZ9p06d5O/vr759+zrV9OnTRx999FGD+gQAAC2PywHovffe05YtWzRkyJBf9IMDAgI0cOBAFRYWKiUlRdIPKziFhYWaNGlSndskJCSosLBQGRkZjrGCggIlJCQ49hkbG6sjR444bXf06FH16NHjF/ULAABaDpcDUIcOHXT99de75YdnZmYqLS1NMTExiouLU25urqqqqjR+/HhJ0rhx49StWzfl5ORIktLT0zV06FAtWLBAycnJWrNmjfbt26dly5Y59jl16lSlpqbqV7/6le666y5t3bpV77zzjnbs2OGWngEAQPPn8kXQc+bM0axZs5yeB9ZQqampevnllzVr1ixFR0frwIED2rp1q+NC57KyMp06dcpRP3jwYK1evVrLli1TVFSU3n77bW3cuFH9+vVz1Dz44IPKz8/Xiy++qMjISL3++uv6r//6L91xxx2/uF8AANAyuHwfoP79++vrr7+WYRgKDw9Xq1atnD4vLS11a4Pe4Mp9BBqC+wABAOB+rvz9dvkU2I/X6wAAADRXLgeg7OxsT/QBAADQaBp0I8Rz587p9ddf17Rp03T27FlJP5z6OnHihFubAwAA8ASXV4AOHjyoxMREhYSE6LvvvtPjjz+u66+/XuvXr1dZWZlWrVrliT4BAADcxuUVoMzMTD3yyCP66quvnJ6wfv/992vXrl1ubQ4AAMATXA5Ae/fu1ZNPPllrvFu3bld9hhcAAEBT4nIACgwMlM1mqzV+9OhR3XDDDW5pCgAAwJNcDkAPPPCAZs+erStXrkiSfHx8VFZWpueff14PPfSQ2xsEAABwN5cD0IIFC3ThwgV17txZ33//vYYOHapbbrlF7dq109y5cz3RIwAAgFu5/C2wkJAQFRQU6KOPPtLBgwd14cIFDRgwQImJiZ7oDwAAwO1cfhSGGTSXR2HUB4/LAACYhdsfhbFo0aJ6//Bnn3223rUAAADeUK8A9Morrzi9P336tC5evKj27dtL+uHO0G3atFHnzp0JQAAAoMmr10XQ3377reM1d+5cRUdH6/Dhwzp79qzOnj2rw4cPa8CAAZozZ46n+wUAAPjFXP4W2MyZM7V48WL17t3bMda7d2+98sormjFjhlubAwAA8ASXA9CpU6f0j3/8o9Z4TU2NysvL3dIUAACAJ7kcgIYNG6Ynn3xSpaWljrGSkhI9/fTTfBUeAAA0Cy4HoOXLl8tisSgmJkaBgYEKDAxUXFycQkND9frrr3uiRwAAALdy+UaIN9xwg95991199dVXOnz4sCQpIiJCt956q9ubAwAA8ASXA9CPevXqpV69ermzFwAAgEbh8ikwAACA5o4ABAAATIcABAAATIcABAAATMflABQeHq7Zs2errKzME/0AAAB4nMsBKCMjQ+vXr9dNN92k4cOHa82aNaqurvZEbwAAAB7RoAB04MABFRcXq0+fPpo8ebK6dOmiSZMmOd0dGgAAoKlq8DVAAwYM0KJFi3Ty5EllZ2fr9ddfV2xsrKKjo7V8+XIZhuHOPgEAANymwTdCvHLlijZs2KAVK1aooKBAgwYN0oQJE/S3v/1N06dP17Zt27R69Wp39goAAOAWLgeg0tJSrVixQm+99ZZ8fX01btw4vfLKK4qIiHDUPPjgg4qNjXVrowAAAO7icgCKjY3V8OHDtXTpUqWkpKhVq1a1anr27KlRo0a5pUEAAAB3czkAffPNN+rRo8c1a9q2basVK1Y0uCkAAABPcvkiaF9fX/3tb39zvC8uLlZGRoaWLVvm1sYAAAA8xeUA9Nvf/lbbt2+XJFmtVg0fPlzFxcX6wx/+oNmzZ7u9QQAAAHdzOQAdOnRIcXFxkqS//OUv6tevn3bv3q0333xTK1eudHd/AAAAbudyALpy5YoCAwMlSdu2bdMDDzwgSYqIiNCpU6fc2x0AAIAHuByAbrvtNuXn5+uvf/2rCgoKdO+990qSTp48qY4dO7q9QQAAAHdzOQDNnz9ff/7zn3XnnXdq9OjRioqKkiRt2rTJcWoMAACgKXP5a/B33nmnzpw5I5vNpg4dOjjGn3jiCbVp08atzQEAAHhCgx6F4efn5xR+JCk8PNwd/QAAAHhcvQLQgAEDVFhYqA4dOqh///7y8fG5ai1PhG9awrO2/GzNd/OSG6ETAACajnoFoBEjRji++ZWSkuLJfgAAADyuXgEoOzu7zn8DAAA0Rw26BkiSSkpKdPjwYUk/fDW+f//+bmsKAADAk1wOQBUVFRo1apR27Nih9u3bS5LOnTunu+66S2vWrNENN9zg7h4BAADcyuX7AE2ePFnnz5/X559/rrNnz+rs2bM6dOiQbDabnn32WU/0CAAA4FYurwBt3bpV27ZtU58+fRxjffv2VV5enu655x63NgcAAOAJLq8A2e12tWrVqtZ4q1atZLfb3dIUAACAJ7kcgO6++26lp6fr5MmTjrETJ05oypQpGjZsmFubAwAA8ASXA9CSJUtks9kUHh6um2++WTfffLN69uwpm82mxYsXe6JHAAAAt3L5GqCwsDCVlpZq27Zt+vLLLyVJffr0UWJiotubAwAA8IQG3QfIx8dHw4cP1/Dhw93dDwAAgMc1KAAVFhaqsLBQFRUVtS58Xr58uVsaQ+Opz/PCJJ4ZBgBoOVwOQC+88IJmz56tmJgYdenS5ZoPRgUAAGiKXA5A+fn5WrlypcaOHeuJfgAAADzO5W+BXb58WYMHD/ZELwAAAI3C5QD02GOPafXq1Z7oBQAAoFG4HIAuXbqkhQsXaujQoZo8ebIyMzOdXg2Rl5en8PBwBQUFKT4+XsXFxdesX7dunSIiIhQUFKTIyEi9++67V6196qmn5OPjo9zc3Ab1BgAAWh6XrwE6ePCgoqOjJUmHDh1y+qwhF0SvXbtWmZmZys/PV3x8vHJzc5WUlKQjR46oc+fOtep3796t0aNHKycnR//yL/+i1atXKyUlRaWlperXr59T7YYNG/TJJ5+oa9euLvcFAABaLh/DMAxvNhAfH6/Y2FgtWbJE0g/PGgsLC9PkyZOVlZVVqz41NVVVVVXavHmzY2zQoEGKjo5Wfn6+Y+zEiROKj4/X+++/r+TkZGVkZCgjI6NePdlsNoWEhKiyslLBwcG/7ADrUN+vnTc1fA0eANCUufL32+VTYO50+fJllZSUON1F2tfXV4mJiSoqKqpzm6Kiolp3nU5KSnKqt9vtGjt2rKZOnarbbrvtZ/uorq6WzWZzegEAgJbLqwHozJkzqqmpUWhoqNN4aGiorFZrndtYrdafrZ8/f778/f317LPP1quPnJwchYSEOF5hYWEuHgkAAGhOvBqAPKGkpESvvvqqVq5cWe9rkqZNm6bKykrH6/jx4x7uEgAAeJNXA1CnTp3k5+en8vJyp/Hy8nJZLJY6t7FYLNes/+tf/6qKigp1795d/v7+8vf317Fjx/S73/1O4eHhde4zMDBQwcHBTi8AANByeTUABQQEaODAgSosLHSM2e12FRYWKiEhoc5tEhISnOolqaCgwFE/duxYHTx4UAcOHHC8unbtqqlTp+r999/33MEAAIBmo0EPQ3WnzMxMpaWlKSYmRnFxccrNzVVVVZXGjx8vSRo3bpy6deumnJwcSVJ6erqGDh2qBQsWKDk5WWvWrNG+ffu0bNkySVLHjh3VsWNHp5/RqlUrWSwW9e7du3EPDgAANEleD0Cpqak6ffq0Zs2aJavVqujoaG3dutVxoXNZWZl8fX9aqBo8eLBWr16tGTNmaPr06erVq5c2btxY6x5AAAAAV+P1+wA1RdwHqG7cBwgA0JQ1m/sAAQAAeAMBCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmI6/txtA8xGeteVna76bl9wInQAA8MuwAgQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEyHAAQAAEzH39sNoGUJz9ryszXfzUtuhE4AALg6VoAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpNIkAlJeXp/DwcAUFBSk+Pl7FxcXXrF+3bp0iIiIUFBSkyMhIvfvuu47Prly5oueff16RkZFq27atunbtqnHjxunkyZOePgwAANBMeD0ArV27VpmZmcrOzlZpaamioqKUlJSkioqKOut3796t0aNHa8KECdq/f79SUlKUkpKiQ4cOSZIuXryo0tJSzZw5U6WlpVq/fr2OHDmiBx54oDEPCwAANGE+hmEY3mwgPj5esbGxWrJkiSTJbrcrLCxMkydPVlZWVq361NRUVVVVafPmzY6xQYMGKTo6Wvn5+XX+jL179youLk7Hjh1T9+7df7Ynm82mkJAQVVZWKjg4uIFHdnXhWVvcvs/m5Lt5yd5uAQDQArny99urK0CXL19WSUmJEhMTHWO+vr5KTExUUVFRndsUFRU51UtSUlLSVeslqbKyUj4+Pmrfvn2dn1dXV8tmszm9AABAy+XVAHTmzBnV1NQoNDTUaTw0NFRWq7XObaxWq0v1ly5d0vPPP6/Ro0dfNQ3m5OQoJCTE8QoLC2vA0QAAgObC39sNeNKVK1f0r//6rzIMQ0uXLr1q3bRp05SZmel4b7PZCEEeVJ9TgJwmAwB4klcDUKdOneTn56fy8nKn8fLyclksljq3sVgs9ar/MfwcO3ZMH3744TXPBQYGBiowMLCBRwEAAJobr54CCwgI0MCBA1VYWOgYs9vtKiwsVEJCQp3bJCQkONVLUkFBgVP9j+Hnq6++0rZt29SxY0fPHAAAAGiWvH4KLDMzU2lpaYqJiVFcXJxyc3NVVVWl8ePHS5LGjRunbt26KScnR5KUnp6uoUOHasGCBUpOTtaaNWu0b98+LVu2TNIP4ec3v/mNSktLtXnzZtXU1DiuD7r++usVEBDgnQOFSzhNBgDwJK8HoNTUVJ0+fVqzZs2S1WpVdHS0tm7d6rjQuaysTL6+Py1UDR48WKtXr9aMGTM0ffp09erVSxs3blS/fv0kSSdOnNCmTZskSdHR0U4/a/v27brzzjsb5bgAAEDT5fX7ADVF3AeoeWAFCADwz5rNfYAAAAC8gQAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMx9/bDQANFZ61pV51381L9nAnAIDmhhUgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOgQgAABgOjwKAy1efR6ZweMyAMBcWAECAACmwwoQIFaJAMBsWAECAACmQwACAACmQwACAACmQwACAACmQwACAACmQwACAACmQwACAACmw32AADfifkIA0DywAgQAAEyHFSCgnuqzugMAaB5YAQIAAKbDChDQRHE9EQB4DitAAADAdAhAAADAdDgFBjRjnCYDgIZhBQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOX4MHWji+Kg8AtRGAABCSAJgOp8AAAIDpsAIENLL6rLaAVSkAnsUKEAAAMB1WgADUS31Xrpraqkxz7RuAZzWJAJSXl6eXXnpJVqtVUVFRWrx4seLi4q5av27dOs2cOVPfffedevXqpfnz5+v+++93fG4YhrKzs/Xaa6/p3LlzGjJkiJYuXapevXo1xuEAptaYp644nQigobx+Cmzt2rXKzMxUdna2SktLFRUVpaSkJFVUVNRZv3v3bo0ePVoTJkzQ/v37lZKSopSUFB06dMhR8+KLL2rRokXKz8/Xnj171LZtWyUlJenSpUuNdVgAAKAJ83oAWrhwoR5//HGNHz9effv2VX5+vtq0aaPly5fXWf/qq6/q3nvv1dSpU9WnTx/NmTNHAwYM0JIlSyT9sPqTm5urGTNmaMSIEbr99tu1atUqnTx5Uhs3bmzEIwMAAE2VV0+BXb58WSUlJZo2bZpjzNfXV4mJiSoqKqpzm6KiImVmZjqNJSUlOcLNt99+K6vVqsTERMfnISEhio+PV1FRkUaNGlVrn9XV1aqurna8r6yslCTZbLYGH9u12KsvemS/QHNRn9+txv498dTvO4DG8+PvsWEYP1vr1QB05swZ1dTUKDQ01Gk8NDRUX375ZZ3bWK3WOuutVqvj8x/Hrlbzf+Xk5OiFF16oNR4WFla/AwHgkpBcb3dQW1PsCUDDnD9/XiEhIdesaRIXQXvbtGnTnFaV7Ha7zp49q44dO8rHx8etP8tmsyksLEzHjx9XcHCwW/eNnzDPjYN5bhzMc+NgnhuHJ+fZMAydP39eXbt2/dlarwagTp06yc/PT+Xl5U7j5eXlslgsdW5jsViuWf/jf8vLy9WlSxenmujo6Dr3GRgYqMDAQKex9u3bu3IoLgsODuYXrBEwz42DeW4czHPjYJ4bh6fm+edWfn7k1YugAwICNHDgQBUWFjrG7Ha7CgsLlZCQUOc2CQkJTvWSVFBQ4Kjv2bOnLBaLU43NZtOePXuuuk8AAGAuXj8FlpmZqbS0NMXExCguLk65ubmqqqrS+PHjJUnjxo1Tt27dlJOTI0lKT0/X0KFDtWDBAiUnJ2vNmjXat2+fli1bJkny8fFRRkaG/vjHP6pXr17q2bOnZs6cqa5duyolJcVbhwkAAJoQrweg1NRUnT59WrNmzZLValV0dLS2bt3quIi5rKxMvr4/LVQNHjxYq1ev1owZMzR9+nT16tVLGzduVL9+/Rw1v//971VVVaUnnnhC586d0x133KGtW7cqKCio0Y/v/woMDFR2dnatU25wL+a5cTDPjYN5bhzMc+NoKvPsY9Tnu2IAAAAtiNdvhAgAANDYCEAAAMB0CEAAAMB0CEAAAMB0CECNKC8vT+Hh4QoKClJ8fLyKi4u93VKzlpOTo9jYWLVr106dO3dWSkqKjhw54lRz6dIlTZw4UR07dtR1112nhx56qNaNNOGaefPmOW438SPm2T1OnDihf/u3f1PHjh3VunVrRUZGat++fY7PDcPQrFmz1KVLF7Vu3VqJiYn66quvvNhx81NTU6OZM2eqZ8+eat26tW6++WbNmTPH6dlRzHPD7Nq1S7/+9a/VtWtX+fj41HoAeX3m9ezZsxozZoyCg4PVvn17TZgwQRcuXPBIvwSgRrJ27VplZmYqOztbpaWlioqKUlJSkioqKrzdWrO1c+dOTZw4UZ988okKCgp05coV3XPPPaqqqnLUTJkyRe+8847WrVunnTt36uTJkxo5cqQXu27e9u7dqz//+c+6/fbbncaZ51/u73//u4YMGaJWrVrpvffe0xdffKEFCxaoQ4cOjpoXX3xRixYtUn5+vvbs2aO2bdsqKSlJly5d8mLnzcv8+fO1dOlSLVmyRIcPH9b8+fP14osvavHixY4a5rlhqqqqFBUVpby8vDo/r8+8jhkzRp9//rkKCgq0efNm7dq1S0888YRnGjbQKOLi4oyJEyc63tfU1Bhdu3Y1cnJyvNhVy1JRUWFIMnbu3GkYhmGcO3fOaNWqlbFu3TpHzeHDhw1JRlFRkbfabLbOnz9v9OrVyygoKDCGDh1qpKenG4bBPLvL888/b9xxxx1X/dxutxsWi8V46aWXHGPnzp0zAgMDjbfeeqsxWmwRkpOTjUcffdRpbOTIkcaYMWMMw2Ce3UWSsWHDBsf7+szrF198YUgy9u7d66h57733DB8fH+PEiRNu75EVoEZw+fJllZSUKDEx0THm6+urxMREFRUVebGzlqWyslKSdP3110uSSkpKdOXKFad5j4iIUPfu3Zn3Bpg4caKSk5Od5lNint1l06ZNiomJ0cMPP6zOnTurf//+eu211xyff/vtt7JarU7zHBISovj4eObZBYMHD1ZhYaGOHj0qSfr000/10Ucf6b777pPEPHtKfea1qKhI7du3V0xMjKMmMTFRvr6+2rNnj9t78vqdoM3gzJkzqqmpcdzd+kehoaH68ssvvdRVy2K325WRkaEhQ4Y47gputVoVEBBQ68G2oaGhslqtXuiy+VqzZo1KS0u1d+/eWp8xz+7xzTffaOnSpcrMzNT06dO1d+9ePfvsswoICFBaWppjLuv6/wjzXH9ZWVmy2WyKiIiQn5+fampqNHfuXI0ZM0aSmGcPqc+8Wq1Wde7c2elzf39/XX/99R6ZewIQWoSJEyfq0KFD+uijj7zdSotz/Phxpaenq6CgoEk8TqalstvtiomJ0Z/+9CdJUv/+/XXo0CHl5+crLS3Ny921HH/5y1/05ptvavXq1brtttt04MABZWRkqGvXrsyzyXAKrBF06tRJfn5+tb4VU15eLovF4qWuWo5JkyZp8+bN2r59u2688UbHuMVi0eXLl3Xu3DmneubdNSUlJaqoqNCAAQPk7+8vf39/7dy5U4sWLZK/v79CQ0OZZzfo0qWL+vbt6zTWp08flZWVSZJjLvn/yC8zdepUZWVladSoUYqMjNTYsWM1ZcoUxwO3mWfPqM+8WiyWWl8M+sc//qGzZ896ZO4JQI0gICBAAwcOVGFhoWPMbrersLBQCQkJXuyseTMMQ5MmTdKGDRv04YcfqmfPnk6fDxw4UK1atXKa9yNHjqisrIx5d8GwYcP02Wef6cCBA45XTEyMxowZ4/g38/zLDRkypNZtHI4ePaoePXpIknr27CmLxeI0zzabTXv27GGeXXDx4kWnB2xLkp+fn+x2uyTm2VPqM68JCQk6d+6cSkpKHDUffvih7Ha74uPj3d+U2y+rRp3WrFljBAYGGitXrjS++OIL44knnjDat29vWK1Wb7fWbD399NNGSEiIsWPHDuPUqVOO18WLFx01Tz31lNG9e3fjww8/NPbt22ckJCQYCQkJXuy6Zfjnb4EZBvPsDsXFxYa/v78xd+5c46uvvjLefPNNo02bNsYbb7zhqJk3b57Rvn1747//+7+NgwcPGiNGjDB69uxpfP/9917svHlJS0szunXrZmzevNn49ttvjfXr1xudOnUyfv/73ztqmOeGOX/+vLF//35j//79hiRj4cKFxv79+41jx44ZhlG/eb333nuN/v37G3v27DE++ugjo1evXsbo0aM90i8BqBEtXrzY6N69uxEQEGDExcUZn3zyibdbatYk1flasWKFo+b77783nnnmGaNDhw5GmzZtjAcffNA4deqU95puIf5vAGKe3eOdd94x+vXrZwQGBhoRERHGsmXLnD632+3GzJkzjdDQUCMwMNAYNmyYceTIES912zzZbDYjPT3d6N69uxEUFGTcdNNNxh/+8AejurraUcM8N8z27dvr/H9yWlqaYRj1m9f//d//NUaPHm1cd911RnBwsDF+/Hjj/PnzHunXxzD+6faXAAAAJsA1QAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQADgJuHh4crNzfV2GwDqgQAEAL/Q5cuXvd0CABcRgAA0KW+//bYiIyPVunVrdezYUYmJiaqqqtKdd96pjIwMp9qUlBQ98sgjjvfh4eGaM2eORo8erbZt26pbt27Ky8tz2sbHx0dLly7Vfffdp9atW+umm27S22+/7VTz2Wef6e6773b08MQTT+jChQuOzx955BGlpKRo7ty56tq1q3r37q0777xTx44d05QpU+Tj4yMfHx+3zw0A9yEAAWgyTp06pdGjR+vRRx/V4cOHtWPHDo0cOVKuPLP5pZdeUlRUlPbv36+srCylp6eroKDAqWbmzJl66KGH9Omnn2rMmDEaNWqUDh8+LEmqqqpSUlKSOnTooL1792rdunXatm2bJk2a5LSPwsJCHTlyRAUFBdq8ebPWr1+vG2+8UbNnz9apU6d06tSpXz4hADzG39sNAMCPTp06pX/84x8aOXKkevToIUmKjIx0aR9DhgxRVlaWJOnWW2/Vxx9/rFdeeUXDhw931Dz88MN67LHHJElz5sxRQUGBFi9erH//93/X6tWrdenSJa1atUpt27aVJC1ZskS//vWvNX/+fIWGhkqS2rZtq9dff10BAQGO/fr5+aldu3ayWCwNnwQAjYIVIABNRlRUlIYNG6bIyEg9/PDDeu211/T3v//dpX0kJCTUev/j6k59ag4fPqyoqChH+JF+CFV2u11HjhxxjEVGRjqFHwDNCwEIQJPh5+engoICvffee+rbt68WL16s3r1769tvv5Wvr2+tU2FXrlzxUqdyCkgAmh8CEIAmxcfHR0OGDNELL7yg/fv3KyAgQBs2bNANN9zgdF1NTU2NDh06VGv7Tz75pNb7Pn361LumT58++vTTT1VVVeX4/OOPP5avr6969+59zd4DAgJUU1NTvwMF4FUEIABNxp49e/SnP/1J+/btU1lZmdavX6/Tp0+rT58+uvvuu7VlyxZt2bJFX375pZ5++mmdO3eu1j4+/vhjvfjiizp69Kjy8vK0bt06paenO9WsW7dOy5cv19GjR5Wdna3i4mLHRc5jxoxRUFCQ0tLSdOjQIW3fvl2TJ0/W2LFjHdf/XE14eLh27dqlEydO6MyZM26bFwDux0XQAJqM4OBg7dq1S7m5ubLZbOrRo4cWLFig++67T1euXNGnn36qcePGyd/fX1OmTNFdd91Vax+/+93vtG/fPr3wwgsKDg7WwoULlZSU5FTzwgsvaM2aNXrmmWfUpUsXvfXWW+rbt68kqU2bNnr//feVnp6u2NhYtWnTRg899JAWLlz4s/3Pnj1bTz75pG6++WZVV1e79O01AI3Lx+A3FEALER4eroyMjFr3C/pnPj4+2rBhg1JSUhqtLwBND6fAAACA6RCAAACA6XAKDAAAmA4rQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHT+PxDTjGH++9WxAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def make_distribution_counts(edges, scale):\n", " bin_names = [str(i) for i in range(len(edges - 1))]\n", "\n", " return (\n", " input_space >>\n", " dp.t.then_find_bin(edges=edges) >> # bin the data\n", " dp.t.then_index(bin_names, \"0\") >> # can be omitted. Just set TIA=\"usize\", categories=list(range(num_bins)) on next line:\n", " dp.t.then_count_by_categories(categories=bin_names, null_category=False) >>\n", " dp.m.then_base_discrete_laplace(scale)\n", " )\n", "\n", "edges = np.linspace(*bounds, num=50)\n", "counts = make_distribution_counts(edges, scale=1.)(data)\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.hist(range(len(edges)), edges, weights=counts, density=True)\n", "plt.xlabel(\"support\")\n", "plt.ylabel(\"noisy density\");" ] } ], "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.9.13" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "3220da548452ac41acb293d0d6efded0f046fab635503eb911c05f743e930f34" } } }, "nbformat": 4, "nbformat_minor": 2 }