{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Triple collocation\n", "\n", "Triple collocation is a method to estimate the unknown error of three timeseries with the same underlying geophysical variable [Stoffelen1998]. It is based on the following error model:\n", "\n", "$$x = \\alpha_x + \\beta_x\\Theta + \\varepsilon_x$$\n", "$$y = \\alpha_y + \\beta_y\\Theta + \\varepsilon_y$$\n", "$$z = \\alpha_z + \\beta_z\\Theta + \\varepsilon_z$$\n", "\n", "in which $\\Theta$ is the true value of the geophysical variable e.g. soil moisture. $\\alpha$ and $\\beta$ are additive and multiplicative biases of the data and $\\varepsilon$ is the zero mean random noise which we want to estimate.\n", "\n", "Estimation of the triple collocation error $\\varepsilon$ is commonly done using one of two approaches:\n", "\n", "1. Scaling/calibrating the datasets to a reference dataset (removing $\\alpha$ and $\\beta$) and calculating the triple collocation error based on these datasets.\n", "2. Estimation of the triple collocation error based on the covariances between the datasets. This also yields (linear) scaling parameters ($\\beta$) which can be used if scaling of the datasets is desired.\n", "\n", "The scaling approaches used in approach 1 are not ideal for e.g. data assimilation. Under the assumption that assimilated observations should have orthogonal errors, triple collocation based scaling parameters are ideal [Yilmaz2013].\n", "\n", "Approach 2 is recommended for scaling if three datasets are available.\n", "\n", "In this example we will show how `pytesmo` can be used to estimate random noise with three datasets using triple collocation.\n", "\n", "We will demonstrate this with an example using synthetic data that follows the error model from above." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABVbklEQVR4nO2dd5icVfXHP3f69t5L6qaRhISEUEIJoXcB+QmoSNEISEchVClCIiBdiBRFEEFFQJQmogmBEEglfdOT7b3N7E597++Pd7J1tk3Zej/Ps09m571z75kk853znnvuOUJKiUKhUChGPobBNkChUCgUA4MSfIVCoRglKMFXKBSKUYISfIVCoRglKMFXKBSKUYISfIVCoRglhCz4Qog8IcT/hBDbhRBbhRA3BRizQAjRIITY6P+5L9R1FQqFQtE/TGGYwwvcJqVcL4SIA9YJIT6VUm7rNG6llPKcMKynUCgUiiAIWfCllGVAmf9xkxBiO5ADdBb8fpOamirHjh0b6jQKhUIxali3bl21lDIt0LVwePitCCHGArOBrwNcPkYI8S1QCvxcSrm1t/nGjh3L2rVrw2miQqFQjGiEEAe6uxY2wRdCxAJ/B26WUjZ2urweGCOltAshzgLeAwq6mWcRsAggPz8/XOYpFArFqCcsWTpCCDO62L8hpXyn83UpZaOU0u5//CFgFkKkBppLSvmilHKulHJuWlrAuxKFQqFQBEE4snQE8AqwXUr5RDdjMv3jEELM869bE+raCoVCoeg74QjpzAd+CGwWQmz0P3cXkA8gpVwGfBe4VgjhBVqAS6Qq06lQKIYQHo+H4uJinE7nYJvSJ2w2G7m5uZjN5j6/JhxZOl8AopcxzwHPhbqWQqFQRIri4mLi4uIYO3Ys/oDEkEVKSU1NDcXFxYwbN67Pr1MnbRUKhQJwOp2kpKQMebEHEEKQkpLS77sRJfgKhULhZziI/SGCsVUJ/gjHq3nRpBbwmiY1mj3NSCmxu+0AlDvK8WiegTRRoVAMEGE9eKUYenxR8gUx5hiOzDwSAKfXicVowSAM7KjdQWVzJWPix3CgseNZjRNyT8AglD+gUIwk1Cd6hOHTfK2PD3ntDo+DZk8zhbWFrC5bzbaabXg0D/XOegAaXA1d5vm8+HMa3Y3sqd9DYW3hgNiuUCgiixL8EcT+hv2sLFmJ0+ukqKmItRVtZSnWVqylzFEGQHVLNV+WfIlbcwNQ76oPON/6ivUUNRVR5iijwlGBR/PQ7Glu/aJQKBThY82aNcycOROn04nD4eCwww5jy5YtYV1DhXRGEPsb9wOwumx1l2vdxfH7yvba7V2em5I8BYvRQrwlHpNB/VdSjBx2VjTR5AzvXlaczcykjLhurx955JGcd9553HPPPbS0tPCDH/yA6dOnh9UG9SkdISwvWj7ga+6o3dH6eEHeAgDK7GVEm6NJsCYMuD0KxXDnvvvu48gjj8Rms/HMM8+EfX4l+CMAr+YdbBNYXrQck8HUasuhLwCFYjjSkyceSWpra7Hb7Xg8HpxOJzExMWGdX8Xwhzlun5vK5srBNgMYGl88CsVwZtGiRTz00EN8//vf54477gj7/MrDH4a4fW5Wla4iOzabUnvpYJujUCjCwGuvvYbJZOKyyy7D5/Nx7LHH8t///peFCxeGbQ0l+MOErTVbafG0kGRLoqipCGBIi71H82A29L2ok0Ix2rn88su5/PLLATAajXz9daA+UqGhQjrDhKrmKuwee6vYD3W+LPkSt89Ng6uB4qbiwTZHoVCgPHxFBFlVuqr1cW5c7iBaolAoQHn4w4JDJ2aHM+sr1mN321FtEBSKwUN5+EMYKSUrilcMthlhodHdyNqKtaTYUpiRNmOwzVEoRiVK8IcQq0pWkReXR0pUCi6fC6d3eHTe6Q81zho8Pg9mo9rQVSgGGiX4Qwi35mZPwx72NOwZbFMiypelXzI/Z77K4lEoBphwNDHPE0L8TwixXQixVQhxU4AxQgjxjBBitxBikxDiiFDXVQxvviz5ckTewSgUQ5lwbNp6gduklFOBo4GfCSGmdRpzJlDg/1kEvBCGdRXDnNVlq2nxtgy2GQrFqCFkwZdSlkkp1/sfNwHbgZxOw84HXpM6q4FEIURWqGuPJEart3uoZLNCMdq59957efrpp1t/v/vuu8NeQC2sMXwhxFhgNtD5iFgO0P7EULH/uS6fdiHEIvS7APLz88Np3pCl3lnPxqqNg23GoHCw8SA+zUdWTBaxltjBNkeh0KncDs7G8M5pi4f0qd1evvrqq7nwwgu56aab0DSNt956i2+++SasJoRN8IUQscDfgZullJ3/pgJ12w2YkC2lfBF4EWDu3LkjNmnb6XWyr2EfVS1VIdeqH+6U2EsosZeoCpuKUc3YsWNJSUlhw4YNVFRUMHv2bFJSUsK6RlgEXwhhRhf7N6SU7wQYUgzktfs9Fxi6hWAGgJ11O6l11g62GUOK5UXLyYrJYnLy5ME2RTHa6cETjyQ//vGPefXVVykvL+eqq64K+/zhyNIRwCvAdinlE90Mex+43J+tczTQIKVUwVtFFw61U1QoRiMXXHABH3/8MWvWrOH0008P+/zh8PDnAz8ENgshNvqfuwvIB5BSLgM+BM4CdgPNwJVhWHfYUmovVd59D2yv3U5GTMZgm6FQDDgWi4WTTjqJxMREjEZj2OcPWfCllF8QOEbffowEfhbqWsOdFUUrmJA4gd31uwfblCFPqb2UnXU7OTrraGwm22Cbo1AMCJqmsXr1av72t79FZH5VPG2AkFIikUrs+8jOup2AXophQ+UGla+vGPFs27aNiRMncvLJJ1NQUBCRNVRphQFC3+pQ9Jeq5ioaXA3sb9jP1JTB2UhTKAaCadOmsXfv3oiuoQR/gNhVt2vQ1j4u5zjcPjfflHfM6Z2UNIlkWzKa1LAYLTg8DnbW7cThcQBgNpgZnzCeqpYqcmJz2Fy9ecBtr3fVD/iaCsVIRQl+hNlbv5cGdwMNroZBWX9ayjRMBhMmg/5PLRAckXEETq+TtOi0DmMTrAkcmXkkUkqqWqpIj04HICtWPxQ9L3Nely8NhUIxfFCCH2EONh0clHVnp89GkxpJtqTW5+ZkzMFitGA1WomzxHX7WiFEq9i3J9oczeSkyRTWFQ54A/WK5goMwqBy9BWKEFCCPwIRCOIscRhExz35nkS+r2TFZpESlYLFaCE1KpXC2kJcPlfI8/aFMkeZEnyFIgSU4EeIZk8zHs0zIGulRqUyOWkyJoMJt+bGarRGdD2L0QJAsi2ZY7KPweFxsKV6C4nWxIgXQ2vxthBlioroGgrFSEUJfoQYqFj3/Oz5HbpHRVrsAxFjjuGorKMAmJw8meVFyyO21tdlXzM5aXLrvoJCoeg7Kg8/AgxEzviM1BksyFswKlsFFtYVDrYJCkVEWLZsGbNmzWLWrFmMGzeOk046KazzKw8/Anxd1rk6dHg5NvvY1rDKUOTwtMP5turbwTZDoQia3XW7sXvsYZ0z1hzLxKSJPY655ppruOaaa/B4PCxcuJBbb701rDYoDz/M6FUkIsfkpMlDWuwBkmxJzMmYE9E1yuxleDXvgG0YKxQDyU033cTChQs599xzwzqv8vDDRLOnGU1qlDvKIzK/yWDCq3mHTew6zhLHgrwFNLga2FC5IezzF9YVtoZ2VB19RbjpzROPJK+++ioHDhzgueeeC/vcSvDDRKQ3aY/MPLL1BOxwIsGawIK8BRHdyFUoRgrr1q3j8ccfZ+XKlRgM4Q/AqJDOMMBqtGI1Wkm2JQ+2KUFjMVhItCYOthkKxZDmueeeo7a2lpNOOolZs2bx4x//OKzzKw8/DNQ56yIyb7wlnqkpU1vLIgxnjs05FtD3OFYUrwjr3BsqNzA7fXZY51QoBoM//OEPEZ1/+CvJILOvYR8HGg+EfV6B4LDUwwYlrz6SRKJqaIOrAbfPPeQ3sxWKwUYJfohEQuwBTsw7MSLzDgXSo9Np9jSHNe1tVekqQK8hlGBNCNu8CsVIIiwxfCHE74UQlUKILd1cXyCEaBBCbPT/3BeOdUciObE5TE+dPthmRJRpKdOYmzm3S7XOcBCJjCCFYqQQLg//VeA54LUexqyUUp4TpvUGHSklnxd/HtY5DcJAQVJkOt0MRWxG1bpQoRhIwuLhSyk/B0ZVV2635kYSvkNWM1JntNajGS3kx+dHZN5InYVQKIY7A5mWeYwQ4lshxEdCiMMGcN2w4tW82N12vir9KmxzzkybSUpUyojboO0Ns8Hc+p7DeXhqR+2OsM2lUIwkBmrTdj0wRkppF0KcBbwHBIxdCCEWAYsA8vMj4wGGwuqy1Xg1b9jmG5cwbljn14fK3My5eHx6GelZabPYWLUxLPMWNRWRE5vTpSeAQjGaGZBPg5SyUUpp9z/+EDALIVK7GfuilHKulHJuWlr4N/VCJZxiD5AbmxvW+YYbZoOZaHM0AIm2RLJiwlM6Yk/9HjZWbgzLXArFSGFABF8IkSn8CdhCiHn+dWsGYu2hzKy0WRgNxsE2Y0gRzo5Wje5Gap2jamtJMQIoKirivPPOo6CggAkTJnDTTTfhdrvDMne40jLfBL4CJgshioUQVwshrhFCXOMf8l1gixDiW+AZ4BIZ6bKSEWBLdcCs06BItCaSaEsM23wjieNyjgvbXJuqNoVtLoUi0kgpueCCC7jgggvYtWsXO3fuxG63c/fdd4dl/rDE8KWUl/Zy/Tn0tM1hyfKi5aRHp1PdUh22OcPRX3akYjKYwlpwrdxRTmZMZljmUigiyWeffUZ0dDRXXnklAEajkSeffJJx48bxwAMPEB0dHdL86qRtH6lsrgzLPCaDiRmpM5Tg94Gjso4KSzOZHbU7lOAr+sXNN9/Mxo0bwzrnrFmzeOqpp3ocs23bNubMmcNZZ51FaWkpAOeddx75+fns3r2bmTNnhmSDEvxecHqdYZknJzZnVB2qCgdRpqiwiX6ts3ZUZ0MphgdSSoQQfPjhhx2ef//998NSh0oJfi+Eq8590GLvaQFzVFhsAEDT4MCXkD4VYvyJUl432CsgMS9864SJKFN43vumqk2qUYqiz/TmiUeK6dOn8/e//73Dc42NjRQVFTFhwoSQ51dJyr2gSS3kOYIub9xSB3uXQ31RyDYAIKU+p9sOxWugejcUfgR7PoOKLfpjAGcjeJzgbAB3MzSWwp7/6q8fBMIl1NtrtodlHoUiUixcuJCWlhZee02vUuPz+bjlllu46qqrQo7fgxL8iFOQVBB81onb3+GqpY/19jUfVO/SvXgArwtc7SpSlm+G4nZ3LDW7us5R+JF+B7BvBRxYpf9ZsVWfq2a3/kUwCKTYUkKeo6K5Imx7MQpFJBBC8O677/L2229TUFBAQUEBMTExPPzww2GZX4V0uqHF2xKW2HFObE7wL+7No3bUgDUWTP6SDLv+rf9Zswc61/lJmQiNJf1Yu92dzaHDZjW79Z9JZ8DOj/VQU8YMiAldjHtjRtoMfJqPlSUrQ5pnW8020qPTw2SVQhF+cnNzef/99yMyt/LwuyEcYm8UET5UVfyNHvLxujpdCPBFUbM7fOvu/Fj/09Oi2+BsgModULs3gC3hI1yH1A6VclAoRhvKww+A2xf6qbawp156nLqgxmaAEOAvR4DU9Pj6YHJgVdvjqkLInAlx/hIJEWjEHCpfln7J1OSpZMRkDLYpCsWAogQ/AOE4nZkSFeYwR8VWcFRCvb/D1lAuCla+Sf8ByDsKrPFgDM9/NZPBFJZ6RttrtyORKj9f0YFDaZHDgWCKFQxh1Rh4DjYeZHnR8pBb74VFRKSEhmL9cXO1LvYdroeePTQgFH0Nuz9t+72lHorXBZ3xc1TWUcxMC+3wySF21O7Ao6nwjkLHZrNRU1MTlJAONFJKampqsNn610RIefjt2NewLyzzFCSG4YBV3X5w1uuPIxgXHzDKN0PyBChZBz63noFkje33NGaDmWRbMgvyFrC3fi8Hmw6GZFZtS60K7SgAfbO0uLiYqqqqwTalT9hsNnJz+1dtVwl+O8LVwSrkzcWaPVC9Myy2DBkaitvuWABcjaB5ICop6CnHJ44PWfC3125Xgq8AwGw2M27cuME2I6KokI6fOmcfc90jjaaNPLEPRNm3cHB1yNMcnnZ4yHMUNYXpYJtCMcRRgu/H5RsCYZP6g7Drk8G2YmBpKIam8o4HxPpBki0p5NaQe+r3hPR6hWK4oAQf2N+wP+Q+qDHmGADSovvZpUvT9FIGhR/pmTijjfLNULoB9q+E5uCalczLnBeyGXvr94Y8h0Ix1Bn1MfyixiL2N+4PaY7JSZNJjUrFbDT374WaNvo8+p4oan/YTcDEk6EPf6fhOJB1sOkg4xPHhzyPQjGUGfUe/p6G0G7nT8g9gazYrP6LfX3RoIq9pkm8mkZVk4v91Y4u173aYKd9Sn0vo+5An0bPz5kf8oobKzeqloiKEc2o9vBDzcGOs8Rh6OsBKGeDnl5pS4CyTXpu/QAhpWRnhZ0Wjw+ryYDL21XMNxbVB3yt1WRgQpqePmkxDbB/UO/PwEka0+tQs6GfX7iBlnPVU19Vr8ooK0YsYRF8IcTvgXOASinl9ADXBfA0cBbQDFwhpVwfjrVDof5QnnsQRJuimZMxp+dBFdt0gXc7oHbgNgadHh81DjdeTaPO0fFLLZDY94TLq7GtrLH195QYC26fRl5SNBaTgfpmN1EWI1bTADRj93mhcptey7+/d1T9wOVzhbwRrFAMRcLl4b+K3rP2tW6unwkU+H+OAl7w/zmobK0JfpM0NTq15wENxW1lECLMjvImnB7fgKxV49DrDLX/EgCYlZcYuUULP4LEMWAv1++SnPUw7oQuw9Kj0zEZTMSZ4yisKwx6ua9Kv1JevmJEEq4m5p8LIcb2MOR84DWpn1leLYRIFEJkSSnLwrF+fymxl4R8qjbBktDzAHdzSPP3Rq3DjSYlxXUtQb3e5XbT0GinxelCSommSYxGA3GxMcTFRmM29e+/xsaiepJjLOQk6ke9jeEumtb+y9Pddc8BYFrKtNbHEsnOulFwnkGh6AcDFcPPAdqfbin2Pzcogr+rLkDjj37SY3E0twOaIvPWHG4vpfUtOFy9e/Ren489+4vZsmM3+4pKOVBUxoHiMqpr62lx9nzuIDY6iuzMdHKz08nPzmRKwVgOmzyenMz0botL1Trc1PrvAKZmxVHR6CI70YZpECpmZsdmhyT4Da4GEqy9fKkrFMOMgRL8QAoRsI6BEGIRsAggPz8/kjYFTa+x+32fR2Rdj09jV0X3B5SklOzeX8TKrzfy9brNbN6xu1XYbVYL+TlZTJ4whhOOOYKk+DgSE+KwWa0YDAYMBoHH48XuaKbR7qC2vpGSskr2HShh5dcb8Hj0CpUJcbHMnTWNY+fO5Og5M8jOCHzuYHtZE6B/CWTEW8mIs2EwDGwVQovBglsLrtT1wcaDzEibEWaLFIrBZaAEvxho3yE7FygNNFBK+SLwIsDcuXPDVrau2dOMEILqELNjDk87vPs691JGJG7fXQaNvqRk+659fPCfL/jvF2sor6oBYNL4MZx3+onMnDqRGVMLyMlMwxCkp+3xeNm9v4htO/eyadsuvt6whc9WftO6zuknHcMZC44hOzOw+Fc0umjxaIxLiQ5f6dkDX+nF16JTIT4r4JAEawJVLcEVwqpx1uDxefqfbqtQDGFEuEqB+mP4/+omS+ds4Hr0LJ2jgGeklL0ej5w7d65cu3ZtWOxbXrQ85DlOyD2h+zTMvSvAE764vcvrY1eFHa8W+N+nvrGJ9z5azj8//Zy9B0owm03MP/Jwjj/qCI6bN4v01OCLkvWGlJJ9B0v5cs1GPv38azZv17tpzZ4+mYvPPZVTjp+H2RzYl5iRk4Ax3J5+wekBG61srd4atOAfQm3eKoYbQoh1Usq5Aa+FQ/CFEG8CC4BUoAL4JWAGkFIu86dlPgecgZ6WeaWUslclH0qCnxqVyvTUTt9lriY9v94SE5ZCYIeQUvJtcUPAa3v2F/Pmex/zwX++wOlyc/i0SZxz6vGcduJRxMf1v9xwOCgpq+STFat576P/UVRaQUpSAhecdRKXnn86yUld4+B5ydGkxFjCa8S4E8ES3eGp6pZqtlRvCWna7NhsJiVNCmkOhWIgibjgR4qhJPgBPb3Cj0KaszM+TbK5JLDQF+7ez7LX/87yVeuwWsycdfJxXHrB6RSM698+x9SsOEwGgUCgITEZDOyqbMLh8jE+LYY4qwm3T7KzopF4mxmH24e7j7n7mqbx1brN/PX9T1n59QasFjMXnX0yl198TsA7jmirkTHJ0eHL4T/UWrGTt9/gamBD5Yagp82Py1dlFxTDhlEv+JrU+Lw4tI3USAu+pkk2BRD73fuKWPb63/ls5TfExkTzg4vO5P/OO5WkhPhu57KaDGQl2jAbDHg0DbPRQIwl9O2aAzUOHC4f8VFmqu09Z/nsLyrl92++z4effYHBaOC7Z5/MT35wQUC7k6LNpMVZiQ6Dja1MPrPDryuKVoTU70CFdhTDhVEv+IW1hZQ5gk+TnJsxl1iLP1wipV7Ot6E4bOURAm3K1tY38vyrf+Pdj/5LlM3G9y88kx9cdCZxsTHdziOA1DgrOYlRoRkUl6UfbvL0nOPf3ZdUe0rKKnnlzX/wj0+WEx0VxY8vO59LvnM6VkvXkM5h2fGYjWFK4TSYwRyl99Q1GHFrHlaVrur9dd2QG5vLxKSJ4bFNoYggo1bwpZSsKF4Rkg3zMucRbfbHhjUf7Pp3SPO150CNg7rmjqUPPF4vf3v/U5a9/neam5187/zT+MkPLiAxvpvMIGB6Tnz3ue6JY9oyh7IO10W8fYOVqCTIP1r/AotKApNNb5C+b0Wb4BtMEJsBjSVdpj/0/0dKaGjxcKA28Mb1ngPFPPXin/nim41kZ6Zx+3WXc+IxXdNbJ2XEhtfTB0geD0ljWbH3Y2RU8Ln183Pmh6Vmj0IRSUat4Hs0D1+WfBmSDR1u5Vvq4eBXIc13iEBe/dbCPdz/+Ivs3l/EMXNmcNu1P2TCmMA9KydlxLK/ppnsBBuJ0d1sgAoDjD9JF2wkHCoj7LKD5gVrPAih/3TGXgWVW2HsCR1j4u5mvYF68RrwOru8zOPT2Fra2OX5Q6xev5nHX3idPfuLWXjckdzxsx+RnprcYcyY5Gjio0zhO61rsoLRSkVTEdvjU/1/H8GhQjuKoc6oFfz1FetpdHcvPr3RwaMLk3cfaGO2xenihT++zRvvfEhqchKLb7iCBcfMCZizHjCtMXUSpExo21OITQd7JRSc1iby4cbn1YXfYNTd+92fdhlSUt9CVVPXWL/H4+W1tz/gpT+9g8lk4oarvsfF554S8JxA2MI8BjNoHvanTWS/o+udSl9Rgq8Y6vQk+CO2PHKdsy4ksQd/yV2PUw931IRejsHl1djeqejY+s07uP/x31FUWsFFZ5/MTT+5lLiYjumFSTFm0mKtWE3GNrFPyIXYTN2Lj+lU5iGnl5PA4cDY6b+OJRYS88Aco3vU1TvJoYr0OCu1DjdlDW13A2aziasvPZ/TTjyaR57+PUufe5X/rVrLAz//KRlpHd+Lw+0lMSp8KZxjfRr7Q3h9cVMxuXGB77oUiqHOiPXwixqLgm9u4nWDwcSCjDlwIPiNvvZ09na9Ph8vvv4Or7z5HtkZafzytkXMPXxal9d18OhjMyDniO4XcTaC0QJmW1hsDplOWUw1DjdFnWL8Ukre+fC/PL7sT5hNRu668SrOOOnYDmPGpcaQEBW+2HntmKPYVLUp6NfPTp+t6uwohiw9efgjsuNVk7sptE5WFZs52mcIi9i7vRobi+o7iH1JWSVX3/ogL73xLueccjxvLVvSRezHp8UwNSu+Y/gmvesXQgds8UNH7AOQkpjIzJyOQimE4KKzT+Yvy5YwNi+bOx95jruWPIejuS1DaF+1o8sXRSgkm+NJtCb2moXUHaHk9CsUg8mI9PAPNB4IqfzxMfZGrGHIxvBpGptLOoZw/r1iNQ898RIAd998dRdvFtpl3Yw7EYpWQ3wOxKRBdHKXsUOaxlL9jsPr1P+MSYOdHwOBN629Ph+v/PkfvPinv5Ofk8lj997MxHF5HcbMzEkIvQibycpWrZmqmkLIOEzPTOonUaYo5mTMwRTCBrBCEQlGnYfvk8E3A0m0JoZF7Juc3g5i7/X5eOJ3b3DHr55h/Nhc/vK7pQHFfkZOgi725mi9VMCEhZA2efiJPUB8NsSk+vcb0jtkAxVkxDI9J57xaW3nCkxGIz/94YUs+/VdNNmb+eEN9/LPTzsemNtU0kB9c3AVMFvxurDZ/WcoNG9QU7R4W0LK61couqOoqYgWb3B3n70xIgX/YOPB/r9I05hQW8yM2uAzONqzp6qtjHFtXQPX3rGE19/+gO+dfxovP35vh8qS0VYjh+cmMHPGLIzj5uvdnMafGBY7hioxFhMmg4F4m5nDczuGeY6cdRhvvvAIh02ewH2PLuOhJ1/C7W47r7C/pplmtxetm8JyfWFcVBrTY/KgamfA9NK+oMnBbvSuGGl4NA976vfwddnXEZl/RAp+v/F5oHIb2dYkjH1tSt4Nm0vqO4QrtuzYzWU/u5vN23fx4O3XsPj6KzpUksxLjmZSehyi4FQMGVP1HriW7k/TDnvyjtLTSA89Ro/j50zreLeTlpLEskfv4spLzuOdD//HNXc8Qm1dWzrrzgp7r6d8e8IgDKRa4gAJFcG3ulQowkWtsxZnkM5HX1GCD1C+iSgpQxJ7p8fHxqJ6fO2cvo//t4qrb30Io9HIq08/wLmntvVhFUBmgq2tauRoqbsenayfGRh7vP44YzqkTyMtewyz8hKZntNWa8dkNHLj1Zew9O4b2LZzL9+//l527u1499bkDC4kc4gj4sbpD2r3BvX65jCWxFaMXuxuO5uqNrGuYl1E11GC7+89Oyd+XNBTSCnZUd7U4feX3niXOx95jsOmTOBPzz3ElIljW69Pnj6Hw8ekknn46XrmzaQzgl572GL11yZKzIOkMfrj6FRMNj2uH2ttuws6fcExvPLEL/H5fFxx0y/535drWq/tqbKzvzpwj9u+0FpQraVOT8ftJ9+Uf8OKotDKdygU3iD3kvrL6Bb82n1QtZ2jEwowieBOpNpd3g61691uD/c9toznX/0bZ598HMuW3tmhQuTUrHiikrKh4FR9UzZpTODSBqORvCNh7HHE28xMzIgnztYm+odNHs+fnvsVE8bmcuv9T/L62x+0Xqtv8fTYFawnog1WABJN0eAIrllKKFU4FQqAiuaKAVlndAt+Sy0ApiBDOc1uL7sr2zZnGxrtXHfnUv716Uquufy7PHTHtVgsbaGaiWmxWE2GkGq5jHiEgPELYOKpTEiLJSPe2nopPTWJlx6/l1OOn8cTv3uDx5e9jqa1xdCCEX2zwciCpGmMsaWBvVwvGREEQzm9WTF00aTGjtodIVXz7Q+jT/C9bv1EqqMGIwamRGcH5d3bXV52tmsoXlFVw5W3PMCm7bt45M6f8dMfXthaC8cg9HTLWJsJ8o8Z0oejhgTmKL1gW/o0shKiOhzWslkt/PqeG7n0O6fzxt8/4s5HnuuQwbOxqD6o7J3WL/3yb4MyeUXxCsrsA/OhVYwcSuwllDvKB2y9sAi+EOIMIUShEGK3EGJxgOsLhBANQoiN/p/7wrFuv/C69YJiFZv1ujj1+zk+aQqZ1sSgpmvv2e8vKuWKmx+gsrqW55cs5syF81uvTcqIZWZuIkabv7xxVHDrjUqSxoA5CoNBMC2rLSxmMBj4xXWXc8uiy/j3itVcd+dSmuxtcfxNJQ399rjjTFF6miZAyTr9p3xLv1I2C+sKcXiC309QjD48mqf3QWEkZMEXQhiB3wJnAtOAS4UQgWoArJRSzvL/PBjquv2mYjM0FIU8TbPb2yF0sLVwL1fe8gBut5uXf3NvhxIJBYdqu1tiIG8e5Pbat13Rmdx5kJCLxWJhcmZbTwAhBJdffA6P3Pkzvt22k6tueZDq2vrW698W91/0Uy1xGGi3n+Jz9TtlU+XmK/pDUGeGQiAcHv48YLeUcq+U0g28BZwfhnnDh6NrZ6qpMTn9nkbTZIcwzur1m1n0i18RHWXj90/+skMmztj0BGLiUyB9qn6QymTtWtVS0TuWaMicAQWnEpU6hsNzE0hu1wD9zIXzefbh2ykpr+SqWx6gtKJt47W7RvA9EWcK0C2spR60vgm5T/OpeL6iT+yu2z3ga4ZD8HOA9q5zsf+5zhwjhPhWCPGREOKwMKzbdzwdb8vNwki6ufuesIHwaVqHgz6frfyGG+95jOyMNP7w5P2Myc1qvTYzJ4HEySfAmGMhaWxIpivakTEdIQT5ydFMTIttffroI2bwwq/vpK6hiatveZADxW2x9P5u5E6PyWNW7JiOT9bugbINem2gXthYtZEt1Vv6taZi9OHRPBTbiwd83XAIfqCcws4uznpgjJTycOBZ4L1uJxNikRBirRBibVVVcGlyrZRthJrd4OiY8jQ/cXLA5iLd0bkI2ifLv+KOXz3DlIljeeWJe0lPTWq9Nmt8NoYxR7flmSvChxB6QbmJpxJrMzEpo+3v+PBpk3jp8Xtwut1cfeuD7NrXdqu8saieZnffsm/MBiOJ5m5OOjf1bVO2xlnTp3GK0UukauX0RjgEvxhoX9IwF+jgCkkpG6WUdv/jDwGzECI10GRSyhellHOllHPT0tICDekbmqZ3qXJ2vK3PtwVctltc3o5i/8F/vuCuJc8x87BJvLD0TuLj2kRnxtzjYexxw7PQ2XDBEq03X5l8JtGHnd2hDs+UiWN55Tf3YTAY+Mltv2JrYdvp2Z0VduyuvqdczojND/pshkLRE1JK1lesH5S1wyH4a4ACIcQ4IYQFuAR4v/0AIUSm8LvUQoh5/nUj5wZpGpQFTq/Ls/Y9ju7y+jp0qPrHx8u599EXmDNzKr99+HZiovV4b3pqCrNOuRRjYm7kWgoqumIwIMadwKwFF9KQPAuA8WNy+P2T9xEbE8VPb3+YTdvaOpXtrrTj9vYtFp9ijiXJFN31Qsk6PaYP+intpsAHZpYXLWd7zXbcvhAreypGHKFU8w2VkAVfSukFrgc+AbYDf5VSbhVCXCOEuMY/7LvAFiHEt8AzwCUykjtbPhfQ9YN9bMIkzH0U5BaPj+1lbeUS3vnwv9z/mxc5avZ0nn7oF0RF2XDbUokyG8lOVZuxg4Y1FkxWTjxiautTuVkZvPLEfSQnxvOzO5eyeXvb5ti2Ti0me2JCdGbgC7V7dKeiajs0dh+HrWiuUCWUFR2oc9axvWb7oK0fljx8KeWHUspJUsoJUsqH/c8tk1Iu8z9+Tkp5mJTycCnl0VLKyH4KAnhdR8VPxNLHE652p5fCdrVx/vr+pzz05MscN28WTz10G1E2/fRnwpgZTC6YBBm9dKJSDAhRKW29ZjPSUnjp8XtISoznusVLOoh+577C3WEzmDksJpdEU4CYfpnqeqXoP99Wfdv7Hs+hcHQEGJknbVu6/oVGGXtvhO3xaZQ3tLC7XS37v/7zU5Y8+wdOPGYOv/nlLVgt+jzJ2eOYnJcJ2bP0k6GKQWdydgqz8hKJMut3cRlpKbz42D0kJsRx3eIlbNmhi77L33ayL6RZ4pkVN6b3gQpFL9jd9u4vuh36D+jORNnGiNgwMgU/SLaWNlLe2NZ79r2Pl7PkmT9wwtFH8Ni9N7XWxUmIMpM/7ejBMlPRHSb9zmvy7ONITdHDbJnpKbz0+L0kJsRx7eKlbC1s63Xcn+ydHjkU0+8Gh8dBcdPAp+AphgZSSlYWr2RtRQ/tWqt26D8lqjzygODydryF+vCzL3jwiZc4du5MHrv3ptamJePnX8S4Yy4YDBMVvZE8HrIOh/hscmef1vp0ZnoKLz5+DwlxMVxzx5IOot/+IF1PHJ1Q0P3F2j3gbNIb6QQ45LemfA2763dTE+DOUzHy0aQ2qBu17VGC76f9Bu2nn3/NvY++wNzDp/L4L29p9ezjxs4hPsqiMnGGKkLofXT9zDr+XOpSjwQgKz2Vl35zLwlxMVx359IOjVQ2FtV3+cLvjK23Psc1O6F8E9Qf6Lauviq7MDrptXx2feglX/rKiBf86TF5zIuf2OOY9vHcFV+t465HnmPm1AKeevDnRNmsuKIySDn8LCZMnBRhaxVhxRrLgsMnkDL7PJzRWWSlp/LiY/cQZbNy7R2PsL+o7bjI9rImWjw9i36yKbTDdP057KcY/jR7mtnXsK/70sctDXoIx1E5YDaNaMFfkDSNVEsc0d1s2FbbXR3EftWab/nFQ08zeeIYnn34dqKjbFRnLWDWUSeSl5YQcA7F0EYIQV5KDPZEPW0zOzONZb++C4Br7niE0vK209ztM7MCYeyrYA9Q9yLF0GZL9RYONB5gT/2ewANqh2ctnWGJx6dRXNd2vHnNxq3cev8TjM/P4fkli4mNiaY+dS4GowGrSYVwhjvHTGg7KzE2L5vnl95Jc4uTa+54hMrqutZrJfXdH3lPMcd1e60DVdshwDETVWNndKEFOAsE6CU6Irw52x0jUvCPTZjEnLiee9RuLW3Lxd68fTc33fs4OVnpPO8vl+CKymBMbjYLp2RE2lzFABBjNTFl9nxaonNwRWUwecIYfvvIYmrqGrh28SPUNej/H6qaXLi6OY2baU1kfsJkTkycGvB6B0rXg6vrHUN1S7U6fTuaqS/qUxG+SDEiBd9iMAUuc+unfRhnz/5ibrj7UVKTE1n267tITtSraLqtyYxL7aaIlmJYYkvJ55hjT6Ap6TBqMo5jxtSJPP3Qzykpq+Rndy6lyaE3tN9e1oivm3LIZoOx77H46p1dntpSvYVVpat6zslWDHsONB7A2bl5js8zoPH6QIxIwe8Jj6/tg1xaXsW1i5dgNpt4YemdpKXoVS8To8wcP3fWIFmoiDSnTMtA+vd15h4+jcfuu5ld+4q48e5HaWnRP6SbSxo5WNvc7Rx9rslUvRM8LW0bdG59TtUZa2SiSY1VJavY17BPf8Lr0junlazTs7gGmVEl+NvKGltDOTV1DVxzxyM4XS6eX7qYnKz01nFjUgIUzVKMKE6emk596lzsCZM4/qjZPHLn9Wzavotb73+ytUdurcPdwUFoz4ToDI7oJWwI6GGdym1tG3R+z77XVD3FsMTpdeLW2oXs7JX+2l5958j48RwZPz7MlumMCsHXNMnGovrWSolNjmauv+vXVNXU8cyvbidz9hk0Js2gJXkaM064ADH2uEG2WBFphBDMmzYOZ0wuNZkncNT5V3HfLT9h9frN3PvoC/j8Qr+1tJG65sAx9/gewoa9saN2B2vK1wT9esXQxO7pHKrr/xd7jNFGjNEWHoM60bdqYsOc9hUSnS43N9/7OLv3FfHUQz8n9+Sf4PIfpDph2qEN2sj8ZSuGFtEWE7PzE9lwsB53VAbnn7GAhiY7T774Z+LjYrjrxqsQQnCgppl4mwmjIQz+UUMRuBohJh0V1BlZlDvK2VG7o+0Jnwcc/WvidGxCZM/6jHgPf2NRPV5N/5b1+nzc8atn2LClkIfuuJb5Rx7eemp24ZT0nqZRjFBSYq0cV6A3xXFFZXD5xedwxf+dy9v/+owX/vh267jNJY00OcOUX+9sgJpd4HWqcgsjiA5iD/2O2R+fOKXPFX2DZcQL/iE0TeOB37zI56vXs/iGKzjjpGOpzZgP6Jt4BoM6BTlasZmNnDItA7vfu7rxx5dw/hkLeOmNd/nzux+3jttT1bVr1tSYQO2b+0jFVjZXBm7UoxjG+LxB5dkbReTleMSGdBqdHvZW6TfNUkqe+N0b/OvTlVz7o++ycNHD1GoeNKNVefaKVk6Yms0Kw0IAbr0nn8amW3js+ddISojjzIW6c7C70s6kjFiiLfpHJ8OSQIYlge2OEircDd3O3S2alwZXAwlWdZJ7WOPzQOUOsMQElXo5PSav90FhYER6+F5NaxV7gN+/+Q/eeOcjLrvgDC679hcgDGhGKydPTVeevaIVs9FAdqK+EetOKmDJXdcz9/Bp3PfoMr74ZmPruJ0VdjSt42Zc0J5+xWY2VKzHo3mCNVsxyDR7mqk8uAo8jqDEfowtlVRLH09xh8iIFPwt7ZqO/+2f/+G5P/yVs08+jtuu+QGu6CwAZuYmqGJWii5My44nPkqvjOnKnMuTD9zKxHF5/OLBp9i4te0gVWFFz3V3+oWznurmrmWVFUOfHRUb+aZ0Fdt6qnXfC9nWpNbHPk3D6fF1mw4cKmERfCHEGUKIQiHEbiHE4gDXhRDiGf/1TUKII8Kxbm/8e8Vqljz7B044ejY3L/kdTalHII0WZucnkh6vMnEUgZk3LpmTpqTjtqXiLDiHB19+n7SUJG6851F279NL2bq8GnuqOqbgHZ84JbgFa/dSuOXPSGff++0qhgDNtZTv+wzKgt+HmR6Th7Vd6e3NJY3sKG/qUPolnIQs+EIII/Bb4ExgGnCpEKJzk9czgQL/zyLghVDX7Y1Va77l7qW/Zdb0yfz6npsgJgWPLZl545NJibVGennFMMdoEBzKwkxKS2fpy3/DZrVy7eIllJTpt+1NTm+HujtGYWBGbH7Qa64ofFuPA5duaD2Rqxi6VO/5T0ivz7YmdQjllDZ0X7gvXITDw58H7JZS7pVSuoG3gPM7jTkfeE3qrAYShRBZYVg7IJu27eK2B59ifH4uTz94Gzarfox+SlYc8bZeGlkoFH6OnZDa+jh9zCSeX7IYt8fDtYuXUFOnb9BuL2vs8EFNMYdQM7+ljtLy9dBUDvtWBD+PYkDY4gitcck4W1rr4xaPj0p/e9WWFifFZRUhzd0d4RD8HKD9Oy/2P9ffMWGhpqaGG+95jLTkRB777UvExcYgERiNgtwkVTJB0XcOpWtOztS9sMRjfsDTD/2Cqpo6rr/r163F1iobXdS3tJ3G7a3hTk/sbC5jed02WlRFzSGLlJJ1+/8X0hzz4idi9ufc213e1l4MHq+XXzz0NFfe/AB2e/gL7IVD8APtfHY+T9yXMfpAIRYJIdYKIdZWVfXvlBpASkoK1998Ky8svZPovBk0pMwmbeapLJiU1vuLFYoA5CX7HQWDkfyTLuex+25m974ibr73cZwuXZj3Vzcj/TXwo40WZseNDWnNKk8jsrIwpDkUkeHLHW/TVL4+pDlsfrGXUrK7Uhd2TdP45WPL+HLNt1z7o4uJjQ2tw1ogwiH4xUD7JNJcoHPB576MAUBK+aKUcq6Ucm5aWnAiveCS67DOuQSEwJqQRm56isrIUYTEKf6yG5opmuPmzeLB269lw5ZCFj/8LF6f3hrx2+KGVtFPMEWTb+tjRc0A7G2pZH/paij8KHTjFaGzdwUUfoS24wO8DQd7H98Dc+LGYRAGmt1evi3WQ4NSSh574XU++u8qbrj6Er5z7hnhsLoL4RD8NUCBEGKcEMICXAK832nM+8Dl/mydo4EGKWU3jR7Dy9Hjg//QKRTtOST6XlMsZy48lsXXX8GKr9bx4G9eQvPXzz/0AQYYZwvtUN8BZzUlzlpd9Lupz6+IMC114LKDRw/fuWXPfY97w4Bo7dWxs6ItZPPSG+/y1nuf8IOLzuKKSy6gPnVuSOt0R8gnbaWUXiHE9cAngBH4vZRyqxDiGv/1ZcCHwFnAbqAZuDLUdXtiVn4iheVNzBmT1PtghaIfnDItg//IuQipcfJPjqOuoYllr71NfHwMt/30Bwgh2FhUT35yNMkxFo5NmMSqhq6NUPrKrpZycmzJsOsTiE7VC6+NORbMwVfqVPSDg6tbH3qlj9UNu0Ka7tCmfvsmTH/956e88Me3OffUE7jyrt9Qa4ucboWltIKU8kN0UW//3LJ2jyXws3Cs1RdSY62kTlSpl4rIcPK0TD7brqdmXnj7s1R4Y3njz6+SFB/H1Zd9B4CDtc0kx1jCUgxLSqmHJA8dzqraoXue+cco4Y8kjW1BCE1qfFEf+p5KvCkaR7t6TP9esZqlz77KCUcfwc8eewPNrGcRHjsxMpGJEXnSVqGIJEIIFkxOa318zZ2PcPbJx/HcH/7K3/7Zlpvd3osLhXJ3p3mayvVOSnuXw77PVbgnnHhd4GzUe8+WbWx9Oqg6SZ0Yb8ugutLALv8m7VdrN7U7K3QjRr/Yj0+Laa3VFG5GbPE0hSKSmIwG5k9M5cvd1RgMBn7580U02h0sefYPJMTHctqJRwO66BckZ7PLHXzj6sLmMrKs3dzmux16mCcqMej5Fe3Y89+ATxc2h7blOCEqg5oqI4dyRzZv381tDzxJ/oQpPPTMiwhzW1esSPbSVh6+QhEkURYjM3L1KpcN+adxxzNvMXPmDO5e+lu+WttWC72qFjLNoVXDrOzJwywNLUVQ0T2a1Pi8bnvI85icbaG3PQeKueHuR0lNTmTJ797Emj6BpqRpJMVYOGVaRkQzCpXgKxQhkNGuJpM1NoH7X/gL4/NzufWBJ9m0rW2Dz1kXmte2zVHS/UWvSy/JIFWf3KBxVEPhx12e/rx+B1oI/YellMQ0pVPRpJ/XKKus5md3LsVkNvPwK/8gOU3P5JozJmlAkkyU4CsUITIpo60eSmxCEk898WvSkhO54Z5H2bO/uPWawR1aIoGzpxLKdfv0htkHv9aFq6k8pLVGFWXfQvEaguk/2xulDU7MQo+c19Y3ct3iJTha3Dzy0ttk5Y9tHZcUYwn72oFQgq9QhEh+SnRrjj6AeeKJvLD0TqwWC9cuXkJpuX5i3OqMx+kJPo97dcMuaro0yW5H1Q5oqQUklPWvvd6ooXp322G2lno9x74x8P7K3pb+17ZvT4vHR47Q/180OZq5/q6llFVU89DzbzB+ymGt446ZMHBnhZTgKxRh4uSp6VhMBhAC2+yLeOKZZ3G53VxzxyPU1DVgNZjJdOXhCyGrZrP9IB6tmy8NT7sKmyEeEBqx1PjDbIUfwcGvYP/KgMM22w9y0Bl8jwKHy4vPbiPRGEtLi5Mb736UXfuKeOy+m5k+5ygA4qPMnDItgxjrwOXOKMFXKMKEEKLVW5MGM3lTZvPMr35BVW09P7tzaWuxtQM1Ldid3qCF/8uGQjSpUjEjSY93Ur3gcHmpaHSRa07D5XZzy/1PsGn7LpbcdQNHnHJR67h545LDYWq/UIKvUIQRs9HQGt6RBjOHT5vE47+8Tc/MuOtRmlucAFQ2uThQE3z98z4dAir8CGr36T+jvTxD4Uew/4teh2lSY5ujuNdx3XGgxkFFo4sCSw4er5c7fvUsX6/fwuJ77mPhicfR4i+q1z4EOJAowVcoIsDJU9PxWuJpSJnFtDOuZMldN7ClcA833/c4Fk/bLXxDS3C9bDUk5a566jyOngdW7dB/ADpv+jbXQnfhoZGCvaotRu/qvS3lVw27qHQH123K5fFxqDNhNDbue3QZK75axw33LOGk711LbeZxABw/KbWHWSKLEnyFIgIIIThlWgYeazI+Sxyz/u92bn/kWdZ+u53fLnkLg0/Pta6xu2lxBye6O5pL+dZ+oO8v2PNfPf1Q0/SOWkVfQ/mmkef51+3XN2ObyqFkbZ9bEDb73HiC3Ptodnspqdfv3qSUPPz0K3z8v1XcePUlnHdpW+mwgoxYrCZjUGuEA3XSVqGIIHnJ0RTVNoMwsPCcC/G2NPLY/Yux/NrMNYsvpsXgpqzBSX5KFCbDAPhfxWv0P9P8/XebysHrhvyjIr92pKjbDzHpYPH3Lajs/0GpWo+dTfbgyh5rmqS8QT8pO8GcxYsvvcO7H/2PH1/2HS696hoO7QacODkNs3FwfWzl4SsUEWRyZlyHeO1pF1/BzXfczeer1/Pi42/j89fSP1jTgscXnKf9VcMu3Jq394HtORTmAX8q5zDF59UF/sAq/fcg71aCFXuvprG/Rt+Ml1Lyp9c/5M/vfMwlF57DdVdcjOZvUD4pI27QxR6U4CsUA8L4tLaTtmdffgM/u+UXLP9iHX96+l9ofqEvqm3BFUSevkvzsKphJ66eDmb1RuFH4Il8E+3w4z8spXn0JiW7Pun3DN829SMs1o4mp5eD/o13KSX//PNyXv7ze5z13e9z5f2/w5Ewiea4cUxIjyU/ZWi0V1WCr1AMAOPTYjscnf/Oj2/jypvu5PPl63nt+X+2NlApqXfS7O6nt+6nzFWPN5T8+0OhEGejHuv3OHUPerjQ/hxCH5BS4pMadd5eNr4D4HB5qWpqK3j2jzf/xz/+spxzzv8ON/3yMQxGI87YPGbmJUW0GFp/UTF8hWKASIqxEB9lptGfmXPZoptwu5y8sexJhBD88LpzMBgMlDe4yEkUWM3929zb76xiv1M/1Xtk/HhijLZeXtEJe0XglorJ4yFtcv/mCieaBvX7IXEsGAzg8+hfSjW7IGlc0NMWNpd1LT3dB7w+jYrGNrF/763P+edbKzjvjBO59a778Pj3YuaNTybeZg7avkigBF+hGEDmjUvG7dX4fKcuzD+6/na8Xg9/efk5fD4fV1x/PgajgZJ6J+nxVmKDPIW5pnEvC5Kmhcfo2r0QkwZGsx4rN9lg/In6tZY6sCVCJHtGNxyEqkJoKIaUiR2zblrqgp42GLGvanLR5Gy76/nHm//jn28u56zzL+S+6y6iPkpPuZydnzjkxB5CFHwhRDLwF2AssB/4Pylll38BIcR+oAnwAV4pZWQaNioUwwCLycDcsUms3V+HEIKrb76bOGcxL//pPXw+jatu+g5Go5HKRhfEQaxtCPhlRd/QGi/3NEPpBojPgZJ1+hfBxFP6No+Uev1+SywYAtzBNJZCbEbHa16/N+129DnFsicavM0Y6P8XVEldCy5v26bw+28t5/03l3Pad77HTQ89SZ3fsz9qfDJxQ1DsIXQPfzHwmZRyqRBisf/3O7oZe5KUMvjiFArFCCIx2sK4tBj2VTkQQnDJ7c/gShjH6799Es2ncfUtF2IyGalsctHk9JCVONitDDtVkmwqb6vI6eu0Wexp0VsvOqr1uwGDCeoPQu0eXejd/kTF2AzImqWHaUDvMlWxRX88+UxoqohIrf8NTfv7/RqvprWKvZSS9974L//66+ccf+7Z3PbQkxj87yEzwTZkxR5CF/zzgQX+x38EltO94CsUinZMSIslNdbKmn21SKOF7/30Rpo8Rbz34tv4fBqLbrsIk9lEi0djb5WjQ6ZPX9jfUsXYqLQIWd+Jwo90AY9K6pjyCWCNazvl6m5Xo8ZeoRcwyz9GF/1DYn9oviFCk7Ntg1ZKyVsvf8x//rmao889hdsefKpV7JNizEzPCa3RTaQJNUsnQ0pZBuD/M72bcRL4txBinRBiUYhrKhQjhoQoM9l+791ksHDq987gkqvPYN2qbbzw6N/wuNu8571Vjn6lbe53VrG8bhs7HMG3V+wX9oquYg89lzRwNeqplNW7I2eXH4fPya7m/vUJcHl8rWKv+TT++Nz7/Oefq7nsu+fx4JLXibHovRAMBpidF/kGJqHSq4cvhPgPkBng0t39WGe+lLJUCJEOfCqE2CGl/Lyb9RYBiwDy8/P7sYRCMTyZlh3P2NRoVu2uYVL6mWR9PxpvfC5vP/kyTz34Z66/6xKiovXmKSX1zn57+uXueqbEZEfC9PBRs6v3MSFg9zpZ27S3z+OdHh+l/lIJAF6vj5effIc1K7dw4aULueUnV1Hv36gemxrNxPS47qYaUggZQls0IUQhsEBKWSaEyAKWSyl7zN8SQtwP2KWUj/c2/9y5c+XatWuDtk+hGE54fRrLC/XsnTpXBf947/e8+fBz5EzI49Z7LiY+MbbD+P4I//GJUzCK0XHsxic1DjqrGWNLRSBYUd+/Ugstbh9lDW1iXxc1ntfv/hXfflPIeddcyuILTsNsSaE+fR7j02IYnxbbw2wDjxBiXXeJMaH+D3gf+JH/8Y+AfwRYPEYIEXfoMXAasKXzOIVitGNqV1o5yZrB3NMXcPXSOyk/UMojd/+JsobgP64r6wOEWkYoRc4aDjir2dtS2W+xBzqIvcPewvM/f4Rvvynkkhsv46QfXAxxBTQmz+CIMUlDTux7I1TBXwqcKoTYBZzq/x0hRLYQ4kP/mAzgCyHEt8A3wAdSyq7dghUKBdDWI/ewhKOZduwcrnnqfuz1jTx261OUHGxru7e3yoHd2feTsGsa94Td1qGIhp5N02vp6E54ffrm+CFqqxp45K4/cWDbLi5/4DaO+t7FmAxmHAkTmT4uk+QB6kMbTkIK6UQaFdJRjEY0TbKhqJ46h5taVwWlLXsp3XOA393yAF63h6uW3M7s/I6f276Gd9It8WRaEkk2Dy/PtD/saa6gyFVDjNGKw+fqdfy+KkeX9uXF+yt46oE/0eL0ctWSxRTMmUFB3CwsBhunTMtARPKgWYhEMqSjUCjCjMEgmDMmiSPGJJFs1UM82RPGcOOyJcQmJfDCzQ+y6n8dDyCVNfSt8FmluzHoypDDBemX776IvdurdRH7HZv2sXTxKwDc8MIjFMyZAYDVGMWs/KQhLfa9oQRfoRiiJMdYOGVaBhNidcFJyc7gpt8tZdzMKbzy5Dv89e0NHLpDb3Hr4QhN69sd+8amA/0vqTyEcGteqt1d0z2llH0K5UgpqXW4Ka7r+EX5xX828MT9r5OUEs/i31xH9oQxrdemZseTFmcN3fhBRAm+QjHEOfOwcRTEzQIgOj6Wnz5xH/POXsgnr73H88/8u0Ou/v6a5j7V1a/3OljVsJONTQfwDEPh/9Z+gC2OInztmrkXOWtYUb8dh9azZy+lpKTOSX1z29+b5tP4yysf84dn3mPCrOnc8MLDWCfMab2eEWclZ9BPO4eOEnyFYohjMho447A8TEI/NmMym7nkzus5+6c/YP1nq1j6y79Q2c7ZLaptoaLR2c1sHan3OviyYSdVQfZxHSxafG6gLXxT5KxhT0tFr6+TUrKvuhl3uy/FZoeT3/z6ff79j684+dyjWPSbe7ElpSEN+t/3+PgJnD/lxAi8i4FHCb5CMQwwGoxcf8w5JJr1aoxCCE65/CKu+NXtlO4t4pEbnmTDQUFznF4u2OHysbfKgcPVN+99q6OYYmctDd7+1ZQfbJp9bjSp9UnsHS4v+6o7vr/ykmoe+cVL7Fyzie/feDEX3XQ1xnY9Z/OSozl72nSizUOjgUmoqCwdhWIY4fX5eHnNZzi1NuEq23uQP9z1a2pKyznv+itYcOFpxDUUtl6PsRpJj7P2ebNxZmz+kM/i+bxuO123W7un82EqgK9W7+X1J9/EbDbxo4fvZOIR01uvWQ1RFGSZGRufz8SkiWGzeyDoKUtHCb5CMQz5/YYPWhupALTYHfz5V8+wZeU3HHHq8fzfz39Cmmtfh9f0NXXTiIEpMdmkWeLDanM4WV63rU/jXB4fJe1LJJjjaIzK5f3n/sjKtz9g3LQJXP7Az0nMbqsekxZr5fszz8Tj82AymIZdVo5Ky1QoRhjjU2MoyGjzwqNiY7jykTs4a9H32fDZlzx+1e1srrDRmDKTxmQ9y2dvlaNPG7o+NLY6imn0Bk71tHudAx7z39tSyRf1+l1Lrcfey2gdTZMdxB6grF7y7HV3s/LtDzjxe+dy3fNLOoj99Jx4xqXqFS/NRvOwE/veUB6+QjEMWV+xnkZ3Ix6fxtbSjuK7Z+M2/vTgkzRW13HmTy5l4WXfwaw5iel00jY3KQqLqXufb0ZsPintQjtuzYvFYGr1rsPWUasXmn1uvmnsezXNykYndldbVVF74mQ0aeCbT1byzlMv6/0H7ryew086psPrJqbFEmszcWLuicNa6JWHr1CMMI7IOAIAs9HA9Jx4DO30acKsafzi1SeZceJRfLDsT7xw8/1UVbflpmumKLzmWIrrWnr0+DfbD9LsP7zU5G1hVcNOyl31EXk/ndGkhkfTRbuvYn+oNEJ7sQdobHTy+3t/w5uPPEvupHH84o9PdhH7WXmJrZ3FhrPY94by8BWKYYqUkhXFK1p/31hU3+X6Nx98xrtPvYIEzll0GcddeAbCpHdksjlKsTirMRogM97WbdP0iVGZlLvrsfucpFviqfSHc45JKMBqiEx3p68bdtOiuTkqfiJf90HwK5uc2J0dhd5nsLDm6928+dzbtNgdnL3o+5z4vXMxGNveZ2KUmfzkaAz+b8zc2Nxht0nbGbVpq1CMUDw+D1+WfgnoMeu9VXbs7o7CV1dexV8fe4EdqzcwdsYULll8HRlj8wCIr9nUOi4/OQqTsX83/TFGK3PixmEIsvSyR/NiEsYOXrVb87KqYWefXt++2BlAS2weUfYiDrYk8rdnXmfbl2vJKRjHZffcSPbEsR3GTkiNIS6q4xfW+ITx5McP7z4cSvAVihHM8qLlHX4PFNeXUrL2kxW89/QruJqdnPB/53DaFReT5j6A0Dr1pKX3+H575saNJ8poocxVT7olHouhb51TPZqPLxsKybOmMCE6o+399CEDp9bupr6lo92a0UZ9zDhWvv0hH7/8JgBn/uRSjv/uOR1y66Gr2KdHp5NsSyYjemgXRusLSvAVihFMZ8GHwKIP0FRbz7+Wvc43H/yX2KQEzl10KcecOpeY5pKAc2cl2DAaRI/in2qOo9qjH/WNNdqYHpuHrVOox6N5+aZxD5mWRGJNNuo9DtIs8R0KucUYrD2WRah1uDuUQ2iJzSfKfhBH/AR8BgubVq7lny+8TnVxGdPmz+WiW39CcmbHrqtGA0zNisdkaHs/U5OnkhGTwUhBCb5CMYLxaB6QcLDpIEVNRR2ueX0aWwII/8Htu3j3qVfYv6WQ7IljOeuKC5gxfw4YLRh9LQjNg8ndgMmfAhlrNRJrM2E2GjD3IewzPSaPRHM0JmFEkxrVnia2OQJ/qfRE+wNT0mDucDfSmDIToXnZu3U3/3z+NfZt2k7G2DzOv/5HTDn6iC6e+vTs+IAhq+GeldMZJfgKxSghkLcPehngbWVdwzwb/vMFH7/yJlVFZeROHs/pV13CYfPndhDAmMbdGD1tJ3ujrUakJomLMmN3ecnoxyne3nB7NewubwdPHnRxBxCaD1tzCS22LHZvKuTfr/6VXes2E5uUwJk/vpSjzjmlS/gG9CycQGTFZDE5uceurMMOJfgKxSihuqWabTXb0GTgdMtah5uDtR3ryfi8PtZ9soJPXv0rtaUVZIzN44SLz2buGQuw2NrKAZvdjfiEGc1kw+htxmeKIrZhF1bc5CXrtWY0TSKRGA1dPWmvptHU4iXJ3ynK7vRS2eQiI96K3eXF4fLhikrD2lKFZrBg0Nytrz0k+D6vl00rvmbl2/9i36YdxKUksfCy8znm/NOxRtm6rNk5C6c9seZYDk8/HHOEMo0GCyX4CsUooztPH6Cotpkah7vL8z6vl/WfrmTFX/9Fyc69RMfHMu+shcw9YwE5BeMCTyYlUY5iXFHpxNYXdrjkiJ/Q5bBXIA6JeXuE5kNIH0avA4PPTVmD4OsPPuOr9z+lqaaO5OwMFnzvPI4692Qs1q416rMTbcRaTdhMxoBiPydjDnGWuF5tG45ETPCFEBcD9wNTgXlSyoDqLIQ4A3gaMAIvSymX9mV+JfgKRXD0JPjQFjrp7O2DHurZ++12Pv/bv9j6xRp8Xi9ZE8Yw5/QTmX7ckaTn5wQM4cTVbe8SY2+f9umxJmF21eGIn4DFVYtmMOEzxeLtRnjrq2rYtPwr1n+6kgNbdyKEYMrRRzD/wjOYetTsDvn07ZmSGYetmzMFALPTZ5NgTej2+nAnkoI/FdCA3wE/DyT4QggjsBO9yXkxsAa4VErZa+6VEnyFIjh6E/xDSCn5trih2+v2+kY2fvYFaz9ewYFtem58SnYGU4+Zw4TZhzH2sEkkpuslm9E0BJK4uq1AR8GXBhNNvZRisNc3cnDbLnat28SOrzdQvk/fgM6eOJbZpxzHEaccT3JWesDXRpmNTM7s3WPPi8tjQuKEXscNZ3oS/L4lzHaDlHK7f4Gehs0Ddksp9/rHvgWcD/St3J1CoQiaWHMs9h6KjQkhiDIbafH4iI8yd6jACRCbGM9xF53FcRedRV15Fdu+Wse2r9bx9b/+wxd//xCAxPQUsieOJS0vm7S8bLLivSTEmMCVgcHVjMEg0CwxNDfacTqaaW6y46hvpLa8kpqSCqqKyygu3ENtWSUAJouZ8TOnceSZJzHt2Llkjsvr8T3aTEYK0vtWznmki31vhCT4fSQHaJ8rVgwcNQDrKhSjlmOzj0UisRqtvXr7E9Nj8Gm05tp3LtFwiKTMNOZfcAbzLzgDr8dD6e797N9SyIGtOynfV8SudZvxuLruDfSE0WQiOSuNvKkTmX/BGeRNmciYwyZ12CwOaEu0mTEpfSv3fIhJSZP6NX4k0qvgCyH+A2QGuHS3lPIffVgjkPvfbRxJCLEIWASQnz+8jzgrFIOFxWhpfXx01tEcaDxAmaMs4FijwUD79PTJmXEUljeREW+lojHwQSiT2Uz+1ALypxbAxfpzmqbRUFlDQ3UN9vomHA2NOB3NIEFKDRBExcUQFRtDdHwsyVnpJKaldBuLD0SsxYjd7SM/ue8dqCYmTmR3/W6yY7P7/JqRSq+CL6U8JcQ1ioH292S5QGkP670IvAh6DD/EtRWKUY/NZKMgqaBbwe9MlNnYmrceZzVT2+yi1tG1/EJnDAYDSZlpJGWmhWJuF1JiLKTHW7EGyK/vjSRbErlxueTG5YbVpuHKQIR01gAFQohxQAlwCXDZAKyrUCj8GISBBXkLqGyuZFtN37fPYm0mYm0m8pP1372axpaSgWl+kp8cTXKMpfeBPZAdo7z69oQk+EKIC4BngTTgAyHERinl6UKIbPT0y7OklF4hxPXAJ+hpmb+XUm4N2XKFQtFvkqxJIb3eZDC0ev/Nbi91Dg9VdheHZccHrN3TV6wmA1mJNgSChKjwHYRKiw7v3cZwJ9QsnXeBdwM8Xwqc1e73D4EPQ1lLoVCEjtnYJqbjE8azt2Fv0HNFW0xEW0zkJEUBbeUL7E4vbp+GxWjwZ/+YsJqMVDe58ElJRnzbidg6h5toqzGocE1vpEUpse/MQIR0FArFEGJ66nSiTdFEm6OxGq1sr90e1vkPdY7q/Dg1rmvmTVKIIZtApEenMy1lYNovDjdUi0OFYpSRGpVKtFnPchkpZYETrYkA5MTmMDV56uAaM4RRHr5CMcopSCpgV92uwTYjJGamzcQnfZiEaUSVOg43ysNXKEY5WTFZpEalkhfX84nWoYxBGDAbzErse0F5+ArFKMcgDExPnU6Dq6FLA5Whzvyc+RiU39pnlOArFIphyZj4MSOuln2kUV+NCoUCaCuCmBqVypTkKYNsTc/EmGMYl9BNjX5FtygPX6FQABBviWdS0iTSotMwG8yYDCa2VG8ZbLMAPdUyzhJHlCmKmpaaYb3fMJgowVcoFK20LzCWGpXK0VlHY/fYiTZF8035N4NmU/tKl6lRqYNix0hAhXQUCkW32Ey21rz943OOHxQbVOgmfCgPX6FQ9AmjwcjMtJlsqtqE1WjF5QtcOjlY5mXOA8Dtc9PobmRvw14mJU1SG7NhRAm+QqHoM8m2ZBbkLQDa2ihmx2bj1bxkRGewuXpzn+dKj06nslnvcjUnY07r6d9oczSJtkTy41U/jHCjQjoKhSIoDsXVJyRMYFrKNFKiUlqvHfpSOESSrWOVzinJU5iWMo35OfOZlDSJuG4amSvCi/LwFQpFUGTHZnfpInV8zvFIf0O743OOxyd9eDUvQIdN30NfDmaDWXWiGkCU4CsUirBhNBg7PDZixGK00OxpBvTWiwWJBSouP0gowVcoFBEn2hzN2PixZMRkEGWKGmxzRi1K8BUKxYAwNmHsYJsw6lGbtgqFQjFKCEnwhRAXCyG2CiE0IcTcHsbtF0JsFkJsFEKsDWVNhUKhUARHqCGdLcCFwO/6MPYkKWV1iOspFAqFIkhCbWK+HVBNBxQKhWIYMFAxfAn8WwixTgixqKeBQohFQoi1Qoi1VVVVA2SeQqFQjHx69fCFEP8BMgNcultK+Y8+rjNfSlkqhEgHPhVC7JBSfh5ooJTyReBFgLlz58o+zq9QKBSKXuhV8KWUp4S6iJSy1P9npRDiXWAeEFDwFQqFQhEZIh7SEULECCHiDj0GTkPf7FUoFArFACKkDD5qIoS4AHgWSAPqgY1SytOFENnAy1LKs4QQ44F3/S8xAX+WUj7cx/mrgANBmpcKDOesoOFuPwz/9zDc7Yfh/x6U/f1njJQyLdCFkAR/KCOEWCul7PZswFBnuNsPw/89DHf7Yfi/B2V/eFEnbRUKhWKUoARfoVAoRgkjWfBfHGwDQmS42w/D/z0Md/th+L8HZX8YGbExfIVCoVB0ZCR7+AqFQqFox4gTfCHEGUKIQiHEbiHE4sG2p78IIX4vhKgUQgzLswpCiDwhxP+EENv9lVRvGmyb+osQwiaE+EYI8a3/PTww2DYFgxDCKITYIIT412DbEgzDvcquECJRCPG2EGKH//NwzKDbNJJCOkIII7ATOBUoBtYAl0optw2qYf1ACHECYAdek1JOH2x7+osQIgvIklKu9x+4Wwd8Z5j9GwggRkppF0KYgS+Am6SUqwfZtH4hhLgVmAvESynPGWx7+osQYj8wd7hW2RVC/BFYKaV8WQhhAaKllPWDadNI8/DnAbullHullG7gLeD8QbapX/hrDNUOth3BIqUsk1Ku9z9uArYDOYNrVf+QOnb/r2b/z7DyjIQQucDZwMuDbctoRAgRD5wAvAIgpXQPttjDyBP8HKCo3e/FDDOxGUkIIcYCs4GvB9mUfuMPh2wEKoFPpZTD7T08BdwOaINsRyj0ucruEGQ8UAX8wR9We9lfWmZQGWmCH6gw/7DyzEYKQohY4O/AzVLKxsG2p79IKX1SyllALjBPCDFswmtCiHOASinlusG2JUTmSymPAM4EfuYPdw4XTMARwAtSytmAAxj0PcWRJvjFQF6733OB0kGyZdTij3v/HXhDSvnOYNsTCv7b8OXAGYNrSb+YD5znj4G/BSwUQvxpcE3qP+2r7KLX45o3uBb1i2KguN2d4dvoXwCDykgT/DVAgRBinH+T5BLg/UG2aVTh3/B8BdgupXxisO0JBiFEmhAi0f84CjgF2DGoRvUDKeWdUspcKeVY9M/Af6WUPxhks/rFcK+yK6UsB4qEEJP9T50MDHriQqg9bYcUUkqvEOJ64BPACPxeSrl1kM3qF0KIN4EFQKoQohj4pZTylcG1ql/MB34IbPbHwAHuklJ+OHgm9Zss4I/+rC8D8Fcp5bBMbRzGZADv+tunHqqy+/HgmtRvbgDe8Dufe4ErB9mekZWWqVAoFIruGWkhHYVCoVB0gxJ8hUKhGCUowVcoFIpRghJ8hUKhGCUowVcoFIpRghJ8hUKhGCUowVcoFIpRghJ8hUKhGCX8P7ui8PSQZdgZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "# number of observations\n", "n = 1000000\n", "# x coordinates for initializing the sine curve\n", "coord = np.linspace(0, 2*np.pi, n)\n", "signal = np.sin(coord)\n", "\n", "# error i.e. epsilon of the three synthetic time series\n", "sig_err_x = 0.02\n", "sig_err_y = 0.07\n", "sig_err_z = 0.04\n", "err_x = np.random.normal(0, sig_err_x, n)\n", "err_y = np.random.normal(0, sig_err_y, n)\n", "err_z = np.random.normal(0, sig_err_z, n)\n", "\n", "# additive and multiplicative biases\n", "# they are assumed to be zero for dataset x\n", "alpha_y = 0.2\n", "alpha_z = 0.5\n", "\n", "beta_y = 0.9\n", "beta_z = 1.6\n", "\n", "x = signal + err_x\n", "# here we assume errors that are already scaled\n", "y = alpha_y + beta_y * (signal + err_y) \n", "z = alpha_z + beta_z * (signal + err_z)\n", "\n", "plt.plot(coord, x, alpha=0.3, label='x')\n", "plt.plot(coord, y, alpha=0.3, label='y')\n", "plt.plot(coord, z, alpha=0.3, label='z')\n", "plt.plot(coord, signal, 'k', label=r'$\\Theta$')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In approach 2 we can estimate the triple collocation errors, the scaling parameter $\\beta$ and the signal to noise ratio directly from the covariances of the dataset. For a general overview and how Apporach 1 and 2 are related please see [Gruber2015].\n", "\n", "Estimation of the error variances from the covariances of the datasets (e.g. $\\sigma_{XY}$ for the covariance between $x$ and $y$) is done using the following formula:\n", "\n", "$\\sigma_{\\varepsilon_x}^2 = \\sigma_{X}^2 - \\frac{\\sigma_{XY}\\sigma_{XZ}}{\\sigma_{YZ}}$\n", "\n", "$\\sigma_{\\varepsilon_y}^2 = \\sigma_{Y}^2 - \\frac{\\sigma_{YX}\\sigma_{YZ}}{\\sigma_{XZ}}$\n", "\n", "$\\sigma_{\\varepsilon_z}^2 = \\sigma_{Z}^2 - \\frac{\\sigma_{ZY}\\sigma_{ZX}}{\\sigma_{YX}}$\n", "\n", "$\\beta$ can also be estimated from the covariances:\n", "\n", "$\\beta_x = 1 \\quad \\quad \\quad \\beta_y = \\frac{\\sigma_{XZ}}{\\sigma_{YZ}} \\quad \\quad \\quad \\beta_z=\\frac{\\sigma_{XY}}{\\sigma_{ZY}}$\n", "\n", "The signal to noise ratio (SNR) is also calculated from the variances and covariances:\n", "\n", "$\\text{SNR}_X[dB] = -10\\log\\left(\\frac{\\sigma_{X}^2\\sigma_{YZ}}{\\sigma_{XY}\\sigma_{XZ}}-1\\right)$\n", "\n", "$\\text{SNR}_Y[dB] = -10\\log\\left(\\frac{\\sigma_{Y}^2\\sigma_{XZ}}{\\sigma_{YX}\\sigma_{YZ}}-1\\right)$\n", "\n", "$\\text{SNR}_Z[dB] = -10\\log\\left(\\frac{\\sigma_{Z}^2\\sigma_{XY}}{\\sigma_{ZX}\\sigma_{ZY}}-1\\right)$\n", "\n", "It is given in dB to make it symmetric around zero. If the value is zero it means that the signal variance and the noise variance are equal. +3dB means that the signal variance is twice as high as the noise variance.\n", "\n", "This apporach is implemented in `pytesmo.metrics.tcol_metrics`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error of x: 0.0200, true: 0.0200\n", "Error of y: 0.0701, true: 0.0700\n", "Error of z: 0.0400, true: 0.0400\n" ] } ], "source": [ "from pytesmo.metrics import tcol_metrics\n", "\n", "snr, err, beta = tcol_metrics(x, y, z)\n", "print(f\"Error of x: {err[0]:.4f}, true: {sig_err_x:.4f}\") \n", "print(f\"Error of y: {err[1]:.4f}, true: {sig_err_y:.4f}\")\n", "print(f\"Error of z: {err[2]:.4f}, true: {sig_err_z:.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimation works very well in this example.\n", "\n", "We can now also check if $\\beta_y$ and $\\beta_z$ were correctly estimated.\n", "\n", "The function gives us the inverse values of $beta$. We can use these values directly to scale our datasets." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "scaling parameter for y estimated: 0.90, true:0.90\n", "scaling parameter for z estimated: 1.60, true:1.60\n" ] } ], "source": [ "print(\"scaling parameter for y estimated: {:.2f}, true:{:.2f}\".format(1/beta[1], beta_y))\n", "print(\"scaling parameter for z estimated: {:.2f}, true:{:.2f}\".format(1/beta[2], beta_z))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABgyklEQVR4nO2dd3hcxdWH39mm1ap3yZZkuchF7rjhgm0MoZeQAIEQwJhQDZjei4HQewumxhBCL4EAhkDAHxiDCy64yEWu6n13pe3lfn/saoslN+2qz/s8erR779yZWZVz554553eEoihIJBKJpPej6uoJSCQSiaRzkAZfIpFI+gjS4EskEkkfQRp8iUQi6SNIgy+RSCR9BE1XT+BApKenKwUFBV09DYlEIukx/Prrr3WKomS0da5bG/yCggJWr17d1dOQSCSSHoMQYs/+zkmXjkQikfQRpMGXSCSSPoI0+BKJRNJH6NY+fIlE0j1wuVyUlZVht9u7eioSP3q9ntzcXLRa7SFfIw2+RCI5KGVlZSQkJFBQUIAQoqun0+dRFIX6+nrKysoYOHDgIV8XFZeOEOJ1IUSNEGLjfs7PFkKYhBDr/F93R2NciUTSOdjtdtLS0qSx7yYIIUhLSzvsJ65orfAXA88Dbx6gzY+KopwSpfEkEkknI41996I9v4+orPAVRfkBaIhGX5JOoKkK3M6unoVEIulkOjNKZ6oQYr0QYokQYuT+GgkhLhVCrBZCrK6tre3E6XVTTOXg9UavP5cdKtb6vhQFvB4wV8LWJeBxR28ciUTS7egsg78GGKAoyljgOeDf+2uoKMrLiqJMVBRlYkZGm9nBfYemaqj6DepLoten4r952BqgeiNs/y9UrvMdK/nG993r9d0MWjDu9d0QZLEciaRH0ykGX1EUs6Iozf7XXwJaIUR6Z4zdo/G6fN/d7Q+Fs7qsOD1+9421AXb9X+Ccq3EPdc4mvIqXGqcJr+Jl065vaCr+N5St8hn5rUugdqt/Pp52z0MiiYRVq1YxZswY7HY7FouFkSNHsnFjmzEikgPQKWGZQohsoFpRFEUIMRnfjaa+M8bu0RzuitrRDHuWQ8EM0BnAaWHl3u8QMYnMypsFpSuodppI0cShU2nYaCnF5LaSo0um0mkMdFMLTBYqDOoYvIoXFK9vZVC6ArJHgT4pih9S0tPYVt1Ek90V1T4T9FqGZiXs9/ykSZM47bTTuPPOO7HZbPzlL39h1KhRUZ1DXyAqBl8I8Q4wG0gXQpQB9wBaAEVRFgFnAlcIIdyADThHkcV0o45r5/fUu5rJrt0CzdXssdWCvRYlNoVyQzbbGzcH2iZrDFg9DgBqXU2t+lpp3kGmLpEapxmA2SlF4DD7bijDToTSlZCUC4n9OufDSfo8d999N5MmTUKv1/Pss8929XR6JFEx+IqinHuQ88/jC9uURJOtSyB1MGQMBWBd8x4sHgdbrBXh7WyNbN/+edgho9saeO1W2nbVtBh7gKWNm8nRJdPotsAWG4NcTlKaa9BWrvfdACR9hgOtxDuShoYGmpubcblc2O124uLiumQePRmZadvTadgBqQPZ0bQXi3/F3lEE3D7GXWwOOd6/No/CpnrQxGDvNw61So1Wdejp3hLJoXDppZdy//33s2vXLm655Raef16uIQ8XafB7CooCLivo4oLvW9jxHaWqrourL9/xNU5tAvFqPbsqfiImazRTBx7XZfOR9D7efPNNNBoNf/7zn/F4PEybNo3vvvuOOXPmdPXUehTS4PcUGnZC3TZQaSEmHuJ8IasexUujswkspV06vVpXU2AvwGHc69tAVmlAq+/SeUl6BxdccAEXXHABAGq1mhUrVnTxjHom0uD3FGyNvu9eFy5LHT/Vb2SA0LLHXte182oLhwl2/+h7PeRYUEv3jkTSHZAGvycgBBtMO7FZqknSGAK+9P3WMesGOLwuYlRaKPkWkvMhy59crSjgdsiVv0TSBcgCKN0Vj8uXCQtgKqO+qRSr1xkWL9+d+dm0HavHyU5bDetKfwieqN0CO7+XWj4SSRcgV/jdlYZdXT2DiFlpDkpCeBUvKqGC5hrfAY8TNLoumplE0jeRBr+bY/M4Wd/cnZ03h8YPK55Cr01gdGwGcWrpzpFIugLp0ulG1Nvqsbv9apYNO3B4Xawwl2D3RjeNvauwu5pYZd7pe2Nr9Gn7SCSSTkMa/G7EhroN/FL5C8tK/4+ljZv52bS9q6fUIWy1VPj2J0plaJ0keuzevfuw9XUWL15MRUXFwRtGmblz5/Lhhx8ecvv2fLa2kAa/u+Fo2q/UQW+h0mlkuXEbUk5J0tV0lcHvKqTB705Y6nzJVX0Ap+Lm/4zFuEylPveOp3e4rSQdw1133cUzzzwTeH/HHXe0KaDmdru58MILGTNmDGeeeSZWq08z6tdff2XWrFlMmDCB448/nsrKSj788ENWr17Neeedx7hx47DZbNx3331MmjSJUaNGcemll7a5KPnggw8YNWoUY8eOZebMmQB4PB5uvPFGRo8ezZgxY3juuecADqm/tubWcnzs2LFMnTqVF154IfIfIiC68ypr4sSJyurVq7t6Gp2CvbmGXzb+s6un0enk69MZFJvpeyNF2LotxcXFjBgxwvemphjs5gNfcLjoEyFzxH5P7969mz/84Q+sWbMGr9dLYWEhK1euJC0tLazNwIEDWbZsGdOnT2fevHkUFRWxYMECZs2axaeffkpGRgbvvfceX3/9Na+//jqzZ8/m8ccfZ+LEiYBPoC01NRWA888/n7PPPptTTz01bC6jR4/mq6++on///hiNRpKTk3nxxRf59ttvee+999BoNIF+9tff3LlzOeWUUzj99NP3O7eWG8esWbO46aabWLJkSasaAGG/Fz9CiF8VRZnY1s9RRul0MYrbRf3eH9lYs7arp9Il7LXXYXRZyNenk1610ZegJYtlS/ahoKCAtLQ01q5dS3V1NePHjw8z9i3k5eUxffp0AP7yl7/w7LPPcsIJJ7Bx40Z+97vfAb7VeE5OTpvjfP/99zz66KNYrVYaGhoYOXJkK4M/ffp05s6dy9lnn80f/vAHAL799lsuv/xyNBqfSW0x8gfrb+vWrW3OzWQyYTQamTVrFuC7WSxZsqTdP78WpMHvIhrtjayvXQ/Vm8Ft6+rpdClmj42NllKmavRojXt88fraWBgwA9TyT7TbcYCVeEfy17/+lcWLF1NVVcW8efPabCP2WSwIIVAUhZEjR/Lzzz8fsH+73c6VV17J6tWrycvLY+HChdjtravNLVq0iBUrVvDFF18wbtw41q1bh6IorcY+lP72Nzej0diqv2ggffhdRKl5LzTu7vPGPpSfTdv5wbjF5+d02cBu7OopSboRZ5xxBl999RWrVq3i+OOPb7PN3r17A8bznXfeYcaMGQwbNoza2trAcZfLxaZNmwBISEigqckn+tdijNPT02lubt5vFM2OHTuYMmUK9913H+np6ZSWlnLcccexaNEi3G434HMNHUp/+5tbcnIySUlJLFu2DIB//etfh/8DawO5fOoCXB4XDeUrfRWkJK342bSdaclDfQXc1TFQML2rpyTpBuh0Oo4++miSk5NRq9VtthkxYgRvvPEGl112GYWFhVxxxRXodDo+/PBDrrnmGkwmE263m2uvvZaRI0cyd+5cLr/8cmJjY/n555+55JJLGD16NAUFBUyaNKnNMW666Sa2b9+Ooigcc8wxjB07llGjRrFt2zbGjBmDVqvlkksu4aqrrjpofwea2z/+8Q/mzZuHwWDY7w3ucJGbtp2FolBtrSYpJplfKn+B8l+7ekbdmulJw6hyGsmNSUUMP6mrp9PnaWtzsLPxer0cccQRfPDBBxQWFnbpXLoLctO2m6Ls/YXiyuXo4nOgubKrp9Pt+cm0FQC9SouruYLsuGyfb1/SJ9m8eTOnnHIKZ5xxhjT2ESANfkdi3AtCDU4Lwu+Pdkpjf1hUO03UNW7D5rYxOD4PKtdDzlgpvNbHKCoqYufOnV09jR6PNPgdRXMNVG8KvN1jq+2yqQw39EMlBJst5WHHY1U6hsf1w+y2kaKJo8ljY6s1/IakExqcipsMbSK1rs7fc6hzNUFTJc7YDN8mt7XO991fuF0ikRw60uB3FDYjADus1TS6LTR7Wod3dQbDDf3IjkkGYLOlHAGMjR+A1eukX0wKAEkaAwDxGj05MSl4FS/ljkaf/zwkNKzZbWd1UxessswVKPZmiMv1vW/YIQ2+RNIOpMHvYEod9V0y7viEAjyKl1RtfODYEQkDiVFpiFFpSSZuv9eqhIo8feuklniNnuGGfmyxVpCvT2OvvfM+W43TjMO7m/EJBZ02pkTS25AGvxeRpDFgcltJ1yaQqI5tlbiRqImNeIzsmGQydImohYpMbRKVTiPljs6ROTa5rZ0yjkTSW5FhDx2B3URp9Ro2NZd1+FACGBmXy6zkEYxPKGB2ShGj4vM6JEuvBbU/WiZeo6fQkM30pKEMis1kXPwA4tQxHTYuQIm1qteriUraR0FBAXV1dYfcfunSpSxfvrwDZ9Q2ixcv5qqrrjqsaw73s+0PucLvCPYsZ0dz+cHbRcgAfToDW4TH2klZoxWz3U1RTmK7+9CqNOTr0wGYpB3M0sbNEc3pQJQ5Gmjy2BnfYSNI+gpLly4lPj6eadOmdfVUOg25wu8ANjTv7dD+B+ozmJ1SdFjG3uJwU2224/F6qW1ysK7UyLpSI3XNTpxuL9VmOxUmGx6vF4fbw7pSI7VNDprs7g78JO3D5LbKIuh9jEWLFjFu3DjGjRvHwIEDOfroo9ts99hjjzF58mQmT55MSYmvpnJtbS1//OMfmTRpEpMmTeKnn35i9+7dLFq0iKeeeopx48bx448/8p///IcpU6Ywfvx4jj32WKqrq1v1v2nTJiZPnsy4ceMYM2YM27f7ihS9+eabjBkzhrFjx3L++ecDHFJ/bc0NoL6+nuOOO47x48dz2WWXRa12hFzhR5vGPdS7mjt0iAGxGYfcVlEUmhxudtZaAKg0tR0t1HK8xuwIHCs3+nR+NCpBskFHbsqh7QHMTilih7W6Qzesbbt/IFYTCwnZkDa4w8aRtKaksYTmKP+Nx2vjGZIyZL/nL7/8ci6//HJcLhdz5szh+uuvb7NdYmIiK1eu5M033+Taa6/l888/Z8GCBVx33XXMmDGDvXv3cvzxx1NcXMzll19OfHw8N954IwCNjY388ssvCCF49dVXefTRR3niiSfC+l+0aBELFizgvPPOw+l04vF42LRpEw888AA//fQT6enpNDT49rRmzJhx0P72N7d7772XGTNmcPfdd/PFF1/w8ssvR/LjDSANfpRwOi14S1fQYKnq0HFmpxQdtI3d5aGu2YnD7YnKCt3tVahrdlDX7CBBr6HJ7mZ4dgJ6rZraJgdxMWoMuvA/pcGGLHL1qey113fIpu6KuvW+F1UwO+2GqPcv6Z4sWLCAOXPmtJIsbuHcc88NfL/uuusAn3Tx5s1BN6PZbA4IpoVSVlbGn/70JyorK3E6nQwcOLBVm6lTp/LAAw9QVlbGH/7wBwoLC/nuu+8488wzSU/3uTVbpJEPpb/9ze2HH37g448/BuDkk08mJSXlkH4+ByMqBl8I8TpwClCjKEqrwovCt4P4DHASYAXmKoqyJhpjdwu8HpavWdShQ0xKHESjq+0olSqTjRitGpfHi83lodHScdWjWm4gW6rC/2HG5SW3ahuj0lJoyKZfTHKweHlH4HGBWttx/UvCONBKvCNZvHgxe/bs4fnnn99vm9BghZbXXq+Xn3/+mdjYAz+hXn311Vx//fWcdtppLF26lIULF7Zq8+c//5kpU6bwxRdfcPzxx/Pqq6+2KY18qP0daG7dWR55MXDCAc6fCBT6vy4FXozSuF2L2wHWBl+Jvg5Er9ISp9aTq09t83yV2cGeeisVRnuHGvsDYbbvf9w4tR6NUJOgjjwstE1KvgWvjNzpzfz66688/vjjvPXWW6hU+zdb7733XuD71KlTATjuuOPCbhLr1q0DwqWRAUwmE/379wfgjTfeaLP/nTt3MmjQIK655hpOO+00fvvtN4455hjef/996ut9LswWl86h9Le/uc2cOTMgibxkyRIaG6NjY6KywlcU5QchRMEBmpwOvKn4dh5+EUIkCyFyFEXp2cIye35Ccdk7zFedooljWFw/1G3cl4srzTjc3nb33WyxUl3XgNVqx+V24/V6McTqSYiPIy0lCUOs/rD6a9kjGNkvEUUBnSZ8zjOShwG+PYX/Mxa3e95tsbRxM7O3I0sk9mKef/55GhoaApu1EydO5NVXX23VzuFwMGXKFLxeL++88w4Azz77LPPnz2fMmDG43W5mzpzJokWLOPXUUznzzDP59NNPee6551i4cCFnnXUW/fv358gjj2TXrl2t+n/vvfd466230Gq1ZGdnc/fdd5Oamsodd9zBrFmzUKvVjB8/nsWLFx9Sf/ub2z333MO5557LEUccwaxZs8jPz4/KzzFq8sh+g//5flw6nwMPK4qyzP/+f8AtiqK00j4WQlyK7ymA/Pz8CXv27InK/DqErUtYYSrB5o1+xEiGNpHhcf0CMe+hmGwudtVZDrmv8soaVqzdxOZtO9m6Yze79lZgsR648Ep6ajIDcnMYXljA+JHDGDdqGGkpSYc85pj+SZgdLpJjW4ucdUTY5uDYLGJTBpHef5J073QA3UEeWdKa7iqP3JYzqs07jaIoLwMvg08PvyMnFQ06wtgDjIzPbXXM61X4rdx00GsVRWHrjj188e0yfvhlDXvLfRvJCfEGhg8p4NTjZpKdkUpmeirxcQa0Wg0qocJqs9NksVBT10hpeRW7Siv44LNv+ddHvlqao4cP4ejpE/ndrCnk5mQdcA4t84zVOhiWnRD+2fyaOJss0UtM22GrBls1+Q3bGZQ6DPKPjFrfEklvobMMfhmQF/I+F6jopLE7jI4Iv5ycOLjNVT2A2XFg/7zVZuffXy3l4y+/Y8fuMjQaNVPGj+JPpx/H1AmjKcjrd8gbQR61HrXHjtPporhkFyvXbuL7n1bz7Gvv8uxr73LkEaP54ynHMGvqEWg1+/8zsrl8Mf0tUT0AGTpfkleuO5WyKEfw7LXXMciW6ffpCziAv1ci6Wt0lsH/DLhKCPEuMAUw9Wj/vbUBd01xVBOsdELDAH06hjakCdaVGg94bYPRzL8++pIPPv+WpmYro4cP4fZr5nHcrCNJSow/4LUAprRxJNWvA8CSMAhnTBqJjRt889JpGVs0lLFFQ7nkvDOoqK7l829+5JMvv+em+56mX3YG8845jVN/NxOdbv+ulEqTnYHp4YJtKdr4qBt8ALvXhX77f31vpF9fIgkQFR++EOIdYDaQDlQD9wBaAEVRFvnDMp/HF8ljBS5qy3+/L92yxGHpStyWGpYbt+Ft2yt12CSqYzkiMTxGV1EUrC4PO2qa8e5nGKvNzlsffckb73+O3eFgzvRJnH/myYwpClYEssXlEut3nTRkTkURatKqlwXOuzXxGDMno3Y141XHoKh8Rjulejlqv6SzgsAWn0eMrRa1x+f793i8/LhiLa+/8ykbtpSQlZHKFRecyanHzTxgFMW+4Zst/vwjkwr5xbT9EH5ah0YgX0GoIXs0JOZEre++iPThd0+6xIevKMq5BzmvAPOjMVaX4mhibfWvUVVt7B+TSqEhu9XxrdXN2F1thxoqisIX3y7j6Vfepr7RxJwZk7jqorMZmN+/VVu3NpH67KMQihev/+mhKbkIlceB3ZCDInxuFo82/EnAo4lD7bHTkDU9cJ01cQg6ex0oHhIbNzF72gRmTT2CFWs28vfFH7DwiZd5/z/fcvP8Cxhb1LZefb3FSVpccCN3YsIgLF4HepWWIbHZlNiik7hW7TSRqU1E4IHaYhAqqNsGBTOgA4XlJJLujMy0PRy8nqhL9Bb4Rcf2ZX/Gfm95FQ888xor125i9PAhPLHwOsYWDcWcOgYafgtra0kcgsN/Mwl9SHAYslGpIF6nQQiB2RbcG8hMjKHG7KAppQiN0xww9i04/fO12+vR26pwx6QwZtYpLB4/kq++X87Tr7zD3AUL+ePJc7j2kj8TH2cIu760wYrL4yU70Rf2Ga/RE4/vda4+lVJHPQ5v5LkExZZyqjUmxiT4w9mqNoDX5UvSkuURJX0UafAPA68zupu0o+Pz0arCfwX1FielDa1vKoqi8M6/v+bZV99Bq9Fw29UXceYpx6BSqbAbcgKGGEARaoTiwW4Id2OMzk0iK7F1fH1JTTNVJjtTB6ehVgkURaGmyYFKZBCjVbFyZwPZSXqqQnR4mlOKaPa7TYTXjaF5DycdM4Ojp01k0T8/4q2PvuSnVb9xz/WXcOSE0WHjVZnsVJnsDM6IJ0Ef/vmnJhVGLWyzwd1MtdNEliYyRVGJpLcgDf4hYHVZ0boc/Fz8ftT6TNbEkbaPG6W0wUq9pXWYZ0OjiXsef4llK9dx1JTx3HntxWSm+7JuFaGiOTnch1efM4tBGXH002sw29wYdGr6Je8/y3VIZjxDMoNzEUKE3RiOLfKFYLYYfK1GhSsk6UsJuWnpY2O57tLzOGbGZBY+/hJX3PoQ5/3xRBZcfC5abfif247aZhJjtQxMM4RFD81MHs4K846orfSzdEkyNl8iQRr8Q2JlxS9g3M1+UgfaRbYuPInJ7fW2aexXrt3E7Q89T1OzlVuvmsuZvz8FW2IBZpUOZ0zrMoQQNNAAmQltNmkXM4f6VDp1GhXNDjd7661U+BU1LYlDcOmScftDLvtnH8U7Qwbw9Cvv8K+PlvDbpu08fOfV9MsKV/o021ysLzORnaQPuHlUQsXUpELWmHdh9hw4QexQWGPexRHJ/j0FczmkthaxkvQMSktLmT9/PsXFxXi9Xk455RQee+wxdDrppjsUZJDywajeBJVro6KXE6vSMSouj9kpRYHC4gAby01sLDeHtVUUhbc/XsKVtz5EYkI8b71wP38460805szAHpeLMzYTVGrfF5Bs0JIzeg55Ew4kaRQZOo0qIJkQH6OhqF8i2Uk+I22Lzw8YewBFpSVGp+OW+Rfy6F0L2LW3nHMuv43lq9a32XdVG7LN4xIGkBgF/R2zx4bH45d9rt0ScX+SrkFRFM444wzOOOMMtm/fzrZt22hubuaOO+7o6qn1GKTBbwuvBxQF7CZ2Vq6KSpdTkwqZkjSEdF1wye3yeFlXasS9T9ylw+nknsde4rEX/8nMqUfwz+fuo3BgPi5d21WpZg3LYGJBKllZWaRFSUb1UBnVP4lji7KYMzyTUf3bll743cwpvP3ig2RnpnP1nY/y7r+/brPdpgpTWKEHlVC1CldtLz8apaHv6fzvf//DYDBw0UUXAaBWq3nqqad4/fXXsVplveNDQbp02mL7fyE+E5pr2GuPjjBajKq1D3lThbnVMaO5iWvveoL1m7dx2fl/4NK//CEQ126Ly2vVHkCr7vr7tkolyE7S4/J42VrVWms8r18Wrz/7AHc+8DSPvPAGu/ZWcNP8C9Co1YE2Lo9CpdlOv6TwVf34hALWNu2OeI51zqawG66kfVx77bUBVcdoMW7cOJ5++ukDttm8eTMTJkzgpJNOoqLCl6h/2mmnkZ+fT0lJCWPGjInqnHoj0uDvj+aaqHV1ZFJhq2NtZc9W1tQx/9aHKa+q5dG7FnDkaRehMvoiVmxxeShqn5/yyMFpxMdosDjcxGrVrfrpSvJSDWQl6lm9uwGvSotHE4fWacSjicOePZEn7hE89/q7vPH+51RU+z5nrD4Y+lljdpAeFxOmtpmkMbQ11GGz0VLKeFUBSXYzaPRgLoPUQVHpW9LxtOjOf/nll2HHP/vssw7Rju+NSIO/HxxeFz9HIfNzVvKIVn+Mdc2OVu2279rL/NsewWZ38MJDtzJx7AjqDNl4VTrc2viAsT9mRGagv7iY7vnr02lUTBuSTkni8eyus/rcYwBCoFaruPaSP5PXL4sHnnmd+bc9zDP330hCfFB2YXOl78knOVbLAH8Ez+yUIkxua8Qr/bVNu5m95ycwpIO1DmJTITY5oj77GgdbiXcUo0aN4qOPPgo7ZjabKS0tZfBgWebyUOh6X0A3JRrGHlpXrbE63ZQ1hkeebCguYd519wHw+lN3M3HsiEAGrEufiqLWMSw7gRmF6T1qJTMkM4Hx+cm+zNZ95v3Hk4/h4TuuZsOWEv56w9+ob2ytAmq0ubA4gglo0VrpN7vt4G0p/djtBVklfubMmYPNZuPNN98EwOPxcN111zFv3jwMhuj8bfR2pMFvA5sncsnjMfH5TEsKlxdwebxsqw5P3lq/eRtX3PoQyYnxLH5mIfkjJmKNL6A+ZxYA/VNiOWZEJnmphoDaZE8iLT6Go4amM2VQ62pdx806kmfuu5G9FVXMu+5eaupaR0JFS68olNVNO9lrLIl6v5KORQjBJ598wocffkhhYSGFhYXExcXxwAMPdPXUegzS4LegKHi2fIFl8yesMEduDFK18ehCEpLsLk+rTdp1m7Yx/7ZHSE1O5NUn7sIw8gSaUkdhTQz6lQsz43vUqr4tYjRqEvRa32o/BIc+g2mTxvLiQ7dS12Dkspv/Rl2DMaxNSyWtFpKjtMrfaathlXmHr0Slp2vKQkoOn9zcXD777DO2b9/Ozp07efbZZ2UM/mEgDX4LDjMbm0ujUmx7SGy4GFqT3d2q6Pe6jVuZf9vDpKUk8eoTd6EdeQrukGSsAWkGJhakoOkGETjRIi0+hlFTjg3IQLhifKv+caOG8eQTj1Jd28ClNz1Awz7undAN7rHxA5ieNJQB+9EgOhwsHgfNVb9B6cqI+5JIegK9x5pEgtsBe5bT6D70soEHIrTYeLPDzY7acDfOlpLdXH3Ho2SkpfDK43eRmZ6KVxMMRSxIj6MwK4FkQ+9buWgSMigomoQzJhVHbFDjZvCss3n2bzdTWV3LZTc/SINxn6ehUmMgSkOr0jAwNpMjEiKP0d9jrwNH6/BYiaQ3Ig0+wJ7lUfHbt0VJTbix31NWyfzbHiY+zsCiR24nM711otSgfQqF9DYyU5Mxp41DUWkxpY3DlDoOgIljR/DM/TdRWlHF1Xc8itUWnn27b3nHRE0sWhHZvkaty4zV44hqGG5vJVr1ryXRoT2/D2nwAdx2LJ7WoZKR8lu5Mex9dW09V9z6EArw4iO3kTh4Ii5tIvVZMwJtji3KQqXq2T77Q2FiQQozh2YwZGABrpAnosnjR/LIndewtWQ3N977FC6XO3BOUWBDuRFPSGby9ORhEc9lpXkHlP8acT+9Gb1eT319vTT63QRFUaivr0evb61+eyC6ZyB3Z2Kpo8ZpYrOlPKJuVAi8KMSqfG6YHbXNeIOCkhjNTVx528OYmyw88do7xI+ejFWowO+WmFGYjtXZtgZ+b6TFXdU/OZZYrZo1exqpz5qByutk1lS467q/svCJl7n7sUU8cOuVgWxjjxd211sYnBFU95ycONhntCOgymEk22kBXe9+umovubm5lJWVUVtb29VTkfjR6/Xk5uYe1jV92+C7bFSWLmdrhMa+X0wKuTFpqIRA4y9A3mQPrkwdTifX3f0EZRU1vPDQLQweMRpPSKHyWcMy0KpVPTLsMhqk+itgKWodHr8ExeknzKbeaOa5194lNTmRG684PxCtFPqzBdqsA3y4bLFWkG0qg4zInxh6I1qtloEDpcpoT6dvG/ydS9lqjCzBalrS0LDwS0VRwqJKFEXh3ideZt2mbTx65zVMHFtEaLT5iH6J3UILp6sZmpXAtuomEIKm5CJcumQu+pNCfYORtz/5iv45mfz5jKAS6LpSI3mphkC5xGlJQ1lu2hbRHFbWb2SAIZmsuKyDN5ZIeiB92uA7ve6DNzoAGdqEVsZ+fVn4xuKiNz9iyXfLuWren/jdrCMB8IYIqWUmRL467Q3kpxnITzOwp97C9mpfWKsQghsu/wuVNXU8seif5PXL4qgp4wPXlDZYAwZfp4r8T9lq3ENxXJo0+JJeS59eWkbit0/XJjAyPqhe2exwtzL2n3/zIy+/9TGnnzCbeeecBkBjxhQUtY5ji7I4tihLru73YUBa0IfujElDpVLxwC1XMnTQAG578HlKdpWGtXe4o7nvoUDddqot1VHsUyLpPvQ9a+Oy4TKVYan6DWMEcff76rrs3CfWfs2GLdz75MtMGlfEHdfMC/ifPdq4sIpUktZMG5JGXqoBc9pYAGJj9Tx93w3E6mO45q7HwhKziiubsDh8T2qj4vIYGz+AyYlD2j+4w0xxQ3FE85dIuit9z+DvXMpPxe+xavc3EXWTEVKMpNnhJrSGSVVNPTfe9zT9szN5/O5r0Wo1ATG00bltFwmRBDHoNAzL9unWu7S+n3NWRhpP33cjjUYz1y18EqczKIewvaYZRVFI1yWQoo3DoNYxfR8do8PC1oi3dKUvNj801Eoi6eH0HYNvLIWtS6LSVYomDr3fD99odYYlVzmcTm649ymcThdP3Xs9iQm+8EFj+kSaEwvDioNLDsyMwnQchpzA+5HDBnHfzVfw2+btPPr3N8La7utO00bi02/Yya66zb7Y/OoN7e9HIulm9BmD765aj1vxYHJHVgotW5fM2IQBgE/9ck99sD9FUXjo2X+wedtO7r/lCgbm9wd8xt6jjWPSEePb7FPSNnqtmrGjR1OXMxtj2hGAr1ziReecxkdffMcnS74Pa7+rLtxFl6Ftf3WrUoe/0pm5ot19SCTdjb4RpdNUzTLjVgSRqZ9PTSoMK1W4r/rlB59/y6df/x+XnHcGU489DbvixaONx61L5Ojhmaj7QAZttEnQa5kzIpvvtqiwxeURayll/tyz2bxtJw8/t5ihg/IZOcxX/MJkc+HyeKO2Eb60cTOzU4rA4wZ13/hXkfRu+sYK3+lzuUSaFB5q7PctYrJ24xYee+FNZkwex+UX/BFnTBrNKUXY4vM5tihLGvsIUKkE/VNisfvdO2q1ioduv4rUlCRuvPfpMKG1TRVmnG6f3z19P0XfD4fNljIoXeET2JOyApIeTt8w+FFgdkpR2PvQMoX1jSZuvv9ZcrLSefC2+ahUKlx+qeOjh2ciiZwROYl4tEE5hZSkRJ6451oajGZue/A53J5geGZLicQsXRIzk4czPqGg3ePWOM1sadgCO76D2i3t7kci6Q70boPvtIDHjbNmc1S7rTAFV/cej5c7Hn6BpmYLTyy8joT4OBoyj0RR65g2JE2u7KPI1MFp1GXPoi5nNgBFQwdxx4J5rFy7iZfeDK91arT61E9VQkWSxkBcBPILVU4jHsUL5sgkOCSSriYqBl8IcYIQYqsQokQIcWsb52cLIUxCiHX+r7ujMe5B2fUD1q2fR5xyH1qqsLjSTI05uLp//d1PWbFmIzfPn0vhwHyM6RPx+mP0DTrp940mcTEapg3NBKGiIWs6jRmTOe34Wfz+hNm89s6n/LImGFGzu96KyxMMqYxUO/9H4xZZGUvS44nY4Ash1MALwIlAEXCuEKKojaY/Kooyzv91X6TjHhSvB4/ijVhF8ajk4YG0fa9XweEOGpHV64tZ9OaHnDhnGqeedir12TNx+/3Gx4yQrpyOwKDTcGxRFl51DB7/jfXm+RcyML8fdz7897ASiaGb6mqhimiVD9Dgaj54I4mkGxONFf5koERRlJ2KojiBd4HTo9BvZHhcEWvc58akog5RtQwtwNFgNHP7Q8+T1y+b2667ElPmZBT/jeHYoqweX4e2uzOjMB2ECmP6BGwD5nDfAw9hsdq44+EX8ISs7Gubgn8DY+MHRDTmb817ZSKWpEcTDYPfHwgVOCnzH9uXqUKI9UKIJUKIkVEYd/943Fi3f8Wapl0RdTPEEKxNG+q393q93PXI3zGZm3nkzmuIScoInJMr+85Br1UzLDsBty4JrzqGrPHHcfP8C1m5dhP/ePezQLtyoy2gt6NTaZiZPDyica1b/hOMzbeboT6yJ0iJpDOJhsFvaym7b/zaGmCAoihjgeeAf++3MyEuFUKsFkKsbnexBbcds9t28HaH2p3XG+a3X/z+5yxf/RsLrlvAsMHBVWNBepxc2XcieakG8lKDmka/P2E2J86ZxotvfsCaDcGImuLKJhotwU3cwbHt1zJaaS7BXr7a92bPT1AX2f6QRNKZRMPglwF5Ie9zgbD0REVRzIqiNPtffwlohRDpbXWmKMrLiqJMVBRlYkZGRltNDorVZWWLNbIMyTHx+YDPb7+xPCTOe+tOXlz8AcfNOpKTz7oAAJcuiTF5SQzJjG+zL0nH0aK5A+DRJXLHgovpn53JbQ8+j8kc9LnvaQhmROfp0yIa8xfTdl8EmETSw4iGwV8FFAohBgohdMA5wGehDYQQ2cK/9BVCTPaPWx+FsduktrEkousnJw4h1R/zvS1EJ8dms3PHwy+QnprMHQvm4dEl0pgxGZE2iMwEqZHTVbSoj7q1CcQZYnnkzmtoaDTxwDOvhdVgDS1ME4nsAsCPa17GrfSdkpSS3kHEBl9RFDdwFfA1UAy8ryjKJiHE5UKIy/3NzgQ2CiHWA88C5ygdWA25qrL9BakztAkY1L6iGiabC7sr+E/9xEv/Ym95FffdfDmG1By8Gj0ebTzThrTvSUQSPSYNTEXxexdHFA7kirln8c0PK/ji22Vh7VqklOPVkd2gPXhZZtwaUR8SSVts2fMDpsr1HdJ3VALF/W6aL/c5tijk9fPA89EY61CweZ3tum50fD6pGl8BDrfXGybGtXT5r3z0xf+48OxTmDRuJHUZk3zXSLnjbkFSrBb8+yeWhMFceNYpLFu5joefX8wRo4fTL9t3U95e08y4vGTy9emkauP5NcKNfYkkKjT5iu644lKpqlxFFTA7Z2zUh+ndmbaHSarGt+mqKOF++7oGI/c++TLDBg/gygvPChwfkhkv5Y67ESqVr+aAotJgypnG326+AgHc+cjfw0I1nW4vQggSNLHkxqRGPrC1IfI+JH0ab/lq1m/9hN2bP+7QcaTB9zPUkBOIsNlRG1zZtxQht9nsPHDbfMwFx1OXczTgi8qRdB+OOGISI0Ydgd2Qg1djoF92BrdefRFrN27ljQ8+D7Rr0dqB8EI27WGVeQd1ZSsj6kMiMbttNLotlDeXdeg40uADGdpE+sWk4PZ6abQ4aXYEi5t/8J9vWbZyHQsu+TODB+T6DgrB1MGRRXpIoo9WqyUmexhT/XsqitBw0jHTOW7Wkby4+AOKtwfdNy0buFoRmVfT4nGwsWGTb5XvdoBJ6u1IDh+H4j54oyggDT6+1T34CmiEhu/tLq3gyZfeYtqksZxz+nGY0nwFTPLTDMTFSJ2c7kpcjIbCrHjqs6bTkDOL26+ZR2pKEnc89AI2ezCfYk+9hRihZVz8AEbF5R2gx0OgdAXs+hGqfgOXPcJPIOlL1JjLKLZ0zkKhzxv8VE08Gr98gsURjMjxeLzc/dgiYmJ0LLzhUowZk3HFpJAWr2NoVmQhfZKOZ0BaHKjUKCoN7qEncd9Nl7OrtIIXFr8faNNodVHWaCVZG4fW7/9vLz82bqHB0eh/J3XzJYdGlaWKzZvf67Tx+rTBn5k8nDEJ+QghaLKHP1K99dEXbCgu4dar5qIecRIenc/Ij89P6YqpStrBjEJfbp+i0jDkpCv502m/4+2PvwrLwm20+hQwkzQGBsW2XxbDg5ffmvdi9bQvQkzSB2muxbb9v506ZJ82+KoQYbQdtcEEqx17yvj74g+ZM30SE8+7E68mFoAJA6Sx70notcFVu6LScM1fz6VfVjoLH38Jmy3odqltcuDxKuTr01G1qRRy6HhkMpbkUClfzR57O+Vj2kmfNfihPttQRUW3x8M9jy3CEBvD7QvmBSJ39Fo1KXG6Tp+nJDKOGZEZMPyGWD0Lb7yc0opqnns9+BhdbrSx1793kxZhBq7RbcVpqYuoD0nfYIUpMkWA9tAnDX6/mBTS/S4aRVEoNwaF1t784As2bd3JLdf8FaXo94HjLe4BSc9CCMH0IcGIqoljR3Du74/nnX9/zer1xYHjJpuLumYHI+L6MSVxSLvH22GrZuXWjw7eUNJ3sRlxVhe3O0E0EvqkwW+JyrE43JSEuHJKdpWy6M0POfaoycw+/pRA5qaUPO7ZtDyluXTJAFw970/k9cti4eOLsIa4dsoabaiEilh1ZE9ybsUjk7EkbaMosPdntpX92CXD9zmDH1rqbntNcyAyx+V2c/dji4iPM3DbNfPCZI6l5HHPZ/KgVEzpR1DXbw6Wgcez8MbLqKiu45lX3glr5/RXNJueNIwUTfsT6+p2fItXunYk+1K2Cq/ipc7V1CXD9ymDn6NLJtG/ARuaXAWw+L3/ULx9ly9mOzkRt9+XO3uYFEbrDSTqtQxuka9Wqck//kr+/IcTeP8/37BizcZAu5YsXK1KTWYEWbgbLaX8sOEN1td2jAiWpGfibK7GFMVaHYdLnzH4kxIHMSyuX+B9SYjs8bade3n5rY85fvZUjj1qMgC2+HyOGJCCRt1nfkS9noHpcWQmBuvazp97NgNyc1j4xMtYrMF/whYh15yYlIjr4DbaG4NvZEJW30VRoHYby03bWN+8p8um0WesmUEV/MfdG5JN6/Z4WPj4IhLj47jlqgvDrkmVUTm9jjG5yYHXsfoYFt54GdW19Tz32ruB4+vLgrWLJyYMiqwsorHU989uroCd30vffl/DaYEd38G2r7DWFh+8fQfTZwx+qB++wRLcHf/XR19SvH03t159ESlJwUf4o4fLjdreysyhPjddY8ZkCmaewzmnH8d7n33D2o0hCVn+vxEhRFi+xmFjqaFp8yeU7vk/35ODw3zwayS9B+NeFJed5cZtrDR3fhjmvvQZg9+CK0Qmd295FS++8SFHT58YcOU0JReRP7AQtUpu1PZWdBoVY3KT8Gjjccckc/4NfyMnK537nnwFh9Nn6Pc0WNkdUg9halJhu8f7tWkXOywV1LuaD95Y0rtQFDx4cXaSONrB6BMGX6/SAmBzedhU4VthKYrC355+Fa1Gw6X3v0JD9lE0ZE7DYcgmdfDErpyupBPIDKljoE3N45r7n2d3aSWvvPVJ4LjR5gq8jvH/DUWCIjV2+iQdWNzvsOnVBj9Ll8SkxEFMSBiI0+1la1UwFOrfXy1l1brN/PXGhaRn5aCodXg1eqYMikJBDEmPIFQqY+L02Zx23EwWv/cftpbsDhwPrYObHmEWrkCAx3XwhpJeQWVzJUvLfuAn07aunkqAXm3wR8T1J06tR406rOhFbX0jT770LyaMGcGJZ50fOH5sURYJ+shXcpKeQYI+XOL6qgXXkJyUwL1PvozbE9TEqTL5InhUEeZjmD02qC+BrUugsesiNSSdw57GbeCyHLxhJ9IrDf7g2CwS1LGB96Ea9wAPP78Yp9PFXdf9FZXK9yNIjJWGvq+hUas4tigr8D4xMYlbrppL8fbdvPVhsERzldmB16tQGJuNRqgpiuvfrvH22kMSsWo2t3vekp6BqOn6qJx96ZUGP0+fxoREX0Zts92NKcQX+78fV/LdslVcdsEfSRl9XOC4VMLsu0ws8P3u7Yb+HHvUZOZMn8SiNz9kT1lloM2eBitalYYZycPI1LW/cP3Sxs0sbdyMy+uGuq6P2pB0DE3OJmzdbHUPvdTghxKqldPUbOHh5xczbPAAzj/zJJyxvvC8OcMzZVROHybZoGPWsAxc+lSMmVO4YuEzaLVa7n/qVbxeX1SXyeai2mwPbMBFuonb4LZA/Xbp0++FlDRs49eNb3f1NNqk1xv8UJ56+W0ajCaueehlTPnB1b1KGvs+j9afUe3RxpPcbxDzr72BX38r5uMvvw+0qTTZMfsL5UxOHExejKxrLAnH6rJSVvIVONsfgjsqLi/ykpv7odcafI/XGxZhsWrdJj5Z8j3nn3kyQ0eODRwP9eFK+jZzQpLtjjn3SiaPG8nTr7xNdW194PiuOguKoqAWKgYb2v+3U2wpZ6WpBFxdp6siiSLWBti6BHPxp+CN7KktXZcQkG+PNr3W4O9tCP4j2R1O7n/qVfL6ZXHOgvsDx48cLFdokiChT3pCCO66/q94PB4efPb1sFjqSnN0NHGsXifs+SkqfUm6mOZqtlur2GKtiKibaUlDozShtumVBt/j9YZt1L705keUVlRz57V/RRcf3JyNj9G0dbmkDzMtpFhK0oijuXLuWfzwy1q+Xvpz4HiNP2onGmyzVoKt8eANJd2bxt2UOyLTSZqdUoRO1bE2qVca/A3lwZj74u27+OeHX/D7E2YzefxIvP7NthH92i99K+m9GHSagIyyw5DNyX+9jaLhQ3n0729iNAcT9+r9WjuTEweTpo1v93gVjkZwWg/eUNI9cTvA23PqGPdKg9+C2+PhvidfISU5kbl3PktdvzmgUpNs0NI/OfbgHUj6JAPT4wIrfa8hnQUPLqKpycITi94KtGkpi2lQxzA6Pj+isohLK3+Ghl09ynBI8O2/7PgOtv8Xj+I9ePsDMD1pGOCr01FS00yjtWPKH/Zqg//PD75gS8lubrlqLnEpwQ25iQVSPkFyYAy64KP1oGFFnHfBBXz+zY8sXxUsaBIaFBBRWUTjHsyVa6Gu+6TgSw6B+hJMbitLGzfzo3HLwdvvh4kJg9Cq1ICvTkezw82e+o556ouKwRdCnCCE2CqEKBFC3NrGeSGEeNZ//jchxBHRGPdA7Cmr5KV/fsSc6ZMCSpggi5FLDh29Vh14ffY1CxmY14+/PfNaWB1cjze4sosklG5N0y6fe6B+h092Qcbnd39MZaxt2h1xN/Ean5Bfk73jFTUjNvhCCDXwAnAiUAScK4Qo2qfZiUCh/+tS4MVIxz0QPiXM19Bqtdxy1dzA8ewkfdg/sURyIMbmBTNqdboY7rr+r1RW1/HCP94PHN9Qbg5E8KTrEugf0/6nx6V7/+db5ddshpJv2z9xSadgcke2Ck/WxDE7JWgqd9R2vHx2NFb4k4ESRVF2KoriBN4FTt+nzenAm4qPX4BkIUROFMZuk0+WfM/q9ZuZf/XVZKYHo3JG9W9/Sryk75Gg14Y9EQ6cfR5nn/o73vn31/y2eXvgeGiFrEJDdkRjVjgau5WcrmT/RLq6H2oImsBKUzCMvK7BSEV1bUR9749oGPz+QGnI+zL/scNtExUaGxt56uW3mTi2iJPOOBsABUE/uUkraQd6rZrpQ3xG361L4uqL/0Rmegr3PfUKLlfwEdziDL6emDCo3eNts1byf8ZirJ6O2bSTRI5X8bJ019cR9TE5cQgG/76Pw+2l2uwInHvouX9w3vw7sViir8UTDYPfli7BvkuUQ2njayjEpUKI1UKI1bW1h3+XS0lJ4e47b+XOay/GGdcPU9p4MscdR5EMw5S0k1hd0A3oGHw8dyy4mB27y3j93c8Cx7dXBx/H4zV6JicOjmjMHbYq3N7uUSVJEkRRFH4o+wGqN0bUT2yIFlNxiHT79z+t4rtlq/jLH08iLi4uojHaIhoGvwwI3a3KBfZNNzuUNgAoivKyoigTFUWZmJGR0a4JjT/tUuImnwdC4IpJoX+6jMqRREaL7IKi0nLUlPGccPQ0Xn37E3bsKQu0cYds4BrUMQw39Gv3ePWuZtauX4xXbt52K7yKFyp/i6iP6UnDAjW2rSFPhk0WKw89t5jCQfmcf/ZpEY2xP6Jh8FcBhUKIgUIIHXAO8Nk+bT4DLvBH6xwJmBRFqdy3o45gdK7020siR6USDMoIrrhuuuJ84mJjue/JV/D46yRvDEn4A8iOSY5oTIvDxA/F7x+8oaTjqd0KO75nV+myiLRyYlW6QAim1elmW8iT4bOvvkN9o5G7r/srWk3HBJdEbPAVRXEDVwFfA8XA+4qibBJCXC6EuNzf7EtgJ1ACvAJcGem4h0pWSO1SiSQSBmUEM2pTU5K48Yrz+W3zdt7/zzeB46Gx+UBYFEa7aK4Cb2RJPZIo0LCTX+o3UFa5OqJuUkOyskON/ZoNW/jw8/9x7hknMGr4EMypYyIaZ39ERbhBUZQv8Rn10GOLQl4rwPxojHUoGGLUWB0exucnd9aQkj7CMSMyWWYtQuMyM+f0P7Pku5947rV3mT1tAjmZvs1dh9tDTBRXaObif5M44veg6tV5kt0TRzOoNBhdFuwRqmCqEOTpfRncnhAtJofTyX1PvkK/7AzmX3iWTxGgg+iVf0HTBqdzbFEWafExXT0VSS9DCEH/AYOxJA3FljCAO669GIAHnn4tEE5ZXNkUkF4AyIqgQhb4k7LqtvqSsayRCXRJDpPdP8LO79lhq4m4qxnJw9D7N2s3lAdDeV/917/ZU1bJnQsuxp4/C4C4DhJ27JUGXyLpSFpcO4pah27cWVx2xRX8tGo9S74LSh3XNgXD7Ea0swZuKI01m3zJWKUrwFh68AskUcGreFnftIcmT2R1C2Ylj0AlfOY2NM9i2869LH7vP5zyu6MYesp8PFrfPlFHeSekwZdI2sHUlloKQsUJF93E6BFDeOzvb9JgDG7c2l1BMbRIffnrm/cE39gaYM9ynxSDpOOwNvCDcQuN7sji4ScnDglE5SiKEkjU83i83PfkyyQkxHHD5X9B8a/+kwzaDlMEkAZfImkHoY/carWae66/lGarjcdffDNwfEtVU9hqLl0bWRWjQCq/uQLsJjCVgczK7TCUvb9EpZ+WBCuP1xuWlf3Ov79i09ad3HzlBSQnBv82xuclR2XctpAGXyJpJ6HlMQcX5HLxuaez5Lvl/LhibeD4+jITTrcvymZYBHH54Evl32wJxv1Ttw22fSVX+tHEv0+i1G5nTRSE0cYnFARel4Xs65RX1vDC4g+YMvNYjp89NewajbrjzLI0+BJJBIQW0pl3zukMGtCfB595HYs1+M/dIoqlVakjdu3UOM2U2uvDD8q6uNHDv0/yfzs+i9hvPy5+AEkaA+CL3Gq0+KJ8FEXhgWdeQyUE19z9CMbMI2lKLkKIjq+xLQ2+RBIBLYV06vrNoTlvNrffciPVdQ089/p7gTYOtzcsCzdSdtiqqXM2HbyhpF1YPdF5Ykr2b8B6vArFlcHf1+ff/sjPv27gmovPITOnPx5tHA5DNrOHZe6vq6ghDb5EEiW86hiGTDmec04/jvc/+4Z1m4IFTfbNwo2UjZZS3Ip/U3jvzwduLDlkqp0mVpp3RNzPrOQRgG81HxqC2dBo4okX32Js0VBOOu/ywPFji7JQq9qSHIsu0uBLJBFSkB6UXHDFpHDVvD+RlZHGfU++gtMZTNZp2cAdE58flXGXGbeytHGz740sjxgx7qZqii3lUemrJSqnrjlc9fSxF/+J1W7njpuuxpFYEJWxDgdp8CWSCBmSGR/0vQqBIVbPnddezK695bz69r8D7VoiNFK18cxMHo5ORDG5Zvt/ZTH09qIoWHcvY9mmtw7e9iAka+IYHXJDD03A+3HFWr76fjkXn/t7Uib9MXC8M5V8pcGXSKLEGL9QX3NiIdMnjeXkY2bwj3c/Y/uuvYE2LWXsVELFEEN0NujWNu3G5fWAsxkq1kHJ/6CpOip99wn2/szKqhVR6WpcwgDS/Ho5Znvw6a7ZYuXBZ15ncEEup1/7KPiTsIZmJXRqrQ5p8CWSKJGZqKcg3YA9LheAG674C/HxBu5/8tWAomZoGbtMXRI5uuSIxzW5rfxk2orDXA5NleBxQuX6g1/YF3E7WslTrKpZE5WuNSI8WWpnbTBh6+lX3qamvoF7rr8UTYwvcmdMXhL5aYaojH2oSIMvkUSRjAQ9CIE5ZSQpSYncfOWFbNhSwrufBiskhRa8GGzICuirRMrWipXYvS6fZrukbfb85JOnAFymcszGPViiFJWTFqKEGaqa+suaDXz0xXecf+bJjB4xJHA8M6HzlXw7RqFHIumjJMX66uAu2w4moeGEoxWWfPcTz//jfWZPnUD/nEwcbi8VJhvpcTHoNGrGJRTwi2n7wTs/CA3u5kA/ww39iKy6bi+lJUlt6xLWmXdEzdgPic0mO8bn0qsIqU9rsdq494lXKMjL4fILzgwcn1iQ0qqPzkCu8CWSKKPXqtGofVEaQghuu+YiVELwtxBFzRqzI7DS1wl1K3dApGyxtllQThJCtIz9uPgB5OpT0Qg1JTVN1ITUp33m1Xeorq1n4Y2XoY/RBY4nG3RtddXhSIMvkXQAs4YGy3PmZKZzzcXn8MuaDXz+7Y+B4y0qOCqhYkbyMCYnDkHS8Sxt3BwMZ40CLQlWAM2OYHjsyrWb+OA/33LeH05kbNFQwJegN2VQ15VclQZfIukAhBDMCsmcPOvUYxk3cihPvPgWdQ3GwPHfyoytL44Sdpu/b0cz1EeeTNSjMVeCI/rZyUUh0td7GoKbtFabnXuffJn8/tmce9PjNCcW4lX56nMk6KOzZ9MepMGXSDoKQxrW+AIsCYNRqVTcff0l2OyOMNeOVwGvv/qRRkT33/GX9a/x04qnqdv2pU9orS8XRK9ch3fXDzS6IpM6bmFMfD5ZuiQy/cVtbK6gVg746tNWVtex8MZL0SWkYo/PoyF7OjMK06MyfnuRBl8i6SiEQN9vOPa4frh0ySRNPJOr5p3N//38K//5Juja+c2feq9TaRifUECcKnqV2lyKh42WPlwwxeMGtwObx8kPxi3hdQUiIFUbH1bYZmtV8Olh9frNvPfZN5z7++MZP2p44PjMoRkdpnN/qEiDL5F0IOPzkjm6qD+m9CPwamI57w8ncsTo4Tz2whtU1QRVLxutvhT8JI2BiYmDSNHE7a/LduFWPFCz2WcAAdxO2LoEGndHdZzugMfrYZdpF16vB2X7f1m6+nlWmEui1n+o5DGEJ1jZbHbufeJl8vplcdVFZ4e102m63tx2/Qwkkl6MEAJViCiW25DJvTddhsfrZeETL+H1q2juqbdi81fIEkIwNmFAVOexzLjVVzilYYevaIrbHzpo3NvriqiUbv6IPZW/UrLl3/yfsTjq/Seqg5mxiqLsk2D1DmWVNdx9/aXExupxxvgqo80M2cTvSqTBl0g6gZH9fXop5pRR5OZkcf1lf2HFmo18+Pn/Am22VjUF/PkQVFyMFksbN0PDTl/RlBYj77RA6cqojtPpOJrDxOO8bisY92Bqroz6UGna+IAwmsPtCatgtXzVet7/zzec98cTmTjW97vzqPUMzozvFqt7kAZfIukUcpJiyU7Sg1DhiM3mjyfPYdrEMTz1ytvsLa8KtPstREpXCNHKfRApSxs3s81aGS6pbGvY/wXdHa8Hdv8IlesCh1ruZRZvdCuBzUweHiaMFqpxbzI3c8/jLzFoQH+unvenwHFDjJqCTpZPOBDS4EskncSo/r6IjqbkEQghuOeGS9Fq1Nzz2KKA1g5AvSUoqdtSMSmaVDgacXh7ScROi4yEtTF4iI5xUalCoqiqTOHVsB587nWMpiYeuHU+Tfm/w5Q2DoCxg/MDTwTdAWnwJZJOJC/VAEJgTJ9AwoBx3DJ/Lus2beOtj74ItCltsIYVP5+VPIJMXXQldH82bWeHtRpziy9/7wqIYlWujqbUXMrS0qXBn5PXhcdSS8Vvb1Pm6PgnlqqQbNol3y3nv0t/4bLz/8jwIQUgBK6YVIZPPQmS8zp8LoeDNPgSSScyLDuBGYXpuHVJWJKGMOuMC5kzfRIvLP4gTEY51DcshKAwNvrKOKWOetY07eJn03Yclmpo7jmSyrvMuwDw+Kt+mdxWftzwps9d1QHk69MCr0OF0apr63noudcZPWIIc885FUvCIBAqxuQmoY9L6pC5RII0+BJJJxMai21PGMB1d95PYnwctz/4AnZH0J0TKqWsVWmYnDgkqjH6LTi8Ln42bcdk3I3NbfMVRe+ArNSOoMXgr23a3WFjzE4pYlBsFl6vEmbsvV4v9zz+Ei63h7/dciUatRqv2qeAmZnY+UqYh4I0+BJJFxCokAXo88dx702XU7K7lKdfeTtwvMnuxu4KRp8Y1DrGRTlcM5S1FctZUbaMpm1LYPcy2PE91BR3X1eP183Plb9EVRdnX0Jdabvqw7N03/vsG1as2cj1l/2F/P7BJ7DZw7pHCGZbSIMvkXQReanBDdnpk8Zy3h9O5L1P/8uPK9YGjm+pasLhDhp9rUrD7JSijptU1W/82rSLKocR3HZfYtb2rw92VcfissOuH3xPHv733sa9viIvNVs6ZMjZKUVMTBjEMEO/wLGWamUA23ft5emX32bG5HEc/+cr8fgT5fQxWjTq7mtWI5qZECJVCPGNEGK7/3ubIs9CiN1CiA1CiHVCiNWRjCmR9BaGZSegDknKuubicxg6aAB3P7aI2vpg1ElxZVPYJm5nsMVagckdUiPX7t9TaNxNY8UaSps6Ua7BXM7mxq0s3fgvrJYa9m7+ACz+/YYoSRyH0uI2i9foUfsjc6rM9sB5m93BrQ88R0K8gYU3XoYzvj/G9AlYEgYzoWh4m312FyK9Fd0K/E9RlELgf/73++NoRVHGKYoyMcIxJZJeQ+jjv06n5cHb52O3O7jnsUWBLFwI38QFyItJo6OpdQYrc2HcS9XWz/FWb2L93u/ZYdzhE2Mr+bZVycBoU2Wro8ZpBoeJlRv+yU5bTYeONzFxUNj7eouTKlPQ4D+x6J/s3FPO3265Ev2ACQAoKg39Bo9E1Y1X9xC5wT8deMP/+g3g9xH2J5H0KYQQ5KbGYo0vwJowkMEDcrn+sr/w868bePuTr8Laris1BjJxBxuyiFV1bBGNMkcD1U4TbsVDde0mtjRuZbe91ney/Few1vuMfulhFgB3O/a/L2DcC14PNreNKosvIW2LaWcEn+LwCY2bVxSF0obgk843P6zgoy++Y+6fTuXICaOx+ROxUuJ03SrBan9EWuIwS1GUSgBFUSqFEJn7aacA/xVCKMBLiqK8HOG4EkmvYXh2ItWG8WwoM6G3lHHmKcewfPV6nnn1HSaOLfLFdvvZWdfMkMwEAKYkDWFTcym1ro6LqCm2lIe9rw5Z9e/d+Q1JGsOhJYd5vVDvL+PY4Dfg+UdCbIgX2FIH1ZvAWMoKtQvcTuqcK6Cp86p35cYEi5N4vUpY5nNFdS33P/UKo4YP5sq5Z2FKHRc4N2FA15QsPFwOusIXQnwrhNjYxtfphzHOdEVRjgBOBOYLIWYeYLxLhRCrhRCra2trD2MIiaTnkpWoJ0arwpowCCEEd19/CanJSdx8/zM0WYIrzNCKSgAj4zs3sSc0Q3enrYa1TbuxeXyhpHa3HcXr9a/SvbDjO58Br9/h2/ht2Bk09oCy52c2Vf2KueUm4lfuXF+3wfe+eiN1jdFTuTwUsvz69lanO8zYuz0ebn/wBbxehYduvwpVTAIuve/mEBpx1d0RkWwGCSG2ArP9q/scYKmiKMMOcs1CoFlRlMcP1v/EiROV1avlHq+k77Bs1RriTVuxG3IoXvk9l9zwN2ZPn8hjdy1olaI/Li8ZoEPDEg+FUXF5lDgbsWu05KaNYIjFGN5Aa0BxWthjryNVG49GqHB63ehUWlZaKyBrBAm6BAbU7KDO1USV09jWMB1GiiaOoYYcYtU+F5nV6WZbdXNYm+dff4/X3vmUB2+7ihPnTMOjMdCYeSRHDEghNa5r6tPuDyHEr/vbK43Uh/8ZcKH/9YXAp20MHieESGh5DRwHbIxwXImkVxJaEWn8qOFc89dz+N+PK3lnH39+KOMTChgTIurV2Wy0lGJ3NYOtkbKy5Sxt3Eyz2463RefGZaXe1cxuey1rmnax0ryDdc17WGkuAbcVLLU0OZvYaCntdGMPvjKFLcYeoKwxXCfnh1/W8No7n/L7E2Zz4pxpgeMj+iV2O2N/MCL14T8MvC+EuBjYC5wFIIToB7yqKMpJQBbwiX91ogHeVhRl/3+9EklfRmsgQa/BoonDozFw/pkns2bDVp56+W1GDR/CmKLCQFNFURBCBHzos1OKuny138LqpsPYaDXu9X11ASmaOLSqcDNodQbdZuWVNdz5yN8ZPqSAW66ai0cdi9rjuyH0T46lpxGRS6ejkS4dSV+krKKCLY0qhOJGb63EU76Oc6+8A6/HyzuLHiQ5MSHQdmS/RLQhoYCVjka2dpCeTG9Dr9IyIWFgwOB7vAoby00BrU2H08ncBQspq67nnRfuo39ONsaMSaTUrmTs4P6IgfvdiuxSOtKlI5FIokxWVjYxOjWKSostPh9vwUwevXMB9UYTdz7897D4/E0V5rCkrJyYlI7NxO0lTE4czJFJhWhVGhosTkw2FxtCjD3AI8+/wZaS3dz6wFPk5mThVfsSssbmJiHoPpLHh4M0+BJJN0OrVnFUYQZZfgEud0wyAyfM5sYrzuenVetZ9OZHYe33TcqSHJjJiUMwqIMidHsbrOyqC9fJ+fSrpXyy5HsuPvd0ps46JnhCpfJtnuu6f8x9W0Tqw5dIJB3E8JwEqv0p/dbEIRxz8T0Ub9vFK//6hMJB+fxu5pRAW6PVSbKhZ20gdgXjEwowhGzQuttIANu4pYQHn/0Hk8eN5IoLz6LRn+Bmi89n5sgCsBjA0PGZzh2BXOFLJN0UrVoVFuOtaPTcdvVFjCkq5O7HFrF1x57Aud311jC9l3x9Gtm65KiXSOzp7JsktrHcHPa+pq6B6+55kvS0ZB6+82rMWVNApaau3xyGDh/p0z5KyAJ1z1wrS4MvkXRzpg0JriadGSN58AGffv519zxBoylosKpM9kB5xEGxWQyP60eSxtApujs9gSOTghFOdpcnTNsewO5wct09T2K12Xn6vhtISUrEowtukGcmdE+N+8NBGnyJpJtj0GkYnuMzPPa4XGKGzubJhddR32Di5vufxeUOyvaWNlhxusPdFIMNWRyV3L1VHDuKFE0cs1OKmJU8Ar1KGzi+vSZcjkJRFBY+/hLF23fxwK3zKRwYntfQk7JpD4Q0+BJJDyA3xYAq5L915LDB3H39Jaxev5mHnv1HWKTO5kozFqc77Hq1UDE5cTAF+u5bnKMjSNXGA+GCaMWVZjz7uO5ff+dTvl76M/MvOpvZ0yaEndNpeo+Z7D2fRCLp5Rw9LFyb8ORjZ/DXP/+eT5Z8z2tvhye5b69ubqWhb1DHkBOT3NHT7Fbk6cPdWeVGG459noCWfLec5//xPifOmca8c04LO5eTrOfIQb3HJSYNvkTSQxBCMLHAp8rYmD4JY/oErpx7FicfO4MXFr/PF98uC2u/vswUJu0LEON3a6gQDI7tHW6K/RFanhCgrNFKbVN4wZRV6zZx92MvcsTo4dxzw6UIIWjMmAz4KpKN7JfUq1b4PXOrWSLpoyQbdBRmxbPdX/AJlZZ7rr+UmroGFj7xEhlpKUwePzLQvt7iJMWgI14f/Fcfn1CAXqUlRqUlQa1nXfMeehMD9OkMjA1/GrK5PNQ1O8OObd+1l+sXPkV+/2yeuvd6YnS+8EuPNp64lCxSBxXS2+g9ty6JpI8wIC2OFpd0ffYMtFoNT9xzHQNyc7hh4ZNs2xmuS1NS28zekJV+ksYQWOkna+M6bd4dyeDYLFQIRsXltTL2DreXrVXhm7TVtfVcdfujxOpjeP7BW4hPCurgJ8ZqKZwwB5HcudLTnYE0+BJJD2TOcL9REyoa0yeREB/Hcw/cjMEQy5W3PsSesnA9nQaLs80kI4CJCYPaPN6TyNOnMTNlBOkhYZQA22qaKK4Mj7VvNJm58raHsVhtPPfAzeRkptOQPQOAFIOWyQNT6a1Igy+R9ECEEAGj79ElUJczm8QhR7LokdvwKgqX3/wgFdXhBYQ2lptbbViCr1j3uPgBXSqx3F7y9ekUxfVv85yiKFj3KRjT1Gzhylsfpryyhqfvu4FhgweEnR+Q1jueePaHNPgSSQ9FpRLMaimCLlTY4vNJmPxnXnzoViw2G5ff/CC19Y1h1xRXmllXasRoDfdnJ2vjUIueZQ50QsOg2Ewy/VWqWrA43awrNbbSGLJYbcy//RFKdpfyxMLrmTg2XGRu2ODBkHHA+k09np71G5ZIJGFo1SrG5IUYPJWanEmn8vwDt1DXYOSKWx+ivrG1uNrueisVxvBCHy1hnEkaA6ma+A6dd6RohZpJiYPbPLd9n2pVADa7gwV3Pc7mrTt59M5rmD5pbNj5OcMziR04GVJ7vnvrQEiDL5H0cDIT9GGVspyxGYwpKuSZ+2+kvLKGS264n5q6hlbX1TQ5cLiDLo8kjYEcXTLDDf0Yk5Df7SQZkjQGkjU+l8vwuP5oVeqw87Y25BLAt7K/6vZHWLNhC/fceRMTTpmLIjSY0sYBkJqRjUrVM+WODxdZAEUi6SV8u7k67H16xXes2bCFq+94lLSUJBY9ejv9slpn2o7ISUCrUrVp9KweB0a3lTh1DGubdnfU1A9I/5hUCg3ZB23XlrE3mZu56o5HKN62i7/deiVzTjgFs9/QA4yJN5LZrwA0Ma2u7anIAigSSR9g9rBwY16XczT5v7uMxxa9gdHcxMXX38/e8qpW1xVXNvFbedua+gZ1DP1iUkjSGJie1DX+7QPJQVicbipMtjaNfUOjiUtu+htbd+zhsXuu5YSjp0FI4ZJh2Qlk5g/rVcb+YEiDL5H0EjRqFceMyCQ13q/3LgSo1AydOIuXHrsHu8PB3GsXsqG4pM3rmx3uNo+3oFWpmZgwCDUq+sWkRHv6TEsayrSkoUxKHMRof8TQuPgBrVw3LXi9Cturm6kxO1qd21texUXX3UtpeRXP3n8TR0+biEuXgiVxSKBNRkLfMfQtSJeORNIL2de9E9u8l7rNP3LV7Y9Q12Dkwdvmc/T0SW1em52kJz1eh0Z14PXgcuM2nIqbkXG5AMSr9aww+24mg2Oz2GGrPtDljInPp9ZpptJpZEriEGLVh17AxeXxsqnC3Oa5dZu2cd3djwPw1H03MmDWeShqXaBEIfQe9cu2kC4diaSPEbqJC2CLy2VAbg5vPHsfhQPzueHep3n74yWtBNbAp6u/b2GQthifUMBAfQYZukQydInEqnUMic0mX59Onj6NYYYcdMIn6TDc0C/s2mlJQ0nVxjPUkMORSYWHZOydbi/1FicWh3u/xv6b//uFy256gMSEeN549j7GjRyKR5cQZuxnDetbiqGhSC0diaQXoteqmVGYzrLtdb4D/hj71ORE/v70Q9zx0PM89uI/KS7Zze3XzCNW39q9sa7USGZiDAatus3yibFqHQNiw41nrj6YpZoTk0KO3/Vj9YS7XTT++Qgh0Asth8K26ibc3rY9Em6Phxf+8T6L3/sP40YO5an7biA5MTzrVqMWzN5HcbSvIQ2+RNJL0WvVAdfFt5urMaZPwKvS4tUYuPOZSXz85I289M+P2LZjL4/fcy15/Vq7OVr84yN0amI0bfvSDwWdymdqhhlyAjeBQ8HrVdhVb8Hp9u7X2Dc0mrjtwedZuW4Tfzx5DjddeUFACK2F/imxjMhJbPP6voT04UskfQCPV+H7LTXhBxUvWz59mjsefgGvV+GOBfP8kSwHZnT/JJweLzq1Cofbg0EX/XWjx6vg9LQWPduXlWt98sZGUxO3XzOP046fBYA1fgB2Qz9A4NXoe7XPfl8O5MOXBl8i6UO0FatfXlnDbQ89z4biEo6fPZXbr7mIxIRDz7Qd2S8RrfrwtwMVRcHu9hKr9T05uDxeqs0O6ppbR93si93h5PnX3+NfHy9hQG4OD99xNcOHFNCUNJwE0xZMqeNw+d1LUwenERfTd5wZctNWIpEAcPTw1j7s/jmZvP7UPcyfezb/+3ElZ116K98tW9Xmhm5bbN7PBuq+uDxeKk1BOYfqJgdbq5qwONxUGG1sqjAfkrFfvb6YP195B//6eAl/Ov043nnxQYYPKQDAEdePhswjA8Z+WHZCnzL2B0Ou8CWSPkZpgzXgKkmv+B4IqYe7bScLn3iZ7Tv3Mm3SWG6+8gIG5Oa0a5wYjYq0eB02lwcUsDm92N0e4mLUxMVo2oyfPxB1DUaefvltvvjfMnKy0rnr2r8ydeIYHPoMYuy12OJysSQNBWBSQSpJhkPbDO5tSJeORCIJw2R1sanShKPJiN5WRaylNHDO7fHw/qf/5e9vfIjT5eLMk49h3rmnk56a3CVzbbJYeevDL/nXx0twulxcePYpzDvndEgdRKy1nOakYT5/vb8qzMCMOAZndG/xt45EGnyJRNIKRVHYUG6ixuwgveK7Vudr6xv5++IP+M9/f0Cr1XDWqb/j3DOOJyczvY3eok+jycyHn/+Ptz76EnOThWOPmsxV8/4UeOKozz4KQ9NuLImDQajolxxLfpqB+D7uwpEGXyKR7JdvN1cjvC4AlJYi524rhqbd6G1V7C2v4uW3PmbJdz8BMHvaRM4+9XdMHFuEuh2btQdCURQ2FJfw8ZffseS75ThdLo6aMp4rLjyTYcOGovLPE6Cu35zA67xUA8OyE9rqss/RYQZfCHEWsBAYAUxWFKVN6yyEOAF4BlADryqK8vCh9C8NvkTS8ewbuRNAURCKh7SqHwCoqKrlg/98yydLvsfU1ExaShLHHjWZWVMnMG7kUGJj9e0a3+l0sWFLCT+tWs/XS3+moqoWvT6GU46dwTm/P57BA3IDxl3tspBSuwJj2hG4Y5IDfcweloEmyjefnkpHGvwRgBd4CbixLYMvhFAD24DfAWXAKuBcRVE2H6x/afAlko6nweKkweJkSGY8Vqeb5SX1YedVbitqt52khnUAWDxqli/7ia//7xeWrViLw+lCo1Yzcthghg8ZwKCCXAb0zyE1JZGUpER0Wg2KAh6Ph0ZTE/WNJsqratixu4xtO/ewobgEu8OJWqViyoTRHD97KkdPn0hCnCEwh9DVfAtHDk7r8+6btjiQwY/op6UoSrF/gAM1mwyUKIqy09/2XeB04KAGXyKRdDypcTpS43yZqQadhgFpBvbUWwPnvRoDXo0BU9o4vCodHm084885ivHnqhCV69i68ntW/1bMmt+28Pm3y7BYbfsbKgx9jI5BA3L5/YlHM3n8SCaOGYEmcyix1vKDXjt9SDqxuvZn/vZVOuP22B8oDXlfBkzZX2MhxKXApQD5+T2vqLJE0tNJjdOxp95KXqqB0oag4XfFBHVyWrR5VOlDmDapgamTJyAUN4qiUFPXQGlFNbUWL011lbjcbgSgUqlITkogLSWJrIw0cnMyUalU1OUcjVA8OIQKh1BhSR6G2tmEztlInDlcylmtEni8CjEa6b5pDwc1+EKIb4G2ys3coSjKp4cwRlvL//36kRRFeRl4GXwunUPoXyKRRJG0+BhmFKaj16ox212YrK79tvVo42nMmIxHE0d65fcIIcjKSCMrI43GjCmk1K44+IBCoIhwU+TRJWDTJeDQZ6DyutBr1Rw5KBW1SqAo9JmShNHmoAZfUZRjIxyjDMgLeZ8LVETYp0Qi6UD0frmD8XnJ2N1easx2YnVqNrUhm+zRth3z7tHG7XNEAArW+HwMzXsBsMYf+Cl+wpD+rRKoDuxBlhyIznDprAIKhRADgXLgHODPnTCuRCKJEI1aRbxaRbw/kSknKZb/FVfTVqyHKXUsKq+TBGNxm33VZ00HQFHrsIZUntoffU0DpzOI6KcphDgDeA7IAL4QQqxTFOV4IUQ/fOGXJymK4hZCXAV8jS8s83VFUTZFPHOJRNIlHDPCpzxZUtPE7jqfj39QRhw7a33n2zL4XpUO5SBFTuJiNOg0gvF5KdJl00FEGqXzCfBJG8crgJNC3n8JfBnJWBKJpHsxJDOBAWlx1DQ56J8cS36qgQaLk4o2HLYpBh0N/tdxMRoyEmLIS43lx2116DQqpg9JRy2NfIcjn5ckEkm70apV9E+OBXzun8xEPfETZuN2OhmdkQkiEbVKBfFZJCQl4/EqZCUGE7SOHJxGjEYljX0nIQ2+RCKJKoaUEHXN7FGgUkNCDumq1nHzMnGqc5E/bYlE0nGkDOjqGUhCkNkLEolE0keQBl8ikUj6CNLgSyQSSR9BGnyJRCLpI0iDL5FIJH0EafAlEomkjyANvkQikfQRpMGXSCSSPkK3LmIuhKgF9rTz8nSgLorT6Wx6+vyh53+Gnj5/6PmfQc7/8BmgKEpGWye6tcGPBCHE6v3VdewJ9PT5Q8//DD19/tDzP4Ocf3SRLh2JRCLpI0iDL5FIJH2E3mzwX+7qCURIT58/9PzP0NPnDz3/M8j5R5Fe68OXSCQSSTi9eYUvkUgkkhCkwZdIJJI+Qq8z+EKIE4QQW4UQJUKIW7t6PoeLEOJ1IUSNEGJjV8+lPQgh8oQQ3wshioUQm4QQC7p6ToeLEEIvhFgphFjv/wz3dvWc2oMQQi2EWCuE+Lyr59IehBC7hRAbhBDrhBCru3o+h4sQIlkI8aEQYov//2Fql8+pN/nwhRBqYBvwO6AMWAWcqyjK5i6d2GEghJgJNANvKooyqqvnc7gIIXKAHEVR1gghEoBfgd/3sN+BAOIURWkWQmiBZcACRVF+6eKpHRZCiOuBiUCioiindPV8DhchxG5goqIoPTLxSgjxBvCjoiivCiF0gEFRFGNXzqm3rfAnAyWKouxUFMUJvAuc3sVzOiwURfkBaOjqebQXRVEqFUVZ43/dBBQD/bt2VoeH4qPZ/1br/+pRKyMhRC5wMvBqV8+lLyKESARmAq8BKIri7GpjD73P4PcHSkPel9HDjE1vQghRAIwHVnTxVA4bvztkHVADfKMoSk/7DE8DNwPeLp5HJCjAf4UQvwohLu3qyRwmg4Ba4B9+t9qrQoi4rp5UbzP4oo1jPWpl1lsQQsQDHwHXKopi7ur5HC6KongURRkH5AKThRA9xr0mhDgFqFEU5deunkuETFcU5QjgRGC+393ZU9AARwAvKooyHrAAXb6n2NsMfhmQF/I+F6joorn0Wfx+74+AfymK8nFXzycS/I/hS4ETunYmh8V04DS/D/xdYI4Q4q2undLhoyhKhf97DfAJPpdtT6EMKAt5MvwQ3w2gS+ltBn8VUCiEGOjfJDkH+KyL59Sn8G94vgYUK4ryZFfPpz0IITKEEMn+17HAscCWLp3UYaAoym2KouQqilKA73/gO0VR/tLF0zoshBBx/k1//K6Q44AeE7mmKEoVUCqEGOY/dAzQ5YELmq6eQDRRFMUthLgK+BpQA68rirKpi6d1WAgh3gFmA+lCiDLgHkVRXuvaWR0W04HzgQ1+HzjA7YqifNl1UzpscoA3/FFfKuB9RVF6ZGhjDyYL+MS3fkADvK0oylddO6XD5mrgX/7F507goi6eT+8Ky5RIJBLJ/ultLh2JRCKR7Adp8CUSiaSPIA2+RCKR9BGkwZdIJJI+gjT4EolE0keQBl8ikUj6CNLgSyQSSR/h/wEp8ozIbg8xpQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "y_beta_scaled = y * beta[1]\n", "z_beta_scaled = z * beta[2]\n", "plt.plot(coord, x, alpha=0.3, label='x')\n", "plt.plot(coord, y_beta_scaled, alpha=0.3, label='y beta scaled')\n", "plt.plot(coord, z_beta_scaled, alpha=0.3, label='z beta scaled')\n", "plt.plot(coord, signal, 'k', label=r'$\\Theta$')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The datasets still have different mean values i.e. different $\\alpha$ values. These can be removed by subtracting the mean of the dataset." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABjLUlEQVR4nO2dd3xc1Zn3v2d6US+2ZcuyXOSOC7gbsKmmmhIINbQkwFLTN+9mQ9jkzWbfZLMpQOIQAoalQyBAqCFgijE2NtUFN1m2ZMmyeps+c94/7mhm7swdFVuy2vl+Pvpo7jln7j2j8txzn/M8v0dIKVEoFArF8Mc00BNQKBQKxbFBGXyFQqEYISiDr1AoFCMEZfAVCoVihKAMvkKhUIwQLAM9ga4oKCiQpaWlAz0NhUKhGDJs2bKlXkpZaNQ3qA1+aWkpmzdvHuhpKBQKxZBBCLE/XZ9y6SgUCsUIQRl8hUKhGCEog69QKBQjhEHtw1coFEOLYDBIVVUVPp9voKcy7HE4HBQXF2O1Wnv8HmXwFQpFn1FVVUVmZialpaUIIQZ6OsMWKSUNDQ1UVVUxceLEHr9PuXQUCkWf4fP5yM/PV8a+nxFCkJ+f3+snKWXwFQpFn6KM/bHhSH7OyuCPRNoOQSgw0LNQKBTHGGXwBzstByES6bvzBX1Q/Yn2JSVEwtBaAztfhXCo766jUCgGHcrgD2baauHQ59Cwp+/OKaM3D28j1G6F3W9Azada255/aN8jEe1m0EnzAe2GoIrlKBRDGmXwBzORoPY91Echbp5G2PdO/LilKnVMUwXsfh2qPtKM/M5XoW5ndD7hvpmHQtFPfPTRR8yZMwefz0dHRwezZs1i69atAz2tQYMKyxzM9HZF7W+H/R9A6Ylgc0GgA0J+cOVp/ZUbuz/H4R3ad09Dal/lRhgzGxzZvZuXYkSyq7aNNl+wT8+Z6bAydXRm2v6FCxeyevVq/v3f/x2v18vVV1/N7Nmz+3QOQxll8IcTFe9p3+u+hPbaeHtmEeRNOvLzRqK+fX+rdkOZdjZUboLsYsgae+TnVSj6gbvuuouFCxficDj4/e9/P9DTGVQogz+U2fkq5E2Gwqn69kRjD9BWo331FfV7tCcATwPUfKbdABSKJLpaifcnjY2NtLe3EwwG8fl8uN3uAZnHYET58Ic6jXsh3LePzd3SsFt/HPRC+Ttw4MNjOw+FwoAbb7yRn/3sZ1x11VX867/+60BPZ1ChVvhDBSkh6AGbO37cyd63oHjhwMwLoHyd9j3o0Z46xi+O7xsoFMeQRx55BIvFwpVXXkk4HGbZsmW89dZbnHrqqQM9tUGBMvhDhcZyqN8FJivYM8CdUNBGRnq2IXusaNwHZhuYLGB1DPRsFCOIa665hmuuuQYAs9nMxo2D6P9iEKBcOkMFb5P2PRLUXnube/S25mAHnrBxVq0vEuSgr5GQDHPAV08wEubdph3UBlqObq4dh7UN5PK3j727SaFQpEWt8IcCRpoZHYd79NZP27VqZytyZhBBUuVroNCWjcts44u2A3RE/DSHOqgLtlHu1c65o+MgVmEmz5pBMBJGIrGZjvBPZc+bkFMCo2dpx1JqoaJq5a9QHHOUwR+shINaJiwYJ0gZ4A0HqA20UOrU3D1fdlTH+vb56jjgq4+9BjCh3Ujqgm0p5/q8/QBWYSYotWSrlbkzj+xzgJap22nw677UkrsmnwYW25GfU6FQ9Bpl8Acrjft6/ZaP2/YRlGEqogY9kU5jn0iErhO7Oo09wLqm7dhNVvzR7N8JjgKK7Lk4TD0vvgBAe/TJJBxQBl+hOMYogz9M2NK6T2eg+4NOYw+w31fP/uhNxG2yszB7MgCtIS92kwV78o1g56tgdWmRPAqFYkBQBn8wUv2JJmHcQ6SUtIW9R3SpJk8AXzBMUbbziN4P0BHx827TDvKsGdRH3UOGLqBEY+9t0lb5KnxToThmqCidwUgPjL2MxuH7I0EO+pu6HR8MR2jzhZBS4gmEONjspbyug6aOIN5AhHZfiDafJqEQCkeoaOjAEwjhD/bsqSGCjBn7HlG7dXCFkioUQEZGxjG5znXXXcezzz7b4/EVFRV9ogmkVviDDU+jYbOUEonEJEz4I0E2tOzGabLhjXRfyMQfDHOwWVPcrEtjkw+3+aP9/ljboRbtdZbDgtViItvZS3+9QqEYVKgV/mCipSrtqndT617ebf6Sj1r2sqFFkzZIZ+xDkQjt/hCNHQHK6zpixv5IafWFaGjXztXiDVLV5CUU1nT12/0hguHUAi3tPZV0bjsUde+oeH3F0fPjH/+Y3/3ud7HjH/3oR4YCahdeeCEnnHACs2bN4v7779f1ffe73+X444/ntNNOo64uNQDimWeeYfbs2cydO5eTTz4ZgHA4zPe+9z2OO+445syZwz333APAT3/6UxYuXMjs2bO58cYbY0/miWzZsoUVK1ZwwgknsGrVKmpqamLtc+fOZenSpdx3331H/kNJQBhNoNcnEeJB4DzgsJQy5blDaMUXfwecA3iA66SUH3d33gULFsjNmzcf9fyGBL4WTYkyDeuatqfta/UGsVpMRCKScERS335syhdazYJgWPv7mVSYKlA1P7OUfd7DNIc8PQvrVCJsQ54dO3YwY8YM7eDwDvC19u0FHFkwakba7oqKCi6++GI+/vhjIpEIZWVlbNq0ifz8fN24xsZG8vLy8Hq9LFy4kHfeeSdWfP3RRx/lqquu4qc//SmHDx/m3nvv1b33uOOO47XXXmPcuHE0NzeTk5PDH//4R958802eeuopLBZL7Pyd3wG+9rWv8dWvfpXzzz+f6667jvPOO48LLriAFStW8MILL1BYWMhTTz3F66+/zoMPPhi7caxYsYLvf//7vPrqqyna/rqfdxQhxBYp5QKjn09fuXTWAvcCj6TpPxsoi34tBv4Y/a4Ih7SqVskKl73gWBn4ZDqNPWhuI7vVrOv/pK0i9tobDuA0dxOGeWirFq+vimArjpDS0lLy8/P55JNPqK2tZf78+SnGHuD3v/89zz//PACVlZXs3r2b/Px8TCYTl112GQBXX301F198ccp7ly9fznXXXcdXv/rVWP+bb77JzTffjMWimdROI//222/zy1/+Eo/HQ2NjI7NmzeL888+PnWvnzp1s3bqVM844A9CeFIqKimhpaaG5uZkVK1YA2s3i1VdfPeqfT58YfCnlu0KI0i6GXAA8IrXHiQ+FEDlCiCIpZR9q9g5RKj8Ef/rNzkP+Zr70VKe017R4CYWlzuj2Bikl7W0efJ4A4VCYSCSCw2nH6XbgcNoQvTS6B5t9mASU5DmRgNmk9xZubNXKNM7PLCXb4jI+SUul9gVgdcKEE8GstpmGLF2sxPuTb3zjG6xdu5ZDhw5xww03pPSvW7eON998kw0bNuByuVi5ciU+n7EL0uj/YM2aNWzcuJGXX36ZefPm8emnnyKlTBnr8/m45ZZb2Lx5M+PHj+fuu+9OuY6UklmzZrFhwwZde3Nzc6//B3vCsfpvGgdUJhxXRdtGrsGXUotUSTD21f4msi1O3GYHwUiYoAwZGntfMIw30PPC5p4OH7u2VlCxp5oD5TXUVNbT2NBKKGhctNxmtzJ6bD6jx+YzYXIRZTNLKJ0yFqut603biISKBi08dGKBC38ogsNg1d8j907QC75mcBf06DMqFJ1cdNFF3HXXXQSDQR5//PGU/paWFnJzc3G5XHz55Zd8+GFc1jsSifDss89y+eWX8/jjj3PiiSemvH/v3r0sXryYxYsX89JLL1FZWcmZZ57JmjVrWLlyZcylY4ouegoKCmhvb+fZZ5/lkksu0Z1r2rRp1NXVsWHDBpYuXUowGGTXrl3MmjWL7Oxs3n//fU488UQee+yxPvnZHCuDb3SrMlyaCiFuBG4EKCkp6c85DSxVm8Gjz37d5en6/heJSCoaepa41FjXwofvfs5nm3ZRvrOKSCSCMAmKiguZUDaW+Uunk5uXhTPDgcViRgiB3xfA0+GjuaGV2upGDuytYfP6bQBYrBZmzZ/M8UtnMG/RNDIy06zSo+yr1+bptpsZnaXXzfmgeRdLsqdgEt3EDBz6HMx2KF3eo8+sUADYbDZOOeUUcnJyMJvNKf1nnXUWa9asYc6cOUybNo0lS5bE+txuN9u2beOEE04gOzubp556KuX93//+99m9ezdSSk477TTmzp3L7Nmz2bVrF3PmzMFqtfLNb36T2267jW9+85scd9xxlJaWsnBhqoS5zWbj2Wef5Y477qClpYVQKMS3vvUtZs2axUMPPcQNN9yAy+Vi1apVffKz6ZNNW4CoS+fvaTZt/wSsk1I+ET3eCazszqUzrDZtO3/OnY9pO1P9cV1tzIK2sq/uIuImHA6zef123n19Czu3ViClZMKUscw+fgqz5k2mtGwsdnvv5AzaWjrYs6OSHV+U88mGL2msb8FitbBg2UxWnLWAspklPXr0HJ/nxGrWG/jl2VOx9kSUTW3mDhmMNhGPNZFIhOOPP55nnnmGsrKyAZ1LfzNQm7bd8SJwmxDiSbTN2pYR57+v3KiFH2aPj/upe/rWRk+Xvnq/L8C7r2/hHy99SMPhZgrH5LL6ipUsXTmHwjFHl8mame1m/pLpzF8ynSu+cTb791Sz/q1P2fD2Z3z4zudMmFzEeZetYN6iabFHWCOaPAFGZepX+utbdrEku6z3ejwKRRq2b9/Oeeedx0UXXTTsjf2R0CcGXwjxBLASKBBCVAE/AawAUso1wCtoIZl70MIyr++L6w56mg+AMEOgI65nn8bYt4b00ghSRkMsOwJpjX0oFOb9Nz/mxSfW0dLUztRZE7jyxnOYs6CsS+N7pAghKC0bR2nZOC659gw+fOdzXntuPff955MUTxjNxdecztyFUw3f2+4LEw57KcrRSzjs9tRwXEYPXXehgFZDt2iuEl5TGDJz5kzKy8sHehqDlj5z6fQHQ9ql034YDm7pdlgwEqLK3xgTIuukttVHhz+9rMHWT/bwxP2vcuhgPVNmlHDJtWdQNvPY73mEw2E2vbuVl556h9rqBo47oYzLvr6KouJCw/FjcxwpG7kAk52jGe9IDZ8DoGAq5E6Ehj1aDV+jwu2KQcFgcOmMJAarS2fk0U1FqtpACzs6DqbtT2fsW5vbefIvr7HxnS8YPS6f2390BXMXTeuzEK7EZKqeYDabWXrKXBaeOIu3XtnEi0++w09u/wNnf+VEzr9sBRar/k+sutnHqEw7GQ59+15vLVkWJy6THasp6YZQv0v76qRxrzL4CsURoAz+AFHrNy4j6A2GqUmzMbvpva387x//TsAXYPXlKznn0pOwWnv+K8xzW8lxpbpCDrf6aPeHmZDvxGzSMnabvUHMJoHFJKht9WM2gYGCQgyL1cKZFyxj6cq5PP3g6/z96Xf5dONOrr/zQkqnjNVfr83P4TY/43IcumStzkStk3KmY+4ugkehUPQaZfAHAE84QGOoPaW9yROgqSNVU8br8fPEn19h/T8/ZdK0Yq6/40LGjjd2mXSSYTfjsJrxBsPYzCZy3el93qOyHIxKODaZBHkJ4ycVan8m5XUd3XwybZP369++mAUnzuaR+17k59/7MxddfSpnXbw8ZV/hYLOPbKeV/Az93N5r/vLoKmwpFApD1DKqP/C1QMuBlOb93jraQz42RbNOEwlHIobGfv/eGn767TV88PZnnH/5Cn74Xzd0a+wnFboZleUgy2lldJajS2PfG0rynJTkOZlU6KYox4Hdkv7PZ+7Cqfzs3ls5YflM/vrIm/z+Z4/T1pp6w2jxBimv64hJM3fijwTZ7TlERPY8wUyhOBpWrlzJsdgzXLt2Lbfddluv3lNaWkp9fWrVut6iDH5/sP+DFPVHKSX7fHVsbkuNIDjU6mN/Q2oBk43vfMEv/vUBgsEQP/j59Vx45amGiSSd5GVYGZVlP/r5p8FiNmGJxtI7rWbG5TqxmdP/CbkynNz0vUu4+uZz2fFZOf9x5xrKdxnX502UZQbY0LKbg/5GDgdaKffWEoz0bzUvhWIkoAz+McJoUzUSiVDZ6MGTtEEbCUd4Zu0b3P/rZymdMo67/ucmps6akPbcnavuHKeNDPux9dIVR69dmu8i15UaTy+E4JRzFvFvv/omZouZX/7bQ2x6b6vBmbS9hOSosdpACwd8DZR7D+sHhwZGME4xuFmzZg3z5s1j3rx5TJw4kVNOOSVlTFeSxY8++ijLli1j9uzZbNq0KeW927ZtY9GiRcybN485c+awe7cmVf7II48wZ84c5s6dy9e+9jUAXnrpJRYvXsz8+fM5/fTTqa1NFUisq6vjK1/5CgsXLmThwoWsX78egIaGBs4880zmz5/PTTfdZCirfCSosMy+pmk/HE7NmN3WXkVdUC8Va+QT9/sC/OlXz/DZR7tYefZCrvjGWSmRLskYSRMnsiBzEhEkH7fpC6MX2/MYa8/DHwniNtvpCPv5vH2/TvNirD2X+kAbExwF7PZ2X4mrsSNAs8dY276tpYP7fvEku7cfYPUVK1l9+cqUG2G6jeXRtmxmuMfFG2xuLcchcwzkT+52XopjQ2KY4J6mPbQHU/eqjoYMawZTcqd0Oy4YDHLqqafygx/8QKdOCaSVLF65ciVlZWX8+c9/5t133+WWW25JkSO+/fbbWbJkCVdddRWBQIBwOEx5eTkXX3wx69evp6CgIHb+pqYmcnJyEELwwAMPsGPHDn7961+zdu1aNm/ezL333suVV17JLbfcwoknnsiBAwdYtWoVO3bs4I477qCgoIC77rqLl19+mfPOO4+6ujoKCvTaUiosc6AIB7Vs2iTlyy/aD9Bg8EdvZOzbWjv4/U8fZ9+eg1x987mccs4iw0sVZTtw2sz4gmFs5q7DMedklJBhiWe42k1WjnOPxxcJUmDLBMAVlS22mSysyJ2JlJKGYHusf6qrCIA8a0ZM9TIdeW4bbpvZsOhKZrab7/7sWh657yVefGIddYeauO72C7BY4m6qxo4gmQ5LitqmLxLUKxIGoj8/f6sy+IoU7rzzTk499dQUYw9dSxZfccUVAJx88sm0trbG9O47Wbp0KT//+c+pqqri4osvpqysjLfeeotLLrkkZow7byZVVVVcdtll1NTUEAgEmDhxYspc3nzzTbZvjy8QW1tbaWtr49133+W5554D4NxzzyU3N7dPfi7K4PcFkTDsedOwy8jYd/hTVSrra5v4zd2P0lDXzC0/vIzjlxgnr5TmuzCZNKNnlMDUyYLMSZiEwGWO+/SXZJdhxoTVZCYDR9r3CiFixj4Rp9nG/MxSPmmrYLprrKGSJ4DdamZSodswecxqtXDDnRcyakwuf3v8bTwdPm7+/qXY7HF30P4GLzaLCZfNHIsWagl5eK/5S07ONfi5hINgVvIMg42erMT7g7Vr17J///6UwiXQvWRx8hNn8vGVV17J4sWLefnll1m1ahUPPPCAoTQyaE8D3/nOd1i9ejXr1q3j7rvvThkTiUTYsGEDTqczpa8/5JGVD/9oCPm1GrTe7ouIdxIMR6ht1W9QVlfW8Z//+hdam9v57k+vMTT2OS4rxbnOmLE3oiSaqTrBUUCGxaEz9gAOkzU1qamXZFtcrMydyRh7DsuzpzLbPT7t2NFZDkZlpm4iCyE4//KVXH3zuXz+0S5++x+P4vXonwgCoQjNniDehCLqEWOBVe1mqzZ1FWhlAf/7v/+bRx991FBepNO4J0oWJ9Kpjvn++++TnZ1Ndna2rr+8vJxJkyZxxx13sHr1aj7//HNOO+00nn76aRoaGgDNZQSaDPO4cZob8uGHHzac75lnnqm7MX366aeA9oTRKYn86quv0tTUcxvTFWqFfzTsX68ZfQN2dlRTH9S7dyIRSWWjPhqn+sBhfvWjtQiT4F//6waKJ4zW9VvNgnE5jrTaOLkWN2WuophbZpJztOG4/sBqslBgy+Qk63RaQ17sJis1/iYq/Q2xMRkOC3ariVBYUtOiN+qnnLMIl9vJX377HL/694f57k+vwZ2hX+kk7zFtbNnD/MxSbMkqm7vfUKqaCu69914aGxtjm7ULFizggQceiPXn5OR0KVmcm5vLsmXLaG1t5cEHH0w5/1NPPcWjjz6K1WplzJgx3HXXXeTl5fGjH/2IFStWYDabmT9/PmvXruXuu+/m0ksvZdy4cSxZsoR9+/alnO/3v/89t956K3PmzCEUCnHyySezZs0afvKTn3DFFVdw/PHHs2LFij6TilebtkeDgcRxJ8lSx6FwhANJxr5qfy3//e8PYzKb+MH/vY4xxanFPjqzX42Y7hrLGHtO7+fdz6STeW7zhVLCLwE++2gXf/jFkxSXjua7P70GV5LRT96UdppsLM42cBfkTICCMuXeGUCUls6xpbebtsql0w+sb96pOw6EDIx9hWbszWYTP/h5qrG3mU2MzrIbGvvZ7vExt8pgJJ3ccabDQkl+qq9y7sKp/MsPL6Oyopbf3J3q3kne4PZG0oRkNu/X3DsHPjTuVyhGOMrg9wNBGfcn+4Nhqpr0xr62uoFf//hhLBYzP/jP6xkzTm/sMx0WivOcuA1i6k0Yb6j2OdZUw9xTlmSXUeYcY9hnMZlSpBQA5i2axs0/uJT9e6v57d2P4vUk7XM0e3XuncZgOxEZMU7I8jZpPv2IytJVKBJRBv9I8DSmrCI94QBtIS97PfrkiuTwxKaGVn591yNEpOR7P7uW0WNTJYHz3eldEoZRKkdLcUL4Z8FUmLDsqE85zpHHipwZzEyMnY+S7TT+fMcvmcGN37+U8l0Hue8XTxJMqLnrC0aoTtgD+Lz9AO82f8n6lp1Gp9J8+rtfP7oPoVAMM5TB7y2Vm+LVqxLY1LqHLW37dBuWh1v1xr691cP/3PUIHW1evn331wx99lrYZeqvZUl2WZcRMWnJScjQnbgCJp+m77dngjsfJiyHKadrMe0OfWQCwgR5k8DadR3bZIQQjLJlc2LOtJS+0Vl2SvKc5CRl5y5YNpPr77iAHZ+V85ffPEckYZXuD0Y43OojEunFvtOuN6B1ZBVXUyjSoaJ0eoO/DTwN3Y9D04ZpT4hB93n9/O6nj3H4UBPfvvvqFMng4lwnwXDEMOwyx+LCYbLisB3BZqQzR9vIjITBGo29L5oLQS9kF0NntIsjS/8+W4Y2ZtIp8fcVTtMKu0TCUPNpj6dgEWbyLBk4zFaq/dqNstNdlee2YbOYOJwQqrrs1Hm0tnTwzENvkJnt5sobz4nFJGs/Uz+jEgqjp4uD1jrDULdDu2nV74LSE+N1hRWKEYYy+L3BwF8clhH2eFIlBxLVH8PhMGt+9Qz79hzk1h9ezvTj9Bl3nQW+bUnqk0W2XKa4RtMj8zTuhNQKW4XTICt6Y0mMXMnS32wMKZqrFXGxJiVoZUSFlDvqoPUgOPPAmasVJemCOZlaWNlUV1FKFE+G3cJh9D77sy5aTktTO2/87QOyczM476srYn3t/rBOzvmd5h3dyykf+gIiQS1JS5VHVIxQlMHvDUFPStN7zV+mtCVHlTz5wGt8sXk3X7vlPOYvma7ry8uwYjVQnJzkHEWxPQ9TTwqBZBfHDTFoq/ZISCuYfqSYrZDRhQxz0RztCyAc6tbgd0eOy5qiwXPpdWfQ2tzO84++RV5hDstOmRvrK6/rYHSWPfak8FnbfvKsGRTb8/olQ1GhGA4og99T/O1aAe0ukFKyr15/U3jzpQ956+VNnHnhMlaetTDlPUbqljZhocSR6t83RJhhzHH6trIzevbevsKc8BmECbrRsJ+fWYoZk04qOs9tI8dppb497gozmUxcf8eFNDe08fA9L1A4OldXt7e21R8rztIU6qAp1EFLyMPsjKQbXcivYvMVCtSmbc+IhKFhd0pzpS/uz/cHwynG/rOPdvLkX15j/pLpXHqt3ghPKnQzqdCNJWmDdnn2VBZldyEGZrHDqJkwdj6UnQlTTks/9lhSOB1KlsLUVVrG65TT0w7NtrjIsDgotufp2k0mofPNA1gsZv7lh18lf1QO9/7nE9QdatT1h5NCL5Ozm2N07vO2pq8jrBg+VFZWsnr1asrKypg8eTJ33nkngYCS1FYGvztqt2khfm16P30wEmavNx6CeTgpg/RAeQ1/+tWzlEwq4pvf+QqmBLeNkZzxGFsOx2WUYDVZsIg0ejeufJh8KuRO0GSBTWbtq5OSpVokzkCQN1HbIO6kByvqQluWYfvYHAdOa/znlZHp4o4fX0kkIvn9zx7H0xGPftrf4CUQ0hv91lBqMRkiUXdRXaoLTjG8kFJy0UUXcdFFF7F792527dpFe3s7P/rRjwZ6agOOMvhGRMIgpVaqsFlfqjAsI3jCAV38dyAUIRiOhwq2tXZwz8+fwJXh4I5/vxK7I75JWGCQdFRsz2O6eyz51oyu5+XI6brfmQO23oVODiTZFldsQ3pFTjy/wGE1U5SjT/waM66AW/71q9RWN/CnXz5DOBzfQE9ObEvW/VeMLP75z3/icrm4/vrrATCbzfzmN7/hwQcfxONJ3YcbSSgfvhG739A2QdsPp3RtbNlDQMYjcLyBsE4ULBwOs+aXz9Da3MH/+a8byMnTZ8VmJSUdlToKKXV2XaM2Rm76qldDArMNwvrH6sXZZfjCAYQQzHIXs60jXgKxONepM+Yz5k7iqpvP5ZH7XuL5R9/ikgQ3WXldh+7Jab+3jmJHPuaebHor+oVvfetbMfXHvmLevHn89re/7XLM9u3bOeGEEzjnnHOortYkvFevXk1JSQl79uxhzpw5fTqnoYQy+OkwMPaAztiHI5EUBchnH36TLz/fx/V3XkhpmT7LdHxeqlzBeEdqpm2MMXPg0Ofa69xSzX8/lDBbwZYJ3kYtrn/CMu1mmoDDZI1p7xTashgXyuOgX/PT2ywmirIdup/xilUL2L+3hlf/+j6lZeNYsCwejlnfHog9Qe3z1VEdaGZpdlnqvHytYHFAa5WWUKYYVnTmZbzyyiu69hdffHHER3Apg98LknVbkqWON77zBW/87QNOPXcRJ542X9dXku9M2aBdlDWl6xVo9jjNyNuzhmbseOfGbacGTg/+2cpcY2gKtuOJCqQ5bWZK811UNMQfxa/45tlUlh/iwd89z9jxhYwdrz0htXqD5LutsX9qfyRIufcwk5yj9BfZvx5cBeCpj+YR5Bzd51QY0t1KvL+YPXs2f/3rX3Vtra2tVFZWMnnyyK6Opp53e0BzsANvkt++1RskMcO/ct8h1t7zAmUzJ3DZ18/Svd8oGifX4o5p2BvSmQHrLhiaxj4RIXqV3bowazJFtpzYsckkcCRs4lqtFv7lh1/FbrPyh188qVPX3FfvIZQQuXPAV29cADrS+aQ2eOXBFUfGqaeeitfr5ZFHHgE0N+u3v/1tbrjhBlyuobPH1R8og98DPm3fn1LLtb497ov2tHu57xdP4s508i//+lVdjdZSAzlg6MKVkzkG8qcc+1j6QYQQgmnusYxJMPpjkzZx8wqyuekHl1Jb3ciDv/ubzqgfaNA/eRklx+Fr7sspKwYRQgief/55nn32WcrKyigrK8PtdvPzn/98oKc24CiXTidSwq7XUpojBklEnkDcjy+lZO29L9JY18IPfnED2bnxSJuipEpVpY5C8qwZZFnSSA+XLNFkCkYKGaOhvTZt9zRXEVZhjgnSjctx6NRHpx83kUuuO4OnH3yd1//2AWddtDzW5w2Ecdq0G28ESTASNi7v6InuL6jErGFFcXExL7744kBPY9DRJyt8IcRZQoidQog9QogfGvSvFEK0CCE+jX7d1RfX7VP8rYbNyYW6D7f6ONQSj7l/+5VNbPlgOxdfczpTpuszPJ0JRcaXZ0+j1FmY3thPXDEyjH3xorgMhDshmzi7OGWoEILJrnjJRrvVzLgcfWLWmRcs5filM3jukTcp3xmP8Klp8emSstLKKNfv0hRQFYoRwFEbfCGEGbgPOBuYCVwhhDBSsnpPSjkv+vXTo71unxLyw/4PDLsOB+I3gkhE6hQw9++t5qm/vM6cBVM584KluveV5ut9hd0WDx9C8fNHhTsfRs3SNk0zi+LtyfIQCSzKipcztFvNuvBLIQTX334BuflZrPnVM3S0x905+5NcO1U+fZZujDQ3e4ViuNEXK/xFwB4pZbmUMgA8CVzQB+c9dqQx9inDEiJFvB4fa/7fM2TmuPn6ty7SuW4mFbp1MsdGevAjGqsDxi/U3CjFi6A4VWMoEaPN7eLc+JOSK8PJTd+/lOaGVtbe84LOn9/hj7vf9ngPUR9II72QJgxX0XsGc53s4cSR/Jz7wuCPAyoTjquibcksFUJ8JoR4VQgxqw+u23eEfClNHWGfTsa3zReKxXNIKXn4vpeoP9zMTd+7lIys+Oo8eWW/MndmeqmE/DIteza5KMlIwp2vd+2kYWXuTJYkxNTbLCZd8ZRJ04r5yjWn8/GGHbz96kex9tpWv65gytaOxD/VBJKlpRVHhMPhoKGhQRn9fkZKSUNDAw6Ho/vBCfTFpq1RvF3yb/tjYIKUsl0IcQ7wN8AgIwaEEDcCNwKUlJQYDelbOuoNmz9qjSs5RiIR6hK0ct55fTMfvbeVr1x7uk69McNu1q3snaY04ZQTlmsbhSYTFEwxHjMSmXxq1L223rA7uTh6ntuGPxjGG9R89WdcsJQdX+zjqQdeY8r08ZRM0lxGFQ0enZTy4UALhdas1CScQAfYUnWOFD2nuLiYqqoq6urqBnoqwx6Hw0FxcereV1f0hcGvAhJ3K4sB3U6nlLI14fUrQog/CCEKpJQp1lZKeT9wP8CCBQv6d5kQ9BqKaSWLb1Uk+IKrK+t48oHXmD1/ii4qJMdlJTdhxbkoawp2U5ofrzBpxl6hx2LX5Bd6QVGOM1Z/wGQy8fVvXcTdd/yRNb98hrt+cxMOp5adnCilvL3jIEU2D9PcRfqTtVRpRWMUR4zVamXixIndD1QMCH1hdT4CyoQQE4UQNuByQBcPJYQYI6LLKSHEouh1e1YrsD8pX6eVLYwSjIRoDLbrxLcOJdSlDQVD/PnXf8XusHHDty7U+e3z3LbYitGEwGW2KR2XI0EITVJi0krD7iJbaiTTxIK4Gy0zy82N37uEwzWNPP2gvoh5YmGamkBTSuUt44dVhWL4cNQWSUoZAm4DXgd2AE9LKbcJIW4WQtwcHXYJsFUI8Rnwe+ByOQidfJ+27+fz9rg6ZiAUwZMQlfO3x9/mQHkN199+Adm5cVG0ZLnj6e5uSgiqmO+uyR4HVuPw1WnuIk7Omc5xGXFXmhCC0VlxnaFps0s56+LlvPP6Fj7dqH+CS5ZSbgklqCc2KZVNxfCmT5agUspXpJRTpZSTpZQ/j7atkVKuib6+V0o5S0o5V0q5RErZs7CYY0xHOO6nD0ciOqXGnVsreO259axYdQLzFsfLFOYbyB2PsmWnv0jpSUNPBG2gcI8ybDYJE/nWDF24pjupctiFV57C+IljeOieF2hpij/FJUsp7/EkJH7JCFR+hEIxXBl5PoegF9pqoV5fwaojHHfdSCl1Mdyedi8P/M9zjCrK0+nkOG0msp29XK3bu9G8V8QpPqHLbpfZppNfSMRitXDjd7+C3xfgoaRQzcSErLawV7/K9xhv4isUw4GRZ/DL10H1x9Cg18ZJXN3XJVWvevRPL9Pc2MY3vnOxrphJUXaarFkj0m3gKrqmm6Ivie6zknynLj5/bMkoLr3uTL7YvJt1CaGayQlZn7RV6E/adgiqNmux+ZGu6/MqFEOJkWPwmyth56uGXRXeOrZ3HCQSieANhHXZtB++8zkb3/mC1VesZNLUeAjU2KQU/86QwRxLmrC+kqUwaoZxnyI92UYpHXqWZ08FwGIyYbOYdPILp567iNnzp/D0g29wqCq+ek/cwAXtbyBG9SfQUafF5td+cZQfQKEYPIwcg1+71bB5t+cQFT7tn/1we0BXbKOxroVH17zMlOnjOeeSE2PtTpsJR4JOTplzDMdnTuSknOnMzTDIHShZqrlyckv75rOMJHJKoGwVjF+cdog16enJbjWT5dTahBBcd8cF2OxW7v+fvxIKxW/mzd5ALCmrwldnnCzUWp3aplAMUUaGwW9Lr8jYWV0J0EXkSCl56J4XiIQjfOM7F2M2awZeCL0rZ0XODMY58rCZLJiFKZ7Mkz0essZpcd2qwMbRYTKBK6/LG+ZJOdN1xwUZ8Y3x3Pwsrrn1fPbvqeblZ96NtTe2B2n0xGWu32neQaXPIFo4HEptUyiGICPD4AfaDZt9kWDstT+or2b17utb2P7pXi697gwKx+TF2icWxF02Ex2F6UumuQuhaI4qodeXZKV37xjlPCS6dk5YNpMlK+bw8tPvsn9vTay91as35nu9BouDyo1aBvDgiyRWKHrFyDD4afiwJR6pk6izXl/bxFMPvc6MuZNYcdaCWHtirPcoWxYTuio+rlb1fY8jq1fD7Va9htEVN55NRpaLB3/3PKFg3ND7gmGdOydlle9vhb1vGWZlKxRDieFt8AMd2uN4/S5dszccoDEYX/V7E1b3kUiEh+55AYDrbr9Al02bGOudbe5CznjiySrWvr8oO1Pz6RtwfOZEFmROYnl2XB4hUcwuI9PFNbeupqqilr8/HXftVDf7dBXMDFf5AK0Hj3LyCsXAMrwN/r53Yc8/Upo3tu6JZdRKKalJWN2/89pmvvx8H5fdsIqCUTmx9uRs2rTJVSVLlQBXf2Iyaz79SadA6Ym6riyLkwyLA6vJHPPpm0xCJ70wb9E0lp06j5efeY+KPfEN2TZfSCelnFywHoBwMLVNoRhCDF+Db/QPm0QgFKHJE/8nrjvUyDNr/8Gs+ZM5+cx40k+iwQDNd59S0MSRDVPOUK6cY4XVAdb0N9ZEn37yPssV3ziLrBw3f/nt8wQTXDu1rfH8i/J0q3yFYggzfA1+0mrMFwnijwQ55G+OtVU1eWmOGvxIJMJDv38Bk0lw3W0X6IxE4usy55hU373FDhOWgVklVx1TTCatDvCkUyA3VaGx0BrXOyrKjm/gujKcXHvbaqoPHObFJ9bp3hOKJlrVBJqNV/kqEUsxhBmeBj8cgvK3dU0ftuxmQ8vuWI3aZBGtt1/ZxM6tFVz29bPIK4y7a0ry4yGYJ+ZMY5wjjxTsvdtMVPQhzlxttT9qetfDbPHYfIA5C6Zy4unzefW59ynflVALN8G9Z1gHd/fr8dh8Xys07D26+SsUx5DhafANKlglkyiiVVvdwLNr3+S4E8o48fT5sfbRWXYs0U3bAmtmauWqWMilktUdjBTZ9VLKBRl2rOb47+qyr59FTl4mD/72eYIB7UkvGNaHXn7QvCs1IavmM+37/vUpAQEKxWBmeBr8bmj2xiMyIpEIa+95AbPVzLW3rda5bxKjcmZnjCeFTp0X5bcfHDj0G+l51gxW5s7UqWqOz4vvx7jcDq677QJqqup56al3Yu2JsgsBGWJzW7z6WbyjI7VNoRjkDE+D79HHUXvDAd1xY3vcv//uG1vYtW0/l92witz8uGsm0ZWTFkeWFimikqsGB2lca8mqmiV58d/t7OOnsPy0ebz61/UcKI8nZHkDcf99R9hPIJKUbbvvXRSKocbwNPiH45WM6gNtbGyNK2MmFrRuamjl2bX/YPqciTpXzqgEVw7AyTkG/mFHtlakw56p6S0oBp4ufg/TXPFyhhazSVfH4Ks3rCIjy8Xae14gHNYMfU2Lj8aO+ELhgxblulEcI1oO9tve0PA0+FHCMsLWjsrYcSgcoaJB0z6XUvLoH/9OOBTh2lv1rpyMBFfOipwZmIxKFU5Y1n8TVxwh0d9hwdTUHiHITVAyzXLEf8cZmS6uuvlc9u+t4Y2/bYi1N3v0kV4HfY0oFP1CW21c8+vQ5/22NzSsDf7hQEvsdTgS4UBjfKN28/rtfLppJxdcdQqjiuKRN4nhe1NdRem1chSDj87cCLNVy3ZOYlZGXN5aCMGY7Hg29IJlMzl+yQxeeOJtDh2MyygfTNjc3+09ZHxdj7oRKI6S6o+1r33v9etlhq3Bbwt52enRfLLJFaza2zw8/qeXmTBlLGesXhJrn1ToxmmLR+KMtacWzGba2TD1rNR2xcCTPwXyyyCr2DDb2SLMrMydGTu2W/R//lfdfC4Wq4WH732RSDTe3h+K0OaL+++bggabtXUG4ZsKxZGQRuixrxi2Br8lFDfwjR36R/OnH3yd9jYv1922OiZ7bOrNT0Kt+gcnJjMUTIn/Mk3G5SdnuDXVTbPJpMuizsnL5LIbVrFr237efWNLrL2uzR/b+/msfT9tIX3FLHzN2io/5Nf8rwrFIGVYGvyWkIc90cfvQChCizdu8Ld/upf1//yUsy5eTsmk+EZeSa4+KmdhlkHkTRdFOBSDkMmnaHIXSYy2ZeM2a+4cIYROBfXE0+czY+4knnnoHzTWx12CiRIcW9r2pV6rcqP2OH7ocwh2nweiUMQI+bsf00cMS4PfnPDYnZhg5fcFePi+lxg9Np/zL1sRa5+Q79SpYi7PnorbrC9hyIRlWhEOxdDBZNbkLgyMvssUj9Jx2y2Mj4ZqCiG49tbziUQiPPrHv8eSrlq8QZ1rZ13Tdl09BQBix0o3X9EL9r51zC41LA3+Pl+dYfvfHnuL+tomrr1tNTZ7/HHfnGDsx9vzU0rmMXFFSlKPYghhYPRzrRm6Y6s5/jdQOCaPi64+lc8+2sWm9+KlMeva/HgCcaOfWE9Boeg17XWw8zVdk5SSD5p3UZsQcNKXDEuDD1q8feLqft/ug/zjpQ9ZcdYCps0ujbUn+nAXZU1hsmu0/kTTzgZbF9r3iqFBkrDdWHturPh5J6UJyXann7eEiVPH8fj9r9DWGn9iPNRy7B6/FcOcg5tJfhoMyQgBGWJHR//sBQ1bg1/f7o8JpIVCYdbe8wLZORlccm18pVeU49CFXbrMtpTzKIYvVpOFyc74DT7RrWcym7ju9gvwevw8+Wf9KixRmsMQX2ufzlMx/AlGwqxr2m4s49GHDEuDH4lI2hMKkr/23PtUVdRy9c3n4XLHffPOhBJ4czJK9CcxW1X45Qggy6LfrE+sg1s8YTTnXnISH77zOV9sibtvGtuDMd/+uqbtqQlZ1R/334QVQx9vM9Tv0TV1hLWNfn/yvlAfMywN/sHmuCunpqqOl558hwXLZzF/SVwiITG1fqqriLwkny6uAhV+Odxwpm6625NCN+1WM257fCFwzqUnMXZ8IY/c9xJeT9ydk7igMEzIUslYCiOkhAMboEG//9MePjauwmFp8DslbiORCA/f+yI2h40rbzw71i/Qp9aPsik9+xFByWJtT6bszFiTw2RlcdYU3RPe6CxHbKVvtVq49vYLaGpo5bn/fTM2pq6tm3/Qyo3K6CtSqfpIdyilpMJbFwsj9wfDNLYHdMEBfcmwNPidvPPaZnZvP8BlN6wiOzde/WhioTvmuy91FKbq3IOKyhnOmMya4Y/iNNtS/gbsCe6+KdPHc+p5i3j7lY/Ys+NArL2iIb6Za+jaqdzYxxNXDHmSlHyr/U1UJEQVHmz20ewN9ltwwLA1+I31LTz78JvMmDuJ5afNi7UXJyRYrciZQWlyucJO8lJL5imGL1kWJzOjGbidJERqcvHVp5FbkMXae1+M1cGNRIhJMIDm2un0xRqiErJGLlJCXaogWkckbtj9we7rcB8tfWLwhRBnCSF2CiH2CCF+aNAvhBC/j/Z/LoQ4vi+umw4pJY+ueZlIOMK1t54fW81bzQJbVD/FYbIqYTSFjlG2bEwJ1ctGZcY3cB1OO9fccj41lXW8/ExcC/9Ak15mIZJcHevwDu2fvbVaK7up3Dwji0CHlli16zVoTJU8rvY3xV4fbO7/BcFRG3whhBm4DzgbmAlcIYSYmTTsbKAs+nUj8MejvW5XfPT+Nj7btJMLrzqFwjHxjbrEakfTXWP7cwqKoUDpiSlyGYmCeU6bOZaBC3DcCWUsWTmHV559n6r9mpRtJAKVjZ7YmJpAM57EgjtNFdo/e+027divQjZHFM0HDKUTpJSsa4rX7UjM4u5o99Le6kl5T1/QFyv8RcAeKWW5lDIAPAlckDTmAuARqfEhkCOEKEo+UV/Q3NzME/e/QumUsZyeoITptMU/6nTXWHKsqWqKAIyZA1njjPsUwwt7piaXUbww1jTZOZqTEgreWM0mXajm5d84C6fLzsP3vEgkrLlzgmFJMPq62t/EplZ9yB0AyRWzFCOD5Ce+KEEZd9+EwhFdEMBTf3mdn9zxBzo6+r6MZl8Y/HFAZcJxVbStt2MAEELcKITYLITYXFdnLJHQFdnZ2Vz8tdO47vYLYkqYmQ4LRdnxldoYe47+TdPOhsmnwqSVkD0Oiub0+rqKIYy7ILaJK4TALEw6107iBm5mlpsrv3kO5buq+Off45uylY1JCpoKRRdsbIkvChLrdGz/rJz1//yEZafMxe1Osyg9CvrC4Bs5wpNvaz0ZozVKeb+UcoGUckFhYZoN1a4mIwQnnXkC4yeOibUVZsbVEBcnFLQG4iF6FrtWslChAE7OnUGBNR7ZNSFBdmHRybOZs2Aqzz36T+pr4z7YxPKZB3z1tId8RGR8UxeAcP8m1igGGb5UTZywjBAmqgKQsOkf8Af53/teYlRRHudfvrJfptMXBr8KGJ9wXAxUH8GYfiE5wcqZLJ9gMgjJVIxMMkbpDhMLn5tNJjIc2t+KEIKr/+VchEnwyB9eimXddpbPBCj3HmZzWzkftiS5dxr2wM5XoWl//3wGxeAh6NVqJSTgCftjlfiC4QgHEgozvfDE2xw+1Mg1t5yvE3fsS/rC4H8ElAkhJgohbMDlwItJY14ErolG6ywBWqSUNX1w7W7Jdmo/uEVZU4wrWCkUMfQPogW2TIrt8U3/xKid/MIcLrnmdLZ9spcNb38Wa69p0bt2AjKN7/7wduN2xfDhwIaUpk2te2OV+BIVAfbvreaNv23gpDOOZ8Zcg1ocfcRRG3wpZQi4DXgd2AE8LaXcJoS4WQhxc3TYK0A5sAf4M3DL0V63J4xKcOUYCqON7dfoUMVQI6ckpWmKa4yuLGKiuurKsxcyZUYJT/7lNVqbtdJ03kCEUDiSch5D6g02dxXDh6TonNaESmmRiKTTmxMOh1l7z4tkZru49Poz6U/6JA5fSvmKlHKqlHKylPLn0bY1Uso10ddSSnlrtP84KeXmvrhuOiYVuplU6CbDYel6YOborvsVIwt3AUxYroVrpkEIEdPaMZlMXHf7avzeAI//+dXYmAONXloTqqzVBVrZ5alJ9ec37FY+/eFIJAIHPowddoZgfpxQKS3R/ffG3zZwoLyGK286F3dG/+4jDttM2y6ZdrYutV6hiOHI0sI1R8/SNS9L0M5PfHIsKi7kvMtW8NF7W/l0U7yYeX17PBZ/W0cV1f4mGoL9W6BaMUjY9w5445v5nRu0RtRWN/DCE28zf8l0Tlg6I9aemP/Rl4wIg59tUQVMFL0kyb1jS6iCJoRgbEJs/tkXL2fchFE8+se/4/XEsyWDSa4dYRSsFlThnMMCT6O2Gb/zVQjF/wYiMkI44clOSkl5XUfs9SP3vYTFYuaqm87VZf4nVmDrS4a1wZ/qKmJ59lTmJmrdq0LkiiMkMUzTkRCbb7FauO72C2huauPZh+OKmj2Kzd+/vk/nqBgg2mtTmiq8dbzb/CUboqUwIxHJvvq4K+f9f3zMl1/s49LrziQ3P67Ym1h5ra8Z1gZ/rD0Xq8mCSSR8TBVrr+gpRXN1h1ah3xNK3MCdNLWY089fwrpXP2LXtnjIZXldRyzWui6YRlYh4fFfMURpqkhpOhRo1g/xxPdrmhvbePqhN5g6awInnakPHkmsvNbXDEuDP9s9nvH2fONOVcZQ0VOyxmpPhC7tb2mKazT51gyOz9SUVIUQ5Lji8dIXXXUqBaNzWXvvCwQD8X/u2lYtWqM20EJIGigiBvpHN0VxDAj5IWKscpnowguFI7QkbOQ/fv8rBAIhrr1ttc7ATyrs++zaRIalwS+wZaYWI590irZRqxKtFL3BlQfjFwFgFiaOyyjRlUXMc8cXEHaHjWtuOZ/agw289NQ7sXZ/MBLLwn2/eSfrmrbj05Wyk9C4L63hUAxSgl5NCXP3GylddYFWvJH4xn1NS9yv//GHO9jywXZWX76CMeMKYu0T+tGV08mwNPiGWB3dj1Eo0lEwVXdYmODPT2TW/MksO3Uerz23nsp98dKH+xv0q/iWUMLxoS+g7kuoT9VLVwxiGlLzKCIywubWcrZ1VMXayus6YlX4PO1eHl3zMuMnjmHVRctjYyYVujH3oyunk5Fj8BWKoyF/su5wpruY2W5NLSTDrn9qvOzrq3BnOFl7zwuEw9qqPVk4KpwUvQFo7oGGvZrsgorPH/y0VKU0+SIh2hOK4HT49ZnWzz7yJq3N7Vx322osFmNvw/GZE2Nuw75GGXyF4ggQQlBg01b5o7Icug3cjEwXV950DhV7qnnzxXgCTmc4HsAuTw3vN3+pP2lbjbbKP7wd9ryJYuhR7dcXuOncvwHYubWCd17bzBmrl1JaFhcLzsvQ6+ZkWZw6t2FfMrwNft7k7scoFD2lZEnaLiGETqhvwfJZzFs0jb899jaHa+JGILGMnYSY8Jpi6PNJWwVVCQbfl/C7DviDPHzfixSOyeXCq06JtWc5LeQ44383ibUY+oPhbfA7QzDF8P6YimOEs2vxvWynlcyonIemqHkeZouJh+97MWbYk8vYfdauVDOHJAaRVbp9GaA64Xf9wuNvU3uwgWtuXY3drhn4XJeVgox41vbK3JmY+9lWDU9LmBnVws8u1sLqJp48sPNRDB+mnK47XJpdxkRHvG5Don5Tbn4Wl1x3Bl9+vo/33/wk1n64NW4ImkMetcofiux7p8vuRPfd3i8ref2FD1hx1gJmJihh5rqPfYj48DT4Y+drIZhCaGF1KtlK0VeY9f5Wu8nKBGchJ+ZMA8BpNVOUILtw8pknMHXWBJ5+8HWaG9sAaPeHafEGYwlZ7zTvML7WvndVqOZgZO/bKU3BhBKWifVpg4EgD/3+b+TmZXHpdWfE2kdn2XXvT3HlmLoRfjxChqfBVyiOMRYRj7hwJsgumEwmrr1tNYFAiMfvfyXW3tAe4GBTXHqh0UhYLdChRe0oBp66nZqhb6rQaeXUBlpY17Sd9S1aSG19u19Xn/aFJ9ZRU1XPtbevxumKLwTc9rhBX5Q1OdWV009PfcrgKxR9xILM+ON6YtTOmHEFrL58BVs+2M7HH8ZX84naap+3HzA+aeNeiPRQX1/RfzSWa4b+sP5pbEfHQd1xqze+ut+3+yCvPb+ek848ntnz46VVE4X33CY7LrN+tQ/AuBP6aOJ6lMFXKHpL0VzImQC5+ljpDIuDMqe2fySE0Bn9VRctZ/zEMTy65mU87fGVfWLVo3VNaapg7X5dGf2Bwt8OQZ9hV32gTXfc5Iln1gaDIR787fPk5GXy1etX6cZ1Cu8VWrNYmG0QSTjtbHCnkYY5SpTBVyh6S9ZYGD0T8lKTY8Y54iURE+VuLRYz191+AW3NHbpiKf5gRLdpu7OjOkl2IUr9Ti0Zy9OY2qfoPyreg/JUnz3A1o7K2OtQOEJTR/z39tKT66iurOPaW1fjcsdX9Ik691NdY1JPWnpSH0w6PcrgKxRHisUOZatS8j0Si6VYzXGjXzplLOd+9SQ2vP0ZH2+IuwYSJXNrAs182LKbtlCStHJThZaMVbkRmitRDBwhGU55GjuQIIVdsaeaV/+6nuWnzee4E8pi7SX5zpjO/dLsMqzJG7NTTgd7Rv9NHGXwFYqjw2SCQr3Ojs1kiUVdjMtxUpKwqjvvqysomVTEI394ibaWeOieN6iPxvmivVIX+aHD2wj7P0ipmaroYwyepqSUvN+8U9cWCMXdbaFgiAd/9zxZOW4u/7relWNJ0Mqxm/TRXkBKBFh/oAy+QtEPmIWJORklmEwCS0L1IovFzNe/fRHeDh+P/OGlmDunptmnq5AVkCHWt+wyllNurQZfi6blomL4+4/KjSlNkRRVJKhKiLZ68cl1HNx/mGtuXY0roT5t4k2/v7Npu0IZfIWin8izZpATLa+ZqJdSPGE0F151Kh9v2MHGd76ItVc2elOSsJJXkzrqd8Gu19RKvy/p3CepT1XCBE0+oZNQOKJLsNqz4wCv/PV9lp82n7kL4099OU5r7KbvNtv1IZhl+qeA/kYZfIWiH5mXWQpAjtOm8+evunAZU6aP57E/vUxTQ7wSVqI/v8eourh9R+c+ScPulK51Tdt1SpiJfnuf189ffvM8+QXZXPHNs3Tvy3XHb/a5lqQCJyaTtlGbVF2tv1AGX6HoC6adrRXZcRemHTI+Lx6maTKbuOFbFxEKhVl7zwu6lb0noPfdNwU7UAwsySGY7T797+jph96grraJG751kS7BamKBKxatlWNxMcFh8Pdhz9Aiv44ByuArFH2F1QG5pSnNk5yjYq+dtvi/3Oix+Vx63Rls/XgP776+JdZ+qEXvoulWYO3AhiObr6LHJIZgAhxOyKb97KNdvPPaZlZduIxps0tj7SV5zpixz7W4mZsxAWtixb3Caf06ZyOUwVco+hJXasJMiaMg5ssvynZSnBvfwFt59kJmzJ3EUw++rpNRTvQNA+z31nV9XaW5c/R0NKQ0BSOpIZiJyXJtrR2svecFbV/m6lNj7bkuq26zfk5GiS4vA5sb8uKZ2ccKZfAVir4k8Z86gU5fPoDNEv+3M5lMXH/HBZjMJv78678SCsUNd317PHNzn68bg7/7DVUM/UiREqo2Q9WmlK71LTsThkn21XXgD0Zix//7h7/T0e7lG9+5GKs1QSk1QQlzefZUvbEHmHBiH3+InqEMvkLR14yaYdhcbI9n4U4qjG/e5RfmcO2t51O+q4qXnlwXa2/1BgklhGrW+Ju6vm6gHao/hT3/hLbaI5r6iOTABuhIvaEml6Csbw/ogjI3rPucLR9s58KrTmH8xHjWbHIx8tQEqzO0zdoBQBl8haKvyZlg2DzJOQp3glCW0xr/91t44myWnzaPl599j51bK2LtBxq9MaO/01PDB81dFDrvqNPKJIYDUPPZ0X2G4UrIn5pQ5WtJHSbDvJdUgjJR9vhwTSOPrXmZspklnHXhct24xGLki7KmkIK5f6SPe4Iy+ApFXyMEFM1LaTYJEwuz4jIMRTn6leAV3zyHwtG5PPA/z9GRILCWGP4XkCE+a0uziducRnFTEWf/+nhCVVttWm2i5PyHxGLkoWCI+//7WUwmwTe+8xVMCb760oTV/dLsMlzmY1/kpCuOyuALIfKEEP8QQuyOfjesASeEqBBCfCGE+FQIsflorqlQDAmyimDcgm6HlebHQzWdLjs3fu8SWpra+N+ELFzQb+I2hTqo8nUjomaUoauIJ6ntfBWqPzbMpk0Og/UEQrpi5M8/9hb7dh/k2ttWUzAqJ9ae47Jiiq7uBWnkEwaYo13h/xD4p5SyDPhn9Dgdp0gp50kpu/8vUChGCCaT0IVqTiwbx4VXncpH729j/T8/1Y2tbo5n4u7xHsIXCRpLLyiOmL2e2pQw2MQw2a2f7OG159azYtUJLFg+K9Yu0DJqO1mRO7Pf53okHK3BvwB4OPr6YeDCozyfQjHsOSFTL6tclO3USS+cddFyph83kcfvf4VDB+tj7b5gRBf//WHL7q6lFxS9wh8JUumPh2ZKKXVPVi1N7fzlN88zdnwhl31Dn007sdCNyWQcoaVj2tl9Nt8j4WgN/mgpZQ1A9PuoNOMk8IYQYosQ4sauTiiEuFEIsVkIsbmurptQNIViMJMmRDPT4mRl7sxYbD5o0gudmMwmvv7ti7BYzaz55TME/HGd9Q5/L1b0nYU7/O2qVGJrDfjb0na3hDxsaInLKQTDEZ3MRSQS4cHfPY/X4+Om71+K3R7/fSUWupngKEi5oTNxhRa5ZTGobHWM6dbgCyHeFEJsNfi6oBfXWS6lPB44G7hVCHFyuoFSyvullAuklAsKC9OnqSsUgx5XPuRPgYKpht2lSWn2iaGaeQXZfOPbF1O57xCP//kV3bjKRg8VDfGV57qm7bSHDKoylb+t+aor3tOE1sIGhVVGCjWfQsX7hl1hGdGJogVCESob9fpEb7ywga0f7+GyG1ZRXDo61j4qy66LsZ/oHEWmRb8Zj82lZWBPPpWBpluDL6U8XUo52+DrBaBWCFEEEP1+OM05qqPfDwPPA4v67iMoFIMUIaCgDHJKwJmnrfQSyLG6GW/XZ+YmGv05C6ZyziUn8d4bH/PBW5/G2oNhmVLxcHNbOZ5wAEUS4VC3aqLJ4ZeJpQoBdm6t4K8Pv8nxS2ew8uyFsfYcl5WMhGLkLtPgisgx4mhdOi8C10ZfXwu8kDxACOEWQmR2vgbOBLYe5XUViqGD2Qoli7WVXhITnYUU2fTBbfaETNwLrzqFabNL+d8//J2q/V0nU33cto+2kNe4RCLA4e2aAQQIBbTVf1NFrz7KkEJK2PMP2PuWYbc/EqQu0JrSnug2a2lq40+/eobCMbnccOeFutV8nltv4HOt/Vutqi84WoP/X8AZQojdwBnRY4QQY4UQnc+ho4H3hRCfAZuAl6WUrx3ldRWKoUmGfpvLJExMcBbohzjiq0az2cyN37sEh8vOmv/3ND5vfLVaXtdBJBIP3QzJMFva9vFhy272eA6lXru1Ghr3aoaws4Ri84HhV0SlfJ1WHKYyVSohkU0te9nWURU7llJS2xp3jYXDYdb88hm8HX5u+eFlKSqYicx2j2ey02AL051uW3NgOCqDL6VskFKeJqUsi35vjLZXSynPib4ul1LOjX7NklL+vC8mrlAMSYrmpzTZhYVCayZuk7ap57aZdf05eZnc9L1LOFTdwCP36ePzKxqM9XOq/Gni9BvLtaIpnecIdHRrGAc9/na9eFzQC4e+0EpBdkEYzS/W4Q/hDYbZ3+DRre6fe+Sf7Nq2n2tuPV/nty/Nd+lW+ouyplBgy8QkDMyp1ZHaNoCoTFuF4lhiMqVonwshmJUxnlkZ4wGwmE1MKnTrpBemz5nIhVeewsZ3v+DNFz/UvT9ZWbNHJEoqd2MYBzWRsLYpXfNpj9+yrmm7TgGzttVPTbOPhIclPt6wg9eeX8/Ksxey9BR9cZLk8MvBlk3bFcrgKxTHmjFzejSsMEsfxnfOJSdx/JIZPPXQ62z7RB9m6Q+mhmuOiMIpnQJnnm6E5dKQ6BLrpPrAYf7y2+eZWDaOy5Pi7ZOF0eZklKSetOxMKI7Gpdizjmhe/YUy+ArFsUYIKFkCufp4bWuSS8BiMumidkwmLT5/7PhC1vzqGWqr40lCB5t9KVWYPmvfz7qm7VT6UnXeUziwkZTQn6FEJFqLdueraYd4wwEOJklSJLvE2ls93PN/n8Bmt3LLDy/TSR7nZ9hiwmhLssuY6ioiz3CjVoA7XytdmDP+yD9TP6AMvkIxEDhzYdR0KJwea7KaLCzJLmOMLUc3dGxO3A/scNq5/UdXIATc+/Mn8Hrim7iH2/x4A6kr/b3eHkglexuhfYhLKhvo4nTSHvKxsXUPu73aZnZyAXKAUCjMH3/5NI31Ldz6b5eTV5it68+OSidYhBmHycpYu4F0WMHUuPSxffBF7SiDr1AMJHkTdUbfYbIy3T2Wyc74JqHDaibHFZdeKByTx80/+CqHDjbwwG+eI5KwMq9p8Rm6KQKRUEpbCp6ojEPQ22VW6mAnGAmzvaOKYCSMPxJkXdN2NreV68Ykl5EEeOqB1/jy831cc+tqpkzXr8wTn7TmZRjLXwNgdabvGwQog69QDDR5E1OaxjvyybEkZN4mxXzPnDuJy76+ik83fslz//tPXV9Fg4dgWO+e+aClCx39TlqqtCSl8nVaVuret+HwjsHr6kmTOVzpr+dwoJVdnhqdXEInta0+Akk/n3WvfcRbr2xi1UXLWH7avC4vm2FJirwxis4ZpAydmSoUI4x5mfqVZHLs92nnLWblWQt49a/vs+7Vj3R9ydIAoEWnBCMhIrILA56YpBTyaYlZu1/v9dz7lKAP9r2rPXl0Hh/+Eva8Cfs/SPs2Tzh1Fd/mC6XoEX2xZTePrXmF404o45JrztD1OawmSvLiq/ZkOQwARs0EW9R9M8iN/+CenUIxwlmYFS90LYRgdELkjhCCK286hzkLpvLon17ms4/0ypnldR20evWr4PUtu9jaXkWv6awK1VQB9Xt6//6jofWgli9Qu01zNZW/DU37tL5gah5CZ4pBR0Rv8Fu9Qera9G3791bzx//3NOMmjOKm71+qK2aSn2FjbI4zVox8Ze5MSp0GBj9nPJQs1fz3GaNT+wcRyuArFIMYt9lBniW++ee2W8hy6jNxb/r+JZRMKmLNL5+hYvdB3fvr2wMpIZuNoXZaQh5dAle3NB+AivWai6ch6iYJB7VVdpqqUX1OR11MAE1KmVJzthNJ6ucKhiO6ovAAdYea+O1/PIY708m3fnI1Tlf8Zmo1CTLtPShF2Ln/YrZA/uS0CqmDBWXwFYrBQP4UyC8z7JqTWcIER1x+oSBDH5/vcNq588dXkpXt5nc/e5y6Q3oDfLDZR0vSSv+Ttgo+a99PfaCHm7MtVeBP0J3Z+Sp4GjSj30V0jCEhf/p9geYD+qxZMJR+2OOt5b3mL3XuKX8kyI6OgylZxpGITHFxtbd6+O1/PEooFObbd3+NnLxMXf/4fFfP9O0N9l8GM8rgKxSDgYIyKJiiCa0ZUJykqplYGhEgOzeTb919NeFwmP/+8SM0NehFwRraAykr+uaQh60dlaxr2k5bKNXn3y3Vn/R8bCQCdTu1r71vafsC3qRkqY56zW1zIJpJHPTBoa3xJ4oEavzaeyPR1Xx9oI0NLbupDegLkh9u9aXE2ns9fn73s8eoP9zMHf9+BWPHx900JlPqXgnADPe41M9UvDC1bZCjDL5CMZhIo51vNZk5MWdazL1jMokUo19UXMi37/4a7a0efv3jh2lr0ceZ76v3EEqzst7Sto9D/uYjM/yJSBldpUc0w167TSu+svt1TcenMSE88sCH8Y1YiCt3dj5J7HsHWiq7vJwnHCAiI2ztMB7XnrRB6/cHuOf/Pk7F7mpu+t4llM3Ub4yX5rt1Ojmg+e5H2/Qx+Vhd4NaL3g0FlMFXKAYj2cUpTRZhZlZGvN1kEroIEtBq4t7x4ytpONzCr+96BE+73oAfaPBS0dBhmKD1paeaLW37eu7mSaR8nbbibzukGfndr2uum+YDmjsoHQc/jr/uSKhwt/PVuGxCF3zcto93k/TsA6EI1c3elMSqYDDEH/7zKXZt2883vn0xxy+doesvykkVOiu25xlfeJD76tOhDL5CMYQwCxMrEwpkW8wmXT1cgGmzS7nl/1xGdWUdv/3pY3g9+mpYkYiWoJXs1+9ka0clezyHerepG/Rqxt5IxMwgkiaGv1W7KXRDREYIJiSPRQw2ZjupavLiC+pvFqGQJnW89ZM9XHvbahavOE7XX5rvxGmNq5SuzJ3J8uypugS44YAy+ArFYMIaddPYMsDm7npslBynjfFJK/3jTijjpu9dQsXug/z6rkfoaE911TS0B9K6eKr8jbSEjA11jb+J7R1HENqZjtptXWrgAGxq3cv6ll3UB9r4ot34BmEklwDayv6P/+9pPt34JVfeeA4nnXG8rt9lN2MyxU1hp7SF1WRJce/Efj9DFGXwFYrBhLsgKqxWqsV2J8guJJK8iWg1p/4rn7BsJrf88DIqyw/xqx+tTfHpg+biiUQihqv5Awmia1JKKrx1hGSYnZ4aDhtUiuovKrx1sSpeWzsqaQi2G447YJBs1umz7zT2p523WNef47QyJkvvypnu1stXx/ZVhAnGdd4slEtHoVD0Bc5czUdstmphf+MXpwwZbctmXJJ/uSQ/Vcdl3uLp3P7vV3LoYD2//LeHaG5M9c9XNHjZV5+6mm8MtceMbZW/kQpfHeXeeNnqXrl8eklLyMPuqFupwleXdlwwHKHNFzJc2Xs9Pn7zk0fZ/lk5199xYYqxz3NbycvQS1YYmvHMIu27xT7kq4Mpg69QDHZceZA3KaW5zDVG52O2mEzRakz6cbOPn8K3fnI1DXUt/OIHD1BdaWxAjTT1K3x1fNiyO6a4We2Ph1Lu8tSkhEH2FZ+0VXDQ38iHremzeqXU4uuTs2cBmhvb+NWP1lK+s5Ibv3sJJ56eWmkswyCxKvkmmoIp6uc3qE88FFAGX6EYChROg0mnpDSPd+jj800mwcSCVN//9OMm8oOfX0cgEOIXP3iAnVsrUsYcbPaxv6GDDn9Ip8CZjppAMzs6DuIJaxmsLSFP2uxXIyIywvvNO1MKiSdWo/KnKcgeiUjDpxKAgwcO8/Pv/5lDBxu47UdXsOik2SljSvNdMcmETtxmu7HBt0QT3XInavsqY4+HMXNTxw0BRH8+lh0tCxYskJs3bx7oaSgUgweDzc1EA9lJIBShqinVp11f28Rvf/oYdTWNXH/nhSxZkb76VqIkcFeUOgpjbpdCaxZTXUVYTfq6vBEZYaenhlyLG4fJSlvYR47FxZa2fbExJY4CqnwNaSNwQpEIgVAEq9lkKA4HsP2zcv7wiyexOWzc+eMrmTB5bMqY8XlO3Z5HsT2PCY7ClDkDMGEZOLJT2wcxQogtUsoFhn3K4CsUQ4im/ZqkQUKxks7ko0/b9xOScbdMMBwxNIwd7V7u+88n2bm1gjMvXMYl156O2Wxg7NAExGxmE06bcX86ZrmLybG4sZrMBCNhmkLtbO842P0bu6Cy0UMwbGyvpJT848UNPPPQPygaX8Cdd11FfmFOyjijm9jJOdONC5ADTDv7aKY8ICiDr1AMN9KEMRqt9iMRmSIvEAqGeOovr/PWK5uYNruUm75/Kdm56Ss0Tch3xsr7HUuklDR2BAnLCO2+1D0G0KQSHr73BT56fxvHL5nB9XdeiMutj7xxWE0UZTtSwiyLbDlMS47KSUQZ/GOHMvgKRRrSGPyWkId93sM0J8XQhyIRDjSkrvY/ePsz/ve+l3BlOPj6ty5i5rzJ3V66NN9FBIkJQUTKFF94X6CpYUrDOSdyoLyG+3/9Vw4drOcrXzudsy5enho7j/HKfpRNcz9ZRBdPL8rgHzuUwVco0tBNopLxSj9ChYEBrdx3iD/96hlqquo5/fwlfOWa07HZjUXcEhGARMtSNR3h6j8UiWCJvldKSasvREOSjLER4XCY1/66nheeXEdGlotvfvtiZsxNjWQCrSaww6o36suzp2I19UD+eJgZ/B58YoVCMeiYsEzTmulUluwBJpOJCflO9icZ/fETx3DXb27m2Uf+wZsvfci2T/dy7a3npwiLJdO5VDzY5GN8fvdhiuFIhHZ/OFYMvMUbpKE9QHGuE38oYhheaUTV/loeue8l9n5ZycITZ3H1zeeRkWV8/VyXNcXYL+vO2JeeGNPdH24og69QDEU6I0dMVkgTumiE2WQi22lN0dGx2a1c+c1zmLdwGg/d8wL/9cMHWXbKXC657swuffsAwYikvK4Du8WEPxQPy8xyWnBazYSiRdU9gRDeQIRwWGK3mmIreaNoIiO8Hj8vPvE2b760EZfbwTe/+xUWn3ycoQtnTLYdu8WEKalvaXYZtu5W9vZMcBVAVlGP5jWUUC4dhWIoE4kY1pytDbQQioTJt2XyoUEhb08gxKEW4xW13xfg5Wfe5bXnP8Bms7DqouWcfv4SXUWoY0kwGOK9f3zM359+l9amdk4+83gu/trpaVf1GQ4zozJTlS/dZjsLs9LsUZhtEM0nGIpunESUD1+hGM74Wros5l0baGGHQUhkJCKpbvERCBknSx2qqufptW/w2aadZGS6OPsryzl51YKUCJj+IuAPsuHtz3j52fdoONxM2cwJfPX6M5k0LVU6GrSyhHkZNtxpShOemDMt/QbttLPj+yLK4A8MyuArFD0kEoHm/VD3pXG3jKToxndipEOj699Vxd8ee4ttn+zF7rCx9JS5nHLOQoon9I90cN2hRt59fQvvvLGFjjYvE8vGceHVpzJr3mRD9w1ApsNCYWbqE0iuxc3czK73IgBl8Ht44kuBu4EZwCIppaF1FkKcBfwOMAMPSCn/qyfnVwZfoeglkTDsfsOwKyTDHA60sstTo2vvzuB3sn9vNf/8+0Y2vruVUDDEuAmjWLh8FvMWT2fchFFHHKkjpaT2YANfbNnNpve2Ur6rCmESzF88ndPPX8LUWRPSGvpOSvKdsWifRKa7xjLGntP9JKadDQe3aMJ1BrpFQ4n+NPgzgAjwJ+B7RgZfCGEGdgFnAFXAR8AVUsrUuLEklMFXKI6Atlqo/tiwKxgJs75lp364L4QnEGJ0lgN/MMzBZp/he2PjWzr46P2tbHpvG3t2HEBKSUami6mzJlAyaQxFJaMYMzafzGw3GVnOWBZvJBKho91La3MH9bVNHNx/mMp9h9i1bX9MxbNkUhGLTprNopNnG2bKJjI6y47LZu7yZpBYLKZLhviqPpF+C8uUUu6IXqCrYYuAPVLK8ujYJ4ELgG4NvkKhOAIy07taOvVinCYbk52j2dpRSabDQqZDMwV2q5kMhzltVitAZrabU89dzKnnLqapoZUdn+/jy8/3sXvbfj7+cEfKeLPZhJTaSj55gZlXmE3ZrAnMmDORGXMmMaqoG7XKKEax9Ym4TXZG2bJSO1wFWvnI2m0wdj5UbQJnz645HDgWYZnjgMQKw1VAqsB3FCHEjcCNACUlJf07M4ViuJLok05iUdYUrMKM1WRmmWUqH7Ts0vW7bRbafWGcNhPeQNfql7n5WSw7ZS7LTtHUI/2+AIcO1nO4ppG2lg7aWjyEQmGE0PIAMrJcZOVkkJufxbgJo454A9jWRXbvDPe41KLjiWQVxUMuR82EzDFHNIehSLcGXwjxJmD0E/mRlPKFHlzDaPmf1o8kpbwfuB80l04Pzq9QKIyYepaWnOVrgcqNsWaXOV70wygm3W23MC5HYLeauxQsM8LusDFh8lhDlcqjZXyeE7MQmqyDSW9WZrmLORRoZqJjFBmWLm4iyd6I3B5s6A4jujX4UsrTj/IaVcD4hONioPooz6lQKLpDCBBmrYBKFwlai7OmsDGp0Ig96i4pynEQDEvafSGEELSmKXzeX4zOsqeEWZoM1pCFtiwKjVw4ibjy05aMHCkcC5fOR0CZEGIicBC4HLjyGFxXoVB0kj8J6nYadjnNNhZnTSEgQ3zSVqHrs5hMWEzgjN4A8t3WtIVH+pqibEe3sszTXWNxmm1djmHCMjDbwXps8gcGM0clcyeEuEgIUQUsBV4WQrwebR8rhHgFQEoZAm4DXgd2AE9LKbcd3bQVCkWvyCntsrtboxlFCMGkQjel+S6ctrj5cFr7RjHTZNKkmCcVurs19llmJ2PsOWRbutHxcWQrYx/laKN0ngeeN2ivBs5JOH4FeOVorqVQKI6CxBh1Rza4R0FDquQCQJ4lA5fZRpW/sYvTCYqynYQiEdp9IXJcNoLhCJ5A2FDt0moWur0At91MjtMaCwHtrd6+02RjXk8SqhQ6lHiaQjFSKFmiacbY3BAKpBj8bIuLiY5CxtpzsZosZJgdfOnpervNYjKR49KeDqxmE9lOEw6LibCU2Cym2BiAdl+IiJRkOePSy0U5Dqxm0a2xz7G4yDQ7sZrMNAc9THKOSl+lSpEWZfAVipGCMzf+2mKDslUpwmsTnIWx12PsOeRa3bSH/bhMtpSN3XTY08THZzhSzY2zi1j6TkocBUxyjko47sEk8iZr8fbdZOiONNQtUqEYqfTAhWI3Wcm3ZuA023qetdrHlDoKej549GztuysPbC6wOvtnUkMUtcJXKBQ9Zml2GXs8teRbM7p19/SGHIubORnjiSCJSIkvEmSnp5qZ7nG9c93kjI8a+9SShgq1wlcoRji9c3nYTVZmZRQzxp5DhlnzrZyQOZFFWZNZnj01Nm5ORvdZ8ouypjDZqclAdBp2izBjM1nIsjhZmDUZt7kH/puMqJRETnQTVxn7tKgVvkIxkpmwFFqroami12+dm1FCY6iDTEvcbTLbPZ6QDJNnzWBuxgTqg214wn5KHAV81r4/Nq6zpqzLnE+xPa9bNUxDckqg+QC4CzRdHOWv7xZl8BWKkYwjW/s6AoNvNVlSNGsKbJmx17lWN7lWbbXtCcdDNbMtLl0hkiMy9gAFU0GYIEttzvYUZfAVCgVMiSqomKMhk4EOaNgLramVso4Eh8mC02RjimsM+daua+SmxWyFcFB/PGpGn8xvpKAMvkKhiBv6TmxuGHOcpia55x9HfXqTMLE4e8qRn6BTr97fDhXvwfi0gruKLlCbtgqFwhghwGyBiSdD8cJ4+0BuitozNOPvGjka9n2JWuErFIqusbm1r+JFWsKWPVOroWsyQe12rZZuX9G5EavoF5TBVygUPcOdH3/dmbSVXawZfCP5ZVcBeOq7P+/Us7RavMKknXf0LE3D39OYtii74shQBl+hUBw5jiwoPRFsGbDrNX3fqBmav707Ol1HuvNGo4cyRkM4VYxNcWQoH75CoTg67JnGYZH25Gic6Ji8SfGmxNdG2FzgzDma2SkSUCt8hULRN4xbAGE/HPrCuH/yqYAEix0Kpx3TqSk01ApfoVD0DRmFmk/fCLNN2/C12I/tnBQ6lMFXKBSKEYJy6SgUir5l3ILUjdZELX7FgKEMvkKh6Fsy4kVUGDUTTGbILBq4+ShiKIOvUCj6j1xVd3YwoXz4CoVCMUJQBl+hUChGCMrgKxQKxQhBGXyFQqEYISiDr1AoFCMEZfAVCoVihKAMvkKhUIwQlMFXKBSKEYKQUg70HNIihKgDjrScTgHQg+oLg5ahPn8Y+p9hqM8fhv5nUPPvPROklIVGHYPa4B8NQojNUsoFAz2PI2Wozx+G/mcY6vOHof8Z1Pz7FuXSUSgUihGCMvgKhUIxQhjOBv/+gZ7AUTLU5w9D/zMM9fnD0P8Mav59yLD14SsUCoVCz3Be4SsUCoUiAWXwFQqFYoQw7Ay+EOIsIcROIcQeIcQPB3o+vUUI8aAQ4rAQYutAz+VIEEKMF0K8LYTYIYTYJoS4c6Dn1FuEEA4hxCYhxGfRz/AfAz2nI0EIYRZCfCKE+PtAz+VIEEJUCCG+EEJ8KoTYPNDz6S1CiBwhxLNCiC+j/w9LB3xOw8mHL4QwA7uAM4Aq4CPgCinl9gGdWC8QQpwMtAOPSClnD/R8eosQoggoklJ+LITIBLYAFw6x34EA3FLKdiGEFXgfuFNK+eEAT61XCCG+AywAsqSU5w30fHqLEKICWCClHJKJV0KIh4H3pJQPCCFsgEtK2TyQcxpuK/xFwB4pZbmUMgA8CVwwwHPqFVLKd4HGgZ7HkSKlrJFSfhx93QbsAMYN7Kx6h9Rojx5ao19DamUkhCgGzgUeGOi5jESEEFnAycBfAKSUgYE29jD8DP44oDLhuIohZmyGE0KIUmA+sHGAp9Jrou6QT4HDwD+klEPtM/wW+AEQGeB5HA0SeEMIsUUIceNAT6aXTALqgIeibrUHhBDugZ7UcDP4wqBtSK3MhgtCiAzgr8C3pJStAz2f3iKlDEsp5wHFwCIhxJBxrwkhzgMOSym3DPRcjpLlUsrjgbOBW6PuzqGCBTge+KOUcj7QAQz4nuJwM/hVwPiE42KgeoDmMmKJ+r3/CjwmpXxuoOdzNEQfw9cBZw3sTHrFcmB11Af+JHCqEOLRgZ1S75FSVke/HwaeR3PZDhWqgKqEJ8Nn0W4AA8pwM/gfAWVCiInRTZLLgRcHeE4jiuiG51+AHVLK/xno+RwJQohCIURO9LUTOB34ckAn1QuklP9HSlkspSxF+x94S0p59QBPq1cIIdzRTX+irpAzgSETuSalPARUCiGmRZtOAwY8cMEy0BPoS6SUISHEbcDrgBl4UEq5bYCn1SuEEE8AK4ECIUQV8BMp5V8Gdla9YjnwNeCLqA8c4N+klK8M3JR6TRHwcDTqywQ8LaUckqGNQ5jRwPPa+gEL8LiU8rWBnVKvuR14LLr4LAeuH+D5DK+wTIVCoVCkZ7i5dBQKhUKRBmXwFQqFYoSgDL5CoVCMEJTBVygUihGCMvgKhUIxQlAGX6FQKEYIyuArFArFCOH/A445casq45d9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "y_ab_scaled = y_beta_scaled - np.mean(y_beta_scaled)\n", "z_ab_scaled = z_beta_scaled - np.mean(z_beta_scaled)\n", "plt.plot(coord, x, alpha=0.3, label='x')\n", "plt.plot(coord, y_ab_scaled, alpha=0.3, label='y ab scaled')\n", "plt.plot(coord, z_ab_scaled, alpha=0.3, label='z ab scaled')\n", "plt.plot(coord, signal, 'k', label=r'$\\Theta$')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This yields scaled/calibrated datasets using triple collocation based scaling which is ideal for e.g. data assimilation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The SNR is nothing else than the fraction of the signal variance to the noise variance in dB\n", "\n", "Let's first print the snr we got from `metrics.tcol_metrics`" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[30.98653541 20.07942605 24.94787017]\n" ] } ], "source": [ "print(snr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's calculate the SNR starting from the variance of the sine signal and the $\\sigma$ values we used for our additive errors. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[30.969095787133575, 20.087734900128062, 24.94849587385395]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[10*np.log10(np.var(signal)/(sig_err_x)**2),\n", "10*np.log10(np.var(signal)/(sig_err_y)**2),\n", "10*np.log10(np.var(signal)/(sig_err_z)**2)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the estimated SNR and the \"real\" SNR of our artifical datasets are very similar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "[Stoffelen1998]: - Stoffelen, A. (1998). Toward the true near-surface wind speed: error modeling\n", " and calibration using triple collocation. Journal of Geophysical Research:\n", " Oceans (1978--2012), 103(C4), 7755–7766.\n", "\n", "[Yilmaz2013]: - Yilmaz, M. T., & Crow, W. T. (2013). The optimality of potential rescaling\n", " approaches in land data assimilation. Journal of Hydrometeorology, 14(2),\n", " 650–660.\n", " \n", "[Gruber2015]:- Gruber, A., Su, C., Zwieback, S., Crow, W., Dorigo, W., Wagner, W. (2015). Recent advances in (soil moisture) triple collocation analysis. International Journal of Applied Earth Observation and Geoinformation, in review" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:pytesmo] *", "language": "python", "name": "conda-env-pytesmo-py" }, "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.9.9" } }, "nbformat": 4, "nbformat_minor": 4 }