{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CS 132 Homework 05 B \n",
    "\n",
    "### Due Thursday August 12th at Midnight (1 minute after 11:59pm) in Gradescope (with grace period of 6 hours);\n",
    "### No late HWs accepted!\n",
    "\n",
    "\n",
    "Please read through the entire notebook, reading the expository material and then do the problems, following the instuctions to try various features; among the exposition of various features of Python there\n",
    "are three problems which will be graded. Each problem is worth 10 points. \n",
    "\n",
    "You will need to complete this notebook and then convert to a PDF file in order to submit to Gradescope. Instructions are given here:\n",
    "\n",
    "https://www.cs.bu.edu/fac/snyder/cs132/HWSubmissionInstructions.html\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are some imports which will be used in the code in the rest of the lab  \n",
    "\n",
    "# Imports used for the code in CS 237\n",
    "\n",
    "import numpy as np                # arrays and functions which operate on arrays\n",
    "\n",
    "import matplotlib.pyplot as plt   # normal plotting\n",
    "\n",
    "\n",
    "# NOTE: You may not use any other libraries than those listed here without explicit permission."
   ]
  },
  {
   "attachments": {
    "Screen%20Shot%202021-08-07%20at%2010.27.41%20AM.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAFVCAYAAADCAXWmAAAK2GlDQ1BJQ0MgUHJvZmlsZQAASImVlwdUU9kWhs+9N52EkkAEpITekU4AKaGHIr2KSkgCCSWGhAAidgZHcCyoiGAZ0VERBUeHImNBLNgGxd4nyKCgjoMFGypzA48wM2+999bba52cb+3ss8/e55671n8BoIRyxOJcWB2APFGBJDbEn5GcksrADwAEqAEqIAHA4UrFrOjoCIDa5Px3e3cLQIr5up0i17///1+NyuNLuQBAaShn8KTcPJQ70fGMK5YUAIAcQP0mRQViBV9FWVOCFojybwrOmuAPCs4YZwx5PCY+NgBlBgAEMocjyQKAbIv6GYXcLDQPWdGDg4gnFKFcirIPV8DhoXwMZdu8vPkKHkTZEo0XA0BBTwcwM/6SM+tv+TOU+TmcLCVP9DVuhEChVJzLWfB/Hs3/trxc2eQe5uggCyShsYr90PO7kzM/XMmijFlRkyzkTdSkYIEsNGGSudKA1EnmcQLDlWtzZ0VMcqYwmK3MU8COn2S+NChukiXzY5V7ZUoCWJPMkYzvi94oSC7LSVD6BXy2Mn+JID5pkguFibMmWZoTFz4VE6D0S2Sxyvr5ohD/qX2Dlb3nSf/Sr5CtXFsgiA9V9s6Zqp8vYk3llCYra+PxA4OmYhKU8eICf+Ve4txoZTw/N0TplxbGKdcWoJdzam208gyzOWHRkwyEIBJwAJehNkkAFPCLCxSNBMwXL5AIswQFDBb6tvEZbBHX3pbh5ODkCIDi3Z24Dm/o4+8kRL845VtZhF71FhQYU75QAwCOBKKPZXjKZ+EHgCoegPP6XJmkcMKHUfxg0aenBjSBDjAAJsAS2AEn4Aa8gB8IAmEgCsSDFDAXrVUA8oAEFIFSsAyUg0qwDmwCtWAH2AX2gYPgMGgDx8ApcA5cAlfBTXAfyMEAeA6GwTswCkEQHqJANEgHMoTMIBvICWJCPlAQFAHFQilQOpQFiSAZVAqtgCqhKqgW2gk1QD9CR6FT0AWoF7oL9UFD0GvoE4zAZFgT1ofN4RkwE2bB4XA8PAfOgvPhErgMXgPXwPXwAbgVPgVfgm/Ccvg5PIIARAWhI0aIHcJEApAoJBXJRCTIYqQCqUbqkSakA+lGriNy5AXyEYPD0DAMjB3GCxOKScBwMfmYxZjVmFrMPkwr5gzmOqYPM4z5iqVg9bA2WE8sG5uMzcIWYcux1dg92BbsWexN7AD2HQ6Ho+MscO64UFwKLhu3ELcatw3XjOvE9eL6cSN4PF4Hb4P3xkfhOfgCfDl+C/4A/iT+Gn4A/4GgQjAkOBGCCakEEWE5oZqwn3CCcI3wlDBKVCeaET2JUUQecQFxLXE3sYN4hThAHCVpkCxI3qR4UjZpGamG1EQ6S3pAeqOiomKs4qESoyJUWapSo3JI5bxKn8pHMpVsTQ4gp5Fl5DXkveRO8l3yGwqFYk7xo6RSCihrKA2U05RHlA+qNFV7VbYqT3WJap1qq+o11ZdqRDUzNZbaXLUStWq1I2pX1F6oE9XN1QPUOeqL1evUj6rfVh/RoGk4akRp5Gms1tivcUFjkIqnmlODqDxqGXUX9TS1n4bQTGgBNC5tBW037SxtQBOnaaHJ1szWrNQ8qNmjOaxF1XLRStQq1qrTOq4lpyN0czqbnktfSz9Mv0X/NE1/Gmsaf9qqaU3Trk17rz1d20+br12h3ax9U/uTDkMnSCdHZ71Om85DXYyutW6MbpHudt2zui+ma073ms6dXjH98PR7erCetV6s3kK9XXqX9Ub0DfRD9MX6W/RP678woBv4GWQbbDQ4YTBkSDP0MRQabjQ8afiMocVgMXIZNYwzjGEjPaNQI5nRTqMeo1FjC+ME4+XGzcYPTUgmTJNMk40mXSbDpoamkaalpo2m98yIZkwzgdlms26z9+YW5knmK83bzActtC3YFiUWjRYPLCmWvpb5lvWWN6xwVkyrHKttVletYWtXa4F1nfUVG9jGzUZos82m1xZr62Ersq23vW1HtmPZFdo12vXZ0+0j7Jfbt9m/nGE6I3XG+hndM746uDrkOux2uO9IdQxzXO7Y4fjaydqJ61TndMOZ4hzsvMS53fmVi40L32W7yx1Xmmuk60rXLtcvbu5uErcmtyF3U/d0963ut5mazGjmauZ5D6yHv8cSj2MeHz3dPAs8D3v+4WXnleO132twpsVM/szdM/u9jb053ju95T4Mn3Sf733kvka+HN9638d+Jn48vz1+T1lWrGzWAdZLfwd/iX+L//sAz4BFAZ2BSGBIYEVgTxA1KCGoNuhRsHFwVnBj8HCIa8jCkM5QbGh46PrQ22x9NpfdwB4Ocw9bFHYmnBweF14b/jjCOkIS0REJR4ZFboh8MMtslmhWWxSIYkdtiHoYbRGdH/1zDC4mOqYu5kmsY2xpbHccLW5e3P64d/H+8Wvj7ydYJsgSuhLVEtMSGxLfJwUmVSXJk2ckL0q+lKKbIkxpT8WnJqbuSR2ZHTR70+yBNNe08rRbcyzmFM+5MFd3bu7c4/PU5nHmHUnHpiel70//zIni1HNGMtgZWzOGuQHczdznPD/eRt4Q35tfxX+a6Z1ZlTmY5Z21IWtI4CuoFrwQBghrha+yQ7N3ZL/PicrZmzOWm5TbnEfIS887KqKKckRn5hvML57fK7YRl4vl+Z75m/KHJeGSPVJIOkfaXqCJiqTLMkvZN7K+Qp/CusIPRYlFR4o1ikXFlxdYL1i14GlJcMkPCzELuQu7So1Kl5X2LWIt2rkYWpyxuGuJyZKyJQNLQ5buW0ZalrPsl+UOy6uWv12RtKKjTL9saVn/NyHfNJarlkvKb6/0WrnjW8y3wm97Vjmv2rLqawWv4mKlQ2V15efV3NUXv3P8rua7sTWZa3rWuq3dvg63TrTu1nrf9fuqNKpKqvo3RG5o3cjYWLHx7aZ5my5Uu1Tv2EzaLNssr4moad9iumXdls+1gtqbdf51zVv1tq7a+n4bb9u17X7bm3bo76jc8el74fd3dobsbK03r6/ehdtVuOvJ7sTd3T8wf2jYo7uncs+XvaK98n2x+840uDc07Nfbv7YRbpQ1Dh1IO3D1YODB9ia7pp3N9ObKQ+CQ7NCzH9N/vHU4/HDXEeaRpp/MftraQmupaIVaF7QOtwna5O0p7b1Hw452dXh1tPxs//PeY0bH6o5rHV97gnSi7MTYyZKTI53izhensk71d83run86+fSNMzFnes6Gnz1/Lvjc6W5W98nz3uePXfC8cPQi82LbJbdLrZddL7f84vpLS49bT+sV9yvtVz2udvTO7D1xzffaqeuB18/dYN+4dHPWzd5bCbfu3E67Lb/DuzN4N/fuq3uF90bvL32AfVDxUP1h9SO9R/W/Wv3aLHeTH+8L7Lv8OO7x/X5u//PfpL99Hih7QnlS/dTwacOg0+CxoeChq89mPxt4Ln4++qL8d43ft760fPnTH35/XB5OHh54JXk19nr1G503e9+6vO0aiR559C7v3ej7ig86H/Z9ZH7s/pT06elo0Wf855ovVl86voZ/fTCWNzYm5kg441IAQQecmQnA672oNk4BgIbqctLsCW09btDE98A4gf/EE/p73NwAaEI1RwyKLHQ+opCz6ExBh0ISxfsB2NlZOf5l0kxnp4lcZFRZYj+Mjb3RBwDfAcAXydjY6LaxsS+70WLvAtCZP6HpFYZDtXyToX7/UPSN4l/BP21C7/+lx3/OQFGBC/jn/CfvGBVP61KWlwAAAJZlWElmTU0AKgAAAAgABQESAAMAAAABAAEAAAEaAAUAAAABAAAASgEbAAUAAAABAAAAUgEoAAMAAAABAAIAAIdpAAQAAAABAAAAWgAAAAAAAACQAAAAAQAAAJAAAAABAAOShgAHAAAAEgAAAISgAgAEAAAAAQAAAWigAwAEAAAAAQAAAVUAAAAAQVNDSUkAAABTY3JlZW5zaG90UbJbhQAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAAnNpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDYuMC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iCiAgICAgICAgICAgIHhtbG5zOnRpZmY9Imh0dHA6Ly9ucy5hZG9iZS5jb20vdGlmZi8xLjAvIj4KICAgICAgICAgPGV4aWY6VXNlckNvbW1lbnQ+U2NyZWVuc2hvdDwvZXhpZjpVc2VyQ29tbWVudD4KICAgICAgICAgPGV4aWY6UGl4ZWxZRGltZW5zaW9uPjc1ODwvZXhpZjpQaXhlbFlEaW1lbnNpb24+CiAgICAgICAgIDxleGlmOlBpeGVsWERpbWVuc2lvbj44MDA8L2V4aWY6UGl4ZWxYRGltZW5zaW9uPgogICAgICAgICA8dGlmZjpPcmllbnRhdGlvbj4xPC90aWZmOk9yaWVudGF0aW9uPgogICAgICAgICA8dGlmZjpSZXNvbHV0aW9uVW5pdD4yPC90aWZmOlJlc29sdXRpb25Vbml0PgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4KbnrttwAAQABJREFUeAHtvQmYHdV1LrrO2JOGbs1oHgBNaACEDAYEOMZmCKMhdpwEHOf6xu8lvr7PcfLZL89T4uH6y712romv7Vw7cTwk19gYg42MwRgwo0AICQlhhGa15nnq6Uxv/atqV+8+Oqf7nO4zVNVZWzpdVbv2sNa/d/21atWuvSM5DqRBEVAEFAFFwHcIRH0nkQqkCCgCioAiIAgoQWtHUAQUAUXApwgoQfu0YeolFjxemUyGyvF8mTzZbHZQsXEeaU36cuoYtOASTla7Tug2lP4liFmzJMAjSPLWDBifVRThhlIftM8apVxx0IT5zRiJRAi/cgLKKDfPYOXbBDAceQYru5xzxfQyuBnZIK/ZL6f8IKXt6+ujN998ky688EJqamqSflPJNg8SFkGQVS3oILTSEDLiAotGowN+iAPhwBrGD/v4pdNpOUaRIChzHlvkwQW8Z88e6u7ulvRIYwc7D8pD6OrqoldffZV27twpx+aPLRPK3r59O504cYJSqRR1dnZST0/PgDqMvKZcvu3IeSOjkcUcm3Soz5YL++YYdaFuyLhlyxbq7e0V8XAe8ZARWwSzj3JN2fn7pm5TPo4Rjh49SocOHZJ8yIPzCPnp7GNTtmkTkwf58sswx6Z+HBcL+XWYY6QHDhs2bKCTJ0/S/v0HvCLsur1I3ak7AkrQdW+C4QtgLipsN2/eTD/4wQ/ou9/9Lj322GN05swZIZxYLEb4GbKMx+NyjFpBTOY8tggg0K9+9at09uxZj7BMPThv50GZCI8//jj94he/EOI1xIHtb3/7W/qXf/kX+s53vkNPPfUUnT59Wso8cuQI/fM//7PcBAwpog7sG1lxzM8AXhziDZEamU39ktbSBemMnIlEQmQEsUEnIx/O79q5i37yk5/QqVOnJM2LL75IO3bskDpRtpHJ1GPkM7IYOZAZ+j3xxK+9vIjLlwtxRi7kRXn4mTbBOeRBQLw5xtYcmzpxXCjk12nnQ3qUtXfvXnruuefo8OFDclyoHI3zBwKFW9kfsqkUQyBgLmZYr5/+9KeprW0ULV26lC2j/UK0sFBh1cJyPHbsmMTDesIFigAShyX71ltv9Vu/zA+wnpEHFi8sLlzUhtiwhYWNGwLKP3z4sNwYZs2aRdOmTfMueKT76U9/KoS78rKVNHPmTJoyZQqNGjVKrHjk/d3vfif1wqpFHbBAX3/9dQKB4xiW5cGDByUN0kIWlAsS3bhxo6Q3hARLGfJu2rRJrEM7HazF5uZmmj59OhnChv6dezvp7rvvpueff16I8ec//zlt27ZNbjTQCzLghgXsQPCI28l4wkUAsjfYQQaQbGfnXtq69S3BGXnxg8zQCfmQDsf79u2jN954Q6xuWN6vvfaayI6nF+RBOmCBGxqOcQNBOuCEdtm0cRMdP34cKgwIBgvovnv3bsEIsqNc4Lhr1y550nnooYeopaVF+grymHwDCtMDXyAQ+ywHX0iiQpSFAC4qWEcgpn/7t3+jG2+8kW699VY677zzaOHChTRu3Dj61a9+Rddccw0tXryYxowZI4+2IIYHH3yQ5s+fL0R33XXXCnH++7//u8S1d7TTj370IyE5EBdIacmSJZ7VjThYvyAdED38mD/72c+kfJTZ3t4ueoDQXnrpJal/2fJlcv4b3/iGEPSkSZPo29/+thAPrGwQMSzDr33ta3KzeOaZZ6ROkN5HPvIRicNNCLIgL6zVrVu3itV+wQUXUEdHBz3wwAN0//33i9yQCYSOpwncrCZPnszOkhx97X9+TUgJNwmEAwcO0Jo1a/jG1iY+WRA8fLMgQzwRXH755bR69WohtQULFtAXv/hFeSpA3ZAB+R9++GEhfbQH2iGbzYksaIPW1lb613/9V5EVuI0eM5oOHzpM73//+ymZTIoMD/70QakPhGtuYCBltAeIfNGiRQRCNTcuyIWbG240U6dOlTJQNwLygZwhE9oY7YN2B8H/1V/9lcgLmd/5zndKGbhhgaiVoAU+X/6J+1IqFWpIBMxFBesQltE999wjeX784x/L4+u9995Lfb199NGPfpQ++MEPCpmDiJAPlhke50Hi73nPXfSXf/mX9Mgjj4glefvtt8vF/yd/8idywX/+85+nW265hcaPHy/kh3QoDwSEe/v5559Pv//7vy9EDCsaZAtiNWTxla98RcgWZAcSBnEb2SEzyPJLX/qSEN+VV17J8ryH/vEf/1F0eNe73iWEA0ID8X/mM5+RekAw0BvWLm44IKvf/OY39Ld/+7c0Y8YMIUS4e0BKc+bMEVz2dO4R69UGFhYpZEDZzz77rDw5QHbIZwL2c0y60Ae4QXcQI2T45S9/SVe8/Qr66QM/lScE3CSBN9wmuDnBegUxQy64FFY/spqWLVtG73jHOwRz3HC+9a1v0cc//nGaN2+etAcIFjde3FSffPJJuSnhaeXmm28W6xntgLrNjRByQjYT8DT0myeeoE988pNys8Y5YIN8uNnhGDcHPAGgPRDs/KYc3foDAXVx+KMdhi0FCA4WIqw6BDyyL1++XFwQiWRCCBQXPKypH/7wh/T000/LRQ/3Bh59YUWBlGCZghxAsCBaEADcArDAEIcAyw2P6EgLK3Xs2LHijsA55M0PIAAQ2l/8xV+ITCDnWDwmaUFysHwnTJgghASyhR4IIBPIh3pBsEgHWZAW7oX/+I//kBsM3ABwx8D6BA6QCwGEA53NMeKikahjtfZzmTx9QM9rr71Wni5gLZuRDciDIDeTqGOZwsKdOHGipAHJApuW5hYpF7qhPtwskA5YwXKFi2T9+vUiz3XXXSdY4+YGkoTb5xOf+ITUfd9994nrCHIj4GkEFvhTbKmjTFjwN9xwg+T5p3/6J2lHlk7S2jeUdevW0eVXXC43EUO8cBnh5oEbMrDEjRqyQn7RzyJ4KVD/+AYBJWjfNEV5guBCBimCuN797nfT17/+dbG4YE2DuEAkIAlYmgjwG+Ox/7LLLpOLFPEgUFhXsPZgGcOHDLLeyf5SEA4IEo/EJoA4YdnBVwuix4UP9wcsQUPiJi0ufJATCBCEgB8sylRfSsjqm9/8ppAPHtlBvH/8x3/MLxsfo5dfflnkQT2Qz8gCeZAfvmfUt3LlSnk8hzsCJI7y8VgPHWEdwpKG6wXkCB2AFfJnM/03EsgMNwLyr1ixgl544QUhfMhs5ICVDEsb+qAu3NQgC/DEFi4m3CCQBgSPJxNgOXv2bCFVECsI+YorriC4YxKJpJQD8kR5aEc8OcBKhk8bAXIBa8j1B3/wB3Jzww0IsuIdA4h97dq1nM4ZQWITNG58Gza8Jm4ZlImbBG5qKA8BOCC9+RkSl5P6x3cIqA/ad01SnkC4wHBRghBAXiAoXNjwPYNsQdS4oI0VivOwVEHUIAKQnbH6QPTYB9nBgkN+XNjnX3A+NSWbhExQFogE/trbbrtNSAgyoE4QG4K56JEfRIl4EAOOYZ2DbFA/rPHjJ47TnXfeKccgv3WvrqNVV6+it7/97UIikAXEhhvR6NGjhcwgoxnLu5RvELO5brgIQEh4kkB6uFSgK1wgwAbWI/IjHaxkBBAscJk7d65gAtlwY4AuSAN3AXC49NJLxeKEDDgHvy1wQT3QCZYo9tEOuEHCj30tW+W48QBruDnwhIC2GD9+nOSFTHhKgIUMXW666SaRGeUhAEPUh7JWrVolTyvQB24PyAzsoROIFiSPLfLAeke74kYBuXDTBV7QEe2ANPZPKtM/vkVAP1TxbdNUXzBYr3j59qlPfaqkygwJlJRYE5WMgOJaMlQNlzAQLg50YNfdVrSBSk5TtATnhJQTkDRGVlin+Nkh/9iktdPAqrr66qu9KDtPofSwvPLTmHRmi8IK7rvv3cw5szXpUa6Jw9bsy3k0fl5+xCOYtLZciM8/NuWZLdKYYMqwj80+tvl5MCLExJt9ibD+mDLz5bCSeLvAFeXY9Zj8SGTvm0yF4kxak8Y+tsuWeOhQBFOTX/JomgHtYrCp1VYt6FohXYN67IsQF71zDTqPvqjePo9jSWPF28eF9pEHwZQzWBqTzqTJPzZlIN6kKRZnn8e+na5Qfjt9/nlzDvEIKMvE5e87KQbKVywt4m25TDpTB7Ymbqh6TL0mr53P5DVxdlqzL3KAWLkLIJ3JY87rNjgIKEEHp61UUkVAEWgwBOo6DhqPf3jBgbfZGJiPfbwkEQug3IaAGwRWYyME0RWKjkxfg/Ng1pgv4GyktgXgjagv9A7I9YvRRXjBjZfd733ve2VruAxqVDL4gqAxFAhkgbf5IGkoWxppgJT5DTa/jc/1dFOEx6RGGDzp4JVEyWdlZbv58+t4giJJ/hoNvueAdOxhwci6ZXm0A3QVfUFeIQ8ZHiYIXaNo3xzeLYzsRux3uILUn8FTGBWDj5YwlYEhamPsVBrruhK0TcIYooTPlIcd0ikiJq2GCI2kKxqUh5Xxlx0N0bSiJA83JP7IKOzE7DUojyVn1vMOg7CDoY72h1DVkrmuBG3fdfABg3xSy19tYXwq7kzYmpccxQCQ8xgBwJ1aLCy2wCtpQcOix1hcbCGvfVMpJlO147P8CXcEnx+DtMJsQcNaxhMV2pZBjchFDAs6hBYldHWfhDLQl4+jPBY7V/LTZLV7XfXKl/bFtRWA/gwOABdgzD4+Uqp2qCtB29cZiC/C5IyAxwd8mYWPBYYapuSRPHzXyAzCqmAw8yKbQf4VLHrYRQEnwQslcGcJbXAJS3RF/5DjEJIzGtDVFbtR6MrtKjelMLcvlOUg/Rl9GgcB0df0SYhczVBfgrY0s8eT4hNlPEKYL76sZDXfxV0Snwj7NliWl29lVMHKQkAMDTY3QnorOhcLKOwo7Tz9WjercxM3VoxvCNqGXSwlt3dm+SWJe2+1k/TvozGZpORRsMS7r2OJ9Rdhuy48i5xPm7skHmk0KAKKgCJQawR8SdACgntHBTnnE2oxkITTh7j7goCN2wTlYt8Mm8EkO5gvAcEQNbZmv1i9Gq8IKAKKQDUQaCjTEEQLUoZFjB/2Qc4ImN0NE51jkhn4nc1NwWyrAb6WqQgoAorAYAj414IeTOphnDPkjJePWDYIljOGyeCFJMgas41hikkzjeUwqtAsioAioAhUFAHfE3RJFixbwvKBShFoDDlj2B6WD9rJ8x1jCkksI4S5e7GqBj6WwRSXqM+kL1KcRisCioAiUBME/EnQ7IpgppSXfymebD3X2+MMv0H8OQGeZ47H8Dr3ZV4kFqc4z8PrjJt1/MkgXszvC+v5r//6r2WECFwaWBoKSyth3DWIGSSuBH0OyBqhCIQGAdCITSWgGvz8GPxJ0EAKiDHpZni1itwpXhWk6CfcDkF7ozgY+QiPn47xhOpC0FZL4CUgyNesiAyyxnp7CJiEvSRrXVLrH0VAEQgiAqADPxNyPqb+JGiXnPH5ZzOvajyigLJcksaqF1g1A4t6Yh/zIZulgEDcJihRGyR0qwiEBwFDznjY3tOZZbcmr385jnglHGesBBjAb4a0Pwna6hODf+rqQCrk6ro4PHIFMbvBjGPGEkB/+Id/KKM0sGhnPMGuEHaHIJg0bhbdKAKKQMgQMLbab59J0//1iRT97sUc3XBHhL78uSQtXRLzHTkDft8T9ODPIw4J42+O0Rdytog5v3+ByDETFX4aFAFFoHEQMNbz6dNZ+uo3ckzOZ+idqw7Qow/OpBXL47RwQZR5AQME/OWPDtc46EHIGV0RBI7hdXgRiC2+UtSgCCgCjYNAOhOhHftyNLt1L03qOMGKx+nYcQwO8CcG4SLoEjCGKwMfp2Ab5bmkNSgCikCjIJDjyfUj9H/fwxOydR2mf39oPCveTDf8Xo5XUO93ifoJDf+7OPyElsqiCCgCgUQAD9dwX2Ca7T+8NUUXdsVoV6SDZi2P0ZUrHXL2m3sDQAeSoAGkHUq998EH7b1EtAvQfUVAEQg9Aubab+s9Qle9vY3escIsEOIwyBAe0rrgEziCLnSX6ydsZ1RHISThc4ZbI51O09q1a+VF4cUXX6yjNwqBpXGKQIgRSJ86Q9TUIho6o8TY1VmqlVdjXALlhDXkzIsZUCePY9y9J8tfAFpvXXNFUJY3s845TIb0hS98QVZDwB3VHv9cY+y1OkVAEagVAq4Vl+W55jNHj1Ccl9hDEKu6CG3USrTB6gkMQRty7unJ0Y9+0kcz3tFDsy7ooW9+u4+/DGTXBU9Kl+03pQfoLHNKMxmf4cVHf/WrX9HSpUvpwgsvlMZRkh4AlR4oAqFGIMcEjakjYmPbXT3z/KU+0z5QBA3stu/I0j1/lCbCCJkxOfqrj6bp1fUYI4Nx0PwiAImKBKzQMmPGDPlYZc+ePbRv3z5valG1pIuAptGKQBgQcB3M6UMHKTZxEkV4iTwn+Nh8ZgEDR9A8nQaHHF23gujaSx1wT52yaNnadRqg/y/mecasdVgE9oUXXqCXXnqJsLyWBkVAEWgMBDIHD8pEatA2l2XDziVuv2pvbiN+lc+Ty+A4bWqErrs5Sk8+4n5kcmGU5s1z7jPwcLgT2nn57B1MljR9+nSZbnT79u00c+ZMWfsQacwbXju97isCikAIEHD9oxme+z3HH6dFR40SpfiZ2/fKBYagDfFiYpP/9T+S9NiNGX7Rl6NVV8Vp8SImaB6lUYxkzTwbEyZMoLvuuksm53//+98vjQMLurm5WV4WFsvv+1ZUARUBRWBIBDLs3ojwHDyx0WOcQdGGVIbMWb8EgSFoG6IF82OEX3/gdQPxtMJLtw8WzBwc8+fP95IZ37OSsweJ7igC4ULAffzOnj1D0TbHeg7KNxEB8kE7zmU8rRhShS8asWbwhhdvItxuZuKxtfdxGsRs4tzksjFxZmuf031FQBEICAIuF2T5/VPm2DFKzJotggfFIPOtBW0APHcLfI2l7HqRcIc0PzltzktbeK4PU5aTpD+NHe/kcIgb+8Y9YuJ1qwgoAsFDINfTTZgRKcZTDktwrWq/a+JLgobVmkqlCHM2yxhmj5CLw4k8EfduWTxVeWdA3H38VYxa0eXhpqkVAb8hkNrTSbEpzscp8sitBD38JhrFb1mx2jaGxSHAkTHYG1chUH5JiCEchazh4UoCvs9k0jwDluO3Gm45mk8RUATqi0D6wD5qvviS+goxjNp9aUGDEGE9Y97mIQPuhJjfubeXojwaQ8bZVdCShosDH7hoUAQUgYAhAB5gfsicPMHrlDZTtLU1YAr4eDY7EHRZAYvKgqA1KAKKgCLACGDMc4TngEjv30/R9nZnBIdL2kEByJcWdFk+X2NBsysiCou76Orfw2+SSrpNhi+F5lQEFIFyEIi4C3JkeQ6eWEdHOVl9k7auBF3Mr1w2ITJJI4+XD6StQRFQBBoXAfedFL4ezJ46SU0LFzpYBIwb/DEOeoR8irHQGkKOwAj7SMjRUfWKIJDrwuQ9EfY/txVJ4e/oilrQMtSN71AHeUKS3/zmNzIS49prr5X5LlY/spraO9rp+uuvp1bXWW9cGblsHsWW6icy6bC1izDx/sZ++NLZug6/FP/ntNsx7DrbuqJljL758f5vtfpLCMzcz7hTPGNlnGewxMvCIIaKWtCGcDGN54YNG+i6664jzH/x4IMP0qjRo2R6z2eeeUZwMmSOA4yS8NwTiCgVTDed+Jqwb/KZLcoKY8An7UZHsw2jnkY3bIf4jD/w6htdoQjve9eDHR94JYsoAB2NnmZbJGlJ0VYZ2RMnKM7Ti1YsmLKxNfsVK/zcgipqQZviQbiY1nM/vz3FCtpHjhyhe++9V2aRe+KJJ0wynt8oKx+BHGCLO4sXfPAbsd5iUfN2yMA3Ssy/keO85hfhO6fcKErJP2QFfkzAiqVZX+4cEcYWS/aU8B2PHxUZWiavfXn+73SMr4cQty0sZrngeQf9mfsx5iwOdftKD+D+bPStVH/mJ3Jgl+LJkSjO/SbBK8VyAEeM+FqBdR7jduFBCcJXUnL1/lSUoM1df9asWfTxv/44bdq4iR577DH5KhCkHfcmyUZfjMj6gOvWraNfrF5N82fPlpUOhFwBQokBRAXgsZQNwAdBy5dCJeYPXDLWF7pGMnG5wdWik9QTI7Qv9I3yjSiLi6KMvlFPuYdVN+sK/bL8WXKUdcVXtGFvX9yUMkZfzM8Mg2OEATc1jHlO7dwhs9dRLEpZ+KIr0HdQRCzBBM0f0TmEz21WxVAVgoa8M2fMpO6ubsIagJiHef369TLN5/jx40UdEHEinqArr7ySxvD38Rs3bqQIj2Nmu5BxxJeDJQZuYAAF6yrS0tIQBC3+Se4k0WSTWFglY1UipH5Jhts0buRywfG4+EiCx8ZX4CLzi37nyOEStPR/1jXKRg10t9vXmC52XH45gUojOnM7c3+OuP0Z+hTSrzS9cuwNc400tqTjEyfyB2y8QKxF/KWV46CaL4cYkGzpRzLcLjA4q9wfK0rQpqPAtQHL+fSp03T33XfLV4E///nPhYjf8573mGTe9KAgbZmUyB23aCxxL+EQO/Koj0dCWM9ocPzCHNgiiESd6VZF55DqalpRdOWLoiHaFn2X2zbCbYyQ374Gk8GaPHBpcN2W0J9L0guEyRhmTp4Uwy0+dZoDFepwQ0nlmMRFtuCcWvBMVQgaLg5MiA//8+jRo0VFrF4CF0cLW7kIYhm5dx/MuSF3JjnDf1yQzWHRrZtO/HTYN/nMtmjGIJ3A/R5dymx5F6NeoCNCqHR1VPL+Gt3kUZ/1xT3JxHmJQrJj68X7eA8jRGLHi6pWPyiqup/TGNnMlpWAjvghnKOvE+38tfLY0da+PH2AoHn8M25uMv0Dzg8od+hyOANnKkDlphz0SSOzVX+ld6tC0JgYv50/rbSDIWoDIM4ZSxlEPiCUagGbdNjaWJr4AYUG9cAoZrash7UbVK1KkttrR1bY6OzFlVRCcBLl62WOzdbTxADhRRTY8XMaI5vZsvjWbgFlrKghEjJhmieOzEEeoDCJR28Y/MxWShuinJLSlFKGJfowd/vt/mEWUCgbSNj+IY05NqQscW5mnNOgCCgCNgINdE1UWNUcT1WcZQs6NmasA6jlf7YRDsJ+VSxom4QNCIXinMcIk0K3ioAioAiMHIH00aM84KCFoq571bOiR150zUuoigVdcy20QkVAEVAEXDdG5ugRJucxjv8ZT+cD3BvBgkkJOljtpdIqAopAIQRcN2kunaZcVxdFzdJWhdIGKE4JOkCNpaIqAorA4Ahk+Qtm+J8TkycPnjAgZ5WgA9JQKqYioAgMjUD27Bn5ejDCH/mEIShBh6EVVQdFoJERMH5m3qYPHKD49BkOGiY+wNgoQQe48VR0RUAR6EdAhtfxtBIx9xsM+YCt/3Qg95SgA9lsKrQioAh4CLijNNI8K2akiecwcYfXmY9WvHQB3FGCDmCjqciKgCJwLgKZ4zz+uaXVmoa3Nl/7nStJ5WKUoCuHpZakCCgCdULAcW90UdydLTP41OwAqQRdpw6l1SoCikAFEHDHP2fPnJY5mhPmBWGAP06xUVGCttHQfUVAEQgWAubrQZ7aONrW6nw1GILRG6YRlKANErpVBBSBwCKA1VPivEhI2IISdNhaVPVRBBoMgVxfH7s3+PNunn9Dguv2CAMMStBhaEXVQRFoRARcIk7t30fR8RMoBhcHQkj8z1BFCRooaFAEFIHgIeASdIanF42N5bmfedks+ThFCTp4bakSKwKKQIgQADnzklbZ3l7K4uvBjnGiXOF554Ort1rQwW07lVwRaFgEzCpM2RMniPp6ZfVuASNE1jP0UYJu2C6uiisCwUXAWMoZnl40Yl4OBledopIrQReFRk8oAoqALxGAewOWMm/T+zopOXeeL8WshFBK0JVAUctQBBSBmiOQ7ekhrKASHTWq5nXXqkIl6FohrfUoAopAZRBwR2+kdu+Wl4PRkEzOXwgcJehCqGicIqAI+B6BzNHDztqDIXsxaAOvBG2jofuKgCLgbwRyWWd4HQ+to2zOG17nb6GHL50S9PCx05yKgCJQYwRyTMoIGQyvi8UoPs4Z/4wXhmEMStBhbFXVSREIIwJMwmaVlMzJE87KKca9YbYh01sJOmQNquooAqFFAFYyEzEm588c5MVhJ09xVA2p9QzllKBD25tVMUUgnAhkeeY6uDSiLS3hVNDSSgnaAkN3FQFFwMcIuG6M9P79FJswMdTjn00rKEEbJHSrCCgC/kXAdW9AQHzeHcXsdQhWvBMRrr9K0OFqT9VGEQg1ApmTJynHkyPF2ttDradRTgnaIKFbRUAR8D0C2VMnKZLN8vzPStC+bywVUBFQBBoAAcuNIe4Nd+7nBtBcR3E0QiOrjopAGBDA2oPZ48coMXWaow5b0mFa3qpQG6mLoxAqGqcIKAK+QyB79gy/FOSxwa3u2oO+k7DyAilBVx5TLVERUAQqiYA7vC61t5MiY8ZQpLnZGb3BS16FPdRVwwjxpNsaFAFFQBEoikD/HBu5M2cp2tYmKWVx2KJ5wnMiXk1VsuwjMkvTmDXEog1w16smplq2ItBQCICf2Y7D3BtZ/sQ76X7ebebkCDsWVSNoELJNxjZR2/thB1j1UwQUgREg4Lo3smw9g5Rjo0e7H6eMoMwAZa2aiwMkfIKnBOzq6iJY0vv27aPjx497FjUwMkRttgHCTUVVBBSBWiHAxl76yGH+vHuCVWNjuEcrbkGDjGE5HzlyhD772c/SHXfcQW3sN3r88cdpNN/93ve+99GUKVP4C82c/IB4mtcVk+A+zsjnm07M0H9xh8U4See/49XGcViD6JunXNj1NY0LtcOuq+jnduaw6wv9XAtZLmAcI9htzOcxe116315qvfKq/vPg53pe5raMjlRV+VtRgjZujUwmQ88++yzt50lNDh06RLt57bB77rmHXn75ZXr++efpzjvv5DbguV3dxgFxi0/J3BS9RitNZ+TN8c/zS5WZv7RafJSKdeW7oCNQ2HWFlo2kL7cn+nHD9GVuXtE1Givan7NnzsjMddFW5wVhf993stTlL+StwbVXUYI2hLtp0yaxiu+99146evQoJZNJmjlzJh07doxeeeUVwRNWdg+vyrtjxw569rnnKMqWd663l88Zy9qw9WDwg+SZnDNpyvY4edHYIH9Gb7CMATqXpwurle3upkg6RVHG1XmbHRZd85vFuYmLvnzTjyYSIWtbW19uZ1zwvHHal1er5sVQz23fvP5gF+Ht+zmNkc3dSn/u4v6cKNCfuf25zft2bKOYrJySY2x4qlG+5osHU37xFI7pPdQ1U6wcliked+TAk3+VSbpiBG0sYrgr1q1bR2vXrhWXRmdnJ01g3xGs6rO8jliclbMD4lP8CNMUY6sB5Oyd7N/zogrumHS8NbuSbsBBwZzBibR0yaFjuboiGr+IdT44So1A0pDqK2pZunm73o6LWf5xISj9nMbIxlvpzyw/oszP6s855pPM6dOUnDmLIrzEFb4mdP2YhZR240z5gySRygY7j3MFyvGiXK7yjocqa3jnB7Ll8MqQXMZ6hmV8++230/XXX08PPPAAzZo1i5cOi9H9999PO3fupBtuuEHSg9ATfHecP38+jRo1ip555hmipmbBfqh7W76YEddilgm88YiYnyB0x7iLJ+ROHn5d2cOR4eGaTUm+QOMN0LZMC7zuHp4WYD02QvvC5yy6svGWr2+an8BxRSemTxfL2XNz1OmaNnwcbWmV668giVdQtooRtJEJBN3R0SG/2267TUi4hVc+ePq3T9N1111Hy5cvl6QgdGMv9/Fd0XFLmFJK3QIuLgcEbX5VfuQoVbKqpjO6ohLsh11nW9+qAuuTwllf/HPIyunjPpGsOmLY7Wv6M7s84WvG7HXYRpKOu8fzzVdHkqFLNfLZMg+da9gpKk7QkMSQ7ezZsz3B7rj9Dm/f7Izc1rXut9auKV+3ikDwEWiwjm3UZVJGyLAFHXcnRzJP6cFv09I1GMzbXnopeSkBpFjIfJcBWeOH4Xf4aVAEFAFFoDgCzNCwTjnkUn2U4YEFyRkziicP+ZmqWNAGM/uOZ++b83iM06AIKAKKgIeAS844zhw9JjPXwT8tIeyuPEfLAX+rYkEPqEEPFAFFQBEoBwGXpPt2bKfY5MnOOHiLuMspKuhplaCD3oIqvyIQNgRc/3Oup7shVu4erPmUoAdDR88pAopATREw04imDx7geZ9bKD7enX+jAd0bAF4JuqbdTytTBBSBUhDInDwlE/NH+GvZRg5K0I3c+qq7IuAzBORrQf66OMOz18UnTvSZdLUXRwm69phrjYqAIpCPAMY/m+F1PNdM9vQpnn9jvJOqQV8QQnklaKcL6F9FQBGoJwLywaTzlUrmxHGKjua1BxvcvYHmUIKuZ6fUuhUBReCcyY/S+/dRfNo0d64LBqhBXxCiayhB6wWiCCgC9UUA1rP7zRom58eMddG2UY5MDezeAABK0E430L+KgCJQVwQchk7xIh+Y7zk2dmxdpfFL5UrQfmkJlUMRaGgEHP9z9sxpZ2pZ83l3Q2OiFnSDN7+qrwj4AAF2Y2CVEkzOnz11muJTpjpCwb3RwP5ngKAWtA/6p4qgCDQ2AkzE/Hl3jpfAo75eio3X4XWmPyhBGyR0qwgoAvVBAGsMsrWc5o9TIq0tshajCNLg1jMwUIKuT5fUWhUBRQAIwI0B6zmbofSBAxTrGOe4NTB3vBK0ErReJYqAIlBHBJiEMVd8tqubJ4DmldvN8Lo6iuSnqtWC9lNrqCyKQCMiwBZ05ugR4gVMKY75nxHUehYYlKAFBv2jCCgCdUMgx8vhnTlDUV5cWohZR294TaEE7UGhO4qAIlBTBJiIsUpp6kwX9e07yKunnOdU73yzUlNR/FpZVdck9KvSKpcioAjUHwHwMCzEZKSPXxaeocT4dvniOxLBGefDlfpLWV8JlKDri7/Wrgg0LAJ4ObhnH1HnMwcp3j2RFvckqLWNKJuLUFT5WfqFEnTDXh6quCJQHwSMi/nNLVn62P+Xps4fH6BdNIv+O5PyH783Q83NMRl9p+8JdRx0fXqo1qoIKAL06GNZWv3jEzRnZTedpFH0oc8xWXc6DmgMg9agBK19QBFQBOqEQF8GHuhOaoqBjTuITvRRJqu+Dbs51MVho6H7ioAiUDME3rUqRa+N2U8/eGEB15mlf/hUgqZNdQaWqXvDaQYl6Jp1R61IEVAEgEDEGatBF43bQ5/857F0a3QmjR0docsub6JRPE+/+/W3gsUIKEFrN1AEFIHaIWDeELKTufuNt2jhOy6mRROTXv05Po/RHRocBPRDFe0JioAiUDMEQMAIfbt3U3TMGMqNn0zps2nK9LLlzPFKzQKP90cJ2oNCdxQBRaCqCMA6xsx1PClSavs2SsyZwxPZ8Zhn9j9ztJJzAfCVoAuAolGKgCJQBQRc6xmrdvP8opSYpBMjDYWyEvRQCOl5RUARGDkCIGeYyRz6duyg5PyFRLHYyMsNeQlK0CFvYFVPEfADAp7vec8eHlGX5XUHp/hBLN/LoATt+yZSARWBgCPg+p6hRXrPLopPnSqLxAZcq5qIrwRdE5i1EkWgcREw1nP60EF+QZil5MxZDhiuT7pxkRlacyXooTHSFIqAIjBcBCzruW/7dopPmkSRpianNB3vPCSqStBDQqQJFAFFYKQIZI4do1xXFyWmTXeKUuu5JEiVoEuCSRMpAopA2QiAhF0ruW83+57PO4+i5ltutZ5LgrNqBA2/U4YHpGf5jS32scWx8UeVJJ0mUgQUgcAjkDlxnDLsf07MmOnootZzyW1atbk48D19zBrnaH9fD5K2j0uWVhMqAopAcBBwreQUD62Lse9ZrGdIr9ZzyW04YoLOt4gN8b755pu0YcMGmjZtGq1YsYJ287f3GzdupIsvvpjm8Ceekg8f348k4E6MxsZ2pGWNRA7NqwhUGgGvP2MngDNUuNdm5vRpyhw5TM2XXiYI5fhJGp97aygNgRER9GCWcCKRoHnz5tGaNWto27ZtdIxfErS3t9MDDzxAf/qnf0rjx4/33B2wtA2xlya2m8q9EyMvL2PWQHdm64JtBGvEUres/hHgxP3XQ0CVd/tlaudOivBCgzGeGAmGVEFyRlqjZiP05zL65bAJ2pDz4cOHxVJu4qEzl112Ga8n1izVz507lzZt2kR9fX1iPS9YsIDuuusu+ta3viXHIGj4pRF6e3tdsnasBadsPmGMB8+a4Dg05IBjJmfXgjYyoSOcmy4vHx+eE/LLPicBR9Q7DTow9DMgyH4hQa24WspsVVuRXaOvafNCbVuRiupYiGkfoyu3rfRliFSsfU2efLHteHvfTmfH2/sjSWPnxT7rku3pYd/zAWpafomcHXBdSxr3D3Q0epqt5LD+FJPTSlL1axMyGFIqJqctTwX2h0XQINYoP6Zg+/3vf5+OHz9O1157rUe4Ri4QL4i7i4fXGOIexW9xQdoIqVSKXlm7ln6+ejXNnz2bct3d0jHRWOZ6NGUV3XJHwOxY6AxRnoBF7tA1Aq+oTNU8gY7f1U2RRJwiybR8NlvN6upeNuub4f4TTadZ50T/hVx3waogAOuKvp+FvokUZTMpbt+Sr4QqCDSCIvmpuJfdnBSPU2z0aMqcOS1Pyedqw+3bDX0T3J/9ry9uMjGWNdt1lnLcJ8XFOgKYhspaNkFDQJAzRmQ8+uij9JWvfIW+/OUv0xVXXEGtra1C0kgDcr700kvp7Nmz9MUvfpEmTJhABw8epL1790o8BItz4y1nn3SOO+ZOHsQeSTaJvjl0SrlbueKjVe1jN1pYnKcrjKQzcmOLNrH1zrLJ3bhY+kLxpjxsi9VV7zSeXKwvf40FsoomkzwpGD+FDKWTn/Wycc3fF515Okpu3wjr6hG0ra+HS35m67iWaaxqy9p1dUUHjJr2ZePmnPa1dbH37crseHu/2mlQPurjazDX10tpHlrXsvJyabsobj58rXrXlycXty/3YbQt2pgJpHh/9vK49dj9AHUjVDlNBOUzb0VgMGAQBHNdNUPZBG0Lgxd/8DPDEoafedGiRTJyAxbyS+x73vzGG3Lu85//POGl4X333UcXXXQRzZk7R4oBQeM3j90hezCJijvqI1LuO4RElEkLFhbUKdRqttTB3xfrmXHDhVDQpxd8FQdoEMEFgScGt38MOBnCA0dffloIYvuCsJiMUwcO8MiNyRSf7EwpKjfXIm0lbYv2hXGFn4+DoeNI3GmfaotaNkHj5QWsZ7zYmzNnDr3vfe+jW265hX7961/TjBkz5EUgzi1ZupRmzJxJo/nxZhIPsVm4cCFdddVVNHHiRGpiS9kOZ86ccSwFEymNbA5466FixWEX8SyPWBmwuvGTuzRvhafxB4ncgF2Pv82OdR7JfJmGZfUw4X3oae7cXryro2yg22B6IVFQ0kB3tqqgM2anzNd3QHtBrwKhlmkKVF9SlMgIXeF/hr5sVSKjp69zzisLJ4vphXiEmqZx5cP1yAZa37at1Lx0mYjhjdyALgicZoBe5tot1L5ODudvXfTiqm2cZZ/lh6Vv9LFlrPB+2QQtPhjXknnllVfobW97G53moTQYpbF161Y6evQoXX311TJKAy8CEZCnra1Nfua4/y214+qwmNNpQMnp/kFHKxRMPBrc7CMdjr1g7Vu73ukBGTnWr2lsnWwZZd+OMJrlxeUdOqnyIvMO/ZOGBTOy2ThAQBPvCFv4by3TFJZg6Fgjo+jHB0ZPs0UJ9r4cFynWlFXrNHIzYet5316xhuMTJ4mA3lOeLf+AfUsPkd1WIP+ce1wkyYD+UJU01k3GEq1au2U/T4BsEbr5hd6WLVuoo6NDXvphjPPrr79OM9lqbmlpEV+0+ZIQ6e0vCW1yxjmnTFdxRGgoHQGFrXSsApUyYA0LXnDdE6mdO6hp0WIHbZcvBoXeVtXeHzRTY5wcNkH38KiJMTy2cezYsTKK4xvf+AZhaB1cGSBg/ODqwAtF7GN7znjnYne4xsBetVQEwoOAS8RYLSUS43dLupxVRdq2bIIG0SJguNzf/d3f0RReGWHz5s0y3G79+vW0gxsIwVjacqB/FAFFILwI2Nbzrp0Unz7jXHdMeLWvqmZl+6CNewJfCsK9kebhJpdffjlNnz6d9u/fL64MSGzSVVV6LVwRUATqjgCMMVzvKV4MFiNQEjxYQAKIm+M1DB+BsgnargoNg2FycGsgYMQGgmkwOdA/ioAiEF4EQM54quZt6q23KM7kDJKWoOQ84nYv28Vh12isZLwABCmbrYm30+q+IqAIhBABvu4R0jwhUi6d4sVgz3OUdOOdA/07XARGZEGbSo1fWonZIKJbRaABEAAJu++kUrt2UWLWbIryCC5Y0+raqEz7j8iCrowIWooioAgEGYH0oUOESfnjU6cFWQ1fyq4E7ctmUaEUAZ8jYFnJqb2d4toQ69nnYgdNPCXooLWYyqsI+AiBzIkTlD11kpJz5jpSyWRHOnKjUk2kBF0pJLUcRaCREHBHaPRt30bR0WMoylM5iO/Z9Uk3EhTV1FUJuproatmKQBgRcEdoZHkqYfiek3PniZYYyaWhsggoQVcWTy1NEQg/Aq71jJEbsfYO/rWL9exNihR+BGqmoRJ0zaDWihSB8CCQ5cnSUju3sfU811FKreeqNK4SdFVg1UIVgZAi4BJxqnMPRXk6UVjQEjAPu4aKI6AEXXFItUBFIKQIgJzZvYH1P+HeSM67QBSVBTMGTMQcUv3roJYSdB1A1yoVgSAjkOLl6bBMVXzcOPU9V7khlaCrDLAWrwiEBgG8HORxzuk9u6lp8UWhUcvPiihB+7l1VDZFwCcIOG4Mol4e9xzh1cbjEyY6krkjOnwiZujEUIIOXZOqQopAhRFg37MzhC5H6R3beVKkWU4FOnKjwkCfW5wS9LmYaIwioAhYCJgPUFKdnRRpG0UJMymSWs8WStXZVYKuDq5aqiIQDgRc6xkujj6ekD/Bi0LLFKNqPdekfZWgawKzVqIIBBMBYz2nDx7gF4QZik+eEkxFAiq1EnRAG07FVgSqjoDneyZn3PMF83l4XcKpVt0bVYcfFShB1wRmrUQRCB4CxnpO7dtLua4unvPZtZ7VvVGzxlSCrhnUWpEiECAELOs5jQn5p06lSDLpKKDWc80aUgm6ZlBrRYpA8BBIHz1C2d5eSsye4wiv1nNNG1EJuqZwa2WKQAAQAAm7VnLf1rdkQqRoc7MuBluHplOCrgPoWqUiEAQEMidPUu5sFyV5tW4Jaj3XvNmUoGsOuVaoCPgcAdd6Tu3mCfknTuQlrUbrclZ1ajIl6DoBr9UqAr5EwLWSM6dPU5oJOjF7tohpRnT4UuYQC6UEHeLGVdUUgbIRcK3nNE/IHzvvPIrxgrBYDDbixpddnmYYEQJK0COCTzMrAiFCwLWesRhset8+b0J+Wa1bCbouDa0EXRfYtVJFwIcIuKtWpXbtlClFY2PHyvzPMveGD8VtBJGUoBuhlVVHRWAoBHhkHTsyKJtKU2r/fmpatNjJoZbzUMhV9bwSdFXh1cIVAf8jAM9GTuiZKMMrdcdGtVHMXc7KjIf2vxbhlFAJOpztqlopAiUhAHKORJyXgD1nM3Ty9a0Un+18Nei6pEsqRxNVB4F4dYrVUhUBRcDvCDjkzFLyzlPPRehn33mLUlta6LqPjafbbiVKJCL6frDOjVgVgs4fM4khOiZu4HAd961EnUHQ6hWBRkQgy9ZzjK3nN96M0j1/cYr6XnuDDo5eRv/rLibsp1N0zaqEEnSdO0ZVXBwgYfsHHc2xIWqJ4w6iQRFQBOqEgOvDOHyMaM9rr9Oly9ro0kswpWiadu3Wi7NOrTKg2qpY0JlMhrq7u4WU29rapMKzPLYyHo9TE68I7AXXgB5oVXtndUcRUASqiAA/13LpUZoe30l30AF6cMM1csyvCmn2rKrYblXUJpxFV4Wgd+zYQY899hiBlG+66SYh5p/97Gc0lsdVvuc976GJ/H0/LGljTafTaQdd9BeQtntnLwlyDANCeue/ZC8rf0mV+CiR6JsnTzl45WX1/SH0NY0LYcOuq+jnduZq6YtLhv0b0ViU0idO0rSeN+i/PrCSOp4dQ12H+uiuOxN0+cqY4C7wixz4U4VQqIIgtHGNZKwKQY8ZM4buvvtu2rp1Kz300ENiTd9+++20fft2euqpp+QcrjljOY/myVhkWXdciwheozmHQ/1F3hz/nKXhy88/VPm+O8+6eh8PlImV73QpRaBG0pfbE/242n05EmX7mReC7XnlZUrMOp9WrZpGl96I71LixKPsHEOH28brXubaLKW9yk0j7YsbAgevQufQt3+jLG8NZK04QcMqnjRpkuC6adMmbvCsuDUuueQSsaBhWSOgg6R703T06FF6ee1aSvf08MJnKT7h3N29HiKpi/xhkid0tHRGJhUHaUnHxt2tmh2qiDi1iWZ9+3plAU/CCheMby06Sm10y6vFbV9MGB/FPtbDC2vbQj+3z0LfCCxcWDF4k2f3ZStdHlr9h4OlEfy4QCaYnlfXUbSlmQl6Nl9DKWprdvsSX4YZ3gVvQoQB9ffX4uwNdR6phkiTY32xIO2Q/XmIckSgKqcBv0XYVZvrZb5iV261Q0UJWoR37yrr16+nLVu20Pvf/36xovv6+lifgQqBvDs7O+n1zZtpyrgOymXY1cEAiOsD5RiwsTUBndWL550I9yLOl4ObhDsZLGnnInZ7tck7IJ9bmJtEykOU39MYlYxLCLoagjZ6Gj2wNTjZ+yadrXv+/kjS2Hnt/WIyDJqGT3L7om2lXQu1D8o18WZ/WHXVuRxDnJAd7Qu92QUh7QsFbR2RxhybrYnD1sZ0wL5zEEkkqffN31H6wH5qveZaiuSybOCkKCtLlPKNga89XFY5c7maOsx2sLrKSYO0HHBzEH1K6c9OFufvYHWZdJVKg/K4LOEm3pU+ybhV2ziqKEEbTLZt20Yf//jH6cMf/jDhJWGUgX/00UfpwIEDNGfOHEkGRRNsEcGyHsdfLb344osUaW6Rc8C0rACS4jqiyI9GboAQiSecFZZxzZUNWMAA4vaNNCXZcnFXlA6Y+OWKC/8w1v+LmhW0yy2gWHr3JpDhp9bUzp3UeuXVMlsdXB2xZBO5ToZiuasa7/XnqtYy8sKN8Rht5TaKMX0C0yqGirGZfCrqWs9dvALwXXfdJS8JQda33norHT9+XF4Orlq1StThezTffBxmAYFXJLjlVaQsXxdiMbK162uRRyJcI+ho48P6VkVlXB/8FNvz6ivUtGABxTo6hGA8f7ctQ033q6VwTZWoSmUVs6BBuCYsWbKE8LODsZy9OE5uHhfScFEM507kWgSSFTcy99irI7Q71l27EXS21A1tk9qKsb64HuSKqkT7WmV0vfg8RcdPpOTcec71Ytdbr33IZ9rYkrVe4vip3gqZrsNTyVjQZlt2KcZiRk/uvz+UXYxmUAR8h4Dpz6aPV0DA3rfeoiy/jG+2jacKll8BEbWIPATqStBwi0gwd8884fRQEVAERoYA/Mt4kZU+dIhSb2ym1svfLqMQGudpc2T41Tt3XQm63spr/YpAqBHAy1V+v5Pt7qKedex3XrGCoviyV90IgWl2JejANJUKqgiUgQCGgLkv37tfeoliU6dSYuo0JecyIPRDUiVoP7SCyqAIVBIBWMjuS5meja/JkL2WJUsrWYOWVSMElKBrBLRWowjUDAHXhZHat5fS/CFY87Jl4odW10bNWqBiFSlBVwxKLUgR8AECrmsjc+oU9a5/lZovuZSirex3dl8W+kBCFaEMBJSgywBLkyoCvkZALGf+RJzJuHfDekrMmUvxyZMdkXU4na+brphwStDFkNF4RSBICFgjMzAJEsVi1LRwkaOBdS5IKqms/J5XQVAEFIHwINC3ezdhrg24NiSoayPQjasEHejmU+EVAUbAtZCzp05S34ZXqQXjnZubxdVhhtopTsFEQAk6mO2mUisCDgIuOed4Ot+ul9ZQ8qIlFBs3Xsi5/pMgaSONFAEl6JEiqPkVgXohAHJ2Qw+/FIyOHkNJfjEIi3rY89uYAnXrCwSUoH3RDCqEIlA+AmZu4r4tb1L29GlqvWxlfyE6aqMfiwDvKUEHuPFU9MZFAEPp4MJIHz1CvW9tcV4K4tNu1+XRuMiES3Ml6HC1p2rTCAi45JzldfF6XnqRmpYuo1h7u5JzCNteCTqEjaoqhRgBWMjuJEg969ZRfPpMSs6Y6Sisbo3QNbwSdOiaVBUKLQL2S8HXN8rK0k2LL3LUtc6FVv8GVEwJugEbXVUOMAJsJaf276PMvn083nml+KHV7xzg9hxCdCXoIQDS04qALxBwX/5lz56VeTaSbDlHR41Sv7MvGqd6QihBVw9bLVkRqAwCLjnneEXu7rUvUXzWbHfyfWc5q8pUoqX4EQElaD+2isqkCBRAoJfXFIzwJEjNZhIkXSm5AErhilKCDld7qjZhQwDzO7PfuW/3Lp58fw+1rLzc0VAnQQpbSxfURwm6ICwaqQjUG4GII0AkSjL5/qaN1HzpClm+Slbqdofa1VtKrb+6CChBVxdfLV0RGB4CsJAReNu95gVKXnAhxSdOkmOdBMmBphH+KkE3QiurjsFCwH0pCKG71r5MsY5x1MQELUEtZweHBvmrBN0gDa1q+h8BMzddNJqjaDJJvdu3U+YYT77Pn3JL0I9R/N+IFZYwXuHytDhFQBEYJgKO1zlLKWqiCH+Iktu5hVovv0L8zvoxyjBBDXg2JeiAN6CKH3wEjEejpztLv34yQmtfPk4z97xAV3/gErqgvYMVhG3tvjQMvrqqQRkIKEGXAZYmVQQqjYBQr3Bvlp59MUq33HyGFtAL9DuaRB9om05fvShH7e0RpehKAx+Q8tQHHZCGUjFDioD4lTFiI0qdnad4+wxdfXcLXfb2xfTdr2Xp1Cl3NEdI1Ve1BkdALejB8dGzikD1EPBe+rGddPwoTdm9husaTz/aPJtOvR6ni6/OUnOL2lDVawD/l6yt7/82UgnDiIDreMbagb07dlDP2hfp4psX0X/7nytobFecbrglR//j75M0cYLje9apnsPYCYbWSS3ooTHSFIpAZREwbwW51O71r1Lm4AFquewKmjxuHP2X+Tm699YkJVubqGNSTF8NVhb5wJWmBB24JlOBA42AS865vj7qfpldGnzcuupaira0iFotLRFqmcoPtgnHcra4PNBqq/DDQ0AJeni4aS5FoDwEDNOyryJz/LiQc3zKFGpewh+hwH9hiJuHdWT68MoQH6s4p8qrSFOHCQEl6DC1puriTwQMObN0KZ6RrofdGk2Ll1ByzhxHXus8uFr8zY4B7U99VKqaIaAEXTOotaKGRACTHplFXl/fxFOG7qaWK95O8fETHKsZoAgjNyQ6qvQQCChBDwGQnlYEhoUArGJ8XsLknOvtpe5X1rLvIs3+5uscf7NF3MMqXzM1BAJVIegcd86sO11ijFeAyD82yEb0HbWBQrdhQkDImRXCXM4nTsgyVbFx46ll+cUOYfO1oVOGhqnBq6dLxQkaZIyxnSBmE/KPTbx8v+od6I4iEAIEQM6uyyLV2Um9G9dTcv5CSs6d5yin5ByCRq6dChUlaEPO3d3dtHnzZhrFqw7Pnz+fevkRb+PGjdTW1kYLFiwQAhcV9UVI7Vpaa6o+ApbboodXQEkzQTeveBtPtD/RqRvkrfM5V78dQlRDRQna4AJCXr16NaVSKfrkJz9Jzz77LK1du1ZcHTfffDMtW7ZM9kHoBYNlhRQ8byJNOmyLFGWS6jaACJj2heh+b1+XnHPc57Hydq4vRa1XXU1RNlJkCB10KOdloNHXxgBlaGgoBCr6qTdcGfA9t7e303vf+16aNm0anWAf3IYNG+jP/uzP6Prrr6c1azDfgBOQHiGRSPRb1YgotSO76aQc7JeaD3UEOli6hllno1sQ2pYt4xOlEm4AABzGSURBVMzp09T1zG8p0txCrVevcsgZ/axc+Tm9uTYaok/b+GDfz8HIZ8tcRXkrbkEbq7iFv4xqamoSwgYBjx8/ns6ePStWNfRBB+zjr6k2bdpEv3r8cZoxaRLlurtEVVNGSXpzOblMhrI9PWyouC9fYHWENvBNkHGKxOP9E7mHVldWjNs303WWoukURbgfedaoH3RGP2NijsTilNrbST2vrqPEvHnUvGAhW9C90i+HQ7C5ri7KwmjhVVV8pW+VMM+ySxRt67v2LaAvuCnGsmbPniE8LQ2nfQsUWzSq4gRtaoIiIGCQc4YJ9NChQ3Tw4EFKotO5AS8SZ82aRRcvW86TeR1m0uELEDfQLHf8Um6k4OEorMkoReMZyS9vx3HhlJLfCBKoLVtXae7MjFUkwVji0TqsurrtK7pCX/QPv/g60EcTfPnwtveNzZTasZ1aePWTxHlTKYebCfqg9aJ8yC4GXWGVQT9X19C3r4DC/TmVDk5/RrvKtcf9Ee8TcFzFUFGClrsLd8ozZ87Q008/TY888gjNY4sCv+9973uUTqfptttuE3VM2klsOV966SX0FKdnNpdz5fINLJgsW89RttgbIUQyvCQSLGj+NUKINiXFmkQ7+ynkuD/3bFhPOXZttP3e9RTll+AIQqwjEDSaZMJiQ6ZR2jfC122Eb3Z+a9/BmjCSZK4p5wY8WGGDnKtKj4f7Yvny5bRw4UIa2z6WzuM5B+bOm0st7Js7//zzRRykAUkjnGHXRw6WoAmIN9YEzEM5xklD3chn4rkc5IVFg625q0l+U6C7PaecvPM4HDIN1+NW3y9PfjmDpTGZy0nTr6vUBD0hJ4LIC1xMuRLr/rHiiupV6TRcnlekaS8jkzlRThouw7QvRm4O0BfH/JNqTF0mwtqa3fz+42XmBAPSeAduBbwx/U122c106hSvuP0SRVvbxN+MixX9sH98symj0BaFIEBm97ylF1x10Fk08uKRntOa4ow8hXS329reRxGmADve3q9JGhZa6hThB2lf97wvZAYwltz2NSiYVedPRQnakC6G0y1ZsmSAxEsuGniMk0iPEJc7kWkMOSHxTge2j91op+uiACcCWyu7F2+Sm61Jb44LbYdMk1dXoTIgjC3PgDTmRJlpbLlMESjX2/d2rNqsODu/lcIqgHet9MNOw2UUK8Y7UW4aK70no1vJOXWZCGtrdj2dTIS1NbuSxj5w9y1CSfP0oD3rXqHEnHnUxMNGJQwgZ8SYMoptnWxeugF6cR5zbLaSHPEmH6qwDxDvHtvx9r5k9UkaWy5bDS9+QKRIPkBfL51zqra6A2pbPiND5bcVJWiIZ0jafEmIuChbtebY/oAF5xAcSxqmgYayEfAsqrJzaoZSEYC1hCczDr1v/o5S27dR8/JLKH7eeU4JIG/3vBNRib8NdD3Yqlo3wkqgGPQyKk7QAAQknU/E+cdIl5NHF+xpUAT8hQB4Aj/2W1A0xuTMJI35NDKnTlLLlVdRbMxYNwGnqZE15S+EVJpaIFAVgq6F4FqHIlBNBMC5Du8yOfOwxrNrXqRISxu1rbrGeQko7M0SKDlXsxkavmwl6IbvAgqAjQDcbeDcQ4ci9NI6oqa+fXRh3ys0dcn5lOA5NSSAnJWYbdh0v0oIOI61KhWuxSoCgUGASReDJ+CeO3YyQl/4chfdcuM6etdtz9Lf/u+FtKfJIeesjNGvzQuiwGCnglYNAbWgqwatFhwIBGxXhfBumvau3U3/9pU36Na3RelA8mr64WPj6a4NKZo7O8HvAjE8VA3oQLRtCIRUgg5BI6oKw0DAJmbOjs92+/bto9j+HZTZlaGTNI8ePjqLprbiIbOPPx5zPoIy2YZRo2ZRBMpGQAm6bMg0Q6ARAMPi5w6LAzFjHo307l3yrVNi/vk0e/5U+vujEfrU33TRPsrSn38kTsuX2vObBxoBFT5ACChBB6ixVNQRIGATM/uZhZg791B6zx7+ZDdKibnzKMGzL2Jel3au5iP/OUPXX9vEq1RFaO6cCE2Z4ryu0XeDI2gDzVo2AkrQZUOmGYKFAF7+8cgMWMwgZp64K7VrJ6WYmDH/Q2LuXIpPndb/iTaPd+bUNHZsjN52Wb+m4Hcl5348dK82CChB1wZnraUeCLisGsGLPZ7YKLVnt8w6F2lqpiTPCRPnmee8uTPwtSAYmIkc7wpxaIeKfyhoF677ikARBJSgiwCj0QFGwJi7IFxm2r6dOyjFP8xAlly0mBKTJns+aOdzQdY1j4HzDgMMhooeZASUoIPceip7PwKGlBEDVwYsZvYxp7ZtlRVOkguZmHlWRc9PgfRuWmdH/yoC/kNACdp/baISlYOAIWbXQey5MthiJp4/Orn4IkpMVmIuB1JN6x8ElKD90xYqSTkI5FnAznC5vWIxY+GHJC87hdVNvJBH5F687igCPkZACdrHjdPoooFT4YiQ93eY8ZMP2HnhwGIsZl5WDeOYUzt3ygok4mO2pwFFaqR10zuZ9a8iEAwElKCD0U4NKyWWnIy28p+YS8zupPS53l7qwzhmJudILEbJ+QsoMdWymMHq+qavYftNWBRXgg5LS4ZID1AxU7IYvfv2Z6lzax81j07QBbx4SUuij3p37Kb0PiZmXrwzOe8C5wMTo7+Y22B1nQfMQKLb4CKgBB3ctgun5OzXkFnl+MvqrVsz9KkvROj/fBd0fZL+474jdPOSvRRLtlDyggt5VAavaGJcF0rM4ewPDa6VEnSDdwDfqO8SbI4JN+pOe/HkC1Em5/10/eLtdOp3R+m+j4yht728kOasOI9nyABps4+a88nalmox+6YpVZDKIaAEXTkstaSyEWBnhvgz+l0SOX7pl+3ppiz7lpNrOinORDyGJ8c4zuOY12yaStlJbVJLLs0Zmci9LwHLrlszKAL+R0AJ2v9tFGIJmZj5P4bIpQ8douzpU5Q+cIAJupdaZk2li9+7iBasaaEHngMEzfSZv0/S5A58gx1hYnbyhhgcVU0RYANFgyJQLQTyxip71bBbQkj56FGeH2MX5Xp6+KVejKJtbdS8ZAlF2ztkZMYSzvDDfyV6fd0JGt0Ro7dfGadRo3mgHRvP6tHw0NSdECOgBB3ixq25aoaQUXGBsceZ48cpc/IkZY4eoczhwxQZPVpe9EXHjKb4hIlOHktojHleelGEls5PED4+ccztc5JZOXRXEQgXAkrQ4WrP2msDUsavACFjrHK2u5tS+/ZS7tQpsZqJXRNxHn3RtHARRZubB5rChuDNyIwcUzR7MjK9XDyfiyaVnGvfwFpjPRFQgq4n+gGqG9yJgRYIkQiTpfmiD74GQ6h8LsvuisyxY5Q9xZby4UM8aVGGouPGUYznw4hPnEjRUaOkDO+PKTSvHDnvupnF3awuZw8y3WkcBJSgG6eth62pM+E9D5pwh78xI3NZ+HFg5saoi8yRI/yCb79jJcPabRvF82Esohi7MSKwlE0A09tBnck2GrqvCAxAQAl6ABx6IAi4Vq1QKS8BhRETPLCCOjuJUhmiqTyd8pjMEeo7foL9yUcpd+Y0W8ajKcaWcnTMGMefnEe8WMnErGqiKCsCikBpCChBl4ZTuFMZqxZbEKtLrq6NTD1nU/STn6bo//n0SVrYspvuuf0s/f41ceoYl6TYxEkUX7iQ/cktvIQUXuRZAeXhxy4QzJehQRFQBMpDQAm6PLyCm9qQMDSwfMaikDl2t5kTxynb1S1f6cVOH6Od6w7S/7mni8bzqMzcoon0oS+dR79eNY1+77IxA/Gw6zD1mLIHptQjRUARKAEBJegSQPJTErgdhAdhmBYSzCZJmxztfeRDOv7B9ZA+cpgy/IFILpMWHzLGKMv31pwnMWUcnZ12CT3CLwWXr4xRbEIL0eYonUw5XQfFYDichPw6nFj9qwgoAsNEQAl6mMDVKxtIOdaEP44EoEYMQZMAgixCktmuLsqePSOknDlzltL791GO45A+2tpKsfHjKcouCrgpYuMnSJyUzQVPG0109z199OPvsQOaZ8GgWVm6cK7lsihSpyOU/lUEFIHhIqAEPVzkqpnPEG5+HUyEW97K0pZNaRrdHqWllySoY6xlR8Mi5rX4ELJnmYQPHaRcbw9bxlnKnmFyZgsZ9m6khWeDmz2HomPHCp9Hmppk6k7JaP9xGXrKJKIvfjpG77yaqLsnR1dfGaPFixyCVm62AdN9RaCyCISCoEFRrg1ZWXSqXRqIOJ+MwXhFWG/jZqI7P5ihrWtgyabpK1/tpv/0B13U1pSldE+Kskd4PouTJ7hILpdHX0RHj3Je3iWjlJg+g0dXTCiukZHF1M1bs4tM58+Lya94AXpGEVAEKo1AKAg6y3yEkWH8/URhv2ylUXPLQ724M4DI+smMIxBfKAhxmtuJm6k/48AcrBCm0hQ/MZ+JZvromYcPMzkfoz/6vTO071CMHvybNF0/PksLLoxTX4r9xZMmUtOMmez+4KFxiaS4KQYWykf5NwSTYKASJtbb2tmwnzeKzkunO4qAIlA5BIJN0MwUMn8wc53MIcxbkEcxzqsUbGKhcmFRc0cAJ7NRK/WCucDBhYInmJMAL+jgF8ZkQfIxCOTnF3QZ/ixaXBJ4WccB9SUSEWo6BuZPUHdvgk7sjFL3KCbkFdN4mFsbtTKXFyLNXBbWNlvDKGgIEkaSYsETHfjqYiXFYNJ4RaCiCASboA3xMCQpdq8mWlwiGi5EYHcEED9vDM/a+yA5mSCez2cxVwSfZIOVx/kioxuQ3yoLL+gyJ9j1wNNpwkeM/LCOszxXBQgZx06ZKCxGkdY2ivPCp9GWVtwF5L1cvCVOq8Y00U27YvTT+3n2Nybqb367mabPcep0SJPzs1nvyct5I2b2eyPbSLcGlJGWo/kVAUVgSAQCS9DgP1h1R44S/eIXfbT1zV5atCRGN9/cRGN5eK6cL6Y+MhYKJh6EaZ239yWazeXfbeqh1Y91UXd3iq5fFaOVy2KUOt1DaZ6xzZkYqE/kE55m0xafO2M6TfnsmSPxVV0CIydGj3HijfkLGYwclgzYvWAh0Tf+e5Z2/ucctY3J0YLFOWpp5k+tszwSA0IiH0g5L58eKgKKQDARCCRBG3Lu7c3R93+Yoo99FK4Afsbn37e+k6MPfZAJcRgshYl+MPoh13WWP9ToIqzu4TC9U1g2w66GSJo69/TQ1+/ron96LE1s49J6njrob76eoGUreYxwy1hKzJ3L1i/vuyES53VBmouMlDCJ8rdQUtjdOsHkO3NGlGaex1/scZnMyJIkKuxspdNdRUARCAUC9SVoPIuPIDCX0i+fQCFZuv7GPnr8l2n67S9P0x/dlKTW1gixi9exYtmtAOKFvxfES93841nWmAH555CvOAbgJwCZ4gUbCDVpzW+JpFEm6PY22rt9ApNzhKYvbaVJY5vpJ88Q3TQqQZetaIKHwbFmOXlJIZ+ETaYCljREQMj18R7/jzBPD+M+5BSifxUBRcD3CNSVoI0vd7gswzxK58+J0ON0hLK7eAwaF7RybIyymxLUwz5hkCUILMd+WIz9jbBPN9bBU1/OmCHkayxxp5XYNRDnTPhYI47fudDARkcYu4D/jOqlztdy1DnRoc2JE10nNA6lUice6SUUM+mLxZt81taQsalCjk2klU53FQFFIBwInMtCFdYLL8tAxAf4U+LVq1dTkj+KuOnGG2kcz3yW5RdlMX4pluUPKQaEgcw54JR3wGl4mC99+INxOn10FL305Dz68z+fSDd+rING87CzigTIYQUM2gD5LphH9MD3IvT/filDbx4h+uqnYrTqSgxtgK7IwH/MFocjCUNhYdh6JHX4Na+nu9UOXpxfhR6mXPl6GZXz4+WpD3VIBytSmcnswzTn6FNEhYLRftLLyFJQ0IpFVojJhpbnscceo/b2dkqzu+GJJ56gu+++m4eFOeO1WvlT4wGrMw9hVXokyNUuXRqhr31tMp0+MpbGTmom/jiuciFPDhyiWWC533Zbkq68MimekgnjiPi+I6E/y2AXRxki9hfYnwk+ZxNvtv1nw7NndIPryduvEK5+Q8nox3LJqB77pfEAWUvR38dpLD09tUruzz7Qy8hv90lPkcrvVJWgjfV8lv2/+/buo3s/cC/18tCyH/3oR7JNso8X53bu3EkZjPkVpzGTIO6ypQQkY8A62jPU0cYTFvMXc7kMuxrkLs3nTDEF2xWRUoC15V0vE/YLB7hOMLRu8iS3YD7Opo31XDhP5WL5xaB8zs114yLGxywF9atcjXUrCc3DF6/oC11xIZfaN+om9EgqRkPy2xC0L+sKoyWXy3u6HEnxvszL7YsxsujEAejPQi3s/sylma9w7VU5VJWgjexwZQB/sQyYUIWA+TjDhPzqq6/Srx5/nM6fMZ1ffjHJMgLOeU7g8Sd2EFAIb8whdvhOluYGznT3Uaw1SrE4X8iCIhLmBRNv8psk5liSo0w+YeIKpWFfBw/oYNeMUz5uphhJgeIlIM9QdQ0njSsL1vqTmxkuXsGWT5i6IYCbTuK8fU5QTK9KphE5hqirzDRZnk8kKroy4Pm4Sh+w9Pd0cXEosy4Px2GXY9dboC1EngJpjF58CvOnoHNlESckwMJAHqML7zptyREF29StRM75NY0rF6si7YuLacj+7OplwBhM90qnQXn8X7gpwwTNK1jgQzPBH+1RpVATgh7Nyx7B5/zaa69RHw9dmzp1KjUlm8TdcdVVV9F4niNi/bp1PBSNh6YJ6FnGAj1yiIAOzHfdKAMFYpahbfgwQ6yOEvIPUfxQp93XgkMlq8p5eZGJUSaGoKtSS50LddsXb3sxokZG1SAu5EG+Kk0meJFc9pvVqC/XFVJwNV7OB6A/C0GDc9AnMZCgyv2xqgQNixkKYfuud72LHn74YYJb47bbbpP+IMry3hisW4fHG5AzhwhM0lKCm17yiU/IzVdq/lLq8GMax2R3JANuYQ2mfaGj0dPEhVVn6MW6etdA2Puyq6/Xvqad/dq+pv9BzurbgLxERpUDyBlh9uzZPMriz7njRYSkEWfO9WAuCvtOhH03H9IVDW46yYt9k89si2YM+AmjK9QIs65GN1vfgDddUfGNrm6bimHj7pd0LRQtOAAn7Pa1cfCj6EY+kbn6AladoG0VmsxQBzfSEDSG2g0IpZAzMph02Np3MxM/oNAQHdi6hkitc1QJezvaCufrao7N1k4btv1G6c/DaLeaPh/DKrAtZbNvtsOQX7MoAiFFIPy+dq/hGkhVT+cSd2pqQRuL2ZNN75weFLqjCCgCikA+AjW1oPMrt0dqnEPe+YkHOW5onm+ER+BB2j6sp+xrI6w6enrZF7D2Zw8W7NTUgh5Qs33ADYQx0fiBqEt2ecBRz29TMR7R/KRYxIc4ZAUnftMP332Yh9mhDbk/iL7QGfqGvG2hMq4D6BwVffExhM1gSBGuEKT+DG6K8/A6tFHJPDWC5vIFQcNaGDVqlMzLMSxdMOQFn/ZVenL6YQlT/UyxNp6ExFgafh+WVAE4YjyPNhYyaJQwUN+6PuTWBPIg9uc27pOYoqLaoa4Ebe5A2K7jD1Vmz57NBjE+by3dAkbKKKeP8ZdIGSZoLIEV9hDj1blzPD42w1iFX1t+zMN0sawrfo0QYqxvjnVtlPaNwRrl6zYI+oKbEvxRzSlelu7111+nd7zjHVXtknUlaDNZ0nheWWT58uW0bds2saJLJWh8Qo6heyd4OakNGzbQpZdeKnc1TMg0Ep92VREfYeHA7OWXX6Zp06bRrFmzeEWXbrmpjbBYX2Y3F8OTTz5JF154oeiccpcI86XAIxAKusrMjtynn3/+eTFWYLBg7pqw9mWj15o1a6Qvz+AFj3t6/N+fcQ2iXa644grCV9IIRpcRdIGCWetO0DCWx4wZQ3fccUdBAQeLRKcGMCDoFJPyLbfcUpPHjsFkqsU53ICWLFlCF1xwAbugs6ElaIPlwYMH6dprr6V58+aZqFBvYZ2tXLmSFizAxOPhD/hQbcWKFdK+Qe3PxtisdGvV/ZmR+XXEAZbHlMmTQ09UBig8cTTzGocI1bpzm7r8sJ0yZYo8KflBllrIMGniRK99a1FfveuYyPqaj9gaoT+Xg3eErdDSHb7llFxm2uGIgTzmzoVHX/iGEIZTVpni1i05OjAsaLxJRoDFEeZODd0wJS1eyiCEuW2hH/TF47MhrEbQN8jXbrWvPd8QNDpnOQEdF+CgcdevXy9ujosvvpgm8Mx4YQ779++XWQFxM1q2bBnBmjZYhE1voxe2L7zwAk2fPp1mzpwZSn1tXXfs2EGvb3qdZs+ZTYsXL/aMkLC1L/TBzXfjxo10+vRpWrRokbxnMFiEUd9ydaq7i6NcgfPTv/HGG/TMM8/Q5s2b6dFHHxXCRho0chgDOjIs6M7OTtHXWNBh1Nfo9NJLL/HqNVfKi+Awty1027NnDz300EPU3NIsbg60bxiDaVu88P71r38toyLuv/9+6uJFnWF4hVXvctsykASNxjWPFhjqgjml7733XiGtY8eOlYtBYNJDb7wYvPnmm+myyy6jQ4cOyYD5wChQhqC4QOG+On78GO3atYuXNfsadXR0lFFCcJLa/Xn79u1CUtAVP+PKCo425UmKscR4hwS3zqRJk86Z6bK80sKXOpAEbZoBliQaFg2Mjowt4sIccGM6cuQIwdp45zvfKX53Y0WHRW8Qlnm38PDDP6fDhw+LnlgazZzDNozh+PHjQlK7d++mH/zgB2JZQs+w6otrFkPVsMXXeVjQQ0M/AoEmaJAy3gDDksSjITqxeZnUr2I49qAbyBlDCr//ve+L/xn+SQTzNBEOTfu1wPsFrL4DKwuPwm+99ZaM++5PEY49tJ8h4LG86vH8+fPpmmuukRsThtyFMZg+i/dHCxcupOuvv16Wv8MNCsHgEUbdy9Ep9lkO5WTwQ1rTobHF3ffpp5+mp556it797nfLBw2Q0XQAP8hbKRmg04svvkj/8A//QOeff748PWAIGm5UhsArVVc9yzFtB6sKY5/x8hdxcO+AvIz7o54yVqNu6IhvAvBV7S9/+UsZC40PuIADgsGlGnXXukzTX/Gye/Xq1fKiH65KvPg2T09h0ne4+AZ2FIetMKwMuDqw7qHpzPb5MO3jrTcG9uNxEJ27vb09VBdusbYyj75YMs1c3MXSBj0eL4LRzujP0DfcIUfHjh2Xl/sYkRQ2Y2OkbRcKgh4pCJrfvwiAjPEzVlU55FxOWv8i0FiSaZsNbO9QEDQa1YSwPxbZukLnsOtr2rWRtnYbN0L7Npq+5fTlUBB0OQpr2sZAAC4gjKltaWmRx2Yz0qURCK8xWrgxtAz0KI7GaCLVshwEjDWGoXmf+cxnCEPzEOAiUXIWKPRPgBBQgg5QY6mopSMAC/rVV1/1PuQ5evSovHgrvQRNqQjUHwEl6Pq3gUpQJQQwZzZGBWAagM997nO0d+/eKtWkxSoC1UGgrvNBV0clLVURcBAAOWOeB8xbcuedd8rYccVGEQgSAmpBB6m1VNaSEcB4eEz0/+EPf1gm+8eE//BDGx91yQVpQkWgjggoQdcRfK26eghg1AZcHB/72MfoueeeU/dG9aDWkquIgBJ0FcHVouuHAF4SbtmyhT7wgQ/Iepcf+tCHZFY8jORQK7p+7aI1l4eA+qDLw0tTBwQBfAL/WZ5mBhP8Y/1GLBGGT6gRdLhdQBpRxST9UEU7gSKgCCgCPkVAXRw+bRgVa+QIwA9t3Bn2/shL1hIUgdogoBZ0bXDWWhQBRUARKBsBtaDLhkwzKAKKgCJQGwSUoGuDs9aiCCgCikDZCChBlw2ZZlAEFAFFoDYIKEHXBmetRRFQBBSBshFQgi4bMs2gCCgCikBtEPj/Af5p1gg8tUwiAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Three\n",
    "\n",
    "The Fibonacci Sequence is the sequence starting with $f_0 = f_1=1$, and where each term $f_k$ for $k>1$\n",
    "is the sum of the two preceding terms:  $f_k = f_{k-2} + f_{k-1}$. \n",
    "\n",
    "\n",
    "Here is a plot of the first 10 Fibonacci numbers:\n",
    "\n",
    "\n",
    "![Screen%20Shot%202021-08-07%20at%2010.27.41%20AM.png](attachment:Screen%20Shot%202021-08-07%20at%2010.27.41%20AM.png)\n",
    "\n",
    "Also on the plot is an exponential curve of the type $cr^k$\n",
    "for certain specific $c$ and $r.$ Notice how the curve fits the Fibonacci numbers better as $k$ gets larger. The\n",
    "curve suggests that the growth rate of $f_k$ is close to an exponential, for large $k.$\n",
    "\n",
    "In this problem we will determine the asymptotic growth rate of $f_k$, i.e., we will find $r$ \n",
    "such that as $k$ approaches infinity, $f_k$ approaches $cr_k$. Note that for asymptotic analysis, we do\n",
    "not need to know the specific $c$ involved, just the $r$. However, once we have $r$, it is easy to find an estimate\n",
    "of $c$ from the last data points. \n",
    "\n",
    "This system can be written as a linear dynamical system, where the state vector can be written as\n",
    "\n",
    "$${\\bf x}_k = \\begin{bmatrix} f_{k+1} \\\\ f_{k} \\end{bmatrix}$$ \n",
    "\n",
    "Therefore, our initial state is ${\\bf x}_0 =\n",
    "\\begin{bmatrix} 1 \\\\ 1 \\end{bmatrix}.$\n",
    "\n",
    "(A) Give the matrix $A$ such that $A{\\bf x}_k = {\\bf x}_{k+1}$, for $k\\ge 0$; thus, the $k^{th}$ Fibonacci\n",
    "number is found in the second row of $A^k x_0 = x_k$.\n",
    "\n",
    "Calculate `A @ x0` and `A @ A @ x0` to verify that your array is correct. \n",
    "\n",
    "(B) Create a list `fibs` (using $A$ and $x_0$) of the first 31 Fibonacci numbers (from $f_0$ to $f_{30}$), and then and print them out in the following form:\n",
    "\n",
    "    f0 : 1\n",
    "    f1 : 1\n",
    "    f2 : 2\n",
    "    f3 : 3\n",
    "    f4 : 5\n",
    "    f5 : 8\n",
    "    f6 : 13\n",
    "    ...etc...\n",
    "\n",
    "(C) Find the characteristic equation for $A$ and find the eigenvalues by hand; both will be real, and the larger of the two is the growth rate $r$ that we seek.  Determine $r$, showing all work. \n",
    "\n",
    "\n",
    "(D)  Determine the constant $c$ by fitting the last data points, i.e., by solving:\n",
    "\n",
    "$$ f_{30} \\approx  c\\cdot r^{30}.$$\n",
    "\n",
    "(E) Plug your values into the code below to recreate the graph above. \n",
    "\n",
    "Note:  The value $r$, the larger of the two eigenvalues, is in fact the *ratio* between successive Fibonacci numbers, and approaches the <a href='https://en.wikipedia.org/wiki/Golden_ratio'>Golden Ratio</a>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution (A)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution (B)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution (C):**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution D\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# (E)  --  verify your results for D\n",
    "\n",
    "\n",
    "fs = fibs[:10]\n",
    "xs = list(range(len(fs)))\n",
    "\n",
    "rs = [ c*(r ** x) for x in xs]\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.title('Graph of Fibonacci Numbers vs  $cr^{k}$')\n",
    "plt.grid(color='r',alpha=0.1)       # alpha sets the transparency, 0 = invisible and 1 = normal           \n",
    "plt.plot(xs,rs,color='r',lw=0.5,label='$cr^{k}$')\n",
    "plt.scatter(xs,fs,color='b',marker='o',label='$f_k$')\n",
    "plt.text(0,45,\"c = \"+str(np.round(c,6)))\n",
    "plt.text(0,40,\"r = \"+str(np.round(r,6)))\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"k\",fontsize=14)\n",
    "plt.ylabel(\"$f_k$\",rotation=0,fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Four\n",
    "\n",
    "In this problem we will consider the Fibonacci sequence from another point of view,\n",
    "using a diagonalization of the matrix A to do the calculation in the eigenspace. \n",
    "\n",
    "(A) For the array $A$ you used in Problem Three, calculate the eigenvalues and eignevectors using\n",
    "the function `np.linalg.eig(...)`, which will produce normalized eigenvectors.  \n",
    "\n",
    "(B) Define the matrix $B$ as the eigenvectors calculated in Part A, and calculate the inverse using `np.linalg.inv(...)`, and calculate the diagonal matrix $L$ which has the eigenvalues on the diagonal (a nice way to do this is with the function `np.diag(...)`). Print these out, and then confirm that $A$ and $B L B^{-1}$\n",
    "are the same matrix, and that $L$ and $B^{-1} A B$ are the same matrix. \n",
    "\n",
    "(C) Now create a matrix $L30$ which has the values $\\lambda_1^{30}$ and $\\lambda_2^{30}$ on the diagonal, and\n",
    "confirm that $B^{-1} L B x_0$ indeed has $f_{30}$ on its second row. In addition, show that the ratio\n",
    "of the two rows in your result (which are $f_{31}$ and $f_{30}$) is approximately equal to $r$ calculated above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution (A)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution (B)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution (C)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Five\n",
    "\n",
    "In this problem we will define a number of functions which we have studied in the lectures on analytic geometry.\n",
    "Fill in the templates as instructed to complete the code and satisfy the tests. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "# Part (A)\n",
    "\n",
    "# Implement the dot product using @ or sum and * on 1D arrays\n",
    "# You may need to use .T as well!\n",
    "# Be sure to return a float, not an array!\n",
    "\n",
    "\n",
    "def dot_v(u,v):\n",
    "    return 0               # just to get it to compile, substitute your code!\n",
    "\n",
    "u = np.array([[1],[2]])\n",
    "v = np.array([[-1],[3]])\n",
    "w = np.array([[9],[3]])\n",
    "\n",
    "print( dot_v(u,v) )    # should be 5\n",
    "print( dot_v(u,w) )    # should be 15\n",
    "print( dot_v(v,w) )    # should be 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "# Part (B)\n",
    "\n",
    "# Calculate the length of v using dot_v\n",
    "\n",
    "def length(u):\n",
    "    return 0               # just to get it to compile, substitute your code!\n",
    "\n",
    "print( np.round(length(u),4) )    # should be 2.2361\n",
    "print( np.round(length(v),4) )    # should be 3.1623\n",
    "print( np.round(length(w),4) )    # should be 9.4868"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 \n",
      " length =  0\n",
      "\n",
      " 0 \n",
      " length =  0\n"
     ]
    }
   ],
   "source": [
    "# Part (C)\n",
    "\n",
    "# Normalize a vector to make it of unit length\n",
    "\n",
    "def normalize(v):\n",
    "    return 0               # just to get it to compile, substitute your code!\n",
    "\n",
    "u1 = normalize(u)                         # should be [[0.4472136],[0.89442719]]                         \n",
    "print(u1,'\\n','length = ',length(u1))     # and length ~ 0.0             \n",
    "v1 = normalize(v)\n",
    "print('\\n',v1,'\\n','length = ',length(v1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "# Part (D)\n",
    "\n",
    "#  Calculate the angle between two vectors\n",
    "#  Don't forget about np.arccos(...)\n",
    "\n",
    "def angle(u,v):\n",
    "    return 0               # just to get it to compile, substitute your code!\n",
    "\n",
    "\n",
    "print( np.round(angle(u,v),4) )    # should be 0.7854\n",
    "print( np.round(angle(u,10*u),4) ) # should be 0.0\n",
    "print( np.round(angle(w,v),4) )    # should be 1.5708 = pi/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "False\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "# Part (E)\n",
    "\n",
    "# Test if u and v are orthogonal and return True or False; use np.isclose\n",
    "\n",
    "def orthogonal(u,v):\n",
    "    return False               # just to get it to compile, substitute your code!\n",
    "\n",
    "print(orthogonal(u,v))      # should be False\n",
    "print(orthogonal(w,v))      # should be True\n",
    "print(orthogonal(u,w))      # should be False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 0 0]\n",
      " [0 3 0]\n",
      " [0 0 1]]\n",
      "True \n",
      "\n",
      "[[ 3.  -1.  -0.5]\n",
      " [ 1.   2.  -2. ]\n",
      " [ 1.   1.   3.5]]\n",
      "True \n",
      "\n",
      "[[2 2 0]\n",
      " [0 3 1]\n",
      " [0 1 1]]\n",
      "True \n",
      "\n",
      "[[ 0.12309149  0.10071122 -0.36927447]\n",
      " [ 0.49236596  0.03357041  0.73854895]\n",
      " [ 0.86164044 -0.03357041 -0.36927447]]\n",
      "True \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Part (F)\n",
    "\n",
    "# Test if all the columns of B are pairwise orthogonal\n",
    "# Test each against all the *other* columns, since a vector is never orthogonal to itself\n",
    "\n",
    "\n",
    "def orthogonal_set(B):\n",
    "    # Your code here\n",
    "    return True\n",
    "\n",
    "B = np.array( [[2,0,0],[0,3,0],[0,0,1]] )\n",
    "print(B)\n",
    "print(orthogonal_set(B),'\\n')              # should be True\n",
    "\n",
    "B1 = np.array( [[3,-1,-1/2],[1,2,-2],[1,1,7/2]] )\n",
    "print(B1)\n",
    "print(orthogonal_set(B1),'\\n')              # should be True\n",
    "\n",
    "B2 = np.array( [[2,2,0],[0,3,1],[0,1,1]] )\n",
    "print(B2)\n",
    "print(orthogonal_set(B2),'\\n')              # should be False\n",
    "\n",
    "B3 = np.array( [[ 0.12309149,  0.10071122, -0.36927447],\n",
    " [ 0.49236596,  0.03357041,  0.73854895],\n",
    " [ 0.86164044, -0.03357041, -0.36927447]])\n",
    "print(B3)\n",
    "print(orthogonal_set(B3),'\\n')              # should be True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Six\n",
    "\n",
    "Now you will write a function which will calculate the orthogonal decomposition of\n",
    "a vector $u$ with respect to a subspace with a basis $B$. \n",
    "\n",
    "Formula and discussion are from lecture and the online textbook.  Complete the template\n",
    "and verify that all tests are satisfied. \n",
    "\n",
    "Solution is shown at the bottom of the notebook, your results should be approximately the same. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0 \n",
      "\n",
      "0\n",
      "0 \n",
      "\n",
      "0 \n",
      "\n",
      "0\n",
      "0 \n",
      "\n",
      "0 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# u is a vector and B is a basis for a subspace\n",
    "# return u_B^perp and u_B = projection of u onto B\n",
    "\n",
    "\n",
    "def orthogonal_decomposition(u,B):\n",
    "    u_proj = 0      # just to get it to compile, your code here\n",
    "    # Your code here      \n",
    "    u_perp = 0      # just to get it to compile, your code here\n",
    "    return (u_perp,u_proj)\n",
    "\n",
    "u = np.array( [[7],[6]])\n",
    "B = np.array( [[4],[2]] )      # projection onto a line\n",
    "\n",
    "OD1 = orthogonal_decomposition(u,B)\n",
    "print(OD1[0])\n",
    "print(OD1[1])\n",
    "print(OD1[0]+OD1[1],'\\n')      # should print out u\n",
    "\n",
    "u = np.array( [[1],[1],[1]])\n",
    "B = np.array( [[1,0],[0,1],[0,0]] )           # xy-plane in R^3\n",
    "\n",
    "OD2 = orthogonal_decomposition(u,B)\n",
    "print(OD2[0])\n",
    "print(OD2[1],'\\n')\n",
    "print(OD2[0]+OD2[1],'\\n')      # should print out u\n",
    "\n",
    "u = np.array( [[2],[-4],[1]])\n",
    "B = np.array( [[2,-1],[-1,1/2],[3,-2]] )      # some plane in R^3\n",
    "\n",
    "OD3 = orthogonal_decomposition(u,B)\n",
    "print(OD3[0])\n",
    "print(OD3[1],'\\n')\n",
    "print(OD3[0]+OD3[1],'\\n')      # should print out u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Seven\n",
    "\n",
    "Now you will write a function which will calculate a orthogonal or an orthonormal basis\n",
    "from a set of vectors.  It should work on *any* set of vectors, as suggested in the online\n",
    "textbook. \n",
    "\n",
    "Formula and discussion are from lecture and the online textbook.  Complete the template\n",
    "and verify that all tests are satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n",
      "[[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]] \n",
      "\n",
      "[[1 2 3]\n",
      " [4 5 6]\n",
      " [7 8 9]]\n",
      "[[1 2 3]\n",
      " [4 5 6]\n",
      " [7 8 9]] \n",
      "\n",
      "[[3 1]\n",
      " [6 2]\n",
      " [0 2]]\n",
      "[[3 1]\n",
      " [6 2]\n",
      " [0 2]] \n",
      "\n",
      "[[1 0 0]\n",
      " [1 1 0]\n",
      " [1 1 1]\n",
      " [1 1 1]]\n",
      "[[1 0 0]\n",
      " [1 1 0]\n",
      " [1 1 1]\n",
      " [1 1 1]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def GramSchmidt(B,normal=False):\n",
    "    \n",
    "    # Your code here\n",
    "    \n",
    "    return B     # just to get it compile, your code here\n",
    "    \n",
    "    \n",
    "B = np.array( [[1,0,0],[0,1,0],[0,0,1]] ) \n",
    "print(B)\n",
    "OB = GramSchmidt(B) \n",
    "print(np.around(OB,10),'\\n')\n",
    "\n",
    "B = np.array( [[1,2,3],[4,5,6],[7,8,9]] ) \n",
    "print(B)\n",
    "OB = GramSchmidt(B) \n",
    "print(np.around(OB,10),'\\n')\n",
    "\n",
    "B = np.array( [[3,1],[6,2],[0,2]] ) \n",
    "print(B)\n",
    "OB = GramSchmidt(B) \n",
    "print(np.around(OB,10),'\\n')\n",
    "\n",
    "B = np.array( [[1,0,0],[1,1,0],[1,1,1],[1,1,1]] ) \n",
    "print(B)\n",
    "OB = GramSchmidt(B) \n",
    "print(np.around(OB,10),'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution for problem 6 should be:\n",
    "\n",
    "    [[-1.]\n",
    "     [ 2.]]\n",
    "    [[8.]\n",
    "     [4.]]\n",
    "    [[7.]\n",
    "     [6.]] \n",
    "\n",
    "    [[0.]\n",
    "     [0.]\n",
    "     [1.]]\n",
    "    [[1.]\n",
    "     [1.]\n",
    "     [0.]] \n",
    "\n",
    "    [[1.]\n",
    "     [1.]\n",
    "     [1.]] \n",
    "\n",
    "    [[-0.71428571]\n",
    "     [-2.64285714]\n",
    "     [-3.64285714]]\n",
    "    [[ 2.71428571]\n",
    "     [-1.35714286]\n",
    "     [ 4.64285714]] \n",
    "\n",
    "    [[ 2.]\n",
    "     [-4.]\n",
    "     [ 1.]]\n",
    "\n",
    "\n",
    "Solution for problem 7 should be:\n",
    " \n",
    " \n",
    "    [[1 0 0]\n",
    "     [0 1 0]\n",
    "     [0 0 1]]\n",
    "    [[1]\n",
    "     [0]\n",
    "     [0]] \n",
    "\n",
    "    [[1 2 3]\n",
    "     [4 5 6]\n",
    "     [7 8 9]]\n",
    "    [[1]\n",
    "     [4]\n",
    "     [7]] \n",
    "\n",
    "    [[3 1]\n",
    "     [6 2]\n",
    "     [0 2]]\n",
    "    [[3]\n",
    "     [6]\n",
    "     [0]] \n",
    "\n",
    "    [[1 0 0]\n",
    "     [1 1 0]\n",
    "     [1 1 1]\n",
    "     [1 1 1]]\n",
    "    [[1]\n",
    "     [1]\n",
    "     [1]\n",
    "     [1]]"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}