{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "##
NYU CSCI-UA 9473 Introduction to ML
\n", " ###
Unsupervised Learning (Part II)
\n", "\n", "
This week we will cover the EM algorithm, spectral clustering as well as some of the linear latent variable models
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise 1. White matter, Gray matter and cerebrospinal fluid\n", "\n", "For this exercise, you will need to install nibabel. You can do this by typing the line 'pip install nibabel' in your terminal. \n", "\n", "We consider the MR image shown below. We want to segment this image into the white matter, gray matter and cerebrospinal fluid. \n", "\n", " - Start by plotting the historgram of all pixel values (you can do this with 'numpy.histogram'). How does the distribution look like ? \n", "\n", "- Use your histogram to initialize the E-M algorithm (means and variance). The plot the resulting distributions on top ofthe histogram. Does your algorithm correctly capture the three clusters?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACDMElEQVR4nO39a4ykWXrfB/5PRGREZsY173Xr6p6ebrIlDqTZtciFIFlLW2uubNnLlbHSkh8swiJkCRDhL/ogUmvYggUCtFcXLLC72qVgQTJgiSIs0yIEwRItwBdB5IocUzRnRj3dPd3VdcvKe2RcMzIz4t0Pmb8T/ziVNTVTVdlVWfU+QCIz4/K+55z3PLf/czkhyzLllFNOOX0nKrzsAeSUU06vPuWCIqeccnoq5YIip5xyeirlgiKnnHJ6KuWCIqeccnoq5YIip5xyeipdmqAIIfyhEMK3QgifhBB++rLuk1NOOV0+hcvIowghFCV9JOnfkHRf0q9L+vEsy775wm+WU045XTpdlkXxQ5I+ybLs0yzLjiX9gqQfvaR75ZRTTpdMpUu67k1J9+z/+5L+d0/6cAghTw/NKaeXT7tZlq1d9MZlCYpwwWszwiCE8B9I+g8u6f455ZTT906fP+mNyxIU9yW9Zf/fkvTQP5Bl2c9L+nkptyhyyulVp8vCKH5d0vshhC+FEMqSfkzSL1/SvXLKKadLpkuxKLIsOw0h/JSkfySpKOlvZFn2jcu4V0455XT5dCnh0e95ELnrkVNOrwJ9Lcuy33PRG3lmZk455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FS6rOa6OeX0XVOj0VChUNBkMpEkZVmmyWSiLMtUKBSUZZnG47FCOGvuzt8hBI3HY52cnLzM4b8RlAuKnL5wCiGoXC5LkorFotbX11UqlXR6eipJOj09jYJibm5OJycnGo1GKpXOtutoNFKhUFAIQcfHxxoMBppMJlF4nJyc6FVo8fg6US4ocvpCCSHx1ltnpznUajWVSiUdHx9LkgqFM2+4VCqpWCzOCIpKpaIQgnq9niqVik5PT3V6eqpms6mjoyPNzc2pUChoe3tbo9Fo5p5YKzk9G+WCIqcvjJrNphYWFpRlmUajkarVamT2LMtULBYVQlClUlGtVlMIIX5+OBzq5OREpVJJi4uLKpfLOj091cnJiU5OTrSwsKB+v68sy9RoNBRC0NzcnEqlkprNpn7rt37rZU//StMzC4oQwluS/ktJ1yRNJP18lmX/jxDCX5D0JyXtnH/0z2dZ9g+fd6A5XS0ql8taW1vTwsKCTk9PNR6PJSniDePxWFmW6fj4WJPJRIuLi1pYWFAIQaenp9E6kM4sglKppJOTE00mE83Pz0uS5ufno/CZTCY6PDxUu92esSYGg4F2d3e1sbGh+fl5VatVlctlDQYDffTRR1/8wlxReh6L4lTSn82y7H8JIdQlfS2E8Cvn7/3VLMv+0vMPL6dXlVqtlubn59Xr9ZRlmebn5yMDj8fjiDEcHx/r+Ph4Bpjkb9wF6czVmJubm8EpHOAcj8cqFAoqlUqaTCaqVCqan59XpVJRsVjUycmJKpWK6vV6vM5kMomfLxQK8ZrD4VCdTkeLi4sRI7l///6MgMlplp5ZUGRZtilp8/zvbgjhX+rscOKcXmNaXl5WsVjU4uKiSqVSZGbwBEmROYvForIsi9hCqVRSpVLR3Nyc5ufnNZlM4nfK5XLEI8bjcbwekY4syyJYeXR0JEkR0MRCwTIplUo6OjqK0RAE03A4jBbN3NycGo2GTk9PFUJQtVpVpVLRycmJhsPhS1jZV5teCEYRQnhH0v9G0v9P0u+T9FMhhD8u6Td0ZnUcXPCd/JDiK0AwLG5AvV5XqVRSCEFZlsXoBRYFkQlJkWHn5+dVLBa1sLCgarWqWq2mhYUFHR8fR6AxyzKdnJzo6Ogo4hVYKNKZUJibm1OWZer3+xqPx3EcYBenp6cqFotxDFghhUJBlUolCpv5+XmNRiOdnp5qOBxqMploYWFB4/E4vl4oFFQsFjUYDL7A1X516blPCgsh1CT9j5J+Nsuy/yaEsCFpV2enl/9FSdezLPsTT7lGHst6hahUKsXw4sLCgprNZmRKmBFNXywWValUdPv2ba2urmpubm5GsHS73WgtVKtV1et11Wo1zc3NaTweR/xiMBio2+2q2+1KUnQrcD2KxWIUWrgvx8fHUbiMRiN1Oh1J0nA41OHhoYbDoQqFgubn59VoNDQajdRoNDSZTNTr9TQej3V0dKRer6ejo6PoJiGEKpWKtra23qQ8jSeeFPZcFkUIYU7S35P0X2VZ9t9IUpZlW/b+X5f0D57nHjl98XTr1q0INkqKJj9uRL1eV6vVUqvVUrVaVaFQUKvV0vr6uqrVanQniGiEEKJFUKlUosvCdcmFKJfLyrJMp6enqtfrks4EBFZBoVDQwsJCtF4QGqPRSL1eT+12O7oux8fH6na7Ojg4UK/X02g0igIA4bSwsBAtB9yO4+NjjUYj9ft9tdttNRoNdTqdN0lYXEjPE/UIkv4LSf8yy7K/Yq9fP8cvJOmPSPr68w0xp8um5eVlra6uajKZREYDOCRUuby8rKWlJTUaDc3NzalYLEo6EwaDwUCdTke1Wk21Wk3ValXz8/MRv3ChwA8WCkzt/49GI9VqNRWLRZXL5XidEIIajYaKxWJ0V3h9fn5ey8vLKhQKMRIiKQKXDx8+1L1799Rut9Xv93V6eqputxstCFwdXJmDgwMNh0MNh0NVKpUYWen1ei/zUb00eh6L4vdJ+vck/XYI4V+cv/bnJf14COGrOnM97kj6U89xj5xeIMFMaG7wAZi0VCpFM71SqajVakUBQfiyWCzq9PRUx8fHOjo60nA41Hg8VqVSifkNWB9YAP4/P3Nzc5pMJhqNRhHUJNx5fHysWq0Wx42mBxMplUrxtzSNmJBjgWsjSYuLi2o2m1pbW9ONGzd0cHCgTqejTqcTBdzOzk7M7pyfn494Be7PaDTS0dFRfK1QKLxx2MXzRD3+qaRwwVt5zsQrQCEEXb9+PTK/pKjd6/V6TI8OIUTMoVQqqV6vq9FoqFqtanl5Wc1mU4uLi1FrI1xOT081Go00HA4joFmv1yOgWS6XValUouCQNCMoyIvAomBs5XI5uhfSWcSi1+tFLAJhUKlU4vuSokWwuLgY3R/uS7IWQqnRaGh9fV3FYlFHR0d69OiRHj16pH6/HxO4sGbK5bKOjo5UqVR0fHwcw64hBLVaLR0cHLwRQiPPzHxN6Nq1azEcyc/y8rIajYaazaYkxQKqxcXFmIsAcEkyU61WixgE369UKpGJJ5NJdFFKpVLU5lmWqVarRfwADAG3A6wBgUROhfv+CK+FhYUofDyfAqFSq9UiWIr7wXwAU/37fJf08YWFBbVareie4Hrt7e2p3W6r3W4rhBAFAJYS4KkXqGH5vO7CIhcUV5DQ1tIUaHznnXcu1LLr6+taWlpSsVicyYg8OTmZ0eyDwSCGJqVpuJN7wYQAlA4ykmaNL4/GB0sgnRrgECaH4Vy4EcrEIuA7/M2c+TxzApzEBUG4EBHhfeZFiBaLq1arqdVqaW9vTzs7O9re3o7AJhaSjyXLspi0VavVoiXyulIuKK4IYS2EELS+vq4bN25EpoLh8a+dwZeWlmIEgdqIdrutg4ODqB1JVpIUcxi4DsCluxO4H7gDJEC5NYHLgUvh1gWgJQwnKVoFzoyMB2vGE6j4mxBruVxWtVrV4uJiTObCNRoOhxqNRtGVWFhYiFYLWAlWlmM1k8lEDx8+nKlEzbIsumiEdQuFgpaXl7W9vR0xnNeNckFxBSiEoPfff38mYYjMQlKpCT0C8hUKBS0uLqper0dLgwjF4eGh+v2+CoWCqtVqZBzcj+PjYw2HQy0uLkYN7VgGbgMMRaQCQeLvMxZJMxWcMDh+P4VhjN8zMcmhgDE9rImrUqvVoouAhdTtdtXpdNTv9yUpRm2wDhgnc2I+zWZT9Xpd9XpdzWZT+/v76na76vV60ZUpl8tqtVoqlUoaDAY6Pj7W+vq63n//ff3qr/5qrG15XSgXFC+Bfvfv/t0xFRkffG5uLoJp1WpV6+vrWlxcjMg8PjXgIgwMETVwfAAMAu0MwbDz8/NqtVqq1+sKIcTcAbIfSYqCmYfD4UxhFQyGZeCCgNc885LX3Y2QZl0lXA7PvVhcXIyp1YPBIK7ZwsJCdIsAMlnH0WgULSbcLLdaiF4ApJ6cnOj09FSlUkmtVkuVSkVra2taX1/X7u6udnZ2IujZ6XRiSJeUcSIjH330kb7yla/o448/fq1wi1xQvECq1+va2NiIqcvj8XgmjEfCDzkIbFZM8VqtFjME9/f3IyD51ltvxfwED1EeHR2p2+1G05kQJkIEoBFtDDOQN1CtVtVsNtVsNmPC0cLCgnZ3d9Xr9XR8fBzdGSIai4uLyrIsRhZwVdDI0pTxYU6KrXBhfA0kRUHiwoTf7gaQ31GtVtVoNKKb4QVpvV5P/X5fg8FA/X4/Zl22Wq3oeoCBYHmARVCfgvBB0DJ/gGGEhqebY0n1+309evRI6+vr2t/fj9miV51yQfECCBMU0x9wDL+ZIqnRaBRBP48KSIoZiWxcTGAiD+5aSIq5C2h80H4YEWYEVIRZwSQqlYoajUZMqfbcitPT05iYNBgMVCqV1O/31Wq1VCgUIj7gzA6jYw1AjI+IBiZ/auVgkXjVJ4KIz2IRVKvViEfA0DDpYDCIuMTR0VH8fL1ej8KUuo/RaBSzQo+OjqIFUy6XIyaEa0Mux3g8jgVkW1tb8XmxFtVqNVpguGW4PleZckHxDFQsFmPNwGQyiTgBCTpof36ILGAS12o1NRqNWCwlKdYuHB4eajQaxdqIVqsVE4hgNJiKzlCDwSBaDl4TQS6At4cDuHS/3jUyKdUIC6/UXFxcjFgC64CLgQDyMKg0xSLw2R2URbh4yNVrLbBQ+CzWDJYT1lOxWIzrMBqNotshKQqUer0+kzGKq4DlAWMjYD2r1HM/iI5Q3LazsxMxH0lx3EdHRzG3g2pV6liuIuWC4nsgNM38/Lxu3rwZ0e1KpRJNeE9jxoeG+brdrvr9ftxshOfw4ymDTjcU5q/jAYQ5FxcX42ewTogSoB09wgCO4NERBEUabZCkTqcTowrM360GfhNR4Z5YAwgBXB9wj5S8SS5rx9p4KzvW3yMyWGuEQsEkeA6OYdC7AqsD4eLj9uiMC0FJMXN1dXVVq6ur2tra0ubmph49eqR2ux1xDs/DQEgNh8MrGxHJBcV3QWjxd999V8vLyzo+Pla73ZYkra6uamNjQ+vr61HTYAq7xiUMd3JyMpO1iKUA4o+fLCmazEQtpGmRVBpVQHN7rgLM5uE9r7OAiWFCIinkBqyursaoARmdWDhEDqSzFGuyJ8ne9KQqSTP5G2hqHxcWF4lNfN4ZF2bF3UBIkEqOReIp4YRBFxcXoyWBUCG6I02Tqlg7xkEGKngT1snc3JyWlpZiqPrevXv66KOP9ODBg4ijMG9copWVFe3u7l7JiEguKL4L+r2/9/dGX5u+jLVaTUtLS7px44ZWV1djNAIBgJBAOEhnm/H4+HiGsdGKMA1hOWnK1J6yjBAhJIlGh0G4Dr0ewEfIg6jVatF/h2E8GcqjJfPz81paWoruBtgGWAjjxuzHZXAwk4YxRFGIODim4QlXgLSUfINVYIm4e8XaMh+iHJ74BWOXSqU4DkKdFH6trq5Ga4N1SMFW8B6EJMJ2fn5e9Xo9Whi//du/rTt37sSOW9KZAKY2ZmVlRe12OzYTviqUC4rvQKVSSe+9916MLHikoF6vR6wCRoaBvJQaDSop5gGg7TCDMf/ZUDAc5j7uA8lHHiFgY3viFGE/18jgBDDv5uZmHIs0BRNhLsx7rov14FYM1gggrId6PWqAVib/wNvt83kKscipcJDWMRcPyZLz4XgD92b8uGzdbjcCjwijfr8fBQrrgEDgWYBNYGVgSaEMEKKLi4u6fft2BLDv3Lmjvb29uG8QMP1+PwLVWZZpd3f3Jezs751yQXEBra2taWlpaUarLS8va3l5OeYdoKFgHM8y9CQhwmeSZsxwPseGDyFEgeNhVI8SeBYjBP7A5kfbsYH5Plrck5dIIILZvMpTmk0V9/Rpz+TkPe7voUzmmQqV4+PjmSgLzEUZN9oXJk3XgXEiUBByXhhGijbduvv9fgQdcQUkRWHmIDCRJe4pzYZ2UQY+DsZVLpd1+/ZtLS4uanNzU9vb2xEUxgLxnqLLy8va399/0Vv4hVMuKBJqNBpRC6FZrl+/ruvXr2t5eVn1ej1uciwCtCAby4E3SbEaEXOe73hI0LW9hxGlWUbz6zowiVuD0HF3Bg3o2IGkmMUIdsB1XXNjsaQ5Erg8JycnkUGGw2EUOo5FeP8K5sv9eB8G4vsAjAhqnwPCFQAUxkutCcKgZKPu7e3FqJKkGLlIhQSC1NeAdXTLgvsyJwQGOE6r1dLq6qoePHighw8f6uDgYCavhVBwtVp95UOouaAwAtEGeKzVatrY2NDt27e1tramWq0WNQnMTvcjmMtz/WFMfsrlsur1egzHYYIjcKjIdBDUE5q4ppvDaFk0taSZ/AyI63pxE23rwVD8fA3HWNCwCAtnHuZKR27cARjYG+q6W+GCFbMffAFgEoGCC8C8Pa/CLRfPY0GojEYjDQYD9Xo99Xq92HeC/hkksPGssMiwVBCMCPUQQgRScWN4dkRbcLuWl5e1tramZrMZ19AzRgFXG41GXLNXlXJBcU6FQkFra2uam5tTvV7XysqK1tfXde3aNS0vL8dohpdyz8/Px7Rr97fR0l6rQEQDcI7GKWkhVL1en8mVQFNixXiOAsAhWonGsIzTCUbAOqCik0az/GDCA7hybwBCt0hciHkuBmNKy8YREgiT9FowOdiOJ2shEDx1nO9iWbG2CKbxeBwLwlhH0tW5lwOYaRQJwNgTqrAopCnTI8zIxQDzAJsgca7Vaun+/fva29uLvTXonNVoNHR4ePjKhk9zQaEzJtrY2FCxWNSNGzf0zjvv6Pr169HPRet4fQWl1ouLi+p0OjExCXNZmvZX4Pg7BAUay4uMPDeCTZviE86YXiAmnUU5AARhGO8riaDwIi+sEDpid7vdGOsn7OhaPQVXYT4yFxGGWCKY4h6mxV3AKvNcCCwGmsQwFqwR1h28RToTEoPBQIPBILqF1Iu4W0eCG0LDXRl3sVhLKk9JAScZiwgMlphXqGJNpM+vVqvpy1/+sm7evKlvf/vb+ta3vqW7d+/q4OAgzgsFcXh4+EpGRN54QVEsFmPxz3vvvRexCDIq0aiecIRGPD4+Vr1e18nJiQ4PD2NlI4AemwWQrVgsxrg9mkSa+v3SFOB0QsNhavMaGpWwIQzMWKnHQBDhWqBR8ef9EGAE33A4nMEsCINyDx+rm9zkU3huBozW7XYjSAk+Qvg2tQpcU3N/B0oZFxmsHqlBKPGMvNuVNBUECHu3mngfQdXpdGbcHqwyhKkLNCpTEdA8WwRdsVjUe++9p+XlZd28eVOffPKJPv/8cx0eHsamP1gwrxq9sYLi5s2bWlhYULfbVbPZ1Pvvv6+bN29qeXk5IuHetMVDY2hPB9ncNHWN7z0OiA4QTkO7pnUSnjQFkzj6joWBAOr1etHFcbANZkyjDlwn1XyEBsEuYAYiIszFE7ukaVYliUkIU0+aqtVqkQE8yclzKTwEjFnO90mIggEdWCSZyjEHXBmsJXeHCKV6vQ2/EchEcRBQ4B64gF6b4yFu1tT7dni6PZW3FBBev35dv/Ebv6H9/X3VajXdvn1be3t72tzc1KtEb5ygwMQjIWlpaUkffPCBbty4oZWVlbh5XDOkITrvoiRNtaA0zZyUphrXmXhubk6j0ShqMtB1zxNIE5acmbxEG6bya0uKeIKT+/NofHeHPOeBa8JYHplgXp5fAfBInoLXpJDSTAUmAg5tzPuMjfVwV8+L6Dyt2/NRsGgc3/GEM9bGXSvvcOXFdghzcJRUULAGuDAAowhZruuEAHKrxHNt/tk/+2dReKVA9KtAz3uuxx1JXUljSadZlv2eEMKypL8r6R2ddeH+Y9kFJ4V9kQTYhk/PRiqXy7p165beeeedmCoNpoAL4ea3NI06wFhEPNLPpRaCI+hpDoIzf0ppbobXWLBR6XTtqdu8luZgsDGdEfh8GqFhnJ4h6q6R13F4qnMa3YChPWSaFqnxd7omMKiHLt26gjndNWDcPCe3phBEJGMxV6wuhB5/u0XpgopnjruCoEgbA0MIBwBa3scNeeedd/T5559HQbuysqK9vb3vea9fFr0Ii+Jfy7LM08t+WtI/ybLs50IIP33+/597Afd5JsI8J0SFeUth15e//GUtLy/PRCfIS0C4oGVhdjYlviqM46nILihgVBjCMzW97sF9dTfDpWnI0xkVBuYQHsbJ2KVpnYoTEYXBYBC1vof1cKuc4dDw7oI4c3vHqlLprJHL/Px8xBhwEZgTn3eBcZFbwzzBVS4K+8K8Hmrm2QOM4o6gILznBKCmd9zi2bgl4S4mgqtWq8WQKa6PA7ee+4JiQbhiqc3Nzel3/a7fpdFopIcPH8Y2Ba+boEjpRyX98Pnff0vS/6CXJChCODswBkFA0dbS0lLM22+1WtHsJ6PPN4czOtdks8AQboq7hnQ0X5r68ggSzwGgQjPLsqiZ2WTkd7jZLU3b72OWO4CIG8HYcGMYA4i+M6ekGTPeoy0pM/NZ1o5MVU7r4gBgCr0QvqyfV2x6ElYaAvY18ggFY/Dn4p/zpCmPxHBPGNoBXQc5HaBNFYekCBbjKuCiICTcPSTTFHcJvIJnNTc3p9u3b8eGQQ76uvX5Mul5BUUm6R+Hs7ND/79Zlv28pI3s/KSwLMs2QwjrzzvIZ6WbN2/G1u6NRkNf+cpXYlybwiPArWazGcub3RR1pvNIA5ZKv9+Pm9B9bpiIUKCk+H6WZfHgGZiqXq9rMpnEaIO7LYeHhzNJVK5FiTp4MhRMwmbndahQOKt92NnZib0uOAELpvfNLk3Tzhkj2pbeGqD+Ozs7sSEt3b+zLJupD4GRAUAhXufz3NfTtFPrjKxSmBhrAkuK0CVWFAqAzxPdIhSMVTccDqPS8Ga8jKlUKqnRaMRxezjW3QyOO/S8FNwWBMvNmzfV6/X0/vvva29vTx999JEmk4l+x+/4HfrGN77xIlnimel5BcXvy7Ls4bkw+JUQwoff7RfDJZ9mfuvWrXi47sLCgj744ANdu3YtbpKFhYWYjCNNTXs2qpdFV6vVx9D5NEHHzWVMbccI2IBkCQ4Gg9hLIoQQD9Lx8B/MjnZPoynORC7IHFj0BCHPuKzX69rZ2ZlhQtq2FYvFx/pilsvlWEhFRITcEUzvt99+WwsLCzo4OFCWnR3cw+HAaEiujwVEhiPJS9yPfBQAQMdosJQ8Td1TwiHC1icnJ9rf34/RE7Aoj6JQjIfAz7Is9hehQY4Dylhz7nbiCuKSkm1JjoqkGLL2UHsIQffv31epVNL6+ro6nY4+++yzV6pg7LkERZZlD89/b4cQfknSD0naCufnj4YQrkvafsJ3f17Sz0uXc5p5o9GIm4kj5TgmT5rGyqWpFnYfGOaTFEuxXYCgVSCPPkh6LHbvQmM0Gs20i5dm0fbxeByZEIZ1TQk5cCppJpEKITgej2Mo1pOvFhYWdO3atahtYVJMZEkxtwLXBqFK4RbmOH53qVSKbh3uzsHBga5du/aYZeB4CwAfwoQxsk5p9If1dByItWEtPFEMJiUhC9cDLIa18uMLqBVhjo5ROc7kY/KwsKSZPIvd3d2I0XiaPgciD4dDbW1tqVar6datWzo5OdEnn3yi9957T59++ulLd0Ge55DiqqRClmXd879/RNJ/KumXJf2EpJ87//33X8RAn2F8sfck2IQ0dRkIpcE0MC0bDW2OMIAJMctB89mobBK0NgLAhQXoPyCah9Hwz7E6OMYuBfBSRmNOmOJ8hg3qJeL47AitZrMZ8yQQCtJZMhit3Dwc6CAtjOV9IrgHnwFQdXAPJvLnxFp5VAQ3DesCC475Olbivjxz8/AoERvcUHABrBrPGXGwGwvE0/Ylzaw15ILEBRdYxtraWvyMF4R52jyJaWtra6rX6xGrWFtb097e3ktN734ei2JD0i+dL0pJ0t/Osuy/CyH8uqRfDCH8pKS7kv7o8w/ze6OFhYWY4Ua/CE999vRmmJG267gIhMc8FwBzkYftLgeuC1qKa7nmow6g2WxGFN83P4wMEJrG7b1fJfdycNWBPQ/18TkYmvE4CEkExIFRBKdrMywptLXXlXBPrBcvbee7Xv7OeB0H6na78T2uzThd+Kaux0UWh4d0yYXgeXkOA2N2y8xxCV+D1Ep04ZxSsVicydXA0vMQLOP0PdhqtdRoNNRut/W1r30t9jelruhl0PMcUvyppN99wet7kv7g8wzqeWllZSWi6uRCEBZz1wBpj5bxCAdaA5McjYs2g5kqlUoMBfK/VyB68hWbQZpqd7cMeL1Wq0XgzRmAIwCc+bCMpFlB4ZoVTcQ1ERh0BUcYHR4eRgDOwUQ2NN9lvLS7T/tBIoDTHA6P1riGRgiMx2Pt7OzE127dujUDDDNOvx5C290qhBGhY/aAl+/zXZLf2CMe9qSAzMfK80jDtG5NpNEnT2ZjfLhGuHTz82cnzWP5NJtN/eAP/qAePXqkb37zm7EDOs/oi6bXMjOz3+/r5s2b+uCDD2InKjaIm9Sk/kpTVJ+H6BrE/WBpVntwOI1n+5EV6UKCH1wg97vR0pREUxdBxyaySQHIwFa8otL9deaH6e8ZmHyGazGmVqullZUVhXDWFBam9PoHmIBxLS0tqVqtRqFD/cpoNFK73Y51GKy5CwonD7/SahDtT84DzxDLyCNQrLELVZiZZ0HiFVYN65dlWbQ6PCs07fDlkQ23thi/C1SvW2GMS0tLMbJEAZlHgiaTiZaWlrSxsaGVlZVo6f3Ij/yIhsOhdnZ2tLS0pPn5eT18+PB75Ijnp9dSUMzNzemDDz7QysrKDLjFAxwOh9Hv9bMycDU8muHalM0kTZkfSwIzmFRlvu++Of53o9GYsTbY1BQPecoxY0FLUqfAe5SVe2VraukcHBxoPB7PdKvC3aBlP1pvbm5O+/v7ESdB+Lm2xjz2+ojV1dW4Nl55yZmeAKGpkPDojTR7ODJHJnoExnM+mAvrynz9WbNOac0L88HqBEdxpeAALOThbp4dP2RVevQGF5SkLy/+8g7gJIaF8xDvZDKJ7fP+8B/+w/q1X/s1ff755xqPx7px48YXLixeK0FRKBS0srIScwM4vIZDWXhopVIpHi4Ls6O9CG1hkoNppAlIuBZpnQGCxzdUCnIhePjbtbwXchGC5LM0h8U1QEB0Op2YIcg90HieSi0pnggG0s+xAYCs3W5XCwsLMYwrnVlN1Wp1ptSbngu4WY6f4DIRokVowTwIF9YG5i2VzloOrqysqNVqxWMVEWKQu2wICS8/9+QtrBx3tRAAaaapR40QYJ6h6taKpJgLghWKYE0tHlcYKBT2CrgMFoO7KQjM69ev6/f//t+vVqulb3zjG3r06JGuXbumR48eXQYbXUivnaBoNpsx1IbA4GGzuQCx6HmJ5PdMw9Tfvwi04j3MY/d/HRlno3joDGbyTZ9iFSSFsYn99CvmgI/rmIBraI9IOEbDmjiwSW+JXq8XGY4fBCvRGkBBN9VhHKoucVMQwo4R+L1ZJ/zz9fX1mOvi32HN/RqsjTTttpUKZgSHR4n8mqnQcgzAn7fvAf+s59zwvLCAwJMYJxaRg6TD4VDXrl2L4VvWuVKpqN/vq1KpaGVlRbdv31a321W73Va73daNGzf06NGjLyR0+loJCh4QoUGy73hQbCz3N/H3wSsuQtMlzWxsj2TwGc+hcJ+W9zFN0XTStFBImj1Sj2ukY/PMwFTzEWZjbF6Dwn2IoqDB+/3+DDiKlTIcDqO1wjxB69NU5hDCjPAAcMXV4UTxWq0WP+8uixfbkapOXwcsPIQDlOY0uHDmh/VnbRzjSCM5PA/cQ3cL07Czk6fZO8DLZ7FIEJIXuSwIEM8MdiuV9QwhRLe2Vqup0+lEYfRF0GsjKEqlUtyMngQkTTW/+7f+YHElnLFTtD7937VhioxfhIpDvpmdsV2Du1+MRcQPc2N+WBUIoRTUYz6u1bk2FabSNFuSdOZqtTojrEIIMcPSmROcBxzDIyue0OSVn255uXleqVRmkqJSgcTaso4uVL1BD/PhmbsLgUDx6A3X5Zn4/6n14+vL9QhtU8PhZfyst1/Ho2oAvqenpzEcyz14JkdHR/GQo8nkrG1et9uNfTBQEpdJr42goN1Yt9vV2traTPYdmosQKeRmI7n+foCPNBu5SC0K11xsMt/40izYJc2evekWiCfhYK63Wq2ZdGPuh3vEZvTqVu7N+47Ae1/LEEI8XBkBgqAALPWO4SQrcdCPpLjGgKRsfAcMPeuSe5DI5nkaaFbW3i0a6CJrizVE4Ls7hKAAX3IrK42aeEEa92AOnh7vz9stD9YD1wvh5dfyKAjzBafAakzdQY92Oeaztramzc3NmIxFHsxl0WsjKPDp6/W63nrrLZXL5XgCFBuAdvJpPL5YLMaQFaCnp1fjA6duC+SCwwFO37hu3TgIKCkySrfbVQghHgYsKWpXcBQ+j1BB0ziDsrmohISh/bAb1oVQJB2kCoXCzPkiCCJcuUePHkVmoG4Bhma++N9o8WKxGE/loo7ENTcZpEQfyGs4PT2dWT+uj+Xk6+wp9R5x8rEBQPu9Xbg6kzoGwprzHYQO3+PZcD+yLN2qSveL40IoMwQCQibNSGU8ZHoOBgN1u13dvHlT7XZb29sXVku8EHotBMXq6qpWVlZULBb15S9/Wevr6/FgWCwDmEN6HG8AgILBECpk6BFCZGM52OUmsDSbueevoeVTIM0JgYX/SbIQOQC4R25JjEYjHRwcRE3jgCWaOT1iD+vCQ4LOKERQwHAIO3LoL+ApWYe4QqllAkOTbdnpdGYiKRzRCMLPGRto1xQ4lmZ7eDioifmP1ejMyPh4doSipSnjX4RfuVvmAoZxYB16CB1FRL4LrqFbkIzX8SmEL3uRSBxREI+UUDzW7/f127/921G4XCa9FoKCB1Wv12N2W6PRiEyJlqJ/ojcm4ftoLjYSWuOirs0pWs+m8Y0mzYJu0hRtT90atAkb2ysMPTTrpjqbw8N/3ojFQ4b4+6nGTlO/IZiIMaLVMJFZI3o5EPXw/IWL5uo4B8J1MplEgI7GL4RweV/SDKMhHBg7jMbzw13y/cG6AaTyHYQNwoC94e5b6pKk65UCk55b470oXCg75jIej9XtdiPoCt7DNebm5uKh0YVCIYLQN27c0GeffRaF7WXSlRcUS0tLsRU7AJwXZvHDQ8SkddAwxR4AujDl0RAepeCzvmlTSgHFVIg4UFconJ2AjWnp9Q38uEbCvEaDoSUdY8H05YfvcO2LUHgPEftaSooCJzXZEXCMjTnxPpmPoPj42lmWxRPbERBuDfgz8Psh5H3dHb/xz6fr7eP05y4panWEDhmh3CcFuVk3rsfaerTNQ9YuTPgOzX0ODw8VQlCtVpvZZ6w7/Ubn5ubi0YjValUbGxu6e/euKpWKGo1GbBXwounKCwrvSlUoFOLpWzAxmpkHhKb1JiVPCpe5WeqItWuDJ7kRTq4N/bq+iVO0201rrCIXZkRBKGiCSby2AF/XrSrul4K0kFsZfAbmS60qtyD4fMqgCB8sG0BKAFOK9rB2fDypwPBog7sK6fidKVnTJ0WC+IyDrlhruFN+Had03zhwzD0Qyjwnx1U8vI1C8rAwzwJrkb4oHiEDzKQfx2Qyie7di6QrLShgDDIHpalkdxPYC77Qlmx+z0nwze1aydFqyDfnRQwCudXgvq9vXB8P1oQLEr+OMzIWBEzqYVTXhIzxIquJNfPIA8Cua0S3ZjzkJ01dJ7e2HEvhO4yJPg9obbd2nHF8Hbh2mg7v7zM3D226MPf0e18TSTMZrCcnJ/F8Eq8z4fMuCLieWxs8C36Tnu3hcM8hATx3q8qtQrfSHADt9XpqtVpaXl6Op86trq7mgiKlmzdvan19Pab6AkC6S+H4gJvsPDDfmCnj8JqXKEOOU/iG889cJDSe9L4DdKnZetG1uJdrOwdL8cE9pJh+N2U2mIkkJdaxUDgrrmo0GjH70MfEfdwqcByFNWeMZIByjYssA8bsAhiL0cFn1+oUcXm6N/f0Xhcwqgs3F9qEqDmHw90+t0xdIbnlmt7bI1xuvfjcsDw8wSzLsjhfQGVcI1LGG42GNjY21G63dXh4eGkh0istKHgI7scTnaBbEtLaT5lKffTU54bcp+QBSrMM78zkWj91U/y3I+WuDVOmwcrxTZYCYswxFXp8P52TX8/nyW/Wx3MPQO9pnOtM79eGGbEk+J8KXh9DSm7B4a6k+ATRC3edfD2wklh/yLU5cyJ07IIS7T8ej9VqtXRwcKDDw8OZc2eZG2sIwwNC4va5u8MzSC0EDwNjCbIWCAfWnqQr1kBSVF5Es9J5v0i60oKC3ImDgwOdnJzEA3xIAKIQiTAnmyBlVGdWtIKDZenGlGbNV2kKprnmuehz/llpulm4NviDa0oXUm5Ku8AgmuNYChvY/XW3mlJ3yufs1gZ5J2hg968R1KwrzAqe4a3zEWSuvblXagH6WvM9nncaZXITP8uyaMmAUfEc0mfo400ZlffSWhUsEV8vBCIhakBoz0ZFQIFV8PlOpxPzZ9jTWEwewWHNPGu1XC5HS47I0dHRkd5++219/vnn35l5vke6soKCmn2Kh1i4FM3HRCON+CJUHNMyZfSU6Z0uAu1ce7gPDznwlzKqM44zQHpvF0wIlVTz+jVTjAQhyDVcCHEdZyi3mFxrO/DoQsWFHdYcWtjnzvz8XggEjxR4G/2joyMdHh5GnAOLkgQxCKHh5jtjdtwmdQURLIyFCA/Kh/dSoci4OV6Q0nVJM+6IWzxYNXx+MpkmzzEOB6ddmXEtadrZnUbAdPx+0XRlBQUbEESd5KSLzGqkNCZeCnaiIT1/wJnsIk3ngsA/5wCpux4Osl507dTnTgUX9+EaDrg6vsF8/Tv8dvzions6wOka3vMK+DsVYinOgAWQpjmnz8XdNn4YP8ln3W5XvV4vltnz/DnYCUbyebr7lrqAjN/xFBfyaZp5ihvxPHENuaYfC+Brz9r5+rhLO5lM1O12Yzib7FnPifE1x2oCR0otjFxQGAHWIfFpq54KAj57UTKNNO0vidBxDec+Jq9BzuwpueUgzaLrT/IjXRu7IPL7uJ/PdSHfdH4atpv5RB6ceVI3i7H474uiCy4YnBG5p2eOkh/CtTy6kYKVzgSkJe/u7urw8HDG9ZEUa2Ge5Fql1pALp/S5+Bgcf0FwoUiwlNgv/lMoTKtncUnSa/mzJfrTbrdjFzDCoI6L+F6VphEa32vpur5ourKCwheD8mTAnYs0rzTVAjzsVMv4Z/m8v+6b0TeRb1Rpallwrka6uVIB45vbx5Myp28K/x6b0QWJZ3oC3noNBGNm07pL9p0EJlrrIq2ZukAI8/F4/JjW93G4uyGdmdMcJPTw4UPt7+/HorKVlZUI8hFKRgAiRDxilb7uCoN19UiYC2VeIwHL8RDe9+Q+QpyOFSEs+R6EoOAwKkKxlJqTZZsKGLdMwWwAa1MA9UXS87Tr/36dHUYMvSvpP5bUkvQnJe2cv/7nsyz7h896nycR50qQRszC+qZwDVEoTA8qZpN7EpM0Cyy6VeAbX3r8lPHz9XjM5OWHzQhI6RvPLYB0w7oPT0KNJ065hcScMPd7vV6sF2FTe1SEeTAn/P0UVPT8BbRq6k7xPsClp4kD7l2UY+H34fqDwUA7Ozv65JNP9PDhQ7Xb7VjYBhPR1r7ZbMZWeT6Xi9xCgGyUhTO6RzRYP+bBOvua4ZLwjMDGvCI2xViGw2HM9XFBXalUdPPmTa2urkZXAnfarScXPh5ZYqzMg7qgF03P04X7W5K+KkkhhKKkB5J+SdK/L+mvZln2l17EAJ9EIQQtLy9reXk5nkTuBVGuoYh6wLCSZhjU/e/UcuB9vy/fB7RziyK9vzMCjJRqKP8cjOVai43l4Bc/jMXdLQd5YX60jqQZDezzCmHawdotKXdf2KTOVG6lIRwRwi6k3FLxSA9AZb/f18HBge7fv6+7d+/GBjqtViuG/ygwu3btmtbX16Nbw1jR7jxLBCnCG6yKe2NBEaEhFOnZu9T74P8DbLJurFfaGEialur7M3fl4srGy/VT65Y19F4k3A8lwl54kkv8PPSiXI8/KOnbWZZ9fhmDdCoWi1peXtbJyUmMdviG97oIBAD5FWiy1FR3YCvFITxUmfrm0uPp2OnDdYQcn9X7FTig5iXlLhg8auBdn1wwuQXgvrRvRj7jJvCT5uSgJf+nYV2uj4ZOIzkpCJeOh3b//X5fnU5HBwcH2t7e1ubmZjwGgCa+WBshBC0tLcUO4C4sL3LRvErUAWzHFRAG9Hyg6S3WAM16PI+C9blIe/t6Ox7BOmIF8H8amnXLJxXMXgmLVePu1GW4HdKLExQ/Junv2P8/FUL445J+Q9KfzbLs4AXdJ2opb2/mjOy1Dm6q85BT5nKtd9Eie1pv6qOmGsEp9XF5yDzYo6OjmbNCuJ4DkR6jdxeAjs4uLFIAyzWdM7xvYhdu7gr4WqRaEErRexfMjkE4ZuQWXqFQiEzYbre1v7+vvb29WDJP9iaClRRlXI61tbXYvt7dSZ5z6hamc/Bxp0yIheMHHLvwTSMfvvdcQLhr5evM+vG3YzRc06MajuUgEPih2bInEr6SYGYIoSzp/yTpZ85f+muS/qLOTjr/i5L+sqQ/ccH3nvmQ4jQjMNUWHoOGER0ISjeMdLHJxkODqd1MBXBKtXGqSdytcdOdTZhlWfS/aQPnmt03FuNjTPju3NfdmJRJ3DIgDMrGRIP52vj6StPNnGa0phaU39PX1+fgOE6v19P+/r4ODg7U7/djla4/N8ZXLpe1vLysjY0NNZvNGOW6CExOhVVKKWP7fByM9jZz7layP/ycFT/6gev5XmKP+X3SdUQ5DAYDdTqdmRR4ngOChCMRuDbP55XCKIz+TUn/S5ZlW5LEb0kKIfx1Sf/goi9lz3hIMWYcvrQTm8WLo7wlu3/G/3a8wTcb5m63243ZfgBN0tTXf5IlgrXCuH1jHR0daWdnJ3anWlpamkkhdgZDwLhpjdYDyPN6AncbGB8bDevIE4FS6yF1Idi8CGG3Tux5Xng9ZwAELWPsdrs6ODhQu92O/SCJHjBvBxsbjYbefvttXbt2LTahdUGIQHZB4S6jj9ndShcoHiVywen4gLsGWDunp6ex+a1XnLowSIWWKyc+i8XY7Xa1v7+vLMtijgSfQZDs7e3F4xxSN+9F04sQFD8uczvC+Unm5//+EUlffwH3iBRCmAlBOUO4metmcGqG+waWZjc2m4v2byT8pMKEtnsppuF/e5Whm/39fj+a1Lu7uzENnQgOjOJujm8412gcZpzOV5pGH9wUddDPNxeb/6KEtUKhENvXeWk060FolTFf5NIAIoYQ4jkfW1tb2t7ejnU5KfgKBgXWsL6+rhs3bsw03U2tCQ/zumDw0Cj+vQuUEMJMFmmWZbG/Ca+R9eiFh2h+T+jzxD3G6O6D40nujkhn0bxer6dOpxNzKwCifZ9iUZTLZwcl8fz9cy+SnktQhBAWJf0bkv6UvfyfhxC+qjPX407y3gsjmN0tBelxE8+ZJNW40uxxey7ROQCHpqXe18F7R6Saw4WWm+KMiySb5eXlOLYHDx5oZ2dHCwsLWlpaiidEcQ0HYGGaYrE40+PAhaEzuwtU1sA3s0dKXPAeHR3F+bgARvAwJ3frPNycMgZWDyBtv9+P6+HukjM2iXSEQZeWlmb6nboV4Oa3Rw94DWuG+wEuMkfWEbcyxa9oKjwajbSzs6MQQgRZAZqxQIvFYsyuTF1ehKDP059ZuVzW0tKSFhcXI4gKPuMRm2q1quvXr8d1Bf96JQVFlmUDSSvJa//ec43o6feMmsl7BV60OVlUNzNPTk5i6zBMTMx3l9ScpIWPSESCWLVret9QrtVdYDEOtDP/M4bPPvssAnoOCLqmR1iBM+CqIDRgXE8y8h+EDGE+NKw0NfOded1fpzEua8icUi3q1tNF+AX39OazKQBI+nKxWNS1a9die8NSqRST2DwcyBgYE2vGs4EJWWvW3oUIY/R6Dt4rlUpqNBoaDofa29vT3bt3tbe3p7ffflvXr1+PGh3BTbsDVxiem8EaujuJ8CX0imXJaXcA33yfdSJZi8OhnpQy/7x0JTMzAQCJvR8dHUVgi02SmtEwF/UCPDgeCu4M5byEyxBKzmwIDWk2K9LNbj7r4BU/xeK0E/VoNNLq6qqOj4+1u7sbhRMMjJBi07pP626Sa1c2vPvKEEKFBrkhTA8ScgHDeqWAHHN2s535elKVPytpttU+8wOT4H2YkpPEFxYWtLKyEntXpGCgR3M8OuBzduAb14bIAYKaZ+3RBNY1hBBD8Ow3LEC0P6nkXsGb7gNcH9aMnBPWDXcWa4Pnwm/cVcLjp6enMye/c94rLtuLpisnKFigmzdvxq7Onl+Pieb5BtK0HyKCJQ1TOdjIRuJ/zHJMT9dKMJsLhRTsY0OiOaRpH0TpjAE6nU50dzi4NmXEi8xYLAuvg3AQ1TEKvlcqldTv96MLwxjcreCe/E6Z37U373s0xbEQ/7zPAzzBs1dpkNxoNGZOIPPwJ9/3tfY5p6AhbsVFURC0O3kTLiz8rBDWb35+PiaANRqNmWeeHvHgCiLFJbiury0gJeuYZpJiaRUK09CyNG3Yw2cPDl5YNkKkKycoJpOJ2u22vvSlL8WWZX4wCklNaAGYHNwBUy31iR2fQCsQIUixDmk2fdrBxovMWZjHmY2NVa/XFUKIAg+BllbCpn4nrpW7W34fxyYYu2s6gFA/0uAiS8BdiHQersHTaANr6wzjuQoOSsPEi4uLWl1dVavVeuzQZccZXAD7b193XnfhypjTOaAA6P3AbzfjSZBbXl7We++9FwU/pehgFu7m+Pq5UGP9L8r9cLfMQ6M03XXc46JQ+mRyVon6ounKCQpJEaEm4QQB4N2F8Fth+MFgEE0zmAqNyuc9mYViHhY/zbXnfWm2pyZaVHr8NCq0pjMrBxURDaGnQdocJWVm30hsOncbuD9mqTMQm7zb7UbcxTGJNCOQTc5c3EpgbKlF40KH8fb7fQ2Hw1gzMjc3Fw/hLZfLMyeYc02Y1SMrKYDM9V37uxCDef3zLsjdhRuNRnE/4eIi8Mjj4CgIXBaugRvra4lwoGwd14KsT67tOBdC/Pj4OLpgWFY8n1KppKWlpWgFXQYu4XQlBYU0rWdwMz/1qR2X6HQ6kQHx+YvFYsQ2pKmfjGVB9idMICnmQbj14g1xGQtawn1h8AE0sI/92rVrOjk50YcffhhDswCtfjq4Cxr377m3l5B3u90YXvPCpsnkrAmrN3uRpsfquSXiZj1CiTmg4Vgf5uyCsN1uR6ZxIcE5IOBD4BJeu8GaQ/SJ8PXl2lgrjUYjPnsHOrEqHfymDICqTVwQ0rjJUVheXo77xLNowZgYB7iQrzN71TEYrBYH1dmTRF3YVwhRWjli+fAaTYDdhbwMupKCYjKZ6J/+03+qGzduRHM9BTBhwuFwqMPDw4hLDAaDeOYmoUoeHiYfzHt6enbmJgg8Jj2bzYufpNncCbLn0DQewWBDEH1A67/zzjtqt9va29uLCUjHx8dx8zvmglaSFJmVFmowBvN36wDXhugNm84tK7S850owLxg1y7JoUbmAcBPY27iRE8BpYySyNZvN2J3Jk8AAdVlXrukZrmhkwrk852q1Gq08gOfhcBj3jwttFAunyJVKZ13dERZEEyTNdLnyLtk8D8dpvHgL69VdIIQY53NkWRZPSiOhcDweR2FeKBSim40l2m634/Ongjq15F4UXUlBAbFZ8BXTngeYd2S6dbtd1Wq1uLn8cCB3LbBEPNzJNb1yk2KzNEkJjdHr9aKwIYTmKbaeQkx47fu+7/v0W7/1W9rd3VW73Y7gFYcEe0aqWyRpvgTr4wDvaDSKbgg9H1ZWVrS0tBTzCpi/4wzu7rAOrBEmOOPnfRgB92I8HscEKgTm4uKiWq1WbL7ruBGM7HOBMVkHBx83NzdnwF+0MS4cp3V5yJxEpvF4HPcPwqzT6ahararX60WBj8uD0E9BUd+TJEwxh6WlpVjB7NGpLMvU7XajQPDOVuwT8k7A4qRp5/NarTazTpdFV1pQkC13cHCgbrcbMyUxId1PZfNWq9WZ3pqYrLgbfggM1YlpjkSpVIpneXroDi0qKeZhECXxrEPu6wg+m7her2tlZUXdbjf69DRMhYG8xJx74w5I027Y/X4/CgWug7m9vb2twWAQ12IymUT/2kEy5u2RCdYLAcWzQNM5aIcF02g0IkBJD8t6vT7Td8FxD8cj/DevY/2xB46Pj7WxsaHRaDQzBtbJgcbxeBzD4B798vTtRqMR188BS0x8z4twbCGEEKNXFJQRDs+yLP4mMQxcjegaY3bchOMGcZ1QaFhNgLEXlTW8KLrSggKTdm9vT3t7e1pdXY0mXqfTiX+TKcdmdx8TxgbMxBx1X5HNw0Yql8uq1WrRGsEER5vMzc1F39VDb2iNi6IMAGKnp6daXl6OFgWZgHNzczFlF82MhUNOCcQ42u12dBEQgmjQ09Ozc0SJvzvGAqVRB4Sb94V0K8AtHgeA5+bm1Gw2VS6XZwA8NKxbcqlA8Ot7aBjrAm3rGaAe2vSoxUXRl263G4WJn9TlNSd83zW84xNeTeyM7BYdiouzYl3xgHU5tuE4GSFTz/HwFHTHaC6LrrSg4MEcHBzowYMHarVaeuutt1QoFKKfDXODR3g4CVyCBeY9JDUPw9OeQb95DwbkesPhMAJPLu3ZDNLjWINbKkdHR6rVavFcCbTIw4cPYyOXxcXFqDUXFhZUrVZn8j/Q4p1OZybMi9AE4Gw2m3F8MMVkMplpweZhTr++uzke5sMVBLkPIcTKWNbI3T5PgYZRPFzIGrFePAuEDgx1dHSkXq8XBSpuD9d0PAJGzLIsZj46o/I3+8a/g6U2Ho9jMZi3YPQwtUekTk9PZ04M4z4Ics/WRbBg9bCmvM+1XVhzn8uiKy0o2Kjdbld37tzReDzWxsaGWq3WDJjoITPflLgGnmkJsu73YCN55iIgIYzm4TWEAp2WwEo84ceBORiGjTY/P6/V1dWIvBcKhdjQZWFhYaaEHusmhBDDvw7MusUkTc9q9ZCw4xi4GXyP8XEQDZrxonoSmNbXpVgsanV1daYTk5+q7gwMEUr0qJEDqm7B8Hqn04m+/cLCQkxaq9VqsQW+W0wIOoQz1iPv8bzcyvE9xPN2SwO3FyHj1amuCNhnKBxKBgCQUTjsMa7tYDnXZLw+r8ugKy8oMLsODw91cnKizc1NDYfDqMnwH0HUe71eNHvZrJh4RCKoFl1dXZ3JTcCEh9gs0nQDcS1aoLExXBsjlChXd9N7fn4+huXW19d1fHys+/fv6+bNm7p79260Evj+eDyOfq8DsikWAojm7lG/349oP4zlroNjOAg7jzZgFUwmk6jFEVRs3GazqV6vF0FfNCKWhW94vx/YEIwI/uGmd7/fV7/f197enjY3N3V0dKQvfelLcX64oIwLABLr0F1TBCPPmWfIAUCsKX/7MwMDYp0Q1F4G4GFzgFT22s7OjprNZlwfj4ggNHd2drS8vKxWqzXj+jHmVGi9aLqygqJWq0mSVldXoyZA2x4eHmowGMSTstE6ntrt4JmkGDrd39+PdRBYBvSM8CgAm8w3GMyEWetZoCRHwSieYShNIypsDKwUGG13dzeCdWxGckI87OoaB/cE0Oz09DQeoEMxEy3lMGXpkXFycjJzBIIDedeuXYtjPTo6ijktMICDnZjunU4nztXBRte8bHZcD+bAfBxDKhQKsY1eqVTSO++8o3K5rBs3bkSXypv74JagMNzK29zcnBkTwpA9gYAFj0BZOIDN52hG9OGHH0bwnEpQcCEEMtjW+vp6LE1wFxJ3BYD5wYMH2tra0jvvvBP3IhgN4Oll0ZUVFJVKRe+++67efvvtaOKjNenMjCaRFNF1wnloTAegut1ufDh0M8ZikGaLoaQpmOobyTVkr9ebES6kAcN0/p1KpRJTbxnn6emptra2tLu7G3MBcJPSeDla0hOCTk9PY1hWUozAvPXWW/r+7//+iEXw3sHBgR4+fKjJZBKFRLPZVKFQULvd1vHxsebn57W7uxtdsSzL4uZ34JckN3phYlYvLCzEhsguKB0fwtrgGUnTIjisDSwgrKmbN2/GKNV4fNaFvN/vK4Sg+/fvx67tHmrkvp5bkgLM7CuwKP7meTo2AZNT6Ed0B3fLFQvCrlqtRmsM19F/cNNarVZUGh999JHefffdKPgRmAjNy6ArJyiKxWLs2XD9+vUZ7UaorFarRbMcJoVJDg4OND8/PxP2lKa5D7glWCdokYvMzd3dXXU6HV2/fn0m5OphTAecnDE8U1KaWhSYnmhQrAlqILgWQsx9V3xeR+EpPPP8iFu3bml5eTniAGAhnuVH5qCkaGmMRiPV6/WZdnWSoqvmvSXAM7waEqzHMY7UZXJGdd/bIzDgMf1+X+12W48ePVKv11Oz2VSn04n+PbgPjYEKhUIMabtQpc0AgpvohLtGWHC4luxFFJBbieAUCBRC6czJLSbmT3iT62LhIXwODg6iu3RycqKtra34bABzsZ4ug66coADpR+PBNDAmoU02LD9oBhjfq0s9jp2Gx8i+dASeaAfdmdbX1+PDxRWQzpJ+SEZKw35pKJDXGYuH8RBobsl4+A3N5/iJpLhpYFCyDqlVwApAWIZwlrjmyVkwHbgPjOV5KClW49dCYOPGeZKSW1tp+NOfG5/lb8xy3IeTkxN1Oh212231ej1JU6sM94V8l2q1Gl1S9hNCizmDgXg+CeA3rg9CCEsHIe9uBYcWkQXsGbMInjSBzsOmCChcJc/hYY945zEU1GXQlRMUPDBCnd1ud+ag4uXl5cc2mTQtNfY272x8HrQ3pEkTm5zB06QXj0Dg62L2E/qCKQnbObgIs7ApXAiUy2Vdv35dkiK+cNH8pKkmYg6Y2ouLi9GEdSuC8F+tVtPa2loM/T148CBqaVoBLiwsqNFoxP4QuEbgGfjVREWojKV2g1CiJz2xtg4Yu7vB56RpJqWkmdT4arWq5eVlSdLm5qbm5uYi/kJ4nDCtW5Hcz+/r1iVZt/45F2q4nJ6khWIZDodxn3jGJ8IKYcO68RwQhm6R+j7B8gIPGo1Gj+WipC7pi6IrJyiIfVM70Gg0YuMQNhUb0usw2CBUZ/IgSIDBbHM/m3Cbx/oB6ubm5rS2tjYT7vOUboArag8wb2Eo1yK4BM440lko8+23347hX4qIvLrU8wvcssIVoDqUArP19fWY48D96fvA/+Sh7O7uxs7Y4Ab37t2LQg8Qr1icHszD5xqNRiyoouALoJdoBjhQGn50UNHfY10RoI1GI+YaFIvFWP6N4PB8DzQyjOcArZ/bgfXIc2XPYJ1hzbD2KBMS8Xi+6bhDCKrVanEtWTtcs1qtFpWgWxrl8tmhxR5+JqrjLjHv566HpoAdJicAD5LXAUt8RIgNLM32AmDBCYexiRAsgG74vNK05+ba2lr8m4fFJq7ValGjoFUce0hzERz/wAcOIcyAVIBuLmjY3I6UuzvgoB3VhoeHh5HJOAx4PB7HOT58+FCFwtlpXKw7DM8cYQ43s/f392MCUqvViqAl6+np3x71cbzG3Q6EvYf9WG/C2Y1GIzaQuX37tjqdTmTMNLvSE94YC/49QoO9w3MDD0lzFS7K5UAgEfnAzWH+Dmoj4Obmzjqwk6vCOR3sEYrluC9uIa4g4XaA6TQn5UXRlRIUmG8wPczsQJ4nG3lkQ1IMUfI6mpsN5P4zklpSrOrzVG2sEgq/GJ8X67i14c1yYXT8dtc+pOw+evRIm5ub2t/f1/HxsQ4ODuLnOXPVD7ItFouxvoS5MhaYmyMCDg8P1Ww2NRqNtL+/L0kReFtfX9e7774bha6b+MViMWpBGIPWawgIxghW4yY2jMJYPa+B9fdoh1tbMJq7edKZ1XXt2rX4XVLFsVokPVb4xz0AXPks91pcXIxt9730HgwDzc4YcX29Hd36+vpML01CmQhBlAlKxF1GhIFnsbJXAc09L4c5YtFcBl0pQcGi+cZiM8OU1FQAIrqmxzRDILCR2Kyeis0mxQ/0iARaEUECtsA93Iz1yIpv+pRgqE6no52dHT169Cia/pjG3kuS1G1PpgL5ZjNSnIS5C6NwaAw9IZrNZtT+DvIydrALSXGNmDtrga+NRk6BOjqMFYvFGIkhQoWf7W4ixNxcsHsJPJEqrMeLIhH8eJ4LbponU1GuDdaEJQlDcoYHxWSUzKMMOP+lWDyrQEVZOCBOTkS9Xtfi4uLMnsalJd8HLMbBeUnRgmB/sn8v64Bi6bsQFCGEvyHp35a0nWXZV85fW9bZSebv6Kwl/x/Lzo8NDCH8jKSflDSW9B9mWfaPXtRgPcTmYUTfDCGEiMp7YpJHK9iY0rTWAg2IRIbpB4NBrK3A9WFjekjMw6f+sBzpd0AMBoQ8cWd/f1+Hh4cz8X6So5rNZiyhBl9xi8nNYnpxoAk9UxCNjGW0vLystbW16P4g6DzagCb0UK8XeqVrgVvCvLrdrsbjs9ocgEK0JbgH5Ovk68jzZ80AAh3153tptMmjYB5pghkdwAYsJGmOH2mqvel2RY4E6wTTsleYm1tUXIe5QuxF5srz9L3CZzyD1cPil0HfjUXxNyX9PyX9l/baT0v6J1mW/VwI4afP//9zIYTfqbNzSH9A0g1J/30I4fuyLHsh1SoICkeeWWTXEpxu7mm1hJMuyl/AlPWsSq93ADziITtImoZOubaHbbmPM4KPAeb0JCU/U7LRaOjatWtaWlqKmggk3TESR72zLIufZcPTiAVB4Q1sOXwIAYv/Lc32ovT8B9YDa8GjMJjcWDmEMRFemOue9JWGI9P1coHr64+gQHN7KBaBALP5+jBOhIRbm17161YiVhy9NLAKWDOPgrg1xBz8vuzj1PrycbmQ8j3kPVG88ZCv2YukpwqKLMv+pxDCO8nLPyrph8///luS/gdJf+789V/Ismwk6bMQwieSfkjSr76IwTpqPRqNYp0+bgIxZsxZ0rBhFGm2CtDBrUqlEje1P3RMQU/AcZMck9EtCt+Ubv6mfjIbBRPy8PBQW1tbOjg4iOE3Tyyjf4P3BpWm1Y4Q4yiXy7GylHoEuiRx2BDCgiy/lZWVmZoPTOM0Q5VwMm4elgnr6SYz/j9Cud1ua2dnJ2Z0+jq5tejWkT+7NKxNVKDT6cT78eww/f07zrx+HSwi8AGEBEC5R5GYP2Fb8i9YAw+hcy9nYrcq2b/SVGn43L04ECzOLQi/7ksTFE+gjez82MAsyzZDCOvnr9+U9Gv2ufvnrz1G4RkOKfaY9mQy0cHBQSwTlhSBJEmq1+uq1+szJr0j3S7xC4WzJC5i/Two7oWlgRmJBAf95/MAUp7O7BrMcx+c8Ug53tzc1Pb2dkS0a7WaVlZWdP369egLez4EGwxtmBKhNTS7Ww5zc3Oq1+vRz0abekQlxV1SPziEEJmH07+dMdwNIAWdZxXCWU/PTqcTI0EnJyczUSxfN2kabXBgtFAoxLWRFJ8/wkmaulc8MwQgGI0LD0/q8twLD5l7zgtKiPsiaGjLx/pgefg+wNXFBXOhyLUGg0FM0wYA9XVlnfj/Zboe3wtdlO1xYbwme4ZDih10Ojk5iS3uyVFAA/T7/QgoIX0rlUqM5YMvuFZAaziaTDgKjQuDkenHpgPzYEOQkJRaDWxQ9zlxOba2tnTnzp3YjHZpaUkbGxva2NiIOQluSruZ7/kTbGqEqdeAANYtLS3NJHeB68DwCEbfjDCFm8Cg/FgWYBU+LiwbxiVNw6FuCfiJ8b6GTszdrTjXzK1Wa+bZIpTIZITB3PqTLtbCCHsHjLkXzOuJZgg6OomhoE5Opg18PQ8EVxAFRI2R73X2Kun0CGLHIsAymK83MHqR9KyCYiucH0YcQrguafv89fuS3rLP3ZL08HkGmBILk+ICad6Eg1JsUkxMUr/d5MQycJ+UzYQP6WCVCwrHLsAP0FIuLFLyBjNUrUqKmoPMRoDUVPD4mrj57iCkJzexIQHoCoVCNKN9w7uvnd6De+OecVgzzJpmizro6dqPRCcyGd2aa7VaceM74CjNnhLvz57xYVn4Z/k+YU13q1gn/mZPpOOH8akv8e7iDt4WCgWtrq7ORDk8j4K1wxIhR4NzVR305PruEvV6vSiYGZuDmen6vyh6VkHxy5J+QtLPnf/++/b63w4h/BWdgZnvS/rnzztICCEhKQJKDiyyGRyIch+YDeZaRZqatGQv8jmYyYuZuE96yK2Donze/Ufu576oRzk4wp4cDNyC1Bd1S8XBK3dzIPfJy+XyzGb2bFLHONBYjhNI0wxQzytg45M4xAbnh+tBLri5lmvU09PTeIIZwjxF8j3ikQoJnoH7+S7I3e3gGs6UbmW4S8kcidy02211Op3HlESpVIrhVWd4X09/9t1uN+a08Ax83fzagMZgTQ6E+rpfFn034dG/ozPgcjWEcF/Sf6IzAfGLIYSflHRX0h89n9w3Qgi/KOmbkk4l/ZkXFfGQFN2NpaUlSZrxZZ3SDMXzecQN4ovqmw2GkqZt3WAqD9/x8MncY1M4LnHRg3PtRAJPu93W9va22u12rDFoNpsx7OZYilsLzhipb5r6uc5snvzFWrnmQrjwOpsXrcv8XaNjgbEGqe/sbgJ5DuPxWN1uN/ZQgFFgIOks+QxgmoiIW2e+nswvzSPwiEcq9HzMjoX4d2FqqmuJSlF8hhsKfkEoPX32/r8LHrJ3Hdj0uflz57m4e3FRxOQy6LuJevz4E976g0/4/M9K+tnnGdSTCMbinAsebvqQPeGHjcIGdhAx3VSpCe9NQ1KwyTWFa2TXHD5GN1OJchwcHGhzc1Obm5vq9/uaTCYz1Y2eaJMKIe7lQiLdZOlmwzriOgC7qYmcApo+byyT9B6+rqmQdEGB5qWI7+DgYCbPhdyLfr+vbrerlZWVeMRgmoiVAp2MPX1PUgyPw3DpM0r/9t8ICu+sTsMe8Ad+s1cQPhflb3iGJyXoF1lgvi+J4Hmo3/eSWxaXQVcqM9MJ35awlDQbp2YBeXjOTEjuNCqB+fwkE99DntJsBiHXZhy+WdNw1/Hxsfb29vTpp5/q3r17arfbcSwelfBEJLci3Dpi7MzbsQwwCt88MCrAIeNGOKK1uL+b637fdGMiIHFN3F1zF4nxk524v78/c9QjP1mWaXd3V4eHh1pfX9fa2lo8jMktAJ87mpc1d1eTZ+X4DnNFyPk4030EyMvewZJIcSTcOebhVon3M51MpkV7vva+nvxGUGDVIbi8NaH3/rgMunKCAq00mUwbroBFSFMAC+AKYcF30s2Qbv40y9PJH6Q0a5GkoKXXA3APNl2/39edO3d07969GPsP4Qyp39jY0PLy8kwLODYUKLsLCRdcRDWoHiVlHeZxJiGZLMVR2GzuQoDos/YwFhYIlpuDuKwB93OBy/tra2txHltbW+p2uzPl26PRKPbD3NjY0Hvvvadbt27FEC1CyQWwR5UcjPa+IG6ZudDw58ezBe9yYQIATGjZo02TySSmcnutDWAkRX7gX+Tn8Mw8Mxj32UF41gW8hDESbclrPYyQ6IPBQJ1OZ8afxyrA4mDTpZl3F2EU/J8KhIvAotR8vciN4Tcbl3DiwcFBbC3Hxi0Wz7pVr6+vx5Oz0JwO+jlw5vcpl8szQJqj9272OvbwJP88FUSeuISJz3w8jyPFSlgvQrEAswifYrGolZWVOJ7hcBjL1emhQUu7hw8fRhzgrbfeiriFd6tGgDAvxigpvpcCnKkwY++4wmEtPaKQNhSCsWHawWAQrTZ39ei0RgIXgp37M4YUU8NC49kzNywMF2aXQVdSUEwmk5hf7yaYNEX/vSQ8LWPmNw+DB+1aNzXzXdswhvSBpmN0kI1oDFqv2WxGrUgE58aNG7px48ZM5y5pFntI/e8sy6KQ8I3vgszHwNwdlWdN+E4qSPg+QpjrsmndAkn9cpjArQzWx6stKU0/Pp4e/1goFCKznZ6e9Q8l72RjYyMeO4DgcvKeHWTXevYoY3F3FXJ31MFE5u/JVay395mAcRESCEjcFQeisU4uwijcPWaOCKeTk5OI24Dz5ILC6OTkRNvb23r//ffjJiM5yn1ugEjfvNIUP3Az2MNWqd+dvp6CfGm0gYeZhuk8jXlxcVHXr1/X0tKSJpNJTEOmvwYMz7z8Xqm/zwZE48OIgI4puZBxAeECKLUwmEfqtnGvNLmMa6dYRXo9X2uYipJ1NLJX8o7H43hI78nJSUw/h4mxBhl/6gY59sReYP7+v4/fAWOPrHgtEDiBJ6nhJnj5ePo8sMQ8gSzdXy6cPWu0UJieR4ugcMX0ounKCYrJZKJ2u/1Y2naKxHsatedOeEZiavL532k04yLtC3E9t1qcARzMZFweWkWz8CNNwUKvWPScEcaAW+CCxYFXDwumiH4ajrsIc0mFiVslbon5NVywsv5utvOeV2nCMFT+Yi2iJVlHxyVWV1fjmaZgOmlNCvkw0uPgc2ohXrTX0rXwcUiK2I2DnG6FpIly/vzcGr1ISKeCzsfkIVFvqefJYi+SrpygSIlOQalEdnMR5nT/nM3r54sSPoNSlyJ9yJ5bcJEmTi0Bjx4QfnUrAaZyxuMHAeJj8nF7Z2jPm7goz8HdC+6Tzk2aVooSBZEUz6+ACVygFQrT7l18BgwDHx/N56FYaZrWPRgM4vfR0Kmw2N7ejj0jvNQ+hGkrfdfibmmkmJRnbzJ3BxLTNcGaweUFFwJIpnzeXdZ0H/jaS3rskKCUHE9hDSiKLJVKMbSeluq/SLqyggIg5+joKDZ3WVxcjP4+WZswDloAyeu5At6s5aIfaVa6OyjIe6mLwPVhZqwZEoe8fZkj7J5a7Sel831pig2kJjD+LpqNcyVg5MlkEpOD0ISpNeDAo7t0LmiYj29qohVuabDW1DGQP4AQQZiVy2UtLS3F1gCeTk79gme0djod7e3tqdlsamlpKQoo1tErQd3K8ufnAK+kx8bNPRk/60waNQ2NSCADmPTO7SlAzr5grI5zucWSPgs+x/4A+C0UCjN9Pk5Pz9oYXgZdSUGRZVmsRDw9PVW73db+/n5MhU1j5qk5inYhJRbfGNPeTXyiDdLj5qmb8E5sglSYkKDjFY5ObArCvl6pynUdoEMYSJqZM4TQ8TEgcNiELgy5R4qDsEFx19xCYVP7hnbshgpKchEI2Tpm481wHPxcWlpSu93W4eGhut1uLI9PmYoxMGYHE7m+p4m7knBgljl470qvjPWOYZVKZcbtwRoEM/Lrw+QuuNN95Dkvjkk45sK1yM4sl8sxqxfs5rLoSgoKaXpWZql0dirVzs6OlpaW4nmRnm/AA5SmjWZSDTmZTGKloCdneT5FynAX+Z2YpjxUNgeNWgC6sG7Q7ICvpPNCFyU2eX9EtJ8XQ2GC0yYOoeDmuLe8d4EGw3i1pRdnOXbjvjRM5uSM3+12Z4Qw13crhP4ZCOdGo6Hl5WX1er2oDLrdbjync2lpKYZRnTmxJlzoeQIUYyMK4paFu6cICH78mEf6nniavYetfX3cynSXhHHgNhAGd4wDDId94K6kpJneJZfVgVu6woKCTVatVnVycqKDgwPdv39ftVotHu7Cw5Nm28A7VuBNahwrSB9mCmbxwFPt4Ildk8lkpsWaC5cUDOU+3vuCze/RBUkR/T89nR6Z6H02uH+aC8C9EFBECqTZaBDVpKkLxuZPTXeEFxaPz9Nb5YEXObnwTXENfP56va5Wq6WVlZVYtenuGXOCgVIrjtcx+Z3RUszAzX5AchgYXIkKYVwO1jbFq1hvXKNUwbgFhmuMkvHxMwZ3gebn57W/v692u60sy2K+yWXRlRQUWZZpb28v9pgI4ezYODL4YHQXEp6OnGo+mIMNwaZIGfpJY0mBK2nakRstnib3POnarvn5360JxoXFwKb00BmdpbEc+C4MAabh+ATfT3si8BsAOI2gMJ4UTGZdnYnconFgz9eQa/EMPLxIU2GYifCiA8A+Nv/flYPjPakbw2t+DXdF6TlCLxPeS++ZAqTujvGccIt9HTwng3tL00pb+nYUi0Vtb2+r1+tFIXNZvSikKyooJMXzOAk1cr7k3t5ebGTj5Ki4x/AJTXq8P/W10Viu+X1jpei2M4e7Mqm5e5Fbc1FBkQOqfJb7X2TKeugs1WBopYvyP8bj8cxBStyL73B8HeuJIHCsx99DU2JZuFVzkaBIsR0HTXHPpGmpu5vpniTnzwXLysflbkb6XKWpq0ouiruYvl9SV+NJGIE/Qwe63fV0a5cepICijq94f9KHDx/OHAFwmXRlBYU0CwKBVbTbbS0vL89EMQqF6UlRbGrqJ7wC0E+zdt/cN0CqgVIXgs2N1Jdmi8Qwf9kYLoic4TFXad2elpdzr8lkmqXK3CaTSQy9pviJXz9F3tnEHsL1CAAVjISdYcLhcDjT98MFCPf1k7YcUHYh4evIb7cCYXju7We0sh/4vAObzBsh7W5Flk3L8N31Yk0d53CrwoUi68kYHMhEaHvKPQIKgeefpZydRkZYzZJm6jk4JZ7vXiaQKV1xQZHGjYmrp2Etaaqxj46O1Gw2Zyo0K5WKqtVqPNfCzX3p8egETOZ5BumGIcZ+EdiZ9lbAx5emNQRoY3cVXDgwf0kxlu7gpK9RCqx56NPTh8FS0pCopCgkAJBZg8Fg8Jh7hOUEE5Nclroy3FOa5lF4Hguvu7AH8/AoB+vu2BHXBQh0jGg0Gml3d1enp6daW1t7LFnLXRqeWWoVOG7DPGFk1sH3BMIrBcoBTYvFs45f3hiH63Q6ndiwhut89tln6nQ6qlar8UiGy6QrLShSs5xmtN4ZClOOKtP5+XnV6/UoZLAoQLAdIWfxU6Z214LogZva/hmIz3ihloc2j4+Po+9LRMPNadwF97NPTk5myo8RQC4c2dSpNQG5VYHlwjEHfNZPL6dxrAs/kp8KhUK0pLw7t3cK8+8x9xTXwPrz9/i+CzzWw2ss0N6caDY3d3ZkH/1VMdXr9XrU3h4t8AhP6p5BCHPCvNSsgBeQ9MU1PLGO5CjHqxAuNC8ikuMnwXe73fhch8OhHj16FNfhi6ArLSg++eST2MJeOnvICAFPnyURCd/OGW5ubi6WCkuPVzp6cY9HCBx4c9MVLczmc9NSmmr44+NjdTqdKAAQaI1GY6YXhJcaA2ixORwQA/3nfszHMyJZI0/3dVMct8KjNh45kBQbB8EkWE58D/eNObL2rIszdak0PeELTYmVgnBKsRhfT8bkZj9rSTiae8DQuJmUicOkjiWxFmlBm1sQWEXkXNCty1sBAHh6yJVwOc+P/BJ6VHhmrUc8wIcODw91//59HR4eqlKpaG9vLwr1y6QrLSho1kF+AunFmLoIC7SPa4gQznIWaIbi4JqTa2NPOU4xhvTzbFjHFKRpC3Y2lvvAaEcHGx2M4/tgEZjU3BeLgms6eOfjwGJxgYlFgqb01x2XAUAGkEMIOPbAeAAfMa8dnc+yLGY2IgwACo+OjuLBwynQyLWfVCl5EaBMngGZrigCMmXTegzG7ACmh1TdPXJ8wAUkCiQ9pxUrB4uQQ5kcJCUlnKxULEdySR49ehSVkQPXl0lXWlBI0s7Ozkz15MnJiRYXF2e0JsLDTUwAqxT0c4ZyyyEFAfnbtZlv6hR45D1SzofDofr9fhwnNQJYDR4mZDMiqEgAQuvBbN44xnMpXNCMx+O4JvjJMIOnVKdMCnPMzc3N9D6A4RBavo6SYpIS1aDeBo58BLQmwhtrAwGSCroUm3Atz3ghbwbMcQzps3NhdVHo3DESPs977DGEHmMkVZ01p08FVtjJyUlsYiNJq6urMa0fVy6EEHGJYrGofr+v7e1tdTqdmO7+RQgJ6TUQFPfu3VO1WtX6+tkZRJ1OR/Pz8+p2uzNmuJvHkmYO03FTW5pujjSW7dEBByXTcKM0Tf12n3o8HkdfOcuymBZMcpE020fBzyDBhXKf2A+xReuQCOTdnjySgTuQMrTjK7gVHpng+w6sOjaDue1ALOCcd2BC20rTKAyfw41ZWVmJzyOtmXGLjXVOw8F+bAMCzDN53TJIi9kcyPUxIrBdCAMGgy15ejoCMnW3uNdgMIip13Q0Y0+i+FgvcnLa7XY81Z5+o18UPeshxf93Sf+OpGNJ35b072dZ1g5nRw/+S0nfOv/6r2VZ9qcvY+BOhIzm5+f1+eef6wd+4AfihsGUhgH48cpSR9aladang1Y8eBjDQ57nazKj8bwQis9TK5CCk2jEo6MjhRBiBMOzLb1OAY3E8YArKyvxAB7HPGAKT+pB07l5yzhwB5irh0ABhGmQw5gxo6UpUxHrB4fh9HQYHu2LNsbHB8W/detWfHbu6nk9hGNGnifhBVd8r1AoRIsNVxTmBUsoFovR0mGsLtB4jTEhWJgDz1hSfBbeVInaHdYHIbiysqKlpSWFEGIxI+ND4BWLRe3s7Ghvby8qgO1tjtL5YuhZDyn+FUk/k2XZaQjhP5P0Mzo7e1SSvp1l2Vdf5CCfRnfu3FGn09EHH3yg7e1tfetb39Lt27dnEOhisaiNjY1oenv5uVeU8ttNVn673woDYjlI0zCp5wl4uLTdbms4HEatwxgAtBwPcLDOy6y5D2MaDAaam5ubCe1yHQQV43JwtFQqRYuKVGmYlutCznTkdaDROarR1w0M5ujoSIeHh9G18sa4WEcAjxx9SIq0J84xL6wZ3CmYGqtwfn4+Wg7OyNIUw3HrAYHKmhDSTTEmyF0eL7hjfyE0ASJhfNYbHC3LMq2trWl5eXkGMGa+9NYkbNput/Xhhx/q4cOHMWJzmenaF9EzHVKcZdk/tn9/TdL/5QWP63umg4MD/eZv/qbef/99ffzxx6pUKlpaWooaxE11JLb79WhVd0H84GI+jwbySEBanepuCCZkr9eLZ1jQFFWaIt/0UkSIuM/N5kV4eDgTlwXBwbUqlYr6/f5jKcBuqXgiEmFQBNrh4WGcj6P0i4uLqlarsaSfNcGaQUDgTmRZFou4vM4DawfmpyCP0KYkHR4eRgASoeZNcrknTMhaEjFhrby03SNHhKVT7e1umGNTWBmsG/vEIzP8zxzZV26V0kMjzc+RNGOJ4KZ9/etf1/379yMgzLP5IulFYBR/QtLftf+/FEL4TUkdSf9RlmX/8wu4x1PJ8wd2d3fV6/Vm+lL6yeagywBHPDTMWR4gm8o1kPvlnhHIJnEN5+FUz4twV4AcDjQLn/NQXa/XizkeDvKVy+Xoz/N5r4MA7GSs7nv78YUp+DcYDLS3tydJUdsDnsFQ5HvAvJRle3tCsjkJR/pYYILT09NYui0palByCHzMKYPhRkJcwwuwsNxSoJL19v4RvpewGi9iZklx3m5x+D70cDBuMGCtR0kAPsFqqFglXLq5uamPPvoonqdbKBSiwvki6bkERQjh/6azE8H+q/OXNiXdzrJsL4Twr0j6b0MIP5BlWeeC737Pp5k/jcbjsTY3N9VoNCI45mYqGo5NAACFVvVsP8ceUiQ83UCOb7gvflH0I8uyaD3A+L55wEtgElwgNhO+vcfpuTdCBq3jwF8a4cAl8XwArBFvN8cYmfdwOFS1Wp0J4/b7fQ0Gg/ibdSfZCDAZC4T7TSaTx6I/WDysEQzngO9Fa+7rwvPhXghM5uzYg7uIKd7kGIk/V7dOfK+4UPOIB8ICy80FM4qFOWNNnJ6eHa/46aefan9/fyYi9jLomQVFCOEndAZy/sHsfOZZlo0kjc7//loI4duSvk/Sb6Tfz57hNPOn0WQy0c7Ojur1ujqdTmR+NEyn03msSUoainOmOB9b3CBYCWk0RJoNo7HJ+C7+MVqZccHkCA0ANndlsBwwW9lYHtrF30cwEHrzPg1ewIR5DyOdnp7OWAWg9w56Alz2+/3oz2PSA9J2Op2YQVgqnfWLwPLwFHVcH1xC8knclel2u1EgSYonyvMcsNo8vyG1qmBuwpL88KxhXD8rlGfOM/Xn4MITwZuGZf2Zp/vFcR/2CPvQXQ5cob29PX3++ecRq2DNXwY9k6AIIfwhnYGX//ssywb2+pqk/SzLxiGEd3V2SPGnL2Sk39v44glUPCQ2CEAcAJNn313kd6YAFowcQogdtfi+pxy7KwEz0y3JLYYQpg13YUiiBPjsuE1ses/ek6ZMxP1hYDQ5ggqBiHDyCAFpwmgssB2YolA4S5GHqf2EtuPjYx0eHurw8FC9Xk/FYlHNZjO2I+Te3gjGIyQe3iSy0m63NRqN1Gq14vyJAOHCuJtFIhM4jOM6+PvMA+uJOQLkMl/WNcUiPKrlFqK/x+cdAMdyoskM+4F6I+bmLf+73a4+//xz7e7uxtwg34tfND3rIcU/I6ki6VfOpS9h0D8g6T8NIZxKGkv601mW7V/S2J9Ic3NzarfbevDggebn57W0tBRPjl5dXVWz2Yy4APiFlwvz28N/YAWU/xLz9iQn36C+aciBIPRYrVZjdAGtiIYktRiNj+lO5ydvbMKmBVyF4TzZDNCRyMF4PI6aDWHomxlzmvAu1/Z50YgGq4miJUnxTFEOMSIsORgMYhai17ik2bLMCQun0WhEhnazm2fHdXz+4BZ8hmIy3M9OpxNPTW82mzo5OYkHCrnAkB4/z9VdG8acjp/7unXX6XS0ubmpbrcb+4O6AvFivpOTE+3s7Oijjz6aiXykLvAXSc96SPF/8YTP/j1Jf+95B/W89OGHH+rGjRu6f/++xuOxbty4oUajoVu3bmltbS12TULjSZoxa3lgxNfZWP1+P5rMaLn04B0eJiY/AkSaJnnhw3tPTFwRGJ3NjkWBtUEiE5+DQdCKWAto0IWFhRhRADh0EBbQjP8dKPXS9vn5eR0eHmoymWh7ezviI2hsMmSXlpZUrVZnBCkJRuAYw+Ewlk8vLi7GyBLWBwKSqA7PCAwCLAbCfeSH1P7hcKhOp6OFhQUNh8NYZclYyG9YW1vTeDyOwK0LL3dJYHwXcG5VeEiV7/FDC4S9vT0Vi8Xolnl+C89td3dXH3/8sfb29uIYXzZd+czMJxHNbKj/uH79+syZCzw8D4V6Ci5CAksCN8bfw/9Ey3rvAypY0TghhHgC2NHRkWq1mqTZIwE9NOrMgXVAQ1VM1V6vFz8Dc3oOAczHOHA1JMUqUZjOcRGKkDwxDUGBYONIREDGwWCgxcXF2IrQXQ7O2vTiJ8dAGo3GTNRjOBxqYWFBrVZr5lr48r6mvobSNEzMs8+yTIeHhzFUnD47xuVgpUdnnOERFFiILhSk2f4TrDPW3MLCgtbW1qIA9GtioYzHZ0cRfPLJJ7p7924UmLmguETa2trS9evXNR6Ptb+/r06no42NDUmKiLKbr5jMPEgkPRoQjAAzl+uwyQCa3Oz3fhmAV7gg9KqAWfBB03g/FgW+ONckjEn9gms6ruOp2w78uXBintwP4A/NmlbhEvLESkBA3bp1a+b8V7QzwBwRlSzLopmP0PUjBBjHwsKCFhcXZ7Jp3Q3wubgl6AAi4dqjo6MZ8FOaVteyjjwfr2wFIHamdqHqxWTSFJfy9G1pil/V6/UIolLbgtDE4vn000/14MGDKBwcAH2Z9NoKin6/HzVguVzWp59+qkqlMqPJ2cTORL4h0TaerOQhODYEoUIXFJQFpwBUsVhUrVabCaUhQNg4uBOeL4BpisVCLgKan9bynuLs+R6e1wETYEUB0GLtIIzSkCP/A8RNJhPV63VVq9WZ3BTcEqwD1pG6C6wE7gEOhMXj2Iy7FK6pIWdS3Ef+pwR8OByq0WjMhFoZr1tV/CB0U0ASIeFCwEFl1tqFiWNG7L92ux2tT55jp9PR3bt39fnnn+vg4CDiVbmg+AJoZ2cnMsn+/r6KxaLeeustNZvNmf4TPExPwXa0HA0G0/A9UHg0vmt2xys8Y48MSMcFnLFTK4IxeRMdxoMv7rF6hAMM4aE7aTaacRGCXyyeNbHFREdju7/NHObn52O3MNwMr6XB96b2Q5qeui4pmvg+RoQjgoZxSNMTtaTH+5SmUSquhaXnGZ24AzwzD6XyHQd//Qel4K6ar6Mn6UnTPArP0GQuhK2Pj4+1tbWlb3/729rZ2YkJVWky18uk11pQSNL+/r4ODw9VrVbV7Xb1la98RTdv3oy1EdQWeKmwNNVIPHg2kucedLvdiOy7+evXkB4Pq/nxgGxyL1LCEknBOQA68hzcLcCS4Lv89vCd51ygfRkbhHXg+Q9gMFhMYB/z8/Px88wJ1wOLDReF6BJrgxCQpjUynuTE3LgWwobkKc99wNrzyBWvAWSSkctauPvFs/DvMh4EiQsdd+l83B4W5btYa6nbwrOUziqeP/74Y21tbUWBGkJQu91+cYzwnPTaCwpJMSw2Pz+vb3zjG9ra2tKNGzd0/fp1Xbt2LTa+cWbFl4RBMLvZANQycH33591sTcN//pvvYg3A9NwbZJ4NiGYkTwGhQ+s0STOWigNqaR2JA3l83q0MrBqYi9PLYDrXjjQA4v4IB8roDw8PH4vMYIUg1ABhEXBe2YplQkIY4WLG5vgBf/NdBEyr1ZqJ4nDfYrE4swYXJU+liXBu0XlCngsWf27uwvC5wWCgR48exYjR9va2RqNRjIp5PdCrQG+EoJDOHu7Ozo7effdd7e/vR2uBUGmpVIqJToBrvlnc5fAKTHotAOTBkJ64w2ZBCHnDGdeivjFJbQYH4T7ZeSYfTEyVogOnvpEdXHULxLW2nwvKeBxwJGvw8PAwRoq81Z7ncTBOrgfGg6CFgUOYnhrmJr2kmQYuCAka33A9F2ie6+KuGIKSKs7UenDz3zW+p6e7Bcn4HBdhjR1AZv0d+6Da00HUXq+n3d3d2BKRPfiysi+/E70xgkI6Y6B79+7p+vXrGgwG2tzcjAg9vTdhZjSvZwQiOEDra7XaTDQERkAIpHn5CAtp9tg7NJsLDTYnJipMTMQE5iRSwwZE0MBchBMJ0dEz0vEHfhAo5XI55l64yc3RfouLi2o0Gmo2m6pWq7HCU5o14SVF4BPBQWapp2sjEGBawEzvfkVnsG63O5Pr4BEaLBWS4mBihB7JcQCIzIvnKk17WiBIuUcKoqbPL50/3/EkuWJx2qXq7t27unfvnvr9fhSah4eHcX+8avRGCQpJMYcev7fVaunu3btaW1vTyspK1LiOdEMepUAgkAiFuewbMgXYuJ7XiDiQ6ZoI7TganR1vTzUs93V8xUFLNiY5C51OR1mWxQxJB/O4P5qZnIg0U9WjIoXCWSo3pj+RFrQzQtRT1InKMB/u7cAnbgTWAVaZ14ZQrs89HIQmksI6MDfP53BQmTWivwbr75aFA8p+XQd23V1xckGPoBiNRnr06JHu3LmjR48excJFDw1fdK1Xgd44QSEpJiplWaZHjx4phBBrMNhQbkp7rNxDpG55wBCe1p0+dDdjpdmEoTTfwV0cz/dwHxuBIU0P1vFoDaenYZYzHq/74JppsplHEqRphmgIQa1WayYhyxORnDlhMIrhGBum9XA4jJ8HaEXoeVTK60HQ+ICKbtG56+XCL40e8By4rltY6ecciOR7rEuKObmFkboxo9FIOzs7+vjjj2MIFBcKa+tVpjdSUEDdblefffaZisViPG1bOtMGuCH0cYQ50cZYFZidMC6+MMKD31z3otBaqVSKCUaOpHNPDjXycJ3f11OyuQ+ZmxxiS76IND04yUOTLgy5Hq/hisA09Xo9ukCehwGm4W6YpJk6F/AXj7aACSGksdAQ1H79ZrMZXQe0cQhnhV0OIrK2kiLewHg8AsW9EBQuMJzh08iXC3fI8yJ87MPhUPv7+/r444/1rW99SwcHB9GVLJfLX3i3qmehN1pQSGcNXL75zW/qy1/+sj766CO1221du3ZNzWYz+uzkFXgoEmHhGyyNejg4ibYheuBFRYBYjUYjhmj9xDMiGoB53AvtjPZMEfpOpxMtJdB9IhdpZy9pttu2a3DuR+QA90dSrD2hlZ+nZpNwRUNfN8cRBIC1WCjStGEQTEZ6NwIH7AXgkrXAfXHwFmEAuZVEpIW/UwGXWgUeuZBmw9uM1zN0Sfja2dnRhx9+qE8//VTdbjdex6NOrzq98YJCOts8n3zyiW7fvq1Op6P79++rVqtpeXlZ6+vrEQ+AufA3MdPJdKTuA2GBZeGgFinF3rClXC5reXk5AqukdHuIs1KpxB4bMDQ5Ec4QDj4Oh8MYFvZcD2cwN5ulaRpyqjE5YMjDvf49adphis94mNW1/WQyPUYPpgaQRdt7pGMymcwArMViMXbFgjkRLrhHhJOpjWFcnljmOE2ayMXc3HrBokgzSV2QEG3Cirh3756+/e1v6/79+9E9I0v1KlEuKIzu37+vEIKWl5djQlW73Va329VgMND6+no8JSuNcsCoMDGmuWtrj7dPJpPYbLdWq6nRaERk35Fv3A2EDEyXZloivADovMkJWII07cvh2IpjFf46TItbBXM4uYtFkZgLGo8kAfx6Ex4EBmvoQC/ulmdnOh5EBS+vM3YENe4I53hiQTn4jFZ3oJp5pxmjKIJisfhYEhXuCdWfDx8+1IMHD7S5uam9vT2NRiMtLi6q3W4/Fg27CpQLCiO0d7vdVq/Xi0Lg0aNHGg6H2tzcjL0sqtWqWq1WDKs6fuBuAEwKA7kvTriS08rYvNIsMOamsudfeMQELUb2IgJraWkpJmR53QTf4/qerQjjO2N4puWTojYwi+ePwPwIGn4QFuPxOApA3CzPt+B7ngzmgDACycftxwQuLi5Kmk2K4zu9Xm8mxTsVAsyBvUG+BG4DQhcXqtfrxa5Um5ubarfb6vf7yrKzQrhOp3MlhYSUC4oLCVO2UChoa2tLxWJRN27cUKfTidYChw6tr6/HLkzus4JTjMezxwM4poHF4ZWaru3BCSAEizdZAahESNDzcjweq1qtRu2Lqc0PmtUFhTTbmMUjKW5tuCXiEQbXsAhDLC9Kq1mPNLMR894zMt1Sk6baXpq27GcMnmDm+S6FwlldCpaHC5ZyuTyDl2AROd7kVoXnvLhLSd+Lra0t3b17NyoWz9HguVxVygXFd6DJ5KwHp6RYFo6GLJVK2tnZ0f7+vlZWVmYsDWmaDwCw6JoEQeEbn42INeChP8cJSHAi14F0ckBSEPfJZBJ7YnBPsABPnXZQlR+Eh1sVjMkjBR65YCx8xl0FBAlCAvcA1+iiaIUDkozHw5ee+k5iGKeuebQJvAdBB7H+/MDICHIiKS6oWS8EEBGcdrutra0t3bt3T/fu3YvYiOMXVyGy8Z0oFxTfJd25c0eS1Gg0tLS0FNviHRwcaHFxUaurq7F2pNFoaDweP3aqFkzj2skxCPz1J2X8eZEUDEoYFI2MoIBJ/FQuzxqVpp2lpWmnJsxrmM8FD+PElSEcipZG0HEtQEKPCvn/jJ+xMA5ncl8HQrFYIAgaLBTW13Mk0ihMmgeRRqd8PJ49igDy9e33+9rd3dXm5qYePnyo7e3tmAk7mUxmzha96pQLiu+ROp2OOp2OCoWCbt++rcFgoE6no3a7HU+brtfrunbtmpaXl2O6dXr4LCar51VQE+DWhqdP1+v1mZwJrpH2fCgUCjFKww+MArl7kwKuae4AAoMO2/V6PYZyaXnHmBlfGpKEHE84OTlRp9NRo9GI7yFcHCdwyyXFJlxgSopZoghPyuGXlpZmckakKRBMbw2sCq5JdzS36ojGDIdD7ezs6P79+9rd3Y2VoFiar2qG5bNSLiiekSaTSbQyJGllZUWTyVnfi1KppNXVVb399ttqNpuq1+vxx6saISpFEShsVscN0KZoUfz51D8nwxEQFUHlJrUnbPFdTx13nALh0e/3Y51HapmQ9AUTu9BIk80IM8LMHiHCWsKNAAgkbEw3LwdIwSo8+uGHIVN8xTWk2QhUrVZTp9OJOJH3+MCtwVKiEvXevXt68OCBdnZ2Yh9Ralru3bt3ibvu5VEuKF4QtdvteG7I22+/rf39/ZhIVKvV4g9p3p6qXCqVYnUm6c/OiEQwvELUczkIq0qKlaTcJ/1O6nMfHR3NAKBuYWBJjEajCORyTQdtuRZ4iJv4DlY6XsPYyfJ01wALiuQkEqM8agR4TMSBrtrMCWHrc/AOYJ5ePj8/P1Os5kAybgSVng8ePNDnn3+ura2t2BsEDMKzNF83etbTzP+CpD8paef8Y38+y7J/eP7ez0j6SZ216/8Psyz7R5cw7leO0Gzj8VgPHjyIzM2ZF/jJ1IIA2JEPICk20yGlGSaEien5KU3Tk7FOPNIgaSbygs8vTX1vtKSHBxk/7zvmMZlMopBIO1M5TpHmIGC5pHUrJycnM8cVeoMc8B/6kTImKnbTCBHZoYCRHunxUnbGjMBC+HgNB2uHK4G1t7+/r+3tbW1ubkZXAzfJ0+NfV3rW08wl6a9mWfaX/IUQwu+U9GOSfkDSDUn/fQjh+7Isu5rB42egyWQSD5Hd3t6ewQtgNkA9STP9KzY2NmYayrrLgKb0zthoVMKiCAZnIM+RcHcmjWg4LgB+4K3n6PdItMWtDgcoU60K87lrw3jT1HjyK8AASD+XFDtUIQCfFC1K0+s9J8PdH++DiqvhmAyuDy7X/v5+xKFIV/fuVa87PdNp5t+BflTSL2RnRwt+FkL4RNIPSfrVZx/i1aWHDx9KkpaXl2MkhI7Q+MJYHuVyWZ1OJ7a7JwRKN6f5+fnYdJU6FKICaH6vp8AlkabAH0zp6cde25AmFTniTxcrwFgvx4fp3VXyXAbPfSCUyxmmkCdL4UoNh8MoRL2GBeHmVgDNgHCPPLEKwBIh5ZW0CBkEJZ9pt9va2dnR3t6eOp3OTI4Kc6MD2ZtAz4NR/FQI4Y/r7FzRP5tl2YGkm5J+zT5z//y1xyhcwiHFryrt7+9rf/+7PzDtK1/5ykwb/izLotuyurqqmzdv6saNG1pdXVWr1dLy8rL6/X4EAvkOzOSmeuoKuOb18m0+i4AhJyOt4yDPgbCgWy4uaMrlctTEjBMmlabChTn7iWKebk0rPMijFZ7nkha6TSaTmeMMWQO6huFudLtd7e/v6+7du9rb24sdp7x69ODg4Pk3xRWjZxUUf03SX5SUnf/+y5L+hKSL0JwL40TZJRxS/LrQ17/+9Zn/Qwj6wR/8QWVZpvv37+vOnTuxjuL27dv6/u//fq2srESG9w7WaD3PGfBEK4jQLMzu7etoost3cT/IFQBjkabhS7cmHIvBIpAU70FjFw+PejGWpBnBiaBxd8GzWMFrPHrT7XYjjuGNcKhO7ff72tvbi2nX7mIwH+o43kR6JkGRZdkWf4cQ/rqkf3D+731Jb9lHb0l6+Myjy0nSmVb82te+Fv+/fv26VldXtb29rcFgoL29Pd26dUutVkvNZjOmlAOmUmtA70i0M0xEVIBwpacd813Sm8FP+K43rvGDiyVF10SSarVadAuy7Kzj1mAwmKm2pJozyzLV6/WYdYrAKBQKMdKQRlH4PhWkYDz0AwXExALBtblz5462t7djMR6WFAIWvINcmTeVnvU08+tZlm2e//tHJKECf1nS3w4h/BWdgZnvS/rnzz3KnGYsgO3t7di5an5+PpbGHxwczPTT9JoVMkqXlpbiwb8wnpdpb2xs6OjoSO12OzKauzMeJXAXxYFCPgdQiMlOtEeapnxjceBmELHA3eL7WCoAj/TF9PDtycmJ1tbW4nUolsNy6fV6Ojg4UL/fj1GQBw8eaH9/f6bs3etBJEUX5HVLovpe6FlPM//hEMJXdeZW3JH0pyQpy7JvhBB+UdI3JZ1K+jNvUsTjiyI0bAhB29vb0c1YX1+P1ZKkGANOekPcZrMZ6x/8BPO5uTnV63UdHByo1+upVqtFjAGcAr9emlZkSprpzeAChUxTQpRe+SqduR/lcjmO0/GISqUyU+FJly0HNgeDgbrdbmR+rJ/BYBBdMCIYe3t72t/fj8BnpVKJJeAksHkK+Xg8jgdUX9WqzxdFL/Q08/PP/6ykn32eQeX03VGWZTMg6enp2dkZhP3m5uaiyd7pdLS7u6tKpaJqtapqtTrTY4FOXp4xmYb/sEKc0T186xWV0pkg6fV6MwVSRE08WoHFkpaM49bQxYuu54CnuEtYE9JZ4hsFaLhQhDUPDw9jc14sL/IhuK/P11sJvumUZ2a+RrS9vR3/LhQK+uCDD2bwBj+3pFKpRG1PRKVcLuvBgwfRsqjX6zGXoN/vq9frqdVqzSRwITw85Iu5XygU1O/342HK3NcLyDhlXJo2pQHroHy+2+3OuCbkW9AlrNfrRdyC+gxCrEdHR9rf348uDoITIejRGzARQMw3MbrxJAqvgrTMox5fDM3NzUWsolwua3d3d6bwihDm4uKijo+PIwBKD4y1tTV9+ctfjjkTWBckfxFNIJy4uLgYu26VSqWYD0LyEx3EaN4D0AkO8ejRI+3t7cUMTXAOXALwC0+jxjXx0CwWCeN0/ITXCZ0i2N5Q+lqWZb/nojdyQZGTpLMakRs3bujk5ERLS0va29uLzAlA6cVmRFVCCGo0GjFz0lOjSRjDFfJSeVwOekiAJ+Ay0HPSO1qTfbqwsBCFh2eaSrPl+IRecSu8elaaukbSmcty1fpYXgLlgiKn70yY5mjiyWQSgVGsi2KxqGq1GlOlYWzCiOAHXI/fMDRMTKblRVmN5FPgdvj5Hqenp9EtQtDgNgwGg5kDhcBReJ/6kXa7PVO5yxjelAzLp9ATBUWOUeQkaXrYjTQNxfrJ2jAvqeEwpadwe60I10P4UG/heICkmf6bng0K4MhpWp4WTnalh2E9auJCyhPIsExyofC9Uy4ocnoieZs8muKAQ3g/DM/GxHWQHj9RHRCSPhO8DiEMvAHORb0mr3LvyatKuaDI6amUZdljAJ+nV0vT5jouKDxM6aXrr0t7uDeJckGR0wuhi0688vM0r8qJWDldTIWnfySnnHJ60ykXFDnllNNTKRcUOeWU01MpFxQ55ZTTUykXFDnllNNTKRcUOeWU01MpFxQ55ZTTUykXFDnllNNTKRcUOeWU01MpFxQ55ZTTUykXFDnllNNTKRcUOeWU01PpqYIihPA3QgjbIYSv22t/N4TwL85/7oQQ/sX56++EEIb23v/nEseeU045fUH0TIcUZ1n2f+XvEMJflnRon/92lmVffUHjyymnnF4Beq5DisNZK6E/Julff8HjyimnnF4hel6M4l+VtJVl2cf22pdCCL8ZQvgfQwj/6nNeP6eccnoF6Hkb1/y4pL9j/29Kup1l2V4I4V+R9N+GEH4gy7JO+sU36TTznHK66vTMFkUIoSTp35X0d3kty7JRlmV7539/TdK3JX3fRd/Psuznsyz7PU/q+ptTTjm9OvQ8rsf/QdKHWZbd54UQwloIoXj+97s6O6T40+cbYk455fSy6bsJj/4dSb8q6ftDCPdDCD95/taPadbtkKQ/IOl/DSH8lqT/WtKfzrJsXznllNOVpvwAoJxyygl64gFAeWZmTjnl9FTKBUVOOeX0VMoFRU455fRUygVFTjnl9FTKBUVOOeX0VMoFRU455fRUelXOHt2V1D///SbQqt6Mub4p85Rej7m+/aQ3Xok8CkkKIfzGm5LO/abM9U2Zp/T6zzV3PXLKKaenUi4ocsopp6fSqyQofv5lD+ALpDdlrm/KPKXXfK6vDEaRU045vbr0KlkUOeWU0ytKL11QhBD+UAjhWyGET0IIP/2yx/Oi6bxL+W+fdyX/jfPXlkMIvxJC+Pj899LLHuez0BM6tD9xbiGEnzl/zt8KIfwfX86on42eMNe/EEJ4YF3n/y1778rO9SJ6qYLivMnN/0vSvynpd0r68RDC73yZY7ok+teyLPuqhc9+WtI/ybLsfUn/5Pz/q0h/U9IfSl67cG7nz/XHJP3A+Xf+3zQ5uiL0N/X4XCXpr54/269mWfYPpddiro/Ry7YofkjSJ1mWfZpl2bGkX5D0oy95TF8E/aikv3X+99+S9H9+eUN5dsqy7H+SlDYmetLcflTSL5y3S/xM0ic6e/5Xgp4w1yfRlZ7rRfSyBcVNSffs//vnr71OlEn6xyGEr503FJakjSzLNiXp/Pf6Sxvdi6cnze11fdY/FUL4X89dE9ys126uL1tQhAtee93CML8vy7L/rc7cqz8TQvgDL3tAL4lex2f91yR9WdJXddaB/i+fv/7azfVlC4r7kt6y/29JeviSxnIplGXZw/Pf25J+SWcm6FYI4boknf/efnkjfOH0pLm9ds86y7KtLMvGWZZNJP11Td2L126uL1tQ/Lqk90MIXwohlHUGAP3ySx7TC6MQQjWEUOdvST8i6es6m+NPnH/sJyT9/ZczwkuhJ83tlyX9WAihEkL4ks46tP/zlzC+F0YIxHP6Izp7ttJrONeXWj2aZdlpCOGnJP0jSUVJfyPLsm+8zDG9YNqQ9EtnJy+qJOlvZ1n234UQfl3SL553NL8r6Y++xDE+M513aP9hSashhPuS/hNJP6cL5pZl2TdCCL8o6ZuSTiX9mSzLxi9l4M9AT5jrD4cQvqozt+KOpD8lXf25XkR5ZmZOOeX0VHrZrkdOOeV0BSgXFDnllNNTKRcUOeWU01MpFxQ55ZTTUykXFDnllNNTKRcUOeWU01MpFxQ55ZTTUykXFDnllNNT6f8P9Cy77R8Ci1cAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import nibabel as nib\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "img = nib.load('T1.nii')\n", "plt.imshow(img.dataobj[:,80,:], cmap = 'gray')\n", "\n", "np_array = np.asarray(img.dataobj[:,80,:])\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#### Exercise 2. Spectral clustering\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2. Image compression \n", "\n", "\n", "#### 2.1. plotting the eigenvalues\n", "\n", "Consider the simple image given below. We would like to compress this image by means of principal component factorization. For this purpose, we will see the columns as our feature vectors. Build the covariance matrix associated to the columns of the image (don't forget to center them first) and plot the eigenvalues associated to this covariance matrix. What do you notice? " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXQElEQVR4nO3de5RcZZnv8e9T1dWV7k46107nRkIGIySBECAmZGTOQlTIMGsGWLOciZeZuOQQBFmjs/RI0HFGjgdBFERmBMEJTkSFUZHLwRwlIJelIhAwEELuIYTcutO5kJBLX5/zR+3WNnTSla6q3rXr/X3WqlW73t5V9by90r/svWvXfszdEZFwpeIuQETipRAQCZxCQCRwCgGRwCkERAKnEBAJXMlCwMzmmdlaM9tgZotK9T4iUhgrxXkCZpYG1gEfBLYCLwAfdvfXiv5mIlKQUm0JzAY2uPsmd28D7gcuKdF7iUgBqkr0uuOBN3s83grMOdbK1Zb1QdSVqBQRATjA3hZ3bzh6vFQhYL2M/cl+h5ktBBYCDKKWOfb+EpUiIgCP+0/f6G28VLsDW4GTejyeAGzvuYK73+3us9x9VoZsicoQkb6UKgReAKaY2WQzqwbmA4+U6L1EpAAl2R1w9w4zuwb4JZAG7nH3VaV4LxEpTKmOCeDuS4GlpXp9ESkOnTEoEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAKuryYmW0GDgCdQIe7zzKzEcB/AycDm4G/c/e9hZUpIqVSjC2B97n7THefFT1eBDzh7lOAJ6LHIlKmSnGh0UuA86PlJcBTwLUleJ8TlqqtJdX4jgYsleFIKx07dsZdhSRQoSHgwGNm5sBd7n430OjuOwDcfYeZjS60yGJJNTaw84Pj4i6jJOqaO6l5SCEgJ67QEHivu2+P/tCXmdmafJ94dBsyEYlHQccE3H17dN8MPEiuG3GTmY0FiO6bj/FctSETKQP9DgEzqzOzId3LwIXAq+TajS2IVlsAPFxokSJSOoXsDjQCD5pZ9+v8yN1/YWYvAD82s8uBLcCHCi9TREql3yHg7puAM3sZ3w2oz7hIQuiMQZHAKQREAqcQEAmcQkAkcAoBkcApBEQCpxAQCZxCQCRwCgGRwCkERAKnEBAJnEJAJHAKAZHAKQREAqcQEAmcQkAkcAoBkcApBEQCpxAQCVyfIWBm95hZs5m92mNshJktM7P10f3wHj+7zsw2mNlaM7uoVIWLSHHksyXwX8C8o8Z67TdoZtOA+cD06Dl3mFm6aNWKSNH1GQLu/gyw56jhS8j1GSS6v7TH+P3u3ururwMbyDUkEZEy1d9jAn/SbxDo7jc4Hnizx3pbo7F3MLOFZrbczJa309rPMkSkUMU+MGi9jHlvK6oNmUh56G8IHKvf4FbgpB7rTQC29788ESm1/obAsfoNPgLMN7OsmU0GpgDPF1aiiJRSn23IzOw+4HxglJltBf4NuIle+g26+yoz+zHwGtABfMrdO0tUu4gUQZ8h4O4fPsaPeu036O43ADcUUpSIDBydMSgSOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4PrbhuzLZrbNzFZEt4t7/ExtyEQSpL9tyAC+6e4zo9tSUBsykSTK50Kjz5jZyXm+3h/akAGvm1l3G7Jn+19iER1ppa65Mi9+nN3dHncJklB9hsBxXGNm/wgsBz7r7nvJtRz7XY91jtuGDFgIMIjaAsrIX8eOndQ8tHNA3kskKfp7YPBO4BRgJrADuCUaVxsykYTpVwi4e5O7d7p7F/Bd/th5WG3IRBKmXyHQ3YcwchnQ/cmB2pCJJEx/25Cdb2YzyW3qbwauBLUhE0kic+91l31A1dsIn2O9djUTkSJ53H/6orvPOnpcZwyKBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgErh82pCdZGZPmtlqM1tlZp+OxkeY2TIzWx/dD+/xHLUiE0mIfLYEOsg1F5kKnAt8Kmo3tgh4wt2nAE9Ej9WKTCRh+gwBd9/h7i9FyweA1eS6Cl0CLIlWWwJcGi3/oRWZu78OdLciE5EydELHBKKehGcBzwGN7r4DckEBjI5WGw+82eNpvbYiM7OFZrbczJa309qP0kWkGPIOATMbDDwAfMbd9x9v1V7G3nFdc7UhEykPeYWAmWXIBcAP3f1n0XBTdyei6L45GlcrMpEEyefTAQMWA6vd/dYeP3oEWBAtLwAe7jGuVmQiCZFPa/L3Av8ArDSzFdHYF4CbgB+b2eXAFuBDoFZkIknTZwi4+6/pfT8foNfeYe5+A3BDAXWJyADRGYMigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCEjlMGPzDXNJDx/e97ryBwoBqSjXXLqUzd8dT7q+Pu5SEkMhIBVn5dzvs/meiaSGDIm7lERQCEjFSVuKl+cu4Y3vTVIQ5EEhIBUpY2kFQZ4UAlKxMpbmN3Pu5uAFU+Mupazlc1ER6UPVmEa6Gkfkta61d9K5ej34Oy67KCUwPF3Lv926mK90fILsz1+Iu5yypBDoh6rx42if2MCWeXUAtI5rp3H83ryee7gtw5GV52JdRt02Z8xj2+hqbqHr0KFSlhy099d0kr59MV/mcgVBLxQCeUrV1tJ63jRazqimdfbb1NcdZljq4Am/Tk11OzXn7AKga5ax868HcWDlmdRtNcb+ZD2dLS3aSiiB82u6+NLt9/CVrk+Q/X8Kgp76DAEzOwn4PjAG6ALudvdvmdmXgSuAXdGqX3D3pdFzrgMuBzqBf3L3X5ag9tIzo2rCeJo/cBJ7znDqT9lHbaaD2iK9fMqcVNoZMXMXXWcam84fQ+2v3sXQ19vJ/uoVvL2tSO8kkNsiaP/373Pz1f9A5rHlcZdTNvLZEuhuQ/aSmQ0BXjSzZdHPvunu3+i58lFtyMYBj5vZuxN3sVEz9i44l33zDjKifhcNJX67lDnDBh+GvznMnrYMBy4+m6k3b6Fjm67WXkzzalvhjnu5+ZMfI/P4i3GXUxYKaUN2LIluQ2aZajrfdzbrvzWb9sv2MqJ+4PfVa6rbGT2lhTU3NbL7irlUnTxxwGuoZPNqW9n1qcOkaou1TZdshbQhA7jGzF4xs3t6dCXOqw1ZObJslp1XzeKtz77N6HftJpvpiLWehpEHqLpsF699aTTpd58Say2V5qXZ97Lm9mmk6uriLiV2hbQhuxM4BZgJ7ABu6V61l6e/40hXufUitLOms/F/n031RbtIp7riLudPNI7dx5ovDmX3FXOxrFq2FUPG0rx+8X+y5vapwQdBv9uQuXuTu3e6exfwXf64yZ9XG7Jy6kWYmjmNDZ+rZuQZu0hZeR6ZH92wn9SlLexceA6WqY67nIqxYd7duSAYNCjuUmLT7zZk3X0II5cBr0bLiWpDlppxGus+l2XU8ANxl9KnlDmZi1rY+clZkErHXU5FSFuKNfPuZM1/nB5suOazJdDdhuwCM1sR3S4GbjazlWb2CvA+4J8h14YM6G5D9gvKuA1ZasZprPt8DQ0jyj8AuqVTXaQvbKH56jkKgiLJWoZV8+5g7bfPDDIICmlDtvQ4zyn7NmSpmdNY97lsogKgWybdRftFLTTZHMZ8Z7nOJyiC2lQ1j154O/88+yrsNyviLmdABfkFIjtnOuv/VzIDoFsm3ZXbNbhqlg4WFsn06ho+sngpXefNjLuUARVcCFg2y8a/r0/EMYC+pFNdZC5sgWnviruUivHx+mYuX/wwXX9xVtylDJigQsAy1ey88hyGT2+Ju5SiSae6WPtPg3QeQRHNH7KXj//nI/h7Z8ZdyoAIKgQ6zjudzIUtZfsxYH81Nr7FuitGgx2rg7ycqI8O2c1HFi+Fc2fEXUrJhRMCZmy6LFN2JwIVS/2pe/C5lf8PdiB9vL6Zv/3e49is0+MupaTCCIHoy0BDJ70VdyUlU13VyYaPDqJq/Li4S6koC4duZ+0VNVhV5X7rPogQqJo4gX3zDsb+XYBSGz15N03zJmm3oMhW/NXtbLjxPRX7KUwQIdB8wYRYvg0Yh30XHMGqwzvhpZSGpmpY9ZHb2fDVsyoyCCo+BFK1teyZUVkHAo+npraVjrnT4y6j4mQtw+r5367IIKj4EGg9bxr1f7Yv7jIGTF22jS0XZYP+QkypZCzNqvn/zob/c3ZFnbJd8SHQckZ1xR8LOFr1qfuxOl0woxSyluHFD9/KphsTc52cPlV0CFSNH0fr7LfjLmPA1WbbaL7s1LjLqFhDUzXc9rffIz3t3XGXUhQVHQLtExuorzsSdxkDLmXOwQlWcfuu5eSvao/wnvtWkTr9tLhLKVhFh8CWeXUVe3JQX2rPaSE9tjHuMira9Q2reO+PVpCakewgqOgQECm1fxm1hvN/+GKig6BiQ6BqTCOt49rjLiM26ZTz1qyxfa8oBbt25Hr+4ge/T+yuQcWGQFfjiLxbg1WilDktZ1TOx1jl7guj1jL7hytJn5q8r3VXbAiIDLTrG1axYUGp29QUXz4XGh1kZs+b2ctmtsrMro/GR5jZMjNbH90P7/Gc68xsg5mtNbOLSjkBkXLy9Me+ztbr/jxRJxPlsyXQClzg7meS6zEwz8zOBRYBT7j7FOCJ6PHRbcjmAXeYWXJ+IyIFGFs1mKev/jpbF81JzDcP82lD5u7efcZNJro5uXZjS6LxJcCl0XKi25CJFGpUuo7fXn0LW66dnYggyLf5SNrMVgDNwDJ3fw5odPcdkOtXCIyOVk9sGzKRYhmaquHZq27hzc/PLvtdg7xCIOo0NJNcN6HZZna8S60ksg2ZSLENTdXw66u+wdZr58RdynGd0KcD7r4PeIrcvn5Tdxei6L45Wi1xbchESmV4upYrPra0rDtL5/PpQIOZDYuWa4APAGvItRtbEK22AHg4Wk5UGzKRUvvM8M2MvG9f2QZBPlsCY4Eno3ZjL5A7JvAocBPwQTNbD3wwepyoNmQiA+X7k55hzP17qJo8Ke5S3iGfNmSvAO/oxODuu4H3H+M5sbchs/ZODrdlqKkO+NTh8L5AWdYWT/w1V/73XLbOP5mOTZvjLucPKvaMwc7V6zny6rC4y4hNe2eKST/ZEXcZcpS7JjzLhPubymqLoGJDAHesM+yr7lqH9sLK0V0TnmXMfbupOmlC3KUAlRwCQN02p8vDDII9LUPwg4fjLkOOYfHEX/Pmh8rjQGFFh8CYx7bR2RVmCAxZmaVz1664y5DjePDTN9N89Z/HXUZlh0BXcwsHVo6Mu4wBd6gtw7inK7fbUqU4JTOYn117M7uumhtrw5jKDoFDh6jbasHtEhw6lCW1cWvcZUgeJmcG89B1X2fXJ8+NLQgqOgQAxv5kPfsPhnMN/i43RjxWQ+c+bQkkxcSqwTwSYxBUfAh0trRQ+6vBcZcxYPa+Vceoh16Luww5QROqclsELQvPHfD3rvgQwJ2hmzs43JaJu5KBsbWGrrcPxl2F9MPEqsHM/Z8vkW4c3ffKRVT5IQBkn3iZA28MjbuMkjvUluHknx/BO8LquFRJ/mP8cxz+QQ1VYwbucvFBhIC3tzH15i3s2j0k7lJKquaBYaSe/n3cZUiBnpz+MIfuHTRgWwRBhABAx7btDH+mcg8QNu0cxqjHX4+7DCmSJ6c/TPuPsgMSBMGEAEDjsm007RwWdxlF196ZouHpDB07dsZdihTRL04bmCAIKgQ6Nm9h6k17ad5VH3cpRdPlRurBkQy/V5dsqDRpS/0xCEaV7qS3oEIAoHPdRob/JlsxJxDtaq6n4ecboUtfFqpEaUvx6Gk/Y/fFpeuAHFwIADQseYm2XzbQ2ZXs6TftHMbUG/fR2dTc98qSWFnLcNf1t7H/I6U5hyDZfwX95K2tjLlzOe2/HEV7ZzJ/BU3bhzH1xj10rtsYdykyAGZms9x5w7d466PFD4Jk/gUUgbe30fjt5+hKYBA0bR/GtK+20Ll+U9ylyACamc1y9w23sf/DxQ2CQtqQfdnMtpnZiuh2cY/nJKMNWVcno+94js7HRiVm16BpxzCmfm1PWV2eSgbOjOpB3HXjbRyYX7wgKKQNGcA33X1mdFsKCWxD1tXJmO8sp73Mg6DLjaamoUz92l5tAQRuRvUgvnPTbbz9d8UJgkLakB1L4tqQeXsbY+56kZFfytDcXH6nF7d3pvAHRjF10Zs6BiBALghqr9xGur7wj7vzapQW/U/+IvAu4Nvu/pyZ/SVwjZn9I7Ac+Ky77yXXcux3PZ6eiDZk3toKv1/FaV89hXVXjKb+1D1UV8X/sVvTzmE0PJ1h+L3P06mPAY/Pnf97xfk8WFv+/f+KwRy8Y3fBr5PXbyvqGzAzakLyYNSG7E7gK+S2Cr4C3AJ8ghNoQwYsBBhEbX9qL4nOdRs55fOb8Lkz2PDRQYyeXPgvuT8OtWWoeWAY0x5/XWcCngD77csE8n1RALqK8BonFJnuvs/MngLmufs3usfN7LvAo9HDvNuQAXcD1NuI4+1eDDx37LcvM/WNcTTNm8S+C45QU9tKXbatpG/b5cbet+pgaw0n//wIqaefRd8HlFLrMwTMrAFojwKguw3Z18xsbHdXYuAy4NVo+RHgR2Z2KzCOBLch69i2nZH37GDUD6rpmDudLRcNpfrU/dRm20hZ8XLrUFuGQ4eyjHishikPvUbX2wf1dWAZMPlsCYwFlkTHBVLAj939UTO718xmktvU3wxcCbk2ZGbW3Yasg6S3IXPHW1tJP/USp/xuEDa4jubL3s3+CUbd2S0ApFN+QqHQfV7CnpYhDFmZZdzTb5Ha9Aade/eS3F+UJJW5x78lXm8jfI712tGsbFk2S3ps7sIP+94zlt2n5/cpaPoITPrpTqy9Az98RKf8yoB53H/6orvPOno8jMOoJeCtrXRs3gLA4M1bGPyT/J+r/+2lnJTv2TEiMiAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgyuKS42a2CzgItMRdSwmMQvNKmkqd2yR3bzh6sCxCAMDMlvd2TfSk07ySp5Ln1hvtDogETiEgErhyCoG74y6gRDSv5Knkub1D2RwTEJF4lNOWgIjEIPYQMLN5ZrbWzDaY2aK46zlRZnaPmTWb2as9xkaY2TIzWx/dD+/xs+uiua41s4viqbpvZnaSmT1pZqvNbJWZfToaT/TczGyQmT1vZi9H87o+Gk/0vAri7rHdgDSwEfgzoBp4GZgWZ039mMP/AM4GXu0xdjOwKFpeBHwtWp4WzTELTI7mno57DseY11jg7Gh5CLAuqj/RcwMMGBwtZ4DngHOTPq9CbnFvCcwGNrj7JndvA+4HLom5phPi7s8Ae44avgRYEi0vAS7tMX6/u7e6++vABnK/g7Lj7jvc/aVo+QCwGhhPwufmOW9HDzPRzUn4vAoRdwiMB97s8XhrNJZ0je6+A3J/TMDoaDyR8zWzk4GzyP2vmfi5mVnazFYAzcAyd6+IefVX3CFgvYxV8scViZuvmQ0GHgA+4+77j7dqL2NlOTd373T3mcAEYLaZnX6c1RMzr/6KOwS2Aif1eDwB2B5TLcXUZGZjAaL75mg8UfM1swy5APihu/8sGq6IuQG4+z7gKWAeFTSvExV3CLwATDGzyWZWDcwHHom5pmJ4BFgQLS8AHu4xPt/MsmY2GZgCPB9DfX0yMwMWA6vd/dYeP0r03MyswcyGRcs1wAeANSR8XgWJ+8gkcDG5I88bgS/GXU8/6r8P2AG0k/tf43JgJPAEsD66H9Fj/S9Gc10L/GXc9R9nXueR2+x9BVgR3S5O+tyAGcDvo3m9CvxrNJ7oeRVy0xmDIoGLe3dARGKmEBAJnEJAJHAKAZHAKQREAqcQEAmcQkAkcAoBkcD9f7UKlzOtYUDhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from PIL import Image\n", "from matplotlib import pyplot as plt\n", "\n", "img = Image.open('shapes.png')\n", "data1 = np.array(img)\n", "\n", "plt.imshow(data1)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2. projection onto the principal vectors\n", "\n", "Given your result from the previous question, compress the image by computing its principal component factorization for increasing numbers of components. Compare the reconstructions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# put your code here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise 3: Face recognition\n", "\n", "In this exercise, we consider the 'Labeled Faces in the Wild' dataset from UM Amherst. We would like to train a Maximum margin classifier on this dataset (scikit learn 'svm') with a Gaussian kernel. In its original form, the dimensionality of the data is too high (if we used the raw images as our feature vectors, each vector would encode 75$\\times$56 features). Although it can sometimes be a blessing, high dimensionality can also be a curse. Working in a high dimensional space means that the feature vectors will be very distant from each other which could result in a meaningless classifier. To avoid this, we will start by computing a principal component decomposition of our set of images. The SVM classifier will then be applied on the PCA representation of the training images. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAADnCAYAAAC+PNwRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABYmElEQVR4nO292W4ryZK066QkTpLWULsKG/2W/TT9hg30RRdqWJI4SCLPRcFSH00Wkcml6v+ggOUAwSkzBg8Pc3OPyMzZ6XSqH/JD/qky//+7AT/kh3xEfhjwD/lHyw8D/iH/aPlhwD/kHy0/DPiH/KPluvfnf/7nf56qqubzeV1dXdVsNqvr6+tar9d1dXVVi8WiVqtVXV1d1Xw+r+vrv4o7nU51PB7reDzWn3/+Wd++fauXl5f69u1bPT4+1ul0Gl6z2azm83nNZrNSXfpMeX19raoa6tJ5+sxsyvF4HO24zpvNZkPfrq6uhvLV5/l8PhxfVUO9Xgb7oVdLpB/pgJ/Vz+PxWK+vr2e6GuuXZ5S8DSwrHd/KSLFutldjPNYm7yffeYx0TluT/Nd//VdU6A8ENkmGN/W3H/L/XroITOHsTCjjA9oa9LGBb5V5dXUV605oyHNbKCmUlee4vr6u+XxeNzc3ZwjA41lOerE9+i69OVIRXYVkQmD99vLyMvzniK3yiNiObJL5fF7H4/FMN/J+aqO+X7IucMkklo6EtGp3qyx917i3pGvAUxuYjE4K8ZeUyeOnuF8ahv6/vr5+5/79s4yTn2WoOn+1Wg2Gu1wuh2MWi8UZnUgGm97dTafXy8vLYGyvr6+D8b28vAzfn5+f6/X1tV5fX+v5+Xk4Vob9+vpah8OhTqdTPT8/D8ck1351dXVm3NTxGLUgbeJESOX27EPl0DY4kRxgND49mYzAx+Nx6MRYQ6cc47N9zE23UM7RVMZLQ725uanZbDZ8lnHKaNfrdS0Wi7q+vq7FYlGLxaKurq6GY8mNvf7U9sQZpUNyXaGtDFKGSOPUf4fD4Z0xPz8/n5VL40rGSEPUea7vFgK3xku/OTB53a3y+VsCgzGbmozAHlzReLzSMTdEl+ZlpEbzN6Irg8jlclnr9XowTqEnjVaoS6qgYNSPEZI7hZgiRDMaLQMYlnd1dTX8fn19PSDoYrGIqOsIrN/2+/1g2IfDYThG7lrIrnZJ75pIavvr6+soWLGPrTH333uTwz9PMd6qEQNmJ6aia8s1tTo4VoZ+Iz2Qsd3e3taXL1/q5uam7u7u6suXL0OWRJmSm5ubWiwWAwIn4yQ9GDPUFrr6MTRgGZCj8WKxaFIMZiRaWQtHYxnzdrutp6enen19rf1+X7vdbjhG1IW0hW2WMb++vkY9fM/emUSnkiR6NjaRJlGIS6jB1OP93MQvW8fI4G5ubmq5XNZisajNZlO3t7d1c3NTm82mNpvNO3QlH26l7sjN3BCr3rh96zd9TykqoZ4Q2AfWP19fX8d2MIDTMc/PzwOS6z8htQeD9H49g2xRvO8x4jH5XvuZjMAqeGxGVOUOXmrUbMP19XUtl8sh4Lq7u6vr6+v69OlT/etf/6rFYlG3t7f1+fPn4VjxWAZ6NN7Upp7rTzSAx1AcOYm6zoddX47qfkwy9KoaaJQmx83NTR2Px9rtdgO12O12AxofDocBvfl5NpvVy8vL0EZvh4tz4e8VB4GpcpEBV30/ukqS220hr9y8DHSxWNTd3V398ssvtVwu68uXL/Xzzz/Xcrms1WpVm83mLHptcWy1Ixkk20g368n7lnFXvTfUlOJiXdSNo3FL/+4xFIBWVa3X6+F8Gufj4+NALZ6enmq/39fLy0ttt9szIyd9GVuoaAWyHxEZ85RJcXEabepvl4qX4ZkFBWOr1arW6/XwvtlsarFYDKirwWTe+JJ2tox1qgG3zmud63W3DDgF0L3fPAAXX1aaTkY9n88HivH8/DzovGraiuZHZEqQN9aG7zLgVs428aQx4u6BIvmtOO3t7W3961//GuiDEPj29nagE8wwtHg02+ApLaKuPnuudsyA0zKwGzPpRc+A+WoFNslgOeldr6vVashyrNfrwZi32+2Qsvv27Vvtdrt6fn6up6enIRedcsutdnmeX33rgR9Rt+eBXLoGzL0NGoAUcLVmfsvFpk4wQJvP57Varerr16+1Xq/r/v6+/uM//qPW63Xd3d3VTz/9NKTKVqtVdxEkzfCEjMzJen42LTZIGDT5REj0YQyBU+aB/aJxOiXiuDH3rYktiuVjutvt6uXlpQ6HQ/3222/1+PhY+/2+5vN5bbfbIUXn7fV2qdyUE04eNgXBTuvGjHjyQkavMR8VclUZowZB/JbvMl5RhSmBpSQFTGOuPxkfz686X8pNxlxVZ5MgxQIqxwex6txAKC13T306GieQYW79+fm5TqdT3dzcDJNV3mjMqP5u2xiTi1biJD0k9ff0coVKycrhKlj7+eef6/7+fkBdcV+5Qk9/tdrnwZUbpC/fOuo6nUhGpmOqztHYP3v9Sc9sb0LglKvmhOBiDxGYK4rSt3SnrM3Xr1/r9va2DodDLZfLIXPx+++/D0FfQmMfAy4xj8VNrWxQ61xK14Av4SK9Buq7jJUug4pVDne9Xtfnz5/rl19+qc+fP9d6vR44sHK7LePV92TEpApunPoszufcmMu9voGGE8H5sy8Z83PSbYodEtd05CXS08hlqAxuZ7PZABT8XFVDKu75+bk2m03t9/t6eHioqqrHx8dhhU/L2ETkZGwyZLWR/Wl5Dxrw37KQ4YVfIi1+zO/skBYbRCG0kiblMz2WyvS29l5jGYJEJVq0IqFyQukxCjFmwAIAT8MR0Xme61oc9eXlZQCPRFWOx+MQ8HEcjsfj2Tm+q8zHN/HmS2zoQxyYyvHfW1yIVCFtcfTvQojFYlHr9bo+ffpUd3d39enTp4E+aLeYkCTRBf/NuamQ0JdS9TsROBkaUZmJfqcKVXWGtERgLuU6hXBjcwNP0T6PpTHNZrNhKVj1zGZ/LVAIAF5eXoYVStECLYBokiyXyyEYPBwOtV6v6+npqU6nv/LLevnGe45Lsg/aQlXeCOTj25KuASfeS2ifYsQ9Y1bjlcddrVZ1f39/Zrx3d3cDAjj6ekc9kk0unBthEoXQZ/afVMHPTSm4VGdVvduLQJfv0bzqTUiaYo0ep9bvHA8tPStHzJVObmziMfv9vm5ubmq/359NhF5mRfUK+TWZ9GoZaM9LUSatxKnSS6C/RxdUthCYLoqLEswyjFEFtm8qTeh9ZjluwG7MbsA9akHX6x6OA9xzz72VxaQblunnyDg5sTQu5K/ylPKW6i+PYf0fpQ5TZdJCRm/7Y2qs/mOGoSrnJxeLxRD53t7e1r///e/6/PnzsGixWq3OENzraxmto6HogVOIXhBHFNd5Ly8v7zaXt4xzCoVItIoG3DJMjUsyZo4BA6FUjzzp4XCoq6urgVpcX1/Xy8tLXV9f1+vr6xmFOx7/2j+h4O75+XnY9ZaQuIey3p9WX1oyyYCpVKcCbASP12z3HfYKHpbL5fD68uVLffr0qTabTf3000/1+fPngRMrUqZrbRkwDcg5K905Dbhl5FXnG8WVG+VyrPbj6lhSDhowubZvQk8GrP4w95pcqcCBmYlEL4SqadxkkOLHqlefuZOP5T8/Pw/ceL/fD/WkbZocqx7fZbunykVZiFR4mvVjs08KIW0gdUjZBpaXlMM0UsuFt15OC2h8/r//zja1qAtX17zdpD/UnUf5SY9EV/Jcfyf4pHFhPS8vL8O5+kyk9/ETBZTRq04GlNLD/wWNmGzAKZOQXpRkyOJRt7e3dX9/X+v1un755Zf6+vVrLZfLAYmZv+T5Lb6aENUNYAx15RrTPloP3JwDMyhzNG5NJupGSMj/0+qX61kejWPEz6Rwqof6VED3+vo60AkFeNpvrHSaylHMIo/y/Pw8HOPZHvZJXoxtdBByqjiGyN+9lNziKcng9JmBmzjuZrOpT58+DbRhs9l0ea/zzWS0jsaOtDxGhspdWvrMFTpHYJXdCuIcjb0tPkCO3KyT/Xc6xXLTWBAVicIsT5NAx8gDVtXZFTAKqkXrtIVVE10rdUJun4QaI27mH5Ox4yblgVvImgpPqCtj5wwW/xV9YMqmV09CYHffjPB7lCHRB+dv3ganDanfPV1KH2681PFY6khlJD15uUx3tvgn+8LMB9N85NFCSpXZ20ylc1O7W0mAS2jGJAROA9UzXjeEqjpLkq9Wq/rpp5/q559/rtVqVZ8/f667u7szAyZyUdw9tVCXhkbk5N5XBla+T5YITUNX+S1EV52JKlBSrOCUx41e785zWzpXHUzT+SIIOaqOoV5FJZg5qnpbIk730GD7lWnhOQ4MiUbov/8TCkHpzf5kwNqULg6sq4mXy+UZSrv7ZQf5ai08OOImtJUB07BTMMd2JDRPA9JCZxqu66o1qDw2ZSv03qMS1J/GoxW30IhoyBJOMo0ZL0BI/VY7mL5reRo34p5MohCt31uDpo5pZorTylCXy2VtNpvhs2Yxy1JHvC4P0FJglD7TYHUu+a0f5wirsnpZjaSnxMH9GO9vD7lbE8PLkXhqk1kLbqwiUup/BwQtXOhYGfbpdBpy+iqLlKO1IKN2JV4/NZ02eSnZf3OX50rQzJTRXl1d1ZcvX+ru7q6Wy2X9+9//rq9fvw5BAbfftaJ2/aZoNhmE2sj2cSmX+VwP4lR2iqC9TEfmlmFyQqhc/c9jL/nMunxysW59pgt3FNS4CWEdmWVY2oVWVWdZCZ3HKzxEM5QnT2lE2YcmAT1T0lFLPhTE8bujhivq+vq67u7u6vPnz7VcLocMhDqiOpIB03i4UYYG7EFNepHXJrpAatFD21Yg5/Xrsxu/69D1NjYeXm8P1WnAVef3klN9CrTSXgWVL0SV52J+WPx4tVoNgLHb7YYMR0tSajbZVk8mB3FTIb2q3nVMSW/dr4G3e2KUfAl3neq6U/62VSa/q5xkvC1jbr37efyfx6QB9H5VtbMULQPvlUmkbfFj33hEekA0Xa1Wgy615C5vRz1PaRfb1pPJCEwFkz+lSpX8vrq6qtvb29psNnV9fV0//fRTffr0abjy4ubm5sy10j27UfmignuHxN8SPXD+TGX7MTQEDS55dGsCJd7bM37q2nXp//P81gC3qI1vvNF/ft83obBzXrp8xTb6Tzn8qhpim+fn5/rtt9/O8uzeh16fW2sNlO9G4LHKmfSWsfKWT0q7eKDgA9xCyCkDyImQDKCFxil4bFECn0heR3qpv5IxRPX6vU0pM5Emh8crzDRQSP30H99JJ1Q2Fzik8+vr63p8fDybBC6p7wmUWvJdHNi3NzIgULpM2yF1/wbRBhmuGt9DK0dJpw8ebHBA3DidTqQX/2f/W64uUQj2i23RMa2J77+7zlPe2D/TUJPB0yCJxikDQL2yD9wQpd1rLE9IrcBdt/ySZ9PmJx1LxGcfpqBv1QVptES4fRWGl8Trng5aJhadWC6XQ9k0UDcuLSroP+78Z+DG6JkTI6XLuAgifsZcMF+p764bUZQeJ05BosvUwWJQ5fU5otNdO1XwjBGzEzQs1TGfz4ebnlRV7Xa7d/siNO78jWOrDNDT09NQT8+Le3Dfku9eyEgIQaP269p4RWzV+0vM3a1zwLkngcdw8hCVW+UlN+4o7e631V+JzmmhXstoL5W068zpiPeZ//H/RD184UAGzjKZheDk96tzZAOn02m4d0dVDd63R/+SZ+/J5LtTOgK7AsSBrq+va7PZDJcCaU+vjJduzg3UUZnkn5eu6DzejZG7rqhk3rnc9wOzLvcELeFyM8tz45A4JWlxa//sem5x4BRDtOrUcdKX6mB/3ZjThBEaX11d1X6/H4I9UUQBl0BGVzxrz4t7ANrSpXIRAtMF+RY+Zhbu7+/r69evZ5fJO1LSKImw5LpSDo1G51a97Vd1OkND1H1zZcy+j0ID4kasAfR33kiaBkzpZQ/c4Pgfj2dcId17+e5pnK542TIaUgWNncqlwepdYMAAjUiqhQ15Xm0L0H9KnQrMxKu1uT8tK7eWul0mG3CCckdgzj6hsT/nILnzXnCVsg/JDdKI0nIzjbNFVVL5/u5l0RDZFm9Xi7r4fzzfjai1jdInh7ffJxjpV2qbUwJHTE0ATnZN6Kq3+xrTNuQx9d/3Iq7Lh/YDM9rUfcuEuroJie9zkDj3dAQUQvcMOKFWC92TwonErSDOB7/qLegRkqRjJDQq1ekI2BNPNbr+SHk4mXy8XJzX6rjZ7G1/MH/jeWoX8+e6tEgpNl6eL7qhPTDyivJkKrfVzp5MvjOPZhCzEHwoipaJF4tF3d/fD1dVcKY7+iU0Y+bBOWvif77wITevtXsiE/tF1EoG78qjcnnFrhtV0puO8T64OMLpRTfuNCEFnz1u2UJ6GhBjC/7P82X8uhhURup3AaIXXq1W9enTp9rv97Xdbuvx8XFAd9Ib/TYWwFV9MAvBXfrMOPjl8FSAK7KVJXCjbbl1pwCOujQgR5KxenSct9+NmWiW9NQ7348d+6w2JRoifXqevicJ+bysVLfQVf1OYMR+k2aSDvmYXEotLjJgNUYIpOVD3UFduV8+X83PT8uwyQidz7UM3K84ljtT4OZcjnSG9MUXSCREQ55HfaTfE0oycGzFFDQ+Bsvijp4lcfTVMSxTkhagkl7YL9bl2y+pN6cT+sw8ve7pfHNzU9++fTvb/ONezzNdLbnoigymYJQe0eIE76KjtJmf78jpKOw8OBmt82YuTIg27Pf7gV+1BpDIz03xSXk8hpOQ/WN/PGdN3p28Aevl5nBlB9zwxFF9YSd5Mz+X18j1jMQnvn5j1kft0ELH4XAY2n04HAavrDIUL+12u2Fjl+4pwbpYz9+WhfBChcK8O7pz3sSdenTAFddrSyuLkAI+tYUu0cvoufge6rpMyQxQn+k73e/YALb+d72yDWNltgyX56bx9N/SGJB26nKllmfo9U9y8UIGN6nrsiDtd+CjXN1Vp2ifNwhhtiHNfnaEOWPmZLnA4G6U6FiVL32XpKDPEdoVq/75SpPvsU10RnXwnblXShpw1ek3ZFF5rkt9ZoaoFat4+xyUpBcG3AqihcTM2qhe0on9fn92kULVG02lJ09y8TMyuEGHBiwu7OvXbsCM9tON9mjEbINPCCXC+aRKJcadJyba4FxRZTt6+oAq2k4o7e5PxztqqU2tScCsBev1zzQ655Bsu09OcuuqtyVePzYFg2ojL9Bk8Hw6nYabY1fVQCu0yUvpNRnw9fV1PT091W63G/osTznFgKffmx8d8hSJo64k0YEWr51CHXpl9Or0c/nu5/TO73FFGZF7EW+jo7uX79wvtXfsXH0fy0b0aMAUcc/k/WxROk1yGSm3AThgjfHgSdfEsRJeDqSb7+kJ7y3FT+WqPQX2jJZCxfQmSsp0sK7kYvUfEZJ56tY1b8kjcelb70Rxld+S1C/Vw9VPttf7ySBXvzPXz5iBE6SVWvMbIIruMCOlNgiFq+rsviD7/f6sPb0JX3WhATPzICPW/WS1ScM7xg6SSrQMOEnLeEkzpFgpnaiYaINzYBdmBDgxuTrIzAc5YCvCT7nPZCDqG1fh3BD9N7a75Q01yarqnV44AchVpc8eCqbMkLgvX87BteVS98XTeWzzmFx0bzQZMRuV8r0ul7ilKeVMdS88r0cTUuqG/WXdPjFYRoun8hj+z8nh/WAbVE/aUskgke10SeW7Ph0YpBMGoaqfIKGy9LunG937uc78RRmznUkGTAKuzeqiDzTkFkf8HmEe1V8cSLrKNBAtxKfQbTKg0Tt3w72+vp4tUzNAIYK5y/TBIyK1EJOf2SfRI9UjY9E+hGQoemdQlIBHkyPtheBvKUhmX2QTROKUPdE2TN7YhhRmDKAmp9H8uW1KpfHJNzr2EqOd4pr40jlEMOUTKYm+cDOPI60ULC6mczkJfG+yt4FXWnOS+QaiqnOjTRPP6USiOSyXx9HAHGll8LykKJXN8x3p/TyOoXiv9gzrLpaccK573z/sMUNPJiMwFyxat/73ipP78AY5UrQM2gcxKZZunIOpPrB+pwlUbuoDz+GWQKdVvTiAOnJuzfrS+eTDEuac2Z6xwMclpdD0nf9p9Y/ITSMmhWB84XGH83SnpVx+/hACqyDtMNNThHSpvPb8Ot9z4+UScc+VS5xPpsCNx8monGZ4ioeoQxRQbtKNie3Vseqz9KPsCydqawsnE/oCAorzRi5mJORmPralpxYV4XhxvB0oiO6J5zOD4p5YWQVSHZXD/TLKRtzc3NRutxvdtUeZ9JAXXW1xe3t79rjXNPCurF7qzBWsjidUauWKpRjNYndPzkETAvuNtMm/aYA6ltxf53IlUBkKTlofcBmwT/qqekdz3Hh1jqiTjmF/k45bRtEzZKcL1B8DNxqnjFeTWXSCL5btTwzVRqyxvdZVEw1Y/I75ut4Vo8loWhRCDefMTQYo9+WczgMnns8J5APC47zt3i5yYlIoBh3k4Mm7JF15hsPb4XlrCdvTardzVjc8nuPluLHqM2kC++TcWMZXdb4qSmpBe6An9cD0QxRCO4n0xPjPnz8P69c9I3b0be1zkBJorFVvSCPDFcJJCey4Zjo/q06dT4VMca1qF72LUEI3avGFm+PxeJbDZB2cLM69q85vV8p9Ib4w4nzXvUVqN/9jPfw9GZaL16vz2AcaI4NFLiczZnBvwsvRvM0tmcSBeVcduc8xgj2Gwj3q4G6LCknneACX0kkyZh/wXnsYJOk3Xu9H4aRjH3Sei3sdD3TdcN1FtwzYy2b9qsczFE4VkoejHlrIr3J5RQe3xiqDI+rjAZ1ASJve/1YKwafS9IyXym4hbhrklriR+j5jGuoU1z3FtU9plx/TQi5NmpQ2S3TIj0mfnTcnfurtbLWPbexNut54+zFsl1MIekJP1aU+jskkCiG3qdxvSgGx4Wxk6mBScDqGXEsuR0ia6nPl9ZTBLIVTAReij09Kp0v+v1Db6+nxfgamrluXXoqypWPFE6kdrhdKb9yoO3lAUQVtfBf6Mk9Mrl71to1S7exlq6omUgjfuN7LYfrgfkTSoDoC0YUS0bycJDSQqvfBSEt6Lt5danpAeUIbf3HS9PrRo2U8Jgn72jNc/42ZhySeB+biSuLdx+PbiqpPolZ7JJMohLs3iSv5Ut7bcsPJ7aUJ42Wkengele7vrlhODk+6e17TEZiDozZQvHxPLdLjkGu7vvg4q1Z60mkW++q6cjRMXLfnqQhsPQNme5IHTotKLZmEwEJdj27doNKAJB6cuJsrNXEyrgrNZrMzmsLMgyvFFzCo2Krzy/39Flan0+mM98steqDB3C+vsUv0YzY7v/dCyqzwveq9tzidzm8fwBtKt4Iy5rRbkz1tIE+bclxotER16UvlJMPmOdLtVO89ekVGL7igOBpNQVz/vTcre661xwPdNToV8YFOEzG5fPFIr7fVd5/4nLStyeoTx+mNl+NpQnmPFq/3svjuOkzbL9MYqT8tQBrzzERgn7RJJgVxCYGpCDbUvwsZ9VsvM9F6TwpqdZ5tojLdM/iAsE5SJ62++dZRIrZ2q/GmLJxoSWc6l+111PZgKC0g0ICZ8XA0ZH80KZPBeR0+eaU7jrHOo0fhpGLq7OXlZVj2T/aQFqR60jVg3zHPFSg3FkfgFJU7T6NQYURI5nYTEjh3dKXzGN/Q7YNApJMyq2q4YZ3cG9u02+2Gh5roAtMx9yf37xeeOv1QG7QSmlb/Utu5gpUCQU9FOs1geTwvjSP3efhSviaTJo3s53T667arrV162pn24Ys6HYla6TOXnqvgq+XeWm3hscmNqa2pPY7E+t3PTYEeXRonrryLr54lT5Voj09a91hV5xc5qi2kP61gsZVJcBpCHY1JK2BXOz29ST3JoH0SSN+keVPRt+qCPLAvYPggJyNJv3knqBy6XVdYUmYLcRNCe50px8wy+C4aoXMULOkGKro5B9Ez5aqTvqredqnJeHV1NYNIohn5MB9jlYLiFDjp5cGZ6zKBi3tF1Tebzc7u4i4d6Lfj8TjsDfa2tsaWAWBLuga8Xq+r6i8qwUuiJRwYdytUEpcS3cC94Uk4e/VdSktbAJ3jJt7NqwS4NMwJyuOJjtry9/T0VI+Pj8P9KZQR0MSX4aVgTRSCdEXGq+vrNCG4n8DpGzl8CnocTHT82P3jUvzhRsXJ6Jt35GHVNxmz9OyBntuA6vmQAfvGYiqp5XISwiRUnCItd6f/WvW74bXOZ5rH9zZ4u4WAvBkLP6eJ0qrbj3Pk1P+8GXfPQNUfvzRd/6fUmRsQ3XnSu+u7R+lIC5IdjOmJk+hDQZyvjrgkhRB5q3LS3js8xhfZUVdCynAw2vYlaP0vb+KIwuCKhilkORwO9fT0NNAGBXFENHJp6cONQmiroItRuj8Tbz6fDwjc8jROqcYG3z3g2ORI40NQY3mMBfgUI+aqqWf1MyH9hxBYqJTu+cAgwgc5cS536y3puXA3XufTjgJULkWRrrvC19fXgdOyP6+vbxdyyoC1iLDb7d5NztnsbaGitaNKZQo5VZffH9l5MXWhMrx8prT4ubcwQnHDJU3xiSIKlzwvwer19fUsC6GJKl1okjNl2/KMlO+6vWrrv4SU6b+/Q6aU44OVkCNNCmYU9L2V76W38UFMgaEkuVGf8AQDtZsGwja48PIlD6wUU1A3U3ScwCcF3zw3xQAJjHo0oiddA9YzvdbrdR0Oh+FpM1rqc0NlJxNKuvQyDu4CW/9z7ymRyet3YyY9UAbh5eVloAd04cfjcUBaoa5QRC6RQo/jCMxytQGeCwB84AzRNQ1kevwYsxUqmwhMz8TfW65auhOV4bj4amECLtIJ0S55Oq0tKLhjCtPb3pKuAW+326qq2mw2td/vh2v3yelaL0c0dUgyNsOYs0zcTpxWnU5ZB3ftDEDpph8eHgYj3m637/YUcMFChuWcnyivwW5lSDw7Q5RKBsxUn/NoegGm1CRcwUt0QoskSb/Ok2X4opRMq6rfVflWAuqPjpMBv76+DnflSRfdfsiA1QgGFVNcd0qPJBlzD71jWhmKZMTkpikIYvCUsgqtO2k633MEmUK5WvSFv1EHKY3JcUn8VeJ6EEfX0q5nmZJ37WWGGG+k43oAk+qbYh9dA9YtL3XLeK1jJ7R1qpCoQysy5nFJ8T4TNVNpbHzqjbIDbJ+QQ0EBg6aHh4eBIvBp9bz7O3mvu21XdtINB1S/8WZ4zruFWHTdKQ/MVTvqim7eUYxtm81mg/dhEKV2+jlcam4hpC9Bsy0qRzqtqrOAzoHyQxxYBkyOmDhdC/Va3NclNdjPdUW5C39+fh4WGHa7XW2323eTYLFYDOUeDoehT4+Pj2fZBHI3GazTBh7H9gnd3IiTAcsgdJPAqnpnwDIsZhBaF9My7dkyeLaBuuGKmeeTObZKh7nxCp1VRlr6JvJLv+qz0mjyCGrnhygEo93kNlsu/BJx11eVN/akc8glmSnwiUblOWd02qCX0wahpfc/IbDaxjih1ffkQhnoMUcsY3DUrXp/A2xOKA+Ge2NE40vtpBHqd++7UwlvQ6/f/H+KdA1YqPD8/Fzb7XbYE+HGnBqTXlPcrG+EScYrw9rtdvXnn3/W6+vr8NwxobHazlyiXJfOFeput9szhJWRK2Kma080iIGb6tG7jND7xbr8pihEZv2XEJdl8iYrQmxfqElCJFaZjtLuTVmHvrsHdW/gYy+9zGazmBOe6sG7BiyOIn652+1quVxOMlo1lEGUo20y4KQMH3h1frfb1dPTUx0Oh8GAPXXFu0YyMhZtED3gow/0rsyDFH46nc4MghtruOGGgRENQQPruiLC07BlwC2PRzeuTIKntjiBW3wyjZ+jIcfVsx1OVbyNTmMYA2jy87Pz4J5MelKnKnO36hW4UbuCUtn+n3fWyyZqMaDqddgNgLTDl28VrAn5aMBqX9X59WMy7BZi8JyWntk/eoHkVsm73VXzQlzXKQ1+LDjibatatKVnYASfsYnDfrPcMbpTNWLAmmn7/b6+fftWp9Pp7MkzrlSmn+h6heQtQ9V3IVpK5TCY0g4wPRyEQQCT6mqX0wOht6L83W43GKsQXb87x/X7YwjhWR/vFsS+ORrRSDmZiMoJxTngOl7IO5vNhvs4X11dDbfCVRDLR8G6u/d6hOrM4xOdeVEpy2QfW9SHdkIb4URx2pVkUh5Y6/9VNTycTqtyDFR8NnFg+L8bKjvdmrF0qzI4bWGUISSl6f/T6XRGFbgJR+UcDof69u3b2X9aRVJblclQWq7q/bItJyoNi+I0K6XRiEo8T79LB37MYrGo/X4/ZAzUXuo6XUGh9lKYWfDFCQV01EFvDB391RcZsHuxHnpLJj9mq8Vz+R+V3jo21dGasWliOIpz9abVRnL2RB3E8RX88Upfd9Hk6zQwum83qESLegNDA5F4BoL820V9Jb+cz+dnga3qYcDnbWVZSVJM0zNcL1/eiuX5+I7J5EuK6MZJuiXOI9PqHRXBAWjNWkbkmqXkorpWTUuSrJ+XuPMRtELkx8fHYYHmzz//rO12O+SEPUhxSsBFHVEXR2b1mQisAUr7ETiIHEDuzmL+NLlxckgdo36J7shY/Snyul8vd+mNcdeqt0uJWHeaeCrPx5k0QoGwyvg/QWCibGoIqQOPT2W3KIQieQ/cVL7Oubq6Gu5kyJUcRvQM9oSwyjCIinz79m1IwQmJ3fBUH6NnIh37xEnLPiVkTkJDr3q7pEmBlY4hgrFOBtOasLPZ7OyqDj7jT8f78/7GFhHcw00x3l6OmVSzx6Epk7IQU6iA87aUjUgRcOKGdPk+AeiqOchsA7MU3AElo+VqHZE9RcKqMwVhqpMbivRbS2ccSI8FqOukG9eTdOp1kN+mO8FX1TtPqgWT4/E4BKZCa/Yt0TTSI/4m3RAEemPNMqjXnkwK4lgJFaxGqYOa7VxIkJKrzg3B+ZtHuKICnJkanOVyOaCL2sn9Grvdrh4eHgaqIFrw+PhYT09Pw+9cyNAONA6AGxr3CaQ2a0LxIX/6jQMpdJfBKFdLL1L1dkFByspogLlJXO3Wc9dms9lAsarqLPPy8PBQ2+12MHK1Qc/+053oeeUKx1rjSRTWixTEjdhjJ3kV3uno+fn5LNfek0lLyR48UdgBX8r1BqRsQwoEqt5uiOEGLIVrELkf16mCMhYy2oeHhzMDVgTv3No9BNGI7pWeQW7aqVQqc8y1Eq1SoMaxYBCr8pbLZa3X62GrJBefuFCj/SJqg4yIe79lZHq5ATptYSZBk5VG7OdyQvJd+hyTyY/Z8oq5dNpzmfrfUyQc/KQU7xB5JAefhqwg7enpabjsR7lirmwxGHWaktCObeotyyY3mSarfufqnVMjD8yEVJ4npX5Vj4xWMQLzuTJQoa4bjLyaOLHKXq1WMe5RnU4FelzaDd/pYrKXllz8pE41nJ1ODXCaoUY512H5NFjf/ab6mVNVQCaj/e233+pwONQff/xRf/zxx7udZrqPAz+rHz4x6PLULnFFfziLjMMN1WkI9SGXfTq9PXeDxsmJq+8UGr1Tnc1mU/f393V9fV2LxWLIkAhFn5+fzxaBFLhK38pGaNP5YrGou7u7dw93Z//VJl9I4vjyeI4/d9wx957KcbkYgVW5BtQblGYNjVgd5Wxl54jmRBzPUug3HSe6sN/vh30RzDYcj2/7a0lzdD4niSvOuXFP0sA5ChMtGRPwBicSN3z95ueL4ihG0EMoSSGWy+WgT+fOngbThBTlEKVgmyjMTI0hJ72w+puC/5S9cplkwKyQaEtu52jayzTw+PRfj5LQyPf7fT09PQ3vQhXlfavOg775fH62wuYbaFi36nIU5uDRo/juLyKVBySiIgziiF7kms4fieb+md5JE1PUw+mIAj3ueaAwhSjOrGNSkJbGlS//399pwMx4jMmkhQxVUPUWXBH2+T8V3UIs7zSVzg5xFYlcS4b3+PhYv//+e2232/r111/r119/PctGKCW0Xq/rdDqd5YFVDssjysq1k3s78rD/Mlq5X7lu7j/QwOhYcVPd9Ui8k1tWXUgbnD5wOZdL4PQoKlvPvNaCzH6/H/RO49I4KBXp2RU+YJJt9KDNgY6AIb1wwUqU50NZCBfnup6no/GSEyahQTjXc/TlsZylog16kdcSbWQoEvJFn2Sqd0oOkn0mDfAcrwegzo09H8zjvd06Nxmxxyl0z2yvDE+rb+62GePQiHW+04XkhV036X8KPevYsZTJV2T4phwamCuwhbwJdWnkbkjkwGqL9v6+vLzU77//Xr/99ttZ1kEzV3lQ3vdMZdLQdLwrTC6e6Ob5YLZdlIEPRBey8nG0kpubm8GLLRaLoe26QpeTaDZ7e9QA20WdcgJweZdjxXGUp1A/0s4zBrcEJ3osUh6fAKJHnEDJq9Boubk/TWKXyQsZnuN9eXkZVmxSpJ2QOSGyBzg6nyjLjj4+Pg7Zhv/5n/+p//7v/x54sDq+Xq/r7u7ujIsyu1F1vlolY6EIuZnkd24rY656owFXV1e12WyGB0Jqn4H6qvN072Ve3Hh1dVXb7fZM7zIC3lgxcUr+fjqdhkvYfTurXldXVwO1YvaD5WqsFQBLd8xYEIU1CWQvnnJsITQNX+lQgktPLlpK5oxO0iLu/t+YtIKoqhrSPrpCZLvdDoNOHi5jYopKRuu8OtEFd//efhqjUIg8mEZP9BKy6Bh9ZxbBNwPxf5cWoiUv6FSC/J3BomiCfueV2kRz3xPuY5W+s/7UF1KeNLFcJnFgzSihhfaaVp1fsuMG3KISrQ4kN8Qo+HQ61Xa7HXaPPTw8nG1WEafj/YypCLnNqqrVajUgSrpshy5WZdLA3LBubm5qs9nU9fV13d7e1mazGeiEEJhL37zkSBuSTqfT2cOuGeVPmfgt/bYCUB7jmQ4aOB9oLmR2dCSSy3uSYqb3JLwGkItlLZmEwFq9qaohWNIAemCSOO0UvuvG67xJKPDw8FD/+7//Ww8PD/Xnn38O1EEPY5zP58Navg+iuDHX3xmsqC1EyarzoI8cmK5diwc3Nzd1f39f9/f3Q30MInmXSVGX1Wo16EQTS3r1m0a3xmnMwD233PKG9EqanM/Pz0MOmbcv0PEy2FYmSr95oOp90FgoIzIlkJuchaBx+Vq/KyB9bynLG5poA+ulJ/BNPgxwGMlTiYzCWze1dneq8lPmIJUp6sD9t1VvuVp5DKcUNBzJ2NL1VBlD4HQcA11mNUirUps4bpe0kbr3rERLJu0HZiJbuVSlYdQhdd6jVjaOnfZjfIJwxUxLn7q8X4a7Wq3ql19+GfKayjxwpjtv1+8MOGjonv5zI03ftbF+s9nUYrGozWYzPNHeua6XIyE9EaozVy0dcZJRWmWm2GNKPCIjErpTb9wXoTLYblGq2extO2fLJkg16Q3H4i3JZASWy9AlN1rNIU8Rb2KOMAlncEJZGq+CNu0e40321ut1/fzzz2c8tKrO8sFUlgZGfFltVYRN3qlj9a7PjqgaoOVyWbe3t8P7ZrM525wzn8/f6Ur6IucUb2YqTa6b9Ip5burV20jjk3AiSfdpjKg3z2/72MqAGdDqc2s/MsdetkBj9itjkky+IoOdSbyxKiPwJS5EQjRk+os8S4oRp6XrVRBAA0h940RS21O/icb+O4M5T7kpSvfov9WWVrs8GGJ5nHSJbjgis44pQVWrnJQSSx6mhbzefv/vb0FgIg0bQg6sdErVOWfytIt3wI3fO8k8MNHn5uamvnz5MiDw7e3tgKSq8+bm5mxDfaqbk4PCRQ2nEnoXNRBtuLm5qdvb27q/vx/24sqdMjJ3r0Q9kP9Sn+6teAUJI38FW6rLvYbE6U/KJvBYF42xZ0c0cZ3HM43IBY5WDEUDTvdedpmEwHQdKpwGQJ6UiDgb1nr3DpxOpyFl5wb8008/VdVfqbDNZlPz+XygDYpi6Q08MKQiaUCOmERc6kIIu1wu6+7urhaLRX369OmdAbNPjjLOZVk+U5RqozImcqtusOTvMsoUcDkqOmqPrXy5TsiBPdD1Y9hWN2AfI4HXWC74or0QLjQ2d0c9DjxWZsp4kItx0wzpgwdeVGJqIyW5Wv7Hc33RwbMPRL5emrEnfixRjxuo/Opwz6R4VoDUjnx/iqQxTlQhTRbaSmvyJoD5WwzYG+azxSP8xHmcKjhPVmO1n3e73da3b9/q6elpULYCnLu7u3dRt9CILrLq7XEIdKtSkEfAzld5nNyjgrTr6+u6u7urL1++1GKxqPv7+7PFCwU0olhsl8pMARYHPXkAlafgU7cUUNBL96y+Mm3Hye0G2Au6+XKkncJ7k91IuHdDL+3l/tsR2GdJMs50PL+nCaFBkAE/PT0NBsxgbb1e1+fPn4cNMNr7m1yb3L23kYjuSOT90zlCVqbL7u/v69OnT8MVC1qJExJXvV2YqUkgXkddOc9OBiyjYf5aXF8LTeoDjYArfm5waqOn37xdnMDUhX5nGS0bcCGF4Jh4+pR7YZJMNmBG23xXA9KA+Od0LI2FKTQGih7o+EAo28CB14JBK1hxzu5cudX3RBtIHdwzuC56NIa/O31wt11V7+iK+pcmRBL3Qv57GjfXs7f5eyV59b+NQgh5tDlc+U5tGeSMZGd7L0ddvUQfdFd4umBXntBVG3RWq9VZIlyRrAZWm91loHS5vs7Pfuhyc70rcOOeB+068+2O9A5uBJx8miB6UnsrpafxqHrbo3A6nc42xXPrq4xDbfEyaEDJ8D1NyCxMq30ubhfenuSJCWI9mbQXYjZ7u9eA7naoTSpEnoQ2nKlEQw2SFiq0x0LvvHKCWQ6iP6NzXQXhCtDxLW4rd6X+yohZj1YdRR9ub2+HFTdNYl2D5qtO5N3qQwrs6CX8Kgc3aJ/EVX/xyGTAvqfa6VPP2HxiUR+JMyfkTvZAg/V38mC2vSWT02h+RSpnU8v19oTHstN0w7rZB4Mgns820Dg4UcboChXGYJRIIwPm/uCpgcsUScbsC0LJYBJgkCZx0lSdP7k+Bd8cA71737yPHv+04qEpNnJJRkTSNWBtA9QG8fV6PQQwy+VyULQayKVTd8OtTpGziqLMZrPBleqKCgVxXJZ0ZGZwIxTlfSCE8IrglWf0h6loWVS0YLVaDdkGeSAZtW8g8gnTGzRSiBT0JQrgBsXJzf25erl7ZlmpXWMxgLeL3ixRyKSLsVW2S4Bg0rOStVVRLlPukqmPdK+EqQ2i8jVpZHzutoQOPtMdtWnoHtnKaLVAoisMdL7oiK6q2Gw27+iCu9WWe2z117MBOt5XAkV/WsZLA3ajdApHWuIURe9OpzjBkkHSxRO5XUgPaLwt1J1qxJMvq/fNNTJaNfh4PJ4hcEIQlueNJYrQbXMQ9FIZnM2kAuLW5FG8lN4XSJiikpHKy5A6yOgcTTz61+fUb6csPUOnh2MZU6gEP6eAjb95cMcFj9aY8TyNN8tJCEz6Rl7PCXapdA2YKz7axqjdaBps7brixZPshKOCIwQjXB2vFBhzmBK5XRrm6XQ6owTctcab++m+Ecw6XF293R9Xm9H1XVRBqOzZhRba0cjUVk0abkpKKUI3RpWn71pilRFKb9ShZ4WS/nwSpgnlga+L0xKfOBxvoTsBx20lTdgxmbwfWJe5ELWqzp8Z0YtCEyeT0EVp1Yi7udhhj1ClFN4Lgvd/SHdfd44oY10ul/X58+fBcPWcCQ+yWp4lBVxsf8pvj1ENvqv/TC+Sf6Z4wNvprjvxa580qZ3uTSgtdKfXSwbc00VLLnpKkT4z0tdlPCTzOjcFFKmRbiApCKBSfEBSUMHfiVC8RZJoA7m9kNh3aKU0lt8QxY02IZwbr5+X9J905bEAU2fyZC0w8UzEmOEkNHVgaSF3r39efut7TyYZsNwvZ+NsNquffvppiMwZfAhtiDzeOEcNXTk8n79t/E7KpcuiayYa003LSIXyzHKIHog2MF2mfjhysv3z+fzsAsRksLz8iV5CfJz1JOPmOLBejYV0q+Pn8/kwSZ2buyfTOdQp6/Ex1cTotZcege1ySYEe0T9RqiST7w/Mz1qfXq1Ww0pXmtVERIl3UA1nB5jPdSRJCuoplIpgJkMri+S9zP2qnwwKPeLmSl7yODQQTgLSCOk0oZr6lwaZARGFKTleAeIyVg/HxSX1OVGJJGxPoki935NcvJmHl7do6beqarlcDlfQ0p0RgZ0GsKFECybeq85dHpGWQRxdNI2V3oFcm1TB757DsonA4p0eLLWQggPrxut80Ce70yEK9UiUVJlXV1fD1c+eudG7A0xC+zHXP/afH5NAS/+n4HWKTEJgpkien5/r4eGhXl5ehuXU7XY78GHlYdMt+91oJQwauDTKqFUTgtkG3iGGtEFUgcIcL4NE/ac+yuUrAFQbdIWHJofohq9OclBkqCyPbefl4zre03ykJtKf2k76wKs0dB5XUD3A5Lk99EwTKx2TPvv3FrL6JCW9TJdJUSZxYFbKiF/pqqoaIv3ZbPbukvcWp0mNJoXgcaQIjoxOH4g6pAW6W6S3hZNE2xJPp/Mn1hOphOREX9XnAY9TBrY7oWDLPVOIqNQVebq3Se1KGYIpSKv/p6Bx7xhHWurJj/swB06FcdWJyKRNOBwkb6i/t+qiy3chwqSdSzRgBWt+ibfq8XQX2y7j1aTUedrQrqVm7v3t0YhWoNlCQ/aLHoh1OK92j0XjdjRLxuaTz/vRk5SNGBPGAEk+zIHZGO4DVtD28vJSDw8PQ85Ut6HnzGlteCGysKHMu3I3Evkpl4D5qFlybQ2YsiR0p4kDqi5RCD7SVhdtKuD7+vXr2d4MonzqJ2kEH8CY+DsN1h//5V6Jhs3HzvK429vbwZjF+3sA4rpO9pACdn73Me2V4ee510iX5FMuQmA3RhmRFCi0UgaB519SD42/hcQeHHHgiMIp6OL5TA26m5eR8XwuMyfq4EJDYNkpeHNDJp3w5VeWrf+dYs1ms+GJQwKdnudTmakPbqw+Fq1zx8T5/aUyGYGrzpd+FSw54XYD7LlVlut1emBA4+LAet1VdRZYeZaAyMUBJ/LyLu9c9fKXGyfz2N5O0gZODjde32Tvx6ou6oZBnOrW7jpmjXpXjRA5W8bEcdUESePWSoux7Y7cPuYa0xa9kEy+vSoNVYn+ZMCCfD8nRZ6kB6269TkFbVXvXU5VDcHafP525TLdIimDsgMyWl4Not/9kVXiveSo5MtugHxxv4YoBNFVHo2GR6pEb0EjTPREY8Sgk2k36o4LR4liuJHrKpjWmPG8KSDGOjixxtJpF+WBSSW88F4+NMmYK/NjPDL3czzQdLqT6ndj9pU91tmajO72ucjg/XBX7CtxydN4kOYGzJx50hGpEutIxjHGjf2Yll4905BkCl/+sAGnNAs3ihDhFO37NV3JVbmiWrwqBS9ECwpRhvle39RCdGSmQQhM+sANQXosK3ev6b/T6TTsE+YEUpaGE4R1Em1JM2SwQmkFlI54nKgMeJ0OcJOTbvrCfdctI3GqoPI4adRfD8g50bnN1ukPx5XSGmeXyQZMt+NXIsi96kU06NEDooQ31AMeKoXBnQxURuuBG+siisk4mc/WZ1IIZSF0fwohn/ZuyIA5aRTkVdVQLieKjHO3253xcRrby8tfj4JVO5+ens5W19RnLn27nqre7u18dXU1LJzIoHR3UV6kKaHx0fCoU9EoGjq9ldroK6saD5XZWqyYshp30XPi2Dl312x0y22nMhOB1/cep/LO+YQiv/NZ31pY8FwtKYYMUPdikLHxaZ8azKo6GzSWxSBSA9qiEAzImDdmv0iF3CVTR+y36iSAjI2ZoywNtYWSYxmaMZliRxdlIdQoPbqUwZunozwn6/zRo1BHS3cvDA519S3LU3qr6vx+Cb4hh4ERc7OkFMmgtaFf/+smJTJgtlfXzbGvDNDc8DwjogBOi0PMQrS4uMZBOuclQB7sCTl1GVUy+Kr2AoOPjeqWHlrxRyu+0cRnO8b4sWTUgGmYVec30JCociIIMwwpwU/pGa/KZ1rMgxgZuDIPbL8MligqbkoD9hcR8HA4DAs2Nzc3tdvtaj6f1+FwGGiA6EHVXxfBap+09hpX1bBMzfY4+jJLoYwIqRQXaahLXh9HCkW3rQki78EVTx47Fsy5AWtM6B1TdqhXDtF9LHVGGaUQHmylBum7u8nWcS3pKY4TgG6aA5qumuV7ys2myD9F8kIJpsuIjLryQ20lF5f36GVS6L38M3Xn9M1pmxsQN/2wHnJQH6+xwKk1PinoGxOCJN+ntmPyFRmMOP3mHaqYuU2JBpLfPQnunDZNAg4IuS43q3MjuiMP6YAvGTudEAKrbVygUT+FwLqlK+mEvs/n89rtdvX09HTmhTyHnVJkVW83JjwbMFzCRSMm6uqYqrcnwVMP8hicaETRJP47g3p5Q6J5z4i9n+k/1tGTyU8pkgE7CkqRx+NxGFgqNmUjPCDwdImjn9qg/zSwpA2kNr61kIGQX1ZP7usbeJz+qD3qk3iqEFjB3eFwqMfHxyE7o/Sa6IQmMNNLjB1owMNAIdvi+y74+2z2toFJulYd3nadp5VGHT81ANf4MRvB3/nu4unRqXW7XPSUInd9rkTyNP0/dTbR9ftnvTuP9kWVMa7lFMFduiO/1+P9oatm7lffdZwQioOWKAR1QXQk5yVVqqp3hsflYpZFXbJOcvKUd51C+1T/pef+HTL5kiJmFSSkE0Keqjq7BF0LG1X9PRAsk5OAkS3PO51OZ3cI0kDpWBopaYNQMm2G9zSV+qGN+8py6BhHYAV0SeSN+LBDn0DH4/HsSuN0JYlPWp9A7o3I7bUQoiBOE8uzFyy/Fcv4mPD3ZEOU7+HYLbnYgIkwLQOWcYkP0oB7WQj/j4ZIRNFnPp9Y5XsQ4IGar4j5/gpSG5XPRwnQC9CA9/v9sDBA8ayGeyTfZaayfVsq4w4aLfciO+2pqneIrTEkTSEf9nGhLj1u4Tv3WNBTa2wvySy4F+vJxRSCQqW1ouqpkmiD15M+J+rAsmi8KQvhG2RYPw3CXyqbK4++0Yh1tNCKrt//52Ypehn/3NKZ65PI6m1MGYn0LukFfB+R1kRpyeTdaHpPg51miitkiviM94jU6+HVFkRbHS/j8b0NDOKYkZBLrXoLivhSpM8HEeolCsGA8XT6K2OhoJOejIbi96FwdPM8ulOIBCAST6mJ0nFlT9SC5XEc2JbWeE015tTO3iT8kAFTGIh4IOcpMB2v99Zs7bkIT7WlnKeMiohCVyUKw3yt77Glgavt5JEpbaftlRRtHGdGgostqlvGxkUDrizqN/aTv1VlOuG6p851HK9GISVk8M0UnXh7S3ops96kSr+3JsgYjfjQU4pcWhwqzaSUcvH/nRd6dqBVdtX7rYvkwBzEMfrQe3lfyU39WM/M0GMkKtDql0/kJE4fUplsT9KbPrNMF3osF57byrz8HTLp3mjkYxwANog3fqZ7dO7oF0D6uwd8LCu9uO9B5ZBSkCJwq2K654NQXWhF6uBXMzjyCZ2FtNKF8qzqR3poSUK5FrpNdds0Yo0Zd5apj2yDEFe6SOPjv7HtU6gk4w96DAbiY32jfNdjtliR3sWv1CkpJK3a9QaFBusInIyXCEc3TJRN6bK0+yxxy2S4nhNWG8mfySt5CZZ7Buf2LQ/g/0t0bsvVs52eV3aaR0/lSM12tdC91Qb3cg6A3qdLAsTJFCK5TJcUKacswVgDXUnpfyqY+Uo3YA5KK6Ccz+dnK400Xn5PE5jtSfnZMZojYepI/WdfWmW08qzJ/Tul8Umoc0kNnCZ4u5K09NQ7x8d7avA/6f7A3mk3SgYrvhLkCJzSJHT7bry9ThBRqAS1iXRB6NsyZua1PevgFIKiurmfgNsD0wKEI5C7U7XH6VNL6A39xXZywcLRl3uPfczpnVzPrYnSGj/vt3uQVkaiJZMpRMuN0SC5scYv+Gyd4/VMneFOJarO3WUr9+v8y8tLE7E1cSnMDav/XL2UW3Zd+GTSfykrkSQFXW4k7KfzYS//dDoNtxLoBWk83j8TkNKxzDtTL98jk4K41m/uOtNAuLtvldsz1hTVp0mRXj3Dnbo6lPrikzAZT688TggGTS2EH3O9/u7pxFRGokfuFRLYjOkw9SF5hSSJsvRkchbC34lUQkDfzOMGzt/GBkbH0oWmyNxdmwzVt09ytcwDimSgSQdjA0m+3bp8RwjIa9l822TSSWvwE8pxn0XrXOmu6u2h5O6pfE+FXx/XagMnpPqcyvfzEyCNyXflgVtImFCOx3sjvQPp/9YkagWHbAMNqIfEXl6PKrTanupt9Y0Tj4ObykzoOWVge1yYk0m/qS28dIl1tfSR2uheupV9kBB1dbyn1VoymUJM4WOOiJTkcseU5AjsdbQMuOocEYmMPLaXmah6f1EkJ4D3hQjjm4O8LT13mlJa3jfXkX8mAvdoUg/VUxunTJx0bOLDrXMu5cKTDNgDMv7HAWXgM0bMe8jQ4s/MdrgH8EF3Q0r/+++eNmNbST/YXhkIt2P6tk2nMN5nGimP5QNp6DWSflJMUJWfptQ6l/0UEqZMQbKTlqHLKKdMXI0zdT5m0N+NwKngFFz1xI0h1ddqVwvlfcCnoIxH6q22Uqlqc9pTQJ6djqlqPwzQ0dvL9L5zSZjBtLfZ+zeFJqkNaXzGxHU+5g143CX1XLQS5xX5LCGX8rST04UWd0rLwh7E9drlwYu3ly+iobdJZSk3ymezEVVUn1+yr808ulGKjqEhept1rlCY90PztkpowL7vNx1LfTIo6xnY91CH76UekqlG3DXgRPIp7v6kFF3FwOOqztfFU9Aio+AigHcoobwjTY/D0U1q4FSW0wtufZRhabFDLlYGxu2U/hxmXfjJm+5xoInQMnxdWu9t4rnUC/eiMHvgCM09GvQQ3K/hOqR7T0aewGFqitLFvcuHKMQUYaeTkSUj8oFw1CNat6jKWJum/NejDj6Y5L80bO6p9QtDhaj+me2gETlndgR2I+FE1uJD1dte5qp6txI4ppPk+i+lhGP/j5WX9oa05OKFDDaGg0zEoht3l55yoxwQrWA5+tJ1sw1CxV6esYXYrnDRBa6c8U42Hpxyn/F+v6/Hx8c6nd5uL0A6obamW5KS9/KWq+m6OTfgNGakdrz7jvQrBHbKkXTmVJAo7ONJmpO8BNuXxA33b0PgxGOTy5DytYXQAwrvuBsz3V3V+XVWaXWGXJavlvIoKapXWbqJHgO2dOWH7gXx+vpa2+12eCazrv7wIC4hXSvVpv/8PPab//GyJV6Y6Tpgf5Ixp438pJI+nmy3Zzx6+k+BqE+YnteQfDcCU1oclP+3uJOfn3ZIjdXNMhKyOspw4nh9zpOlSN2ohLyXBsw7XDLo40RlW5yeuOe4RAdst/rl/abxkd4RnT3mSQG065iTjmOYxmhKZmlqJkryXVmI1EE1lHwwcWD+xs9COyIoEbDFUWmMNJak/BSMUuHMf+o/CQdQvzvX7dGWpFMdT1eeBrClA+nBzyNyJcRk39V/cmiiINE8UTE3Xteb665Fe7zf7o1bMhmBWzDv7lcGfDgc3vGd3pY8Bm+OkKzfB0zHOmck51PmoKUMH1RNJJbXSmOlrIWifdcd+8z9x1Xvb1DiE8+Na0zIVcmvnWOrr+y76uD9ltPiDg3XJ0fVm2E6KrdswunKGAhUXXhRp1/ImMRnpbvC5F4SfWA5iYPz3cuTXOqOqs55Jzk1DdgXDCTk7CqLG3WIUum+Z260bsBTVjgp1Kc+p8no1MPb47/xPO8DU24Em54Rf0QuvqxeM02dY+cdVVzSPmHnbuqou2s/3hXASL7qPO859qwxz/nSkyg4kvJ1NyBeqezn+V6EVubFpaU7R12nQC4+XrqUXu1VvpmXPXHC+JXYPeQlehNwOA7qm8bWPZ4b9iXAM9mA3R2pcWoAjdeDNjVKqEPXzgFnp+TayJ10vAeLREtu6aQBn06nePWGBkH94sMFlQrTpfS6meD9/f3ZXXu8XeLFjuLJzaZBUp/Tqp1TL34mnXL9SKdawND4UaQvGq/f9cg/0yaYIvQJpzFJmZUWN066cbkIgadKMjIZYHIxrMMDQj/GleiUhdKbySkYcWP27AF5tQxYd8dkGxQUcSLRIFPg5e2SwbnBpwmgPsjwmQJk21l3i0c7fWH5fNdn6r8lpIdT6U+vjZSLL+psNUAD5MEDDZj0o6oGNEjZAw+M6F7dG9DdO5orgCMCu2HTgBNP1c39bm5uarPZDI/UXS6X7xD4dDqdrb6xTbqBifTl+vPfZYzkrS0P1JrE1B2N0wNNesWxgI0eymmTyvLx7EmaMPy9J5Nu7ucrJB5g+Dk0PA0ggwjPRrDzrIs7+5MS6J65KsZUHHmfB0haYaOh6bsGdDabDQ/2Xi6XtdlshjtV6hnRmrQyIt4fo7X/wb1J8liObqQlnMTSBc9ldsMDTjcYT5k552Y7OL787tSyZ7Rjkyy1sSXffUVGT3zGUvG8YNAHsTXrOTD6j8f0sh5qr6Ou/9ZygY5eU/RCukT33dKVxwtjLzd2fudETW3j516GIbUljVmL0lB6KMog79Jzqy581CzvgtObHXQz2gtARNH6vG/nIzdUEJc6SSrRQlDn3lweFUrqc899sc+qR6tvogRErcTbW8g5ZohczWO/iOhp/4FH9hTnx9Sn2kHvqXcGauL31HPiwa0FkNSmpP8poHHRkzrpalqFUwnH49vDT4hEygzw4ScMWsiDvX6+05WNGTBfjLJ7gQUHQOXp2XC8z4R4tvdHbU8ZibS7zM/logNpWcuAh0FtpA09N9uiZeLGNFQacWpTy8OorjH9ts77Wwx4SurHf0/IytUfHZPeeb7K1HcaNnlZQjFvU5rVHtx4XxI6JHpU9X6BI9GbxBl71IAGwz5Lnzq+RQuSPtm/RAek5xZ94GRzXUuYW07SolVjlMKla8C64zo7kxrSq0xuR6LB0POGdQyRQcepQ5rp7pIc0ZJrU5nu5nlbVqGnB4AeVLA/DGZms7eb+Knt9EQqkxvdxzhlQj0atY4ZBhI7yjwQ9kUaCYNPUj7SFAbJRF5vF3XFIJDZDt8p6OLjO8WIuwaszde+RCpJPEriLlGdJA9O/E2/MXo+Ho9DJoEdoyv2KNkRXAringIirG7INwXJWKfOUX9S+8gf+ZhaR8TEl91Y+Zu3M3mMHmLqXLX96urqLH+tOIEg5PplnOKTZwxNxzziFJnMgd1dUhEuLWpACkBUTRTF3VWiBM4hW8fzvEQnnFun7AQpjPriqaa0lJ7QypGX/XcDcQPmxLxkcaDn7tkH9pMrockDjEnS4yXiueokXQPe7XZDQdrkLVc5lt5wzkSXquBOD00hd9SLv7lrcXTXZCAySlQ+b9KnPtGAhcDsl77P5/Oz+/qSfvixqlPtZJt4J3jmhyU8nvQj6VaiOqkTz7M7t+6NGYNggo1PPtZJPZOusV8tD+si3WnZXo/pbckkDqxb5zufmoLA/MxVNe4ZZqqKPEkI4GjHehzhnD5Une/Z8NSZK54Iq7rJeRnJO+r7nR3ZRhoky1RdErpqLgJJfHK2JNEG10sS1qexca/BMtyIEwXoecVekKeNUz2ZRCGkbB+wqeK0wylEMgoOLl1lKqtFI6SIFnXQd+9vr+0tGsDPRODEa1lmyxgkieO2jDDRiZT35fEp4KXRT0mTpf8voTYp+zOVC3cNWA+vvrm5OVuqVcd7s4fvaij50+FwGDZMV73tO6ALpQtP/NTpgwc+HBQPSojGFDdYfmYwlm501wrKfLmb/ezVm4Jk5+zuqtUOBas0AsYcrmtHSXkpH+MWhWj9nvTiwgnmd0f6EAeWO3PDneKKkniAwBU5rmKpbFKOltJ8TT5x4BZVaM1w8j+163g8f0SAI1OLa6qf0iMXIdzNJ2H7Wku/Y3sGWn1ku70NY0gvIZBJR/y91S4idKuesRu1VE3kwH6Hc59hKZr3YMkVQiR7fX09ewSUOujHT3EpiZ9OoT5TjIn9Tee3jJnvdOmesvPzp2QZ/Ho6D6ISF3We3+qH+sDyk7HTiNV23fyFfU468zZ6wP0hBH58fKyqv5S03++HJWDlDKVcHwgdp2iS2yZ1rJZk9T/dGfmvymDwxRlPBEkozN8kjr4tN5rcdhqUhMqsn9kGucUWsvYCpPS7/0/K5RmexP3d+DyT48u97KcL6Zm3jffW6PVDtsMLBnoyiQMTgVtIRneReFqa/cwyuPHqGJ1D4/cgyc9NwdWUIKRnGJ6l8PNaLjkhT0o1UY/Oh9mGFvVg29x4vR5HeO+L/kvl+7ZN6mhMxy7Jy8jYeWVITyZxYOUvffXMlcTGpM/+XAYOfmtWuyG2ZjiPTwasSTKfv13+lBTKOjzrkQxEfWS/vH1pENKl9N52Tlp6nRQks97WtYdJEiozVpER89IiorMDFdE/6aLXBg+2W5OcMmkhY7vd1na7rfl8XqvV6swd+TVUVecbzanMxOdkvDKqlAAn/XDU6JWrsp1GtAZMbU8TJQWCNDRmNRjM3NzcvBukFoVh/SrXRW6a7dI720buKP0kL6bPvrlfaKv/ucxOYRCuK6cd+dPkSX2Qrdzc3Jw9KaonkxBYkTQvVqQS3L16lJ+Ci+T+e7OtRwN6bt3pBQcuHd+S1B/3MDqO/UhoxN969Cfx3lZw53u1HbHTJHcEd32kvvvx7oFp0ET+FPxTmDZzFO/JpA3tbrSaHbqokZUK9dJgpbKdM0oBHiHTrSWEkhJpCOTWPG5MepzM3TN5PK+6Zl1CETdu9YX98snG9pCGccL7BJkyoccmLMdHOmAf9BtRl4Yn8aA0eRUaPV9sR0u6Bux7UNUZRYe6Tkw3vpMBq6FKkzlCuQFySVnZC7k0DY4HdBJ3oWw7J4d7DQ90ktBgeQWyDJiTRfde0ID4BGZbJaQNieJ4gMV3vyCzxzkTnWIOXeLekXED9cXnSUsncvtpErUCWp8IHry1DJ4y+Zq4ZABquF+hkch8a/D0OeUkPRImcrg7a/FKz8cmDtiTXr8cbWnAngdP+uQ5yciTzlptTpPRA9kkYwjXqsdRcwwMeuOg9xS4fQiBJZptmnG6Mne9Xtd6vR5QU7NGjfCbPp9OpyEwdPfrKTWV58EIjZ3G5RyTwU66CiShEfe2iiJxV9RyuXyH3h5Usi2twXQ07ckY/XEq5v95brz38sux2BchI4Mroq7Ahtzb25XKla6JvC1PkmSyAWtn0Hq9rru7u7q/v6/NZlOfPn0aEJj3UdDihBTDW+5T0dxzqo4rivclZrpCKkCGzHZw11dSIl2qu3AN2HK5HIxXl9I7TeHtp1qegPWmPbYaUNbvxjhlMqRFEBowl7W5f8T3kjhKEsB4KwFNcqcL7Fcrx8/HEuseG5wQKvNvoxA+Y/yh2IyCX15earFYDJfQ80pg50hO1lP067/ru5/jBjQWPPLFIJKBC1FHKTGKZzcSFXC6kORSV57Od1pGHuv68u+ty32cMvYCNkprRdHLrap3FHQqvaua+LR65nkVzOiuNLpnmKdxmFecz99uEC102e12w14LX30TSvjslhF4EOfHMNAiVWA6kDei5jEMKoS+fFH4OAIiF9tI983v7q7TZEp8n+JI55douRHzukHuLuR+Z9VHGqQJ7EbL9jGTNBZ/0JaoawKiZyNaMvlZyUQhDSY5sBpSVUOGQnsmhMZq9MvLSz0+PtZ2u63j8e1Khfl8/m4HWgqSWgZMxfus1kDpng66XwVpAKmI6JLow2azeWfAz8/PZ9eMMU/uvNN/kwFrgBNPpDFyPPw/GaeXWfX+gTQtA+YVJ9KhtrvSgBnvpPITXWBMw8yDDHixWAzpWNXV8qAuFyFwCmBS7o4GJMSezWZDQ/U7d7gxMJsSNScOyQFoReNV5xdGOnJ6X+l9fFXIc7Gsz4Mm/e6omwIv/0zES/W0kDt9J/q36AX1yM89NGz1KY1lK6vj9jWFRnQN+P7+vqqq7u7uBoj3dWpPK7GzMtTVanV2zdnLy8uQSz4ej7Xdbs/uQiNDkgJoZPotfU4Ziao39/r8/Dw8x0I34BMi6TyhgnZDyevQsJnU9wieHsVRiS9yzkQ9XNJgMpfLAK1lwNyQ1dpc7xyeMYpAxp98pP9bl5vROBUILxaLMwqhz60FjZZ0Dfjr169V9ZcBK0r0ZHMi8xrk4/E4IK46J+RdLpfDHW4eHx/PHpYiuqFB6gUfek9owYEkbdDDWUQhyH1luKvVavhMFKYblQEru0JXLjrh3NT74f0hDx2TjxqwI3HPi0mUcUmez70NbaGqzjguV3GlZwbOY0GipGvA/rC8nhtxVGQHnJuqbA0sH8slVJgaiSbFe708tpXW4bm9unrHJkNs6SUZsLv01C+vJ02CFn3wgLJF01poz0lLO0jBdau8BHqX5n4pXQO+vb2tqhoCtRSBJoThQDCwEHopcpcyZrNZLZfLIe0mOtHid3z3CVJ1PuHSbO7REUdsiR4jq5ylJht5rtDX0S0NbIraq+qdLpO4Uaqc1uVKqY0cI6dFLWNivjttSU0TggGggnvm2OnZnA9PMequAW82m6qqwZ2qIa5AD7x8S6D+Y3Cn9JtmsziZUm5yy5wAUnALQZwDS+GJS7F9/E0DI7pR9deEkAErMlf5XAzhdzduCRcpOMklLYNnGzmhW+krjoUbsB+TMjnM3hBovG0ea+g3GqBog9Kv4sAppmp5nZZMvrmfBsV5kwwwuUu+s1FE4FZGIyEsfyPNaKVcev9J0gynYQlRec8HXsvXcuHJmHzytVx/QjKfcD2q4PrnOCbpGfBUt96iY55ZIIjRYFtj9CEO/Mcff1RVDQHKcrmsb9++1dPT0zB71ut17EwKVHyQ6MLUMd4aNBk6B5dbG6kYHsNr8jjTU1qs6vzyqaq/qMNisajX19dhGVWf2Qedo3dfQFF/qCNP5fWoUis1lX73sp1atGiDPlO/LQNzw+Z3UkUudumGjsox87rJlJ79MAJ/+/atqt5WnA6Hw7AAoYZxpz45ZgtpXNk6T8YpAxblSIplvpmpLTfgVlTr6KCMCT2LKMB8Pq/9fj+0TbrwVULySjfgpIOU9Hc6kQy15dloyE4VEuVqGW3S+dQgS8dzn7AyOqINMmDfR9yKV8ZkEoUQssh9pos8W1G5u8nkxqY0lkplnQkhkjtMwYHXS36q/pOjV70htAcZ7vr9MyN3Grrryr/3KIGLxx3pv5b4gowb75i4/ukZ6R0TpUuT428xYA3afr+vh4eHur6+rm/fvg2ftRDRQ1enD/rNj0/Cc7iwoXcicYvrJn5NxVadewDWoSyEvjMwFIIoz628L9st9+3GWHW+4T712WnFmAE7ukucg+t9DHVdr2MiKlBVg3e+uroaaIP2O6R871gs05NJV2Sw04+Pj/X09FTL5fLd1kiXZMTe6Z7x6p1G5/U5j0sI4osvdFXixh7907voswyeK0euq8RDHWlZvv/ubdB/6XNLemCRuGvSGXXq53lWRe/czsoVN9IGR2b3jqz3wxzY3baiaKWYtKrlO+kT9+Ogtgza6+bLN81PaXcqx3/roQxdP9NI4piiVLPZ7B3nlTh6jrWZkzRN2LG+p/7Tg0ypP6Exy1K7eKzHIp7bTRRiCpUYk64BM8OgwTkcDvXnn3/W6fTXxvUvX77U8/Pzu4ecJIRpBXdp4KUUBQJySZ7C87veuGLULgVhvMzd03WuPAZ0+u/l5aV2u92AKLpvMiem74Xg6qK3NYn2T9NYpDPqk5KQ0j0c6/QgOh2TeLHK9FSYsgzcyTebzYadi9IXr3lMwbjLGAp3DViDTYN5fn6u7XY7JPf3+/1gZEIpReE+aG6kUoYbNZXFTeWqw9fzW8ZBVPBFjbQ8rjJojEqZsW/a26zJ5Er2bERy+bwAkue10HOKJG+SYgLV43rzbEMPiXmO+qNFHm7U8cuOnEKk8vXbhymEz2Z1WhRit9vVdrutqrfnLLgBU1qzvkUxdJzzJ7pyDrLTBP89KSVF+i1dJOrh2QW9t8pKQSRdfGtBwP9LhsnxSpy1F+lfkrqScExIHfzVS8MlWndJOyYhMO8m/vr6Wk9PT3U8Hmu5XNavv/46LDXzqUOtVSo1lgp35dF9SlHaDceofLvdDsu9LaWk3f10oxJu9GYA5pmKNBikGiqL/R+UjYQ9797ICUyvJKNl/tlRU0JkoxDVPYPh/HXMcPw4Bmvc56ALYOWZuXvPy2O5apN7z55M2tBOHng6nQaj2e129fj4WK+vr8M2SQU2rd35Q8UYzLRhPiEwDVjZgzSD9T0hcGuy6HNV+/4FrTJ0jl5+IamXRSMW1VDZXLUjBxYv5mf21RHY29bLRiTe2xOOG6/U4P7edKWx6mI9Xnf63m3LaGs/KJdwuB/yQy6V2Q8D+yH/ZPk/R+Af8kP+L+WHAf+Qf7T8MOAf8o+WHwb8Q/7R8sOAf8g/Wn4Y8A/5R8v/B5HhhKveyp7aAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "\n", "from sklearn.datasets import fetch_lfw_people\n", "from sklearn.model_selection import train_test_split\n", "\n", "lfw_people = fetch_lfw_people(min_faces_per_person=70, resize=.6)\n", "\n", "n_samples, h, w = lfw_people.images.shape\n", "\n", "images = lfw_people.data\n", "targets = lfw_people.target\n", "\n", "training_images, test_images, train_targets, test_targets = \\\n", "train_test_split(images,targets, test_size = .20, random_state=10)\n", "\n", "plt.imshow(training_images[0,:].reshape(h,w), cmap='gray')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 3.1. PCA decomposition\n", "\n", "- Extract the first 50, 100, 200 and 300 principal components from the image dataset and display the corresponding reconstruction for one of the faces. \n", "\n", "\n", "- Display the plot of the reconstruction error (deviation between the true image and its PCA reconstruction for K=1 to h*w)\n", "\n", "\n", "- Finally, plot each of the principal images for K=1,...10" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from sklearn.decomposition import PCA\n", "\n", "# put your code here \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 3.2. \n", "\n", "Now that we have a PCA decomposition, we can learn our Max Margin classifier. We will do this with the svc model of scikit learn. finding the optimal parameters is a little tricky because of the dimension. One approach is to carry out cross validation over a parameter grid. To spare the search for the parameter values, this step is implemented below. \n", "\n", "- Complete the code below with the training of the Max Margin classifier (remember to split your dataset between a training and validation part).\n", "\n", "\n", "- Once you have trained the SVC, apply it on the test set and display the predictions vs actual targets. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sklearn.svm import SVC\n", "from sklearn.model_selection import GridSearchCV\n", "\n", "param_grid = {'C': [1e3, 5e3, 1e4, 5e4, 1e5],\n", " 'gamma': [0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1], }\n", "\n", "clf = GridSearchCV(SVC(kernel='rbf', class_weight='balanced'), param_grid)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### II. Independent component analysis: The coktail party " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this exercise, we will use another approach to dimensionality reduction, known as independent component analysis (ICA). ICA is particularly useful in speech or more generally, source separation. In the classical version of this problem, known as \"The coktail party problem\", one is interested in recovering two distinct signals from their mixing. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Drawing\"\n", "\n", "image credit: [The conversation](https://en.wikipedia.org/wiki/The_Conversation) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercise II.4.a.__ Using the FastICA transform from scikit-learn, recover the two speeches from the mixed1.wav and mixed1.wav files which are given on github. \n", "\n", "(Hint: start by storing the two signals 'samples1' and 'samples1' into a single matrix, then pass this matrix as an input to the FastICA method of Scikit-learn) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(__doc__)\n", "\n", "import os\n", "import wave\n", "import pylab\n", "import matplotlib\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import signal\n", "from scipy.io import wavfile\n", "\n", "from sklearn.decomposition import FastICA, PCA\n", "\n", "###############################################################################\n", "\n", "\n", "# read data from wav files\n", "sample_rate1, samples1 = wavfile.read('mixed1.wav')\n", "sample_rate2, samples2 = wavfile.read('mixed2.wav')\n", "\n", "fig, axs = plt.subplots(2)\n", "axs[0].plot(samples1)\n", "axs[1].plot(samples2)\n", "plt.show()\n", "\n", "# add your code here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercise II.4.b.__ Once you have recovered the original signals, plot them as time series. Then use the lines below to store them in new .wav file. Then play them and compare them with their mixings." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# write data to wav files\n", "\n", "scaled1 = np.int16(S[:,0]/np.max(np.abs(S[:,0])) * 32767)\n", "wavfile.write('recovered-1DemoJune28.wav', sample_rate1, scaled1)\n", "\n", "scaled2 = np.int16(S[:,1]/np.max(np.abs(S[:,1])) * 32767)\n", "wavfile.write('recovered-2DemoJune28.wav', sample_rate2, scaled2)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }