{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "線形代数-固有値(LAEigen)\n", "
\n", "
\n", "
\n", "file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/laeigenvectors\n", "
\n", "https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python\n", "
\n", "cc by Shigeto R. Nishitani 2017 \n", "
\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 固有値\n", "\n", "$A$を対称正方行列，$x$をベクトルとしたときに，\n", "$$\n", "Ax = \\lambda x\n", "$$\n", "の解，$\\lambda$(lambda)を固有値，$x$を固有ベクトルという．$x$がゼロベクトルではない意味のある解は特性方程式det$(A-\\lambda E)=0$が成り立つときにのみ得られる．\n", "\n", "この特性方程式をpythonで特には，sympyってのを使わないといけない．一応示しておくと，\n", "\n", "$$\n", "\\begin{align}\n", " A\\, &= \\, \\left( \\begin {array}{cc} 3&2/3\\\\2/3&2\\end {array} \\right) \\\\\n", " \\det( A\\, - t E) &= \n", "\\, \\det \\left( \\begin {array}{cc} 3-t&2/3\\\\2/3& 2-t \\end {array} \\right) \\\\\n", "& = t^2 -5t + 6 - 4/9\n", "\\\\\n", "& = t^2 -5t + 50/9 \\\\\n", "t (\\lambda) &= 10/3,\\,5/3\n", "\\end{align}\n", "$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "⎡t 0⎤\n", "⎢ ⎥\n", "⎣0 t⎦\n", " 2 \n", "t - 5⋅t + 5.55555555555556\n" ] }, { "data": { "text/plain": [ "[1.66666666666667, 3.33333333333333]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sympy import *\n", "t = symbols('t')\n", "A = Matrix([[3, 2/3],[2/3,2]]\n", ")\n", "xx = diag(t,t)\n", "pprint(xx)\n", "eq1 = simplify((A-xx).det())\n", "pprint(eq1)\n", "solve(eq1,t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "固有値を求めるコマンドEigenvectorsを適用すると，固有値と固有ベクトルが求まる．ここで，固有ベクトルは行列の列(Column)ベクトルに入っている．" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([3.33333333+0.j, 1.66666667+0.j])\n", "array([[ 0.89442719, -0.4472136 ],\n", " [ 0.4472136 , 0.89442719]])\n", "array([0.89442719, 0.4472136 ])\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":4: DeprecationWarning: scipy.array is deprecated and will be removed in SciPy 2.0.0, use numpy.array instead\n", " A = scipy.array([[3, 2/3],[2/3,2]])\n" ] } ], "source": [ "import pprint\n", "import scipy.linalg # SciPy Linear Algebra Library\n", "\n", "A = scipy.array([[3, 2/3],[2/3,2]])\n", "l, V = scipy.linalg.eig(A)\n", "pprint.pprint(l)\n", "pprint.pprint(V)\n", "pprint.pprint(V[:,0]) # columnの取り出し方" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "固有ベクトルは規格化されている．" ] }, { "attachments": { "image.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAHKCAYAAADfFkj4AAAAAXNSR0IArs4c6QAAAERlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAAB4qADAAQAAAABAAABygAAAADZGuYSAABAAElEQVR4AeydB7gUxdKGW685K4pZMWDOOYsBBVFREEUUA6iYEwbEhAFERVEUEUURUTABRswC5pxzQExXEbNe4+/dv9/y9mHPORtmdyf07FQ9z57dM9PT4euZqQ5VX82Ss2JUqkJg6tSpZqWVVqrqWt8uOu+888wZZ5zhW7W0PopAZAicfvrpZsCAAZHlH1fGCy+8sJk8ebJZZ5114ipSywkZgVlDzk+zSykCZ555pjn//PNTWnuttiKQXQS+//57s+OOO5o33ngjuyCkvOWqiEPowBVWWMGwsJDWz3zzzWf+9a9/GZRx//79Q0BEs1AE0oMA93xan93bbrvNLLHEEmbGjBlmhx12UGWcntuuUU1VETeCI5v/zDnnnGbo0KGijFmeVmWczftAW50+BGabbTazySabmA4dOqgyTl/3NdRYFXEDFNn+sddee5lRo0apMs72baCtTyECrGaNGzfO7LLLLqqMU9h/VFkVcUo7Lopq77fffuaGG25oUMb1YMgSBU6apyLgGwKsao0fP76RMn7zzTd9q6bWpwgCqoiLAJPVw/vvv3+DMq4Xq9Ks9qW2O1sIqDJOb3+rIk5v30VWc1XGkUGrGSsCkSLglHH79u3N119/LQZcOjOOFPJQMldFHAqM9ZdJU2V8wQUX1F8jtUWKQB0igDKeMGGCUWWcns5VRZyevoq9pijjkSNHyp5x3759jSrj2LtAC1QEqkKg0Mz4rbfeqiovvSh6BFQRR49xqkvo3r17I2U8cODAVLdHK68IZAWBueaaSwy43Mx4++23N6qM/ex9VcR+9otXtcpXxqeddppRZexV92hlFIGiCDhl3K5dO9kzVmVcFKpET8yWaOlauDcIXHTRRWbuuecuWZ/ddtvN3HXXXQZlvMoqq5hOnTqVTK8nFQFFIHoE3nnnHdOvX7+SBW244Ybmww8/lA8EIJ9++qlp0aJFyWv0ZHwI6Iw4Pqy9LemUU04pq4Sp/Lrrrmu23HJLacdTTz3lbXu0YopAVhBYffXVTdeuXc0ss8xS8jP77LMbeALmn39+8+uvv5rp06dnBaJUtFNnxKnopmgriSIOKgsssIB58skngybXdIqAIhAhAquttpo5++yzA5dw++23m7fffjtwek0YDwI6I44HZy1FEVAEFAFFQBEoiIDOiAvCkv6DLFUVEw1BXQwZPa4I+INA/jOsz6w//RJFTVQRR4FqkzzzH6j8U/kPV5A0+dcG+Z2fv0tfrBx3Xr8VAUVgJgLFnpf8ZytImpk5BvtFnk3LyP8/WC6aKi0IqCKOoafcA9T04covmjSFzhc6ln9dJb/DzKuScjWtIpBWBHx5dou9H9KKq9a7MQKqiBvj4d1/TR9AlGkpcS8O910qrZ5TBBSB6BCo9tmlRvr8RtcvPuasitjHXilRJ31AS4CjpxQBjxGo5dktNwD3uNlatQAIqNV0AJDqIYkuS9dDL2obsoiAU8K1KPIs4pamNuuMOE29ZevqHspi1daHtRgyelwRSBaBap5dd40+18n2XdSlqyKOGuGQ86/mgXQPc8hV0ewUAUWgAgQqfXZ5biu9poLqaFKPENClaY86o1BVwnoY9YEuhK4eUwSiQ6CWZ7fQtRxTqU8EdEYcQ7/mP0D5v/OVozvuvl218tO4Y/qtCCgC8SCQ/zzm/85/Lt1x9+1qlp/GHQvy7fJx30Gu0TTpRkAVcQz9F+SBDJKm2qpGmXe1ddLrFIE0IBDk2QmSppK2hp1fJWVr2mQQ0KXpZHDXUhUBRUARUAQUAUFAFbHeCIqAIqAIKAKKQIIIqCJOEHwtWhFQBBQBRUARUEWs94AioAgoAoqAIpAgAqqIEwRfi1YEFAFFQBFQBFQR6z2gCCgCioAioAgkiIAq4gTB16IVAUVAEVAEFAFVxHoPKAKKgCKgCCgCCSKgijhB8LVoRUARUAQUAUVAFbHeA4qAIqAIKAKKQIIIqCJOEHwtWhFQBBQBRUARUEWs94AioAgoAoqAIpAgAqqIEwRfi1YEFAFFQBFQBDT6kt4DioAikDoEPv/8c3PDDTeYxRdf3LRs2bLhw//zzTdf6tqjFc42AqqIs93/Vbd+4sSJpk+fPmaxxRarOg+9UBGoFoEll1zSXHPNNeazzz5rlsU888wjijlfSef/5p599913zc0339zs2qwcGD9+vFljjTWy0lzv26mK2Psu8rOCvMh4Ga699tpmo402MksttZT87745t8QSS5jZZtNbzM8eTHet/vWvf5kePXqYc845p1lDfv31VzNt2jT5NDuZ8QN///23IHDmmWeaq666ynTp0kWe4ZVXXtnw4fmddVbdsYz7NtG3ZNyI10l5Cy20kPnhhx/Mq6++Kp9CzeKBZvaRr5zdb7633XZbQz4qikA1CKCIzz//fOOUSyV5cG9uvPHG5rnnnqvkstSnZQCDzDLLLObLL780Q4YMadSmueee26y44oqilJ1yXmmlleT/5ZZbzrjrG12k/9SMgCrimiHMZgbdu3c3o0ePFmVcDIH//ve/Zvr06fJ55ZVXGiVbffXVTdNjjRLoP4pAGQRQDDvvvLNhm6QSYZXm/vvvN5MmTcqcInY4HXXUUebKK690/zZ8//bbb+att96ST8PB//1ACTvlzEz6wAMPbJpE/68SAV2DqBK4rF82++yzm2uvvbYqGBiNX3311WbOOees6nq9SBFwCBxyyCHuZ6Bv7r1Ro0aZHXfcMVD6ek3EasLyyy9fUfNYeXjvvffMI488YjbYYIOKrtXEpRFQRVwaHz1bAoG99trL7LTTTiVSFD7Vs2dPs8022xQ+qUcVgQoQ2HXXXcUWIeglF110kenWrVvQ5HWbjkFw//79q2rfueeeK/vKVV2sFxVEQBVxQVj0YFAErrjiiopmtliv8jJUUQTCQICVmbZt2wbK6oQTTjAnnXRSoLRZSLTvvvuaDTfcsKKmbrXVVqZ3794VXaOJyyOgirg8RpqiBAKrrLJKRS+3yy67zCy88MIlctRTikAwBFgm7dy5s7npppvKXrDPPvuYQYMGlU2XpQQYrF188cWBm4x/Nr7barAVGLLACVURB4ZKExZDoG/fvqZVq1bFTjccb9eunenatWvD//pDEagGge+++84cc8wxZq211jL4w84xxxxiRFQsr+233172hdUtpzlC2223nWF5P4gwkMGCWiV8BFQR14Aphh+I+64hq1RfCoECM91ygm/xX3/9VS6ZnlcECiKAFT5bIVhLY/GL8RDW+7BsnXfeeQWvWXfddUVZq2FgQXjk4IUXXljW33/zzTc3hx12WPFM9ExNCKgirgk+vdgh0LFjR9OhQwf3b7NvZiMjR440W2+9tfnqq6+andcDikApBB5//HGz3nrrmWOPPdb85z//ESPBN954w9x4441m0UUXNXvuuad85+eBVTCuTQsuuGD+Yf3dBAEYtrCiLiW4Gt5+++2lkui5GhBQRVwDeHppYwQgB4AQoKnwAn3iiSfMsssuK36bGIg8++yzTZPp/4pAMwQgjcHKGfIXFC/Gfrfccot58MEHzZprrtmQnhnvAQcc0PB/ixYtzAMPPCBkMg0H9UdRBPr161eUo3uHHXYwv//+u2wrQaCiEj4CqojDxzSzOcLIc+qppzZqP4YdcAJvscUW5p133jFt2rQx//73v+XFOnbs2EZp9R9FIB8BFC73FPcJAzzoLD/55BOD4VUhcT7FpL377rvNaqutViiZHiuAANtGhSzKsenAb/jSSy+VLTioMY888siq2MwKFKuH/oeAKmK9FUJFAEUM+44TGHygEkTmnXde8/DDD4uhzZ9//ikzHVwhqqEodPnrd/0h8OmnnxoM+3Cv+f7778U9idnwWWedVdJVDrY2Zs4obgZ+KpUhwLOIQnYCDe3QoUPlX1y/7rjjDhkQDRs2zOyxxx7m//7v/1xS/a4RAVXENQKolzdGYK655mrgr2UpuulSFvSCLGGPGTNGDEQYaUNT+MsvvzTOSP/LJAIEIlh11VVl6ZmZ7SWXXGIeeuihwNa6EyZMMNgrqFSOAO5JLogGBqjXXXedWWSRRRoyYh+egTTH7r33XmEn+/HHHxvO64/qEVBFXD12emURBNq3by/GM1i4zj///AVTMdt59NFHxcCGb4y4WLJWySYCBCDA2I8VFPYjYW2bZiMonXjiiRUBoj7qFcHVLDFGW+y99+rVS1YlmibYcsstzZNPPmmWXnppM2XKFINV+owZM5om0/8rREAVcYWAafJgCOD4X25mAs0lRluESySK0yabbGLefvvtYAVoqrpB4M477zQrrLCCWDgvsMACYoyFhW7Lli3rpo1paQg2HfDADx48uGiV2QJ46qmnDGQ+7NmzHaCD6KJwBTqhijgQTJqoUgR4oQYRCALefPNNmRF/8cUXhhE3rioq9Y8AkX4w/GHJ848//pAZGAZ9xYyx6h8RP1oIjSVbTKUE1zBmxrg+0WcMoj/77LNSl+i5EgioIi4Bjp6KBwFcTQhJ16lTJwmryJ4xe30q9YsAy864sWH4g+vR5ZdfLjNiDIRU0oEAscZxS8QYk0E0K1z0q0rlCKgirhwzvSICBFgSYzny6KOPlj1C4p1iLKJSfwhg1bz22mvLTKp169bmtddeE6KOrDPUpbGnMdzCvWmdddYRJcwyNcvVKpUhoIq4Mrw0dYQIwL6FgRd0hbg0HXrooSX3qiKsimYdAQLQm0KTCEEHVvL77befefnll8VKOoLiNMuYEGAbipkxLmO4nsFfrcvUlYGvirgyvDR1DAicccYZopCZIWE126dPnxhK1SKiRABaUxiarr32WsOAi5i2RE3CZUYl/QigjCdPnmw222wz8/HHH0tfYwmvEgwBVcTBcNJUMSPAEvXNN98sL21I6U8//fSYa6DFhYUAPMUbbbSRzJqWWWYZ8/TTTxsYmlTqCwFiQ99///2y9//BBx/IzBiKUpXyCKgiLo+RpkgIAej12E/kAR8wYIBGf0moH2opln1/rHAx5mHpEve0TTfdtJYs9VqPEVhooYXE6wFfZOJFY3j5888/e1xjP6qmitiPftBaFEFg7733NrfddpvMjFnW1GXqIkB5eJjBE65Iv/76q+nZs6dYxhcjePGw+lqlKhEgLCoMXPCEP//880LUonSYpcFURVwaHz3rAQLw2o4bN05mxixTs4es4i8CvHQxtHPbCRjfjRgxwswxxxz+VlprFioCcFZjwAVZD9/EjSaetEphBFQRF8ZFj3qGAMqYiDosU/fv31+Wqj2rolbHIsDsd5dddhHFy8xo/PjxOnDK6J2BTzihKIkHTSStk08+OaNIlG+2KuLyGGkKTxAgIs/o0aMNPsfMtkaNGuVJzbQaIPDtt9+KtSzLkgRsgKSFAZRKdhGAi5rBGKQtBHi58sorswtGiZarIi4Bjp7yDwH2HNkrxrXp4IMPNjfeeKN/lcxgjaZPny7GWHCHwxsNSQe0hyqKwPbbb29Gjhwpz+zxxx8vkZsUlcYIqCJujIf+lwIEUMAQ0+dyOTECuueee1JQ6/qt4tSpU0UJv//++2attdYSkg4Ys1QUAYcA0db69esnRD38hl9eZSYCqohnYqG/UoQADE34omIYxCwZ31SV+BHALYkQlihj3JII2IELi4oi0BSBs846y+CSCKvabrvtZr7++uumSTL7vyrizHZ9+hsOOxOzY6L47L777ubdd99Nf6NS1ALwhvCfEHj4CLMnrPGAU9SBCVSVJerNN99ceKl5ZqGyVTFGFbHeBalGgP1iHmgMhTDmUiafeLrTKWFoDKGudAZa8ZSupaQVAcIrTpw40cCw9txzz0mwj7S2Jcx6qyIOE03NK3YEXNQmQuoR9QVlzAxZJToEPvzwQ1G+LDG2adNGjG9wVVJRBIIgwNYFltQo5auuukoIe4JcV89pVBHXc+9mpG0QRTDKxlqXUTZRfVSiQYCoOltuuaUsRxNlB9zLBZGPpiaaa5oRYEuDWNTIgQceaOAjz7KoIs5y79dR21u2bClKgdH2hAkTDIYhKuEigHENy9B8wx+NtTr+wiqKQDUIHHTQQebwww9viD+e5W0lVcTV3EH/uwb3GRV/EFhttdVkmWu22WYz559/vrn11lv9qVzKa/Ljjz/Ksj9RddZee21z3333mXnnnTflrdLqJ43AZZddJpG5PvroI/F+SLo+SZWvijgp5LXcSBBo27atueSSS8THGIvqN954I5JyspTpH3/8YTp27CjLh6uuuqp55JFHDPFnVRSBWhGAcYsIXaysPPTQQ5ll3lJFXOudpNd7h8Cxxx5r9t9/fzHaQoFgUa1SHQIQ9bOEOGXKFHFN4mXJNoCKIhAWAq1atTI33XSTZHfiiSeat956K6ysU5OPKuLUdJVWtBIEiPYDwcTHH39s9tprL5khV3K9pv0HAcJOQtgPcT9+wsstt5xCUwQBrPV1u6oIOGUOd+rUyRx55JHmr7/+Ml26dJHgIWUuqavTqojrqju1MQ4Blrzuvfdemb1NnjzZnH322e6UfgdEgMHMxRdfLLGgcTeBwF+lOAIssc4666zqG1scopJnBg0aZNZcc03zzjvvZC9Skx3BqVSJgDUwwForkx+7hFQlavFe9uijj+asr3HOviBzDz74YLyFp7i0p556KmeN3uTevu6661LckuJV79u3byafXd5Zdvm3ODAJnnn11VdzdhAt/WJDKCZYk3iL1hlxyTGankw7AkR+gQqTvU78i+FGVimNALzRsJXB433KKaeYHj16lL5AzyoCISHAqgseDwj33XfffRdSzp5nE6/e19LCQGCxxRaTEaN9YYaRXSbyaN++vWC27bbb5qyCyUSbq2nkzz//nLPuSYLVrrvumrNcwNVkk8lr3D129913l23/gAEDBGNWaqxfdtn0WUrAPWcDiQg+3bp1y0TTdUbs+UCpafUwPpoxY4bsfcIkpRIMgVGjRpmlllpKrH/diDvYldlKRYg6XL5WX311c/PNN8ueZ7YQqL61rCQgQZ7Lt99+W9KyUsNevMpMBNhnJ844fupjxowxdmAz82Sd/lJFnLKOff7556XGGnS9so6zqwjyUM8yyyyy9KVhE5vjd+mll4qBGz6dvPzUV7g5RsWOoFCnTZtmuL+CKOLXXnutIavzzjvPQGihMhMBXJouuOACOXDEEUfUfTAXVcQz+z4Vv1588UWpJ1ytKpUhYJelzamnnip7n/gZ//TTT5VlUMepGeCxH4wiYRay8sor13Frw28aHNwQnyyxxBJlGcdIlx+yE7ennXbayRDJSmUmAkcddZRQqRJm8+STT555og5/qSJOWae+8MILUuONNtooZTX3o7oYboEdS/zHH3+8H5VKuBbffPON6dy5s8SGRRnvscceCdcofcW/9957UmmYx8oJy9L4y+YLy9oaxjMfERuj1y5RX3/99QbKWmu5b5544onGCeroP1XEKepMlr9efvllqTFkFSqVIzD77LObsWPHSsQggpRnYf+pFErWEsZABfr555/L7EP3z0uhVfxcJYo4f1k6P8fXX3/d7LbbbhrGMw+U1q1bmzPOOEOIUnr16mX+/PPPvLP181MVcYr68v333zfWqlXYjVq0aJGimvtVVZZdBw4cKJViiZoZYVZlyJAhsi/M/cSSNLMPlcoRcMZXBB4pJ8UUMdc9+eSTZu+99242Yy6XZz2fh90NXCH6gEe+HkUVcYp69aWXXpLabrjhhimqtZ9VhY96s802k4HNcccd52clI64V1tG85BBoLJdddtmIS6zf7N98801p3FprrVW2kaUUMRfDCNezZ0+ly/wfkrDkudjF/fv3N59++mlZjNOWQBVxinrMLUurIq690zBKwj0nSy4S+ag5gpPff//dsOS344475p/W3xUiEKYipujRo0cbAiCo/INAmzZtJEzif/7zH9O7d++6g0UVcYq69JVXXpHarr/++imqtb9VXXHFFY0lVpAKEqD8l19+8beyIdfstNNOE3/hlVZayQwePDjk3LOV3SeffCLuNbjIYTVdSpjNBWWLIlYvM0CVfxBgVjzffPOZO+64w8AfX0+iijhFvekU8QYbbJCiWvtdVZaoN998c3EdwbUpC/LMM8/IXpvl4JZVAfyGVapHoJLnstyydNNaYKg0fPjwpocz+f/CCy9sGEAieDxYBq66wUEVcUq6EnebH374QUbc5UbdKWmSN9W85pprzBxzzGGuvvrqBqt0byoXckWwOu3evXuDq5Ja39cOsNsyCrJSVakipnb40zILVDGyXA/ZBzji2lQvooo4JT1po5JITddbb72U1Dg91cTABsIA9k0hmifYQb1Kv379hMWJcHPqqhROLzuSnSArVZUq4uWXX94cffTRZZe8w2mJ/7nMNddc5sILL5SKEtqUPeN6EFXEKelFVcTRdtTpp59u2DPmRTl06NBoC0sodwyKXHxh+I0hTFCpDQH8sB3tbJDVhXKKmO2CnXfe2bA/jLsitJn83mqrrWqraB1d3aVLFwPWMJHVjTtTJkJb1EEjO3bsKNFILBlFHbTGzybcc889grHlWM79+OOPflayyloR0ca+vKR91l2rylz0sqYIfPDBB4Kp3S5qeqrZ/0S2ItqS1YmNPnbWm9tzzz3lmGXmanadHmiOgOWKF7zmn3/+uohepUPilIwUYd1BiNepEg0CNuyfxOGFg7refIvZB3/uuefMMssso0vSId4+LngIPunlBL9ttj9gd9tuu+1kdYJVCma9t956q1gEMwv+6quvymWV+fMYWPK8QnDkgkOkGRRVxCnoPW42HlasW6F8U4kOASIQQSBAGDa35BhdafHkbOPdNvikXn755fLCj6fk+i8FJiwkyNIxCnjcuHESbOSxxx4zJ510kmGvHuHclltuKSQeU6ZMkWP6pzQC2DiwvcJWEoEh0iyqiFPQe4yk7eKMWWONNZSCMOL+wq/2hBNOkJkLRjLgnnbp27ev8BdDitCpU6e0N8er+j/11FNSnyCKmGAj4I/BUSGhfxBVxAJD2T+sDkIHiidA2v2tVRGX7e7kE6CIkSD0ecnXNv01wHALcgYiXaXdbQSLXoJb4J6l/qjh3pvTp08X/mNWUIJYTJcrnTCdiCrickjNPH/WWWcZDNyIzkQoyrSKKuIU9Jyjz1t77bVTUNv0VxH2HufaA81gmiO+EFSdfUkIEFZZZZX0d45HLYDdiRUTZrIsLdcqzJjnmWceQwCJLAciqQTH1VdfXWbFxHh2gVwqud6XtKqIfemJEvVwilhnxCVACvnUIYccYogtS3hARzgfchGRZ8d+JDPiRRZZRELJRV5gxgqYNGmStBjDqzAEZY4RElLPsXfDwCo/j3POOcfAHQ/BBy5NaRRVxCnotbfeektq6Qw7UlDl1FcRI5BBgwZJO9h/wmAuTcIs3lF2wqdt3TzSVP1U1PXhhx+Wem6//fah1XfrrbeWvB5//PHQ8qz3jDBgZe+dACZp5U1XRez5XYrF64wZM8yCCy5oll56ac9rW1/Vwz0CIxzwTxtxAO5KH330kczqmd2rhIuA9R82U6dONYsuuqgJMxqaU8TOGjvcWtdvbo6DmtWrNAZvUUXs+b3pAo5jMc3yi0q8CODOhDDSDho1J94aNi8N2j9nRQodIMYsKuEi8OCDD0qGbdu2DZWhDMYolqhh0kujQgkX5eC5MRiiL8Asjcx4qoiD93UiKd955x0pF6MElfgR2HjjjYVyEJIPt1Qdfy0qK3HIkCFCCgHJxO67717ZxZo6EAITJ06UdO3atQuUPmgi4mMTPAK+82effTboZZrOIoBfNsL9nzYDS1XE0nX+/smfEftby/qu2bnnniurEVdccYWQMfjc2vwBA7NiXUUJv7dYccBQCzuC9u3bh14AxB6I81EOvYA6zZAZMSuHkHvAVJYmUUXseW+9++67UkOdESfXUZtssonZZZddZNnLRX5JrjalS2awwBI6PqlhGhGVLjVbZzHSwjCIZWT8zcMWp4gdfWbY+ddrfgw682fFaWqnKmLPe8spYlxpVJJD4Mwzz5TCr7zySvP9998nV5ESJbMc5/a0CRGnEg0CEyZMkIyjWvZ3Lkxwg+MDrhIcgX333VcGR7jtpWkgo4o4eB/HnpIlsC+++EIo8QiGrZIcAsx+WPpi6Rdl7KOwN8ZsmL3hsHxbfWxnknX666+/zL333itViIoudKmlljI87zYCmHGui0m2OU1lQx966KGHSpXTZLSlitjju+y9994T5p6VV15ZLV896Cc4mxGWf3/99VcPajSzCjAL6Wx4Jh5R/WJvmMEOPv1RMpW5WfEzzzwTVVPqNt9evXrJ+xJCG1wP0yCqiD3uJUKiIbos7UcnQWXIzJiHGxYfn4RoUbAKYXFLYHmVaBC47bbbJOO99tormgL+l6sLq6iW05XDvNxyy4lNB4NTeNbTIKqIPe4lp4ijHHl73Hwvq+bYqvAr/vvvv72oI3zHjlHolFNOUUvpiHqFPXi3P7zPPvtEVMo/2boZsSri6mA+/PDD5UKCQaQhgpoq4ur6OZarVBHHAnNFhXTs2NGwVQCr0l133VXRtVElvv/++yUKEDOBqGdqUbUhDfk+9NBDsixN+L2ovRgog/1OtqfYK1apDAFWhZZZZhnDOzQNvN2qiCvr31hTf/jhh1IeL34VPxDAd/S4446TyrhZaNI1u/zyy6UKxxxzjMarjrAzRo8eLbljmRu1ELZyvfXWE6tpwnGqVIYAbHIHHXSQXOTbNlKhlqgiLoSKJ8fgs0UgNVfxB4HDDjvMLLTQQgY+4FdeeSXRijFjwq8VRibllI6uK5iV3n333ULisd9++0VXUF7O2CMgzz//fN5R/RkUAayn8S3GaAsPFJ9FFbGnvYOvKtaZRM1ZfPHFPa1lNqvFbKVHjx7S+KRdma666irZA0M5MDhQiQYBmJog8cAtjCXPOAR6VURnxNWhzVYN5CjwT48fP766TGK6ShVxTEBXWgyRc5CVVlqp0ks1fQwIHHHEETI7YrmSAVMSggvVqFGjpOijjz46iSpkpky3vOmWO+NouCri2lHu3r27ZHLTTTfVnlmEOagijhDcWrJWRVwLetFfy749BiEQPDhlGH2pjUu4/fbbxZAHC9u111678Un9LzQE3nzzTQPLFSsOnTt3Di3fchmxJbXAAgsIqc/06dPLJdfzBRDAup1oVo899pghpKyvoorY057BKhdZccUVPa2hVgviAOTaa69NxEWCchHHJCT/6J/QEXA4d+vWzcw999yh518sQ/Y33axY94mLoVT6OHHciZBFNKs77rijdOIEz6oiThD8UkU7RbzCCiuUSqbnEkSgQ4cOZumllxbXobhdJHDLIDoPNgQoCJVoEGD5H7IUBCO9uIU4u8hLL70Ud9F1U17Xrl2lLT5HZFJF7Ont5hSxzog97SBbrdlmm80cfPDBUsG4GXyccth7773NnHPO6S9IKa/Z2LFjzQ8//GC22GILg29v3KKKuHbECc6BTzYDV1+X+FUR197PkeQwbdo0yVdnxJHAG1qmGO+whMh+7c8//xxavqUyIiKPU/wHHHBAqaR6rkYEnFU8xnlJiFPEL7/8chLF10WZ8803nwRsgQnvzjvv9LJNqog97BZumM8++0xe8Jjgq/iLAFbt22yzjfgpoozjEPyXCX7OvbH11lvHUWQmy5gyZYp59dVXzRJLLGG6dOmSCAasiLHPSX9/9dVXidShHgp1RnaOotS3Nqki9q1HbH0IfYg1Li8AllRU/EbgwAMPlArG5SJxyy23SHn4DjMbV4kGAcdYhlFeUsv/9C+BPJCkyWOiQTmeXHfddVeJyDR58uTYVq4qaZkq4krQiintp59+KiW1atUqphK1mFoQYLSNNS0zKFYyohSsP50idkYoUZaX1byhl4VJC6rEpJalHfZQXSLMzlWqQ6BFixayz09EpgcffLC6TCK8ShVxhOBWm/Unn3wil+qydLUIxnsdvp677bab8AI7JRlVDVD2sK7hY7rOOutEVUzm8yW2M1tE2AAkzWznFLHOiGu7LXlGkfvuu6+2jCK4WhVxBKDWmqWbEasirhXJ+K53gQCidpFwVH1Rh+GLDzn/SiLe9IgRI4Q5rXfv3olX0C1Nv/baa4nXJc0VaN++vVSfGbFvoRFVEXt4Z7nlTVXEHnZOkSpBGoBRDf6eLmpWkaRVH8Za2gWm79SpU9X56IWlEWA2jI0GM6iowx2Wrsk/Z1dbbTUDvzn3le/BC4K0J6k0a621lll22WXNl19+aXwb1KgiTuquKFGuU8RxkcuXqIqeCogARnXEKkaisp5+8cUXzTfffCPW0m6WFLB6miwgAiz7E0gD6du3b8Crok2GEmZAwEAMuk2V6hGAlhYhtrRPoorYp974X12cImb0ppIeBKJ2kbj33nsFDAgKVKJBYMiQIeann34yO+20k9lkk02iKaSKXB2ZiG8zuSqakuglO+64o5T/6KOPJlqPpoWrIm6KiAf/476E6IzYg86ooAq8vCEPYObq9vkruLxsUmdkArWmSvgIMBu+7LLLJOPzzjsv/AJqyNEF9XjjjTdqyEUvJYwlLmH44mNB7YuoIvalJ/5XD2Kefvvtt7IntNhii3lWO61OKQRYnsYgBEOQu+66q1TSis9BzYfVLP6svExUwkdg8ODBQmfZtm1br2bDtFQVcTj93bJlS8ESDnEiavkiqoh96Yn/1QMGHV7kSy65pFhtelY9rU4ZBNyysVtGLpM88OmHH35Y7ovtt98+MXKJwJVNYUIspd1s+JxzzvGuBU4Rv/76697VLW0V2nbbbaXKcQdqKYWTKuJS6CRwzi1LE9VHJX0IMCOGBCJsBh9nXMLyt0r4CPTv318Yl2BgIr6zb7LUUkuZRRZZRHzIsfpVqR4BRwv7+OOPV59JyFeqIg4Z0Fqzcw8ZM2KV9CEAgw8v8j///NMwiw1LJk2aJFk5Y5Ow8tV8jLChYSk966yzmgEDBngLyZprril1e+utt7ytYxoqBjc88vTTT4slug91VkXsQy/k1cEpYkbAKulEIJ84IIwW4D/6+eefG2wG3Ms4jHw1j38QOPXUU8VvGO5utwTsIzau71UR19Y7MKURTOOXX37xxh1MFXFtfRr61U4R64w4dGhjyxByDyQsTluWuRH2tjTIg0AR2p/nn39euLvhCj///PNDyzeKjNZYYw3J9u23344i+0zludlmm0l7n332WS/arYrYi26YWQkX6ozISyrpRACyjUUXXdTAGf7ee+/V3AhcLRBnZFJzhpqBIIBR5JFHHilGcMcff7wQpfgMjVPE77zzjs/VTEXdnCL2xXJaFbFnt40qYs86pIrqMGvFBQZ57LHHqsih8SXsZSFbbLFF4xP6X00IjB07VihJF1poIdOnT5+a8orjYke3qexataO90UYbSSZQ0vogqoh96IW8OuAviiQd8SWvSvqzCgRwM0JqZfDBp/yDDz4QohCNtlRFRxS5BM7mk08+Wc7CLU0ELd8FuxH4zCEewd1KpXoEYCqbbbbZDMv8v/32W/UZhXSlKuKQgAwrG1XEYSGZbD5OEWPtXEukF7eHxQieF4dKOAjAnIXPPrgeeOCB4WQaQy4EgEDefffdGEqr3yLmmWce4e8muIcPbGWqiD261yB1Z6TL0qayannUMVVUBatMuMK/++47GXVXkYVcAl0msvHGG8u3/qkdAayOmQXjruTclmrPNZ4cVBGHh7Pj7/aBJEUVcXj9WnNOLDn93//9n2HPiogrKulGIAzigBdeeEFAcHta6UYk+dqzOnHooYeKu1KvXr1SN8BZddVVBcQwjACT741ka+C2elQRJ9sP3pVOiDtEZ8PedU1VFXKK+Kmnnqrqei569dVX5doNN9yw6jz0wpkIjBgxwjzzzDNigzFo0KCZJ1Lyyyni999/PyU19reaThHr0rS/fZRIzZwBBq4vKulHwFk58+KvRhiYQXk6//zzmxVWWKGaLPSaPATw0T/llFPkyOWXX27YJ0ybrLLKKlJlVcS195wjSPHBL1uXpmvvz9BycDNiVcShQZpoRjzoKNGpU6caZ4RXSYXckhlsT+xnqtSGwCGHHCLRleCT3meffWrLLKGrV155ZbkXuKfYxlKpHgH4/Hk+v/76a4l4V31OtV+pT3ftGIaWA64qiCri0CBNNCOCP7jg8tUQBzgqQ59pFxMFuILC8RmeOHGiuIENHz68giv9SkqoTeKUY+0LYYxK9QhgFOuM35ImSVFFXH0/hn6lmxETOEClPhBwRlbO+rmSVrmXg3tZVHKtpp2JAKsRhx12mBwg1GHaedxbt24tbcG/XKU2BNxSP3zuSYoq4iTRb1I2ri6IKuImwKT43zAUsWNUSjEMiVUdK+mePXsKwT9Rd/iddmF5GklaeaQdR+rvC5aqiD26m5wiJu6oSn0g4BSxc0OqpFXOIEdnxJWg1jjtNddcY+677z6J5TtmzJjGJ1P630orrSQ1/+ijj1LaAn+q7RRx0liqIvbnnhDyB6qjitijTqmxKssvv7z4hbPt4HjEg2T566+/Gqx855xzTtkTDHKNpmmMAAZNvXv3loPDhg0zGOfUg/iiPOoBy1atWkkzpk2blmhzVBEnCn/jwiH0QBZeeOHGJ/S/1CKAQYhj8HnttdcCtwMlwrIqLwqMvlQqQ+CPP/4we+21l4FTunv37mbvvfeuLAOPU8PahnCPqNSGAANlJGnDN1XEtfVjqFf/8MMPkh/MWir1g4AjDnDuSEFa9vHHH0sy99INco2mmYnAqaeeal555RXZAxw6dOjME3Xwy90T3CO18JjXARQ1N4FVktlnn11Wqxi8JSWqiJNCvkC5qogLgFIHh9Zaay1pRSXEAZ999plcs9xyy9UBAvE2YcKECWbIkCFCE3v33XeLr2i8NYi2NHxfYd9j+6KS7Y5oa5XO3PHPX3LJJWVAQxCQpEQVcVLIFyj3xx9/lKOEOlOpHwRcQHfnFxykZZ9++qkkI3CESnAEWK7t0aOHvFgvuugiibAT/Or0pPRlbzM9iBWvqXNnU0VcHKPMnGGJ6aeffpL2qiKur26vhkrv888/FxBUEQe/F4gr27lzZ2HP6tSpkzn22GODX5yylKqIw+swZ8Snijg8TFObE8tMUNbBf6txZ1PbjQUrjvEdbGkYDgV92B0l5hJLLFEwTz3YHIHDDz9cgmTAPHX99ddLONHmqerjiCri8PqxZcuWkhlUl0mJLk0nhXyTct1seIEFFmhyRv+tBwQcg09QNiS396eKOFjvX3HFFebGG28UCsv777/f1PuqklPESVv7Busdv1OpIva7f2Kt3c8//yzlYYihUn8IVEpL6GbEiy++eP2BEXKLHn30UXPiiSdKrqNGjTLOOC7kYrzKzhnxqSKuvVtc2FkX/a72HCvPQWfElWMWyRWqiCOB1ZtM811OylUKewF8yvFBVnKX0mhBxMB+MNs6ffv2ld+lr6iPs87/1Rn11UerkmmFe8Ycj0MStVBFnATqBcr85Zdf5KjOiAuAUweHXDxh5x9cqklsU6BYuBfwcVQpjAAvzvbt24uR4+67727OO++8wgnr8KibEasirr1zHW+DKuLasUx9Dk4RzzfffKlvizagOQJOEQeh0nOc48qw1hxHd4QwgDBnvfvuu4YwkTfffHOmYjazB449Ce+NJBWI6480f7vnzPE4JNEWnREngXqBMp0innfeeQuc1UNpR8C5ITmijlLtcdsUarhXGCWW7g888EDz2GOPSUhDgjpkcQCLdTgS5J4qjKQeBQG3CumeuyRQUUWcBOoFysS1BVFFXACcOjgEaQAsPgRyYNm5lLhBWRaVSylc3LkzzjjDjB07VgJiTJw40bhBjjuflW/X7i+++CIrTY6kne45c89dJIWUyVQVcRmA4joNGQGCH7FK/SHAXi8W0H///XdZWkL3QnAj9fpDo/oWjRgxwgwYMED2zu+8886GgBrV55jeKx0RhSN/SW9Lkq25U8RuMpREbVQRJ4F6gTIh9EBUERcAp04OwWmLOB/hYs36/fff5dRcc81VLEkmjxNPuFevXtJ2CDvatWuXSRxco50i1hmxQ6S6b/ecueeuulxqu0oVcW34hXa1U8Rzzz13aHlqRn4h4HyCnY9wsdq5KDDEIlb5B4FHHnnEHHDAAea///2vufjii83++++feWhUEYdzC7jnzD134eRaWS6qiCvDK7LUbjSmijgyiBPP2CnicjNi90KYY445Eq+zDxWAsGO33XaTZf0+ffqYk046yYdqJV4HH4IVJA5CCBWAUpiY39husHWUhKgiTgL1AmU6ReyWSQok0UMpR8Ax+Hz77bclW8KsD+HlkHV5/vnnza677mp4Po488kjZH846Jq79ThFjAKhSGwKO37+cIWVtpRS/WhVxcWxiPaOKOFa4EynMMfgEVcRYWWdZnnzySbPjjjuKEia04ZVXXlnXgRwq7Wtnc6CKuFLkmqd3z5obBDdPEe2RbD/p0WJbUe5uOdLtV1R0sSZOBQItWrSQepZTxPjJIlBcZlWeffZZYc3Ct7N79+7muuuuyzQehe4DghWgQOBITmpJtVC90nhMFXEaey2COjtFrPuCEYDrSZaOSu/HH38sWaOkXwolKxfDSYg62rZtK6xRGGiNHDkyhlLTVwTLqWx3oISTDOGXPuSa1zjpwa/OiJv3SSJH/vzzTylXZ8SJwB9LoY4py4W8LFao45eGxjFrAkEHe8L4UjMTvuGGG3SvvMRN4AwAy1nil8hCT1kE3IpCUnYZqog9uQ3dS9e9hD2pllYjRASCKuKkDUdCbHJFWd1+++1mzz33NJDbYJiFEs7y8nwQ8FQRB0GpfBpVxOUxykQKVcT1381BGXzcqojbrqh/ZIwZNmyY6dq1q2Fl6IgjjjBDhw7NVBCHavtYFXG1yM28DgMtrKXZEnKD4Jln4/mlM+J4cC5biirishClPoHzEXfkLcUa5NjVyqUrdn2ajrM3d8opp8gMmN8DBw40V111VZqakGhdMdhCdI+4+m5wA143AK4+p+qvnJWln3r4HHLIIcK6Uz0UyV7p/NeSGpEl23o/Sx83blyoz0br1q2loR999FHJfHfaaSdJ9/DDD5dMF/S55QUDL7NvwuCzZ8+ewpTFbGTUqFHm1FNP9a2aXtfH+aZjOZ126dSpUyj3e9DnwqVzA1+2RNyxWr7XWGONsjS2TfuqbmbEuDfAQ5uUH1hTYCv9P+k9ikrrq+nTgwDLvfvss4+56667vKk0bkk777yzWEQTcezuu+8W4yxvKpiSitSTInaWyymBvmg133nnHbP99ttXpIwbFDEgpPWz6qqrGhipiMySVmWsM+Ki93XiJzp37pzaZ2PatGkSQB5lvPfee3uhjKnT5ptvbiZNmiQRqd544w3ToUOHxPs5jRVwivibb75JY/UL1nn8+PGpfd5OP/10w3ZBpcq4QREXRCQlBwkXN3jwYIlclFZl7Gbyzoc0JdBrNVOAwMILLyz7sD4o42eeecZssskm5q233jJrrbWWgbhjhRVWSAGKflbRkcTUkyL2E+lgtWKZu1u3bhKesxJlXBeKGIg22mgjc++99zYo48MPPzxVy9SqiIPd6JqqOgQuvPDCRsqYpeC4Zfjw4aZNmzbCBLXDDjsYeKRbtWoVdzXqqrxFF11U2lOOra2uGu15Y1DGRAtbZ511ZGbMvV7Oz7tuFDF9s91225l77rlHlPG1115rUMZp2XdQRez501UH1UMZn3zyyeIi1KVLF9mXjaNZWKWyZcTzyKyc6EkPPvigcVbkcdShXstwM2JVxH71MAMkooahjN9++23ZMy6ljOtKEdMVbJKnVRn7dStpbeoRgYsuuigUZTx58mTjXO5K4fT555+bTTfd1FxzzTXiozlmzBixkk6KwahUXdN4jm0HLHx/+OGHVK0AphHrSutciTKuO0UMWPnKmBdAGmbGbubOQ6WiCESJAMqYWSmz02pmxrjK7LvvvhINqVQ9oatcf/31zWuvvSb7wCxFc51KeAgwoFlwwQWForEch3l4pWpOQRFAGbtl6lIz47pUxICEMmYfjPV6lDFsPU7ZBQVR0ykC9YrAxRdf3EgZs4oUVA477DBxzejfv78p9PLnOWMZHM5ojIjat28vy3MoZZXwEXDhNb/77rvwM9cca0YAy3aU8dprry3PQaE947pVxKBHg50yxlDEZ2XsZsI6WKj5vtcMAiKQr4z32msv2dIpdyn++o4chH1JFG6+TJ061WyxxRamT58+EqwBpqz77rtP3Avz0+nv8BBwUb1YnlbxEwGUMXvGKGM8BtBN+Wxoda2I6RKnjDEM8VkZqyL28wGq91qhjHv37i3L1OWUMYxgJ5xwQiNILr/8cvPvf/9bjvEbViFcknjxPPHEE8KU5e7tRhfqP6EhoIo4NCgjzaipMmbV1injulfEINtUGRPZxbeZp3tZ+VavSO9MzdwLBAYNGtRIGeMG2FRgfiM2MIxY+QIfNrNfyEKOP/54g4U0+87vvfee2WyzzfKT6u+IEMBgC9EZcUQAh5itU8b40DMzdsp4thDLSDQr9oGbvkCccnMV4wUxduxYc/XVV8uI/dxzz3WnEv92VqSO6jLxCmkF6gYBXtDnnHNOyfZAigPbFYQbe+yxh3nppZeElMBdxBLz008/7f5t9D169Gj5n+hSrDpBaKASHwIYayGF9uvjq4WW5BBgJahfv37u34Lf7dq1k9kwyhjbibqYER966KFmqaWWatTgpkqYkyuuuKIBAOTNN9+Ub1/+OEYt50/sS720HulGgJc0M1VWWsp9CDax3HLLiQVu/vPx4osvllXkGAzBJKRKOP77xcW5VkUcP/ZNS9x6661lBRb9U+rDoBWjR4StnbqYEbsGNQWl0P8TJkwIZJRS6Nooj+mMOEp0s5s3+4flRuf56LAPfNNNNzUcYum5e/fuZX2Gsdj95JNPzDLLLNNwrf6IBwGdEceDc5BSUMR8gsr5558vSetiRhy00T6nc+EPXfAHn+uqdcsOAsQKfvfddwM1WEMYBoIp9ERsKyBN9+9DL0gzjAyBVCjiplP8yNBIMGNVxAmCn+Ki3bMRRRMeeOABc9VVVwXO+qmnnvIiulPgCtdJQrc0rYo42g51zxrfYUsgRZxfgfzf+ZXJP57/Oz9NNb/JC8nf36omH9+vmX322aWKQWgDfW9L1uqXf7/n/87HIf94/u/8NJX+Jp+orOx5qffo0aPi/Pv27St7zJW2RdNXj0CWZsT5z07+73z08o/n/85PU+lv8snXQfwfpgRSxK4CFOx+N30BuP/zz/O7lgq7a13eYTbct7xUEfvWI8Hr4+55rnC/m96z7v/88/x293jw0mamdHnOPBLeL2bCX375ZcUZQuM3cuTIiq/TC6pHAMMf5Jdffqk+k5Rc6Z4fqut+N30O3P/55/nt67NGWwIpYhJWK00BAIxSn6blOFCbHq+3/+eYYw5pEvy/KtlBoNbnI2ykfv/9d8nyjTfeqDprjMMw8lKJB4EsKeJaEPXtWctvi1hN1zJSyM8syO9KFWt+3Sq9tlB98vMrdD6pY6qIk0Ler3LDuMdradErr7zScDkveEjrYaXjM9dcc8nH/XbfHHe/3fdXX30l7oINmemPyBCYd955Je8szIjDBLGaZ83pj2quLVV3UcSlEvhwzjUaENxvH+oVZh3mnHNOyQ5mIhVFICkENthgA4ML0yWXXGJOPPHEpKoRe7m8V/DnZDmeQQTUg0SZ4vP9998LaxXf//nPf2QJmD10Vq94Xn/77TcJQQgHAIQ8cALwwSURI0yebT78hgWLAQ6DF1iW+J893iWXXFL+xx8bX254EVwwh3JgOEWsqxDlkKr9vNM/Yeui2BWxG1EUg8Q11J3P/5/fYQPgykn62ylitzSYdH20/GQQqPT5CLuW7j5s2bJl2Fknnh+B2T/++GNxx/rwww/NZ599JoOOTz/91HzxxRehxfNFIfPBFRFFjfJ2QplBBbsRSIhQzMsuu6xp3bq1WWmllSSkJAMmRwJEhDlEFXFQZP9JV8uzFrYuil0R5yvWymCr79SMkBGdEdd3P5drnT4f5RAqfx4FSOzjF154QRj0UL6PP/54SVISXspLLLGEzESXXnppWZJnxrr44otLvF+IUVq0aCFhVZmB4jLEdhIDF5bj82fB9CEzYz54QTBzRknybLN8DOUos2i+3Sz7888/l9k3x6gvM3N+w9nNp5CgoAmy0apVKzkNs5ZTEIXS67HGCPj0rEWuiGudweZfX24E0xjmdP3nFHG1M2IeQpiNGNkTA5YHmWU1GI/4cP6nn36SFwEvBUbpvCTciyIfLZbQGI3zkuHjXjy8cHgh8WEUziidlxXMPozcWU5jic35ROfnqb8LI5B/fxdOoUdLIYBCw38ZxQvxCAqX56CQcF9uvPHGMqtkZskHZbb88svLt5thFrq2kmP0KWXxcSsMQZeZ88thQPH+++9Le1guh5eYbQOW0KEdZSbPxwlKnTIJtUc7mTVvuummwhnumPtc2ix+h/msha2LAini/ELzf+ePKNxx9+06Oj+NOxb0243u8vOsJT/KrfX6oHWvNB1KDuHFUkioNw8hS1uvvfaaYWmNB5KXDw9rmMEiGMHzyV9SK1SnQsfoKxQ0LzfoDldffXVZUuN/XgphvewKlZ3Usfz7M/93/r3mjrtvV9f8NO5Y0O/8vPJ/15Jn0LKTSsdAlRjjxHZFAcNv3ZSfHaVD5Cci3KCUmDXyYXabJkGpuro3rbd7H/BOYBBy1llnNWzb8X7g44RBPs/elltuKR84xck7jZJ/n+f/zr/n3XH37dqZn8YdC/odZl6FygzUG0EaECRNoQqUOxZVvuXKjft8viIm4DrWqzxMRMPhYXv99debvXDy68gsFcXHrJQlNvb4mK1i9epmsW4mywyXD7NePs5i2+XHSJyZMktpbo8L4xQUM8tlzLiZXTPjZkmNkTiDAurJ7LvpSN3ly/cKK6wgypmII+uuu65Zc8015WWTnyZtv4Pco0HSVNruKPKstA5Rp2cp984775QZ4F133WWmTZvWqEju3Q033NBssskmZqONNpLPqquuKkqpUcI6+wfFsPLKK8tn5513lqAcDMYZqLz88ssG97Mnn3xSlDQD9SlTpsgHGBgMg1nbtm1NmzZtzHbbbZcaxRzkng+SptLbIYo88+sQSBHnX6C/w0WAh4dlJpbUkGuvvVYsVguVwkxzlVVWMbxoMNxAqbEkzKwzzBEuLzc+KOtKhdkJS4Mso2Ecw8z9gw8+EEtUBhXsf/GZOHFiQ9aM2LfaaitZTttiiy0kRqczQGlIFOIPHipeTuCo4hcCPA/MdO+77z4JxcjMN1+4L3fbbTeZ8TLDYwk2zHs/v6w0/eYZYqDMIJpwlnxcMBwGM5MnTxZcGdBDXcr+OZ8BAwaIdXf79u3ludtll11Cfy4YtGOFzsBBpTACooij1vaFi87mUUaszz77rDwYDz/8sMx48/HnpmXEygifPR5mjauttposKaVhn4e6M0Dg01RoJ4YnzJxRysz4Gb0zm37kkUfk465hSRHlvO2228oLgtl9WDJ+/HhzxhlnyAComsFGWPXQfP5BAEWB4mXGe//998uqi8MGxbvDDjsYllOZvfE8NF0mdGmz/M1eNIqYFaym9zTuUrvuuqt8wIjVLhTzc889J6sNPIPEcueDyxp75x07dpQPg51a3zvQnhLfmji9Omgqfpfm7Cn7jsyG2JewtHfPPfeMrcGTJk3KDRo0KGeXkHL2gZHywZyPfank7Iw2Z1828r9dLsrZWWVsdfOhILu8nbPLjzkbEiy3zTbb5OxyeTOM1llnnVzv3r1z9mWRs7OmqqttZww5u+8m+ffs2bNkPnfccYek69y5c8l09XRy//33lzaPQCoXXgAAQABJREFUHj060mbZ7Yzc9ddfn7MzsJxVto36m/456aSTcnZwlrlnoVrQrZGkYGgHtRVnYe1BcuPGjcvR93b7qlFf2K2tnI33nnvooYdyVoFXnLddEWt4ns8+++yS1++xxx5SNu/orIjTA/b7H4WQlYbHoYh52duYrrkjjzwyZ10eGt3Y4L3eeuvlbLD2nI2NnLP7wQK9XaqVdO3atctKVxRtJwMRO0rPXXDBBTk7I87ZWXYjDBnM2BF7btiwYTlrAV40n0InbrzxxkZ53XrrrYWSyTFVxEWhqeqEnfnmbr755pydnTW8nHke7IxL+nnw4ME5u21RVd5Zv8gaQ8p9bffQa4KCd5fdW86dfPLJOTszbvSsoJSPOOKInN1GCzxA2nfffRvysLPh3NNPP120fqqI7cOQFYlKETNavPrqq3PceHa/puHm40XD/wcddFBuzJgxOetGVBBqu2wj19jl2ILns3wQbFlVOPbYY3PWyKsRtqwogNmFF15YFFuHHSP/pi8XZgDFXv6qiB1y1X+zesHMtnv37jm7ZNrQd7yUWf255pprctbor/oC9EpBwO6/CrbW9iFURKzRV85aZOfs9lhD3/FOs25fObu9k7PbTEXLYzDddBDN82e33wpeo4pYFXHBG6PcQV4y99xzjyiIpkvOdi9LbmDrWhBo9GitpOVGZ7asUhoBBjMs9TOz4oXOi4EPStnuaeUuvvjinN0va5bJ8OHDG9K6a/i2BmIFl91UETeDMPAB68+eYynSGhg2YO76Z+jQoap8AyMZLKHbbrFGksEuqCKVtenI9enTR5Rw/vPDM3fDDTc0W51iKy4/nft98MEHFyxdFbEq4oI3RrGD1spZRoPWCKLRjWZ99XKXXnpp2ZlZoXzZS+FGZcSoEhwBZrkjRozIdenSpdFeI8udHTp0yFkr9BzLbdY/O2fduxr1l3sx8H3mmWc2K1QVcTNISh5gYGp9fHPWqlkGRQ5fa9mf69evX8nZU8mM9WRZBLChAO9XX321bNpaE7B1ZI29cihUy5Pd8EwxKOYYK0zW2r3huLsP8r95tpqKKmJVxE3viWb/W2vEHDMqSyfX6Aazfrs5S5BfdLmlWUZFDlhXH8nXWgcXSaGHyyGAUr7uuutEEaCI3YPPywIjMPd/oW/S83LJF1XE+WgU/83eL8uXlvaxAWO2Y7p165Z77LHHAq0IFc9dzwRBwHpYCPbWOjlI8tDS0Pc8c9ZdqqHveb4K2cbkP3ecZ9UkX7KsiGe14KiUQABXIzvKExrHXr16GWsMIdFT7AxKSDeglMTknwgqtQjctQjuSyrVIQA5SY8ePYR5CQISa/wjLmD8dn7axXK2szljrUaFkKRYGj3eGIG3335bfFWhbzz33HPl3rUDVfGD57mwhllCFmGXpBtfqP+FjoCdjUqedvUn9LxLZYirFM+cNcIy3A9HHXWURJaClKiUcN7azWCcVCpZZs6pIi7Q1XZmJS8RaOFwjLf7H+Kfx41jZ0ni58iLx+7nFri6ukM45OML6NisqstFr3IIwFRmLdPFf/G0005zh0t+wxJmXZpKptGTRvxBIdWAQhICGvxSIYLAD9guS8rAtBpuZcW2egTsio5czIAyKYFY6LLLLhPe7iB1gEfh8ssvD5K07tOoIs7rYgIjDBw40DCqZ3YEhyuEAtYNSUgnRo4caaxPac0O7nlFNvpJ8AQEGkmVcBAAS2vNHjgz6BStW1Tg9FlJyMwFwgeYz+wyv/xmFmb9fYXqFEKO3XffPStweNdOiHQQmO2SFOsbLqx1QevAIBk6zqyLUlzaO4CXNaTpvLAZ3SNwIDOjQiEzW41DUMRQwTEgSBtBfRz4VFMGAe7hwq5ELHGI2XrrrSu5pG7TooCty58sPUOPiMBffvTRR5vjjjtOeMzrtvEpapgPM2Kiup1zzjkVoQbToJv0VHRhnSXOtCKGE9n6nxprdduggCFCZ5TPd9x7WwRlQHRGHM5TRlCKapa+iIBlLa5ldSScmqQzF0t2Yqy7ithF0ALi9PJsWKalZjSK6Wxh/dTavauS3HMdMmSIrBxWiioDPOxvsiyZVcTWzF8CKLAfzLKOdX8R/mHrBpDY/bDwwgtL2ZXO4BKrsOcFM8jCUKsaIYIUWxJZFesDbKZOnSrNJ6IX3NwY5WDHoOIfAkkrYqKu8bxVK6NGjZKIcdVen/brMquIMSrh5j3mmGNkWZrltqRFFXF4PcAAi1i1BArAkrSSDzNijObcyoTbrgivdv7mBG4IShhLfuv/ayyjWWR2Ef4ika6aJa2IeZ9iMc1zw3KzMzot9tul41njOWNZm225rEpmFTGjfFxaLIGGN33vFDGjS5XaEMDIDgOiaoVINiyX4YKDW1RWBNyQrl27GgxvXJzsrLRf21kdAsQ35lOtsAKFISBub1mUzFpN45rkkxLm5rNO7nIPqiJO/lHEP9JG6Eq+IgnVgD1yVcIJgV9FsW5v2M2Mq8gi0Ussd7XElk60EgkWnllFnCDmRYt2vpeqiItCpCcUAUWgAAJpV8QFmpSpQ6qIPepuNyP+5ptvPKqVVkURUAR8R0AVse89VLp+qohL4xPrWWcwpoo4Vti1MEUg9Qg4Ri3nT5z6BmWsAaqIPepwp4jL8bR6VGWtiiKgCHiAgGPUcgxbHlRJq1ABAqqIKwAr6qROEUNEoaIIKAKKQFAEdEYcFCk/06ki9qhfWrZsKbXJsj+dR92hVVEEUoOAi7rkojClpuJaUUFAFbFHNwKhFHEZwYcVR3cVRUARUASCIKCKOAhK/qZRRexZ37hZ8VdffeVZzbQ6ioAi4CsCqoh97Zlg9VJFHAyn2FK5qEvTp0+PrUwtSBFQBNKNgKMmdcxo6W5N9mqvitizPldF7FmHaHUUgRQg4PjQs0THmoJuCVxFVcSBoYonoVPEujQdD95aiiJQDwjojDjdvaiK2LP+W3LJJaVGqog96xitjiLgMQJEM0I0TKXHnVSiaqqIS4CTxCmniL/88sskitcyFQFFIIUIqCJOYaflVVkVcR4YPvxURexDL2gdFIF0IUDcX2SuueZKV8W1toKAKmLPbgTiJCNZjcvpWXdodRQB7xHAUAtmLQy1lGva++4qWEFVxAVhSe6gU8T//ve/k6uElqwIKAKpQcCR/+hsODVd1qyiqoibQZLsgSWWWEJGtdBcOpeEZGuUzdJ5uY0ZMyZzjXcDwJ9++ilzbU9rg92yNKx8aZVPP/3UvPDCC2mtfs31nq3mHDSDUBGAKxZlzNI0BlvLLbdcqPlnMTMMWX744Qf5fP/9981+Nz328ccfm2nTpslyX9bwci/Do446ypx88slmwQUXlH1HrHH5QBhR6HfTc4cddphZddVVswZfIu399ddfpdx55pknkfJrKfTnn382/fv3N4MGDcrk8+awU0XskPDoe5lllhFF/Pnnn6sirrJfoPzbZJNNzHvvvWfci6rSrHixVXttpWX5kr5NmzbmnnvukerQ9mrav9JKK5mBAwf60qS6rwfc9EiaFDHvtgEDBpjrrrvO4APN3jaTDmbGWRRdmvaw11HECDerSnUIsLJw4oknVqVIKJFlvssuu6y6wlN8FTNgZJ111qm6Feedd54YDlWdgV5YEQJusDTvvPNWdF0SiT/44ANzyCGHGAZrw4YNMwyY9913X/PWW2+ZDTbYIIkqeVGmKmIvuqFxJZZddlk5kNXRYWM0qv9vv/32M1tssUVVGVxxxRVmkUUWqeraeriIpelqtkXWW289s88++9QDBKlpwy+//CJ19VkRv/3226ZTp06yXcEsGPuXbt26yRYQthhZ38ZQRezh4+ZegJ999pmHtUtPlWaZZRZz+eWXm1lnrew2b9u2renRo0d6GhpBTVnmvPHGGyt2h7ngggsqxjuC6mcqS7c0Pd9883nV7lwuZ66++mqz7bbbmjXXXNNMmDBBbAwOPfRQ8/7775ubb77ZuEmHVxVPoDKVvaESqGAWi3Q3p86Ia+/9jTbayBx00EGBM+Jlds011xiUeNaFF+gpp5wSGIZtttnGtGvXLnB6TRgOAm5G7IsiZrY7fPhwgyvmEUccYR5//HGz0EILmVNPPdVgCMnztfLKK4fT+DrJRRWxhx25/PLLS60++eQTD2uXviphFOL2PsvVnrStWrUqlywz58855xzDYCaIvPjii+bYY4817777bpDkmiYkBLA8Ruaff/6QcqwuGyzue/bsKc/a4YcfbuDLxwOEWfGMGTPEgM8xB1ZXQv1epYrYw75VRRxupxDR6swzzyyb6VZbbWXYG1WZiQBsTTfddJMpt//IfjpGQ+ytr7766gYsR44cadyy6cwc9VfYCDif7yQUMW6BQ4cONeuvv754KVx//fUGv+YDDzxQrO/xS+/Vq5fBeFKlOAKqiItjk9iZli1biivCd999Z9xDllhl6qBgMGR5rJTASjRixAjd3ywAEoY0l156aYEz/xzC9eSJJ54wH330kbyAeek+9dRTss/ODIj99smTJxv2DFXCRyDuGTGWzhMnTjRdu3Y19O/RRx9tXn31VbGUR+ni7XHDDTeYXXfdVbd4Ana3KuKAQMWdzC2PsqeiUj0C+BFjKHL33XeXJMTv169f5i03Qdkpy6Z75BB0dOzYsWBHdO/e3ayxxhpmxRVXlBcwAx9mw1isoyT4vd1228mSP/uEvLRVwkPgxx9/lMyCbr9UUzL3xTPPPCNbD0svvbTp0KGDufXWW2X2u/3228tv9qpZhnY0vZWW0/Seq/T6NKfPrCJ2LxxfO2+FFVaQqsHwpFIdAozaIfVghL7uuuuaN998U14gTXNjD7R3795ND+v/TRC49tprZQaUf5iVBPaR8wUfbAzkmBV/+OGHsi3AdgvGhxdddJEsYzLLPuOMM1Qp5wNX5W+3arbAAgtUmUPhy5zyPemkk2QQxcCKrQfod+k/LOR5th599FGz9957i0V04Zz0aDkEMquIywGT9HmniKdOnZp0VVJZ/pVXXml22203WdrHr/Wll14SEgGWWKFpdMIeKH6NuoflECn+vdhii8nsNn/mglWsc7crdCXEDeeee65Yy7J8zR48+eC+ArUhe4vc65CvTJkyRQgeCuWjx4ojEOaMGDrYBx54wBx55JGGmS/K95JLLpFBFERDKGWeJQzy+vTpI2mK10zPBEVAFXFQpGJOxzIfooq4MuDZv+Ilcswxx8gyK0pg7NixDf6wq6yyijnuuOMaMj3ttNNqYpFqyCgjP3beeWfBluYyA+vbt2+glqO8MeBigASH+iOPPCKuLRjSseozePBgA70mRl/sPeLDjNWtSnkEMJhCcBGqRvDOwKVojz32MAsvvLBp3769sF7RT7hSnnDCCbK6wYrGxRdfnGkGrGrwDXKNmrIFQSmBNKqIKwcdC90dd9zRPPvssxKYYNy4cQWXolkSHT16tFl00UXN6aefXnlBGb/iwgsvNI899pgsR4JhpYJx1w477CAfFDP9ddddd5kHH3zQvPbaa7LfyP4jynvttdcWRqaNN95Y+jZ/NaPScus1PUFLEJRoEMGqmcEpM9uHH35YVifyr9twww3F0Gr33XeXFYv8FZD8dPo7PARUEYeHZag5saSHYImqUh6Bb7/9VpTuc889JwqWF3sxektmcigTDIz0xV4e26Yp2Be+5ZZbQiFlgPWMfuJDnzALvuOOO0RB3Hfffeb111+XD3VAITAjp9+23HJLg5FQtbPApm1K8/9OERfDgpks+/Xs5eLr+8YbbzQY5dFutmdQupCx7LnnnqZFixZphiOVdVdF7Gm3MSPmxYPV9N9//92wtOppdROtFiEjYYFi0II7xZNPPikWvKUqdcABB5Q6refKIIAlehQCAQTuMHy479lXRoFMmjRJZs7sX/Jx7lTsL2OQB6sXdUKho1iyJE4Rs6zPHu+9995rCK7w/PPPy747bpD5gj0EePFhYAN+ldLA5uenv2tHQBVx7RhGkgNcv7gBoGTgnG6lbE8FcWYPnWVO9hmJGMRLWtl7CkKVuoMsYbNvzAdBMbO3zFI2M7ynn35aBqoMVlnKdoIVPAZka621ltl0003l2VlttdXqStn897//FXci7nv2iBm0E3CjEC0ubk3sz7OKwBI/KwmqeN3d4se3KmI/+qFgLeBjRREzulVF3BwifITZE8aFAncKrG6LLc81v1qPpA0BFDMzOD4IyoioPsz8+Lz88svygWqTz/jx4xuaiKJi7xOyHLZ9+DD75r5BafsYaQtaSBQrzz/vAVZ8CBfIwLwpvwCuRqRlNYAZLkp5s802k9+tW7cWRd0Ahv7wDgFVxN51ycwKYeGLcuFBJCKQykwEcH+BJALLTpalCWafBMXfzBrpr7gRYFbHrJePi5aFQoIw5J133hG/cfaYGbChxFDOxQQlzzI3qykoZVx3cLNiNsk3AzzuL+wLoPtkxYq98jnnnFPsDLie+vBhgMDsnW+Wit0HY0JILyA5wfeXmSy2DSwto3S5l7/55hsDLSSWzFxfTBhYMIigvizb843fPApYJX0IqCL2uM8YySIoYpWZCICHU8Iss2Ftq37AM/HJ8i8UFL7JfPIFBf3KK6/IbBKljB+sU3jcTyhLyEf4+CIoepbUmbEze2cPnJk8W1YoXNrKEj3Lzq3s1pUqYV96rvJ6qCKuHLPYrmBGjDCiV/kHAWYNKGFeonwzE1YlrHdHOQRQWhtssEFRH1hmqdxTLPlCkMEMlQ/HYZJi9sqMlnPMan/77Tehd0SB//nnn41mwcyK3QyZGTMflCqMY4QqZFbN7JpZNrNt/ucbwgxckJjdMjsnj3LCjBqpxo2sXN56Pj4EGhQxN2qWhCUrHiYeDl+FpSeEZVgVI64tWHqyX7b55pub+++/X15yUWODP3LWng8IHphV4iqUBUEx8ry5Zy4tbWawgNSDIsa/GenUqZN8Z+nPrFlqbH5bsbZlOYewXezL+Ci4MDHbm2YtIxl5Z1l4SHfaaSdZOsQqFutoZhoq0SBAtCqWQuEQxlJZxU8EnCJmRp1W4V0MsQ5bTJkVu3cSu9x+++3EQ8ttvfXWsZdNgXZPKGeXiKQO1MPOdnL2pZOzvLk5yziTs1a4idSrUKF2eVrqaZ3wC53OxDGrhHPWIEtwAA/rF1n37XbttYY4sbfVvhhzlgY0Z2eJDc8IuPPcqviFgOXolj4aNGiQXxUrURtrhJazxCI5y3CXs8xpDfcY72Lue2vIVuLq8E/ZLQmpg7U4Dz/zgDkmMiPG4RyBlD8JwS0o3+/QYiVuAcOGDTP77ruv7NWQBkvMG2xcTUZsSQnGGgjGJVkVwuxhPY4lK7OzoFR+acbLkVLAnR23sD952WWXiWETEXZgWmJ7pEuXLrJcPWHChEbMTHHXT8ubiQD71whuWT4L++jMeOGBh78af+bzzz9fWL5cvbGJYaWLbYI4ZdSoUVLcX3/9JfYAcZbdUFZAhR1aMkZD1hhBRiDWJy60fKvJyBJANBqNWVCK/m8NKXJWSeesss5Rb9oRh9j4rVInG7wgjuK8K8O13xq25Oy+vnf1i6pC1l1N+t2+vKIqInC+9iWas/STOWtU1PB82LCSOfvSDJyHJowGAXef+NgX1i0rd/PNN+fs9kaje6fQe9YauOWs7Uc0IJXI1W755ewAv+G+ZlU0CWFkG6tYUndpNIotabHh7xo6oNDNUeqYJQeIpfo2qLrUkUFA1sT1j90nz9k9y0w13xKVSL9bUn5v2s1LywZpyFn3mYbnhiVrS6bhTR2zVhE3mbCuWV403dqz5IYMGZKzbHc5lnpLvUPzz3Xo0CGR+tvZcKM6UmfrpRJ7XWJXxOxl0AEHH3xw7I1tWuCvv/6as9aGjToi/+Yo9tsu3eXsklDT7CL53wYxkPoxA8mSWJaknA3IIG23Aemz1HRpq/WPlrZbnmXv2s5zY8Ph5eyStdQRG4v99tsvZ0kovKtrvVfIhpGUPmCfM2mxTGZSl2LvzVLHsc2JW1jVdAOZ/LrZwBdxVyUXuyK2+8LSWTYMXeyNLVSgjUdb0c1jmXZyb775ZqGsIjlm/RjFmMz6IebsfmEkZfiWqaXwy7kXTO/evX2rXiz1wZCRl4PdG4+lvGoKYcl64MCBDTMf6wqYO/vss3PWLbCa7PSaChGwe5o5lnStv7E37wYbB7yi9yn3OEaB1ke7wtbXnvyhhx4qWlcbOKb2AirIIVZFbGnfctbQRhrvy+jZ8rPmWPrMHxEV+00667taAbzhJF1++eWlftZgK5wMPc6FwYZ1T5L2sv+VlcFH0y6xfMGCASsivgteBp07d254hhhEqYV19L3GnirvKsu6FX1hAUtglmmNXBvuhWLv0vzjBx54YMDcw01mOcuL1tPyFMRmB0SrYrWatjNJ4VXFfxfaNh8EC7499tgjUFUGDx4sMTsDJQ4xEVy6CITv9S5HHXWUcAJzj9iXeSB2oXrExJEbwMjku2DNTgxhQhZCszh9+nSxsCbGLZzJKtEgAMscAhOXLwLxDWQw++yzT+Aq2W2NwGnDSkhMZjsjLprdM888YyDyiUtiVcSEL0MIx+WTHHvssWWrg4IgRmoS4hQxA5l6lptuuskMHz5c6ACJnAPhflbF7sNK0wkukBaB85gISAxYYayDfhQ3wH79+pUMYJCW9vlWTyg5EZ8UMfWBmtNuPRoGYuWEusMXH7cQz9pOREsWa7ctDS5NcUisitgus0mbCNPlk9j9uJKE6TA4wfySlDhFzCiuXgWy/cMPP1yad8UVVzQj7a/XdhdrF9F5EHiI0yTMiI4//niJdsSsCD/oc845RwITEBFJJTwEnCImCIRvgh/80KFDyw6m4W0IwqkdZvvAbcyYMWWz5J109dVXl00XRoJYFTExQxGCdfsmxWbFjOyhl2zfvr2ELkui3pZ9RoqtV0UMvlApQqrPy/uQQw5JAmavyiS4AJI2RexAZKZzyy23CIkDioKlPrirrdeEzo4dSDV+w7mOsDXgm0DLC2kH93EpKtoklqWtC54E6giC2XnnnSdtCJK2pjThbn8Xzw2XB4yd+EBZ6Jtg6dnUlcnyt+Zwo7HsVrKpj/FQXEQe+fiAF/5tWEeCY72Jo+mzId5ip7fzEUvuRftQ5+wLzMfqVVwnG60od8ABBzQYxtiA9Tn8TVVqQ8AZRdk92doyCvnql156SQzIuIdtxCvxy3UGmBxzHxtQJOSSy2eHF4ozGHb1KPeNZ03UEpvVNNyiNBhOZ1+lb9++DTcJL0FrfCJV5aVhKeTknJ05J1J9uzwt5dcbecJjjz0mLhj4otZb26q9UbDk51mB8aeexAaubyADgaUL1iWV6hFwVr/33Xdf9ZmEfCXPM66W3L+44DneaBtYpxmvdP/+/UMuvXx2kI2UU7xNz+OWxzMZpcSmiEeMGCEA4Pjvq+C/6lyZYFzJFxuAW2YodFISBBPgRtngWC9il61yzjUrqxSehfqSVRj62logFzqd6mP4i1ovhYaXITNlZikqlSPAjJL7BLZCH8RGsmvwKbf88Dn8zPPFWnnnYGKjzgy8P/744/zTkf/GFdJGtGu496hH0M9BBx0Uaf1i2yN+++23bZuNcYZH8o9nfwjMbVlVDNZy9gXRqHZbbLFFw8Y9FtRYh8YplllLirMPXZzFRlqWXZIW9xaM98Bc5R8EHJF/mkPbFetLrMAJGmEHs2beeec1N954o+wl1rtHQDE8ajlu/bflclwwk5aLLrpIguRgZUygEPoVg618sf7OErTFDr7FeA8XxTiF+65UAB8MDfHUwL4Ba3/LumXsNopYdX/77bdiwxJZfSNV83mZ77LLLjL6sGDkHfXvp7WUK7kPbC1CpR1wZcdFcwlKkP/bmyCx0JFh9xQcyoyKobEkLKXKTASYWdDXSREdzKxJtL8gqLFKRNrKVpCSgATH+4cffhDcCOeapMA/zn3K/cpqorUyLlsdyJys33nZdGEnePXVV3MsnVs3Wgkgw7sealDqQv1Z6k9KYluadksC1oUhqbaGUi7LLdZfUjquTZs2JZV2KAX+LxNLkiBlsrcGQ1maBeMdFzhgwIABaW5KJHW3lprS13EYiUTSgAoyZamaZT9ehHygSGy6pFlBdplJSiQy8ErC4MmBjBImuAP1sCsdOZ8ClLg6Bvm2rkzShq5duwZJHkmaWJam8SW0m93G8qIaq5Btv6VXWG6BbIL4n5MnTzbW4CCWxlAebgr4l3700UexlBlVIXY/2ODLxzbFySefHFUxqc3XLTn66JYSNqgsVdsIY+YGG/fbro4YfMht5CljLcfDLqqu8nOMZUkxFOKehDuaDUoiS7mwqtFvaZTvvvtOqp1knPNYFDFKGGXMi4WHLe3C3h0sUAwsICuwhlyxNIkbH7Ehz2IpL4pC2AuEeQknfl7AdjkrimJSnScvOSTuPTQpNKE/dnnT2FCXBjsNvu1Mr1HQ+ISq5W2xvFMR9lvjlhdffNFgs4LdjzW+MvBDWDeluKsRWnkzZsyQvJK0yYhFEbsXiw0fGBp4SWdkfYrNKaecIgMMu6Rh7HJr5FVyN3vchmJhNqxnz55CGweLlvUtDDPrusnLWpNKW+rpeQnSORD9WB9UYbnjnYERXyk+4CB51msa906Ne0YM//s222wjK3N2T9VYt1QZPKUZZ2ccyapjUhKLIrZuQdK+uG+aqEG1+5ti8clSot3birq4hlFnWhWxjTkqo2cbncdccsklkeOVxgIsYUxDoIQszYhdX/EyhIVr//33NwS+sEae5qqrrnKn9ft/CDhFHOdg7cILLxTmO7YNunXrZu69997UMr/l30iZUcQ+c6Lmd0ilvzF3h9ycfS72uKyxQqVZVJR+ww03lPTMGqzFQEXXJp0Y+spTTz1VqsEAphTtXdJ1TbJ8ZsNQfuKSwn2VRSHiFM8VgzUGJrgLcs+ozETAKeI4Bmvcj7hz9unTx/DOo18sGUvdbCs5qtAkObtjmRH7GK5r5i1d269VV13VMFJEuFndxn9tuRa+mv0zfPEowy1fFk7p31E4hlkZYTnaWsn6V0FPavTee+9JTbivsi74mTsjLoKusI+MrYmKaXj+o1bERAEjOhIDI/y+LYuXoV/qSXyYKMaiiN3UnyXJehRG7Ntuu6356quvTO/evSNtottXxWAiLQIuKGIEwnWM3FQKI2B9a+WEKuJ/8GFwy0oTqwOQRKAUsm5RjeeEpYwUTBiYRyXWRcpYSmLz9NNPS/x4jFLbtWsXVXGJ5MuKCxNFZvp1PyOGlQRp0aJFImBHXSidCFMQkZoYwU+ZMiWyIp0ixkgiLUI8WozZLLWhl5G3fMLRRdjiBajyDwIYB+EqaIOyGNxk2De2tJiZhcethrE/zLsnCmH/l+hJLIHzjWW0Y/eLoryk8mRZGjYwBjRJbpfFMjVxy7X1qoi5iVq3bm3OPvtsuZ8OPfRQ2eeL4ubioUB4MNIgxPS0TFGyn3TZZZelocqJ1tEpYuj1VGYiwH3PzAxljFLeaaedjGWXmpkgQ78cTWMUnAzYnrAf37FjRwkVaDmjxZ2sXlcznT92Em5g+bdsLIrYxVZdaKGF8suuu9/snTCTsZSNDfvGYTfSKWKWpi3DVtjZh54fftaMONnfS/pmD71xIWfIMtlbb70ls5x6nH3UCheD3eeee84wE7Q0hQYXwiwqYwa3iA0bWiukja63jGamU6dOhv14ZODAgbIdgPFcvUqcRm+lMIxFETsfWzb761lg3cLVguUiiM8dQ1KYbcbpnJEwRhS8tH0WBiSWPk6WfM4880yfq+pF3ZgN47LDC9ZyCHtRJ98qwb0P4QfPAYPRLM6MHbNemIoYgpDNN9/c3HnnneKSZENWNng5+HYPhFmfqAY1ldYxFkWM6wqShZcL+1n42PFCjcpwC+IDhNmBzwKVJbM8SDx0Nly+p3BLQ5wdQPkrspkC7wF86YmQg60EBlwMTLMiTnnQ/jCEAT0rbWDKqgOrDZB1ZEHCxrJazGJRxM7KEWOmLAhLOsz+b7vttkiUZRoUMTf4LbfcIpSmzn84C31fSxudAZ4q4vIoooyxk2CZGspXDLjce6b81elOwUoTgtKsVYYOHSqc0Xi2bLfddrLKsPrqq9eabWquDxPLWhoduSJmRsReJi4r8AtnQXhJQH+JHHfccaGTbxAjE2Hk6qsQnxSfT4w96o1RLSrMMUZCXP9GVU695AtJP0EHeN7wVNh9991TYTdRC/4MNtjyYhusllUm9oP32Wcfc/TRR8tzyvsKNzEb3a2W6qXqWgzTbDRAqfNqq62WbN0jiemUlymhsmwLJe5s3uG6/0moPxtgWto+bty4UNsLptaAImcHN7nvv/8+1LzDyMz6DTfUz1p4hpFl3edhDRpzdqCaIy6vfdnWfXvDbKDdM83ZPWN51vbbb7/YQpOG2YageVk7AmmnDbYQ9JJm6ex+cM4uRUs+1j8798gjjzRLk4UDlmBIMLC0qok3N/IZsbPszcps2A2rWJp2BkpYIdqedqdq/iaCFQEgWG3wcZ8Y0g72yHGBiJMLt2ZgE8wAfmWeFfq1nq1Uo4AYA64HHnjAzD///EK9ePzxx0dRjBd51sq8BikHVLlsg4AbqzA2prAXbYu7Es7YlUhfSUvkitixKKE0siaHHXaYKCLYkgibGKZg4Yi45cww864lL3hpHUn/SSedVEtWmboW31gEYz+VyhFgADNhwgSx0B8yZIipV5/1999/X8Ah/GAlwkTg4osvNm3atDGE/cPA7bXXXqtLko6guDiffeKiJy2RK2I3E3Yz46QbHGf5tB1WKeT8888Pdf9qyy23lHx9U8SQwUPgQgi7LbbYQuqof8ojMGnSJEmEwYxKdQgws8NAEPdBPBZwxak3qWZGDAtZly5dxG6F9zArdOwHZ8GLpVT/MxBBvCDPiXpx3M6Ec/bBkA+/syaWzCJnrRtlL8Jy5YbWfPZh7T2UY4/HGkWFlm+tGa2//vqht7XWOvl+vSWlkP1+a4CTw7ZApTYE+vfv3/Bs2OD1tWXm2dXWY0LaZn2pA9XMrsbl7PaQXIP9wT333BPouiwkWnvttQUXu72XeHPZu4xcrNuSNNj6+kVelo8FjBw5UtpvWbdCNSSxfoSSr/U/9aLZdp9T6mNpCHMYlKkEQwBjPgZVdlk62AWaqiwCXbt2FUytRXXORtcpmz4tCRZccEFpl3U3KltluzqVs7NeSW9nfTnrqlP2mqwksNwWudlmm00+PuilyJemmfo7/+EsOd3TbifWklPiy2IcAJl6WLLVVltJVk8++WRYWdaUD4EvkB49eoj/cE2ZZejiBx98UFqbFRKFOLp21KhRZuuttxZXnz333DPUbaE46l+oDML1QReM2xbMYsUEexx473nvwGrIN2QxYRGAFCs3TcdZlsa9Ekpip5+SrH8sithRWzqGrSQbnETZ+Py5GJ4YTIQlPiliQrNBYML+3CGHHBJWEzORjxuc1VuIuSQ7D88C9ojxYcez4PDDD0+yOqGU7XxeSxkXEZkJwh9Y7bC+Hz58uBiK2tlfKHWol0xcGFksyH2QWBSxC/aQRYJ218lEZMK9gjBuUMmFIU4Rk2fScvvtt8vom7jMYTD+JN2euMp/9dVXDTMdYqHa/fW4is1EOYsssoi56667JG7viBEjJFRpmhtu97ul+sWYrxjQ4YqDkrEcBkL4g+eGSnMEHBmSYylsniLeI7EoYpZSEEs+EW/rPCqNVQGUMYJ7RRhC8HjiaFrDLePcGsLIt5o8Ro8eLZcRZUklOAJ33323JN5tt91kNSH4lZoyCALrrbeeufrqqyXpMcccIy47Qa7zMY1TxE39XmHJYsUNZjH894n7TVqN4FW8F/HbR3xhsYtVEbu4xMXhqe8z0Mnh0oSLhQsNWWuL2QdDiEiTlBC5hfKtBbfp3LlzUtVIZbnOxYaXqEo0CECzetRRR0mMcO7PsJ69aGpbPNc333xTTrKv6QSOaFwZBw8eLDG/hw0bJv7UbhXSpdPvmQhMnz7dsISP+1apZf6ZV0T/KxZF7IJKA0CWBZap9u3bywvBGTbViocjgIBrNykZO3asMIcxEmf5XSUYAgTGIGCBtYTNLLtRMKRqT3XJJZcIoxQhBFHMaRTHBOUU8fjx4w0cySxFsxfOLK8e9sKj7hs3aWEA4winoi6zXP6xKmKWULMubs8GRWxdBWqGA6YcJElFjJEWAom8SnAE2FdHmA1bH8/gF2rKihEAX/aLGfRYX9rU7Rd/8cUXsrXHNh9734QWZXbPdh8DYGwNfDE8qrhzYr7A2dS41cSYiy9YXCyKmH1MRBWxkXBtSy+9tOzpuhuiYM8EPMjo2PrtGktgblzA8ICXhpKMJR6Mz3jBqftNZZDeeuutcsFee+1V2YWauioEeO7cfvEJJ5xgoJ5Ni7hlaQwhCZN5/fXXy2wOPnuoPZ0dTlrak2Q93aQlc4rYhcH75JNPksTfi7LZIz7ooIOkLpZpq+Y64S6EpTLiaBJrzrSCDO644w5JrbO6CkCzSXmx4svI7EbdlirDrpbUlujDYFCIKyUrOPiSpkEcLzJbGdw7cE2zJI2bkkpwBNhTB0vsWXyxmKb2scyIXdxMjHpUTMMeFUuTYQQzd/zESShilvsQlsdUgiPggoDAAYzPq0p8CFxxxRXibfD6668LB3x8JVdXEm6fjn/AUuaK4Rl1V3e3yvHkHcmWIK6fPm0HxaqIp02bFsq+aOXw+3UFbkeYzUOCwX5VreIUMUHS45RvvvlGDEQgDthpp53iLDrVZUG879y99t9//1S3JY2Vx6CQbQEMdQjGglLzVWDNwwWLmRxCEBnCjPqkRHzFrlC9bOxlOexb6MdYFPECCyxgbPBlA8Ul5AUqxrBEhrh9wlowwcEfB36s0h37Ti35Bb124sSJEhOZkGpZj+QSFDPSQWnJc8DyoouiVcn1mrZ2BPA2IG4xgyIGQ/ji+iTQVLL/y7YTW3psQTFw0NCitfWSo5P1beIQiyIGOhc/M2niidq6MbyrMdDhwbr//vsNYcpqER5SlCES56yYYOwILlkqwRG47rrrJPHBBx+sJB7BYQs9JbNhjJ/YMxw4cGDo+VebIW5thBClfo4yluXUlVZayTi64GrzzvJ17K1j1MqkxTeyk9gVsYunmeUbgrZjwclsiD1iZpa1StyKmBE7gwhEjY2C996XX34p2xEMwpSFLDhuUaSE7B83QpTdOeecE+tqUrH2MEhj7xd+bN4R+Lw6616WqFWqR8C9r/DuoM99ktgUsXNCd9Z/PoGQVF2ICoM4msNa6rHjjjvK5ZMnT47FEhTrTYxIMMTTqC7Be44XPwY3nTp1kpF58Cs1ZRQIsPSLbz8DyyOOOCIxGxbsLfA8IGAKEZN4N0BTycwYH2FEFXFtd4Czx+nQoUNtGUVwdWyK2AZhluqrIp7Zi47WkJEaL+daBBcxlv9Rji6ySC35lbvWWWi7AUC59HreyAAJi13kyCOPVEg8QeCCCy4wsP/hX0r4xLiFYA28H1EUUFPaOMIG1ixsaxAGvYhaSQsMVf359ttvzdNPPy1Gbj7yHcSuiPGdDINRqqre8Owi9nwgcIcdh5ukVnFK8eGHH641q7LXO0XsLLbLXqAJDD7XzHwYMClu/twQkGEMGjRIKnTKKafExkWNDzMGWQT8gOwIAzJIebp169YADu9KVcQNcFT9g8EOhnlt2rTxkoY3NkUMuxah3nDZwRhB5R8EnKGTs+arBRdnCfjQQw/Vkk3Za1nGcwMHp/zLXqQJhJgfGNTy1b+bActplqlnzJhhzj777MgryN4vbowYZOH+d+mllwohDwQv+TJ16lQZGPDudAyF+ef1dzAExo0bJwnddmCwq2JMZUdcsYldioVcOTdmzJjYyvS9IKuABRPLE1tzVW1Umdzss8+es0HAc/yOSuz2gtTZ7g9HVUTd5WvpTAUzS0eas258dde+emiQXa3LWeY7+Vij0kiaZN2kcpZeM2eNheR+sIFgcnYvuGhZ1r1R0u26665F0+iJ0ghYr5ScHezkrIFkzq48lE6c0NnYZsSMLTbeeGMZYjz//PPyrX+MjMIZEb/00kuGfYxahD0liEJY8nrsscdqyarktU899ZScx5BEJRgCF154oSRkbxhrXRX/EFhnnXVMr169ZAkTH+OwhfcebjOELITq9qyzzjJ4kcADUEx4LyAa0KEYQuWPYwxLnGbYtFwkwPJXxZsiVkW8+eabS+vcizzepvpZGgw5DpcwgkA4QwTn4xtFq91AytU7ijLqKU+Ym+677z7xASU4vYq/CODGRAATDCjDsrWALKRv377irgjhDqELcU+iLLuCVRIMVcQl4Ql0kvjviM/R4WJVxMzW7LKpmONDuq7yDwIuaIOLClILLs6nN0pF7F4ORIFRKY8AVrl2xUtcU4iUpeIvAvQPShM58cQTazYsxdCKGS/3AMZCuEhhsLrBBhuUBYF7xnlAbLLJJmXTa4LmCOBFwrsQveN1lLO4l8Tt8rTsedjRZtxFe1ueXUYWTMLYJ7aGVDn2Ie0tWXLvqVowLAGJ7EGzD20HU9Vmk5nrbDB32ZuygR1yltYyM+1Oc0OtK2HOugPKM2RdiapqCs9Jnz595FnhWbSz4Jw1cKwoL/apudYSe1R0nSaeicDw4cMFQ2vIOvOgh79inREzXsF8HIlyD1MKSNEfRruM2FjCrDUaE4wxzhLbMcmECQU0cexBs7xGKDGV0ggQpg4r80MPPVQJPEpD5c1ZnkVnOc13pT7+cEMz44U2k75nLxj+hEq3cl544QXBxNnWeANQiiqSluAqsStiF/XCRcFIUZ9GVlX4Y9daay154N2yby2F7bLLLnJ5FIrYEbI4gpZa6lnv14IVoS4xxnPLnfXe5nppH/SjDDZxtbzhhhsCNYvtNoy8VlxxRaHLhCMANz/2glHulYoq4koRa5yesLvYIxGQBiY7nyV2RYzTOvFX2fuo1UrYZ2ArrZsLUu0MoSq9Pj89/sTMjPFVDHsvXhVxPtKlf9ulyQbqRPxAVdKDAFbN/WzIQWTAgAFlZ8UYWkLje/nll0swFwZeL7/8ck3B5zHoQnR/WGCo+A90snYV2nTu3Nn7YBmxK2JcNzBOAqAwSCwq7h1PL3DGG45XtpZqQgpAQAmsNcPGGP5bxHGH11LPer6WARXBPHAp09lwOnu6S5cuwnxHHPVi1JffffedIYoWEwyWpHGBsnYBpn///jXFDObZ5V1AcBBdmq78/sEwzq1k9OjRo/IMYr4idkVM+9zSKS4dKv8g4Ajdw1DE5OiIzcPG2IWxhBVIpTACDDKPPvpoOQmLllpKF8bJ96MoQSgoERQrL/d8wS0Gi2he+Eww2BNma8mFfM1PW+lv3gP4vrI8jjuVSmUIMAH5/PPPpS9c9KrKcog5dRIGZB988IFYslmC85w1/EmiCt6VCduSY/WxD2DN9bOGX4KxdWDP2RdIzfmRAfWijlhMwxCkUhgBayDSYO1qI+kUTqRHU4EA7yerWKU/x44dK3W2fOE5rHDtq1o+1gA1Z/niQ22PXeKWvO1sO9R8s5KZ5e8W/CyRTiqanMiMmLB5LG3i46XW0/+MvBhRY+TBqNsOVGoejmFM1apVKzN9+nTjjD5qzfTjjz+W+pFvOSKCWstK6/V//PGHYW8YYW9RA7mntSf/qTd7xSeffLL8Q39ecsklEvoTPneCRYwYMULeYURNClOeffZZya5SS+sw65DWvDDSYlsIsiS2DdIgiShigGEDHXFk3PJPxv+wDIW8++67oSBBVBeEyCNhCHtlCIpYpTACWMh+8cUX4r7SvXv3won0aKoQwIIauwsMFdlqwABy3333lee0Z8+ekQSZf+aZZwQjSJBUKkPA+g7LhAEds9hii1V2cUKpE1PEzpwcRYxfqopp4JwNWxHDtRqGqCIujSL750TRYW9x2LBhkbygS9dAz4aNAEZT+IITqhTBd559YBu4xrRs2TLs4iQ/S/xieNbYG1ajyMogZl8da2nkqKOOquziBFMnpoghP8fQgfis6lP8zx3Akj1CTNIwBOt0HmaIQlhWrlU+++wzycKyDtWaVV1ef/jhhxuWprHSVJeT9HexZdWSZWhCFeIHzHYMhDtRG0+52TAujQzqVIIjYPfxJZQl9LtpCkqTaC+7ANjc8CpG9ojBIQylST74azvu6TBmxV9++SXZSlxp+aF/GhCAwWfSpEmyhIn1rEp6EeD5w9KWGMU2bJ7MSrFiZoBlLX/EVzjK1rmgOGlSJFHiETRv+mbQoEGS/Nhjjw16mRfpElXE++23n4AwYcIEY2NGegFIkpWwsUml+LAUMZl17NhR8rzzzjvlu5Y/ThEvueSStWRTd9dCTOOWwQhx16JFi7prYxYahKHkRRddJC4vTz75pLidjRw5UvaGYckihCVy4403hk6Uk4+vU8RwAagER4BoWfAcQJ7jc6Slgi1K2rbbLp+Kmbld10+6KomXb5c1JWA47kFhuRxZy/QcAQfIc8aMGTW10fo6S1/ZPbKa8qm3i7t27Sq42BlMvTUtM+257bbbcssuu6z0o10Oztm4xDk7wGrWfhvTVtJcc801zc6FcQB3N7sELs8rAe1VgiPQtm1b6Ru7IhX8Ik9SstSSqNjRpYBnzfQTrYcvhVsrP8HDzj5Dq5JdnpY87ei+pjzdi8oaktSUTz1dPH78eMHW8tnmpk6dWk9Ny0RbrN1DDj9gO0uRzzLLLJOzrGhF237TTTdJOqLIRSGPPvqo5B9GJLYo6udrnpYyWXDjOUzjACbRpWmm6CwhYPyAgQJxOrMujpMYy8mwZI899pCsrNKoKUvo/BBdev0HRixp7cxJ/mFf2G0t/HNW//qMwE8//SSuSPTZ5MmT5Z7G7QWaylKUkrjE4D+Mbz5GkGELnNVIKtigwm58DfkR7xkh3jNBHtImiStiDIoOOuggwe2qq65KG36h13eJJZaQPL/++uvQ8kYRY33JHkq1e/G4ceA/ieVoGm/00MDMywiyALvcb+yMqmGPOO+0/vQQAcISsu8LRSvkHHaGZ7B2x/XssMMOK2ulTCQtZ2R6/fXXh95CVcSVQ0qELCYZ9M0JJ5xQeQY+XOHDUoP1m5W9UcuEkmNPM8tiH3JZYmEJLEyxI2zJ1/LjVpUttH72fs3Z2XBV19fbRdg0gAc0rZbJp96aV5ftsS/sHEvK9BsfawyVY2m6UsFGguutH3HOxiqu9PKi6aGNtX7Kknet9hxFC6nDE5ZcRTCzxnSpbV3iM2J7Q8vo1G60iw8my0NZFhh8kLBDRNbKZOZm0jobNhJrlrizCMQddu9cfusfPxFgmwdXpNatW8uSst0HNrhMMvvkd6VCpDSsqFm1CjO6GaFhLee85K2BQoL1yjvvvGOsoZ3QWTpq2WBX+pXKC0UMJCeeeKIgc8UVV5SN/ekXhOHWhv0nBB7uMAUmM2IUw8HKw16puLjGWVfEEDpYK2lZpof6kN8qfiLAfQ7lKNGQULxwD59xxhlCTcnyMs9DtYJiR8LkQJgyZYrkCRGPSjAEYD3D7Qyq0TQPiP+/vXMBuqqs/v/z//mf6QJeSzEdBUGQAkXNIpMIhQYRxsuIjRqSFIioCIWheCE1JRXpQsJAIaJgSVAioKgx/PA6oCQJSCIoVARjDRbZZKV1fuuz7Dnt93Au+5yzzzn7stbMfs95z9ln7+f5Pnvv9az1rPVdsVHEFLPv2bOnlq6CHSWr4hWdV3xR4cBFClMPx12xYkXRw1LIfNy4cVrHuHAH1ogRHmZZFvAhSIc1xnvuuSfLUMS27+KfdFOnTlUKyptvvjk/aYIZ7lvf+lYkhTiGDx+uinzZsmVlJ7YEf0G7GEYIGkNMESsMFf8Q3Is1TJzRDTfcUHH/OO8QG0XM7BRCdYQIVIIqsiheEUs+YeTdHzp0qB5z0aJF+WPDHETQCgXNJWVCa6tyYReKV8RZrrpE3Vl4bMEHjnQ/VoVY2f+tQ4DsC2p7T5w4URUwVIe4oBm7KN29Rx55pNKYcp/6iS380FwX119/vRs4cKAWHKAt0GNWEllrdhB58Bw0RVwJrfe/x7uBnoBMx2ebhPtlDPeK0+o2wQqSTqAL7yTYZ1HI9ZXLJCeR5JF3X2bnGhQnCiQnLrXckCFDlDiA8/lN1r6KnlceZroPQV9ZFOo7Q7QAThItm0UIYt3nLVu25CQ7IH8dS/ZBbuHChQ1pszz8NWdcrGI9H+ciiNHfQ8FXWRIK1Ybnn39efy/8+6H2z/pOMmlRvPbff/+crNUnHo7KUzW5qpolWFvXXnutphN885vf1FKJWSM997NnZshRC7WJqRjDq6cXLTwHM/1iIld6sY8z8RkpSpSUZExGjRqVmBqnWRgc8oF5ZgjTlVpHPENIYcEFXcyzUwsmjD/rt1RdYiOHOBjDgVeplBCEGkbgKUdOP/30MLtneh+eRT6mCC9qUkodlh20uE0loHmE3UYanZO14rg1r+HtIb2IvgvRSSTnktq4uTvvvDOHpctxK22lLPGsWsR4aTytIexvst4XybjYQepDQGIdclIVKdeuXTu9pvfbb7+cTJJyUTLS+Rbu3LlTU5Uq3TvFvpdKav4wZV/POOMM7cfixYvL7mdf5tTTAdbCeZ+DEjQN0nKKy2IgzpkzRy/Kzp07R5qnV+xccftM1m+177KeW1fTnnzyydygQYNyPKCKPSBKfSZBD0XPKyT4ehyUUpZk5MiR2m/xFOQkDSZLXY9lX8nblRTHnKwJ5q9rllgkjaWh7YV6stp7ScqahmqTROLnhIwiB8c1+fompREAq06dOunYN4rvu/TZG/dNbIK1RDHkBaYtqOeEu9fde++9+c+z8MYHqdXrkoc2VDhzNbS/GtxK5VX6IC0ftFXNMZO6L5V4ZFKojD2PPvqos6pTrRtJeQRqIFSPHj2UVpTcYGo+b9261RG53L1794Y2TixWR13iaiSsW5ogLSKrqdFu9LHlEf7Od77jCIoDK8pSpkViqYhl5qnlyACZtWJPJpEW0Mv1IypFzEOKda1qowlLrRH7tCUK32dBiH6FIIAJEbmi3PgmrUGAtB5ygYn6h4oSYg7SVtasWePE6mxao1iL9mVFw5yUlMwwIta27jZgwIAwu2d2H7jdPac05UbRE2mRWCpiwIUJisLYBBZhmWRF3nvvPe2qD9qqp99YD6RuiIs/9GEqKeKwOZGhTxjDHXnwE8yGFSbr6w4yFJPmI4BHB+VEABN8wngkYDJ75ZVX3AUXXFAXIUctvSG16P777w+l/Ll/wwZerVy5UpvTv3//WpqVmd9Q0IF0MXRDWGwTA07jvN71H1miEzXdhvWTsEEP9Z+1tUfwHMasTUYlrG0KWUp+TU0uzpLvZeJT9LRw8vI7AunSLKQpfehDH9K+Xn311Wnuamz7Joo2N3jw4Pw1Cv8yAYesD8dBhEgizwld6l6CxzqMUPOYtWdR3DkC0EyKIyBR5Xo9EJxHGmbaJLYWMTMZkvEvueQSXT/x4ep8nmbx9JPy8Imsm1gSuKlh1ionpHuUSgWQfD39KekiaRXWG1nXg8aSilW4v0yah4AEKik5Awx7rMnLQ1fJMSQSWgk6ovASRdEbyG+wzMtJNW5pKBr79u3rorzny7Utad+xHIY1jMCgdfTRRyetC5XbG/eZBdbcAQccoLMhuTnj3ty62ydrINpXWZ+s+1iFB6Bgtri/8paGXB1t3hONWErkYaHeCSI7ITRIm5CiAjEDmMhD1NKUmjjA27dvzwlvdz4qGQ+YUInmJD+3ia2o/lRSPrHN/RO8nyDoCCOkXPE7YRMMs3sm9xGaUsUIshPSW9MosUxfKgRaLBMdiI4dO6Ymb6ywj/5/FDA35pQpU/xHkb6SB3veeefpOYIPDt5XcqdR8o/93nrrrUjb1OqD4XYn1cRjkJbcxFbjWun8khWRE7L+PHBA5YcAACLlSURBVGMZ7tnRo0fXVJqw0rka8T33UrCsor+fpHBLTmI9Qp1SrDu97oTnPdT+WdtJgvPUbS/r8znx6qW2+7F2TcuFrTJ27FjlQYZAnSjqNAuRgYivwhR1X4l+JuJULJB9Dl0qdcnv6Ll6cSGmRUiD8cFABLd5l2ha+hfHfuBqZsmpS5cump5IgCLVcxiLWbNm1VSasBX95F4SAo59OKxJdQoT0UsJP6ll7Tp06KD82K3oQ5zPSQbJiBEjnL8+cN+nVRKhiLmoH3jgASVP/973vqdpC2kdELE2tWuNUsQcnLU24bR2THCCUipi2u+TNkVMNZ5+/fppRC4FL4gwJ//apDEIEF8wadIkLVe3YMECVVZUMSIliXztUvEJjWlNNEdlvZK+BPP+w+YP+1rGrCcTkW3SFoEZM2ZoIQxSMKmmlWZJhCJmACjEDa8ogQ0kcqc1jQZeW6TRDyVu/OnTp7vPfOYzej7+VFLEzNyRcty6ukMC/kAWwwybAC2UMKQKjZz8JACShjVRMh7UsoGsgspqXHso4M2bN4dOB2pY4yI4MJWWZB0zfyQs/TDy+OOP62783qQtApB2MGlDZs6c6WRZrO0OafsvSU531mQ8Z/KECROS1PTQbZU6t7pmRApHM0RI7PV8cl3nfvzjH5c9pUQu6r4/+MEPyu4X9y9fffXVPEUi/NGyHBD3JieyfZs2bcpdfPHF+SAs0nSENS8nk59E9qdcowlm9MF+kmtcblf9jlQlT2uZhupBFTtcxQ5g6fnduX6yIIkI1goOxPr16/XGJnpXEuGDXyX+PcoXhZiGTVxJsRwPqVebO+SQQxRjSjoyuTOJFgF4n4MlCblXiYpuREGGeloudYNTca+J9yxHbnNahGcHz0D4xMmzzoIkThEzKFLiTAeKCzBNA5UmRcyNFDdlLIxZ+Wo9KAoI5E2iQ0DiODT1y08kIUa56qqrchKQFN1JIjxSWhQxePMshIwm6fKrX/0qJ0Fw+nxfsWJF0rsTuv2JVMSkBsjapg6W1IlNTV6rV8S435MsQgyilWR4QNx9992x6IoQMORdpJSYxP1lEg0Cy5cvzwnJhd6PjDnF2idOnBj7PGCviG+//fZogGjBUSQDIu8ST7oyFjKjnGQu6HV0xRVXtADN1p0yMcFacoPnhSjqhQsXOnExauUVKnKYxAcBorLvuecejSQlwG7atGktaxwpEHJTKzMP7wmqkTrXbaJcW9a4BJ+YKlxE3hNEKWUInVhjToh3HPci6W1wdPvgvgR3M/ZN516jwItQgjoCPeGr3rhxY+zbXayBPCvEGHFC3JH6KOl9+t+6OUD9Z166dKmyPUmJvhxrx0mXtFjEcANDiiElLFtqGUu6TI5atXLRazvmzZuX9Euk5e0HUynCkg92A1v4fydPnpy49fY0WMQ///nPlaCHWAfPz51Ey1jysfU+5VlOkF/WJJEWsZ9NiFtaU5qEDN4NGjRICQH8d/baWgRIUSHNTIp35y3jZnouIH/BUhC3qdZ4XbVqVVESk9ailJyzk98uRTAUS3E7670GJ/T8+fPdn//8Z3fLLbc4XyozOb1KT0vBntKdZ511Vt4yFoWWiA6SqgTBC8IzAmKdrEmiFTGDJVSQ6o4ht5VydeQZm8QHARiTvDKWlDO90WptnSc7qfR72LFOPPFEJ2lKWjAeV93nP//5Sj+z74sgIOltWhLysMMOc5K25pj0whwlgTTqjh42bJgSxBT5qX3UZARQxmIhJ0oZU9CBkpYUWoHYRIL7moxaPE6XeEXMGsmiRYuUrWft2rVKHCBujXiga61QBAqVcS1VjYREXydcrPOWEiZhN910kzv77LPVSqNuKTVtqT5lUh0Cjz32mLKOUQFN8suVhENyOp1EtToK2Z955pnGBlUdpE3ZO6iMJT9Z75k4W8bjx49369at05rpxP1kVtLiiyfsvX379rrOkNQoyDStERerrUqtZXJK5WbLiQsq9KVHiponx3/ooYeK/o5KPRSt4NjiFs+x/idKu+i+9mFxBFj//f73v58TZijFESxJQRJ+9xzfpU3StEZcODasGctynY6jeDNy4hUq3KXl/0uwn7YPYpOsF71IZPpSqSvokUce0Qc9D+JKLFGljtHKz9OuiMG2WmWMMhULN68YSO0qTD2iPCYPGxSHWL+5/5Ui4ibhERAKSi076MuNgiPKWHjdc6SUpFXSrIgZs0JlHKcgKPFUKbMY15pYwmm9xEL3K1WKmF7PnTtXH8jksooLLTQQcdgxC4oYnGXNOG8ZU+KynGA5c7MGN0k/0p+QT37rrbeqBcz31Fqm5rJJZQSY4AjXsU5ymLh6fGEbW7JkyT6TncpHTN4eaVfEjAikNUHLOA7KGEpZqfSm11zW8oVL3SWpU8R09Otf/7oOMjN8IZYv1ffYfZ4VRQzwYZSxrPnnmFB5JeFfKRC+Zs2aHK98xj6weBVayrEb4Bg0CAsXZjpq5no8cT9TF5jlnSxJFhQx44kyljV9HW88R61UxrQFfneuvd69e4eu25z26/L/CyCpE2FzclR8EVe1Rss+/fTTGj2buo7GuEOyTu9EQZZtIUQQy5Ytc1/72tdc165dlZTA/4CUmAsvvNBBHFEo1HGVm5lJpNau5RhESZuURkC4iLWKjRQkcESqIlTbIkp11KhRmpZU+tf2TZwR4H4gfayYcI8gVBjbtm2bblRcI2WIaljNFuoLC9+769SpkxMeiFB1m5vdxlacL/FR08VAI4eVSGqfU0d90DQVsy/W5zh9dsMNN1RUwrT35JNPdjIr1qYTFR2UkSNHuu3btwc/avOeB4x4PrSMoSnhNtDk/yHCnPtAgth0okIaGRObAQMGaEF76jFfd911LXkg5xtpb+pCABYqoWzVSSn3RHDjwDwL2YQoQ9PQhH7UCdmOe/PNN+s6by0/Fk+Mk2BLJ0G1OgEnJc7kfQRSaRHTNS48cuoo/C5uTEcahgTxuGOOOcbGvsEIiMsv9BkYJ8YnKNBjQk5QSVDkEnFZabfMfY/lM2vWLKWg9LnX4n7WiQuTF6hhTdKBQPfu3ZW2NWxvmJhRB7rZAq2sRN+rBcx7yGBM/otAKi1i3z1y6p588km1CGBaQhmXs7L87+y1dQhIGoOypYVpATNsI3B5HynczVgbkG1069ZN+b1RwrgkJVLdSYCMu+2220wJh7mwbJ9IEYDVDpc01rpE4is3eaQnSMHBUq2IGR9cMbAAffazn3U8mPr27eu2bNmSgqFrfhe8myv4GmUrJFdV3Wx+DbPSsRlHZtdZFtYHsXJZ773ooovU6yO52uqGfPHFF5UsATc/k1KTZCLA/ZZUgQDmvPPO07gErtOsMmdVGr/UK2IAQBlLqoYGbu3cuVPXJePMNlNp0IKKMPg++Lvg58H3wX1qeR9cg+J9lDJ69GgNJqnmmFm0ioUsxUlBDZ1cUv0IpjIhPXG9evVyUu5RWcUWLFigHqBqsLR990UgeO8E3wf3DH4efB/cp9b3HC+psmvXLq0BwASbwEvJbEhqVxre7kwoYlBEGUPbN3DgQLd3714n+ZLuqaeeajjAjTiBV4Yc278vVIr+/+D3vI/rjY1LGtdqtfLaa68pBWO1v0va/owd0f9YvQceeKDDyiX6lPeXX365Wr5YH7wnGMYkGgT8/cPR/Hteg+L/D37P+7jea8G2N+r97t279RlLDQCetRQHwVNjUhyBTCEj5fk0pYloatJjIBlHOWdJCh8QwRl8sfeNxsanJ+GxqFWwioXco9afx/p31JhlbZd1XwpXMFkhGpr3pCJhdWAFsxZsEi8E6r3XuB85RtJEctUdmSpvvPGGVkCjCAs1AUxKI5ApRQwMrJWRd8paBUqAUopZdplwo5fbCi+doLIu/K6W/1nHDAqz5v32209vXCKq2chHZtzYiJImApiNiZXUwnXMvsNEWQfPE+f3Qk2oSnbo0KGuQ4cOWsiCSOgjjjjCkRq2detWt3r1ajd8+HDFIM59sbb9F4Fy91mhwk2qEsbAIQ5HyInccccd55544gn1Rv4XBXtXDIFMTlN42FPSjdD/sWPHOuqr4uKcPXu2uU+KXSX/+azYw4KvCj8vc4h9voJcAK8ElX0efPDBfb7PygdYub/4xS80v3fevHl5C5+JB+U9qe0sFJ5ZgcP6mUAEfF14SmdC0EO66EEHHZTAnjS/yZlUxB7mK6+80h1++OFqWcyZM0eVMXl2WUo0Z+ZdTsopWb6r9Ptyx+Y7rF+kU6dO+pq1P5SAo8wgZBsEYSFgyroawWsEuXiMsoZN2vpb6V7x9xr7+fdJwYCALOG0Vk4A7mVKZVr50fCjl2lFDEzUrJVKMxrERTAMTDVYJpBFZEGSdsOnYUxY95Vyg0o5SX6vFyHCd5dddpkbNmyYEc94UFL0Ws29Vqi0/f/VHKNZ0KGEmThu2LBB68KzbHLUUUc16/SpOE/mFTGjCEUiXLysyT333HOaFkJKyJgxY1IxyMFORDnbjvJYwTam8T1rZyjf5cuXa4Sz7yOuuy9/+cs6EcSiMEkPArXeH4XKttbjNANJuBmkoIQqYSxh3NEdO3ZsxqlTdQ5TxP8ZTlzUMMBMmDDBQbEo5bnUvUI+ZtxoFLkxvQTfB29g/7l/9fsH9/GfhX2N8lhhz5nk/UiTmz59ugaswKXtsSfI7Itf/KIGCkJ2UIhrkvuctrYHxyb43o8l/fWf+1ePQXAf/1k1r8Hj+ff1HrOa81faF89Ov379lDITQhk8imYJV0Kt+PemiAO4EJ1LEFefPn20Ig2RuOvXr9cgIoKK4iJhbsYw+1Tbn0Ycs9o2xH1/LF8mbzyUFi9enFe+RHyfc845GniF58XWfeM+ku+3L8w1H2afWnrbqOPW0pbC30AZTPonQa4s561cuVKj+gv3s//DIWCKuAhOVDMhL/OCCy5wkCRAjzl+/Hh31113WT5cEbyy/hGVvVC+uOVIjfMPUCZ2gwcPVusX5Wu5lFm/UtLRf2hVyRP+/e9/r25oiJEOPfTQdHSuRb0wRVwC+GOPPda98MILbty4cUqYwJoxD1oiXJkBmmQbAUoIonyJDmXzguULYQyTOQIBTfl6ZOw1DQgQQ3P22Wcrbz+kMtR8h93NpD4EMkfoUQ1ckEnMnDlTFTDlE7GOTzrpJHfrrbcWLVhf7NhY0VYhqBgyyfoMK3ft2rVu0qRJyuF89NFHO8o9ooQhF2HNl7KbpCDxiiI2JZysMbbWlkfg4YcfVkuYAK0hQ4aoO9qUcHnMwn5rijgEUj4ggQAuKgNRV5MIQbh+ywnFt3lYk45iyrgcUvH8DoYr6PnI54XVijiBO+64w0FYQLTzpZdeqmQkFFpfuHChVpmxtd94jqW1qj4E8AiyvPLOO+9oNsmSJUtsolkfpG1+ba7pNnCU/ofI6RkzZqjlQ64nQQqnnXaaPox5OBcjAfnpT3+qChh+YFyWc+fONeau0hDH4htI6kkxYqOWNQ8eL+T5wnJ17rnnat6kWbweGXtNKwJwuEN8BOEMkdvwut94441p7W7L+mWKuEroWRchWIEC17gp77vvPofLZvLkycpfjTvbS7CaEAT9BO9Ao+lTEfx+9to6BHA5Q+DC2hfWL1WgfLAVrYLYBcULJzn55iaGQFYQgGyGJRbuD5ZfoF5lCcYkegRMEdeAKVzVFLnGOho1apSulfA/68F33323FmXfsWPHPq7rH/3oR2oZk1tqyrgG4CP6CdVhmBiR28skylNLcnjGlmArUo3I8f3oRz8a0VntMIZAchDA4wdRx/bt27WsJulJvXv3Tk4HEtZSU8R1DBjrxMwWV6xY4a655hpNbGc9mFzk448/vo1l5U8DWQhuahS2SXMQYH0epUsQFZYv/M5Bq5da1V/60pc0GnTAgAFa8ak5LbOzGALxQ4D1X9jeoK7EI8T/RtTR2HEyRRwBvlATDhw4UNOcKFNHdC1bKZk2bZoq49tvv73ULvZ5nQgwk6eqExMl0s54qHhhbfeMM87QCFBczpRrMzEEso4Ak1OeTVSj4z3BWaRrBpfbso5Ro/pvijgiZHFpEtRAuTrWjuEVLidTpkxRZczaskn9CBDNTk4j3gnWeqHfC8onPvEJZQJCAeN6tujmIDr2PusI/OEPf1CvEC5oJqoEoEL3a9IcBEwRR4wzQQ24OsMIaVC4qa+99towu9s+AQRY10XprlmzRi1fAuiCAtMPbmYYgPBWkH5kYggYAvsiwLINecEEZ5EXDDsc1ZRMmoeAKeIGYP2Tn/wk9FGvu+46VcZQaJqURmDPnj1alAM6PfJ4UcBBYRZPgBV0pP3793cnnHCCBcQFAbL3hkABAsRO3HbbbZqSxHsyQnh2WR3hAqCa8K8p4ohBhhbz9ddfr+qoRFyT2tRPKpmYvI8Aa7zPPvuse+aZZ7SAwpYtW9pAA159+/ZVdzP53JBtsDyQZKFgBFYJG2vabER4s+EBePfdd3X75z//qWt4PuAMNzsbmLCeh1emffv2WjUM7wwbn1FhDE+BRewn+SqJpu1cWzxvuMe4b6ZOnaquaLs2osG32qOYIq4WsQr7V2MN+0PxQGV9GcWSVdm9e7eWDCTIjZSiIJEGmKBYsHbBCIuXVIqkPDSwNjZt2uRgWtu6davbuXOnrmH/9re/dRCIQBkIgX6zBKVMFCwkNB/5yEeUJY4ydlhCXbp00f/NKmrWaLTmPFROevvtt3XJZv78+Rq82JqW2FlBwBRxhNfBv//9bwebVq1C6TzEWzq1HieJv4MYJSgwmVH9qlevXrpeRTWsOAdYMfbkXuIRoSAE7ymhyQSD6kxhhPW5Qw45ROkzDzjgAF2vo3YxW7t27dTaxeLF8mUS4jcUPQxIWMwErUHNCe0mr1g+WNj8zyttYZJD+9hKCechmpziJ127dtVCJwS89ezZU9tS6nf2ebwR4FpBUMKsC88Tkg4mYyatRcAUcYT4s365a9euikfkIdtJcpApJOE3Ho6rVq1SWsWkWHoVO1rFDjwMqOpCVPMpp5ziunfvXsWvm7srSu/xxx9X65ZCIChfXOelJlBMIFBqWKGMO69YnFih3l3M/80Yd9qICxzrnEhZlDNLKVjkRJq/+uqrDjIa3mPFsxUKfWANHqYxNiZJFMEwiT8CfjJrVJXxGitTxBGOh3dLEwndsWPHvJL1yta/lpqBktPao0ePCFuUnENRWCGOedVYuhs3blSWNDwWWLnFlC5KtHPnzg6rESsSy7Fbt2555RuXkaCdBx98sG7lJjusQ7/yyitu27ZtqpiZcBCZ7hU1ynrp0qX5bsFAxjo926mnnqqBP/6hn9/J3sQGAVgBTeKDgCniCMdixIgRyjmNpdMM6ybCptuh/oMAFiOlDYnMRvHi5QhSYLIbCgYrEGX76U9/Wq3CtLlscX9T8pMtKODDRGTDhg2qmKlARsAPLm9fLIP9CQACIzwcBNORRsZyg4khYAjsi4Ap4n0xqfkTLAGTxiCAO/iNN95QKzPqMxBABY0fCoWKS6ytBgVXMkFijC8KlxzLpEdoB/tXzXsmmFjShdY0AWjgR1oZ7myWWV588UXdOD6/I0qXgDvyusHSqldVg7ztm2YETBGneXRT1DdYflDEkA3UKyh1FAXc0w8++KAGMgWPyXpunz59NDobpUEAlUl5BCgReeGFF+rGngSO4VFAGUM1CmkEVKNsLEEQJ4GVTMDQ4MGDNa2q/BnsW0MgvQiYIk7v2KamZzzIKaSBW5QHOlZVtUIuLkp84cKFurbpo0c5DhHJBIph6VK0IywzWrVtyNL+RF2TZsYGaQ1r7Ux+GMsnnnhCC6T87Gc/c2x4F7CQKTfJ2iVr7SaGQJYQMEWcpdFOYF/Js2XtHSWMXH/99W716tX6vtIflC8uZ6xeopxRBl5IiyI4Dncp65i2pu+RacwryhbKUTaE7AK4wVlXRkFTFYvtG9/4hlb8IXWNWrgEOJoYAmlHwBRx2kc4wf1D+aKESbPxQvAU67gQEhQTn1p07733quXrlS+KgHVeLC7qDFu6TTH0mvcZ3N9jxozRjdQ9CnU89NBDOnF66aWXHBuTLohbKFGJ29tqQzdvfOxMzUXgf5p7OjubIRAeASpYYckWCqUmvYXsv2P9eNSoUepmxtLFEmYf1nqpAY0FhhIfN26cKWEPWkxeYfqi5N7ixYuVzhNLmf+ZPBH8NXbsWGWAwnXNuLL+bGIIpAkBU8RpGs0U9YW8VdYWi8m6des00Iro5tmzZ+vaItSMc+bMceS/EtFLGTcoJOGqhj60Q4cOxQ5ln8UMAdKmWK9ftGiRBtGRm3/mmWcqcxgKGm8GgWFULCPH2cQQSAMC5ppOwyimrA+s7V588cVK11iqa5dffrnS9EHpiEAPSaDVpZdeqsxcpX5nnycHAfKOfSQ2ecoPPPCAY8lh8+bN7q677tKN9f2rrrpKlbcRiCRnbK2lbREwi7gtHvZfDBCgGlVhfeHCZvFgRglTug2XJoUTcEFDj2mSPgRYH+a6gO0LEpGvfOUrugxBoBfr/kzEJk2atE8qWvqQsB6lEQFTxGkc1QT3iTVA3M1hhIArArfOP//8zBJshMEpbftAo4llDD82cQTkfcN+xnIE9LE33XST27t3b9q6bf1JMQKmiFM8uEnrGg/WkSNHhm42a8A//OEPQ+9vO6YLAazgq6++WvmvyU0eNGiQBnJR7B4SFiZoxBOYGAJxR8AUcdxHKCPtI8J5+PDhbs+ePVX1GJamQi7oqg5gO6cCAdLZIAthyYIgPyKuYU771Kc+paQiMHqZGAJxRcAUcVxHJmPtgqSD9b5qBcIPWLdMDAEQwBL+9re/rbnnBHHBZ811RVAX7F0rV640oAyB2CFgijh2Q5KtBlG4HiEAp1YhgpYauyaGgEeAUo9M0LCQcVUT7EVOMvzWlKokrc3EEIgLApa+FJeRyGg7goQd5JC2b9/e8coGX3G5LbjPCy+8UJJtK6PQWrcFAdaRIYAZP368KuapU6dqRD4sa8ZpbZdIXBAwRRyXkchoO6jn+/LLL2veMJzQJoZAIxCgsAdrx6RATZs2zd1yyy1azYtzFbK0NeL8dkxDoBwCpojLoWPfNRwByuQhMCh97GMfczfeeKM76KCDGn5eO0GyEdixY4dauOSSw6bGa+H7Yv/zGcsYvJoYAnFBwBRxXEYio+2g3i90llQ/wlIhL5ToaXJBO3XqlFFUrNuVECCHnMCrDRs2VNq15PeHHnqo++Mf/2iVt0oiZF80CwEL1moW0naeoghQIB6BshJ6StyEc+fOdV27dnWXXHKJMikV/aF9mGkESE+aMmVKzRhQ05oiISaGQBwQMEUch1GwNqg7+r777nO/+93v1CKmfOGCBQtcz549lUf42WefNZQMgTYIDB482H3uc59r81mYfygA8vDDD4fZ1fYxBJqCgCnipsBsJwmLAOvE999/vwbSkAdKibxly5bpA5eqOzNmzHC+xnDYY9p+6UWAnOFqhCUQ6DEPO+ywan5m+xoCDUXAFHFD4bWD14pAx44dNRjnN7/5jZs8ebJyCEOBiXI+/PDDVSETpGOSbQROO+00N2TIkNAgXHHFFQ5L2sQQiBMCpojjNBrWln0QIKCGVJM333xT6w1jFRNgg0JmfZl15HrIQPY5oX2QGARef/11d80117innnoqVJsh8iCP2MQQiBsCpojjNiLWnqIIUGv2q1/9qqPQAwT/I0aM0LQV1pEJvGEt+bvf/a4yKRU9gH2YCgTeeecdR7756aefrgF9RNq//fbbypxVroOQv/A7ljpMDIG4IWCKOG4jYu0piwBrfBD8E1m9e/durUGLq5o6tZA1HHnkkVp155FHHtFKPGUPZl8mAgEi6Z977jk3evRod8QRR7hhw4a51cJNTuQ071977TVHPvoHPvCBkv2hOMiJJ55Y8nv7whBoJQKmiFuJvp27LgSIfiWFhUhromBR0O+++65W3Tn33HOVIGTMmDHu6aeftgCvupBuzY9//etfK8FLly5dXJ8+fbTkJWQcvXv31prVvJ8/f75axuSco6iLCQUfmKSZGAJxRcAIPeI6Mtau0AhQYQfFy8ZaMu7qefPmuU2bNrlZs2bphiU1dOhQ3QjwwZoyiR8CeDYWL17sFi1a1CaHHE8H1i9LEscdd1zRhsMpTQocrmovVGMiCt/G2yNir3FEwJ5GcRwVa1PNCGAlT5gwwW3cuFFZlyZNmuSwlnbt2uWmT5/uIPsnRQoyh+XLl7u///3vNZ/Lflg/Av/617/U7Txx4kS1bFnrv/nmm1UJo0Qvu+wyRy1hYgPuuOOOkkqYlpCSVGj5zp492xHgZ2IIxBkBs4jjPDrWtroQOP744x0b64Pr1q3LW1rbt2/XCOw5c+a4D37wg1o4/qyzztIi8hSSN2ksAiwfLF26VJcTHn300TYlLClXiGcD7wUuZapvVSNMwmbOnKmR9TC1cRwTQyDuCJgijvsIWfvqRoAALxQs25133qmWMoqAjSAflAEbgiI4//zzHVWhLrroIouyrRt9p+v2q1atcpSqXLJkiVu/fn2bikesAZ9zzjm6sWxAhHytsv/++2vZQ2oR4wExMQSSgIAp4iSMkrUxUgROOOEEx0alpz/96U+OCOvHHntMlQSF5HFnspEu1aNHD1XKWGcoiWOOOSbStqTxYH/5y1/UnUykM+v0pJsF2dBQtIMGDcp7Irp16xYpDFdeeaUysaGUTQyBJCBgijgJo2RtbBgCBx98sMOFyUaaDC7sFStWOLitSZEheIiNICCEEo2sM2MxEzQEqxPu7awKbua1a9eqwn3mmWfU6oVoI1jjl0ApIp379++v+b/9+vVzBNg1Sjj2ySef3KjD23ENgcgRaNzdEHlT7YCGQGMRCLqwORPKBJfqL3/5Sw0oQtFgQXu3tm8NEb0oZVzfBBsRHHbqqafW5WL1x47LK1hs2bLFsb6OO59gOOhHfT3pYDshzaD/eBAgW/nCF76QKiyCfbX3hkAUCJgijgJFO0YqEUAxY8WxecHaY63z+eef17Vm6DXhwGZDaQeFmrko5Y9//OOOdVCid4899ljH50R3x03IyyU6mY00MJQt/eU9noKglevbjpsZooxPfvKT7pRTTlFPAQFyli7kEbJXQ6AyAqaIK2NkexgCeQRQqGwEcnl5+eWX3bZt29Q9+9JLL6nygozCKzUIRQoFRYVCJpWKtBvYweDVJmWHDRf4gQce6FjnbNeunWvfvr0yR8Eexcbv2VCErL+ykQpEIYx//OMfmpYFHeRf//pXx5otryhaeLr37t2rr0wesPBJ7cLSfe+99wqbmf+fSUnnzp3V8vfR6HA3n3TSSY7vTAwBQ6B2BEwR146d/dIQUAR69erl2Ii29oL1iIKGhhOrkqAlLEuUMy7dPXv2uB07dujmf9PqV/iYu3btqhME725n0sGEAYVrVm6rR8jOn1YETBGndWStXy1FACsRl20pfmOsVxQy1uhbb72lrm0itrFa/SuWLCxRWLN/+9vf1MrF2mULWsFB6xhlisXMOi1BZFjTWNcf/vCH9ZXgNN4fddRRaoFjdRMJjmXOexNDwBBoPgKmiJuPecUzbt682dx9FVFK9g64lHH1spm0FgGoMdlMDIFWIWAUl61C3s5rCBgChoAhYAgIAv9P1rJyhoQhYAgYAoaAIWAItAYBs4hbg7ud1RAwBAwBQ8AQUARMEduFYAgYAoaAIWAItBABU8QtBN9ObQgYAoaAIWAImCK2a8AQMAQMAUPAEGghAqaIWwi+ndoQMAQMAUPAEDBFbNeAIWAIGAKGgCHQQgRMEbcQfDu1IWAIGAKGgCFgitiuAUPAEDAEDAFDoIUImCJuIfh2akPAEDAEDAFDwBSxXQOGgCFgCBgChkALETBF3ELw7dSGgCFgCBgChoApYrsGDAFDwBAwBAyBFiJgiriF4NupDQFDwBAwBAyB/wNFmsIjlP6+nAAAAABJRU5ErkJggg==" } }, "cell_type": "markdown", "metadata": {}, "source": [ "# Googleのページランク\n", "\n", ">多くの良質なページからリンクされているページはやはり良質なページである\n", "\n", "Googleのpage rankは上のような非常に単純な仮定から成り立っている．ページランクを実際に求めよう．つぎのようなリンクが張られたページ(有向グラフ)を考える．\n", "\n", "![image.png](attachment:image.png)\n", "\n", "\n", "計算手順は以下の通り\n", "1. リンクを再現する**隣接行列**を作る．ページに番号をつけて，その間が結ばれているi-j要素を1，そうでない要素を0とする．\n", "1. 隣接行列を**転置**する\n", "1. 列ベクトルの総和が1となるように**規格化**する．\n", "1. こうして得られた**推移確率行列**の最大固有値に属する固有ベクトルを求め，適当に規格化する．\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 課題 \n", "\n", "1. 上記手順を参考にして，pythonでページランクを求めよ．\n", "1. このような問題ではすべての固有値・固有ベクトルを求める必要はなく，最大の固有値を示す固有ベクトルを求めるだけでよい．初期ベクトルを適当に決めて，何度も推移確率行列を掛ける反復法でページランクを求めよ．\n", "
\n", "
隣接行列
\n", "
\n", "\n", "$$\n", "{\\it A1}\\, := \\, \\left[ \\begin {array}{c|c|c|c|c|c|c|c} \n", "&1&2&3&4&5&6&7\\\\\n", "1&0&1&1&1&1&0&1\\\\\n", "2&1&0&0&0&0&0&0\\\\\n", "3& & & & & & & \\\\\n", "4& & & & & & & \\\\\n", "5& & & & & & & \\\\\n", "6& & & & & & & \\\\\n", "7& & & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
転置行列
\n", "
\n", "\n", "$$\n", "{Transpose}({\\it A1})\\, := \\, \\left[ \\begin {array}{c|c|c|c|c|c|c} \n", "0 &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "1& & & & & & \\\\\n", "1& & & & & & \\\\\n", "1& & & & & & \\\\\n", "1& & & & & & \\\\\n", "0& & & & & & \\\\\n", "1& & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
規格化
\n", "
\n", "\n", "$$\n", "\\left[ \\begin {array}{c|c|c|c|c|c|c} \n", "\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
遷移
\n", "
\n", "\n", "$$\n", "\\left( \\begin {array}{ccccccc} \n", "0 &1 &1/2 &0 &1/4 &1/2 &0 \\\\\n", "1/5 &0 &1/2 &1/3 &0 &0 &0 \\\\\n", "1/5 &0 &0 &1/3 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &0 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &1/3 &0 &1/2 &1 \\\\\n", "0 &0 &0 &0 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &0 &0 &0 &0 \n", "\\end {array} \\right) \n", "\\left( \\begin {array}{c} \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7 \n", "\\end {array} \\right) \\, = \\, \n", "\\left( \\begin {array}{ccccccc} \n", "\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \n", "\\end {array} \\right)\n", "\\, = \\, \\left( \\begin {array}{c} \n", "0.32\\\\ \n", "0.15\\\\ \n", "0.11\\\\ \n", "0.06\\\\ \n", "0.29\\\\ \n", "0.04\\\\ \n", "0.03 \n", "\\end {array} \\right) \n", "$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 1. 0.5 0. 0.25 0.5\n", " 0. ]\n", " [0.2 0. 0.5 0.33333333 0. 0.\n", " 0. ]\n", " [0.2 0. 0. 0.33333333 0.25 0.\n", " 0. ]\n", " [0.2 0. 0. 0. 0.25 0.\n", " 0. ]\n", " [0.2 0. 0. 0.33333333 0. 0.5\n", " 1. ]\n", " [0. 0. 0. 0. 0.25 0.\n", " 0. ]\n", " [0.2 0. 0. 0. 0. 0.\n", " 0. ]]\n", "array([0.29404762, 0.14166667, 0.15833333, 0.13690476, 0.13214286,\n", " 0.07261905, 0.06428571])\n", "array([0.69945653+0.j, 0.38286042+0.j, 0.32395882+0.j, 0.24296911+0.j,\n", " 0.41231122+0.j, 0.1030778 +0.j, 0.13989131+0.j])\n" ] } ], "source": [ "from pprint import pprint\n", "from numpy import array, zeros, diagflat, dot, transpose\n", "from scipy.linalg import eig\n", "\n", "A = array([[0,1,1,1,1,0,1],\n", " [1,0,0,0,0,0,0],\n", " [1,1,0,0,0,0,0],\n", " [0,1,1,0,1,0,0],\n", " [1,0,1,1,0,1,0],\n", " [1,0,0,0,1,0,0],\n", " [0,0,0,0,1,0,0]])\n", "\n", "diag = []\n", "for i in range(0,7):\n", " tmp = 0.0\n", " for j in range(0,7):\n", " tmp += A[i,j]\n", " diag.append(1.0/tmp)\n", "\n", "D = diagflat(diag)\n", "tA = dot(transpose(A),D)\n", "print(tA)\n", "x = array([1/7,1/7,1/7,1/7,1/7,1/7,1/7])\n", "pprint(dot(tA,dot(tA,x)))\n", "\n", "l, V = eig(tA)\n", "v0 = V[:,0]\n", "pprint(v0)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 0.200, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.500, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.333, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.000, 0.500, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000]])\n", "array([[ 0.000, 1.000, 0.500, 0.000, 0.250, 0.500, 0.000],\n", " [ 0.200, 0.000, 0.500, 0.333, 0.000, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.333, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.333, 0.000, 0.500, 1.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000]])\n" ] } ], "source": [ "import numpy as np\n", "from pprint import pprint\n", "\n", "np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n", "#pprint(A)\n", "pprint(D)\n", "#pprint(transpose(A))\n", "pprint(dot(transpose(A),D))\n", "x = array([1,0,0,0,0,0,0])\n", "#pprint(x)\n", "#pprint(dot(tA,dot(tA,dot(tA,x))))\n", "l, V = eig(tA)\n", "v0 = V[:,0]\n", "#pprint(v0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 累乗(べき乗)法により最大固有値が求まる原理\n", "\n", "累乗(べき乗)法は，最大固有値とその固有ベクトルを効率的に見つける算法である．すこし，固有値について復習しておく．正方行列$A$に対して，\n", "$$\n", "A x = \\lambda x\n", "$$\n", "の解$\\lambda$を固有値，$x$を固有ベクトルという．$\\lambda$は，\n", "$$\n", "\\det( A - \\lambda E) =0\n", "$$\n", "として求まる永年方程式の解である．\n", "\n", "では，なぜ適当な初期ベクトル$x_0$から始めて，反復\n", "$$\n", "x_{k+1} = A x_k\n", "$$\n", "を繰り返すと，$A$の絶対値最大の固有値に属する固有ベクトルに近づいていくのかを見ておこう．\n", "\n", "すべての固有値がお互いに異なる場合を考える．今，行列の固有値を絶対値の大きなもの順に並べて，$|\\lambda_1|>|\\lambda_2|>\\cdots>|\\lambda_n|$とし，対応する長さを1に規格化した固有ベクトルを$x_1, x_2, \\ldots, x_n$とする．初期ベクトルは固有ベクトルの線形結合で表わせて，\n", "$$\n", "X_0 = c_1x_1+c_2x_2+\\cdots+c_nx_n\n", "$$\n", "となるとする．これに行列$A$を$N$回掛けると，\n", "$$\n", "A^N X_0 = c_1 \\lambda_1^N x_1+\n", "c_2 \\lambda_2^N x_2+\\cdots+\n", "c_n \\lambda_n^N x_n\n", "$$\n", "となる．これを変形すると，\n", "$$\n", "A^NX_0 = X_{N}\n", "= c_1 \\lambda_1^N \\left\\{ x_1+\n", "\\frac{c_2}{c_1}\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^N x_2+\\cdots+\n", "\\frac{c_n}{c_1}\\left(\\frac{\\lambda_n}{\\lambda_1}\\right)^N x_n \\right\\}\n", "$$\n", "となる．$|\\lambda_1|>|\\lambda_i|(i\\ge2)$だから括弧の中は$x_1$だけが生き残る．\n", "\n", "こうして最大固有値に属する固有ベクトルが，反復計算を繰り返すだけで求められる．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Jacobi回転による固有値の求め方\n", "\n", "固有値を求める手法として，永年方程式を解くというやり方は回りくどすぎる．少し古めかしいが非対角要素を0にする回転行列を反復的に作用させるJacobi(ヤコビ)法を紹介する．現在認められている最適の方策は，ハウスホルダー(Householder)変換で行列を単純な三重対角化行列に変形してから，反復法で解を追い込んでいくやり方である．Jacobi法は，Householder法ほど万能ではないが，10次程度までの行列には今でも役に立つ．\n", "\n", "$$\n", "B := \\left[{\n", "\\begin{array}{ccccccc}\n", "\\cdot & \\cdot & 0 & 0 & 0 & 0\\\\\n", "\\cdot & \\cdot & \\cdot & 0 & 0 & 0\\\\\n", "0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n", "0 & 0 & \\cdot & \\cdot & \\cdot & 0\n", "\\end{array}}\n", "\\right] \n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapleでみる回転行列 \n", "\n", "行列の軸回転の復習をする．対称行列$B$に回転行列$U$を作用すると\n", "$$\n", "B.U = \n", "\\left(\n", "\\begin{array}{cc}\n", "{a_{1\\,1}} & {a_{1\\,2}}\\\\\n", "{a_{2\\,1}(={a_{1\\,2}})} & {a_{2\\,2}}\n", "\\end{array}\n", "\\right)\n", "\\left( \n", "\\begin{array}{cc}\n", "\\cos(\\theta) & -\\sin(\\theta)\\\\\n", "\\sin(\\theta) & \\cos(\\theta)\n", "\\end{array}\n", "\\right) \n", "$$\n", "\n", "となる．回転行列を4x4の行列に\n", "$$\n", "U^t B U\n", "$$\n", "と作用させたときの各要素の様子を以下に示した．\n", "```maple\n", "> restart:\n", "> n:=4:\n", "> with(LinearAlgebra):\n", "> B:=Matrix(n,n,shape=symmetric,symbol=a);\n", "```\n", "\n", "$$\n", "B := \\left[{\n", "\\begin{array}{cccc}\n", "{a_{1, \\,1}} & {a_{1, \\,2}} & {a_{1, \\,3}} & {a_{1, \\,4}} \\\\\n", "{a_{1, \\,2}} & {a_{2, \\,2}} & {a_{2, \\,3}} & {a_{2, \\,4}} \\\\\n", "{a_{1, \\,3}} & {a_{2, \\,3}} & {a_{3, \\,3}} & {a_{3, \\,4}} \\\\\n", "{a_{1, \\,4}} & {a_{2, \\,4}} & {a_{3, \\,4}} & {a_{4, \\,4}}\n", "\\end{array}}\n", "\\right] \n", "$$\n", "\n", "c = cos($\\theta$)\n", "\n", "s = sin($\\theta$)\n", "\n", "```maple\n", "> U:=Matrix(n,n,[[c,-s,0,0],[s,c,0,0],[0,0,1,0],[0,0,0,1]]);\n", "\n", "#U:=Matrix(n,n,[[c,-s],[s,c]]);\n", "```\n", "$$\n", "U := \\left[ \n", "{\\begin{array}{ccrr}\n", "c & - s & 0 & 0 \\\\\n", "s & c & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}}\n", "\\right] \n", "$$\n", "\n", "```maple\n", ">TT:=Transpose(U).B.U;\n", "```\n", "\n", "$$\n", "\\mathit{TT} := \\\\\n", "{\\begin{array}{c}\n", "\\left[ \\right. \n", "(c\\,{a_{1, \\,1}} + s\\,{a_{1, \\,2}})\\,c + (c\\,{a_{1, \\,2}} + s\n", "\\,{a_{2, \\,2}})\\,s\\,, \\, - (c\\,{a_{1, \\,1}} + s\\,{a_{1, \\,2}})\\,s\n", "+ (c\\,{a_{1, \\,2}} + s\\,{a_{2, \\,2}})\\,c\\,, \\\\\n", "c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\\,, \\,c\\,{a_{1, \\,4}} + s\\,{a_{2\n", ", \\,4}} \\left. \\right] \\\\\n", "\\left[ \\right. \n", "( - s\\,{a_{1, \\,1}} + c\\,{a_{1, \\,2}})\\,c + ( - s\\,{a_{1, \\,2\n", "}} + c\\,{a_{2, \\,2}})\\,s\\,, \\, - ( - s\\,{a_{1, \\,1}} + c\\,{a_{1, \n", "\\,2}})\\,s + ( - s\\,{a_{1, \\,2}} + c\\,{a_{2, \\,2}})\\,c\\,, \\\\\n", "- s\\,{a_{1, \\,3}} + c\\,{a_{2, \\,3}}\\,, \\, - s\\,{a_{1, \\,4}} + c\n", "\\,{a_{2, \\,4}} \\left. \\right] \\\\\n", "\\left[ c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\\,, \\, - s\\,{a_{1, \n", "\\,3}} + c\\,{a_{2, \\,3}}\\,, \\,{a_{3, \\,3}}\\,, \\,{a_{3, \\,4}} \n", "\\right] \\\\\n", "\\left[ c\\,{a_{1, \\,4}} + s\\,{a_{2, \\,4}}\\,, \\, - s\\,{a_{1, \n", "\\,4}} + c\\,{a_{2, \\,4}}\\,, \\,{a_{3, \\,4}}\\,, \\,{a_{4, \\,4}} \n", "\\right] \n", "\\end{array}}\n", "$$\n", "\n", "```maple\n", ">expand(TT[1,1]);\n", "expand(TT[2,2]);\n", "expand(TT[1,2]);\n", "expand(TT[2,1]);\n", "```\n", "$$\n", "c^{2}\\,{a_{1, \\,1}} + 2\\,c\\,s\\,{a_{1, \\,2}} + s^{2}\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "s^{2}\\,{a_{1, \\,1}} - 2\\,c\\,s\\,{a_{1, \\,2}} + c^{2}\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "- s\\,c\\,{a_{1, \\,1}} - s^{2}\\,{a_{1, \\,2}} + c^{2}\\,{a_{1, \\,2}}\n", "+ c\\,s\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "- s\\,c\\,{a_{1, \\,1}} - s^{2}\\,{a_{1, \\,2}} + c^{2}\\,{a_{1, \\,2}}\n", "+ c\\,s\\,{a_{2, \\,2}}\n", "$$\n", "この非対角要素を0にする$\\theta$は以下のように求まる．\n", "\n", "|　　　　　　　　 |\n", "|:----|\n", "\n", "このとき注目している$i,j=1,2$以外の要素も変化する．\n", "```maple\n", ">expand(TT[3,1]);\n", "expand(TT[3,2]);\n", "```\n", "$$\n", "c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\n", "$$\n", "$$\n", "- s\\,{a_{1, \\,3}} + c\\,{a_{2, \\,3}}\n", "$$\n", "これによって一旦0になった要素も値を持つが，何度も繰り返すことによって，徐々に0へ近づいていく．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## C 言語での実装\n", "\n", "lapack_codes/c_versions/Jacobi2.c にコードを入れている．\n", "```bash\n", "$ gcc Jacobi2.c\n", "$ a.out < input.txt\n", "```\n", "\n", "で固有値を求めていく過程が表示される．\n", "\n", "``` c\n", "#include \n", "#include \n", "\n", "#define M 10\n", "void PrintMatrix(double a[M][M], int n);\n", "\n", "int main(void){\n", " double a[M][M],v[M][M];\n", " double eps=0.0001,div,r,t,s,c,apj,aqj,aip,aiq,vip,viq;\n", " int i,j,n,iter,count,iterMax=1000000,p,q;\n", " \n", " scanf(\"%d\",&n);\n", " for(i=1;i<=n;i++){\n", " for(j=1;j<=n;j++) scanf(\"%lf\",&a[i][j]);\n", " }\n", " PrintMatrix(a,n);\n", " for(i=1;i<=n;i++){\n", " for(j=1;j<=n;j++) v[i][j]=0.; \n", " v[i][i]=1.;\n", " }\n", " \n", " for(iter=1;iter<=iterMax;iter++){\n", " count=0;\n", " for(p=1;p<=n-1;p++){\n", " for(q=p+1;q<=n;q++){\n", "\tif(fabs(a[p][q]) 1000 [dim] 2.5200 [sec] #BOB\n", "> 1000 [dim] 0.4700 [sec] #LAPACK\n", "```\n", "となった．2006年に初めてこの計算に用いたPCはMacBook(2GHz, Intel Core Duo)であるが，\n", "この計算での0.47秒は1.4GFLOPに相当する．\n", "07年のMacBook(2GHz, Intel Core 2 Duo)ではさらに速くなって\n", "```maple\n", "bob% gcc -O3 bob.c -o bob\n", "bob% ./bob\n", "1000\n", " 1000 [dim] 1.7543 [sec] #BOB\n", "bob% gcc -O3 lapack.c -llapack -lblas -o lapack\n", "bob% ./lapack\n", "1000\n", " 1000 [dim] 0.1893 [sec] #LAPACK\n", "```\n", "で，3.5GFLOPSが出ている．今(2016年)は，MacBookAir(2.2GHz, Intel Core i7)で...\n", "top500.orgが毎年２回High Performance Computerのランクを発表している．\n", "今は，Top1は100PFlopsであるが，\n", "初回の1994年6月の500位は0.4GFlopsで，今のlaptopがはるかに凌いでいる．\n", "まさにlaptopスパコンの時代なんですよ．\n", "\n", "ライブラリーは世界中の計算機屋さんがよってたかって検証しているので，バグがほとんど無く，また，高速である．\n", "初学者はライブラリーを使うべきである．ただし，下のサンプルプログラムの行列生成の違いのように，ブラックボックス化すると思わぬ間違い（ここではFortranとCでの行列の並び順の違いが原因)をしでかすことがあるので，プログラムに組み込む前に必ず小さい次元(サンプルコード)で検証しておくこと\\footnote{少し前(2002年ごろ)GotoBLASが開発されて，性能が10%ほども上がった}．\n", "\n", "添付のコードはちょっと長いが時間があればフォローせよ．コンパイルは，OSXでは\n", "```maple\n", "> gcc -O3 -UPRINT lapack.c -llapack -lblas\n", "```\n", "とすればできる．linuxではLAPACK, BLASがインストールされていれば，\n", "```maple\n", "> #include \n", "```\n", "をコメントアウトして，\n", "```maple\n", "> gcc -O3 -DPRINT lapack.c -L/usr/local/lib64 -llapack -lblas -lg2c\n", "```\n", "などとすればコンパイルできるはず．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## python謹製 lapack利用逆行列" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 [dim], 0.02294 [sec] # python\n", "2000 [dim], 0.08076 [sec] # python\n", "4000 [dim], 0.52680 [sec] # python\n" ] } ], "source": [ "import numpy as np\n", "import scipy.linalg as la\n", "import time\n", "\n", "sizes = [1000, 2000, 4000]\n", "for n in sizes:\n", " A = np.random.random((n, n))\n", " b = np.random.random((n))\n", " start = time.time()\n", " la.lapack.dgesv(A, b)\n", " print('%s [dim], %10.5f [sec] # python' % (n,time.time()-start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "|size|elapsed_time[sec]|\n", "|:----|----:|\n", "|100 |0.004878|\n", "|1000 |0.085786|\n", "|10000|46.031718|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 計算速度競争\n", "\n", "\n", "Mapleと，C, pythonで競争させてみました．Cはそのままです．Mapleは\n", "``` octave\n", "with(LinearAlgebra):\n", "data:=[1000,2000,4000];\n", "for n in data do\n", " A:=RandomMatrix(n,n,generator = 0..1.0):\n", " b:=RandomVector(n,generator = 0..1.0):\n", " st:=time();\n", " LUDecomposition():\n", " print(time()-st);\n", "end:\n", "```\n", "\n", "です．pythonのcodeは，lapackにあるdgesvを指定して\n", "呼び出すようにしています．scipy.laのsolveでは何を使っているかはよくわからないので．\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "結果は，次の通りです．\n", "\n", "|size|Maple[sec]|C[sec]|python[sec]|\n", "|:----|-----:|-----:|------:|\n", "|1000 | 0.339|0.0589|0.04101\n", "|2000 | 1.384|0.3826|0.12682\n", "|4000 |45.482|3.0111|0.81373" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pythonの圧勝ですね．絶対値はあまり意味がありません．MacBook Air(13-inch, Early 2015/2.2GHz Core i7, ElCapitan 10.11.6)ですが，softのversionにもよります．\n", "\n", "ただし，どんなrandomを生成しているかで結果は大きく変わるので，ちょっと怪しいです．\n", "たとえば，MapleのRandomMatrixでgeneratorを指定しないと\n", "size=1000でも9.773[sec]ととんでもないぐらい時間がかかります．\n", "\n", "MapleではさらにRandomMatrixの生成にも時間がかかります．\n", "MapleではNAGのライブラリを使っているので．．．" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 課題\n", "\n", "1. 次の行列$A$　の固有値をpythonで求めよ．\n", "また，対角化行列$P$ を求めて，対角化せよ．\n", "$$ A\\, = \\, \\left[ \\begin {array}{ccc} 4&5&5\\\\\n", "-4&-5&-7 \\\\\n", "4&4&6\n", "\\end {array} \\right]\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.000 9.000 2.000]\n", "[[-0.905 0.229 -0.707]\n", " [-0.302 0.688 0.000]\n", " [ 0.302 -0.688 0.707]]\n" ] } ], "source": [ "import numpy as np\n", "np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n", "\n", "A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])\n", "l, P = np.linalg.eig( A )\n", "print(l)\n", "print(P)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. -0. 0.]\n", " [-0. 9. 0.]\n", " [ 0. 0. 2.]]\n" ] } ], "source": [ "from pprint import pprint\n", "dA = np.dot(np.linalg.inv(P),np.dot(A,P))\n", "np.set_printoptions(precision=3, suppress=True)\n", "print(format(dA))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 行列のpprint設定イロイロ" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "'[[ 1. -0. 0.]\\n [-0. 9. 0.]\\n [ 0. 0. 2.]]'\n" ] } ], "source": [ "from pprint import pprint\n", "dA = np.dot(np.linalg.inv(P),np.dot(A,P))\n", "pprint(format(dA))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1., -0., 0.],\n", " [-0., 9., 0.],\n", " [ 0., 0., 2.]])\n" ] } ], "source": [ "np.set_printoptions(precision=3)\n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1., -0., 0.],\n", " [-0., 9., 0.],\n", " [ 0., 0., 2.]])\n" ] } ], "source": [ "np.set_printoptions(precision=3, suppress=True)\n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1.000, -0.000, 0.000],\n", " [-0.000, 9.000, 0.000],\n", " [ 0.000, 0.000, 2.000]])\n" ] } ], "source": [ "np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": { "height": "12px", "width": "252px" }, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": "block", "toc_window_display": false }, "vscode": { "interpreter": { "hash": "f3f87633aac09da3bda522f97956bee375b5501d1579e6458804e567301cb62a" } } }, "nbformat": 4, "nbformat_minor": 2 }