Simulating datasets Shannon theorem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['seaborn-dark', 'seaborn-darkgrid', 'seaborn-ticks', 'fivethirtyeight', 'seaborn-whitegrid', 'classic', '_classic_test', 'fast', 'seaborn-talk', 'seaborn-dark-palette', 'seaborn-bright', 'seaborn-pastel', 'grayscale', 'seaborn-notebook', 'ggplot', 'seaborn-colorblind', 'seaborn-muted', 'seaborn', 'Solarize_Light2', 'seaborn-paper', 'bmh', 'tableau-colorblind10', 'seaborn-white', 'dark_background', 'seaborn-poster', 'seaborn-deep']\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import langevin\n",
    "import lmfit as lm\n",
    "from lmfit.models import ExponentialModel\n",
    "from scipy.stats import pearsonr\n",
    "from scipy.optimize import minimize_scalar\n",
    "\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = 18\n",
    "BIGGER_SIZE = 20\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "SEED = 35010732 # from random.org\n",
    "np.random.seed(SEED)\n",
    "\n",
    "print(plt.style.available)\n",
    "plt.style.use('seaborn-white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate A and B from the dataset\n",
    "def OUanalytic1(data):\n",
    "    N = data.size\n",
    "    data1sq = data[0]**2\n",
    "    dataNsq = data[-1]**2\n",
    "    datasq = np.sum(data[1:-1]**2)\n",
    "    datacorr = np.sum(data[0:-1]*data[1:])\n",
    "    coef = [(N-1)*datasq,\n",
    "       (2.0-N)*datacorr,\n",
    "       -data1sq-(N+1)*datasq-dataNsq,\n",
    "       N*datacorr]\n",
    "    B=np.roots(coef)[-1]\n",
    "    Q=(data1sq+dataNsq)/(1-B**2)\n",
    "    Q=Q+datasq*(1+B**2)/(1-B**2)\n",
    "    Q=Q-datacorr*2*B/(1-B**2)\n",
    "    A = Q/N\n",
    "    P2A = -N/2/A**2\n",
    "    Btmp = (N-1)*(1+B**2)/(1-B**2)**2\n",
    "    tmp = (2+6*B**2)*(data1sq+dataNsq) + (4+12*B**2)*datasq - (12*B+4*B**3)*datacorr\n",
    "    P2B = Btmp - tmp/A/2/(1-B**2)**3\n",
    "    PAB = (N-1)*B/A/(1-B**2)\n",
    "    dA = np.sqrt(-P2B/(P2A*P2B-PAB**2))\n",
    "    dB = np.sqrt(-P2A/(P2A*P2B-PAB**2))\n",
    "    return A,dA,B,dB\n",
    "\n",
    "def OUresult1(data,deltat):\n",
    "    A, dA, B ,dB = OUanalytic1(data)\n",
    "    tau = -deltat/np.log(B)\n",
    "    dtau = deltat*dB/B/np.log(B)**2\n",
    "    return A,dA,tau,dtau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate A and B from the dataset\n",
    "def OUanalytic2(data):\n",
    "    N = data.size\n",
    "    data1sq = data[0]**2\n",
    "    dataNsq = data[-1]**2\n",
    "    datasq = np.sum(data[1:-1]**2)\n",
    "    datacorr = np.sum(data[0:-1]*data[1:])\n",
    "    coef = [(N-1)*datasq,\n",
    "       (2.0-N)*datacorr,\n",
    "       -data1sq-(N+1)*datasq-dataNsq,\n",
    "       N*datacorr]\n",
    "    B=np.roots(coef)[-1]\n",
    "    Q=(data1sq+dataNsq)/(1-B**2)\n",
    "    Q=Q+datasq*(1+B**2)/(1-B**2)\n",
    "    Q=Q-datacorr*2*B/(1-B**2)\n",
    "    A = Q/N\n",
    "    P2A = -N/A**2/2\n",
    "    Btmp = B**2*(1+2*N)\n",
    "    tmp = (1+Btmp)*(data1sq+dataNsq) + (2*Btmp + N + 1 -B**4*(N-1))*datasq - 2*B*(1+B**2+2*N)*datacorr\n",
    "    P2B = -tmp/((1-B**2)**2*(data1sq+dataNsq + (1+B**2)*datasq - 2*B*datacorr))\n",
    "    PAB = (N-1)*B/A/(1-B**2)\n",
    "    dA = np.sqrt(-P2B/(P2A*P2B-PAB**2))\n",
    "    dB = np.sqrt(-P2A/(P2A*P2B-PAB**2))\n",
    "    return A,dA,B,dB\n",
    "\n",
    "def OUresult2(data,deltat):\n",
    "    A, dA, B ,dB = OUanalytic2(data)\n",
    "    tau = -deltat/np.log(B)\n",
    "    dtau = deltat*dB/B/np.log(B)**2\n",
    "    return A,dA,tau,dtau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "A,D = 1.0,1.0\n",
    "delta_t_list=np.linspace(0.01,2,200)\n",
    "#M=1000 # number of data sets\n",
    "N=50000 # length of data set\n",
    "#print(np.exp(-delta_t*D/A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us first calculate the standard deviation for each point of the acf\n",
    "# all acfs are stored in acf_all for later analysis\n",
    "\n",
    "result_Bayes_array = None\n",
    "for delta_t in delta_t_list:\n",
    "    data = langevin.time_series(A=A, D=D, delta_t=delta_t, N=N)\n",
    "    \n",
    "    results = OUresult1(data,delta_t)\n",
    "    if result_Bayes_array is None:\n",
    "        result_Bayes_array = results\n",
    "    else:\n",
    "        result_Bayes_array = np.vstack((result_Bayes_array, results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a22f336a0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAD6CAYAAACccFGwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNXdx/HPZLLvKxAICUnAgyggKqAo1rWCu1Zbl1qx7hW12mrd2mJ9tFZrfdyXPrUuWKgFiyguCG6oiKAsLnjYl7AnBAgJ2ef5405iGBIIYSYzmfm+Xy9fI+fe3PtLvObLuefcc10ejwcRERF/iwp2ASIiEp4UMCIiEhAKGBERCQgFjIiIBIQCRkREAiI62AWEAmNMHDAU2AA0BLkcEZGuwg3kAnOttTW+GxUwjqHArGAXISLSRY0EPvFtVMA4NgC88sor9OjRI9i1iIh0CRs3buSSSy4B7+9QXwoYRwNAjx49yMvLC3YtIiJdTatDCxrkFxGRgFDAiIhIQChgREQkIBQwIiISEAoYEREJCAWMiIgEhALmAFXXNXDiXz9kzoqyYJciIhJSFDAHaFdtAytKK1m8YUewSxERCSkKmAOUEOsGoKpOS5iJiLQUEU/yG2NuAa4APMA84JrWFmbriLhoJ6OraxUwIiIthX0PxhgzHLgcGAYMxFn980Z/Hd/lcpEQ42aXejAiIruJhB5MOTDWWlsJYIxZBOT78wQJsQoYERFfYR8w1tolwBIAY0wucBMwxp/nSIhxs6u20Z+HFBHp8sIiYIwxFwJP+DRvstYe0mKfAuBt4Glr7Qx/nj8h1k21ejAiIrsJi4Cx1k4EJra13RhzGPAG8KC19nF/n19jMCIiewqLgNkbY0wO8A7OOMykQJzDuUWmgBERaandAWOMOQm4HxgEbAZeAP5krfXbb1ZjzFnAK9balFa2XQXcBuQBC4BbrLWz23HYm4FU4G5jzN3etmnW2rv8VDbxsW527Krz1+FERMJCuwLGGHMMzvjFv4A7gCOAe4FG4B5/FGKMGQGMB1ytbPsF8AzwJ2AucAPwrjFmsLV25d6Oa629E7jTHzW2JSEmik3b1YMREWmpvT2YB4Dp1tox3j+/b4zJAk7AJ2CMMRnAcGvtOy3aegFF1tpZvgc2xsThzOy6F6gEYn22u3CC5Tlr7T3etvcAi9M78dszLR2lMRgRkT3t80FL7xjGMcBzLduttbdba49v5UuuB6YaY0a3+PoZwOPesPA1GqdXdCvQ2gB8X6AAmNri3HXANGDUvurvDHoORkRkT+15kn8gzm2rSmPMG8aYamPMZmPMOGNMa1//Z5wZW5ONMWcD7wIpwDnWWk8r+88FCq21j+Es5eLrIO/nMp/2FUCxMcbdju8hoOJj3FoqRkTER3sCJsf7+RLwPU6P4yngbpxex268g/4XAbOAKUBv4BRr7arWDm6tXWet3baX86d6Pyt82iu89Se143sIqET1YERE9tCeMZgY7+e71tqmQPnAGJONMzPrr/uYSeZu53na0nRbzbd309Qe9EfoE2Lc1Dd6qGtoJMYd9su7iYi0S3t+G+70fr7j0/4ekAz0adnovWU1ATgWOANYCUw3xvTtYI3bvZ++U5eTccKlsoPH9Zv4GOcunXoxIiI/aE/ANI19xPq0N/VsfHsWd+IEy3nW2mnAqTgLTk5pY5B/X5Z6P4t82osA28a4TqdqeieMxmFERH7QnoD5DlgHXODTfjqwHljl0/4kcKa19l0Aa20pcApwXQfDYCmwFjinqcEYE+M9/8wOHM/vEtSDERHZwz7HRqy1jcaYO4EXjTFPA5OAk4HLcEKj0Wf/rcB0n7Z1OCG136y1HmPMA8ATxphy4FNgLJANPNKRY/pbU8BUqQcjItKsXSPS1tqXgItxxlWmAecD11prnw1gbS3P/xTOjLVLcQIuHTjVWruiM86/L/Gx6sGIiPhq9+wua+0EnMH7gLHWjgPGtbHtYeDhQJ6/oxJjNAYjIuJLc2r9IEE9GBGRPShg/ECD/CIie1LA+EHzczC6RSYi0kwB4wfNz8GoByMi0kwB4we6RSYisicFjB/8cIss6MuiiYiEDAWMH7ijXMRGR1FVVx/sUkREQoYCxk8SY/VOGBGRlhQwfqLXJouI7E4B4ydOwGgMRkSkiQLGT+Jj3HoORkSkBQWMnyTEuvUcjIhICwoYP9EYjIjI7hQwfqJbZCIiu1PA+IlukYmI7E4B4yeJMW690VJEpAUFjJ8kxLqprNWT/CIiTRQwfpKRGEtFdT11DXoWRkQEFDB+k5kUA8C2qrogVyIiEhoUMH6SkRQLQHlVbZArEREJDQoYP8lMdAJma6UCRkQEFDB+09SD2aYejIgIoIDxm8ykph6MxmBEREAB4zfpic4gv8ZgREQc0cEuoDMYY24BrgA8wDzgGmttjT/PERftJinWrTEYERGvsO/BGGOGA5cDw4CBgBu4MRDnykiKpVwBIyICREYPphwYa62tBDDGLALyA3GizKRYtuoWmYgIEAEBY61dAiwBMMbkAjcBYwJxroxE9WBERJqERcAYYy4EnvBp3mStPaTFPgXA28DT1toZgagjMymWFaU7A3FoEZEuJywCxlo7EZjY1nZjzGHAG8CD1trHA1WH04PRNGUREQiTgNkbY0wO8A7OOMykQJ4rMymGnTX11NQ3EBftDuSpRERC3n4HjDEmDlgAzLHWjvFnMcaYs4BXrLUprWy7CrgNyPOe/xZr7ex2HPZmIBW42xhzt7dtmrX2Lj+V3eyHp/nr6J6qgBGRyNaRHswfgf7AHH8WYowZAYwHXK1s+wXwDPAnYC5wA/CuMWawtXbl3o5rrb0TuNOftbal5Xpk3VPjO+OUIiIha7+egzHGDMF5hqR0L/tkGGNG+bT1MsaMbGP/OGPMbcAHwB5v7DLGuHCC5Tlr7T3W2reAs7w13Lw/9Qda84rKmkkmItL+gDHGRAPPAw8B6/ay6/XAVGPMaO/X5QAzgMe9YeFrNHAHcCvQ2gB8X6AAmNrUYK2tA6YBo1rZP2ia1yPTszAiIvvVg/kdEAv8eR/7/RlnxtZkY8zZwLtACnCOtdbTyv5zgUJr7WM4S7n4Osj7ucynfQVQbIwJmcGOjET1YEREmrQrYIwx/YG7gCuttXv97WmtbQAuAmYBU4DewCnW2lVt7L/OWrttL4dM9X5W+LRX4NSftM9voJOkJeitliIiTfYZMMaYKOAfwD/aOWvLl5sDmw7ddFvNt3fT1N54AMf2q9joKOKio9hZs8dQkohIxGlPD+YGnDGQPxhjor1jMQCuFv/ezHvLagJwLHAGsBKYbozp28Eat3s/facuJ+OES2UHjxsQKfEx7KhWwIiItCdgzgV6AVuBOu8/g4FfAHXGmD4++9+JEyznWWunAafiLDg5pY1B/n1Z6v0s8mkvAmwb4zpBkxofTUW1bpGJiLTn1tU17Nl7eAVnAcl7gPU+257EeQhzOoC1ttQYcwpQ1MEwWAqsBc4BpgMYY2KA03FmkoWU5PhoKtSDERHZd8BYa61vmzFmF1BmrZ3Xyv5b8QZBi7Z17H1q897O7zHGPAA8YYwpBz4FxgLZwCMdOWYgpagHIyICdJEXjllrn8J5TuZSYBKQDpxqrV0R1MJakRIXox6MiAgdnN1lrT3M34V4jzsOGNfGtoeBhwNxXn9KiY/WLDIREbpID6YrSYlXD0ZEBBQwfpfs7cE0NIbU5DYRkU6ngPGz1HjnrqNuk4lIpFPA+FmKN2A0k0xEIp0Cxs9S4p31yNSDEZFIp4Dxsx96MAoYEYlsChg/S47TLTIREVDA+F3TLTL1YEQk0ilg/KxpFplWVBaRSKeA8bMfejC6RSYikU0B42fxMVFER7nYqR6MiES4A3nTZJdhjLkFuALnrZjzgGustTWBOJfL5dKS/SIiREAPxhgzHLgcGAYMxHmF842BPKeW7BcRiYweTDkw1lpbCWCMWQTkB/KEWrJfRCQCAsZauwTn7ZsYY3KBm4AxgTxnim6RiYiER8AYYy4EnvBp3mStPaTFPgXA28DT1toZgawnJT6Gddt2BfIUIiIhLywCxlo7EZjY1nZjzGHAG8CD1trHA12PxmBERMIkYPbGGJMDvIMzDjOpM86pW2QiIu0MGGNMLPAH4FIgG5gD/NZa+5U/izHGnAW8Yq1NaWXbVcBtQB6wALjFWju7HYe9GUgF7jbG3O1tm2atvctPZe8hLSGGiuo66hsaiXaH/UQ9EZFWtbcH8whOuPwOWI4zzfcDY8wga+1qfxRijBkBjAdcrWz7BfAM8CdgLnAD8K4xZrC1duXejmutvRO40x81tlev9AQaPbBhezW9MxM789QiIiFjn3+9NsakAVcB46y1T1trpwMXADE4oeO7f4YxZpRPWy9jzMg2jh9njLkN+ADY476SMcaFEyzPWWvvsda+BZwFlOL0TkJOfpYTKqvLqoJciYhI8LTn/k0lMBz4Z4u2Opyn4uNa2f96YKoxZjQ0j4HMAB73hoWv0cAdwK1AawPwfYECYGpTg7W2DpgGjGpl/6AryEoCYPXWyiBXIiISPPu8RWatrQfmAxhjonB+2d+DEzDjW/mSPwNDgMnGmIuAPwIpwGhrraeV/ecChdbabcaYca1sP8j7ucynfQVQbIxxW2sb9vV9dKYeqfHEuqNYs1U9GBGJXPs7Av17nF/slwJ/sdZa3x28v+wvAmYBU4DewCnW2lWtHdBau85au20v50z1flb4tFfg1J+0P99AZ3BHucjLTGCNbpGJSATb34D5L3A8Tg/mD8aYe9vxNW4ObDp00201395PU3vjARw7YPIzEzUGIyIRbb8Cxlq7yFr7kbV2HPAYcKsxJqblPsYYNzABOBY4A1gJTDfG9O1gjdu9n75Tl5NxwiUkBzoKMhNZs7UKj6e1u4IiIuGvPbPIehhjLjfG+P6Cn48zyJ/l034nTrCcZ62dBpyKs+DklDYG+fdlqfezyKe9CLBtjOsEXX5WEjtr6imv0hP9IhKZ2tODSQeeB873af8xsNn7T0tPAmdaa98FsNaWAqcA13UwDJYCa4Fzmhq8vabTgZkdOF6nKMhsmqockh0sEZGAa88ssu+NMZOBh71P9K8AzsMZ6P+ltbbRZ/+twHSftnXAuo4UaK31GGMeAJ4wxpQDnwJjcVYUeKQjx+wMBd5nYdZsrWJIfkaQqxER6XztHYP5BfB3nOdV3gSOAi6w1v5zr1/lJ9bap3Cek7kUmITTqzrVWruiM87fEU1P8K8q1UC/iESmds3ustZW4SwT87tAFuOdPDCujW0PAw8H8vz+FB/jpiArkSWbfGdXi4hEBq3EGEAH90jluw07gl2GiEhQKGACaEDPVFaVVbKzRkv3i0jkUcAE0IDcVDwesBvVixGRyKOACaABPZ1Vbr5br4ARkcijgAmg3LR40hJi+G6DBvpFJPIoYALI5XIxIFcD/SISmRQwATagZyrfb9hBQ2NIrmgjIhIwCpgA69stmZr6RjbuqA52KSIinUoBE2B5GQkAlOjlYyISYRQwAdYr3Rsw5buCXImISOdSwARYTwWMiEQoBUyAxce46ZYSR0m5bpGJSGRRwHSCvIwE1m1TD0ZEIosCphPkZSTqFpmIRBwFTCfIy0hg/bZdbNpRzfOfrKRRz8SISARo1/tg5MD0ykigvtHDH1//lne+3Uh+ZiInD+ge7LJERAJKPZhOkJfhvN3ynW83AvDi7FXBK0ZEpJMoYDpB08OWAMMKM5m1tJRlm3cGsSIRkcBTwHSCpoctc1LieOKiIcS6o5jwxZogVyUiElgag+kE8TFuDs9P56SDu9MtNZ7BvdP4umR7sMsSEQkoBUwnee1XxzT/e99uybzzzcYgViMiEni6RRYExTnJlFfVUbazJtiliIgEjAImCPp1TwHQQL+IhDUFTBD07ZYMwFIFjIiEMQVMEPRMiycx1q0ejIiEtYgY5DfG3AJcAXiAecA11tqgDYC4XC6Kc5JZvkUBIyLhK+x7MMaY4cDlwDBgIOAGbgxqUTi3ydSDEZFwFvYBA5QDY621ldZaD7AIyA9yTfTtlsyG7dVUVNcFuxQRkYAI+4Cx1i6x1n4EYIzJBW4CXg9uVTCwVxoAny4rDXIlIiKBERZjMMaYC4EnfJo3WWsPabFPAfA28LS1dkZn1teaY/pm0zMtnlfmrGHUobnBLkdExO/CImCstROBiW1tN8YcBrwBPGitfbzTCtsLd5SLnw3N55EZS1hdVklBVlKwSxIR8auwv0VmjMkB3gFuDpVwafKzob1xR7n4lxa+FJEw1K4ejDHGjTN2cRXOAPlq4CngSe/AuV8YY84CXrHWprSy7SrgNiAPWADcYq2d3Y7D3gykAncbY+72tk2z1t7lp7I7rEdaPD8e0J2JX6zl1ycdREKsO9gliYj4TXtvkf0euB24F/gcGAn8L5AIPOiPQowxI4DxgKuVbb8AngH+BMwFbgDeNcYMttau3NtxrbV3Anf6o8ZA+OWxhbz9zUYmf1XCz48qCHY5IiJ+s89bZMaYKOAW4CFr7X3W2pnW2nHAs8BvW9k/wxgzyqetlzFmZBvHjzPG3AZ8ANS3st2FEyzPWWvvsda+BZwFlOL0Trq0IwsyGJSXxvOfrqSx0W+dQRGRoGvPGEwa8BLwmk+7BXKMMb6j09cDU40xo6F5DGQG8Lg3LHyNBu4AbgVaGyPpCxQAU5tPbG0dMA0Y1cr+XYrL5WLMiD6s2FLJgpJtwS5HRMRv9hkw1tpya+1Ya+18n01nAiXW2kqf9j/jzNiabIw5G3gXSAHOaWO8Zi5QaK19DGcpF18HeT+X+bSvAIq940Nd2hEFGQAs26Qn+0UkfHRomrIx5krgZFpZcsVa22CMuQgnZKbg3Mo6zlq7qrVjWWvX7eN0qd7PCp/2CpyATAJ2tLv4EJSXkUisO0prk4lIWNnvacrGmEtwBtwnsefDja1xc2DP2zTdVvPt3TS1Nx7AsUOCO8pFYXaSAkZEwsp+BYwx5mbgZeBN4JLWbnl5b1lNAI4FzgBWAtONMX07WGPTy+t9py4n44SL7y26LqkoJ4kVW8LiWxERAfYjYIwx9wN/wwmY8621tW3seidOsJxnrZ0GnIqz4OSUNgb592Wp97PIp70IsP58DieYinOSWb21itp6p0O2qrSSxRu69J0/EYlw7QoYY8xNODO9HgXGWGv3mE7cwpPAmdbadwGstaXAKcB1HQyDpcBa4JwW9cQApwMzO3C8kFTcLYmGRg9rtjq9mN+//g03TPCdVyEi0nXsc2zEuwLxX4Cvcdb7Gm6MabnLvJaBY63dCkxvuYN3IH9fg/mtstZ6jDEPAE8YY8qBT4GxQDbwSEeOGYqKsp3XKC/fUklxTjKLSrazs6ae2vpGYqPDfkUfEQlD7Rl8PxWIw3lZV2tLs+TgzBQLGGvtU8aYBJzlam7GWSrmVGvtikCetzMV5TiPEy3fspP129PYvst5T8yqskoO6r7HyjkiIiFvnwFjrX0BeCHglTjnGgeMa2Pbw8DDnVFHMKTEx9A9NY4VWyr5dt325vZlm3cqYESkSwqL5frDhemRypery+mVnoDLOx1Cr1UWka5KN/dDyJmDcllZWsl/5q2lMDuJXukJLFXAiEgXpYAJIacNzCUp1s367dUc0jONvt2S1YMRkS5LARNCkuKiOWNQTwAO6ZlK35xkVmzZyd/eW8Itry4IcnUiIvtHARNiLh6eT3SUi2GFmfTtlkxNfSOPzVzK6wvWU1W7t8ePRERCiwb5Q8zg3uksGvdjEmOj8Xic51LTE2PYVlXHt+t3MLRPZpArFBFpH/VgQlBirJP7g/PS+fXJ/Rh/xXAAFq7V+2JEpOtQDyaERbuj+PXJzutweqUnsEABIyJdiHowXcTg3mkKGBHpUhQwXcTgvHRKyndRurMm2KWIiLSLAqaLGNw7HYAvV5cHuRIRCRffrd/BX975nne+2RiQ42sMpos4rHc6PVLjuffN7xjWJ5OMpNhglyQiXdCG7bt4fcF6psxfx/cbK4iOcnHbKLPvL+wABUwXER/j5umfH87Pnv2cGybM5/kxQ7WMv4i0y7aqWqZ/u4n/zl/H5yvL8Hjg8Px07j37EE4f1JPMAP2FVQHThQzJz+C+cw/l1kmLuHXSQq79UTFx0VEU5SQHuzQRCTHVdQ289fUGpi5czydLS6lv9NAnK5GbTurHuUN6UZCVFPAaFDBdzAVH9mZzRQ0PvWt5fcF6EmLczLv7ZJLi9J9SJNI1NnqYt7q8OVi2VtbSOzOBK0YWcvrAXAb2SsPl6sib6ztGv5W6oF8dX0z/HiksKtnOozOX8sXKrZzQv1uwyxKRIPB4PMxfu42pC9bz9jcb2LSjhrjoKE46uBs/H17A0cVZnRoqLSlguiCXy8VJB3fnmL7ZPP3RcmYtLVXAiESQ+oZGPlqyhW/W7eC9xRv5Zt0OYqOjOMHkcPqgnpzUv1tI3NUIfgXSYfExboYXZvLJsi1U1zWwuqwK00NvvxQJV5t2VPPaV+sY//lq1m3bhcsFpnsK9517KGcf1ovkEAiVlkKrGtlvx/bN5s9vf8/Pnp3NwpLtTLr2aI7UgpgiYaO6roEZizcx6csSPl6yhUYPDCvM5I9nDmBkvxwSYt3BLrFNCpgu7th+2fA2LCzZTlpCDHf99xvevPFYYtyawizSVVXXNTD9u01M/3YjHy/Zwo7qenLT4vnV8X35yRF5FGYHfgaYPyhguriDe6Qy6pAeHNM3i57pCVzx4jyufflLfnVCMUcUqCcj0lV4PB5mLt7MlAXr+GjJFiqq68lJiePHh/Tg7MN6MqI4G3dUcAbrO0oB08VFRbl45tIjmv/8m1MO4u+zVnDBM5t5+6bjNCYjEsK2Vtby4DvfU1ZZy67aBj5ZVkp2chynHtKDc4f04uiiLKK6WKi0pIAJMzec1I+Lh+cz4oH3eWn2Ku47d2CwSxKRFmrqG3h59mq+XredT5eVsX1XLfmZiWzfVc/dpx/MmBF9iA6TW9wKmDCUlRzHWYN78tpX67htVH/SEmKCXZJIxFu2uYLPlpcx/vPVLNm0k7yMBAb0TOX2Uf0Z0DM12OUFhAImTF02og//+bKEFz5dxU0n9wt2OSIRw+PxsGF7NfUNHnpnJvDl6nKe/GAZH9gtAORnJvL8mCM5sX/3IFcaeAqYMHVorzROPrg7j8xYwpad1fzprEO79L1ckVDX0Ohhwhdr+N8ZSyjdWQtAVlIsZZW1ZCbFcuuphrMG9yQvIyFoT9Z3NgVMGHv654dz/1uL+eenqzjp4O6cYPS0v4i/1NQ3MGfFVqJcLgb1TuO68V/y6bIyhhdmctPJPcHj4YtV5RzWO52LhvUmMTbyft1G3nccQWLcUdw+uj+TvyzhjQXrGVGcxcdLSjnB5ITNIKJIMLz2VQnjpn7Ljup6ABJi3NQ1NPLAeQP52dDezT2US4/uE8Qqg08BE+biot2cNjCXNxauJzHOzfjP1/DbHx/E2BM1LiOyPz5esoVxU7+lrrGRtVt3MaxPJtceX8TWyjr+O7+Eq0YWcbzuEuxGARMBzjqsJxPnrmX852tIiYvmsZnLGD0wl2K9R0akTbX1jfx1uuXfc9cytE8mHy/ZQu/MBAblpXPJ8AKuPLaw+U7A+UfkBbna0KSAiQBHFWaRmxaPC3jlqqM4+4lPGDf1W16+YniwSxMJOa/MWc2T7y+jorqeipp6jjsohy9WlnFwbgov/nIY6Yl6XXl7KWAiQFSUixcuH0ZcdBR9spMYe2Jf7n/rexas3UZDYyM7qus1AUAiVl1DI0s37eS7DTuYtXQLry9Yz7A+mZgeKRxvcjjp4O7UNTTiAo1d7icFTIRouWTMxcMLePKD5dw95WuWbtpJXUMjz48ZqvvHEtYWb9jB6rIqeqTFc1jvdL7fuIPHZi7lQ7uFqtoGAOJjorj6uCJ+N6r/but+afHYjlHARKDkuGjGjOjDozOXUpSdRGx0FDf8az7/vX4Efbtp7TIJDztr6pm6YD1HFGTwybJS7n3zu+ZtV40s5NV5JbhccM6QXhxVlMWA3FQKs5O63IKSoUwBE6GuGFlIdV0DPz+qgKgoF2c/8QlXvDiP164bQVVtA73SE/RgpnQp32/cQX2Dh0N7pfHVmnJu/vcCVpdVNW8ffWgPrj+hL09/uJy/z1pJz7R4Jl59NPlZiUGsOrwpYCJUanwMd5x2cPOfn730SC567nOG3jeDRg/8zzmH8vOjCoJYoUj7VNXW89QHy3n6o+U0NHoYXpjJF6u2kpsaz/NjjuTrkh1U1dVz648N0e4oHr9oCCf278ZRxVn0Sk8IdvlhTQEjABxRkMGTlxzOrKVb+GRZKf+as0YBIyGnsdHT3LP2eDz889NVPPHBMrZW1vKTw/PISo7l5dmruezoPvzmxweREh+zx5pfUVEufqJpxZ1CASPNThnQnVMGdOfl2av4/evf8s267STFRXPvm9+xqrSSaTeODOnXs0p4mL+mnLe/2cg1xxVRXd/IgjXb6JOdyD9mreS9xZt4fsxQDs5N5fbJi3hz0QZG9svmppP6Nb8q/PZR/XV7N0QoYGQPZw3uxb3TFvObVxeysrSSqCiormvktfklXDJcvRoJnPLKWq4d/yWbdtQw4Ys1VNU20NDoASA6ykVmUixXvDCXlPgY1m/fxe2j+3PNcUW7LR6pcAkdChjZQ1piDGcMzOW1+es4b0gvbh/dn1++OJfnP1nJRUPz9T+wBMTmimpun/w1WytrefyiIbz9zQZy0xI4bWAPlm+pZHBeOomxbi54ZjZx0VFMuvZovRY8xClgpFX3nH0IvzqhuHna8pXHFvHrfy/goyVbOKG/87zMtqpapn+3iQuOyIuY5cfFf178bBULS7YxZkQfXp23lolfrKW+0cMfzxzAmYN7cubgns37tgySD357PDFulx567AIUMNKqlPgYUuJ/eBPmaQNzeehdy1+nW447KAd3lItnPlrBMx8tp2+3ZA7PzwhitRLq1pRVMX9tOUcUZJCZFMszHy7nsfeXEeWC175aR5QLLh6ezy+PKaRoH2vkaRyw61DASLvERkfxu9H9uXHCfCZ/WcIFR+bxxsL1ALwsTUNgAAAPwElEQVS/eLMCRnazqrSSWctKmb+6nK/WlLOqxfMoMW4XdQ0efnJ4Hreeanhj4XqOLs7i0F5pQaxYAkEBI+125qBcXvh0JQ+++z2pCTGs27aL2OgoZizexJF9Mrh98tf859qjqW/0cPeUr3nw/MF6zqCLW75lJwvXbuPcIb32ehu0vqGRaHcUC9Zu4+kPlzH9u014PJCdHMeQ/HR+flQBR/bJ5NNlpVRU1zOyXzYjirNwuVxcdVxRJ35H0pkUMNJuLpeL+84dyHlPfcb1//qKuOgorvlRMY/NXMptkxaxuaKGv89awfZddXy6rIyHp1v+9tPDgl22HIDbJi3iy9XlfLxkC1eOLKJXegIZSbGsLqtkwhdrcbngs+VlLFy7jYzEGMqr6khLiGHsCX356ZG993g98GG904P43UhnU8DIfjk4N5XHLhrC1S/P48T+3Tn7sJ48NnMpmytqGNgrjX/PXUtdQyNZSbH8d/46rjmueLeFNtuyuaKa7KQ4zVALAVsra9laWcv2XXV8ubqc4YWZTFmwnikL1hMfE8X95w7kkRlLWFe+Cw/Qv0cq1x1fTHllLX27JXPhsHyS4/SrRRQw0gGnDOjOa9eNoCAricykWA7tlUq/bilcf0IxJ//tY2Kjo5hw9VH85KnPuOCZz+iZnkB2chz9e6Rw8fD8PQZxt1fV8aMHP+TKkYX85scmSN9V5CmvrOVv7y3hta9KOGNQT24dZahraOSnz85m/bZq8jMTSUuI4fkxQ1ldVsXa8iqeeH8Zt7y6kLjoKF771TEMzkvTDEJpkwJGOmRIi0H9Kb86BpfLhTvKxVUjC8lKjuOg7ik8c+kRTF2wnrLKWkp31vDi7FX83ycrGdkvm+uOL2ZEcTYAs1eUsauugf+btZLLRvQhOzmu1XNOXbieWLeLUw/poV9q+8nj8VBSvovemc7Cjl+XbOeql+axZWcNx/XLZvJXJfzny7Ukx0Xj8cCPDsrh/e83c/0JxSTFRTOgZyoDeqYyojiLe9/8jlMG9NDtLtknBYwcsJbPI9x1+oDmfz+mbzbH9M1u/vPmimpenbuWV+as4bLnv2DydSMYlJfOZ8tLiYuOoqa+gUfeW8LvRvcnOTYal4vmICnbWcNvXl1AXYOHkw/uzhMXDyE+pmtOV21o9LB2axV9spMO+FiNjR7Kq2qJi3G3eluquq6B+Bg3L3++mj+8/i1nDMqlV3oCL85eRVZSHK9ffwyH9kpj2eadTF2wDrupgquPK+bw/HTmrS5ncN7uIZISH8OD5w8+4LolMihgpNN0S4ln7In9+PlRBZz+2Cdc/6+vePOGkXy2vIyjirLonhrHK3PW8MqcNQDEuqMoyknixpP6UVJeRV2DhzEj+vDCZ6uY/FXoLFtTVVtPbX1ju1+l+8ep3zD+8zVMvPoojirKanO/iuo6/vTGd9Q1NPLXCwYzce5a1myt4k7vKthfrNzK2H99xeaKGlLio5nsfdXCe99t5IYT+/HS7FU8PH0JD10wmEfeW0J+ZiLvfruRhkYPpw/qyR/OGEBOitNb7NstmVt8bk8O7aOn5OXAKGCk06UnxvL4xUP46TOzufblL1m2eScXHJHH5ccUcmL/bqwuq6K6rpHK2no+slu4+d8LyEyK5YiCDP545gC+XF3OP2at5Lh+OYyfs5prjysmIymWHdV1jP98NSnxMVy6nytBz1lRRmFOEt1S4lvd3nIV3yYNjR6e/GAZ//hkJdFRLt688Vh6pMZT3+jZ4w2Imyuqmbl4M5U19Yz/fA3uKBfjpn7Lc5ceyewVpZx3eN5uX7N5RzUXPDubtVuraPTAum27mLuqHIDjD8qhvKqOX/97Pr0zErnmR8U8/eFyfvnCXMp21rKrroG5K8v5ck05bpeLGyfMx+WC8VcOJzPJCcHcNE0fl8BTwEhQHJ6fwe2j+/M/0xYDzu202OgoRh2au9t+Vx9XxGmPzmLD9mpuOeWg5ucmbpwwn9GPzmJnTT2LN1RwyfB8bv3PQnZU1wPQr1vyXnsH4ATE9l11PPvxcp79aAUnH9yd/7vsSOCHMYvuqfE889Fynvt4BfeecwjnDvlhmfcJX6zhb+8t4cT+3ZizooyrX/oSDx7Kdtby+thjmsNqa2UtFz33Ocu3VAIwKC+NK44t5KaJC/jRXz/A44FZS0t59MIhzW9TvP+txWzYVs3Eq49m5uJNPPvxCkYUZ7GytJK7p3xDybZdDM5L5x+XDSUtMYZDe6Zy8f/NoW9OMqcNzOWRGUvolZ7A+CuH86tXvmJ4YSaH9NSDjNK5FDASNFccW8gXK7eyqGQ7A3JTW90nOzmOZy49gglz1nDGIGdtqtMO7cFDmQnsqm3kkqPyefajFXy8ZAsDe6XxxzMH8Jv/LOS2SYt4+6aRJMVFs3F7NU99uIzvN1ZwdFEWFw7rTWVNA5c9/wXrtu0CIC8jgQ/tZkp31pCdHMer89byu8lfEx3lor7RQ/fUOG55dSGfLC0jKzmW4/rl8NfplqOKMvnHZUcy7esNjP3XfHLT4imvquXmfy/gpV8OZ2HJNu7+7zeUlO/i+TFHEut2c3BuCplJsXy+ogyPB3JS4nj8/WUsLNlGbmoChxdkMGXBeq4/oZhhhZkcWZDBEQUZHNM3mzcWruf2174mPzORv//iSNISneV8hhdlMf3m4+iRGk9SXDT9c1Po1y2Zwuwk3rrx2M75Dyriw+XxeIJdQ9AZY/oAK2fOnElenl5E1JkaGj1U1taT2mLds/bYXFFNrDuKtIQY/vKOpaK6jt+fMYD4GDdzVpRx4d8/55CeqZzUvzvPfryc+gYPB3VPYfHGHSTEuImPcRPlcnHd8cX07ZZMz7R4TnnkY35/xgAuGtabHz30Id1S4prf1X76oFx++5+FfLKslKqaBmobGnFHuZh247H07+GE47frt1OYncTUBU4INMlOjuWh8wc3LxLamlfmrObTZaWsLqvi2/U7yE2LZ+ZvfkRi7O5/B6xvaOTvs1Yy6tAeFPphkoDIgSgpKeGkk04CKLTWrvLdroBBAROOZi7exE0TF7Czpp5Rh/TgztMOJj8rkbVbq7jnje9YvGEHL1w+lH7df3gI9OwnPqGytoGhfTKY8MVaJl/X+nLwFdV1vDqvhLSEGM5v5c2IHo+HNxZtYMWWnaQlxPDTI3uTtB8PHi7esIPkuOjmKcUioUoBAxhjbgGuADzAPOAaa21Ni+19UMCEnTVlVWzYvovhrYzFeDyePZ6lmfDFGu7w9jzOGtyTxy4a0il1inRV+wqYsB+DMcYMBy4HhgFVwEvAjcBDwaxLAi8/K5H8rNZ7Aa09qHnh0N4cVZSFx+OhT5ZuP4kcqLAPGKAcGGutrQQwxiwC8oNbkoQil8ulcQ0RPwr7gLHWLgGWABhjcoGbgDHBrElEJBKETcAYYy4EnvBp3mStPcS7vQB4G3jaWjujs+sTEYk0YRMw1tqJwMTWthljDgPeAB601j7eqYWJiESosAmYthhjcoB3cMZhJgW7HhGRSBFyAWOMOQt4xVqb4tN+FXAbkAcsAG6x1s5uxyFvBlKBu40xd3vbpllr7/Jj2SIi4iOkAsYYMwIYD7h82n8BPAP8CZgL3AC8a4wZbK1dubdjWmvvBO4MTMUiItKWkAgYY0wczuyue4FKILbFNhdOsDxnrb3H2/YeYHF6Jzf6oQQ3wMaNG/1wKBGRyNDid2arL2cKiYABRgN3ALcCWcBvWmzrCxQAU5sarLV1xphpwCg/nT8X4JJLLvHT4UREIkousNy3MVQCZi7OUgPbjDHjfLYd5P1c5tO+Aig2xrittQ1+OP9IYANwoMcSEYkUbpxwmdvaxpAIGGvtur1sblrHvcKnvQKIApKAHQd4/hrgkwM5hohIhNqj59Ikqq0NIaRpwN93Vc6m9sZOrEVERNopJHow+7Dd+5kCbGrRnowTLpWBPPn+To82xhwKPAoMB7YCT+I84Bn+y1bToZ/Xm8DprWxKsdbuDEyVoamtKfqt7BfR11iT/fh5Rew1Zoxx40ygugpnDcbVwFPAk21dL/68vrpCD2ap97PIp70IsIH8n6rF9OjxwE+AbTjTowvb2L8bMAOnt/VT4DngPnaftBC29vfn5TUI52I+2uefqsBWG1ramqLfyn4RfY01ae/PyyuSr7HfA/fj/KzOAl4F/hdnQtUe/H19dYUezFJgLXAOMB3AGBOD8zeSaYE6aQenR1+P8zM9y1pbBbzlnYJ9hzHmUWttXaDqDbaO/LyMMelAb+Ada+3nnVhuyNjbFP02ROw1Bvv/84rka8wYEwXcAjxkrb3P2zzTu7rJb4EHW/kyv15fId+D8fZQHgCuNcbcZ4w5DXgdyAYeCeCpW50ejRNqbU2PPhmY6f0P02QKkAkMDVCdoaIjP69B3s9FgS0tpLWcot+edfIi+RqD/f95RfI1lobz/qvXfNotkGOMae3dFH69vkI+YACstU/hXFCXApOAdOBUa+2KAJ52n9Oj2/ia1vZvebxw1ZGf1yCgBvgfY0yZMabKGPMfY0yPQBYaYpqm6D/GnhNZWhPJ1xjs/88rYq8xa225tXastXa+z6YzgZKmd2T58Ov1FXK3yKy144BxrbQ/DDzciaV0ZHp0ahv7tzxeuOrIz2sQEOfd51yccbX/Ad43xgxp+VrrcLWPKfqtieRrrCM/r4i/xloyxlyJ00tpawUUv15fIRcwIaQj06NdrezfJNynU3fk5/U3YIK19gPvnz82xiwGPscZYHzZ71V2fZF8jXWErjEvY8wlOJNwJrHnu7Oa+PX66hK3yIKk5fTolvY2PXp7K/untNgWzvb752Wt/b7F//hNbXNwZp8NDkSRYSCSr7H9pmvMYYy5GSdM3wQu2cvsW79eXwqYtnVkevTSNvYHZ2AtnO33z8sYc6Ex5jifNhfOLY3SgFTZ9UXyNbbfdI2BMeZ+nJ7cy8D51travezu1+tLAdO2ltOjgd2mR89s42tmAif7zM44ByjDeegwnHXk53Ud8Kh3OmWT04AE4OMA1dnVRfI11hERfY0ZY27CmXX3KDDGWlu/jy/x6/WlMZg2WGs9xpgHgCeMMeXAp8BYWkyPNsYUAzkt5tc/hfOumreMMQ/hdMHvAG7fx98aurwO/rzuB94Gxhtj/okzS+VeYLK19rPO/h5Cka6x/aNr7AfGmFzgL8DXOK+TH26MabnLPJxHCwJ2fakHsxftmB79e2B2i/034MzQiPbufzVwl7X2r51Zd7B04Of1Ls7TxX1x5trfBTzv/Xpx6BrbP7rGfnAqzq3AgTg/E99/0gnw9eXyeCJq+SIREekk6sGIiEhAKGBERCQgFDAiIhIQChgREQkIBYyIiASEAkZERAJCASMiIgGhgBERkYBQwIiISED8P7wgZhMDA9CBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(delta_t_list,result_Bayes_array.T[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us first calculate the standard deviation for each point of the acf\n",
    "# all acfs are stored in acf_all for later analysis\n",
    "\n",
    "result_Bayes_array2 = None\n",
    "for delta_t in delta_t_list:\n",
    "    data = langevin.time_series(A=A, D=D, delta_t=delta_t, N=N)\n",
    "    \n",
    "    results = OUresult2(data,delta_t)\n",
    "    if result_Bayes_array2 is None:\n",
    "        result_Bayes_array2 = results\n",
    "    else:\n",
    "        result_Bayes_array2 = np.vstack((result_Bayes_array2, results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11b990860>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD6CAYAAACoCZCsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4HNW9xvHvqq16s4olS7LkNi64N2wDhpheHVoAB0ggQCAEknDDJeGSSwgBUiChcyFAiCGhkwAGE2KqDW4Yd3tcZUm2ZMlW72X3/nEkW5YLlryrlTTv53n8LJpdaX9aRvPOOXPOGZfX60VERJwnKNAFiIhIYCgAREQcSgEgIuJQCgAREYdSAIiIOFRIoAs4WpZluYHJQCHQEuByRER6i2AgDVhm23ZD+yd6TQBgDv6fB7oIEZFe6kRgYfsNvSkACgFeeukl+vfvH+haRER6haKiIubMmQOtx9D2elMAtAD079+fjIyMQNciItLbHNR1rovAIiIOpQAQEXEoBYCIiEMpAEREHEoBICLiUAoAERGHckQA/H1JHlc+uyTQZYiI9CiOCIC80lqWbC8NdBkiIj2KIwIgNiKExmYP9U1aQkhEpI0zAiA8FIDK+qYAVyIi0nM4IwAiWgOgTgEgItLGGQEQbpY8qqhrDnAlIiI9hzMCIEJdQCIiHTkjAMLVBSQi0pEzAiDCdAFV1qsLSESkjTMCQC0AEZGDOCIAwkODCQsJ0jUAEZF2HBEAYFoBlRoFJCKyj3MCICJELQARkXacEwDhoboGICLSjnMCICJUo4BERNpxTgCEh1ClFoCIyD7OCYCIUF0DEBFpxzkB0DoKyOv1BroUEZEewTkBEBFCY4uHhmZPoEsREekRnBMAmg0sInIA5wSAVgQVETmAcwJA9wQQETmAYwIgTi0AEZEDOCYAdFtIEZEDOScA9t0YXl1AIiLgoACIjwwlLDiIvL01gS5FRKRHcEwAhAYHMTojjq92lAW6FBGRHsExAQAwcWACa3dWUt/UEuhSREQCznEB0NjiYe3OikCXIiIScI4LAEDdQCIiOCwAkqLdZPeLZLkCQETEWQEAMHFgIisUACIizguAoanR7K1ppLpB8wFExNkcFwApMW4AiivrA1yJiEhgOS4AUmPDASiuaghwJSIigeW4AGhrAexWC0BEHM6BAWBaACVqAYiIwzkuAGIjQnCHBKkLSEQcz3EB4HK5SIl1qwtIRBzPcQEAphuouFItABFxNkcGQGqsm+IqtQBExNkcGQBqAYiIODQAkmPcVDU0U9eoZaFFxLkcGQD7J4OpG0hEnMuRAbB/Mpi6gUTEuZwZALGt6wGpBSAiDubIAEhtnQ1cWK4AEBHncmQAxEeGkpEQwR/+bfPkJ1sDXY6ISEA4MgBcLhdv3jidmcOS+d38jXydpxvEiIjz9LgAsCwrzrKstZZlDfHn+6TEhvOn74wjxh3Cswu3+/OtRER6pB4VAJZlzQAWAcO64/2i3SFcNiWT99cWsbO8rjveUkSkx+hRAQDc0PpvV3e94dXTswF4eWled72liEiP0KMCwLbtq2zbXtSd75mREMm4zHgWbdnTnW8rIhJwId39hpZlXQY81mHzbtu2R3V3LW0mZyfy7MJt1De1EB4aHKgyRES6VbcHgG3bLwMvd/f7Hsnk7ASe+tTL13nlTBvcL9DliIh0ix7VBRQokwYm4nLB8tzSQJciItJtutwCsCzrfOAl27ZjOmy/DrgdyABWAj+zbfvLY6rSz+IiQ7FSY1iqABARB+lSC8CyrOnAi4Crw/argKdan7sIKAc+sCwrpzM/37btbNu2t3Sltq6alJ3AVzvK+Mvn2ygoq+3OtxYRCYhOBYBlWW7Lsm4HPgaaOzznAu4BnrZt+9e2bb8HnA/sAX7qo3r95uzj0gh2ubh33gbufXdDoMsREfG7zrYAzgJ+AfwceLTDc0OAgcDbbRts224C5gFnHkON3WL6kCRW3306s8els3xHKV6vN9AliYj4VWcDYBmQY9v2I0DHI2Tb7N2OXTfbgMGWZfX48ZUul4tJ2YnsqW4kr1TdQCLSt3XqIrBt2zuP8HRs62NVh+1VmKCJAio7836BMCk7AYCvdpQxsF9UgKsREfEfXw4Dbbsg3LFl0Lbd48P38puhKTHEuENYvkMrhIpI3+bLAKhofYzpsD0ac/Cv8eF7+U1wkIvxAxNYoQAQkT7OlwGwufVxUIftgwDbtu1ec1V1YlYC9u4qymoaA12KiIjf+DoA8oHZbRssywoFzgEW+PB9/O7UkSm4gN/MWx/oUkRE/MZnAdB6hv8A8EPLsn5rWdbZwL+AJOBPvnqf7jAqPY6bvzWUN1fsZN7qwkCXIyLiFz5dC8i27ScwcwSuBF4H4oEzbNve5sv36Q63fGsIVmoMf1nY60oXETkqXV4LyLbtu4G7D7H9QeDBrpfUM4QEBzFrRApPf7aN2sZmIsO6feFUERG/0mqgRzAlJ5Fmj5evNCJIRPogBcARTMpOJDjIxZJtWiVURPoeBcARRLtDOC49liXb9wa6FBERn1MAfIOpg/qxKr+CP35gM39tUaDLERHxGQXANzhpaDKNLR4e+3gLP3zxK15emhfokkREfEIB8A1OGJrER7fNZNX/ns7MYcn84q01rMovD3RZIiLHTAFwFAYlRxMXEcqjV4wn2OXifXUFiUgfoADohNjwUCZlJ/CJXRzoUkREjpkCoJNOtlLYWFRFYUVdoEsRETkmCoBOOtlKBuBTuyTAlYiIHBsFQCdZqTGkxYXznw3qBhKR3k0B0Ekul4vZ4wfwnw27WbxNE8REpPdSAHTBj781hMzECO54YzU7y3UtQER6JwVAF0SGhfDAhWPIK61lxgMfMfvxRZobICK9jjMCYM3r8MYPfPojZwxJ4sOfzeTOs0ewq7yO2U8sYsGG3T59DxERf3LGIvdFq2H923ChF1wun/3YwcnRDE6O5rIpmZz9yOc8/dk2Zo1I9dnPFxFnqmtsYe2uClbll/N1fjlFFfU8MWcCqbHhPn0fZwRARAK0NEBTHYRF+vzHx4SHcvmULH4/32ZLcTVDUqJ9/h4i0jc1t3jYWlLDqoJyVuaXsyq/nI1FVbR4vAAMiI9gak4iEWHBPn9vZwRAeLx5rC/3SwAAXDIxkz99uIm/L8njV+eN9Mt7iEjvt6u8jmW5pSzdXsqanRVsLKqisdkDQEx4COMy47lx5mDGZcYzNjOe5Bi332pxRgBEtAZAXRnEpvvlLZJj3Jx5XBp/X7qDydkJrC+sZNueGh69bDxBQb7rdhKR3qOitom1uypYsm0vX+eXs6W4msKKegBi3CGMyYzj6mkDGZUex3ED4hiUFNWtxwuHBECCeazz70idX507ktw9Ndz40op92+ZMyWL6kCS/vq+IBF5xZT1rd1Wwbmcl63ZVsnZXBQVlZph4kAtGpMUyNSeRMRnxTMlJZERaLMEBPjl0RgCEt2sB+FFyjJuXrz+eRz/awoSseG57bRWvryigqLKeFxfvYO61U4lyh1BV30RMeKhfaxER/6iqb2JbSQ1bS6rZUlzNul3mgL+numHfa3KSohibGc8VU7MYlR7H+Kx4Ynvg37wzAqCtBVDv/7H6Ue4Q7jhrOADn2iW89XUB760ppL7Jw2vL82n2eLn//Y38dvZxXDYly+/1iMix+TqvjDdX7GRzcRXbSmoortp/oA8JcjEkJZqTrWRGpccyKj2OEWkxveYEzyEB0D0tgI4unjiAfyzNIz4ylJykaJ75fDvltY2EBQdxx5tr+HRTCScOTea8sWm9ZocR6au8Xi+7KxtYs7OCtTsrsIuqyN1bw8aiKqLCgrH6xzBzWDKDkqMZnBzFoORoshIjCQvpvdOpnBEA7lhwBfv9GkBHE7IS+O7xWZw2sj+1Dc3c+NIKwoKDeO/WE3hlWT7vrCrk/bVF3PfeBn513kgunZTZrfWJOJHX66WpxUtji4d3Vu1i0ZY9lNY0sml3FXuqGwHTZ5+dFEVWYiQXT8zgsilZRLv73uGy7/1Gh+JyQXhct3QBHfi2Lu6dPRqAFo+XCVnxfGt4CkNSYrjznJH88uwRrNlZwV3/Wsf9723g/LHphIf6fqyviNMVV9Xz2vIC/rVyJzv21tLQOuwSzDj71Fg3J1spHJcey+iMOEakxRIZ1vcPj33/N2wTkdDtXUDtBQe5ePOmGQdsc7lcjMmI5/YzLOb8ZQnvrSnkwgkZAapQpHcrrWkk2h3CRxuLeeGLXHaW1xES5CIl1s3y3DKaPV6m5CRyspVCbHgIXi9MG9yPiQMTcPlwhYDexEEBEN/tXUBHa/rgfgxKjuJvX+5gaEoMW0qq2FvdyLfHD6BftP8mgYj0VvVNLdhFVWwurqa8tpEP1hWxLHf/CV52v0jGZ8XT0OQhv6yW78/I5rIpWQxO1iz99hwUAAlQWxroKg7J5XIxZ+pAfvPues57bOG+7X/5fDuPz5nAxIEJLN62F7uoiqumDXTs2Yo4z/pdlbyyLI/VOytIiwunqcVLfmktW4qraW5dKgEgMzGC204bRrPHS1ZiJBeMSyckuPdenO0uzgmA8Hgo3RboKg5rztQs01yNcTMkJZraxhZueflrrn1hGfNvPYlb/vE1xVUNrNtVwX3fHq2dW/qU2sZmXl2WT+7eWk4flUq0O4TXvypg7uIdhAYHMS4zno1FVYQFB5EWF86sESkclx6H1T+G+Mgw4iNCNeO+C5wTAAG+BvBNwkODuXp69gHbnvruRM555HMufGIRxVUNnD82nVeXF1BZ18zDl4/DHaILxtI7eDxemjwe3CHB7Nhbw57qBoYkx/Dwgs0s2Lib3ZX11Dd5CAsO4q9f5AJm7Mb3pmfzk1nDiIvUMGl/cFAAxEN9BXg8ENQ7zp5HpMVy5fEDeeHLHcwansIjl49nXGY897y7nguf+IJTR6RyzYwc/XFIj+D1enlteQFJMWHMHJbCmysKaGzxMCQ5mtteW0VZTSMTsxNZuLkET7uV2U8dkcqs4amcPbo/I9Nj+XzzHlzA0NQYcpKiAvo79XXOCYDwePB6oKFy/8SwXuBnp1k0ebxcf+IgAK45IYd+0WH836fbeOSjzeSX1fLQpeMAKKyoY8WOcs4e3V/XCcQviqvq2VBYxUlDk1q/bmBPdQOlNY28vCyfeasLAUiPC2dX66JnbV+fMao/i7bu4erp2UwcmMDKvHLOHZvOuMwD/x7PGNW/+34hh3NOALRfDqIXBUBcZCj3fXv0AdsuGDeAC8YN4N531/P8F7ncOmso7pBgLv2/L8kvrWPutVOYPjiJwoo6MhL8s/y19G0VdU08v2g7zy/KpaahmQEJEZw/Np2XluRRWtPI6AFx1DQ2s62kZt/3BLngv88cjsfr5Y2vCnjo0rHkJEWxLLeUiyZkHDSi7dwx/lmZV46egwKg3XIQCdkBLcVXrj9pEH9bvIOfvrKS3ZUNVNQ1kRLj5o8f2GQm5jNvTSGv/3A6EwcmHPL7vV4v2/bUMCgpSi0GB8rdU8M/luWxp6qRQclRDE6OpraxmY/tEv69roiGZg+njUxlaEo0S7eX8uhHWxjeP4Yff2sIc7/cQVK0myuPH0haXASJUWFkJkaQFhcBwI9OGbLvfcZnHXr/k8BzUAB0z5LQ3SklNpw5U7N4flEuEwcm8Mjl49laXM3tb6xmVUEF4aFB3PPOOv5y9WTmLt7BvNW72F3ZQEx4CLedbrEyv4wXF+fx4CVjuWiiJqD1RjUNzfxu/kaGpcZw8cSMQ84kX7ergucX5TIsNZrrThzExqIqHvtoC++tLSTY5SIxKow3VhTse31iVBiXTsrkO5MzOW5A3L7t+aW1pMS6cYcE8/0ZOd3y+4l/OScA2t8VrA/55dkjuOGkwfSPM/cKHZsRx7w1hYzLjCcrMZLbXlvFjAc+osnjYdqgfswclsLK/DL+67VVAES7Q3j84y2cMyaNeasLOWV4ColRYQe9T2V9EwWldYxMj+3W38+J6ptavnFJkPzSWkprGrnn3fV8tcOMbnvw3zbnjkmntLaRjYWVRIQFU1Rh+uhDg100tXh5ZVk+W0tqiHaHcOPMwVw9PZvU2HCqG5rZXlJDcJCL4f1jDjmkMjNR3Yl9jXMCYF8LoOcOBe2K0OCgfQd/gJDgIF64Zgpght59sK6IIJeLn59p7ZsF2eLx8vyi7bhDgkiMcvOjv6/gzD9/Ru7eWtLjwrnr3JGkxIYT7Q4hMiyYstpGbv771+SV1nLXuSO59gSd/XXW1pJq0uLCv3F9mTe+KuD2N1Zz5fEDuXxKFhsKK3lvTSFbiqtp8XqZkJVAVX0T/9lQDJjliJ+cM4GEqDDmLt7BK8vyiYsMZUJWPI3NHob3j2X0gDhmjx/A618V8NzC7dw6a+hBo8ei3SGMzog7XFnSR7m8Xu83v6oHsCwrG9i+YMECMjK60F3R3AD3psLJd5h/ApiQOP3Pn5G7p4ZbZg3l1eX5++5i1F5yjJuRabF8uqkEKzWG1Lhw8vbWMCAhgrNHp3HZ5KxD3t3I0zpbs+2MsqCslkue+pJR6XFcc0I20wf3rLulrcgr460VO7nznBFdWpivtrGZmoYWYiNCcIcEU93QzP3vbeClJXkMSYnmmasmHTS00ev1Ut3QzKr8Cq55YRkpMe4D/h+kxLiZnJOI1+vli6178Xrh2hNyGJEWS05SJENSYva9tqG5hdCgIE2Kkn0KCgqYNWsWQI5t27ntn3NOCyDEbe4HXLYj0JX0KEFBLp7/3mSqG5oZkRbLD07MYd2uSmoazIGsprGZhmYPp49MJSnazbMLt/HF1r2UVDUwIi2WzcXV3PnWWhZu3sMvzjKrmz728RZiwkO4fEomj360hWCXi0cuH8+ItFheWpLH7sp6mlo8XPHMbq6YmsUlEzOIdodQUF7HqLRYUmLDD6ixqKKe+MjQwx6Qj6bL5Gisyi/nqmeXUt3QTE5SFNeckMPe6gYSIsMICnLR0NzCYx9toanFy/UnDeKVZfnkldZy4tAkTh2RyvIdpVz71+XUNbXQLyqMX503kic/2cqm3VV8Z1Im/15fxHmPLuQnpw5lXGY8K/PLeWfVLjbtrqauqQUwK1O+ffMJ7NhbQ+7eGoYkxzAyff+tA1s8Xrxe72FngmtyoHSGc1oAAM+fbeYCXDPfl6U53l8+38a98zbs+3pwchTltU3srWkkIyGChmYPFXVNPH7FBO54YzUTBibw6OXj+dOHm3j682203wWDXDA1px9jMuKIjwxjfWEl767exeSBicz9wRS2FtdQ3dDMsNRoXC4Xf/hgIy8uzmNAfAQzrWROG5nKl1v30uLxcvrIVJ78dCv5pbVcPiWLwop68ktryU6KIrtfFMPTYhibEU9tYzNzF+/g0QVbSIoJIzHKza7yOs4YlcqLi80NfYalxLCnuoFte2r21enxQlRYMDWNLeQkRVFS1UBaXDhXThvIS4vzsHebG4k8+d2JnDQsmYKyWu58ay2fbirZ9/uOyYhjcnYiKTFuUmLdzBiSRErMgQEociyO1AJwVgD88ybY+hHcttGXpQmwLLeU7SU1JMe4OXFoEjUNLXy4YTdnHdefuqYWrn5uKesLK/F64YVrpjBzWDJgLmZuLKqitrGZlJhwFm4p4RO7hE27q2hq8RIZFsxpI1P518pdDEqK2ncAbu/SSRlUNzTz0cZi6ps8hAS5cLmgqcVLtDuEQclRrC6oICwkiMyECPLL6mhsXQ8+PjKU6vpmmlsD49cXjCJ3Ty2XP7MYgIsmZBAcxL415G+dNZTYiFBe+CKXSydlcvygRD6xS3hg/kbqGlt47YfTSI+PoLaxmecWbueU4SmMSt/ft+71elmyvZS6xhaGpETrwqr4nQKgzad/gI/vhTuLIDTCl+XJNyipauDip74gOMjFf3468xv7qJtbPDR7vAQHuQgNDuLxj7fw8ILNXH/iICYMjGdbSQ2NLR6mZCcyKTsRMJOXlm4vZUJWPM0eLws2FHOylUxaXHjrRdgIotwhtHi8ZtZ0Xjmf2iWkxLo5bWQqE9qNV7/77XWkxLq5cebgo5oj4fF4afZ4e/XtAaVvUgC0Wf0avPkDuGkJpAz3ZXlyFGobm2lo8pBwiGGmR6O5xaNVUEU66UgB4Ky/prYZwGW5gazCsSLDQrp88Ad08BfxMWf9RSW2jl8v2x7YOkREegBnBUBkPwiLVgtARASnBYDLZbqBStUCEBFxVgCACQC1AEREHBwAvWT0k4iIvzgvAGIHQEtDn1sUTkSks5wXADGp5rGqKLB1iIgEmAMDIM08VhUGtg4RkQBzXgBEt7YAqncHtg4RkQBzXgDE9DePagGIiMM5LwDCosAdC1VqAYiIszkvAMC0AtQCEBGHc24A6BqAiDicMwMgWi0AERFnBkBMf3MNQLOBRcTBnBsAmg0sIg7n3AAAXQcQEUdzZgBEay6AiIgzA2DfZDCtByQizuXwAFALQEScy5kBEBYFsRlQtDbQlYiIBIwzAwAg63jI+1JDQUXEsZwdAFWFUJ4X6EpERALCwQEwzTzmLQ5sHSIiAeLcAEgZAe440w0kIuJAzg2AoGDInKIWgIg4lnMDAGDgNCjZoHsDiIgjOTsAhp1lHu15ga1DRCQAnB0AKSMgcTCsfzvQlYiIdDtnB4DLBSPPh9zPobY00NWIiHQrZwcAwIjzwNMMm+YHuhIRkW6lAEifYFYH3fpRoCsREelWCgCXCzImQcHyQFciItKtFABgAqBsO9TsDXQlIiLdRgEAMGCSedz5VWDrEBHpRiGBLqBHSB8PriDYuRx2rYCkYXDchYGuSkTErxQAAO5oSB4By56F2j0QHAaJgyB9XKArExHxG3UBtRkwwRz8+4+GqGR4/fvQVBfoqkRE/EYB0GbwKRASDhc8Aec8BKXbYNsnga5KRMRv1AXUZtSFMOxMc7vIZAtCI2HLArDOCnRlIiJ+oRZAG5fLHPwBQtyQcxJsXRDYmkRE/EgBcDiDZ5luoNJtga5ERMQvFACHM2SWedyiVoCI9E0KgMNJHAQJOfD1XGhpCnQ1IiI+pwA4HJcLTr0bClfBR/cGuhoREZ9TABzJqNkw8Xuw6M9aLVRE+hwFwDc5435IHg5v/RCqSwJdjYiIzygAvklYJFz8HNRXwAe/CHQ1IiI+owA4GqmjYNwVYM/XBWER6TMUAEdr0CnQWKUlo0Wkz1AAHK2cE82S0Vs/Nl+3NMHnD+pm8iLSaykAjlZEgrlvQNsCcZs+gAX3mHkCIiK9kAKgMwadDAXLoL4S1v/LbNNMYRHppRQAnTH4W+Btga/+Cvb74AqGvC+hsSbQlYmIdJoCoDMGzjAh8OFd5oLw8TdCSyPkLoKKAvB6A12hiMhRUwB0hssF5z0CYTEQHgcn3wEhEfD+z+FPo2D5s4GuUETkqCkAOis+E74zFy54HNwxkH0ClOWaUFjxt0BXJyJy1HRHsK4YfMr+/z7nQdP9U7Qa5t8BxRshZXjgahMROUpqARyrhIGQPQOOu9hcFF79cqArEpG+xk/XF9UC8JXoZBhyqhkhNOgUGDQz0BWJSG/UVA+715pVB3augF0roLIQbvgU+g326VspAHzpjN/CPy6Hv10Ax98Ep/wSgkIgNDzQlYlIT+RpgT2bYefy/Qf83evA07rmWHQqDJgI46+E+Cyfv70CwJeShpqU/vddsPhxWPwE4IVZv4ITbwt0dSISSB4P7N0ChSvNjaYKV8GulWZIOYA7FtLHwfSbIX2COfDHppvRh36iAPC1sCg49yEYc6lZLmL7Z/DFozDlBnBHB7o6EekOLc3mQF+02nTf7Pq6dRWBcvN8sNusMjz2O+ZAP2Ai9BsKQd17WVYB4C9Zx5t/+Uvh2dPgs99DVREMnG7uMiYifUfNHihYbg7yBctMV07bmT0uc1OpkedDxhRzsE8aBsGBP/wGvoK+LnMKZE2HRQ+br1e/AuV5MPkHpnknIr1LWS4UroaiNeZibdFaqMgzz7mCzZn9mEvMHKGMKRDTH4JDA1ry4SgAusOZ95tZwjN+Ap/9wSwj/fmDJgTOeTDQ1YnIoTTVQ10p2O+ZZeDdsabfvnided4VBP2GQMYkmHwtZEw2ffhhUYGtuxMUAN0hfRyc/6j579lPwtQbYPnzsOwvkHm8OVsQkcCqKzMjcfZuhfVvw46F+59LyDHrfsUOgDMfMC375BHmlrG9mAKgu7lc5r4C54yGEhve/YlpUk75gbnngIj4T1MdbP4QCpaafvvoFHNBdtvH5uDv9ZjXJQ6CE//LdN8MmGj+Zv04GidQFACBEhxibjb/zi3w8b2w6X249kPTrIQ+ubOJdKvijebAXrrN9MtXl8CXj0J9hTnoRyVBdbFZ4j19Apz0c9Nv328IxKQ54m9QARBIcQPgu2/Ayn/AP39org9smg+44IpXzNmJiBxZZaGZOBUWDQsfMhOpGqrNWT6YFXub68x/W+fAlOsg+0RzEtbSBM31ZmFHB1IA9ARjL4OVL8En90NoFOCFZ0+H780zISEi+9VXmpV3y7ZDeT5s+dB03QSFmMf+Y0xLetb/mr+tmDTT3QoHL9QYHNpjR+h0BwVAT+Bywbl/Njeamfnf5qxk7rfh5Svg1LthyVMw4SoYfk6gKxXxr4YqaG6AyH7m76J0O+QvgT2boHIXVO40I3HqKyAi0UyuPOGnEN3fdPVMuNJ093SkFXoPSQHQUyQNgcv/sf/ri54x6wrNnW3ObDbNN+uBnHm/Y5ur0sfUlsKa16Gq0By0lz4D+YvNc+kTIHUkrPy7Oat3BZsz+dh0sM6GKdfDgAmBrb8PUAD0VNZZZuhoWS5M+5FZTmLRn2Hbp6ZbKDoVLnq2R8wmFDmi5kbTH7/6FdiywIyXDwqGDe9CS4PprvF6zD59yv+YfXrZc+ba2JTrYcLVZp0tB3fV+IuOHj3ZhCv3//ep/wtDT4f/3G3GI6//p5lfMOMnZs2Rje+Zqeb9RwesXHGQ/GUQHmuWNNi71XTXeL1mv+w3GEbONjdKWvwkfD0XGiohNBIGnQy5C81F24lXm1Zt4iCzZk7/MfvXy5p2MzTVmluvit8oAHqTgdPg2g/MH9qrV8HH98GKuVC61Ty/8iW44TMzvE0//is8AAAMe0lEQVTEH5ob4cNfwZInzdfh8fsXOGsvaZhZ5tgVBMddCCMvgJyZJjQ8LeY1QcH7Xz9w+oHfHxwKwTr4+5sCoDdyucwSEsUbIDYNZtxizqJevBjeuBbmvL6/ubxxnllHXC0D6ayitWYAQmIORCaZs/u8xebMfOoPzX61e53p0gkOM9uHnWmWTlj9Kpz0X2bwQsd17Nsf+CWgFAC9VXQK/Hj5gdvOeRDevtmEwEXPQfF6eOW7EJcJP/4KcJkzsm5eclZ6oIZqM/O1bLtZnLB0uzmh8LZAykgzQSp/sZkw1TaGPiHHHNCHnmbufnc4U28w/6THUwD0JROuNH2tH/wS6i6ExmpzZla+Az79Pax704yiuPwVeO175o/9Oy9CaESgK5euqq+AN683B+0p15nlRNr+f+7daoZONjeY60QtrXeZyl1ohla23XXKFQzxmWbJYlzmtTFp5vrSjFugscaM2Ok/2hGzY51EAdDXTPuRGSb6/h3QVAPnP2YmzXz2e3MRbu8WeGKqOesDeO37Zh2isGiz2uHhRltUF5v1zoef3b2/j9N5WiDvS3PhP/tE8/+mrtyMqOk/Gr54DDb/29x8aOFD5ntyTjIXVBc/sX9tGwBcgNd837SbzOuShkFM+pFHk0UkQFyGP39LCRAFQF804SpzsMhdCOPmmDO7+f8NZ/8R1rxmDgwn/MwsdPX+7WYdojahkaZPN/sEmHTN/gvK7/wE7Hlwxv3m4HEoNXthxQswZFbrbEydLXbKni3gaTaTlurKYflzsOT/oLrIPB+RCGljzRr0NSX7v+/M35nRNds+htq9ZqXZ7Z/B2CvMTNigEOh/nAn55oZev4Kl+I7L6/UGuoajYllWNrB9wYIFZGTobKTLPC2m73fAJHMtoHS7ObtvrILaMnM3o7wvzAXAmDS45K9mVuZjkyAy0XQFDD3dnBW6gswwvX6DzVjtV6/aHyajL4ELn+k5IVBXZkasHE09DVXwwZ2QMgKOv/HIry2xzRn48HPMbNTVL5v1ZmJSzWe97i2z+mRjdetomG9D7mfmukxcpmmZRSWbA/vcC02rLSHbLHHgbYHBs2D8dyEkHDa+a67rhMfByb8069I3VMH0Ww78vRqqTPdP+rhj+sikbygoKGDWrFkAObZt57Z/rkcEgGVZPwOuBbzAcuAG27YbOrwmGwVA9ylcDa9eacZyp4w0B7ofL4dPfrd/Kr7XYx4bq8xBqywXTrnTHOwWPWwOTFWFpk/6jPtg93pzHaJojemqOtqlLVqaTMulYDmc8duDr1m0NJuZ0iHhpvXR/mD41V9h4Z9MbcPOgpP/2wTgoFNMcLWp2m3WYqorMyNb9m422ydfZ5YQjsuAGbfuP3tuaYZ3b4WvXzRfh0ZBVD/TtZYyEk7/Dcz/hVnCICrFXIupLIC4LHP3qKBQ8/4lG/fXkDjI3C40d5Hpphlxng7icsyOFAAB7wKyLGsq8H1gClAL/A24BfhDIOtyvLQxcP2n8M8bzbC+8Vea4XyzHz/4tWvfhH/dDANPMGuou1zmgPvFI6ZLqbke1v0LGirMyoxRSfDyHJj1Kzj+JqjdAxveMWfSYdHmbHjYGbDqZfOz68vNEEMwE4VOu8ec5b7yXRNUQSFQU2yeHzCpdTnfVHNxc+FDkDnVnJUvfXp/CyUsBk6725y1531p5lM01+/v677yn+b9lz0D7jhT+1fPQ/xAMxO7odosQjbtZhhzKSy4xxz8z7jPTNZ78SJzhn/JCzDifBOWnz5gfs/zHjYH+U0fwMXPm5bU2tdN91p8pgkakW4Q8BaAZVnDgDTbtj9t/frnQJZt2z/u8Lps1ALofl6vObseOP3IszJrS83BPjTcfN1QbbpDhp9nzqY/ecCc8U+4CnDBm9eZLg13rBm5BJBkmT7w0m2YxqALRpxrDqQ5M83rV75krmWsewt2fGH6uJvrYdS3zQ0+lj5tWiBVRebC6cgLWpfMCDVhsWsFpI6G939uWgJgzs6HnAqn/casydTG44GSDaaugqXm2kl9heleqSo0C/Ud6mC96QPTBz/z9iN/Zh6PhuSK3/X4LqA2lmWlAcuA79m2/Z8Oz2WjAOg7vF7I/Ry+fsmMPBp5gXkEqNhpWgPp4w/sAqkrh6dOgIp8c9Y8+ykY+51D//zGWjNDOmXkoScetTSbg3tLoxkJ09kF9prq94edSA/WI7qALMu6DHisw+bdtm2Pan1+IPA+8GTHg7/0QS6XGYaYc9LBz8UNgEnfP3h7RDzctNgsCRweb7p5Dics8sizn4NDjm12tA7+0gd0WwDYtv0y8PKhnrMsaxzwDvB727Yf7a6apBdyR0OyFegqRPqEnnAROBmYD9xs2/brga5HRMQpOh0AlmWdD7xk23ZMh+3XAbcDGcBK4Ge2bX95FD/yp0As8D+WZf1P67Z5tm3f2dnaRETk6HUqACzLmg68iJlT3n77VcBTwD2Yi7g/Bj6wLGusbdvbj/Qzbdv+JfDLztQhIiLH7qgCwLIsN3Ar8BugBghr95wLc+B/2rbtX7du+xCwMWf3t/i4ZhER8YGjHYR8FvAL4OdAx4u0Q4CBwNttG2zbbgLmAWf6oEYREfGDow2AZZgxpI9gZui0N6z1cUuH7duAwZZl6e4PIiI90FF1Adm2vfMIT8e2PlZ12F6FCZgooLLzpR0kGKCoqMgHP0pExBnaHTMPOhn3xTDQtgvCHVsGbds9+EYawJw5c3z040REHCUN2Np+gy8CoKL1MQbY3W57NObgX+OD9wDTDXUiUAi0+Ohnioj0dcGYg/+yjk/4IgBa181lEAdeBxgE2LZt+2SxodbloRf64meJiDjM1kNt9MVShJuBfGB22wbLskKBc4AFPvj5IiLiB8fcArBt22tZ1gPAY5ZllQGLgJuBJOBPx/rzv0lnZyBblnUc8DAwFSgFHsesQdRzlkX1sy58Zu9iAr2jGNu2q/1TZc9zuFnwh3id4/exNp34zBy7j7WOlLwVuA7IAnYATwCPH26f8dU+5pPFyG3bfgIzR+BK4HUgHjjDtu1tvvj5h9NuBvKLwEVAOWYGcs5hXp8C/AdzwfpS4Gngt8Bt/qyzJ+nsZ9ZqDGZnm9bhX61/q+05DjcL/hCvc/w+1uZoP7NWTt7H7gLuw3xW5wOvAn/GHFMP4st9rEfdD6AzWmcgbwfet237xtZtoZgZyO/atn3QDGTLsn4N/Ahzw5na1m2/AW4C+rdOYOuzuviZxQNlwFm2bc/vznp7gkPNgrdtO/oIr3f0PgZd+swcu49ZlhWEOQl72Lbtu9ptfxy4xLbtlEN8j8/2sd58O6KuzEA+FVjQ9qG1+ieQCEz2U509SVc+szGtj6v9W1qPdaRZ8Ifi9H0MOv+ZOXkfi8PcBvfNDtttINmyrKhDfI/P9rHeHABdmYE87DCvb//z+rKufGZjgAbgXsuy9lqWVWtZ1muWZfX3Z6E9yJFmwR+K0/cx6Pxn5th9zLbtMtu2b7Zt++sOT50HFNi2fahh9D7bx3pzABzNDORDfc+hXt/+5/VlXfnMxgDu1td8G9PMnAZ81NrU79Ns295p23Z5J77F6ftYVz4zR+9jHVmW9QPMWf7vD/MSn+1jAb8hzDHoygxk1yFe38ZXM5Z7sq58Zg8B/7Bt++PWrz+zLGsDsBhzAWquz6vs3Zy+j3WF9rFWlmXNwQzSeJ2Db6Hbxmf7WG9uAbSfgdzekWYgVxzi9THtnuvrOv2Z2ba9sd0fZtu2JZgLV2P9UWQv5/R9rNO0jxmWZf0UE3bvAnOOMKTTZ/tYbw6A9jOQ2zvSDOTNh3k9mIsufV2nPzPLsi6zLOukDttcmCb7Hr9U2bs5fR/rNO1jYFnWfZiW0FzgYtu2G4/wcp/tY709ADo7A3kBcGqHK+uzgb2YCVF9XVc+sxuBh1uHq7U5G4gAPvNTnb2Z0/exrnD0PmZZ1q2YUVMPA9+zbbv5G77FZ/tYr70GcDQzkC3LGgwk27a9uPXbnsDcrvI9y7L+gGle/gK44xsSt0/o4md2H/A+8KJlWc9jRhn8BnjDtu0vuvt36Gm0j3We9rH9LMtKA34HrAFeBqZaltX+JcsxQ7f9so/15hbA0cxAvgv4st3rCzFX10NaX389cKdt23/szroDqQuf2QeY2YlDMGON7wSea/1+0T7WFdrH9jsD09U1GvOZdPwXjx/3sV47E1hERI5Nr24BiIhI1ykAREQcSgEgIuJQCgAREYdSAIiIOJQCQETEoRQAIiIOpQAQEXEoBYCIiEP9P+IxYtmCAMISAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(delta_t_list,4*result_Bayes_array2.T[3])\n",
    "plt.semilogy(delta_t_list,result_Bayes_array.T[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dt_tau(r):\n",
    "    return np.exp(r)/r*np.sqrt(1-np.exp(-2*r))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11b9a7da0>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD6CAYAAACrklzBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XHW9//HXTCaZrJO12dOk6XJKd7pQ1kIpCAVBZFOsil5Fr4Jy9YpXrj8XRFQQ8aKCC4p6gSuyiYhsUspS6AKUbrQ9LW26JGnTNM2eZp/fH99Jm4akNOkkZ3Lm/Xw88pjmzJJPJqfv+Z7v+X6/xxMMBhEREffxOl2AiIgMDwW8iIhLKeBFRFxKAS8i4lIKeBERl/I5XUAPy7L8wDxgL9DlcDkiIqNFDJAHvGnbdlvvOyIm4DHh/prTRYiIjFJnAct7b4ikgN8L8NBDD5Gbm+t0LSIio8K+fftYsmQJhDK0t0gK+C6A3NxcCgsLna5FRGS0eV/Xtk6yioi4lAJeRMSlFPAiIi6lgBcRcSkFvIiISyngRURcyh0Bv/LX8JdPOF2FiEhEcUfAH9wBu99wugoRkYjijoD3p0BrA+jqVCIih7kk4AMQ7IKOQ05XIiISMVwS8Cnmtq3B2TpERCKIOwI+PtXctjU6W4eISARxR8D3tOBb1YIXEenhkoAPmFt10YiIHOaSgFcfvIhIX+4I+PieFrz64EVEergj4NUHLyLyPi4JeLXgRUT6ckfAe2MgNkl98CIivbgj4MH0wyvgRUQOc0/A96xHIyIigKsCPqA+eBGRXlwU8CnqohER6cU9AR+vFryISG/uCXj1wYuIHMVFAZ+qFryISC/uCfj4ALQ3QneX05WIiEQE9wR8z3IF7U3O1iEiEiFcFPCh5QrUDy8iArgq4HuWDFY/vIgIuCng43XRDxGR3twT8FpRUkTkKO4L+NZ6Z+sQEYkQLgr4not+KOBFRMBNAZ+UBR4vNFU5XYmISERwT8DHxEJyDtRXOF2JiEhEcE/AAwQKoEEBLyICrgv4fAW8iEiIuwI+tdB00QSDTlciIuI4dwV8IB86mjWSRkQE1wV8gblVN42IiFsDvtLZOkREIoC7Aj41FPD15c7WISISAdwV8Mm5ZrKTWvAiIi4L+BifCXn1wYuIuCzgwYykUReNiIgLAz61QF00IiK4MeBT8hXwIiK4MeCTssxkp45DTlciIuIoFwb8GHPbfMDZOkREHObCgM8yt83VztYhIuIw9wV8YijgW2qcrUNExGHuC/jDLXh10YhIdHNvwLco4EUkurkv4P0B8MaqD15Eop77At7jMa34ZvXBi0h0c1/Agwl4ddGISJRzZ8AnZqmLRkSinjsDPilLo2hEJOq5NODHaBy8iEQ9dwZ8Yia0N2k9GhGJau4MeE12EhFxa8CHFhzTSBoRiWLuDPie9Wg0Fl5Eopg7A17LFYiIuDzgm/Y7W4eIiIPcGfD+gOmmOWA7XYmIiGPcGfAeD+TNgL3rnK5ERMQx7gx4gLyZsH8zdLY5XYmIiCPcHfDdnbB/k9OViIg4whUB39UdpLWj6+iNeTPNrbppRCRKuSLg735xKx/77YqjN6aPA38q7F3vTFEiIg5zRcA3tHay40Dz0Rs9Hsidrha8iEQtVwR8IN5HU1sn3d3Bo+/ImwlVG6Gz3ZnCREQc5IqAT4mPJRiE5vbOo+8oPh06W6H8TWcKExFxkEsC3gdAY2ufgC85Ezxe2LHMgapERJzlkoCPBfoJ+IQ0KJgDO14e+aJERBzmkoDvacF3vP/O0oVQ8TYcqhvhqkREnOWygO98/52l50CwG3YuH9GaRESc5pKAN100Df214AvngS8Bdr0+wlWJiDjLFQEfCLXgG/prwfviIGeKGS4pIhJF3BHwCT0nWftpwQNkT4GqdyEY7P9+EREXckXA+31eYmM8/ffBA+RMg5YaXQBERKKKKwLe4/GQEh87cAs+Z4q53f/uyBUlIuIwn9MFhEtKvG/gFnz2VHNb9S40VkH2SZA/a+SKExFxQHQEfFImJOfC+r/Cvg0QKIDrV4M/eWSLFBEZQa7oogFI8R+jiwYgZ6oJd38AGirglZ+MXHEiIg5wT8AfqwUPJuABFn4bZn8aVtwLjftGpjgREQe4KOBjjx3w068ywT73szDnMxDsgl1vjFh9IiIjbcQD3rKsVMuyNlqWNSGcr5sS7+t/JmuPvBlw6S/B54fcGRCbCHtWhbMEEZGIMqIBb1nWGcDrwKRwv/aAF/3oT0ysWWVy98pwlyEiEjFGugX/xdBXZbhfeMCLfgykaL456drWFO5SREQiwogGvG3bn7Zte1hW/TrmipL9GXua6YeveHs4yhERcZyrTrLCIAK+aB7gUTeNiLiWiwL+GBf96E98qumHf+cBaG8exspERJzhwoA/zhY8wIduhfo98Oqdw1SViIhzBh3wlmVdallWYz/br7Msa5tlWYcsy1phWdZp4Snx+Bzzoh8DKT4dZl4Db/wSqrcOU2UiIs4YVMBblnU68CDg6bP908BvQvddAdQBz1uWNa6/17Ftu8S27feGVPEAAgnHuOjHsZz/AzMm/plvaL14EXGV4wp4y7L8lmV9E1gGdPa5zwP8APidbdu32Lb9DHApcAD4WpjrHVBGYhxeD1Q3tg3uicnZsOg7UPYKvPvE8BQnIuKA423BLwZuBm4CftnnvglAMfBUzwbbtjuAfwIXhqHG4+KL8ZKV7Gd/Q+vgnzz332DMZFj1u/AXJiLikOMN+DeBcbZt/wLo24/RMyu1b5fLDmC8ZVkxJ1DfoOQE4qkaSsB7Y2DCeVD5DnS2h78wEREHHFfA27ZdYdt23QB3B0K3fU+8NoZeP2mItQ1aTsBPVcMgu2h6FM6DrjYzu1VExAXCMUyy54Rr35Z9z/buMPyM45IdiGd/4xBa8ABFp5jb8tXmNhiEHS9Dd1dYahMRGWnhCPj60G1Kn+3JmHAfsVlEOSnxHGhqp6NrCJ8pgXwIFMKeUMCXvQr/+xHY/NSxnyciEqHCEfDbQrelfbaXArZt2yM29jAn4AeGMJKmR9E8KH/L/Pu9F82t1owXkVEqXAG/B7isZ4NlWbHAxcDSMLz+ccsJxAMM7UQrQOEpUL8bGvbC9mVmm9aqEZFR6oQvum3bdtCyrJ8Av7Isqxaz3vsNQBbw8xN9/cHIDrXgh3yitfQcc/vK7VC1ARLSoWojtDWCv28PlIhIZAvLWjS2bd+LGSP/KeAxIA24wLbtHeF4/ePV04If8onWnCkw9aPw9h/N96fdAMFu023TOcQPDRERhwy6BW/b9veB7/ez/WfAz068pKHLSIzD5/UMvYsGYNF3YfPTpsU+73Pw0g9h2W1QuRY+9gBYi8NXsIjIMDrhLppI4vV6yE45gbHwABmlsPh2IGi6aHKmQvmb5r71jyjgRWTUcFXAgxkLf0IteDAt9x5nfs1Mfmqqgi3PmJmuvrgTe30RkRHgmvXge+QE/Ow/kRZ8X9OvhPNvgZMugbZ62K1hkyIyOrgw4OPZW38o/C9ceg744sF+LvyvLSLRrWsQ17EYBNcFfFF6Ig2tndS1hHnRsLgkE/Lr/wp714X3tUUkutSXw9q/wJPXw//MgB/lm21h5ro++JIss7ZZ2YFmTh4b5r7y82+FBz4Kf7wILv8dTL44vK8vIu4TDMLBHbBnFex6HXYuh9qd5r6EdCg+A876TwgUhP1Huy7gx4UCfmdNMyePTQ/vi4+ZBJ9/ER6+Bh7+BFgXg8djTsQWzg3vzxKR0au+HHatMEue7FhmBmkAxKdC8Zkw/9+h5EzIngre4etIcV3Aj81IxOuBsgMtw/MDAnnw2efgX9+FLU9Dy0Hoaocljw7PzxORyNdQaZY32fEy7F4B9XvM9oQMmLDIXP+56FRzYaFhDPS+XBfwcT4vBekJ7DwwjItYxsbDRXeYr6W3wvK7zB84kD98P1NEIkN3N+xbB1v+CRVvm+6Wg6FJ+0ljTJfLaTfA2PmQO8NcUMghrgt4gJLMJHbWjNAqxSd/El67E168BfatN8MpF/73yPxsERl+PX3oZa/AjlfMUuKHDoLHC7nTIWeauexn6UIzMdLj+eDXHCGuDPhxWUn87Z0KgsEgnuF+szPGQclZsP5h8CWYhcrikuHUL0OMK99eEXfraIW9a80M9j2rzVpUjZXmvkABTLoQSs+GCedDUqaztX4AVyZQSWYSja2dHGxuJzPZP/w/8LzvwzsPmpb701+Df33HtOov+JFp4YtI5Nq/Bba9AIkZpnX+7t/MeTWAtGIoOQOK5psWeub4iGqhfxBXBnzvkTQjEvCFc4+MornqT+bk64p74JmbzNj51MLhr0FEjk9XhwnxLf80o1t2r+TwFUfjUmD2tTB+oblOc3K2o6WeKFcG/JGx8C3MKc4Y2R8eE2uWHM6fDffMh+e+BVc/MKo+9UVcJRiE/Zvg3SdNP3r1FmitN5foTC2ABd8wfegdhyA5B/zJTlccNq4M+ML0BOJivGzb3+hcEenFcPZNsPQHcN9C+PD/QP4s5+oRiQblb8Hq+8xkoqYqSCsyJ0Nr3jO3hfNg6uVmVdgJ54/okEUnuDLgY2O8TMxJZlNlg7OFnPE1SMkzI2we/xx8eZVpyXu8atGLnKidr4P9DFTbpo+8egtsfMx0s0y6wAxbri0zLfNTv2xGuI3yLpfBcmXAA0zJC/DSlv0jM5JmIF4vzPqEGVXzyKdg9e9g7f9BfMBMjIpLcqYukdGkq9PcdneaOSdV70LzAdizEmL85mj5vX9BTByc/S04/Suu6mY5Ea4N+Kn5AR59u5zqxjayQ5fyc8zkD0POdHj+ZrNDdnfAI9fCNX8xffYicrR9G80Eov2bYeU9ZnJRfCrU74Ysy0weuuDHMPezEJsADXvNUXFKrtOVRxTXBvyU/FQA3q1scD7gvV44//vw+HXw0d+YWa9P/we8eiec/U3TXzj2NF1IRKJXezPU7oK6XeZk6PqHj9w38UMmuGu2wyU/hwnnvf/5gbyRq3UUcW3AT85LAWDT3gYWTo6AfrcJ58FN24+c1Nm9Al79qVm7Ys9KM7356gcifuKEyJDtXW/Gm3e2wdTLzPmpl35o/g8c3H7kcTFxZnXFKR8xrfb0EqcqHvVcG/CB+FjGZiTybmW906Uc0fuM/eLboew1qHwHTvkivP0n+P25cM1fYYxlHqMTsTLatBw0S+CCWUnxjV+Y/vKcabDxcQh2mfte+xkkpEFrgzkhOusacz3ktBIzmSghzbFfwU1cG/Bg+uEdH0kzkIR0+NwLZsZc5niYcTX85RozpDImzpw4+vxS9dFL5Os4BNv+BSvvNUemebOguwuqNpip/RmlZmLRjI/BBbeZ57z4PdO/fvFdkDfD2fpdzPUB/+zGfdS3dJCaGIFBmVZ05N+Fc+G6l2DZbeY/zKYnTav+lOugab/pl7QWH/0ckeFy4D1IyQF/irnQvNcH7U3mimZ5M6HoFGhrgtfvNsHe3mQmDp1xI7y31LzGZb+GaVeac0vB4NFHpJf+0pnfK8q4OuB7ZrGu2V0bGf3wHyStyJyEDQbhTxfDyz8x06g3/d2MvNnwKPzb866fnCEOW/cwPPklM5684GRz4QqP14xcaW8CT4xpjW97AVoOwJTLYM61ULLALLB3/g/e/5rqbnSEqwN+VlEaPq+HN3ceHB0B38PjgQ/dCvctMoe+8z5vJmgsvQXe+oNp1ffo2zISGYyDZWZCkNdnWt72M+bIseRMSMoyY87nfMbc39EM068y8znW/R9MvAAW3ARF85z+LWQArg74hLgYphak8tbOWqdLGbyCOXD9arNWRlySCfKyV+DF78O4BeZEbGc7/OYMsC6C829xumKJZF0dJqzBrGFeX26W0dj4mDkfFJdixpjHJsHMa+DDd5nx5f0pPsOcTNWIr4jn6oAHmFeczgMrd9HW2YXf59yVVYZkzKQj//Z44CP3wO8Wwl8+bvrr330SDmw14+rP+rpZQMmfYv7DtjWZVlesw3MAZOR1d5kLPL+3FA7Ypj/94PYjS+B6fWZWqNdn+szry+FQHVzwQ7PWue8DVmD1eBTuo4TrA35uSTq/X17GxooG5hSH+SLcIy21ED72oOmf/+PFJtBTx5qW1wvfgQ2PmREL1z4Fv18EsYnwuX9BXKLTlctQ1VeYv3fxGbDwZnNJuJ4A3rfBzOBsbzKjV7o6zDDErc+bhbY8MeaCNJkTYeL5kH+yCf9968w656ULIWuCs7+fDCvXB3zPida3dh4c/QEP5jqPSx6BRz9jAn7JY2bC1Jo/m0khVRvgt2cfuejvP//TXIgkPmAOwwc6QdvWCHV7IGfKiP0qElJfAQSPXDegq8OcWC8+w1xApqHSzOxc+6C5f8pHIHuKOQnfs455bJI5WutsN2uZT/mIGV/uT3n/z5tx1Uj8VhIBXB/wY1L8TMhOZvl7B/ji2eOdLic8xp8LX3jZjLCZcJ45ZP7nN+DqP8PrvzD9qmf8hzkEf+1Oc0IMzPdFp5rW3KQLIXvykdd8+utmIspnnzUfIv3p7jaH+pkTdGI3HCrfgWU/MifSwfxdSs40XW+Va8Aba0ZPXfBj83d+70VoqDDL4W76O0y7AuZ/ycyVyJmqORPyPq4PeICzJ43hgZW7ONTeRULcKOuHH0hGqfkC85//xrXm3x++y5yEnflxExBjT4XGvWbGYNM+2P6ymWTy4vfMld8/9ENoa4DNT5nD+8c/B1feD4mZZlhcYuaRVuDyu+ClW82V4hffDsWnO/Krh8X+LWZkUuJxXBCmu3vgI5/uLrMoVltj/yHbVH1kZcMHrzT94NOvhF2vw6anzPu74CZz/7qHzdDD+FS45Bem26WzFeZ/0fwtes7JzPms6WaZerk+aOWYoiLgF0wawx+Wl7GqrIZzrFE0XHIo4lPNmOQeE89//2PqK0zLfsWvzMm1gpNNkFx0Jzx3M/yhz3PSiuHin8Frd5nRPc3V8NdPwQ1vmqMCb8yRpY+7u805gbTikQmf7m6zls+e1aZFe6yJYN1dgAc2/x0e/7zpm/7cC+YDrL9aGypNN8j6R8zaKGMmwfPfNuO8p11ulpp49pvmakFgjo4W3gwv326OgkrPgYeXmKsEZU0yoZ5aaJ6TmGmWtV3wDfM3Azj329BcY/rY/clH/x17y5qgvnM5Lp5gMOh0DQBYllUClC1dupTCwvBew7S1o4uZt7zAJ+aP5XuXTA3ra49awSC8cju8/GMT0pkT4MsroW63CaxDdRDsNmG++j5oKDdHBNevMh8Gv10ABXPNRRYS0uBjD4XGUP/ZPPaMG+G8Wz445FvrjwRp+Zum66nnogydbaZl3Nlqup7iA3Dm14+cNO7qMMs7vBfq4kjJg9O/CtuXmrkD1uIjP2fbi/Dkv5vX62qH7KlQvdmEfONeM+z0o781JyUBGqvMieqmKjP1vny12R4TZ96HOdeaGZxpY013WFeHudh6Vzv4U6EttAZS5kRzErRxLyz8tvmgqHkPMsabSUEiJ6i8vJxFixYBjLNte2fv+6JiD4uPjWF+aSavbK12upTI4fHAOd8yQyqf/aaZzOLxmDVw0ouPfuy0K+DRa2HyxWbdHDBB++odUHiKCazfnGG2j19kll14/W7AY07wNlTCOw+YvuXETNNtNOUy0+Wz5WnzPF+CCcdX7oCL7zRdG/eda44G4Eh/9Lq/QkaJCdb2FhPu591iTkg+8imz5n5sEmxfBou+Y9ZCsZ8xPzt7ipmB6Y2Bs//LzAx+4bvmpOT2ZXDPKZCQYbq+Wuuhpca08PNmmck9jXth9qfNBLSV95rXuuTuI+PF82eZ/vTTv2I+ZLY8Y7qywCwJfdKlpqunZzE5kWEWFS14gD8sL+PWpzfx8jfOOXxRbgmprzCzGQfTpdLdBRVrTJfNwe2w/Odw8idNv3x3N/zjqybUEzNNUHq8psuio9W01Ls7TGv4tOtNoI4/13QBrfkzfPJxs1Z++Zumf7qrA05eYmZdLv85dLSYPvS2etNi/9CtpqbGKnMEkj0ZHv4ElL1qtidkmOef898DDxmt3QVv3meOXPatNwthXXm/ucxbXxVrzJHLzGvUBy6OO1YLPmoCvry2hTNvX8Z/XTiZL53jktE0kW77S6Z7J2+WCf/UArO9cR9s/odpdfceltlYBb+YZQIc4LLfmGVk+9PZbkI2d3r/Idvdba7H2dVuukMGezGV7i7T0heJcFHfRQNQmJ7IjMJUntu4VwE/Usafa776Ssk9ej2dw9tzzMSs2jJIzj32Gie+uGMvM+v1HulOGgqFu7hAVC1LeOG0XNaV11Ne2+J0KTKQ3GmmW0QLWImcsKgK+MXTzHUbn9u4z+FKRESGX1QF/LisJKbmB/j72kqnSxERGXZRFfAAV8wuZENFPfa+RqdLEREZVlEX8B+ZlY/P6+HxNeVOlyIiMqyiLuAzk/0snJzNE2sq6OzqdrocEZFhE3UBD3D13CIONLXx4uYqp0sRERk2URnw507OpiAtgT+/scvpUkREhk1UBnyM18MnTy1mxY4atlbpZKuIuFNUBjzAx+YVEefz8sfXdzpdiojIsIjagM9IiuOK2QU8vqac/Y2tTpcjIhJ2URvwAF9cMJ7Orm7+sLzM6VJERMIuqgO+JCuJi6bn8dDK3dS3dDhdjohIWEV1wANcv3ACTW2d/O617U6XIiISVlEf8CflBbhkZj73L99JdWOb0+WIiIRN1Ac8wNfPn0R7Vze/emmb06WIiISNAh6zyuTH5xXx4KrdvLdf4+JFxB0U8CFfP38SiXEx3Pr0ZqdLEREJCwV8SGaynxsXTeSVrdW6IIiIuIICvpdrTy9han6A//fkRupa2p0uR0TkhCjge4mN8XLHlTOoa2nnB09vcrocEZETooDvY2p+Kl8+ZzxPrKlg2Zb9TpcjIjJkCvh+3HDuRKycFG5+YgP1hzTDVURGJwV8P+J8pqumuqmNbz62jmAw6HRJIiKDpoAfwMyiNG5ePJnn363ivtd2OF2OiMigKeCP4XNnjuOi6bnc/pzNqh01TpcjIjIoCvhj8Hg83H7FDIozErnhL+9Q1aB140Vk9FDAf4CU+Fh+/ck5NLd18tk/vkljq066isjooIA/DlZuCvcumc3Wqkb+/cG3ae/sdrokERnl9je08ve1Fdz8xHo++ftVw3JNCl/YX9GlzrGy+ckVM/jGo+v45mPruOvqWXi9HqfLEpFRoqapjZU7DrJixwFWbK9he3UzACnxPhZMGkNMTPjzRAE/CFfOKaSqoZWfPm+TEh/LLZdOVciLSL/qWtpZueMgK3fUsGJ7DXaVWak2KS6GeeMyuHpuEaeNz2Rqfioxw5QjCvhB+vI542k41MFvX91BZ3eQ2y6bppAXEQ42t7O67CCry0yob97XQDAI8bFe5pVkcOmsfE4bn8n0glRiY0amd1wBP0gej4dvLZ6ML8bDPcu209XdzY8vnzFsn8AiEpmqGlpZVXaQVTtqWF12kG37mwDw+7zMHpvO186bxGnjM5lZmEacz5nTnQr4IfB4PHzjQxY+r5e7l26jqa2Tu66eRXxsjNOlicgwCAaDlNceYlXZQVaX1bCq7CC7aloASPb7mFOczmUnF3BqaQbTClLx+yIjCxTwQ+TxePja+ZNIifdx2zOb2Ve/kvs+PZfMZL/TpYnICQoGg2yvbg51uZhA31tv5sGkJcYyrySDT51azPxxmZyUl4JvhLpcBksBf4I+f1YphekJ3PjwWj567xvc/5l5TMhOdrosERmEjq5uNlU28NauWt7eZfrRDzSZa0JkJfuZX5rB/HEZzB+XycTs5FFz3k0BHwYXTsvj4S/Ec93/vsVHfrWcO66cycUz8pwuS0QGUN/SwZrdtby16yBv7axlXXkdrR1mfktBWgILJo7hlHEZnDIug3FZSXg8oyPQ+1LAh8nJY9P5x1fO5PqH1nD9/63h7V3juPmiySN2tlxE+hcMBtlV03K4df7WztrDJ0RjvB6m5gf4+LyxzC1JZ25xBrmp8Q5XHD4K+DDKS03g4S+cxo+f3cz9r5fx1q6D3HX1LHXZiIygts4uNlY0HA7zNbtrD3e3pMSbE6KXzsxnTkk6s4rSSIxzbwy69zdzSJzPy/cumcr8cRnc/MQGLv7Fa/zXhZP5zOklo6bfTmS0CAaDVNa3snZ3He/srmXtnjrWV9QfXk6kODORBZPGMKfYtM5HU/95OCjgh8mF0/KYXZzOzY9v4AdPb+K5jfu47aPTmJiT4nRpIqNWU1sn6/fU8c6eOtaGvqob2wDTuJqWH+Da04qZU5zO7OJ0slPc090yFAr4YZSdEs/vr53LY2+Xc9szm1l892tct6CUr547kYS4yBgnKxKpOru62VrVFApy0zrftr+JnguslWYlcdaELGaNTWNWURqTcwOOTSiKVAr4YebxeLhqbhHnTs7mx89u4dcvb+fJdyq46QKLy2YVRNXhoshAeiYSrS+vZ31FHWt317Ghop6W9i7AjD2fVZTGRdPzmFVkAj0tMc7hqiOfAn6EZCb7ufOqmVw1p5Af/nMzX39kHfe/XsbNi0/i9PGZo3YYlshg9Q7zDRX1bKwwtz0XuI+N8TAlP5Wr5xYdDvPizET9HxkCBfwIm1+ayd+vP4On1lVyx3NbWPL7VcwtTueriyZy1sQs7cTiKscT5lZuChdNz2VaQSrTC1KxclMiZqr/aKeAd4DX6+Gykwu4cFouj761h1+/vJ1P37+amUVp3LhoAgutbAW9jDpd3UF21jSzeW8DGysa2FhRz8bKeupaFOZOUcA7KD42hk+dVsLH5o3l8TXl3LPsPf7tT28xfkwS155ewuWzC0n2608kkaeprZMtexvYvLeBTXsb2by3AXtfI4c6TJ95bIyHSTkpLJ6mMHeS0iMCxPm8XHPKWK6cU8g/1lXy5zd28t2/v8sdz9lcOaeQT546lgnZGl4pI6+ni2Xz3gY2h4J8876GwyspAqQmxHJSXgofP6WIk/ICTMkLMDEnWWEeARTwESQ2xsvlswu5fHYha/fU8ec3dvLQql386Y2dzCpK4/LZBVwyI5/0JI0ekPAKBoNUN7WxraqJrVWNbK1qYltVI3ZVI42tnQBT60gEAAAKAElEQVR4PFCckcjU/ABXzi7kpLwAJ+UHyE+NV5dihPIEewaVOsyyrBKgbOnSpRQWFjpdTsSobmzjyXcqeHxNOVv2NRIb42Ghlc2HZ+az0BpDSnys0yXKKHOgqY2tVY2Hw3xbVRNb9zce7isH0yqflJPMxJwUpuQFOCkvwOTcFJLUZRhxysvLWbRoEcA427Z39r5Pf60INybFz3ULSrluQSmbKhv42zvlPLm2khc2VREX4+WMCZlcMDWX86bkkKW16CWkqztIZd0htlc3UXagmR3VzSbM9zdxsLn98OMC8b7DfeUTs1OYlJPCpJxkxqT41Sp3AQX8KDIlP8CU/Cl8a/FJrNldy/Mb9/H8pn0se2IDnr9tYHpBKmdNzOKsiWOYPTZds/qiQG1zOzsONLGjupkdB5opq25mx4Emdta0HF6PBcwiWxOzk/nQlBwmhkJ8Uk4K2QpyV1PAj0IxXg/zSjKYV5LBty8+ic17G3lxcxWvbavmN6/s4J5l20mMi2FeSQZzi9OZW5LBrKI0LY8wCnV3B9nf2Mbugy3sPtjCntDXzhoT6L27VXxeD2MzEynNSuYcK5vSrCRKxyQzLiuJrOQ4BXkUUsCPch6PJ9SyD/DVRRNpaO1gxfYaXt1azeqyg/xsazVg/vNPzQ8wpziDk8emMSU/QElmki4WHgGa2jopr21hd82RED8c6LWHjmqJezyQF4inKCORxdPyGD8miXGhIC9KT4jYS8eJMxTwLhOIj+WCqblcMDUXgLqWdnPlmp21vLWrlodW7eL+18sASIiNYXKeOYk2JT/A5NwA48ckaY2PMGrt6GJffSuV9YfYW9dKZd0hKutb2dvzff2hw6NUeqT4fRRlJDIxO4VFJ+VQlJHI2IxEitITKEhP0PBDOW4KeJdLS4zj3Mk5nDs5B4D2zm627W9kU2UDm/Y2sKmygafWVfLQqt2Hn5OeGHv40L50TBLjMpMoSE8gPy2BzCQd6geDQZraOqlubKO6sY39odvqpqO/39/QSk2vE5o9MpLiyE+LZ2xmIqeWZpCXlkBBWgJjQ0Gelhgb9e+xhIcCPsrE+bxMzU9lan7q4W09k1nsfY1mxEXopN0rW6t57O3yo57v93nJT0sgPy2evNQEspL9ZCXHkZkcR0aSn8ykOLKS/WQkxY2Kk7wdXd00t3VSf6iD2pYO6lraqWvpoDZ0W9fSTm3oe/OYdqob2w5fv7M3n9fDmBQ/Y1L8FKTFM6solbxU88GYnxpPXloCeanxxMeqBS4jQwEveDweijISKcpIfN99ja0d7KppMV0Loe6FitC/l287QE1zGx1d/c+lSIqLITneR7LfR3J8LCl+H0n+GJL9sST7Y4iPjSHO5yU2xnvUbVyM56hr2faeqtH7JwWDQdq7umnv7PUV+r6t178PtXfR1NZJc1unuW3vpLnNbOvdv92fQLyP9KQ40hLjSE+MozQriaxkP9kBE+RjkuMZk+InO8VPakKsln+WiKKAl2NKiY9lWkEq0wpS+70/GAzS2NZJTVM7NU1tHGhq52Cz+XdtS8fhUG0MBWx1Y5v5vrXjcAiHc66dxwNxoQ8Kv89LXIyXRL+PJL+PZH8M6UmJJIc+aJL8PpLjzH2BhFjSE2NJS4wjLTGW9MQ4UhNidRJaRjUFvJwQj8dDID6WQHws47KSBv38YDBIV7dpiXd0Bmnr6qKjK0hHZze9u6E9HPmm9/a4UIjH+cyXz+tR/7VIiAJeHOXxePDFeMzwvjgALb0gEi6RfxZMRESGRAEvIuJSCngREZdSwIuIuJQCXkTEpRTwIiIuFUnDJGMA9u3b53QdIiKjRq/MfN8aGJEU8HkAS5YscboOEZHRKA/Y3ntDJAX8m8BZwF6gy+FaRERGixhMuL/Z946Iuei2iIiEl06yioi4VCR10fTLsqzrgG8ChcBa4Ou2ba84xuOnAXcD84GDwD3AHbZtR82hyhDes6eBi/u5K8W27abhqTLyWJZ1KfCQbdspH/C4qN/HegziPYvafcyyrBjgRuA6YCywC7gXuGegfSZc+1hEt+Aty/o08BvgQeAKoA543rKscQM8Pht4EbNs+NXA74DbgP8ckYIjwGDfs5AZmJ3ptD5fLcNbbeSwLOt0zHt2zKUotY8dcbzvWUg072PfAX6Eea8uBR4B/ge4qb8Hh3Mfi9g+eMuyPEAZ8Kxt218KbYsFbOBp27a/2s9zbgGuB8batt0S2nYr8GUg17btjr7PcZMhvmdpQC2w2Lbt50ay3khgWZYf07q6FWgG4mzbTj7G46N6H4MhvWdRu49ZluXFNLLutm37O7223wNcZdt2dj/PCds+Fskt+AlAMfBUz4bQL/ZP4MIBnnMesLTnTQl5EsgA5g1TnZFkKO/ZjNDt+uEtLWItBm7GtKZ+eRyPj/Z9DAb/nkXzPpYK/C/wRJ/tNjDGsqz+LqIQtn0skgN+Uuj2vT7bdwDjQ/1a/T2nv8f3fj03G8p7NgNoA35oWVaNZVktlmU9allW7nAWGkHeBMbZtv0Ljr4i4ECifR+Dwb9nUbuP2bZda9v2DbZtv9PnrkuActu2m/t5Wtj2sUgO+EDotrHP9kZM3f198gUGeHzv13OzobxnMwB/6DEfxRwGnga8FDoUdzXbtits264bxFOifR8bynsW1ftYX5ZlfR7TSr9jgIeEbR+L5FE0PSdu+rYQerb3d7VkTz+P73Hsqyu7w1Des7uAv9i2vSz0/auWZW0GVmJO8DwQ9ipHt2jfx4ZC+1iIZVlLMIMgHgN+NcDDwraPRXILvj5023f4VTLml+zv0Ka+n8en9LrP7Qb9ntm2vaXXf7yebaswJ4ZmDkeRo1y072ODpn3MsCzra5gPs6eBJccY8hi2fSySA35b6La0z/ZSwB7gzdk2wOPBnNRwu0G/Z5ZlfdyyrAV9tnkwh9QHhqXK0S3a97FB0z4GlmX9CHMk8wBwpW3b7cd4eNj2sUgP+D3AZT0bQkP+LgaWDvCcpcB5fc5MXwbUYCb8uN1Q3rMvAXeHhnP1uAhIAF4dpjpHs2jfx4Yiqvcxy7JuxIw6uhv4jG3bnR/wlLDtYxHbB2/bdtCyrJ8Av7IsqxZ4HbgByAJ+DmBZ1nhgjG3bK0NPuxf4CvCMZVk/xRz+3Qx86wM+MV1hiO/Zj4BngQcty/oj5iz9rcDjtm2/MdK/Q6TRPjZ42seOsCwrD7gd2AA8DMy3LKv3Q97CDG0eln0sklvw2LZ9L2as7acwJyXSgAts2+4ZMvQdYEWvx+/FnJ32hR7/BeDbtm3fOZJ1O2kI79nzmNl1EzBjbb8N3B96vmgfGwrtY0dcgOmKmo55T/p+pTGM+1jEzmQVEZETE9EteBERGToFvIiISyngRURcSgEvIuJSCngREZdSwIuIuJQCXkTEpRTwIiIupYAXEXGp/w8/TCo8YSE2fAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(delta_t_list, dt_tau(delta_t_list))\n",
    "plt.semilogy(delta_t_list,500*result_Bayes_array.T[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/hstrey/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \n",
      "/Users/hstrey/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "     fun: 2.485267512660051\n",
       "    nfev: 16\n",
       "     nit: 12\n",
       " success: True\n",
       "       x: 0.7968121301739481"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimize_scalar(dt_tau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate A and B from the dataset\n",
    "def OUlargeN(N,A,B):\n",
    "    data1sq = A\n",
    "    dataNsq = A\n",
    "    datasq = (N-2)*A\n",
    "    datacorr = (N-1)*A*B\n",
    "    P2A = -N/2/A**2\n",
    "    Btmp = (N-1)*(1+B**2)/(1-B**2)**2\n",
    "    tmp = (2+6*B**2)*(data1sq+dataNsq) + (4+12*B**2)*datasq - (12*B+4*B**3)*datacorr\n",
    "    P2B = Btmp - tmp/A/2/(1-B**2)**3\n",
    "    PAB = (N-1)*B/A/(1-B**2)\n",
    "    dA = np.sqrt(-P2B/(P2A*P2B-PAB**2))\n",
    "    dB = np.sqrt(-P2A/(P2A*P2B-PAB**2))\n",
    "    return dA,dB\n",
    "\n",
    "def OUlargeN2(deltat,N,A,tau):\n",
    "    dA,dB = OUlargeN(N,A,np.exp(-deltat/tau))\n",
    "    dtau = np.exp(deltat/tau)*tau/deltat*dB\n",
    "    return dA,dtau\n",
    "\n",
    "def OUlargeNdtau(deltat,N,A,tau):\n",
    "    dA,dB = OUlargeN(N,A,np.exp(-deltat/tau))\n",
    "    dtau = np.exp(deltat/tau)*tau/deltat*dB\n",
    "    return dtau\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "dA,dtau = OUlargeN2(delta_t_list,10000,1,1)\n",
    "dtau2,dA2 = OUlargeN2(delta_t_list,10000,2,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEQCAYAAACTEVJTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XlYVdX6wPHvARTEAVRUnHFiOYXDdR5JTTPNodRMcyjzV9lNr5lmOZRps1ZWdh2uXTO910qtq9fMElPT1Chn0uWsOI+oiILA/v2x4VxmOXAO+wDv53l4YK+19z4vx+15WXutvZbNMAyEEEIIZ/GwOgAhhBAFiyQWIYQQTiWJRQghhFNJYhFCCOFUkliEEEI4lZfVAVhNKeUNNAfOAQkWhyOEEPmFJ1ARCNdax6asKPSJBTOp/GJ1EEIIkU+1B7akLJDEYrZUWLp0KYGBgVbHItzMwIEDAVi2bJnFkQjhXs6fP8/gwYMh6TM0JUksSbe/AgMDqVKlitWxCDfj5WX+F5FrQ4hMpetCkM57IYQQTiWJRQghhFNJYhFCCOFUkliEEEI4lSQWIYQQTiWJJReuXLnCgQMHrA5DCCHciiSWXBg4cCANGza0OgwhhHArklhyoUaNGhiGQUKCzAQjhBDJJLHkQnBwMIZhcOvWLatDEUIItyGJJReKFi0KwIULFyyORAgh3Ickllw4deoUALt377Y4EiGEcB+SWHKhXbt21KlTR+aREkKIFGQSylzo06cPffr0sToMIYRwK9JiyYW4uDh+++03zp1LN2u0EEIUWpJYcuHo0aO0bNmSSZMmWR2KEKIAWLduHUopIiMjHTru4sWL1K1blz/++MNFkTlGEksuVKhQAYDKlStbHIkQoiDYv38/fn5+VK1a1aHj1q9fT5kyZWjSpImLInOMJJZc8Pf3x2az4eEhb6MQIvf2799P/fr1HT4uLCyMTp06uc1nkXTe54KHhwe+vr4cP37c6lCEEPlIfHw8c+fOZeXKlVy+fJnmzZszY8YMIiIi6NevHwBz587lww8/zPQcgwYN4rXXXiM6OpodO3bw6aef2uuye6yrSGLJpdjYWLZs2WJ1GEKIfMIwDMaNG8f27dsZNWoUwcHBbNq0ieHDh3P9+nX7/IM9evSgVatWAKxdu5ZFixaxePFivL29gf/dgt+4cSNFihShTZs29tfI7rGu4h7tpnysa9eu9O3b1+owhChUQkNDWbRoEQB3794lNDSUJUuWABATE0NoaChfffUVANevXyc0NJSVK1cCcPnyZUJDQ1m9ejUA58+fJzQ0lB9++AGAyMhIQkNDWb9+PQDHjh0jNDSUTZs2AaC1zlXsy5cvZ926dcybN49hw4bRunVrJk6ciL+/PwANGjQAoGrVqjRu3JjGjRtz+/ZtKleuTMuWLe1l5cqVA8z+lfbt29tnAnHkWFeRFksurVmzxuoQhBD5yPz58+ncuTONGzdOVV6zZk2OHDlCtWrV0h1z6NAhlFLpyuPi4vjll1+yvK2V2bGuJIkllyIiIrh06RKhoaFWhyJEobFx40b7z0WKFEm17evrm2rbz88v1XZAQECq7cDAwFTbVatWTbVds2bNVNu5+ZCOjIzk1KlTjB49Ol3dxYsXqV+/PjabLVW5YRgcOnSIYcOGpTtm27ZtxMbGZvr5k9WxriS3wnKpf//+PPjgg1aHIYTIB5InrA0ICEhVfunSJcLDwzNc3+nMmTPcunUrw4QWFhZG8+bNKVWqVIavl9WxriSJJZfuu+8+SpYsaXUYQoh8oHTp0gAcOXIkVfknn3xCbGxshkONL168CKTvcDcMgw0bNtClS5dMXy+zY12t0NwKU0r5AVuBPlrrI/faP7uqVq1KTEyMs04nhCjAgoKCqF27NnPmzKF48eIEBgayZs0ae19tcsd9Sr6+voD5VH5CQoK9b2b37t1cunSJzp07Z/p6mR3raoWixaKUaouZVIKdfe7ExERiYmJksS8hxD15enry8ccfU7t2baZNm8aECRPw8PDg8ccfx9fXl6CgoHTHKKUYOHAgX3/9NePHj7eXr1+/noYNGxIYGJjp62V2rKsVlhbLM0lfS5194jNnzgDmkMT77rvP2acXQhQwtWrVsg+NTunll1/OcH+bzca0adOYNm1aqvKwsDB69+6d5WtldqyrFYrEorUeCrkbzZGZfv36cePGjXSdcUII4UrJz924owKTWJRSA4FP0xRf0Fqnv2npRP3796d///6ufAkhhMhXCkxi0VovA5bl9etGRUXxww8/0KZNmwwfbBJCiMKmUHTeu1JYWBiPP/44c+bMsToUIYRwC27XYlFK9QKWaq1LpikfCUwAqgC7gRe11tssCDGV4GBzoFleP4AkhBDuyq1aLEqpNsASwJamfCgwN6nuUSAKWKeUquHI+bXWQc58hgWwT+YWGxvrzNMKIUS+5RYtFqWUNzAGmA7cAoqmqLMBbwDztdbTksp+AjQwFkg/6U4eSp6R1F2WBBVCCKu5S4ulO/AKMB74JE1dbaA6sCq5QGt9F1gDWD5Jl4+PDx4eHmzevNnqUIQQwi24S2IJB2porT8GjDR1yU/Lp72FdQyopZTydHVw9zJy5Mg8fapVCCHcmVvcCtNan8miOnnazptpym9iJsbiwA1XxJVdc+fOtfLlhRDCrbhLiyUryR35aVsyyeWJeRhLhn7//Xfeeecdq8MQQgi3kB8Sy/Wk72nnpi+BmVQsn/3xpZde4pVXXuHGDUsbTkII4RbyQ2I5nPS9ZprymoDWWqdtyeS55PUQDMPyUIQQ+di6detQShEZGZmu7uuvv7bPVpyZixcvUrdu3QxHqWbneGfJL4klEuiTXKCUKgL0AMKsCiql6tWrA3D+/HmLIxFC5Gf79+/Hz8+PqlWrpiq/desWH3/8MQBa60z/iF2/fj1lypShSZMmOTreWdw+sSS1SN4BnlVKvamUegj4DxAAfGhpcEmSlwWVIcdCiNzYv39/hqtIzp8/n0uXLtGxY0diYmI4ffp0hseHhYXRqVMnPDw8cnS8s7h9YgHQWn+G+YzLEGA54A9001ofszSwJMlP32/ZssXiSIQQ+UF8fDyffvopnTp1IiQkhBEjRnDu3DkiIiLSJZZz586xaNEiunTpwqBBgwCz1ZFWdHQ0O3bsSLdUcXaPdya3GG6cktb6deD1DMpnAbPyOp7saNCgAS+99BKDBw+2OhQhhJszDINx48axfft2Ro0aRXBwMJs2bWL48OFcv36dhg0bptp/5syZJCQkMH78eIoUKQLAwYMH0yWQjRs3UqRIEdq0aZOj453J7RJLfuTn58f7779vdRhCFBqhoaHpygYMGMCoUaOIiYnhoYceSlc/fPhwhg8fzuXLl+nXr1+6+ueee47HHnuMyMhIhgwZkq5+3LhxPPzww2itczXp7PLly1m3bh3Lli2zr0HfunVrdu3aBaRe937v3r2sWbOGoUOHEhQUhGEYFC9enEOHDqU77/r162nfvj1FixbN0fHOlC9uheUH8+fP580337Q6DCGEm5s/fz6dO3e2J5VkNWvWpESJEqnWdXrnnXfw8/Pj+eefB8ylhmvXrp3uVlZcXBy//PJLulZIdo93NmmxOMnEiROJj49n0qRJVociRIG3cePGTOt8fX2zrA8ICMiyvmrVqlnW56a1EhkZyalTpxg9Ov3cuRcvXqR+/frYbOaz32vXruWPP/5gzJgx2Gw2+3Ny1atX57///S+3b9+mWLFiAGzbto3Y2NhULTlHjnc2SSxO0q5dO44fP251GEIIN3bhwgXATG4pXbp0ifDwcHs/bVxcHLNmmV3Ks2fPZvbs2enOdfjwYUJCQgBzNFjz5s3tI1QdPd7ZJLE4SeXKlfn111+tDkMI4cZKly4NwJEjR2jdurW9/JNPPiE2NtY+Imzx4sVERkYyadKkdC2ko0ePMm3aNLTWhISEYBgGGzZs4LnnnrPv48jxriCJxUl8fX25cuUKe/bsoVGjRlaHI4RwQ0FBQdSuXZs5c+ZQvHhxAgMDWbNmDWvWrAHMjvurV68yd+5c2rdvz9ChQ9Odo169evbEALB7924uXbpE586dARw+3hWk895JypcvD8D27dstjkQI4a48PT35+OOPqV27NtOmTWPChAl4eHjw+OOP4+vrS1BQEJ9++imxsbFMmTIlw3OUKlWKwMBAe2JYv349DRs2JDAwEMDh413BVtjnt1JKBQHHw8LCqFKlSo7Pc+LECb7//nv69etnTzIi/0vuDM2qM1cIKz344IP07t071a2wvHD69OnkVlINrfWJlHVyK8xJgoKCGDVqlNVhCCEKmR9++MHqENKRW2FONGXKFJ566imrwxBCCEtJi8WJ5s6dK1PnCyEKPWmxOFHHjh0pW7as1WEIIYSlsp1YlFJdlVJ+rgwmv6tVqxbHjx8nISHB6lCEEMIyjrRY/g1MdFUgBUFAQAB3797lq6++sjoUIYSwjCOJxQacdVUgBUFwcDBgTpUghBCFlSOd95OAqUqpKOBH4KI7rDfvTjp06MCqVato27at1aEIIYRlHEkso4FSwKLkggxm+TS01oV2pFnp0qV5+OGHrQ5DCCEs5UgS2J70JbIwY8YMFi5ciNY61YI7QghRWGQ7sWitn3RlIAXF+vXruXDhAteuXaNChQpWhyOEEHnO4dtWSql6QC+gOhAHnALWaFcvSZZPdOvWjU2bNuHj42N1KEIIYQmHHpBUSr0D7APeBp7F7HeZCUQopd5zfnj5z3333QfA/v37LY5ECCGs4cgDkk8DE4DvgdaAP1AGaAOsBsYppYa5Isj8JDmxyISUQojCypFbYX8FftZa90pTvh3oq5QKS9rnC2cFlx9Vq1YNb29vvL29rQ5FCCEs4citMAWszKJ+JVAvd+Hkfzabjd27d7Np0yarQxFCCEs4klhuAoFZ1FcEbucunIKhbt26+Pj4EBMTY3UoQgiR5xxJLOuAF5RS6RZ0V0o1Bl4AfnJWYPnZ4cOH8fPz44UXXrA6FCGEyHOOTunSDfhdKbUOSB5eXBfoCkQBk50bXv50/fp1bt68KVPoCyEKpWy3WLTWp4AWwAqgAzA26asD8C3QUmt9zBVB5jchISF4e3uTmJhodShCCJHnst1iUUp1BXZorQcqpTyAAMwZjy9preUTNIWiRYvStGlTtmzZwpUrV6TlIoQoVHK0HovWOlFrfVFrfUGSSsYaNWrEjh07+Oyzz6wORQgh8pSsx+IiHTp0oESJEjRv3tzqUIQQIk/JeiwuMmDAAAYOHIjNZrM6FCGEyFOyHouLeHp6kpCQwI4dOyhfvjy1a9e2OiQhhMgTsh6LC3366aeMHTuWcePG8f7771sdjhBC5AlHEsu/MUeFXXdVMAVNkSJFMAyDxx57zOpQhBAiz+RoVJjInvvvvx+A3bt3WxyJEELkHRkV5kJ169YlMDCQhQsX8uOPP1odjhBC5AlHEssk4FWl1BClVAWllAx3ugebzUZoaCjh4eF88UWhXk1ACFGIyKgwF+vbty8eHh4sWLDA6lCEECJPyKgwFxswYAADBgywOgwhhMgz2U4sWusnXRlIQWYYBpMmTSIuLo6ZM2daHY4QQriUw7etlFIBwANANeAr4BYQoLU+4OTYCow333yTd999ly5dumAYhjyNL4Qo0BzpvEcpNQ44BSwF3gJqAm2A/UqpOdKhn7HQ0FASExN5+umnJakIIQq8bCcWpdQg4H3Mte37Yw4/BtiJuR7Ls5irSIo0WrVqRenSpVm9ejV37tyxOhwhhHApR1osLwE/aa2fADYmF2qtI7XW/YDvgZHODa9g8PLyonfv3qxYsYKAgACOHZP10IQQBZcjiaUe8J8s6ldj3hoTGejfvz8xMTG0adMGDw+H7kAKIUS+4kjn/U3AP4v66kB07sIpuLp06cJ3331Ht27d8PHxsTocIYRwGUf+dP4BGKWUKp+2QikVAjwP/OSswAqaokWL0rt3b3x8fDh69Cg7d+60OiQhhHAJRxLLK5gd9hHAfMAAnlFKrQTCgbvAVKdHWIDcunWL119/ndDQUP72t79ZHY4QQrhEthOL1voM0AxYA3TCTDL9ga6YfS8ttdbSK50Fb29v5s6dS82aNVm2bJnV4QghhEs49ICk1vocMDzpeZUAwBO4pLVOcEVwBY2XlxeDBw/mk08+oWjRolaHI4QQLpGj4Ulaa0NrfUlrfV6SimOGDRvG3bt3mTNnDgMGDODPP/+0OiQhhHAqGfeax0JCQmjSpAnffPMNW7dulcQihChwJLFY4Pnnn6dOnTrs27ePfv36WR2OEEI4lSQWC4wYMYJvv/2WMmXKALBv3z6LIxJCCOeRxGKhI0eOsGDBAkJCQti6davV4QghhFPItPkWOXXqFMHBwUyePJkPP/yQpk2bWh2SEEI4hUOJJWna/OmAD+YDkuFASWClUmou8FetteH0KAugatWq0b17d+bNm8fJkydlmhchRIEh0+ZbaNy4cVy8eJElS5awa9cuOnTowIULF6wOSwghckWmzbfQ/fffT5MmTXjvvffw8vLi/PnzREZGWh2WEELkikybbyGbzcbUqVOJjIzkzp07HDhwgGbNmlkdlhBC5IojiUWmzXeB3r17c/z4cZo3b46npyeJiYksXLiQmzdvWh2aEELkiEybbzGbzUZgYCCGYXD+/Hn27NnDyJEjWbx4sdWhCSFEjjgyKuwV4DfMafM3879p8/8K9ABuINPm59iYMWP4z3/+w6FDh9i2bRstWrSwOiQhhMgRmTbfTfTs2ZNTp06xcOFCWrZsic1m4+LFixw5csTq0IQQwiEybb6beOCBB2jfvj1vvPEGTzzxBCVLlqRbt254enoSHh6OzWa790mEEMINZDuxKKUOA0uBpVrrw8All0VVCNlsNj744ANatGjBtGnTmDVrFh999BH+/v6SVIQQ+YojnfdngMnAQaVUuFJqjFIq0EVxFUrNmjVjxIgRLF++nNu3b9OxY0caNWoEyESVQgjniY+Pd+n5HeljCQWqAGOBOOADIFIp9ZNSaphSqqRrQixc3nvvPfbv30+xYsXsZevXr6dRo0asWLHCwsiEEPlVykQyduxYmjdv7tLXc2h246QVIz/WWrcFamCOFCsFfA5cUEp97YIYC5XSpUtTsmRJ4uLi0FoDEBoayrvvvkuPHj0sjk4IkR+kTCQff/wxgYGBxMXFAeadkZ49e5KYmOiy18/xtPla61PAJ8AMzOlcfIBHnRRXoTd48GC6dOnC9evX8fLyYvz48fj4+BATE8O2bdusDk8I4UYSExNJSDDHUK1atYrSpUtz/PhxwFy1dtiwYcTExADmZ8v06dPx8HDdqikOn1kp5a2U6quU+hdmB/53QDAwDVBOjq/QGj9+PGfPnmXcuHGpyidOnEiXLl04f/68RZEJIaxmGAZ3794FYO/evZQrV44ff/wRgHr16jFkyBD7vqGhocyaNQt//6wmTnEuR2Y37pMimawA7se8BdZKa6201tO01vLQhZO0aNGCl19+mYULF/L999/by19//XX+9a9/ERgo4yaEKCwMw+DOnTsAREVFUblyZebMmQNA7dq16dOnD+XLm5Oi1KlTh88++4waNWpYFq8jLZaVmE/YrwS6AZW11n/TWoe7JDLBa6+9RoMGDRg5ciTXrl0DoEyZMvTu3RuAbdu28eWXX1oZohDCBQzDsP+fNwyDhg0bMnbsWAD8/f3p168f9erVA8DX15eFCxfyl7/8xbJ403LkAcmBwCqt9R1XBSNS8/b25osvvuC5557j6tWrlC5dOlX9Rx99xK5du+jfv78sFCZEPpaYmMjZs2epUqUKAF27diUhIYENGzZgs9kYOnRoqhbIxx9/bFWo2ZLtxKK1lhFfFvjLX/7Cjh07MnxI8ssvv+TSpUv4+PiQkJBAfHw83t7eFkQphHBEQkICWmvq168PwMiRI/nhhx84ffo0NpuNJ5980t4ZD/Dyyy9bFWqOZJpYlFIJwBCt9b+SthMxJ57MiqG1dmiaGHFvNpuNmzdvMnr0aMaNG0fDhg0BKFq0KJUrVwZg8uTJbNy4kbCwMHx9fa0MVwiRRlxcHDt37qRZs2Z4eXnx9ttvM3XqVK5evYq/vz9DhgyhQ4cOJCYm4unpyaBBg6wOOVeySgKLgaNptmU9e4vcunWLtWvXsn37dsLDwylRokSq+mbNmhEbGytJRQg3EB0dzcaNG2nZsiXlypVj5cqVPP7444SHh9OsWTMeffRRatWqRZEiRQBz5FaBYhhGof4KDg4OCg4ONiIjIw13t2HDBsPDw8MYPHiwkZiYmOl+x44dM1577TUjLi4uD6MrmDp27Gh07NjR6jCEm4uOjjYWLVpk7N+/3zAMw9i5c6cBGEuXLjUMwzAuXrxoLF++3IiKirIyTKeKjIw0goODjeDg4CAjzeeqI8ONNyilOmdR/7BSKsIp2U5k6P7772fatGksXbqUmTNnZrrfihUr+PDDDzl37lweRidE4REfH8+sWbNYu3YtYHa+P/XUU3z77beA+VDixo0b6dOnDwDlypXj0Ucfxc/Pz7KY85LNMDK+u6WU8sWcGj/ZCWAMGa977wFMAIZrrfPVvRilVBBwPCwszD4iw50lJiby+OOPs3HjRg4dOpTphXr69Gn77/Ptt9/Ss2dPe7NbZF/yLYqNGzdaGoew3htvvEHp0qV54YUXMAyDatWq0adPHz755BMADh8+TK1atVz6RLs7OX36NJ07dwaoobU+kbIuqz6W4sBuIPmTywA+SvrKiA1ZmtjlPDw8WLRoERcuXMjyr5/kpBIeHs4jjzzCnDlzGDVqVF6FKUS+FB8fj5eX+bE4ceJEzp49a18mfNu2bVSoUAEwB9T8+eeflCz5v7l369Spk/cBu6lME4vW+pJSajDQAjNpTAW+BfZmsHsC5hP5y1wRpEitWLFiBAUFYRgGs2bN4rHHHqNq1aoZ7tu8eXPWrl1Lp06dAHP6h4CAACpVqpSXIQvhdhITEzl58qT9+ZCXX36Z5cuXc/SoOWapWLFiqQbDfP/996mG/adMKiK1LIcGa63XAmsBlFLVgbla6x15EZi4txMnTjB9+nQWLlzIL7/8QkBAQIb7Pfjgg4A5UGPEiBHExsayZ88eWUBMFCpXr15l69atdO/eHS8vL2bMmMG0adO4ceMGxYsXp2XLlhQpUsTeannttddSHS//X7LPkfVYnrxXUlFKlct9SCK7atSowerVqzlx4gTdu3fn5s2bWe5vs9lYtmwZf//737HZbMTHx/Ppp58SHR2dRxELkXdOnDjBtGnT7INY1q5dS69evThw4AAAffr04R//+Id9/0ceeYQZM2bYb4WJnHPoHVRKDcGcGr8EqZOSF1ASaAAUdVp04p46dOjA119/Td++fenduzdr1qxJtUhYWrVq1aJWrVoAbNiwgRdeeIGqVava5x8TIr9JTEzEw8ODEydOMGbMGMaPH0+7du24dOkS06ZNo2XLllSsWJFu3bqxZcsWe19ISEgIISEhFkdfMDky3HgCsAhzAsr7gI5ALcw+mHaYU+e79wQ2BdTDDz/MokWL2Lp1K7/99lu2j+vatSs7d+6kV69eAMybN49x48a5fNlSIXLCMAyOHz/OqVOnALhy5QpVqlRh7ty5APj5+XHo0CH75I1NmjThxo0b9lvBAQEBtG3bVubVywOOjIt7EtgDlAdaY3bo3485aux5zIW+tjs7QJE9TzzxBEeOHKFjx46A+Z8wO5o0aWK/d3z48GF2795tvxVw4sSJbJ9HCFdYtGgRq1evBsyWSYMGDfjggw8Ac6bvhx56iNq1awPm6qsHDhzg4YcfBsDLyyvdDBUibziSWIKAxVrrm1rrY8A1oL3WOkFr/XfgK+BvLohRZFPyyLDly5fzwAMP3LPPJa2ZM2faFwuKjo6mcePGjB8/3ulxCpHMMIxU1+m4ceNSXXOzZs1i0aJFAHh6erJ06VL+7//+DzD7DOfPn0/Xrl3zNGZxb44klrtAyk+qw0DKG5Q/Y94OExaLj49n48aNhIaGcvHiRYeO9fT0BMy/9mbPnm2fDC8yMpIHH3yQ/fv3Oz1eUXjs3LmTr7/+30Tp/fr1s7eyAW7fvm1f0Arg559/5ptvvrFv9+3b1z4jsHBfjiSWA0CbFNsaaJZi2x+QOdvdwMCBA1m1ahUHDhygbdu2HDt2zOFz+Pj4MGzYMJo2bQrAsWPHOHDggH1c/86dO/n888+5ffu2U2MX+d/Vq1ftPy9dujTVwJCFCxcycuRI+y3WQYMG8cwzz9jrP/vsM/uT7GD2ixSWJ9kLEkf+xf4JPKmUWqKUKg6sAtorpV5TSg0AxmL2wQg38NBDD7FhwwauXr1KmzZtuHTpUq7O17FjR06cOEHNmjUB+Prrrxk9erS9fvv27fzxxx+5eg2R/xw4cIC3336bmJgYwLx1VbZsWW7cuAGYt1QvX75sr3/11Vc5ePCg/fhHH300VWIRBYMjz7HMBd4CemLeFluJ+aT9a0nffYH8tRpNAdeqVSu2bNnCqFGjKFcu948YpXxA7O2332b//v32oc1Tpkxh2LBh9vqwsLBUHyAif7p58yZhYWH2kVYbNmygZs2a9n/bvXv38uqrr9pbxZ06dWLWrFn2FskzzzzD1q1b7S3dypUrU7FiRXnYsIBzqI2ptZ4MBGit47TWhtZ6EBAKPAIEa623uSBGkQv16tVj6tSpgPkhMHXqVBITE3N9XpvNRlBQkH176dKlfPHFF/btp59+msmTJ9u3Fy9eTESETH7tbuLi4vjjjz+4cOECAEePHqVr16788ssvAOzbt48uXbrw66+/AlChQgVatmxpP75Xr17cuHHDvvhckyZNePHFFwvNLL4iYw7fvNRax6fZ3qy1/k5r7VgvschzK1asYPr06fTp08fhEWP3Ur58ef7yl7/Yt9evX8/06dMBiImJYcSIEfz73/8GzGVZhw8fzoYNG5wag0gvISGB3377jePHjwNw48YN+vfvz3fffQfAxYsXadasmX3b19eX69ev2/vOQkJCCAsLo127dgA0aNCAf//739StWxcw59OSObNEWlktTZyT//WG1jrTNVuEtV5//XUAW+eIAAAZPElEQVTKly/PmDFjaNmyJV999RX33XefS14r+el+MD+sIiMj7bdHzp8/T1hYGPfffz8Ap06dolmzZixYsIDevXtz/fp1tmzZQuvWrSlTpoxL4svPDMMgPj7evgzCd999R2BgIK1atcIwDDp16kSvXr0YO3YsAG3atGHixInMmDEDX19fIiIi7Le2KlWqxMqVK2nRogUAFStWZMeO/83cVKJECfsEpkJkV1YtlppADQe/aroyWJE7NpuN559/nh9//JGrV6/SokWLVB8irhQYGEjFihUB8z57ZGQkQ4cOBcwPyt69e1O9enXAHHHWs2dP+2CAHTt20Lp1a/bt2weY60CsXLmS69evA+Zf5fn5Qc7ExER75zbAr7/+ytatW+3bU6ZM4cMPP7Rv169fnxEjRti3R48ezd///nfA/DcuV66cvU/D09OTNWvW8PTTTwPmMPI///yTJ598EjCXYejbty+VK1d23S8oCp+0S0oWtq/8tDSxM50/f9546aWX3HL54ujoaGPr1q3G9evXDcMwjK1btxqdOnUyjh07ZhiGYSxZssQAjIMHDxqGYRiLFy82vL297fVr1641BgwYYFy9etUwDMMIDw83PvjgA+P27duGYRjGkSNHjB9//NGIj483DMN8LyIiIuzLPUdFRRmnT582DMNcmrht27bGmTNn7PGdP3/e0Frbtw8ePGhs3brVvr1lyxZjxYoV9u1//vOfxttvv23fnjhxojFy5Ej79v3332906NDBvt22bVujU6dO9u0ePXoYTz/9tH37o48+Mr755hv79qFDh4wrV67c830VwpmyWpo4px/GlYKDg1sGBwf7BQcHFw0ODvbIyXnc4auwJpaULl68aLRv39749ddfrQ4lW27cuGHs2rXLiI2NNQzDTBwTJkwwoqOjDcMwE0/dunXtiemdd94xACMmJsYwDMOYPn26ARh37941DMMwpkyZYgD2xDJhwgTDx8fHMAwzsVSuXNkoVaqU/fWfffZZo1y5cvbtp556yqhSpYp9+4knnjBq1Khh3x4yZIjRokUL+/Yrr7xijBgxwr69dOlSY/HixfbtiIgI48iRI7l5i4RwuawSS6ZLE2dEKdUWc6LJxklFD2D203wOvKi1/jqzY91Vflua2BX27t1Lr169iIyM5MUXX+SNN97Icobk/CYhIYHo6GhKlSqFzWbj9OnTnDhxwt4hHRERwZ9//kn//v0B+O2334iIiODJJ58kNDSUGzduMHr0aIYPHw6Yq3KePn2avn37AuazHFFRUbRu3Row+5ASEhLk9pIo0LJamjjbiUUp1RzYDEQCqzHnBXsAuAF8A1QFHk5aHCzfkMRiunnzJuPHj2fevHkopfjnP/9p/6AszGTNeyEyllVicWS48QzgONAIeBtzdmO01r8nlR0AXnVCvMICJUuWZO7cufz000/cvn2b9957z+qQhBD5lCOJpTXwT631bSBVM0drfQOYDzR0YmzCAl26dGHfvn3Mnz8fgEOHDrFixQocuWUqhCjcHH1AMjaLOp8cnE+4oVKlStmngJk9ezb9+vXjoYce4siRIxZHJoTIDxxJBDuAQRlVJE1K+TQQ7oyghPuYPXs2H330EVu3bqVhw4ZMmjTJPsGgEEJkxJHEMhVoopTaBAzDvB3WUik1GnNW45rAm84PUVjJy8uLMWPGoLWmf//+vPXWW8yaNcvqsIQQbsyR2Y23Yc5sXAWYidl5/ybwEVAMeExr/bMrghTWq1ixIl9++SXh4eH2qUI2b97MqlWrpP9FCJFKthOLUqqs1vonoDbmAl+PYd4aawNU11qvdE2Iwp00a9YMf39/wLxN1rt3b1q3bs2PP/4oCUYIAWQxCWUGdimlFmitpwM7k75EIbZs2TIWLVrEjBkz6NatG23btuWdd96xP3gohCicHOljKQecd1UgIv8pUqQII0eO5NChQ3z22WecOHGCP//8E4C7d+86Zd0XIUT+40hiWQr8X9KT6kLYeXt789xzz3HkyBH7rLnz5s0jJCSEL774gri4OIsjFELkJUduhSUCdYHDSqkjwEUgIc0+sh5LIebj42P/uVq1anh4eDB8+HAmT57M2LFjGTlypCwKJUQh4EiL5QHgMnAa82HIash6LCITvXr1Ys+ePaxdu5Y6deowbtw4HnvsMavDEkLkgWy3WLTWNVwZiCh4bDYbDz74IA8++CDh4eEkJJgN3AsXLvDEE0/wzDPP0Lt3b/tKiEKIgkGmYBF5onnz5rRq1QqAo0ePcvjwYfr3709QUBCvv/66fU12IUT+J4lF5Lk2bdpw9OhRVq9eTUhICG+88QZ16tTh0qVLVocmhHACRzrvhXAaT09PevbsSc+ePTl58iSbN2+2T3z52GOPUbx4cYYPH067du3w8JC/f4TIT+R/rLBc9erVGTJkCGAule3v788333xDx44dqV69Oi+++CJ79+61OEohRHZJYhFuxWazMW/ePM6fP8/SpUtp2rQpc+bM4YcffgDg1q1b7Nq1S6aPEcKNya0w4ZaKFy/OoEGDGDRoEFFRUfbyVatWMWjQIKpXr06vXr3o3bs3HTp0kJFlQrgRabEIt+fv72+f+LJr16784x//ICQkhAULFtClSxfKlSvHuXPnAKQlI4QbkBaLyFfKli3LiBEjGDFiBDExMaxfv56tW7cSGBgIwFNPPcWpU6fo2rUrXbt2pVGjRtL5L0Qek/9xIt/y9fWlV69evPvuu9hsNgD7sOWJEyfStGlTAgMDefXVVy2OVIjCRRKLKFBeffVV9u7dy5kzZ/jiiy/o2rWrfQ6z+Ph42rZty5gxY1i5cqU8NyOEi0hiEQVSpUqVGDp0KEuWLGHq1KkAXL58mWLFirFgwQIeffRRypcvT4MGDVi9ejUg/TNCOIv0sYhCIzAwkPXr1xMXF8fvv//O5s2b2bRpE6VKlQJgw4YNPPvss7Rq1YpWrVrRsmVLDMOw32YTQmSPJBZR6BQtWpQ2bdrQpk0bJk6caC8vVqwYDRs2ZP369SxZsgQwn6tp3rw5AKdOnQKgatWqkmyEyILcChMiSZs2bfj22285e/YsJ0+e5KuvvqJy5cr2PpqZM2dSvXp1KlWqRI8ePZg8eTLLly+XW2hCpCEtFiHSsNlsVKtWjWrVqvHZZ5/Zy5999lmUUuzYsYNdu3axbt06KlSoQL9+/QCYNm0a165do0mTJjRu3Jh69epRtGhRq34NISwjiUWIbKpfvz7169fn+eefB+D27ducPn3aXn/gwAFWr15NTEwMAF5eXjz66KMsW7YMgJ9++olKlSpRp04dSTiiQJPEIkQOFStWjDp16ti3ly1bRkJCAocPH2bXrl3s37+fihUrAuaIs0ceeYTo6Gi8vLwIDg6mQYMG9O/fn/79+wMQFRVln2FAiPxMEosQTuTp6UndunWpW7duqnLDMNi8eTN//vknERERREREsHPnTurVq0f//v2JioqidOnSlC1bluDgYOrUqUNwcDA9evSgcePGFv02QuSMJBYh8oCHhwdNmjShSZMmqcqTO/5tNhszZ87k0KFDHD58mLCwMBYvXkxAQACNGzdm3759dOnShRo1ahAUFGT/3r17d6pVq2bFryREpiSxCGGh5GHLfn5+jBs3LlXdrVu37D/7+PjQq1cvTpw4we+//87KlSu5e/cua9asoVq1aqxdu5ann36aoKAggoKCqFKlClWqVGHAgAFUqFCBuLg4PD098fT0zNPfTxROkliEcFPFixe3/1ynTh0WLFhg305ISODs2bOUKVMGgICAALp168bx48f59ddfOXv2LHFxcXTs2JEKFSqwcOFCRo8eTaVKlahcubI98UyaNImyZcty9uxZbt68SYUKFfDz85PndESuSGIRIh/y9PSkatWq9u3mzZvbH+QESExM5PLly5QuXRqApk2bMmHCBE6fPs2ZM2fYu3cva9asYcqUKQDMmTOHt956CwBvb28qVKhAYGAgmzZtwsfHh3Xr1nH48GECAwPtdeXLl8fPzy8Pf2uRX0hiEaIA8vDwoHz58vbtli1b0rJly1T7pHywc/DgwdSvX5/z589z4cIFLly4wNWrV/H29gbgX//6F4sXL051fKlSpbh+/ToAEyZMYMeOHZQtW9b+Vb16dZ577jkAtNaA2bLy9/eXW3IFnCQWIQqplLe7kp/Rycznn3/O+++/nyrx3L59215fokQJbDYbhw8fZvv27Vy+fJlatWrZE8vIkSP55Zdf7K9bunRp2rZty6pVqwB4+eWXuXTpEv7+/vj5+eHv709wcDDdu3cH4ODBg3h7e+Pv70+pUqUkMbk5SSxCiHvy9PSkfPnyqVpBKU2dOtU+izSYraE7d+7Yt99++21OnDjBlStXuHz5MleuXKFSpUr2+oiICPbs2UNUVBTR0dEAdO/e3Z5YunTpwpkzZ+z7lyhRgkGDBjFv3jwAevfubS8vWbIkJUqUoF27dvTp0wcwnzEqXrx4qvry5cvb+6iEc0liEUI4nc1mo1ixYvbttm3b0rZt20z3/+9//2v/OT4+nhs3bpCQkGAvmzt3LpcvXyYqKorr168TFRVFw4YN7fU3btzg2rVr3Lx5k+joaG7evMmdO3fo06cPd+/e5fHHH0/3muPHj+e9994jKiqKChUq4Ovri6+vL8WKFcPX15cxY8YwYsQILl++zLPPPpuqrlixYvTq1YtWrVpx7do1Vq1alarO19eXOnXqEBAQQGxsLJcvX8bb29v+VaRIkQI9QEISixDCrXh5eaVrSfTs2TPLY37++ed0Zcl9SJ6enhw8eJDo6Gh70omOjiY4ONj+emPHjiUmJoaYmBhu375NTEyMfWDC7du3OXDgQKq6mJgYqlSpQqtWrTh58iTDhw9P9/pffPEFQ4cOJTw8nPbt26eqs9lsrFy5kj59+vDzzz8zdOjQVInHx8eH2bNn07x5c3799Vc++uijVHXe3t6MGTOGatWqsWfPHtatW0eRIkUoWrSo/XufPn3w9/fn2LFjHDhwwF6XXN+oUSOXTS0kiUUIUSAltwg8PDxQSmW6X4kSJXjnnXcyra9atSoRERHpypMTV/369Tl27Fi6xJTcoqpVqxbz58/nzp07xMbG2r+SYypTpgxdunRJVRcbG0uRIkUAuHbtGvv27UtXP3jwYKpVq8b27dt5+eWX08XXunVr/P39+fbbb3nppZfS1Z85cybV7UhnshX2Kb+VUkHA8bCwMKpUqWJ1OMLNhIaGArBx40ZL4xAiMwkJCdy5c4e7d+8SFxfH3bt3uXv3LpUrV6ZIkSJcuHCBkydPpqqPi4vjgQcesI/6y4nTp0/TuXNngBpa6xMp66TFIoQQ+Zinp2eqh2nTqlChAhUqVMjDiGShLyGEEE4miUUIIYRTSWIRQgjhVJJYhBBCOJUkFiGEEE4liUUIIYRTSWIRQgjhVJJYhBBCOJUkFiGEEE4liUUIIYRTSWIRQgjhVJJYhBBCOJUkFiGEEE4liUUIIYRTSWIRQgjhVJJYhBBCOJUkFiGEEE5VoFeQVEq9CIwADOB34Bmtday1UQkhRMFWYFssSqmWwJNAC+A+wBMYbWlQQghRCBTkFss14K9a61sASqm9QDVrQxJCiIKvwCYWrfUh4BCAUqoiMAYYbmVMQghRGOT7xKKUGgh8mqb4gta6QVJ9dWAt8Het9fq8jk8IIQqbfJ9YtNbLgGUZ1SmlGgOrgfe01p/kaWBCCFFI5fvEkhmlVDngB8x+luVWxyOEEIWF2yQWpVQvYKnWumSa8pHABKAKsBt4UWu9LRunHAuUAiYrpSYnla3RWk9yYthCCCHScIvEopRqAywBbGnKhwJzgTeAcOAFYJ1SqpHW+nhW59Ravwq8mo2X9wQ4f/58DiIXBV18fDwAp0+ftjgSIdxLis9Mz7R1NsMw8jaaFJRS3pijtaYDt4CiWusSSXU24DiwVmv9XFJZEUAD/9VaO+WZFKVUO+AXZ5xLCCEKofZa6y0pC6xusXQHXgHGA2WBcSnqagPVgVXJBVrru0qpNcCDTowhHGgPnAMSnHheIYQoyDyBipifoalYnVjCgRpa6yil1Otp6oKTvh9JU34MqKWU8tRa5zoRJE3xsuWeOwohhEjraEaFliYWrfWZLKpLJX2/mab8JuZUNMWBG66ISwghRM5Z3WLJSnJHftpOoOTyRFe+uKOj0ZRSDYHZQEvgKjAH8/kZ6zqx8lgO3rP/Aj0yqCqptY52TZTuJ7MRkRnsV+ivsWQOvGeF9hpTSnli9mGPxJzO6iTwGTAns2vGWdeYO09CeT3pe9oLpwRmUrnlqhdOMRptCfAoEIU5Gq1GJvuXB9ZjJsEBwHzgTVL3GRVojr5nSUIwL+LWab5iXBut+8hsRGQG+xX6ayxZdt+zJIX5GpsCvIX5XvUCvgY+wuzTTseZ15g7t1gOJ32vSep+lpqAdtVfaUmj0d4A5mutpyWV/YQ5Gm0sGc+Q/Dzme9lLax0DfJ804u0VpdRsrfVdV8TqLnLyniml/IGqwA9a6+15GK5byGhE5D0OKdTXGDj+nhXma0wp5QG8CLyvtX4zqTgs6cHxl4D3MjjMadeYO7dYDgORQJ/kgqThxj2AMBe+boaj0YCsRqN1AcKS/jGSfQeUAZq7KE53kpP3LCTp+17Xhua2Uo6IzM50Q4X9GgPH37PCfI35AYuBlWnKNVBOKVU8g2Ocdo25bWJJapG8AzyrlHpTKfUQ8B8gAPjQhS99z9FomRyT0f4pz1eQ5eQ9CwFigRlKqStKqRil1DdKqUBXBupGkkdEfkz6fsSMFPZrDBx/zwrtNaa1vqa1/qvWeleaqoeB08nLiaThtGvMbRMLgNb6M8y/ToYAywF/oJvW+liWB+ZOdkajZXRMRvunPF9BlpP3LATwTtqnLzAK8973hqTmd4GmtT6jtY5y4JDCfo3l5D0r1NdYWkqppzFbJRndBgMnXmNu08eitX4deD2D8lnArDwMJSej0WwZ7J/MpaPX3ERO3rMPgH9rrX9O2t6slDoAbMfsOPzS6VHmb4X9GssJucaSKKUGYw6uWU76ZUaSOe0ac+sWi0VyMhrtegb7l0xRV9A5/J5prQ+m+A+fXLYDczRZI1cEmc8V9mvMYXKNmZRSYzGT6H+BwVkMfHLaNSaJJb2Uo9FSymo02uFM9gezs6ygc/g9U0oNVEp1SFNmw7x1cdklUeZvhf0ac5hcY6CUeguz5fYl0E9rHZfF7k67xiSxpJeT0WhhQJc0Iy36AFcwHxQs6HLynj0HzE4aFpnsIaAYsNlFceZnhf0ay4lCfY0ppcZgjqKbDQzXWsff4xCnXWNu08fiLrTWhlLqHeBTpdQ1YCvwV1KMRlNK1QLKpRgb/xnmlP7fK6Xex2xmvwJMvMdfCAVCDt+ztzCXjF6ilPon5qiT6cAKrfWvef07uBu5xhwn19j/KKUqAu8C+zBX2G2plEq5y++Yjwi45BqTFksGsjEabQqwLcX+5zBHW3gl7f9/wCSt9cy8jNtKOXjP1mE+DVwbc6z8JODzpOOFXGM5IdfY/3TDvOV3H+Z7kvbLHxdeY5auxyKEEKLgkRaLEEIIp5LEIoQQwqkksQghhHAqSSxCCCGcShKLEEIIp5LEIoQQwqkksQghhHAqSSxCCCGcSqZ0ESIPKKVKA+cwn4YeorVeksW+TwIlkxa0yqh+DuYcTlVctUS3ELkhLRYh8sZgzKQSC4y4x77vAQ9kVJE0O29v4D+SVIS7ksQiRN4YgbnM69+BjkkTJqajlKqNOXnn9ozqMdcer4w595UQbkluhQnhYkqppkBjYBrwH+BvwJPA5DT7fYfZGgFznfYZyT9rrack/dwXc9Gln5OO6Qb8cI8Q/tBaN8vt7yFEdkliEcL1km99LdZaH1NK7QOGKaWmaq1TLvk6H/AEemKuJRKdVL4txT59gTVa67tJ2ydJPVvve8BdzJl8kx1DiDwksxsL4UJKKR/gLBChtW6fVPYS8D7wkNZ6bZr9/wO01lqXz+Bc9YA/gf5a6+UZ1HthLgP9ndb6Maf/MkJkk/SxCOFajwClgcUpypYCCcBTGezfFNiVybn6YHb+Z3brqx5QFNiTo0iFcBJJLEK41gjgDvB1ckHSgkrrgV5KqYDk8qSfqwA7MzlXX+AnrXV0JvWNk77LUsXCUpJYhHARpVQN4H5gldb6eprqxZitiydSlP0l6Xu6xKKUqgI0I+vRYI2SvkuLRVhKEosQrvMkYCP1bbBk3wI3Sf1MS5Ok7xm1WPoABrAqi9drAERprc84HqoQziOJRQgXUEp5AMOBC8C6tPVa69uY64o3VEo1TyqumfT9VAan7Ats0VpfyuJlKwIXcxqzEM4iw42FcI2uQFXM/o6XlFIZ7eOd9H0EEM7/hgV/rJTahtnB/y/AH+gATLjHa54DOiulJgN7tNarc/UbCJFDMtxYCBdQSn0D9Mvm7tcxWxs2YC7QHfPp+1Na6+pKqaHAF0BNrfXxLF6zKfA50BBYp7XukYtfQYgck8QihJtTSn0L1NBaN77nzkK4AbkVJoT72wYssDoIIbJLWixCCCGcSkaFCSGEcCpJLEIIIZxKEosQQginksQihBDCqSSxCCGEcCpJLEIIIZxKEosQQgin+n+FxsYQIK2chQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(delta_t_list,dtau,\"k:\",label=r'$d\\tau/\\tau$')\n",
    "plt.semilogy(delta_t_list,dA,\"k--\",label=r'$dA/A$')\n",
    "plt.axvline(x=0.7968,color='k')\n",
    "plt.yticks([1e-1,1e-2])\n",
    "plt.xlabel(r'$\\Delta t/\\tau$')\n",
    "plt.ylabel('relative error')\n",
    "plt.legend()\n",
    "plt.savefig(\"optimaldt.png\",format='png',dpi=300,bbox_inches='tight',facecolor=\"white\",backgroundcolor=\"white\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: 0.007859125189659028\n",
       "    nfev: 15\n",
       "     nit: 11\n",
       " success: True\n",
       "       x: 0.7968052720351725"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimize_scalar(OUlargeNdtau,bracket=(0.1,0.7,2.0),args=(100000,1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a241a5ef0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD6CAYAAACoCZCsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4HOW58OHfNmlVVr1bkotkj9x7xw7YgMEm1AQ4cSAFCAQChARygJATSIAQEpJAwMnh5JB89CRAQjOEgyEGV9xwwfK4S5as3utKW74/3pUsybItySuNpHnu69pr7NnR7uPx6H3mrWPx+/0IIYQwH6vRAQghhDCGJAAhhDApSQBCCGFSkgCEEMKkJAEIIYRJ2Y0OoKc0TQsFZgNFgNfgcIQQYqiwAanAFl3X3R3fGDIJAFX4f2p0EEIIMUQtAtZ13DGUEkARwEsvvURKSorRsQghxJBQXFzMypUrIVCGdjSUEoAXICUlhfT0dKNjEUKIoeakpnPpBBZCCJOSBCCEECYlCUAIIUxKEoAQQpiUJAAhhDApSQBCCGFS5kgAW/8ML1xhdBRCCDGomCMBVB2Bo+uNjkIIIQYVcySA0CjwusHjPvOxQghhEuZIAM5otW2uNTYOIYQYRMyRAEKj1NYtCUAIIdqYIwE4AwmgucbYOIQQYhAxRwKQGoAQQpzEHAmgvQYgCUAIIdqYIwFIDUAIIU5ijgQgNQAhhDiJORJAiEttpQYghBDtzJEAbHZwREgNQAghOjBHAgDVDOSWYaBCCNHGPAkgNEpqAEII0YF5EoAzSvoAhBCiA/MkAKkBCCFEJ+ZJAFIDEEKITsyTAKQGIIQQndiNDqArTdOigfXA5bquHwzaB0sNQAghOhlUNQBN0xaiCv9xQf/w0GjwNIOnJegfLYQQQ9GgSgDAzYHX8aB/slPWAxJCiI4GVROQruvXA2iaFvwPD+3wTICIhOB/vhBCDDEDngA0TbsWeLrL7hJd1yf26xdLDUAIIToZ8ASg6/qrwKsD/b0nagCSAIQQAgZfH0D/kRqAEEJ00ucagKZplwIv6bru6rL/JuBHQDrwOfADXdc3nlWUwSA1ACGE6KRPNQBN0xYALwKWLvuvB/4YeO8qoBr4l6Zpo3vz+bqujwrqHAAAZ7TaSg1ACCGAXtYANE0LBe4Efg40ACEd3rMAPwOe1XX9ocC+/wN04C7gjiDF3DfOGLDYoKHM0DCEEGKw6G0N4GLgPuAe4Pdd3ssGRgJvte3Qdb0VeBe46CxiDA6rFSKToa7Y6EiEEGJQ6G0C2AKM1nX9KcDf5b222btdm24OA1maptn6EF9wuVIkAQghRECvmoB0XS88zduBXlbquuyvQyWaCMDYBnhXKlQdNTQEIYQYLII5DLStQ7hrzaBtvy+I39U3rmSoKzI6CiGEGBSCmQDaHrjr6rI/ElX4NwTxu/rGlQpNleBxGx2JEEIYLpgJ4EBgO6bL/jGArut615rBwHOlqK30AwghRNATwDHg8rYdmqY5gBXAmiB+T9+5UtW2vsTYOIQQYhAI2lpAuq77NU17DHha07Qq1Lr+3wMSgN8G63vOSmSy2ko/gBBCBHctIF3XV6HmCFwHvAbEAMt0XT8czO/ps7YagDQBCSFE32sAuq4/CDzYzf4ngCf6HlI/Co8Hq11qAEIIgZlWA4XAbGCZDCaEEGC2BAAyG1gIIQIkAQghhEmZNAEE/5nzQggx1JgvAUSNUA+GlwfDCCFMznwJID5LbSsHx8hUIYQwivkSQFxbAjhkbBxCCGEwEyaAwNMppQYghDA58yWAkAg1I7hCEoAQwtzMlwBANQNJE5AQwuRMmgBGSxOQEML0zJkA4rOgoUyGggohTM2cCUBGAgkhhEkTgMwFEEIIkyaA2MBQ0AqpAQghzMucCSAkHGJHQckXRkcihBCGMWcCAEiZDMW7jY5CCCEMY+IEMEX1AbjrjY5ECCEMYd4EkDwJ8EPpXqMjEUIIQ5g3AaRMVtviXcbGIYQQBjFFAnhr53G+/+qOzjuj08EZI/0AQgjTMkUC2FNYw3t7ujwG0mKRjmAhhKmZIgG4Qu24PT5aPL7Ob6RMgZK94PUYE5gQQhjIFAkg0mkHoMHdpaBPmw6eJiiV+QBCCPMxRwIIVQmgrrlLAsicq7b5mwc4IiGEMJ4pEoDL6QCgzt3a+Y3oDHClwbFNBkQlhBDGMkkCUDWA+q41AItF1QKkBiCEMCFTJIBTNgEBZMyD2gKoKRjgqIQQwlimSADtNYCuncDQoR9AmoGEEOZiigTQNgqorrsEkDwZHBGQt2GAoxJCCGOZIgG4QgOdwM2tJ79ps8OohXD44wGOSgghjGWKBOB0WLFbLSd3ArfJWqpWBq08MrCBCSGEgUyRACwWC5FOe/d9AADZS9X20EcDF5QQQhjMFAkA1EigbkcBAcRnQ3SmJAAhhKmYJgG4nI5TJwCLBbLOg8NrwdtNP4EQQgxD5kkAoXbqu84E7mjshdBSB3nrBy4oIYQwkGkSQKTzNE1AAFlLwBEOuW8PXFBCCGEg0yQA1+k6gQFCwiH7fMh9B3y+Ux8nhBDDhGkSQGSo/dTDQNtMuAzqi6Fgy8AEJYQQBjJPAjhTExCofgBbCOx9c2CCEkIIA5kmAUQ5HbR4fbg93lMf5IyC7Atgz2vylDAhxLBnmgTQtiLoGZuBpv0H1JfI0hBCiGHPdAmgR81AYbGw85UBiEoIIYxjmgRw2iWhO7KHwqSrYN+70FQ9AJEJIYQxTJMA2peEPlMNAGD618HTDDtf7eeohBDCOKZJAKddErqrtOkwYiZs+RP4/f0cmRBCGMM8CaCnTUBtZt8EFQfg8L/7LyghhDCQaRJAZG8TwMQrIDweNv93P0YlhBDGMU0CiA5zYLFARX1Lz37A4YTZN8L+96B0X/8GJ4QQBjBNAnDYrCREhlJS29zzH5pzM9jDYP2T/ReYEEIYxDQJACA12klRTS8SQEQ8zPwG7P4bVOf3X2BCCGEAUyWAlCgnRTVNvfuhBbeDxQqf/Kp/ghJCCIOYKgH0ugYAEJ0OM78FO16CikP9E5gQQhjAXAkgJoy6Zk/PRwK1WfRDtUrox4/2T2BCCGEAcyWAaCcAxb2tBbiSYf6tapXQgq39EJkQQgw8UyWAlKg+JgCAc+6CyGR4/z6ZHSyEGBZMlQBSo8MAON7bjmCAUBcs+QkUfAafvxzkyIQQYuCZKgEkR4cCfawBAExbCRnz4F/3Q11JECMTQoiBZ6oEEGq3kRAZ0vuRQG2sVrjsaWhtgvfuCW5wQggxwEyVAABSop0U96UJqE3CWDj3XvXc4L1vBS8wIYQYYOZLAFFhfa8BtFlwO6RMgXd/CI2VwQlMCCEGmOkSQJ8mg3Vlc8Blz0BTFbx5m4wKEkIMSaZLAEmuUGqaWnF7vGf3QalT4MKfg74aNjwVnOCEEGIAmS4BxEWGAFDV0IMng53J3Ftg/KXw4UOQt+HsP08IIQaQ6RJAfIRKABUN7rP/MItFNQXFjoK/fwvqS8/+M4UQYoCYLgHERai5AJUNPXwwzJk4o+Dq56G5Bl6+BloagvO5QgjRz0yYAFQNIGgJACBlEnzlf6Hoc3jtBvD2crE5IYToqrUZjnwCHz2iWhgaKoL+Ffagf+Ig194E1NNHQ/ZUzgq4+HFYfTe89yNY8YRqIhJCiJ7wtsLxHXBkrSr48zeD162eR5IxDwj+aEPTJYDoMAc2qyW4NYA2c26CmgJY/zv1QPnz7pckIITons8HJbtVYX/kEzWQpKVevZc8WT2TfPRiGDkfnNH9EoLpEoDVaiE23EFFfyQAgKU/hcZy+ORx9XdJAkIIUAV+6RdwdD0c/RSOroPmavVewjiYem2gwD9HPY52AJguAYDqB6gMxiig7lit8OXfBx4j+Tj4fbDkAUkCQpiNzwvFu1SBn7de3eG3FfgxmTD+Ehj9JRi1CKJSDQnRxAmgn2oAoJLAJU8CFvj019DaCBc+ovYLIYYnbysU7VR39nnrIX8TuGvVe3FjYPyXYdQ5MHIhxGQYG2uAKRNAfEQoucW1/fslVitc8jtwhMGmVapv4Mpn1d+FEEOfpwWOb+9Q4G+G1sAw8IRxMOmqQIG/AKLSjI31FEyZAPq9BtDGaoWLHoPoDPjgx/B8CVz7yoC17wkhgqilAQq3qTv7o+vg2GfgCawsnDQBpn0NRi1Ud/iRScbG2kOmTQDVja14vD7stn5ulrFYYMH3IHoEvHEzPHsuXPM8pE3v3+8VQpyduhI4tkkV+PmbVHu+zwNY1Nyfmd9UBX7mgiF7U2fKBBDfth5QYyuJrtCB+dKJV0B0JvztevjfZbD8VzDzGwPz3UKI0/P5oHx/5wK/6oh6z+6EEbNg4Z2QOR/SZ0NYjLHxBokpE0DH2cADlgAA0mfCzZ/A6zfA23dA/ka4+Jf9NsZXCHEKrc1q0lVbgX9ss1reHSA8ATLnwewbVIGfMgXsIcbG209MnQDUgnCugf3yiHj4+uuw9nE1TPToerjij6oqKYToH7VFULhVtdsf26wKf2+gHzB+LORcogr9zPlqxI5Jhm2bMgHEB3tBuN6y2uC8+yB7KbzxHfjLCph3q5o0FhppTExCDBetTVC0Cwq2BF5bobZAvWd1qP63uTerwj5jLkQkGBuvgUyZAPplQbi+yJgDt6yDDx6ATc9A7luw/NegXWRsXEIMFX6/aqsv2HqiwC/eHeisRfW7ZcyB9NtU233KZHA4jY15EDFlAogNd2CxQFldP80G7o3QSPjy72DKNfDO9+GVa1R19IKfQXyW0dEJMbg0Vavmm7YCv3ArNAZWyXREwIgZ6pnd6bNVx60r2dh4BzlTJgC7zcro+Aj04jqjQzlh5Hy4+VPY+DR88mvYP1dVUxffDWGxRkcnxMBrrlVDL4/vOPGqPHzi/QQNxl0M6bNUgZ+YAzZTFml9ZtqzNSEtih351UaH0Zk9BBb9QE0o+ehh2PgM7HhB9Q/MvWXYDD0T4iTuetV007GwrzhI+xLI0RmQNg2mf1214afNkN+HIDBtApg0Ipp3dhVR3dhCTPggG+LlSoHLnlY1gH8/Bv/+BWxcBfNugXnflRqBGNpaGqFkDxz//ERhX66rhRMBXGmqkJ9ytdqmToPIRGNjHqZMmwAmpkUBsPd4LQuyB+kogJTJcO1LakTDJ4/D2l/Cpj+oCWSzb1TPIhZiMKsrUXf2JbvVtni3urNvK+wjklQhP+GywJ39NHUDJAaEiROAmny153jN4E0AbVKnwDUvQvEe+PQJVRvY8DRoy2Hud9SSsiYZtywGKZ9XFezFuzu/GkpPHBOdqW5qJl6ptiNmgCtVrl0DmTYBxEWEkBrt5Ivj/bwqaDClTIKv/hlqCmHrc7Dtz6C/qzrDpn1NjSQyaF1xYSL1ZVCWC6W5UPKFas4p2XtiYTSrA5JyYOwFqqBPmQzJE6XpchAybQIAVQsYUgmgTfQIWPoTWHwP7Hkdtv0FPvwprHkIxpynkoG2HELCjY5UDGVNVVC670Rh3/ZqLD9xTFgsJE+CWd8+UdgnjBu2SycMNyZPAFF8tK+EereHyNAheCocTpi+Ur3KD8LOV2DXX9VaQ/YwGHs+jL8Uxi2T9YbEqTXXQPmBEwV8W4FfV3TimJBISBoP2sVq6eOkHLWNTJYmnCFsCJZ6wTN3dBxP+uGzIxUsyRniE0YSslWt4LwfQ9462Psm5L4DuW+rKvmYL0H2BWr5ifhs+aU1G68HqvNUQV9xILA9qLYd2+ntYZA4Dsacq8bVtxX20RlyzQxDpk4AM0bGEmq3su7AMEgAbaxW9WDp0Yvh4l+pmZK5b8G+d+H9/1THRGdC9hL1S545X0ZdDBc+n7prrzqqJkxVHFA1w4oDUHkEfK0njg2PV4ugjbtQbRPGqgI/dpRaq0qYgqkTgNNhY/aoONYfLD/zwUOR1arWQcmYAxc+rAqBQ2vg4EewO9B3AOqXvm1hrMz5qg1Xnl88OLU0qjv5qqPq/7PqqFoLp+ooVOWBt8PyJrYQtbJlwjjIWXGioI/PhvA4g/4BYjAxdQIAWJidwC/f30dpXTNJrmG+SFTcaIi7Uc0h8Laqafb5m9RzCQ5+qPoQAEKjT3TopU5R66EnamBzGBv/cOf3q/b4mgKoLYSaY2rEV00BVOerQr6+uPPPhLggbpT6/xm3DGJHq4QeN1rV9GRpBHEapr86zslO4JfAxkMVXDZthNHhDBybA0bMVK/5t6nCp/Iw5G1QD7ou2qVqCG1D+2whqpBpv4scq/od4rMhdICfqTAUeT1q9ExdMdSXqoK8rrhDYV+gXi31nX/OalcPFI/OVJ36saMChXygoA+Pk7Z50WemTwAT0qKIctrZfKTSXAmgK4tFrT4anwUzrlP72ib3FO1StYXSveqh2F/8g/Y1WgAiUyAmQ3UURqdDTKbaRqermZ7h8cPvTtTvB3ctNFZCUyU0VgW2FYECvkS96gLbxvITs187ikiEqBEqkY45Tw3xjU6HqMD5i0ySNnnRb4bZb2Xv2awWclKj2Fc0BOcD9DerTd31J2ow5asn9rc2d+hkDHQw1uRD0eew750TT1pqZ1F3qhGJgVeC2obFqdpD+ysqsI1UW3uYqqnYQ8EWenb9En4/eNyqjdzTEti6VawtDeCuU3ff7vrAtq7zvuaaQAEfKPCbqk6sOX/SebOrxOdKVgX6iBlquKQrWW0jU1TBHpksa9MLQ5k+AQDkpLh4Y3shfr8fi1Snz8zhhOQJ6tWVzwcNZYEmjWPqzw3lgW3gz8V71J+be7kaq8UWSAYO1STVxt+hNtKxZuL3q74OT3PnETC9+b7QSNXOHhajJj0l5ahtWJxKat1tw2KlE10MCZIAAC3FRb3bQ0FVExlxMnv2rFit6k7XlQzpM09/rM974u67/VWr7rqba1XB7W1RL09g234H39Kl7bvDnzvut4Wol92pZqfaQgNJJOTENiRCTXRqq32EBLZ2p7Svi2FNEgCqBgCgF9dJAhhIVhs4o9RLCDHgJAEA45IDCaCkjhC7lbQYJ9lJMrJFCDG8SUMl4HI6GBETxoe5JXz7L1v4xnNbaGrxGh2WEEL0K0kAATkpLnbkV2O1WCisbuLpjw8YHZIQQvQrSQABWqAf4NvnjOaK6SN49pPDlNY2GxyVEEL0H0kAARdMSGb+mHhu+dIYbjhnNK1eP5uOVBodlhBC9BtJAAHTM2N55TvziAkPISfFRXiIjW1HJQEIIYYvSQDdsNusTMuIYVt+ldGhCCFEv5EEcAqzRsaSW1RHg/sU0/2FEGKIkwRwCjNGxuL1+dl5rJfLFQghxBAhCeAUpmfGYrHA1jxpBhJCDE+SAE4hOszB+JQo3tl1HI+3m2V8hRBiiJMEcBp3LM1mf0k9z2/MMzoUIYQIOkkAp7FsYgqLxibw2//bT2mdTAoTQgwvkgBOw2Kx8NClE2n2eHnsvX1GhyOEEEElCeAMxiRGcuOiMbyxvZCtMjFMCDGMSALogduXZJMa7eR3H8oCcUKI4UMSQA+Eh9hZNjGF7flVMiJICDFsSALooemZMTS2eNlfUm90KEIIERSSAHpoekYsANtlfSAhxDAhCaCHMuLCSIgMYUe+WhqisLqJyT/9F5sPVxgcmRBC9I0kgB6yWCxMy4hlxzFVA/hoXyl1bg8f7Ss1ODIhhOgbSQC9MD0zhsNlDVQ3trDuQBkA22StICHEECUJoBdmjlT9AB/sLWHDoQosFthVWIPbIw+QF0IMPZIAemHOqDgmpEbx0FtfUNfsYcXkVFo8PvYU1hodmhBC9JokgF6wWi3ce3EODS1eLBa4Y+lYAN7fU8SPXtvJ4TIZIiqEGDrsRgcw1Cwel8h5WiINbi/jkl1kxoXzP58eASDMYeOhyyYZHKEQQvSM1AD64NnrZ/HCjXMAOFdLJDbcwcS0KD7MLcXv9xscnRBC9IwkgD5w2KyE2m0APLBiApvuX8p180ZSWN2EXlJncHRCCNEzkgDOUohdJYMlOUkArMmVeQFCiKFBEkCQJEU5mZIezXt7imTBOCHEWfP7/Rwqq+fFTXk8/M5e6ppbg/4d0gkcRF+bk8m9b+zmOy9s4+mvTSc8RE6vEKJn/H4/+ZWNbDxUwcbDFWw8VEFpnRuAzLhwblo8BpfTEdTvlBIqiK6dk4nX7+cn/9zD5c+s57GrpmC3WhiX7MLpsBkdnhBikDlS3sAn+8vYcrSSLUcrKalVBX5CZCjzs+KZPyae+VnxjIoPx2KxBP37JQEE2cq5I8mMC+eOV3Zw5aoNANx4zmgeuGSCwZEJIYzk9/vJq2jksyOV7QX+0YpGAFKjncwdHc/sUbHMz4onKzGyXwr8riQB9INFYxN5//uLWX+wnNe3F/D69gJ+dFEOIXbpchHCLJpavOwrriW/spFNhytZd7CMY5VNAMSEO5g1Mo5vLhjFkpxkMuPDDYlREkA/SY5ycuWMdOIiQvjmn7fwYW4JoXYrY5Nchv1nCyH6T1sb/rqD5Xyyv4xP9pfT1KrWCYsMtTM/K56bFo1h/hh1h2+19v8d/plIAuhni8YmkhLl5N7Xd1Hb7CElyskbty4gLSbM6NCEEGfB7/dzrLKJbfmVbD5cybqD5RRUqTv81GgnV80cwaKxiWTGhZOdFInDNvhaACQB9DOb1cI1szN4cs0Brps3kn/uKOTrf9rM9fNHcsX0dKLDg9urL4ToPyW1zbz5eSFbjlaxI7+K8voWAFxOOwuy4rl58RgWZicwOiFiQNrwz5YkgAFw+5JsrpwxgpHxEVw8OYUf/2MPD769l7d3FfH3m+cPiqqgEOIEv99PSa2bsjo3f92az/t7SogJd3C0vAGPz8/ohAgWj0tkRmYsMzJj0VJc2Ibg77EkgAFgt1kZGR8BwIKsBD6++1xe+Syf+97YzRs7CvnKzHSDIxRCeLw+9hXXsbeolpc35/P5MfX4V4fNwoUTU2j1+Dh3XCLXzR/Z/vs81EkCMMg1szL429Zj/GJ1LhmxYcwdE290SEKYQnVjC9vyqiipdZMa7STEbmX17iLe31NMRYNq0smMC+f+5TlkxIYzLTOG1Ojh2WcnCcAgVquFX1w5mW8+t4Vrnt3ELV/K4t6Lc4wOS4hhxefzU1TbTGqUk/IGN7/+l84/Pz9Oi6fzci1hDhtLxydxwYRkJqRGMSYxckg26fSWJAAD5aRE8fHd5/Kfr+/iT58eZuXcTErr3FgsMCMz1ujwhBhyPF4ftc0eqhtbeObjQ3zwRTF1bg/JUaE0tXhp9vj4ysx0Lp82grQYJwVVTdQ3e1iQHW/KpVvM9y8eZMJCbNy3PIf39xTzw7/tVO2OFnj+23OYJ81CQpyWz+fng70luD1eYsJD+K8395AXmF3rdFi5bOoIclJdbDpcgd8P916cw5jEyPafT48195wcSQCDQGp0GNfOyeD5jXnkpLjw+Pzc9PxW3rxtYaeLVQgBNY2trFp7kEOlDRwpr+dQWUP7exlxYTywYjxWi4UVU1JJjnIC8K2Fo40Kd1CTBDBI3Ll0LE6HjRvPGU2rz8+Kpz7lh3/fyW+unsaqjw/y9XkjmZoRY3SYQgyoVq+P3KJavjiullTIr2hk4+EKqhtbGJfsIj4ilDvPH0dGbBgHSutZPjmVyFAp1npKztQgER8Zyv3Lx7f//aFLJ3Lnq59zwW/W4vH5eWvncX57zTSWT041MEohgqu0tpmKhhbGJEaw6uNDvLeniOZWH8snpzItI4YH/rmH8nq1QqbDZiE9NpxZI2O58/yxTEyL7vRZ06XfrNckAQxSl05N49MD5Rwqq+e/LpnAw+/mcutL27n7wnHMyIwlKszBpBHRZ/4gIQaJ4ppmPswtYX5WPFaLhf9ee4jXtxfQ6vUT5rDR1OrlnOwEHDYLf1x7CICcFBc//fIEpmXEkBYTZoqROQNJEsAgZbFY+PVXp7b//aUb5/Kj13bx6w/2A+pRlB/ffS4jZE0hMYj5fH7W7i/j7Z3HeWdXES2Bp+VZLOrZ2tfMzmBiWjSfHank4kkpXDgxBYDNhyvYVVDD9QtGtj9/WwSfJIAhwumw8eS107hi+gjcHh93vLqDJz7Q+cWVk1m9u4jXtxXytbmZ0kQk+p3P5+cvG46SFhPG4nEJbD5cicUCXp+f5zfmMS45knuW5bAtr4pHV+eyu7CGKKedq2enc/WsDD49UE6Lx8fKeZkkuVQn7X/Myez0HXPHxMvkyAEgCWAIsVgsnBd4+Py3Fozi2U8Ps1Yvo6KhhRCbld2FNcwZHUdCZKjBkYrhyufzc98bu/nr1mOAupP3+0+8HxcRwtr9Zby+vZDKhhZSopz85uqpXDIlrf15GFPSZTDDYCEJYIi69dxsPtZLyYyL4BsLRpIS5WT5U5/y8Dt7+d2109uPO1haR1SYo/1OS4ieamzx8Pq2ArISI0mKcvLS5jw+zC3hWGUTty/JJiclil0F1SzMTsDpsFHX3MqisYl8sLeY59Yd4fvnj+WrMzMIC5EmnMFKEsAQFR3u4IO7vtRp33fPzeapNQc4Z2wiX5mZTk1TK1es2sCU9GheunGeQZGKwczv91PR0ML+kjp2HqvB5/czJT2a49VNrPr3ofZJVQAhNiuLxibwwws0LpuWhiUw1r6rS6akccmUtIH8Z4g+kgQwjNyxJJutRyu5/x+7GZ0QzroDFdQ1e1h/sIIDJXUcLm8gJszB3DHxfHG8Br8fGUk0DByrbCTRFYrT0bM77YOl9by3u4h/7y/jQEkdtc2ebo/LjAvnhRvmUFHfQlmdm8unjyDRJc2Lw4kkgGHEbrPyzNdmcPmq9az802bsVitzR8ex41g1d/99J7sKa4gMsfP8DXO4/rnP8Pn8vHPHIkYnDI+lbc1oT2ENlz+znrSYMO5YOpYkVygzRsYS7rDx8mf5HKtspLHFy5ajlbR4fPiBI+Vq5uy0jBgunZZGVmIkWYmRTA7cDOw5XkNaTBij4yPkWRXDnMXfsQdnENM0bRRwZM2aNaSny/r5p1Ne7+aWF7axLb+KN29byPMb83htWwETUqM4XF6P1+fHYbMSYreSGh3GzYvH4HLaiQl3MDEtuts7yaPlDXx6sJyvz80cEk86Gi5avT7W5Jbg9vhYNjEFp8Mlq3A9AAAMrElEQVRGaW0zv//oIPOz4nnm44OU1LqJDrO3L4mQGu1ES3Hxb72MELsVu9XCzJFq7khzi5dFYxO4aFIqKdHSL2QGBQUFLF26FGC0rutHO74nNYBhKCEylJdvmkdhdROjEyK4Y0kIfj/858Uab+44ziOrc3nw0vGkRju5+YVtfP+vn7f/rNNhZWFWAkvGJ3Hl9PT2DryH383lw9wS8Pu5bv6obr+3we1h9e4iluQkES8jkXqtpqkVqwVcTgc+n5+/bzvGEx/sp7ROzYSNCXcwe1Qc2/OqqGho4YVNeQD8YeUMzp+QzP6SOkpr3Ty6Ope1+8t4YMV4bjhntCRscUpSAzAZv9/PobIGshLVM0trGlspb3BT1+yhrM7N+oPlrNmnRnrMGhnL/35zNh6vj7mPrsFus+Dzw4+WaSS6QrFaLMSGh5CVFEFqdBgPvvUFf9lwlFC7lXsvzhnSC3C9u6uI9NiwM66/5Pf70UvqGJukHglYVucmITKkvdCtqHez7mA5tc0eluYkkRYTRovHh8NmwWKxsOVoJbHhIcSEO7j09+uoaWrlggnJ7DhWTV5FI7NGxnLreVk47Tb+vq2AnQXVxIQ5ePTKyewqqKGqoYXvLB7TqZBv8fgoqW0mI87cK10K5XQ1AEkA4iR+v593dxdx118/Z2ySiyU5STz98UFeuWke976xq9PIEFAPvr9z6VieWnOAZRNTqHN72HCwnLe+dw5WK9gsFsYmuwBVIG4+Usl5WlKvhgd6vD7K6t29fjJTTVMr/9xRyGdHK/mP2ZmcMzbhlMcWVDXS4Pby/p5ifvvhflyhdv5x20KcDiuuUAfR4Y5Ox+88Vs1P3/qCz49Vc66WyNT0GJ766ADfWTyGey/K4YVNefzqXzp1gU7WhMgQblw0hj+uPcS4ZBcXTkjmkdW5hNisjIqP4GhFA0vHJ7FWL2PGyFiunpXBJVNS5Q5enBVJAKJP1u4v4+YXttLc6mPyiGjevv0cfD4/lY0tVDe24PNDZUMLT390kHUHy4kIsfHve87DbrVwwW/Xtt8Rh9it/OyySXx6oJzVu4vw+vwsGpvAn74xq9M0/wa3B6fD1r7eS4vHx86Caoprmnnm44McKK3nzdsWto9cKqhq5HBZAz6/n6c/OojDZuUXV05mVKBTu8Ht4ao/bGBfcR0RITZavD4unTqCtfvL+NqcDO66YBwWiwW/X81sfeTdXDw+9fuwfHIKmw9X0tzqpaHFS2y4g8eumsLMkbHEhDk4WtHIFavWExFi56JJKby4KQ+Pz09mXDj5lY0syIpnw6EKFo1N4O4LNSwW+O6L2ymsbmJiWhSHyxra177x42f9wQqe+OpUrpLnQ4sgG/QJQNO0HwA3AH5gK3CzruvuLseMQhLAgNtytJJbX9rO/ctzuGJ69+e9udXLL1bnMi0zpv2Yd3cVcfsr27l6VgZfHK9ld2ENToeV6+ePIj4ihF+8t4/F4xJ55PJJHCqr5+XN+azZV0qo3cqMzFi+PDWVv2zII7eoFoARMWHUuz1My4jh/317DnkVDVzy1Drq3OruOjkqlMYWL82tXlxOR/uzXncV1PDsdTOZNSqOm57fyva8KnJSXewprGXF5FQy48P5KLcUvaSO88cn8eWpavz6JVPS+PxYNY+/v49FYxN4d3dxeywup53QwKzWf962kPTYcLblVVJQ1cSyiSlcsWoDuUW13LNM49Zzs9rv4EvrmtlypIqLJqVwpLyeD/aW8O2Fo3HYrByrbGxPXEIE06BOAJqmzQX+BMwDGoHngV26rv+qy3GjkARgCL/f36dmiAa3h4hQO/VuD3/dcoyLJqW0L1738uZ8fvbOFzS3qsXB4iNCuHz6CLw+Px/mllBQ1UR8RAg/XjGerMRItBQXz288yqOr9/HAivG8sb2QwuomnvjqVFq9Ps7VkqhpauXP649Q7/agF9exPb+KBy+dyPWBTmuvz0+920OU084v39d5bv0RWjw+pmbEsHJOJl+ZmX7KYY/NrappqKaplb3Ha9lVWMPDl09i5siTlyCuqHdzrKqJafL8BjEIDPYEMA5I1XV9beDv9wCZuq7f3uW4UUgCGFYKq5t4ZXM+OakuLpyQ0r5WjNfnZ3t+FWMSIjqNJmpu9bL0ibUUVjdht1r4n+tnta+N1B23x3vGlST7mtyEGCoG9TBQXdf3A/sBNE1LBe4EvmlkTGJgjIgJ4+5l2kn7bVYLs0fFnbTf6bCx+o5FlNU3Ex0WcsZZqT1ZRlgKf2FmA5YANE27Fni6y+4SXdcnBt4fCbwH/EHX9Q8HKi4xtESHnzwaRwjRNwOWAHRdfxV4tbv3NE2bBrwNPK7r+u8HKiYhhDAzw5uANE1LBN4Hvqfr+mtGxyOEEGbR6wSgadqlwEu6rru67L8J+BGQDnwO/EDX9Y09+Mi7gCjgAU3THgjse1fX9R/3NjYhhBA916sEoGnaAuBFwNJl//XAH4GfAVuA24F/aZo2Vdf1I6f7TF3X7wfu700cQgghzl6PEoCmaaGo0Tk/BxqAkA7vWVAF/7O6rj8U2Pd/gI66u78jyDELIYQIAmsPj7sYuA+4B+jaSZsNjATeatuh63or8C5wURBiFEII0Q962gS0BTWJoFrTtAe7vDcusD3YZf9hIEvTNJuu696ziLGNDaC4uDgIHyWEEObQocw8aWJMjxKAruuFp3k7KrCt67K/DlXDiABqe/I9Z5AKsHLlyiB8lBBCmE4qcKjjjmAMA23rEO66pkTbfl8QvgNULWQRUAQEo0YhhBBmYEMV/lu6vhGMBFAT2LqAkg77I1GFf0MQvoPA6qDrgvFZQghhMoe629nTTuDTORDYjumyfwyg67pu/HrTQgghThKMGsAB4BhwOfABgKZpDmAFaiRQv+rtBDRN0yYBTwJzgUrgGdQSFKZJVH04Z++g/j+7cum6Xt8/UQ4+p5oE2c1xpr/G2vTinJn2GtM0zYYaZn8TkAnkAauAZ051zQTrGjvrGkDgCx8DbtE07RFN05YDbwIJwG/P9vNPp8MEtBeBq4Bq1AS0bh9Gq2laEvAhqr/iauBZ4BHgh/0Z52DS23MWMAV1sc3v8mo8zc8MK6eaBNnNcaa/xtr09JwFmPka+wnwKOpcXQr8Dfgdatj9SYJ5jQVlLSBd11dpmhaGymJ3oe4ql+m6fjgYn9+dPk5Auw31b75U1/VGYHVgktt9mqY9GZi/MGz15ZxpmhYDZADv67q+aQDDHRRONwnyFEx9jUHvz5mZrzFN06zAD4Bf6br+SGD3msAaaXcDj3fzY0G7xnpdA9B1/UFd1yO72f+EruuZuq6H67q+oIfrAJ2NvkxAOx9YEzhpbf4JxAGz+ynOwaQv52xKYLurf0MbtE43CbI7Zr/GoPfnzMzXWDTqKYhvdNmvA4mapnX3nNCgXWPB6AQ2yhknoJ3iZ7o7vuPnDWd9OWdTADfwsKZpFZqmNWqa9ndN01L6M9BBpG0S5FOcPNS5O2a/xqD358y015iu61W6rn9P1/UdXd76MlCg63p3oyiDdo0N5QTQkwlo3f1Md8d3/LzhrC/nbAoQGjjmCuBWVNvsR4Fq57Cm63qhruvVvfgRs19jfTlnpr7GutI07UbUXX53zT8QxGvM8OcBnIW+TECzdHN8m2BNWBvM+nLOfgO8ouv6x4G/f6JpWi6wCdUB9ULQoxzazH6N9YVcYwGapq1EDdJ4jZOfoNgmaNfYUK4BdJyA1tHpJqDVdHO8q8N7w12vz5mu6/s6/GK27duMGj00tT+CHOLMfo31mlxjiqZpd6GS3TvAytMM6QzaNTaUE0BfJqAdOMXxoDpdhrtenzNN067VNG1xl30WVJW9vF+iHNrMfo31mlxjoGnao6ia0AvAV3RdbznN4UG7xoZ6AmibgAZ0moC25hQ/swY4v0vP+uVABWro6nDXl3P2XeDJwHC1NsuBMOCTfopzKDP7NdYXpr7GNE27EzVq6kngm7que87wI0G7xoZsH4Cu635N0x4DntY0rQpYD3yPDhPQNE3LAhI7jC1ehXpa2WpN036Fql7eB9x7how7LPTxnD0KvAe8qGnan1GjDH4OvK7r+oaB/jcMNnKN9Z5cYydompYK/BLYDbwKzNU0reMhW1FDt/vlGhvKNQB0XV+FGmt8HarTJIbOE9B+AmzscHwRqnfdHjj+O8CPdV3/9UDGbaQ+nLN/oWYnZqPGGv8YeC7w80Kusb6Qa+yEZaimrsmoc9L1FUM/XmMWv990y5MIIYRgiNcAhBBC9J0kACGEMClJAEIIYVKSAIQQwqQkAQghhElJAhBCCJOSBCCEECYlCUAIIUxKEoAQQpjU/wfMyoZmTC3CtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(delta_t_list,result_Bayes_array.T[3])\n",
    "plt.semilogy(delta_t_list,dtau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}