Created
June 10, 2017 06:17
-
-
Save ruoyu0088/03f8e4eb51e0c8e743bfe5480a21c7ae to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"%matplotlib inline" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# 拉格朗日多项式插值" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"`scipy.interpolate.lagrange(x, y)`返回经过`(x, y)`所有点的多项式系数。" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFkCAYAAACuFXjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xt8jvXjx/HX5545bJbz6SvZCE3lsEnkNHLWnCk5JodC\nK9L3FykKpYNTzoeEYkkISQ45xnLY5FsshPGNEBVDGPfn9wft6zg226572/v5eNx/3Nd9Hd7bPOy9\n6/pcn8tYaxERERG5FZfTAURERMSzqSyIiIhIglQWREREJEEqCyIiIpIglQURERFJkMqCiIiIJEhl\nQURERBKksiAiIiIJUlkQERGRBKksiIiISIISVRaMMc8ZY7YbY05eeW00xtS/zTYhxphIY8w5Y8xu\nY0zHu4ssIiIiqSmxZxb+C/wfEAQEA6uAhcaYwJutbIzxB74CvgXKAqOBqcaYOknMKyIiIqnM3O2D\npIwxJ4C+1tqPb/LZu0ADa22Zq5aFAzmstQ3v6sAiIiKSKpI8ZsEY4zLGPAX4ABG3WK0SsPK6ZcuA\nykk9roiIiKSuTIndwBjzEJfLQVYgFmhmrf35FqsXBI5et+wocI8xJou19vwtjpEHqAfEAOcSm1FE\nRCQDywr4A8ustSeSY4eJLgvAz1wef5ADaAnMNMZUT6AwJEU9YFYy7k9ERCSjaQvMTo4dJbosWGsv\nAvuuvN1mjKkIvAg8f5PVjwAFrltWADh1q7MKV8QAfPrppwQG3nTspKQxvXv3ZuTIkU7HkGSin2f6\nop9n+hIdHU27du3gyu/S5JCUMwvXcwFZbvFZBNDgumV1ufUYh3+cAwgMDCQoKOju0olHyJEjh36W\n6Yh+numLfp7pVrJdxk9UWTDGvA0sBQ4Cflw+xVGDywUAY8w7wL+stf/MpTAR6HnlrohpwONcvnSh\nOyFERETSiMSeWcgPzAAKASeB/wB1rbWrrnxeECjyz8rW2hhjTCNgJBAG/Ao8a629/g4JERER8VCJ\nKgvW2i63+fyZmyxbx+UJnERERCQN0rMhJFW0adPG6QiSjPTzTF/085TbUVmQVKH/jNIX/TzTF/08\n5XaS424IERFJYQcPHuT48eNOxxAPkDdvXu67775UPabKgoiIhzt48CCBgYGcPXvW6SjiAXx8fIiO\njk7VwqCyICLi4Y4fP87Zs2c1UZ3ET7h0/PhxlQUREbmRJqoTp2iAo4iIiCRIZUFEREQSpLIgIiIi\nCVJZEBERkQSpLIiIiCTS2rVrcblcrFu3zukoqUJlQUREHDVjxgxcLtdNX/3793c024QJE5gxY8ZN\nPzPGpHIa5+jWSRERcZwxhsGDB+Pv73/N8oceesiZQFeMHz+efPny0bFjx2uW16hRg7///pvMmTM7\nlCx1qSyIiKQz1toU/as3pfZfv379NDWPREYpCqDLECIi6UJsbCxhYQMJCKhNkSJNCQioTVjYQGJj\nY9PE/hNy4MABXC4XM2fOvOEzl8vFW2+9Ff9+0KBBuFwu9u7dS6dOnciVKxc5c+akc+fOnDt37obt\nP/30Ux599FF8fX3JnTs3NWrUYOXKlQAEBASwY8cO1qxZE39ZpFatWsCtxyzMnTuXChUq4OPjQ758\n+Wjfvj2HDx++Zp1OnTrh5+fH4cOHadq0KX5+fuTPn59XXnkFa+1df79Sgs4siIikcbGxsVSu3ILo\n6D643YMAA1jGjVvGqlUtiIiYh5+fn8fu/x8nT57kxIkT1yzLkydPovbxzxmP1q1bU6xYMYYNG0ZU\nVBRTp06lQIECvPPOO/Hrvvnmm7z55ptUqVKFwYMHkzlzZjZt2sSqVauoXbs2o0ePplevXvj5+TFg\nwACstRQoUOCGY/1j+vTpdO7cmUcffZRhw4Zx9OhRRo0axcaNG9m2bRv33HNP/HZut5t69epRqVIl\nhg8fzsqVKxkxYgT3338/3bt3T9TXnCqstR73AoIAGxkZaUVEMrrIyEib0P+JL7zwhnW5llqwN7xc\nrq9tWNjAuzp+Su9/+vTp1hhzw8vlcllrrY2JibHGGDtjxowbtjXG2DfffDP+/aBBg6wxxnbt2vWa\n9Zo3b27z5csX//6XX36xXl5etmXLlglme+ihh2zNmjVvWL5mzRrrcrns2rVrrbXWxsXF2QIFCtiy\nZcva8+fPx6+3ZMkSa4yxgwYNil/WqVMn63K57NChQ6/ZZ1BQkH3kkUcSzHO7fwtXrwME2WT6vazL\nECIiadzixRtwu+vd9DO3uz6LFm3w6P3D5b+2J0yYwMqVK+NfK1asSPK+rv/rvFq1apw4cYLTp08D\nsGDBAqy1vPHGG3edHWDr1q0cO3aMHj16XDOWoWHDhjzwwAMsWbLkhm1ulnHfvn3Jkie56TKEiEga\nZq0lLs6Xy5cGbsYQF+eT5EGJKb3/qz3yyCPJNsDx+icy5sqVC4A///yT7Nmzs2/fPlwuV7I9xfPA\ngQMYYyhZsuQNnz3wwANs2HBtocqaNesNl1hy5crFn3/+mSx5kpvOLIiIpGHGGLy9z3D5rPPNWLy9\nzyT5F3lK7/9OM9yM2+2+5TZeXl43XW49ZADhrfJ5KpUFEZE0LjS0Ci7Xspt+5nJ9Q+PGVT16/7fz\nz1mBv/7665rlBw4cSPI+ixcvjtvtZufOnQmud6clqGjRolhr2bVr1w2f7dq1i6JFiyYpp6dQWRAR\nSeOGDu1LYOAIXK6l/O8MgMXlWkpg4EiGDHnZo/d/O35+fuTNm/eG2xTHjRuX5DMaTZs2xRjDW2+9\nleDZBl9f3xtKys1UqFCB/PnzM3HiROLi4uKXL126lOjoaJ544okk5fQUGrMgIpLG+fn5ERExjwED\nhrNo0Qji4nzw9j5L48ZVGDLk7m9rTOn9w+0vD3Tp0oVhw4bRtWtXKlSowLp169izZ0+SLysUL16c\n1157jSFDhlCtWjWaN29OlixZ2LJlC4ULF2bo0KEABAcHM3HiRIYOHcr9999P/vz5qVmz5g2ZM2XK\nxLvvvkvnzp2pXr06bdq04ciRI3z44YcUK1aMl156KUk5PYXKgohIOuDn58fo0YMYPTplZlhM6f3f\nbn9vvPEGx48f54svvmDu3Lk0bNiQpUuXkj9//iRnefPNNylWrBhjxoxhwIAB+Pj4UKZMGTp06HDN\ncQ8ePMj7779PbGwsNWrUiC8L1x+3Y8eO+Pr6MmzYMF599VV8fX1p0aIFw4YNi59j4XZfr6c+b8J4\nymCPqxljgoDIyMjINDX1p4hISoiKiiI4OBj9nyh38m/hn3WAYGttVHIcV2MWREREJEEqCyIiIpIg\nlQURERFJkMqCiIiIJEhlQURERBKksiAiIiIJUlkQERGRBKksiIiISIJUFkRERCRBKgsiIiKSIJUF\nERERSZDKgoiIpBsHDhzA5XIxc+ZMp6OkKyoLIiLiqBkzZuByueJf2bJlo1SpUrzwwgscO3bM6XiC\nHlEtIiIewBjD4MGD8ff359y5c3z33XdMmDCBpUuX8tNPP5E1a1anI2ZoKgsiIuIR6tevH//Y5c6d\nO5M7d25GjhzJwoULefLJJx1Ol7HpMoSIiHikWrVqYa1l//79AOzfv59WrVqRJ08efH19qVy5Ml9/\n/XWC+5g+fToul4vt27ff8Nnbb79NpkyZ+O233wAICQmhTJkyREdHU7NmTXx9fbn33nt5//33b9j2\n999/59lnn6VgwYJky5aNcuXK3TBO4p/xEyNGjGD8+PEUL14cX19f6tWrx6FDhwAYPHgwRYoUwcfH\nh6ZNm/LXX38l6XuV0lQWRETEI/3yyy8A5M2bl2PHjlG5cmVWrFhBr169ePvttzl//jyNGzdm4cKF\nt9xHy5YtyZYtG7Nmzbrhs9mzZ1OrVi0KFSoEXL4U8scff9CgQQPKly/PiBEjCAwM5NVXX2XZsmXx\n2507d44aNWowa9Ys2rdvzwcffEDOnDnp1KkTY8aMueE4n376KRMmTCAsLIy+ffuydu1aWrVqxYAB\nA1i+fDmvvvoq3bt3Z/HixfTt2/duv20pw1rrcS8gCLCRkZFWRCSji4yMtOn5/8Tp06dbl8tlV61a\nZY8fP25//fVX+9lnn9m8efPa7Nmz28OHD9uXXnrJulwuu3HjxvjtTp8+bYsVK2aLFSsWvywmJsYa\nY+yMGTPilz399NP23nvvveaYUVFR1hhjZ86cGb8sJCTEulwuO2vWrPhlFy5csIUKFbKtWrWKXzZq\n1CjrcrlseHh4/LKLFy/axx57zN5zzz329OnT12QpUKCAjY2NjV+3f//+1hhjy5cvby9dunRNzqxZ\ns9oLFy7c8nt1J/8W/lkHCLLJ9HtZYxZERNKRs3Fn+fn4zyl+nAfyPoCPt0+y7c9ay+OPPx7/3hiD\nv78/4eHhFCpUiKVLl1KxYkUqV64cv46vry/dunWjf//+7Ny5k9KlS9903x06dOCzzz5j9erV1KxZ\nE4BZs2bh4+ND8+bNr1k3e/bsPP300/Hvvb29qVixIvv27YtftnTpUgoWLMhTTz0Vv8zLy4uwsDCe\nfvpp1q5dS8OGDeM/a926NdmzZ49//+ijjwLQvn17XC7XNcs/++wzDh06hL+//x1931KLyoKISDry\n8/GfCZ4cnOLHiewWSVChoGTbnzGG8ePHU6JECTJlykSBAgUoVapU/OcHDhygUqVKN2wXGBgY//mt\nykKdOnUoWLAgs2bNombNmlhr+eyzz2jatCm+vr7XrHvvvffesH2uXLn48ccfr8lSokSJm2ax1nLg\nwIFrlhcpUuSa9zly5Ljpsf5Z/ueff6osiIhIynkg7wNEdotMleMkt0ceeST+bojk5HK5ePrpp5k6\ndSrjx49n/fr1HD58mHbt2t2wrpeX1033YS9fIk+SW+0zJY6VUlQWRETSER9vn2T9i99TFC1alF27\ndt2wPDo6Ov7zhHTo0IERI0awePFivv76a/Lnz0/dunWTnOXqMw2JzZIW6W4IERHxeA0bNmTz5s1s\n2rQpftmZM2eYPHkyAQEBt7wE8Y+HH36Yhx9+mClTpjBv3jzatGlzzXiBxGY5cuQIc+bMiV926dIl\nxowZg5+fHzVq1EjSfj2ZziyIiIjjbnfq/dVXXyU8PJz69esTFhZG7ty5mT59OgcOHGD+/Pl3dIwO\nHTrQt29fjDG0bds2yVm7devGpEmT6NSpE1u3bsXf35+5c+cSERHB6NGjbxgHkRieeAkCdGZBREQ8\ngDEmwc/z589PREQEdevWZezYsfTv35+sWbPy1Vdf0bhx4zvaV9u2bfHy8qJUqVJUqFAhUTmuXp41\na1bWrl1L27ZtmTlzJn379uWvv/5i+vTp9OrV64btbrbPOzmOJzGe2GKMMUFAZGRkZIoMdhERSUui\noqIIDg5G/yfenRMnTlCoUCEGDRpE//79nY6TJHfyb+GfdYBga21Uchw3UWcWjDH9jDGbjTGnjDFH\njTELjDElb7NNDWOM+7rXJWNM/ruLLk7xxIIpInI7H3/8MW63+6Z3QUjCEjtmoRowBth6Zdt3gOXG\nmEBr7d8JbGeBkkBs/AJr9dzRNCQ2NpbXXvuAxYs3EBfni7f3GUJDqzB0aF/8/PycjicickurV69m\nx44dvP322zRr1oz77rvP6UhpTqLKgrW24dXvjTGdgGNAMPDdbTb/3Vp7KlHpxCPExsZSuXILoqP7\n4HYPAgxgGTduGatWtSAiYp4Kg4h4rLfeeouIiAiqVq3Khx9+6HScNOluBzjm5PJZgz9us54BfjDG\nHDbGLDfGPHaXx5VU9NprH1wpCvW5/KMEMLjd9YmO7s2AAcOdjCcikqDVq1dz7tw5Vq5cGf/QKEmc\nJN86aS4P2RwFfGet3ZnAqr8B3bl86SIL0BVYY4ypaK39IanHl+RjreXXU78SfTyaw7GHOXL6CEdP\nH+XImSP88fcfrDu/BXeHNeAaDK6L4PaGczngfA7c5/2YfmgZ92/Kg39O//iXXxadaRARSS/uZp6F\n8UBpoEpCK1lrdwO7r1r0vTGmONAb6JjQtr17946fK/sfbdq0oU2bNkkKLJeLwe4Tu1kTs4ao36L4\n6fef+OnYT5w6/78rRLmy5qJg9oIUyF6APNnyYM76wJn7wJ0JrBe44iDrSch+BPLs4kz233llxSuc\nv3Q+fh/5fPJRpkAZyhcsT/lC5SlXsByl8pTCy3Xz6U1FRCTxwsPDCQ8Pv2bZyZMnk/04SSoLxpix\nQEOgmrX2tyTsYjO3KRkAI0eO1G1CyeBw7GGW7F7CmgNrWL1/Nb+d/g0v48VD+R/iofwPEVoylAfz\nPUjpfKW59557yZIpyzXbB/xfbWJiZvK/SxBXsxTxr8Pefcs5evoo+//aT8xfMew5sYftR7fzRfQX\nfBDxAQB+mf2oVrQaIUVDCPEPoXyh8mRyaV4wEZGkutkf0FfdOplsEv0/9ZWi0ASoYa09mMTjluPy\n5QlJIUdPH+WLnV8wZ8ccvjv4HcYYggsF075Me2oG1KRKkSp3fKkgNLQK48YtuzJm4Vou1zc0blwV\nl3FRyK8QhfwK8ViRa4ek/Pn3n/xw5Ae+//V71hxYw6C1gzi78ix+mf14vNjjNCnVhEYlGpHPN1+y\nfO0iIpK8ElUWjDHjgTZAY+CMMabAlY9OWmvPXVnnbaCwtbbjlfcvAvuBHUBWLo9ZqAnUSZavQOLF\nXYpjwc8LmBQ5iTUxa3AZF3WK1WFak2k0KdWEXNlyJWm/Q4f2ZdWqFkRH26sGOVpcrm8IDBzJkCHz\nEtw+V7Zc1AyoSc2AmvSr1o8Lly4QeTiSVftXsWTPEjov7AzAY0Ueo3GpxrR+sDX+Of2TlFUkPfvn\nQUWScTn1byBRMzgaY9xcvvvhes9Ya2deWedjoKi1ttaV968A3YB/AWeB/wBvWmvXJXAczeCYCEdP\nH2VK1BQmbp3IodhDVLuvGh3KdqDZA83I45MnWY4RGxvLgAHDWbRoA3FxPnh7n6Vx4yoMGfLyXd82\nefT0UZbsWcKiXYtYvnc5f1/8mypFqtD24ba0erAVeX3yJsvXIJJWHTx4kMDAQM6ePet0FPEAPj4+\nREdH33K+iJSYwVHTPadhe//Yy5D1Q5j942y8jBftyrSjV8VelClQJkWPa61NsfnLT184zZc/f8ns\nH2ezfO9yjDE0LNGQrkFdaXB/Aw2QlAzr4MGDHD9+3OkY4gHy5s2b4MRSKgsCwIG/DjB43WCm/zCd\n/L75ebnyy3Qu3znJlxk81bEzx/h8x+dM2zaNbUe2ce8999KlfBc6l+9MkRxFnI4nIuKRVBYyuN9i\nf2PIuiFMiZpCzqw56Ve1H89VeI5s3tmcjpbiIg9HMjlyMrN/ms3ZuLM0LtWYlx59iepFq3vsU9pE\nRJzg+IOkxBlxl+IYETGCUmNLEf5TOG/VfIt9L+6jd+XeGaIoAAT/K5hJoZM43Ocw4xuOZ/eJ3YTM\nCCF4cjCfbP+EC5cuOB1RRCTdUlnwcGtj1lJ+UnleWfEKHcp2YG/YXl6t+irZM2d3Opoj/LL40b1C\nd356/ieWtVtGgewF6PBlB/xH+fPehveumVxKRESSh8qCh/r9zO+0m9+OkBkh3JPlHrZ23crYhmPT\n3biEpDLGULd4XZa2XcqOHjtoVKIRA1YNoOiooryx+g2On9VAMBGR5KKy4IEW/ryQhyY8xDe/fMO0\nxtP4rvN3lC9U3ulYHqt0vtJMaTyFfS/u45lyzzA8YjhFRxWlz7I+HDl9xOl4IiJpnsqCBzl1/hSd\nF3am6ZymPFr4UXb02MEz5Z/BZfRjuhP33nMvI+qN4MBLB+hTqQ8fbfuIYqOL8cryV/j9zO+33M4T\nB/mKiHgS/RbyEGti1lBmQhm+2PkF0xpPY+FTCymQvcDtN5Qb5PXJy+Bag4l5MYaXK7/MxMiJBIwO\noN/Kfpw4ewK4PMlUWNhAAgJqU6RIUwICahMWNpDY2FiH04uIeB7dOukwt3UzdN1QBq4ZSPWi1Zne\ndLqmOk5mJ86e4IONHzBm8xi8XF68WOFFvnh5A7t+egW3ux7/m756GYGBI4iImHfXs1KKiDhFt06m\nM3/8/QdPzH6CgWsGMihkEKs6rlJRSAF5fPLwTu132PfiPjqV7cTQ74YSXWcb7qD94Lp4ZS2D212f\n6OjeDBgw3NG8IiKeRmXBIVsPbyVoUhCbDm1iadulvFHjDY1NSGH5ffMzusFoCs2rBPuegEY9oWdp\nKP0F/zzyxO2uz6JFG5wNKiLiYfTbKRVcf6ln2rZpVJlWhfy++YnqFkW9++s5lCzjsdbCn3lhwUyY\n+AOcKAmtW0HnqnDv94AhLs5Hgx5FRK6SqEdUy52LjY3ltdc+YPHiDcTF+eLtfYYnQh8jU71YRm0d\nRbegbnzY4EOyZMridNQMxRiDt/cZwMLRMjB7CRRbCXX7QpfK8NOT8PNxTSEtInIVnVlIAbGxsVSu\n3IJx4yoTE7OCQ4cWEnNoIWOPrWLUllG8E/IOE5+YqKLgkNDQKrhcy/63YF9tmBQJX34M963ktxab\neHXlq8Se150RIiKgspAiXnvtA6Kj++B21wcMZD8CnWpC8R8wcwby2/xz+svVQUOH9iUwcAQu11L+\nGauAdeH6TwECV5bl35Vf4cNNH1JybEk+3vYxbut2NK+IiNNUFlLA4sUbrtySB+TbCV0eBb9DMG09\n9ueBGkDnMD8/PyIi5tGr1yb8/etSuHAT/P3r0qvXJjZ99yXv1HuHXb12UdO/Jp0XdabilIpsOKif\nmYhkXCoLycxaS1ycL2Cg8GZ4phqcywlTNsOR8mgAnWfw8/Nj9OhB7N+/gv/+90v271/B6NGD4udX\nKJKjCLNbzOa7Z77DGEPVj6vSfkF7fov9zeHkIiKpT2UhmcUPoAv4Fjo8DscfgOlrIbbwlTUs3t5n\ndBnCgyT0s6hyXxU2ddnE1NCpLPtlGSXHluT9De/rkdgikqGoLKSA0s3zQNv6cLAqfLL88pmFK1yu\nb2jcuKqD6SSxXMbFs0HPsvuF3XQu15l+3/ajzIQyLN+73OloIiKpQmUhmc34YQbf5PiCe37Lg5nz\nPMT5XPnE4nItJTBwJEOGvOxoRkmanFlzMrrBaLZ130bB7AWp92k9Wn7ekoMnDzodTUQkRaksJKMp\nkVPotLATz5Z/lgPDf+aFHlE3DKDTcwfSvocLPMzqjquZ3Xw2G/+7kcBxgbyz/h3OXzzvdDQRkRSh\nB0klkymRU+j2VTd6PtKTMQ3GXHMd3FqrMQrp1Knzp3hzzZuM3jSa4rmLM67hOGoXq+10LBHJwPQg\nKQ81NWrqLYsCJDyATtK2e7Lcw/B6w/nhuR8omL0gdT6pw1NfPMXh2MNORxMRSTYqC3fpo6iP6Lq4\nKz0q9LhpUZCM4aH8D7Gm4xpmNp3J6pjVPDD2AUZGjOSi++LtNxYR8XAqC3fh420f03VxV56v8Dxj\nG45VUcjgjDG0L9ueXb120aFsB15e/jLBk4PZ+N+NTkcTEbkrKgtJ9PmOz3l20bN0C+6moiDXyJk1\nJ2MbjmVL1y1k8cpClWlV6LqoKyfOnnA6mohIkqgsJMHSPUtpN78dbcu0ZXyj8biMvo1yo+B/BRPx\nbATjG45n7s65lBpbimnbpulZEyKS5ui3XCKtP7CeFp+3oEGJBkxrPE1FQRLk5fLi+UeeZ1evXTQs\n0ZBnFz1L9Y+r8+PRH52OJiJyx/SbLhEiD0fSaHYjKt1biTkt5+Dt5e10JEkjCmQvwMxmM1ndcTUn\n/j5B+UnleWX5K5y+cNrpaCIit6WycId+Pv4z9WfVJzBfIAufWkjWTFmdjiRpUIh/CNuf287gmoMZ\nu2UspceVZkH0Aj1YTEQ8msrCHTgce5h6n9Yjv29+lrZdil8WzcAoSZfZKzP9qvVjZ4+dlClQhuaf\nNyc0PJT9f+53OpqIyE2pLNzGqfOnaDirIZfcl/im7Tfkzpbb6UiSTgTkCmBxm8XMbz2f7Ue38+D4\nB3l7/dt6oqWIeByVhQRcuHSB5nOaE/NXDEvbLqVIjiJOR5J0xhhDs8BmRPeMpucjPXlj9RuUnViW\n1ftXOx1NRCSeysItuK2bzgs7s/7ger586kseLvCw05EkHcueOTvv132fbd23kSdbHmrNrEW7+e04\ncvqI09FERFQWbqXfyn7M/nE2nzT7hBD/EKfjSAbxcIGHWffMOj5u8jHL9i6j1NhSjN08lkvuS05H\nE5EMTGXhJiZsmcB7G99jRL0RtH6wtdNxJINxGRedynViV69dtHmoDWFLw6g4tSKbft3kdDQRyaBU\nFq6zdM9Sei3tRVjFMF6q9JLTcSQDy50tNxOfmEjEsxFYa6n8UWW6Le6maaNFJNWpLFxl+5HttP6i\nNY1KNGJEvRFOxxEB4NF7H2VL1y2MaTCGz3d8TsmxJZkSOeWm00ZrvgYRSQkqC1ccjj3ME+FPUCJ3\nCWa3mI2Xy8vpSCLxvFxe9KzYk129dhFaMpRuX3Wj8keV2Xp4K7GxsYSFDSQgoDZFijQlIKA2YWED\niY2NdTq2iKQTKgvAmQtnCA0PxVrL4jaLyZ45u9ORRG6qQPYCTG86nfXPrOfcxXNUnFKRoj0DGTvt\nQWJiVnDo0EJiYlYwblxlKlduocIgIskiw5cFt3XTdn5bdp/YzZKnl1D4nsJORxK5rar3VSWyWyTV\nztTjz8J/YXs+D8GTwVwCDG53faKjezNgwHCno4pIOpDhy0L/b/uzePdiPmvxGWULlnU6jsgdy+TK\nxMG5cTDmF9jVGEKfg66Pwr3fA+B212fRog0OpxSR9CBDl4UZP8zg3Q3v8kGdD2hUspHTcUQSxVpL\nXJwvnCkICz+GqRsBC10qQ9NOkP0ocXE+GvQoInctw5aFDQc30O2rbjxb/lndIilpkjEGb+8zwJUy\n8GtlmLIZFk2GEkvghZKcLbeTOHecozlFJO3LkGUh5q8Yms1pRqV7KzG+0XiMMU5HEkmS0NAquFzL\n/rfAekFUVxizG36oyV8V9vLwhIf5es/XzoUUkTQvw5WF2POxhIaH4pfFj3mt55HZK7PTkUSSbOjQ\nvgQGjsDlWkr8GQYsrgvf8+DBv9nYYSOF/QrTaHYjGs5qyM/Hf3YyroikURmqLFxyX6Lt/LYcPHmQ\nxW0Wk9fHZ4MRAAAgAElEQVQnr9ORRO6Kn58fERHz6NVrE/7+dSlcuAn+/nXp1WsTERHzqFSsEt92\n+Jb5refz8/GfeXjCw/RZ1oe/zv3ldHQRSUOMJw5+MsYEAZGRkZEEBQUl2377rezHexvf46s2X9Gg\nRINk26+Ip7DW3vKy2rmL5xj1/SiGrBtCNu9svBnyJt2Cu5HJlSmVU4pISoqKiiI4OBgg2FoblRz7\nzDBnFmb9ZxbDNgzjvdrvqShIupXQ+JusmbLyatVX2fPCHkJLhtLr616UnViWb375JhUTikhalCHK\nwqZfN/HsomfpVK4TfSr3cTqOiKMK+RViWpNpbO22lXw++WgwqwENZjVgx7EdTkcTEQ+V7svCoVOH\naDanGUGFgpjYaKLufBC5IqhQEKs7rmZ+6/nsObGHMhPL0H1xd46cPuJ0NBHxMIkqC8aYfsaYzcaY\nU8aYo8aYBcaYknewXYgxJtIYc84Ys9sY0zHpke/c33F/03ROU7xcXsx/cj5ZMmVJjcOKpBnGGJoF\nNmNnz50MrzucuTvnUmJMCYasG8LZuLNOxxMRD5HYMwvVgDHAo0BtwBtYbozJdqsNjDH+wFfAt0BZ\nYDQw1RhTJwl575i1ls6LOrPj2A4WPrWQgtkLpuThRNK0zF6ZeanSS+wN20u3oG68tfYtSowpwUdR\nH3HJfcnpeCLisESVBWttQ2vtJ9baaGvtj0An4D4gOIHNngf2WWv/ba3dZa0dB3wB9E5q6Dvx9vq3\n+eynz5jRdAZBhZLvjgqR9CxXtlwMrzecn3v9TPWi1emyuAtlJ5blq91fadpokQzsbscs5OTyTDB/\nJLBOJWDldcuWAZXv8ti3tCB6AQNWD2BgjYG0erBVSh1GJN0qlqsY4S3C2dxlM/l98xMaHkrIjBA2\n/brJ6Wgi4oAklwVzeaTgKOA7a+3OBFYtCBy9btlR4B5jTLIPIvjP0f/QfkF7WpZuyRs13kju3Ytk\nKI8UfoRvO3zL109/zZ9//0mljyrRfE5zon+PdjqaiKSiu5mNZTxQGqiSTFlu0Lt3b3LkyHHNsjZt\n2tCmTZubrn/szDEahzemRJ4STG8yHZdJ9zd7iKQ4YwwNSjSgbvG6zP5xNq+vfp2HJjzEM+WeYWCN\ngRTJUcTpiCIZVnh4OOHh4dcsO3nyZLIfJ0kzOBpjxgKhQDVr7cHbrLsWiLTW9rlqWSdgpLU21y22\nSfQMjucvnqf2J7XZc2IPm7tu5r4c993hVyMiiXH+4nkmbp3IkPVDiD0fS89HevJq1VfJ55vP6Wgi\ngofM4HilKDQBat6uKFwRATx+3bK6V5YnC2st3b/qzpZDW1jw5AIVBZEUlCVTFl6s9CL7wvbRr2o/\npkRNodiHxXhj9RucPJf8f9GIiPMSO8/CeKAt8DRwxhhT4Mor61XrvG2MmXHVZhOBYsaYd40xpYwx\nPYCWwIhkyA/AuxveZcb2GXzU+CMqF0mxcZMichW/LH4MDBnI/hf383yF53l/4/sEjA7g3e/e5cyF\nM07HE5FklNgzC88B9wBrgMNXvVpftU4hIP4iprU2BmjE5XkZfuDyLZPPWmuvv0MiSRZEL6Dft/14\nvfrrtC3TNjl2KSKJkMcnD+/VeY+9YXt56qGneH316xT7sBgjI0byd9zfTscTkWSQpp86ue23bVT9\nuCqNSjTis5afaUCjiAeI+SuGwWsHM2P7DApkL0D/qv3pEtRFM6iKpBKPGLPgKQ7HHiY0PJTS+Uoz\nvanufBDxFP45/fmoyUf83OtnHg94nLBvwigxpgQTt07k/MXzTscTkSRIk79hT184TWh4KACLnlqE\nj7ePw4lE5Hr3576fmc1msqPHDqoVrUaPJT1UGkTSqDRXFi66L/LkF0+y58Qeljy9hEJ+hZyOJCIJ\neCDvA8xqPosdPXZQ9b6q8aVhwpYJKg0iaUSaKgvWWnou6cnyvcv5ovUXlC1Y1ulIInKHAvMFMrvF\nbH7q8RNV7qtCz697UvzD4ozdPJZzF8/ddBtPHFMlkhGlqbLw7oZ3mRw1mclPTKZu8bpOxxGRJCid\nrzThLcLZ2XMntQJq8eI3LxIwOoCRESM5c+EMsbGxhIUNJCCgNkWKNCUgoDZhYQOJjY11OrpIhpVm\n7oaY/eNs2s5vy8AaAxkUMsjRfCKSfPac2MM7373DzO0zyZU1F67NOfl9yTDsueaAASwu1zICA0cQ\nETEPPz8/pyOLeLQMezfEqv2reGbhM3Qs25GBNQY6HUdEklGJPCWY1mQae17YQ8E/i3IsMAb7YhcI\nGQTZ/gAMbnd9oqN7M2DAcKfjimRIHl8WNh/aTJPPmhDiH8Lk0MlcftiliKQ3AbkCOP1ZThi9H37o\nBFXeh5eKQp1XIPtvuN31WbRog9MxRTKku3nqZIqr16Yjp5/aS9l7yzC/9Xwye2V2OpKIpBBrLXFx\nvhB7LywbCev7Q6VRUHEsPDoGtnXm718ur6c/GkRSl0efWTj+yBHOHSvIyfHZcJ93Ox1HRFKQMQZv\n7zPAlXFUZ/PBqqEw8iCsfQNKz+Xok9/S8cuO7Di2w9GsIhmNR5cFLmaFTyLY/Z9/61qlSAYQGloF\nl2vZtQvP54D1/TGjJ1P973qsiVnDQxMeoslnTfj+1++dCSqSwXh2WVgyEc4U0LVKkQxi6NC+BAaO\nwOVaSvwZBiwu11JKl5jAV69/zi9hv/Bxk4/ZfWI3lT+qTM0ZNVn2yzLNySCSgjy7LJz+Z3ZGQ1yc\nj/4zEEnn/Pz8iIiYR69em/D3r0vhwk3w969Lr16b4m+bzOyVmU7lOrGjxw7mtZ7H6QunqT+rPkGT\ng5jz0xwuuS85/WWIpDsePc8CRAJBgMXfvw779yfLU61FJI24k8GM1lpWx6xm2HfDWLFvBcVzFafv\nY33pVK4TWTNlTaWkIp4jw86z4HJ9Q+PGVZ2OISKp7E7uejDGUCugFsvbL2dr160EFQqix5Ie+I/y\n55317/DXub9SIalI+ubhZeHytcrAwJEMGfKy02FExMMF/yuYz1t9zq5eu2hSqgmD1g7ivpH38cry\nVzgce9jpeCJplkeXhUKFelxzrVJE5E6UyFOCSaGTiHkxhh6P9GBy1GQCRgfQZVEXdh3f5XQ8kTTH\no8csXP1sCBGRpDp57iSTIicx8vuRHD19lGaBzfi/Kv9HxcIVnY4mkuwy7JgFEZG7kSNrDv5d5d/s\nf3E/k56YxI9Hf+TRqY9Sa0Ytlu9drjutRG5DZUFEMoysmbLSNbgr0T2jmdtqLqfOn6Lep/UInhzM\n5zs+122XIregsiAiGY6Xy4uWpVuypesWVrZfSR6fPDz5xZOUGluKyZGTOXfxnNMRRTyKyoKIZFjG\nGB4v9jgr2q9gS9ctlC9Unue+eo5io4vx/ob3iT0f63REEY+gsiAiAlT4VwXmtppLdM9oGpZoyGur\nXuO+Uffx+qrXOX72uNPxRBylsiAicpVSeUsxtfFU9r24j2fKPcOI70dQdFRR+izrw6FTh5yOJ+II\nlQURkZu49557GVFvBAdeOsDLlV/m4x8+JmB0AN0Wd2PvH3udjieSqlQWREQSkNcnL2/VfIsDLx1g\nSK0hLNy1kFJjS9FhQQeif492Op5IqlBZEBG5A/dkuYd/V/k3MS/GMLLeSFbtX8WD4x+k9dzWbD+y\n3el4IilKZUFEJBGyeWfjhUdfYG/YXiY9MYmth7dSblI5mn7WlKjfEp4sT5M/SVqlsiAikgRZMmWh\na3BXdr+wm+lNprPj9x0ETw4mNDyUrYe3xq8XGxtLWNhAAgJqU6RIUwICahMWNpDYWN2WKWmHyoKI\nyF3I5MpEx3Idie4ZzSfNPmHPiT08MuURGs1uxNo9a6lcuQXjxlUmJmYFhw4tJCZmBePGVaZy5RYq\nDJJmqCyIiCSDTK5MtCvTjh09djCr+Sz2/rGXkNkh7CgTizt/AcBcWdPgdtcnOro3AwYMdzKyyB1T\nWRARSUZeLi+efvhpdvTYQd51pSHvCXguCJ5sDvl/jF/P7a7PokUbHEwqcudUFkREUoDLuMiy634Y\ntxMWTIcC2+H5stDiaci9BzDExflo0KOkCSoLIiIpwBiDt/cZcHvB9o4w9mf4aiIUXQe9AqHxs5ic\nJzDG3H5nIg5TWRARSSGhoVVwuZZdfuP2hshu8OEvsPwDKDWPwy2/p8+yPnr2hHg8lQURkRQydGhf\nAgNH4HItBa5cbriYBdfmUjywLIj+VfoxNWoqxT8szpB1Qzhz4YyjeUVuRWVBRCSF+Pn5ERExj169\nNuHvX5fChZvg71+XXr02sfm7hQyuM5h9L+6jc7nODF43mOIfFmfClgnEXYpzOrrINYwnDq4xxgQB\nkZGRkQQFBTkdR0QkWVhrbzlGIeavGAauGcgn2z+hZJ6SvFfnPUJLhmpMgyRaVFQUwcHBAMHW2oSn\nFb1DOrMgIpJKEvrF75/TnxlNZxDVPYoiOYrQ5LMmhMwIYcuhLamYUOTmVBZERDxIuYLlWN5uOUvb\nLuWPv/+g4tSKtJvfjl9P/ep0NMnAVBZERDyMMYb699fnh+4/MCV0Civ2raDkmJK8ueZNzsaddTqe\nZEAqCyIiHsrL5UWXoC7seWEPL1R8gbe/e5tSY0sx+8fZmsxJUpXKgoiIh7snyz28W+dddvbYScXC\nFWk7vy3Vp1fnhyM/OB1NMgiVBRGRNKJ47uLMaz2Ple1XcuLsCYInB9Pr61788fcfTkeTdE5lQUQk\njXm82ONsf247H9T5gJnbZ1JyTEmmRE7Bbd1OR5N0SmVBRCQN8vbypnfl3ux+YTeNSjai21fdqDqt\nKtuPbHc6mqRDKgsiImlYwewFmdF0Bms7reXk+ZMETw6mz7I+xJ6PdTqapCMqCyIi6UD1otXZ1n0b\nQ2oNYeLWiQSOC2RB9AKnY0k6obIgIpJOZPbKzKtVX2Vnz52UL1Se5p83p/mc5hyOPex0NEnjVBZE\nRNIZ/5z+LHpqEXNazmHjfzcSOC6QiVsnagCkJJnKgohIOmSMofWDrYnuGU2r0q14fsnz1Jheg13H\ndzkdTdKgRJcFY0w1Y8wiY8whY4zbGNP4NuvXuLLe1a9Lxpj8SY8tIiJ3Ile2XExtPJXVHVdz5PQR\nyk0qx/CNw7nkvuR0NElDknJmwRf4AegB3Ol8oxYoARS88ipkrT2WhGOLiEgShPiHsP257Txf4Xle\nWfEKVT+uSvTv0U7HkjQi0WXBWvuNtfYNa+1CIDEPWv/dWnvsn1dijysiInfHx9uHEfVG8F3n7/jj\n7z8oP6k87373rs4yyG2l1pgFA/xgjDlsjFlujHkslY4rIiLXeazIY/zQ/QdeqPgC/b7tR/Xp1dn7\nx16nY4kHS42y8BvQHWgBNAf+C6wxxpRLhWOLiMhNZPPOxvt132fdM+s4cvoIZSeWZdLWSXqapdyU\nuZt/GMYYN9DUWrsokdutAQ5Yazve4vMgILJ69erkyJHjms/atGlDmzZtkphYRESud/rCafou78uk\nyEnUv78+HzX+iH/5/cvpWHIHwsPDCQ8Pv2bZyZMnWbduHUCwtTYqOY7jVFl4D6hira1yi8+DgMjI\nyEiCgoKSnE9ERO7c13u+psuiLpy/dJ6PGn9E0weaOh1JkiAqKorg4GBIxrLg1DwL5bh8eUJERDxE\nwxIN+fH5H6letDrN5jSj++LunLlwxulY4gGSMs+CrzGm7FVjDopdeV/kyufvGGNmXLX+i8aYxsaY\n4saYB40xo4CawNhk+QpERCTZ5PHJw/zW85n0xCQ++c8nBE8OZttv225YT2MbMpaknFmoAGwDIrk8\nf8JwIAp488rnBYEiV62f+co6/wHWAA8Dj1tr1yQpsYiIpChjDN2CuxHVPQofbx8enfooIyNGcurU\nKcLCBhIQUJsiRZoSEFCbsLCBxMbqCZfp3V2NWUgpGrMgIuIZLly6QL+V/Rjx/Qj8Dufl9Kxx2DOt\nuHxHvMXlWkZg4AgiIubh5+fndFwhfY1ZEBGRNCCzV2aG1xvOE6eeIjbnBWy3l+G+DVc+Nbjd9YmO\n7s2AAcMdzSkpS2VBRERu66f5v8PEn+CvAOgUAtWGgrn8FEu3uz6LFm1IeAeSpqksiIhIgqy1xMX5\nwqkiMGMVrO8HtV6Hp5+AbCcAQ1ycjwY9pmMqCyIikiBjDN7eZwAL7kywejB8uhQKb4buQVB4E97e\nZzAmMY8LkrREZUFERG4rNLQKLtey/y3YWw8mboPYf0HnqhRtlVlnFtIxlQUREbmtoUP7Ehg4Apdr\nKZfvmgdO3YuZ0Y/ce//FWt+ltJ3fVpM4pVMqCyIiclt+fn5ERMyjV69N+PvXpXDhJvj71+WFHlHE\nTPyJOS3nsGjXIh6b9pieYJkOaZ4FERFJNGvtDWMUfjr2E83mNOP42eOEtwin/v31HUqXsWmeBRER\n8Qg3G8z4UP6H2NJ1C48VeYyGsxry9vq3NY4hnVBZEBGRZJMza04Wt1nM69Vf57VVr9FqbiuNY0gH\nVBZERCRZuYyLN2u+yYInF7Bs7zKqTKvCgb8OOB1L7oLKgoiIpIimDzRlY+eNnDp/ikemPML6A+ud\njiRJpLIgIiIp5uECD7O562YezP8gtWbWYnLkZKcjSRKoLIiISIrK65OX5e2W0y2oG92/6k7vb3pz\nyX3J6ViSCJmcDiAiIumft5c34xqNo3S+0rz4zYvs+WMP4S3C8cuix1qnBTqzICIiqaZnxZ4seXoJ\n6w+u18DHNERlQUREUlW9++uxsfNGYi/E8ujUR9n06yanI8ltqCyIiEiqezD/g2zuspniuYsTMiOE\n+dHznY4kCVBZEBERR+Tzzce3Hb6lSakmtPy8JSMjRmrGRw+lAY4iIuKYrJmyMrvFbPxz+tNneR/2\n/bmPUfVH4eXycjqaXEVlQUREHOUyLobVHkZAzgB6fN2DAycPEN4iHN/Mvk5Hkyt0GUJERDxC9wrd\nWdxmMav2r6LmjJocO3PM6UhyhcqCiIh4jIYlGrLumXUcOHmAKtOqsPePvU5HElQWRETEwwQVCiLi\n2QgMhsemPcbWw1udjpThqSyIiIjHKZarGBs6b8A/pz8h00P45pdvnI6UoaksiIiIR8rnm49VHVZR\nM6AmoeGhfLL9E6cjZVgqCyIi4rF8M/uy4MkFdCjTgQ5fdmBkxEinI2VIunVSREQ8WiZXJqY2nko+\n33z0Wd6HY2eO8fbjb2OMcTpahqGyICIiHs8Yw7Daw8jnk4++K/ry+9nfmfjERDK59GssNei7LCIi\nacbLj71MPt98dF7YmRN/nyC8RThZM2V1Ola6pzELIiKSpnQo24GFTy3km1++odHsRsSej3U6Urqn\nsiAiImlOo5KNWN5uOVsPb6X2J7U5cfaE05HSNZUFERFJk6oVrcbqjqvZ9+c+akyvweHYw05HSrdU\nFkREJM0KKhTE+mfWc/L8SapOq6rpoVOIyoKIiKRpD+R9gA2dN+Dt5U21j6ux8/edTkdKd1QWREQk\nzbsvx32s67SOvD55qf5xdaJ+i3I6UrqisiAiIulCgewFWNNpDcVzF6fmjJpsOLjB6UjphsqCiIik\nG7mz5WZl+5WUK1iOup/WZeW+lU5HShdUFkREJF3xy+LH0rZLqV60Oo1mN2LxrsVOR0rzVBZERCTd\n8fH2YeFTC2lUohHNP2/OFzu/cDpSmqayICIi6VJmr8zMaTmHlqVb8uQXTzLrP7OcjpRm6dkQIiKS\nbnl7efNps0/Jmikr7Re059zFczwb9KzTsdIclQUREUnXvFxefNT4I7JlykaXxV04d/EcPSv2dDpW\nmqKyICIi6Z7LuBjXcBxZM2Wl19JenL90nj6V+8R/bq3FGONgQs+msiAiIhmCMYbhdYeTxSsLLy9/\nmdizsZxY6Gbx4g3Exfni7X2G0NAqDB3aFz8/P6fjehSVBRERyTCMMbz9+NvYS5ZB3w3C/NgOG7MC\nMIBl3LhlrFrVgoiIeSoMV9HdECIikqEYYzi7JAtmVQdsyKdQ63XAAga3uz7R0b0ZMGC40zE9isqC\niIhkOIsXb8Cumw7L34fqQ6HO/3G5MIDbXZ9FizRV9NV0GUJERDIUay1xcb6AgY194VJmaPAiuC7C\nsuGAIS7OR4Mer6KyICIiGYoxBm/vM/xz6YFNYeD2gka9LheGpaPw9j6jonAVlQUREclwQkOrMG7c\nMtzu+pcXbOkJ7kwQ+hy49hFavIqzAT1MoscsGGOqGWMWGWMOGWPcxpjGd7BNiDEm0hhzzhiz2xjT\nMWlxRURE7t7QoX0JDByBy7WUf8YqENkNs+hFqLCE2BoHcVu3oxk9SVIGOPoCPwA9iP8O35oxxh/4\nCvgWKAuMBqYaY+ok4dgiIiJ3zc/Pj4iIefTqtQl//7oULtwEf/+6vFA1JxPqT2DGjzPosqiLCsMV\nib4MYa39BvgGwNzZBZ3ngX3W2n9feb/LGFMV6A2sSOzxRUREkoOfnx+jRw9i9OgbZ3DM7pOdjl92\nxG3dfNT4I7xcXg4mdV5qjFmoBKy8btkyYGQqHFtEROS2rv/bt12ZdngZL9otaMcle4npTaZn6MKQ\nGmWhIHD0umVHgXuMMVmstedTIYOIiEiitHm4DV4uL56e9zSX3JeY2WwmmVwZ874Aj/6qe/fuTY4c\nOa5Z1qZNG9q0aeNQIhERyUhaP9gaL+PFU/Oewm3dfNr8U48qDOHh4YSHh1+z7OTJk8l+HGPtbcco\n3npjY9xAU2vtogTWWQtEWmv7XLWsEzDSWpvrFtsEAZGRkZEEBQUlOZ+IiEhy+PLnL2k9tzVNHmjC\n7Oaz8fbydjrSLUVFRREcHAwQbK2NSo59psZ0zxHA49ctq3tluYiIiMdr+kBT5rWex8KfF/LkF09y\n4dIFpyOlqqTMs+BrjClrjCl3ZVGxK++LXPn8HWPMjKs2mXhlnXeNMaWMMT2AlsCIu04vIiKSSkJL\nhbLgyQUs2bOEVnNbcf5ixhlyl5QzCxWAbUAkl+dZGA5EAW9e+bwgUOSfla21MUAjoDaX52foDTxr\nrb3+DgkRERGP1qhkI7588kuW/bKMFp+3yDCFISnzLKwlgZJhrX3mJsvWAcGJPZaIiIinaVCiAYva\nLKLJZ01oNqcZ85+cT9ZMWZ2OlaL0iGoREZFEqlu8Ll+1+Yo1MWtoHN6Ys3FnnY6UolQWREREkuDx\nYo/zdduv2fjfjTwx+wnOXDjjdKQUo7IgIiKSRCH+IXzT7hu2HN5Cg1kNiD0f63SkFKGyICIicheq\n3leV5e2Ws/3odurPqs+p86ecjpTsVBZERETuUuUilVnRfgU7f99JnU/q8OfffzodKVmpLIiIiCSD\nioUr8m2Hb/nlj194fObjHD973OlIyUZlQUREJJkEFQpiTcc1/HrqV2rOqMnR09c/RzFtUlkQERFJ\nRg8XeJi1ndZy4uwJQmaEcDj2sNOR7prKgoiISDILzBfI2k5rOX3hNNU/rs6Bvw44HemuqCyIiIik\ngBJ5SrCu0zrc1k21j6ux58QepyMlmcqCiIhICgnIFcD6Z9bjm9mX6tOr89Oxn5yOlCQqCyIiIimo\n8D2FWdtpLQV8CxAyPYTIw5FOR0o0lQUREZEUlt83P6s7rub+3PdTa2YtNhzc4HSkRFFZEBERSQW5\nsuViRfsVBBUKou6ndVn2yzKnI90xlQUREZFU4pfFj6+f/ppaAbUIDQ9l7o65Tke6IyoLIiIiqSib\ndzbmt55P6wdb89S8p5gSOcXpSLeVyekAIiIiGY23lzczm80kZ9acdPuqG3+e+5N/V/m307FuSWVB\nRETEAS7jYkyDMeTOlpv/W/l/HD97nHdrv4sxxuloN1BZEBERcYgxhrdqvkXubLnpvaw3R88cZWro\nVLy9vJ2Odg2VBREREYe9VOklCvgWoOOXHTl25hhzW80le+bsTseKpwGOIiIiHqDNw21Y8vQSvjv4\nHbVm1OL3M787HSmeyoKIiIiHqFO8Dms6ruHAyQNUmVaFfX/uczoSoLIgIiLiUYL/FczGzhuxWCpN\nrcSmXzc5HUllQURExNMUz12ciGcjKJGnBCEzQpgfPd/RPCoLIiIiHiivT16+7fAtjUs1puXnLRkR\nMQJrrSNZdDeEiIiIh8qaKSvhLcIJyBnAy8tfZu8fexndYDSZXKn761tlQURExIO5jIthtYdRPFdx\nenzdg91/7GZOyznkzpY79TKk2pFEREQkyboGd2V5u+VE/RZFxSkV2fn7zlQ7tsqCiIhIGlEzoCZb\num4hm3c2Kk2txJLdS1LluCoLIiIiaUixXMXY2HkjNQNqEhoeyjvr38Ft3Sl6TJUFERGRNMYvix8L\nnlxA/2r96b+qP83mNOO/v/+XsLCBPPHEc8l+PJUFERGRNMhlXAypNYRFTy1iXcw6ir9XirFf5Oe3\n3yYk/7GSfY8iIiKSakJLhdL4SHviThXGdu4LpRYl+zFUFkRERNK4dQt3wrT/wH/aQ423kn3/Kgsi\nIiJpmLWWuDhfuJgNFk+Gle8k+zFUFkRERNIwYwze3meAK1NB76ub7MdQWRAREUnjQkOr4HItS7H9\nqyyIiIikcUOH9iUwcAQu11LizzAkI5UFERGRNM7Pz4+IiHn06rWJQoV6JPv+VRZERETSAT8/P0b/\nf3t3GyNXVcdx/PtLBRqa8CTahdDwVIPEaCEItRC21UYIGtCgiVFiSUwqBkgaXiiJMZHwxpRn8AEI\nhCIBNukrbSpYRWxMlJbYUjHYBiPFKGVLgWabFGrM9u+Lc7feTmfOzp3euTPs/j7JfTF3zjnzn5zc\n2f/eh/O//zbWr/c6C2ZmZtYwJwtmZmaW5WTBzMzMspwsmJmZWZaTBTMzM8tysmBmZmZZThbMzMws\ny8mCmZmZZTlZMDMzsywnC2ZmZpblZMEaMTY2NugQrEaez5nF82nT6SlZkHSTpJ2S3pe0SdLFmbZL\nJR1s2SYlfbT3sO2Dxj9GM4vnc2bxfNp0KicLkr4G3A38ELgQ+AuwQdKpmW4BfAwYKbbTIuKt6uGa\nmQxy/XcAAAW+SURBVJlZ03o5s3AL8HBEPBERO4DvAO8B35qm356IeGtq6+FzzczMbAAqJQuSjgEu\nAn43tS8iAngOWJLrCmyTtEvSbyRd2kuwZmZm1rwPVWx/KjAH2N2yfzdwXoc+bwI3AH8GjgNWAhsl\nXRIR2zr0mQuwffv2iuHZsJqYmGDr1q2DDsNq4vmcWTyfM0vpb+fcusZUOjHQZWPpNOANYElEbC7t\nXw2MRkTu7EJ5nI3APyPi+g7vfwN4quvAzMzMrNV1EfF0HQNVPbPwNjAJzG/ZPx8YrzDOi8Blmfc3\nANcBrwMHKoxrZmY2280FziL9La1FpWQhIv4raQuwHFgHIEnF6wcqDHUB6fJEp895B6glGzIzM5uF\n/lTnYFXPLADcAzxeJA0vkp6OOB54HEDSj4DTpy4xSFoF7AReIWU7K4HPAp8/2uDNzMys/yonCxGx\ntlhT4XbS5YdtwJURsadoMgIsKHU5lrQuw+mkRyxfBpZHxB+OJnAzMzNrRqUbHM3MzGz2cW0IMzMz\ny3KyYGZmZllDkSxI+r6kP0raL+ndCv1uL1aFfE/SbyUt7Gec1j1JJ0t6StKEpL2SHpU0b5o+a9oU\nHXumqZjt/6oUiyvaL5O0RdIBSa9KaruGig2Gi//NHJIul7RO0hvF3FzTRZ+jPj6HIlkAjgHWAg92\n20HSrcDNwLeBS4D9pIJWx/YlQqvqaeB80mO1XwRGgYe76Pcs6cbZqaJjX+9XgNZe1WJxks4C1pOW\ngV8E3A88KslPPA0BF/+bceaRHiy4kTRPWXUdn0N1g2OR7dwbEad00XYXcGdE3Fu8PoG07PT1EbG2\nv5FajqSPA38DLoqIl4p9VwK/As6IiLYLeElaA5wYEdc2FqwdQdImYHNErCpeC/gX8EBE3NGm/Wrg\nqoj4VGnfGGkuv9BQ2NZBD/O5FHgeODki9jUarFUi6SDw5YhYl2lTy/E5LGcWKpF0NinbLRe02gds\nJl/QypqxBNg7lSgUniNlwYun6btM0m5JOyT9TNK0iaPVp8dicZ8p3i/bkGlvDXHxP6Om4/MDmSyQ\nEoWgfUGrkebDsRYjwGGnLCNiEniX/Pw8C6wAPgd8D1gKPFP8J2TNyBWL6zR3Ix3anyDpuHrDs4p6\nmc+p4n9fAa4lnYXYKOmCfgVpfVXL8dnLCo5dKVZyvDXTJIDzI+LVfsVg9ep2Tnsdv+Xy0SuS/gr8\nA1gG/L7Xcc2se8Vvcvl3eZOkc0mr9frG1Vmqb8kCcBewZpo2r/U49jjpNNl8Ds+Y5gMvte1hdeh2\nTseBw+6cljQHOIUKBcciYqekt4GFOFloSi/F4sY7tN8XEf+pNzyrqKnifza8ajk++5YsFMWg3unT\n2DsljZPutH8ZDt3guBj4aT8+07qfU0kvACdJurB038JyUoK3uXPPI8Y5A/gwmaJjVq8ei8W9AFzV\nsu+KYr8NUFPF/2yo1XJ8DsU9C5IWSFoEnAnMkbSo2OaV2uyQ9KVSt/uAH0i6WtIngSeAfwO/bDR4\nO0JE7CDdQPOIpIslXQb8GBgrPwlRnlNJ8yTdIWmxpDMlLQd+QTodWluZVevKPcBKSSuKJ1seoqVY\nnKSfl9o/BJwjabWk8yTdCHy1GMcGr9J8Slol6RpJ50r6hKT7SMX/fjKA2K1F8Vu5qHQPyTnF6wXF\n+/05PiNi4Bvp1PZkm2201GYSWNHS7zZgF6lA1QZg4aC/i7dDc3MS8CQwAewFHgGOb2lzaE5JFUl/\nTTpldoB0OeNB4COD/i6zcSM9w/068D7pP5BPl95bAzzf0n4U2FK0/zvwzUF/B2+9zSfw3WIO9wN7\nSE9SjDYds7eOc7kUONjm7+Vj7eaz2HfUx+dQrbNgZmZmw2coLkOYmZnZ8HKyYGZmZllOFszMzCzL\nyYKZmZllOVkwMzOzLCcLZmZmluVkwczMzLKcLJiZmVmWkwUzMzPLcrJgZmZmWU4WzMzMLOt/vc6S\nDOWqHu0AAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<matplotlib.figure.Figure at 0x7fb2dc44aa90>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"import matplotlib.pyplot as plt\n", | |
"import numpy as np\n", | |
"from scipy.interpolate import lagrange\n", | |
"\n", | |
"def f(x):\n", | |
" return (1 + np.exp(-x)) / (1 + x ** 4)\n", | |
"x = np.linspace(-1, 1, 10)\n", | |
"x2 = np.linspace(-1, 1, 100)\n", | |
"y = f(x)\n", | |
"p = lagrange(x, y)\n", | |
"plt.plot(x, f(x), \"o\", label=\"Function\")\n", | |
"plt.plot(x2, np.polyval(p, x2), label=\"Polynom\")\n", | |
"plt.legend(loc='upper right');" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"下面是维基百科中关于拉格朗日多项式的说明\n", | |
"\n", | |
"给定 $k + 1$ 个数据点:\n", | |
"\n", | |
"$$(x_0, y_0),\\ldots,(x_j, y_j),\\ldots,(x_k, y_k)$$\n", | |
"\n", | |
"若$x_j$中没有重复值, 则拉格朗日插值的多项式是一个线性组合:\n", | |
"\n", | |
"$$L(x) := \\sum_{j=0}^{k} y_j \\ell_j(x)$$\n", | |
"\n", | |
"其中,\n", | |
"\n", | |
"$$\\ell_j(x) := \\prod_{\\begin{smallmatrix}0\\le m\\le k\\\\ m\\neq j\\end{smallmatrix}} \\frac{x-x_m}{x_j-x_m} = \\frac{(x-x_0)}{(x_j-x_0)} \\cdots \\frac{(x-x_{j-1})}{(x_j-x_{j-1})} \\frac{(x-x_{j+1})}{(x_j-x_{j+1})} \\cdots \\frac{(x-x_k)}{(x_j-x_k)},$$\n", | |
"\n", | |
"下面使用NumPy的广播运算功能实现上述公式。在运算中将会出现三维数组,其各个轴的含义分别为:\n", | |
"\n", | |
"* `i`: 进行插值运算的X轴坐标\n", | |
"* `j`: 与公式中$\\ell_j(x)$对应的轴\n", | |
"* `m`: 与公式中$x_m$对应的轴" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 38, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"xi = np.linspace(-1, 1, 100)[:, None, None]\n", | |
"xj = x[None, :, None]\n", | |
"yj = y[None, :, None]\n", | |
"xm = x[None, None, :]\n", | |
"ell = (xi - xm) / (xj - xm)\n", | |
"ell[~np.isfinite(ell)] = 1.0\n", | |
"ellj = ell.prod(axis=2, keepdims=True)\n", | |
"yi = (yj * ellj).sum(axis=1, keepdims=True)[:, 0, 0]" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"NumPy的广播功能虽然好用,但是需要我们记住每个轴对应的含义,并且在进行聚合运算需要设置`keepdims`参数为`True`。可以使用`xarray`简化该类运算。\n", | |
"\n", | |
"在下面的运算中,我们分别创建沿着`i`, `j`, `m`轴的`DataArray`对象。在进行除法运算之前先通过`where()`方法将分母数组的对象线上的元素设置为`NaN`。在进行累乘运算时,通过`skipna=True`跳过`NaN`元素。" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 65, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFkCAYAAACuFXjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8jvXjx/HX5545L+dTko3QKIetRA7N+dScUyIkh0JE\n9f199a10oDNSlPOpsihkiCFnlsOGijmEqcghiiGMfX5/bPmOLzebe7vube/n43H/seu+Du91te3t\nuj/X5zLWWkRERESux+V0ABEREfFuKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqC\niIiIuKWyICIiIm6pLIiIiIhbKgsiIiLiVorKgjHmaWPMNmPMyaTXemNM0xtsE2KMiTLGnDPG7DbG\ndL21yCIiIpKeUnpl4Vfg/4AgIBhYDswzxgRea2VjjD+wAPgOqAKMAiYaYxqlMq+IiIikM3OrD5Iy\nxhwHXrDWTrnGe+8Czay1lZMtCwPyWWub39KBRUREJF2kesyCMcZljHkMyA1EXme1GsCyq5ZFADVT\ne1wRERFJX9lSuoEx5h4Sy0FOIA5oY63deZ3ViwNHrlp2BLjNGJPDWnv+OscoBDQBYoFzKc0oIiKS\nheUE/IEIa+1xT+wwxWUB2Eni+IN8QHtgujGmrpvCkBpNgC88uD8REZGsphMwwxM7SnFZsNZeBPYl\nfbnFGFMdGAA8c43VDwPFrlpWDDh1vasKSWIBPv/8cwIDrzl2UjKYgQMHMnLkSKdjiIfofGYuOp+Z\nS0xMDJ07d4akv6WekJorC1dzATmu814k0OyqZY25/hiHf5wDCAwMJCgo6NbSiVfIly+fzmUmovOZ\nueh8Zloe+xg/RWXBGPMWsAj4BfAj8RLHQyQWAIwxbwO3W2v/mUthLNA36a6IyUADEj+60J0QIiIi\nGURKrywUBaYBJYCTwA9AY2vt8qT3iwOl/lnZWhtrjGkBjAT6A78BT1lrr75DQkRERLxUisqCtbbH\nDd5/8hrLVpM4gZOIiIhkQHo2hKSLjh07Oh1BPEjnM3PR+ZQbUVmQdKFfRpmLzmfmovMpN6KyICIi\nIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6p\nLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiI\niIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhb\nKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsi\nIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLi\nlsqCiIiIuKWyICIiIm6pLIiIiIhbKgsiIiLilsqCiIiIuKWyICIiIm6lqCwYYwYbYzYaY04ZY44Y\nY+YaY8rfYJuHjDEJV70uGWOK3lp0cYq11ukIIiKSjlJ6ZaEO8DHwANAQ8AWWGGNy3WA7C5QDiie9\nSlhrj6bw2OKguLg4+vcfQkBAQ0qVak1AQEP69x9CXFyc09FERCSNZUvJytba5sm/NsZ0A44CwcDa\nG2x+zFp7KkXpxCvExcVRs2Y7YmIGkZDwGmAAy5gxESxf3o7IyNn4+fk5nFJERNLKrY5ZyE/iVYMT\nN1jPAFuNMYeMMUuMMQ/e4nElHf3nPx8kFYWmJJ5KAENCQlNiYgby8svDnYwnIiJpLEVXFpIzxhjg\nQ2CttXaHm1V/B3oDm4EcQE9gpTGmurV2a2qPL55jreW3U78R80cMh+IOcfj0YY6cPsLhM4c58fcJ\nVp/fREKXleB6E1wXIcEXzuWD8/lIOO/H1IMR3LWhEP75/S+//HLoSoOISGaR6rIAfAJUBGq5W8la\nuxvYnWzR98aYssBAoKu7bQcOHEi+fPmuWNaxY0c6duyYqsCSWAx2H9/NytiVRP8ezU/HfuKnoz9x\n6vx/PyEqkLMAxfMWp1jeYhTKVQhzNjecuRMSsoH1AVc85DwJeQ9DoV2cyXuMF5e+yPlL5y/vo0ju\nIlQuVplqxatRrUQ1qhavSoVCFfBx+TjxbYuIZEphYWGEhYVdsezkyZMeP45Jzch2Y8xoIBSoY639\nJRXbvwfUstZes2gYY4KAqKioKIKCglKcT650KO4QC3cvZOWBlazYv4LfT/+Oj/HhnqL3XH5VKlKJ\nikUqcsdtd5AjW44rtg8IaEhs7FL++xFEchZ//0bs3beEI6ePsP+v/cT+Fcue43vYdmQbWw5vIfav\nWAD8svtRp3QdQkqHEOIfQrUS1cjmupW+KiIiV4uOjiY4OBgg2Fob7Yl9pvg3dVJRaAU8lJqikKQq\niR9PSBo5cvoIX+/4mpnbZ7L2l7UYYwguEcwTlZ+gXkA9apWqddMfFYSG1mLMmIikMQtXcrkW07Jl\nbVzGRQm/EpTwK8GDpa4ckvLn33+y9fBWvv/te1YeWMlrq17j7LKz+GX3o0GZBrSq0IoW5VpQJE8R\nj3zvIiLiWSm6smCM+QToCLTkyo8WTlprzyWt8xZQ0lrbNenrAcB+YDuQk8QxC32BRtbaldc5jq4s\npEL8pXjm7pzLuKhxrIxdicu4aFSmER0qdaBVhVYUyFUgVfv9790QA5MNcrS4XIsJDByZ4rshLly6\nQNShKJbvX87CPQv5/rfvAXiw1IO0rNCSDpU64J/fP1VZRUSyurS4spDSspBA4t0PV3vSWjs9aZ0p\nQGlrbf2kr18EegG3A2eBH4DXrbWr3RxHZSEFjpw+woToCYzdPJaDcQepc2cdulTpQpu721AodyGP\nHCMuLo6XXx5OePg64uNz4+t7lpYtazF06PO3fNvkkdNHWLhnIeG7wlmydwl/X/ybWqVq0eneTjxS\n6REK5y7ske9BRCQrcLwspBeVhZuz98Rehq4ZyowfZ+BjfOhcuTP9qvejcrHKaXpcay2JN8N43ukL\np/lm5zfM+HEGS/YuwRhD83LN6RnUk2Z3NdMASRGRG/CKMQvivAN/HeDN1W8ydetUiuYpylv136J7\nte6p/pghpdKqKADkzZ6XzpU707lyZ46eOcqs7bOYvGUyoWGh3HHbHfSo1oPu1bpTKl+pNMsgIiJX\n0pWFDOT3uN8ZunooE6InkD9nfgbXHszT9z1NLt8bzbad8UUdimJ81Hhm/DSDs/FnaVmhJc898Bx1\nS9dN0/IiIpLRpMWVBT11MgOIvxTPiMgRVBhdgbCfwnij3hvsG7CPgTUHZomiABB8ezDjQsdxaNAh\nPmn+CbuP7yZkWgjB44P5bNtnXLh0wemIIiKZlsqCl1sVu4pq46rx4tIX6VKlC3v77+Xftf9N3ux5\nnY7mCL8cfvS+rzc/PfMTEZ0jKJa3GF2+6YL/h/68t+69KyaXEhERz1BZ8FLHzhyj85zOhEwL4bYc\nt7G552ZGNx+dbuMSvJ0xhsZlG7Oo0yK299lOi3IteHn5y5T+sDSvrniVP87+4XREEZFMQ2XBC83b\nOY97Pr2HxT8vZnLLyaztvpZqJao5HctrVSxSkQktJ7BvwD6erPokwyOHU/rD0gyKGMTh04edjici\nkuGpLHiRU+dP0X1ed1rPbM0DJR9ge5/tPFntSVxGp+lm3HHbHYxoMoIDzx1gUI1BTNoyiTKjyvDi\nkhc5dubYdbfzxkG+IiLeRH+FvMTK2JVU/rQyX+/4msktJzPvsXkUy1vM6VgZUuHchXmz/pvEDojl\n+ZrPMzZqLAGjAhi8bDDHzx4HEieZ6t9/CAEBDSlVqjUBAQ3p338IcXFxDqcXEfE+unXSYQk2gWGr\nhzFk5RDqlq7L1NZTNdWxhx0/e5wP1n/Axxs/xsflw4D7BvD18+vY9dOLJCQ04b/TV0cQGDgixdNX\ni4h4E906mcmc+PsED894mCErh/BayGss77pcRSENFMpdiLcbvs2+AfvoVqUbw9YOI6bRFhKC9oPr\nYtJahoSEpsTEDOTll4c7mldExNuoLDhk86HNBI0LYsPBDSzqtIhXH3pVYxPSWNE8RRnVbBQlZteA\nfQ9Di77QtyJU/Jp/HnmSkNCU8PB1zgYVEfEy+uuUDq7+qGfylsnUmlyLonmKEt0rmiZ3NXEoWdZj\nrYU/C8Pc6TB2KxwvDx0ege614Y7vAUN8fG4NehQRSUbPhkgjcXFx/Oc/HzB//jri4/Pg63uGh0Mf\nJFuTOD7c/CG9gnrxUbOPyJEth9NRsxRjDL6+ZwALRyrDjIVQZhk0fgF61ISfHoWdf2gKaRGRZHRl\nIQ3ExcVRs2Y7xoypSWzsUg4enEfswXmMPrqcDzd9yNshbzP24bEqCg4JDa2FyxXx3wX7GsK4KPhm\nCty5jN/bbeDfy/5N3HndGSEiAioLaeI///mAmJhBJCQ0BQzkPQzd6kHZrZiZQ/h9zjn9y9VBw4a9\nQGDgCFyuRfwzVgHrwvVDMQKXVeFfNV/kow0fUX50eaZsmUKCTXA0r4iI01QW0sD8+euSbskDiuyA\nHg+A30GYvAa7c4gG0DnMz8+PyMjZ9Ou3AX//xpQs2Qp//8b067eBDWu/4e0mb7Or3y7q+deje3h3\nqk+ozrpfdM5EJOtSWfAway3x8XkAAyU3wpN14Fx+mLARDldDA+i8g5+fH6NGvcb+/Uv59ddv2L9/\nKaNGvXZ5foVS+Uoxo90M1j65FmMMtafU5om5T/B73O8OJxcRSX8qCx52eQBdwHfQpQH8cTdMXQVx\nJZPWsPj6ntHHEF7E3bmodWctNvTYwMTQiUT8HEH50eV5f937eiS2iGQpKgtpoGLbQtCpKfxSGz5b\nknhlIYnLtZiWLWs7mE5SymVcPBX0FLuf3U33qt0Z/N1gKn9amSV7lzgdTUQkXagseNi0rdNYnO9r\nbvu9EGbmMxCfO+kdi8u1iMDAkQwd+ryjGSV18ufMz6hmo9jSewvF8xanyedNaD+rPb+c/MXpaCIi\naUplwYMmRE2g27xuPFXtKQ4M38mzfaL/ZwCdnjuQ8d1b7F5WdF3BjLYzWP/regLHBPL2mrc5f/G8\n09FERNKEHiTlIROiJtBrQS/63t+Xj5t9fMXn4NZajVHIpE6dP8XrK19n1IZRlC1YljHNx9CwTEOn\nY4lIFqYHSXmpidETr1sUwP0AOsnYbstxG8ObDGfr01spnrc4jT5rxGNfP8ahuENORxMR8RiVhVs0\nKXoSPef3pM99fa5ZFCRruKfoPazsupLpraezInYFd4++m5GRI7mYcPHGG4uIeDmVhVswZcsUes7v\nyTP3PcPo5qNVFLI4YwxPVHmCXf120aVKF55f8jzB44NZ/+t6p6OJiNwSlYVUmrV9Fk+FP0Wv4F4q\nCnKF/DnzM7r5aDb13EQOnxzUmlyLnuE9OX72uNPRRERSRWUhFRbtWUTnOZ3pVLkTn7T4BJfRf0b5\nX8G3BxP5VCSfNP+Er3Z8RYXRFZi8ZbKeNSEiGY7+yqXQmgNraDerHc3KNWNyy8kqCuKWj8uHZ+5/\nhl39dtG8XHOeCn+KulPq8uORH52OJiJy0/SXLgWiDkXRYkYLatxRg5ntZ+Lr4+t0JMkgiuUtxvQ2\n01nRdQXH/z5OtXHVeHHJi5y+cNrpaCIiN6SycJN2/rGTpl80JbBIIPMem0fObDmdjiQZUIh/CNue\n3sab9d5k9KbRVBxTkbkxc/VgMRHxaioLN+FQ3CGafN6EonmKsqjTIvxyaAZGSb3sPtkZXGcwO/rs\noHKxyrSd1ZbQsFD2/7nf6WgiIteksnADp86fovkXzbmUcInFnRZTMFdBpyNJJhFQIID5Heczp8Mc\nth3ZRqVPKvHWmrf0REsR8ToqC25cuHSBtjPbEvtXLIs6LaJUvlJOR5JMxhhDm8A2xPSNoe/9fXl1\nxatUGVuFFftXOB1NROQylYXrSLAJdJ/XnTW/rOGbx77h3mL3Oh1JMrG82fPyfuP32dJ7C4VyFaL+\n9Pp0ntOZw6cPOx1NRERl4XoGLxvMjB9n8FmbzwjxD3E6jmQR9xa7l9VPrmZKqylE7I2gwugKjN44\nmksJl5yOJiJZmMrCNXy66VPeW/8eI5qMoEOlDk7HkSzGZVx0q9qNXf120fGejvRf1J/qE6uz4bcN\nTkcTkSxKZeEqi/Ysot+ifvSv3p/najzndBzJwgrmKsjYh8cS+VQk1lpqTqpJr/m9NG20iKQ7lYVk\nth3eRoevO9CiXAtGNBnhdBwRAB644wE29dzEx80+Ztb2WZQfXZ4JUROuOW205msQkbSgspDkUNwh\nHg57mHIFyzGj3Qx8XD5ORxK5zMflQ9/qfdnVbxeh5UPptaAXNSfVZPOhzcTFxdG//xACAhpSqlRr\nAgIa0r//EOLi4pyOLSKZhMoCcObCGULDQrHWMr/jfPJmz+t0JJFrKpa3GFNbT2XNk2s4d/Ec1SdU\np3TfQEZPrkRs7FIOHpxHbOxSxoypSc2a7VQYRMQjsnxZSLAJdJrTid3Hd7Pw8YWUvK2k05FEbqj2\nnbWJ6hVFnTNN+LPkX9i+z0DweDCXAENCQlNiYgby8svDnY4qIplAli8LL333EvN3z+fLdl9SpXgV\np+OI3LRsrmz88lU8fPwz7GoJoU9Dzwfgju8BSEhoSnj4OodTikhmkKXLwrSt03h33bt80OgDWpRv\n4XQckRSx1hIfnwfOFId5U2DiesBCj5rQuhvkPUJ8fG4NehSRW5Zly8K6X9bRa0Evnqr2lG6RlAzJ\nGIOv7xkgqQz8VhMmbITw8VBuITxbnrNVdxCfEO9oThHJ+LJkWYj9K5Y2M9tQ444afNLiE4wxTkcS\nSZXQ0Fq4XBH/XWB9ILonfLwbttbjr/v2cu+n9/Ltnm+dCykiGV6WKwtx5+MIDQvFL4cfszvMJrtP\ndqcjiaTasGEvEBg4ApdrEZevMGBxXfieSr/8zfou6ynpV5IWM1rQ/Ivm7Pxjp5NxRSSDylJl4VLC\nJTrN6cQvJ39hfsf5FM5d2OlIIrfEz8+PyMjZ9Ou3AX//xpQs2Qp//8b067eByMjZ1ChTg++6fMec\nDnPY+cdO7v30XgZFDOKvc385HV1EMhDjjYOfjDFBQFRUVBRBQUEe2+/gZYN5b/17LOi4gGblmnls\nvyLewlp73Y/Vzl08x4fff8jQ1UPJ5ZuL10Nep1dwL7K5sqVzShFJS9HR0QQHBwMEW2ujPbHPLHNl\n4YsfvuCdde/wXsP3VBQk03I3/iZntpz8u/a/2fPsHkLLh9Lv235UGVuFxT8vTseEIpIRZYmysOG3\nDTwV/hTdqnZjUM1BTscRcVQJvxJMbjWZzb02UyR3EZp90YxmXzRj+9HtTkcTES+V6cvCwVMHaTOz\nDUElghjbYqzufBBJElQiiBVdVzCnwxz2HN9D5bGV6T2/N4dPH3Y6moh4mRSVBWPMYGPMRmPMKWPM\nEWPMXGNM+ZvYLsQYE2WMOWeM2W2M6Zr6yDfv7/i/aT2zNT4uH+Y8Oocc2XKkx2FFMgxjDG0C27Cj\n7w6GNx7OVzu+otzH5Ri6eihn4886HU9EvERKryzUAT4GHgAaAr7AEmNMruttYIzxBxYA3wFVgFHA\nRGNMo1TkvWnWWrqHd2f70e3Me2wexfMWT8vDiWRo2X2y81yN59jbfy+9gnrxxqo3KPdxOSZFT+JS\nwiWn44mIw1JUFqy1za21n1lrY6y1PwLdgDuBYDebPQPss9b+y1q7y1o7BvgaGJja0DfjrTVv8eVP\nXzKt9TSCSnjujgqRzKxArgIMbzKcnf12Urd0XXrM70GVsVVYsHuBpo0WycJudcxCfhJngjnhZp0a\nwLKrlkUANW/x2Nc1N2YuL694mSEPDeGRSo+k1WFEMq0yBcoQ1i6MjT02UjRPUULDQgmZFsKG3zY4\nHU1EHJDqsmASRwp+CKy11u5ws2px4MhVy44AtxljPD6I4IcjP/DE3CdoX7E9rz70qqd3L5Kl3F/y\nfr7r8h3fPv4tf/79JzUm1aDtzLbEHItxOpqIpKNbmY3lE6AiUMtDWf7HwIEDyZcv3xXLOnbsSMeO\nHa+5/tEzR2kZ1pJyhcoxtdVUXCbT3+whkuaMMTQr14zGZRsz48cZvLLiFe759B6erPokQx4aQql8\npZyOKJJlhYWFERYWdsWykydPevw4qZrB0RgzGggF6lhrf7nBuquAKGvtoGTLugEjrbUFrrNNimdw\nPH/xPA0/a8ie43vY2HMjd+a78ya/GxFJifMXzzN281iGrhlK3Pk4+t7fl3/X/jdF8hRxOpqI4CUz\nOCYVhVZAvRsVhSSRQIOrljVOWu4R1lp6L+jNpoObmPvoXBUFkTSUI1sOBtQYwL7++xhcezAToidQ\n5qMyvLriVU6e8/y/aETEeSmdZ+EToBPwOHDGGFMs6ZUz2TpvGWOmJdtsLFDGGPOuMaaCMaYP0B4Y\n4YH8ALy77l2mbZvGpJaTqFkqzcZNikgyfjn8GBIyhP0D9vPMfc/w/vr3CRgVwLtr3+XMhTNOxxMR\nD0rplYWngduAlcChZK8OydYpAVz+ENNaGwu0IHFehq0k3jL5lLX26jskUmVuzFwGfzeYV+q+QqfK\nnTyxSxFJgUK5C/Feo/fY238vj93zGK+seIUyH5VhZORI/o7/2+l4IuIBGfqpk1t+30LtKbVpUa4F\nX7b/UgMaRbxA7F+xvLnqTaZtm0axvMV4qfZL9AjqoRlURdKJV4xZ8BaH4g4RGhZKxSIVmdpadz6I\neAv//P5MajWJnf120iCgAf0X96fcx+UYu3ks5y+edzqeiKRChvwLe/rCaULDQgEIfyyc3L65HU4k\nIle7q+BdTG8zne19tlOndB36LOyj0iCSQWW4snAx4SKPfv0oe47vYeHjCynhV8LpSCLixt2F7+aL\ntl+wvc92at9Z+3Jp+HTTpyoNIhlEhioL1lr6LuzLkr1L+LrD11QpXsXpSCJykwKLBDKj3Qx+6vMT\nte6sRd9v+1L2o7KM3jiacxfPXXMbbxxTJZIVZaiy8O66dxkfPZ7xD4+ncdnGTscRkVSoWKQiYe3C\n2NF3B/UD6jNg8QACRgUwMnIkZy6cIS4ujv79hxAQ0JBSpVoTENCQ/v2HEBcX53R0kSwrw9wNMePH\nGXSa04khDw3htZDXHM0nIp6z5/ge3l77NtO3TadAzgK4Nubn2MJ3sOfaAgawuFwRBAaOIDJyNn5+\nfk5HFvFqWfZuiOX7l/PkvCfpWqUrQx4a4nQcEfGgcoXKMbnVZPY8u4fif5bmaGAsdkAPCHkNcp0A\nDAkJTYmJGcjLLw93Oq5IluT1ZWHjwY20+rIVIf4hjA8dT+LDLkUkswkoEMDpL/PDqP2wtRvUeh+e\nKw2NXoS8v5OQ0JTw8HVOxxTJkm7lqZNprknHrpx+bC9V7qjMnA5zyO6T3elIIpJGrLXEx+eBuDsg\nYiSseQlqfAjVR8MDH8OW7vz9c+J6+keDSPry6isLf9x/mHNHi3Pyk1wknE9wOo6IpCFjDL6+Z4Ck\ncVRni8DyYTDyF1j1KlT8iiOPfkfXb7qy/eh2R7OKZDVeXRa4mBM+i2T3D//SZ5UiWUBoaC1crogr\nF57PB2tewowaT92/m7AydiX3fHoPrb5sxfe/fe9MUJEsxrvLwsKxcKaYPqsUySKGDXuBwMARuFyL\nuHyFAYvLtYiK5T5lwSuz+Ln/z0xpNYXdx3dTc1JN6k2rR8TPEZqTQSQNeXdZOP3P7IyG+Pjc+mUg\nksn5+fkRGTmbfv024O/fmJIlW+Hv35h+/TZcvm0yu092ulXtxvY+25ndYTanL5ym6RdNCRofxMyf\nZnIp4ZLT34ZIpuPV8yxAFBAEWPz9G7F/v0eeai0iGcTNDGa01rIidgXvrH2HpfuWUrZAWV548AW6\nVe1Gzmw50ympiPfIsvMsuFyLadmyttMxRCSd3cxdD8YY6gfUZ8kTS9jcczNBJYLos7AP/h/68/aa\nt/nr3F/pkFQkc/PyspD4WWVg4EiGDn3e6TAi4uWCbw9m1iOz2NVvF60qtOK1Va9x58g7eXHJixyK\nO+R0PJEMy6vLQokSfa74rFJE5GaUK1SOcaHjiB0QS5/7+zA+ejwBowLoEd6DXX/scjqeSIbj1WMW\nkj8bQkQktU6eO8m4qHGM/H4kR04foU1gG/6v1v9RvWR1p6OJeFyWHbMgInIr8uXMx79q/Yv9A/Yz\n7uFx/HjkRx6Y+AD1p9Vnyd4lutNK5AZUFkQky8iZLSc9g3sS0zeGrx75ilPnT9Hk8yYEjw9m1vZZ\nuu1S5DpUFkQky/Fx+dC+Yns29dzEsieWUSh3IR79+lEqjK7A+KjxnLt4zumIIl5FZUFEsixjDA3K\nNGDpE0vZ1HMT1UpU4+kFT1NmVBneX/c+cefjnI4o4hVUFkREgPtuv4+vHvmKmL4xNC/XnP8s/w93\nfngnryx/hT/O/uF0PBFHqSyIiCRToXAFJracyL4B+3iy6pOM+H4EpT8szaCIQRw8ddDpeCKOUFkQ\nEbmGO267gxFNRnDguQM8X/N5pmydQsCoAHrN78XeE3udjieSrlQWRETcKJy7MG/Ue4MDzx1gaP2h\nzNs1jwqjK9BlbhdijsU4HU8kXagsiIjchNty3Ma/av2L2AGxjGwykuX7l1Ppk0p0+KoD2w5vczqe\nSJpSWRARSYFcvrl49oFn2dt/L+MeHsfmQ5upOq4qrb9sTfTv7ifL0+RPklGpLIiIpEKObDnoGdyT\n3c/uZmqrqWw/tp3g8cGEhoWy+dDmy+vFxcXRv/8QAgIaUqpUawICGtK//xDi4nRbpmQcKgsiIrcg\nmysbXat2JaZvDJ+1+Yw9x/dw/4T7aTGjBav2rKJmzXaMGVOT2NilHDw4j9jYpYwZU5OaNdupMEiG\nobIgIuIB2VzZ6Fy5M9v7bOeLtl+w98ReQmaEsL1yHAlFiwEmaU1DQkJTYmIG8vLLw52MLHLTVBZE\nRDzIx+XD4/c+zvY+2ym8uiIUPg5PB8GjbaHoj5fXS0hoSnj4OgeTitw8lQURkTTgMi5y7LoLxuyA\nuVOh2DZbl0yMAAAZ+0lEQVR4pgq0exwK7gEM8fG5NehRMgSVBRGRNGCMwdf3DCT4wLauMHonLBgL\npVdDv0Bo+RQm/3GMMTfemYjDVBZERNJIaGgtXK6IxC8SfCGqF3z0Myz5ACrM5lD77xkUMUjPnhCv\np7IgIpJGhg17gcDAEbhci4Ckjxsu5sC1sQJ3RwTxUq3BTIyeSNmPyjJ09VDOXDjjaF6R61FZEBFJ\nI35+fkRGzqZfvw34+zemZMlW+Ps3pl+/DWxcO483G73JvgH76F61O2+ufpOyH5Xl002fEn8p3uno\nIlcw3ji4xhgTBERFRUURFBTkdBwREY+w1l53jELsX7EMWTmEz7Z9RvlC5Xmv0XuElg/VmAZJsejo\naIKDgwGCrbXupxW9SbqyICKSTtz94ffP78+01tOI7h1NqXylaPVlK0KmhbDp4KZ0TChybSoLIiJe\npGrxqizpvIRFnRZx4u8TVJ9Ync5zOvPbqd+cjiZZmMqCiIiXMcbQ9K6mbO29lQmhE1i6bynlPy7P\n6ytf52z8WafjSRaksiAi4qV8XD70COrBnmf38Gz1Z3lr7VtUGF2BGT/O0GROkq5UFkREvNxtOW7j\n3UbvsqPPDqqXrE6nOZ2oO7UuWw9vdTqaZBEqCyIiGUTZgmWZ3WE2y55YxvGzxwkeH0y/b/tx4u8T\nTkeTTE5lQUQkg2lQpgHbnt7GB40+YPq26ZT/uDwToiaQYBOcjiaZlMqCiEgG5Ovjy8CaA9n97G5a\nlG9BrwW9qD25NtsOb3M6mmRCKgsiIhlY8bzFmdZ6Gqu6reLk+ZMEjw9mUMQg4s7HOR1NMhGVBRGR\nTKBu6bps6b2FofWHMnbzWALHBDI3Zq7TsSSTUFkQEckksvtk59+1/82OvjuoVqIabWe1pe3MthyK\nO+R0NMngVBZERDIZ//z+hD8Wzsz2M1n/63oCxwQydvNYDYCUVFNZEBHJhIwxdKjUgZi+MTxS8RGe\nWfgMD019iF1/7HI6mmRAKS4Lxpg6xphwY8xBY0yCMablDdZ/KGm95K9LxpiiqY8tIiI3o0CuAkxs\nOZEVXVdw+PRhqo6ryvD1w7mUcMnpaJKBpObKQh5gK9AHuNn5Ri1QDiie9CphrT2aimOLiEgqhPiH\nsO3pbTxz3zO8uPRFak+pTcyxGKdjSQaR4rJgrV1srX3VWjsPSMmD1o9Za4/+80rpcUVE5Nbk9s3N\niCYjWNt9LSf+PkG1cdV4d+27usogN5ReYxYMsNUYc8gYs8QY82A6HVdERK7yYKkH2dp7K89Wf5bB\n3w2m7tS67D2x1+lY4sXSoyz8DvQG2gFtgV+BlcaYqulwbBERuYZcvrl4v/H7rH5yNYdPH6bK2CqM\n2zxOT7OUazK38j+GMSYBaG2tDU/hdiuBA9bartd5PwiIqlu3Lvny5bvivY4dO9KxY8dUJhYRkaud\nvnCaF5a8wLiocTS9qymTWk7idr/bnY4lNyEsLIywsLArlp08eZLVq1cDBFtroz1xHKfKwntALWtt\nreu8HwRERUVFERQUlOp8IiJy877d8y09wntw/tJ5JrWcROu7WzsdSVIhOjqa4OBg8GBZcGqehaok\nfjwhIiJeonm55vz4zI/ULV2XNjPb0Ht+b85cOON0LPECqZlnIY8xpkqyMQdlkr4ulfT+28aYacnW\nH2CMaWmMKWuMqWSM+RCoB4z2yHcgIiIeUyh3IeZ0mMO4h8fx2Q+fETw+mC2/b/mf9TS2IWtJzZWF\n+4AtQBSJ8ycMB6KB15PeLw6USrZ+9qR1fgBWAvcCDay1K1OVWERE0pQxhl7BvYjuHU1u39w8MPEB\nRkaO5NSpU/TvP4SAgIaUKtWagICG9O8/hLg4PeEys7ulMQtpRWMWRES8w4VLFxi8bDAjvh+B36HC\nnP5iDPbMIyTeEW9xuSIIDBxBZORs/Pz8nI4rZK4xCyIikgFk98nO8CbDefjUY8Tlv4Dt9TzcuS7p\nXUNCQlNiYgby8svDHc0paUtlQUREbuinOcdg7E/wVwB0C4E6w8AkPsUyIaEp4eHr3O9AMjSVBRER\ncctaS3x8HjhVCqYthzWDof4r8PjDkOs4YIiPz61Bj5mYyoKIiLhljMHX9wxgISEbrHgTPl8EJTdC\n7yAouQFf3zMYk5LHBUlGorIgIiI3FBpaC5cr4r8L9jaBsVsg7nboXpvSj2TXlYVMTGVBRERuaNiw\nFwgMHIHLtYjEu+aBU3dgpg2m4N7bWZVnEZ3mdNIkTpmUyoKIiNyQn58fkZGz6ddvA/7+jSlZshX+\n/o15tk80sWN/Ymb7mYTvCufByQ/qCZaZkOZZEBGRFLPW/s8YhZ+O/kSbmW344+wfhLULo+ldTR1K\nl7VpngUREfEK1xrMeE/Re9jUcxMPlnqQ5l805601b2kcQyahsiAiIh6TP2d+5neczyt1X+E/y//D\nI189onEMmYDKgoiIeJTLuHi93uvMfXQuEXsjqDW5Fgf+OuB0LLkFKgsiIpImWt/dmvXd13Pq/Cnu\nn3A/aw6scTqSpJLKgoiIpJl7i93Lxp4bqVS0EvWn12d81HinI0kqqCyIiEiaKpy7MEs6L6FXUC96\nL+jNwMUDuZRwyelYkgLZnA4gIiKZn6+PL2NajKFikYoMWDyAPSf2ENYuDL8ceqx1RqArCyIikm76\nVu/LwscXsuaXNRr4mIGoLIiISLpqclcT1ndfT9yFOB6Y+AAbftvgdCS5AZUFERFJd5WKVmJjj42U\nLViWkGkhzImZ43QkcUNlQUREHFEkTxG+6/IdrSq0ov2s9oyMHKkZH72UBjiKiIhjcmbLyYx2M/DP\n78+gJYPY9+c+Pmz6IT4uH6ejSTIqCyIi4iiXcfFOw3cIyB9An2/7cODkAcLahZEnex6no0kSfQwh\nIiJeofd9vZnfcT7L9y+n3rR6HD1z1OlIkkRlQUREvEbzcs1Z/eRqDpw8QK3Jtdh7Yq/TkQSVBRER\n8TJBJYKIfCoSg+HByQ+y+dBmpyNleSoLIiLidcoUKMO67uvwz+9PyNQQFv+82OlIWZrKgoiIeKUi\neYqwvMty6gXUIzQslM+2feZ0pCxLZUFERLxWnux5mPvoXLpU7kKXb7owMnKk05GyJN06KSIiXi2b\nKxsTW06kSJ4iDFoyiKNnjvJWg7cwxjgdLctQWRAREa9njOGdhu9QJHcRXlj6AsfOHmPsw2PJ5tKf\nsfSg/8oiIpJhPP/g8xTJU4Tu87pz/O/jhLULI2e2nE7HyvQ0ZkFERDKULlW6MO+xeSz+eTEtZrQg\n7nyc05EyPZUFERHJcFqUb8GSzkvYfGgzDT9ryPGzx52OlKmpLIiISIZUp3QdVnRdwb4/9/HQ1Ic4\nFHfI6UiZlsqCiIhkWEElgljz5BpOnj9J7cm1NT10GlFZEBGRDO3uwnezrvs6fH18qTOlDjuO7XA6\nUqajsiAiIhnenfnuZHW31RTOXZi6U+oS/Xu005EyFZUFERHJFIrlLcbKbispW7As9abVY90v65yO\nlGmoLIiISKZRMFdBlj2xjKrFq9L488Ys27fM6UiZgsqCiIhkKn45/FjUaRF1S9elxYwWzN813+lI\nGZ7KgoiIZDq5fXMz77F5tCjXgraz2vL1jq+djpShqSyIiEimlN0nOzPbz6R9xfY8+vWjfPHDF05H\nyrD0bAgREcm0fH18+bzN5+TMlpMn5j7BuYvneCroKadjZTgqCyIikqn5uHyY1HISubLlosf8Hpy7\neI6+1fs6HStDUVkQEZFMz2VcjGk+hpzZctJvUT/OXzrPoJqDLr9vrcUY42BC76ayICIiWYIxhuGN\nh5PDJwfPL3meuLNxHJ+XwPz564iPz4Ov7xlCQ2sxbNgL+Pn5OR3Xq6gsiIhIlmGM4a0Gb2EvWV5b\n+xrmx87Y2KWAASxjxkSwfHk7IiNnqzAko7shREQkSzHGcHZhDszyLtiQz6H+K4AFDAkJTYmJGcjL\nLw93OqZXUVkQEZEsZ/78ddjVU2HJ+1B3GDT6PxILAyQkNCU8XFNFJ6ePIUREJEux1hIfnwcwsP4F\nuJQdmg0A10WIGA4Y4uNza9BjMioLIiKSpRhj8PU9wz8fPbChPyT4QIt+iYVh0Yf4+p5RUUhGZUFE\nRLKc0NBajBkTQUJC08QFm/pCQjYIfRpc+wgtW8vZgF4mxWMWjDF1jDHhxpiDxpgEY0zLm9gmxBgT\nZYw5Z4zZbYzpmrq4IiIit27YsBcIDByBy7WIf8YqENULEz4A7ltI3EO/kGATHM3oTVIzwDEPsBXo\nw+X/wtdnjPEHFgDfAVWAUcBEY0yjVBxbRETklvn5+REZOZt+/Tbg79+YkiVb4e/fmGdr5+fTpp8y\n7cdp9AjvocKQJMUfQ1hrFwOLAczNfaDzDLDPWvuvpK93GWNqAwOBpSk9voiIiCf4+fkxatRrjBr1\nvzM45s2dl67fdCXBJjCp5SR8XD4OJnVeeoxZqAEsu2pZBDAyHY4tIiJyQ1f/27dz5c74GB86z+3M\nJXuJqa2mZunCkB5loThw5KplR4DbjDE5rLXn0yGDiIhIinS8tyM+Lh8en/04lxIuMb3NdLK5suZ9\nAV79XQ8cOJB8+fJdsaxjx4507NjRoUQiIpKVdKjUAR/jw2OzHyPBJvB528+9qjCEhYURFhZ2xbKT\nJ096/DjG2huOUbz+xsYkAK2tteFu1lkFRFlrByVb1g0Yaa0tcJ1tgoCoqKgogoKCUp1PRETEE77Z\n+Q0dvupAq7tbMaPtDHx9fJ2OdF3R0dEEBwcDBFtroz2xz/SY7jkSaHDVssZJy0VERLxe67tbM7vD\nbObtnMejXz/KhUsXnI6UrlIzz0IeY0wVY0zVpEVlkr4ulfT+28aYack2GZu0zrvGmArGmD5Ae2DE\nLacXERFJJ6EVQpn76FwW7lnII189wvmLWWfIXWquLNwHbAGiSJxnYTgQDbye9H5xoNQ/K1trY4EW\nQEMS52cYCDxlrb36DgkRERGv1qJ8C7559Bsifo6g3ax2WaYwpGaehVW4KRnW2ievsWw1EJzSY4mI\niHibZuWaEd4xnFZftqLNzDbMeXQOObPldDpWmtIjqkVERFKocdnGLOi4gJWxK2kZ1pKz8WedjpSm\nVBZERERSoUGZBnzb6VvW/7qeh2c8zJkLZ5yOlGZUFkRERFIpxD+ExZ0Xs+nQJpp90Yy483FOR0oT\nKgsiIiK3oPadtVnSeQnbjmyj6RdNOXX+lNORPE5lQURE5BbVLFWTpU8sZcexHTT6rBF//v2n05E8\nSmVBRETEA6qXrM53Xb7j5xM/02B6A/44+4fTkTxGZUFERMRDgkoEsbLrSn479Rv1ptXjyOmrn6OY\nMaksiIiIeNC9xe5lVbdVHD97nJBpIRyKO+R0pFumsiAiIuJhgUUCWdVtFacvnKbulLoc+OuA05Fu\nicqCiIhIGihXqByru60mwSZQZ0od9hzf43SkVFNZEBERSSMBBQJY8+Qa8mTPQ92pdfnp6E9OR0oV\nlQUREZE0VPK2kqzqtopieYoRMjWEqENRTkdKMZUFERGRNFY0T1FWdF3BXQXvov70+qz7ZZ3TkVJE\nZUFERCQdFMhVgKVPLCWoRBCNP29MxM8RTke6aSoLIiIi6cQvhx/fPv4t9QPqExoWylfbv3I60k1R\nWRAREUlHuXxzMafDHDpU6sBjsx9jQtQEpyPdUDanA4iIiGQ1vj6+TG8znfw589NrQS/+PPcn/6r1\nL6djXZfKgoiIiANcxsXHzT6mYK6C/N+y/+OPs3/wbsN3McY4He1/qCyIiIg4xBjDG/XeoGCuggyM\nGMiRM0eYGDoRXx9fp6NdQWVBRETEYc/VeI5ieYrR9ZuuHD1zlK8e+Yq82fM6HesyDXAUERHxAh3v\n7cjCxxey9pe11J9Wn2Nnjjkd6TKVBRERES/RqGwjVnZdyYGTB6g1uRb7/tzndCRAZUFERMSrBN8e\nzPru67FYakyswYbfNjgdSWVBRETE25QtWJbIpyIpV6gcIdNCmBMzx9E8KgsiIiJeqHDuwnzX5Tta\nVmhJ+1ntGRE5AmutI1l0N4SIiIiXypktJ2HtwgjIH8DzS55n74m9jGo2imyu9P3zrbIgIiLixVzG\nxTsN36FsgbL0+bYPu0/sZmb7mRTMVTD9MqTbkURERCTVegb3ZEnnJUT/Hk31CdXZcWxHuh1bZUFE\nRCSDqBdQj009N5HLNxc1JtZg4e6F6XJclQUREZEMpEyBMqzvvp56AfUIDQvl7TVvk2AT0vSYKgsi\nIiIZjF8OP+Y+OpeX6rzES8tfos3MNvx67Ff69x/Cww8/7fHjqSyIiIhkQC7jYmj9oYQ/Fs7q2NWU\nfa8Co78uyu+/f+r5Y3l8jyIiIpJuQiuE0vLwE8SfKont/gJUCPf4MVQWREREMrjV83bA5B/ghyfg\noTc8vn+VBRERkQzMWkt8fB64mAvmj4dlb3v8GCoLIiIiGZgxBl/fM0DSVND7Gnv8GCoLIiIiGVxo\naC1crog027/KgoiISAY3bNgLBAaOwOVaxOUrDB6ksiAiIpLB+fn5ERk5m379NlCiRB+P719lQURE\nJBPw8/Nj1KjXWLBA8yyIiIhIOlNZEBEREbdUFkRERMQtlQURERFxS2VBRERE3FJZEBEREbdUFkRE\nRMQtlQURERFxS2VBRERE3FJZEBEREbdUFiRdhIWFOR1BPEjnM3PR+ZQbSVVZMMb0NcbsN8b8bYz5\n3hhzv5t1HzLGJFz1umSMKZr62JLR6JdR5qLzmbnofMqNpLgsGGMeBYYDQ4BqwDYgwhhT2M1mFigH\nFE96lbDWHk15XBEREUlvqbmyMBAYZ62dbq3dCTwNnAW632C7Y9bao/+8UnFcERERcUCKyoIxxhcI\nBr77Z5m11gLLgJruNgW2GmMOGWOWGGMeTE1YERERSX/ZUrh+YcAHOHLV8iNAhets8zvQG9gM5AB6\nAiuNMdWttVuvs01OgJiYmBTGE2918uRJoqOjnY4hHqLzmbnofGYuyf525vTUPk3ihYGbXNmYEsBB\noKa1dkOy5e8Cda217q4uJN/PSuCAtbbrdd5/HPjipoOJiIjI1TpZa2d4YkcpvbLwB3AJKHbV8mLA\n4RTsZyNQy837EUAnIBY4l4L9ioiIZHU5AX8S/5Z6RIrKgrU23hgTBTQAwgGMMSbp649SsKuqJH48\ncb3jHAc80oZERESyoPWe3FlKrywAjACmJpWGjSTeHZEbmApgjHkbuP2fjxiMMQOA/cB2EttOT6Ae\n0OhWw4uIiEjaS3FZsNbOSppT4Q0SP37YCjSx1h5LWqU4UCrZJtlJnJfhdhJvsfwBaGCtXX0rwUVE\nRCR9pGiAo4iIiGQ9ejaEiIiIuKWyICIiIm55RVkwxrxkjFlnjDljjDmRgu3eSJoV8qwxZqkx5q60\nzCk3zxhTwBjzhTHmpDHmT2PMRGNMnhtsM+UaDx37Nr0yy3+l5GFxSeuHGGOijDHnjDG7jTHXnENF\nnKGH/2Uexpg6xphwY8zBpHPT8ia2ueWfT68oC4AvMAv49GY3MMb8H9AP6AVUB86Q+ECr7GmSUFJq\nBhBI4m21LYC6wLib2G4RiQNn/3noWMe0CijXltKHxRlj/IEFJE4DXwUYBUw0xuiOJy+gh/9lOnlI\nvLGgD4nnyS1P/Xx61QDHpLYz0lpb8CbWPQS8b60dmfT1bSROO93VWjsrbZOKO8aYu4EdQLC1dkvS\nsibAQuAOa+01J/AyxkwB8llr26ZbWPkfxpjvgQ3W2gFJXxvgV+Aja+1711j/XaCZtbZysmVhJJ7L\n5ukUW64jFefzIWA5UMBaeypdw0qKGGMSgNbW2nA363jk59NbriykiDEmgMS2m/yBVqeADbh/oJWk\nj5rAn/8UhSTLSGzBD9xg2xBjzBFjzE5jzCfGmBsWR/GcVD4srkbS+8lFuFlf0oke/id46OczQ5YF\nEouC5doPtCqe/nHkKsWBKy5ZWmsvASdwf34WAV2A+sC/gIeAb5P+JSTpw93D4q537opfZ/3bjDE5\nPBtPUig15/Ofh/+1A9qSeBVipTGmalqFlDTlkZ/P1MzgeFOSZnL8PzerWCDQWrs7rTKIZ93sOU3t\n/q/6+Gi7MeZHYC8QAqxI7X5F5OYl/U5O/nv5e2NMWRJn69XA1SwqzcoC8AEw5Qbr7Evlvg+TeJms\nGFc2pmLAlmtuIZ5ws+f0MHDFyGljjA9QkBQ8cMxau98Y8wdwFyoL6SU1D4s7fJ31T1lrz3s2nqRQ\nej38T7yXR34+06wsJD0M6nga7Xu/MeYwiSPtf4DLAxwfAMakxTHl5s+pMSYSyG+MqZZs3EIDEgve\nhutv+T/7uQMohJuHjolnpfJhcZFAs6uWNU5aLg5Kr4f/iVfzyM+nV4xZMMaUMsZUAUoDPsaYKkmv\nPMnW2WmMaZVssw+Bl40xocaYe4HpwG/AvHQNL//DWruTxAE0E4wx9xtjagEfA2HJ74RIfk6NMXmM\nMe8ZYx4wxpQ2xjQAviHxcqjHHrMqN2UE0NMY0yXpzpaxXPWwOGPMtGTrjwXKGGPeNcZUMMb0Adon\n7Uecl6LzaYwZYIxpaYwpa4ypZIz5kMSH/412ILtcJel3ZZVkY0jKJH1dKun9tPn5tNY6/iLx0val\na7zqJlvnEtDlqu1eAw6R+ICqCOAup78XvS6fm/zA58BJ4E9gApD7qnUun1MSn0i6mMRLZudI/Djj\nU6CI099LVnyReA93LPA3if8CuS/Ze1OA5VetXxeISlp/D/CE09+DXqk7n8CLSefwDHCMxDsp6qZ3\nZr2uey4fAhKu8fdy8rXOZ9KyW/759Kp5FkRERMT7eMXHECIiIuK9VBZERETELZUFERERcUtlQURE\nRNxSWRARERG3VBZERETELZUFERERcUtlQURERNxSWRARERG3VBZERETELZUFERERcev/AfiXzVfq\nL8KLAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<matplotlib.figure.Figure at 0x7fb2cea01b00>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"import xarray as xr\n", | |
"\n", | |
"xi = xr.DataArray(np.linspace(-1, 1, 100), dims=\"i\")\n", | |
"xm = xr.DataArray(x, dims=\"m\")\n", | |
"xj = xr.DataArray(x, dims=\"j\")\n", | |
"yj = xr.DataArray(y, dims=\"j\")\n", | |
"\n", | |
"ell = (xi - xm) / (xj - xm).where(xj.j != xm.m)\n", | |
"ellj = ell.prod(dim=\"m\", skipna=True)\n", | |
"yi = (ellj * yj).sum(dim=\"j\")\n", | |
"\n", | |
"plt.plot(x, y, \"o\")\n", | |
"plt.plot(xi, yi);" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": { | |
"collapsed": true | |
}, | |
"outputs": [], | |
"source": [ | |
"den = x[:, None] - x[None, :]\n", | |
"num = x2[:, None] - x[None, :]\n", | |
"\n", | |
"with np.errstate(divide='ignore', invalid='ignore'):\n", | |
" div = num[:, None, :] / den[None, :, :]\n", | |
"\n", | |
"div[~np.isfinite(div)] = 1\n", | |
"\n", | |
"y2 = (np.product(div, axis=-1) * y).sum(axis=1)\n", | |
"\n", | |
"plt.plot(x, y, \"o\")\n", | |
"plt.plot(x2, y2);" | |
] | |
} | |
], | |
"metadata": { | |
"anaconda-cloud": {}, | |
"kernelspec": { | |
"display_name": "Python [default]", | |
"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.5.2" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 0 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment