{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Wasserstein Barycenters using Mosek" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wasserstein Distance is a way to measure the distance between two probabilty distributions. It allows to summarize, compare, match and reduce the dimensionality of the emprical probability measures to carry out some machine learning fundamentals." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wasserstein distance of order $p$ between two probabilty measures $\\mu$ and $\\upsilon$ in $P(\\Omega)$ is defined as:\n", "
\n", "
\n", "$$W_p(\\mu,\\upsilon) \\overset{\\underset{\\mathrm{def}}{}}{=} \\bigg( \\underset{\\pi \\in \\Pi{(\\mu, \\upsilon)}}{\\mbox{inf}} \\int_{\\Omega^2} D(X_i,Y_j)^p d\\pi(x,y)\\bigg)^{1/p}\n", "$$\n", "
\n", "where $\\Pi(\\mu, \\upsilon)$ is the set of all probability measures on $\\Omega^2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the distributions are discrete $W_p(\\mu,\\upsilon)$ is equilavent to the objective of the following LP model:\n", "
\n", "
\n", "$$\\mbox{minimize} \\quad \\sum_{i=1}\\sum_{j=1} D(X_i,Y_j)^p\\pi_{ij}$$\n", "
\n", "$$\\mbox{st.} \\quad \\sum_{j=1} \\pi_{ij} = \\mu_i , \\quad i = 1,2,..n$$\n", "
\n", "$$\\quad \\sum_{i=1} \\pi_{ij} = \\upsilon_j, \\quad j = 1,2,..m$$\n", "
\n", "$$\\pi_{ij} \\geq 0, \\quad \\forall_{i,j}$$\n", "
\n", "where $D(X_i,Y_j)$ is the distance function on $\\Omega$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are more efficient ways to approximate this metric but LP approach will be applied in order to compare the performance and modeling structure of Fusion, Pyomo and CVXPY." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wasserstein Barycenter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wasserstein barycenter problem involves all Wasserstein distances from one to many measures. Given measures $\\upsilon_1,\\ldots,\\upsilon_N$ we want to find the measure which minimizes the sum of distances to the given measures, that is solve the problem:\n", "$$\\mbox{minimize}_\\mu\\sum_{i=1} \\lambda_i W_p(\\mu,\\upsilon_i)$$\n", "for some fixed system $\\lambda_i$ of weights of distances to specific distributions that satisfies $$\\sum_{i=1}\\lambda_i = 1.$$\n", "For simplicity uniform weights are used in this problem. Then the barycenters problem becomes:\n", "$$\\mbox{minimize}_\\mu \\frac1N \\sum_{i=1}^{N} W_p(\\mu,\\upsilon_i).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this problem, Wasserstein Barycenter of One's are visualized using images with size $28x28$ using $20$ handwriten '1' digits from MNIST database http://yann.lecun.com/exdb/mnist/. Computations are carried out by Intel(R) Xeon(R) CPU E5-2687W v4 @ 3.00GHz processor. Similar experiments are performed by Cuturi and Doucet in http://proceedings.mlr.press/v32/cuturi14.pdf." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABIEAAAILCAYAAAB/zYxFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xu0nVV5L/5nJoQQrnIRiCGA3MQ7SExBWkul9ofUFmyPVHpq1dITa8VLtT1V2zGk9niK91MvRbHwC7ZWy69oRUqPUg7VU7lIQlNuEUEaJJKCgFyFkGTP3x9ZdgSy5srKus5kfj5jZOy932e/73z2Gvlm7/3kXWumnHMAAAAAsH2bNe0GAAAAABg/QyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0IAdhjk5pXRSRPx5RMyOiL/MOZ/d6/N3THPzTrHLMEvCNuvxeDSeyGvTJNaSTeifbEKdZBPqJJtQp36zmXLOAy2QUpodEd+NiJdHxOqIuDYiTs8531w6Z/e0V/6pdOJA68G27pp8eTyU7x/7N0zZhK0jm1An2YQ6ySbUqd9sDvN0sMURcVvO+fac8xMR8cWIOGWI6wGjIZtQJ9mEOskm1Ek2YQyGGQItiIg7N/l4defYk6SUlqSUlqWUlq2LtUMsB/RJNqFOsgl1kk2ok2zCGAwzBOp2m9Fmzy3LOZ+bc16Uc140J+YOsRzQJ9mEOskm1Ek2oU6yCWMwzBBodUQs3OTjAyLiruHaAUZANqFOsgl1kk2ok2zCGAwzBLo2Ig5PKT0zpbRjRLwmIi4eTVvAEGQT6iSbUCfZhDrJJozBwFvE55zXp5TOjIivxcYt+87POd80ss6Agcgm1Ek2oU6yCXWSTRiPgYdAERE550sj4tIR9QKMiGxCnWQT6iSbUCfZhNEb5ulgAAAAAGwjDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABqww7QboA7zvrFfsXbAzg8Ua7eftEuxtuG++4fqCQAApu22jx5brH331/6iWFv0/jOLtX3/4sqhegIYlDuBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAF2B2vI2le8uFi75LBPF2szMVOsHX/qW4u1vc+7qr/GoHGzn/usYu20i64o1l6/+z3F2o9nnijWXvE75d1Kdrrk28Ua0J8dDjm4WLv4/36pWJudyv83t/jdbyrW9rzA91sY1g77l3fK/fypnyrWyj8l2wEMRmHNO19SrC1/xyeKtRPfXP6+Oe/vB/t597FTF4/8mtMw1BAopbQqIh6OiA0RsT7nvGgUTQHDkU2ok2xCnWQT6iSbMHqjuBPo53LO947gOsBoySbUSTahTrIJdZJNGCGvCQQAAADQgGGHQDkivp5SWp5SWjKKhoCRkE2ok2xCnWQT6iSbMGLDPh3s+JzzXSmlfSPispTSd3LO39z0EzphXRIRsVPsPORyQJ9kE+okm1An2YQ6ySaM2FB3AuWc7+q8vScivhwRm71cds753Jzzopzzojkxd5jlgD7JJtRJNqFOsgl1kk0YvYHvBEop7RIRs3LOD3fe/4WIeN/IOmPkvv+Lo38JqEcXpGJt75GvRj9kc9vz4HP3LNb+625rirV1uXzNOWl2sbb658v/Fhx2SfmaDEc2G7JhQ7G0ZsOPi7X9Zs8r1o59y7Ji7da/3alYm3n88WKNjWSTiIgHlu5SrB3TY65w74bHxtANEbLZkjved1yx9rnX/nmxNtPjmvccXR51HPT3/XS1uZ3vfLRY6/FjeXWGeTrYfhHx5ZTST67zNznn/z2SroBhyCbUSTahTrIJdZJNGIOBh0A559sj4oUj7AUYAdmEOskm1Ek2oU6yCeNhi3gAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0YZncwGvGFhxcUa4d89vZibf04mgGAyuXH1xZrX3v0sGLtN3f/QbH2kflXF2s///NvKtZ2uuTbxRq0ZoeDFhZr5x75+R5n7lisHH/FW4u1w+O6ftqC5u1+zL3F2gvL8evp0AvuKtYG/T01L79pwDPr4k4gAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADbBG/HZr56aO6Hv+Lk5YOdL0/ueqXirUj1iwf6JoAsL1ad+SCYu03d//Hga65ev1jxdqcRwbd7Ba2PzvM379YO/JLPyjWjpgz2D7UB144e6DzoDUPvPa4Yu3vnv+hHmfOHWi99bevGui8FrgTCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAFvEb4e+99vdt6r8uXmP9DjLPBCm5aHTHxr5NWdipljb8QF5h2HN2m23Ym3JZy8a+Xqn/ut/K9b2/+frRr4ebKvWHbxfsfbB/S8t1man8lbvh/7t7xRrh/3D1f01Bo275gPnFGvr8ryBrnnyd04t1mbFnQNdswV+EwAAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAW8QBT9l8OWTHyaz4480SxduBZV458PWjN99/y/GLtl3f554Guedu6tcXaPp/YeaBrQmv+/VXlrMxELtZueuLxYu3Ij64u1tb31xY0YfZhzyzW1uXlxdpMzAy03pw37VisbRjoim3Y4p1AKaXzU0r3pJRu3OTYXimly1JKt3be7jneNoGnkk2ok2xCnWQT6iSbMFn9PB1saUSc9JRj74qIy3POh0fE5Z2PgclaGrIJNVoasgk1WhqyCTVaGrIJE7PFIVDO+ZsRcf9TDp8SERd03r8gIk4dcV/AFsgm1Ek2oU6yCXWSTZisQV8Yer+c85qIiM7bfUufmFJaklJallJati7Kz3UHRkI2oU6yCXWSTaiTbMKYjH13sJzzuTnnRTnnRXNi7riXA/okm1An2YQ6ySbUSTZh6ww6BLo7pTQ/IqLz9p7RtQQMQTahTrIJdZJNqJNswpgMukX8xRHxuog4u/P2KyPriL7ssP9+xdrvvfifRrrWgV+aPdLrMVayCXWSzW3R4vI28Of+9idHvtxrVpxRrO1/eXlrXYYim9ugdb+wqFj71ukf7nHmTsXKr1z09mLt0Duv7qctRks2t0G3nPW0kV/z2Ze/sVh71vdvHvl6Lehni/gvRMRVEfGslNLqlNIZsTGML08p3RoRL+98DEyQbEKdZBPqJJtQJ9mEydrinUA559MLpRNH3AuwFWQT6iSbUCfZhDrJJkzW2F8YGgAAAIDpMwQCAAAAaIAhEAAAAEADDIEAAAAAGjDoFvFMWd5tl2JtydNuG+lau1xZvt6Gka4E27cfvum4rsfftc/He5xlVg/Tcvuv7lqsHTu3fN5Mj2vetm5tsbbPJ3buoytg1S/NLtb2njWvWJuJXKw963+tLtbW99cWNOOB13b/mfa6Ez7a46wdi5UrHit/vz3y7EeKtQ2PP95jPUr8dgEAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAt4gEm5A1nXtr1+CzzeKjSea8+Z+TX/NXP/n6xtvDyK0e+HmzL1rzzJV2Pf/vUD/U4q7xF/BFffVOx9qz/WNFvW9C8x/ZNXY/vnMrbwM+K7udERPzxB36rWNv75qv6b4y++M0DAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADDIEAAAAAGmB3sO3QIDsNzUmzx9AJsKnD5/7HxNY69uJ3lPuIaybWB2wLfvOWO7seP27uhh5nlXc5OeGGVxdrB354ebGWe6wG26s0p7yb0CG//L2ux/ecVd4B7Ftryz8HP/sj9xZrG9Y9UaxBi2bvvVex9l9ff1nX4zMxUzznwZlyxubdVz6P0XMnEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAbaI3w712pqvZJ19aWEk7v+t44q1F839VqGy00BrLV9brh358fuKtV6bXsP2avW7X1KsvWbXTxYq5W3g12z4cbG260m3F2u+3cKTPf7zLyzWvnrYp7se7/WT7nu/d0qxtuOt5WwCT7byQ4cUa1/e62tbfb0X/8PvFWtHfOmarb4eg9vinUAppfNTSveklG7c5NhZKaUfpJRWdP6cPN42gaeSTaiTbEKdZBPqJJswWf08HWxpRJzU5fjHcs5Hdf5cOtq2gD4sDdmEGi0N2YQaLQ3ZhBotDdmEidniECjn/M2IuH8CvQBbQTahTrIJdZJNqJNswmQN88LQZ6aUru/cvrdn6ZNSSktSSstSSsvWRY8XsABGRTahTrIJdZJNqJNswhgMOgQ6JyIOjYijImJNRHyk9Ik553NzzotyzovmxNwBlwP6JJtQJ9mEOskm1Ek2YUwGGgLlnO/OOW/IOc9ExGcjYvFo2wIGIZtQJ9mEOskm1Ek2YXwG2iI+pTQ/57ym8+GrIuLGXp/P6N133H4jvd4RX19Srv1oxUjXYnxkc/pe+ubyFpd7zhpsK/iS07/2pmLtiFu+PdK1GI5sTsaGE15UrF34xuJ/IsdM7LjVa5207I3F2oK4aauvx3TI5vTt9Ad3bfU59254rFhLH9inx5l3bPVaTIdsTt/T93twpNd79h/dVqxtGOlKbMkWh0AppS9ExAkRsU9KaXVEvDciTkgpHRUROSJWRUT5JyFgLGQT6iSbUCfZhDrJJkzWFodAOefTuxw+bwy9AFtBNqFOsgl1kk2ok2zCZA2zOxgAAAAA2whDIAAAAIAGGAIBAAAANMAQCAAAAKABA20Rz/Q9cspDI71eerTHX4UZm/bBph78jWOLtT/d7+M9zpy91WvNxEyxtteKrb8ebM/WvGVtsXbEnK3fBv7k75xarC34FdvAQ792OGhhsfbBQy7scWb33B5/xVuLZxz+T8v7bQuaN/tZhxVr3zrqiz3O7H4vyRvuOLF4xob77u+3LcbMnUAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAbYIp6IiNhlla2moV/3vLhcm5NGm6X/ee9Rxdo+n7lqpGvBtmDDCS8q1i485lM9ztz6LeJXf6O8rfWBcedWXw9a9eBn5hRrR8zZ+mweeKGfW2EUbvmdfYq1mZgp1v51bfd7SVa///DiOXPj2v4bY6zcCQQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAt4is2a7fdirUDnvbgSNd6xoevHOn1gNH42od+pljbI66eYCdQhzVvWVusDbLVdETE5Y/t3PX4wR+7oXhOeeNcaFM+7oXF2jee///2ODMVKz97w3/penyXf7DVNIzCd077VLHW6/vcW24+vevxvWRzm+BOIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAA2wRX7GZZx9crF185PkjXWuHgxYWa+vvuHOkawFP9jt3/myxtudXbirWbFFNi15ywKqRX/NDb/yNrsd3eHj5yNeC7dUDz9q5WJuJPNA1H79ov67Hd4nbB7oeMBoPXr931+N7TbgPBrPFO4FSSgtTSleklFamlG5KKb2tc3yvlNJlKaVbO2/3HH+7wE/IJtRJNqFOsgl1kk2YrH6eDrY+It6Zc352RBwbEW9OKT0nIt4VEZfnnA+PiMs7HwOTI5tQJ9mEOskm1Ek2YYK2OATKOa/JOV/Xef/hiFgZEQsi4pSIuKDzaRdExKnjahLYnGxCnWQT6iSbUCfZhMnaqheGTikdHBFHR8Q1EbFfznlNxMbgRsS+hXOWpJSWpZSWrYu1w3ULdCWbUCfZhDrJJtRJNmH8+h4CpZR2jYiLIuLtOeeH+j0v53xuznlRznnRnJg7SI9AD7IJdZJNqJNsQp1kEyajryFQSmlObAzk53POX+ocvjulNL9Tnx8R94ynRaBENqFOsgl1kk2ok2zC5Gxxi/iUUoqI8yJiZc75o5uULo6I10XE2Z23XxlLhw2b9cSGYu2u9eVbHZ+xw9ZPwF96yXeKtc99d3GxNuuaPcp9fPDKre6D/snm9uO+tbsUazMP3z3BThgF2Rzehp97UbH2nv0/3uPMeQOtt9Pt93Y9vn6gq1Er2RzerN12K9Yu+dMP9zhzp2JlzYbHirW9b3i0n7bYxsnmeN3xvuOKtVlxXbH25Uf3KdYO/8s1XY/7vrlt2OIQKCKOj4jXRsQNKaUVnWPviY1hvDCldEZEfD8iXj2eFoEC2YQ6ySbUSTahTrIJE7TFIVDO+V8iIhXKJ462HaBfsgl1kk2ok2xCnWQTJmurdgcDAAAAYNtkCAQAAADQAEMgAAAAgAYYAgEAAAA0oJ/dwZiSmRU3F2uvft8fFGuXntV9i849Zu1YPOcde5W3iH/D4hXF2nGPnVmsAf1518JLi7V3vObNxdpuX7x6HO3A1N11fHk76QN2GGwb+J/5t18r1p62+vsDXRNac8dbn1+s7T3rG8XaTORi7c71O5cXvPr6vvoCesil19zunc0/+tKvF2uH3H7VUC0xXe4EAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBbx26i9zytvy3fsz7y16/GVL//0QGv91D91v15ExBFvWD7QNWFbdsjfP1GsffeUcu2IOTt2PX7M3PJaL3/3/y3Wrv7inPKJsA075Of/feTX/NHD5W2o91hXzi0wXn/8piXF2o6xbIKdwPZpj9vK28D3PO97I26EargTCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgd7Dt0OGv775j1y/Hiwe63hFhBzDY1OwrrivW3vpbZxZrnzz/E12PHzanvD3YX/3zzxRrh8fVxRpsyx74xIHl4sfLpU8/cEixdti7HizW1vfTFBAL339lsXby+1800DXtAAbj9bS/Ku8qHWdPrg/q4U4gAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADbBEPMEI7/J/lxdrbD37JVl/PNvC0aJeLrinWXnnRMQNe9Y4BzwOA7dMrF5S/p+4dPbaWZ5u2xTuBUkoLU0pXpJRWppRuSim9rXP8rJTSD1JKKzp/Th5/u8BPyCbUSTahTrIJdZJNmKx+7gRaHxHvzDlfl1LaLSKWp5Qu69Q+lnP+8PjaA3qQTaiTbEKdZBPqJJswQVscAuWc10TEms77D6eUVkbEgnE3BvQmm1An2YQ6ySbUSTZhsrbqhaFTSgdHxNER8ZMn65+ZUro+pXR+SmnPwjlLUkrLUkrL1sXaoZoFupNNqJNsQp1kE+okmzB+fQ+BUkq7RsRFEfH2nPNDEXFORBwaEUfFxsntR7qdl3M+N+e8KOe8aE7MHUHLwKZkE+okm1An2YQ6ySZMRl9DoJTSnNgYyM/nnL8UEZFzvjvnvCHnPBMRn42IxeNrE+hGNqFOsgl1kk2ok2zC5PSzO1iKiPMiYmXO+aObHJ+/yae9KiJuHH17QIlsQp1kE+okm1An2YTJ6md3sOMj4rURcUNKaUXn2Hsi4vSU0lERkSNiVUS8cSwdAiWyCXWSTaiTbEKdZBMmqJ/dwf4lIlKX0qWjbwfol2xCnWQT6iSbUCfZhMnaqt3BAAAAANg2GQIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAEp5zy5xVL6YUTc0flwn4i4d2KL91ZLL/rYXC29jKKPg3LOTx9FM6Mmm1ukj83V0otsTkctvehjc7X0IpuTV0sfEfX0UksfEfX0IpuTV0sfEfX0oo/NTSybEx0CPWnhlJblnBdNZfGnqKUXfWyull5q6WMSavpaa+lFH5urpZda+piEmr7WWnrRx+Zq6aWWPiahlq+1lj4i6umllj4i6umllj4moZavtZY+IurpRR+bm2Qvng4GAAAA0ABDIAAAAIAGTHMIdO4U136qWnrRx+Zq6aWWPiahpq+1ll70sblaeqmlj0mo6WutpRd9bK6WXmrpYxJq+Vpr6SOinl5q6SOinl5q6WMSavlaa+kjop5e9LG5ifUytdcEAgAAAGByPB0MAAAAoAGGQAAAAAANmMoQKKV0UkrplpTSbSmld02jh04fq1JKN6SUVqSUlk147fNTSveklG7c5NheKaXLUkq3dt7uOaU+zkop/aDzuKxIKZ08gT4WppSuSCmtTCndlFJ6W+f4NB6TUi8Tf1wmTTZls0sfVWSz5VxGyGZnbdl8ch+yWQHZlM0ufcjmlNWSy04vsimb/fYxscdk4q8JlFKaHRHfjYiXR8TqiLg2Ik7POd880UY29rIqIhblnO+dwtovjYhHIuJzOefndY59MCLuzzmf3fkHa8+c8x9OoY+zIuKRnPOHx7n2U/qYHxHzc87XpZR2i4jlEXFqRLw+Jv+YlHo5LSb8uEySbP7n2rL55D6qyGaruYyQzU3Wls0n9yGbUyab/7m2bD65D9mcoppy2elnVcimbPbXx8SyOY07gRZHxG0559tzzk9ExBcj4pQp9DFVOedvRsT9Tzl8SkRc0Hn/gtj4l2EafUxcznlNzvm6zvsPR8TKiFgQ03lMSr1s72QzZLNLH1Vks+FcRshmRMhmlz5kc/pkM2SzSx+yOV1y2SGbm/Uhmx3TGAItiIg7N/l4dUzvH6QcEV9PKS1PKS2ZUg+b2i/nvCZi41+OiNh3ir2cmVK6vnP73thvE9xUSungiDg6Iq6JKT8mT+klYoqPywTIZplsRj3ZbCyXEbLZi2yGbE6RbJbJZsjmlNSUywjZ7EU2p5TNaQyBUpdj09qn/vic84si4hUR8ebOrWpEnBMRh0bEURGxJiI+MqmFU0q7RsRFEfH2nPNDk1q3z16m9rhMiGzWr/lsNpjLCNncFsimbP6EbNZFNtvLZk25jJDNEtmcYjanMQRaHRELN/n4gIi4awp9RM75rs7beyLiy7Hx9sFpurvzHMGfPFfwnmk0kXO+O+e8Iec8ExGfjQk9LimlObExCJ/POX+pc3gqj0m3Xqb1uEyQbJbJZgXZbDSXEbLZi2zK5jTJZplsyua0VJPLCNkskc3pZnMaQ6BrI+LwlNIzU0o7RsRrIuLiSTeRUtql80JMkVLaJSJ+ISJu7H3W2F0cEa/rvP+6iPjKNJr4SQg6XhUTeFxSSikizouIlTnnj25SmvhjUuplGo/LhMlmmWxOOZsN5zJCNnuRTdmcJtksk03ZnJYqchkhm73I5pSzmXOe+J+IODk2vmr79yLij6bUwyER8W+dPzdNuo+I+EJsvM1rXWycWJ8REXtHxOURcWvn7V5T6uOvIuKGiLg+NoZi/gT6+OnYeKvm9RGxovPn5Ck9JqVeJv64TPqPbMpmlz6qyGbLuex8/bIpm0/tQzYr+CObstmlD9mc8p8actnpQzbLfcjmFLM58S3iAQAAAJi8aTwdDAAAAIAJMwQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABhgCAQAAADTAEAgAAACgAYZAAAAAAA0wBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANGCHSS62Y5qbd4pdJrkkVOPxeDSeyGvTtPvoRjZpmWxCnWQT6iSbUKd+sznUECildFJE/HlEzI6Iv8w5n93r83eKXeKn0onDLAnbrGvy5RNbSzahf7IJdZJNqJNsQp36zebATwdLKc2OiE9FxCsi4jkRcXpK6TmDXg8YDdmEOskm1Ek2oU6yCeMxzGsCLY6I23LOt+ecn4iIL0bEKaNpCxiCbEKdZBPqJJtQJ9mEMRhmCLQgIu7c5OPVnWNPklJaklJallJati7WDrEc0CfZhDrJJtRJNqFOsgljMMwQqNsLDuXNDuR8bs55Uc550ZyYO8RyQJ9kE+okm1An2YQ6ySaMwTBDoNURsXCTjw+IiLuGawcYAdmEOskm1Ek2oU6yCWMwzBDo2og4PKX0zJTSjhHxmoi4eDRtAUOQTaiTbEKdZBPqJJswBgNvEZ9zXp9SOjMivhYbt+w7P+d808g6AwYim1An2YQ6ySbUSTZhPAYeAkVE5JwvjYhLR9QLMCKyCXWSTaiTbEKdZBNGb5ingwEAAACwjTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCASQi4KXAAAZnUlEQVQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADdph2AwAAwPZrhwXPKNbWHbxvsfbIATtt9Vrzfrhuq8+JiNjh/ywv1mYf9sxi7dJvfrlY+8WX/HKxtn7V9/trDGDE3AkEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiALeIBpuzpVz6tWPvrg/+5WFv87jcVa3tecNUwLQHAyNz6loOKtW/8+oeKtX1mz9vqtR6ZWVusbYhcrN34xG7F2n6zryxfM5e3sb//Jc8o1na3RTwwJe4EAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAbYHYwteuyUxcXar/3ZPxZrS/ZYVayd+LvlXY3mfeXbffUF24uZnIq1dXlDsfaKd3yzWLv6gjlD9QQtmbVT9919vnv2UcVzTjj2xmLte3/y7GJt7j9e239jsJ3Y9Y5y7SWXvqNYSzt1/x540nNvKp5z6+8dWaz90QWfK9aOnvtosbZzmlus9bLHdx8p1sr7lEF7Hjr92GLtmx/+VLH2ygXHjKOd7d5QQ6CU0qqIeDgiNkTE+pzzolE0BQxHNqFOsgl1kk2ok2zC6I3iTqCfyznfO4LrAKMlm1An2YQ6ySbUSTZhhLwmEAAAAEADhh0C5Yj4ekppeUppSbdPSCktSSktSyktWxdrh1wO6JNsQp1kE+okm1An2YQRG/bpYMfnnO9KKe0bEZellL6Tc37SK5XmnM+NiHMjInZPe3kNNJgM2YQ6ySbUSTahTrIJIzbUnUA557s6b++JiC9HRHkbKWBiZBPqJJtQJ9mEOskmjN7AdwKllHaJiFk554c77/9CRLxvZJ1RjfuPLP816bUN/EzMDHTNBV/pqy0KZLNes591WNfjJ+71LwNdb04qbx8fYYv42shmvVb99xd1Pf7dV3+yeM7sVP5/tOcsekGxtvAf+++LyZDN8Xv6OVeVawNc73s9arNiRbH2P3779cXaXS/ZqVj7t9/9RB9dbW72vQ8Va+sHumJbZLMdz3nbjcVar98p7zvjuGJt7/PK/+60bping+0XEV9OKf3kOn+Tc/7fI+kKGIZsQp1kE+okm1An2YQxGHgIlHO+PSJeOMJegBGQTaiTbEKdZBPqJJswHraIBwAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQgGF2B6MRs477UbkWqdeZxcqCKx4eoiPYNj30vL27Hv/N3X8w4U4AoD2zr7iuWHvs119crPX6efeFHz+zWFuw6sr+GoMG/OAPX1KsXbrwk8XaTI/fKW0DPxh3AgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADDIEAAAAAGmCLeCIi4r4zjivWLnnRh4q1mZhXrH3qgUPLC377hr76Asr++jvl7WwPChmDcbp3w6PF2kFffaBYmxlHM0Bfdth/v2LtYyd8oVibiVysLfiAbeChH48e/kSx1itjM75zjpw7gQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADbBFPBER8cCJjxVr82eXt4GfFalYu+Cck4u1fcN2mrTn3hfMHun1dv+HXUd6PWjV4wvWbfU5P849trNdcfMw7QBjcs/JhxRrPzX3omLt7g3j6Aa2P+tfdkyx9oWXfaZY6/U75ZI7X9ZjxYf7aYuncCcQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABW9wiPqV0fkS8MiLuyTk/r3Nsr4j424g4OCJWRcRpOecfja9NRmLx84ulvz72vGJtJmaKteVry3PE+Zf/sFiz0+bwZHPb87KTr5t2C0yAbG57Lnj5Z6fdAhMgm9y7uPwT6D6z5xVrR1745mLtsLh6qJ6Qze3JHb9dztjRc8u/U17b43fKO999eLE2O/xsPYh+7gRaGhEnPeXYuyLi8pzz4RFxeedjYLKWhmxCjZaGbEKNloZsQo2WhmzCxGxxCJRz/mZE3P+Uw6dExAWd9y+IiFNH3BewBbIJdZJNqJNsQp1kEyZr0NcE2i/nvCYiovN239G1BAxBNqFOsgl1kk2ok2zCmGzxNYGGlVJaEhFLIiJ2ip3HvRzQJ9mEOskm1Ek2oU6yCVtn0DuB7k4pzY+I6Ly9p/SJOedzc86Lcs6L5sTcAZcD+iSbUCfZhDrJJtRJNmFMBh0CXRwRr+u8/7qI+Mpo2gGGJJtQJ9mEOskm1Ek2YUz62SL+CxFxQkTsk1JaHRHvjYizI+LClNIZEfH9iHj1OJtkNI75zL8Vay+em4q1mR6zwtf+3ZnF2iErr+qvMQYim21YuW5dsbbrmnKN6ZFNqJNstiEff1SxduXJHy3WvvX47sXaYb9nG/hxks3tx+8ddXmxNqvH75T/885fLNZmX2Eb+FHb4hAo53x6oXTiiHsBtoJsQp1kE+okm1An2YTJGvTpYAAAAABsQwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0YIu7g7Htue+M47oe/529P1Q8ZybmFWsvvf60Yu2Q/24beNjUDgsPKNaO2PnGrb7eX93fPc8REXO+vmyrrwcA27NHDtipWNtndvnn3VvWjaMb2E4tfn7Xw0v2WFo8ZSZmirXbLz2kWFsQ/9F3W/THnUAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAbYIn479IZ3XtL1+Pwe22LeveGxYm39/7dvj9W+129b0IT7f7q8Rfybn/aVCXYCAACjt+Y967senxWpeM7yteX7TxZ84Mqhe6J/7gQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANAAQyAAAACABtgdbBt13xnHFWtL9vhk1+MzMVM85xXLlxRrzzj/qv4bAwCAKZl53b0DnXfh/Yt7VB8frBnYTr3yoJu6Hp+JXDznN6767WLt0PjXoXuif+4EAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBbx26gnfumBYm1WpGKl5Bl/NnvIjgBg2zQ7ZrofT+Xvm3PG1QywRbMPe2ax9q0XXtjjzNLPyBH/uPwFxdoR8e1+2oLtykOnH1us/cm+n+p6vPx7aMRuV80buidGY4t3AqWUzk8p3ZNSunGTY2ellH6QUlrR+XPyeNsEnko2oU6yCXWSTaiTbMJk9fN0sKURcVKX4x/LOR/V+XPpaNsC+rA0ZBNqtDRkE2q0NGQTarQ0ZBMmZotDoJzzNyPi/gn0AmwF2YQ6ySbUSTahTrIJkzXMC0OfmVK6vnP73p6lT0opLUkpLUspLVsXa4dYDuiTbEKdZBPqJJtQJ9mEMRh0CHRORBwaEUdFxJqI+EjpE3PO5+acF+WcF82JuQMuB/RJNqFOsgl1kk2ok2zCmAw0BMo5351z3pBznomIz0bE4tG2BQxCNqFOsgl1kk2ok2zC+Ay0RXxKaX7OeU3nw1dFxI29Pp8BLX5+sXTJi84p1mai+/Z7L73+tOI5u3/7hv77olqyCXWSzbptKPyf2Ibcfev4iIhfW/naYm1e/PvQPTEZsrltuuvk+QOdNxO5WNvz32YP2g5jIJvT95y3lR/ymej+/fHTDxxWPGf+58rX29B/W4zAFodAKaUvRMQJEbFPSml1RLw3Ik5IKR0VETkiVkXEG8fYI9CFbEKdZBPqJJtQJ9mEydriECjnfHqXw+eNoRdgK8gm1Ek2oU6yCXWSTZisYXYHAwAAAGAbYQgEAAAA0ABDIAAAAIAGGAIBAAAANGCgLeKZjIcO3aVYmz+7+zbwERGzInU9vssH9xi6J6C3B3/1kZFe76uXHFusHRRXjXQtoH9rVuxfrB1ii3gYq4eeNdiG0uc9eGCxtt9f99oOG7ZPOyw8oFg7d+HFxdpM4V6Sz3zuF4vnLHjoyv4bY6zcCQQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAt4iv2sfd/qlib6bFZ5fK13Wd7O/7Hw8VzBttoE9q0/mXHFGt/dcw5Pc6c3fXoj2YeL56x73LpBIBNveB5qwY67+51exRrMw+Xf06G7dXN792/WJuJ3KPW/XfRvW9eP3RPjJ87gQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANMAQCAAAAKABhkAAAAAADbBFfMVePDcVazM95ncXP3h01+MbVt46dE9AxH3PnVusvWDH7tvA9/JHd/1CsTbvK9/e6usB47f/1d23xwVGaPHzux5+ywF/M9DlXrXHdcXaYbcsKNb+5mU/Vayt/8FdA/UCk7L+ZccUa999xWeKtVlR/l108Z+9revxfb96Zf+NMTXuBAIAAABogCEQAAAAQAMMgQAAAAAaYAgEAAAA0ABDIAAAAIAGGAIBAAAANGCLW8SnlBZGxOciYv+ImImIc3POf55S2isi/jYiDo6IVRFxWs75R+NrtT0zkXvUylvTfnXV87oef0bcPHRP1EM2p+f0JZeN9HrfuPwFxdoz46qRrsX4yWYbdrv5vmJtwwT7oH/bYzbT3LnF2obFzynWnv2xG4u1X9lzebF2z4bdirWVj5W3WP/6miOLtYN2Kz/Uf3rAOV2PH7jDvOI50WNb62fPmVOsLZx9Z7H23refVqwd+ge2iB/W9pjNmtz33PK/E71+p1y+tny/yPzLf9j1uO9/24Z+7gRaHxHvzDk/OyKOjYg3p5SeExHviojLc86HR8TlnY+ByZFNqJNsQp1kE+okmzBBWxwC5ZzX5Jyv67z/cESsjIgFEXFKRFzQ+bQLIuLUcTUJbE42oU6yCXWSTaiTbMJkbdVrAqWUDo6IoyPimojYL+e8JmJjcCNi38I5S1JKy1JKy9bF2uG6BbqSTaiTbEKdZBPqJJswfn0PgVJKu0bERRHx9pzzQ/2el3M+N+e8KOe8aE6Un48IDEY2oU6yCXWSTaiTbMJk9DUESinNiY2B/HzO+Uudw3enlOZ36vMj4p7xtAiUyCbUSTahTrIJdZJNmJwtDoFSSikizouIlTnnj25SujgiXtd5/3UR8ZXRtweUyCbUSTahTrIJdZJNmKwtbhEfEcdHxGsj4oaU0orOsfdExNkRcWFK6YyI+H5EvHo8LbZrTppdrK0r7x4f+ZqnjaEbKiSb24n9lpW352SbJJsN2HDLbdNuga233WXz8Ze9oFi77C8/XazN6rGN+kz0+CEz7i9WXrVLufbH+5S3pO+9Xq+t4LfeTU+sL9Ze9/HfL9YO/diVI+2DzWx32Zy0HRYeUKytePdfFGvrcvmekDPff2axtvfKq/prjCptcQiUc/6XiOJ3ihNH2w7QL9mEOskm1Ek2oU6yCZO1VbuDAQAAALBtMgQCAAAAaIAhEAAAAEADDIEAAAAAGmAIBAAAANCAfraIZ4zuO+O4Ym1dXl6szUSPLaV77bQJDO3HMztOuwVghP7X6pd3PX78oV8rnvO9vzmqWDv011cUazBKOy9bVawddfVvFmvp2j2Ktcdf8ONiba89Hi3W7n9g12Jtp3lPFGu5x8+trzzkpq7HnztvdfGcDy49rVjb4/byz88LLv9OsbahWIE63PHrBxZr63L5b3Cv3yn3Ps828NsrdwIBAAAANMAQCAAAAKABhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgi/gpe/zpqVibk2YXa1c/Xr7mgg9cOUxLwBZcc8bRxdrlX7ytWLtn/W5dj+92ywPFc2xLC+N370ef2fX4Y58sb2t9y8+eX6wd9ftnFmvP+LDv0YzOhh/+sFg74FfLtXHYcwzXXFE8fkDxnAUxWMZ8v2WbVv6VsufvlIf83e8Wa4fHNcN0RMXcCQQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANsDvYlO190/pibV0u71PwhgveUqwdOOCuCEB/8rIbi7WPHfbsAa54y+DNAEOb9/ff7nr8tH/+f4rnrPzAEeXrzRu6JQDoXy6XPvGjg4q1I/94ZbFmx7ztlzuBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAAN2OIW8SmlhRHxuYjYPyJmIuLcnPOfp5TOioj/FhE/7Hzqe3LOl46r0e3VTl/tvi1tRMQrv3pMsWYbeGQT6iSb248NDzxYrB3xxmsn2AmjIJtQJ9kc3oIPlH83vOQDe/Y486HRN0P1tjgEioj1EfHOnPN1KaXdImJ5SumyTu1jOecPj689oAfZhDrJJtRJNqFOsgkTtMUhUM55TUSs6bz/cEppZUQsGHdjQG+yCXWSTaiTbEKdZBMma6teEyildHBEHB0R13QOnZlSuj6ldH5Kqdd9ZsAYySbUSTahTrIJdZJNGL++h0AppV0j4qKIeHvO+aGIOCciDo2Io2Lj5PYjhfOWpJSWpZSWrYu1I2gZ2JRsQp1kE+okm1An2YTJ6GsIlFKaExsD+fmc85ciInLOd+ecN+ScZyLisxGxuNu5Oedzc86Lcs6L5sTcUfUNhGxCrWQT6iSbUCfZhMnZ4hAopZQi4ryIWJlz/ugmx+dv8mmviogbR98eUCKbUCfZhDrJJtRJNmGy+tkd7PiIeG1E3JBSWtE59p6IOD2ldFRE5IhYFRFvHEuHQIlsQp1kE+okm1An2YQJ6md3sH+JiNSldOno2wH6JZtQJ9mEOskm1Ek2YbK2ancwAAAAALZNhkAAAAAADTAEAgAAAGiAIRAAAABAAwyBAAAAABpgCAQAAADQAEMgAAAAgAYYAgEAAAA0wBAIAAAAoAGGQAAAAAANMAQCAAAAaIAhEAAAAEADUs55coul9MOIuKPz4T4Rce/EFu+tll70sblaehlFHwflnJ8+imZGTTa3SB+bq6UX2ZyOWnrRx+Zq6UU2J6+WPiLq6aWWPiLq6UU2J6+WPiLq6UUfm5tYNic6BHrSwiktyzkvmsriT1FLL/rYXC291NLHJNT0tdbSiz42V0svtfQxCTV9rbX0oo/N1dJLLX1MQi1fay19RNTTSy19RNTTSy19TEItX2stfUTU04s+NjfJXjwdDAAAAKABhkAAAAAADZjmEOjcKa79VLX0oo/N1dJLLX1MQk1fay296GNztfRSSx+TUNPXWksv+thcLb3U0sck1PK11tJHRD291NJHRD291NLHJNTytdbSR0Q9vehjcxPrZWqvCQQAAADA5Hg6GAAAAEADDIEAAAAAGjCVIVBK6aSU0i0ppdtSSu+aRg+dPlallG5IKa1IKS2b8Nrnp5TuSSnduMmxvVJKl6WUbu283XNKfZyVUvpB53FZkVI6eQJ9LEwpXZFSWplSuiml9LbO8Wk8JqVeJv64TJpsymaXPqrIZsu5jJDNztqy+eQ+ZLMCsimbXfqQzSmrJZedXmRTNvvtY2KPycRfEyilNDsivhsRL4+I1RFxbUScnnO+eaKNbOxlVUQsyjnfO4W1XxoRj0TE53LOz+sc+2BE3J9zPrvzD9aeOec/nEIfZ0XEIznnD49z7af0MT8i5uecr0sp7RYRyyPi1Ih4fUz+MSn1clpM+HGZJNn8z7Vl88l9VJHNVnMZIZubrC2bT+5DNqdMNv9zbdl8ch+yOUU15bLTz6qQTdnsr4+JZXMadwItjojbcs6355yfiIgvRsQpU+hjqnLO34yI+59y+JSIuKDz/gWx8S/DNPqYuJzzmpzzdZ33H46IlRGxIKbzmJR62d7JZshmlz6qyGbDuYyQzYiQzS59yOb0yWbIZpc+ZHO65LJDNjfrQzY7pjEEWhARd27y8eqY3j9IOSK+nlJanlJaMqUeNrVfznlNxMa/HBGx7xR7OTOldH3n9r2x3ya4qZTSwRFxdERcE1N+TJ7SS8QUH5cJkM0y2Yx6stlYLiNksxfZDNmcItksk82QzSmpKZcRstmLbE4pm9MYAqUux6a1T/3xOecXRcQrIuLNnVvViDgnIg6NiKMiYk1EfGRSC6eUdo2IiyLi7Tnnhya1bp+9TO1xmRDZrF/z2WwwlxGyuS2QTdn8Cdmsi2y2l82achkhmyWyOcVsTmMItDoiFm7y8QERcdcU+oic812dt/dExJdj4+2D03R35zmCP3mu4D3/f3t3jCJFEIZh+KtAk800MlTwFsYGm5mZbeAx9g5eQIxEzAQ39wQm7qoYyB5gz7BoGXQJItOTWX9DPw8UM8ww8E/BmxQ9PRVD9N5veu8/e++/krzKpH1prd3JEsLb3vv78XLJnhyapWpfJtLmOm1uoM2ddplo8xhtarOSNtdpU5tVNtNlos012qxts+IQ6FOSx621h621u0meJ7mYPURr7WTciCmttZMkT5N8Pf6p/+4iydl4fpbkQ8UQfyIYnmXCvrTWWpLXSb733l/+9db0PVmbpWJfJtPmOm0Wt7njLhNtHqNNbVbS5jptarPKJrpMtHmMNovb7L1PX0lOs9y1/TrJedEMj5JcjvVt9hxJ3mW5zOs2y4n1iyT3k3xM8mM83iua402SL0muskTxYMIcT7JcqnmV5PNYp0V7sjbL9H2ZvbSpzQNzbKLNPXc5vr82tfnvHNrcwNKmNg/Moc3itYUuxxzaXJ9Dm4VtTv+LeAAAAADmq/g5GAAAAACTOQQCAAAA2AGHQAAAAAA74BAIAAAAYAccAgEAAADsgEMgAAAAgB1wCAQAAACwA78BY010KJtAWJAAAAAASUVORK5CYII=\n", "text/plain": [ "