{
"cells": [
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def log_1(t):\n",
" if t <= 0:\n",
" return -float(\"inf\")\n",
" return np.log(t)\n",
"log = np.vectorize(log_1)\n",
"\n",
"# comme en TP, on fournit des fonctions permettant de vérifier le gradient et/ou la hessienne\n",
"def verifier_gradient(f,gradf,x0):\n",
" N = len(x0)\n",
" gg = np.zeros(N)\n",
" for i in range(N):\n",
" eps = 1e-5\n",
" e = np.zeros(N)\n",
" e[i] = eps\n",
" gg[i] = (f(x0+e) - f(x0-e))/(2*eps)\n",
" print('erreur numérique dans le calcul du gradient: %g (doit être petit)' % np.linalg.norm(gradf(x0)-gg))\n",
" \n",
"def verifier_hessienne(gradf,hessf,x0):\n",
" N = len(x0)\n",
" H = np.zeros((N,N))\n",
" for i in range(N):\n",
" eps = 1e-5\n",
" e = np.zeros(N)\n",
" e[i] = eps\n",
" H[i,:] = (gradf(x0+e) - gradf(x0-e))/(2*eps)\n",
" print('erreur numerique dans le calcul de la hessienne: %g (doit etre petit)' % np.sum((H-hessf(x0)))**2)\n",
" \n",
"# on donne également la fonction pas_armijo/gradient_armijo\n",
"def pas_armijo(f,gradf,x,v):\n",
" t = 1\n",
" m = np.dot(v,gradf(x))\n",
" alpha=0.3\n",
" beta=0.5\n",
" while f(x+t*v) > f(x) + alpha*t*m:\n",
" t = beta*t\n",
" return t\n",
"\n",
"def gradient_armijo(f,gradf,x0,err=1e-4):\n",
" # ne pas hésiter à copier-coller et adapter gradient_optimal...\n",
" x = x0\n",
" G = []\n",
" k = 0 # nombre d'itérations\n",
" maxiter=200\n",
" while(True): \n",
" k = k+1\n",
" if k > maxiter:\n",
" print('erreur: nombre maximum d\\'itérations atteint')\n",
" break\n",
" # calcul de d, t, \n",
" d = -gradf(x)\n",
" G.append(np.linalg.norm(d))\n",
" if np.linalg.norm(d) <= err:\n",
" break\n",
" t = pas_armijo(f,gradf,x,d)\n",
" x = x + t*d\n",
" return x,np.array(G)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Partie 1.2 -- Étude numérique \n",
"\n",
"$$\\newcommand{\\eps}{\\varepsilon}$$\n",
"\n",
"**QN1.** Écrire deux fonctions `f(x)`, `gradf(x)` calculant respectivement $f(x),\\nabla f(x)$. *(On pourra utiliser la fonction `verifier_gradient()` pour vérifier le calcul)*\n",
"\n",
"**(Attention, il faut utiliser la fonction log donnée ci-dessus plutôt que np.log)**"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"erreur numérique dans le calcul du gradient: 3.98877e-10 (doit être petit)\n"
]
}
],
"source": [
"#<\n",
"eps = 1e-3\n",
"def f(x):\n",
" return .5*np.linalg.norm(x - z)**2 - eps*(log(x[0]) + log(x[1]) + log(1-x[0]-x[1]))\n",
"def gradf(x):\n",
" return x - z - eps*(np.array([1/x[0], 1/x[1]]) + 1/(x[0]+x[1]-1)*np.ones(2))\n",
"#>\n",
" \n",
"verifier_gradient(f,gradf,np.random.rand(2)/3)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**QN2**. La fonction `x,G = gradient_armijo(f,gradf,x0)` donnée en préambule implémente la descente de gradient avec rebroussement d'Armijo. Elle retourne la solution calculée `x` et un tableau contenant $(\\|{d^{(k)}}\\|)$. En utilisant cette fonction, tracer sur un même figure et en échelle logarithmique $k\\mapsto \\|{d^{(k)}}\\|$ en partant de $x^{(0)} = (\\frac{1}{4},\\frac{1}{4})$, pour plusieurs choix du paramètre $\\eps \\in \\{.1,.01,.001,.0001\\}$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd4VFX6xz9nenqHdNIgJEBoAZSisFhAAVFXBLtiX9eyuqurq65i7wXXimJHl99KUcGCSBGlCdJCD5DeezL9/P64SUggIYGQzCRzP8+TJ5kzt7wzk7nf+5bzHiGlREVFRUXF89C42gAVFRUVFdegCoCKioqKh6IKgIqKioqHogqAioqKioeiCoCKioqKh6IKgIqKioqHogqAioqKioeiCoCKioqKh6IKgIqKioqHonO1ASciNDRUxsXFudoMFRUVlW7F5s2bi6WUYW1t55YCIISYCkxNSkpi06ZNrjZHRUVFpVshhDjcnu3cMgQkpVwqpbw5ICDA1aaoqKio9FjcUgCEEFOFEO9UVFS42hQVFRWVHotbCoDqAaioqKh0Pm6fA1BRUXFfbDYb2dnZmM1mV5vikZhMJqKjo9Hr9ae0v3Dn9QDS09OlmgRWUXFfMjMz8fPzIyQkBCGEq83xKKSUlJSUUFVVRXx8fLPnhBCbpZTpbR3DLUNAKioq3QOz2axe/F2EEIKQkJAOeV9uKQBqElhFpfugXvxdR0ffe7cUADUJrOIpVK9ahS0nx9VmqHgobikAKiqeQs49f6Nk3jxXm6HSDpYvX05ycjJJSUk888wzLW6zevVqhg0bhk6nY+HChV1s4cnjlgKghoBUPAHpdOKsrcV6JMvVpqi0gcPh4C9/+QvLli1j165dfP755+zateu47WJjY5k/fz5XXHGFC6w8edxSANQQkIonIC0WAGzZ2S62pPvzySefMHLkSIYMGcItt9yCw+HA19eXe++9l2HDhjFx4kSKiooAeO2110hNTSUtLY2ZM2e26/gbNmwgKSmJhIQEDAYDM2fOZPHixcdtFxcXR1paGhqNW15aj8Mt5wGoqHgCzvrqDVtODtLpRHSTi0ZrPLZ0J7tyK0/rMVMj/Xl06oATbpORkcEXX3zBL7/8gl6v5/bbb+fTTz+lpqaGYcOG8eKLL/L444/z2GOPMXfuXJ555hkyMzMxGo2Ul5cDsHLlSu65557jju3t7c26devIyckhJiamcTw6Opr169ef1tfqClQBUFFxEbJeAKTNhr2oCH3v3i62qHuyYsUKNm/ezIgRIwCoq6ujV69eaDQaLr/8cgCuuuoqLrnkEgDS0tK48sormT59OtOnTwdgwoQJbN26tdVztDRfqidUP7mlAKgzgVU8AWfd0fptW1ZWtxeAtu7UOwspJddeey1PP/10s/E5c+Y0e9xwwf7mm29YvXo1S5YsYc6cOezcuZM1a9ac0AOIjo4mK+toriY7O5vIyMhOeDVdi1v6nGoOQMUTkOa6xr+tah7glJk4cSILFy6ksLAQgNLSUg4fPozT6WysxPnss88YO3YsTqeTrKwsJkyYwHPPPUd5eTnV1dWNHsCxP+vWrQNgxIgR7Nu3j8zMTKxWKwsWLGDatGkue82nC7f0AFRUPAFnkxmctmx1LsCpkpqayhNPPMF5552H0+lEr9fzxhtv4OPjw86dOxk+fDgBAQF88cUXOBwOrrrqKioqKpBScs899xAYGNjmOXQ6HXPnzuX888/H4XBwww03MGCA4vE88sgjpKenM23aNDZu3MjFF19MWVkZS5cu5dFHH2Xnzp2d/RacMmovIBUVF1Gzbh1HbpgNQMD06UQ+83Qbe7gfGRkZpKSkuNqMFvH19aW6utrVZnQ6LX0Gai8gFRU3p8EDEN7eWLPVuQAqXU+XCYAQIkUI8ZYQYqEQ4rauOq+KirvirFNyAMakJDUE1Al4wt1/R+mQAAgh3hdCFAohdhwzPkkIsUcIsV8I8QCAlDJDSnkrMANo0zVRUenpNJSBGpOSsBcU4LRaXWyRiqfRUQ9gPjCp6YAQQgu8AUwGUoFZQojU+uemAWuBFR08r4pKt6ehDNSYlARSYs/NdbFFKp5GhwRASrkaKD1meCSwX0p5UEppBRYAF9Vvv0RKORq4srVjCiFuFkJsEkJsapi6raLSE5GWegHoq8x3sWappaAqXUtnlIFGAU0zWtnAKCHEeOASwAh829rOUsp3hBB5wFSDwTC8E+xTUXELGj2AxEQAbDmqAKh0LZ2RBG5pfrSUUv4spbxTSnmLlPKNEx1AnQim4glIcx3CaEQXHo7Q69WmcG5Oe9pBWywWLr/8cpKSkhg1ahSHDh0CoKSkhAkTJuDr68sdd9zRhVafmM4QgGwgpsnjaOCkgptqO2gVT8BZZ0aYTAiNBn1kJFa1EshtaW876Hnz5hEUFMT+/fu55557uP/++wFl8fY5c+bwwgsvdLXpJ6QzBGAj0FcIES+EMAAzgSWdcB4VlW6N02JGYzIBoI+JwZalzgU4VdylHfTixYu59tprAfjzn//MihUrkFLi4+PD2LFjMdV/3u5Ch3IAQojPgfFAqBAiG3hUSjlPCHEH8B2gBd6XUp7UXGgp5VJgaXp6+k0dsU9FxZ2RdU0EIDoK8/btLraogyx7APJP82sIHwSTWw63NOBO7aCbbqfT6QgICKCkpITQ0NCOvAudRocEQEo5q5XxbzlBorct1G6gKp6A02xGeHkBYIiOxlFRgaO6Gq2vr4st6164Uzvo7tY22i2bwakegIonIOvq0BiNAOijlbtGW3Y22v79XWnWqdPGnXpn4U7toBu2i46Oxm63U1FRQXBw8Ol4mZ2CW/YCUpPAKp6A02Jp9AD00dGAujzkqeBO7aCnTZvGhx9+CMDChQv505/+pHoAJ4vqAah4ArKuDq2fHwCGWMUDsB4+7EqTuiXu1A569uzZXH311SQlJREcHMyCBQsajxEXF0dlZSVWq5VFixbx/fffk5qa2mnvS3tQ20GrqLiIAxdciDG5H9EvvwzA3jFj8Z0wnsgnnnCxZe1HbQftenpcO2g1BKTiCTjNdWhMXo2PDfFxWDMPucweFc/DLQVAnQms4glIswVhMjY+NsbHY83MdKFFPQtPuPvvKG4pACoqnoDTbG7uAcTF4ygtxaF6vipdhFsKgBoCUunpSCmVMlCvozNDDfHxAKoXoNJluKUAqCEglZ6OtFpBSsQxOQAAi5oHUOki3FIAVFR6Og2rgWma5AAM0dGg06kegEqXoQqAiooLaFwQvklzMKHXY4iNxZp50FVmqZyAjrSDBnj66adJSkoiOTmZ7777rnH8hhtuoFevXgwcOLCzX8JxuKUAqDkAlZ6OrF8QXuPl1WzcEB+PRfUA3I6OtoPetWsXCxYsYOfOnSxfvpzbb78dh8MBwHXXXcfy5cu79PU04JYCoOYAVHo6LXkAAMb4OGyHjyDrLw4q7cPd20EvXryYmTNnYjQaiY+PJykpiQ0bNgBw1llnuaxfkFu2glBR6ekczQE0FwBDfDzSZsOWk4MhNtYVpp0yz254lt2lu0/rMfsH9+f+kfefcJvu0A46JyeHM844o9n+OTmuXwBIFQAVFRfgPIEAgFIK2t0EwFV0h3bQ7tomWhUAFRUX4KzPATQtA4WjAmDJzMT37LO73K6O0NademfRHdpBt3f/rqbLcgBCiOlCiHeFEIuFEOd11XlVVNyRxhCQV3MPQBcUhDYgQO0JdBJ0h3bQ06ZNY8GCBVgsFjIzM9m3bx8jR47soneodTq6JOT7wBSgUEo5sMn4JOBVlCUh35NSPiOlXAQsEkIEAS8A33fk3O6Eo7qG3Pvuo/e/HlJquVVU2sBptgDHJ4EBDAkJWA+qpaDtpTu0gx4wYAAzZswgNTUVnU7HG2+8gVarBWDWrFn8/PPPFBcXEx0dzWOPPcbs2bM77w1ripTylH+As4BhwI4mY1rgAJAAGIA/gNQmz78IDGvP8YcPHy67AzWbf5e7kvvL0s8+c7Up3YYtBVvklP9NkZWWSleb4hJKP/tM7kruL21FRcc9l/PPB+WesWNdYNXJs2vXLleb0Co+Pj6uNqFLaOkzADbJdlxjOxQCklKuBkqPGR4J7JdSHpRSWoEFwEVC4VlgmZTy99aOKYS4WQixSQixqaFsy91xlJcBYD18xMWWdB+2F2/nUOUhMis8s+bdWddQBup13HOG+DgcRcU41G6WKp1MZ+QAooCsJo+z68f+CpwD/FkIcWtrO0sp35FSpksp08PCwjrBvNOPo6xeALKy2thSpYEys/Ke5dXkudgS1+A0108Ea9IKogGj2hTutKC2g26bzqgCaqm2SUopXwNea9cBhJgKTE1KSjqthnUWDQJgO6Iu59deyi1K/XV+Tb6LLXENss4Mej1Cd/xXUF9fR27LycFr0KCuNk3Fg+gMDyAbiGnyOBrI7YTzuA320noP4EgW0ul0sTXdA08XAKfFfNwcgAb09eWBtpwe/bVRcQM6QwA2An2FEPFCCAMwE1jSCec5/dQUw4Z34STXSW7wAKTFgr2+FE3lxHh6CEjWtS4AWj8/NH5+2HJVAVDpXDokAEKIz4FfgWQhRLYQYraU0g7cAXwHZABfSil3nsxxpat6AW37Er69D0r2n9RuDQIAYD2iJoLbg8d7AGYzwuv4BHAD+shIVQBUOp2OVgHNklJGSCn1UspoKeW8+vFvpZT9pJSJUsonT/a4LusGWlnfm6P05GqwHWVl6Pso0/ZtqgC0C4/3AMx1aIzHJ4AbUAXA/eisdtCtHXfu3LkkJSUhhKC4uLhTXpPaDbQppygA9rIyTCmpoNerpaDtQEpJuaUcndBRai7F4rC42qQux2m2qB5AN6Kz2kGf6Lhjxozhxx9/pE+fPp32utxSAFznAdR/4UoOnNRujrIydKGhGKKi1BBQO6iyVeGQDhIDEwEoqClwsUVdj6yrazUHAIoAOKuqcFRVdaFV3Zfu2g76RMcdOnQocXFxp+cNagW3bAYnpVwKLE1PT7+pS0/cIADHeAA/HfmJGlsNUxOnHreLtNlwVlWhDQ5C3ycWa5YqAG1Rblbi/ykhKewp20NeTR6x/p7V+dJpNqMNab0HvD6qvhIoNxdtcnJXmdUh8p96CkvG6W0HbUzpT/iDD55wm+7eDro9x+0s3FIAXILTAVX18ehjBGDe9nmUW8pbFABH/T+QLigIQ0wsdZs2I6V0i1avHcVRXs6ByRcQ9cor+Iw6fY2rSs3K5PH+wf0Bz0wEO8116FuYBdxA01JQUzcRAFfRndtBO1soG+/Ka4dbCoBLJoLVFIHTDl7BUH4EHDbQ6gHIqspSwhZOB1qNttluDXMAtEFBGGJjcdbU4CgtRRcS0nW2dxKWgwdxlJVRt3XraRWAhgqglOAUwDMTwdJsaXEWcANHBcD1i4a0l7bu1DsL2c3bQbuyTbRb5gBckgRuSADHjQXpUEQAqLZWU2Ypw+60U1h7fI1/QwmoNigYQ30lUE9JBNtylQvz6b4INVQA9fbpTYgp5LR4AA6ng0prZYeP01U4zeYW+wA1oA0JQRiNaiK4HXTndtDtPW5n4ZYC4BIa4v9x45Tf9WGgrKom6lydfdxuDY3gtEGB6GPqS0FPMQ+wLmcdZrv5lPbtDGx5yntyugWgwQMIMgYR7hN+WgTgf/v/x6SFk6i11Xb4WF1BW0lgIQT6iAhVANpB03bQaWlpnHvuueTl5TVrB/3TTz/xyCOPNLaDHjRoEEOHDj2ldtApKSnMmDGjWTvoJUuUua6zZ8+mpKSEpKQkXnrppcayzqbtoCdNmtTYDvpEx33ttdeIjo4mOzubtLQ0brzxxtP+3qkhoAYaBWCs8rslAajKZkT4iGa7NXgAuqAgNAEBoNGckgeQVZnFLT/ewsNnPMyM5Bmn8AKaU1BTwBXfXsF/Jv6H5OBTiyHb8zrJA7CUYdAY8NJ5EeETwcGKjve+312ymypbFZmVmQwIGXAarOxclIlgrQsAqKWgJ8Pll1/eGO9vypw5c44LBa1du/aUznHBBRdwwQUXHDf++OOPN/5tMpn473//2+L+Dz30EA899FC7j3vnnXdy5513npKt7cUtPQCXhYC0BgjrDwa/xlLQBgEQiGZi0IC9VEloagMD0RgM6CMiTqkU9ECFcr7T1R55R8kOCmsL+b2w1c7bbdIYAsrNPa09jsrN5QSaAhFCNHoALSXJToacGkWkDpSfXAmvK5A2GzgcJ/QAQKkEUgVApTNxSwFwCZW54B8JGg0ExzfzAIKMQUT6RrYcAiorR+Pvj9ArCWNDn9ijArD3O/houlJh1AYNF/6WROZUyK5SbD1SeQIx2vgefP+vVp+21XsA0mrFfhpnIpZZyggyBgEQ7hNOrb22w/H7nCpFAPaXn1wbD1fQsCB8S6uBNUUfGYmjpKRxe5WTQ20H3TaqADRQmQv+UcrfwQmNApBdlU2MXwzRvtHkVB8fCnGUlaENOhpD1MfEHm0HsfsbOLgSStu+q28QgCNVpyeB3CAkJxSUrZ/DpvmtNr+z5eWhr1/i8nSGgRo8AIAInwigY6WgUsrGSqKD5e6/lGLDgvCaEySBoUklUK57V0l11HtTOXU6+t6rAtBAZY7iAYAiAOWHwWEnqyqLaL9oov2iG++qm+IoK0MXGNT42BAbi6O8HEdFxdGmcgXb2zz9ocpDgCI4jnZ4DG3RcOFvVVCkhKI9YK06WgHVBEd1Nc7KSrzT04HT25q43FLezAOAjglAcV0xFocFgegWHkBrC8Ify1EBcN8wkMlkoqSkRBUBFyClpKSkBFMbnuSJUJPAoFwMG0JAACGJ4LRjKztAfm0+MX4xmHQmSs2l1Npq8dZ7N+5qLytD37t342NDfByg1NB7F+9VBvO3w4CLT2hCZkUmJq0Js8NMQW0Bkb4dqwVuEIAGQTl2/gIV2crFH6BoNwQ0X8y+IQHsnT6cikWLTqsHUGYpI9B41AOILZTYvlgE9519Ssdr8MwGhQ5iW/G24z4jd6MxBGRsrwC471yAhiqV7rJ8a0/DZDIRHR3d9oat4JYC0OWtIGpLwGFtHgICcvI245ROYvxiMGqVSTvZ1dn0C+rXuKujrAxTSkrj44a/zX9sxrum/kuRv+OEpy8zl1FuKefs6LNZlb2KI1VHOiQAdqedvOo8gk3BlJpLWxaUoj1H/y7cDUnnNHu6If5vSEhEGxzcTACK64pZl7uOaYknX69sd9qptFQSZFI8gBCvEKZshKhty3Hc9G+0p5D4bxCAcdHj2Fa8ze0rgdrrAeh69wat1q09AL1eT3z9EpYq3Q81BARHQyBNQ0BAVpGyjEGMXwzRforKNg0DSSmPywHowsPRBgVh3rpRGfAOUTyAE9AQ/jkr+izlvB1MBOfV5GGXds6MPBNoJQxUlKH81vsoHsAxNMSd9VGR6KOimgnAgt0LeGjtQy2GxNqi0lqJRDYKgEZoSCpS/g3NO09q2YhGcquVC+S4aGUOhztWAtny8zHvUUT3RAvCN0XodOh693JrAVDp3qgCAEfnADQIgG9v0PuQVaZcSGL8YojyVbyDpolgWVuLtFjQBR3NAQghMKWkYN5df4edMg2qcqG2tNXTH6o4BMCoiFEYNAayKjsmAA0X5jGRY4BWKoEKd4NPL4ga1rIA5OWBTocuNPQ4AcgoVcRjd+nJN/5qmAXckAOQVisRhXYA6nacmgDkVOcQbAqmX1A/dBqdWwpAwZNPknXTzUgpkQ0LwrfhAQAYIqNUAVDpNLpMAIQQCUKIeUKIhV11znbT6AHUh4CEgOAEsmry8NJ5EeoVSqAxEB+9T7O7XnuZMqNV20QAAEwDUrFkFSIxQP8pyuAJvIDMikz0Gj3RvtHE+MV0uBKowYMY3ns4Rq2xZY+iKAN69VfmPRTtOa4SyJaXi753b4RW21iP3jAXYHeJcuFvEIKToUEAGqqALAcOoHUo5zbvOHGorDVyqnOI9o1Gr9ET5x/ndgIgpaR2y1bshYXY8/NxmpX1D9rKAYA6F0Clc+nokpDvCyEKhRA7jhmfJITYI4TYL4R4AEBKeVBKObsj5+s0KnNBowOfsKNjwfFkW8uJ8o1CCIEQgmjfaLILt0PJAaSUfPSd0j1QG9S8ra8pNRUcEgtxEDlEGWxDAPr490Gr0RLj33EByK7KRq/R09u7NzF+MRyuPNx8g4YKoLAUCEsGS+VRL6gee24e+gilRFMfFdU4F6C4rpjCOqXnyql4AE3bQACY69sHZ/YW1G7fdtLHA0UAGnIciYGJblcJZM/NxVE/j6Ju2/aT8gB0kZHYCwqRdnun2qjimXTUA5gPTGo6IITQAm8Ak4FUYJYQIrWD5+lcKnLALwKaVsr0SiVLWojxDm8civaNJLtgK3xzL5nFNSxbq1y8tAY7vDoYsjcD9QIAmGtDwCdUOXZB63e3hyoPER+gJNJi/GLIqszqUFldVlUWUb5RiqD4xRzvAVRkgbVa8QB6KUlrWdA8/GLLy0MXqQiAIUrxjGw5OY0X/QifCDJK2ucBNH0tZZZ6D6C+CsickYE0GVmXInDk5WMvKTmp1+pwOsiryWsmALnVuW7VE6hu21FhM2/f1iQH0LYAGOPiwOHAst+9RE2lZ9DRNYFXA8cGt0cC++vv+K3AAuCi9h5TCHGzEGKTEGJTl5WWNZ0DUI8zcQLZOi0x9qM1+dFSQ45W4Dz8C7/vyybAUgOArmgjlB2CPz4HQB8ZjkbvxFxuUHbsPbBVD8DmsJFVlUWcfxwAsX6xmB1miupaeO352+GPL44f/+Iq+O2txodZVVnE+MU0Hi+rKgunbNLKobD+zj0sBcL643TAoTufJO/hRwCQDge2ggL0Ecp7om8UgNxGAbgo6SKK6ooorjvxDOHiumImfDmBZZnLAKisUF5XYwgoIwNTcjIHoxTxPdlEcFFdEXanvTFHkxiQiESSWXl6WmqcDuq2/oEwGjGmplC3bTvORg/gxElgAO/6RURqflnXqTaqeCadkQOIAprecmYDUUKIECHEW8BQIcQ/W9tZSvmOlDJdSpkeFhbW2mYnxO44yb41TecA1FMUHItFoyGm4ugEpejKQiwaDcU4qMj4EX+rIgCagjXKBnuWgZSI8sOYAm2Y85QvOuGDlJCL3XrcqbOqs3BIR6MHEOundBRtMXG77AFYdGvzcE3+DshYCuvfAimRUjYXAP9YLA5L81bWDRVAvfqDTyjFu3tjPlRIxddf46ytVdo+2O1HQ0BNJiRllGQQ7RvNyHBlfYC2wkCf7PqEEnMJSw4sofz//o8zrnuD4VkGjFoj0unEvHs33gMGYKrvgFh3knmAhqR8gwAkBSpzR1yVB7Ds30/eww/jaNKGoG7bNkwDBuA9dBjmHTtw1ireyYkWhW9A37s3hsREatapAqBy+ukMAWhpORsppSyRUt4qpUyUUj7dwjZHD9DBNYE/+edlfHzlSBbcN5W18/7JvlULyNu7mZrKEpyOY2bZNk4Ci2o2nFVfWhiTnwH1yc+ofOVil+PlR2juKoJsNTiEQJO/HgJjoTJbCfUU78UUZMN8pBDpcED4QHDaWqy2aWgB0RgC8lcu3MeFbUoz4fBakM5GTwOAbQuU32WZUJihTFaz1zYKQMPvZscr3A2+4eAVRO3vv1OyXYspXI+sq6Nq5crGpKO+PgSk8fZGGxSELSeHjNIMUkJSGjuMnkgAKq2VfLHnC3RCR+3adeQ98igah5NJW5TnbTk5OKurMfbvT1qfUeQGC2q2nVwe4FgBiPGPcVklkNNsJvvuuyn/70Iq6td1lVYr5p078UpLwyttEM7aWsy7doFWC/X9o9rCZ/RoajdtwmmxdKb5Kh5IZ0wEywZimjyOBrq0jME/6zD9ttaglVXw9X7sLKIcKK9/3qYFuw6cGoFDpwUC4dulaJ/5CZ3egMHLFzQ27rc66GWwUub9Ej5njia69DB4R3KwVypnHtiE3ecCLAYDGhww6VlYcAU7tn3Ms6WbSUnw5uK9Vv712fU4I/wZkhtKwjU3YfALxcvbH6PJB6HRYKrJ4e+VDnx/e53sskzQaLmnzon/yrfJDd9ApaOa3VUHMNgqMdWG4OPtS8Bn8wmt7Y9PVCzGLQvRRI+A7E2w+xuyUs4FaOYBgOJRjAgfQWbBHtb9sR5jYBgBOxYR+fdX0AWZiJxQzJEVCVR+u4yAC5XWtA0eAChhoLqsw2QlZDE9aTr+Bn+ifaPZVbKrxc/Amp3Nt7/NQ1RU8VDMtSS88D7WPlHsDXcycH0+9rIyzLsUT8SUksqI4Ch2RLxD+PY/Go9RuWE9lj17MPYKR9crDFNyMhrv5jN8GwQgwrfeW9HoGW6NxOerlVTk9UMXFoYxMQFdaOgp/z/VbtxIzW/r0QYGog0MRBcWhj4yAl14OBqDoXG7wueex7r/ANrQUCq+WkTwlVdi3rMXabXiNTgNY3L/+uNtQmMytXvpP5/RZ1L28cfU/f47Pmcqcztqt2yhdv16hE4HWh0aH2/FtuBgjP37o/X1PeXX2xb2sjIqv/4GNAKh16Px8kLj54fW3x9DXBy64NbXOm7AWVND5Q8/oO/VC31sH/QR4Qitts39AOq278CY3K/Ze69yanSGAGwE+goh4oEcYCZwRSecp1Uu+XwTNnMtG3/5jm3rfqQ89zDe1ioCnRb8NHa0divSYsHhkDilA4cTrFqB1V6J1iHR24vQOSC0FhwlXuRnzEPo5xM80IgmCjYag5kqSknRlmEzOqnWh+DbbxJl0cO4O/c7nBodsb21gMT7YAE5BYcZuMbAkbBKKrVVGColPhovYn2ioaaUCKsW554/sNYUIyUkanWYcg5Rc6CCCks1EXY7BhsY7EY00oYEilbeRxEgBegjLfiH9yXQ9BVZ0cpFpkEAwr3D0Wl0HCnfj91sZvf1VzLsYB1QB2/9Eyfw8JVa9sYE8Nd+Ns5YtZLKSH80gO4YAajYpdydp8oIqn78kZTg/s08gPe2v0d+TT43acdTOvsvDLZamQfA+5T4a/j59jSK8g4ycF0elUu/xl5WClotxn59GapxsDRCw1k7y7EVFGLOPcKR665D2zSaZzDgM2oUvhPGE3jppWiMRnKqcujl1atxpra9tJRb383Bp9xC7sI1xXuiAAAgAElEQVT7Gnc1pqbgO2YsQbNmNoa0jqV84UJKP/mUmLffamzv4aiuJuvOu3DWr/vQDCHwOfNMgq6+CqSk7LPPCL72WvRRkRQ89TTmvXup26YImtfgwYpg+PnhrKhAexKC5D1iJOh01Kxbh8+ZZ+KsrSX7r3c2VhYdh06H15DB+I4dR+DlM5rNU2lK2ZdfYj1wkF4P3N+iGNlLSih84UVCbrgeY9++jeOFzzxDxeIlLZ9bo8F7+HD8zjuPgGlTW53ZXfDCC5R/vqDxsfD2xnvIYLzS0/E/99xm52tK9apVZN1yKz7jxhE99/V2hdFUWqdDAiCE+BwYD4QKIbKBR6WU84QQdwDfAVrgfSnlSWX2TkcrCL3Jm9ETL2b0xIupttj53+/ZvL82k0MltVx7Zh/+OakvxTt+4sDPHxNasYOHvR6if/9+pCcK/AOK2F26G6dXCOev+YDaomqK11kp2RLA45U6Xpz8B8PjY5nzyx8E+9pZZ0znHAEP+hsorXXwcYWT1N7x7DHlcr01nZpVvyL8tZxxkY7DQ6eRIay8nL+WWkcuPnofEo19eX/rTzDsGjj/SV744VZKinZyYXUJL/gbeTn1NiZ+8zClU14hLyiJvE+vodivLweKSiG/lsSCUtK2SEo3S4I3vkziCEmUnxIS0Wq0ROv8yNr4Dr+9sYy4gzUUj6wjftCFVNmCqQn35bL+TnI3vsm+tEhGry2jcuFXGLx0VOsdNHx99VFRsHIFgzKdhL79NNmlZUyels4PqUeoslaRU53Da7+/RlClk4kffobGz8Tb4zXcFXklEWYji2IPsbx6A15BXpT0CcT41Vfoe/XCmJCAxmjEG5ApifDjXmrWrePQS89Q4Qc/33s2ecWZmHOzGZolOOvATmrWrKF240aiXnqJ3Jrcxgog6XSS+4/78a5x8PBVWqYPu4qLA8+ibtt2ataupeSDD6j45muiPvqQ+cVLOa/PeSQFKTmDiq+/UZLgUlL47LNEvfQSACXvv4+zrIx3/5LAc1d+jKO8HHthIba8fKyZmVQsXkz2bbcDYExOJuzev+GsqaHg+Reo+N9X2EtL0IaFoouIQAiB16CB1Kz79aQuXFpfH7yGDFYSwffeS9lnn+EoLib2ow/xGjQIabfjrK5WbCsupnbjJmrWrqXolVco/eADwv72NwL/fGmzO2zrkSMUzHkCabOh6xVGyOzmFdrO2lqybr0N8/btWA8epM/nnyE0GiwHD1Kx9GuCrrma0FtuQdpsOGvrcFZX4aiopG7LFqp++J6CJ5+k9KOPiHnnbYzHtIow795N+RdfEnjZZfhfeCHWrCNYMnZTu3kzxa/PpfjNt+j9978TdPVVzYRJSknhq6+iDQigZs0asu+8k+jXX1c9gQ7QIQGQUs5qZfxb4NtTPe7pbgbna9RxzZlxzEiP4bnle3j/l0y+31VAQaUTH8O1zBwVQ0hJLYu25PPZegcPXZDCHWcpvXFKMjMJyXkW3XAdlUPOga/38OwHBv6YbKd3jYPdYRreDDrCmnWPsrYuh4dLyhhQVQ3x52JK9qPi//4HWi1xD0zD69B7BK1+lSHAuQZv/pU4iF/MeZxbVgwhSTDpaTD4EBM2gE1lu3nD18B44cvEIzsQRn9CBs8kRO/FwLTJsGsx+MKeKefyWoAP/9m9mlfWVmLaX8QTGQ6qgj/CMPsGRFUusRUFhOw0ErK+gN/O0HBdXBni8qshVqkwGVlTDN8/gzxvGvuivsI/J5fDvRw89M0sXp3wKn2D+qKPikRjtfPwAtAnBuM1bhwsXsLMYsHu8zN4e9s7hAl/XvvRH7vtMP+caSUsPp7BlyoLhQ8/9D2frfqZCksFuWenE/LRb1gPHsTv/PMaP6eoIWNwir3kPfZvtFYL3989iCcuexMhBJkVmTy1/iney/uNm7ZFcO43yylOTCInPIehvYcqn9W771Gzdi0R//43Ub3W8Vre/zhv/A30Gj2a0FtvwbxrF4evuZYdV1/GJzOq+WLPF3w46UNCNmeSe//9eKen4zVkCCXvvkvgZZdhTEqi+P33+a2/4Af/I2RpK4hPSMCYkNBoc9hf76Dqxx+p+uEHQu+4A43BgMZgwG/82VQsXYrG2xuvtMGNFzLTwEHUrPu1zdXAjvsfHjOGotdex5qVRcm77+Fz1jh8Ro5sfF7r59cYsvMdNw7+dg/mPXspmDOH/EcfpXzhQqLnzkXfuxcABc88C3o9PmeeQeFLL+OVlob3CGW1O2m3k33PPZh37iRg+nQqFi2iYskSAqdPp3juXITJROitt7YY6vEdN5awO/9K7ebNZN/xVw7PuoLo//wH72HKZySlpODJp9D6+9PrvnvRBgTgc8aoxv3tJSXkPfwIBU89Re2mTUQ8+QRaPz8Aqr7/AcuuDCKffQan2UL+o4+Sc9fdRL/6CkIVgVPCLVtBdNaKYCa9lkempvLJ7FEEeOm5bnQ8q/4xgYcuTOXda9L5/eFzuXBQBE9+m8FHvx7icEkNd21WvjBGYWdx3/MJ+egzLE4jZy12ElQBw7BS6zCxaP8iJsdN5jJdiHKy0H4YU5Ua+7A7/oLX1U/BvwrgH5lw88+E9j2f/2Ss59XCUm4qKYJL3wODD1BfCuq0glbPg4f3ILZ/CQMvAX192eCQK5U6fms1ycNu5IWzXyAspi93ne/Ps7O17EsLpuillzh81dXk3HAZVyyA6Sslm/oKzupbpGTpw5osE+kTCt4hiL3LCTzvTwBEJw2lzl7Hld9eyTcHv8GUXJ/0HdGb+C+/IOLppzFdMpVL1kny77yH8a+v48V5TjR7DxH/xL950FHEM+ajX8oxUWPQa5SkZ834oQi9HmmxYEo5OkVkWNxoskMBs4VPJmi48tJ/N1444wPieefcd3jurOf4Mt3C+qHeFM+dS9KGXFKPQOErr1D06qv4X3ABgZfP4L70+7A77bz6+6tHP//UVAr/fSPGogqe/cqX6T9Us3XWdLLvvBNjagqlj9/C5gsT0cfGkP/4HApffRVptbDkHOUCtOLIiuP+p4Rej//kyUS99FIzYQi4+BIcJSXYsrLwGjy4cdwrbRDQ9loAx+IzejRISc5dd+OoqCDszrva3MeU3I/Yjz8i8vnnsR44wJFrrsFWUED12l+o/uknQm+9laiXXsIQE0P23/5G7caNVCxdSvbdd1OzajXhjzxMxFNPYhqcRuGLL1L7++9UfruM4KuvbjPO7z18OHELPkcT4M+R66+n6D//wV5aStV331G7cSNhd93ZYnhIFxJC9Btz6fX3+6hasYJDl83AcjAT6XBQ9PprGBIS8J8yhaDLZxD+6CNUr1xJzj/uV4otVE4at+wG2tntoMf2DWX53WcdN27Sa3ll5hCsDiePLN5JsI8B6YzB5hOOs66SZzNCyA2Gr86+k0/3zEdmHiEwKJLqzLu5Z5qVm0dMZf/hWvqWfUy+IYagy8ej9Q8g5OabG14YeAcrP5fNR5M+mz+tfBIGXXZ0xjBHK4L+MuwuIpwL4dAa5aLfQJ/REBQPDhvEjsZLo+Hl8S8zc9F0MrzsbLl1DBNKzqDoxWexW0vRh4fyXVQFhr/cQGLUGCjeA17HxIVH3QY/P4V/1RZK8Cc0Ko4vprzMfavu44E1DzAtcRq//cXARaMvRuOjCFXcE8/wVs4Khm0rxTvAQFBqGoH3TcE/opCxm81QfLQSx0fvw8jwkfyS+wu+oRH4njORqmXLMaX0b9xmaK+hPD5ER0SJA+3l00gNaT5/UAjB5PjJxAfEc6tjNiHFVu5YYocliykRAu8zRhH++GMIIYjxj+Hq1Kt5f8f7XJx0McN6DyO3Opf7Kz9k4jXxXPFxNpNyJZm9JCvTjSw+K4fc1bcC8M9Zkxj67NdYMzP5Yahg8rgbWJW9ih8P/8iNg9q3MLfvuLFoQ0NxFBfjlZZ29H9skPJ3W8tBHotp4EA0/v6Yd+3C79xz8BrYvm6nQggCpk5BHx1F1o03cfjqaxAaDfrYWIKvuxaNwUDUq69y6PLLOXz1NcpOOh2hd/6VoJkzAQj/1784dNkMjtx4ExpfX0Kuv65d5zb06UPc55+T99C/KH7tdUrefgdhMmFMTiZwRuvrXgshCJk9G6+0NLLvvItDl19OwNSpWPcfIOqVlxtDWUGzZuGsM1P43HPk+foQMWdOuxPrKvXI+tpxd/wZPny4dAVmm11e9/56mfbv7+T27HIpN38oLWtek+Oe/Un2uf9r2fehb2VtWZksuPkcWbTsY9nn/q/lO6sOSIvNIa966n3528Mj5X/Xbj/l8zucDvlLzi/S7rBLWVUg5e+fSOl0Nt8ob5uUOVuaDf34+9ty4PyB8v0306ScN0nKp2OkfHu8zK/Kkc9veF7WWGtOfOKifdK58EZZcEm0rH1ggJQ5W6TNYZNzN70i0+YPlAPnD5Tff3GplGWHG3e57Yfb5MD5A+WqrFVHj/PpDCkf9Zfy0QApLUfPuSBjgRw4f6D84dAPsnbbNnnommulo7q6mQlXfH2FHP7xcJlblXtCU3cV75LnvXuGfGx2itz45RvSXl5+3DbV1mo5/ovxcuD8gXLIh0PkiE9GyDM/PVMeqTwi7VVV0mGxyC0FW+Sliy+VD655UC7LXCYfXPOgHDh/oFw/+zK5dchAec5bw2W5uVzO2z5PDpw/sNEup9MpP9n1iTxUcahVGwteeEFmpA2W9qrmr3HvWWfLwzfedMLX1xJZd/xV7uqfIut27znpfaWUsnbrVrl7eLrcldxfVq74qdlzdXv2yIply6V53z7ptFiO2zfnwQflruT+svD1uad0bvP+/TL3kUfl3rHjZM2mTe3ez5KVLQ9Mu0juSu4vD0y/WDodjuO2KXz1Vbkrub/Mf/qZU7KtJwJsku24xgrpxiv5pKeny02bNrnk3FJKzDYnXoajibM1+4q4et4GRsYH8+UtZzaOn/HUCs5MDGFkfDD//N92tBrB9CFRvDhjcEuH7kyj2fnNX4mvKMTbWqOsbzz5eWXC18mQvRm+vFpZJ+Hs+2HbF2ysymRRZF8eOLIfP+mE5MlQmcfqin1s8Q/izuvXI3QGsNXBs/FKWKkiC27+GSKV+G+FpYKnNzzNP0b8g2BTyyGE7UXbKbeUN7Z2PhEZJRl8sOMDHh39KD56nxa3yarMYm3uWgprCyk1lzItcRrDew9v9Zg2h41bf7yVrfmb8a2RXDjiSv4x4h8crjzMlK+mcP+I+7kq9SoW7l3IY78+xnl9zuPF8S827u+UTvaW7aV/cH+cViv2vDwMffo0O0flDz+g9fVtLOlsL5aDB7Hs3Yf/pPNPar+mmPfsoW7rHwTOuOyk7pYd5eWUL1xI0KxZjR5gV+GsqaH4rbfxO//8Fj0fKSUFc+ZQ9tnnxP33v3gNGtil9rkjQojNUsr0NrdzRwFoEgK6ad++fa42pxnzf8kkqZcfY/seLeO7Yf5GDhXXYLE7CfMzEupr5EBRNSvvG+86QztKdREsvF4JP/mGw/T/QNJEZSWxn5+G/SsgOBF8w2DnVzD9TRhyBez9Hj67TJkXsfx+mP4WDGmxVsBtqbBUcNW3V5Fdlc2yS5c1Llt5yZJL8NP78eTYJ7l0yaVYncrM7hWXrWgUtE8zPuWZDc/w5ZQvSQlJafUcKqcXR1UV+8dPwHfCBKJeeN7V5ric9gqARyWBTwfXjYlvdvEHGBDpz8HiGnLK67jn3H6kxwWRWVxDSXU3nrnpGwZXL4I/fwC3rVMu/qAsHXnRG3Dvbrj+G+X58EGw5kVwOmDvcmWRmWFXg9ZwtO1ENyLAGMCHkz/kkws/abz4A5wbey5bCrfw91V/RwjBaxNew+60s/TAUkDxHj7Y8QEAv+X95hLbPRWtnx+Bf/4zlcuWNa5mp9I2bikA3Y3UCH8AhsUGclbfUIb3URKsmw+3MHmoO6HVKdVHPiGtbyMEnPV3KNmveAJ7v4PECUpFU0jfo43nuhnBpuDjlpWc2GciEsmOkh3cP+J+xkWPY3DYYBbuXYiUkqUHl1JQW4CXzosN+RtcZLnnEnzN1SAlpZ984mpTug1uKQAd7QXU1aTHBRMX4s0/L0hBCMGgqAD0WsHmI91cANpL/6nKwjLLH1D6IfWrj1GHJXdLD6A1+gb2JTkomXNiz2F60nQALu17KYcqD7GpYBPzts8jNSSVaYnT+L3gd2xOm4st9iz0UVH4nX8e5V/+F0d1javN6Ra4pQC4cwioJcL8jPz89wmMiFPiwCa9loFRAWw+5CECoNHAuPugpr6Fdd/6yV29UqD8CFh7xpdRCMHnF37OC2e/0JhAPT/ufHz0Pvxr7b84UnWEmwbdxIjwEdTaa5utl7AxfyNLDrTSPkHltBFy/fU4q6qo+N//udqUboFbCkBPIL1PENtyKrDYPWSCysBLlNnMUcPBrz5uHlZffVS0x3V2nWb0Wj3aJgsHeeu9uSD+AnJrcokPiOdPsX8ivbeSe2sIAzmcDh755REeW/eYWy1U0xPxSkvDa9gwyha0sG6GynGoAtBJDO8ThNXuZEdOpatN6Ro0Wrh2KVz+6dGx+tXGWmqD3ZOYkTwDjdBwS9otaISGEK8QkgKT2Ji/EYA1OWvIrs7G6rSyJmeNi63t+fiMGY01MxOn2exqU9weVQA6iWH1ieDfu3si+GTwjwT/ox1ECYpXKoEKe04eoCX6B/dn5YyVXJhwYePYiPARbCncgs1h45OMT+jt3ZtgUzArDh/fTkLl9GJMTAQpsWa6z6pw7opbCkB3SwK3RC8/E7HB3mw+XIbd4WR3fiU55XWuNqtr0eqUSqAe7gEAx01sGxE+gjp7HYsPLGZ93npm9p/JhJgJrMpehcXRjcuDuwGG+p5Mlv2uWRWuO+GWAtDdksCtMbxPED/tLmTQv79n0itruPJdD6wN79W/25aCdoSGPMDzG5/HqDXy575/ZmLsRGrttazPW+9i63o2hrg40GiwHFQFoC3cUgB6CpcNjyY9LojLR8QwJS2CQyW1lNUcvy5wjyYsBSqOgKW67W17EEGmIPoF9aPWXsuUhCkEmgIZFTEKX71vi11FVU4fGoMBQ0wM1gMHXW2K26MKQCcyOimUz246g39PG8DlI5QVunbleUhSuIGGPkTFPacSqL2MDFf69V+RoiyIZ9AaOCv6LFYeWYndaXelaT0eQ1KS6gG0gy4TACGEjxDiQyHEu0KIK9veo2cxIFIJZ+3M7b55jVMirL4S6L/XwzsT4MNpUJ514n16CLMHzeb1P71Ov6B+jWMTYydSZiljS+EWF1rW8zEmJGA9fARpUyfjnYgOCYAQ4n0hRKEQYscx45OEEHuEEPuFEA/UD18CLJRS3gRM68h5uyPBPgYiAkzszPUwDyA4QVlroPdA8A6BQ2th47uutqpLCPUKZXzM+GZjY6PGYtQa+eHwD64xykMwJCaAzYY1yzNuNk6VjnoA84FJTQeEEFrgDWAykArMEkKkAtFAw6fhIbOjmjMg0t/zBECjgcnPwKzP4KqF0G8SbP1MWczGA/HWezMqYhRrc9a62pQejTExEQDLATUMdCI6JABSytVA6THDI4H9UsqDUkorsAC4CMhGEYETnlcIcbMQYpMQYlNRUVFHzHM7UiMDOFhUTZ3VI/VPYdg1SsuIvd+52hKXMTpyNFlVWWRVqnennYUhXikFVRPBJ6YzcgBRHL3TB+XCHwX8D7hUCPEmsLS1naWU70gp06WU6WFhYZ1gnusYEOmPU0JGvod5AU1JOkdZX2DLx662xGWMiRwDwLrcdS62pOei9fVBFxGhJoLboDMEoKVlhqSUskZKeb2U8jYp5actbHP0AD1gIlhLDIhU2kZ7XBioKVqdsnDMvu+hMtfV1riEPv59iPSJ5JfcX1xtSo/GmJCgegBt0BkCkA3ENHkcDXjmN/0YogK9CPDSs8vTKoGOZehVIJ1KLsADEUIwOmo0G/I3qC2jOxFDYgKWgweRTqerTXFbOkMANgJ9hRDxQggDMBNQ++CifPE9MhF8LCGJEDdOCQNZPbM75pjIMdTYathWtM3VpvRYjAmJyLo67OoKYa3S0TLQz4FfgWQhRLYQYraU0g7cAXwHZABfSil3nsxxe0oriJYYEOnP7vwqbA4PvysZdQuUHYLXh8HmD8HhWROjRkaMRCu0ah6gEzEm1VcCHVTDQK3R0SqgWVLKCCmlXkoZLaWcVz/+rZSyn5QyUUr55Mket6fmAECZEGa1OzlQ5FmtEY4jZSpcv0xZY3jpnfDBZPAgV93f4M+g0EGsy1EFoLMwqKWgbeKWrSB6ugcAsNNT1gk4EX1Gw+wf4NzHIXsDZHlWs7zRkaPZWbKTcnO5q03pkeiCgtAGBamJ4BPglgLQkz2AhDBfTHoN//d7Nq+v2MdL3+/xbG9ACEifDTov2P5fV1vTpYyOGo1Esmj/InWlsE7ClNKf2k2bkFK62hS3xC0FoCd7AFqN4MyEENYdKOHFH/by2k/7eXaZ57VLbobRF5Inw85FHjVDeEDIAMK8wnhx84uc+fmZXLb0MrYWbnW1WT0Kv8mTsWZmYt6xo+2NPRC3FICeznvXjmDnY+ez/8nJXDc6jp/3FlFl9pwLX4sMugzqSuHgz662pMvQaXQsmr6I/0z8Dzen3UylpZK7Vt5Ffk2+q03rMfhPmoQwGqn4apGrTXFL3FIAenIICBQvwMeoQ6fVMHVwBFa7kx8zClxtlmtJmgimAI8LA/kb/BkXPY6/DPkLb57zJhaHhbtW3oXZrq5nezrQ+vnhN3Eild98g7R62Foc7cAtBaAnh4COZWhMEJEBJr7+w8NrlXVGSJkGu79R5gbUlMCi25USUQ8hITCBp8c+za6SXcz5bY4atz5NBEy/CEdFBdWrV7vaFLfDLQXAk9BoBBcMimD1viIqaj09DPRnsFbDyifhrTGw9VP46QmPmiMwIXYCtw+5nSUHlvD2trddbU6PwGf0aLRhoZQvUsNAx+KWAtDTQ0DHMmVwJDaH5PtdHh77jRsHvr3h17lg8IEJD0FNIRxc6WrLupRb0m5hasJU3tj6Bl/u+dLV5nR7hE5HwJSpVK9ajb2szNXmuBVuKQCeFAICGBwdQEywF19vU8JABZVmFm3JweH0sBCARgvnPAZn3gE3r4Ixd4NXsMf1DNIIDY+NeYxxUeN44rcn1MVjTgMB0y8Cm43Kb751tSluhVsKgKchhODCQZH8sr+Ye7/8g7HP/sTdX2zlp92Frjat6xkyC85/UikN1RmUsNDub6DOsyZL6TV6Xhz/IoPDBnP/6vv5NfdXV5vUrTElJ2OIj6d6jZoHaIoqAG7ClLQI7E7J19tymTkiFqNOw28HS1xtlusZPBMcFtjlefFbL50XcyfOpY9/H+5aeZfaOK6DeI8YQd3m35EOD16Q6RhUAXATBkYF8MXNZ7DugT8xZ/pAhsYGsj5TFQAih0FoP/hjgastcQkBxgDeOfcdQr1Cue3H29hXts/VJnVbvEek46yuxrJnj6tNcRvcUgA8LQncwKiEEEJ8jcrf8SHszK2kos7DK4OEgMGz4MivUOqZPV3CvMN459x3MGlN3PzDzWRWZLrapG6Jd3o6ALWbNrnYEvfBLQXA05LALTEqIRgpYdOhY5dc9kDSLgcEbHjX1Za4jGi/aN49712c0sns72ZzuPKwq03qdugjItBHR1O7URWABtxSAFRgWGwQBq2G9ZmqABAQBcOuhvVvQ8FJLS3Ro0gITOC9897D7rRzw3c3qIvKnwLe6elqc7gmqALgppj0WgbHBLBeTQQrnPOY0iri67951LoBx9I3qC/vnvcuFoeFG76/gawqVQROBu8R6TjKyrCqi8QAXSgAQogEIcQ8IcTCrjpnd2dUfAg7ciuptnjOTNhW8Q5W1g3I+g3+8Kx5AceSHJzMe+e9R529TvEEVBFoN415gI0bXWyJe9AuARBCvC+EKBRC7DhmfJIQYo8QYr8Q4oETHUNKeVBKObsjxnoaZySE4HBKNQ/QwJArIWYUfP8wrHoOVjwOPz8LdourLety+gf3591z31VF4CTRx8ai69VLzQPU014PYD4wqemAEEILvAFMBlKBWUKIVCHEICHE18f89DqtVnsIw/oEotMINQ/QgEYDU15W/l75JKx9GX5+yqMaxjUlJSSlUQSuX369mhhuB0IINQ/QhHYJgJRyNXDsVWgksL/+zt4KLAAuklJul1JOOean3VNahRA3CyE2CSE2FRUVtfuF9ES8DTrSogPUCWFN6T0A7tsLD5fAI6UQcwb88grYPbPVb0pICvPOm4fVYeX65ddzsFyNbbeF94h07AUF2LKzXW2Ky+lIDiAKaOp3ZtePtYgQIkQI8RYwVAjxz9a2k1K+I6VMl1Kmh4WFdcC8nsGohBC2ZVdw/8JtLN+RT42aDwCtHrQ6ZY7A2X+HypzmeQFrDVR7ThuN5OBk3j//fZzSyfXfXc/esr2uNsmtOTofYLOLLXE9HREA0cJYqz6VlLJESnmrlDJRSvn0CQ/soRPBWuKGMfFMGhjOt9vzuPWTzUx5fS1Wu+dWwRxH4kRltvCal5S20aUH4c0x8NY4j8oNJAUl8cGkD9AJHTd8dwO7Sna52iS3xZCQgDCZ1BnBdEwAsoGYJo+jgdyOmaNyLGF+Rt64Yhi/P3Iuz12aRmZxDUv+UN/mRoSAs/4O5Ydhxb9h3nlQlQfV+UoTOQ8iPiCe+ZPm463z5sbvbuSPoj9cbZJbIrRaDAnxWPbvd7UpLqcjArAR6CuEiBdCGICZwJLTY5bKsei1Gi5Lj6Z/uB/vrj6oJrCakjwZeg+Cda+D3gtuWQ0BsbB5vqst63Ji/GP4cNKHBJoCufn7m9mYr5Y7toQxKUkVANpfBvo58CuQLITIFkLMllLagTuA74AM4Esp5WmZpqm2gmgZIQQ3jUtgT0EVq/Z6doK8GULA5GchdTrM/gHCkmH4NZC5CkoOuNq6LifCN40zGpAAAB5VSURBVIL5k+YT7hPObT/extqcta42ye0wJvXFnp+Po6rK1aa4lPZWAc2SUkZIKfVSymgp5bz68W+llP3q4/pPni6j1BxA60wdHEm4v4l3VqvVHs2IGwMzPgS/cOXxkKtAaOH3j1xrl4vo5d2LDyZ9QEJAAn/96a/qojLHYExKAvB4L8AtW0GoHkDrGHQarh8Tx7oDJWw+XMqafUU8u3y3OlnsWPwjlNDQ1k89tkQ02BTMe+e/x8CQgdy36j4W7fe8NRVaw9hXFQBwUwFQPYATM2tULL5GHZe++StXz9vAmz8f4NUVap/44xh+HdQUwd5lrrbEZfgb/Hn73LcZFT6Kh395mE8zPnW1SW6BPioK4eWFVRUA90P1AE6Mv0nPkxcP5PoxcXxw3Qhmjohh8+EybA61PLQZiX+CgBj47S3w4KS5t96buRPnMjF2Is9seIa3/3jb44sIhEaDMSEByz5VAFS6IRcNieLRqQOY0L8XZ/ULo9bqYHuO6jE1Q6OFMXfBkXWw/0dXW+NSDFoDL5z9AlMTpjJ361ye3/Q8TunZNwxqJZCbCoAaAjo5RsYHA6gtI1pi2LUQFAc/PubRbaQBdBodT4x9giv6X8HHuz7mkV8ewe703Jnlxr5J2AsLcVRWutoUl+GWAqCGgE6OUF8jfXv5sv6gmgg+Dp0B/vQwFGyHHfWdyEsPwuI7oDDDtba5AI3Q8MDIB7h9yO0sPrCYe3++F4vDc2ZMN8WgVgK5pwConDyjEoLZdKgUu5oHOJ4Bl0B4Gvz0hLKs5P+3d+dxVZb5/8dfn8POwQUQFMEFXDDcFcQlG5UWTc2ZqaYsTUuznBydyr5Wzrd+TdP0aGoqK8cZ01abrGn/OqZOri2GopZaipomoiagiCDrgev3x31siFDRw+G+z32u5+NxHnLuDvC+u+F8uJb7uhZcCtteN1YU9UMiwvTe03lgwAOsObSG6Z9Mp6SyxOxYTS6kcxcAvx4HsGQB0F1AFy49MZrTldV8c8R/m7Nn5XDA5Q8by0Usnw3tBhj3CexeDkWHzU5nmpsuuYnHhz7OtmPbuG3lbRwv868uxKC2cUh4uG4BWI3uArpw6UnGOEDmAf/6JW6wThnwizkwdh5MfN9YRVTVwJaXzU5mqjFJY5g3Yh4Hig4wacUkcov9Z4lkcTgI6dSJin3+O4XakgVAu3CxzUJJinHypR4HqJ8IDH/QuDdAxBgY7nqVsZmMn94odsZlCZfx4pUvUlheyC0f30L2Cf9ZJdPfZwLpAmAj6YnRbD5wguoa/57j3WBpU+F0HuzSaxj2ie3DqyNfRUS4dcWtbDnmH2vlh3TuTHV+AdUnT5odxRS6ANjIwKQoiitcfKvHARqmUwZEJsLmRWYnsYTOkZ1ZMmoJ0WHRTFs1jdU5q82O5HUhXY2B4P1jr+HgpMnkP/8Cqrra5FRNx5IFQA8CX5z0xGgAfvm3zxn6lzVMfTWLvFPlJqeyMIcD0qZAzkY4vNXsNJYQFxHHa6Neo1tUN+5Zdw9vZ79tdiSvcqanEztnDs4hQ6guPkXB/PmUbNhgdqwmI1a+JTw1NVVlZWWZHcOnrNl9jG05Jzl4vJTlO44yeXBH/jAmxexY1lVeBM/1hZhuMPnfxviARmlVKbPXz+bTw58yvfd0pveejtj8/42qqmLvsOGE9e1DuxdeMDuOR0Rki1Iq9Xyvs2QLQLt4I7q15t4rk3lufF+u6tGGf23JpbzKf5q0Fyy0BQyfCwc/h28/NDuNZYQHhTNvxDzGdRrHgq8X8MjGR2x/17AEBdFi3DhK1q3HVVBgdpwmoQuAjd2c3p6isiqWbT9qdhRr6z8ZWveAVf8LVWXGscNbYM8qU2OZLcgRxKNDHuX2nrfz7t53uXvt3ZS5ysyO5VUtr7sWXC6KPvSPiQFNWgBE5Jci8qKIfCgiVzbl9/ZHg5Ki6RTjZMmXB388VumqIa9Yjwv8hCMARj4ORTlGEXhrArw4Av75Gyj83ux0phIRZvabydz0uazPXc/UlVM5UW7fqcYhSUmE9e3LyXff9YsVUxtcAETkJRHJE5GddY6PFJFsEdknIvef62sopT5QSt0OTAZuuKjEWoOJCDend+CrQyfZebiIvOJyrl3wBRl/Xa+7hepKvAy6jYHNL8J3a41VRB0BxtIRGjd2u5Fnhj9DdmE2E5dP5NCpQ2ZH8pqW111L5f79lG37yuwoXnchLYBXgJG1D4hIADAfGAWkAONFJEVEeorIsjqP2Fqf+gf352ledm2/BEKDHPx1VTa/mv8FOw4XUVzuYqdeOvrnxjwDVz4Gs76GK/4IKeNg6+tQ4d/7xp6R0T6DRVcuoqiyiAkfT2Bnwc7zf5IPanbVSCQ8nJPvvmN2FK9rcAFQSm0A6rb9BgD7lFL7lVKVwFJgnFJqh1JqTJ1HnhieAD5WStU7705EpolIlohk5efrjc891SI8iLG92rI2O5/K6hpevjUNgK05hSYns6CIWBg8A5ytjOfp06GiCL5609xcFtIntg+vj3qdsMAwblt5G+sOrTM7UqMLiHDSYvRoit57n8Oz76MyJ8fsSF7j6RhAPFC7LZjrPnY2vwMuB64TkTvre4FSaqFSKlUplRoTE+NhPA1gxojOXNc/gfd/O5jhybG0jwpn60H/vPPxgrRLg/hUyPy73+8lUFtii0SWXL2EpBZJzFo7i7d2v2V2pEYXO2cO0VOnUvzJJ3x39WiOL15sdiSv8LQA1Dcx+KwjJ0qp55RS/ZVSdyql/n7WL6pvBGtUHaKdPHV9bxIiwwHo3yGSLTmFfjHI5bGB0+HEd36/o1hdrcJa8dJVLzE0fih/yvwTT2952lY7jAVEOIm99x46rVqJc0AaBfP/Rk2l/daM8rQA5ALtaj1PAI54+DU1L+vXviX5xRXkFtp7Sl+jSBkHzeLg83l+va9wfcKDwnl2+LPckHwDL+98mTkb5thuc5mg2FgiJ06kprSU0k2bzY7T6DwtAJuBLiKSKCLBwI2Af0yg9WF920cCehygQQKCYOi9cPAz+PYDs9NYTqAjkLnpc7m7/92s+H4F01ZN42S5vboXnYMGIaGhlKxda3aURnch00DfBDYCySKSKyJTlFIuYAawEtgFvK2U+sbTUHo/AO/q1qYZ4cEBbMux1y+q1/S/FVr3hJVzocL/ds46HxHhth638eRlT7KzYCcTPp5Azin7DJw6QkNxDhpEydq1tus2vZBZQOOVUnFKqSClVIJSarH7+HKlVFelVCelVKPssafHALwrMMBB74SWbDmoWwANEhAIo5+CU4dhw5PGMVcF7FkJZbqInjEycSSLrlpEUUURE5ZP4Ks8+8yjjxg+jKojR6jYY6/NYyy5FIRuAXhfvw4t2XX0FGWV+oawBmk/EHrfBBvnw4oH4OkU407hVXPNTmYpfWP7suTqJTQLbsaUlVNYcWCF2ZEaRcQvhgHYrhvIkgVAtwC8r1/7SFw1iu25+i/YBrviEQgKhy8XQEIadB0J29+G4mNmJ7OUDs07sOTqJXRv1Z37NtzHoh2LfL7rJKh1LKHdu+sC0BR0C8D7/jsQrAtAg0XEwh3r4Pfb4aalxl3D1VV6Q5l6RIZG8uKVLzIqcRTzts7j4S8epqq6yuxYHokYPpyy7dtxHTf23VY2uDfEkgVA874oZzBJrZz8c9NBfvOPjaT+6T88vnyX2bGsLyoJWrY3Pm7VGZJHGQWgSk+prSskIIQnhj7BHb3u4P1973PnJ3dSVOG7rfqI4cNAKY499hgHJ00mu3cfCt/27Q1zLFkAdBdQ0xjdK46Schc1NYooZzD/zMzRi8RdqEF3QdkJ+FovF1EfEWFG3xk8duljbM3byoTlE3x2NdHQlBSC4uM5tfxjqgsLcUREULzat7fN1DuCaQCs35PPpJc2sXBif67s3sbsOL5DKVg4DCpPw12bjG0mtXplHs1k6qqpzE6dzaTuk8yOc1Gq8vJAKYJat+boQw9zavlyumZ+iQQEmB3tJ/SOYNoFGdwpmsjwIP6vzuYxuYWl1NRY948E04nAoBlwfC/s0juKnUt6XDrJkcmsyVljdpSLFhQbS1Dr1gCEp6VRU1JC+e7dJqe6eJYsALoLqOkFBTgY1TOO1buO/Tg19It9BQz9y1r+b7te3eOcuv/K2FFsxYN66ejzyGifwba8bRSU+f6Wi+Fpxh/YZT7cS2HJAqBnAZljTK84SiurWbM7jwpXNX/4YCdKQdb3+oaxcwoIhDHPQvFRWPtns9NY2oj2I1Ao1h7y/emUQW3aENSuHac3++4aQZYsAJo50hOjiWkWwrLtR/jH+v3sLzhNq4hgtuvNY86vXRqk3mosHX3EfQdswT7I9d2/Dr2ha2RX2jVrx+oc3x48PSM8NZWyrC0+OyVUFwDtRwEO4eoebVi9O48X1u5jdK84ft0vgV1HT1Hp8s0f8CaV8TCEt4J3boO/DYIX+sOiy41CoAHGrKCM9hlkHs2kuNL3u8vC09KoPnmSin2+eY11AdB+YkzvtlS6aggOcPDQmBR6xreg0lXDnmO+/8vqdWEt4eq/GBvJh0Ua20oGBMPG581OZikZ7TNw1bj4NPdTs6N47Mw4QKl7HMBVUMDxxS+hXC4zYzVYoNkB6iMiY4GxnTt3NjuK3+nfPpKMbrGM7hVH6+ah9EowxmF2HC6iR7wekzmv7r+CbmONcQGAE/uNLSWHzzXuJNboFdOLVmGt+CTnE3rF9GL5geUopbij9x1mR7tgQQkJBLZpQ+nmzbQYO5acqbdTsXs3od274xyYbna887JkC0APApvH4RAWT07j1/0SAGgfFU7z0EC25+pxgAYLqPV31aDfQXUlbFpoXh6LcYiD4e2GszpnNaPeG8Xz257nha9eoNxVbna0CyYihKemUro5i9y7ZvzYFVS+c4fJyRrGkgVAsw4RoVdCS3Yc1msGXZRWnaHbaGO5iMrTZqexjBuSbyC1dSqz+s1iVr9ZAPxw+geTU12c8LQ0qgsKKN20ibaPP05QfDxlO3aaHatBdAHQzqtXQguyfyjWy0RcrMEzoawQti0xO4llJEcls/iqxUztOZXeMb0BOHr66Hk+y5qcQ4YgISG0fuB+WowdQ2jPnpTv1AXgJ0TkEhH5u4i8IyLTm+r7ap7rldCCqmpF9g96IPiitE+HdgONfYX1jWI/E+eMA3y3BRCcEE9y1maiJhnLW4T17EHV4cO4Tlh/zaMGFQAReUlE8kRkZ53jI0UkW0T2icj95/oaSqldSqk7gd8A512jQrOOngktAfT9AJ648lE4dQTW/MnsJJbT2tkaQXy2BQAgQUE/fhzaoyeAT7QCGtoCeAUYWfuAiAQA84FRQAowXkRSRKSniCyr84h1f841wGeAPe4C8RNtW4QS7Qxmh9485uK1GwADbofMf8Ah371z1BuCHEHEhMf4dAGoLbR7CohQtsP6A8ENKgBKqQ1A3fbMAGCfUmq/UqoSWAqMU0rtUEqNqfPIc3+dj5RSg4Gbz/a9RGSaiGSJSFZ+fv7FnZXWqESEngkt2J5bhFKK9XvyeeXzAz6/y1OTy3gImreFj34Hrko4fRz2fgKl1u8q8LY4ZxxHS+xRAAIiIghOTKS81kBwTWUl1SXWmwTgyX0A8cChWs9zgbNOfBWRYcCvgRBg+dlep5RaCCwEYzloD/JpjahXfAs+3VvAmOc/45sjpwAYkBhNStvmJifzISHNYPTT8OYN8GwPKHFvJXnJWLjBvweI45xxfHv8W7NjNJqwnj0o+eILlFKICEfunU3loUMkffC+2dF+wpNBYKnn2FnfsJVS65RSM5VSdyil5p/zC+vVQC0ntWMU1TWKsqpq/jD6EgDW7ckzOZUPSh4JQ2ZB277G0hF9J8CuZVCw1+xkpopzxvHD6R+oUfZYciS0R0+q8wtwHTtG6bZtFP/nP1Ts3k3NaWu1AjxpAeQC7Wo9TwD0usE2NbRLK1bdfRmdYiIIcAjvbj3Muux8fjtM3619wa74438/LsmHHe8YM4TGvWBeJpPFRcRRWVPJifITtAprZXYcj4X17AFA2Y4dFL7+39Zdxd69hPXpY1asn/GkBbAZ6CIiiSISDNwIfNQYofSdwNYjInRt3YwAh9HwG54cw5aDhZwq9+2Nvk0XEQN9J8LXS41ZQn7qzFRQu4wDhHTrBoGBnFi0mNJNm4iaPBmA8uw95garo6HTQN8ENgLJIpIrIlOUUi5gBrAS2AW8rZT6pjFC6S4g6xuWHEt1jeLzvb6/sYfpBs8AVQMbz9kzams/FgCbzARyhIYS0qULZV9/TVDbtsTc/XscTicV2dlmR/uJhs4CGq+UilNKBSmlEpRSi93HlyuluiqlOimlHvNuVM1K+rVvSbPQQNZl65laHovsCD2uhS2v+O2MoDZOYx9quxQAgLAeRjdQqxkzcISEEJKcTPkeHywATU13AVlfYICDoV1asX5Pvp4O2hiGzILKElj/F7OTmKJ5cHPCA8N99m7g+rS87loibxpPi3HXABCS3JWK7D2W+n2xZAHQXUC+YVjXWH44Vc5uvUSE59r0gLTbjR3FcjLNTtPkRMS4F8BOLYDevWnz0ENIQAAAocnJ1BQX4zpqnXO0ZAHQLQDf8IvkGADdDdRYLv9/0KIdfHgXVJWZnabJxUXEcaTEvgPhIV2TASi30DiAJQuA5htaNw8lJa45H351mO8LrDW/2SeFRMA1z8Hxvcbm8kW5sH897F9ndrImceZeALsK6doFgAoLzQSyZAHQXUC+Y/LgjuzNK2HYU+uYuDiTrTmFZkfybZ2GQ79b4Ivn4Jnu8No18No4OLTJ7GReF+eMo7CikDKXPVs/ARERBCUkUJ692+woP7JkAdBdQL7jN2nt+OL+EdxzRVd2/1DMHa9v0fsGeOqqP8OwB41lIya8B2FRsOEps1N53ZmZQLZuBSQn6xaAZi+tm4cyM6MLz4/vS35xBUs35ZgdybeFNINhcyBtCnTOgIHTYe9KOLrd7GReZbd7AeoTmtyVyu+/p6bcGttfWrIA6C4g3zQwKZoBHaNYsP473QpoTANuh+Bm8NnTZifxqrYRbQH73A1cn5CuyVBTQ8W+78yOAli0AOguIN81M6MLx05V8K8tuWZHsY+wSKM18M0Htl40LiY8Boc4bN0CCEnuCmCZO4ItWQA03zWkczT92rdkwdp9VLpqOHSilOU7jlKs1wzyzKC7IDAEPnvG7CReE+QIIibMPhvD1Ce4fXskNJSy7dt/dkOYqm76VrMnq4Fq2s+ICDMzujD55c0Menw1x09XAnDPFV2ZmdHF5HQ+LCIW+t8Km/5h/NsuzexEXhHnjONQ8aHzv9BHSUAA4ampnHzrLUozM2k+ejQ15WWUfplJ+a5dxM97luZXXNFkeXQLQGt0v+gaw3X9E+jfIZKHx6bQOTaCjd8dNzuW7xv+IDSPhw/uhMpSs9N4RWqbVLblbeOJTU/YZm+AuuKffYY2j/6RwNhYCubPp/C113GEh+OIiKB4xcomzWLJFoCIjAXGdu6s15r3RSLCU9f3/vH5oRNlvJF5kApXNSGBASYm83GhzY09A14bB2sehZGPm52o0c3oM4MyVxlLdi0hrzSPPw/9MyEBIWbHalQBERFEXn89kddfj+vECRxhYTjCwjgyZw4lGz5F1dQgjqb529ySLQA9CGwv6UlRVLhq2J6rZ3V5LGmYsWbQlwtg6+vw7Uew7Q04aY+ptwGOAOakzWF26mxWHVzFtFXTKKqw789NYFQUjrAwAJyXXkp1YSHl3zTd1piWLACavaR1jAJg0wH/XOq40V3xiLGE9Ecz4O2J8OFv4Z0pYKFVJj0hIkzqPoknL3uSHQU7mPjxRA6XHDY7ltc5Bw8G4PTnnzfZ99QFQPO6KGcwya2b8eV+PQ7QKIKdcMd6mPxvuONTYxG53E1wYL3ZyRrVyMSRLLxiIQVlBdz875v55nij7DdlWYHR0YSkXMLpzz5rsu/ZpAVARJwiskVExjTl99XMl54UxZaDhVRV23Ngr8mFtoCOl0JcLxj4W2jW1pZ7CaS2SWXJqCWEBIRw64pbWX/IXkWuroghl1L61VdUlzTN4ooN3RLyJRHJE5GddY6PFJFsEdknIvc34EvNAd6+mKCab0tPjKa0spqdh+3bn2uawBBjQ5mDn8P3TffXY1NJapnEG6PfILFFIjPXzmTp7qVmR/Ia55Ah4HJRuqlp9oRoaAvgFWBk7QMiEgDMB0YBKcB4EUkRkZ4isqzOI1ZELge+BY41Yn7NRwxINMYBMvU4gHf0nwTOWFu2AgBahbXi5ateZmj8UB7LfIy/Zv3VltNEw/r1RcLDOf1Z04wDNHRP4A1A3d/cAcA+pdR+pVQlsBQYp5TaoZQaU+eRBwwHBgI3AbeLiB5/8CMxzULoFOMkU48DeEdQGAyZaYwD2HRHsfCgcJ4d/iw3Jt/IK9+8wuz1syl3WWNRtcbiCA7GmZZGyedN05Lz5E04Hqh9y16u+1i9lFJzlVK/B/4JvKhU/eVbRKaJSJaIZOXn652m7CQ9KZqs7wuprrHHbBXLSb0NmsXB0a/NTuI1gY5AHkx/kPtS7+OTg58wZdUUjpfZ648K56WXUnUwh8pD3r8j2pMCIPUcO+9vtlLqFaXUsnP894XAI8DW4OBgD+JpVpOeGEVxhYtvj5wyO4o9BTth5jZIn2Z2Eq8SEW7pfgtPD3ua3OJcjpfbrQAMIaxPH6pPen+8zJM7gXOBdrWeJwD23dBT89jApGguvyTW7Bj2FhRmdoImc3mHyxkSP4SwQHudc0hiIh2Xvtkk30vqrkh31heKdASWKaV6uJ8HAnuADOAwsBm4SSnVaJN1U1NTVVZWVmN9OU3TNL8gIluUUqnne11Dp4G+CWwEkkUkV0SmKKVcwAxgJbALeLux3vz1hjCapmne16AuIKXU+LMcXw4sb9REmqZpWpOw5FRMvRicpmma91myAOguIE3TNO+zZAHQLQBN0zTvs2QB0DRN07zPkgVAdwFpmqZ5nyULgO4C0jRN874G3whmBhHJBw5e5Ke3AgoaMY4V6XO0B32O9mClc+yglIo534ssXQA8ISJZDbkTzpfpc7QHfY724IvnaMkuIE3TNM37dAHQNE3zU3YuAAvNDtAE9Dnagz5He/C5c7TtGICmaZp2bnZuAWiapmnnYMsCICIjRSRbRPaJyP1m5/GUiLQTkbUisktEvhGRWe7jUSLyHxHZ6/430uysnhKRABHZJiLL3M8TRSTTfY5viYhPbxMnIi1F5B0R2e2+noPsdh1F5G73z+lOEXlTREJ9/TqKyEsikiciO2sdq/e6ieE59/vPdhHpZ17yc7NdARCRAGA+MApIAcaLSIq5qTzmAu5VSl0CDATucp/T/cBqpVQXYLX7ua+bhbG/xBlPAM+4z7EQmGJKqsYzD1ihlOoG9MY4V9tcRxGJB2YCqe7NowKAG/H96/gKMLLOsbNdt1FAF/djGrCgiTJeMNsVAGAAsE8ptV8pVQksBcaZnMkjSqmjSqmt7o+LMd404jHO61X3y14FfmlOwsYhIgnAaGCR+7kAI4B33C/x6XMUkebAZcBiAKVUpVLqJDa7jhj7jIS5dw0MB47i49dRKbUBOFHn8Nmu2zjgNWX4EmgpInFNk/TC2LEAxAOHaj3PdR+zBffWnH2BTKC1UuooGEUC8PUNd58F/geocT+PBk66d58D37+WSUA+8LK7m2uRiDix0XVUSh0GngJyMN74i4At2Os6nnG26+Yz70F2LABSzzFbTHUSkQjgXeD3SqlTZudpTCIyBshTSm2pfbiel/rytQwE+gELlFJ9gdP4cHdPfdz94OOARKAt4MToEqnLl6/j+fjMz60dC0Au0K7W8wTgiElZGo2IBGG8+b+hlHrPffjYmaal+988s/I1giHANSLyPUa33QiMFkFLd1cC+P61zAVylVKZ7ufvYBQEO13Hy4EDSql8pVQV8B4wGHtdxzPOdt185j3IjgVgM9DFPesgGGMA6iOTM3nE3Re+GNillHq61n/6CJjk/ngS8GFTZ2ssSqkHlFIJSqmOGNdsjVLqZmAtcJ37Zb5+jj8Ah0Qk2X0oA/gWG11HjK6fgSIS7v65PXOOtrmOtZztun0E3OKeDTQQKDrTVWQ5SinbPYCrgT3Ad8Bcs/M0wvlcitGE3A585X5cjdFHvhrY6/43yuysjXS+w4Bl7o+TgE3APuBfQIjZ+Tw8tz5AlvtafgBE2u06Ao8Au4GdwOtAiK9fR+BNjDGNKoy/8Kec7bphdAHNd7//7MCYEWX6OdT30HcCa5qm+Sk7dgFpmqZpDaALgKZpmp/SBUDTNM1P6QKgaZrmp3QB0DRN81O6AGiapvkpXQA0TdP8lC4AmqZpfur/Ax4bHQOQPWNdAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#<\n",
"for eps in [.1,.01,0.001,1e-4]:\n",
" x,G = gradient_armijo(f,gradf,.25*np.ones(2))\n",
" plt.semilogy(G,label='eps=%g' % eps)\n",
"plt.legend()\n",
"#>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**QN3** Écrire une fonction `hessf(x)` calculant $\\mathrm{D}^2 f(x)$. *(On pourra utiliser `verifier_hessienne()` pour vérifier le calcul)*"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"erreur numerique dans le calcul de la hessienne: 1.05158e-21 (doit etre petit)\n"
]
}
],
"source": [
"#<\n",
"def hessf(x):\n",
" return np.eye(2) + eps*(np.diag([1/x[0]**2, 1/x[1]**2]) + 1/((1-x[0]-x[1])**2)*np.ones((2,2)))\n",
"#>\n",
" \n",
"verifier_hessienne(gradf,hessf,np.random.rand(2)/3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**QN4**. Implémenter la méthode de Newton avec rebroussement d'Armijo (appelée *(Newton)* dans le sujet). La fonction Python `x,G = newton_armijo(f,gradf,x0)` retournera la solution calculée `x` et un tableau contenant $(\\|{d^{(k)}}\\|)$. En utilisant cette fonction, tracer en échelle logarithmique $k\\mapsto \\|{d^{(k)}}\\|$ en partant de $x^{(0)} = (\\frac{1}{4},\\frac{1}{4})$ et pour plusieurs choix du paramètre $\\eps \\in \\{.1,.01,.001,.0001\\}$."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xd4VGXa+PHvPS2dhAQSSCahBUgokRIIit11F3dFBBtYd2XFhruvW151f7uurrq6vusWBXfFslgpYgFc164LFkoo0nuRUNIT0pOZeX5/TIIhUgJzMudM5vlc11xhTs48506Auec891NEKYWmaZoWnmxmB6BpmqaZRycBTdO0MKaTgKZpWhjTSUDTNC2M6SSgaZoWxnQS0DRNC2M6CWiapoUxnQQ0TdPCmE4CmqZpYcxhdgDHIiLjgfFxcXG3DBgwwOxwNE3TQsqqVatKlFLd23OuWHnZiNzcXJWfn292GJqmaSFFRFYppXLbc67uDtI0TQtjlkwCIjJeRGZVVlaaHYqmaVqnZskkoJRarJSaFh8fb3YomqZpnZolk4CmaZoWHDoJaJqmhTGdBDRN08KYJZOALgxrmqYFhyWTgC4Mh4C6Clj9MpTtMjsSTdMCYMkZw5qFeRpg5XOw5P+grhwcUXDR/ZB3K9jsZkenadopsuSdgGZBPh+sex1m5ML7v4HU4XDDW9D3PHj/PnhhHBRvMztKTdNOUdCSgIj0FZHnRWRBsK6pGWTXf+HZ8+HNn0JEPFz/pj8B9LsQpsyFSc9C6Xb459mw9C/g9ZgdsaZp7RRQEhCRF0SkSEQ2tDk+TkS2isgOEbkXQCm1Syk1NZDraUFWuBFeuRJeugxqy2DiLLh1CWRe9O05IpBzNdy5AgaOg48fhOcugkMbjt+upmmWEeidwGxgXOsDImIHZgKXAIOAKSIyKMDraAapqG3k/721noVr93PcxQMr98Pbd8I/xkLBCrj4IZieD2dcA7bj/JOJTYarX4KrXoTD+2HWefDpo+Bp7LgfRtO0gAVUGFZKLRGR3m0OjwZ2KKV2AYjIXGACsCmQa2mB23qoimkv57O3tJZXl3/DexsO8fDlQ0iKjfCfUF8Jn/8Vlv0DlA/Omg5n/wKiE9t/kcGXQ59z4b174b+PwebFMGEGpI3omB9K07SAdERNIA3Y1+p5AZAmIkki8k9guIjcd7wXi8g0EckXkfzi4uIOCC88/Wf9QSY+/QW1jV5ev+1M7hmXxcebi/jB35bw4bpv/G/8fx/mTwKDJvg/+X//4VNLAC2iE2HSLJgyD+rK/N1DH/4emuqN/8E0TQtIRwwRlWMcU0qpUuC2k71YKTVLRA4C410u18jTCcDrU9htxwoj/Hh9ir9+uI0Zn+5geEYC/7x+JCldIhnVO5ELBibx9iszGLjgdrAV4+l1Lo4fPASpw4y5+MBxkLEMPvwdfPE32PJvmDATMvKMaV/TtIB1RBIoANJbPXcDB06lAaXUYmBxbm7uLacTwMs/G0fsNwfw2QSf3Y6y21EOOzic4HBgc7mwOV04IiNxREbiiorCHhmNPSIacUVgd0Vgc/rPsdkd2OwOxGHHZrMjLc9ttubv2bHZnYjdht3mQBwObDYHdrsDm8OBzWbHZndgtztxOlw47S4cdid2mwObzeYvrIo097ULIvj/3HxcWr7fco74k5vIyZNcZV0Td89byydbirgmN50/XD6YCEfzWP7dS8j68H7urVlDcWwmP66cyraDo3i8xs3Zp/NLP56oBLjsKRg8ERb9HF74AeTdBhf9DlwxRl5J07TTEPDOYs01gXeUUkOanzuAbcBFwH5gJXCtUmrjKbQ5HhifmZl5y/bt2085pjduHkv3TWXYfGD3gd3b/LX5zw7fKTdpWaolGQj+pAHNiULw+BTgvyuyiQDq24fyAjbE4QKbA6+Ceo8PpRQOu40Ip/3bRNM64bQkoZbndjvYbUhzYsRuR+x2/1ebDRzNxx12/0sPFyDVBxFXJKRkI7HdEZcTR0oPXBnpONPTcaWn43S7sUVGdvBvT9M6p1PZWSygOwERmQOcD3QTkQLg90qp50VkOvA+YAdeOJUEYIQrXvji6ANKgbcR1VBFdVUF5WVlHC4vprK8hNqqMuqqKmisqaCpvgpffTWqsRaa6rB56nGqJiLFQ5zTR5TNi83nQXkbUU2NKKXwISif/63VpwQU+BQoBKX49iEOvHYXXpsTn82BVxx47Q582PHabPiw4RMbPqT5AT7Aq0ApH16fF5/Xi8fXRKO3gUZPI/6rghzJ48qfC5oTu1MpXD6IUAqXT+FCEYEdlzhwRHXDFpuCzebAJjbs2EAJBWX1HKhoINppJ6tnAl2jXNiwYZNvH3ZsiAg2bIjPh/iU/+H1gc/X/FUhPh80H6Plz65+qOgkVPEO1O61EN0NFdGNmi++xFdbe9RfmyM5GWdGOi53Os50N66MDJxu/1d7YmK77oY061CNjTQVFeMpKsRTVEz0yBE4urdrG1ytA+k9hk/A61N8tbOU11ft470Nh2jw+MjqEceVI91MHJZKUiTQVAuN1dBYC4010FTj/9r20XAY6iv8a+7UV/hH4hz582Fa3tCPyeaAyASIjPe3VV2ID0WtCFU2G4dtNqqcERyOSWKfimFLvZ2ayBiSUxNpjHBRJVClvFT5GqlqqqaqsYoGbwNNviY8PvMmdtkRRPmwAzZ7BPH1dpIroEe5onuFIrncR7cyL93LvSQc9h712gaXUJbkojzRRXlSBJVJkRxOjMDr8nf9KbsN5bTjs9vA6UDZbeB0gt2Oz2FDHA6Uw44036nYxIZd7AhCrCuWUT1GMabnGGKcusuqPbzVNf4390OHaCoswlN4iKbCQjyHCvEUFtJUWIi3tPSo1zjT0ui94HUcXbuaFHXndSp3ApZMAoF2B3WEyromFn99gNdXFfD1vgocNuGi7GSuGpnO+QO747AHMNDK54OGyqMTw1FfK7/9sysGuqT6H3Gp0KUnxKVSZYvj7vnr+GhzIVeMcPPIxCFEOk++lo9SCo/PQ5Ov6cjD4/NQUVfHU59u4T8bCuiVFMldF/UhPSmCJm/TUed6lRdf852KQvmf+3x4VfNzX/P3j/Hcp3z4Du/Ht/U/+Bqr8fa7AF9Spv+48uHDd6RtafIQXVxFTFE1scXVxBbX0KW4lriSWrqU1OFoOv0+Pq8NPHbBZwevTWiyKbyiUDZwOiKJdEYRFRGD0xHxbZeXzeav0dhtiM0ONlur74l/HSWbIGI7up4jfHtMpPkcf/fdcc9pPnakRgQ0n/Tt92j9ffm2C6/VOf76ku3o67Z67n99888l+H+eVueIzf9nb3XVUW/unsJCfNXV3/m92hMScKSk4OiRgjO5+WtKCo6UFFRDA/vv/gVRI0aQ8dyziNN52n9/2neFfBJoYfadwPFsK6zi9fx9vLVmPyXVjXSLjWDSiDSuGummf0pc0OPZWVzNtJfy2VNay+9+lM1NZ/U2rKvk0y1F3PPGOspqGvn5Rf25/fx+gSW8Y6ktg/k3wp6l/nkJF/7u+JPSjkH5fHiKS/AcPICvsRE8HlRTE8rjQTW1/nMjyuOBI8+b2ny/CeVpwtfYSGltMYVVhyiuKaSmoQqbgmh7FN0jkujmSiTB1QWbwt/FpXworw+8Xv8EPK8X1dIFphTqSH+h70j/4JFjSn173pG+xGOc0/KAb89t+R5tzmk+plofh1bxqKOv2+b5SdntOLp3x5GSjDOlB46UFJw9UnAkN39NScGRnHzSmk7FW29z8L776HrDDfT4f79p99+3dnIhnwSseCdwLE1eH59tLeb1/H18sqUIj08xLD2Bq3LdjD8jlS6RHf/p5pMthfx8zlqcDhszrx3Bmf2SDL9GRW0j9y/cyKKvD3CGO54nrh5GZnKssRfxNsG7v4JVsyHrUpj4DEQYfI3TdLD6IEv3L2Xp/qUsP7icOk8dLpuLUT1HcU7aOZyTdg4ZXTLMDtMwqnVC8Pn8yaQlYXh92KIi/cV/AxQ++ihlL75Ez0ceIeGKSYa0qXWCJNDCqncCx1JS3cDba/YzP38f2wqriXDYGDekB1fnpnNm3yRsBs9b8PkUMz/dwV8+2sagnl2YdWMuaQlRhl6jrX+vO8hv315PbaOX/x2XxU/O6m3sz6UULH/Gvypp8mCYMgcS0k/+uiBq9DaSX5jP5/s/Z2nBUvYc3gNAry69jiSEkT1GEmGPMDfQEKE8Hr655Rbq8lfR6+WXiBpm0ByVMKeTgImUUqwrqOT1VftYtPYAh+s9pCVEce6A7qQnRpGRGE1612jSE6PpGu08rW6b6gYPv5r/Ne9tPMTE4Wk8Omlou/r/jVBUVc99b6zn4y1F5PVJ5K/XDCPV6OSz/SNY8BNwRMLkVyF9tLHtG2jf4X1H7hJWHlpJg7eBKEcUo3uMZmzaWEaljKJvQl9soldtPx5PeTl7rr4GVV9P7wULcKYkmx1SyAv5JBAq3UEnU9/k5YNNhSxYVcCG/ZWU1Ry9mFpshAN31+bEkBhNetcoMpL8ScLdNZoo13ff2PeU1HDLS/nsKqnhvkuymHp2n6APlVRK8fqqAv6weBORTjvP3DCCkb1OY3mJEyneCq9dA4cP+Nceyrna2PY7QL2nnpWHVvqTQsFSCqoLAEiISGBE8ghGpoxkZI+RDOw6EIdN7+fUWv22beyZPIWIzEx6vfwStgh9JxWIkE8CLULxTuBEqhs87CurZV9ZLd+U1VJQXnfkz/vKa6lvM8Kle1wE6V2jSE+MJiMxmi6RTp76ZDs2mzDz2hGMzexm0k/it6Ooip++mM/+ijoeuXwoV48yuOumtgzm3QB7P4dzfgkX/PaUCsZmUkpRUF3AqsJVRx77qvxLasU4YxiWPIzclFxyU3IZnDQYp12Pjjn8wQfs/9nPib/8cno++kc9DyQAOgmEIKUUxdUN7Curo6C8lm9K/Ynhm7Ja9pXVcbCyDp+C7J5dmHXDSNITo80OGYDK2iamz1nN0u0l/Pis3vz2R9nGjh7yNPoLxqtf9BeMJ80K2eUmCmsKWV20mvxD+awqXMXOyp0ARNojyemew8iUkeSm5DK0+1CiHB1b37Gq4qdmUDJzJim/uY/EG280O5yQFfJJoLN0BxmpyevjUGU9PeIjcRo9RDNAHq+PP767hRe+2M3Zmd2Yce1wEqJdxl1AKf8qpx/8P0gZ7N/NLN5tXPsmKasvY03hGvIL/UlhS9kWFAqHzcGQpCH+7qOUkYxIGRE2k9aUz0fBz35G9aefkfHsLGLOOsvskEJSyCeBFuF0J9AZzM/fx2/f2kDPhEieuzHX+DkT2z+EBTf7C8ZT5oC7Xf/GQ0ZVYxVritYc6T7aWLIRj/KQEJHAH8/+I+e4zzE7xKDwVtewd8pkmoqK6fP6fFwZnWf4bbDoJKCZZtXecm59eRX1TV7+ds0wvjcoxdgLFG2BOdfA4YP+ZalzrjK2fQupbaplbfFansh/gm3l27hl6C3cMeyOsCgqN37zDbuvuhpncjK9587BFhMed0JGOZUkYK1+BS3kjezVlcV3jaVPtxhueTmfmZ/uOP42lqcjOQt++gm4R/k3vv/4If9Epk4o2hnNWaln8eoPX2VS/0k8u/5Zpn04jZK6ErND63CujAzS/vIEDTt3cuDee/0zsLUOoZOAZrie8VG8ftuZjM9J5f/e38rP5q6lrtF78he2V0wS3PAWjLgRlv4ZXr/Rv7BeJxXpiOTBsx7k4bEPs754PVctvoqVh1aaHVaHix07luT//TVVH35EyT/+YXY4nZYlk4CIjBeRWZWVlWaHop2mSKedv08exv+OG8g76w5w1TNfcqCizrgLOFww/kn4waP+HcteGAeVBca1b0ETMifw6o9eJdYZy08/+CnPrX8On+rcn5ATb7qJ+AkTKHlqBlUffWR2OJ2SrgloHe7jzYX8fO7ajptYtv1DeP0n4IqGm9+HxD7Gtm8xNU01PPDlA7y35z3OdZ/LI2MfISEyweywOoyvoYG9199A486d9J43l4j+/c0OyfJ0TUCzlIuyU3jrjrOIjbAzedYy5q/cZ+wF+l8MP/0QPPX+5SY8jSd/TQiLccbw+LmP85u83/DlgS+5+p2rWV+83uywOowtIgL3jKeQmGj23Tkdb0WF2SF1KjoJaEHRPyWOt+8cS16fJP73jXU8uHgjHq+BXRnJ2f7RQgfWwEcPGNeuRYkIU7Km8PIlLyMIN753I69uftXYIryFOFNScD/5JJ6DB9n/i1/6lwTXDBG0JCAiMSLyoog8KyLXBeu6mnUkRLuY/ZNR3Dy2D//6Yg8//tdKKmoN/NSePR5G3wrLZsKWd41r18KGdBvC/PHzGZs6lsdWPMav/vsrqhu/u8FLZxA9fDg9Hvg9NV9+SdGfnzA7nE4joCQgIi+ISJGIbGhzfJyIbBWRHSJyb/PhScACpdQtwGWBXFcLXQ67jfvHD+LxK3JYvruUCTO/YHthlXEX+P5D0PMMePt2qDC428mi4iPiefLCJ7l75N18/M3HTP73ZLaWbTU7rA6RcMUVdL3+espmz6bi7bfNDqdTCPROYDYwrvUBEbEDM4FLgEHAFBEZBLiBlv+VBo4X1ELR1aPSmTttDDUNXiY+/SVf7zOon9cRAVf+C3xeeGOqf7OaMGATGzcPuZnnvv8ctU21XPfudby1/S2zw+oQKff8L9F5eRy6//fUre+8tZBgCSgJKKWWAGVtDo8GdiildimlGoG5wASgAH8iCPi6Wucwslcii6aPJT7KyR2vrqay1qA37KR+MP5vsG85fPqIMW2GiNweucwfP59h3Ydx/5f389vPf0udx8ChuRYgTidpf/srju7dKbhzOk2FRWaHFNI64s04jW8/8YP/zT8NeBO4QkT+ASw+3otFZJqI5ItIfnFxcQeEp1lJakIUM68bQVFVPb98/WvjCptDr4QRN8Hnf4Ud4TW+vFtUN565+BluzbmVRTsXcd2717G7crfZYRnK0bUr7pkz8FZVsefKK6n+4guzQwpZHZEEjrUIuFJK1SilfqKUul0p9erxXqyUmgU8CKx2uQxciVKzrGHpCdx3STYfbS7kuaUGvlmNewySB8Gbt/rXGgojdpud6cOn84/v/YPi2mImvzOZj/d+bHZYhorMyqL3nNewdenCvqk/pfCxP+Fr7NzDgztCRySBAqD17iJu4MCpNKCUWqyUmhYfH29oYJp1/WRsb8YN7sGf3tvCqr3lxjTqivbXB5pq4c1b/HWCMDM2bSyvj3+dzIRMfrXkV51uPkFkVhZ9FrxO12uvpWz2bPZcfQ0NO3aYHVZI6YgksBLoLyJ9RMQFTAYWnUoDetmI8CMi/OnKHHomRHLXa6sprzHoE11yFvzwz7BnKSz5P2PaDDE9Ynrw9PeeJjkqmV/+95dU1HeuyVa2qCh63P873P94Gk9REbuvuJKy117rtHMmjBboENE5wFfAQBEpEJGpSikPMB14H9gMzFdKbQw8VK2zi49y8vS1IympbuQX89fi8xn0n3j4dXDGFPjsMdi9xJg2Q0x8RDx/Of8vlNSVcN/n93XKNYfiLriAvosWEp03msI/PETBbbfjKS01OyzL02sHaZbz8ld7+N3CjdwzLovbz+9nTKMN1TDrfGiogts+h9juxrQbYuZtmcfDyx/mruF3MS1nmtnhdAilFOWvvkbR449ji4sj9dE/EnvuuWaHFVQhv3aQ7g4Kb9eP6cWlOT358wdbWbG77Qjk0xQRC1fNhrpyeOvWTrsHwclcPfBqftT3R8xcO5NlB5eZHU6HEBESr7+O3gtex5GUxL5pt3Lo4Ufw1debHZolWTIJ6MJweBMRHp00lIzEaO6as5qS6gZjGu4xBC55DHZ+DF/+3Zg2Q4yIcP+Y++nTpQ/3LLmHwppCs0PqMJEDBtD79fkk3nQj5a+8wp6rrqJ+a+ecSR0ISyYBfSegxUU6mXHtcMprm7h7noH1gZE/gcET/TuSfbPcmDZDTLQzmr+c/xfqPHX8esmvafJ13lnVtogIUu67j/Rnn8VTUcGeq66m7KWX9E5lrVgyCeg7AQ1gcGo8D4wfzNLtJcz81KBhfyIw/u+QkO7ftL7WoO6mENM3oS8PnvUga4rW8LdVfzM7nA4Xe87Z9F24kJizz6bwj4+yb9qtNBXpmcZg0SSgaS2mjE7n8mGp/PWjbXy506C9dSPj/fMHqgth4Z1g4cERHemSPpcwJWsKL216iY/2dv5Z1Y7ERNwzZ9Djgd9Tm5/P7gmXU/XJJ2aHZTpLJgHdHaS1EBEemTiUPt1i+NmctRRVGVTcSxvhX3F067uw/J/GtBmCfpX7K4Z2G8rvvvgdew/vNTucDicidJ08mT5vLMDRowcFd9zJwQcewFfXudZXOhWWTAK6O0hrLSbCwdPXjaS6oYmfz1mL16j6QN5tMPCH8MHvYP9qY9oMMS67iyfOewK7zc4vPvtFp1ts7ngi+vWj97y5JE69mYq589g96Qoav/nG7LBMYckkoGltDewRxx8mDOGrXaX8/ePtxjQq4t+NLDbFvy1lfXjeefaM7clj5zzG9vLt/HH5H80OJ2hsLhcpv/41Gf96AU9xMYWP/cnskExhySSgu4O0Y7k6N50rRrh56pPtLN1u0Aqz0Ylw5Qv+DWgW/Sxs6wNnp53NtJxpvL3jbd7c/qbZ4QRVzJlnkjT1Zqo/+SQs9yewZBLQ3UHa8Tx0+WD6J8fyP3PXUnjYoPpARh5c9DvY9Dbkv2BMmyHo9jNuZ0zPMTyy7BG2lG0xO5yg6nrDjdgTEih+8imzQwk6SyYBTTueaJeDp68bQV2Tl7teW2PcZvVn/Rz6XQTv3QeHwu/TIPiXn/7TuX8iITKBuz+9m8ONh80OKWjssTEk/XQqNUuXUrt6jdnhBJVOAlrIyUyO45GJQ1ixp4y/fLjNmEZtNpj4DER1hdd/DB6DZimHmMTIRJ447wkO1Rzit5//NqxW4ux67bXYk5IofvJJs0MJKp0EtJA0cbibyaPSefqznXy61aBJP7Hd/RPJSnfA5uNuftfpDUsexi9yf8Gn+z7lxY0vmh1O0Niio+k27RZqly2jZvkKs8MJGksmAV0Y1trjgcsGk9Ujjl/MW8uBCoOGNvb/PiT0gtXh8+Z3LNdnX8/FvS7mb6v/Rv6h8FnJN2HyZBzJyRQ/+WTY3AVZMgnowrDWHpFOO09fN4JGj4+75qyhyYj6gM0GI2707ztQujPw9kKUiPCHs/6AO87Nr5f8mpI6g2ZrW5wtIoKk226lbtUqar740uxwgsKSSUDT2qtv91geuyKHVXvL+b/3DVohcvj1IHZY/ZIx7YWoWFcsfzn/L1Q3VnPPknvw+DxmhxQUCVdeiSO1Z9jcDegkoIW88Wekcv2YDGYt2cW6AgO2TozrAQPGwdrXwNt5V9hsjwFdB/DbMb9lxaEVzFw70+xwgsLmctHt9tupX7eO6s8+MzucDhe0JCAifUXkeRFZEKxrauHjnnFZRLvsvLLMoPVvRt4ENUWw9T/GtBfCJmRO4Ir+V/Dc+uf4777/mh1OUCRcfjnO9HSKn3qq098NtCsJiMgLIlIkIhvaHB8nIltFZIeI3HuiNpRSu5RSUwMJVtOOJy7SyWVnpLLo6wNU1hnw6T3ze9AlDVbNDrytTuC+vPvITszmvs/vo6CqwOxwOpw4nXS78w4aNm2m6sMPzQ6nQ7X3TmA2MK71ARGxAzOBS4BBwBQRGSQiQ0XknTaPZEOj1rRjuC6vF/VNPt5esz/wxmx2GH4D7PwEyjv/6ponE2GP4InznwAFDy9/2OxwgiL+0ktx9elDyVMzOvUmNO1KAkqpJUDb3TdGAzuaP+E3AnOBCUqp9UqpS9s89O4NWocb6o5naFo8ry3/xphb+OHX+7+ueSXwtjqB9Lh0bhh0A1/u/5ID1QfMDqfDicNBtzvvpGH7dg7/p/N2CwZSE0gD9rV6XtB87JhEJElE/gkMF5H7TnDeNBHJF5H84mKDFgnTwsa1eRlsLaxi9TflgTeWkO7vFlrzCnjDY2TMyUzInADA2zveNjmS4Ojyw0uI6J9JyYyZKE/n/DcQSBKQYxw77scvpVSpUuo2pVQ/pdSjJzhvFvAgsNrlcgUQnhaOLjsjldgIB68uN2ht+JE3QdUB2NH5d95qj9TYVMb0HMPbO97G6/OaHU6HE5uNbtPvonH3birfecfscDpEIEmgAEhv9dwNGHKPqCeLaacrJsLBhGGp/HvdQSprDSgQDxgHMcm6QNzKpP6TOFhzkOWHlpsdSlDEXfw9IrKzKZn5NKqp8w0ZDiQJrAT6i0gfEXEBk4FFRgSll43QAnFtXgYNHh9vrDZgFIvdCcOvg+3vw+HO3w/eHhdmXEh8RDxvbX/L7FCCQmw2ut91F0379lG5cKHZ4RiuvUNE5wBfAQNFpEBEpiqlPMB04H1gMzBfKbWx40LVtPYZnBrPGekJvLp8rzEF4hE3gvLBmlcDb6sTcNldXNr3Uj7+5mMq6g2YnBcCYi84n8icHEqe/geqsdHscAzV3tFBU5RSPZVSTqWUWyn1fPPxd5VSA5r7+R8xKijdHaQF6rq8DHYW17Bid9tBbachsS/0OQ/WvASdeKjgqZiYOZEmXxP/3v1vs0MJChHx3w0cOEDFG2+YHY6hLLlshO4O0gI1PieVuEgHr60wsEBc8Q3s+tSY9kLcwMSBDEoaxJvb3+z0M2pbxJw9lqgRIyj55zP4GjrPfhOWTAL6TkALVJTLzqThafxn/SHKagy4fc+6FKISdYG4lUmZk9hWvo1NZZvMDiUoRITuP/sZnsJCKubNMzscw1gyCeg7Ac0I1+b1otHr441VBhSIHREw7FrY+i5U67mPAJf0vYQIe0TYFIgBYsbkEZ2XR8msZ/HVGbSHhcksmQT0nYBmhIE94hjZqytzVhg0g3jETeDz+FcX1eji6sL3en2Pd3e9S72n3uxwgqb7z3+Gt6SE8tc6x78DSyYBTTPKtaMz2FVSw1e7SgNvrPsAyDjLv+tYmPSDn8ykzElUNVXx4d7Ovchaa9EjRhBz9tmUPvsc3uoas8MJmCWTgO4O0ozyo5yexEc5jZ1BXLYL9iw1pr0Ql9sjF3esO2yWkWjR/WeIDA1qAAAgAElEQVR34a2ooPyVl80OJWCWTAK6O0gzSqTTzhUj3Hyw8RAl1QaM6Bg0ASLjYVV470HcwiY2JvafyIpDK9h3eN/JX9BJROXkEHv++ZS+8C+8hw+bHU5ALJkENM1I1+al0+RVvJ5vQIHYGQU518DmRVBrwByETuCyfpdhExtv7QifAjH47wZ8hw9TNju0PxBYMgno7iDNSJnJcYzuk8icFd/g8xnQlz/yx+BthK/nBt5WJ9AjpgdnpZ7Fwp0Lw2JRuRaRgwYRd/HFlL34Ip5yA1atNYklk4DuDtKMdl1eBt+U1fLFzpLAG0sZDGm5ukDcyqT+kyiqLeLLA1+aHUpQdbtrOr7aWspe+JfZoZw2SyYBTTPauCE96Brt5DUjC8TFW2BfeKykeTLnu8+na0TXsOsSihwwgC6XXELZK6/gKTVgBJoJdBLQwkKEw86VI918uKmQosMGjGkfPAlcsbpA3Mxpd3Jpv0v5dN+nlNWHV62k2/TpqIYGSp99zuxQTotOAlrYmDI6A49PMT/fgFEsEbEw9ErY+BbUhcdKmiczKXMSHp+HxTsXmx1KUEX07UP8+PGUz5lDU2HozSa3ZBLQhWGtI/TtHstZ/ZKYs2IfXqMKxJ46WP964G11ApldM8nplsNb298Km0XlWnS78w6Ux0PprFlmh3LKLJkEdGFY6yjX5mWwv6KOJdsN2L86dTj0yPF3CYXZm97xTOw/kZ2VO1lfst7sUILKlZFBwqSJVMyfT9PBg2aHc0osmQQ0raN8f1APusW6jC0QF66HA6uNaS/Ejes9jihHFG9uf9PsUIIuado0VFMTh9973+xQTolOAlpYcTlsXDkynU+2FHGo0oAC8dCrwBmtC8TNYl2xXNzrYt7b8x61TbVmhxNUrvR0XL17U7PsK7NDOSVBTQIicrmIPCsiC0Xk+8G8tqa1mDI6Ha9PMW+lAQXiyHgYPBE2vAEN1YG31wlM6j+JmqYaPtj7gdmhBF30mDzqVuaH1Ib07U4CIvKCiBSJyIY2x8eJyFYR2SEi956oDaXU20qpW4AfA9ecVsSaFqBeSTGc078bc1d+g8drwHaRI26Cxmp/ItAYkTyCXl16hdU+Ay1ixpyJr7aWug0bTn6yRZzKncBsYFzrAyJiB2YClwCDgCkiMkhEhorIO20eya1e+tvm12maKa4dncHByno+22pAgTh9NHTP1ruONRMRLs+8nNVFq9lTucfscIIqOm80ALXLlpkcSfu1OwkopZYAbWeBjAZ2KKV2KaUagbnABKXUeqXUpW0eReL3J+A/SildSdNM871BKXSPizBmD2IRf4H4wGo4FF6jYo5nQr8J2MUedjOIHV27EpGVRc2y0JlJHmhNIA1o3bFa0HzseO4CvgdcKSK3HesEEZkmIvkikl9cbMCnNE07BqfdxjW56Xy2tYj9FQZsE5hzDdgjdIG4Wffo7pyTdg6Ldi7C4/OYHU5QxYwZQ92aNfjqQ2O3tUCTgBzj2HEHTCulnlRKjVRK3aaU+udxzpkFPAisdrlcAYanacc3eXQ6CphnxN1AdCIMugzWzYfG8BoVczwT+0+kpK6EpQXhtQFP9Jg8VGMjdWvWmB1KuwSaBAqA9FbP3cCBANvUk8W0oHB3jea8Ad2Zl7/PuAJxQyVsWhh4W53AOe5zSIpMCrsuoejcXLDbQ6ZLKNAksBLoLyJ9RMQFTAYWBRqUXjZCC5ZrR2dQeLiBj7cYsOZL77MhsZ8uEDdz2pxc1u8ylhQsoaTOgCW8Q4Q9NpaooUNDpjh8KkNE5wBfAQNFpEBEpiqlPMB04H1gMzBfKbWxY0LVNONdmJVMjy6RxswgbikQ71sGRVsCb68TuLz/5XiVl0U7A/5sGFKix+RRt2ED3mrrzx05ldFBU5RSPZVSTqWUWyn1fPPxd5VSA5RS/ZRSjxgRlO4O0oLFYbdx9ah0lmwvZl+ZAX35Z1wLNiesfinwtjqBvvF9GZ48POwWlYsZMwa8XmpXrjQ7lJOy5LIRujtIC6bJo9IRYI4RBeLY7pD1Q/h6DjSFxuiQjjYxcyJ7Du9hbfFas0MJmqjhwxGXi9oQqAtYMgnoOwEtmFITorgwK5n5+QU0GVUgriuDLe8E3lYn8IPePyDaER1Wi8rZIiKIGjGCmhCoC1gyCeg7AS3Yrs3LoKS6gQ83FQbeWN8LICFDF4ibRTujGddnHO/veZ+aphqzwwmamDF5NGzdiqfM2jutWTIJ6DsBLdjOG5BMWkKUMQVim81fG9izVO861mxi5kTqPHW8t/s9s0MJmpgxYwCoXbHC5EhOzJJJQNOCzW4TrhmVzuc7SthTYsCn1XT/GjIc/DrwtjqBM7qfQd/4vry5I3y6hCKHDMEWE0PNV9buErJkEtDdQZoZrhmVjt0mzFlpwN1A6nD/14PhUww9ERFhYuZE1hWvY2fFTrPDCQpxOIjOzbX8fAFLJgHdHaSZIaVLJBdlJfN6fgGNngALxNGJEJ8BB3QSaDG+33gc4girJaajzxxD4969lt5y0pJJQNPMcvnwNMpqGtl08HDgjfXM0XcCrSRFJXFe+nks3rWYJm/obLoSiJa6gJWXkLBkEtDdQZpZzkhPAGB9gQEF3dRhULYL6vW/4xaT+k+irL6M/xb81+xQgiJiwADsXbtaukvIkklAdwdpZkmNjyQpxsW6AgPeuHu21AXWBd5WJ3FW6ll0j+oeNovKic1GdF4eNcuXW3bGtCWTgKaZRUQYkhbP+v0GJIHUYf6vukvoCIfNwYTMCXy+/3MKawyYkxECYsbk4Tl0iMY9e8wO5Zh0EtC0NnLc8Wwvqqau0RtYQzHdoItbF4fbmJg5EZ/ysXjXYrNDCYrovDwAapdbsy6gk4CmtTE0LR6vTxlTHE4dpu8E2sjokkF2YjbLD1rzTdFort69cfToYdnisCWTgC4Ma2bKcRtYHO45DEp3QL0BCaUTyU7KZkvZFsv2kxtJRIjJy6N22TKUz4C1qQxmySSgC8OamVK6RNA9LoJ1RtYFDunicGvZidlUNFRwqOaQ2aEERfSYMXgrKmjYts3sUL7DkklA08wkIuSkxbPekBFCZ/i/6rrAUbISswDYXLbZ5EiCI2aMvy5gxVVFdRLQtGMY6o5nR3E1NQ2ewBqKTYa4VL2GUBsDug5AELaUhccObM6ePXH16kWtBdcRCloSEJFsEfmniCwQkduDdV1NOx057niUQheHO0i0M5o+8X3YXBoedwLg7xKqXbkS1WSt2dLtSgIi8oKIFInIhjbHx4nIVhHZISL3nqgNpdRmpdRtwNVA7umHrGkdb0iavx5lzKSxYVCyHRqqAm+rE8lKzAqb7iCAmDPH4KutpX6jtbZhb++dwGxgXOsDImIHZgKXAIOAKSIySESGisg7bR7Jza+5DPgc+Niwn0DTOkByXCQ9ukQat3wECg6tD7ytTiQ7MZvC2kLK6q296YpRokf7lxe3Wl2gXUlAKbUEaPs3NRrYoZTapZRqBOYCE5RS65VSl7Z5FDW3s0gpdRZwnZE/hKZ1hKHueGNGCPVsHiGki8NHyU7KBmBLaXjUBRyJiUQMHGi5+QKB1ATSgH2tnhc0HzsmETlfRJ4UkWeAd09w3jQRyReR/OLi4gDC07TA5KTFs6u4hqr6APtw41IgrqeuC7QRbiOEwL+qaN3q1fgaGswO5YhAkoAc49hxZ34opT5TSv1MKXWrUmrmCc6bBTwIrHa5XAGEp2mBGer21wU27DdiWelh+k6gjfiIeFJjUsNmhBBA9Jg8VGMjdWvWmB3KEYEkgQIgvdVzN3AgsHD89GQxzQqGNheH1+83qC5Qsg0aqgNvqxPJTsoOqzuB6FGjwG63VF0gkCSwEugvIn1ExAVMBhYZEZReNkKzgqTYCNISogwaIXQGujj8XVmJWew9vJeaJgP2dQ4B9thYooYModZCdYH2DhGdA3wFDBSRAhGZqpTyANOB94HNwHyllLXGPmlagHLc8WwwsjisJ40dJTvRXxzeWrbV5EiCJ3rMGOrWr8dbbY27wvaODpqilOqplHIqpdxKqeebj7+rlBqglOqnlHrEqKB0d5BmFUPS4tlTWktlbYDF4S49ITZFF4fbaBkhFE5dQjFnjgGvl9r8fLNDASy6bITuDtKsIqelOHzAoLsBXRw+Sveo7iRGJobVzOGoYcMQl8syXUKWTAL6TkCziqFGzhxOHQYlW6ExPPq/20NEjiwrHS5skZFEDR9umeKwJZOAvhPQrCIh2kVGYrQxI4R6DgPlg0MbTn5uGMlOzGZnxU4avY1mhxI0MWeOoWHLFjzl5WaHYs0koO8ENCsZ6o437k4AdF2gjazELDzKw/aK7WaHEjTfbjm5wuRILJoENM1KctLiKSivo6wmwE+qcT0hJlnXBdoYlDgIILzqAkOGYIuOpmbZV2aHYs0koLuDNCtpmTm8PtChoiJ6WeljSItLI9YZG1Z1AXE6iR41yhLFYUsmAd0dpFlJy7LSxswXOAOKt0BjbeBtdRI2sTEwcWBYDRMF/3yBxj17aDpk7hablkwCmmYlXSKd9O0WwzqjNp5XPijU8ypby07MZlvZNrw+r9mhBI1Vtpy0ZBLQ3UGa1Qwxas9hXRw+puykbOq99ew5vMfsUIImYuBA7AkJpncJWTIJ6O4gzWpy3PEcqKynuCrAJYC7pEF0N10cbiMcl5UWm43ovDxqli9HqeMuwNzhLJkENM1qhhpVF9DF4WPqG9+XCHtEWI0QAn+XkOfgQZr27jUtBp0ENK0dBqfFI2LgnsNFm6GpLvC2OgmHzUH/hP5hNUII/MVhwNTdxnQS0LR2iI1w0K97rHF7CyivLg63kZXk33jezK6RYHP17o0jJYWa5eYVhy2ZBHRhWLOinDSDZg4f2XPYOrtLWUF2YjZVjVXsr95vdihBIyLEjMmjdtlylM9nSgyWTAK6MKxZ0VB3PEVVDRQerg+soXg3RCXqukAbLXsLhF+X0Jl4y8tp2G7OshmWTAKaZkUty0oHPFT0SHFYbzDTWv+u/bGLPaxGCAHE5I0GoNak+QI6CWhaOw3qGY9NYJ1RO40VbYamAO8qOpFIRyR94vuE3QghZ2oqzl4Z1HwVBklARGJEZJWIXBrM62qaEaJcdvonx7HeiJnDqcPA54EiXRxuLTsxvPYWaBEz5kxqV65EeTxBv3Z79xh+QUSKRGRDm+PjRGSriOwQkXvb0dQ9wPzTCVTTrGCoO571+ysDH8FypDis6wKtZSdlU1xXTEldidmhBFXMmDx8NTXUbwz+h4L23gnMBsa1PiAidmAmcAkwCJgiIoNEZKiIvNPmkSwi3wM2AYUGxq9pQZXjjqekupGDlQF24yRkQFRXXRxu48jM4TDrEmrZX8CMLiFHe05SSi0Rkd5tDo8GdiildgGIyFxgglLqUeA73T0icgEQgz9h1InIu0qpUx4T1dTUREFBAfX1ui/VDJGRkbjdbpxOp9mhmKL1dpOpCVGn35CI3nP4GFqSwJayLZzjPsfkaILHkZhIxMCB1CxfRrfbbg3utQN4bRqwr9XzAiDveCcrpf4fgIj8GCg5XgIQkWnANICMjIzvfL+goIC4uDh69+6NiJx28NqpU0pRWlpKQUEBffr0MTscU2T37ILDJqzfX8G4IT0Cayx1GHw5AzwN4IgwJsAQF+eKIz0uPexGCIG/S6h87jx8DQ3YIoL37yGQwvCx3oFP2lGqlJqtlHrnBN+fBTwIrHa5XN/5fn19PUlJSToBmEBESEpKCuu7sEinnQEpcazffzjwxnoOA1+TnjncRlZiVth1BwFE541BNTRQtya4d4eBJIECIL3VczdwILBw/E42WUwnAPPo372/LrC+oMKA4vAZ/q96vsBRshOzKaguoKqxyuxQgip69Ciw24O+hEQgSWAl0F9E+oiIC5gMLDIiKL1sBLz33nsMHDiQzMxMHnvssWOes2TJEkaMGIHD4WDBggVBjjB8DUmLp7y2iYLyABeA69obIhN0cbiN7KTwnDlsj40lcsjgoO8v0N4honOAr4CBIlIgIlOVUh5gOvA+sBmYr5TS97UG8Hq93HnnnfznP/9h06ZNzJkzh02bNn3nvIyMDGbPns21115rQpThK8fIPYd7nqGLw22E6wghgJi8MdStX4+3uiZo12xXElBKTVFK9VRKOZVSbqXU883H31VKDVBK9VNKPWJUUFZfO+iVV15h9OjRDBs2jFtvvRWv10tsbCy//OUvGTFiBBdddBHFxcUAPPnkkwwaNIicnBwmT57crvZXrFhBZmYmffv2xeVyMXnyZBYuXPid83r37k1OTg42m574HUwDe8ThtIsxi8mlDoOiTeBpDLytTqJbVDe6R3UPuzsBgJgzx4DHQ92q/KBdM5DRQR1GRMYD4zMzM0943oOLN7LpgAEFulYGpXbh9+MHH/f7mzdvZt68eXzxxRc4nU7uuOMOXn31VWpqahgxYgRPPPEEf/jDH3jwwQeZMWMGjz32GLt37yYiIoKKCv9M008//ZS77777O21HR0fz5Zdfsn//ftLTvy23uN1uli83dws67VsRDjtZPboYs6x0z2HgbfQngpatJzWyk7LDcoRQ1PDhiNNJzbLlxJ53XlCuackkoJRaDCzOzc29xexY2vr4449ZtWoVo0aNAqCuro7k5GRsNhvXXHMNANdffz2TJk0CICcnh+uuu47LL7+cyy+/HIALLriAtWuP3wVwrIKjLshay1B3PIu/PoBSKrC/m9Z7DuskcERWYhZf7P+Cek89kY5Is8MJGltkJFHDhwd183lLJoH23gmc6BN7R1FKcdNNN/Hoo48edfyhhx466nnLG8O///1vlixZwqJFi3jooYfYuHEjS5cuPeGdgNvtZt++b6dgFBQUkJqa2gE/jXa6ctLieW35N+wtraV3t5jTb6hrH4iM99cFRhoXX6jLTszGq7xsL9/O0O5DzQ4nqGLOHEPx35/EU16Oo2vXDr+eJTuTrVwTuOiii1iwYAFFRUUAlJWVsXfvXnw+35EROq+99hpnn302Pp+Pffv2ccEFF/D4449TUVFBdXX1kTuBto8vv/wSgFGjRrF9+3Z2795NY2Mjc+fO5bLLLjPtZ9a+a6jRxWE9QugoLSOEwrFLKDpvDNhsNGwOzs9uyTsBKxs0aBAPP/ww3//+9/H5fDidTmbOnElMTAwbN25k5MiRxMfHM2/ePLxeL9dffz2Vlf4Fx+6++24SEhJOeg2Hw8GMGTP4wQ9+gNfr5eabb2bwYP9dz/33309ubi6XXXYZK1euZOLEiZSXl7N48WJ+//vfs9GEBajC0YCUOFwOG+v3VzL+jADv0noOg+X/BG8T2MNzOY62UmNS6eLqEpZJICpnKAOWL8MeFxeU64kV9/Ns1R10y/Y2u+1s3ryZ7OxscwI7gdjYWKqrq80OIyis+ncQbBNmfkGU08bcaWcG1tD6BfDGVLh1KfTMMSa4TmDq+1OpbaplzqVzzA4l5IjIKqVUbnvO1d1BmnaactLi2bD/MD5fgB+kUof7v+ouoaNkJ2azrXwbTb4ms0Pp1CyZBEJRuNwFaN8a6o6nusHD7tIAJ/Z07QMRXfSksTaykrJo9DWyu3K32aF0apZMAnrZCC0UGLbnsM2mi8PHMChxEBB+y0cEmyWTgO4O0kJBZvdYIp02Y2YO9zwDDm3wF4c1AHp16UWUIyosl48IJksmAU0LBQ67jcGp8cbMHE4dDt4GKNafelvYbXb6d+0fliOEgkknAU0LwNC0eDYeOIw30OKw3nP4mLITs9lathXfqW9CqLWTTgIW1Z6lpBsaGrjmmmvIzMwkLy+PPXv2AFBaWsoFF1xAbGws06dPD2LU4SfHHU9to5ddxQEODEjsC644XRdoIzsxm+qmagqqCswOpdOyZBII98Jwe5eSfv755+natSs7duzg7rvv5p577gH8+wA/9NBD/PnPfw526GGn9Z7DATlSHNYbzLSWldS8rLTuEuowlkwCVi8MW2Up6YULF3LTTTcBcOWVV/Lxxx+jlCImJoazzz6byMjwWXjLLH27xxLtsge+fAS0Kg57Am+rk+if0B+HOPQIoQ4U2stG/OdeOLTe2DZ7DIVLjt39AtZaSrr1eQ6Hg/j4eEpLS+nWrVugvwWtnew2YUhqPOsKjCgODwNPHZRshZTgL45oRS67i34J/fQIoQ4U2knABFZaSlovOW0NQ93xvLJsLx6vD4c9gJvr1sVhnQSOyErMYun+pYEv260dU9CSgIicDzwEbATmKqU+C7jRE3xi7yhWWkq65Ty3243H46GyspLExEQjfkztFOS442nw+NheVE12zy6n31BSJrhi/cXh4dcZF2CIy07KZuHOhRTXFZMcnWx2OJ1Oe/cYfkFEikRkQ5vj40Rkq4jsEJF7T9KMAqqBSCBkS/1WWkr6sssu48UXXwRgwYIFXHjhhfqTkglaisOGzBzukaOHibaRndi8rLTuEuoQ7b0TmA3MAF5qOSAidmAmcDH+N/WVIrIIsAOPtnn9zcBSpdR/RSQF+AsQkh91rLSU9NSpU7nhhhvIzMwkMTGRuXPnHmmjd+/eHD58mMbGRt5++20++OADBg0a1GG/l3DWOymGuAgH6/dXcvWo9JO/4ERSh0H+v/zFYbvurQUYmDgQQdhctpnz0oOz5WI4afdS0iLSG3hHKTWk+fmZwANKqR80P78PQCnVNgG0bccFvKaUuvJk18zNzVX5+UdvuGzVZYz1UtLhbcqsZdQ2eVl459jAGvp6Hrw1DW7/ClJ00m4x/q3x9Evox98u+JvZoYSEYC0lnQbsa/W8oPnY8YKaJCLPAC/jv6s43nnTRCRfRPJbhllqmtXluOPZfPAwjZ4AZ7Ye2XNYzxdoLSsxS3cHdZBAksCxOp+Pe1uhlHpTKXWrUuqaExWFlVKzgAeB1S6XK4Dwgitc7gK0YxuSFk+jx8e2wqrAGkrKBGeMnjncRlZiFgdqDlDZEJ4TSDtSIEmgAGjdAeoGDgQWjp/VJ4tpWls5Ru05bLP756ro4vBRwnnP4Y4WSBJYCfQXkT7N/fyTgUVGBBXuy0ZooScjMZoukQ5jlpVOHQaH1oHPG3hbnUTLCKEtpXrmsNHaO0R0DvAVMFBECkRkqlLKA0wH3gc2A/OVUnqXcy0siQg57gRjlpXuOQyaaqFk+8nPDRNdI7vSI6aHvhPoAO1KAkqpKUqpnkopp1LKrZR6vvn4u0qpAUqpfkqpR4wKSncHaaFoqDuerYeqaPAE+An+SHFYdwm1lpWYpZNAB7DkAnK6OyiwpaQBHn30UTIzMxk4cCDvv//+keM333wzycnJDBkypKN/hLCTkxZPk1ex9VCAxeFuA8AZresCbWQnZrOncg+1TbVmh9KpWDIJhPudQKBLSW/atIm5c+eyceNG3nvvPe644w68Xv+n0x//+Me89957Qf15wsVQt1HLSjcXh/WdwFGyE7NRKLaVbzM7lE7FkknA6ncCVl9KeuHChUyePJmIiAj69OlDZmYmK1asAODcc8/V6wt1kLSEKBJjXIEvHwH+usDBdeDTO2q10COEOoYl56UrpRYDi3Nzc2850Xl/WvEnw9cZz0rM4p7R9xz3+6GwlPT+/fsZM2bMUa/fv3//af9OtPYREYakxbPOiL0FUofBimegdAd0HxB4e51ASnQKCREJem8Bg1kyCVhZKCwlrZeYNk9OWjz/+O9O6pu8RDrtp99QzzP8Xw+u1UmgmYiQnZitZw4bzJJJQETGA+MzMzNPeN6JPrF3lFBYSrq9r9eMN9Qdj9en2HTwMCMyup5+Q90GgiPKXxzOudq4AENcVlIWL296mSZvE0670+xwOgVL1gSsXBgOhaWkL7vsMubOnUtDQwO7d+9m+/btjB49Oki/ofB2ZOZwoHUBuwN6DNHF4TayE7Px+DzsrNxpdiidhiXvBKwsFJaSHjx4MFdffTWDBg3C4XAwc+ZM7HZ/18SUKVP47LPPKCkpwe128+CDDzJ16tSO+4WFmR5dIukWG2HMzOGew+Druf7isM2Sn9eCrvXeAlmJWSZH0zm0eynpYGrVHXTL9u1Hz5q06jLGeilprcXNs1eyv7yO9+8+N7CG1rwCC++E6fnQrb8xwYU4n/Jx5mtnMiFzAr/J+43Z4VhWsJaS7jBW7g7StJMZmhbP9qIqahs9gTXUes9hDQCb2BiYOFCPEDKQJZNAKAqXuwDt5HLc8fgUbDpwOLCGumeBI1LXBdrITsxmS9kWfErPoTCCTgKaZrCWPYcDrgvYHZAyRG8w00ZWYhZ1njr2Ht5rdiidgk4Cmmaw5C6RpHSJCHxvAfBPGjv4tZ453ErLzGHdJWQMSyYBqy8boWknMzQtgXUFBiwrPfZ/4I6vQE/2O6JffD+cNqdePsIglkwCujCshbocdzy7Smqoqm8KrKGEdIh36yTQitPuJDMhU88cNoglk4DWcUtJH6/dGTNmkJmZiYhQUlLSYT9XuJg8Kp0lv76A2Ag9FacjDEoaxJayLcdcIkU7NToJWFBHLSV9onbHjh3LRx99RK9evYL6s3ZWyV0iSU+M1ms2dZCsxCwqGioorC00O5SQF7QkICI2EXlERJ4SkZuCdd2OEKpLSZ+o3eHDh9O7d29jfkGa1sFaZgtvKv3uhyPt1LTrXlVEXgAuBYqUUkNaHR8H/B2wA88ppY7db+E3AUgDyoCC0464lUN//CMNm40dIRCRnUWP3xx/JmKoLyXdnnY1zeoGdB2ATWxsKdvChRkXmh1OSGtvh+VsYAbwUssBEbEDM4GL8b+prxSRRfgTwqNtXn8zMBD4Sin1jIgsAD4OLHRzhPJS0r5jDDPU3RVaKIp2RtO7S29dHDZAu5KAUmqJiPRuc3g0sEMptQtAROYCE5RSj+K/aziKiBQAjc1PA9yJ2+9En9g7SqgvJa2XmNY6i6zELFYVrjI7jNCnlGrXA+gNbGj1/Er8XWxoN/YAAAfASURBVEAtz28AZpzg9dHA88BTwJ0nOG8akA/kZ2RkqLY2bdr0nWPBtHHjRpWZmakKCwuVUkqVlpaqPXv2KEDNmTNHKaXUQw89pKZPn668Xq/avXu3UkqpxsZGlZycrMrLy096jaamJtWnTx+1a9cu1dDQoHJyctSGDRu+c96MGTPUrbfeqpRSas6cOeqqq65SSim1YcMGlZOTo+rr69WuXbtUnz59lMfjaVe7vXr1UsXFxSeMz+y/A01TSqnZG2arIbOHqNK6UrNDsRwgX7XzvT2Q8WvH6kc47ngtpVQtcNI1i5VSs0TkIDDe5XKNDCC+DhHqS0kfr90nn3ySxx9/nEOHDpGTk8MPf/hDnnvuuQ76LWpa4FqKw1tKt3BW2lkmRxO62r2UdHN30DuquTAsImcCDyilftD8/D4A5e8OMkRubq7Kz88/6phVlzHWS0lrWnBVNlRy9tyz+Z8R/8PUoXpPjNaCtZT0SqC/iPQRERcwGVgUQHtH6GUjNE07mfiIeH7Y54ekxKSYHUpIa+8Q0TnA+UC35gLv75VSz4vIdOB9/COCXlBKbeywSC0uXO4CNM1K/nTun8wOIeS1d3TQlOMcfxd419CI/O0uBhbn5ubeYnTbmqZp2rcsuWzEybqD2lvH0Iynf/ea1rlYMgmoE6wiGhkZSWlpqX4zMoFSitLSUiIjI80ORdM0g1hyicNWG81/53tut5uCgoIja/NowRUZGYnb7TY7DE3TDNLuIaJmONYQUU3TNO3EgjVEVNM0TQtxlkwCep6ApmlacFgyCZyoMKxpmqYZx9I1AREpBvae5su7AaG2T6KOueOFWrygYw6WUIv5RPH2Ukp1b08jlk4CgRCR/PYWRqxCx9zxQi1e0DEHS6jFbFS8luwO0jRN04JDJwFN07Qw1pmTwCyzAzgNOuaOF2rxgo45WEItZkPi7bQ1AU3TNO3kOvOdgKZpmnYSIZ8ERGSciGwVkR0icu8xvh8hIvOav7+8eYc004hIuoh8KiKbRWSjiPz8GOecLyKVIrK2+XG/GbG2iWmPiKxvjuc7a3mI35PNv+d1IjLCjDibYxnY6ne3VkQOi8j/tDnH9N+xiLwgIkUisqHVsUQR+VBEtjd/7Xqc197UfM52EbnJ5Jj/T0S2NP+9vyUix9xD9WT/hoIc8wMisr/V3/8Pj/PaE76/BDHeea1i3SMia4/z2lP/Hbd3M2IrPvBvZrMT6Au44P+3d3YhWlRhHP89oBGV2JpUm3WR0U1dmIuIfUlgbLqEWxGxIRgphJQXXQQJQkR3BnUTUdAHWUhJ33uhpNRFV2vR0qph5CqCi9sKGVp0Uda/i3MmptmZd8f13Zl5eZ8fLHNmzjPsn2eePc+c55x3X8aAWzI2TwKvx/YQsLtmzb1AX2wvAH7K0XwP4as8a/dxStMJYHGL/gFgL+G7p1cBB+rWnIqRnwn7phvlY2A10AccTl17EdgW29uAHTn3LQKOx2NPbPfUqLkfmBfbO/I0l4mhijU/DzxTInZaji9V6c30vwQ81y4fd/pMYCUwLum4pD+BD4DBjM0gsDO2PwLWmJlVqPF/SJqUNBrbvwFHgCV16Wkjg8C7CowAV5pZb92igDXAMUmz/dDhnCHpa+BM5nI6XncCD+Tceh+wX9IZSb8C+4G1cyY0RZ5mSfsknY+nI0Cj/s1sgZ/LUGZ8aTut9Max6xHg/Xb9vk5PAkuAk6nzCaYPqP/ZxEA9C1xViboZiKWp5cCBnO7bzWzMzPaa2a2VCstHwD4z+87MnsjpL/Ms6mCI4j+YpvkY4BpJkxBeGICrc2ya6muATYQZYR4zxVDVbI0lrLcLym5N9PPdwJSkowX9F+zjTk8CeW/02e1OZWwqx8yuAD4GnpZ0LtM9SihfLANeAT6rWl8Od0rqA9YBT5nZ6kx/4/xsZpcA64EPc7qb6OOyNM7XAGa2HTgP7CowmSmGquQ14CbgNmCSUGLJ0kQ/P0rrWcAF+7jTk8AEcEPq/HrgVJGNmc0DFjK7qWHbMLP5hASwS9In2X5J5yT9Htt7gPlmtrhimVlNp+LxNPApYaqcpsyzqJp1wKikqWxHE30cmUrKaPF4Osemcb6Oi9P3AxsUi9NZSsRQZUiakvS3pH+ANwq0NMrPcfx6CNhdZDMbH3d6EvgWuNnMboxvfUPAcMZmGEh2TzwMfFUUpFUQa3pvAUckvVxgc22ybmFmKwnP6ZfqVE7Tc7mZLUjahIXAwxmzYWBj3CW0CjiblDVqpPCtqWk+TpGO18eAz3NsvgD6zawnljH647VaMLO1wLPAekl/FNiUiaHKyKxXPVigpcz4UiX3Aj9KmsjrnLWP53qlu4KV9AHCDptjwPZ47QVCQAJcSigHjAPfAEtr1nsXYUp5EPg+/gwAW4At0WYr8ANhN8IIcEfNmpdGLWNRV+LntGYDXo3P4RCwombNlxEG9YWpa43yMSFBTQJ/Ed46NxPWq74Ejsbjomi7Angzde+mGNPjwOM1ax4n1M6TeE52410H7GkVQzVqfi/G6UHCwN6b1RzPp40vdeiN199J4jdle9E+9k8MO47jdDGdXg5yHMdxLgJPAo7jOF2MJwHHcZwuxpOA4zhOF+NJwHEcp4vxJOA4jtPFeBJwHMfpYjwJOI7jdDH/AiiLJrmgv49KAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#<\n",
"def newton_armijo(f,gradf,hessf,x0,err=1e-4,maxiter=100):\n",
" # ne pas hésiter à copier-coller et adapter gradient_optimal...\n",
" x = x0\n",
" G = []\n",
" k = 0 # nombre d'itérations\n",
" while(True): \n",
" k = k+1\n",
" if k > maxiter:\n",
" print('erreur: nombre maximum d\\'itérations atteint')\n",
" break\n",
" # calcul de d, t, \n",
" g = gradf(x)\n",
" d = -np.linalg.solve(hessf(x),g)\n",
" # on pourrait aussi utiliser np.linalg.inv (mais ça serait moins efficace)\n",
" # d = -np.dot(np.linalg.inv(hessf(x)),gradf(x))\n",
" G.append(np.linalg.norm(g))\n",
" if np.linalg.norm(g) <= err:\n",
" break\n",
" t = pas_armijo(f,gradf,x,d)\n",
" x = x + t*d\n",
" return x,np.array(G)\n",
"\n",
"for eps in [.1,.01,.001,1e-4]:\n",
" x,G = newton_armijo(f,gradf,hessf,.25*np.ones(2))\n",
" plt.semilogy(G,label='eps=%g' % eps)\n",
"plt.legend()\n",
"#>"
]
},
{
"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.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}