diff --git a/tutorial/ex4_real_experimental_data.ipynb b/tutorial/.ipynb_checkpoints/ex4_experimental_data-checkpoint.ipynb similarity index 98% rename from tutorial/ex4_real_experimental_data.ipynb rename to tutorial/.ipynb_checkpoints/ex4_experimental_data-checkpoint.ipynb index 671d72b..7b24649 100644 --- a/tutorial/ex4_real_experimental_data.ipynb +++ b/tutorial/.ipynb_checkpoints/ex4_experimental_data-checkpoint.ipynb @@ -35,7 +35,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 1) read in the impedance data from the csv file" + "## 1) Read in the impedance data from the csv file\n", + "### IMPORTANT: the frequency value should be sorted ascendingly (from low to high)" ] }, { @@ -44,8 +45,6 @@ "metadata": {}, "outputs": [], "source": [ - "# read the experimental impedance data stored in a csv file\n", - "# IMPORTANT: the frequency value is sorted from low to high in default\n", "Z_data = pd.read_csv('EIS_experiment.csv')\n", "freq_vec, Z_exp = Z_data['freq'].values, Z_data['Z_real'].values+1j*Z_data['Z_imag'].values\n", "\n", @@ -54,7 +53,7 @@ "xi_vec = np.log(freq_vec)\n", "tau = 1/freq_vec\n", "\n", - "# define the frequency range used for prediction, we choose a wider range to better capture the DRT\n", + "# define the frequency range used for prediction, we choose a wider range to better display the DRT\n", "freq_vec_star = np.logspace(-4., 6., num=101, endpoint=True)\n", "xi_vec_star = np.log(freq_vec_star)\n", "\n", @@ -66,7 +65,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2) show the impedance in a Nyquist plot." + "## 2) Show the impedance spectrum as a Nyquist plot" ] }, { @@ -88,7 +87,7 @@ } ], "source": [ - "# Nyquist plot of the impedance\n", + "# Nyquist plot of the EIS spectrum\n", "plt.plot(np.real(Z_exp), -np.imag(Z_exp), \"o\", markersize=10, fillstyle='none', color=\"red\", label=\"experiment\")\n", "plt.plot(np.real(Z_exp[40:80:10]), -np.imag(Z_exp[40:80:10]), 'o', markersize=10, color=\"black\")\n", "\n", @@ -99,7 +98,7 @@ "plt.legend(frameon=False, fontsize = 15)\n", "plt.axis('scaled')\n", "\n", - "# be arful about the axis range, modify this according to real data\n", + "# this depends on the data used - if you wish to use your own data you may need to modify this\n", "plt.xlim(1.42, 1.52)\n", "plt.ylim(-0.001, 0.051)\n", "plt.xticks(np.arange(1.42, 1.521, 0.02))\n", @@ -108,7 +107,7 @@ "plt.xlabel(r'$Z_{\\rm re}/\\Omega$', fontsize = 20)\n", "plt.ylabel(r'$-Z_{\\rm im}/\\Omega$', fontsize = 20)\n", "\n", - "# label the frequency points, should be modified according to real frequency data\n", + "# label the frequencies - if you wish to use your own data you may need to modify this\n", "label_index = range(40,80,10)\n", "move = [[-0.005, 0.008], [-0.005, 0.008], [-0.005, 0.008], [-0.005, 0.01]]\n", "for k, ind in enumerate(label_index):\n", @@ -125,7 +124,7 @@ "metadata": {}, "source": [ "## 3) Compute the optimal hyperparameters\n", - "Note: the intial parameter may be adjust according to different problems" + "### Note: the intial parameters may need to be adjusted according to the specific problem" ] }, { @@ -182,7 +181,7 @@ } ], "source": [ - "# initialize the parameter for global 3D optimization to maximize the marginal log-likelihood as shown in eq (31)\n", + "# initial parameters parameter to maximize the marginal log-likelihood as shown in eq (31)\n", "sigma_n = 1.0E-4\n", "sigma_f = 1.0E-3\n", "ell = 1.0\n", @@ -354,7 +353,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 4e) Plot the imaginary part of the GP-DRT impedance together with the exact one " + "### 4e) Plot the imaginary part of the GP-DRT impedance together with the experimental one" ] }, { diff --git a/tutorial/ex4_experimental_data.ipynb b/tutorial/ex4_experimental_data.ipynb new file mode 100644 index 0000000..7b24649 --- /dev/null +++ b/tutorial/ex4_experimental_data.ipynb @@ -0,0 +1,415 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Gaussian Process Distribution of Relaxation Times. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## In this tutorial we will show use the GP-DRT method to analyze actual experimental data\n", + "\n", + "The impedance data in the csv file named `EIS_experiment.csv`. The file has three columns. The first column is the frequency, the second one the real part of the impedance. The third column is the imaginary part of impedance. To use this tutorial for your own data, we recommend the frequencies go are sorted ascendingly." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from math import sin, cos, pi\n", + "import GP_DRT\n", + "from scipy.optimize import minimize\n", + "import pandas as pd\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1) Read in the impedance data from the csv file\n", + "### IMPORTANT: the frequency value should be sorted ascendingly (from low to high)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "Z_data = pd.read_csv('EIS_experiment.csv')\n", + "freq_vec, Z_exp = Z_data['freq'].values, Z_data['Z_real'].values+1j*Z_data['Z_imag'].values\n", + "\n", + "# define the frequency range\n", + "N_freqs = len(freq_vec)\n", + "xi_vec = np.log(freq_vec)\n", + "tau = 1/freq_vec\n", + "\n", + "# define the frequency range used for prediction, we choose a wider range to better display the DRT\n", + "freq_vec_star = np.logspace(-4., 6., num=101, endpoint=True)\n", + "xi_vec_star = np.log(freq_vec_star)\n", + "\n", + "# finer mesh for plotting only\n", + "freq_vec_plot = np.logspace(-4., 6., num=1001, endpoint=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2) Show the impedance spectrum as a Nyquist plot" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAADmCAYAAADhjzMqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3zcVZn48c+TUmgLwtCrWLlNVH4CLWuagl1d1qVTtCgCuxmBQkEuJguyK5U1pS4LC66wkxXQteyayM3CFmiiQNV2JWllWWmBJhFBFimdaBFQaJsMWAq95fn9cb6TTJK53yd53q9XXjPzvc2ZbzvzfM/5PuccUVWMMcaYQqsqdQGMMcaMDRZwjDHGFIUFHGOMMUVhAccYY0xRWMAxxhhTFAeUugClMnXqVD3mmGNKXQxjjKkYXV1d21V1Wrb7l23AEZE6IALUqGpTOutFpA/oATpUdWmy4x9zzDF0dnbmv+DGGDNKicjWXPYvyyY1L5igqh1AREQCaa4PquqcVMHGGGNM8ZVlwAHm4moqeI81aa73iYi/8MUzxhiTqXINOL5hr6ekuX4y0CsizfEOKiL1ItIpIp3btm3LQzGNMcakq1wDTgQXPDJar6otqhrBNbPVJVhfq6q106Zlfd/LGGNMFso14GxisBbjB9pTrfdqL9Egs6PwRTTGGJOJsgw4qtoG+KPJAF5yACLSnmT9KmISCLxtjDHGlAkZq6NF19bWqqVFG2NM+kSkS1Vrs92/LGs4xhhjRh8LOMYYY4rCAo4xxpiisIBjjDGmKCzgGGOMKQoLOMYYY4rCAo4xpvDCYViyBGbMgHHj3OOSJW65GTMs4BhjCmvtWvj4x2HiRNiwAXbvdo8TJ7rla9eWuoRFFQwGaWlpKfr7luI9hyvb+XCMMaNAOAwXXQSrV8O8eYPLq6vh5pvhzDPh85+Hp55yy8aAhoYG/P7iDmofiUSIRCJFfc94rIZjjCmc5cvhS18aGmxizZsHl18Od9xR3HKVUCAQKHrAWbq0PKYIs4BjjCmclSvhssuSb3P55W67Aunp6aGpqYm2tjaWLl06cKXf1tbGnDlzqK6uJhKJ0NPTg4jQ0NBAT08PHR0dVFdX09DQQEtLC21tbQPr0jl2dXU1TU1NtLS0MGfOHCKRCN3d3cyZM2cgAHR0dDBnzpyBY3R0dAy8R1tbW0bvOfxY0X1XrFhBT08P7e3ttLS00NHRUbBznZKqjsm/OXPmqDGmwKqqVPfuTb7Nnj2q48YVrAh+v3/geTgc1kAgMPC6r69P/X6/9vX1aV9fnzY3Nw/ZNxQKaV1d3Yjt0zl2KBQaeB173ObmZm1sbBzyuqamJuHr1tZWra+vT+s94x2rsbFRQ6GQhkKhkScnQ0Cn5vC7azUcY0zhTJ0KW7cm3+aVV9x2BdDS0kJNzeCEwX6/n9hBe30+H6FQiGAwyKpVq6ivrx9xjNjmL5/Ph9/vp6WlJeWxo9sDQ447efLIqb5i32Py5MlDtvH5fPT29qb1eeIdqxzu3URZ0oAxpnAWLYK77nIJAonceafbrgDC4TCRSGRIM1Jra+uQberq6mhujjtJcFx+v5+wl86d6thz585N65jDg1A0UEVFg0Y6n2f4sXp7e6kuk4QMCzjGmMK56iqX+nzmmfETBzZudAHnqacK8vZz586lu7ubQCCQcJvu7m6WLl1KQ0NDWjf0e3p6CAaDTJ48OeWx8y2dz5PIjh1uXsq2tjbq6kZMiFwU1qRmjCmc6mpYscKlPi9b5tKk9+51j8uWueUrVhQsJbquro7e3t4hzUptbYNzM0YiETo7OwkEAjQ3NxMMBkccI/aGfTS5oL6+PuWxExm+Tz4/TyJ+v78smtashmOMKayFC10N5o474BOfgO3b3T2bRYuK0v+mtbWVW265ZaB5K3oPpKmpiVtuuWXg/kq0xhIMBlm2bNmQeyUdHR1EIhE2bdpEe3t7ymN3dHTw0EMPAa55LPoe3d3dtLa2DmTBTZ48eeB1W1sbfr+f5ubmgdc1NTU0NzfT2dlJS0sL9fX1Cd8z9tjDj3XuuefS29tLS0tL0VOyY9mMn8YYk0BTUxM7duwgFAqVuihlwWb8NMYYUxEs4BhjTBzRZrFoh0yTO7uHY4wxcQQCAbq6ukpdjFHFajjGGGOKwgKOMcaYorCAY4wxpigs4BhjjCkKCzjGGGOKomwDjojUiUhARBozXS8i1kvLGGPKTFkGHBGpA1DVDiAiIoF013vPSzd2gzHGmLjKMuAAc4HoiHk9QE0660XEH7N8BBGpF5FOEenctm1bfktsjDEmqXINOL5hr6ekud6vqgkDjqq2qGqtqtZOmzYt1zIaY4zJQLkGnAgwclq8JOtFJOA1sRljjClD5Tq0zSYGazF+oD2N9b3e/Rsf4BeRGlXtLkZhjTHGpFaWNRxVbcMFjYD3ugNARNoTrVfVbm+7yYxscjPGGFNiNh+OMcaYtNh8OMYYYyqCBRxjjDFFYQHHGGNMUVjAMcYYUxQWcIwxxhSFBRxjjDFFYQHHGGNMUVjAMcYYUxQWcIwxxhSFBRxjjDFFYQHHGGNMUWQ8WrSI+HAToE0G+oAXVfX3Mev/EvhLVb0pb6U0xhhT8dKu4YjIkSLyEPAm8N/ASmAt8DsR+V8Rmett+inghnwX1BhjTGVLq4YjIrOBx4DpuMnPngZ6gQ8CJwCfADaIyNUFKqcxxpgKlzLgiMhE4BFcE9q1wHdUdXfM+gOBxcA3gX8HXixMUY0xxlSydJrU6oFjgKtUtSk22ACo6h5VvQs4BRdsjs97KY0xxlS8dALOXwO/UdWWZBup6lbgHODdfBTMGDNSR0cHHR0dNDQ0EIlESl0cYzKSTsA5AVifzsFU9WXgdODSXApljBmpu7ub1tZWAoEA4IKPMZUknYDzPlyiQFpU9UlV/UH2RTIm/7q7uwkGgyl/pCORCE1NTbS1tdHU1ER3d3fRy5SoDDU1NTQ3NxOJROjt7R0IPMZUinSy1HYAR6Z7QC9TbZaqXpZ1qYzJo+gPek9PT8ptg8Egzc3N+P1+ABYsWEBrays+n69oZUpVhs7OTubOnTtiP2PKXTo1nE3AZ0RkUqoNReR04DbgizmWy5i8CQQCBAIBJk+enHS7SCRCT0/PwA89gN/vT1orSnYfJdm6RGVKpwzRms0tt9yS+MMYU4bSCTgrgGnAd5NtJCKfAVblo1DGlEJnZ+eImozP56O9vT3hPqtWraKlZWQ+TTAYpLe3N69laGpqGngvn89nSQOm4qQMOKr6Q2Ad8EUR+W8ROTm6TkSqRKRGRO4FVgOCG33AmIoTiURG1DimTJmSNHDU19cTiUSGBJ1gMEgoFBpSS8lHGerr6wdqO11dXYRCoYyPb0wppTuWWhB4FJeBtkBE3sUlEkwHxuECzavedp/x/oypONnUShobGwdqH+3t7VkHm1Rl8Pl8A81pljBgKlFaY6mpagQ4DbgC+CUwEfgALmBtBW4GTlTVp71dJP9FNaaw4jVT7dixI+W9H3BBp7m5mblz5+YUbHIpgzHlLu3BO1V1v6o2q2otgwHnYFX1q+p1qvq2t+m9wF/lWjARqRORgIg0prveex0QEWtrMBmrra0dUbuIRCIsWLAg5b7BYJDW1lZ27NhBW1tbScpgTLlLGXBE5NsicqqIDNRavOFs/qiqI0YVUNWtqvo/uRRKROq8Y3UAEREJpFovIjXAAm9ZjYhkf5lpxoyenp6B1GSfz0dtbe2QVOXOzs6UzVcNDQ0DzWihUIhNmzZlHXSyLYMxlSCdezhXAX8H7BCR1cDDQLuq7ilgueYCD3nPe4AaoCPZelVtArq9+Xp6VDV1pwszJnR3d9PR0UFnZyehUIju7m4aG13FONqRsrm5GYDW1lZaWlrw+/309PTw/e9/P2kfnKamJpYuXTqkGS0UCg0sq6mpybhMmZbBmEohqpp8A5HpwNm4cdL+ChgPvAOswQWfNar6p7wWSqQZaFbVbq92s0BVl6azPt72MfvV4wYj5aijjpqzdevWfBbbGGNGNRHp8m6rZCWdtOg3VbVFVRfistIWAz8DzgAeAN4UkZ+IyGUiMi3bggwTwU2HkPF6r0nNF212G7auRVVrVbV22rR8FdUYY0w60k4aAFDVt1V1paoGcZ1BzwYexE1N8H3gdRH5HxH5iogcnUO5NgHRNgQ/MLzn3Yj1IhLyajCQOmAZY4wpsowCTixV3a2qq1X1EmAGMB/4HnAscDvQIyJdWR67DfBHkwW8Wgsi0p5kfbP3ngHAl2o6BWPy4b333mPnzp2lLoYxFSHlPZysDioyF3fP52xVLcsJ2Wpra7Wzs7PUxTAVbsmSJXzwgx/kmmuuKXVRjCm4XO/hpDvSQEZUdROu2evrhTi+MeVgz5493H///Tz99NOpNzbG5C/giMj7gT/HJRYMaapT1f/I1/sYUy5+/OMfc+KJJ+Y0soAxY0leAo6IXAjciRvSpg+IbadTwAKOGXXuvvtuLr3UJrc1Jl35quF8E2gCblLVfXk6pjFl67XXXmPjxo20traWuijGVIyss9SGORS414KNGStWrFhBXV0dkyalnJfQGOPJV8D5L+CzeTqWMWVNVbnnnnusOc2YDOWrSe2rwCMiMh94Htgbu1JVb8rT+xhTck8++SQHHHAAp5xySqmLYkxFyVfAacBNurYd+BAjkwYs4JhR4+677+aSSy4hZgB1Y0wa8hVw/gm4RlVvz9PxjClLO3fu5OGHH+bmm28udVGMqTj5uoczDlidp2MZU7ZaW1s59dRTef/731/qohhTcfIVcO4BLsjTsYwpW9b3xpjs5atJbRJwuYh8GniOkUkDf5+n9zGmZDZv3szmzZs544wzSl0UYypSvmo4HwV+CewB/h8wK+bvxDy9hykTDQ0NpS5CSdx7770sXryY8ePHl7ooxlSkvNRwVPWv8nEcU/6i0yKPNfv27eMHP/gBjz32WOqNw2FYvhxWroTt22HqVFi0CK66CqqrC19YUzYikQgdHR309vbS3t5OKBQa02Pv5auGY5Lo7u4mGAzS0dGRdLtIJEJTUxNtbW00NTXR3d1d9DIlK0MkEsHv9zN58uif2y68fj1XnnQSh1ZVUSXCoePHs6+vjwlvvJF8x7Vr4eMfh4kTYcMG2L3bPU6c6JavXRvnzcKwZAnMmAHjxrnHJUvcclPRVq1aRU9PD/X19SxYsIBQKFTqIpWWqmb1h8tKOzTmecK/bN+jkH9z5szRYmhvb9f29natqanR9vb2pNsGAgENh8NDXvf19RW1TMnKEN02EAjkvUwZ27JF9eqrVadPV62qco9XX+2W57Ktqq658UadBDq+qkpx/cgU0CoRnQS65sYbE5dp6lTVDRvir9+wwa2Pfd81a9yyZcvc8r173eOyZW75mjX5Ow8VpqurS+vq6lJ+b/r6+jQUCmlra6uGQiHt6uoqepnSKUNjY6OGQqGCla0YgE7N4Xc3lxrODgY7eO5I8TdmBQIBAoFAylpBJBKhp6dnSHXb7/cnrRVFIpGs1iUqU7IydHR0EAgEkn6GosmkFpFhjSO8fj11N9zALmBvf/+Qdf2q7ALqbriB8Pr1I8u1fDl86Uswb178cs+bB5dfDnfc4b1ZGC66CFavhptvds1tBxzgHm++2S2/6KLENZ1salMVItoM1dPTk3LbYDBIXV0ddXV1NDY2snTp0qT//wtRpnTK0NPTQ2NjY97LVVFyiVaV/FesGk5UIBBIeqUWrXHEamxs1Pr6+oT7NDc3a3Nz84jldXV1Q2op6ZYpWRm6urq0tbVVW1tb1e/3p7zqLJhMahFZ1DiumD17RM1m+N/4qir98kknjTze9OmpaxZbtqjOmOGeX321q8kkc+21qkuWxD9OprWpCpTqe9PX16d+v3/Isvr6em1tbU26TzbrEpUpnTJUes0mikLXcETkbRFZISIFmR3UOJFIZESNY8qUKfT29ibcp76+nkgkQktLy8CyYDCY9Y3JZGWoqamhrq4OIGmZCi6TWkSmNQ7g/uefH1GzGW5vfz/3PffcyBXbt8PRRycv/1FHue3AJRVcdlny7S+/3G03XBafbTTeK+rs7MTn8w1Z5vP5aG9vT7jPqlWrhnxnooLBYFb/t1OVoa2tjfr6eoCU93FHu3Sa1A7Bdep8UEQsyaCAsvnP3tjYOBB0cgk26Zahrq6Ovr6+/DavZfJDmMmPdBY/6DtVk2w8KO52U6fC1q3Jd3zlFbcdZB6gYmX62UZp81u5X6h1d3ezdOlS5s+fT3V1dVpNhKNZugFEgXOAlZJixEIRuVBEvp1zycYYn883os13x44daWWENTY20tzczNy5c3MKNrmUIWuZ/hBm8iOdxQ/6IWkOyBl3u0WL4K67ku94551uO8g8QMXK5LPleq+ozJXzhVpNTQ3hcJiuri7C4fBATWesSjfgLAeeBYLAihTbVgN/l0uhxqLa2toR/2kjkQgLFixIuW8wGKS1tZUdO3bQ1tZWkjJkJZsfwkx+pLP4Qb9w1izGVyX/WoyvqmLx7NkjV1x1FXz/+7BxY/wdN250AefLX3avMw1QsTL5bNk0v1WIMXuhVqHSDTi9wALgBWCRiNxduCKNHT09PQNVbJ/PR21t7ZAqd2dnZ8qmq4aGhoGrs1AoxKZNm7IOOtmWIWvZ/BBm8iOdxQ/6NbffzvgU93DG9/ez5LbbRq6oroYVK+Dzn4dly1yg3LvXPS5b5pavWDHY+TPTABUrk8+Wy72iqDK9/zMmL9QqWaqsAqAfuN57Pg0XdPYDzQm2vwHYn0smQzH+ipWl1tXVpaFQSH0+nwYCgSHZKsOz0DLtTxAKheJmozU2NibdN1mZitmnIeOsrujrAmapqSbuhzO+qip5P5zYMi5Z4so9bpx7XLIk/meN9sO59lq3fs8e93jttcn74WTy2aqqXP+eZPbscWWNJ9e+QjmIl6UWDodH9BWLfV1TU5My26y+vn7IPo2NjUkz21KVKZsyVCJyzFLLKOB4r2cAL3lBZ3mc7S3gjHXpdkbM9ocwkx/pLH/Qt6xbp18+6SQ9VESrQA8V0S+fdJJuWbcuixOSQiYBKla6ny2bwB67vATp13ahVp6KHnC8ZTOBLV7QuX3YOgs4Y1kmV8O5/hCm+yOd7Q96JUjns+XS3yeTfUfxqAfGKUnA8ZYfBfzWCzqhmOV5CThAHRAAGtNZD/i8ZXWx5Un0ZwGnADK9Gs7lh9CkL5daSroXBYcfXrJmN1M8JQs43rpjgVe8oPMvmqeAEw0c3vN6IJBqvfdY7y0LRZ8n+hstAefVV1/V/fv3l7oYTqYBZIz0li8L2d4rSqfZ88UX3U+J/TuOerkGnJw6cqrqb4HTgD8Cy0TkhlyOF2MuEE2V6gFqUq1X1RZVjfbk8gMjuvSKSL2IdIpI57Zt2/JU1NLp7+/n5JNP5g9/+EOpi+Jkmg2VaVaXyd7ChfDUU66f0yc+4fo5feIT7vVTT7n18aSTfh0KwcEHV1Ta9a5du0pdhDEpnYBzHfCrRCtVdQsu6GwDrgcW56FcvmGvp6S7XkT8QK+qjujS6wWlWlWtnTZtWh6KWVrPPPMMhx12GDNnzizOG6ZKjc2m53y2P4Qmc9XVcNtt8Mc/wr597vG225IH9HTSr1etgi98Ifk2qdKui+iXv/wln/rUp0pdjDEpZcBR1ZtV9dEU27wEzMeNDJ2P2YUiQLJeU8nW16nqmJiS8pFHHuGcc84pzpulMyJAtj3ns/khNMWRTl+hXbtg6dLkxxl+oVHCfj133XUXn/vc5wr+PmakvI2Npqov4O6l9OXhcJsYrMX4geEj8cVdLyJ1qtrkPS+TsfQL55FHHuHss88u/BulOyLAGWdk33PelKd0mj0PO8z9f0gm9kKjhOO6vffeezz44INcfPHFBXsPk1heB+NU1eeAWcD5OR6nDfBHg4aqdgCISHui9d7zkIh0iUhXLu9fCX7zm9+wc+dO5syZU/g3S3dEAJHse86b8pWq2fOSS9K/0CjxuG6PPvooH/vYxzg6VdOvKYxcMg5w46rdkssxSvVX6Vlqt9xyi1555ZX5OViq/hOZ9JfJNhvKVK5Msg1LnAp/+umn68qVKwty7LGAUmapARcCZ+Ua9Ezm8taclk7zRibJAJYEMPZkkm2Yj3HdsvTKK6/Q2dlZnGZoE18u0QrXR+f/cjlGqf4quYbz6quv6uGHH6579uzJ7UDpXplOnpz9iABm7Ehn1INMhjPK86gF3/jGN/SKK67Ien9T+hqOKYHVq1dzxhlnMH78+NwOlO69maOPtmQAk1o62YbpZDLedx/09+c1qaC/v5977rmHSy+9NON9TR7lEq2wGk5JnH766WmPbJtUuvdmpk61EQFMfqS6h7Nli+rEiarnnht/fZb/137+85/rrFmztL+/P6P9zFBYDWdsiUQibNy4kU9/+tOZ7Riv38Obb7or0WSOOgr6+mxEAJMfqfr1fP3r7vGb34y/PstRC+655x4uueQSUkxYbArMAk6FWbt2Laeeeirve9/7MtkpfmLApEnuC5ysiSLaf8KSAUw+pEowaG11Tb3JLl4yTCp4++23efTRR7nwwgvz8AFMLlL01jLlJuPRBWL7PcTeq6muhvp62LHDrX/qqfhf8th7M9E2+nizXRqTrujFyx13uIuW7dvdRU30/9lFFyXff/ioBSk89NBDzJ8/n9EwnFWlE9csl+XOIv3Ab1T1+PwVqThqa2u1s7Oz1MXIyO7du5kxYwabN29m+vTp6e20ZImrjdx888h14bCr+Sxc6L7wwwPJxo3uSjRRMDIm32bMcLXvZP/fwmEXqP74x7QOOW/ePK677jo++9nP5qmQY5eIdKlqbbb7W5NaBVm3bh2zZs1KP9hA8n4P0eaNn/wEWlrs3owpvXQGC/23f4MPfCCtcdhefPFFtm7dmvk9T1MQFnAqSMrOntkkBixc6K4od+2yezOm9FIlFdx2m7s4OuWUtFKm77nnHi666CIOSDXWmymOXFLcsLTootm3b5/OmDFDX3755fgbJJraedIkNxtjsiFlrNOmKSeJhkdqaFAVUb311vj7bdigWw49VK847jh9n4gCKqDnfeQjumXduuJ+hlGKEqdFvwK8luMxTBqefvpppk2bxoc+9KGRK5MNiFhfD5/7XPIBEa3TpikniTIin3kGGhrgq1+Nu9va9nZmv/02d770En/y7k0r8MMtW5g9fz5rb7qpiB/CxJNT0kAlq7SkgcbGRiZMmMBNF1/s0kZXrhzM7pk50zUx/Od/jtzREgPMaJEkoSC8fj2z588n2Tyek4Dn1q2j+rTTClbEctTR0UEoFKK9ffgsL5mzpIExQFV5+OGHOXvatPj9aV56CR56KH5/GksMMKNFkkFkb12yhL1VyX/O9lZVcXuC2tFoFgiUz9RgFnAqwIsvvsjud97hYzfeGL/Z7L334NFHEzebWWKAGQ2SjMN2//PPs7e/P+nue/v7ue+559J+u+7uboLBIB0dHUm3i0QiNDU10dbWRlNTE93d3Wm/R6YSlamYZciFpW5UgIcffpizZ8xAFi6MP9Dm1KkuTTQ65Ee8jpnjx8P06Wn3XTCm7ERTpuP0KduZ5q2BnaquaW7RIpcRl6BmH/1B7+npSXnMYDBIc3Mzfr8fgAULFtDa2orP50uxZ2aSlalYZciV1XDKUHj9eq486SQOraqiSoTrr7uO3z33HOHjE/SvjX4Rkw35YYkBptIlSZk+JM0x0g4RSWvk6UAgQCAQYPLkyUmPF4lE6OnpGfihB/D7/UlrRZFIJKt1icqUTRlKxQJOmVl7003Mnj+fO3/9a/6kiuJyz9f29zN78eL4mTbRL+Lrr8cf8sOmdjajQZJx2C6cOpVUk3WMr6pi8ezZeZ3OurOzc0QtwufzJb1Bv2rVKlpaWkYsDwaD9Pb2FqUMpWIBp4yE16+n7oYb2AUj2qP3AbuAuhtuILx+/dAdo1/Es86CCRMsMcCMXglSpq859dTUAae/nyWxzc1ZjjwdKxKJjKhxTJkyJWngqK+vJxKJDAk6wWCQUCg0pJaSrzK0tbXR09NDS0tL0hpUMVjAKSM5ZdosXAhf+AIcd5wlBpjRLc5Eb9VtbbTdeCOTcDWZWOOrqpgEtN1448iU6DxMZ51NraSxsXEg6OQSbNIpQ11dHeFwmPr6+pLf07GAU0ZyyrTZuBF++ENYtSr5jIvGjFILr7+e59ato37WLA4VoQo4VIT6WbN4bt06Fl5//eDG0WGg5s2DN95IOh5bMj6fb0StYceOHSnv/YALOs3NzcydOzenYJNLGYrNAk458P7zZ5RpY81mxoxQfdppLH/2Wd7q72f/9Om89fLLLH/22aE1m9j5oR58cLBDaRZTWNfW1o6oXUQiERYsWJBy32AwSGtrKzt27KCtrS3t98xnGYrNAk6pxfznPyTNXQ4BazYzJpV4I08PHwaqvd1tl0EiQU9Pz0Bqss/no7a2dkiqcmdnZ8rOlg0NDQPNaKFQiE2bNmUddLItQ0nkMhBbJf+VxeCdW7a4QQo3bFBV1Stmz9bxVVWKGwIq7t/4qir98kknlbjgxlSAYd8vVVW9+mo3wK2qWz51qtsu1rXXateiRRoKhdTn82kgENBQKDSwurGxUevr6wde9/X1aSgU0tbWVg2FQtrV1ZW0WKFQSMPh8IjljY2NSfft6upKWKZMy5Atchy808ZSK6Vhk6PZeFDG5Nnata7Gcvnl7m/ePNeM1t7uugqsWDGyZSDDCd7GklEzlpqIfEZEXhKRLSJybZz1B4nIQ976p0XkGG/5FBH5uYjsFJHlxS53ToZNjlZ92mnZZdoYY+Ibnka9bZtrQkvWDJ3hFNYmfWURcERkHHAHsBA4HjhfRIZ3q78M6FPVDwG3AyFv+XvAPwH/UKTi5kc47CZH+/M/HzJr4cILLhiZaQPxM22MManFplFPnw5PPpk8e/N//xcOOiitGUVNZsoi4AAnA1tUtUdV9wAPAmcN2+Ys4Afe8zZgvoiIqr6jqr/ABZ7KEE0UmDQJHnhgxKyF1bt3D2babNnCWzNmjMy0McZkLtUU1mvXuvmjTjghrRlFTWbKJeDMBH4f8/pVb1ncbVR1H/AWMCWTN9IcjkUAABNFSURBVBGRehHpFJHObdu25VDcHMRmydTXQ0fH4KjP8bJkbAw0Y/In2RTW4TCcf76r3TzwgPtObt3q5p+66y7YsQM++1m45JKcaztf/OIXy3ZE50Iql4ATb+S94dkM6WyTlKq2qGqtqtZOmzYtk13zZ/ly+NKX3M3LeP/5Y4fbsDHQjMmvJOOxcd55sGePu7daXT20v86GDW5dfT288EJOtZ3XX3+dRx99lOOOOy7PH678lUvAeRU4Mub1B4HXE20jIgcAhwGZjylRCtFezTNmwLe/7SZCW7LErYv3n3/BAreNdeY0Jv8STWH9wgtuosKFCxNP2/61r8Err8B//Af8zd/AlCkZ3+e57777qKur4+CDDy7Chy0v5RJwNgEfFpFjReRA4Dxg9bBtVgMXe8/rgPVa7jnd4TCccw58+MPwne/A/v1u+QMPDLYJw8j//Oef7yZLs86cxhRGnPHY2L0bTj3VrY9tiYh11FEu0+3KK+HP/swNmJvBfR5V5e677+bSSy8t4IcrY7l04snnH3AGsBkIA//oLbsJ+Lz3fALQCmwBngH8Mfv+Dlfb2YmrCR2f6v0K3vFzzRrVww9XnThRddUq1b17XQezSZPc8jVrEnc827JFdcaMwpbPGDPU9OmD38XY57HWrVMVcd/deN/TRN9pzy9+8Qs97rjjtL+/P8+FLw6s42d2CtrxMxx2VzqnnTaYDBC1ZIm7+bh2ravB3Hmnu0KKHTZ92bKRy4wxhRXbEXvcOPcdPGDYpMhz54KIa6X4znfgu9+Fqio36250FtF432nPZZddxnHHHUdjY2ORPlR+jZqOn6NGOAzBoGsSW7Vq8H5NtG33qqtcsFm40CUGDB8e3RIFjCmN2CSeqVNdhlqsjRuhqwv++q/dBeXevW674U1qH/5w3CkPdu7cyY9+9CMWL15cpA9Ufizg5FM0q+Wll+DHP3ZXQk8+ObRtN5ol85OfuGC0d6/r1WyjPhtTWrEZbEcdBd/61sgR2cEtX70aVOHII2HmTPjIR1zq9GmnwTXXxB2poK2tjU9+8pMcccQRRf5g5cMCTr6sX++ufPbvd7Wb8893s2++9trI/jULF7orol274JOfdPvYqM/GlF40g23WLGhuHjki+8SJrgVj40Z3wXjKKUM7iFZXw7vvjmyKg7GdLBCVyw2gSv7La9LAmjUuOWDePHezcNo0d3Nx7ly3fM0at92116ouWeKeR284xi4zxpSPNWtcAsC117rv6549qgcdpPrRj7rEgVtvjb/f2WerwpDEgc2bN+v06dN19+7dRSp8YZBj0oDVcHIVzdefMAHuu89d4VxwgRtB4IEH4MAD3c3EcHjo/Zo774RPfcru1xhTruL119m9GzZvdunQ8aZ637gRfvEL9/yOOwYW33vvvVxwwQUceOCBRSp8mcolWlXyX95qONH5NaqqXOqz6tB5OKK1n7lzVV98UXXcOJcmPXHiYHq0MaYyTJ+uesgh7rsbW/PZssW9njpV9a673KOXMr1v3z6dOXOmPv/88yUufO6wGk6JRacYiM1qib35+MQTrq33hRdcJ7L9++Hcc+HTn4ZNm+x+jTGVZNEi2LnT1WSGj1Swezfcfjt8/esuaeCNN6Cqisc+8hE+cMghnHjiiaUufclZwMnV9u2up/LMmTB79uAwF489Bg895P4TXnONSxDYswdqauDll+Hhhy0TzZhKc9VVLvv0uedGjlTg88Hixa6f3be+5X4HOjq45623uOSll+Cmm0pd+pKzgJOr973P1VxOOcWNMvv444M5+eee68ZF27ABJk920xGsWmWBxphKVV3tWi4uvnjo+If33w833ACHHOIyUrdvh0WL2D57Nj/bu5fzb78d/vmfXTbrGGYjDeRi/Xo4/XQYP97VXg49FN57z/1n/NrX4PXX3c3Fo4+G3/wGfvQja0IzptKFw27EgdNPdxeY27e7Ws+0aW7ytjffdEHpqaf495/+lKeeeoqVK1e6C9M9e1zn0QqV60gDJb95X6q/nJIGtmwZTH0ElyIZDLpU6L/9W5cQ4PO5BIEJE9zjunXZv58xprwMT5kG1fvu0y3nnadXjBun7wMV0CrQM489VresW+d+A6qqSl3ynGBjqWUn6xrO2rUu7fm991zN5pln3GgCV13l1i9f7oa2OOss+MIXXBMaQG9lzKRgjElTOOxSn1euhDfeYO3BB1P3zjvsFWFvzO/qAcCBQNvEiSx8913YsqVim9VtLLViio4msGuXCzhvv+3mxfjLv4Tnn3fV6Msvd/1rdu92WWgbNrjtjDGjS8wUB2Gg7p132AVDgg3APmAXUPfuu4TBjWJw991FL245sICTruhc5x/7GBx8MLz4oksE2LPHjZO2eTM8+CAsXQpf+YrLYnntNVcLmjq11KU3xhTQrRMmsDfFNnuB2wH6+11Xin/918IXrMxYwEnH8NEEIhF3dXPRRXD44UPHSYuOJnDUUe5m4p13utx9Y8yodf9776UVcO4D1yfvhBNcltsYq+lYwElHdPa/t95yGWfRTp7R4czBBZo77hgMNK+8AocdZkPXGDMG7Mxku+XL4Q9/cAsuuwwuuSStqalHAws46Rg+msCiRW4o8tgRBfr6XO0nHHa1ngsvdPd6bKoBY0a9Q0TS2w7cKNQzZgwuvP/+MXNfxwJOOrZvdzWbaKCJnagpOsDfQQe57Y4/Hv70J/jVr+CnP7V+N8aMARfOmsX4quQ/p+OBxUcc4ZrmJ0xwfXeqqtwIBXV1rpVklAcdS4tOx4wZLtsMXILA6tXuPs5FF7n/JJdf7nobf/zjLolgwgT4r/+yYGPMGBFev57Z8+ezK8k2k4DngLjtHevWwb33Qmsr/PrXZdsqYmnRxRCvCe2JJ9xYabt2uWFtPvpRd4/HBuU0ZsypPu002m68kUnA+GHNa+NFmITrh5MwjAQCLut13Lgh0xqMNlbDSUc4PFizmTdvaIev7dsHh7T5yU/cFLPGmDEpvH49t3/1q9z3q1+xE3fPZnFVFUv6+xMHm+GmToVt2wpXyBxYDacYYms2y5a5ZaGQG2Hga19zVyU//KEFG2PGuOrTTmP5s8/ylir7gbeA5SLpBxtwF7Fr1xamgCVmASdd8Wb/i53r3JrQjDGxpk93I8Rn04p05pmjMlXaAk4mYoayGJgD47bbyvYGnzGmhBYtcvd4x4/PfN/9++HSS/NfphIrm4AjIp8RkZdEZIuIXBtn/UEi8pC3/mkROSZm3TJv+Usi8uliltsYY+KKTta2e3f89ePGuYzWI46Iv/6JJwpXthIpi4AjIuOAO4CFwPHA+SJy/LDNLgP6VPVDuCGJQt6+xwPnAScAnwH+wzueMcaUTnSytkT273fJRvv2wTHHFK1YpVQWAQc4Gdiiqj2qugd4EDhr2DZnAT/wnrcB80VEvOUPqupuVf0tsMU7njHGlNatt8ZfHq3VVFW5JIF3342/3Si7j1MuAWcm8PuY1696y+Juo6r7cAkgU9LcFwARqReRThHp3FamaYfGmFGkutr1sRnu3Xfh5JPdFPVHHglvvBF//1HWJ6dcAk68gYiGp3Yk2iadfd1C1RZVrVXV2mnTpmVYRGOMycL3vjdyWSTiJm986y34vXe9HG9onJUrC1u2IiuXgPMqcGTM6w8CryfaRkQOAA4DetPc1xhjSqO62s0AnEg0bbq/f+jyI45wzW2jSLkEnE3Ah0XkWBE5EJcEsHrYNquBi73ndcB6b47t1cB5XhbbscCHgWdSveFUmxTNGFMst94avwaTTFNTOU7emFMELJuhbUTkDODbwDjgblX9pojcBHSq6moRmYCbv+hjuJrNeara4+37j8CluNlcr1bVlN10ReTXwHuF+TQVZyo5/kcaRexcDLJz4eTlPMyE978/wf3l4XbCW+/CuwKy1bXilIsJqnpitjuXTcApNhHpzGVMoNHEzsUgOxeD7Fw4dh4G5XouyqVJzRhjzChnAccYY0xRjOWA01LqApQROxeD7FwMsnPh2HkYlNO5GLP3cIwxxhTXWK7hGGOMKaJRF3BE5G4RedNLe0623VwR2S8idd7rPxORjSLygog8JyLnFqfEhZPtuYhZfqiIvCYiywtb0uIRkYCItKexXSjmeY2I1A0/P5Uuy3NR5+1XX9jSFVeqcyEifSLSleBcNBanlMWR6bkQEV/0+xF7fuIZdQEHuBc3anRC3mjSIeBnMYt3ARepanTU6W+LiK9QhSySe8nuXER9A/if/BerdFS1I9U2IhIA/DGLGlS1DfCLiD/BbhUn03PhPe/x9usRkZoCF7Fo0jgXQVWdo6pLwQWbmP0i3rkZFTI9F8AXgMned4RkFyOjLuCo6hO4jqHJ/B3wQ+DNmP02q+rL3vPXvXUVPeBatucCQETmADOAxwpTuvLkBZSemNf1QJeI+FW1KdrZeCwYfi6ATqDVCzR+Ve0uTclKwjfsYmMug+emBxg1wTcNQ86FN0ZlNJnADyQMWKMu4KQiIjOBc4A4I+oNbHMycCAwusYGHybRuRCRKuBW4GulKFeJ+YcFlWrvr1dEmkdBrTcTQ86FqkaAZqAVd07Gksl4/we818P/H0wpcnlKafi5AAYuUHqTXZSNuYCDGz5nqaruj7dSRI7ADaFziar2x9tmFEl0Lq4E1qjq7+PsM2qJSCBBc0LY+7HtAkbVvYtE4p0LrxmpQ1WrY16PCd5VfATXfFYHRHA/vGNOnHMRVaeqDcn2PaCwRStLtcCDbu42pgJniMg+VX1ERA4Ffgpcp6pPlbKQRRL3XADzgL8QkSuBQ4ADRWSnqo6Y+nuU6fXa4n24+zU1uIFloz8sPtwPzVgQ71z4o+30wC24tvtRz2tW7fU++w5v8SYGazl+IGXyxWiQ4FwgInWq2uQ9T3ThNvZqOKp6rKoeo6rH4GYOvdILNgcCDwMrVLW1pIUskkTnQlUvUNWjvOX/gDsnoyLYeFdktbFXZtGMHFXt9r4ok/F+TLwvli96UzimrbriZXougBZvEsMA8IWxci6AVcQkBqhqW0wSSXRZygSMSpHpufCeh7zMta6kxx5tHT9F5AHgU7gr9jeAG4DxAKo6/F7FvcBPvJN2IXAP8ELMJl9U1WeLUOyCyPZcDFv+RaBWVa8qfImNMaPZqAs4xhhjytOYa1IzxhhTGhZwjDHGFIUFHGOMMUVhAccYY0xRWMAxxhhTFBZwjDHGFIUFHGOMMUVhAceYLIjIP4qIpvm30xsQNV/vPcObv+jfY5aJiPyNiDwqIn8QkT0iskNEnhCRq0VkUr7e35hsjcWx1IzJh/8Dbkyy/nDgKtxF3aN5Hgj2LO+4DwOIyOG4IUcCwHbgv4Hf44akOQ24HfiKiJytqr/KYzmMyYiNNGBMnonIFNxgjh8DHgHOVdU9eTz+Wtx8LDMAAdYBpwI/AK5S1Z0x2wrQAHwX6APmqurWfJXFmExYk5oxeSQi04Cf44JNG252xHwGm8NwtZYfe9NKXIYLNh24KTV2xm6vzveA63ATCn4rX2UxJlMWcIzJExF5P/A4MAt4ADhPVffl+W0+i5sc8Efe6y95j9dr8uaK7+CmVjhHRMbkPC6m9CzgGJMH3uyp/wMcD6wALkw0yV+OzgHeAdpF5CDc1MbvAk8n20lV3/O2GYdrjjOm6CzgGJMjETkKF2w+AtxFgWaLFZEJwGeAtV4AmYy7hxNJ8/16vUer4ZiSsIBjTA5E5FhcsKkGvgd8qYBTk5+Om4H1Ye/1W97jdBFJJ+P0g95jX74LZkw6LOAYkyUR+RAu2BwDfFdVr0hxH2X4/gdm+JbnAHtw06CjqruA3+Gayf4sxXsdBJzgvXwxw/c1Ji8s4BiTBRE5DhdsjgRuVdW/T2Ofx0XkP0XkWyKyDXjS67DZKCJhEXlXRJ73Zp8dvu844Exgvaq+FbNqpfeY6v0vxDWlbbS0aFMqFnCMyZCInIALNh8A/lVV/yGD3S/E3Xf5C+Ai4F9wqc1fxiUc3AI0i8hnh+13KjCFwea0qCbgVWCRiJyRoLxHA98E+oGvZlBWY/LKOn4akwERmY3r8zIN+IaqXp/Bvo8Dk1V1tvf6YNzIAKer6v/GbPdt4COqekbMsu8CVwIfUNU3hh3348Bj3sugqv4sZt2HcUHqBOArqvrvGFMiFnCMSZOX+vwcrmnqt7j052TCqnpfzP6PA79V1Uu813OBZ4BdQOwXcTzwO1U9Lmbf33vL/iJm2THAF72XJwMLvefXqOptIvJJ4GfAJNx9m1Xe+sdV9fF0PrMx+WRjqRmTvk8ymFJ8LHBDiu2/C9w3bNk7Mc+jTdpnAq8M225v9IkXmD6IGxMt1jEJyjDbe/wQLtgAfHTYto8nLrYxhWEBx5g0qepDwEN5POT/AbuBo1V1fZLtzvEeh9y/8WopkmgnVb0XuDenEhqTRxZwjCkRVf2TiHwL+JY3yOYTuH42Hwf6VbXF2/Qc4Feq+tsSFdWYvLCAY0xp/RPwBvAPwH8CbwPP4rLPAFDVj5amaMbklyUNGGOMKQrrh2OMMaYoLOAYY4wpCgs4xhhjisICjjHGmKKwgGOMMaYoLOAYY4wpCgs4xhhjisICjjHGmKL4/2Tooim0yq1kAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Nyquist plot of the EIS spectrum\n", + "plt.plot(np.real(Z_exp), -np.imag(Z_exp), \"o\", markersize=10, fillstyle='none', color=\"red\", label=\"experiment\")\n", + "plt.plot(np.real(Z_exp[40:80:10]), -np.imag(Z_exp[40:80:10]), 'o', markersize=10, color=\"black\")\n", + "\n", + "plt.rc('text', usetex=True)\n", + "plt.rc('font', family='serif', size=15)\n", + "plt.rc('xtick', labelsize=15)\n", + "plt.rc('ytick', labelsize=15)\n", + "plt.legend(frameon=False, fontsize = 15)\n", + "plt.axis('scaled')\n", + "\n", + "# this depends on the data used - if you wish to use your own data you may need to modify this\n", + "plt.xlim(1.42, 1.52)\n", + "plt.ylim(-0.001, 0.051)\n", + "plt.xticks(np.arange(1.42, 1.521, 0.02))\n", + "plt.yticks(np.arange(0.00, 0.051, 0.01))\n", + "plt.gca().set_aspect('equal', adjustable='box')\n", + "plt.xlabel(r'$Z_{\\rm re}/\\Omega$', fontsize = 20)\n", + "plt.ylabel(r'$-Z_{\\rm im}/\\Omega$', fontsize = 20)\n", + "\n", + "# label the frequencies - if you wish to use your own data you may need to modify this\n", + "label_index = range(40,80,10)\n", + "move = [[-0.005, 0.008], [-0.005, 0.008], [-0.005, 0.008], [-0.005, 0.01]]\n", + "for k, ind in enumerate(label_index):\n", + " power = int(np.log10(freq_vec[ind]))\n", + " num = freq_vec[ind]/(10**(power))\n", + " plt.annotate(r'${0:.1f}\\times 10^{1}$'.format(num, power), xy=(np.real(Z_exp[ind]), -np.imag(Z_exp[ind])), \n", + " xytext=(np.real(Z_exp[ind])+move[k][0], move[k][1]-np.imag(Z_exp[ind])), \n", + " arrowprops=dict(arrowstyle=\"-\", connectionstyle=\"arc\"))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3) Compute the optimal hyperparameters\n", + "### Note: the intial parameters may need to be adjusted according to the specific problem" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "sigma_n, sigma_f, ell\n", + "0.0373346 0.0045175 1.0000012\n", + "0.0003046 0.3966704 1.0003612\n", + "0.0003029 0.3966881 1.0004665\n", + "0.0003008 0.3965080 1.6176298\n", + "0.0002903 0.3961564 3.0251495\n", + "0.0002962 0.3957943 3.8833941\n", + "0.0003193 0.3955872 3.7220880\n", + "0.0003122 0.3956723 3.6404464\n", + "0.0003056 0.3957285 3.5879491\n", + "0.0003067 0.3957004 3.6022622\n", + "0.0003071 0.3955843 3.6067656\n", + "0.0003091 0.3944528 3.6278175\n", + "0.0003113 0.3922271 3.6506848\n", + "0.0003141 0.3880471 3.6778552\n", + "0.0003179 0.3801509 3.7139034\n", + "0.0003467 0.3043847 3.9826235\n", + "0.0003572 0.2720155 4.0972295\n", + "0.0003689 0.2256385 4.1635160\n", + "0.0003861 0.1250404 4.0463634\n", + "0.0003902 0.0925561 3.9641940\n", + "0.0003909 0.0663511 3.8302884\n", + "0.0003866 0.0398763 3.6085232\n", + "0.0003716 0.0206307 3.2806321\n", + "0.0003422 0.0131601 2.8723553\n", + "0.0002875 0.0126253 2.2688182\n", + "0.0002737 0.0072530 2.0055273\n", + "0.0002861 0.0054548 1.8126150\n", + "0.0002953 0.0062130 1.7550819\n", + "0.0002955 0.0063377 1.7755749\n", + "0.0002953 0.0063697 1.7758190\n", + "0.0002953 0.0063749 1.7756881\n", + "0.0002953 0.0063749 1.7756463\n", + "0.0002953 0.0063749 1.7756444\n", + "0.0002953 0.0063749 1.7756443\n", + "0.0002953 0.0063749 1.7756443\n", + "Warning: Desired error not necessarily achieved due to precision loss.\n", + " Current function value: -577.496686\n", + " Iterations: 35\n", + " Function evaluations: 167\n", + " Gradient evaluations: 155\n" + ] + } + ], + "source": [ + "# initial parameters parameter to maximize the marginal log-likelihood as shown in eq (31)\n", + "sigma_n = 1.0E-4\n", + "sigma_f = 1.0E-3\n", + "ell = 1.0\n", + "\n", + "theta_0 = np.array([sigma_n, sigma_f, ell])\n", + "seq_theta = np.copy(theta_0)\n", + "def print_results(theta):\n", + " global seq_theta\n", + " seq_theta = np.vstack((seq_theta, theta))\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", + "\n", + "# minimize the NMLL $L(\\theta)$ w.r.t sigma_n, sigma_f, ell using the BFGS method as implemented in scipy\n", + "res = minimize(GP_DRT.NMLL_fct, theta_0, args=(Z_exp, xi_vec), method='BFGS', \\\n", + " 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" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# 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)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4b) Factorize the matrices and solve the linear equations" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "\n", + "# Cholesky factorization, L is a lower-triangular matrix\n", + "L = np.linalg.cholesky(K_im_full)\n", + "\n", + "# solve for alpha\n", + "alpha = np.linalg.solve(L, Z_exp.imag)\n", + "alpha = np.linalg.solve(L.T, alpha)\n", + "\n", + "# estimate the gamma of eq (21a), the minus sign, which is not included in L_im_K, refers to eq (65)\n", + "gamma_fct_est = -np.dot(L_im_K.T, alpha)\n", + "\n", + "# covariance matrix\n", + "inv_L = np.linalg.inv(L)\n", + "inv_K_im_full = np.dot(inv_L.T, inv_L)\n", + "\n", + "# estimate the sigma of gamma for eq (21b)\n", + "cov_gamma_fct_est = K - np.dot(L_im_K.T, np.dot(inv_K_im_full, L_im_K))\n", + "sigma_gamma_fct_est = np.sqrt(np.diag(cov_gamma_fct_est))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4c) Predict the imaginary part of the GP-DRT and impedance" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize the imaginary part of impedance vector\n", + "Z_im_vec_star = np.empty_like(xi_vec_star)\n", + "Sigma_Z_im_vec_star = np.empty_like(xi_vec_star)\n", + "\n", + "gamma_vec_star = np.empty_like(xi_vec_star)\n", + "Sigma_gamma_vec_star = np.empty_like(xi_vec_star)\n", + "\n", + "# calculate the imaginary part of impedance at each $\\xi$ point for the plot\n", + "for index, val in enumerate(xi_vec_star):\n", + " xi_star = np.array([val])\n", + "\n", + " # compute matrices shown in eq (18), k_star corresponds to a new point\n", + " k_star = GP_DRT.matrix_K(xi_vec, xi_star, sigma_f, ell)\n", + " L_im_k_star = GP_DRT.matrix_L_im_K(xi_vec, xi_star, sigma_f, ell)\n", + " L2_im_k_star = GP_DRT.matrix_L2_im_K(xi_vec, xi_star, sigma_f, ell)\n", + " k_star_star = GP_DRT.matrix_K(xi_star, xi_star, sigma_f, ell)\n", + " L_im_k_star_star = GP_DRT.matrix_L_im_K(xi_star, xi_star, sigma_f, ell)\n", + " L2_im_k_star_star = GP_DRT.matrix_L2_im_K(xi_star, xi_star, sigma_f, ell)\n", + "\n", + " # compute Z_im_star mean and standard deviation using eq (26)\n", + " Z_im_vec_star[index] = np.dot(L2_im_k_star.T, np.dot(inv_K_im_full, Z_exp.imag))\n", + " Sigma_Z_im_vec_star[index] = L2_im_k_star_star-np.dot(L2_im_k_star.T, np.dot(inv_K_im_full, L2_im_k_star))\n", + " \n", + " # compute Z_im_star mean and standard deviation\n", + " gamma_vec_star[index] = -np.dot(L_im_k_star.T, np.dot(inv_K_im_full, Z_exp.imag))\n", + " Sigma_gamma_vec_star[index] = k_star_star-np.dot(L_im_k_star.T, np.dot(inv_K_im_full, L_im_k_star))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4d) Plot the obtained GP-DRT" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEUCAYAAACVjRnNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXRb130n8O992Bdi4SaSoigKpCRrs2SS3h0vkWy3sePasRQ3cdJMmrGUnJM5045TyZ7OTNPWaSK5c2aSnDaW3Kynk9Sm7LipY9cW40Z27HihaMvULnERRVGiKBILgYf14c4fWIgHbgAJ8GH5fc7hkd7DA3D1BOL37n2/+7uMcw5CCCGkkAlKN4AQQgiZDwUrQgghBY+CFSGEkIJHwYoQQkjBo2BFCCGk4FGwIoQQUvDUSjeg2FRXV/Pm5malm0EIIUXlyJEjVznnNQt9PgWrLDU3N6O7u1vpZhBCSFFhjJ1fzPNpGJAQQkjBo2BFCCGk4FGwIoQQUvAoWBFCCCl4FKwIIYQUPApWhBBCCh4FK0IIIQWPghUhhJCCR8GKEEJIwaNgRQghpOBRsCKEEFLwKFgRQggpeFTIlsyoq6sLnZ2daGlpgc1mAwDs3LkT+/btw+7du9Hf34+9e/fiwIED2L59O+6++264XC709fWhpaUFu3fvnvW1+/v7sX//fuzbt0/23PHxcQDA3r17Zzx2586daGlpAQD09fXh7rvvxvbt2wEAPT092LFjB7Zv346qqip88MEH6OrqwpNPPgkAOHTokOxPQkiR4ZzTTxY/7e3tvNTt3r2bb9++fdr+7du3823btsn2AeBHjhyZdtzOnTvnfZ+ZntvZ2cnb2toyOratrY3v378/+bxDhw7JXsdms8mOz6RNhJD8ANDNF/HdS8OARKarqwsHDx5EZ2fntMcSvZT57Nq1CwcOHFjQ+2/fvh0dHR3Yt2/fvMdu27ZN1gvbtm2b7PHKykrZdnt7+4LaRAhRHgWrfGKsMH6ysGfPHuzatWvGx9ra2uBwOHJxZua0Y8cO7NmzZ97jXC4X2traACCjdnV0dCy6bYQQZVCwIjI9PT3JADCT1J7MbDo7O7Fz584FtyERVPr7+2d83OVy4cCBA5iYmMCzzz4LAHO2OSGTYwghhYkSLEhSIjikD5+lSiRbpOrq6kJ/fz8mJibQ19cHm82WUVCb7z36+/tlPabU9+ns7MT+/ftnbA8hpPRQsCJJicAwMTEh29/f34+uri7Zcan3h7Zt25bTXovL5QIwPWimvk9HRwfa29vhdDpz9r6EkMJFwSqfOFe6BVlra2tDT0+PLBg5HA7s3LkTO3bsQH9/P44cOZLx6yWek5DJcxPHzxUA29ra4HK55h22JISUBrpnRWT27t2L/fv3z/iYw+GYc4hwJp2dnThy5EjyJxPPPffcnPO0Us12X4sQUlooWBGZbdu2Yfv27dixY8e0x2YLDOnDhovR09Mjm8w71/s4HA588MEHACAbpsxHuwhJ8Pl8cLlciEajSjelrNAwIJlm79696Orqwq5du5IVLFwuF5599ll0d3cDmKosASD5Z/o8p9mkP7e9vV1WwSK1B5Y4dvv27cm5X4n36ezsxJ49e3Dw4MHk/bb+/n4cPHgQzz33HFwuF/bs2YOWlpZFZScSAgDBYBBDQ0MYHh4GYww6nQ6NjY2oqamBXq9Xunklj/EivK+ipI6ODp74wiaElIeRkRGcO3cOjDGYzWYIgoBwOAxRFCEIAtrb22EwGJRuZkFjjB3hnC94smPR9awYYzYAOwH0A3AA6OKc92R7LGOsDUCiK3A9gP2c8+ljSYSQsubxeHDmzBlYLBaoVKrkfo1GA6vVCp/Ph+PHj+O6666TPU5yq+iCFYBOALs45/0AwBg7xBjbwTl3ZXnsNs75vvh+G4ABxtjW2QIfIaT8SJKEU6dOwWAwQB2JoOrnP0fFe+8hsHo1xh98EKGVK2EymZJFnFevXg2WZdUYkpmiClbxoOJIBJ+4fsR6SAczPZYx1g/gSQD7AIBz7mKMdcdfh4IVIQQAcP78efhFEc3vv4/6//t/ob10CQBg+d3vUPvjH8Pb0YGrjzwCvnUrhoeHYbFYUFdXp3CrS1NRBSsAHQDSe1AuAHcjLVjNdSznfBdjLD3dzTHD8YSQMuV2u3Hpww/R/tRTMH/00YzHmLu7Ye7uxtjnPgfp8ceTw4VGo3GJW1v6ii113QYgPR95HMBMk3/mPDb1/hRjzBHf/3zOWkoIKVqSJOH0iRO49q//etZAlarmF79A9WuvQRAEDA8PL0ELy0+xBStg5sC02GP3A9g6y30vMMZ2Msa6GWPdY2NjWbw9IaQYjY+Po+pnP0NFb69sP1erMf7QQ/Bt2TLtOY1/+7eoGh3FyMgIAoHAUjW1bBRbsHIh1mNKVYXpPaiMj2WM7QawZ67ECs75Ac55B+e8o6amJvtWE0KKBuccVw4fxqof/Ui237dlC0698AKG/+qvcO4nP8HZn/4UUZ0u+bgqEMCq3buhCYUwMjKy1M0uecUWrLoxvbdkAzDTWuXzHssY247p6eyEkDLmmZhA0//6XxBCoeS+iMWCwaefRmjlyuQ+cfNmXHziCdlz9f39uOZ738Pw8DBCKc8ni1dUwSo+TNcdv8eU0AGgC4jde0o8lsGx2wC4UgKVI/44IaSMBb71LVhOnZLtu/jEE4jMMKoy8eCDmLj/ftm+ql//GpbeXly+fDmv7Sw3xZYNCAA7AOyMp587ADyWcq9pF2K9p11zHRsPTIcApM+JoHXPCSlj/qEhVP/DP8j2ue+6C64//EMAQCAQgCiKYIxBEASYTCZc/Mu/hPHECehTamc2P/88jm3ejPr6emg0miX9N5SqogtW8cC0b5bH9mRybHzuFc3cI4TIBL77XRhSh/9sNgz/j/+BcCSCyclJVFRU4LrrroNOp8PY2BguXLgAbySC4SeeQGtK/Unbm29Cd/YsxhwONDQ0KPFPKTlFF6wIISQfIqII409/Kts3+pWvIGizwevxYN26daitrYUgxO6eNDU1obGxEWfPnsUIgPqNG2E6diz53FWdnejfsAH19fVU1SIHiuqeFSGE5Ivvpz+FLl75HwAkgwETDz4It9uN1tZW1NXVJQNVgiAIaGlpgcFoxPDnPy97rOrQIUj9/fD5fEvS/lJHwYoQQgCof/AD2bbzgQfgYQx2u33OoTy1Wo3169dj9OabEUjJFmSRCJoOHgTNzcwNClaEkLIXeustmNImAI9+9rOIRCJYs2bNtB5VOpPJhDXXXIOBz35Wtr/+5Zdx5cQJWqgxByhYEULKnvR//o9s23Prrbhit6O1tTXjOn/Lli0De/RRBFNS3IVgELUvvgiPx5PT9pYjClaEkPI2MgL9v/2bbNfFhx+G1WpFfX19xi/DGENjSwsubN8u21936BBGac7VolGwIoSUtcg//RNYJJLcDjQ3Y3TzZqxcuXLe4b90ZrMZoc9/HlH1VKK14cIF+A4fRiTlPUj2KFgRQsoaf16+2MKV7duhMxhgt9sX9HrLN23C+E03yfbVHjoEl4tWIFoMClaEkPJ1+jQ0x48nN7kgYPjWW7FixYqse1UJFosFk5/+tGxf7RtvYGRoaFFNLXcUrAghZSv885/Ltr0dHYhUVqK2tnbBr8kYg+0LX0DEZEru0zqdQFcXFbddBApWhJDylTYEePmOO9DQ0ACtVruol7XV1WHik5+U7Vt26BBlBS4CBass+f1+ujoipBScOAFNSnV1rlLhyq23ZpUBOBtBEKD+T/9Jtq/6d7/D2MDAol+7XFGwylI4HMaxY8cQDoeVbgohZBGkX/xCtu1ub4e5uRmmlOG7xai47z7ZnCtVIAD2q19BkqScvH65oWCVJcYYfD4fjh8/Th86QooV59OyAC/ffjsaGxtz9hYanW56osVrr9FQ4AJRsFoAi8UCj8eDkydPUhkVQorR8eNQnzmT3ORqNSZuvx1WqzWnb6P58pdl2/aeHkz09eX0PcoFBasFslqtuHr1Ks6fP690UwghWeLPPSfbdnV0wN7aCrU6t6smmW++GWJqcVtJgvRv/0YXuQtAwWoRrFYrBgYGMJ6yrAAhpPBFDx6UbY/ecQeWLVuW8/dRqVQI3HuvbJ/t8GF4vd6cv1epo2C1CIIgoKKiAidPnoTf71e6OYSQTAwOQpWWBejMwxBggnbHDtl25fvvw0m1ArNGwWqRtFotGGM4efIkJVwQUgx+/WvZpmfTJthWrYJKpcrL2xnvvBOhqqrkttrvh/jrX4Nznpf3K1UUrHLAbDZjcnISQ1ROhZCCJ6VVWB+/8UbU1dXl7f0EtRrBe+6R7bP+9rcQRTFv71mKKFjliMViwfnz52ksmpBCJooQDh+W7XLecgssFkte31b9mc/ItqveeQdupzOv71lqKFjliCAI0Ol0OHPmDGX6EFKo3ngDLBBIbgbq6mC+4Ya8DQEm6O+7D5LBkNzWXb2KybSgSeZGwSqHjEYjPB4PRkZGlG4KIWQG0Zdflm2P33gjavOQBZiOGQwI3HGHbJ/+3/+dKuFkgYJVjlksFvT19VF2ICGFhvNpyRVLMQSYIDz0kGy76u23MTk5uSTvXQooWOWYSqWCWq1GH81SJ6SwHDsGYXg4uSnpdBC2bs37EGCC7jOfAU95L/PAADwffrgk710KKFjlgclkwtjYGK0MSkghSetVua67DtUrVizZ2wvV1fC3t8v2RV9+me5xZ4iCVR4wxmAwGNDf308fREIKRDQ9Zf2mm5ZsCDDp/vtlm7Z33qEU9gxRsMoTg8EAt9tNpZgIKQROJ9i778p2BbZuhV6vX9JmaB98ULZt++gjeKiaRUYoWOWRyWRCX18fVbYgRGldXWApoxy+Vatg37x5yZuh3rgRoeXLk9tCKAT/K68seTuKEQWrPNLpdAgEAhgdHVW6KYSUNf7aa7LtiRtugM1mW/qGMIZIWjUL/X/8B60+ngEKVnlmNpsxMDBA8ykIUQrn4K+/LtvlvummnK0InC1N2oKMle++i0lakHFeFKzyTKPRIBKJYGxsTOmmEFKeTp+GcOFCclPSaqH95CchCMp8/WnuuQdRrTa5bbh8GZ4PPlCkLcWEgtUSMJlMGBgYoHtXhCghvVe1eTOqc7h8fdZMJoRuvlm+75VXKHN4HhSsloBGo0E4HKbeFSEKiP77v8u2nddfj4qKCoVaE8Puu0+2TSns86NgtURMJhMGBwfp6omQpRQMgqUVjA3ddRc0Go1CDYrR/NEfybatH38M98WLCrWmOFCwWiJarRaBQIDmXRGylN5+GyylxxKsqYHlppsUbFCMsGYNQitXTm1HIvCnFdklchSslpDRaMTAwAD1rghZKmn3qyba22FVImV9BtF775VtGymFfU4UrJaQTqeDKIpUM5CQJRJNm1/luekmGI1GhVojp37gAdl25fvvw+N2K9SawkfBaonp9XoMDQ0p3QxCSt/oKISPPkpucsYg3HOPYinr6dRbtyKaUu5JPzYGz+9/r2CLClth/K+VEb1eD7fbDZ/Pp3RTCCltXV2yzck1a2BfvVqhxsxAr0fo1ltlu9irr9JtgllQsFpijDGoVCpcunRJ6aYQUtJ4esp6RwfMZrNCrZlZegq7/d136UJ2FhSsFGAymTAyMkI3UwnJl2h0Wokl8ROfWPIq6/PRpqWwW3p74abbBDOiYKUAQRDAOcfVq1eVbgohpam3F8KVK8nNiMEA/V13KdigmTGHA6GWluS2IEkIpi0SSWIoWCnEZDLh/PnzND5NSD6kZQG62tpQuWyZQo2Z27QU9sOHEQwGFWpN4aJgpRCNRoNgMEhp7ITkwUwllgrtflWCOr0K+3vvwU3fC9NQsFKQTqfDhZRq0ISQHPD5wN5+W7ZL2roVKpVKoQbNTX3XXZAMhuS2bnwcnrfeUrBFhYmClYIMBgOcTif8fr/STSGkdBw+DJaSvORvaIClrU3BBs1Dp0PoE5+Q7VK9/joikYhCDSpMFKwUlEhjv5JyI5gQsjjTVgW+/npYrVaFWpMZIS2FvfLddzE5OalQawoTBSuFGY1GDA8P01pXhORI+vyqyQIqsTSbaSnsx49j4uxZhVpTmChYKUytViMcDsNNNcEIWbyhIQhnziQ3oyoVNPfeC8aYgo2aH1u5EsG1a6e2o1FIL79M2cIpKFgVAJ1Oh4u0lg0hi5c2EdizYQNsTU0KNSY70fQFGd96i6pZpKBgVQAMBgPGx8cp0YKQRYq++qps29nRofiqwJnSfOYzsu3K99+Hc3RUodYUnqILVowxG2NsN2Nse/zPWdN85juWMdbGGOtkjG3Lf8tnxxiDIAhU0YKQxQiHpxWvDdxxB7RarUINyo765psRrqqa2hZFiGnBt5wVXbAC0AngIOf8IOd8H4C9jLHZVlOb9dh4gKoE4FiSVs8jkWhBY9SELNA770DweJKbIasV5ttvV7BBWRIEhO+5R7bL/B//QSMucUUVrOKBxsE570/Z3Q9gWs9ovmM5512c8y4AE3lscsYSFS0o0YKQheFpy8JP3HgjrJWVCrVmYVQPPSTbrn7nHTgnCuIrSnFZByvG2BbG2JZZHnt4tsdypANAeh0SF4C7F3lsQdBoNLh8+bLSzSCkKEXTgpXz5psLtsTSbLSf+hSiKcOW+tFRON98U8EWFY6MgxVj7BuMMQnAEQBHGGMSY+wfGWPJu5ec8xdih7J8TRqyYXpPaByx4bzFHDsnxthOxlg3Y6w7nz0fo9GIK1eu0NIhhGRrcBCqU6eSm1ylAu69t2BWBc4UM5kQuO022T7Db34DURQValHhyOh/kjH2DICvAngCwD3xnycBtAJwMca+nTiWc/4hgHxOasgm2ORkDIBzfoBz3sE578jnTPjEL9YEdfsJyU7ashruDRtQ6SiI29HZS5sgXP3OOxgfH1eoMYVj3mDFGLsOADjnrZzzpznnv4n/7OOc34NYQOhnjD3PGPvPjLF81jVxIdZjSlWFme87ZXNswdDr9TTnipAsSb/6lWx74uabYbFYFGrN4ugffli2bTl5Eld7e8E5V6hFhSGTntVWzvlXZ3uQc+7mnD/LOf8sYtl3HQD25KqBaboxvbdkA3BokccWDL1ej8nJSZoMSEimfD4Ihw/LdnnvvBOGlErmxURYvhz+TZtk+0xdXWU/FJhJsBrI9MXiges3nPOnF9GmuV7fBaCbMZbav+8A0AUAjDFH4rH5ji1kgiBQcVtCMvXGG2ApixUG6upQccMNCjZo8aJpQ4G1v/1t2Q8FZhKsCq3vuQPA9sREXwCPxQMTAOyCvFc367HxCcG7Ee8Jxv9eEEwmE0ZGRmjOFSEZSB8CHL/pJlSmTK4tRtrPf162bf3oI4yV+VCgWukGZCsebPbN8tieLI7tAdAz2+NKSi1ua7fblW4OIYWLc+CVV2S7nLfcgroiS1lPp1m3DuK6dTCePAkAYJzD8vrr8N5+e9GUj8q1THpWNzDGMrpTyRj7ZHyu1XOLbFfZ02q1uHTpktLNIKSwHT0K1chIclPS6aAq4FWBsxFJmyBc+9vflnVJtkyC1X4AnanzqVLFA9Qz8QA1EZ9rtT2XjSxHRqMRY2NjNOeKkDlIzz8v23a2t6N6xQqFWpNb2i98QbZt6e3F1aNHy3btu3mHATnnA4yxFwAMMsa6AHyAWAq4A7HSRRMAdnHO38hrS8tMYv2diYkJ1NXVKdwaQgoQ5+CdnbJdY5/4BBxFmrKeTr9uHbzr18N84gSA2FCgrasLE7feipqaGoVbt/QymhTMOT8A4BEALYjd49kT//sTnPPViUDFGFvFGPsLZJFBSGZnMBhozhUhszl+HOpz55KbUbUawXvvhU6nU7BRuRV+8EHZdt2bb2J4eFih1igr41ok8cKvHZxzIf7TwTl/Nu0wW3zicGuO21mWdDodzbkiZBbSc/Jb4672dlS1ltZXj/bRR2XbFb29CJw9W5bfCTktnBUvtURyiOZcETKz9CHAK7ffXnLZs8Z16zC5caNs37LDh8uy4PWcwYox9h3G2CeXqjFkOppzRcgMTp6E+vTp5CZXqeC+804YjUYFG5V7jDGE07IC619/HSMXLyISiSjUKmXM17PaD+CeeMXxH+R5+Q8yg9Q5V4SQmNAvfiHbdrW1oXL16qKrsp4J3Re/CJ7y7zKcOwfT8eNlV9Fizv9ZzvkA5/wJznkHgAMAvsoY+4Ax9m3GWPNSNJDQOleETDPDEGCpZsgZW1vhvOUW2b7GV1/FhQsXyqqiRTYJFh9yzr/KOb8esfp6+xhjr8UrrZdGrmiBonWuCJkSPXUK2tS1qwQBzjvuKNnKDowxRL78Zdm+qq4uiJcvw+VKX1+2dC2ozxwvVvtZzvm9AJwADsYD12dy2zwCxJIsOOdl1+0nZCaBn/1Mtu3esgVV11xTElUrZmN46CEEq6uT2yq/H41vvYWBgYGy6V0teoCXc/5CfF2rzwKoYoy9zhh7jhIzcstoNOLixYtl88EkZEbRKNRpwerKHXegOuWLvBSZbTaMffrTsn11L78Mt9sNp9OpUKuWVs7uRqasa3UPYisKt1NiRu7odDp4vd6ynF9BSELotdegTZkoH9VqMb51a9EutJgpxhj4n/4pOJtahN144gSqL1xAf39/WWQL5yV1Jp6Y8XRKYsYfU3HbxVOpVBgdHVW6GYQoJvLMM7Lt8TvuQGVra0kPASbYt2yBs6NDtq/+5ZcxOTmJiYmCXgA9J/Ke5xlPzHiCc/5Ivt+r1CXmXJVrIUtS3qJjY9C/+qps38U//EPU1tYq1KKlZTKZMP4ZeVqA/ZVXUBGJlEXvakHBihIplKFSqRCNRsviKoqQdL5nnoEQDie3A42NmGxrK/khwATGGEyf+xxCKVU6VD4fGl98EaIolnylm4VmA77IGHuMMfYNmm+1tHQ6HRW3JWUnKklQ/eQnsn2X77sPtXV1ZTEEmFBdX4/h7fIVmKr/3/+DBcDZs2cRCASUadgSWPAwYDyZ4u8RS6T4BvW2loZer4fL5YIoiko3hZAl4379dRj7+5PbXKXCyN13l80QYIJWq0XosccQTplTpp6cRP2LL4IxhrNnz5ZsxnAu7ll1IbYkyA2MsXPxtPX/TD2u/GCMUaIFKSvRaBThtMQK9223gdfVwWq1KtQq5SxrbcXwww/L9tX88z/Dqlbj6tWrJfvdsOBgxRj7TDzDbyuArngSRWs8keI3iPW4nomnrn+bqlzkjslkwvDwMMIp4/eElCrniROoTkusGP6DP0BjY2NJ1gKcj9Vqxegjj0AymZL71E4nqjo7YbFYcObMmZIcDlxogsV3AOwF8Bjn/EXOuazKajx1/YV4eaavAfgO59yTg/YSTCVaUEULUuqi0SjCf/d3ssSKUF0dxq+/vuyGABMEQUD9unW4mJYZWPPTn0ITiUClUuHUqVMllzW80MuSnQD2ZhqA0oMZWTyj0YihoaGST1cl5e3q8eOoefFF2b6LX/wibNXVMBgMCrVKeTU1Nbjw8MOQUs6BZnwcy374Q5jNZrjdbpw5c6akvh8W04cuzbt4RUKr1UIURXg81GElpSkQCCD07W9DFQwm94VqazG0bRsaGxsVbJnyjEYjzM3NuJzWu6r98Y+hP3sWVqsVly9fRl9fX8kkXCw0WD2BWC1AoiCtVovh4WGlm0FIznHOMfjBB6j/5S9l+0e/9CWojEbYbDaFWlY4Vq5cif4//mOEU5ZGYZEIGv/mb8CiUdjtdgwPD+P8+fMlEbAWOs/qAIAuxti3c9wekgWj0YirV6/C7/cr3RRCcurKlSswPPMMVCmJAuHqagzdcw+WL19eVnOrZmO1WmGsr8fAf/tvsv2m3l5UP/88GGOw2WwYGBjAuXPniv4e1mLmWT0N4ADNr1JOIo390qVLSjeFkJwJBoMYfPddNL70kmz/lS99CWG1umwTK9IxxtDc3IyRm26C+667ZI/Vff/70Fy6BEEQYLfbcfHiRRw9erSoswQXlfcZz/p7cf4jSb6YzWZcuHChqD+EhCRIkoRTJ0+i9X//b6hSJr6HKytx4VOfQlVVFUwpKdvlzm63w2AwYODxxyGZzcn9KlFE8+OPQ/D7wRiD3W5HIBBAd3c3xsbGijLxovwmKZQYQRCgUqno3hUpetFoFKdOnYLml79E1VtvyR678pWvwM8YVq5cqVDrCpMgCFi5ciVcJhMu/df/KnvMeOIEVu7eDUQiAGLzM3U6HY4fP47u7m5MTEwUVdBSK90AsniJScINDQ0wGo1KN4eQrEWjUZw5cwbuc+dww/e+J3vMt2kThh54AJUVFWVTtDYbNTU16Ovrw+UHHoDl8GFYfve75GOWt97C8u98Bxf/8i8BxqDVaqHVahEIBPDxxx9Dr9ejpqYGdrsdJpMJKpUKjDEwxiBJEoLBIEKhEILBIERRTP6kBjm9Xo/q6mpUVFTAZDJBo9Hk5d9JwaoECIIAtVqNoaEhXHPNNUo3h5CsRKPR2Jft5cvY/I//CLXLNfWYRoML3/wmgpEI1jc3K9fIAqZSqbBq1SqcOXMG5/fuRctXvgLjqVPJx6sPHoRkseDy178OxCt+6PV66PV6hMNhXLp0acaRGc45WMpijyqVChqNBmq1WrY/FAphYGAAnHMIgoCmpiY0NDTkPGhRsCoRZrMZly9fRmNjI8wpY9eEFDJRFHHq1Cl4PB443ngD9tdflz0++tWvwllfD5vBQL2qOdTV1eHixYvwSRIGvv99rP6TP4E2JfFq2Y9+BH1fH4aeegrRlCK4Go1m0UFFpVJBr9cDiN1zPH/+PIaGhtDc3IyGhoacZW7SPasSwRiDRqPB4OBgScypIKUtGo3i0qVL6O7uRiAQQPMHH2DFt74lO0Zcvx5XvvQlBINBNDc3y67miZwgCFizZg38fj/C1dXo/4d/QCQlKAGA9fBhrHn0UejPns1bO1QqVSyl3mhEX18fjh49mrOpNRSsSojJZMLY2FjJVl0mxU+SJFy5cgXd3d04ffo0TCYT6np60PTf/ztYyn2QqFaLC9/8JsRQCFartSyrq2fLarWioaEBk5OTCDocGPj+9xFJWagRAHRDQ1jzuc9h+d/9HdRXr+atLWq1Gna7HX6/H/zmuqEAACAASURBVEeOHMlJHVMKViWEMZasuuzz+ZRuDiEAYvc+fD4fLly4gPfeew8nTpwAEEu7tr/7Lpq/8Q0I8Yw1ILZW1fnvfAe+lhYEg0GsXr2aelUZSvRAI5EIxC1bcObnP4e4fr3sGBaJoPr553HN/fej7nvfgyaP8zTNZjN0Oh2OHj0KlUqlXcxrUbAqMWq1GhqNBidOnEAk5QuAkKUkSRLcbjfOnz+P9957D93d3RgYGIBWq4XdboeeczTs3QvH178OIaX2HwAM/c3fwPPJT8Lj8aClpYXuwWZBp9OhpaUFHo8HnHOE6+tx7sc/xvhDD007VhUIYNmPfoR1n/oUHF/7GqyvvQaWh/maWq02kZSxqJtXlGBRgoxGI9xuN/r7++mqlCyZcDgMt9uNK1euYHx8HNFoFIIgwGAwyCqkG48exYq//mvoU1b+TRh+8km47rsPoijCbDajoaFhKf8JJWHZsmVwuVwYGxuD1WoF1+kw/Fd/hcmbb0b9974HXVrmH+McFb//PSp+/3tIRiM8d94J1z33YPKWW8C1i+oM5RQFqxJlsVhw8eJFqFQqNDc3Uy01kheRSAQTExO4fPkynE4ngNiVtNlsli+MGI3C8uabqPnZz2Du6Zn2OlwQMPLnf47xRx5Jzu/ZtGkTfW4XIDXZYnJyEhXxRAv3PffAc9ddqHzhBSw7cACaiYlpz1WJIuyvvAL7K68gYrHAde+9cN5/P8RrrwUUvuhllDmWnbVr1/Lnn39e6WZkJBqNwuPxwGw2Y926dTRhmOSMz+fD6OgoLl68CEmSkvN2ZL34aBTGjz+G9Y03YO3qgm5kZMbXCi5fjqGnnoJ43XXgnMPpdMLhcFC1ikUKBoP48MMPAWDa2l+CKML+8suofOklGOP3EOcSWLkSVz//eUw88AD4AtYRc7vd2Lp164lwOLwh6yfHUbDKUjEFqwSfz4dwOIzW1lZUVlaW9aJ1ZHG8Xi/Onz+PsbExqFSqZNWDBMHjQcW778Ly9tuoePttaObJOJv4oz/Cxb/4C0Tj96VcLhdqa2uxdu3aslyyPte8Xi8+/PBDaDSaWX/v9adPo/Kll2A7dGje/6+I1YrxHTsw9uijkNIyDedCwUoBxRisgNhwjc/nA+ccZrMZdXV1MJvNMBgM0BbQuDQpTF6vF4ODgxgbG4NWq4XJZEr2otRjY8nek7mnByyDpSgmb7wRV778ZXhvuim5z+12w263Y/369TT8l0NerxfHjh1DOByee2K1JMH04YewvfYabIcOySqJTDvUbMbozp24+rnPgWcwqZiClQKKNVilCgaDCAQCycnDGo0mlkZst6OiogIGg4GuagmA2GdlaGgIw8PDsiDFQiFYu7pQ9cILMPX0gGXwPRLVauG6+26MffGLCKSVBXO73bBYLNi4cSPUarqVnmuhUAinT5/G+Pg4rFbrvL/fLBxGxdtvw/7yy7D89reyqQWpgitWYOTxx+G58845X4+ClQJKIVilSy1YCcTS3xsaGlBdXQ2z2UzZhGVIkiRcunQJAwMDAJBMmNCMjqL65z9H5a9+BXU8oWLO1zEa4bn9drjvuguTt92GaNryHokU98rKSqxfvz5vRVBJ7B72hQsXMDg4CLVaLesdz0V95Qqq/+VfUNXZCfXk5IzHOD/1KVx84glIs/TcKFgpoBSDVbpIJAJRFCFJEgwGA5qamlBTU0NXvGXC6XTi7Nmz8Pv9qKioiBUwvXwZtT/8ISpfeglCODzn8wMOBzy33YbJ226D77rrZh0mEkURwWAQra2taGhooN78EhFFEYODgxgdHYVer8848UoQRVQ9/zyW/dM/QeX1Tns8VFuLC9/8Jry33DLtMQpWCiiHYJUqFApBFEWoVCqsWLECdXV10Ol0SjeL5IHf78fAwABGR0dhNBqh1+uhcjpR98wzqHzhhVmHgoBYHT/3tm1wbd2K0DxZfIFAAH6/H2azGddccw1N+lWI2+1GX18f3G538v87E6qJCdT94AeoeuEFWYmshCtf+hIu/Zf/AqRc3FKwUkC5BasESZLg9XrBOceKFSvQ0NCQ8YebFLZwOIyRkREMDg5CpVLFhvwiEVT/4hdYduDAjFfRABCpqIDzgQcw/vDDCDocc75HNBqFKIoIh8OoqKhAU1MTqqqqqDelsMRUgb6+Pvh8PpjN5oyHYvWnTqHpf/5PGGYojDt5ww04/53vQKqsBEDBShHlGqwSJEmCz+dDNBpFfX09Ghoa6Mq4SEmShLGxMfT39yeDiEqlQsWbb2L5009Dd+HCjM8LNjTgyp/+KZz33w8+xwVLJBKB3++HJEkQBAHV1dVoaGiAxWKh+6AFJhqNYmxsDGfjgSfTe9UsFMKyZ55B7U9+Mq2XFaqrw+Df/z38GzdSsFJCuQerhGg0Cq/XC0mSYLVa0dTUBIvFkvFVGec8mZXo9/vh9XqT88Gi0SgikUiyzqFOp4PRaERFRQX0ej0MBkNBpzZLkoRAIIBgMAi/349AIIBwOIxwOAxJkqBSqZILZiZKEWm1Wuh0Ouh0urx/kSeC1ODgIILBYHJ1V+3wMBqefhrWw4dnfF6woQFXHnsME/ffD8zw/5z6fwrEskxra2tRVVUFi8VS0P9nJCYYDCYXwszm99n44Ydo3r0bmrEx2f6oToehb30LQ9dfT8FqqVGwms7v9yMYL0ZqsVhQW1sLo9EIlUqVXCY7HA4jFAohEAjA5XLB4/FAkqTkaqQqlQpqtVq2rDbnHJIkIRqNJoNY4ovcbDYnl9LOZugiH4LBIHw+H9xuNyYmJuD1epPtB2LZlYIgQBCE5H7OOaLRaPLfl9ivUqlQUVEBm80Gi8UCo9GYs3uEfr8fV69exfDwcDJIabVaCH4/an78Y9T+5CcQ4hmhqSSzGaOPPRabUzPDnLxAIJCcClFRUYH6+vrkmkbUgyo+nHNcvXoVJ0+ehE6ny3i4Xz0+jpW7d8N85Mi0x/p37cLaH/6QgtVSomA1u8SVdTAYlC2Jnfh74k+tVgutVrvgK23OeTLwJT6/JpMJVVVVsFqtMJlMeU0CCYVC8Hq9cLvdGBsbSy4uJwgC9Ho9NBrNgr+ko9EoQqFQ8hwCsVp7Npst+W8zGAwZvYckSRBFEZOTk7hy5QrcbjcYYzCZTLHMzmgUtldfRf13vwvtlSvTns8FAeMPP4zLX/ta8t5D+mtHIhFUVFRg+fLlsNlsdB+zhExOTqK3tzdZSCAj4TAavvtd1PzzP097SKNWU7BaShSsCg/nHOFwGIFAANH4uHnqF3wi00mr1WZ9Qz8xpOf3++F2uzE+Pg5RFMEYS1YUn7dXxzkEvx8qjwdgDFGNBlynQ1SvBzII2JFIBKFQCOFwOBnAEkHHaDQme6SCICTPQ6LNqRcIydp9nKPid79D3Q9+MGtdON/mzRh+8slpk3cTiTYA0NDQkKyEQkpTIBDAsWPH4Pf7565+kcb+0ktY8dRTYCkZpBSslthigxULhSB4vVB5vVD5fGDBIFg4HJu7EokAKhU4Y4BajajBAMloRNRsRsRqnfNmNpFLfMGHUoa1GGPJe14GgyF5T0ytVsuGHIPBIERRlA1vAbEluxM9p5kIHg9Mx45Bf/YsdAMD0Pf3QzsyApXbPePcJK5SIVxbi1B9PUINDQi0tMC/di3811wzrSeTLnFfLxKJJIcUOecQBEE2pJpe+bzi7bex7MABmHp7Z3zdcFUVLv3Zn8F5//2yKtuSJGFychKCIKC5uRl1dXU0gbdMhMNhHD9+HJOTk1kFLPP772Pl448nJxKXXbBijNkA7ATQD8ABoItzPn3NgXmOzeZ1UqUHKxYMQj0xAfX4ONQTE9CMj8f+7nTG9k9MQO10QuXxQOVyQbWIxc0ksxnhmhqEa2sRXLECoaYmBJuaEGhpQWj5cqDI0oAFrxf6/n7o+/uhGxiAZmwsec4EUZw6kDFErFZEKish2e0I1dUhuGoVAs3NCK5cmXEVaM558gs+EZgSPbFEDyT9/ll80bgZX0999SrM778Pc3c3jEePQt/fn1HZoUyE6uogbtoE37XXQrz2WvjXrl3wxYrm0iVU/uu/ovJf/xXaWVaFjWo0uPqFL2D0K19JFpUFYudlcnISnHOsWrWKglSZCofD6O3thSiKySVHMqEbGMCqr38duosXyzJYHQKwi3Pen7K9g3M+reriXMdm8zqptpjN/O2WllgwcjpnnYOy1CSjEYHWVvivuQbixo0QN21CcOXKwglg0Sj0Z8/C9NFHMPb2wnT06Kyp0dngjCHY3Az/unUQ162Df8MG+NetQzQPleVVbjdMR47A/MEHML//Pgx9fTl/j9lwtRqBlhaIGzYg4HAg2NSEUFMTwlVVseFEjQaQJKidTmiuXoXu/HmYjhyB6ciRedvp2roVl/7szxBasUK23+fzIRQKYfny5WhqaqLJ4GUuFArh448/RjAYzGroVzUxgcY//3PUHD9ePsEq3hs6wjlvSdm3H8AhzvnBTI8F0JXp66TrYIx35+Rfk3+S2Zy8Mvdt3gxx40ZEs7gqWtybSzCcPQtTdzfM8S9NtcezJG/NBQEBhwP+9evhX7MG/rVrEWhthWSzZbyAHAsGoTt/Hsbjx2E8dgzG3l7oz55dcM8pqtVCslpjrx0OgwWDUMUTM3KBq1QA5zNWFJiN+847MfrVr8Kfdl8qkUBSWVlJy8oTmVAohKNHjyIcDsOUVudxLm6XC1u3bVtUsCq2Ym8dANJ7Pi4AdwNIDzJzHevK4nVyiqvVkMxmSCYToiYTono9uEYT+0n5wmGRCARRhCCKUHm9UDudGS29kErl9cLyzjuwvPNO7L0ZQ8DhgHjttRA3boT/mmsQaG0Fz8EVsyCKMMS/2E0ffQRTT8+sRS/zjUWjMJw7B8O5c7L9ksmEUEMDQvX1kCwWRA0GRA0GsGgUgs8HQRShnpiA7sIFaEZHsw5MnDEEWlshbtwY6/04HAg2NyNcVTXjEB7z+6G9fBnaS5egO38e+jNnYDh1Cvpz5+atvzfttTL8bET1eri2bsXVRx+Ff/16+WPRKCYnJ6FWq7Fx40ZUV1dT6jmR0Wq12LRpE3p6euD3+zNfGy8Hn6NiC1Y2AOlrMY8jds8pm2OzeR0wxnYidn8L7WmPcZUKkcpKROx2RKqqEK6qim1XVcX2VVYiYrNBstsRsVpjVacX8h8XjULtckE9NgbtyAh0Q0Oxn8FB6M+cySgwMM5h6OuDoa8PVb/8ZbL9gVWrEGxuRmjFCgRXrEC4pgaSzYaIzSYbTmPhcHL4Uz0+Dt3587H3HxiAbnAwq6v6xHsHV65EwOFAwOFAqKEhdt4qKyFVVCTPE4tEoHK5kkNc2gsXoB8chG5wENrh4YyDisrng+Hs2RnLwywEFwSI69fDe8MN8HV0wLdpU1Y9V24wILhqFYKrVmEypfgnC4VgOH0axqNHYfr4YxiOH4fu4sWFt5MxiBs3wvnpT8P5B3+AaNpNcs55ckL2ypUr0djYSPelyKz0ej2uvfZa9PT0QKVSLdl6eMUWrABg7jSpzI/N+HU45wcAHACADcuX83NPPRULRHZ7rCT+UtwXEoRY4KusRGDt2vQGQnP5MgynT8Nw4gRMvb0wHDuWWQCTpBl7IfkQqaiAr60Nvi1bYr279esXtER2KsHng+HUKRhOnIDxxAkYjx+HbmgoRy2ezr9mDbzXXw9vRwe87e3TvvhzgWu1EDdtgrhpExLrtqqcThhPnoT+9OnkhYp2eBgqrxeC35+8UIhYrbGLppoa+NeuhbejA74tW2Ztp9/vh9/vR21tLVatWpVxBW5S3sxmM6699lp89NFHyWos+VZswcqFWK8oVRWm95LmOzab15GJmM3wtbVl1NglwxjC9fUI19dPLYIWjcbuuRw9ClP8RzcwkLNstUxELBb42trgbW+H9/rrEVi9OqN5RdmImkzwtbfD1z7V51V5PDCcOAHDmTOxobUzZ6AbGoKQRSYmZwzhZcsQWL06mbAibtiQvO+U0WvEsw8TmYep1SqSq+yq1TOnmaeR7HZM3nKLrAeW8kbJ+SyZrNoKxObPiKIIq9WKtWvXwp7FEuWEAIDNZsO6detw4sQJWK3WvJfTKrZg1Y3pPSIbYkkT2RybzesUJ0FIDjE5H3wwtmtyMpYw0NsLw8mTMJw6Bd3ISM7eMrBqFcQNG2Ip19ddh0BrqyLZiJLFAu9NN8mWTAfnsWHEixehHR2N3aPy+yH4/YBKFZvPZjRCqqhAqLERoYaGrO7lzTQxGUByXpdOp4NWq4VarU6mzEej0WTQ8Hq9yecxxqBWq6HT6TK7YmUsoyCVXvl8y5YtsNlsdF+KLNiyZcsQCoVw7tw52O32vH6WiipYxVPOuxljjkTKOWKJFHsAgDHmiB/XP9ex871OqYpWVEz7Eld5PND190N34UJyaEntdELldkPtcoGl1ooTBERstuQQaLihITbXyeFAoLk5b5mGiVp6ib/HmiJk94vBWHIY1b9pU07alF7yyWw2o7GxEVarFXq9HjqdLqurzVAoBL/fD5/PB4/HA6fTmawWwRhLBrxsqnAkVoEOBoNgjKG2tpYqn5OcamxsRDAYxPDwcF4vfooqWMXtALCTMZaYzPtYytyoXYj1kHZlcOxcj5UNyWKBuGULxC1blGuDJMmqkgOQFYIFkCxwmxBJKeOSOFYQBGi1Wmg0mrwMSYTD4eSSF4wxmM1mrFy5Ejabbare3iIkaiZarVY0NDQAQLKahtfrhdPpTBYAZowhGo0mg3biJ1HJIrXqht1uR1VVFSorK5fsZjgpH4wxOBwOBINBjI+Pw5rFUHlW71NM86wKAdUGXJz0YrcAkgv+VVRUJKuMJ4bMUiuxp7+OJEnJskqJyueTk5Nwu92yoKfRaJJllTK96otEIslK8Yl2GgwGVFdXw263K1bpPVGBPnEOEzUDE/fEEjUAE0urGAwGWuCQLAlJktDb2ztjWaZcrGdVjD0rUoQSdfYAoKKiAg0NDcngtJAq5Yn7Omq1Olnpu6amBsBUQPT7/RBFER6PBx6PB6IoygrBpl+opS7fodfrYTabYbVak+0shF6JIAjJda8IKSQqlQobNmxAb28vvF5vzieTU7AieZM6f6eiogJr165FVVVV3nskiYK1er0edrsdy5cvBzDVKwmFQskCsImf1HqAiV4dISQ7Go0GGzduxNGjR3MesOg3kuRFYomKqqoqrFq1qiBK9lCvhJD8S1S5OHr0KHw+X1ZlmeZCg9kkpyRJgtPphCAI2Lx5MzZu3FgQgYoQsnT0ej02b94MtVqNyRyVXaNgRXImEAjA7XajpaUF7e3teZ93QQgpXImApdfrIaYu+bNAFKzIonHO4Xa7AQDt7e1YsWIFZaARQqDT6bB582bU1dUBwKJSz+meFVmUcDiMyclJLF++HA6HgxITCCEyGo0GmzZtQiQSCS7mdeibhSyYz+dDJBLB+vXrsWzZMqWbQwgpUPFJ+tSzIktLkiRMTk7CbDZjy5YtVKmbEJJ3FKxIVhLzphwOBxoaGvJeaZkQQgAKViRDiergNpsNmzdvztncCUIIyQQFqxInSRICgQDCKcukM8YgCAI0Gs2cVbwjkUjyuRaLBddeey3sdjtl+hFClhwFqxIVjUbhdruhUqlQXV2Nmpoa6HS6ZHVzURThdrvh9XqTVbw557L6eDqdDlVVVaivr6clJQghiqJgVYIkSYLL5YLD4UBTU9OcPaFoNJqslZeo3C0IAgwGQ0EUbiWEEICCVcmJRCJwu91Ys2YNGhsb5z1eEIRk1XJCCClUFKxKiCRJ8Hg82LBhA817IoSUFLpTXkLcbjdWr15NgYoQUnIoWJUIURRly6ETQkgpoWBVAiRJQjAYxJo1ayitnBBSkuibrQRMTk6iubmZ1o0ihJQsClZFLhAIwGAwYMWKFUo3hRBC8oaCVZETRRGtra1Uo48QUtIoWBWxQCCAiooK2O12pZtCCCF5RcGqiImiCIfDQWWQCCElj4JVkQoEAjCZTLDZbEo3hRBC8o6CVZHy+/1YtWoVpaoTQsoCfdMVoWAwCKPRiKqqKqWbQgghS4KCVRESRRHNzc3UqyKElA36tisykUgEGo2GelWEkLJCwarI+Hw+NDU10bwqQkhZoWBVRDjniEajqKmpUbophBCypChYFRGfz4dly5bRYomEkLJDwaqIhMNhWgKEEFKWKFgViUS6usViUbophBCy5ChYFQm/34+mpiYqrUQIKUsUrIqAJEkQBIHS1QkhZYuCVREQRRENDQ3QaDRKN4UQQhRBwaoIRCIRLFu2TOlmEEKIYihYFbhgMAiz2QyTyaR0UwghRDEUrAqc3+9HY2MjJVYQQsoaBasCFo1GwRijxApCSNmjYFXARFFEbW0tJVYQQsoeBasCFg6HUV9fr3QzCCFEcRSsClQ4HIZer0dFRYXSTSGEEMVRsCpQPp8PjY2NtMAiIYSAglVB4pyDc47q6mqlm0IIIQWBglUB8vv9qKyspKVACCEkjoJVAQoGg1i+fLnSzSCEkIJBwarASJIEtVoNm82mdFMIIaRgULAqMKIoor6+HiqVSummEEJIwaBgVWAkSUJtba3SzSCEkIJCwaqAhEIhGAwGmM1mpZtCCCEFRa10A7LBGLMB2AmgH4ADQBfnvGchxzLG2gA8CWA/57wr323PhCiKWL16NRWtJYSQNEUVrAB0AtjFOe8HAMbYIcbYDs65K5tjGWPb4sc4lqbZ8+OcAwDNrSKEkBkUzTBgvKfkSASfuH4A27I9lnPeFe9NTeSxyVkRRRE1NTXQarVKN4UQQgpO0QQrAB0A0ntQLgB3L/LYghAKhahoLSGEzKKYhgFtmN4TGsfMQ3nZHDsvxthOxO5/gTEWufHGGwcX8jpzvwVTBQKB8Ry/br5VA7iqdCMKBJ2LKXQuptC5mLJ2MU8upmAFAJV5OnZOnPMDAA4AAGOsOxAIdOTqtYsZY6ybc07nAnQuUtG5mELnYgpjrHsxz1c0WMV7LO3zHLY3fu/JhViPKVUVZr7vlM2xhBBCCpyiwSreY8lUN6b3lmwADi3yWEIIIQWuaBIs4unp3Yyx1PtOHQC6AIAx5kg8Nt+xi5RNgC11dC6m0LmYQudiCp2LKYs6Fywxv6cYzDXRlzG2F4CNc74rg2PbEEtjfxKxXtghzvm+HLd1f6It5Sh+/rch1sO9G8CetKkEJY0xth3x6RK5/mwVk3L/HMym3L8fgFgHA7HPxgSAnvk+F0UVrIpFfNLxXs75fPfjSlb8fqSNc74vcW+yXH454///Ds75gcSFUbkGrHL+HMyGvh9iGGOdnPMd8Qs7cM4PznV80QwD5gNjrI0x1plS0SKx38YY280Y2x7/sy2L17QhdkVdVMkcuT4XnPMDKV/QLQD6ct3mpbKAc3M3Yp8BoMDn92Ur23NRSp+DdAv5nSnW74f5ZHsu4gHqA8aYg3N+cL5ABRRf6nrOzFNyKZuyTuk6OOddxVTfL4/nIsHBOd+z2HYqYSHnBvJM1AnkcBqFknLwOSnaz0G6RZyLovt+mM8Cf0euR3w+LGNsN+ao85pQtsEqUbyWMSa7wpmnVNPB+FDGTK93gDG2rVCK4mYjH+ci5TV2c8535L7VS2OB5yZ16kQlSuQqeqGfk/gxRf05SLeQc8EYcxXj98N8Fvi5GAcwHq/V2oNY/sCcn4+yDVZzmKtU08F50u0nEuOvABzFGrxSLOZcJLr6icnUxX4u0s11bjoxdZXpQOlPmZjzc1Lin4N0c52L/SX2/TCf+X5HUm8pzHtBV9b3rGYxW6mmeYdyOOc9KWOvpTD0s+BzER+b3gvgN4yxPhRQhfscmfXcxL+AbPHhkbYySK6Y9VyUwecg3Vyfi1L7fpjPfL8jiQvaNgDzDg9Tz2pmi/ogxT+Q894wLBILOhfx8eeWHLel0Mx6blICVClfOaea8VyUyecg3Zy/MyX2/TCfTH5HMkI9q+moVNMUOhezo3Mzhc7FFDoXU3J6LihYTUelmqbQuZgdnZspdC6m0LmYktNzQcEqTZ5LNRUVOhezo3Mzhc7FFDoXU3J9Lsq2ggWbo+QSm6NUUymiczE7OjdT6FxMoXMxZanORdkGK0IIIcWDhgEJIYQUPApWhBBCCh4FK0IIIQWPghUhhJCCR8GKEEJIwaNgRQghpOBRsCKkSDHGHIyxvUq3g5ClQMGKkAKRCD6MsZ0pS0nMZRdSStfEn3uEMcYZY/tTV22Nv+ah+GOds61FRkihoknBhBQIxtgRxBag2wVgG+e8fb7j04+JB6G9nHP7DMe3ATgCwJ7lSs+EKI6WCCGkAMQDiYNz3h9f92nOYp/x47uXpHGEFAAKVoQUhkcQL/A53wrMcbsA7M9riwgpIHTPipDCsA3ZLZ3QUcrFUQlJRz0rQhTEGNuN2Eq6bQDuZoy1A9g/VyCKJ07kZMmJ+HDibwB8G7HK2ECsOvZe0L0tUkAoWBGiIM75vnjA2Mk535Hh03YB2DPH47Z4EEw30/LylQAeiy+1DgBgjB0CsIcCFSkkFKwIUV4Hpno1mbBxzuc63pVYTyhVIiimvxZSemnxbMLKmZ5PiJIoWBGivHYAGd1/is+/6szhe3clelDxFV33xttDSEGhBAtClNcB4IMMj90F4PlcvXHaUF8nYsN/2fTyCFkSFKwIUV4bMkiYiC8Rnh5gciJxjys1bT4+bEhIQaBhQEIUlAgIGaahfxZ5mFsVH/57EinDf/F9lbl+L0IWinpWhCgrm+SKHalZezk00/DfdgATeXgvQhaEelaEKCuj5Ip4T2fO4b94BfZtiKWu7wfQyTnvij+2E7G6gwDwLGPsOc75wfh+B4CJePJGZbxNOzFzqjshiqBCtoQoKF689tvz9Zji95R6EsGHkHJDwYqQJRbvwbg4512MMc45Zxk8Z1qFdULKCd2zImTpPQugLV42ad7Jt1RhnRC6Z0WIEhKlU+Gu9QAAAFRJREFUku7mnM9VNinhEVCFdVLmaBiQkALHGOvMom4gISWJghUhhJCCR/esCCGEFDwKVoQQQgoeBStCCCEFj4IVIYSQgkfBihBCSMGjYEUIIaTg/X+wZ5ra0t+VPgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# plot the DRT and its confidence region\n", + "plt.semilogx(freq_vec_star, gamma_vec_star, linewidth=4, color=\"red\", label=\"GP-DRT\")\n", + "plt.fill_between(freq_vec_star, gamma_vec_star-3*np.sqrt(abs(Sigma_gamma_vec_star)), gamma_vec_star+3*np.sqrt(abs(Sigma_gamma_vec_star)), color=\"0.4\", alpha=0.3)\n", + "plt.rc('text', usetex=True)\n", + "plt.rc('font', family='serif', size=15)\n", + "plt.rc('xtick', labelsize=15)\n", + "plt.rc('ytick', labelsize=15)\n", + "plt.axis([1E-4,1E6,-0.01,0.025])\n", + "plt.yticks(np.arange(-0.01, 0.025, 0.01))\n", + "plt.legend(frameon=False, fontsize = 15)\n", + "plt.xlabel(r'$f/{\\rm Hz}$', fontsize = 20)\n", + "plt.ylabel(r'$\\gamma/\\Omega$', fontsize = 20)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4e) Plot the imaginary part of the GP-DRT impedance together with the experimental one" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEZCAYAAAApEwoTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df3xbV33/8ddxrKRRbNVN6QMY343WLixdaUdtd2P8+gJ1YCl8xzbkFIizbqy1m4QCJTQmjO+gG5A6Jcv6JSW1G36tYaO1CttgycCmY8AeYzRxWaBsjMjtgA1Gm0SJEyet4pzvH/de5epKliVLsq6s9/PxuI9Y9x5J1ze2Pj73fM7nGGstIiIiYdZU6xMQERGZi4KViIiEnoKViIiEnoKViIiEnoKViIiEnoKViIiEnoKViIiEXnOtT6BUxpg2oB+YBNqBcWvtRKltjTGd7j6Aa4Ej1trtVT59ERGZB1Nvk4KNMWPAgLV20ve411qbKqWtMeYYcJ0veFmga7bAJyIitVNXtwHdnlK7F3xck0DPPNr6A1Wbuy8n4ImISO3VVbACuskNKClgdaltAz2otUAiENhERCQk6m3Mqg04Gth3hPNjTyW1Nca04/S0Vltre2d7U2NMP87YFytWrOhatWpV6WcuItLADh48+JS19pL5Pr/eghXAykq1dXtSI8YYjDGjswUsa+0IMALQ3d1tDxw4UMIpiIiIMeY/y3l+vd0GTOH0mPwuJrcHVVJbNxh1GmO2VOIkRUSksuotWB0gt7fUBoyV0tYY02mMSQaOTQIdFTlLERGpqLoKVm56+gF3rMnTDYyDMwblHZujbQpIBF6+m/xBT0REaqwex6x6gX5jjDfR92bfHKsBnN7TwBxtU8aYMTdxAqALGLTWBgOYiIiEQN1NCq41JViIiJTOGHPQWts93+fX1W1AERFpTApWIiISegpWIiISegpWIiISegpWIiISegpWIiISegpWIiISegpWIZRMJtm4cSOxWIympiZisRgbN24kmQxWiKpf4+PjdHV1MTg4WOtTEZE6oGAVMvv37+fqq69mz549TE1NYa1lamqKPXv2cPXVV7N///5an+K8jIyMZD3u6elhYGBgltYiItkUrEIkmUwSj8eZnp4mnU5nHUun00xPTxOPx+uuh5VKpUiltAiziMyfglWI7NixIydIBaXTaXbu3LlAZ1QZutUnIuVSsAqRvXv3FhWs7r///qqdw8jICOPj4yQSicxtukQiQUdHB6tXr870kHp7e+nq6mJiYiIz/rR9+3YSiUTmuZOTk4yPjzM5OcnY2FjmtYOCzylkcnIy8z6Dg4OZ80kkEnR1ddHR0UEqlWJychJjTNZ5dHR0MDAwwMjISNHvJyIhYa3VVsLW1dVlq8UYY4E5t6ampqq8/9DQkD148GDWY8/o6KiNx+OZx2NjY/bYsWOZx8PDw7azszPr8ZYtWzKv43+tYp4zm/b29szXyWTS9vT0ZB4fO3bMtre322PHjtljx47Z4eHhnO/P/z147UWk+oADtozPXvWsQqSlpaWi7UrV3t7OzTffzMjICKlUiv7+/syxeDzO+Ph4pieTSqVoa2vLeb5n5cqVRY1TlfKckZEROjs7s57rr4Df1tbG0NAQvb29PPjgg1nnn+/92traaG9vz0n+EJHwUbAKkb6+PiKRSME2kUiE9evXV+X94/E4W7duZXR0lIsuuihnrGnt2rWZQOb/0PesXJm9MPPRo0fnfM9SnpNMJkmlUoyPj2e20dHRnO+hFO3t7XWXsCLSiBSsQmTz5s1FBavbbrutKu8/Pj5OPB5nbGyMY8eOceDAgawxncHBQYaHhxkfH8/q4RTjyJEjgDO2NF/XXnst4KS9+ze/iYkJBgcHGRoaKmo8anJyko6Ojnmfk4gsDAWrEOno6CCRSBCNRnOCViQSIRqNZpIdqmFsbCzzAd/W1pYTCNrb22lrayuqxxR8XiVS1+PxOEePHs16LX/wS6VSHDhwgJ6eHoaHh+nt7c15DX8A8xIx8t0uFJFwUbAKmTVr1nDo0CH6+/uzKlj09/dz6NAh1qxZU7X37ujoyNxeSyQSXHvttTm3+wYGBli7dm3WvomJCUZHRzPPm5iYYHh4mImJCRKJRCbIjIyMZMa55nrObEZHR9m2bVsmg9Dr4W3fvp3LLrssc0tv5cqVTExM0Nvby8TERNZreO+5bds2xsbGyr5uIlJ9Wta+RI2+rL0XfOrR9u3bOXLkCENDQ7U+FZGGo2XtpeoGBgYy86VKHasSEakEBSuZU29vL5OTk0xMTOTNAqwH4+PjPPDAAyQSibwTk0Uk3HQbsESNfhtQRGQ+dBtQREQWPQUrEREJPQUrEREJPQUrEWlYjbAq92KhYCUiDWmxrsq9WClYiUhD8PeijDFcf/31i25V7sWsudYnIOHkVTTv6OjIlEjq7+9n+/btbNmyhcnJSYaGhhgZGSEej2cWZkwmk3R0dLBly5ZZX3tycpLh4WG2b9+e9Vyv2K2/woS/bX9/f6YuYjKZZPXq1ZlqGl5ppXg8zsUXX8wjjzzC+Pg4W7duBciUVVJ5pca0f/9+4vE46XQ6JzgtAS4HpoH/Bmbc/dPT01xxxRXcdNNNbN68WQWPa62cxbAacavm4othsWXLlqxFCj3xeDxrsUNrrQWyFmz02vX398/5PvmeOzo6mrUgY6G2nZ2dmQUWR0dH7djYWNbrtLW1ZbUv5pxk8Tl8+LCNRqNZC5g+H2w/2IfApsBad3sG7CTYh8FuB/tLkFkU1RhjV6xYYa+88kq7YsUKa4yxra2tdsOGDfbw4cO1/jZDDy2+KJXkFXkNrhMFZHopc/GWjp+PeDxOd3c327dvn7NtT09PVi8sWCU+uFZWV1fXvM5J6pN32++KK65genoagBcCDwNPAMPA7wIX+p4TAS4DXg3cDvwHcLe1PAfnD/tTp07x2GOPcerUKY1xLTAFK8kyODjIwMBA3mOdnZ0LUm6pt7c3Z+HHfFKpVKZWYTHn1d0978nzUmf8yRPpdJom4N3Ad3ACUbGWAbcCSWAIaM3TRmNcC0PBqpqMCcdWgomJiYLFaoupWD46OlrWGlFeUJlt8cRUKsXIyAhHjx7lvvvuAyiqwK6K8DaGZDJJPB7PJE+8EPgGsANYPstzngJ+XuA1o8AW4F+A58/SJp1Os3PnzvmetsxBCRaS4QWH4O0zPy/Zws+ryH706FGSySRtbW1lLcPhvcfk5GRWj8n/PqOjowwPD+c9H2lsO3bsyCRRXA+M4gSboH8Gvgh8GXgUZzArihOMXgb8EXBp4DlXuM97vfscv3Q6zf3338+uXbsq8n1INgUryfACQ3Al4MnJyaxK5e3t7VnjQz09PRXttXgrAQeDpv99uru76erq4tixYxV7X6lvyWSSHTt2sHv3bgBeC3we51ae31HgHcBn87zGNPBv7vYZ4A+B9wPP87V5LvB1II4T6PxOnjxZ3jchs9JtQMnS2dmZs7Jue3s7/f39jI2NMTw8nJPIUEhvby9dXV2ZrRheD69QAOzs7CSVSuWcqzQm/xgVOONSf01uoPpb4EryB6qgNHAv8ALgc4FjLcCXgBuD+1taSjtxKZqCVTVlEmJrvJVgaGiI4eHhvMfa29sL3iLMZ3R0lIMHD2a2YjzwwAMF52n5zTauJY0jOEb1cpzbe8HxqVuBNwI/8+2LRCJEo1E+8YlPEI3mu1kIp4G3AncF9jcDnwDW+NuePq1yTVWiYCVZenp6iMfj9Pb25hybLTAEbxuWY2JiImsyb6H3aW9v55FHHgHIu6BiJc9Lwss/RnUtsA9YEWhzKxAcSVq6dCn9/f0cOnSIt73tbSQSCaLRKJFIJOc9LE6CxTuAc779S3B6Xb/iPk6n00plrxKNWUmOoaEhxsfHGRgYyFSwSKVS3HfffXgLT3qVJYDMv8XeHgw+t6urK6uChb8H5rWNx+OZuV/e+4yOjjI4OEgikciMt01OTpJIJHjggQdIpVIMDg7S0dFRVnaihNvevXtJp9M8C/gCuenl7yY3UEWjUQ4dOpRVlWLNmjUcOnSInTt3cv/992fqBRpjsO4dio/hVLl4ACdQAcRwbi/+OnAEMlUy4vF4znvI/Gml4BJppWCRcGlqagJr+Tuyb8kBvBdnfpQnEokQiURIJBKsWRNsnSuZTGaC18mTJ4lGo1x66aWs/v73+bNz57Lafg0nqcMr5hSJROjv71d2oEsrBYtIQ/IqVFhruZ3cQPVRsgMVkLntV0ygAujo6GDXrl0cP36cmZkZpqam+O53v8ueaJRgjZZXkd2D81LZpTJ0G1BE6o6/MO3LgA8Hjn8T8I96VrqXc/LUKd6OU77pVb79/Thp7V62oVLZK0c9KxGpK/7sv1g6zV+R/Vf3EeAtwFnfvkgkwm233Vaxc2hpaSGNM9cqmPe3i/PzspTKXjl1F6yMMW3GmC3GmLj776yTcQq1NcZ0uvu2GGNGjTHFTx4SkZrxZ/99EvjFwPHfA37ifu2lpicSiYomOvT19RGJRDiCkw5/2nesDSelPdLczPr16yv2ng2vnJLttdiAMaA98Lit1LbAFt/+NuAY0DnX+zfCEiEiYdba2moBuzbPrMIh3zIggN20aVNVlu8ILjvyrjznMuAuK6JlRBw00hIhxpg2nODjn/AzCeT0igq1dXtYmVva1toUcCDf64hIOHgJFVNTU7QBdweO/zNOPT9PU1MTu3btqkrqeEdHR9a8rLtxsgH9PgpcpmVEKqaughXQDaQC+1LA6lLaWmsngOCs1/Y87UUkBILllO4EnuM7/jTwNrLHqao9XuTNy+rv76c1FuNtwJT//YFP43zIahmR8tVbsGrDqUPpdwTIVwOoYFtrbabkgTGm3d3/YMXOVEQqIlhO6WVAcMW1bcC/+x5HIpEFGS/yp7b/5oYN3L5kSdbxVwCbfI+1jMj81VuwgvyBqdy2w8B17u3AHMaYfmPMAWPMgSeffLKEtxeRcvkTKpZCzvymf8cJVn6Vzv4rxt69exmemWFfYP+fAs92v9bcq/mrt2CVwukx+V1Mbg+q6LbGmC3AoHtrMC9r7Yi1ttta233JJZeUftYiMm9eOSVw6vP9SuD4APCM+3W1sv+K4c2pugk47tt/IbA9TzspTb0FqwPk9pbacLL8Sm5rjIkD416gKpQGLyK14X24X0Z2AgU4KeJf9z0utUJFJXljZD8F/jhw7PeAlwfaSWnqKlh5WXvuGJOnGxgHZ+zJO1ZE2x4g5QtU7e5xEQkR78P9TuAC3/6fA7f7Hsdisapl/xXDm3sFcA9wKHD8Hpzit1pGZH7qKli5eoG4N9EXuNk31jQADM7V1g1MY8CYMcYaYyzORHRVqBUJmb6+Pv73kiWsDey/HWdyJCxcQkUhmzdvzgSrGbITKwCudvdpGZH5UdX1EqnqusjC8Japv/8zn+Hh6Wmu9R17BGdJDu/TK9+SH7Xgr1mYTqf5DM4tQM9xYBXnF4AMy3kvBFVdF5FFxz+v6o2BQAXOGlWW2iZU5OOfexWJRNhCbrLFnb7HSmUvnnpWJVLPSqS6kskkV199NdPT0ywHfkB2/b8Ezv39lpYWbrzxRm677bZQBKqgWCzG1NQUtwL/L3CsE3jU1+748eMsdupZicii4p9X9R6yA9XTOIPSkUiEG2+8saYJFXPxshg/DnwvcGxHnnZSmIKViISKN6/quWRnS4HTQ5mkPibXelmMMzhB1+/VwP8JtJPCFKxEJFS8nsYfAyt8+58EPpSnXVj5U9m/7G5+dwHLtYxI0RSsRCRUWlpauBynEoTfB4ATgXZh5k9lB6d3NeM7/svA286e5dOf/rTmXRVBwUpEQqWvr48PG5O1+u9h4D7f4zDMq5pLcBmR7wGfCrT5ILDk1CnNuyqCgpWI1Jy3VlUsFuNbu3ezNpCl/H6qu0x9tfhT2VtaWvi/gP/m5bOA96ElRIqhYCUiNeWfUzU1NZVTQf1Rzq/dE7Z5VcXwlhFZv349RyKRrKK2AO/kfMaj5l3NTvOsSqR5ViKV459TBU6W3MOBNr8JfMUYWltbWb9+fWjnVc3Fm3cVBf4DeJ7v2KdwFo/02i3GeVeaZyUidcs/pwpy16X6GvBwczMbN27k+PHjoZ5XNRcve3Ga3KrsNwJXBtpJNgUrEakZ/1pVv41T789vK5A+ezb0c6qK4c9e/Azwfd+xJuAjedrJeQpWIlIzXi+iiew5VAB/DXwr0K6e+eddzeAEYr/fAl61ZEnosxxrRcFKRGrG60W8hfO3wQDOkb3Q4mLobQTnXf0t8E+BNh+emeHTn/qU5l3loWAlIjXT19fH8uZm7gjs/yznb5PVw5yqYgTnXUFuOamXAj3T05p3lYeClYjUzObNm3lbUxP+lIk0zmRZT73MqSpGcN7VP+H0sPy2Aec07ypH0cHKGPNiY8yLZzn2ptmOiYjMpuN5z+OuWCxr3ydxitXW45yqYvjnXUUiEd5HdhmmK4Dfd7/WvKvz5gxWxpj3GGNmgIPAQWPMjDHm48aYVq+NtfYhp6mZmfWFRESC7r2X5U89lXl4BviwMcRiMfr7+zl06BBr1qyp3flVkZcJ+RhOdqDfHcBy6qO6/EJpLnTQGHMv0AO8F5hwd3e5+1LGmO3W2q0A1tpHjTGmmicrIovIyZPwkY9k7brgXe/iRw3Sk/BnOH4AeCtwgfv4ecA7gCEWRyZkJczaszLGXANgrb3cWnuXtfar7rbdWvtaYCUwaYx50BhzkzHmwoU6aRGpb8lkkr+57jp48snMvjPNzTzxlrfU8KwWlj/D8Sfkrib8XpwP2cWQCVkJhW4DXmetvWW2g9ba49ba+6y1a4FRoJvc5BYRkSz79+/nlVddxSu//e2s/TvPnePKV7+6YTLg/POuAO4EjvmOtwF/1NS0KDIhK6FQsHq82BdxA9dXrbV3VeCcRGSRSiaTxONxNp0+zUW+/ceAoXPnGioDLjjv6hi55aY2nTvHVz/5Sc27onCwUoVbEamoHTt2sPKZZ3hnYP92wCvd2igZcPnmXX0M+LGvzTJg8PRpzbtC86xEZAF461Xt3r2bLWfPZi1X/zOyx2saKQMuOO/qDLlFbn8PWKV5VwWD1a8ZY2IFjmcYY17jzrV6oELnJSKLhH+9qucDA4HjH8KpRO7XSBlwwXlXfwF8z3e8CTJrYDVKrzOfWdezMsZcBtwLxK21U3mOvwZYC1wEbLPWfscYM2OtXVLNE641rWclUrzgelWf4vyEV3AGxn8Zp2qF32Jd06kQb70rgDcAXwwcfx3wFer32lRtPStr7ePAQ8ATxpgH3MnB29yvjwDDwIPW2hustd+Z7wmIyOLlX6/qCiCY1/ZBcgPVYqkFWCp/b/JLwD8Ejt+F84F94sQJYrFYwyVdzLlSsDGmByerstPdNQEMW2vv87W5DIgDA9bay6t0rqGgnpVI8fy9hQTwJt+x7wNX4VRY94tGoxw6dGhRlVgqhv9agfOBezDQ5g+AT7tfRyIRIpEIiUSiLqp8VH2lYGvtuLW221rb5G7d/kDlanMnDi/qQCUipfF6Cy8hO1ABvJ/sQLVYawEWKzjvagIIppl8CIi6X6cbLOmiItmA1tpHK/E6IrK4eNUXtgf2fxv4QmDfYq8FOJfgvCtwAvoZ3+PnAe8OPK9Rki7yBitjzJ1uAoWISMm8VPUzZ87wW8ArAse3+L6ORCJs2rSJXbt2NWSPypNv3tWPgD8PtBsEnu173Cip/rP1rIaB1xpjDhhjdmv5DxEplj9V/Vw6zVDg+JeAf/Q9XkzrVZXLP+8q5i6dsg14ytemBfhw4HmNkOqfN1hZax+31r7XHQwbAW4xxjziZgNeupAnKCL1wyunND09TTqd5g+BVb7jMzgFWkFjVLPx5l0dP36c1tZWTpC9GCU4iRZdvsfnzp1b9NmBxSRYPGqtvcVaey0wDmw3xnzZrbRe1KRhEWkM/lT1FZCzXP2ngMeApUuXNvwYVTG8pIthnOvmaSK3SvtiL8k0Z+r6rE805k04k9EtTir75yt5YmGl1HWR2fnTr/8v8Ce+Y9PAC4D/pn4nti40/6TqHmAscLwP+GxgX1hT/6ueuj4ba+1D7rpWa4GLjTFfcScMKzFDpEF5YyfPITuJAmAnTqDyt5PC/EkX/xiJ8NeB40OQVWcRFm92YNmp6751rV6Lczu6S4kZIo3JS1W/EycRwPMk2enrWlCweP6ki83A075jzwPeF2i/WLMDK1p13U3MuMuXmPFmFbcVaRx9fX28dMkSbgzsvwM44X7dqOWUyuElXTxuDDsCxzYD7YF9i7HnOu8xq0alMSuR2SV/+EOOrlrFtefO16b4LnANTiYghHdMpR7EYjHOTU3xA5xelWcf8PpAu7CNCdZszEpExONNAr7r6quzAhXAO3EClVLVy9fX18czkQiDgf3X4xRnhcXbc61asDLGvKdary0i4eFNAv7cfffxgTNnso49BHzNGGKxmFLVK8AryfRZ4OuBY3cDMRbvJOuKBCtjzM3GmKPGmCPudhRyJq6LyCLjnwQ8ePYsz/UdOwO8B1i+fDkTExMNX06pEvzZgbc2N/OM79gv4FS2OH36NNdcc82imyRcsZ6VtXaltfZid1sJ3FKp1xaRcPImAb8ACP4t/1HgCRZvKnWteNmBrxgY4O6lS7OObQSutZapqalFN0m4IgkWxpjrrLVfDeyLWWtPzPacMt6rDegHJnGSYMattRPzaWuM6QS24kxqHi/m/ZVgIXKeNwn4YeDVvv0/wVkBeNrXLmwD/ovB5GOPYa+6ig7f5/h3gG7Cl9ASlgSLpDHmd40xL/Y2qncbcBRIWGsT1trtwJAblEpq6y4quZLcrE8RKdLJkyf5fbIDFcDtnA9UXjupvI/ecw9vb8r+GH8x8C7f48XSs61UsHovzty07b5tbYVeO8MNNO3W2knf7kmgp9S27qKS48DRSp+nSKO4NBrlo4F9XwY+F9inScDVsXfvXv5+Zian5NKf4vRsYfFMEq5UsBpzVxB+rbfh3H6rtG4gFdiXAlaX2VZE5uFzv/ALXOx7PA1sCLRZrKnUYeD1WN8NHPPtXw58BlgSaFfPKhWsjuXZV400lDZye0JHcG7nldNWREr1la/waz/8YdauDwKPB5ot1lTqMPB6rD8H3hE49us4t2P97epZpYJVh1sL8CZ3u5nqjVmVEmwqEpiMMf1uvcMDTz75ZCVeUqS+TU/DLdkJv4dwitV6NAm4+rwlRAD2Al8IHL8DuAonnb3eU9krFawGgOPARe7WBll3Byol5b6238XkH3cqpW1B1toR9zZn9yWXXFLq00UWlWQyycMveQk8fr4PdQ5IvO51RGMxmpqaNAl4gXiThD23kL2q8FKc24Gk03Wfyt5codcZzJO6XlQqeIkOkNtbaiN3mZdS24pIEfbv3889v/M7fOnpp7P2f7ypiR3f+AaJRELBaQF5k4Tj8TjpdJqfp9PcAiR8ba7BWVvsj9Np0uk08Xg8FKnspapIzyoYqFz5xrHKfZ8UcMAY408378ZZwRhjTLt3bK62IlIcr+5fS0sL666/nnsDgepxYOu5c0xPTxOPx+v6VlM98i8hEolEeAj4y0Cb9wGvdL+u11T2eU0KNsa8xlr7sO/xTcEmQNxa+7oyzy/fe8860dcYMwS0WWsHimjbiZPGvhWnFzbmzsUqSJOCpZHs378/81d7Op3mL4G3+I6fw/kQ/Cf3cSQSob+/n127di34ucr5SdoXAY9BVvmr/8KZg/UUtZmkXe6k4PkGqwPATdba7/geB9etuqGcEwsrBStpFP4l1QHeDPxVoM02chf/U7WK2mlqasL7TF8NfCVwfD/OUiKmqYmZmRkWUrnBatYxK2PMEeA11tp/DR7L84Y3W2sfDTxft9tE6phX9w/gfwEfDxz/DvCBPM9bDHN66lVLSwtTU1OAMzj/EbL/mFiDU1x4uA5T2QuNWV0ETBhjfmeuFwkGqtn2iUj4eWNUu3fvJp1O04zTo7rI1+YMsA5I53n+YpjTU6/8qewAfwx8M9DmI8CvnjpVd6nshYLVCM54zkPGmM3Bg8aYC40xtxtjtrl1AWNVO0sRWRDe2lR79uzJ7LsLeHmg3XuB7+d5vqpV1FYwlX0GZ4zxiK9NM7B3ZoaHdu/m8ssvz0w1CH3wstbm3YDd7r9xnHHUjxdoezvOdfn72doslq2rq8uKLEaHDx+20WjUApntBrA2sO0Ha3xt/Fs0GrWHDx+u9bfS0Pbt22ej0aiNRCKZ/5fX5/l//DrYiO//LhKJ2Gg0avft21eV8wIO2DI+e+dMXbfWJnBSvt9sjPn7fD0oa+1dOCXBVHdPpE75x6gArgD2BNr8J87tv2BalqpVhEcwlR3g7yCn4PArAH/OZjqdDvX0g0LBKjOh1jrp3t3AC3DmLj0/2NhaO4JTxUJE6khwjAqgBWdJev/o09M4t1mCJWBaWlpUrSJkOjo62LVrFxdccEFm31YgOCG2H2fBRr+wzsMqFKyy1nmyzlIbnTh/XE0YY341z3OU0y1SR/KNUTUB9+P0rPxu5fwvuNeT2rdvH1NTU1qyPqT8mZlncdZtCvaZ7iZ7PbKwLilSKFh1BntQ1trj1trVONU88mUKBpfkEJGQSiaTxONxpqens27//Rnw24G2nwbu8z1WT6o+BDMzjwK/BUz59jXjrFL7At++EydOhC7polCwMkBiljGqAZxeZSKQKaiFDEXqRHCMCuCd7ub3Hc7fKopEImzatEk9qToRTGUHJ4tzHU7WnOdinEUzn+3bNzU1Farit4WCVQfwILDHGPOeYNCyTmmiG4C7jDHB+YIiEnJ79+7NCla/i9Or8vsJ8AbgtPtYa1PVl2Aqu+eLwPsD+y7DqXDR6tsXpqSLWYOVtfZxa+1d1tq1OHcActaGCmQKfpnseYMiEkJeQoVX6QDgJTjrIfk/EE7glOb5L5TtV6+8quzRaDQnaG0DPhFofw3OmlhLA/vDkHRRVNV1d6zqiVmOeZmCl+MkC4lISOVLqHgxTmrzcl+7szi/zIfcxxqjql/+VPZYLIYxBgBjDAPAlwLtr8NZA8v49qXTae65557ajmGVM0nLv+GsFfVgpV4vrLyfEv4AABN/SURBVJsmBUu9yjfp9yqwT+WZMPoHvomimzZtqvWpS4UdPnzYbtq0ycZiMRsF+895fgaG80z+LmfiMGVOCp5X1fVGpqrrUq82btzInj17MuNUVwL/AATXvr4D+KD7dTQarcuF+qR4sViMpVNTfBNYFTj2cWBTnufM5+ei3KrrlVrWXkRCzp9QsQpngmgwUN2FE6g0RtU4+vr6OBGJ8Js445N+G4E/z/OcWoxhlRWsjDHbKnUiIlJd3gTRF+P0qJ4dOP7nwBb3a41RNQ4vY/A/gdcAPwscfye5pZpqMXG43J5VT0XOQkSqrqWlhVcD/wg8J3BsF+AlpMdiMc2jaiD+jMHHIxFeA/w80GYzudMaFnrdsnKDlZm7iYjUkpeq/obpafYDwVn+u3FKKYGW+GhU/ozBf8PJCHwq0OY2nHWjvKCx0OuWlRuslJ0hEiJeYIrFYjQ1NRGNRlm1ahWR4WH2zsywLND+o2QPoGvSb+Pyit9u2LCBH0Qi9JBbkuhm4LNABDh9+vSCprIrwUJkkfDPoZqamsJay9nTp7n77FnuPncu55f9PTgL0VmUUCHneWNY/wp5bwm+Gfg8sCSdXtByTApWIotAvqK0lwDj5FkCAlgP7HAfL126VAkVkuEfw/p+JMIrgB8H2rwBp5ZgywKWY1KwElkEgkVpO3GW83hloN0pnKrbe337LrjgAiVUSBb/GNbjkQgvB34YaPNK4JvALwHT09NcccUVVb0tWNakYGPMI9baayt4PqGnScESRrFYLFPrrx8nDX15oM3jwBuB7wb2NzU1MTMzU/VzlPrk/Ww9GxgDrgoc/ymwBvhXnNvJkUiERCKR00vXpGCRBuYvSnsRzkJzw+QGqq/iFPAMBipY+KwuqS9eivr/4PSmvh44/lzgG8Bqzldpv/7662lpaaloT0vBSqRO+RMqXo6z7tSb8rTbCbyO/IvNKVVd5uL/YyYFvBZn7Si/VmAfsMG379SpUxVNwFCwEqkjXk+qpaWF66+/nnPT0/xJOs3XcMYO/KaAPuDdwGw3+ZSqLnMJLuD4NE5GYHCScDNOLcGPAUvcff71sCBn5kRJNClYpE74e1KnTp3ipcCjwHs5/+HgOYCTZPHZWV5LqepSrHwLOFqcqhbvInvFYYC34yw5c6Fvn5v8E6zwVZJyg1Vvmc8XkQKCPanp6WmWptPsxBknCFbJBqcY7UuBw+7j5uZmlixZQjQapampiVgsplR1KVqhBRzvBn4bCBZeeh3wLeCF7mM3WAXrJpekrGBlrX28nOeLSDZ/BQpjDJdffjn33nsvp06dAuAG4N9x/qIN/vL+GOdDYgvOXCrPwMAAP/jBDzh16hQzMzMcP35cqepSEn8qezBgfRF4GfCjwHNWAd/GmZNVCVrPqkRKXZdq2b9/P/F4nHQ6nTVnCuBFOGMBr5rlucM4QepEYH8sFuP48eOVPlVpYMlkkquvvprp6ems/c8GvgD8Rp7nfAD4E8BaO++hIyVYiIRAvgoU4FRH/zjO2NSr8jxvEqckzi3kBipl+kk1zHZb8H+AVwN/kec5d1TgfRWsREIgWIHiQuDDQBInHbg50P4Z4E6cCZr/MMtrKtNPqsV/W9Cf2v40cCPwDuBshd9TwUpkgQQrosdiscykSW8V3wuBrTg9pvcB0Tyvsx/ntuBWYDrPcWX6yULwqrRPTU2xb9++rJ7Wx3AWO3yygu+nMasSacxKipVMJtmxYwd79+7NlEIyxuD/nfPK08Smp3knTi/qwvwvRxJnTaEvFnjPlpYWbrzxRm677TYFKllQyWSSnTt3cv/993PihHNT+hdxKrR348xzKmfMSsGqRApWUoxCyRJ+L8JZT+r3gQtmafNTnMHpT5Cd5ecpVI9NpBY2btzInj17SKfTLMOpVbkBBasFpWAlc5ktW8oTwSmLtBF4RYHXOQ4M4cxlyf9K6klJOM32O6BsQJEKKzS+NJdgsoTnN3Du5f8X8FfMHqieBP4IuBTYRm6g8sak9u3bx9TUlOZMSegUmkg8b9ZabSVsXV1dVha3ffv22Wg0aiORiMWpLGMBG4lEbDQatfv27bPWWnv48GG7YcMG29raao0xdsWKFfbKK6/MtDdgfw3sh8Amwdo5tifAvh3sct97AtYYk/k3FovZTZs22cOHD9f4KonM7fDhw3bTpk02FotZJ9zM/7NXtwFLpNuA9cmf7HDy5ElaWlro6+tj8+bNWb2SuW7heZYtW5ZZA+rs2fNJuitx5kO93t2KKYb2MM5cqr8hN903Fouxfv163eaTulfuelY176nU2xbmnlXwL/3W1la7YcOGBf0rvBLnUOnvo5iekveewTaFNgO2A+wNYHeBPVRE78nbjoK9G+yqAq/f1NQ0r+9XJIyAA7aMz96af/jX2xbWYFXsratan8NcgWi21/DfCvNut61YsWLWYOa9z4oVK4oPPu575NueA/aVYG8CuxPs18CmSghOFuxpsA+CfSPYpUWcTywWq/r/mchCUbBSsLKHDx+20Wh0zg+/FStWFOyh5Askb33rW+26devm7OUUcw7Lli2zy5cvzwlEzc3NdsmSJXbZsmVFB5bg5g+IswW82XpHzwL7IrDXgV0HdovbU/pbsI+CPV5iUPJvKbAPgF0PNlbi97Np06aF+PERWRAKVjUIVvO5TVXqc+Zq7z8+3w91v1I+4PO9Rqm30KqxtYJ9Pthr3MDTC3YA7Fawd4H9BNgvgP062O+D/TnYs2UEonzbDNiDYP8M7GvARgLnWKj35t+i0aiSKGRRUbBa4O3yyy8v+XZbqbfHZvtg89rfcccdRQeWQtvy5cvtlVdeaZcvXz6v5/s/UEsNmnNtF+KMB/062NeDvRHsZrDbwN4H9vNu0HkM7M/APlPhoFPsdgzsw2DvBLvGPe9C35eXyReG27YiC6ncYFV32YDGmDagH6d8Wjswbq2dKLVtKa/jt2TJEnvuXHBtzGzLly+nvb2dJ554IrMOUZHfG/X0/xGJRLjhhhtobW1l9+7dWceWACuAGNDqbjHgIndrc/+92N1Wuv8+y90qNDOjYk4C/wH8wP33EDABPFHk8yORCP39/ezatSuzz1+exstQVOafLFblZgPWY7AaAwastZO+x73W2lQpbUt5ncBrVu2CNePUhYsFtlagxd1acYLAct92AbDUt0VwgsUSnFnfwSXPPRZnSep8/3qbX1Ngi/jeb6l7Hstxiq+GLdgUkgJ+hrPEwc/c7Ufu9mP335+W+R7RaJRDhw4pCEnDKjdYBVceCDW3N9TuBRjXJE6B30SxbY0x48W+TtAKnHpu3od5E85F9LalZAeRKNmBpoXsQHShb8tXYVuKdwo4ChwBjrlfH3MfH/U9fsq3HSV/vb1K8dftU6ASmb+6ClY4xXuDPZ8UsJrcIFOobaqE18myCvhW8ecr83QSJ5g8SXZwOZLnsReMnq7SuSxbtoympibOnj2bVUapubkZay3Lli3jzJkzRKNRLr30Up544gmmp6d1W0+kguotWLXhfC75HcEZcyqlbSmvgzGmH2d8i67SzrdhzeDUtDsBTLnbCZy/CI5x/q+FI2QHnCfdr88s/Cnn8PeKXvjCF2p8SaSG6i1YgTMWX4m2Rb+OtXYEGAHoruKY1QzOB/px919v8z7sT7rbtLuddrczOL2KZ9wt7b6Wt3njUH7Gt3ljUP7HwdLIxn0d/+ulfe+Xds/DO69nyroSuSqRfOIFn8HBQYaGhnKW7/DewxhDa2trTjDatWtXVoKEiCycegtWKZxekd/F5PaS5mpbyutkOQX8C+c/2Gdw6rl5W5rzAcQLJv5Ac5LzQei4e+y4u52c682LFI1G+djHPsatt94653pKYdDc7PwY+mvsFerVeLfbJicn56zhB+QNPuvWrVNPSaSelJP3vtAbToA5Ftg3DMRLaVvK6+R53YrOJ6rklq8q+KZNm2xLS8u8Xq+lpcX29fXZdevWVe18vXlHXmXmpqamkiqLa76SSH2g0SYFA2M4mXze44NAm/t1e+BYobazHpvj/Rcs+BRb7cDbCn3AF1uhYrYP+WpUqKhUlYZygp2ILIxGDFZtwBacHtIWoNN3bAgYLrLtrMfmeP+sD9zZ6t0VCgLFFGv1PnDXrVtX1GsXU0cu+KHe0tJiX/SiF9mWlpY5P+SLrT/o/15nq7ShXo9I42m4YFXrLd8H7nyCQLG9gWKCxELVkSslyPprGKrXIyLlBqu6q2BRa+4H8oIOxu/fv594PJ6TLOFPQlizZk3VzwNUIkhE5qfhyi3VWq1WClaQEJF6pmC1wLSsvYhI6coNVk2VPBkREZFqULASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQa671CZTCGNMG9AOTQDswbq2dmE9bY0wnsBUYttaOV/vcRURk/uoqWAGjwIC1dhLAGDNmjOm11qZKaWuM6XHbtC/MaYuISDnq5jag21Nq94KPaxLoKbWttXbc7U0dreIpi4hIhdRNsAK6gWAPKgWsLrOtiIiEXD3dBmwjtyd0hPy38kppOydjTD/O+BfA08aY783ndSTHs4Cnan0Si4iuZ2XpelbWL5fz5HoKVgArq9S2IGvtCDACYIw5YK3trtRrNzJdy8rS9awsXc/KMsYcKOf5NQ1Wbo+la45mQ+7YUwqnx+R3MfnHnUppKyIiIVfTYOX2WIp1gNzeUhswVmZbEREJubpJsHDT0w8YY/zjTt3AOIAxpt07NlfbMpUSYKUwXcvK0vWsLF3Pysq6nm7WdtGMtbayp1NFhSb6GmOGgDZr7UARbTtx0ti34vTCxqy12yt0jl4qfS8wOMscMCmCey0HrbXK4pwHXb/K0u92ZRljjuEMzYx7n9sF29dTsAo7NwgOWGsHjDHDOEEwUevzqmfGmDF92M6frl9l6He78owxPaVUD6qb24DVYIzpNMaM+v5i8va3GWO2GGPi7r+dxbyetXbC/WFuwxkza5gyTpW+lpJN17eySr2ejfy7XYx5/ny2BYZqCqq31PWKmaPkUillnfLpBh6pwGnWhSpfy4an61tZZV7PhvrdLkYZ13MlcNQdwtk2189swwYrr/tpjMlKZ5+jVFPCTbfP93ojvq/HfYVyByt97mFTzWsp87++C3eG9aWc69lov9vFmO/19H7PjTGPUMT1bNhgVUChUk2JQh+kxpgtQMptkwI6qnaW9WHe11KKUvD6Lvzp1L1Zr6d7u0q/26UpdD0BVvqu55waesxqFrOVaiqmIsYIMOl2i7vQX17lXEuMMXGg3RjTX2qaa4MoeH11/UpW6Hrqd7t0s15PNznlqPszuhrYNteLqWeV37xKNbn3XL2BVw3AOuZd9sr9gVYPobBZr6+u37zkvZ763Z63uX4+ocifUfWscqlUU+XoWlaXrm9l6XpWVkWvp4JVLpVqqhxdy+rS9a0sXc/Kquj1VLAKqHKppoaia1ldur6VpetZWZW+ng1bwaJQyaVCpZokl65lden6VpauZ2Ut1PVs2GAlIiL1Q7cBRUQk9BSsREQk9BSsREQk9BSsREQk9BSsREQk9BSsREQk9BSsROqUMabdXQtIZNFTsBIJCS/4uFXS40U8ZQBf6Rr3uQeNMdYYM+xftdV9zTH32Ohsa4mJhJUmBYuEhDHmINCLE4R6rLVdc7UPtnGD0JC19qI87TuBg8BFWklY6o2WCBEJATeQtFtrJ40xSeYo9um2P7AgJycSAgpWIuFwA26BzyJXUB4Ahqt6RiIhojErkXDoobSlE7pVYFUaiXpWIjVkjNkCdACdwGpjTBcwXCgQuYkTFVm2wr2d+FWcZcUn3d3twBAa25IQUbASqSFr7XY3YPRba3uLfNoAMFjgeJsbBIM68uxbCdzsW2IcY8wYMKhAJWGiYCVSe92c79UUo81aW6h9yltPyM8LisHXwtdLc7MJV+Z7vkgtKViJ1F4XUNT4kzv/arSC7z3u9aDcFV2H3PMRCRUlWIjUXjfwSJFtB4AHK/XGgVt9ozi3/0rp5YksCAUrkdrrpIiECXeJ8GCAqQhvjMufNu/eNhQJBd0GFKkhLyAUmYa+lirMrXJv/23Fd/vP3bey0u8lMl/qWYnUVinJFb3+rL0Kynf7Lw4crcJ7icyLelYitVVUcoXb0yl4+8+twN6Dk7o+DIxaa8fdY/04dQcB7jPGPGCtTbj724GjbvLGSvec+smf6i5SEypkK1JDbvHabXP1mNwxpQkv+Ig0GgUrkQXm9mBS1tpxY4y11poinpNTYV2kkWjMSmTh3Qd0umWT5px8qwrrIhqzEqkFr1TSamttobJJnhtQhXVpcLoNKBJyxpjREuoGiixKClYiIhJ6GrMSEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQU7ASEZHQ+/+tIycktsayEQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.semilogx(freq_vec, -Z_exp.imag, \"o\", markersize=10, color=\"black\", label=\"synth exp\")\n", + "plt.semilogx(freq_vec_star, -Z_im_vec_star, linewidth=4, color=\"red\", label=\"GP-DRT\")\n", + "plt.fill_between(freq_vec_star, -Z_im_vec_star-3*np.sqrt(abs(Sigma_Z_im_vec_star)), -Z_im_vec_star+\\\n", + " 3*np.sqrt(abs(Sigma_Z_im_vec_star)), alpha=0.3)\n", + "plt.rc('text', usetex=True)\n", + "plt.rc('font', family='serif', size=15)\n", + "plt.rc('xtick', labelsize=15)\n", + "plt.rc('ytick', labelsize=15)\n", + "plt.axis([1E-3,1E5,-0.01,0.03])\n", + "plt.legend(frameon=False, fontsize = 15)\n", + "plt.xlabel(r'$f/{\\rm Hz}$', fontsize = 20)\n", + "plt.ylabel(r'$-Z_{\\rm im}/\\Omega$', fontsize = 20)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}