#!/usr/bin/env python # coding: utf-8 # # Estimating the Date of COVID-19 Changes # # https://nbviewer.jupyter.org/github/jramkiss/jramkiss.github.io/blob/master/_posts/notebooks/covid19-changes.ipynb # In[ ]: import pandas as pd import numpy as np import seaborn as sns; sns.set() import matplotlib.pyplot as plt import matplotlib.dates as mdates from sklearn.linear_model import LinearRegression from scipy import stats import statsmodels.api as sm import pylab # from google.colab import files # from io import StringIO # uploaded = files.upload() url = 'https://raw.githubusercontent.com/assemzh/ProbProg-COVID-19/master/israel.csv' data = pd.read_csv(url) data.Date = pd.to_datetime(data.Date) # for fancy python printing from IPython.display import Markdown, display def printmd(string): display(Markdown(string)) import warnings warnings.filterwarnings('ignore') import matplotlib as mpl mpl.rcParams['figure.dpi'] = 250 # In[ ]: data.tail() # ## Create country # # In[ ]: # function to make the time series of confirmed and daily confirmed cases for a specific country def create_country (country, end_date, state = False) : if state : df = data.loc[data["Province/State"] == country, ["Province/State", "Date", "Confirmed", "Deaths", "Recovered"]] else : df = data.loc[data["Country/Region"] == country, ["Country/Region", "Date", "Confirmed", "Deaths", "Recovered"]] df.columns = ["country", "date", "confirmed", "deaths", "recovered"] # group by country and date, sum(confirmed, deaths, recovered). do this because countries have multiple cities df = df.groupby(['country','date'])['confirmed', 'deaths', 'recovered'].sum().reset_index() # convert date string to datetime df.date = pd.to_datetime(df.date) df = df.sort_values(by = "date") df = df[df.date <= end_date] # make new confirmed cases every day: cases_shifted = np.array([0] + list(df.confirmed[:-1])) daily_confirmed = np.array(df.confirmed) - cases_shifted df["daily_confirmed"] = daily_confirmed fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 6)) ax = [ax] sns.lineplot(x = df.date, y = df.daily_confirmed, ax = ax[0]) ax[0].set(ylabel='Daily Confirmed Cases') ax[0].axvline(pd.to_datetime('2020-12-20'), linestyle = '--', linewidth = 1.5, label = "Vaccination start: Dec 20, 2020" , color = "red") ax[0].xaxis.get_label().set_fontsize(22) ax[0].yaxis.get_label().set_fontsize(22) x = df.date ax[0].set_xticks(x[::170]) # ax[0].xaxis.set_major_locator(mdates.MonthLocator(interval=5)) #to get a tick every month ax[0].title.set_fontsize(20) ax[0].tick_params(labelsize=22) myFmt = mdates.DateFormatter('%b %-d, %Y') ax[0].xaxis.set_major_formatter(myFmt) ylabels = ['{}'.format(round(x)) for x in ax[0].get_yticks()/1000] ax[0].set_yticklabels(ylabels) ax[0].set(ylabel='', xlabel=''); ax[0].legend(loc = "bottom right", fontsize=22) sns.set_style("ticks") plt.tight_layout() sns.despine() plt.savefig('/content/sample_data/israel_daily.pdf') print(df.tail()) return df def summary(samples): site_stats = {} for k, v in samples.items(): site_stats[k] = { "mean": torch.mean(v, 0), "std": torch.std(v, 0), "5%": v.kthvalue(int(len(v) * 0.05), dim=0)[0], "95%": v.kthvalue(int(len(v) * 0.95), dim=0)[0], } return site_stats # In[ ]: cad = create_country("Israel", end_date = "2021-04-26") # ![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAroAAACBCAYAAAAv3xleAAAgAElEQVR4AezdV6xm1RUH9slTpLwkT3nJQ6RIkRIlipSHRFGkYFpohgDGFIMRYMCY3rvpvYPpvXdsejGmGUwvpmMw2BgMmDpDZ4CZE/02rDt7Nqd83517Zy4ze0l3zveds8vaa6/y32vv882splKVQJVAlUCVQJVAlUCVQJVAlcBSKIFZS+GY6pCqBKoEqgSqBKoEqgSqBKoEqgSaHyXQve6665rddtut+dWvftW8/fbby+Q0Pv/8882BBx7YbLXVVs1jjz32Axm8/Le/NYceemh6/sQTT/zg+ZK68Y9//KM5/vjjm1//+tfNGWecsUhs/Oc//2nOPvvs5pNPPhm5nW+++aY599xzm1dffXXkOouz4EMPPTQxr3ffffdIXc/UuR6J+UkW+uyzz5rXX3+9+eKLL3pbePfdd5t///vfzfz583vL1YfjS2Bx2tJU+o3xR7r01piOOZyMXx6S8I/Rx/3xj39MOGWLLbZojVHLgk5fcsklzfbbb98ccMABQ1M8rc+nHOiut956zaqrrtocddRRzeGHH94st9xyzU477dSccMIJzUYbbZS+C1CLSvfff/9Ybc2bN28wKI7K00cffdQcfPDBzQorrJB42H333ZtHH3101OpTVu6dd95J/d97772tbcZzspppZJFy0EEHLRJbN998cxr/M888M3I7b775Zqpz0UUXDdaZSp0Z7Cwr8PHHHzfLL798c8MNN2R3F/4I4OEvaCbPdfDYd/3666+br776qq9Ievb0X//abL311mmxdMUVVyTfsv/++6eAkle+8cYbU7lTTz21Of/885ttt922OfPMMxv9oMsuu6xZbbXVki7suOOOzSOPPJLuz507tzn22GObn/70p80qq6zS0BNtrLnmmqnsZptt1lhoo2uuuaZZZ5110n2O/Mgjj2z22muvZu+9914i/iAxtRj/GceWpoqtqfAbU8XL0tBO2xyOaotd45+MXy7bavO9P0Yfd9tttyX/MHv27HKIE9/H0ekLL7ww+Sa46uc//3nrAp6f83z11VdvTjnllIl+FseHL7/8csLHRn/8IrC/JGlagG4EjbfeeisJ/Po//CGNEUBce+21UyZmUQdthWcyRwXNf/7zn1P5HBwsKg+A/Morr9yqbIva9ij1OSQy6AK6lM7zmQh0gYFFBbrm8l//+tcoolqozBtvvLEQSFzoYfZlOnQma773I5DVBXQBYfP65JNPTrQxk+d6gsmeDyeffHKzww479JRoGhlucnnhhRcmys379ttmjz32SFmDuKmtX/ziFw1/EwTA7rLLLgkYR3b3L3/5S5KjYFSSNu66666J27feemsqm/ft4fXXX5/uy84EycqvscYaiz3IRP+L8zqqLU0VT1PhN6aKl6WlnXIOR7HFvrFP1i/nbbb53h+jj5MA46v7gO64Oi3eb7nllqldC/+S9txzz/RscYNcfGy++eYpsZDzBHgvdUCXkQSVQNd92RD3F5VmAtA95phjUlZoUccy2frLOtCdrNxGrdfmbEetu6jlKtBdWIKff/55yp4effTRCz9omubll19ufvOb36T7jvT85Cc/aW6//fYflANSBR2gFQnI6667bgLAZeFdd911ocWQbUh19ZXTLbfcku7nQNfzyGp5Pgrh98477/xB0Xvuuac56aSTUla5Laj9oMJSfmNcUFCKIxY55f36fYEEFhXoLmhp8p/afG8Fut/JE/ZxZE/SEAbJiY+zW8VXnX766fmjxfJ5mQG6uTTbgG7+/KabbkrbfM6ScmBlsMjLytycc845zT777JO2CB0XMJmR0fXc5B933HEpa6Oc4Iguvvjixpaj8jvvvHP6c8YPXXXVVY2g5ryrIMqYRqU2oCub4wyqM8Tbbbddc+WVV6ZUvvZtUQiAzqZuuOGGTQSuO+64ozniiCMmzmaW50e7eBwV6GrfURLbr7Z583PNsl6e4Rl/Z5111kLDtz1rKxeQ8JyRtdHQGNSz0pQFx4PxR0ZX5oys9KGc87uMGC+2s/EWvDv/hcjOXJpX8iKLUsZrrbVWyvb9/e9/T3U++OCDpD8xDzGOtjH26Ywz0fo+7LDD0phi4WaLHO+O6QAuG2+8cZp//RgHOdJL29u23XMCpGQmzYVxrLjiiq0Z3QcffDBlPumyceCDLCII9M21/rp0KeeFvNgE2T777LPNfvvtlzKUrmzKYtURJUcIQoejvvIxz2wg301gc2zYVhY5bbrppqmabX/t2WozHjIoKYBm7BaVz2OBTcZk8+mnn5ZF0nfA1pmxIMcSHBN5//3341bSJ3zmFP2PCnT5I8caAoDnbeWfBaMNNtgg8VwepzEWsjS3dOxnP/tZ8h95fZ/J3Jgcz2Bf6JVXXkkAXpbc8yE7//DDD5tDDjkk6a5+83P9Fg/miIz5PLy02VKu//SQjdMhoCUn/ocNOP5Fz/qOffX5jWizL47885//THyzKbL47W9/G9WSH+zio82PT1T8/gNAYX61ya/T51/+8pdNuUPQNd777rsvzRs7Cr/HN7ZR1/z0xb0h/trmsMsWh/QneC79svuj6kW00eV7p8rHkZl5Zy/sK2yfb33qqacSGzEOxy7pkPjCV5APXRKv490QuAUeoetPP/108pd0Yc6cOUm3+aM8ozuKTocs2q58EJmeeOKJyS/jLUgW13gs9kug6z7eA2/xN7HDnc/RuLFL396F4LslaPgzn6+99trEVmR0+W4yWn/99X8Q27rww0Ky/etfk2xhAT7KzuaoNOVHF/KO+4CuLT9ndgOMcoYCHYNuIyCU0lFSBEBSoAC6gq+zwUgZTpbBxHeCVF5/oRgmAsCyrYm0HwEz3Rj4pwS62rVlGu1deumljTM1yKSstNJKacvB5Dny4DlQgS/KjwQqzj+oj0f9qTt0dIExmwt/gq+XsYI4VosMFPMlMCIGBXwhBqGdEth4NjQG58CMN4InwCNoB9DVxmmnnZbkAwByaJdffnkaG8eifvCegwGOyfhjOzlkLHDgiS5ZYeb9cDh40T7qGiMdatMZ82RbWpBAF1xwQepDpsgfAEfGDzzwQFrscICInI0RGT9dj3kDogD7aPOll15K42o7umAe4qyaYEwH3Isg0DfXfbqUGMv+Cdly6IJcZErJ1jwCkvSUQwsic47OFVmUOAf78MMPp+9XX311yjb4QgbkRGbGwJaABp+9IFMSnTXXfYthdYAZsu0i88NRBr344oupXbwFcczlkZhxga628ELXIphE+23X8AfxjJ0Zr+MVQRYpeA+fGfddBSflw4+49/vf/76JY2N9dk7mfHFkuskCYDM35G1ObW8jCzVBBpW2FPpvMQHE0Un+LOZZHffYPv1CxskvBthIN7//ZxS/MRRH8Bv+4bnnnks+XvN9fJBHlx/P+fMZyDUnfKZ5BrrNQ+h8Xz/qS3p414PfEyva9KVvfvri3ij8lXOorzZb7NOfUibhO0Luo+hF3kaX7yVLsp0KH/fee+8lvfayVJD4b5xBgGLoKYCYxxE4gb8yNsQX4w2IjM98Znl0YRSdjv67rgF06Zw+2RriNyVfUAl0JVpygKksX87XoZijcWNXqvx9fbojuaNNnwOr+c4nWwziQwKLzZtPNIQfQp58iRhpEa9+W3wMfsrrEgG6HIJsW65kBE3IIficUcoh6wI8BHHoJjmArvOKMmlBwBulC5I1VN4EIDzoT/AIAqLcG5VKoIsnExAZWS8M5dkKbQPoSCDFCwfPgGLF57lAgIZ4VN+YAjClStk/4RjCEDwCXPKMlqCWK4ygRiERXqxMY+VqXBHwsm56x6AcBZV9zKncgsSDvoP0aWz51q/5zIEVvVAmHKq6ZJxnSwHRck7JN4Bu3xhLndE+hyd4BnkRDg8BjMxlLBwARH/AgjI5CHFeVAYNyWbKFuREFvm85M8Yu/byX9sYmushXcrb95ld6QPoDnLOylwG0RP6HsBUkMgXacqRB3CJyEX9cICyvfHZTowMSxeph582QJTXkU3g0LvIwk07HG4QkBeZV3KyuCppMkDX/OordiHKNvPv+eLLfcBesMqzFhy89iL45vXxDWDGYsozgDT8XZ+de3GPvkVZANlizndjyOeUbACNoNyW3DPfsmVBAdjDPgBvbeMXudq9wENJQ35D3aE4UtpR+Mo+Pob8eM6neFWeP6TvFjmorx/Pc79HJ/m0kvrmZyjuDfGnr3IO22yxT39Kftv88pBelG20+d6p9nH8FXCL+OlNNtkkLTRiIRn4Qdxjd5EAUj78Yyxo4jvdgWViQV4C3SGdLuXQ9j2Arn74ugDn/AK5oRLoWhDlNq4M3eO/I75PJnaV/JFhuRtGB2Vyg+xwkGckvvowkDoh2zwWmTc7QqPSEgG6kTnMARiGMW+7s6QIMuoFUSjCIoScKBhnz+kDJUGl4YSwOSTO298222yT3mgMJxx1u64l0JXJNdFA+b777vuD7CfAlQPrvF3ZR1tNAKHsHhricVSgm28fA1SlYzZeSidrLrsgS4Js2QjAVmMWEX6mqY/axqA8IMHAc2oDupGRj3LmN98GtIp2CD+ozaGWMg4DjzquuWPvG2OpM+pyLLJdoTMytgJ3/PIDZ5EDYXXizVsOKepZ+QIEMgvGGdm04LMM0HHftQ/ods31kC7l7fscziUH5+zTQimI/eJddhZ5C7jcEbHQiCxVyIEM7bZEPXXbgmv04+pYkr5yfvLn8RmPZNdFkXHPn0e2mIwAh9D/vEz4IEEmp64zusrQjRi7LIY5jz/HN3Iqga5zucYrAAcJCO4Bom2kD/pv8WDhXAacLjt3VIRO5xSLFwsHR4eCBNdYnLiX25Lvpf6/9tprC80b/ZEBDTtwBVbtbpU05DdGiSP6IzPzTv/Ct/fxMeTHcz7bgKTt49hV6OtHO8BG6ffy9n3um58o2xX3hvhTv5zDLlvs0p/gIa5tfnlIL6JuXNt8bwDdqfJx9MEiC7C1aGf77NBZeeAvFo3xImq+6CQLcT70ts1fGksJdId0Osbfdw2gq0xkR/FrYRu7ySXQFX9zsKluLELtdKByjtwLn90Wu1Kl4p8uoJtjDn6fTebzqJku/NAmWzspkVAqWGj9ukSAbqSqy4HKqvipn5KsikunXwJdTk+m0rkVZ2OAjciYaa80nMiydZ35K3lo+14CXWX0/bvf/S45czzHORXPShDmHuMV9GRNgC5AWDk0xONUAF2GKMgJ7trjdPPsitWsBYOVHyeQZ08TkwNjUEYgy7Pt7s0UoIuXrjGWOqOsIC14dFGbswBM6IK5Linm2EtHOU010I1+RtX3NucyBHQBfrqfE4AqIARwsivDTsmD04+ztF3BNdqKYyQcbx+xf23nQSnKx66RceQUspHlx0db3ckAXdvfkSl29g+oib8yK1sCXTZoHHQwCOh1TxttZLfFc23zO5FFVbbPzoFNuwptRN/LRWpergRJpf6XQJcvkVAYhYb8xihxhN7JRpoLsons1xAffX48570NSFpgmE/6NtTPKEC3b36G4t4Qf8ZSzmGbLfbpTy4Pn5cU0A07HtXHWQzyTXBIxCfz5V2Eu/70p4mMI1uiO+GrjJEvV9fZVtTmL90vge6QTqfGBv4BdM8777xUKmzALmbue0ugS+/LXdXgLc4ll7arg77Y1cbmZIBuHwbSR5tsfxRA1xYNJQFgg6yQOAdGVpItWormPF1QCXRlgvNtbdtrcV5FnQAtseIJReUEc+KcRqUS6MqKxhEE/ViJm/igNqBrRQhcBggSbDkeNMTjogJd8yDzkJ97NQdxdCEPsrKI/nOKtvnpG4NxxKH/kIPrTAG6fWMsdQbfgEuZsc11ps1ZhEMps4Hq0RO2kB+30M8oQFe7QXSFjeSLxzx7P6RL0U5c25zLENAlG4EiJ3YZNpDL2jYU3Yt79CpAYV4/PgOfftvWubI2ctbRGC1ayDN0OC8bZ9ra/iMOGQfHdMpMa9QfF+hGX+WuVbRXXkugSy+AM75A1lgwsgNgjmViusgiQhCIwK3ckJ2TnXbzHbNon7xl6ruoBEml/pdAVwYKeM6zwtrObSj6GvIbQ3HEkZHIVmkzzs/yZX18DPnx4M+1DUjK4gao6OtH/VGAbt/8DMW9If7wUM5haYtD+pPLw+epBLoRr7U71T5Om3Z0+awArOxVUgd2iMV5+G9AOohesZnwX23+UtmoG7hgSKej/b4r35IfdeO38OxMcFAJdONXGuK5q7htDHEcrLRdZYL/ttiVtxWf+Xo6n1Opg2VGdwg/tMl2RgHdYLBtux4IzM9+EaTJsUIpicKVKX9ZC5MU5aXlY7UOAHoJKc/oRvl8+90ZHNnMOKPC+Qp0jEtgkV2O9kuefFc2Vu6+W/nl24UOUXN6QY4klIDGVq/7saVm+8Gxi6A+HvGZG1vUiaszwp7HuTT3rfpiCyNWtLb3kTcngQQBAZFnfg6Uw23LqA2NQaC2RRQZJmOV1bFyDpJRJssggQ/vOWCRrYsAolxksPJgVsrYnOTtqke+jAv1jbFNZ/BjLLleeIlEQEecZn6G0T36a7vemaKYZ84lfudQxowexv/yBdRxXH4Noo04ihLMDc21dvp0qezH7gL554tL9pAvHumVMpHp8EKUrG6ccTNuoCb0iw3bngqiA2GPXvwgo5BPlMmvsbPDceYkSxyByn3b/kBivu0vSAKBbUej1Imf5MnfA8j7iC28PJh4DqiQAV0JArDYWMlnPC+vxkyn4sXVeM4+2RZ9Ejy052hRBM0ol1/NAd1wDRqyc7rL94ZNqGe83iAPmec2Zks3fGZuS+qV+s9Ocj3iA/CXHw8RTMPnBM+uo/iNvjiib3MeOkUv6aDvfXwM+fGcR3NiIRB98MmyduEn+/rRTun38rbjc9/8DMW9If70Uc5haYtD+hN8xrXNLw/pRdSNa5vvnWofpy+JLvppgYyAeomnOJ/rHjwBwMXL7VHPAkF51OYv3Q+gGH5vFJ22uydplx/tSp18/4+jFfGSs1vedcqPHkUCLD9GZk7ElPycsYRNYCbtlHPk3lDs+p6liYuFcR4jPKBPgTl8Fy/IPJIAQ/ihTbb4LuPsBBMtH6bt6AJAJwNoQIKf858RyGOwUDkgyNhlc2J11MJnOpvnvJhfSQAMCI/DFHxtp1Mgb3gbPDAnS0lhI8VPaTg5Rg1gyV4AFCabwuLFFlEAGJkx7eUrp+BLAJUxif9ZSaaLQgt6xsH5MiArY5OEBAeyMIY8i+wFE5neqCdAUUhKh+cuHjltstMmUJ8DEv15Tk6ey1aREflyyoJlAAPAWqbM+C1IKBAelWXsxga8yzTnL7eELFyHxiAIUGYBwOKD7IxXRk/bjA+PeGUkDJW++M5A8O5XOcwT3tUBFM2jMnRCFiFkDDBxlI8//niqr4w5QcC6oK6MOesbY5vOGEtkaWSlzVP8/qngBtjJxgKxwFUQvbLy9mflbO7CkQkkHBtwZs7MDd3yXSaxjegq/bZYpKujzHWXLpXtOzfMLsgNUABG2ahx0RXHCIzHHCpDd7RtcYJ3OsMWtCHQRvaOE2SvfANQFzqof2ectU0+xtJF9JK9kA1doSfOAZcZQeBaf/wEH4BXekO32sgvWVgkcewlOcOrP2MV8AAhxF74FPfxRLfZCF1rO+JTtuu7rUFBSxvswXhygB516B2dzbcn41l+pXN8XK57nvfZuef8Z+gT3Yoz4+RhgYk39sY/xNhKWwr9t7A0x3wsGzE2OkHGSEadL5BNs/gzL2005DfUETS74giQya+ICRaNZBfBtY+PPj9e8km/yY0uAuvsUezJ9bFrvLnfozt5xrDsp2t+huLeEH/lHOq3zRaH9Cf4bfPLo+pFtOFa+l5AdCp9XPTlXD6fk5OYkCeHPDM3/AmfRbd9DqyQ+0uxKBa9eObv6D87cBxmFJ02ZzBA2w6LZ2KMNr3QDF8oF7vj+BZLPKeX+c6OOOEdF/5LPLZDhCcUczRu7Mrl5jNsg3eZa/5Rn/wWrCY+i3Xkiz/4kAz78INkaR6L2DT7FCvwGuMu+Si/TxvQLTvq+i5AYp4CDBHnIejGdkbuTNQF7iLT4Lu2y8DGEMu+1LOdVZKf9slBafm8/B7tMtK2YFWWj++AAL6CLAgCBMW9Lh7j+aJc8Roy1Q5ZGEuMx/d8kdLW1yhjMKZ4A72cu7Y2F8e9UcbYpjMAQATucfi0si/nNurrh86iUeRTzlu0M3SdTl3SN5AlMES2I/gha+Mit5B7PHNlq5H5yO+3fRagQpfansc95TjLtv6iTFzHsdmosziubAtYB7Yje97Xbyn3KFvqS9h5PCcj8mqTlTm14zNVFHrQtrAo+xjFb3TFEWOhV1260sZHjF+/QzoBeADT/KesWRlvYixt/cSzUa9d89MX90blr+ShzRaH9KdsYyq+t/neUdodx8flsU/bxh46UPbFd3XpUlm27/uQTuN/Osi4+MOIM+P00Re78naMzd84NAp+GKe9suwSB7olQzPluwyh7EVX0JgpfFY+qgSqBJZeCTgyAEjxReMGj6VXKjNnZICkHbWZSjOdv5kqt8rX0iWBCnQ75tNZvRpYOoRTb1cJVAlMuwQcEXButm9Le9qZqB10SsALlY7EOKrhCMFksmSdjU/Bg5nO3xQMsTZRJTCSBCrQHUlMtVCVQJVAlUCVQJXAAgk4X+gXMOIvPza3oNSS+zTT+Vtykqk9L2sSqEB3WZvxOt4qgSqBKoEqgSqBKoEqgWVEAhXoLiMTXYdZJVAlUCVQJVAlUCVQJbCsSaAC3WVtxut4qwSqBKoEqgSqBKoEqgSWEQlUoLuMTHQdZpVAlUCVQJVAlUCVQJXAsiaBCnSXtRmv460SqBKoEqgSqBKoEqgSWEYkUIHuMjLRdZhVAlUCVQJVAlUCVQJVAsuaBCrQXdZmvI63SqBKoEqgSqBKoEqgSmAZkUAFusvIRNdhVglUCVQJVAlUCVQJVAksaxKoQHdZm/E63iqBKoEqgSqBKoEqgSqBZUQCMxbofvTRR4t1Cv7zn/80Z599duN/t3n++eebAw88sNlqq62axx57rHn//feb4447rvn1r3/dXHLJJRN85XUmbi5lH0pZxPBuuumm5s9//nN8Hel67733NnvuuWezxRZbNJ9//vlIdZaWQtOlKw899NCErt59990jiYttnXzyyUmfL7vsspHqjFJo3rx5zeWXX97ss88+zQEHHNC8+eabEzY1Sv2pKHPdddc1u+22W/OrX/2qefvtt6eiySlto5TR119/3dn+dM1TZ4fZg+nS16yLafn48t/+1hx66KHJdz/xxBPT0odG/Xe/Z555ZrP99ts3559/fvrf0abSlqaN8cXQ8OzZsxdDL+N38c033zTnnntu8+qrr45feRFrzCR7Whbj8IwCunPnzk2BcbvttmtuuOGGRVSt8arffPPNzXLLLdc888wzqeI777yTvlMKJEBttNFGzYknnjjRcFln4sEUfXj22WebnXbaKfGx4oorNocffnjz1ltvTVHrozdTykLNn/3sZ80uu+wyeiPfl7ztttvSeGbaf5c59kDGrDCduiLoLr/88r0288UXXyQdDrbp84YbbticcsopcWuRr+ecc04K+q+99lpz1FFHJX5ym5pMB/jE+zh0//33Jx17/fXXx6m2WMqWMiqB7pdfftnk96Z6nq6++upm7bXXTvI56KCDEij8+c9/3uy6667NP/7xjwkZjKOveN5///2bn/zkJ0kP/Ze4Qebu4IMPblZeeeWGD9OnP2XVseBaddVVm1/84hfJtwKOdGZR/H/4K3owDgFA++67b+qfbRx77LFpYf7b3/62uetPf1qoKTHqqaeeasSHiy66qDnssMOaVVZZpRHDZiLN+/bb5phjjknzYC5OOOGENLZ11lknzY9kDhp3Lks/zubo0/z582ecGCy86Zb5Wtw0jj1NNW9tPnRZi8MzBuha8Wy99dYJzAnci5sow7/+9a+JbgUbRhFA14Pf/OY3CwHdss5E5Sn+wHHIhC4papMFx/jZZ5+NzdLDDz+c5Fo6yLEb+pFVmG5dEWS7wAF7ostPPvnkQlKzQzFVQFdgA1gef/zxiT6mYsx2DfCurVFJVk+dmQZ022RUjmnzzTdPi4X8/lTOk3avv/76JJ8AtuwY+Nl4440nup3M3F1zzTWp3fXXX/8HvuGKK65ojj766NS+BftZZ5010dd6662XdiXc0O/OO+/cqcsTlXo+AGvmf1ygG00Cf3b0gk4//fTU3p133pluhX7JEAbpUwyb6WTRsdZaa02wGWMB3HMadS7zOj7LcpO9RcBMpDfeeGMsXzJVY5iMPU1V320+dFmLw9MGdO+7774URK2KrR6GVnhW8lacbWQLSiai7c/qfTqoDdyVQHc6+m1rUybZdvCSojZZTJaXZc3AJiuncestaaBr4SPAvfDCC+Oy3lu+zUn3VmiaJoL3TAO6o8hocQDdW265Jc1VAF3ytHNkVyAHb0NyLp+/9NJLKUu70korpYx+/vwvf/lLc+mll6Zbstr50aUc6CpAh+666668+lifFxXolvy8/PLLSV52KZDsrjH+GMkiIwe6xiC+mPt8MTnqXOYyUH+//fZrVl999R/Mf16ufl68EmjzoctaHJ4WoOs866mnntrYLuG01lhjjXTGtWt6H3jggeRIBKivvvoqnX/Ky2qPIXJAsgXxJwvx4Ycf5kWTg7Q6BUq1Jxtiq04GQdvHH398s+aaayaHHCtwW22yCJttttnE+Z02cJcD3bY6GHHcQPaVQ3FWMM8qOMuFH9tGd9xxR9o61idF7KM2oOtsprHow3ivvPLKtOUpE+F8oq2SM844I/WBV3PhDDJZCmqAcx5s9H/VVVelLUxtyL4IGKiUhcWL7LtrUBs/8Sy/hoHRC0cfzIWsoj6CXnzxxbQVKLCQT5w9vf4Pf2i23Xbb9GfRgwQd4ydzZ8PMMb6Mb6+99mpkktoon4vbb789ZUb0jpcAACAASURBVLPIEA21ASgccsghadvV1muAg656pa7QEQs7Y6En6JVXXkny2GGHHZIOdbWVCjdN88c//rHZY489UkAxX7aF2zK6Dz74YKNNIJRe0HP8ILqIfzJyFMX3MhPjjLo6tmbJuO3ojHOxbEMf9EJ585zblPnt0k28yCCZN+3gk+1efPHFaf61qy1/bbsIdBt4MudHHnlks/vuuydecqALZOrf+c299947bV/mgb2t/5B1l13E87gaI7/nfLK+jjjiiHTm3/M2GeVjeffdd9P4LFg22GCD9Pnaa69NTU/VPAWfJdAlBz41do1KfY16QzIECNkVeZmz3PdFsiLayq8lsMyf9ekf388O+VO6E2dyA+jysXxI+PtRz2yX/PBVxuP9DMmYX/7ylyke0Ud+MvwSvUX8KtumhzKowJ/vbAJ12baFhzp83tNPP53qbbrpps2cOXNSvfjHImDcGBd1S6Brt4fvyLP5yk5mLh955JHkg8gI2DXONtKnuEUmdI7uk5mjLO7ZcRiKV3zuaaedluZXHPz9738/0VWbXnzwwQfJP0RsVFgyjh/mswB7drbuuusmvzPRWNOks/5iibl0xOfRRx/NHy/0uc2PlPbU5gstPvjzv//97wu15x0Z/oSfJtfc/47il7p86ChxuEtPF2LwR/JlyoFubJMGECEHTh9Q7Tq/ZCKdieKUnNXiVASjCATAGXCQE+Ux+W3EAKy4BXBOzwsy2owXZCiL7an8rA5AoAynhiij731HF8o6lFSgCmUFpFdbbbUJByerHcbEyPDG8eClL+NdAl28kVPIU6bkwgsvTHyTv7FvueWW6cydxYDnAq2tZcSJkDcjCBJUOc9oExgO+ZayEBi32WabicDYx0+0H9cwMCD8008/bUKGEaS0zegDdArKxqAPREcsnAJcuufMX7xgwFmaf0R/ONx8DtODpknyNhdkb6HF8QoyqK8N/GyyySYTTtz8xvm9vnoxztAv4Jp+AXRBnLWgOcQDkGvxxnkjTlpbbUCXPOPMosUIObqHjF+AY0uOknDiHH8Q3sg6+rngggsaGcdSV+mTII0HdhpzVY65SzcFVU4e4U3ww5N2gRXtCoTRbvAXV2CDviqP6I46AXQ5bEAu5EN3BDfbrKirf8/67CJVzv7h5+hiEPsB3MmrS0ZRVhnjMx/4SvP0/Ximap6irwC6FpCApMUBUG5+gsq5G5KhegGOzKH5oKPxUvFkgG6f/pEPv3jrrbcmlsUbAJQcA+jGwiz8vReRRqEc6GqPj2YHgDLdoUfhk8IPOd7gHmJLAFvYmfrmNfS3z0+wUbprQRafgZ2SJhPjtJEDXXLCt5hVgrfJzKW26RDbNYZIUJS8++4IgTgVft49PoQNoaF4pS9/CK/whXjbpxcSIeYIHggSh1ZYYYV0jzzMbe473JMEYA/I2PBt0VdSnx8p7Sl8IT8ELwDnfGvuP8jD3JAVspCP5NKofqnLhw7FYf316Wk59pn+fcqBLsGarOeee25i7JTXqjEMfeLB9x9MsJXre++9l+7Y5qJ8MalxNiqvB5S0KZsylJWCBHE8lJeDD+JIrMiD4qB6ABG8qpODpDyjq15Zh2IACzmRhUAV5LvVYVA4hfx8cDyLawl0BQDGFuDOCx+5owLewoFo11icz5T5CLKQIAPEIatjhRhkEeAeapOFupEBGuIn2nQNA8uDKhCSByFgNhY22jYPZI3++c9/pu90BAk0ceQlwJY6QbLGsj5tZC7IAQnI/vraEPSs+OOXN0JuHGdfPe2XuqKutgKUK0PfyXqoLYFY5jAn+h5ALr/vM6BKhkBNTvQyFjPuO2LE7mKxw2451yAvamqnTVfL8anTdq9NN+kq+49z2/Q6nLvslT7JpY30IchZrASFzgTQFTTJJ2+DrNiQPrv6j/ntsovozxXP+Ay9dU//7tF51CaP9CD7xyJKdjqnqZwn7QbQBXIskPk1L2HldlPyOiRD7QY48hnosMgM/ZkM0O3TvxtvvHGhObVIBEbNcQBdyYQgQD5fxMX9tiv7AkxlCiUMyCfigvJ2ywLURv38nv5DH8O/h+0N2XbojLngb5Rvo8nEOO0Ah/TefP/0pz9NC77Y2cz7GXcuJRX4U4RvPj18a95u/ll58g2SDdUv6otXZORlxrA1u2HGBHP06YV2gdYc6JZypD/a0g6SeKBXfAFyhWXiebr5/T9dfsTj0p7c4wvzHUeJhIi5notbOZ6QqBNzx/FL2mnzoUNxeEhPtftjoikHuuXggTCTlwf0sowgVL4UwylxJqFgeR3OIwdt+TOfKW9kL+OZgCOQBwnwuZGVikjh1RkH6HppLAcO+qK8+TgYdzh/z72hrp88yASPcS2BLiAiIArwHBZ55ETeeXDOn1FgcwFkxepZpgIPtrgZkz8ZW46Q/NtkkQPdIX7y/sPAAtR4ZjupnH/BQjaPYeMNwA0SeLwJjWyPxmo73iS14IhxCFj5wiLacC3nwr2+NgQEvOTnByOj01dPu6V+uec4jbmyOJRtCIDT15bFIB4ik6UdNFmgm8s9jhDF3AhWsmQhS4tLTj9+meT7rtOlbXxt99p009yyEcCIXf/73/+eaLrNSU88/P4IB3nkW3ol0BXcbQfnFABEcOzqf8gu8vbMBz7yBRzbYaNxNrVNHnkbPncB3amaJ30E0M1BFHumQ3Gv5HVIhtrNwZHvocdA4GSAbp/+WQzTzZzCFgPo5kcnLAz9rCGSIOEX8r/c/wK6kQ3WVkk5qI1nbffwYwEXflbZkEmXj8qBbrTddp1MjNOOeYwzuvgC2gDFksadS+PiH84777z0x6dbNEdGv2zfd8kZNiMGokhY5GXb4lXshuXHFmPu+/RCu21At8QKjrpEosgCyffwga7kFzad89rlR5Qp7cm90hfGYjLadITJYivIAkKsGMcvqdvmQ4fi8JCeBk8/luu0AV0OglEBS1YtJqmLOBYKmlO8vVmeT1LGKqdtFRr1J+MESkVsA3dDGV0g4He/+12wka7OxQp2YYgluJoM0NUwueiLEXIWcZ7Ps9KA3AMELB78PJq6wHZkOjkTbThj1UZtssiBrjp9/ORtDhmYsvRFlp8jDNCSZxH9hq9VN1Bx0kknTSyGZHWMoy045TzE53Iu3O9rI+SU7wxEW331lCn1yz1ZS/wC6uYvxtjXVvBwzz33RNfpOh1Al27F9v5CnbV8aRtf27023dScjKjAa16B3sigtTnpvHsZFzLMA2roDNCAgIo4Pxl1I8jGmeS2/kPWXXYRbbmaP3w4jhNED9m+c6uoTR5RNq6TAbrjzJN+2oBuyCN2VkpeR5EhuQM6OakHSLDZ/DhbXob/z8FgPOsbF8ChXhuRu7noArqOwYgR+V8Omrr4ib7aQG3bPX7M2GOnUv0+2/Z8cQJdPps9WgTy8TmNO5fOIXu3Qj1/dtzMAfDWRRaCZC3xYvs+z5L2xauQYb64jT769EKZcYEuvZT0GZXa/Ii6pT25V/rCEujSfxiqpHH8krptPnQoDoeMR42lJY8z7fu0AV2GY+tIAKDMeUaiFALgmh8j8JyCAY4lCUxxpKF8Ft+XFNAFhEuH7eyn4BVUgqvJAF3ZkfhRbtlUi4S8j9KA9I2PXMYyFl4wQhGQy0AUi5MhoDvET4zddcjAZNg4SC9jIIDH9zzLZ5sMsJMBzc9ZR7CO7a/oN8YR3+NazoX7fW3YnbBl1uZ8+uppt83RuW/xITud71D0tWW+gad8y0s7owBd7eZkSzy3yzKjS5/z3Qd1u2TZNr62e226yREHOWbhP2pxLh+FkzbuNrIlTD8E2SCB1r0AuhabzovmJCOjjOxeV/9DdpG3F3MmCAWFLkf7bfKIsnFlx/lPb7k/lfOkvTagGxmcWNiUvA7JULu2ksuFPj9lziU8Sv8SY+4Cln3652yl+WsDO0NAN/rtunbxE+XbQG15D7hlk5EZlOgAvENPunzU4gS6xoMncozdpBjjOHPJhiJpEvVdAWhnUPvIwoh+OEYT7+Qo3xev4uxyuaulXp9eeD4u0GWLAKdMak5tfjDsXLnSj5X2pEzpC0ugS3Z2iUsaxy+p2+ZDh+LwkJ6WPM3079MGdPOBx0HsCLSyUUBXKLYAa5sjztxa6VGCcLp5W85ctjm3vIy3NmWFgiglY+aMgmQ28yxPZNcEJySwqpMrL8ebbyWUdWypAOcyBohzYyT5qhYQBjKDrGT1kwfqeBZXW3jxopR7Fg+5Y2L4tliCBPUSCHE6sisIcJUxzZ2TDK2twNi2ZtheUiGHNlmoyxmhIX5Soe//CV2QTQgyD3Hkw+KIPKyMkRe9fI8zm1EHuAf4bOMEkTdZeTuWDiE6lYO5KOtazoV7Q22oQzdji5oOa3+oXqkrwQedMQ7XoKG2AG1zBXgjvMiEOuPWRnjUR67/ynGkMuJBEfRibMrb2oyXK5VzpjO2GqOea9v42u616Sa9jDOM2hKsAC8U+pIvdNKD7/8hK0d78sVH1Am+8UE+caZPVQA+7KGv/z67yPlgU0Bq/oInYCewxtvnbfLI2/DZnMQCNJ5N5TxpUyKATcURFLZOHuY6Muklr0My1K6XJPMzhcF/zEcX0GVP5YJK3T79o4MWnfkWMp2xQGYXxpcfOwPAy+MrwV95xU/oRvnM99hRyp+V92QW9Uc/kWQA2xmybdvfeO+LB9qbTIxTj28sE0h8Gv+Ax6Bx5tIuQB7jog22YCxhh3E/v/Lf5jF+Yzme9cUrtsauxNHw8/rg//r0QtuO7OU6U8pRGfKJ+GmXjWyUC4JjzHdJfX6ktCd1S18ojua4JXarApOoI/stRo/ql9QJ+8t9aNzrisNDelqOfaZ/n3KgK6jmAYUACJjCh/IAsAJPDlgpH4dOWZ2lAigjQIQQHVAvjzjEs7jqW339CRiMghH5zng4cj/nxVCAa4ARGLL1ooxfHmA0gKPvACGn4793pISAKydQ1rFiA6htU3rjWDDhvI01VoMcsfraIQvjd2BfP+qUvwnsbVuA0nPG5rM62ga6yULwkP3iIBFDUN75njywyOL4hQDOgdOXCbWtFluNwA1+yUWGkaMmB6C9lIWzs84pGYezwH38xLy4Aq8ymPjTl+1m8qILgovABFSbP/wbr4UG8Oot5jwLYp7JrCQ8OxfnT8YW77GgysvGXMi6AKpWyUF9bVip65ecBAgOJ5xFV702XYm+9Mv55v171tWWZ7JkFop+eQPAI0Nz6zvdbCPzab4BEY6avZG7eTSfdJzumxsOm8wEkcho0FPjbXsxNB+febJwze9pl3106SYZkilb9KJIvJRnHPgAZMmIjeY+I8Ypg0tf9EMeskPsZccdd5wAbkC8M/mCiUWVc9sxb339d9lF9J1fZXP5MPPBvn02jyiXR8gorxuf2ax5UYZtTeU86QNfskTmmV2xMXpD/vEfiuS8xtyp2ydDcwsgaFdb5iQnc5r7I8/oCT+gjvkKnxP1hvSPTwudVleGj7+iA9rkR/h7PsSY+Xtz00XACJ+hLt/Gr+ZHGtTjdyzMlVEWsMrv0SVg23PvPFhIaod+xi/jdNm2LHDIg9yBkTaaTIwTg+wa8Xd4I684m2uM5o7e8ZfsY5S5FFsBZ0CVfUYswbO5EG/1xbe3ZV9jbHxS+Z7JULwyr2Rqjvkli+MAvW16oS9l8Erv+cBSjhZI5gvPfE7Yg1+P4CfNp/HmSaYYg2uXH2mzp/CFeDHP/qOdkFdgHGBTQs67C3wfvxwL0XH8UulDZXOH4rDxdOlpPuYfy+cpB7ocKUXJwUWkySOjSziR9cwFBeQAQ/k5t/w5BypgznQCWgS9EqhPFd9h0GScn00cap/MI2OrLGPhrHJSBpgbhybLT18fflsx2nVtG2effDnvXAf7+up61tcG2QVQKuv31SvL+j7ZcXCg+EBtW2llX2TIxsYlDrdchI3bRl/5mGd6F1nqsryxRrnyme/Gb0Ed42uTh/q2hkNm0U6029f/OHZBVn3vEES/XVd6Oxndne55wm+XDLvGEvfxVi7m4tnQtW9c+JEVjDkcamsmPR/XT8wU3hdlLtvGwGbb7HWUeMXO2vznVOsF/ti1sXdR6GCfH+mq23ef3YiHbTSOXxryoW3tu/dj1dN8PFMOdAFVq8E4qyYDIyva9kZlzsjQZ9sq5ZbeUJ36vEqgSqBKoEqgSqBKoEqgSmDZlcCUA12itMK2ReT4ga37eLN5UcQsUxxvpS9KO7VulUCVQJVAlUCVQJVAlUCVwLIhgWkBusuG6OooqwSqBKoEqgSqBKoEqgSqBGayBCrQncmzU3mrEqgSqBKoEqgSqBKoEqgSmLQEKtCdtOhqxSqBKoEqgSqBKoEqgSqBKoGZLIEKdGfy7FTeqgSqBKoEqgSqBKoEqgSqBCYtgQp0Jy26WrFKoEqgSqBKoEqgSqBKoEpgJkugAt2ZPDuVtyqBKoEqgSqBKoEqgSqBKoFJS6AC3UmLrlasEqgSqBKoEqgSqBKoEqgSmMkSqEB3Js9O5a1KoEqgSqBKoEqgSqBKoEpg0hKoQHfSoqsVqwSqBKoEqgSqBKoEqgSqBGayBCrQncmzU3mrEqgSqBKoEqgSqBKoEqgSmLQElhqge9NNNzV//vOfJy2IUSp+8803zbnnntu8+uqroxRfJsr85z//ac4+++zmk08+mfLx/vOf/2yOOOKIZuutt27uv//+KW+/Njg9Enj44YebffbZp/nVr37VvPfee9PTyQxq9eW//a059NBDm6222qp54oknOjmbTlvp7LQ+aJZWfZw3b16vfc2dO3eRZv/jjz9uzjzzzGb77bdvzj///N62Pvroo+bkk09ufv3rXzeXXXZZb9ll8eE555zT3HvvvT+6oT///PPNgQcemHzbY489Nmn+lzR2WmqA7s9+9rNml112mfREtFX84osvGs4k6M0332yWW2655qKLLopby/z15ptvTjJ55plnpkUWAPRPfvKT5pZbbpmW9md6o48++miz+eabJxlvu+22zbHHHptA5K677trccMMNzddffz0jh/DAAw8knt96660Zyd+oTB1++OFJ//bff/8UyFddddXmF7/4RXPiiScmAMAfmId33nknjbdvQTbdtjLqmGZyudLnThWvi0sfp4v/Ug5//OMfm/3226+59NJLkz8ANHMCNk899dT81tift9tuu+app55KAK0t5pVjFSs33HDD5pRTThm7r5la4csvv5wSH7v++usn3z1Tx9nHV/i2cYB6qRtLGjstNUD3/fffbz777LO++RrrmdWsIPbkk08uVO+NN95YCPwu9HAZ/MK5/etf/5rWka+00krNrbfeOq19zOTGZbbp4vV/+MMEm/Rws802a7bZZptmUTM3E41O4QeZADwvDUD3rLPOmpDMeuutlzIcbtD9nXfeOQFdAdF4+4Du4rCVCUZ/hB+6fO5UDGVx6ON08l/KYKONNmqee+65dPuggw5qrr322okidhV22mmnZt63307cG/eDXQr6LBPXRl1jldFdmoCuJMNQNrtNPuW9vfbaq7nuuuvK2z+K75IpdGFUoNulG0sSOy01QHeqNaZrsqa6n9resASWdaALLHI0OdAlNQuM5ZdfvpF1nGm0OIBFjPm+++5rTjrppLTNKgNVBmeO+sYbb2wAVluIDz74YFQdvCr/+eefT5TLga6bL7zwQnPXXXc1owDdiUbqh1YJTKfPXRz6OBn+LYyuvvrqH8gD0LzgggvSsbCrrrqq+fDDDyfKOALDH8yePTvdO++881J215d333232WKLLRYqP1FxjA93/elPDb/bRV1jrUC3XWKO903Xrmd7j1N3d6qA7tRxNH5LUwp0Kf9uu+3W2GLdc889k9EJErb93Lv++uvTKtOkH3fccWnV6SxfHkgEqdNOO6056qijmuOPP775/e9/PzEqxn7IIYc0J5xwQtoGiPNwtnOd43RFAp9zRe699NJL6dzQuuuu21x88cUTbfnw4osvNocddljqS3bs7rvvTs8Fwh122CE5E+cMZW2c/41zh7YggyiBLaIDDjggZXqcKY3zqqPyEW3FlVPWnjNPxxxzzERWzAq9S3aCra2m3/zmNw0nyeGsvfbaKbh/9dVXSZZrrrlmmguOEj300ENpnJtuumnz+OOPp8/kdPrppy+UDbjjjjvSWdk4qxNnlJ/+61+TbMgu7mnXyu93v/td0gX8kJ+/l19+udl9991TJlJdW294Mm90p4s4XFt0+jcmOhZjiDqCge18ZY4++ugEPDwjsyuuuCLplDOUMh14oUOoSwc8s5InR2XJwLacsQ6dBXdefO+9907nNl3/8Y9/pL5cJzP+LqCr0d/+9rcJ7AJayFyTJ12VRTD2nBwBYUN0y1+Awq56f//735NcjfvZZ59Nc7bGGmukK7uVpQD+2Jo5DQpgQS/ZknnDV559fvvttxOPBx98cOrDMQ2U203oMh/SRsYA6Bs/4E+X2b7sKcKjrDd9dM+uDz3gYyZDJdCNNgLo0hO+K2zNGFGbreAFkD7yyCOTnuG9jXI9vP3225uNN964ufLKK1PRLhnaumZ79IO/0vYvf/nL5rbbbkv1+C22wr95fsYZZyT9jjnsaldlfoNvZod8TvDimXN87ItfFQMio5+Poc2W2nxu8HLNNdck3TEe/NKJNuryU8qOoo92BcnEeWt2a7s+9Igd6Z/tIACVv2EXxtjHfxuv7tHXVVZZJdlqXsaxBEfxYqes/K7sBhtskMbkM/sBls0p3WenQ9QXt8QcumIRbS7505z6xspfsknyMgbfLT5z6tKRvEx87vJXffw/8sgjyc/Tef6N3+U/6A5/hIb00YLB2M0PWfssa57XC1u88MILF8I9jpHwB3ww3GM30lyGjx7F7ufPn5/mla+ihzDGz3/+86SDME1OXfFmlPg+pNf6Iesyo9tla226sTixUy6X/POUAl0NS08DJrnzo2gcCBIYnXNDQAjwkANQQSsyVIARYwNqCNt2TWxhM2zGSCE4I8GMYw3iuFdYYYXm8ssvTwrmHJ3Jev3117/re968Zq211prgk0GtvPLKqR/txbkUyqJv96ygldFmkHHZOgpi5BwivtAQH1EvrmTFuMgRMZYA8EOyE7zJnvIyKnwaM9DsjAyHvM466yx0xlibygh2xihDtfrqq6fsmP4BHc8juJgbDj6IYnuuHrIYwUOcGSML5xnJEJGn8tr54IMPmldeeSWVNz9dpD3O3JgAYsDJOcnYmuOABNwAUcCseUBAk2CEgB46YUHDARpvlw4obw45auBfO/rHNxnG/KaGs38s5uhpLN4YOXlGRmYy4+8DuhwvecY2JoceII5D1XdsOdHxTTbZJIFhLK+22mqNzA3qqxdzbBFnXoEG56aBW4vNTz/9NOmEYBAUwMLckDO9JjegDpGlQKhtBNSYZ2ADhd3QZYGE3WknJ23wD2QQZEFEHvQKXXLJJcmelA3ycpIydHtcGgK6Ae7C1ry8GhRyDFvhw+LFHXPVpVehh547awpgWjANyRDIBbgBH7ITEI3b+BFbIvMtt9wyzQ8Zk19fu+yY7YWtKW9+EB9hEcSukYykrV/8xxi6bAl/bT5XDNhjjz1Se8qQbwDgdPP7f4b81JA+Wug5Rxl+iH+g317GCuL3AJcgYJc8xZQu/qNsedWf9vTBxwfxG+YsZBr3+SE2GkR3yEI8lMEFzCzI+fNRqC9uGTs5RDz0Pae+seLTQswcSfgYI/Ad1KcjUSauff6qj3/1YQpzE3GHrvMVMb9D+ui5usZCB3wWb6JeaYuBeyRcggBSi6aSRrV7NsZHArv8Iru0iIafIr4OxZuh+I63Pr323NjJMuJIn6116cbiwE6lnPPvUw50NQ74cZ5BVuQcFnLmNbJpvkfmy2cOQwCNFWmAIEHc1iMAGIor08upxncKlQNdCq18kHKcunaCKE30xQBNJkCIOGvfyzcNKV4AXVlMZaIN9YzBvQgmo/AR/LjKGFC8ICA1xtUnO+XLvjgavHAYQeSUA5II+pGFVo5hc3IMi4GZz9gms4Ahg6A4ZB7B27wAV0EAJ0ccFPLJV6WAKSDbReYtFjjK6Mu4AFaGxenkDsYK2D1EfjLHQQJvfoasTwfUMXaZ0SAOXN+RbYn7rngBnIGrII4RLxEwJzP+PqBLn/EjuwBM+hyLEjx4QZNO4QPICN5CbnS5r5422ng2p7EgVcZCyTxFUAxggfcgfQk2dMoxDPaLD+S64oorTthnrstAQdhltBVXPMTOjnvqkUHoFz5zfVdGe8rIpo5LQ0DXgijIwiAP8qWt8H34jwUb4BSfo4240kPlkSDnb0iGdC63PXUBEQvFILoZSQk6zU/2tUu3zHPs4HjpJDLxwEcO2mzVknPYypAttflcvMnMhX/SbyQBYgyuQ35qSB8Bgjy+aJMuGWv0LeGQA92IGewDtfGfHrT8oy3+3JzmMgvw/PTTTy9Uy1yynVhgRH+yl2yH/fPt9AKQssDs2iUbJW6xZzGgi7rGSr8iyaCuHQTJheB7SEeivz5/NQr/d955Z9K9wAfaNb/mNGhIH5WTGCj9RJstKmtsQGmQmMS3ljSO3ZNnHq8kFeikxb15H4o3Q/Edb0N6TYbsOIDukK116cZ0Y6dSzvn3aQG6HB/BvPbaa6kvWyElUQCZJ8GXgiDGql5kv9yLwAmQyODmFM/cawO6kTmOOlbK4dTjHjAnDQ9M6tuLP2iUyQK+1MkdCuXjkGQ6EMMahY9U+PstKduCQQy+DH5tsuvqC3+xXakMJ5QvQtoMIZxtnvEyJ7JPtn9sQweVwRvo16fsDJKNcYYsKEBTDsZkniKIR7n8yrBzoEvGnKdsmS1W/QneQK0/2f2f/vSnyRFY2cvUIHLk7HJg5H6XDnjGqeWBiE7rL+c/Nd40KWPuWQ52PAPkZeHQZMbfB3QjgymQm2f9y/yELGQkAPU41xdbd3gJ++mr18WzMQFyQcasb5lJ07UqHAAAIABJREFU1AYsBG9lyF9d9hh8unLafXYTfXVd2YpFogx+2Iysddsiik3GNnRXe233h4Au2wmy3ZgDzdJWQu70k56G7KJ+fi310LMhGbYBXUEzX4iWi8ShdgEWwZ+P23fffRfKrhoHHx1zKutsMRNnE8sxlLbU5nP5Y4ALzxIk//73v3Ox/OBzl58a0kcLjtj5iUZjURu7JUOAoI3/aCu/inORsS2BLpDKRmJxEPWAprb7nvPTfJ5Fh90kIFVyia9uo1Hi1qIA3RyYxa9dxGJhSEeC3z5/NQr/owLdId/eBXTzesGzRYc5CnBLX9poHLsvga72+BTJm4gLffFmlPg+pNcl0I0xddlalx3kQHeUORwXOwVfbddpAbqAiIAAyDLCPItqcmQ5bGfPmTMngYgIOLKBFCXPAgXTnKc2u2gyQNe2qAwfYB6r83Awo0yWVTN+rbKCgCZBILZTx50swT/PlEW7rn2y87ytL/wxrKBRgK4stnqCg/HYbgYgBB0ON7Kl2iyDt3uyI8C0wARkxDaLZ5MBej8Aut9+m1a19ItTwSsn00aMlH5ZTPlz/jinPh1Qbig4523Flk4Odjx3fGPHHXdMRScz/nBo5ctoGuSk6Jt5CvvxuaSQU57djzJ99ZRp43kyQDcWQbbwzYXg3EVtutxV1n07DnyEAJDrG10N2Ud9gNLOUVcginJt16kEutoHBPhDOgyk5L4k77/UQ8+GZNgGdGWnAEeLAtQGdIfa5bedeear8M0PIt/12UXlGEYButqSwcMTP4D32D3K+xnyU21AN9dHi8MSGEbCJs6YDgGCrpiR88mHWiCIkagEuvSeTKPPqMtvuS+blhMAyddaAJBL7LaZX4uDAF15nVHi1nQB3SEdCT77/NUo/C8JoGtOyd+xqxL3xLjiOqrdtwFduwoWkaPEmzagm8d3/AzpdQl0h2ytyw5yoDvKHLbFAPpTJitDpn3XaQG6OpTF40QZaJ6pkNnKtxJlUr24gOL8Yp69C+YpD0NvA8HKjAt0rdK1F1tEsR0eGYOYrNiWCz7yyQpHmDuTaDcyZ+NOljNbDp23UZ/slG/ryxjHBbrKAwKCrgyb7FeAJ9vP4Uz12QZ0AQ4veeWZ7hhPG2gayuja0s51wpEW4+KM8QXoyYTkFIHcqlOQtICxZZ1TzFWXDig7FJzz9rSPlxyQcn6Cs5cv0WTG3wV0jV0wi63C0Mc4JhS8kYVsjzltW0T11evieQjohmxze41sFRvzCwicVmRec159btPlKFNe6QagQudKwqesbuiD55GpK3WmrNv2fSqBbvgI/ThqYS7ze3n/pR56NiTDNqArC5wDujag29eubHwcY5LdtdMm64Us6MpMVy73cgxdQDf3ubk8+GRbw2FLuXyG/NSQPnrJN9+p0nbslgS4tHA3xqBIjrBp1BUzojx7lTDIF2Il0A3dLIM5PaYfOZEtef/lL39JtyUY8qMVkji574+6Ye99cWtUoJvPlfZLYFZmdId0JHgcxV/18e9YmxgRRya023Z0IdfXUh/VodvsIadSj/NnFpJ2qeGVHPfkZXKdHrL7Up7akRW3WzRKvGkDunl8196QXpdAd8jWuuxgurFTLuPy87QBXVvKAmv5xqbtoThUT4CMMTK6vhMG8BsrXqsW2zCUUHuEHGTCAqRoI7aHPXdIG8DIyTaaDB6KrVbZAuSlHIYR578oKdDC4HOixMEDfhmCbccgwdMYAlQN8RH14gok4YNjDpK5tHLvk52yZV8cYQDCaEsmPQ90YQgAaxDjj21pIEoAiPmw+iODIPLK+Y1D/wye7Bi1bEqQjIbyOSihD3EOOcrlV/OQ/2C5DCQZc4bIIseWaWyPAU+y0Jyc7Kpgbt4BIueMYo6HdEDbQHLOW6yic/5zXgX+/Iw10Elv1UOTGT+QTmb0IIh+4U02kB4iRxE4QVn0mC9BOrYSlSeLWIDQcc+G6rXxLEsaC1R9kyseIyMZ2bLYtlbGAkg2CxlTWhRcf3367h9B0wtTqNTliULZB/PMvxiXLLGFcvwF2I9jOPxHEPt1tCXkYHzk2DWnUc+V/PLgGM/oovHHOTb3ZT3z7fDSVmRwLcSCZLhjoR334lrqoftDMgR0LZpDF9iD4yE5+GHb4ROjr752ZWLyM4vkHb6CvVuUhq5rzxZ9yL4cQ2lLbT6Xb8jfkwAgcv6D5yE/NaSP5kbGWLkg8xyxyj3+Iz/Ha7zmPMbbxn+05Ypvvpe84o8fs3jy3XY9n62MxX+QeQNyY5cw7vOJ+cuO5s3iEZlzvLb9z4SjxC12SB5d1DVWiRo/9RcU9he2NqQjUc+1y1+Nwn8kB8L/WIzwN7BE0JA+Kmc8uZ8LvvKYEO25sl860XZcM8qNY/eArvgWZHFlXsw1Goo3Q/FdG0N6Tf+MKQD6kK116cZ0YCc+bpSduWkDugQo4Fuh5mTrVIaFoggEjFXGMM5xyqz6OQ/nUCgixxaOWl1ltStomiDkjC0HDtjaWues1Dc5lFQQYny+UzIvAZg8ZfQFJJtEIMGKOIKkPqI/QRgvQItysZK1qmQMnJCA4XM4vlH4yGXjM9ABjHJstmE52tiq65Nd2RdnAAQYs2CrDb8AwLE63xrBKgzBIsEREwZK7pQVyZ4K7lbijMr5MoamDGPDoz786kGAZdk4/Ktnrj1Xl9MVPHwHLNRXFujglPNMaMjF3GsHsLKgkOmQ1Yk5Uo4T1a6xCRDmLeZAfX0HP3SEwxN0hnRAGYFDHXPLecrA4F+GJc4hB6+ugB4eBC5gjdzCQUxm/ACc88f6xIf+rcA5P1mLPGOmf+P28o4/mQd8xFxaaeObnMyfNmxD99UreebIjct8mTe2p0/jxKOxmw82xr7MszmWgXBeO7KB+sQ/uzU+4Dx0MtdltplnbhKz3/8TmX39ln/GFyS4sUs/08Q/sIcIgMqQMZ3wvItkpkJ3lc11zJvyMuV4MGa2Zs6BTLbGN1hwlLaCJ3pNtwDCEsgEL6GHZG5hErsrnnfJ0DNBgP/iAwEXfpO/DZ2xcMIzH1hmt7vaNZfhC9SRXbVQQGwtssH01BzYQkYxhiFbKn0uHaWzdMNukqAWC7vU8Pf/9Pkp+j+KPgJljgEA84K5s+1hH7oxh0Ao/THf5o0u2KGMnYuS/5xHZcg1/+PL/LkXi0Tl2I35sXCyTY0f/ipI/BFDIzbGffIRP/hx8bWLhuKWcdINPiQWKmVb5VgtIMUGNs0vWDiKC9oRx8xDn46U7ff5qz7+ox38iR18iGy3RJfElF2JUfWRjhsTXMB3RL02W4x+6WvE7LiXX0e1e3UAXQDRzrgEj+98UVBfvFFmKL4r06fXEoFiiDmUlDSnQ7amzVI3pgs7kYdfgQmfFnIpr9MKdBlmGwMCQ2TfMCQwls7L6jayojnTDEW2uDTwvMw4n/0sS7Tlmm8racf33MF0tQ304HmqSDDDW0mjyK6s0/c9DAEA4WDD2eZ1ZM4YQ5CFQ4CnuBdXL6xxwLGCd5+ztkDIA3SUH7qakwCteACm23RKO2TDOeaEF0EyXrxS16JHsA4a0oEoN87V+AH50K1x6k5FWZnCrjnCWx7A8/766uXlhj6zG/Nl/OasnJeobz7YTsxP3J/qKz4Asi4bpTuLm/AU418UPYk2ShnG0QX+S8ay9LFD421rN/ikW6WvjPbw0bYIjOdD19znRn/0J3Zwuur3+alR9VF/sn+5/yr7E39ClmRUUs5/+az8Xh5dyJ/jge8b12+q1+bH87bj86LGrXHGGn26jqMjff5qiH8+KPxg21zlPHV9Vj/a6CoT9+lP29GaeO6qTNhW6Hf+PP8MyFncmk/xuas8GbXFm1Hie/Q3pNdRzrXP1qLcqLoxNIfRXtuVHoUttj2Pe9MKdKOTep25EghDYChTQVZXVp453XPPPSmLmN9bHJ85OSvRMiMoI2qVX6lKYGmWAKArC1Np5kqgD+jOXK4rZ6UEHIOziLZzUu5il2XH+Q7oyuZPlqY6vk+WjyVdrwLdJT0DS7B/KzhbXMCglzG6zgeOw6IzWLZCbWPZvtWurbf8nO447S1qWeOzDWfL1dY0kOucqNVmpSqBpVUCXnKxa+GogCNPU7WQXVrltaTGxU/6q/TjloDjYI5Fdh0/mszoHFFwPMz7J44djpuRno74PplxzIQ6FejOhFlYQjzYMrD6jD+GMRVkq83ZHueUrHLHNdCp4CFvQ5B3iN9ffkY0L1M/VwksTRKw5R127ZofFVuaxlnHUiUwEyTgaEO5c7iofOX26/O4cXS64vuijmtJ1K9Ad0lIvfZZJVAlUCVQJVAlUCVQJVAlMO0SqEB32kVcO6gSqBKoEqgSqBKoEqgSqBJYEhKoQHdJSL32WSVQJVAlUCVQJVAlUCVQJTDtEqhAd9pFXDuoEqgSqBKoEqgSqBKoEqgSWBISqEB3SUi99lklUCVQJVAlUCVQJVAlUCUw7RKoQHfaRVw7qBKoEqgSqBKoEqgSqBKoElgSEqhAd0lIvfZZJVAlUCVQJVAlUCVQJVAlMO0SqEB32kVcO6gSqBKoEqgSqBKoEqgSqBJYEhKoQHdJSL32WSVQJVAlUCVQJVAlUCVQJTDtEljmge69997b7Lnnns0WW2zRfP7559Mu8Ml04H8X+/TTTzuret5FX3/9dev/qPJjGLcx+W97d9tttzQ/i+t/d/I/3Bx//PGN/2f8jDPO6BLtjLrv/1j330ROJfnffs4555xmu+22a3baaaeRm/6x6NbIAxqh4EMPPdQceOCBzVZbbdXcfffdgzVmso5N1f8eOJPHODhBi1jgkksuabbffvvmgAMOWMSWls7q9b9gX/R5ffCdec1P/vBFM+vcT5tb/vntojc4iRb+8fH85qMv50+i5oIqr83prj/326bpfrqgjaFPUw50BTn/N/Nyyy2X/v7973//gIeHH344PfvJT37SHHLIIc277777gzLTdQPw++qrrxZq/rbbbkv8LC4gtVDnPV9uvvnm5vDDD2+uvfba5oQTTmj22Wefxn/tGfTyyy8nkH7llVem67PPPhuP0vWZZ55JztaY22imjrvkNficqgBctt/1/Ve/+lVz0EEHdT2eUff9X+u77LLLtPB0yimnNOuss05n2z8mm+ocxBQ88F9NL7/88s0NN9wwcmszTcfmzZvX0KU2vz3yoIqCM22MBXtT8tV/e1762SOPPDIt0Kekg6WkETHNAuC00077gbyWkiFO2zC+nd80n8xdGPb9fc68ZtZBHzWXv/zNtPXb1vDJf53brHDDF81hj81tfn7bl83/dNXnzWP/mddWtPPeI+/Ma/6HKz5vfvvw3HS9682Fwfqf3vi2+S8u/LQBdheVphzoYshqbdNNN21WWmml5oILLvgBj7JlG2644RJxAieffHKzww47LMRTAO+ZBHSfeOKJBL5ff/31CV4PPfTQhWR29NFHN5deeml6Lrjuv//+E2U/+OCDVPa9996buFd+mInjLnn0/dFHH02yWNxAd++99/7RAN3333+/kYGdDrKQ6gO6Pxabmg7ZlG2ussoqYwHdmaZjjzzySLK1888/vxzapL/PtDFOeiA9FTfffPOmlNmZZ565kL/uqb7UP5o7d25K2myzzTaNTH+l8SVw2d++SaB2XoZ1P/t6/mIHurLHwPWzHyxgJGWWz/pkrEH99OYvmj0f/C7peMxTc5v/7ZoFO+r//nR+M+usT5rXP1nQx1iNF4WnBejqw3EAwOsXv/jFQl1a9R588MFpK3Trrbde6Nni+DKTgvJ9993XnHTSSQ2H+NRTTzXffLNgVXbhhRemhUJ+ZMFqeOWVV544imCxIIOOAOM11lgjfZ737bdJvk8++WSvSCvQ7RVPsywE6H4JfPe0At1RpPRdmcUNdOfPn5pAECM86qijUsZto402iluLfF0W7KgC3X41Oe+889JOrx2DLrJYv+KKK5qzzz67Ef/y3cuuOqPen3X07Cb9HfJRAmmj1ptJ5WYK0N3lz181sw6ds9CRhROfmtvMOmz2QscMvvx2fnP8U3Obre7+svnNPV811/x9Ab4h11mnfNxc/OJ3u80JPB89O4n7m3lN899c8llz6+tTkMr9fgKnDejutddeCYQ5wvDCCy9M6AtgdtdddzU777xzUwJdIPjUU09N55qcdzviiCOayLIChbY81Hn5b39L5yfzDOZVV13V7LrrrumcnEynraSSnJdab731mtVXXz31z6mjAHz4tP275pprNrZr862oF198sTnssMMadTbbbLOJc3jK4NX2HCDqTOdaa63V7LHHHr2GCnA7loDPf/3rXykDbnzhCB577LGUWQGCg7S97777xtfG1lie0Y1nZMhhDNGijPv6P/yh2XbbbdPf1Vdfnbq6609/Suc5LXJkXx0ROfbYY9ORC/rQx9Mdd9yR5jvOOb766qsT7OcZ3TZ5A/xP//Wvqfzbb7/d6Mtiij6oG9Q1h/GcXuHdvNAt7cbRBTplvM6qPv3006mK+XDPYiUnGYvdd9896Ymy++23X5rfOXPmpGJ9uvr8888n/acfxxxzTPPWW29NNE0n2A09xGc8I2N24Wr73JlmfCnjWJCz58bj3vXXX5/a62rLQ3PnSBEZGr8sTFdGd1Fsahz9uOaaa9L4fvOb3yRbM1dBk5Wnc81AmJ0S18g0telYl007Q87W+QW6u+KKK/ZmdPt0zHj03eUDPf/nP//Z0A192Zn67W9/G2Jo+nR/olDPB76WTzEmfvu5555bqLSMpXPrjlGxV/bBF5Znw4fGuFCjTZPa4utL2wd02LA+HMtiRxbzrnT6uuuuS/6c7of9R9vKhy2zh/vvvz89GnVuJQv4K1vszl2ze7Zn7DmxL/ctcDbYYIP02VEzFBldWXJjWH/99X+gG118Rh+T1fuo79qnU+HXdtxxx+bxxx9P1cKv8Seoy1fk+nD77bc3G2+8cWNRXBIASz7sVHxrO5/L36+77roTiRu2yOdEIqdsc9zvloNXvfxNM2sA6J73wtfNrAs+bWad/2njDOys8z5tZh07J4G1z7+Z36xz65cJMMs+lmdLT3hqbvM/X/l5Ojfr+vT73y1CXf/7yz9rZp35SXPnG982/9vVnzezTv+4efeL757f+Nq3zX99yWfpKMD/eMVnzcuzf7gY2O0vX6X6MqnK+pv91fwmMrpnPPt1s9rN3/H2v1/zefNK1sbbn81PzxLvZ3+SxhLy2x5oPe/T5me3fdloA/icdcYnDVDdRfjFx6/uXnAEdPO7v2r+l6sXZGQ/njs/tbvxH79sZKDx6izxpnctqPP/3vRFs9eDc1M3xz41d6L+L+/8stn/oe/ud/Ew7v1pBboC2WqrrZYcczDGkdrGaAO6HF0AC+U5dIEtshYMcIUVVkgOyIpPdpPhcCxemNEu4pTULYnBAw/a9DkyqAH4tC+D+uCDDyZHL0uK9CHQhRHfcsstqW9tIODCMQ0Ol3P+8MMPGyv8fCyp4Pf/ALfO8uVbXQCr4PLKK6+kUvrkaNwjK8GH08+dBBDFCQOYwIgXYiwiRn0BYlHHLbgLPCFHjBtzgFTgSqBAttUtMNocF5kZZ4AXQNNYg3Kg617Ie8stt0zghB6QH7k6X2j+kOBnXjjaoTl8880305zGnONXWzGHdNlxHOekg7Rpzr/44ou4NXE1F8YERMVnILZPVzl3AeGNN95I7egrgg3ZkLUjKciRIDrGNvBh/gV2pL5xh766p232hfraMpfO2J977rmprH/w0QV0J2tT2h1VP5xFByaRsRpnAJvJyhPgl7WMF1CBNfrJdlHoWJ9Ns8m11157Yk5eeumlNOddZ3SHdEy/Qz7QXETiABCNlwT7dD8NaIR/8M12tEUPS1BH1wBdgETiQTm2SjfCR48yxpyVIdsPXwwI831syLsd5oWt8td8BR8ZpE38uyIv64pD/B0aZW6NDzhFdETc8ZJh7us8M242AOQBtj4Dych3OiUBwn+cddZZyS7JDQ3xOVm9T41n//TpVPg1C72gmGd+rc9XRDnz/8ADD6QFtgV+Sb///e9TvBMLvPQt9vExsbBU3hyGLkd98tN2xPW4P5nrba9/28w6Ynaz7X1fDWZ0gbFZh36UwBYgue9Dc1Od//Paz5sXZ89r/jZ7XjPr+DnNrg8sAG1HPzE3gUQAD6Xs61Gzm7c+++77ta9+d+zAFn98vu/f874D00fPbmzVo50AzzM/WSgz6r4sZwKiB33U6EO2FAXQddYVX8GbcQYBvoA3Ss8P+qiJ87BaSWD+hDnNpX/7uvn06/nNSjd+kcb3XQ/RyoKrrte85cskE4D7zGe/bv67Sz9rAOqgPQDzI2an9uLedd/LIM7yAv3A7zkWF+d9muRy/otfN+Q81TStQBezHKFgwPg5mACgJdAFjgCD/IUq51PdCwfFEXNgiIFyqgIfY7BaDLrooos6A/Nxxx2XMltR1jUAH/6CrL7zYA/UBW+MH1/6D8JDnrEERNzrInIJUKWMsWlTsAziCAQW7Xgmk1eCKk7zL3/5S3LGygsAxnHrrbcm8Jef8Y1247qo45Zdwpf+kSBgIYHw4hlZBcmWyxSWBIgClHEGF0ADMoNKoOs+mQSQkxEXYGSZgUE6gVxl2G688cb0vW8OzYesfE7llitwpL3QE2PTZxuF7iojIJDHkK6STQ7wLVgCvApWOcj2oiH5GjsCqKOs7+RpIRAkMyRwor62LOIEohij8uTcp8uTsalx9EP/Fhmxu8NXAPOTlad6Fq7ejA8yR8YouAb53mfTdodkX3Pin7qA7pCOjeIDy/Zj4Tik+zmPXZ/pVwA5vNrZYlc50Ss7JkEWHLkeDo0x6sV1yPbDjnK/CCzpJwiQtLAL3i2gcjtSDt98adDQ3LI9meMgi0q7fF20ySabpF8oyZ/TpQDL7su4k1X4/SE+J6v3OQ+j6BR95ddixwkAj126Pl+hH3INkGwhkidigg8xf9VVV234LIQnRxrZdOACchG3c4rz4jLK45BsY/r7fivcS0+zjprdyC4iWd0+cl4USAv68MvvzsGe8tcFWUaAFchDgN+sY+Y0wF1QApDHz5nIejrPiidZYs/++t53oPD/uO6L9CJX1PMClnLPf7gANMazlG0+6KOFXs4KoAukBv2/gOqFC36l6aBH5jYypkHG5qhB0P/9+y+a//F7IOzeHRYFB33UPJedwY2ycZWlTgD5+DmprBfT8hflnK8N+UQdGXHtbnPvgt12/F/98jcJvKc2z/2kef/L+c0pz3ydzu/m54Cjnclcpx3oBkgBqgCFyEaUQBcwo+x5gBWQBN18e57B5BTOwzYe5+RPduunP/3pBODJy/cF5QikygM9pWMDKm3ZARR4BfSCOM4cbNtWC1AeZbqugizHzClGNsAKWZaA47CiJQPO3EqYXEqSgQRy8WQ7CuAke4BRMGmjALqLMm7Z8dg+tT0Y2dT4pQTOPObFePIgWfIkm3bZZZcl2VscBYUOBRB2v5S3e7I+gnP05wrQhP4o0zWHsnt58FS2BLrkZA5kKJBAlssu3fz+nwjQOdAf0lVbn17UDKIXoQ8WXuY+xiZzAtRH8CiBbsjstddeS83FAsSXvraU+/nPfx4spOuiAN1cPrlNjaMfdFrWXnZMljF+EWCy8nTkg/3K2uUko2SbO6jUsdymveSpDX4rpxKI5s+GdGwUH0jH9Qu0kWH4glF0P+el/Ewm9MtZSn+OQekntvyjPGCTL7jol3Kh50NjjHbKa5ftt9mReTLeIPOIh3gZk/5GQiXKSDzEDqB7fXPr+cUXX5zsxGc2aF4DoEab+bUL6ALlQfjLZTrE52T1PvpzHUWnZMXJhn4jmeeIw32+QtlSH1IDxT/8aC4HjyWRyMKuhAW4z+UCkU65H4mKotnOr8BU0EsfzUtHD/6/W7/8QZY0ypTXBHSPXAB0Pdfm755ZAA43vOPLZtbZ37185aiB5znYVOe/vOSzlOn0OYDuw+8sHLtnnTQnHWOQxfSXjjgcPbsBeEvqA7r5ry7Y+gc0c3J84OZ/ftt4AQzQl4EOAnT9akLQU+99B3SB8jZKR0BO/ThlcL/4Zn46fuDM7qzTPk6gX51ZR85u/u8//HCn0/3/p+X+nHS04ZPmmffnJx7XuOXL5g+vfdPMOmFO88YUvJA27UCXkwBaZKtk1IJKoGsLklLnL18BJfkWP0MogW5khqz+RqHJAl2ZHat6ACIMMLJp+h1ynF28AW/Ai+xBvhrmzMutQyCbjALc5G0Cm/fcc0/KwAAEkW0WlLqyjqMA3aFxO+cI/HGMzqpG4AUG8WoOh0gZR1qAfc7dgoE8gwK0DQFd2QeLnC7qGwtAXMq7BLradW7VgkL2qFwI5f22BeghXQXSS7AdbXqWZxvjflxLoGseZBxtF8rQ5MGiry0ylGnJaTqA7jj6gRdZILzRNUHZgnmy8iQPulmCOIs2ZxWD+mw6+mZzOfUB3SEdG8UH0jvZNhkxY7CQREO6n/PY9lk2zVEMvi3+jD/ajzolsCmB7tAYo524Dtl+mx0NAV0LwN/97nfRRbp6wUksiaxv39yqIJMtZpGrv9NPP32h9sovkwG6o/A5Gb3PeRtFp5Q3VsfgyCdfJPT5CvVKfcj7js98mvPcOUW21o6EOEqX86N8ysaxFUfzxqEAut7Yn3XinHT29ZvvE6TxrK+9cYHu4+9+B3SvKH7ia9b5nzb/1UXfZVY7ge4xs5tf3T0cI/E7WaB7/T++SWAauHXkAdj0kljQuEDXi2LO9OYUxyoCoDvaEWOPcgnMHvRRY9FR0v9+7efNhS9+nbLVsu+OiKD/9arPmyMzUF7WG/X7tAHdfJvoxBNPTAFKti6oBLoBZgSQIKs9BuDcKWoDugGGY6sl6sqGtRGgK6DlNAT4go94CUmQxVdklrQ15Djz/uKz1bZA4iWpkgSM3OF4LjObyyPqAHBxFjYCT2w5agOIbKOpGLcsheDuPxbcwcuLAAAgAElEQVTIt55iPmO7PPpvmxcZVwuYAMWXX375SEcX8gy69mUiOObIguZ9Ds1hvLwVdVzbgK4XNsyBjE+8vJHXic9tAXpIVwHoMpsa7dHZPJPmfi7LEuh6LjNHLwXpyHa539eWebRQCkCg/ChAd1ybGkc/wv7x4oyygMyOJytPx56AnnwBaGEAQGs3qM+m7bJog+3l1Ad0h3QsZNLlA501pcdBFpl0kUz6dD/K912Bx1yflOU7LCzyrHwJbMLfAMdoaIwlD0O232ZHQ0CXLgKnOcnUA6NBfXOrjAyzNoAw+jJE2jYHOVmY5pnMMqM7xOdk9T7nYUinoqxsNV2S+c4z132+Qt1SH6K9/MpXsq2ISZ7Z3dWfRSdbomf8bU7ew1GmjCF5mbbPwOx/Pp+fwB0AFSBX2ekAumlL/uCPFgJk6TjDYbMnQF0n0D3/0wTm8nG0I5cFQFcWNSiOLnRldBNvR81e6DyxX0fwO7hB4wJdxzRktHOSdSVbZ2yRbDZAvYDTpvnj90ciHKXIyYtn8ZJaZJPjDPJGd3zZeGltUWlagC5AZrUawMXZVkEh30L38pgMaRAj4CwYRRDwaus9HA3jYDAlCfK23cIhAzrx0ltZVtC3HROZR89j5RhnlNyzzRpAM7bGrK6RFSYDjJeG3JO1zoOelXQbr8riz/ECjlTfXlaKvzBqDsRLHzmfgD4gl2d+/YSYTGiAO+XVE4AQB3/nnXemz+U/UzFubcrUm19byUGAEjnLgMYYzH9bFpSc0znu749keAHLGILCWceLQu6X8nZPUEoA5vtfFnBPXYBgaA6dTXVOLbL0eJYdzs/pac9943IMowQGngfJTNORchHTp6uAlzo5mLGAodfOIeJPYAgSCGKeZWTybXdlzIcXd+haTn1thQzJI4jNAL9dNBmbGkc/LAbzc3peGrVtjyYrTzqbn+Nkd2SVy7fUsdKmZar4nTg3b1dDwHYeuo2GdGzIB+It/2UWL0k5KkAnY974yKDQfd+BrnixMZ7HFbCJ40dxz5Xfpo/5ePis/Cw4nnI9Hxpj3r7PQ7bfZkey7t5XCJIVxEPsBvoFAPEnXjSka/xmbM2rNzS3sv3AMJ8uIaGP3N9H33G1SM15cp9d5Gd08YfPODIzxOdk9T54ch3SqSjLl8m62inI/Vqfr1C31IdoL7/y22JhnCf3THwDoiM2iAvknS+wHc/Kd1jgAUdq8kV73k987gOzfc+ifnqx7LAFRxeANfVycLieowvZ8YBVb/qi+W8vXfA75ulc8EEfNbK9yJa8Nu5/e+HjAOlnuQ75aKH/bMHPdwF8Jfl5Lm04jhHkXKx7l3wPMN13FtkLbSgB/oM/Sv+xg++vfzy/mXXwR41fiAj6v5zRvWIB75GhLnmN8oCxIwX57/k6AzzrmNkTL6QB3vjKx+EXFvzE23vf/9qE9vyEGLl52Q5pU9tRz0tuZz+34Pzxd6XG/3fKga7tdVufDNqVw2A4kXH0qwIcAkCiDEOJ/8JWJoPDsIXBwfgcgYfTtTpWR1DPsx4CjGwXUMxhOwoQ9UqR2PZ3fhdvAhXwKnhoVxtApP4FLIbHOK049e0Mpb6ttIEd2QsB0tjUdw94lOnj4NzLj2sEL2TgWdsfeSCgW0DBJ1kIUvj01nGQc4IAeQ58PQOIOQTjsEoOhx/1XKdi3NGeuZFlKckc4N8fx+ZcXZuTIg+y5vjISwAgf7IA+jg7sjIWcgl5m48yk+/NaNlwZ7aBbBlKNDSHHK6gq645dowBPwBetBHjk/WQde4i80KX8GxBRyeC+nSVk7f7oU8AW6CLM+34i6wdOZBNLGDYHL4FkzLLDczFLxQED31tKaM9PNAh82HuzAdQWeqa8pOxKfVG1Q/90i/zQO5sIbJDk5Un0MFX0Eng0FznGbTQsT6bdpTGzhRwwJewN2/3++4YQEmj6FifDwRm+SF+Afi0PR+gSV9duu+ZHbYSxLjvhTxzbX6NIeyTTVs40WHP9GVxATDSMz7JuV66qIz5eeeddxJwGdWO9N9n+7K5YUd8nx00cyVrzofTUzpk7vBgPumDeGMu8MQXagPQj2TAKHNrrug+2fBNxixmxQIL7znxQ+QkJrBBcY/uqKM/ukKP8WlHAt9DfE5W73O+fO7TqbwsXeDbcurzFaEP5gNQjcRWXj8+O6ogkWV3yHzw6RH3lVGXvPkcMZcf5Lvy/zWV//PrF95f6SMAq++vr65fIwBg1feClf+VK51r1eYZnzT3vTXvu19UOHFOOuvq92GR/wZX9tiLYAkon/9pOmrgmSMUtuATTxd82lyd/Z6szK/fmQUA/TIC0NkF7NLWv5//OmFO+okyYBCPqd2zvuNNRjWd+z3ko8ZPh6H0CwlHz078HfDId/8xw6yTP078OXecfl/4sNnpVyb8KgM+tOmIgv95rSQ/iwYYk4dsrKMIXn7zKxI5HfLo3GbWuZ80Zz33daNf4DuONiiXjpac/ckEOI668VNj+MdL/JJFPJ/MdcqB7mSYKOtwmLkRlM+7vgN08fNLXWXcFyTz7GBf2fwZo2P4yLUt6Oflp+IzPmUFySNfaQ+1bYx9/yvaUP38+Sjjjqx7Xi8+G0ME0LhXXgWhPOMvSzZUp2wjvpMTHcqzA/FsaCz6DN3rkregKaAuCvXpKqefO/i8H2MytlEJwO8aR19b9AewoOdd9XMeJmtT2hjSj7A5th3Z07xvnycrT/MIQEYfZbujfKe3oQ+jyGoUHevygfgk69DRkj/9t+m+uVZvcdEoYwxeptL2o01XdgTk9fmmvHz+2WLCwip8CLl6BwGo7iJj9jcudfEZOjlZvS/56NKpKGceYrxxL659viLKDF2Nx+Il9/NlHbKwCMh3V/MyFg2x65nfzz8n4NcDdvOyU/1ZttKvFQCw49DX85pWUNnWhmMC47bvp7/yIw9+zizPyLb1M3TPT6fJ+v7j4/kLHVHI6+FTNrv8zeG8TNtnC4yp+l/RtD8jgW7bwOu9KoElKQGZVdllgTP/FYclyVPtu0qgSmDqJQCIybyy9Zzs8tjtqbRkJGBnzG7BEMGYfX9D9evzpU8CFegufXNaRzQNEnDcxhad6+LMjE3DUGqTVQJVAgMS8DN7joU5+uAYCpBrW31x7OINsLbMPjYPo+yYLLMCqgPvlEAFup2iqQ+qBBZIwNad84Rxzm/Bk/qpSqBKYGmUgOMo8VNr+U8bLo1jrWOqEliaJVCB7tI8u3VsVQJVAlUCVQJVAlUCVQLLsAQq0F2GJ78OvUqgSqBKoEqgSqBKoEpgaZZABbpL8+zWsVUJVAlUCVQJVAlUCVQJLMMSqEB3GZ78OvQqgSqBKoEqgSqBKoEqgaVZAhXoLs2zW8dWJVAlUCVQJVAlUCVQJbAMS6AC3WV48uvQqwSqBKoEqgSqBKoEqgSWZglUoLs0z24dW5VAlUCVQJVAlUCVQJXAMiyBCnSX4cmvQ68SqBKoEqgSqBKoEqgSWJolUIHu0jy7dWxVAlUCVQJVAlUCVQJVAsuwBKYc6F5yySXN9ttv3xxwwAEzWqz+O8Hddtut2WKLLZpPPvlksfDq/04//vjjm1//+tfNGWecMVafX3zxRXPaaac122yzTXPOOeeMVXdpKvz0X//aXHbZZTNmSA899FBz4IEHNltttVVz9913zxi+ppORe++9t9lzzz2T7Xz++eedXeVztSTsrZOx+mAkCfxYfPlIg2kpxB/7r3354/vuu6+lxGi3ZsJ/C/zwww83++yzT/pvi997773RGJ9hpfzvk+eee27z6quvzgjOPvvssxRrt9tuu2annXbq5Cn3c52FlvCD6667LuEd/63122+/nbiZbnnDKWS3//77L+HRN82UA10jOvLII1MQXOKjG2Dgtttua5Zbbrlmcf/3jpTtoIMOGuCu/fGWW27ZHH744e0Pl4G7hx12WLPKKqs0c+fOnbLRfvnll83XX3896fb8V6HLL798c8MNNwy2oa9DDjmkWWGFFZLuGU8bCb50c+ONN25uvPHGtiJL9F7YTt8isZyrqLO47W2JCupH1HmbHfxYfPm4Yp4zZ06z/vrrN59++mlz/vnnN/fff/+4TTQ333xzSupIQCyK/xi7444KDzzwQPIZb731VkeJmX37zTffTPxfdNFFU8bovHnzGkmiRaFTTjmlWWeddTqbKP1cZ8El/ICOiymvv/564mQy8iZLMh2V+I9NN9101OLTVm5agO6ZZ575owC6jz76aJr4xR14995770kD3d13332ZBrqC8X/+858pNYjNN988BbtFaRT4HgXoRh8nnnhiY9GinjHlJCOzySabJN2cTADO25quzzJInGYf0C3naknZ23TJYGlrt80Ofiy+fNy5uPLKK1O2adx6yltkSzbYXZMVnin0/PPPJ5v8sQJdcnzjjTfGAlJDsv/zn/+cZDIOOCvbpCt9QLf0c2X9mfL95b/9LckigC6+xpG3hA6f/+STT448pLPPPrsC3ZGlNU0Fl1TgrUB3miZ0ks22BfhxmxoX6F511VXN1VdfnRyHbf2cLr/88sYfp/Lggw/mj2bM51GAbsnskrK3ko/6vV0CbXawtALdY489ttljjz3aBZHdvf3225s777wzu9M05513XvPLX/5ySgHZQh1M8svSAHQnOfTOaosD6HZ2PsMetAHdcVisQLeQVjjHRx55pNlss83SFlGZ7brjjjuaI444YuJ8Y5zLOfroo5ttt9222XHHHZvHH388tew8q3ucE3LGZK+99moOPvjgZtddd20E0Daa9+23jRXFcccdl87YOMOUnynMA6+tJ2ctHSuwJaXPDTfcsHH+Bj322GPNzjvv3NimcD4xXzV3jSV4omDqyAI4r6JdRxcojnPCxua5s15Wh/h079Zbb40mJq4yusZ9+umnJ7nit1xh9fHaJTtn1Jyt3nrrrZuXXnopnVtbd911m4svvnii7/zD9X/4Q+IRnwAbuutPf0pZEmORJR+SC8fsLPfJJ5/cHHPMMQvJ9JZbbklb/J75c54o+jRm1DZna621Vgpgf//731MZ/5DrUUcdlc5Hq3vWWWelZ++++26aUyB1gw02SJ+vvfba9Oyrr75K+mYu6NoVV1wx0Z4PwKlAqV16s+KKK46V0SUz82Zbh7xyIhN61wZ0++a2S96yTvSGLWp3v/32a9Zcc800PjqIyMRcsrvQJ1tbYbeyzPvuu28qc8899zQBdF944YVml112Se3Z4ost3HKu9JHbW4y3bzxRpu3aZ9tDujyq3d10002NRemhhx6arpG9G1Xv8G1Bw0fREb6tzN7nY+vqb1F8Yi4LfsiRmPLMXJ8djOLLR53DvjkjB/ydeuqpyU96d4K/9UfHUJfvymUYn5999tkJn8vHxs6Ic5falKFjAz47v1gS/8onsMF8K/39999PuzDmVZZwlPO5H374YfJlzgOLYU888UTqrk8e+byN4o81GED3rrvuanbYYYdm7bXXTv3lx7yG/Jp2yE4s4A/imNwrr7yS7Fy7nnf51JBj15iDz9Lvf/DBBxNnjMVf5EgJfSU3vk3c5MOA15y6/J7Ypbw5DF0y/6hPZ8Uux8vEWXFa5r4ro1v6OXGHveuXnPjaNdZYI11hD7q23nrrpTgb2IJPYJNik5gPe9BN/WqPP2YPjtqEP47x99kE/XJOVgxzhEAMIAsZ3TZ5a/PFF19MGEdcM4Z470TCxdyrj0/yDP77eMgzuuHHnHd++umn0xAC25100kkxpGm5TtvRhdVXXz0BRoYliK600koTTt7kERjHhhgT5UDKC/4CTND8+fOTwjsfQil+9rOfTWS6CFvbHFBJlGrVVVdNt006Q8mBWxl4BUBt2VIW1FZeeeXm0ksvTXxSVsqBLrjggkb2A199Y1EWWNBOODeGhv84o2vrQJ8cZxDHJri2EWUlH8FJwHU+DMjSDyLTLl6HZEfpnB2VTVSWUYVhtPHCGPQFhAYZl0XLkFzIF8A0fsQYYyED5Nq6pwtotdVWSyDaZwGIPINizgB0fXKw5ibkqxwnEvpkgWJMnJD5I0PnYAV0n+lJ1CFbZM7os5ewEJAriIQ+mC9tlk4oFe74J4CuF36c742AaQxetuMItJlndPvmdkjeXprTHlvDt8BF74Jn4ybDsENs03MLh5AJXkImAXTpjHOOnmk/9Fz9cq5Ke+sbT4fYJm4P2faQLg/Z3fXXX99stNFGEwtjwZUO0C80it5ZNHkZI4AGB2/R1kZ9/S2qTwxZWKxdeOGFyX7+/3buJkeKJIkC8NH6AGzZcAMOwSk4BhK34AKIFQuEkJA4RYvRl9OPsfEOd48sqgJm2kyqisz48Z/nZs/MzT2yLuWu7IBdrLj8njHcjRm7jY2pl6O92eSff97F+2wBtzgSW51wCJ0lyhTI4NJq87eLw7/4gJx++/btzV7ZgeCD7cruZhKU+3JUPj1K0oLdux/mOzwybmf5OIEuvTO+2iRAE+REcGFseOS13OMok82ejW9E3wV2ZMaprq36vOJ9ASa89ZfASKAr4SLw55Pwlz65Rla8h7cEwfohyNQustJZvsz4eCkuwjfNAl33jDwXLpTIw+vG5Y8//rjxK37El3hWsBixZxwfJ7iVWKA3gluYJ46CT/qx8+cCS5wT/rYFAxbZujDiTWfwvfsIP5z6XPv27dvteb5EG5zbtaEGuuExeEaUwe/87D7qlDc7Plmgi6QiIv7qBAWlAhFAE8AK/iKcr+DN4BPKnKyhQRdcAYg4uvfohR0zIbPBSLIz+T46XucpdAb6y5cvtwGVjamD8/79+1t/XN/1hWEmA5l6taMGYsr31n5E0DsjTuSszAiF44jSz1Vbd9jBnYOIKJvxHWHrns+fP99wePfu3e0RJCEzS3a4cDQ1qDLDl9lEYIhNAEiMrzEJ+ZntM74qrteMqwCtEhPdSUDnOX1M1sB3QXX9JQvYjyQva6nNxIxc5qmKMmsd9drR5wS6CbyTSebgYXcU6K7Gdoc3ctMnQXkECcpaRDhe2GaCIYOLoLOygvzZMkmgK+CLIOXqIMaxGu1t1Z+UOTvubPuMLs/sjs4h/OigNtBLOmV8Iiu9i97WCazMYNXLlHOmPv15KCdWLIxtJsWpP8fRDpxPwJl7Ri6/ZwxXYxabS3CaJefo14670j5HQVjlFufwraApIpEwrqTkWj2yh3CP8yYqkid8ADGpf/HixS35kICiPo87cUOCE8Ei/+X7Cg9l1HHzfcfHCXTraqO2C8bxQzCuwWvltdpuOomHExS7JhGRfqw4ddXnGe+nbnFAxdu4WVGLSGzhMb6X7HhPdtv9abdnVjoruINX9M794oEju7014Pv3my+pPumIa02Kqt/GjXxrTRKNfGKFjK+JsFt9STZ0ZRPu1Q8vKEaMu+cT6Do/4i15JQlEcn/4QoLE87LhkVUb3FMDXd9N6PFY8FWHMp5anizQNbARM0cAZfko52VHZK8EgjJkETMeisPxEhnhAGOWJK0vMMofpyTzOhMGzmAZLiWPjI7XecpWnZNzHLhZXuoTbCKrkJ17Zn0xK6sK7t4x0LXFAz7aSZLZvH0Z/o2BrstmhpaayKqtO+wQazLgqRbWCfxzrh5fvnz5/dWrV7dTxrdmIJ2c4WJZ0E+tRQQSHIXsCywQVKSSwRg8uWccM3pTA3b3IG6zabNcWWsGFxkdfH4dgMPMmMv6Ilwzbe1LhiZlPDTQ9TwMOWFtRDTkKNBdjW3aMcM75FudnHFLptvzCWLYhYynyYPlquhjzQwl0K0vo7Fj5BwZx2q0tzP9SVmz48y2z+jyzO4y+ZBFqmJiYBk8stK7YAm/6JBszbNnz27jnDIcz9T3M5x4hEWtP59HO3BeoLvi8oeM4dGYyXqzIbZH+IVwmu877ro99Ne/58+f/y1zbvLLp7Ax8tBAl71UPJRlcocTPnz48FcL/nNgz3xHlcpnzh/h4fzRuK34+CjQDY/gvhWv1fblM79Lx/GypFRNBrhnxqmrPs94P3WOgZdAtyaZPn36dMO68phnZ7x3FOiudFaihv5UeWigW9uIO+hwBLfQGbFRZOQTOqutETGQZ5JYWtmEVUf31kmP9ji3CnTVJUsrE24i6H4JLXIU6K7a4Jkx0OUvBPgmfAS/VB9yO/kE/35JoAtITtOsG4gCS4NcxcxPlhMp1OU+gSqHcUYMMrL0hrvsMINJVs7zo+N1blQ255BLzeTUund9EYQn25rnxkAXYTBwS1UyGrMMquePAl04JiO8ausOu3uJVXvsK6S4jNA+mziSHS7aOU4AlIf0GZeZ9ZGMwZN7xjEbA13jzNkIbs3sBfMV49HBM0Jt0IdR0j77VKs8JNDNMr9srvq0OwF+HFSdOKzGdof3mUBXf6zEvH79+pZVMekw21avz/bMRR4j0F31J/XMjjvbPqPLM7vLcug4MTchsYc5stK76Ilgeidn63soJx5hcdSm0Q7cswt07xnD3ZjJJJlQSkqYVGZvrnbsuKv2RxKCDlfhcGW4EmQ+NNDFWbbAVcmEKdsu6jWTnJqVq9d2eByNG7xniYejQBemuAWPrHittiufTXbzLI5KFtX1Faeu+qz9R7yfOu8NdHe8dxTornSWno0/h/W7Brorm8Dbxi5b4uB7JtCV3LCFyPjm/oz7UaC7aoM6x0DXOauItiuwxZoYiQ48xfGXBLqyr4KNBBKWKih4FUGAgTKrSUDgulkmRR2XibJnp5Yh+1L3wZih1N8tPRvocnB1VqmO1Lfriw394xLZGOgqz2xZxlmwW2d5tT8+61MlCu2wdSH9WrV1h929xKo92irI0/760sYOF4o+zpyVZ6+O5fLax4rBvYGupVr41LbJ7CgnwsHDJhK9+PjxY07djrCWeeIw61YJF+8NdC2LZ1lJNkKZAqfs5zwKdFdju8P7bKArmIVHJpeIUttkf2uG4jEC3VV//gv4gy872z6ry0d2R2f0uS6pCYrpjRdbI6tAF7cpI1uu8kx4I98dz9b3UE48wqLWn8+jHTi/C3TvGcPdmLEpE1yrJqPsuKver00ccBV161/koYGudzzoQV0KN4Hmq7LtInU44nPXBLWj7PA4GrdVoCujPNZF/5z7+vXrj8TOEa+Nbct3iSK2X5M1O05d9XnG+6nv3kB3x3sJdMOr6lnpLD7gLzIhcv/vGuiubML2AuPu5bJIAtdZRjf6w/cQE83oju8JdPnHyKoN7jkKdG2FUy5byra4lPdUxycJdG3Ornt0LbvpWJYCOVFbFZL9szQqyKvCIZg52/9UnYPZxc0JlaVnwMssjqINlp8JYjJTOcroCjQi2jUGMYIi+0oqkXmpwzLKri+I27OZFemzjLS3MasgIhhlj2u9Vj/L3tpXFZE5gkeCkFVbd9ghbCReRXZkxKNe99lSlTZk/6ZzO1wy46zLfTI5ljE4KUFEtqsIpjPzSwa5tmEcM9mH9MOSm7ZlG4yX+Hyv+oJ8M1FQLpKzZGTmGR21FyxtEITLEGcDvXbKar958+ZHswQIWfL/cbJ8cA1JR4xrDaJCBvUN49XY7vDOfupKfGxjnIR5SY0eBi/t0zYOqNphXrjIPnr32bqgHZFxrDKBiL2t+iPAlj1Nhjtl5riz7bO6PLM7Oq3fEcGBCVjlgJXeeU4wRU+yNGdybiWrOt2Uf6a+h3LiERaptx5HO3Btx+WrMaxl+7wbM1sC2JwAz9IrZx2sdtxV6/IiFd7KL+ywZwFi1Wl+QKC5EraPu3F9hO7ilpq9tbwueApX5F5HPoLeVFu3hUAwscPjaNxWfJzsbd1SJ7tqrz3Z8Vptdz7DEl86RnacuurziveVLw6oWPEFlaOy+hEe2/FeeIqdR1Y6Gz2rK4psQPA7k5HnjrgWl1UfQ3/wrNgoMvKJpENNAMLdM1lpSlvpSSSxkLEet00Gi8phFe9sp8gvYPkVJfXlhXF+mC7AL7Jqg3sEwuPWHXbCv1q9qT5FjOaXjGr5qednj48e6BoEnQCIwMXgIHyAWV4HMmcqkEEOCJ4RCRQodc1mynrJ6I7iJy9sCbD/DTGahR0JZfW2LUOxlCWrJ5PsdxANJuXTLhlWDlt7fbePqGZiDExmLu7Vzvy24q4vnmWM2otczYz1m+GM7baPpy7XHfVJmymD/a0cgrLSFvev2ur6DDskydHoP4MUxNmK4DuD8eLETDyr7VV2uDBEW0rgIOgXdKXvZo7KY+Swpj/GRz2CKW1yTuCQMaNzDFmAyIG4h24RQaW9kTITtsmoyxgniDLW9E/2Pfuz6aklLH8cmb1I0U06baXAJEzQi5Dome8cNLHndpykOa8f9IADVSf9ZeDaEvKm98hBH7QzRLYa2xXeZvBWJJQnQ4Oc6A5MZKJr5lIbjbcxiCAeYxVhO8pRnnIFpTDQH3aNxMex0r7R3lb9sVXC+PsFlCNZ2fa9unxkdxwQfTHu8Gdn0RftOaN3JkDw0Q9lCTqqk6n92tWXe+/lxIoFvTMxnsloB2e4fDWGYz2rMXOvSS/H649u8iGSE1n5m3HXWA/nSR+NKz03BiaeWQW0/xAfC1hxcyastRzL/LaJ0XEcZW9/loFtVZAdlulTNr2mrzPRb37HxIkOZH//Co86bmf5GEfjcLap316Y5SvwVWTFa7mnHmFvPDIGubbj1FmfV7wvE2xSgMsFbCYEJijGSbJFVpz/NSbG1q8ArHgPX/vDZ/rAzyhjp7P0w5ibIPAh+B+38TnRgeAw8pw20YlwrQAbf+BZOq1sY4BP3IMX8ETlE323esN/wMM1/WC/njEhTWJrZRPu4T/88gs/JWBnU/SVrx3xNqmkP55RF74zFvxiVgHob3RZO8msDfhDWepUVxVxXX3p0DV6Cqddsq+Wc/bzowe6ZytGOrJkEWSTQCLn3MMwjgSZUarZ9TxjVp+MinOUqi475b4zR3Wpc5QzfdG3kGGdxaQsxlczejk/HqPg+sRpzfoya6vyzmI31r36boT8arMAAAK/SURBVDlrlDO4IFBZ1iMxVjVjeHTP2XMIKtkhzwjkYB4xPqP+uSaDc3TeNfqrjWQcU/jPxub2wE/8m43tGbzPVFtxcj+c1PlUMusPbBHrTB7Dtnd2pw0mBlVXZu2ZndfOOnGY3ef8rr7H4MRV/TM7WD3j2mwMx+dmY+ZtehkeWEfwpSAnDtX5e7gLt+DII25KHQ890gcTyOrDVmW534rXqEczPFZlza7hOO1RhzflVzq34rWx/Bl+O06d9Vn5K94f6999P8N7waWWtdJZ3C1I1YeR22sZv8PnlU24ph/h9DN94Y+jp45jcD+OOwxWbTjCyMQzvrNeFweeaWN95sznXxbonmncP+EeywUI3fIHsm9pBBqB77dlZ1mFvPH72Ji03T02oj9XngzhOKnhcAW6Nfj9uVr66UagEfhVCMgiW7Ux+azbU65oTwe6V6C8qMOypqUVS20tjUAj8G8EkOJTBjhtd7+Xpvl9Z8ukltstmVrWtOXMEnhLI9AI/O8jYAuOrQmOT7XaOUOpA90ZMhedt1S42jd3UTO6mkbgH4VA293vN9yWSU1uZH0svV/tDH8/RLpFjcD/DwK2ithTbavJ1dKB7tWId32NQCPQCDQCjUAj0Ag0Apcg0IHuJTB3JY1AI9AINAKNQCPQCDQCVyPQge7ViHd9jUAj0Ag0Ao1AI9AINAKXINCB7iUwdyWNQCPQCDQCjUAj0Ag0Alcj0IHu1Yh3fY1AI9AINAKNQCPQCDQClyDQge4lMHcljUAj0Ag0Ao1AI9AINAJXI9CB7tWId32NQCPQCDQCjUAj0Ag0Apcg0IHuJTB3JY1AI9AINAKNQCPQCDQCVyPQge7ViHd9jUAj0Ag0Ao1AI9AINAKXINCB7iUwdyWNQCPQCDQCjUAj0Ag0Alcj0IHu1Yh3fY1AI9AINAKNQCPQCDQClyDwLxqgnO4Df9LNAAAAAElFTkSuQmCC) # https://www.dw.com/en/israels-clever-coronavirus-vaccination-strategy/a-56586888#:~:text=To%20date%2C%20Israel%20has%20administered,40%25%20of%20the%20country's%20population. # In[ ]: cad_start = "2020-12-01" # 13 confirmed cases cad = cad[cad.date >= cad_start].reset_index(drop = True) cad["days_since_start"] = np.arange(cad.shape[0]) + 1 # ## Data for Regression # In[ ]: # variable for data to easily swap it out: country_ = "Israel" reg_data = cad.copy() # In[ ]: reg_data.head() # In[ ]: reg_data.shape # ## Change Point Estimation in Pyro # In[ ]: get_ipython().system('pip install pyro-ppl') get_ipython().system('pip install numpyro') # In[ ]: import torch import pyro import pyro.distributions as dist from torch import nn from pyro.nn import PyroModule, PyroSample from pyro.infer import MCMC, NUTS, HMC from pyro.infer.autoguide import AutoGuide, AutoDiagonalNormal from pyro.infer import SVI, Trace_ELBO from pyro.infer import Predictive # # First method # In[ ]: # we should be able to have an empirical estimate for the mean of the prior for the 2nd regression bias term # this will be something like b = log(max(daily_confirmed)) # might be able to have 1 regression model but change the data so that we have new terms for (tau < t) # like an interaction term class COVID_change(PyroModule): def __init__(self, in_features, out_features, b1_mu, b2_mu): super().__init__() self.linear1 = PyroModule[nn.Linear](in_features, out_features, bias = False) self.linear1.weight = PyroSample(dist.Normal(0.5, 0.25).expand([1, 1]).to_event(1)) self.linear1.bias = PyroSample(dist.Normal(b1_mu, 1.)) # could possibly have stronger priors for the 2nd regression line, because we wont have as much data self.linear2 = PyroModule[nn.Linear](in_features, out_features, bias = False) self.linear2.weight = PyroSample(dist.Normal(0., 0.25).expand([1, 1])) #.to_event(1)) self.linear2.bias = PyroSample(dist.Normal(b2_mu, b2_mu/4)) def forward(self, x, y=None): tau = pyro.sample("tau", dist.Beta(4, 3)) sigma = pyro.sample("sigma", dist.Uniform(0., 3.)) # fit lm's to data based on tau sep = int(np.ceil(tau.detach().numpy() * len(x))) mean1 = self.linear1(x[:sep]).squeeze(-1) mean2 = self.linear2(x[sep:]).squeeze(-1) mean = torch.cat((mean1, mean2)) obs = pyro.sample("obs", dist.StudentT(2, mean, sigma), obs=y) return mean # In[ ]: tensor_data = torch.tensor(reg_data[["confirmed", "days_since_start", "daily_confirmed"]].values, dtype=torch.float) x_data = tensor_data[:, 1].unsqueeze_(1) y_data = np.log(tensor_data[:, 0]) y_data_daily = np.log(tensor_data[:, 2]) # prior hyper params # take log of the average of the 1st quartile to get the prior mean for the bias of the 2nd regression line q1 = np.quantile(y_data, q = 0.25) bias_1_mean = np.mean(y_data.numpy()[y_data <= q1]) print("Prior mean for Bias 1: ", bias_1_mean) # take log of the average of the 4th quartile to get the prior mean for the bias of the 2nd regression line q4 = np.quantile(y_data, q = 0.75) bias_2_mean = np.mean(y_data.numpy()[y_data >= q4]) print("Prior mean for Bias 2: ", bias_2_mean) # In[ ]: model = COVID_change(1, 1, b1_mu = bias_1_mean, b2_mu = bias_2_mean) # need more than 400 samples/chain if we want to use a flat prior on b_2 and w_2 num_samples = 400 # mcmc nuts_kernel = NUTS(model) mcmc = MCMC(nuts_kernel, num_samples=num_samples, warmup_steps = 200, num_chains = 1) mcmc.run(x_data, y_data) # In[ ]: # Save the model: import dill # with open('israel_new.pkl', 'wb') as f: # dill.dump(mcmc, f) with open('israel_new.pkl', 'rb') as f: mcmc = dill.load(f) samples = mcmc.get_samples() # In[ ]: # extract individual posteriors weight_1_post = samples["linear1.weight"].detach().numpy() weight_2_post = samples["linear2.weight"].detach().numpy() bias_1_post = samples["linear1.bias"].detach().numpy() bias_2_post = samples["linear2.bias"].detach().numpy() tau_post = samples["tau"].detach().numpy() sigma_post = samples["sigma"].detach().numpy() # build likelihood distribution: tau_days = list(map(int, np.ceil(tau_post * len(x_data)))) mean_ = torch.zeros(len(tau_days), len(x_data)) obs_ = torch.zeros(len(tau_days), len(x_data)) for i in range(len(tau_days)) : mean_[i, :] = torch.cat((x_data[:tau_days[i]] * weight_1_post[i] + bias_1_post[i], x_data[tau_days[i]:] * weight_2_post[i] + bias_2_post[i])).reshape(len(x_data)) obs_[i, :] = dist.Normal(mean_[i, :], sigma_post[i]).sample() samples["_RETURN"] = mean_ samples["obs"] = obs_ pred_summary = summary(samples) mu = pred_summary["_RETURN"] # mean y = pred_summary["obs"] # samples from likelihood: mu + sigma y_shift = np.exp(y["mean"]) - np.exp(torch.cat((y["mean"][0:1], y["mean"][:-1]))) print(y_shift) predictions = pd.DataFrame({ "days_since_start": x_data[:, 0], "mu_mean": mu["mean"], # mean of likelihood "mu_perc_5": mu["5%"], "mu_perc_95": mu["95%"], "y_mean": y["mean"], # mean of likelihood + noise "y_perc_5": y["5%"], "y_perc_95": y["95%"], "true_confirmed": y_data, "true_daily_confirmed": y_data_daily, "y_daily_mean": y_shift }) w1_ = pred_summary["linear1.weight"] w2_ = pred_summary["linear2.weight"] b1_ = pred_summary["linear1.bias"] b2_ = pred_summary["linear2.bias"] tau_ = pred_summary["tau"] sigma_ = pred_summary["sigma"] ind = int(np.ceil(tau_["mean"] * len(x_data))) # In[ ]: mcmc.summary() diag = mcmc.diagnostics() # In[ ]: print(ind) print(reg_data.date[ind]) sns.distplot(weight_1_post, kde_kws = {"label": "Weight posterior before CP"}, # color = "red", norm_hist = True, kde = True) plt.axvline(x = w1_["mean"], linestyle = '--',label = "Mean weight before CP" , # color = "red" ) sns.distplot(weight_2_post, kde_kws = {"label": "Weight posterior after CP"}, color = "red", norm_hist = True, kde = True) plt.axvline(x = w2_["mean"], linestyle = '--',label = "Mean weight after CP" , color = "red") legend = plt.legend(loc='upper right') legend.get_frame().set_alpha(1) sns.set_style("ticks") plt.tight_layout() sns.despine() plt.savefig('/content/sample_data/israel_weights.pdf') # In[ ]: print(w1_["mean"]) print(w2_["mean"]) # In[ ]: 1 - 14/120 # In[ ]: start_date_ = str(reg_data.date[0]).split(' ')[0] change_date_ = str(reg_data.date[ind]).split(' ')[0] print("Date of change for {}: {}".format(country_, change_date_)) import seaborn as sns # plot data: fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 5)) ax = [ax] # log regression model ax[0].scatter(y = np.exp(y_data[:ind]), x = x_data[:ind], s = 15); ax[0].scatter(y = np.exp(y_data[ind:]), x = x_data[ind:], s = 15, color = "red"); ax[0].plot(predictions["days_since_start"], np.exp(predictions["y_mean"]), color = "green", label = "Fitted line by MCMC-NUTS model") ax[0].axvline(19, linestyle = '--', linewidth = 1.5, label = "Start of Vaccination: Dec 20, 2020" , color = "red") ax[0].axvline(ind, linestyle = '--', linewidth = 1.5, label = "Date of Change: Feb 17, 2021", color = "black") ax[0].fill_between(predictions["days_since_start"], np.exp(predictions["y_perc_5"]), np.exp(predictions["y_perc_95"]), alpha = 0.25, label = "90% CI of predictions", color = "teal"); ax[0].fill_betweenx([0, 1], tau_["5%"] * len(x_data), tau_["95%"] * len(x_data), alpha = 0.25, label = "90% CI of changing point", color = "lightcoral", transform=ax[0].get_xaxis_transform()); ax[0].set(ylabel = "Total Cases",) # xlabel = "Days since %s" % start_date_, # title = "Confirmed Cases in China") / ax[0].legend(loc = "lower right", fontsize=12.8) ax[0].set_ylim([100000,1000000]) ax[0].xaxis.get_label().set_fontsize(16) ax[0].yaxis.get_label().set_fontsize(16) ax[0].title.set_fontsize(20) ax[0].tick_params(labelsize=16) plt.xticks(ticks=[19,46,78,121], labels=["Dec 20", "Jan 15", "Feb 17", "Apr 1"], fontsize=15) ax[0].set_yscale('log') plt.setp(ax[0].get_xticklabels(), rotation=0, horizontalalignment='center') print(reg_data.columns) myFmt = mdates.DateFormatter('%m-%d') sns.set_style("ticks") sns.despine() plt.tight_layout() ax[0].figure.savefig('/content/sample_data/israel_cp.pdf')