diff --git a/tutorial/ex1_simple_ZARC_model.ipynb b/tutorial/ex1_simple_ZARC_model.ipynb index e602278..ddc4146 100644 --- a/tutorial/ex1_simple_ZARC_model.ipynb +++ b/tutorial/ex1_simple_ZARC_model.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Gaussian Process Distribution of Relaxation Times. " + "# Gaussian Process Distribution of Relaxation Times" ] }, { @@ -65,21 +65,7 @@ "\\mathbf \\Sigma_{\\gamma| Z^{\\rm exp}_{\\rm im}} &= \\mathbf K- \\mathcal L_{\\rm im} \\mathbf K \\left(\\mathcal L^2_{\\rm im} \\mathbf K + \\sigma_n^2 \\mathbf I \\right)^{-1}\\mathcal L_{\\rm im} \\mathbf K^\\top\n", "\\end{align}$$\n", "\n", - "The above formulas depend on 1) the kernel, $k(\\xi, \\xi^\\prime)$; 2) the noise level, $\\sigma_n$; and 3) the experimental data, $\\mathbf Z^{\\rm exp}_{\\rm im}$ (at the log-frequencies $\\mathbf \\xi$). \n", - "\n", - "To optimize these hyperparameters, we maximize the marginal likelihood, i.e. the probability, $p(\\mathbf Z^{\\rm exp}_{\\rm im}|\\mathbf \\xi, \\mathbf \\theta)$ of measuring the data, $\\mathbf Z^{\\rm exp}_{\\rm im}$. The maximizing $\\mathbf \\theta$ is the one that would have most likely resulted in the measured experimental data. We note that $\\mathbf Z^{\\rm exp}_{\\rm im}\\vert \\mathbf{\\theta}\\sim \\mathcal N\\left(0, \\mathcal L^2_{\\rm im} + \\sigma_n^2 \\mathbf I \\right)$. Therefore, we obtain the following marginal log-likelihood (MLL)\n", - "\n", - "$$\n", - "\\log p(\\mathbf Z^{\\rm exp}_{\\rm im}|\\mathbf \\xi, \\mathbf \\theta)= - \\frac{1}{2} {\\mathbf Z^{\\rm exp}_{\\rm im}}^\\top \\left(\\mathcal L^2_{\\rm im} + \\sigma_n^2 \\mathbf I \\right)^{-1} \\mathbf Z^{\\rm exp}_{\\rm im} -\\frac{1}{2} \\log \\left| \\mathcal L^2_{\\rm im} + \\sigma_n^2 \\mathbf I \\right| - \\frac{N}{2} \\log 2\\pi$$\n", - "\n", - "In practice, we minimize another function $L(\\mathbf \\theta)$ called negative (and shifted) MLL (NMLL):\n", - "\n", - "$$L(\\mathbf \\theta) = - \\log p(\\mathbf Z^{\\rm exp}_{\\rm im}|\\mathbf \\xi, \\mathbf \\theta) - \\frac{N}{2} \\log 2\\pi$$\n", - "\n", - "and the optimial hyperparameters $\\mathbf \\theta$ are obtained by solving the following problem with gradient based approaches\n", - "$$\\mathbf \\theta = \\arg\\min_{\\mathbf \\theta^\\prime}L(\\mathbf \\theta^\\prime)$$\n", - "\n", - "In the developed GP-DRT model we solve this minimization problem using the `optimize` approach provided by `scipy` package" + "The above formulas depend on 1) the kernel, $k(\\xi, \\xi^\\prime)$; 2) the noise level, $\\sigma_n$; and 3) the experimental data, $\\mathbf Z^{\\rm exp}_{\\rm im}$ (at the log-frequencies $\\mathbf \\xi$). " ] }, { @@ -129,6 +115,7 @@ "tau = 1/freq_vec\n", "\n", "# define the frequency range used for prediction\n", + "# note: we could have used other values\n", "freq_vec_star = np.logspace(-4., 4., num=81, endpoint=True)\n", "xi_vec_star = np.log(freq_vec_star)\n", "\n", @@ -142,12 +129,12 @@ "Z_exact = R_inf+1./(1./R_ct+C*(1j*2.*pi*freq_vec)**phi)\n", "gamma_fct = (R_ct)/(2.*pi)*sin((1.-phi)*pi)/(np.cosh(phi*np.log(tau/tau_0))-cos((1.-phi)*pi))\n", "\n", - "# finer mesh for plotting only\n", + "# we will use a finer mesh for plotting the results\n", "freq_vec_plot = np.logspace(-4., 4., num=10*(N_freqs-1), endpoint=True)\n", "tau_plot = 1/freq_vec_plot\n", "gamma_fct_plot = (R_ct)/(2.*pi)*sin((1.-phi)*pi)/(np.cosh(phi*np.log(tau_plot/tau_0))-cos((1.-phi)*pi)) # for plotting only\n", "\n", - "# random and adding noise to the analytically computed impedance\n", + "# we will add noise to the impedance computed analytically\n", "rng = np.random.seed(214975)\n", "sigma_n_exp = 1.\n", "Z_exp = Z_exact + sigma_n_exp*(np.random.normal(0, 1, N_freqs)+1j*np.random.normal(0, 1, N_freqs))" @@ -157,7 +144,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### show the synthetic impedance in the Nyquist plot, as shown in Figure 7 (a)" + "## 2) Show the synthetic impedance in the Nyquist plot - this is similar to Figure 7 (a)" ] }, { @@ -216,7 +203,34 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2) running the GP-DRT model to get the optimal hyperparameters" + "## 3) Obtain the optimal hyperparameters of the GP-DRT model by minimizing the negative marginal log likelihood (NMLL)\n", + "\n", + "We constrain the kernel to be a squared exponential, _i.e._\n", + "\n", + "$$\n", + "k(\\xi, \\xi^\\prime) = \\sigma_f^2 \\exp\\left(-\\frac{1}{2 \\ell^2}\\left(\\xi-\\xi^\\prime\\right)^2 \\right)\n", + "$$\n", + "\n", + "and modify its two parameters, $\\sigma_f$ and $\\ell$ as well as the noise level $\\sigma_n$. Therefore, the vector of hyperparameters of the GP-DRT is assumed to be $\\boldsymbol \\theta = \\begin{pmatrix} \\sigma_n, \\sigma_f, \\ell \\end{pmatrix}^\\top$.\n", + "\n", + "Following the same derivation from the article we can write that\n", + "\n", + "$$\n", + "\\log p(\\mathbf Z^{\\rm exp}_{\\rm im}|\\boldsymbol \\xi, \\boldsymbol \\theta)= - \\frac{1}{2} {\\mathbf Z^{\\rm exp}_{\\rm im}}^\\top \\left(\\mathcal L^2_{\\rm im} \\mathbf K +\\sigma_n^2\\mathbf I \\right)^{-1} \\mathbf Z^{\\rm exp}_{\\rm im} -\\frac{1}{2} \\log \\left| \\mathcal L^2_{\\rm im} \\mathbf K+\\sigma_n^2\\mathbf I \\right| - \\frac{N}{2} \\log 2\\pi\n", + "$$\n", + "\n", + "We will call $L(\\boldsymbol \\theta)$ the negative (and shifted) MLL (NMLL):\n", + "$$\n", + "L(\\boldsymbol \\theta) = - \\log p(\\mathbf Z^{\\rm exp}_{\\rm im}|\\boldsymbol \\xi, \\boldsymbol \\theta) - \\frac{N}{2} \\log 2\\pi\n", + "$$\n", + "\n", + "the experimental evidence is maximized for\n", + "\n", + "$$\n", + "\\boldsymbol \\theta = \\arg \\min_{\\boldsymbol \\theta^\\prime}L(\\boldsymbol \\theta^\\prime)\n", + "$$\n", + "\n", + "The above minimization problem is solved using the `optimize` method provided by `scipy` package" ] }, { @@ -228,16 +242,16 @@ "name": "stdout", "output_type": "stream", "text": [ - "sigma_n, sigma_f, ell\n", - "0.8903599197476706 5.001415082296882 1.0120874912591977\n", - "0.8136354402155858 5.003533687387537 1.029191223973519\n", - "0.8291862644351911 5.035786747187512 1.2588675436350178\n", - "0.8303934003754703 5.083236934181329 1.2117785066002995\n", - "0.8304435946349606 5.206077455873398 1.2283658609053703\n", - "0.8305174795222053 5.387448882835493 1.2524148478215245\n", - "0.8305140530808819 5.406904963014705 1.2546655487001486\n", - "0.8305265249777404 5.407074309971135 1.2546848761160714\n", - "0.8305262553195482 5.407086124627953 1.2546867263395471\n", + "sigma_n, sigma_f, ell\n", + "0.8903599 5.0014151 1.0120875\n", + "0.8136354 5.0035337 1.0291912\n", + "0.8291863 5.0357867 1.2588675\n", + "0.8303934 5.0832369 1.2117785\n", + "0.8304436 5.2060775 1.2283659\n", + "0.8305175 5.3874489 1.2524148\n", + "0.8305141 5.4069050 1.2546655\n", + "0.8305265 5.4070743 1.2546849\n", + "0.8305263 5.4070861 1.2546867\n", "Optimization terminated successfully.\n", " Current function value: 53.657989\n", " Iterations: 9\n", @@ -258,15 +272,33 @@ "def print_results(theta):\n", " global seq_theta\n", " seq_theta = np.vstack((seq_theta, theta))\n", - " print(theta[0], theta[1], theta[2])\n", + " print('{0:.7f} {1:.7f} {2:.7f}'.format(theta[0], theta[1], theta[2]))\n", " \n", "GP_DRT.NMLL_fct(theta_0, Z_exp, xi_vec)\n", "GP_DRT.grad_NMLL_fct(theta_0, Z_exp, xi_vec)\n", - "print('sigma_n, sigma_f, ell')\n", + "print('sigma_n, sigma_f, ell')\n", "\n", "# minimize the NMLL L(\\theta) w.r.t sigma_n, sigma_f, ell using the Newton-CG method as implemented in scipy\n", "res = minimize(GP_DRT.NMLL_fct, theta_0, args=(Z_exp, xi_vec), method='Newton-CG', \\\n", - " jac=GP_DRT.grad_NMLL_fct, callback=print_results, options={'disp': True})" + " jac=GP_DRT.grad_NMLL_fct, callback=print_results, options={'disp': True})\n", + "\n", + "# collect the optimized parameters\n", + "sigma_n, sigma_f, ell = res.x" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4) Core of the GP-DRT" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4a) Compute matrices\n", + "Once we have identified the optimized parameters we can compute $\\mathbf K$, $\\mathcal L_{\\rm im} \\mathbf K$, and $\\mathcal L^2_{\\rm im} \\mathbf K$, which are given in equation `(18)` in the article" ] }, { @@ -275,18 +307,42 @@ "metadata": {}, "outputs": [], "source": [ - "# collect the optimized parameters\n", - "sigma_n, sigma_f, ell = res.x\n", - "\n", - "# calculate the matrices shown in eq (18)\n", "K = GP_DRT.matrix_K(xi_vec, xi_vec, sigma_f, ell)\n", "L_im_K = GP_DRT.matrix_L_im_K(xi_vec, xi_vec, sigma_f, ell)\n", "L2_im_K = GP_DRT.matrix_L2_im_K(xi_vec, xi_vec, sigma_f, ell)\n", "Sigma = (sigma_n**2)*np.eye(N_freqs)\n", "\n", "# the matrix $\\mathcal L^2_{\\rm im} \\mathbf K + \\sigma_n^2 \\mathbf I$ whose inverse is needed\n", - "K_im_full = L2_im_K + Sigma\n", + "K_im_full = L2_im_K + Sigma" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4b) Factorize the matrices and solve the linear equations\n", + "We are computing\n", + "$$\n", + "\\boldsymbol{\\gamma}|\\mathbf Z^{\\rm exp}_{\\rm im}\\sim \\mathcal N\\left( \\boldsymbol \\mu_{\\gamma|Z^{\\rm exp}_{\\rm im}}, \\boldsymbol \\Sigma_{\\gamma| Z^{\\rm exp}_{\\rm im}}\\right)\n", + "$$\n", + "\n", + "using \n", + "$$\n", + "\\begin{align}\n", + "\\boldsymbol \\mu_{\\gamma|Z^{\\rm exp}_{\\rm im}} &= \\mathcal L_{\\rm im} \\mathbf K\\left(\\mathcal L^2_{\\rm im} \\mathbf K+\\sigma_n^2\\mathbf I\\right)^{-1}\\mathbf Z^{\\rm exp}_{\\rm im} \\\\\n", + "\\boldsymbol \\Sigma_{\\gamma| Z^{\\rm exp}_{\\rm im}} &= \\mathbf K-\\mathcal L_{\\rm im} \\mathbf K\\left(\\mathcal L^2_{\\rm im} \\mathbf K+\\sigma_n^2\\mathbf I\\right)^{-1}\\mathcal L_{\\rm im} \\mathbf K^\\top\n", + "\\end{align}\n", + "$$\n", "\n", + "The key ingredient is to do Cholesky factorization of $\\mathcal L^2_{\\rm im} \\mathbf K+\\sigma_n^2\\mathbf I$" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ "# Cholesky factorization, L is a lower-triangular matrix\n", "L = np.linalg.cholesky(K_im_full)\n", "\n", @@ -311,12 +367,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### plot the DRT together with the analytical DRT" + "### 4c) Plot the obtained DRT against the analytical DRT" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -352,12 +408,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 3) calculate the impedance" + "### 4d) Predict the $\\gamma$ and the imaginary part of the GP-DRT impedance\n", + "\n", + "This part is explained in Section `2.3.3` of the main article " ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -393,27 +451,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 4) plot the imaginary part of impedance together with the exact one" + "### 4e) Plot the imaginary part of the GP-DRT impedance together with the exact one and the synthetic experiment" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEZCAYAAACq1zMoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3zbdb348dcnTXpJL3QdIMptpAiIOKHtUI8HL6xDp4IC6bgNAWEtbEc9ON2Y53i8OzrPRD0TaRkqMtHRAP4EmdoOUVSEtQXGRZClgwETdmmz9bKtt/fvj2+SJWmSpmnSJs37+Xh8H22/+X6TT79Nv+98bu+PERGUUkqpVLNNdwGUUkrNTBpglFJKpYUGGKWUUmmhAUYppVRaaIBRSimVFhpglFJKpYUGGKWUUmlhn+4CRDLGVAG1/h/nAU0i0uZ/bAUwG9gIVAALRGTltBRUKaVUXBkXYIBaEVkDYIwpB7YbY+aLSKf/8Xr/1gYsmaYyKqWUGkdGNZH5ay+rAj+LiA9o53CNxicis/xbnf9xpZRSGSijAoy/llIXsdsFhAUSY0yVMcY1ZQVTSik1YSaTc5H5g0gHcJKI+Iwx9UA3VvNYLTAvVh+M/9h6gOLi4urTTjttikqtlFIzQ0dHxx4ROSrZ8zM9wLQCK0P6XyIf9wINgUEAsdTU1Eh7e3s6iqiUUjOWMaZDRGqSPT+jmshC+UeMhQUXfx9NqE5gwZQWTCmlVEIycRQZxhg30BYILiGBZTMwK+TQcsA7xcVTSimVgIyrwRhjarFGiwWCiwuo8f8cOSzZBdwzxUVUSimVgIyqwfiDSav/+9CHqv1fu/xNZz7/Ph2qrJRSGSqjAoyIdAEmzuOdWP0uSimlMlzGNZEppZSaGTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFIqLTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFJToLm5ebqLMOU0wCilVJr5fD58vtzLaqUBRik1acaEb7E0N4cfV18f+9jq6vBjOzpSX+6psnJl1HURZ7yMykWmlFLp1tXVhcfjweVysWXLFlatWkV5eTkej4fVq1fj8/no6Oigu7ubyspK3G43jY2NuFwuOjs76erqoru7m46ODhobGykvLw8+b1NTE/PmzaO8vJyKigqqqqpoa2ujq6uLrq4uysvLcblc1NbWTvNVmCIiMuO36upqUUqlD4RvsTQ1hR+3ZEnsY6uqwo9tb09NWV0uV/B7r9crtbW1wZ97enrE5XJJT0+P9PT0SFNT05hzOzo6RESkpaVF6uvrwx7r6ekREZGOjg6pqqoKPtbY2CiNjY2p+QWmENAuk7j3ag1GKZUzmpubqao6vDCuy+UidDn18vJyGhsbqauro66ujvqINryOjo5gjcXlctHV1QWAx+OhvLw8+FhVVRWbN29O96+T8TTAKKUmTSSx4+rr4/e7hEpHn4vX68Xn89HW1hbc19LSEnaM2+2mqakp5nOsXLmSyspKfD4f3d3dgNU8VlFREXZcINjkMg0wSqmcMW/ePDo7O+P2gXR2drJy5UoaGhqora3F5XIB1kiw6urqYC2ms7OTjRs3AnDKKacEv49l7969gFXbcbvdKfqNMpuOIlNK5Qy32013d3fYkGGPxxP83ufz0d7eTm1tLU1NTdTV1QUfa29vD2sGCzSPdXV1ccIJJ4Tti3xel8ulw5SVUmqma2lpYfXq1Xg8HjweT7BPZs2aNZx00kl4vV4AKioq6OzspK6uLljrqampobm5mba2NqqqqqipqQl+v3nzZhobG8c8LxwObM3NzTnVdGYk0cbTLFZTUyOhHXlKKaXGZ4zpEJGaZM/XGoxSSqm00ACjlFIqLTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFIqLTJuJr8xpgoITLOdBzSJSJv/sXKgHugCXECbiHROS0GVUkrFlXEBBqgVkTUQDCjbjTHz/YGkBWgQkS7/463GmDoRyb0pskopleEyqonMX3tZFfjZHzjagVp/sHEFgotfF4drO0oppTJIRgUYfy2lLmK3C/ABNf6voXzAgikomlJKTUhbWxvV1dU5u5olZFiAAQj0twAYY1xABXAPUA50Rxy+1//4GMaYemNMuzGmfffu3ekqrlIqCV6vl6VLl1JWVobNZqOsrIylS5cG84Blo+bm5rCfa2traWhomKbSZIaMCzARmoD5IX0sUYNJNCLSLCI1IlJz1FFHpad0SqkJ27RpE3PnzmX9+vX09vYiIvT29rJ+/Xrmzp3Lpk2bpruIE+bz+XIyW/J4MjbAGGNWACtDRon5sGoxoWYztlajlMpQXq8Xt9vNwMAAQ0NDYY8NDQ0xMDCA2+3OuppMLjeDxZORAcYY4yZkCLK/87+dsTWYcqB1iounlErS2rVrxwSWSENDQ9xyyy1pK0Mg3b7H4wk2YXk8HiorK1mwYEGwJlJXV0d1dTWdnZ3B/pQ1a9YE0/E3NDTQ1dVFW1sbXV1dtLa2Bp87UuQ58XR1dQVfZ+XKlcHyeDweqqurg6tpdnV1YYyhrq4uWI7KykoaGhpobm5O+PXSSkQyasMaFVYb8rMLqPd/34o1kizwWAdQPt5zVldXi1Jq+pWWlgow7lZWVpaW129sbJSOjo6wnwNaWlrE7XYHf25tbZWenp7gz01NTVJbWxv284oVK4LPE/pcocdUVVVFPScWl8sV/N7r9Ya9Zk9Pj7hcLunp6ZGenh5pamoa8/uF/g6B45MFtMsk7ucZVYPxd+q3Aq3GGDHGCODFqr2ANcLMbYxx+5vQlojOgVEqa/T19aX0uIlyuVwsWbKE5uZmfD4f9fX1wcfcbjdtbW3BGoPP5xuzOFjozxUVFQn1uwSWXE7knObm5rCFylwuF6FrWZWXl9PY2EhdXR333HNPWPmjvV55eTkul2vMAISpklEBRkS6RMRE2Tr9j/tEZI2IePxfdRa/UlmkpKQkpcdNlNvtZtWqVbS0tDBr1qwxfSeLFi0KBp/QG3VARUV4K3139/hdwBM5x+v14vP5aGtrC24tLS1jfoeJcLlc09anlVEBRik1sy1evBiHwxH3GIfDwZVXXpmW129ra8PtdtPa2kpPTw/t7e1hfRQrV66kqakpuAzyROzduxew+kqSNW/ePMAa4hy6hers7GTlypU0NjYm1L/S1dVFZWVl0mWaDA0wKqfMxPkX2WT58uUJBZgbb7wxLa/f2toavCmXl5ePuXm7XC7Ky8sTqplEnpeKYcput5vu7u6w5woNWD6fj/b2dmpra2lqaqKuLnJeOmFBJzAYIFpT2lTQAKNyxkycf5FtKisr8Xg8OJ3OMYHG4XDgdDqDI7rS9fqBpiePx8O8efPGNIU1NDSwaNGisH2dnZ20tLQEz+vs7KSpqYnOzk48Hk8wMDQ3Nwf7acY7J5aWlhZWr14dHHkWqEmtWbOGk046KfhhqKKigs7OTurq6ujsDO8tCLzm6tWraW2dxoG2kxkhkC2bjiJT27ZtE6fTOe7opeLiYrnhhhtk27Zt013kGW3btm2ybNkyKSsrE5vNJmVlZbJs2bKMuO4tLS3TXYSkNTY2jjtKbSKYSaPIlEqXROZfAPT392uNZgpUVlaybt069u3bx8jICPv27WPdunXT1lfQ0NAQnM8y0b4XFZsGGJUTNmzYkFCAgeyeUa6SE5is2NnZGXX0WDZoa2tj48aNeDyeqJM9p4OxakEzW01NjYSOJVe5x2azMdH3usPhoL6+nnXr1qWpVEplNmNMh4jUJHu+1mBUTkhmXsXQ0BB33XVXGkqjVG7QAKNmvL5Dw5z7iYvJs098Adf9+/frMGalkqQBRs0okfNcSkpKeNcZZ7Dp/l8xMjyc1HNqp79SyZn4RzqlMtSmTZtwu90MDQ0FO/T7+/vp3/Zi2HHHAOcDRwDPA08Dr8d53sDzud1utm7dOm0jnZTKNlqDUTNCvHVGAmqAu4BXgGbgu8BvgdeAPcBDwIfivEa608grNdNogFEzQrx5Lu8F/gpsARYD+VGOmQ0sBDYDqwAT5Rjt9FdqYjTAqBkh1jyXjwB/BP4tweexAd8B7gVKozyerjTySs1E2gejspbX62Xt2rVs2LCB3t7eMY8vAH4NFEY595/An4HTgbOLnNgPDIQ9fiHwDv/XF0L2pyuNvJp6gVT4lZWVwfxh9fX1rFmzhhUrVtDV1UVjYyPNzc243e7gapder5fKykpWrFgR87m7urpoampizZo1YecGMi43NjZGPba+vj7Yx+f1elmwYEEwPX8g75jb7Wb27Nls2bKFtrY2Vq1aBRDMOTatucciTSbPTLZsmots5nnooYfE6XSKw+GImlNsPsgAiERsvwNZCGL8xzlLSqX1mZ3y1Q/Wyu4ox+8BmeM/1m53yLJly6b7V1cpsGLFirCVHwPcbnfYCpIiIkDYKpiB4+rr68d9nWjntrS0hK1yGe/Yqqqq4KqVLS0t0traGvY85eXlYccnUqaJQHORqVwzXof+ucADQFHE/mXAR4FNWBEjz26n9nw32Gy4bvoW788vIHIFu9nAr7Cq+naHnaWf/Vyqfx01xQKZhiMX8gKCtYHxBNa9T4bb7aampoY1a9aMe2xtbW1YbSdyeYHIxcyqq6uTKlO6aIBRWSdeh/6pwG8YG1w+C9wasc9ud3DxpxsAeNsJc7jqhz/h3IJCfm7Cu/jfA9ycZ+crt6zHdsQxKfgNZhBjMmObgJUrV9LQ0BD1saqqqinJRVZXVzdmNc1ofD5fMPlmIuWqqUk6q0taaIBRWSdWh74N+ClQHLH/80BoNrE8u52CwiL+5/vredsJc4L7zz5nPj/89SPce8lVPJgX3j25fGSYOVsdfOOreQwMJjdhU2WGzs7OuBmTQ2sMsbS0tExqEa9AIIi1IqXP56O5uZnu7m5uv/12gISyPGdaJmjt5FdZJ9ZIrs8D74vYt8qRz2/muCh6/VUOHhigqLiEj124iOVf+AKnnXIyvQeH6Ts0TP+hYfYfHOJtJ8zhs1+5GcdnV3Dw4loK39gZfK7aW5exnKf59SdKufx87ezPRoEbemTTUqhAh3+oQCr/7u5uvF4v5eXlCQWi8V6jq6srrGYS+jotLS00NTVFLU+20ACjsk5JScmYUWMnA9+OOO7BPDvzn3yF+cZgs8HRpYUcX+HkiKLDKymWOw/Pitl/cIhnXtvHgcERhsoreHbNrVRffRFmdBSAo9nNBhbz9R/8ho/XDoU9j8oOgZt55JLIXV1dYSnuXS5XWH9HbW1tSmsHgSWRIwNd6OvU1NRQXV1NT09Pyl53qmkTmco6ixcvDltu1wB3EN7v0gN4zr8Y/MHl3ceVc8axR8QNCmWFDubNqaCixAo6vur30rX0i2HHzOdh3rX5Z/ypYyDaU+SeMePupmmbgKqqqjFLDLtcLurr62ltbaWpqWlMZ3o8dXV1VFdXB7dEBGpS8YJWVVUVPp9vTFmziQYYlXWWL18eFmCWAR+IOOZLDgfnNtyIzQZzjytndklBQs+db7dx1vHlzDnSCcD2+s/Tffb7w45ZU/p1Sgr2sqfv0GR+DTVNGhsbaWpqivqYy+WK23wWTUtLCx0dHcEtERs3bow7jyZUrH6abKABRmWdyspKfvUrD8Y4mYOdmyMe/53NhuuHP+W4OXOYe1w5RyYYXAKMMZx8dCmnvKUU8vJ47uZ1DOU7g48X9+7mxDtvw7tLZ/Vno9raWtxuN3V1dWMei3Uzj2xSm4zOzs6wCZLxXsflcrFlyxaAqKtUprJc6aB9MCorlZV9FJGn+Q4fo5iXgvv7HQ767riX9847m3cdO/HgEuqE2U66BwbZuuMQO08+mUue3xp87C23/i8f/sk6Xjt4kJKSEhYvXszy5cs103KWaGxspK2tjYaGhuBMfp/Px+23305g9dvADHsg+DXRprPIc6urq8Nm8ofWdALHut3u4NycwOu0tLSwcuVKPB5PsP+oq6sLj8fDxo0b8fl8rFy5ksrKykmNakuXjFsy2RhThZVvsElE2kL2r8Ca97YRqAAWiMj4A8nRJZNnotd9B+i8+3HOX/bhsP3Pf30tO91X8I63lXFseeRsmIn7zYO/5ZJFdRQODvLSyAhHhjz2f0Bg2qXD4cDhcODxeFi4cOGkX1epTDCjlkw2xtRiBY9YM4rqsRLeNgCrp6pcKrOICNt39/PelvC3QO+p72TnRZdRUZKfkuDi9Xq57JJFHDxwAN/ICN+MePx6Dr9Rh4aGGBgYwO126+qXSvllVIARkTZ/rSVaw6JPRGb5tzoR8U11+VRmeHP/IfLbn+CoR/4Qtt/7uZXkOfI4/a1lKXmdyIwBtwHbQx53AN+KOEfXjFHqsIwKMIkwxlQZY9Kfy0FlrFf29nPyD8JrL74za9jzwQWcfFQJhY68lLxOZMaAQeC/Io65DAgdaKprxih1WFYFGGOMG+gCqowxyU+jVVmru38Q+yN/pOLxv4Tt937uJmaVFHB8hTPGmRMXLWPAr4AnI/Z9LYHzlMpFWRNgRKRZRDwi4hMRD+D299lEZYypN8a0G2Pad+/ePYUlVekyOAjbdvaNqb3sfd8H2P++c1LWNBYQbe0XASJHlpwPnDLOeUrloqwJMP7RZaE6sdaUisofkGpEpOaoo45Kb+HUlPjZXSN8v+rPHLE1fGaz93M3UXlUCUX5qWkaC4jMGBDQirX8cqgb/V8dDgdXXnllSsuhVLbKigDjDy6bI3aXAzpcJ4f84Psv8Z994WnWHz/2eLYfdSTHzZr8qLFIkRkDQq2N+PkqrDH0DoeDG2+8McoZSuWerAgwItIJLInY7QLumYbiqGmw7kcPMvvZMzmbXWH7r39jJ9dc8CF+//vfpfw1Kysr8Xg8OJ3OMYHmXmBHyM9FwGftdn664Zc62VIpv4QDjDHmTGPMmTEeuzjWYxPhHyG2AqgBVvq/D+gyxqzw9600ATpUOUd4vV6Wf+ESvkh47q8HgKdGRjiQxvknCxcuZOvWrdTX11NWVobNZqOkpITKd5zOrY78sGNXlJTyzrnvSXkZlMpa462pDHwRGInYbgVKI447CxiZzPrN6dqqq6snuBK1yiTXX3+DnJ6XNyaH7jlWn7sA4nA4ZNmyZVNWppGRUXms0ytDxSVhZXr+W9+TQ0MjU1YOpdIJaJdJ3Hvj1mCMMbdhTVi+CTjPv63CWn7DZ4wJDucRkSexMqcrlVIbfrGBz4+MhO17Ang05Oepnn9isxlOdL2V1y++Imz/8T9r4g3fgSkrh1KZLGaAMcacBSAiJ4vId0Vks39bIyLnYaV06TLG3GOMuc4Yc8RUFVrlluLeXj4dse9/oxw31fNPjikrZNfV9Yjt8L9RybYX2f+b305pOZTKVPFqMPNF5PpYD4rIPhG5XUQWAS34+01SXUCV23wDg3ze4aAwZN924L4ox071/BNjDEfPPZU3z/tE2H7n2tvo6R+c0rIolYniBZjtcR4L4w82m0Xkuykok1JBr72+l2Um/G16C1ZHYKjpmn9ybHkRtx8RPizZ9fwf+dqVi4ODAsrKyli6dKkmwVQ5J16Ayaw8/irn7Osb4eCPfk7Z4OHRYz3AT6IcO13zT+x5Nk6+4mwe5+yw/cf+2kNvby8iQm9vL+vXr2fu3Lls2rRpysuo1HTJinkwKjfdefchjvjB+rB9TcbQH/Kzw+HA6XTi8Ximbf7Jpz5SgOfI68L2XSVC6MwZTeevclG8AHO2MSah5E7GmHP9c2E2pqhcKseNjgreH27m7WwL7hsydl668JKwpqf6+nq2bt06rYt8FeXn4bj2QvabwyHlLcAnoxyr6fxVLom5oqUx5iSsJTDcItIb5fFzgUXALGC1iDxljBkRkdQmhEoBXdEy+7Q/e5A33uXmExwekfXPD9RReN/POGF26jImp4pvYIiNRxTTMHw4vX8b0ZPllZWVsW/fvikrm1LJmuyKlvZYD4jIdmPMvcDLxpg2rPx+s7FStNRiLQrWICIPJ/viSsWy4y8v8CkeCtvXu+waXOWFMc6YXuVOB7cODxGaKa0WqGRswjxN569yRdw+GBFpBi7B+j9ZgzUMuRK4SUTeHgguxpiTjDFfYgIjz5SKpffgENWPN2ELGWfy6nFVlH3w37HnZV63odfrZenSpWwFHot4LDKBHmg6f5U7YtZgAsRawni8KlK5f4iyDlNWk/ba63s5+b5fhu3rXXY1rhQuJpYqmzZtwu12B1e+bAbeF/L4NcBXgEDDmabzV7kkJR8H/WlilJq0oZFRzN1349h/uI9icFYFUrcoZUshp4rX68XtdjMwMBAMMBuB0AysRwMXhvys6fxVLokaYIwxN/s78ZWaMl6vl89c18Dg174Utv/F887nhGMrpqlUsa1duzYYWAIOAJEZ0erJjOHUSk21WDWYJuA8/5LDP05FKn6l4tm0aRNz587l1bt+ytzR0eD+EeDiX/+Kv/yxbfoKF8OGDRvGBBiwmslCzQf+44KLpn04tVJTLWqAEZHtInKTf3haM3C9MWaLMWa1MWbOVBZQzXyhTU3XR2RN/n/AS4cOZeQExVijwZ4F/hax77NHHovL5Up7mZTKJOP2wYjIkyJyvYjMwxrav8YY83t/BuWEJmIqFU+gqekY4OKIx9b5v2biBMV4o8Ei09kU33kb5ZqbTOWYCXXy+xNaLhKRj2ClhfL4g81F6SmemskCw3t//OMfMzQ0xBIIS6/yPPBH//dTvd5LIhYvXjxmKeWAeyAspc3RBwd4b1+f5iZTOSXpUWQicq9/XZhFwGxjzB+MMRt1cIBKRKDPZf16K9eYHcImKYK1bGqo/fv3Z9Qn/+XLl8cMML1Ya1iE+kzI95qbTOWCSQ9TDlkX5jyslS+rdXCAiifa8N7zgWNDjulj7GgsIKM++VdWVuLxeHA6nWMCjcPh4Oe28H+vT2HlVQqViU1/SqVKSqdF+wcHfDdkcMClmgBTRYo2vHdpxDEbgP1Rzs20T/4LFy5k69at1NfXj0nC+URhYUiqTigALo84PxOb/pRKlZjJLmcSTXaZWcrKyujtPZw/9VTghYhj3g1sjfMcDoeD+vp61q1bF+eo6WWz2VglwrdD9nUC1VGOGxmJXEJNqek32WSXmZfYSc14kcN7I9fl/gvxgwtkxyf/kpIS7gRGQ/ZVYQXPyOOUmonSFmCMMV9M13Or7BZ6Q3UCV0c8Htm5H0umZyVevHgxuxwO/hCx/5qQ7zU3mZrJUhJgjDFLjDHdxpi9/q0baEzFc6uZJ3R47+VAechju4B7E3yeTP/kHxhlFjknZjGQ7/9ec5OpmSxlNRgRqRCR2f6tgrEtH0oB4cN7Izv31wODCTxHNnzyD4wye8hexN6Q/bOBi2x5mptMzXipCjBdUfYlNXrMGFNljGkxxtRG7C83xqwwxrj9X6uSKqmadpWVlbS0tPCB/HzOCtk/CtyRZ6ewsJCCgoK4z5Etn/wXLlxI26NbuZt3he1ffvQxbOl4UnOTqRktVQHGa4y5yBhzZmAjiSYyf1CpwFo1M1IL4BERj4isARqNMeVRjlNZoOacc/n6W+aH7XsIOx/6zLU8++yz3H///THnl2TbJ//3vvdkfJetD9tXvetfOA/N/BGcKrelKsDcBHwZa9XLwLZook8iIm3+Bc66Q/f7A4lLREJrSl1Yq9KqLPTK09t5/6utYft2X3MP65t+TGVlZdz5JdmYlfhzd5xF72lnBH82o6OYn/98GkukVPqlKsC0ikiNiJwX2LCWwUiVGsLXccL/84IUvoaaIv2Hhim/6w4cDAf3vWQ7hfd9tRZjTHBfZWUl69atY9++fYyMjLBv3z7WrVuXNTWXUEcUOdiz6IqwfUe13E1336FpKpFS6ZeqANMTZV8qp1mXE1GrAfZiNaepLLPjDR8n3Bc+h+WF86/FdWzxNJVoajgWX8GoIz/4s/PVl/H9fvM0lkip9EpVgKn05x67zr8tIfXDlCcUTIwx9f6caO27d+9OcVFUsgaHRxn1eCjYsyu4b9hZTOXNl5Fvn9nzfo8+6Vh2z/9o2L6iDXcyODwa4wylsluq/qMbgH1YufxmYdU4ZqfoucFqDovs0J/N2FpNkIg0+5vtao466qgUFkUly+v18pklDfSuWBa2f9v8hRx7wjHTVKqpk2+30X95+NDqit8+yBuv7YpxhlLZLVUBZqV/BczvBjZgSYqeG6CdsTWYcqA1yrEqAwXS8//zrp/y3tHwT+yX/OEB/vanzFsSOR1KPvExfGVHB38uGBpg9alzMmoZAqVSJSUBRkSiNSRH65dJ9vl9QLsxJnT4cg3WCpsqw4Wm518akdTxD8DWDF0SOR3aH3uEH+0Pr3gvHjyUUcsQKJUqSQWYyEXFQvpeQvtgmpJ43ipjzAqs4LHS/31AHeAOTLQElvgDj8pwgfT8RwKXRjz2f/6vubAuitfrpa6ujp+GjJ4DeD/gyrBlCJRKhWRrMGsiFhO7nsP9L0n3wYhIp4isEZFZIrLAP6Ey8JjP/5jH/7UzybKrKbZhwwaGhoZoAApD9ncBD/m/z4bsyJMVCLRe4JGIx671f82FQKtyR8z1YIwxe4FzReTpcZ/EmLNE5Mnx9k0XXQ9metlsNvJFeBkI7cr/IrA24riZvC5K6Do4VwKh0yx3A8dh5WErKytj3759U19ApSKkcz2YWUCnMebC8Z4kWiDJlOCipl9JSQlXEB5c9gO3RzluJgtdXqCF8E7Ko7CWVI48TqlsFi/ANAOrgHuNMcsjHzTGHGGM+ZIxZrU/D1lZ2kqpstrHPllH5BvodsKXRM6G7MiTFRpADwKRDYL1UY5TKpvFCzDi7wNZBHzXGHNrxIP7/EOSVwGVQI8x5ndpLKvKQn2Hhql/+zs5PWTfMPCDiOOyJTvyZISugwNja3DzgdPy7DM+0KrcMW4nv4h4sEZ1XWqM+V20mop/3ssNaG4wFWH77n7OfDB8CbGNwKv+77MxO3KyQtfBAXgW+FvEMdcZWPrZz01puZRKl3gBJjix0T9iqwZ4O9Z8lBMjDxaRZqzZ/CoHeb1eli5dGpb5eEnD9ey4/34qtoTfRm8tcmZ9duRkBBYgC12GIHIs/1WjxbyxWzNPqBlCRKJuwJYo+47Amj2/F3h3lMf/EOv5ppuXfx8AACAASURBVHOrrq4WlT4PPfSQOJ1OcTgcAgQ3u90hG215IhDc9p79fnnxjf3TXeRptW3bNlm2bJmUlZWJE6Qn5PoIyJp/3yAjI6PTXUylBGiXSdx749VgqiJrKmL1uywAPEQfYaYTH3NM6Cz9oaGhsMeOHR7iotHwYcc7rrmB42c5p7KIGSd0GYK9g8M8/77wrEpVf/sZz28/OE2lUyp14gUYA3hi9Lk0YI0w80SMMIuZfFLNTIHJg9F8AbCH/LyjrJydp7+Dovy8KSlbNih05HHEV64K2zd/tI2/3v3iNJVIqdSJF2AqgXuA9caYL0YGGrFGmF1ClBFmKncEZulHeitjV5z7Rm8vdQvP0XxbEd5yTg1db50Xtu8j22+np39wmkqkVGrEDDAisl2sYciLsEZUjlmPRcJHmP0ea3KmyiGxJgWuJDwtzKvAXTKi+baimF2cT3/Dp8P2HXvv3bz2uq5jpLJbQrnI/H0vL8d4LDDC7GTAnbqiqWwQbVLgMYytvXwHKw0KaL6tSMYYSq5exOARhz+fOfb7+PMVdWGj8jSlv8o2qUrX3wVUA/eOd6yaWSInD4JVeykK+flV4CchP+dCYsuJeuvbZrOzbnHYvg+1P0Zvby8iQm9vr6b0V1knZWvUipXteFGqnk9lh8jJg8dgLW8aKrT2EqD5tsK9vuNlGrvfDEvk/06gNuTnIU3pr7LMzF4EXaVd6ORBu90+bu0lQPNtHRZY7XP9A/ePaQL4fJTjtYlRZYtJBRhjzOpUFURlr4ULF/LbR/7OFZ+4OKHaSy4ktkxU6Dyi4aGhMTnaPoHVuRlKmxhVtphsDaZ2/EPUTHdoeAQpewvfKClLqPaSC4ktExU5j+gxYEvEMZ+Ncp42MapsMNkAY1JSCpV1QnOPFeU7+NyZJ3HMhjvCjrnZZgurveRSYstERZtHFFmLuQaInO2sTYwqG0w2wERfDlPNaME+g/Xrg6OcvnlwgHxGg8fsP/IY3qxbHDbMNpcSWyYqWk3kHuBfIT+XYgWZAG1iVNnCPv4hSh0W2mcQ8CEOr8YYcMc7/pPG/1uG56g7p7J4WaekpCS4jHLAEPBj4Bsh+24EbvU/pk2MKlvoKDI1IZF9BjYgcjzTYxgePuYFTqjI7aSWiYg2jwjgNqxVLwNOBAJ1lgMHDnDWWWfpxEuV8YyVkTnJk43ZIiLzxj9yetXU1Eh7e/t0F2NGKCsrC/vEfS2wPuKY9wDPl5bSu38/Kj6v18vcuXPDaoQB/wf8R+ixwKlAID+1w+HA4XDg8Xi02VGlhTGmQ0Rqkj1fazBqQkL7DEqAb0U8/gvgCWCgv38KS5W9oi1CFrA2Ly9skEQlcHnIzzrxUmU6DTBqQkJHL63CmrkfcMC/L/I4Fd/ChQvZunUr9fX1YYMizrnqWv5w/JywY/+Lsf+0OvFSZSoNMGpCFl12OXl2O2cAX4x47LtYc190lNPEhS5CNjIywr59+7ij+ces3LMrLH3MqUBdxLk68VJlKg0wKmGDw6Oct+haCvLs/BTID3lsJ7DG/72OckoNR56Nfxw8wIaI/f/N2AloOvFSZaKsm2hpjFlhjGk0xlQZY2qNMY1TXYZcJCI8u3MfFW89gfVnXERkr9/ngEN2u06kTLGSkhK+w+GOfYAzgMi1yrVJUmWiyQaYyNr6VKkHNmMl7tV8aFPAu7uP7r5Bhre8xEUdnrDH7gE2lZTy6Wuu1YmUKbZ48WK22x1sjNj/Pxz+59UmSZWpJhVgRGR7qgoyAT4RmeXf6kTENw1lyCm79h/k5T0DMDTMCZ+7kYKQsU27OZLiB5/l0ede4Y7m27TmkmLLly8nP9/BtyP2vxu42v+9NkmqTJW1fTD+JjLXdJdjpus/NMxz/9rPzh0vs/3yT3H6/vD5RHd/+EsUuo7klLeUTlMJZ7bAMOauIie/NOEt0t8Gjiwo5Ifr79LArjJSVgYYY4wb6AKqYvXBGGPqjTHtxpj23bt1bfNkHBwa4ckdPh57pI3bPvlBrnw+PLjch2Hl377Kts6/UFKgWYfSZeHChTy2pYPN57s5ELL/GODhTy5izpnvZ9/AUKzTlZo2k5rJnwmMMV6gQUTaYh2jM/kn7tDwCB0v97Btm5cvf+pD/O3QQU4IeXwv1oqLbwJOp5OtW7fqp+g069zRQ/l3vomr6fCcl1FHPn978FFsLhfvcc0mz6YJzlXq5NxMfmNMVcSuTmDBdJRlphocHqXzFR8DgyPc99Nb+fngobDgArAUK7iATvSbKq4ji3nl2v/g4NGHp7fahgZ5+9pvMjA4wh8e6+emmyDLPzOqGSSrAow/uGyO2F2OlaZJpcDQyChP7uih/5A1ve+D9/+ScyPuWLdgjRwLnqMT/aZEuTOfsqNn4f3PL4ftf8sfHuSVH7Vz8XlOGhvhB5ELyig1TbIqwIhIJ7AkYreL8PudStLg8ChPveqj96AVXAp++SA3RiyG9SdgRZRzdaLf1DjpyGL+db6bfWecGba/6tb/ZnDAWo/nS18S/va36SidUuGyKsD4dfknW9YbY5oAHaqcAv2Hhtnycnews7j3wYc461vXhx3zGrAIwlKXBOhEv6lRUZxPRVkh/7zpm2H7z+IpbuJmAIaHDS0t2k6mpl/WDf3x12I6p7scM0lP/yBPv+ZjeMS6KXXd9RPqbv4yoSHjEHAxsCvK+TrRb2qd8pZS7j/yKAZOOIkP7jg8Fe2rfJXWvDN475fPpf4GGLvQslJTKxtrMCqFdvoO8OSrPcHg0vfXR/jUzV9mVsRxn8NKwx+NTvSbWo8+3ErDhedy6es7ggMtAByM8rPRizj2qPt4vecAj7Y/y9KlS8MyNOsiZWpKiciM36qrqyWXbdu2TW644QYpLS0VY4wUFxfL6aefLkXOYjHGiLO4RM6/9Gq596f3ye4ip4g1ECm4fQOEKJvD4RCn0ykPPfTQdP+KOWPbtm3idDqDf4MLIv5WAvL9vDz5wje+JwWFRWJ3OPRvppIGtMsk7r3TfvOfii2XA8xDDz0kTqdTHBE3msjt5Lw86TJmzM3qe3HOWbZsmWzbtm26f8WccsMNN4z5W94RJcjMt9ni/r2dTqf+7dS4JhtgtIlsBvB6vVGbQh5++GHcbjcDAwMMDcWe6f0e4K8jI5wk4R3DTcAXYpxjs9lYt26dTq6cYhs2bBjzt/xP4OWI4+4YHWV2nOfRuUtqKmT9TP5EzOSZ/Js2bcLtdjM0NBR243E4HIyOWsNWR0ZGYp2OG/g5UBSxfwNwFTAa47yysjL27ds3iZKrZNhsNqL9z34A+CPhnaqPArUQtuxyKP0bqvHk3Ex+dZjX641ZQxkaGmJkZCRucFkBtDA2uPwCK1NvrOCio8amT6zh4H8Gvhex7xzgjjjPpXOXVLppgMlia9eujdv0FUsFcC8QLUvoN4HFhC9wFUlHjU2fxYsX43A4oj72ZeCRyOOBr8Z4rtHR0UmPKovVPKsj1RSgnfzZrLS0NG5HbrStFuS1KJ3CgyBXjXOujkCafpGjyCK3WSAvRPn7XpmGv2msAST6Ppk50E7+7JLKT3wTaeIoANYCrcCxEY/5gI8Ad8Y5v6ysjPr6el2xcpoF1odxOp1RazI9wMeAyAUq1gOfjPJ8Q0NDDAwM4Ha7J/QeHK95NpnnVDPQZKJTtmyZUoNJ9Se+RGswF4B4o3yqFZAnQN4+zvk2my1NV0Qla9u2bbJs2bKYf7N/AzkQ8bceBrk2Tk1m2bJlCb9+tOHSk33OTBY5l6y0tFRuuOGGGT/UG50Hkx0BZrymDSYwN+HA4LBs390nF15xjeTZ7TGf7xSQTTECywjIt0DsCQSosrKyKbhCKhnGmJh/t0Ux/varYgWE/DLp6x9N6HUT/XAzE947udwUONkAo01kUySRDvmhoSG+9rWvRW1C++dL23hz/0Ge3NHDX7ftYduuPj61uAG7fWwzyduAHwLPAB+N8jo7jOHS407gq3n2qIkrQ+mIscwWL8noPcC1jB2w8R3g+4wd4TM0uJ+3vvUGfnL3c4yOStzXTbR5NttHqmlT4CRNJjply5YJNZiJdMhHflKy2x1SWFQk377tF9L63Bth27dv+4UUFBZJnt0ubwP5QZSmkdAmknV5eXJ3S6u0PveG3Lnp71JQWJSSWpWaHok0VV0AMhDl/fBHkBPGHH+4RlRcUipLGq6P+vfPlRpMtjcFTrZpD20iy44AE68pI9GtoLBI7tz09zFB5te3/VI2nXxqzMASuJmcAZJnt8sFl10TPHfN7b/M2er/TJBI06stL08+aMuTnijvCx/I4njvu8D7NuLmlCs33mwOpKlo2tMAkyUBJpkhxZFbWHB4Zqc8+aOfy55/+2DMoCIgL2O1xYfVSkpKpfW5N6T95W45MDgc7DAuKysTm80mZWVlmmcsi8S7kRQWFckXvmklvnwXyOsx3icbQY5J4D1otzukqMgpt97WnLI+xak2kRtvoh8MM20gTCr6fPsODmmASWTLhACTyCe+RLZTncWybdmXZODY4+MGlldA6kEcUZ7DGCMv7+mb7kuiUijWh4SnnvuHPPzCm8Gm1OPy8uR3Md4z/SCrQcoTeB8WFBbJl759ixQWOcVuT/wT8nSPxprojTdbazCTqWGOjIzKS2/2yl9f2q0BJpFtOgPMgcFheXP/AXn48aelsCj+GzvWVgxyGcjvsUZ/xQssL8cJLIGtNMP+GVR67e07JA//4025c9Pf5YLLrhFAlhG9X0ZAukFuAimL8x7Ky7Nq04HndJZYAaPIWSyVp5wmzmJrKYiSUqsf56WXXkrZaKzJBKmJ3niztSkw2cDYd3BIHu/aK63PvSF/0QCTfICJ9ia9/PLL5YorrkjqjTs4PCI9/Yfk1e5+eeFf+6X95W555MVdMTvkw/7QUarhRSAXg9wT50YQuv3V2GTVO94lBXl5WffPoNLPNzAYfD86i0sEkNNA2uO8p/pA1oOcHeO9FGhqHe/9nWe3S35+gTjyC8atOTz57D+k7+CQHBgclsHhERkeCR82PdkgNdEbbyqnF0yliTbtHRwalq7dffLwP94M/j1TEWByIptydXWNdHQczqYcKwNxNA6HA4fDgcfjYeHChYyOCn2Dw/QdHKbvkLX1Hxrm0FCs1JDhdu54mXt/3kTbAx4G+noBMMYgIhwFfBy4ADgPKB7nuQ5iDUX9AdYa0vn5BWAMg4cOxjzH6XSydetWTbOfg/oPDfPkDh9r/udLPOTZwMjwMHZgCfAV4K1xzt0K/BJ4EHjWv88Ywx+e/Rdgva/rL/wwhw4emFQZ7XYHC91X4L7qet52wpzgfpsN/vXqK1z3yQ9x8EDs1ygqcvLQn/5OZWUl+XYbjjwbBXYbhY488mwmZjbqSDabLZgoNl7G8tB7QyYpKyujt7d33ONKS8v4y/OvsKfvEKMRt7Ci/Dz+/e1HaTbl8fQdGua3f3mSa+uvp6SkhI997GPjrpESEBjrftHFbu79YzuP/HMXT3R18/zO/ezYO0B332DCwQXgbSfM4bP/vZoft7TiLCjkfcBXRPgr8AbwU+BC4geXJ4AbsG4IV2EFF4CR0RHOeu85FBQWkWe3h53jcDhwOp14PB4NLjmquMBOzZxZXLlkWXD+1DDwY+BkYBVWqplo5gKrseZWdWHNs3IXFODYayWl8dx5G8PDE0+8Gml4eIiHPBuov/DDPPHo5uD+0VHY+JMfj/s/Ozg0yC233MKLb/TyzGv76Hylh8e8e/njC7t45MVdFDnH+9hmCZ1ftHDhQrZu3Up9fX3Y/LRMTp0ULylqgN1u58OfuJhd+8cGl1TJiRrMcXNcsufNNxgeHmJkeLyphdHl2e18vO5KPvvfq5MvyMgIpf98nvL2v9NzVzPveP1Vjkjw1Jexait3cfgTZDTOklJ++uvNbPrVHdx3zy/p6+ujpKSEK6+8khtvvFGDi2JoZJT1d9/HjfWfHvM/UQ58DrgOOD7B5+s/0cX9r+/g78PDPIP1/uxOQTkLCotovv+PgBXAHvjVzxI+9/xLrx5TCwL44TdvCtbeYrE7HFz9mWtpuvVWbDaTTNGnTaCF5fkXXuLcf5vHgQMDMY8NXN/IaxSQihpMTgQYfzvjpJ+n0FnMggvq2PyAhwMD/RQ5i5l/vjvqGxkgr3c/Rzz7FGVbOyl/8gnKn9yCvW/8amvAa1jrtWwEHk/wHGMMIyMjGJNd/xhq6v39qee5+bv/S+tvxjbX5mFlgVgCfALIm+Bz/wv4B1ZtZ7t/e8W//w0gkYa0PLudmvd/mKce/8uEPxzm2e3Y7Q7+5/vrOfuc+cH9iTTlBW68x82ZgzPfTlmhg9LCw18zJegMDo+y9fkX+P73b+HXLRsZ6O8Luye99oqXb/zndWOuXaxrE0kDTIL8HV4pkWe3R/1jfWPNrZwxPMKrG27HubWDs4aHOZ2Jt0F2AL8BHgCeTKJ8ukqhmoiDQyM8t3M/Pf2DMfsHjxHhAuB8YD5QmILX3Qe8idUkF7r1hWz9WIHoAFZ/44EoW7//2Fi9jtE+pT/x6Oakb7zGWE2NJQV2SgutryWFdgrsEw3BifF6vaxdu5YNGzbQ19dHcXEJH7uwjouuuh7vtpfG/T2OO7Ey+Dc90N9HUXEJtee7ufjTDTFrLjt3vIznztvY/ICHgf4+RCTpiKoBJgklwJlAlX87CzgdsMc7KYbdWCn0A9vrkyiXw+Ggvr6edevWTeJZVC56Y99Bunb3MTB4OHNZaMAJ3Jw+vvCT3PCOuVR6X6T8qS2UvPActjirpk6VYaxAsw+rea4b2AvsAd6w2Sg54yz2FBbS9syTvHBgACly8pZjj2fXztc4eGAgoRtvPPY8gzPfjjM/j+ICOwX2w4MLAoMN4hER/vHiS3zve99j46/upr+vj/yCQob9fU4jIyG1N2OscX7jGK8JLJpowVcDzDgmE2BKgOqI7e0kPzpiD9byto9irT74NNZ4wVTQEWJqMkSE130H2L6nP+GBK7aBfo549inMnzfz4p1NnD46wjuA/PQWddJeBbzG4LXZOPHCS6n4yAX0vf00Bo882rqBp0GezQQ3A4yIMCowKsLfH2nj61FqI5N6vSj9xjt3vMy9P/sxnQ94KBro54iiIt7/gVrO++gnsQ0P883/+jy+wUP0Y9UO95BjAcYYUw7UYzXvuoA2Eekc55yEf8lTgX8H3uPf3snE258DRoHnsfpPHgf+ArxA6gJKQCYPl1TZZ2RU2Ok7wOu+A/QdTPxmF/j0y9AgJ42MUAmcBFQag8sYao6fwxH9fRTs3Y09Q+87e7CGZD9jdzA0732c/JlllLz3HGuc9DhCm5YS6aMNPS8VQ7wDCrHuY28HTs3P5/pPXUrRazsY3f4S8q/XmU3i9zRD7gWYVqBBRLpCfq4TEV/sc2wS67Z+OnAu8AH/9pZJlM0LtGP1o3T4v9+fwHnR+nUCVeCRcZof8vPzWbJkiY4QU2mx78AQr/cc4M3eg4yMjH+viNasNqbpaXSUR+9+nvvWfILSkYPMAmZhjWArxmo1CGyFQFGczQmUYq3Ymi4D9kJ2HlvNzuPPZODd76b042dw4Pg5YUHniUc38/XPX8fI8FBYc1Yi/TqJjGyLxoEVSN7l396JdT9zkbr5JzkVYPy1lw4RqQzZ1wS0iogn9nl5YtUnrDfxAqwlgj8CHJdEOUaBF7HmnwS2p7CWHk7GBZddE/yHLC4p5fLLr+CSRW7OP/98BgZiDzPUJjE1VUZHhY5n/8H3vvc9fuPZOKFP6LHE6my35R3+fD2aYP+OAysglQMVwGz/16Oxlgg/zv/1BP822S75YWcxfae+k953nMGrRx/Dqlv/l62Dg3HvAYVOJwsuWMR7PrCMb36hGptNyLPDQF85o6OxR5fagTlYASSwnQGcRvqbInMtwNQCjSJSHbKvESgXkYZY5xXaTpLr5HUuZJgPIhPujP8HsIXDNZOnsNonU6GktJRtr+2iojh/TEdgNs4gVjNTrPei3e7A7rDzlVviD3mNJbLGE7NjO4XysW7YbwdOwbpZB2oAzkk+95tYHz5fxhqw8zqwE2vAQT9wKC+Pg3l2BgZvp4D3k88g+byDUqxgGNiOBSqxaiMnktwAonh8/jIdBAaBIayBEgVY16DYvx1JbgUYN1bz2IKQfSuAeSJSF3FsPVZfDdVQ3U5i+m0lPG4v55HBnfydUZ7AGpmSDomM+vJ6vdxyyy3cddddOmlSTQuv18vcuXPHrU0/8tgWjj72RCuF0sFhBoZGeHX79oT7JSbSFxFoerrkuv9g4/p1Kekct2Hd0M/i8AjRaqzaULbpAv5pszHn4xdhO+0MVn73a7yMFQD3YAWURE0mwEx7IsqJbIAbq4ksdN8KoCXeedVxEvoNFxXJ7g/Mlxe/9FX5+8bfSdvTrwWTAaZ7y8QkeUpFSjajcKzElNYKrU753h2/kj//c1cwweL5l149NhFslM3ucAQzOQdWZr3gsmskzz755TCibSdiJZ5dV/Zu+ccJZ8tgWfm4yWenansV5LcgN2MtHHcWiDP0HlNcIudferUUjpOwM94mk7hnZ1sNphZokvA+mHGbyGqMkdAaTF/lKez5QC17//3D+KrORvLDuwjPO+OtCSXES5Y2calskmjixNBJvonWegJ9iKOjQnn5EQm9TklpGU9ue51REUZGrW14VHhlu5dLP3pO3GSYkxHssL/lds45+TRKX3iO0heeYettt3DyyAinkL7BBv/CGoH6XMj2LDFyx0XMk8mz24P3s0T7tELJJGowqW7aS7d2rL67UOVYcxTj8p1Zw+75C9k1fyEHTnTFPbbIWcxAf1/ypfTLz89n0aJFiAgPPPCANnGprNTXl9j/Quhxa9euHTcx5dDQELfccgvr1q3DZjMJv85Afx8nH10yZv+Zx1dz3733JpwpfaI33pHhYUaGh/nGjUv45o/u4tG//tGa7e4/18bhvp3jgLdh9aW8DTiCw/0axVid54Mh2wFgV8S2Hdhus9EFDOUXBJsOA+l8Qn+PPFsegjA0OBgWXALlni5ZVYOBqMOUO4D5EmeY8mmnvVPW3fdwwq+R7LDBUDrCS80UydRgpuqcaKL1W15wwQVRP+h94oJPcdGFn+RAnJpWJJvNhjE2MOm/eV9w2TVhQ7xjDQPv27+PP/3+N3HLk2S5R0Uk6UF32RhgJjzR8tQzzpQf3fP7hF9j185XuPaCD3MwXibSggJsNhvDw8M6wkvNaEuXLmX9+vVxawSRA1aSWXclmddJhYmsDzWVQtfbGc8nzz45oVYXZ3EJF15yGf+v5Vdx09Hk2e04HA4OHjjwkoicktxvkIXrwYiIT0TWiIjH/zVucElUXp7hmCMKOeuEci6bP4/77vXgdDrHrKkQWFfl/vvv55lnnsmqNSKUSsby5cvHXVvE4XBw4403Bn8OXU8lntDjknmdVAhd7yWTFBWPfw1tNnAW5HFgILGJEwcPDLDhjmZ69+9ndHSUAwP9PNb5NHWLr8ZZUooxBmdJKR+vu5I7H/gTJDZXPLbJjBDIlo2Q0RSBkSeB7U8v7pLXewbGLM0qYi2XumzZMikrKxObzSZlZWWybNkyHfmlcs5ElypO9cizRJdEnqxEl1QebysoLIq+PHqCW57dLhdcdo20PveGPPzCm/LXbbul85VueeFf++WVPf2ya/9B6T80JKOjoxMqd2Ap6GgGh0fkyR09KV0yedpv/lOxhf7RCgqL5Nu3/UJan3tDnnnNJ4eGRsZ7zymlZGIfuCazlv10frBLJDAmshmbTTqfeV6uua5B7Ek8X5HTKU8/94IMDSd2f0o2oEfj3dUrbc9rgJlwgAlshUVF8sTTz417oZVSyZvu2kgyEgmME60pTOQ5k702kwno0ezuPSiPd+2ddIDJuj6YVBgZHubO5lunuxhKzWjZuJZ9ZWUlHk/s/te8vDzy8uIPqnI4HFx55ZUJPWdg5VljzKSuzXjldjqdeDyehEe1HllSwJnHl0+oDFFNJjply8YE2yKVUrktVjPd5s2bM7rpL9WvwSRrMFk3TDkZ0daDCR0eqZRSicqlJLTGmA4RqUn2/JxsIoPEh1EqpVSobGz6my45WYPRteuVUmp8WoNJQjomaymllAqXbckuJyW0jVRzhCmlVHrlTA1G20iVUmpq5UQNprq6mvb2RNe0VEoplQo5U4NRSik1tTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFIqLTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFIqLTTAKKWUSgsNMEoppdJCA4xSSqm00ACjlFIqLTTAKKWUSousStdvjFkBzAY2AhXAAhFZOb2lUkopFU1WBRi/ev/WBiyZ5rIopZSKIdsCjE9EZk13IZRSSo0vK/tgjDFVxhjXdJdDKaVUbNlWg8EY48ZqHqs1xjTE6oMxxgSa0gAOGWOenaoy5oAjgT3TXYgZQq9laun1TK1TJ3OyEZFUFWTKGWO8QIOItI1zXLuI1ExRsWY8vZ6po9cytfR6ptZkr+e01mD8tYzqcQ5rFJEu//FVItIZ8lgnsACrRqOUUiqDTGuAEZHmRI81xlQBm4HQTv5ywJvqcimllJq8rOnk99dcIoclu4B7Ejg94UCmEqLXM3X0WqaWXs/UmtT1zKo+GH8tphbwYTWtNUU0mU32+ZtEpCFVz5eLjDHlWH+jCqzmy5WBJk6VGP9Ali6gVkTWTHd5spW+F9Mn0XtlVo0i8weTlAWUUMaYWkA7BydvEVAuImuMMQArAQ3aCfK/DytExGOMwRizQoNM0vS9mAYTuVdmTRNZgH8OTIv/lwzdX26MWWGMcfu/Vk3gOcuxPjF2p7q8mS7V11NEmkNuiJXkeB9ZEtd3AdZ7Eaya+oKpLG8mm+i11PdifMn870/0XplVNZiQCxFtkmUL1pDlwIizVmNMnYj4EnjqGhFp83/KyRlpvJ4BrlzOFZfM9cUauBLQjdW8NwJk9gAABFFJREFUk/NS8F7N6fdipElczwndK7MqwATmuxhjwqKnP6q6ItpXu7DaXz3+4dDRnq/ZGFM73jyamSod1zPkOVaISF3qS509kry+Pg4HmQpysFYdTbLvVf8xOf9ejJTM9TTG+CZ6r8yqABNHDdY/ZqhA84JnnOHQ3f5OVQBXLgecEJO5noFO6mb/93o9x4p3fVs4/KnSBbROYbmyUdz3qr4XJyze9Wya6L0y6/pgYihn7Ce9vSTQvCAinSLi8f+ozRGWpK+nv722Edjsz7SgOePGinl9/f+w5f4mjCrt4B9XzGup78WkxHtvTvheOVNqMDDJ4OC/cJ5xD8wdSV1P/0i/yhSXZSaKeX1Dgop+2k5M1Gup78Wkxf3fn8i9cqbUYELbrQNmo+3XydLrmV56fVNHr2VqpfR6zpQA087YqFuOtl8nS69neun1TR29lqmV0us5IwKMf/hce8QaMTVoE0NS9Hqml17f1NFrmVqpvp7ZmipmFVakbQ20V/uH19VjDalzAW2pTCMzE+n1TC+9vqmj1zK1pup6ZlWAUUoplT1mRBOZUkqpzKMBRimlVFpogFFKKZUWGmCUUkqlhQYYpZRSaaEBRimlVFpogFFqihljXMaYxukuh1LppgFGqUkKBAxjTH1IOvN4GghJveE/t8MYI8aYptAVBv3P2ep/rCXWWjxKZSKdaKnUJBljOoA6rMBRKyLV4x0feYw/cDSKyKwox1cBHcCsCa4oqtS0mknp+pWacv6bv0tEuvxrjsRNCug/vn1KCqfUNNMAo9TkXII/EeB4K336NQBNaS2RUhlC+2CUmpxaJpbKvEYTMapcoTUYpZJgjFmBtVpiFbDAGFMNNMULHv7O+5Skkfc3tW0GVmNlvQUr820j2lejMoQGGKWSICJr/Df5ehGpS/C0BmBlnMfL/YErUrRlfyuAJSFrpGOMaQVWanBRmUIDjFLJq+Fw7SER5SIS73hfYE2OUIFAFvlchNSG/KPQKqKdr9R00QCjVPKqgYT6U/zzY1pS+NptgZqKf/XBRn95lMoY2smvVPJqgC0JHtsA3JOqF45oBmvBahqbSG1KqbTTAKNU8qpIoNPevwRtZFBIiUCfTegQaX+TmlLTTpvIlEpC4Cae4JDjRaRh7ou/aWwVIU1j/n0VqX4tpZKhNRilkjORDv660NFeKRStacwNdKfhtZSaMK3BKJWchDr4/TWKuE1j/szKtVjDlJuAFhFp8z9Wj5XnDOB2Y8xGEfH497uAbv8Aggp/meqJPqxZqSmnyS6VSoI/weXq8Wom/j6SzkDAUCqXaIBRKkH+moJPRNqMMSIiJoFzxmROVipXaB+MUom7Hajyp3wZd0KjZk5WuU77YJRKXCDNywIRiZfyJeASNHOyymHaRKZUmhhjWiaQp0ypGUcDjFJKqbTQPhillFJpoQFGKaVUWmiAUUoplRYaYJRSSqWFBhillFJpoQFGKaVUWvx/IvVYjnM0O3QAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "plt.semilogx(freq_vec_star, -np.imag(Z_exact), \":\", linewidth=4, color=\"blue\", label=\"exact\")\n", "plt.semilogx(freq_vec, -Z_exp.imag, \"o\", markersize=10, color=\"black\", label=\"synth exp\")\n",