Newton's Divided Difference Polynomial Interpolation.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"Newton's Divided Difference Polynomial Interpolation.ipynb","provenance":[],"collapsed_sections":[]},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"markdown","metadata":{"id":"czoA0EnWua5v"},"source":["## CSE330 Lab: Polynomial Interpolation using Newton's Divided Difference form\n","---"]},{"cell_type":"markdown","metadata":{"id":"8dQJpP15jzP9"},"source":["###1. Importing necessary libraries\n","\n"]},{"cell_type":"code","metadata":{"id":"iLsFZK0JjoYk","executionInfo":{"status":"ok","timestamp":1601827901905,"user_tz":-360,"elapsed":1211,"user":{"displayName":"SUHAIL HAQUE RAFI","photoUrl":"","userId":"10412558163487613106"}}},"source":["import numpy as np\n","import matplotlib.pyplot as plt"],"execution_count":1,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"7RrXSL8lj3YR"},"source":["###2. Writing the `Newtons_Divided_Differences` class.\n","\n","Newton form of a $n$ degree polynomial:\n","\n","$$p_n(x) = \\sum_{k=0}^{n} a_kn_k(x)$$\n","where the basis,\n","$$n_k(x) = \\prod_{j=0}^{k-1}(x-x_j)$$\n","$$ n_0(x)=1$$\n","\n","and the coefficients, $$a_k = f[x_0, x_1, ..., x_k]$$\n","\n","where the notation $f[x_0, x_1,...,x_k]$ denotes the divided difference.\n","\n","\n","Expanding the Newton form:\n","\n","$$p(x) = f [x_0] + (x-x_0) f[x_0,x_1] + (x-x_0) (x-x_1) f[x_0,x_1,x_2] + . . . + (x-x_0) (x-x_1) . . . (x-x_{k-1}) f[x_0, x_1, . . ., x_k]$$\n"]},{"cell_type":"markdown","metadata":{"id":"qDM6A5gZo7_T"},"source":["\n","### Task: \n","1. The `calc_div_diff(x,y)` function which takes as input `x` and `y` and calculates all the divided differences.\n","\n","2. The `__call__` function which takes as input `x` and calculates `y` using all the difference coefficients. `x` can be a single value or a numpy. In this case it is a numpy array.\n"]},{"cell_type":"code","metadata":{"id":"H6hA0gsKjUYq","executionInfo":{"status":"ok","timestamp":1601827901909,"user_tz":-360,"elapsed":1203,"user":{"displayName":"SUHAIL HAQUE RAFI","photoUrl":"","userId":"10412558163487613106"}}},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","class Newtons_Divided_Differences:\n","  def __init__(self, differences):\n","    self.differences = differences\n","    \n","  def __call__(self, x):\n","    \"\"\"\n","    this function is for calculating y from given x using all the difference coefficients\n","    x can be a single value or a numpy\n","    the formula being used:\n","    f(x) = f [x0] + (x-x0) f[x0,x1] + (x-x0) (x-x1) f[x0,x1,x2] + . . . + (x-x0) (x-x1) . . . (x-xk-1) f[x0, x1, . . ., xk]\n","\n","    \"\"\"\n","\n","    n=len(data_x) - 1 \n","    res =self.differences[n]\n","    for k in range(1,n+1):\n","        res=self.differences[n-k]+(x-data_x[n-k])*res\n","\n","    return res"],"execution_count":2,"outputs":[]},{"cell_type":"code","metadata":{"id":"NTND3V4xxkdM","executionInfo":{"status":"ok","timestamp":1601827901912,"user_tz":-360,"elapsed":1198,"user":{"displayName":"SUHAIL HAQUE RAFI","photoUrl":"","userId":"10412558163487613106"}}},"source":["def calc_div_diff(x,y):\n","  assert(len(x)==len(y))\n","\n","  lng = len(x)\n","\n","  x=np.copy(x) \n","  b=np.copy(y)\n","\n","  for k in range(1,lng):\n","    b[k:lng]=(b[k:lng]-b[k-1])/(x[k:lng]-x[k-1])  \n","  \n","  return b "],"execution_count":3,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"GYpfaqvRuSS8"},"source":["###3. Creating and plotting the polynomial"]},{"cell_type":"markdown","metadata":{"id":"uH7EhC-fuTQD"},"source":["* `data_x` and `data_y` are the coordinates of the given nodes.\n","\n","* `differences` is a list which contains the divided differences as each of its elements: $f[x_0], f[x_0,x_1], f[x_0,x_1,x_2], ...$\n","\n","* `obj` is an object of type `Newtons_Divided_Differences`. Creating the object runs the constructor of the class where the `difference` are stored in `self.differences`.\n","\n","* `X` contains $x_i$ values through which we want to plot our polynomial.\n","\n","* Calling the object using `obj(X)` executes the `__call__` function of the class which returns a numpy array containing the corresponding $y_i$ values, and storing them in variable `F`.\n","\n","* Using `plt.plot(X,F)` we plot the $(x_i, y_i)$ pairs of the polynomial."]},{"cell_type":"code","metadata":{"id":"YCICBXzRjrNB","executionInfo":{"status":"ok","timestamp":1601827902662,"user_tz":-360,"elapsed":1921,"user":{"displayName":"SUHAIL HAQUE RAFI","photoUrl":"","userId":"10412558163487613106"}},"outputId":"3b509e1c-47bb-4d24-dbbb-447c40fd7054","colab":{"base_uri":"https://localhost:8080/","height":265}},"source":["data_x = np.array([-3.,-2.,-1.,0.,1.,3.,4.])\n","data_y = np.array([-60.,-80.,6.,1.,45.,30.,16.])\n","differences = calc_div_diff(list(data_x), list(data_y))\n","obj = Newtons_Divided_Differences(list(differences))\n","\n","#generating 50 points from -3 to 4 in order to create a smooth line\n","X = np.linspace(-3, 4, 50, endpoint=True)\n","F = obj(X)\n","plt.plot(X,F)\n","plt.plot(data_x, data_y, 'ro')\n","plt.show()"],"execution_count":4,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3yV5f3/8dcnmwxCICGbDEIiK0SJQcCB2wKC4jatu9aqrf22v29rq11aWlutba3jK1pbK1G0VgVEUcCNMiUEElYY2SRhhCSEzHP9/shBowJJyDnnPuPzfDzyIOc+437Xwjt3rvu6r1uMMSillPItflYHUEop5Xpa/kop5YO0/JVSygdp+SullA/S8ldKKR8UYHWAvoiOjjapqalWx1BKKY+yfv36fcaYmGM95xHln5qayrp166yOoZRSHkVEyo73nA77KKWUD9LyV0opH6Tlr5RSPkjLXymlfJCWv1JK+SAtf6WU8kFa/kop5YM8Yp6/Usp97Wtuo6S6kZb2Llo7ujjS0fXF92FB/lx+WhKRgwKtjqm+RstfKXVS6ppaefrDXcxfVUZbp+24r3vk3e18Z3IKt0xNIyYi2IUJ1Ylo+Sul+qW+qY2nP9zJ/NVldHQZLj81kSsnJhEREkBoUACDAv27v4L8Ka1r5okPSvm/D3fy3Ce7uS5vBN89O53EIYOs/p/h88QT7uSVm5trdHkHpay1v7mNpz/axb8/20N7p43LT03iB+dlkBod1ut7d9U389QHO3l9QxUAV5yWxP0zRxMRosNBziQi640xucd8TstfKdWbDeUHufX5dTS0tHNZTiI/OH8UaX0o/a+rajjCvA93UrC6nJzkITx/Sx5hwToA4SwnKn+d7aOUOqEVW2q57plVRIQEsPRHZ/PoNTknVfwAiUMG8dvZ43jsulPZUNHArc+v5Uh7l4MTq77Q8ldKHdeCNeV899/ryIyN4L/fn0JmbIRDPnf6+HgevXoCq3cf4PYX1tHaoT8AXE3LXyn1DcYY/rp8O/e+tomzM2N46btnEB3u2Jk6s3MS+eMV2Xy8Yx93FXxO+wlmDCnHc0j5i8hzIlInIpt7bBsqIstEZIf9zyj7dhGRx0SkVESKROQ0R2RQSjlGZ5eNn7+2ib8u38GVE5N45oZcp43LX52bzO8uG8eKrXX88KUNdHbpDwBXcdSR/7+AS7627V5ghTFmFLDC/hjgW8Ao+9ftwFMOyqCUGqD2Tht3zF/PgrUV3H1uBg9fmU2gv3MHCL59Rgq/mjmGpcV7+fErG+myuf8kFG/gkB/nxpiPRCT1a5tnA9Ps3z8PfAD8zL7936Z7mtEqERkiIvHGmBpHZFFKnRxjDPe9vonlW+p4YPZYbpic6rJ933JmGu1dNh56eyupw0L58UVZLtu3r3Lmj/TYHoW+F4i1f58IVPR4XaV921eIyO0isk5E1tXX1zsxplIK4NmPd/Of9ZX88LwMlxb/UXecM5I5pyby1Ic72VHb5PL9+xqXnPC1H+X363c5Y8w8Y0yuMSY3JuaY9x9WSjnIii21/P7tLUwfH8ePLsi0LMd9M0YTHhzAz1/bhE2Hf5zKmeVfKyLxAPY/6+zbq4DkHq9Lsm9TSllg294mfvjSBsYmDObPV+Xg5yeWZRkWHswvpo9mXdlBFqyt6P0N6qQ5s/wXATfav78RWNhj+w32WT9nAId0vF8pa+xvbuPW59cSFhzAMzfkMijI3+pIXDkxicnpw/jD21uoa2q1Oo7XctRUz5eAz4AsEakUkVuBh4ALRWQHcIH9McBbwC6gFHgGuNMRGZRS/dPW2cUd89dT39TGMzfkEh/pHoutiQhzLx9HW6eNBxaXWB3Hazlqts91x3nq/GO81gB3OWK/SqmT0z2zZzNr9xzk79edyoTkIVZH+or0mHDuPjeDR5dt54rT6jj3lOFWR/I6eoWvUj5o/upyXl1fyT3nj+LSCQlWxzmmO84Zyajh4dz/xmZa2jutjuN1tPyV8jGldU387s0Szs6M4Z7zR1kd57iCAvz4/ZzxVDUc4S/Ltlsdx+to+SvlQ9o7bdyzoJCw4AAeuTLb0pk9fXF66lCuyxvBcyv3sLnqkNVxvIqWv1I+5NFl2ymubuShOeMZPjjE6jh9cu8lpxAVGsTvlujJX0fS8lfKR3y2cz9Pf7ST6/KSuWhsnNVx+iwyNJA7zkln1a4DrC87aHUcr6Hlr5QPONTSwU9eKSR1WBi/nDnG6jj9dl3eCIaEBvLUB6VWR/EaWv5KeTljDPcv3ExdUxt/vSaH0CDPu21iWHAAN01JZfmWOrbt1XV/HEHLXykv90ZhFYs3VvOjC0a53Xz+/rhpSiqhQf569O8gWv5KebGKAy386o1iTk+N4vvTMqyOMyBDQoPInzSCxUU1VBxosTqOx9PyV8pLGWP46atFGODRq3Pwd/NpnX1x21np+Ivw9Ec7rY7i8bT8lfJSL6+t4LNd+/nF9NEkDw21Oo5DxA4O4YqJibyyrlIXfRsgLX+lvNDeQ63MXbKFyenDuC4vufc3eJDvnT2Szi4bz32yx+ooHk3LXykvY4zh/jc202Gz8Yc54xHx/OGenlKjw5g+Pp75q8o4dKTD6jgeS8tfKS/zZlENy7fU8pMLs0iNDrM6jlN8f9pImts6eeGzPVZH8Vha/kp5kQOH2/nNomImJEVy89RUq+M4zdiESKZlxfDcyj0cae+yOo5H0vJXyos8+GYJja0d/PHKbAL8vfuf953TMjhwuJ2X15ZbHcUjefffDqV8yPtb63h9QxV3TsvglLjBVsdxury0oeSmRPHMx7vp0pu995uWv1JeoKm1g/te30RmbDh3njvS6jguc9PUVKoajrCydJ/VUTyOlr9SXuBPS7dR09jKH6/IJjjA+puwu8qFY2KJCg3k5bUVVkfxOFr+Snm49WUHeGFVGTdPSePUEVFWx3Gp4AB/5pyWxLsle9nf3GZ1HI+i5a+UB2vvtPHz1zaROGQQP7ko0+o4lrjm9GQ6ugyvb6iyOopH0fJXyoPN+2gn22ubefCysYQFe95SzY6QGRvBaSOGsGBtBcboid++0vJXykPtqm/msfdKmZEdz3mnxFodx1LXnj6C0rpmPi/XO331lZa/Uh7IGMMvXt9EcIAfv77U8+7M5WgzsuMJC/JnwRo98dtXWv5KeaD/rK9k1a4D/GL6aIZHeMaN2J0pLDiAWTkJvFlUQ1OrrvfTF1r+SnmYfc1tzF2yhbzUoVyT610rdg7ENaeP4EhHF4s31lgdxSNo+SvlYX73Zgkt7Z38fs44/LzgBi2OMiEpklPiIligyz30iZa/Uh7kw+31vFFYzZ3TMsgYHmF1HLciIlxzejJFlYcorj5kdRy3p+WvlIdoae/kvtc3kR4T5lNLOPTH5acmEhTgxyt6xW+vtPyV8hCPvLOdyoNHeGiOby3h0B9DQoP41rg4Xt9QRWuHLvV8Ilr+SnmAz8sP8s9Pd/OdM1LISxtqdRy3ds3pyTS2drJ0816ro7g1p5e/iOwRkU0iUigi6+zbhorIMhHZYf/TtxYkUaof2jq7+NmrRcQPDuGnl2RZHcftnZE2jJRhoXritxeuOvI/1xiTY4zJtT++F1hhjBkFrLA/VkodwxPv72RHXTNz54wnIiTQ6jhuz89PuDo3mVW7DlC+v8XqOG7LqmGf2cDz9u+fBy6zKIdSbm1LTSNPvl/KnFMTOTdruNVxPMbsnAQAFhdVW5zEfbmi/A3wroisF5Hb7dtijTFHr8TYC3xjYRIRuV1E1onIuvr6ehfEVMq9dHbZ+Nl/i4gcFMgvZ+oSDv2RFBVKbkoUiwq1/I/HFeV/pjHmNOBbwF0icnbPJ033MnzfWIrPGDPPGJNrjMmNiYlxQUyl3ERBAaSm4h8YwJP3zeGZgG1EhQVZncrjzMpJYFttE9v2NlkdxS05vfyNMVX2P+uA14E8oFZE4gHsf9Y5O4dSHqGgAG6/HcrKEGNIaqzn1N/9rHu76pfp4+Px9xMWbdR1/o/FqeUvImEiEnH0e+AiYDOwCLjR/rIbgYXOzKGUx7jvPmj56klKaWnp3q76JTo8mCkjh7F4Y42u838Mzj7yjwU+EZGNwBpgiTFmKfAQcKGI7AAusD9WSpUfZ3ri8barE5o1IYHyAy0UVjRYHcXtOLX8jTG7jDET7F9jjTFz7dv3G2PON8aMMsZcYIw54MwcSnmKjsSkYz8xYoRrg3iJi8fFERTgx0I98fsNeoWvUm6ivdPG3869iSOBwV99IjQU5s61JpSHGxwSyLlZMSzZVEOXTYd+etLyV8pN/P29HTyeMIntDz4KKSkg0v3nvHmQn291PI81a0Ii9U1trNq13+oobsU37/islJtZX3aQJ94v5aqJSUy4agb87E6rI3mN80cPJyzIn0WF1UzNiLY6jtvQI3+lLHa4rZMfv1JIwpBB/Ervx+twIYH+XDw2jrc319DWqSt9HqXlr5TFfrekhPIDLTx6dY6u3eMkl+Yk0NjayUfb91kdxW1o+StloWUltby0poI7zhmpSzU70ZkZ0USFBrKwUC/4OkrLXymL1De1ce9/ixgTP5j/uSDT6jheLdDfj+nj41m+pZbDbZ1Wx3ELWv5KWaCzy8YPX9pAc1snf702h6AA/afobLMmJHBR4QokLRX8/CA11aeXzdDZPkpZ4E/vbOOzXft55KoJZMbqjdhd4fRPl5L9zuMM6mjr3lBW1r2OEvjkVFo93FDKxd4sqmbeR7u4YXIKV048zhW9yuH87r/vy+I/yofXTdLyV8qFttc28dNXi5iYEsX9M3Rap0vpuklfoeWvlIs0tnbwvRfWExoUwJP5p+k4v6sdb30kH103Sf/2KeUCNpvhxy9vpOJAC0/mn0bs4BCrI/meuXO710nqyYfXTdLyV8oFnvyglOVbarlvxmidz2+V/HyYN4/2xGRsCIfjEn163SSd7aOUk72/tY4/L9vOZTkJ3DQl1eo4vi0/n8Drr+fcRz4geWgoL+RPsjqRZfTIXyknWrfnAHcWfM7ouMH8YU42ImJ1JJ8nIszIjufTnfvZ39zW+xu8lJa/Uk6yueoQN/9rLfGRITx/Sx6DgvytjqTspo+Pp8tmeKe41uooltHyV8oJSuuaufG5NUQEB/DCbZOIiQju/U3KZcbEDyYtOowlm3z3Dl9a/ko5WOXBFr7zj9WIwPzbJpE4ZJDVkdTXiAgzxsfzmQ8P/Wj5K+VAdU2tfPvZ1Rxu6+Tft0wiPSbc6kjqOKaPj8dmYGnxXqujHFtBQff6Q05ah0jLXykHaWhp54Z/rKG2sY1/3pzHmITBVkdSJzA6PoL06DCWFNVYHeWbCgq61x0qKwNjvlyHyIE/ALT8lXKAvYdayX92NbvqD/PMDblMTImyOpLqxdFZP6t27Wefuw393Hdf97pDPTl4HSItf6UGqLCigVmPf8KefYd5+oaJnDlK7xPrKWZk24d+NrvZ0I8L1iHS8ldqABYWVnHN058RFODHa3dO5dys4VZHUv2QFRtBeowbDv24YB0iLX+lToLNZnjknW3cs6CQCclDWHjXVLLidF1+TyMizBwfz+rd+6lvcqOhn7lz6Qr52iwxB69DpOWvVD8dbuvk+wXrefz9Uq7JTWb+rZMYFq7z+D3VjOwE95v1k5/P8v95kMrBMRgRSElx+DpEWv5K9cP6sgPMefJTlpXU8suZY3joivG6NLOHy4wNZ2RMGEuK3OuCrzfGTOPb97+M2GywZ4/DF6DThd3UySko6J55UF7ePQ45dy41M+fweVkDG8oP0tjagSCIQPdyNt1r2gwJDWRkTDgZw7v/wUWEBFr6P6Ovahtbeejtrby+oYq4wSE8d9PpTNPxfa/QPesngcff20FdUyvDI9xjue2SmkbGOnG6sJa/6r+jc5CPTkUrK+PIzbfyh4s3sGjsuQQH+BEVGoTBYAwYuqcqg+HQkQ46uswXHxU3OISM4eFkxUWQmxLF6WlDiXajIZS2zi6e+2QPf39vB51dhrvOHcmd0zIIC9Z/Ot5kZnY8j63YwdLNe7lhcqrVcWhs7aBsfwtX5yY7bR/6N1j1W/vP7iXoa3OQB3W08fu1L3Hrk/czOn7wcYdCOrpslB9oobSumdK6ZnbWNVNa38z8VWX845PdAKTHhJGXOpS8tKGcnjqUpKhBLl8N0xjDii11/G5JCXv2t3DhmFjunzGalGFhLs2hXCMzNoLM2HAWb6x2i/LfWtMEdK9B5Cxa/qrPdtQ28ed3t/NkVdUxnw+vq2FC8pATfkagvx8jY8IZGRPOxWO/3N7eaWNT1SHW7jnA2t0HeGtTDQvWVgCQEBlCXtpQJqUPIy9tKOnRYU75YWCMYUtNE4s2VrN4YzVVDUdIjwnj+VvyOCczxuH7U+5ldk4iD7+zjcqDLSRFhfb+BicqqT4E4NSrxC0rfxG5BPgb4A88a4x5yKos6sQqDrTwl+XbeWNDFaFBATQPj2dw3TFOjg1gDnJQgB8TU6KYmBLFHeeMxGYzbKttYs3uA6zZc4BPSvfzRmH3PqPDg8hLG8q4xEhOiYsgK24wCZEhx/+BcIzzEz1Pnu3ed5jFG6tZtLGa0rpm/P2EMzOi+clFmczMTtATuj5i1oQEHn5nG4s2VnPntAxLsxRXNxIdHsRwJ64Ga0n5i4g/8ARwIVAJrBWRRcaYEivyqON7YVUZDywuxk+E285K545zRjI4809fHfMHx89B9hNGxw9mdPxgbpySijGG3fsOd/8wsP9AeGvTl1PzIkICyIqNIDMugoTIEEKDAggL9idj+SJyHvgp/q1Hul9YVkb7rbcxf+VuXht9DuX7W2hs7QQgL3UoD142junj4nTqpg9KHhrKxJQoFhVaX/4lNY2Mjh/s1OFOq47884BSY8wuABFZAMwGtPzdhM1m+OPSrTz90S7OO2U4v798PHGR9lkQR4+aT3A07WgiQnpMOOkx4Vyb1/0bRmNrB9v3NrF1bxPb7F9vbqz+oswBPnnqgS+L3y6orZXpLz3Oh3++iJzkIWTEhHPR2DgSdOllnzc7J4FfLSxm695GTomzZmG+9k4b22ubuOXMNKfux6ryTwQqejyuBL5yM00RuR24HWCEAy9pVr1r7ejiJ//ZyJKiGr5zRgq/mTUWf7+vHYHk51t+4+vBIYHkpg4lN/WrN0Rv77RxpL2Lw+2dxP9p3zHfG3eojudvyXNFTOVBpo+P57eLS1hYWM0pl1hT/qV1zXR0GcYmRDp1P247mGmMmWeMyTXG5MbE6Mk2Vzl4uJ1vP7uaJUU1/GL6KTww+xjF7+aCAvyIDA0kYcggxAVrpCjvER0ezFmjollUWI3NZnp/gxOU1DQCzp3pA9aVfxXQcwJrkn2bslDZ/sPMeepTiqoO8cT1p3H72SM9/4bjc+d2n4/oycHnJ5R3mZ2TQFXDET4vP2jJ/ourDzEo0J+0aOdOK7aq/NcCo0QkTUSCgGuBRRZlUcCmykNc/uSnNLS08+Jtk5iRHW91JMfIz+9eEyUlpftSYyeskaK8y4Vj4ggJ9GNhoTXLPZRUN3JKfITTf+O2pPyNMZ3A3cA7wBbgFWNMsRVZFNQ1tnLr82sZFOjPf78/5Rtj6B4vP797bRQnrZGivEt4cAAXjI5lyaYaOrpsLt23MYaSmkanD/mAhWP+xpi3jDGZxpiRxhj9Hdwi7Z02vl/wOU2tnfzjply956xSwGU5iRw43M4nO449YcBZKg8eoam10+kne8GNT/gq13jgzWLWlx3k4auyLZvappS7OTszhshBgSwsdO2pyOJq+8leF9z/Wcvfh728tpz5q8r53jnpzMxOsDqOUm4jKMCP6ePjebeklpb2zt7f4CAl1Yfwk+47jDmblr+PKqxo4JdvFHPWqGh+evEpVsdRyu3Mzkmgpb2LZSW1LttnSU0jI2PCGRTk7/R9afn7oLqmVu54YT2xkcH8/bpTPW4ev1KukJc6lPjIEBa5cNZPcXWjS4Z8QMvf57R32rir4HMajrTz9LdzGRIaZHUkpdySn58wa0ICH26v5+Dhdqfv78DhdmoOtTr1Bi49afn7mD8u3craPQf54xXZLjvCUMpTzcpJoNNmeGtzjdP3teWLK3udP9MHtPx9ysaKBp5buZsbJqcwOyfR6jhKub0x8YPJGB7Oq+srnb6vYhes4d+Tlr+P6LIZ7n9jMzHhwfzvxVlWx1HKI4gI+ZNGsKG8gY0VDU7dV0l1I/GRIQwNc81QrJa/j3hxTTmbqg5x/8wxHnPTdKXcwZUTkwgL8uf5T/c4dT/F1a65svcoLX8fUN/Uxp+WbmVqxjAu9ZY1e5RykYiQQK7KTWZxUTV1Ta1O2UdrRxc765tdeh5Oy98H/OHtLbR2dPHA7HGev0qnUha4YXIKHV2Gl1ZX9P7ik7BtbxM2g8tm+oCWv9dbtWs/r31exe1npzNS1+1R6qSkx4QzLSuG+avLaO90/GJvXyzr4KKZPqDl79U6umz88o3NJA4ZxN3njrI6jlIe7aYpqdQ3tfG2E6Z9ltQcIiI4gKQo191KVMvfiz33yW521DXz21ljXXK5uFLe7OxRMaRHh/HPlXsc/tnF1Y2MThiMnwuvttfy91LVDUf46/IdXDA6lgvGxFodRymP5+cn3DgllcKKBjY48C5fBw63U1zdSHai64Z8QMvfaz2wuASD4deXjrE6ilJe44qJSYQHBzh02udLa8pp77Rx9enJvb/YgbT8vdCG8oMsLd7LndMySB4a2vsblFJ9Eh4cwJUTk1iyqYa6xoFP++zosvHCZ2WcmRFNpguWce5Jy98L/W3FDqJCA7nlzDSroyjldW6ckkpHl6FgdfmAP+vtzXvZ29jKzVNTBx6sn7T8vUxhRQMfbKvntrPSCQ8OsDqOUl4nLTqMc7NiKFhdPuBpn/9auZvUYaGcmzXcQen6Tsvfy/xt+XaGhAZy45RUq6Mo5bVumprGvuY23tp08tM+N1Y08Hl5AzdOSXXpLJ+jtPy9SGFFA+9vq+e7etSvlFOdlRFNekwYz36yiy6bOanP+OfK3V+cQ7CClr8XeWzFDj3qV8oF/PyEH5yXweaqRv7vw539fn9dYytLNtVwVW6SZQstavl7iY0VDby3tU6P+pVykctyEpmZHc+jy7bzeT/n/c9fVUanzXDj5FTnhOsDLX8v8Tf7Uf8Nk1OsjqKUTxAR5l4+nvjIEH740gYaWzv69L7Wji4KVpdzXtZwUqPDnJzy+LT8vUBR5ZdH/bpWv1KuEzkokL9deyo1h1r5xWubMKb38f83i2rYf7idm6daOxVby98L/G25HvUrZZWJKVH8+MJM3iyq4T+93O7RGMM/V+4mMzacqRnDXJTw2LT8PVxRZQMrttZx25lpetSvlEXuOGckk9OH8euFxeysbz7u69buOUhxdSM3TUmz/N4aWv4e7rEVO4gcpDN8lLKSv5/wl2tyCAn04wcvbqCts+uYr/vnyt0MCQ3k8lMTXZzwm3RaiAcrqW5k+ZY6fnJhph71K2WxuMgQHr5yArf9ex2/XljM5JHDqDx4hPL9LVQctH8dOMId54x0iyXWtfw92HMrdxMa5M8NFk4XU0p96YIxsdw4OYXnPytjwdruWz5GhwczYuggThsRxVUTk7nVTdbc0vL3UPua21hUWM01pycTGapH/Uq5i19dOpZLJyQQOSiQpKhQtzjKPxanjfmLyG9EpEpECu1f03s893MRKRWRbSJysbMyeLOCVeW0d9m4yYLVAJVSx+fvJ+SmDmVUbITbFj84/8j/L8aYR3puEJExwLXAWCABWC4imcaYY58hUd/Q3mlj/uoypmXF6E3ZlVInxYrZPrOBBcaYNmPMbqAUyLMgh8dasqma+qY2yy8SUUp5LmeX/90iUiQiz4lIlH1bIlDR4zWV9m1fISK3i8g6EVlXX1/v5JiewxjDc5/sIWN4OGePirY6jlLKQw2o/EVkuYhsPsbXbOApYCSQA9QAf+7PZxtj5hljco0xuTExMQOJ6VXWlx1kU9UhbpqSavlFIkopzzWgMX9jzAV9eZ2IPAO8aX9YBfS8U3GSfZvqg3+u3MPgkADmnGb9RSJKKc/lzNk+8T0eXg5stn+/CLhWRIJFJA0YBaxxVg5vUtVwhKXFe7lu0ghCg3SWrlLq5DmzQf4kIjmAAfYA3wMwxhSLyCtACdAJ3KUzffrm35/tAdCLupRSA+a0I39jzHeMMeONMdnGmFnGmJoez801xow0xmQZY952VgYKCiA1Ffz8uv8sKHDarpytpb2Tl1aXc/HYWBKHDLI6jlLKw3nv2EFBAdx+O7S0dD8uK+t+DJCfb12uk/Ta51U0tnZyi07vVEo5gPeu6nnffV8W/1EtLd3bPYzN1r0G+PjESCamRPX+BqWU6oX3ln95ef+2u7FPSvexs/4wt5yp0zuVUo7hveU/YkT/trux+avKGBYWxPTx8b2/WCml+sB7y3/uXAgN/eq20NDu7R5k76FWVmyt46rcZIID3HeRKKWUZ/He8s/Ph3nzaIiJx4ZgRoyAefM87mTvy2sr6LIZrs/zvN9YlFLuy3tn+wDk57Ms62z+99Ui3vvJOaR72AqYnV02Fqwt56xR0YwYFtr7G5RSqo+898jfLisuAoDttU0WJ+m/D7bVU3OolfxJKVZHUUp5Ga8v/4zh4YjAtr3NVkfptxfXlDM8IpjzRw+3OopSyst4ffmHBgUwYmioxx35Vx5s4f1tdVxzejKB/l7/f5NSysV8olUyYyPY5mHl//LaCgS4Vk/0KqWcwCfKPys2gt37DtPW6Rnrx3V02ViwtoJpWcN1HR+llFP4RPlnxkXQZTPsqj9sdZQ+WbGllvqmNvIn6VG/Uso5fKL8s2I9a8ZPwepyEiJDmJalJ3qVUs7hE+WfFh1GgJ+wba/7l3/Z/sN8vGMf1+aNwN9P1/FRSjmHT5R/UIAf6TFh7n3kb7/3wIiYCFY+dTM37l5pdSKllBfz7it8e8iMjWBjZYPVMY6tx70HBEhsrId77oJBgR63HIVSyjP4xJE/dI/7Vxw4wuG2TqujfJMX3XtAKeUZfKb8M+3LPOyoc8Mrfb3o3gNKKc/gM+X/xYwfdzzp60X3HlBKeQafKf/koaGEBPq555W+c+2bRJYAAAqaSURBVOfSHhzy1W0eeO8BpZTn8Jny9/cTRg2PcMsZP61XX8uvZtzD/ug4EIGUFI+894BSynP4TPmDfY0fNxz2WVJUw4JRZ7Ft1Saw2WDPHi1+pZRT+VT5Z8WFU9fUxsHD7VZH+YoX15STHh3G5PRhVkdRSvkInyr/TDdc5mHr3kbWlx3k+kkjENErepVSruFT5e+Od/V6cXU5QQF+XHFaktVRlFI+xKfKP25wCBEhAW4z46elvZPXP69ixvh4osKCrI6jlPIhPlX+IkJWbATba93jQq/FG6tpauvkel26WSnlYj5V/tB9pe/22iaMMVZHoWB1OZmx4eSmRFkdRSnlY3yu/LNiI2ho6aC+qc3SHJsqD1FUeYjr8/REr1LK9Xyu/I/O+LF63P/FNWWEBPpxuZ7oVUpZYEDlLyJXiUixiNhEJPdrz/1cREpFZJuIXNxj+yX2baUicu9A9n8yMmPDASy92KuptYOFhdXMmpBA5KBAy3IopXzXQI/8NwNzgI96bhSRMcC1wFjgEuBJEfEXEX/gCeBbwBjgOvtrXWZYeDDR4UGWTvd8o7CalvYurp+UYlkGpZRvG9DNXIwxW4BjjVnPBhYYY9qA3SJSCuTZnys1xuyyv2+B/bUlA8nRX5mxEWyzaMaPMYb5n5UxJn4wE5IiLcmglFLOGvNPBCp6PK60bzve9m8QkdtFZJ2IrKuvr3douMzYCHbUNmGzuX7Gz8c79rGttombpqbqiV6llGV6LX8RWS4im4/xNduZwYwx84wxucaY3JiYGId+dlZcBC3tXVQ1HHHo5/bFMx/vIiYimNk5CS7ft1JKHdXrsI8x5oKT+NwqILnH4yT7Nk6w3WW+mPGzt4nkoaEu2+/WvY18vGMf/++iTIID/F22X6WU+jpnDfssAq4VkWARSQNGAWuAtcAoEUkTkSC6TwovclKG4/pixo+LT/o++/FuQgL9yNcTvUopiw10quflIlIJTAaWiMg7AMaYYuAVuk/kLgXuMsZ0GWM6gbuBd4AtwCv217pUREggyUMHsaG8wWX7rGtsZWFhFVfnJus6Pkopyw10ts/rwOvHeW4u8I37EBpj3gLeGsh+HeH8U2J5cU05Ta0dRIQ4f67985/todNmuGVqmtP3pZRSvfG5K3yPunRCPO2dNpaV1Dp9Xy3tncxfVc5FY2JJjQ5z+v6UUqo3Plv+pyZHkRAZwpKiGqfv69X1lRw60sF3z0p3+r6UUqovfLb8/fyEGdnxfLSjnkMtHU7bT5fN8I9PdpOTPISJunqnUspN+Gz5A8zITqCjy/BOyV6n7WNZSS1l+1v47lnpelGXUspt+HT5T0iKJHnoIN504tDPsx/vIilqEBePjXXaPpRSqr98uvxFhBnjE1hZuo8Dh9sd/vkbyg+yruwgt0xNI8Dfp/9TK6XcjM830szseLpshqWbHT/08+zHu4kICeDq05N7f7FSSrmQz5f/2ITBpEWHsWRTtUM/t7j6EG9vruH6SSMIDx7Q5RRKKeVwPl/+IsLM7Hg+27nfYbd2tNkMv3xjM1GhQdx5ToZDPlMppRzJ58sfYGZ2AjYDSzc75sTvq+sr+by8gXu/dQqRoXqnLqWU+9Hyp3uht4zh4Sx2wKyfhpZ2Hlq6ldyUKK7Q+/MqpdyUlj9fDv2s3XOA2sbWAX3Ww+9so6GlnQdmj8PPT+f1K6Xck5a/3czsBIxhQMs9bKxo4MU15dw4JZUxCYMdmE4ppRxLy98uY3g4p8RFsGTTyZV/l83wy4WbiQ4P5n8uzHRwOqWUciwt/x4unZDA+rKDVJ/E7R0XrC2nqPIQ988YzWAXLBGtlFIDoeXfw8zseKD/Qz/7m9v409JtnJE+lFkT9N68Sin3p+XfQ8qwMMYnRvLCqjIqD7b0+X1/XLqVw22dPDh7nC7eppTyCFr+X/OL6aM5eLidy55Yybo9B0742s4uG4+/t4NX1lVy65lpjLLfGF4ppdydlv/XTB45jNfvmkJ4cADXPbOK/6yrOObrttc2MeepT3nk3e3MyI7nngtGuTipUkqdPC3/Y8gYHsEbd00lL20o//tqEXOXlNBlM0D30f5TH+xk5mOfUHnwCE9cfxpPXH8aoUG6fo9SynNoYx3HkNAg/nVzHg++WcIzH+9mR10z95w/it8sLmFjRQPfGhfHg5eNIzo82OqoSinVb1r+JxDo78cDs8eRFRfBrxcW88G2eoaEBvLYdadyaXa8ntxVSnksLf8+yJ+UwsiYcN4p3sv3p41keESI1ZGUUmpAtPz76Iz0YZyRPszqGEop5RB6wlcppXyQlr9SSvkgLX+llPJBWv5KKeWDtPyVUsoHafkrpZQP0vJXSikfpOWvlFI+SIwxVmfolYjUA2UD+IhoYJ+D4jibJ2UFz8rrSVnBs/J6UlbwrLwDyZpijIk51hMeUf4DJSLrjDG5VufoC0/KCp6V15Oygmfl9aSs4Fl5nZVVh32UUsoHafkrpZQP8pXyn2d1gH7wpKzgWXk9KSt4Vl5PygqeldcpWX1izF8ppdRX+cqRv1JKqR60/JVSygf5RPmLyIMiUiQihSLyrogkWJ3pRETkYRHZas/8uogMsTrT8YjIVSJSLCI2EXHbqXMicomIbBORUhG51+o8JyIiz4lInYhstjpLb0QkWUTeF5ES+9+De6zOdDwiEiIia0Rkoz3rb63O1BsR8ReRDSLypqM/2yfKH3jYGJNtjMkB3gR+ZXWgXiwDxhljsoHtwM8tznMim4E5wEdWBzkeEfEHngC+BYwBrhORMdamOqF/AZdYHaKPOoGfGGPGAGcAd7nxf9s24DxjzAQgB7hERM6wOFNv7gG2OOODfaL8jTGNPR6GAW59ltsY864xptP+cBWQZGWeEzHGbDHGbLM6Ry/ygFJjzC5jTDuwAJhtcabjMsZ8BBywOkdfGGNqjDGf279voruoEq1NdWymW7P9YaD9y227QESSgBnAs874fJ8ofwARmSsiFUA+7n/k39MtwNtWh/BwiUBFj8eVuGlBeTIRSQVOBVZbm+T47MMohUAdsMwY47ZZgb8CPwVszvhwryl/EVkuIpuP8TUbwBhznzEmGSgA7rY2be957a+5j+5fqwusS9q3rMq3iUg48F/gR1/7TdutGGO67MO/SUCeiIyzOtOxiMhMoM4Ys95Z+whw1ge7mjHmgj6+tAB4C/i1E+P0qre8InITMBM431h8MUY//tu6qyogucfjJPs25QAiEkh38RcYY16zOk9fGGMaROR9us+tuOOJ9anALBGZDoQAg0VkvjHm247agdcc+Z+IiIzq8XA2sNWqLH0hIpfQ/eveLGNMi9V5vMBaYJSIpIlIEHAtsMjiTF5BRAT4B7DFGPOo1XlORERijs6cE5FBwIW4aRcYY35ujEkyxqTS/ff1PUcWP/hI+QMP2YcpioCL6D6D7s4eByKAZfbpqf9ndaDjEZHLRaQSmAwsEZF3rM70dfaT53cD79B9QvIVY0yxtamOT0ReAj4DskSkUkRutTrTCUwFvgOcZ/+7Wmg/WnVH8cD79h5YS/eYv8OnUHoKXd5BKaV8kK8c+SullOpBy18ppXyQlr9SSvkgLX+llPJBWv5KKeWDtPyVUsoHafkrpZQP+v8CEeE+E7AO8AAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]}]}