12_monitor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "# hide\n",
    "# default_exp monitor\n",
    "from nbdev import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monitor\n",
    "\n",
    "> Monitoring functions for HECSS sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "from numpy import sqrt, loadtxt, array, linspace, histogram\n",
    "from numpy import median, abs, convolve, ones, arange, cumsum\n",
    "from IPython.display import clear_output\n",
    "import subprocess\n",
    "from time import sleep\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.pyplot import plot, figure, subplot, legend, show, sca, title\n",
    "from matplotlib.pyplot import hist, semilogx, semilogy, axvspan, axhspan\n",
    "from matplotlib.pyplot import xlabel, ylabel, xticks, xlim, ylim, axhline, axvline\n",
    "from scipy import stats\n",
    "import sys\n",
    "from spglib import find_primitive, get_symmetry_dataset\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "from ase.data import chemical_symbols\n",
    "from ase import units as un\n",
    "THz = 1e12 * un._hplanck * un.J # THz in eV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def plot_band_set(bnd, units=THz, lbl=None, **kwargs):\n",
    "    if lbl is None:\n",
    "        lbl=''\n",
    "    kwa = {k:v for k, v in kwargs.items() if k not in ('color',)}\n",
    "    plt=plot(bnd[0], un.invcm * bnd[1] / units, label=lbl, **kwargs)\n",
    "    for b in bnd[2:]:\n",
    "        plot(bnd[0], un.invcm * b / units, color=plt[0].get_color(), **kwa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def plot_bands(bnd, kpnts, units=THz, decorate=True, lbl=None, **kwargs):\n",
    "    plot_band_set(bnd, units, lbl, **kwargs)\n",
    "\n",
    "    lbls, pnts = kpnts\n",
    "\n",
    "    if decorate:\n",
    "        xticks(pnts, lbls)\n",
    "        xlim(min(pnts), max(pnts))\n",
    "        axhline(0,ls=':', lw=1, alpha=0.5)\n",
    "        for p in sorted(pnts)[1:-1]:\n",
    "            axvline(p, ls=':', lw=1, alpha=0.5)\n",
    "        xlabel('Wave vector')\n",
    "        ylabel('Frequency (THz)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def plot_bands_file(fn, units=THz, decorate=True, lbl=None, **kwargs):\n",
    "    bnd = loadtxt(fn).T\n",
    "\n",
    "    with open(fn) as f:\n",
    "        p_lbl = [l if l!='G' else '$\\\\Gamma$' for l in f.readline().split()[1:]]\n",
    "        p_pnt = [float(v) for v in f.readline().split()[1:]]\n",
    "    kpnts = (p_lbl, p_pnt)\n",
    "\n",
    "    if lbl is None:\n",
    "        lbl=fn\n",
    "\n",
    "    plot_bands(bnd, kpnts, units, decorate, lbl, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def run_alamode(d='phon', prefix='cryst', kpath='cryst', dfset='DFSET', sc='../sc/CONTCAR',\n",
    "                o=1, n=0, c2=10, born=None, charge=None):\n",
    "    fit_cmd = f'/home/jochym/Projects/alamode-tools/devel/make-gen.py opt -p {prefix} -n {sc} -f {dfset} -o {o} --c2 {c2} -d {n}'.split()\n",
    "    b = ''\n",
    "    if charge is None:\n",
    "        charge = prefix\n",
    "    if born is not None:\n",
    "        b = f'-b {born} -c {charge}'\n",
    "    phon_cmd = f'/home/jochym/Projects/alamode-tools/devel/make-gen.py phon -p {prefix} -n {sc} {b} -k {kpath}.path'.split()\n",
    "    alm_cmd = f'/home/jochym/public/bin/alm {prefix}_fit.in'.split()\n",
    "    anph_cmd = f'/home/jochym/public/bin/anphon {prefix}_phon.in'.split()\n",
    "\n",
    "    with open(f'{d}/{prefix}_fit.in', 'w') as ff:\n",
    "        fit = subprocess.run(fit_cmd, cwd=d, stdout=ff, stderr=subprocess.PIPE)\n",
    "\n",
    "    with open(f'{d}/{prefix}_phon.in', 'w') as ff:\n",
    "        phon = subprocess.run(phon_cmd, cwd=d, stdout=ff, stderr=subprocess.PIPE)\n",
    "\n",
    "    alm = subprocess.run(alm_cmd, cwd=d, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "    anph = subprocess.run(anph_cmd, cwd=d, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "\n",
    "    for p, l in zip((fit, phon, alm, anph), ('fit', 'phon', 'alm', 'anphon')):\n",
    "        if p.stdout is not None:\n",
    "            with open(f'{d}/{prefix}_{l}.log', 'wt') as lf:\n",
    "                lf.write(p.stdout.decode())\n",
    "        if p.stderr is not None and len(p.stderr) > 0:\n",
    "            with open(f'{d}/{prefix}_{l}.err', 'wt') as lf:\n",
    "                lf.write(p.stderr.decode())\n",
    "\n",
    "    return all([r.returncode==0 for r in  (fit, phon, alm, anph)]), fit, phon, alm, anph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def get_dfset_len(fn='phon/DFSET'):\n",
    "    try :\n",
    "        with open(fn) as dfset:\n",
    "            return len([l for l in dfset if 'set:' in l])\n",
    "    except FileNotFoundError:\n",
    "            return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def show_dc_conv(bl, kpnts, max_plots=4):\n",
    "    prev_n = sorted(bl.keys())[-1]\n",
    "    plot_bands(bl[prev_n], kpnts, lbl=f'{prev_n}', color='C3')\n",
    "    y_lims = ylim()\n",
    "    alpha = 1\n",
    "    plotted = 1\n",
    "    for n in reversed(sorted(bl.keys())):\n",
    "        if n > prev_n*0.75 :\n",
    "            continue\n",
    "        alpha *= 0.8\n",
    "        plot_band_set(bl[n], lbl=f'{n}', alpha=alpha, color='C0', ls='--', lw=1)\n",
    "        plotted += 1\n",
    "        prev_n = n\n",
    "        if plotted > max_plots :\n",
    "            break\n",
    "    ylim(y_lims)\n",
    "    legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def build_bnd_lst(directory='phon', dfset='DFSET', prefix='cryst', kpath='crast', sc='../sc/CONTCAR',\n",
    "                  order=1, cutoff=10, born=None, charge=None, verbose=False):\n",
    "    N = get_dfset_len(f'{directory}/{dfset}')\n",
    "    bl = {}\n",
    "    for n in range(1,N+1):\n",
    "        if verbose :\n",
    "            print(f'Using first {n:3} samples', end='\\n')\n",
    "        run_alamode(d=directory, prefix=prefix, dfset=dfset, kpath=kpath, sc=sc,\n",
    "                    o=order, n=n, c2=cutoff, born=born, charge=charge)\n",
    "        bl[n]=loadtxt(f'{directory}/{prefix}.bands').T\n",
    "    if verbose :\n",
    "        print()\n",
    "    return bl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def build_omega(bl, kpnts):\n",
    "    omega={}\n",
    "    eps=1e-3\n",
    "    kp_dict = {k:v for k,v in zip(*kpnts)}\n",
    "    for k,v in kp_dict.items():\n",
    "        omega[k] = array([[n] + list(bnd[1:,abs(bnd[0]-v)<eps][:,0]) for n, bnd in sorted(bl.items())]).T\n",
    "        omega[k][1:] = omega[k][1:,-1][:,None] - omega[k][1:]\n",
    "    return omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def plot_omega(omega):\n",
    "    for k, o in omega.items():\n",
    "        if len(o[0])<2 :\n",
    "            return\n",
    "        l = k\n",
    "        if k == 'G':\n",
    "            l = '$\\\\Gamma$'\n",
    "        semilogy(o[0, :-1], (un.invcm * o[1:,:-1].std(axis=0))/THz, '-', label=l)\n",
    "\n",
    "    legend()\n",
    "    plt.gca().set_xscale('function', functions=(lambda x: x**(1/2), lambda x: x**2))\n",
    "    rng = 10*un.invcm*median([o[1:].std(axis=0) for o in omega.values()])/THz\n",
    "    if rng > 1e-3:\n",
    "        ylim(None, rng)\n",
    "    #axhline(0, ls=':', lw=1)\n",
    "    ylabel('Frequency convergence (THz)')\n",
    "    xlabel('Number of samples')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def monitor_phonons(directory='phon', dfset='DFSET', prefix='cryst', kpath='cryst', sc='../sc/CONTCAR',\n",
    "                    order=1, cutoff=10, born=None, charge=None, k_list=None, \n",
    "                    fig_out=None, once=False):\n",
    "\n",
    "    def update_fig(fig, bnd_lst, kpnts, k_lst):\n",
    "        if fig is not None:\n",
    "            plt.close(fig)\n",
    "        fig = figure(figsize=(14,5))\n",
    "        (dcplt, omplt) = fig.subplots(1, 2)\n",
    "        sca(dcplt)\n",
    "        show_dc_conv(bnd_lst, kpnts)\n",
    "        sca(omplt)\n",
    "        if N>1:\n",
    "            if k_lst is None:\n",
    "                plot_omega(build_omega(bnd_lst, kpnts))\n",
    "            else :\n",
    "                plot_omega(build_omega(bnd_lst,\n",
    "                                        ([l for l in kpnts[0] if l in k_lst],\n",
    "                                         [v for l, v in zip(*kpnts) if l in k_lst])))\n",
    "        show()\n",
    "        clear_output(wait=True)\n",
    "        return fig\n",
    "\n",
    "    bnd_lst = {}\n",
    "\n",
    "    if get_dfset_len(f'{directory}/{dfset}') < 1:\n",
    "        print('Waiting for the first sample.', end='')\n",
    "        sys.stdout.flush()\n",
    "        while get_dfset_len(f'{directory}/{dfset}') < 1:\n",
    "           sleep(15)\n",
    "           print('.', end='')\n",
    "           sys.stdout.flush()\n",
    "        print('done.', end='')\n",
    "    print('Calculating the plots.',)\n",
    "    sys.stdout.flush()\n",
    "    clear_output(wait=True)\n",
    "\n",
    "    N = get_dfset_len(f'{directory}/{dfset}')\n",
    "    run_alamode(d=directory, dfset=dfset, prefix=prefix, kpath=kpath, sc=sc,\n",
    "                o=order, n=N, c2=cutoff, born=born, charge=charge)\n",
    "    bnd_lst[N] = loadtxt(f'{directory}/{prefix}.bands').T\n",
    "    prev_N = N\n",
    "\n",
    "    with open(f'{directory}/{prefix}.bands') as f:\n",
    "        p_lbl = [v if v!='G' else '$\\\\Gamma$' for v in f.readline().split()[1:]]\n",
    "        p_pnt = [float(v) for v in f.readline().split()[1:]]\n",
    "    kpnts = (p_lbl, p_pnt)\n",
    "\n",
    "    fig = update_fig(None, bnd_lst, kpnts, k_list)\n",
    "    if fig_out is not None :\n",
    "        fig_out.append(fig)\n",
    "\n",
    "    while True :\n",
    "        N = get_dfset_len(f'{directory}/{dfset}')\n",
    "        if N > prev_N:\n",
    "            r = run_alamode(d=directory, dfset=dfset, prefix=prefix, kpath=kpath, sc=sc,\n",
    "                            o=order, n=N, c2=cutoff, born=born, charge=charge)\n",
    "            if r[0]:\n",
    "                bnd_lst[N] = loadtxt(f'{directory}/{prefix}.bands').T\n",
    "                fig = update_fig(fig, bnd_lst, kpnts, k_list)\n",
    "                if fig_out is not None :\n",
    "                    fig_out[-1]=fig\n",
    "                prev_N = N\n",
    "        else :\n",
    "            SN = N//2\n",
    "            all_done = True\n",
    "            while SN > 0:\n",
    "                for NN in range(N, 1, -SN):\n",
    "                    if NN not in bnd_lst:\n",
    "                        all_done = False\n",
    "                        r = run_alamode(d=directory, dfset=dfset, prefix=prefix, kpath=kpath, sc=sc,\n",
    "                                        o=order, n=NN, c2=cutoff, born=born, charge=charge)\n",
    "                        if r[0]:\n",
    "                            bnd_lst[NN] = loadtxt(f'{directory}/{prefix}.bands').T\n",
    "                            fig = update_fig(fig, bnd_lst, kpnts, k_list)\n",
    "                            if fig_out is not None :\n",
    "                                fig_out[-1]=fig\n",
    "                    if get_dfset_len(f'{directory}/{dfset}') > prev_N:\n",
    "                        SN = 0\n",
    "                        all_done = False\n",
    "                        break\n",
    "                SN = SN//2\n",
    "            if all_done:\n",
    "                if once :\n",
    "                    break\n",
    "                sleep(30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def load_dfset(base_dir='phon', dfsetfn='DFSET'):\n",
    "    '''\n",
    "    Load contents of the DFSET flie and return dfset array\n",
    "    '''\n",
    "    N = get_dfset_len(f'{base_dir}/{dfsetfn}')\n",
    "    dfset = loadtxt(f'{base_dir}/{dfsetfn}').reshape(N,-1,6)\n",
    "    nat=dfset.shape[1]\n",
    "    \n",
    "    confs = []\n",
    "    sets = []\n",
    "    es = []\n",
    "    with open(f'{base_dir}/{dfsetfn}') as dff:\n",
    "        for l in dff:\n",
    "            if 'set:' not in l:\n",
    "                continue\n",
    "            s, _, c, _, e = l.split()[2:7]\n",
    "            s = int(s)\n",
    "            c = int(c)\n",
    "            e = float(e)\n",
    "            confs.append((s,c, \n",
    "                          dfset[s-1,:,0:3]*un.Bohr, \n",
    "                          dfset[s-1,:,3:]*un.Ry/un.Bohr, e))\n",
    "            \n",
    "    return confs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# hide\n",
    "confs = load_dfset(base_dir='example/VASP_3C-SiC_calculated/1x1x1/T_300K/', \n",
    "                   dfsetfn='DFSET.dat')\n",
    "len(confs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def plot_stats(confs, T=None, sqrN=False, show=True, plotchi2=False):\n",
    "    '''\n",
    "    Plot monitoring histograms for the configuration list in confs.\n",
    "    If len(confs)<3 this function is silent.\n",
    "\n",
    "    confs - configuration list\n",
    "    T     - target temperature in Kelvin\n",
    "    show  - call show() fuction at the end (default:True)\n",
    "    '''\n",
    "\n",
    "    if len(confs) < 3:\n",
    "        return\n",
    "\n",
    "    #E0 = Vasp2(restart=base_dir+'/../calc/').get_potential_energy()\n",
    "    #es = [(Vasp2(restart=d).get_potential_energy()-E0)/nat\n",
    "    #          for d in sorted(glob(base_dir+'/../calc/T_600.0K/smpl/0*/'))]\n",
    "    \n",
    "    es = array([_[-1] for _ in confs])\n",
    "\n",
    "    if T is None:\n",
    "        T = 2*es.mean()/3/un.kB\n",
    "    \n",
    "    nat = confs[0][-3].shape[0]\n",
    "    \n",
    "    E_goal = 3*T*un.kB/2\n",
    "    Es = sqrt(3/2)*un.kB*T/sqrt(nat)\n",
    "    e = linspace(E_goal - 3*Es, E_goal + 3*Es, 200)\n",
    "    n = len(es)\n",
    "\n",
    "    plt.hist(es, bins='auto', density=False, label=f'{n} samples', alpha=0.5, rwidth=0.4, zorder=0)\n",
    "    h = histogram(es, bins='auto', density=False)\n",
    "    de = (h[1][-1]-h[1][0])/len(h[0])\n",
    "    N = len(es)\n",
    "    if sqrN :\n",
    "        plt.errorbar((h[1][:-1]+h[1][1:])/2, h[0],\n",
    "                     yerr=sqrt(h[0]), fmt='+', color='C0', alpha=0.66,\n",
    "                     capsize=6, label='$\\\\pm\\\\sqrt{n}$')\n",
    "        # plt.errorbar((h[1][:-1]+h[1][1:])/2, h[0],\n",
    "        #              yerr=2*sqrt(h[0]), fmt='+', color='C0', alpha=0.33,\n",
    "        #              capsize=4, label='$2/\\\\sqrt{N}$')\n",
    "\n",
    "    plt.axvline(E_goal, ls='--', color='C2', label='Target energy')\n",
    "    pdf = N*de*stats.norm.pdf(e, E_goal, Es)\n",
    "    plt.fill_between(e,  (pdf-sqrt(pdf)).clip(min=0), pdf+sqrt(pdf), label='$\\\\sigma, 2\\\\sigma, 3\\\\sigma$', color='C1', alpha=0.1, zorder=9)\n",
    "    plt.fill_between(e,  (pdf-2*sqrt(pdf)).clip(min=0), pdf+2*sqrt(pdf), color='C1', alpha=0.1, zorder=9)\n",
    "    plt.fill_between(e,  (pdf-3*sqrt(pdf)).clip(min=0), pdf+3*sqrt(pdf), color='C1', alpha=0.1, zorder=9)\n",
    "    plt.plot(e, pdf, '--', color='C1', label='Target normal dist.')\n",
    "    fit = stats.norm.fit(es)\n",
    "    plt.plot(e,  N*de*stats.norm.pdf(e, *fit), '--', color='C3', label='Fitted normal dist.', zorder=10)\n",
    "    if plotchi2 :\n",
    "        fit = stats.chi2.fit(es, f0=3*nat)\n",
    "        plt.plot(e,  stats.chi2.pdf(e, *fit), '--', color='C4', label='Fitted $\\\\chi^2$ dist.', zorder=10)\n",
    "    plt.xlabel('Potential energy (meV/at)')\n",
    "    plt.ylabel('Samples')\n",
    "    plt.xlim(E_goal-3*Es,E_goal+3*Es)\n",
    "    plt.legend(loc='upper left', bbox_to_anchor=(0.7,0.5,0.5,0.5))\n",
    "    if show :\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEGCAYAAADMsSqUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACTRElEQVR4nOydZ5Qb1dmAnztFZXu1vV6XtTHFxg3jRjGY3nsntIQEEkJCII3kIwlJSEIJCZACoZcQWmihl4ADmGqDsQ0GjHH32t7e1DX3+3FHW7ySVtqVtnmec3R2pRnN3JE0971vF1JKHBwcHBwcBgJtoAfg4ODg4LDz4gghBwcHB4cBwxFCDg4ODg4DhiOEHBwcHBwGDEcIOTg4ODgMGMZADyAVysrKZFVV1UAPw8GhC+ua1wFQVVA1oONwcEjE0qVLa6WU5QM9jmQMCSFUVVXFkiVLBnoYDg5d+PqLXwfgniPvGeCRODjERwixfqDH0BOOOc7BwcHBYcAYEpqQg8Ng5KLpFw30EBwchjyOEHJw6CX7jN5noIfg4DDkcYSQg0Mv+az+MwD2KNljgEcyeAmHw2zatIlAIDDQQ9kpeeWVV6Z9/PHH6wZwCBawMhKJfHPvvffeHm8HRwg5OPSS696/DnACE5KxadMm8vPzqaqqQggx0MPZ6YhGo5GpU6fWDtT5LcsSNTU1U7Zu3XoncHy8fZzABAcHh6wRCAQoLS11BNBOiqZpsry8vAmYmnCffhyPg4PDTogjgHZuNE2TJJE1jhBycHBwcBgwsuYTEkJ4gDcAt32ef0spfyWEuBr4FlBj7/pzKeXz2RqHg0PGkRKsKLplIZAQCQICNF09HBLy51e+yOjxLj9stx73+cY3vsGzzz7LiBEjWLlyZfvr9fX1nHHGGaxbt46qqioeffRRiouLAfjDH/7AXXfdha7r3HLLLRxxxBEZHXe6XH311eTl5fGjH/1oQMeRDbKpCQWBg6WUM4CZwJFCiPn2tj9LKWfaD0cAOQx+IkEINEHTZqhfAw1fccWkU7l8l1OhcSM0roe6NVC/Fpq3QrAVopGBHrUDcMEFF/Diiy92e/3aa6/lkEMOYfXq1RxyyCFce+21AHz66ac8/PDDfPLJJ7z44otccsklRKPR/h72TkPWhJBUtNpPTfvhtHF1GDpYFgSaoWEDNG6AthqwImDkgCufKaNmM6ViDrhywZUH7nzQXRAJQMtWqP8Kmqsh7Ffak8OAcMABB1BSUtLt9aeffprzzz8fgPPPP5+nnnqq/fUzzzwTt9vNhAkTmDRpEu+//36391955ZVMmTKF6dOnt2sozzzzDPPmzWOvvfbi0EMPZdu2bYDSZM4//3wOP/xwqqqqeOKJJ/jJT37CtGnTOPLIIwmHw4AqUfbTn/6UuXPnMnfuXL788stu512zZg1HHnkke++9NwsWLOCzz1SqwGOPPcbUqVOZMWMGBxxwQN8/uH4iqz4hIYQuhFgGbAdekVK+Z2+6VAixXAhxtxCiOMF7LxJCLBFCLKmpqYm3i4NDdmgXPuugRU0iuPLAzFVCxna0r6z7hJV1n3R9r9DAcCvB5M5XGlTjRvUI+/v3OhySsm3bNioqKgCoqKhg+3aVxrJ582bGjh3bvt+YMWPYvHlzl/fW19fz5JNP8sknn7B8+XKuuuoqAPbff3/effddPvroI84880yuv/769vesWbOG5557jqeffppzzjmHgw46iBUrVuD1ennuuefa9ysoKOD999/n0ksv5Qc/+EG3cV900UX85S9/YenSpfzxj3/kkksuAeA3v/kNL730Eh9//DH/+c9/MvMh9QNZFUJSyqiUciYwBpgrhJgK3ArsgjLRVQM3Jnjv7VLK2VLK2eXlg7oIrMNwIuxXWk/LNiVw3Hmgm3F3vWvFndy14s7kxzPcShhJCY2blIYUDWdh4A6ZQsbRWneM8CsoKMDj8fDNb36TJ554gpycHEDlRR1xxBFMmzaNG264gU8+6VikHHXUUZimybRp04hGoxx55JEATJs2jXXr1rXvd9ZZZ7X/feedd7qct7W1lbfffpvTTjuNmTNncvHFF1NdXQ3AfvvtxwUXXMAdd9wxpMyH/RIdJ6VsBBYBR0opt9nCyQLuAOb2xxgcHJJiWdBWqzQWIZTwyWSQgWELtJBPaViB5swd26FXjBw5sn0Cr66uZsSIEYDSfDZu3Ni+36ZNmxg9enSX9xqGwfvvv88pp5zCU0891S5Qvve973HppZeyYsUK/vGPf3SpFOF2uwHQNA3TNNsFm6ZpRCId/sPOAm9H4WdZFkVFRSxbtqz9sWrVKgBuu+02rrnmGjZu3MjMmTOpq6vr2wfUT2RNCAkhyoUQRfb/XuBQ4DMhREWn3U4CVsZ5u4ND/xEzmfkbldlNd2XvXKYXDK/SiFq2gjV0VqzDjeOPP5777rsPgPvuu48TTjih/fWHH36YYDDI2rVrWb16NXPndl0rt7a20tTUxNFHH81NN93EsmXLAGhqaqKysrL9mL3hkUceaf+7zz5d6xMWFBQwYcIEHnvsMUBpbR9//DGgzH3z5s3jN7/5DWVlZV0E6WAmm2V7KoD7hBA6Stg9KqV8VgjxgBBiJipIYR1wcRbH4OCQnGCLCh7QbT9Of6DpykQXaoVICPJHKU1pJyCVkOpMc9ZZZ7Fo0SJqa2sZM2YMv/71r7nwwgu58sorOf3007nrrrsYN25c+8S+5557cvrppzNlyhQMw+Bvf/sbut5VK25paeGEE04gEAggpeTPf/4zoAIQTjvtNCorK5k/fz5r165Ne7zBYJB58+ZhWRYPPfRQt+0PPvgg3/nOd7jmmmsIh8OceeaZzJgxgx//+MesXr0aKSWHHHIIM2bM6GIOHKyIePbPwcbs2bOl09TOIeP4GlTEm5nTK9Pb5YsuB+DPC//c+zFEgirirrBSaUnDjFWrVjF58uSBHsaQIdbAs6ysLCPHW7lypW/q1KmrMnKwPvDxxx+XzZgxoyreNqeAqcPOh5TgqwNfvTK/9bKszHdnfrfvYzHcENWgaSPkj1Z+IweHnQhHCDnsXEgJrTUQaOyTAAKYVDQpM2PSTRXa3bwZCkYrU53DTknnKLmdBad2nMPOQ0wABZvURN/HwppLti1lybalmRmbpqs8pOZq5adycNhJcDQhh52DmAku2KQ0oAzw4Kp/AjB75N4ZOZ4SRDnQvAUKKh3TnMNOgaMJOewc+BqUD8jspwi43hLTiFq2OBUWHHYKHCHkMPwJNIGvts8+oH5D00H3QNMmu0K3g8PwxRFCDsObkE+V4HHlDg0BFEM3QXMp05xTjbvX1NXVMXPmTGbOnMmoUaOorKxsfx4KhTJ6rsbGRv7+979n9Jg7A44Qchi+RIIq4sz0quizoYbhUr6slmpVVsghbUpLS9vL23z729/m8ssvb3/uciVOEO5cRidVBkII9Wacg40heGc6OKSAFVWRZpoLtOzE31y+9+VcvvflWTl2O6ZXCdO22uyeZyfijjvuYM6cOcyYMYNTTjkFn88HqL5DV1xxBQcddBA//elPWbNmDfPnz2fOnDn88pe/JC+vI1DkhhtuYM6cOUyfPp1f/epXgGrtsGbNGmbOnMmPf/zjbuf95z//ydy5c9sLj8aKjObl5fF///d/zJgxg/nz57e3f6ipqeGUU05hzpw5zJkzh8WLFwOqKsNFF13E4YcfznnnnUdNTQ2HHXYYs2bN4uKLL2b8+PHU1tbyi1/8gvvvv7/9x/+9732v8pprrhmRpY+11zjRcQ7DDyntFgyWqtOWJcblj8vasbvgylV5TaYbPIX9c84s8fUXv97ttSOqjuDMPc7EH/FzyauXdNt+wqQTOHHSiTQEGrhi0RVdtt1z5D1pj+Hkk0/mW9/6FgBXXXUVd911F9/73vcA+OKLL3j11VfRdZ1jjz2Wyy67jLPOOovbbrut/f0vv/wyq1ev5v3330dKyfHHH88bb7zBtddey8qVK9vryHVm1apVPPLIIyxevBjTNLnkkkt48MEHOe+882hra2P+/Pn87ne/4yc/+Ql33HEHV111FZdddhmXX345+++/Pxs2bOCII45oL1a6dOlS3nrrLbxeL5deeikHH3wwP/vZz3jxxRe5/fbbAbjwwgs58sgjjeuvv55oNMpTTz1V/MEHHwx49YQdcYSQw/DD3wjhtoyFYifi7S1vA7Dv6H2zeh5ACaKWrarGnenJ/vmGMStXruSqq66isbGR1tbWLq27TzvttPY6ce+88057o7uzzz67vXHdyy+/zMsvv8xee+0FqGKmq1evZty4xIuS//73vyxdupQ5c+YA4Pf726t2u1wujj32WAD23ntvXnnlFQBeffVVPv300/ZjNDc309KicsiOP/54vF61wHrrrbd48sknATjyyCPbW5RXVVVRVFTE4sWLvdXV1eaee+7pGzVq1KCrmOsIIYfhRcin6sFlWQABPPaFKnjZL0JIaGB4VKBC8fjMtpnoR5JpLl7Dm3R7sae4V5rPjlxwwQU89dRTzJgxg3vvvZdFixa1b8vN7TmEX0rJz372My6+uGvt5WTVDqSUnH/++fzhD3/otq1zWwdd19v9PJZl8c4777QLm850Hmey+p8nn3xy+M477yzbvn27+fWvf31Q9nZwfEIOw4doRDnxTe/QioRLFd0FAmjd7rQL7wMtLS1UVFQQDod58MEHE+43f/58Hn/8cQAefvjh9tePOOII7r77blpbWwHViXX79u3k5+e3ayo7csghh/Dvf/+7vXtrfX0969evTzrOww8/nL/+9a/tz+OZ+UB1c3300UcBpaU1NDS0bzv00EOjr7/+euHHH3+ce8oppzQlPeEA4Qghh+GBlGpyRmQtEGFQYHhVWZ+g0xSvt/z2t79l3rx5HHbYYeyxxx4J97vpppv405/+xNy5c6murqawUPnjDj/8cM4++2z22Wcfpk2bxqmnnkpLSwulpaXst99+TJ06tVtgwpQpU7jmmms4/PDDmT59Oocddlh7Q71E3HLLLSxZsoTp06czZcqULn6pzvzqV7/i5ZdfZtasWbzwwgtUVFSQn6/qD7pcLvbdd9/m448/vt4wBud94bRycBge+Bv7zQwXIyOtHHqDtJTPq2i8qsI9iBnKrRx8Ph9erxchBA8//DAPPfQQTz/99EAPqxvBYBBd1zEMg3feeYfvfOc77VrT8uXLfaeffrp47LHH1kybNm3AMp+dVg4Ow5tIENq2D/6SPJlCaKC5VQRg4RjQHINGNli6dCmXXnopUkqKioq4++67B3pIcdmwYQOnn346lmXhcrm44447APj000856qijvEcfffT2gRRAPeFoQg5DG8tS5W2k7PfupNv9yr4/wjtAqRehNvAWQ27pwJw/BYayJjQccJraOThkm0ATREP915q7EwMmfGKYOaoyuCvXCdt2GLI4erzD0CUc6GjP3e9IXl//Xxat/69qzy0HoKyOEMon1LLVKevjMGRxNCGHoYllqWg4w90/4dhWGCIhiPghFAAZYfEqFba7MH+i2kcIFUZtelVOj+ECkeV8Ht2lzHL+hkFtlnNwSETWhJAQwgO8Abjt8/xbSvkrIUQJ8AhQBawDTpdSNiQ6joNDXPrDDCejKvk10AzRIKCpJFHNAM0kHAsFd9mamERpRYFmkI0d29wFtrksS8LSMcs5DGGyaY4LAgdLKWcAM4EjhRDzgSuB/0opdwX+az93cEidSBB8WTTDySj4m6BhI7TZSeauXHB5lXaTKBpNALqhBIErR40vElLmssZNEGpFSaoMEzPLtW53zHIOQ46sCSGpaLWfmvZDAicA99mv3wecmK0xOAxDYkmpmisLZjipkkAbNkKgwRYmXiVYeoNACS1Xjgqrbt0OTVsgEsjoqAFllosElYbo4DCEyKpPSAihA0uBScDfpJTvCSFGSimrAaSU1UKIuCFGQoiLgIuApIUBHXYygi1qEs90Umo0BK016q/hyXzujW6oRzSs6r+5CyCnOLM+I1eO3UE2t9/D1VOmcWNmhbDhgaKxGTvcxo0bOe+889i6dSuapnHRRRdx2WWXpX2cQCDAAQccQDAYJBKJcOqpp/LrX/86Y+McTmRVCEkpo8BMIUQR8KQQYmoa770duB1UnlB2RugwpIhG7GCETLZnkOBvBn+tSgB1pW7iu2L2D9M/nW6qR6hNPfJGqECGTCA0EIbqPVQ4OjPHzDSZXkCEWnveJw0Mw+DGG29k1qxZtLS0sPfee3PYYYcxZcqUtI7jdrt57bXXyMvLIxwOs//++3PUUUcxf/78jI53ONAvIdpSykZgEXAksE0IUQFg/93eH2NwGAb46pQJLlMVpK0ItGwHf72qtmCYab290JVPoSu/d+c2PSrAoaVa+Z8y5SsyPWpiDsYvpLmz8vHHH3PAAQcwZcoUNE1DCNHejK4zFRUVzJo1C4D8/HwmT57M5s2b0zoGgBCivQleOBwmHA63V8pO5zg7A9mMjisHwlLKRiGEFzgUuA74D3A+cK39d/AVY3IYfIT9yt/h7uWkvyNRO2BAyrS0n868vnERAAeNXdi7MeiGEkT+eqUh5JVlxjxnepVp0cwZsi0fMkkgEOCMM87g/vvvZ+7cufziF78gEAhw9dVXJ33funXr+Oijj5g3b16vjhGNRtl777358ssv+e53v9vr4wx3sqkJVQCvCyGWAx8Ar0gpn0UJn8OEEKuBw+znDg6JsSylsWSqWGfYD02blfmqDyHNizYuYtHGRX0bi0AJwUgAmreqfKS+ohmApYq6OvDqq68ya9Ys5s6dC8D06dOpr69v10zi0drayimnnMJNN91EQUFBr46h6zrLli1j06ZNvP/++6xcubJXxxnuZE0TklIuB/aK83odcEi2zuswDInlBLkz4EsItqhip4Z3cGkJpgfCQRU9VzBa+Y36gpED/jr1mQ3yStvZZuXKlUybNq39+YcffthucotHOBzmlFNO4Wtf+xonn3xyr47RmaKiIhYuXMiLL75IJBLp9XGGK07ZHofBTTRiR3xlICco2GwHNgxSM5XpVua45i1K6PYFIVQYe1vNTt8Ar7S0lOXLlwPwxRdf8MQTT3DmmWcCqtlczOcDqkvphRdeyOTJk7niiitSOka849TU1NDY2AioVt6vvvoqe+yxR4/H2RlxyvY4DG589XbUVx/XS8FmFTXmyh3cXVcNE6IojahwtMr/6fWx3ErzC/sGpMBrXAxPZiPajJ7NqWeddRb/+c9/mDp1KmVlZTz00EOUlpZiWRZffvklJSUl7fsuXryYBx54gGnTpjFz5kwAfv/73yc8BhD3ONXV1Zx//vlEo1Esy+L000/n2GOPpbW1NeFxdlYcIeQweAkHINDY92CEmAAycwa3AIoRM8VlQhCZXqX9FY0fHH2HMpjTkyp5eXk888wz3V7/9NNPOeWUU/B6O0Lk999/fxK1t4l3jETHmT59Oh999FHKY9mZcfoJOQxOpFSJjX3tExRqVZNwFgSQP6r6hHn1LPlcImFVQqivPqJQG+SUquTYfsbpJzSwDIV+QoNgaeTgEIdQqypD0xcBFPZD6zYVhJAFDciru7MngECZ5oRmt2qI9P44plf51aIZiLxzcMgwjhByGHxYlp3n0oeK0LE8IMObNTPUi+te5sV1L2fl2O0YLtWrqGVb73sWxXxqPqdYvcPgwxFCDoOPQJOacLVeuiytMDRXKxNWFqPg3tnyNu9seTtrx2/H9Khr8tXS68oKhlf518JZKJ7q4NAHnMAEh8FF1J5se9umQdqJrYi+59oMJkwvb69az2vrVuM3CtpfXlPTBsAu5d2j3w7eYwSH7DFSPYk13PPVQkHl0AjQcNgpcISQw+DC19CHkGypJlkrnLmioIOIfXcfx74T7aKndhHQO99cC8A3F0zo+QCDMWTbYafHMcc5DB4iQWUy6m2VbH8zBFuHpQAC7P5Edsh1JNi7Yxgep/mdw6DCEUIOg4e2WmUy6o2pKOxXZWqy1W11sKBp6jNq3abCt9NFN1UVitDOUWW7rq6OmTNnMnPmTEaNGkVlZWX781Coj1UpdqCxsZG///3vGT1mJqmsrJxWXV2d1Pp1yy23lJ533nnjAK6//vryv/71rwkzaZ999tn8V155pc8qtWOOcxgchHwqn6U3ialW2C7H4+lXX8ev9726387VBd1UBU9ba1UeVbrXbHrt6hF5g7N8UQYpLS1l2bJlAFx99dXk5eXxox/9qMf3RSIRDCO96TEmhC655JLeDDXj4+krP/nJT2qSbX/ttdfy8/LyoocddlhbX87jaEIOA4+UalLsVaFN+71ksM/QUMDwQNiHGe1FCRxNV5/5TtoK/I477mDOnDnMmDGDU045BZ/PB8AFF1zAFVdcwUEHHcRPf/pT1qxZw/z585kzZw6//OUv2/sDAdxwww3MmTOH6dOnt/cCuvLKK1mzZg0zZ87kxz/+cZdzrlu3jsmTJ/Otb32LPffck8MPPxy/3w/AsmXLmD9/PtOnT+ekk06ioUGF0i9cuJCf//znHHjggdx8880sXLiQyy+/nAMOOIDJkyfzwQcfcPLJJ7Prrrty1VVXtZ/rxBNPZO+992bPPffkkUce6VFy3XzzzaVVVVVT58yZs/vbb7/dfpFXXHHF6F/+8pcjAa655poRu+yyy5677bbblGOPPXbi559/7rr//vvLb7vttpF77LHHlBdffLHX1YUdTchh4Am1qpV9b7Qgf6MKO85EgdM0eXqNKr9ywi7H9fu5ATBz8IQ3Ymm9SOg1vapJoDu/f6MI7zmm+2t7nghzv6W04QdP67595tmw19egrQ4ePa/rtq8/l/YQTj75ZL71rW8BcNVVV3HXXXfxve99D1BFRV999VV0XefYY4/lsssu46yzzuK2225rf//LL7/M6tWref/995FScvzxx/PGG29w7bXXsnLlynbNa0dWr17NQw89xB133MHpp5/O448/zjnnnMN5553HX/7yFw488EB++ctf8utf/5qbbroJUNrV//73P0CVDXK5XLzxxhvcfPPNnHDCCSxdupSSkhJ22WUXLr/8ckpLS7n77rspKSnB7/czdepU87vf/a4+atSouLbb9evXm9dee+3opUuXriopKYnuu+++u0+dOtW343633HLLqPXr16/wer2ytrZWLysri5533nk1eXl50d/85jfb0v4SOuFoQg4DS3tiai+CCcJ+8DcMmB/ow21L+XDb0gE5NwACLOHCG6pL3z8Ui0DcCXsOrVy5kgULFjBt2jQefPBBPvnkk/Ztp512GrquNOp33nmH005TQvHss89u3+fll1/m5ZdfZq+99mLWrFl89tlnrF69usfzTpgwob0o6t577826detoamqisbGRAw88EIDzzz+fN954o/09Z5xxRpdjHH/88QBMmzaNPffck4qKCtxuNxMnTmTjxo0A3HLLLcyYMYP58+ezbds28cknnyTM+n7jjTdy58+f3zJ69OiIx+ORJ598cn28/XbffXf/SSedNOHvf/97iWmaGa315mhCDgNLqKV3ialWRPmBdLeKGttJkZqBsKJKq8kdkd6bDa/q6uop7Ft5pHRIprm4cpJvzy3tleazIxdccAFPPfUUM2bM4N5772XRokUdp8jt2c8upeRnP/sZF198cZfX161bl/R9bneHuVnX9XZzXDJ2HE/sGJqmdTmepmlEIhEWLVrEq6++yjvvvENOTg6zZ8+2/H5/UmUjlYZ6r7/++uoXXngh/6mnniq6/vrrR69evXplj29KEUcTchg4rKhd3TpdLUiqyVOgWmTv5FiaW4WmB9OMeIv1HPLVZWdgg5SWlhYqKioIh8M8+OCDCfebP38+jz/+OAAPP/xw++tHHHEEd999N62tyh+3efNmtm/fTn5+Pi0t6X0HhYWFFBcX8+abbwLwwAMPtGtFvaGpqYni4mJycnL47LPPWLFiRdI5/oADDmh7991387du3aoHg0Hx5JNPdqtyG41GWbNmjeu4445r+fvf/76ppaVFb2pq0vPz86MtLS19dsQ6d7DDwBFoBmT6ianBVvXoz4RLK6qqOUTDEA1CNExx0IfArvYtNNWQTjdVczrNpf7vr2g9w454Mzzp+XhMu79PONC3Wn1DiN/+9rfMmzeP8ePHM23atISC46abbuKcc87hxhtv5JhjjqGwsBCAww8/nFWrVrHPPvsAqj3DP//5T3bZZRf2228/pk6dylFHHcUNN9yQ0njuu+8+vv3tb+Pz+Zg4cSL33HNPr6/tyCOP5LbbbmP69OnsvvvuTJs2LWlC2Pjx48M//elPt8yfP39yeXl5ePr06b5oNNrlRxuJRMTZZ589oaWlRZdSiosvvnhbWVlZ9JRTTmk89dRTd3nhhReKbrrppg11dXX6Bx98kHvTTTdtSWfMTisHh4EhGoGGtXaLhTSEUDQMTZvUZJvt/jjSUn6nYKsSPIgOX4rQuWnZX0DAD/a6TJV0k1H1sDr5Z8wcJSyzNMH/6z3lBzh73ljVHtwwIX8UadkoIyHQdSgck/HxDeVWDj6fD6/XixCChx9+mIceeoinn356oIeVFkOhlYOjCTkMDP7elOexy/IIPbsCKJbMGWylPfQ7TvDDD2Zd1vFEoMaFDjEDhUQJrzafuk5vYfpCNx1Mt4ow8zerc6WK4VKmvJBvQKIMBytLly7l0ksvRUpJUVERd99990APaViSNSEkhBgL3A+MAizgdinlzUKIq4FvAbFEqJ9LKZ/P1jgcBiGREAQawEwztSDQkt1wbCuizhFqUQLF8PYt6EGgqhvoqChAXyOIJhUIkK0246ZXVY5wedPryGq4bf/cWKe4qc2CBQv4+OOPB3oYw55sakIR4IdSyg+FEPnAUiHEK/a2P0sp/5jFczsMZvwNIIz0JrtoWDnQs1UXLtCi6tYJDYyclITPM1+pSK3jJsbJfdkRTQPNo4SRv1G1HPeWZN5MFws2aK1RrcFTlaK6S2l+TnFTh34mazYNKWW1lPJD+/8WYBVQma3zOQwRIkGVqZ+WMLHNcFqagisVomFo3moXTvUojSDFU3xe/xmf13+W3vk0zRY8GrRtV1XDe9usLhGGqZr6+ZvTfJ9bCa8h4Cd2GD70S4i2EKIK2At4z37pUiHEciHE3UKIuI3vhRAXCSGWCCGW1NQkLWHkMJRoqwMtzQz9oG2Gy3QuS6AFWqoBqYRif5qhdEP5h8I+1YCvt1WxExEzy6XT0ls3VR2+dEO9HRz6QNaFkBAiD3gc+IGUshm4FdgFmAlUAzfGe5+U8nYp5Wwp5ezy8vJsD9OhPwgHVDhwOiYoK9zLXKIkSEuZ9gINyu8zkM3vDLfyP7VuU0IxUwihhH1bDWl1YzU86vN2Wj049BNZFUJCCBMlgB6UUj4BIKXcJqWMSikt4A5gbjbH4DCI8NWm5yxHQlt9+v6jZEQj0LLVzovJGRxOeN2w2283KOGYKXOY4VIaVjqajWaoEPNh1OpB1/X29g0zZ85k3bp17LvvvoCqcvCvf/2rfd9ly5bx/PPpx0ktXLiQwZhG8rOf/cx1zz33xLU2xfj8889du+66654Ab7zxRs4FF1wwNtm+t912W0kmx5g1ISRULYi7gFVSyj91er2i024nARkr/+AwiAn7VQhwOpWyQ20QblOhx5kgElQCCNHLit1dyTXzyE03wi8RQiihGPJDy7auuUZ9wfQqE6iVhlku1uohU2MYYLxeL8uWLWt/VFVV8fbbbwOZE0LZIBKJ9Ps5DzjgAN+99967MdH21atXux955JGMCqFsRsftB5wLrBBCLLNf+zlwlhBiJspGsA64ON6bHYYRUiqHt57GxG9F1OTZ2y6rOxLy2cENroyV+vnOjCz8dE2PEpat21QtuL6OVdgJtm31kD8ytfdouoptDTRDTtJF9JAlLy+P1tZWrrzySlatWsXMmTM566yz+Nvf/obf7+ett97iZz/7Gcceeyzf+973WLFiBZFIhKuvvpoTTjgBv9/P17/+dT799FMmT56csA5cVVUV559/Ps888wzhcJjHHnuMPfbYg/r6er7xjW/w1VdfkZOTw+2338706dO5+uqr2bJlC+vWraOsrIzddtuNtWvXUl1dzRdffMGf/vQn3n33XV544QUqKyt55plnME2T3/zmNzzzzDP4/X723Xdf/vGPfyStCffmm2/mfPOb36zyer3WvHnz2vuBPPvss/k33njjyNdff/3L5557Lu+HP/zhOFD15d5+++3P/u///q/yq6++8uyxxx5TzjrrrNpf/epX2/v6XWRNCEkp3yJ+nNHgWGY49B9hf/qtGmLVnTORlBpqU2YuvR+qLGQCw62i21q3Qd6IvvusTI8SwmFf6hXHTY/6zDwFGe3TtP7c87q9ln/UkZScfTaW38/Gi7oL9sKTTqLo5JOINDSw+fuXddk2/oH7ezyn3+9vr149YcIEnnzyyfZt1157LX/84x959tlnARg5ciRLlizhr3/9KwA///nPOfjgg7n77rtpbGxk7ty5HHroofzjH/8gJyeH5cuXs3z5cmbNmpXw/GVlZXz44Yf8/e9/549//CN33nknv/rVr9hrr7146qmneO211zjvvPPaW0AsXbqUt956C6/Xy9VXX82aNWt4/fXX+fTTT9lnn314/PHHuf766znppJN47rnnOPHEE7n00kv55S9/CcC5557Ls88+y3HHJW4xcuGFF1b9+c9/3nDMMce0XnzxxXFLZdx4442jbrnllvWHH354W1NTk5aTk2P97ne/2xwTUj1+8CkyBO5IhyFNbxrWRQIqjyYTOTShVjWZGt6MC6AnVj/JE6uf7HnH3qC7lAbTsi29CLdEtIdfp2hiE5paQgbSDPMehHQ2x3UWQKnw8ssvc+211zJz5kwWLlxIIBBgw4YNvPHGG5xzzjkATJ8+nenTpyc8xsknnwx0tG8AeOuttzj33HMBOPjgg6mrq6OpSTUZPP744/F6OywARx11FKZpMm3aNKLRKEceeSSg2jnEjvf6668zb948pk2bxmuvvdalPcWO1NXV6S0tLfoxxxzTCvCNb3wjbgXb+fPnt/7oRz8ae80114yora3VTTM7ATxO2R6H7BL2KfOSO0XfibT7C6UVwJCAUCv46u3KB5kPQPiq6auMH7MLMQ2oZZsypfVFI9J0pV0FmlSSbCoYXmXCdOdnzISZTHPRvN6k243i4pQ0n0wipeTxxx9n991377YtlRYI0NF+Qdf1dj9PvJqdseMla99gmmb7frH2DYFAgEsuuYQlS5YwduxYrr76agKBQNJrSmXsv//977eeeOKJTU8//XThvvvuO/nFF1/8IoXLTRtHE3LIHjFfUDpaUKBFOcT7aoIKtWVVAPUbuqkESOv2vmtEhleZOVPNSRKa+uyGcRvwHdsv7Pj8iCOO4C9/+Uu70Pjoo48AOOCAA9rbQKxcuZLly5endd7O71+0aBFlZWUUFBT06hpiAqesrIzW1lb+/e9/J92/rKwsmpeXF33ppZfyAO699964q5JPPvnEPXfuXP/vfve7rdOmTWtbuXKlp7CwMNra2po5+yyOEHLIJqFWFZWVqkCJhlWfoL6a4UK+DhPcUBZAMXRTCYTW7SpgYwfW17eldhyBCsxoqyXl3CEjR32W0f6P1OoPpk+fjmEYzJgxgz//+c8cdNBBfPrpp8ycOZNHHnmEX/ziF4TDYaZPn87UqVP5xS9+AcB3vvMdWltbmT59Otdffz1z56aXaXL11VezZMkSpk+fzpVXXsl9993X62soKiriW9/6FtOmTePEE09kzpw5Pb7nrrvuWvf9739/3MyZM/fwer1xfwzXX3/9iF133XXP3XfffYrX67VOPfXUprlz5/oNw5C77777lF//+tcj1q1bZx544IGTej14nFYODtlCSmhYp/JOUu2a2rJNrdJTDMl+88taFn9Z2+W16tpGiq1GykvykXRdsM2ZUMK8qsxFl/5xicqz/tHsH2bsmEmJhtTfvJHtwQL/em8jb35Zw61fS+wY70bIB7llqQeKhH2q6GpuWZoDHtqtHIYDTisHh52XoG1WM1LUasI+lROURvHMBZPKWDCp08QYDfP0Wx8htZGcOCthvl3GKPL0c/iybiefttYoH1FvtTzDo8LfTW9qCwTDq4rOegoHtrqEw7DEEUIOmcey0ouIk1ZHV9BenzMCrduRQkeKjJqsE/LNqd/ol/PEePzDzQDoVpANzV+xqtFFY0BFu33nn0sBGFuSw7gSJcjPmT8u/oE029fjb0xNu4nlGvkbIc8poeWQWRwh5JB5QrYWlGq9t0CLElxGL4WHtKC1FhDIdIujDkGimptxBRajij18uE2yvt7HfpPSFA6GR4XBu/NTWyzEygp5CjNfSNZhp8YJTHDILJbVYepJhWhYVXvuS4HStnqQkX6fHB/+/FEe/vzRfj1njKjw4I604LZ8vTtAukEKQqgafv6G3p3PwSEBjhByyCzBZqWZpJpl74sVKO3l+QJNEPH3zZTXSza1bGRTS8IyW9lFCCKah3yrmTJPL4OLYn2Hgq097wvqMw40Zb7thMNOjSOEHDKHFbVDo1MNRvD3rUBpyKcmxUzVlxtqCI2wMBnrakWTvQyhNjxqIZBKJYVYewhHG3LIII4QcsgcgWYVmp2KFiQt1esmnaKmnYmG7dYQnt5rUcMAS+hIoeEJNwC96AGkaYBMPSHV9Kjv2dGGHDKEE5jgkBliWlC6wQiuXgihmADTzKFRkDTLRISJJiO4Iy0EjcL0DxCrpODKS61ckmYqv1/h6PTPtQPPr6jm+RXVXV77YpsyD+42snupp6OnVXD0tIpurzsMXRwh5JAZAs1KIxEpCAUrVhmhl2Y0f4PSuAY4SmtEzogBPX9nopoHM9pGVDOJaClWyo4hUILFVwf5o+hRtTQ9qhpGONDn6hbxhMrNr64G4LJDd+3TsfvCOeec015WJxWGQtL/YCWlZaQQ4nohRIEQwhRC/FcIUSuEOCfbg3MYIkQjyjSWqm/G39hRlyxdAi2qLlwGmtKlhZQdHU9rPoMl93JeUxPnNTXCu7fCh/erJNIBJCrceMKN6LIXNeYMl914MMUSQDGhlSW+2Jadzq6LFi3iggsuSLpPdXU1e+yxB1LKlB8OvSdVW8bhUspm4FhgE7Ab8OOsjcphaBFo6kho7IlYu+nerKAjIQjU918gQqAJvnwVFl0Hj5wDtZ+p1/118OXLsGUZVK+ETR/Aquc6nPufPw//+R68+3fY8G5mWjGkgtCQGHhDDQh60RU1VklBpuBbMtx259v4zdz6yurtKUbsZYGHH36Ys88+u8trJ510EldddRULFixg1KhRvPrqqwM0uuFHqua4WAbg0cBDUsr6VMuYOwxzohHbtJZKuR27t5DWCzOatOzOqO7sFyVtq4V3/grbVqjnnlKoOkD5TADG7Qfj9uP+Tx8A4Lwp5yKjFuFttYTXrSa3qAzyK9j+9Ef4ty/Dsu5Dyy/EGLMbrrGjKD/raHVJUQuhZ9anZWkmuhXEFemFJqHptqBvAW8KviXdrbShwrg90QYV8+bNIxgM0traSn19fXuTu+uuu47Zs2ezZMkSjjjiCADWr1/PxIkTu7x/5cqV7Lfffrz55ps88cQTPPjggxx66KH9fRnDklSF0DNCiM8AP3CJEKIcSNywwmHnIdCUumkt2ArRYFr14dpp9wNlyY0ZDUPzZiiuUt1Egy0w9XQYOxdKJnTZ9b119Xywtp4N/rWMWR/itXv/RvFXa/EEg4Q8Ht78zkUgTmc343VK9Y0UenyIQJDAmk1EW/zQuhXyRrHuxzcioxb5+8ygYP+9cI9Jsf12T5eiuTGjPqWppPtZmx6l6blzlMktGYZLfU5hf9+SjW1iQQqxwIRL//UhkJlghPfeew9Q5rh7772Xe++9t33b0qVL+dGPfsS8efPYsmVLtwZ1Pp+PpqYmLr/8cgAikQhFRUV9Go9DBynd0VLKK4UQ1wHNUsqoEMIHnJDdoTkMetpbL6Qw0cloR3+fdAm1qpygDEx0cdmwGD78pxJyJ/5dRYgdc2PC3edVqWrcz/0+yMR3mtHzw9RP3pWGqnHsc+CeTN9lrGoadlCcCvdttfDUJTByMvnTdqf18+3U/ut5ah98DnfVaMpOP4KCBWlUxE6AJVx2S3N3eg3phAChK6Gfm0LgRaxja9HYPmuoMWFz86ureW5FNX89u++fQyrsvffenHbaaTz22GNs3bqVyy7r2kL8k08+Ye+990bXVerB8uXLmTp1ar+MbWcg1cCEHOC7wK32S6OB2dkalMMQIZ0Ag/YcojTNT9Ew+BqyUxGhYR289HN4809qMp13cUK/lpSS1qWfsv5nN+P75EsA1u+dzzvnj2LXB35P9cnH0TRrBt5J45J3rXTnw6xzoXEDZa7HqTo1n0l3/JyRF52KcLuwgqpdg+UPEmnqvXNeCl1VyG6r7QioSBXTo7TWSArGDt2ltNss+YYyzcKFC7toQTHOOeccHnjgAfx+f7fmcitXrmw334ESQsnaeTukR6pLpHuApcC+9vNNwGPAs4neIIQYC9wPjEJl0d0upbxZCFECPAJUAeuA06WUTgr2UCMSUgUtzRTadkfDamWdkt+oE1IqR7lmZN4PVL8WXviJajs+9zsw6eCEAijw1Sa23/MUbcs+xygvJtqmJlxfqYmv1EToaRReNdww+TjY5WBY9k/4/HnMDe9TcvwtlBx3INJSQQENL75F7UMvUnbGERQffyCa2YvCrLpLhVEHm1Xh0XTQXJ1ygXr47HW7Bp3Zd20oxq4jUmwHnyIxn9COXHfddRxxxBFMnDiRaDTKwoULu+2zYsUK5s2b1/585cqVjiaUQVIVQrtIKc8QQpwFIKX0i54jEyLAD6WUHwoh8oGlQohXgAuA/0oprxVCXAlcCfy0l+N3GCj8DXbNtxQmHX+9MvGkOz8Fm1WLhr52Wu1MNKQmzZIJMOMs2PVQcCduq7z1tkdpeP4t9DwvI791CsVH7Y8w1W0zJr8PPYtcuTD3YphwENR+3h5yHrut8mbviW/ll2y/92kaXlrMyG+cBK7S9M8Tq/dmeNILazdM5VMKtvbc/E53qf3Cvt75++Kw28gUG+6lSMwnlIzrrruOfffdt9vrf/rTn7o8/+qrrzI2LofUhVBICOHFLrcrhNgFSFq3Q0pZDVTb/7cIIVYBlShf0kJ7t/uARThCaGgRCaqJzZXCajUS6J2DPHYOI83Ey2RsWAwf3ANH/EH1xZl6co9vcVWOpOSkgyk77XD0vK5jOXP30/s+pvLd1ANUNN6Se2DBFbjHjmHsLy6m9cNVbLvzCTb97g5GLNiH7YcemN7xO1fLLqhILYw+huFVfjxXjlpEJN3X9g0V56SlDcWrmBAjFpjQmWxWTNh///2zclyH5KQqhH4FvAiMFUI8COyH0mhSQghRBewFvAeMtAUUUspqIcTgSTt3SA1fvYqc6nGysc1pqZSC6fI2qyOUOxPWHSsKyx6EVU9D6W5Ja9tFW31s/cdj5M2ZSuEBe1NyXJqTfl/H6auH534M+/8Axs4jb9Zkcm+5krrHX+WrWAM6S4KWxgejGxCJKu01Jw1tStMgIu2Q7aIezmHaASQpaE6dcMrwOKS0LJJSvgKcjBI8DwGzpZSLUnmvECIPeBz4gZ3wmhJCiIuEEEuEEEtqagY2E92hE+FA6smmwVbb/JWmP8NvF9NMJ6orEaE2eO23SgBNOhwO+y3klMTdNfDVJtZefgPNb35IpL7ngp53rrybO1fe3fcxxqiYCcf8CYrGwRvXwydPACAMnbIzjiAwepTa7dkXGfnCqxBNIyE1llwaSrP/kGm39rZSSLg1PGrxYPWikKrDTkvSu1wIsWOMZExvHieEGCel7K4vd32/iRJAD0opn7Bf3iaEqLC1oApge7z3SilvB24HmD17tlMXY7Dgq09Ns2kPyU7TnxPyq86sZobMcCsfh22fwLzvwKTEyYWNr77L1lsfRc/PZfy1PyBnjwkJ921/TyAL8TQ5JXD4b+HtvyjtrWCMylWKISWWaVL67hI8W7ex9sBDCOak+FnpdtsGw516v6dYPcBUQrY1Qy1SQq0q16p9yDJ5xKDDsMayLEGSEu89LTUTJ0so/9DBiTbagQt3AauklJ09e/8Bzgeutf8+3cMYHAYLYX/q5pZAM5BmSLYVVUEMvW3vYLOpsdNqf8aZMHZ+h98lDr6VX1J984PkTN+Nyh9fgFGUWad42ugu2P8K2LhPVwEEIATbjjoUf2UFo59+gYP//RjvHnEkjEocXNGOpoElVMh7Xlnq44mFbLsLel5UmF6lDbnyQNPweDzU1dVRWlrqCKKdEMuyRE1NTSGwMtE+SYWQlPKgPpx/P+BcYIUQYpn92s9RwudRIcSFwAbgtD6cw6E/aatNTUD0NiQ70Kj+prpKT0BtbQ28dC8cdKWaDBMIoNgKPWfqJMb84mLy9p6S8TI6vUYIGGdHajVXw0f3Y+Se3F4hu3n6noTKSxn54OPs+/xzfLXHJUh3Chqq4VZRbOkGi+gphmxrOkQsFdnoLWLMmDFs2rQJx6Q+MGzdutWIRqNprDgyjgWsjEQi30y0Q0pGdyGEB7gE2B+lAb0J3CalTJjNJqV8i8S/1kNSOa/DICLkU5pQKlqQvxctu0P2xNhHM1xJ25eUtq0DV4My7SWI4AvXNbH5hnsY+c2T8U4aR/7cQZz30bQRNi1hgXszi8dfQkhX30GgYhSvn3Iq+Q0NlKcigGJ0Mcul6HfTTfUdpRKybXhUtQZ3PqZpMmFCz6ZNh+wwZcqUFVLKQV1YIFXP7/1AC/AX+/lZwAM4WszOgbQLj6aSZ9KbkOxYVJienv/o8Q83t/+/udFHXV0NpW3rCGPy8/CFWC9VM7q4kcoiJdhOmVUJQHDzdjZc9VeirW1Em1NsXRCHiYUTe94pE4ydCwt/RsGi61iw7hbeHP99QoYSBCGvlzqvl3Kg+L2lGC2t1BxyQPLIRU0DSwNfY3pmuVgr8J5CtjXdjqprhpzi1I/vsFOSqhDaXUo5o9Pz14UQH2djQA6DkLBP5e24e8oL6mVItq9BTZp96JI63VXNfqG/02gW8lsu5MCJ8YuBBtdXs/6qv4AlGf+HH+Cd1PuE05N3PanX702bylksHv8d9l1/K/uv/xtvVF1KRO/6fbi311CyZBl6IMjWow9LfjzD3dE/KNUFQzoh26ZXVT1352cmytFh2JLqXf+REGJ+7IkQYh6wODtDchhUSKmSEFPRgnoTkh3yQcTX5yZ1bWYp2/N258HCC4kkWFuFttSw/ue3IDTB+D98v08CaCCozdmVd8Z9i4CRjxTdr3HrsUdQu988Sj74kIpnX2ovAZQQ3a00GyuS+iBSDdkWmnoEeg51d9i5SXWJMg84TwixwX4+DlglhFgBSCmlU81vuBJsURNOT9URehOS3UszXGe84XoCZiEBs4h3x15E29ZmCjzxJ0hzRAn5+82k9MSDcY0u7/U5Y9z68T8A+M6Mi/t8rFSpyZ1MTe4egEC3QphWkLBmC3Ah2H7YQqSmUf7mO1T/9SEqvnd24qi03pjl0gnZNrzKP+gpSD9XzGGnIVUhdGRWR+EwOLHsygWptF/oTUh2oLFPZrjcUA0Hrr2JTYV7sXzUqe2v53u7Tnju7TVEcnMRhk7FJWf06lzxaAsPVPdPlXax7/pbmRoRPJZ3bqdNQvmEdI0R40b0HBbdm2i5VEO2Y912fQ2Q7xRGcYhPqv2E1gshioGxnd/TU7KqwxAn2KxK6PQUMt2bkOy+RsP5G9lvw99ARlhb0r3oZAxXTR3j732IQMVIWHBF7841SNjc6GNLQ0fLBBmdwqmRZzi8+XEeC5xIl3DEqim4Z46hFAhvr8coL04skHRPR7uMVMPjUw3ZNryq2rq3sM8mV4fhSaoh2r9FlexZg13ElB6SVR2GOFbULs+fghbkSzMk27LrmPU2KTUShNd/jyfcwltV36PFNTrubmZ9I+PvewiEYOtRPTjqhwCVRTntkX6Kw1hZE2Hv2hfIyRnBpyOP77L/vKoSQlu2s/YHN1B8/IGMOOfY+AeOJbEGGlOvLZdqyLYQqgagr14VUHVw2IFUzXGno9o5hLI5GIdBRKCpw/6fjLAfwmmGZAfsEoK9TUp95y/QsIYPxnyLem9V3F30llbG3/8wIhJh/de/Rqgsfr24oc7n5UfiiTQzuf5V/GYJa0u6VoI2K8opWLAXdY+8hO71UHpKgtJFhlsJFMOrQrBTIdWQbdNjtwEvzmxbDodhQapCaCVQRII6bw7DjGhYJRv2aF7rRUh2OACh5vSrKXRmj2Ng1EyqWyYn3KXi+VcwWttYd8FZBEf2PQghHruX7JGV46aHYHnFqUihU5O3a/etQjDqkjOx/EG23/s0RkkhhQfNiX8ow2ObVT2ptXzQNIja+UDeHvKBdDe02W3AHRw6kaoQ+gMqTHslnfoISSmPT/wWhyFLqm272yPnUlw5S0sJt95Gw7VshfxRUD5ZPTolq+5I9TGH466pJTAmvqkuExw38Zguz7vUrOtHJHqnwAyJGW0lrHeYyISuMfryc4k0tbDllgdxjRmJd9dx3Q+k6Spc29+UepKpYYdsu/KSR8AZLrvVgy/134vDTkGqQug+4DpgBUmqoToMAyIhFVbr6qE0i4xCW5oh2cFYBF0vzHDVy1RLhv1+AFUL4u9jScav+pQNu+1ONK8AX15munymyuZOQQMDxdStTzK6dQWvT+gahCFMgzE/+yb1T72OZ0Jl4gMYbruKuTc105lA+QN99ZAfP0G4HT3W+G5c5tu1OwxZUo2NrZVS3iKlfF1K+b/YI6sjcxgYfHV2M7keJomYzyjV8OpIUJlteqMFtWyFN2+EwnFQmcCUBJQveou9Fy2isp/aL9/00c3c9NHN/XKuVNlSMBNvuIF5G+6GaNckVD0vh/JzjkEYOpGmFsKJeiZpLrUQkSmuN007zDvcgxDWTbBCSoN2cLBJVRNaKoT4A6oNQ2dznBOiPZwI+9UE0VOBymjI9h2k0N4bVNUFX53djTXdMfng9d8robjwpwlX5wXLP6H8f4tZt8cebJo0iWx6a2I167Y2+/C1lnPlV8tpDqgJ/8rHlwMwutjbrWZdf1CfM5EPR5/N3C3/hA/ugPnf6baPtCw2/PLvCCEY/4fL0Lw7RCnqhqoBGGwGT1FqJ9ZN5R8sqiR5yLanS6sHB4dUfwV7AfOB36N6DN0I/DFbg3IYAGJFSlMJm27XllI8drBVJb72Jmt+8c3QshkW/AjyRsXdxbtxM6Offp628WP56IADs2PqkaiAjUgQ3QqgWwFcVoTCnK3MH59HVaGGmxDzx+cyf3w+YwrSrJ+XQTYWzmVV6WGw5lX47Llu24WmMeLcYwms3cym6+9GxuvQqnuU5hpJMSBWN5U/qSctRzOUhuWU83GwSTVZtS99hRyGAuEUWzWEfSrCLVXncjQMgXoweumfqZwDo/eCUQkqQ0UtKh9/hkh+PpvOOAnZkkYdtJ6IRkBGlIAWQglow03ALEaiUW8nX/pdpTTpOo34CJhFaDKMEVXCClCmSM3s15X/pyOOYXJBCPLi5+bkzd6TUd8+ja1/f4Tt9/6HkRfuUIxVoMbsq+PNbQaL19R12by+TgVhjC/t+B0IGcWwvmCvaVM5eEqSgBDT297qwSlu6pDyL0AIcQywJ9BuD5FS/iYbg3LoZywLWrf3HGQg7TI+6YRkt6WpNcWIhtR5dk3ckhsAXWPjGSchDZ1obg60NKd5oh2QEqJBda26W5mjdE+XyTKiNaoh2gEWUWESFS48ORDRvICXkA5goVsRNdmGWiEcVfk0ei8+j7TRYN/vd7ouq1vYdfFR+xPcUE39U6/h2XUchQfs3fUQugnhAAvG5rFg164Gzn+9txGAs+ftEHIdDoC7B206FnkZaITcgey35jAYSGlpJoS4DTgD+B7q9jkNGJ/FcTn0J6EWtervyVwWaLHNaimuXYKtKoQ7XTNc8yZ48mLY/FHCXWTUIu/zL9VpKkYSKu/jZGZF1QQaDYG7EApGq2gvV17C6x3pmsxIV0euUqF3R+GsEdVcqoBnwWjIG6l8WhGf8rlISb/w+XPwyq+6BSoAjLzwZEpOOoTc6QnanxseVXYn2kPV7C77N6nPMel+dnHTVM19DsOWVO0D+0opzwMapJS/BvZB1ZFzGOpYURU221N5HiusJo1UM96jEbV/OiHcoATBouvV/8WJ1zk1Dz7HuH/9m5x1GxLukxJWVJkhZRRySpSw8OSn1HF0knchk7wLUz+X4bLPUakEXSRoC6PeDz8lPIVQ8yl8eG+3TcLQGfmNEzGK8pGRKJEdm/wJ1Gfhq0/tXDEzXlstSS9MCBCx/Rx2ZlIVQrHYS58QYjQQAZyevcMBf6M9cfSQu+NvSC2BtctxjfSDBN67VQUi7H+5mrDj0Pz2Muoee5mGWTPwVcVJukwFiRJ4VqRD+Lhy+id/RdOVoCuoUKa6iC+7GsH4/WH34+CLF2DtGwl323TtXWz89a1YoR20Ht2lNJtQilXDDZcSsKEeutaaHttMOfD5VQ4DR6pC6FkhRBFwPbAUWAs8nK1BOfQTkRD468DoIcggElCmtVS1oFAbRPxqMkqHz5+H9W/B9LMSBiIEN1RTfdM/8exexdZjelmUNBpSE7/HNrulU/euE4ubbmVx0629GwPYwqgQ8keriT7sV5pZNph1jqoy8c6t0LAu7i6Fh8wj8MV6tv7tYeSOpkLDrrSdagM8w638gbKH6zHsBNb+Mk06DDqSCiEhxBwhxCgp5W+llI1AHqpqwmPAn/thfA7ZxFenTCJJV/9p1oezomqy6k2F7NbtMHpvmHpK/JGEI2z6/Z0It4sxV34DaaQZWSWlmuiFDvkVShsZDJn7uqGayuWUqUk+Euz5PemiGSrM3Z0Ptavj7lKwzwzKzj6Kptfep+H5N7tubO8N1Jji+XT1eQd6CBTRXSoQxElg3WnpSRP6BxACEEIcAFxrv9YE3J7sjUKIu4UQ2+16c7HXrhZCbBZCLLMfR/dt+A69JpaY2pN2E2xRTulUgwv60qhu7wvgwJ8kFAzCNCg/7zjG/PQbmGUp1jaz0ayw0s48RSrgYDB2+nR5VW0806NC4Xtqz50u3iI44S+wa2INsuyMI8mbM5Vtdz6B/7O1XTcabqVBhlKskRdrBd5TUIPhVb6hbGmBDoOanmYKXUoZ80ieAdwupXxcSvkLYFIP772X+B1Z/yylnGk/nk9vuA4ZQUplAumpyZgVUfXhUjXDhQPKFJdO8zIp4d1boVZFuiUKCAjXNQJQsO9Mcqb29NPrim7ZmkVM+xnMaLrq6ZNTCtFA6lFpqRL7bjZ/BF91r7wlNI3RV5xLzp67IMw434WuKm0LmYJZrr2uXF3y/TTdSWDdielRCAkhYr/EQ4DXOm1LaguRUr4BpBhS49CvBFuUCaQnE5u/0TbDpGCyaq+QnaYZ7tMnVWZ//ZqEu7Qu+YQ13/w1bcs+S+/Y0kK3/IR1D23ussGp/STClat8RQjlk8s0nz+rhH8c/5Cel8P4a76HZxcVANvFP2RruO5IikEK7XXletCezBz1+3FCtnc6ehJCDwH/E0I8jYqQexNACDEJZZLrDZcKIZbb5rqENhUhxEVCiCVCiCU1NTW9PJVDN6yo6uti9BCSHasdlqoWFGgi7QrZNZ/Bsgdh7HzY7Yi4u4S317P5xvtxjR2Jd/LElA8tZBRdBgkYhQSNIlKPwUmd0e4ZjHbPyPhx29ENyB+hzFphX2ad9/t+Xwm6N/+YMDpNWhZbb32Umgee6brBcGNEfWhWir4rw2MHHyQxtwmhfHV+Z926s5H0zpRS/g74Icq0tr/sWBJpqMTVdLkV2AWYCVSjatAlOvftUsrZUsrZ5eXZaUq2U+JvpGdhkWYwQsR2LKdTITvYAm/+GXLKYZ9L4o8iHGHTdXeDZTHmygvR3KmNR5cRNBnBb5YR0bPXzmGCZ18mePbN2vEBFQyQUwqeYnQZQPQUbZYqnkI44ApoqYb3/pHg1BoyGqXusVdoeXd5l22W5sITbkit0naqQQqmV+3jhGzvVPQYXiSlfDfOa1/05mRSym2x/4UQdwDP9uY4Dr0kElQh2T1Vv44FI7h60Jag9xWyP3tOrXqP+H3CLqvb7n6KwBfrqfzZhbhGd1+IbG70sWWHHj4GYUDw6gaIiq7K+viyHBZMylyZmIhU/hpD9IOZz5NPwCzGE27AkiYyhWTaHhmxJ0w/E5Y/DLsshIqZ3XYZedGpBNZsZMtN/2TCn37c/j1IoaPJiNKAe+qqCh1BCq7c5IubWMh20djBEbnokHX6tXqgEKJCSlltPz0J1Tbcob9oq+k5JNsKKy2opwoKMYLNKoorVbNdjGmnweiZUNa9JTUoP4RekEPJiQdTsO/MuPtUFuW0t0sAFYBgCR2/qxgZ56edSQEE8F7znQDsV9i9XUI2iGhe/KaON1yPtCSWlgHhN/UUKBoXVwABaC6TMT/7Jmt/cB2b/nAnVTf8EM2jhIglXHaEZU7PwSidCqKSP4qEKxbdpYJbgs1KW3MY9mStrK8Q4iHgHWB3IcQmIcSFwPVCiBVCiOXAQcDl2Tq/ww6E2lRobU/CIp3KCJGQWgmnU5qnaaN6j6ZDeeKuP0IIys86mhHfODGlw+pWgKjmwucqjSuAhgtRzYXfVYpAolsZiJwTGoydp/5v2aoiHHfAHFHC6B+eT3hbPYG1mzq9V3QIllT8VYarI4IyGaZXaUNxat05DD+ydrdKKc+K8/Jd2TqfQxIsC1q29azdhP0QbEutTYOUqp5YOma4sB8W/UGZA4++Pv5QAyE2XvMPyk47nNwZuyNSEIYqAi6HoFFAFtdVg4aoMPG5SskJ1aNbIVUkta8EW+CFH6nWGftd1m1z3t5TmHTX1ej5O5hOdRPCwdQ1F8PToWmLBH7J2CLIV68CMxyGNcP/jnVQmoe0khfllJYy16UajBBsVXlE6YQ9v3cbtG6DOV+PPwQp2XrrI/iWr0ZGU3B4S4lh+QnruQSNQnamn7MlDNrcJVgIdCsDYc3ufNjjOFj3Bqx5Pe4uen4uUkoaXlyMu3prxwbDrX5jqYRXx5KY/Y3J9zO8KvE5jmbmMLzYee7anZVICHw1ym6fjEBz6m0aomFV3j8dM9yXr9p14c5UNczi0PTKuzS99j5lZxxB3qz4+7QjwZABgnqerQHtfE5siYHfXYIltMwIommnqWCF929X7TTiYPkD1D78IpWPPIkWsAVEu7+nPkWznN3uIVn+kxB2kMJ2p67cMMcRQsMZKVMLRoiG7DYNKQYjtKUZDde0Ed6/C0ZOS1gXLvDVJrbe9ii5M3en7Myjkh9PAhEfQT2PkJHPQAmgsZ7ZjPXMHpBzx5AY+F0ZEkRCg/1+oCb/N/4YtyeQnuOl8qdfx2xspuLJ5zoSWdvbe6eQxCpQGndbXfIQ7/a6cn1sVOgwqHGE0HAm1KqcwEmDEaSq25WqUAm0pN+ozlMEE/ZTE5yI/5Nr/O976Pm5jP7h+Qg9yc/SFkB4CgZUAAGMc89hnHvOgJ0/hkTPnCDKKVE+obLEQSM5kyey/bCDyF/1BfVPdzLdGR7Vyj2VUkO6qYRcoIfCpWaO0oacIIVhiyOEhitWVN28PWk3wRZlFkml7UJvzHDSUv6GfS5VBTQTMPKbJ1N14w8xinqo7Rbxq+N5ihhoE1zIaiVkpVi+JssoQVScGR/R6L1g/reT+gcb9p3D+nFVbL/vP4S321UOBKqVe1uK0XKm3V01mdASml2922l+N1wZvrGsOzv+BlRlhCRfcSwnqKcSPjHa6tTxdjDtvbeung/Wdi23sqnBz4zoChbK93iu5OsE9I4E2TkTSphXpRrWNb2xFO9u43GNKuu5MnbYr9ptp5Ic2Q980PIA0H95Qj0R8xHlBusz4yNq2gCL/6oqK+SN6rpNCBbvfyD77+rFHNGp+aBu2lGWrT0Xi42V6vHVJs8dMnOUz9JdkFrkpsOQwhFCw5FwQDmJXT1URmirV6vMVNouxMxwcTSreVUdQiXGS+8uY+GaZ3GXVfHNg6fFLRPk/2wtW/50PwUL9qbyh+clP384oCagnIEXQBvq29hY76MxT030izvVNlz8Zfc6h5mu1JAMiYHPVUJOqC69Fhzx0D3QuoX6F//AGxMuJyrM9ioVzYEIoHHN+iisX85ukVYKdxkLmsYpe1Uqs5zp6fn8plvlrwVblYabCMOjIiuLxveuTYjDoMURQsMNy1JmOMOdPBgh1KqKYqaysmw3w6WoMUXDzNl0D1II2O/yuAIo0tzGpuvuwSwtYtTFpyY/XiSgJiFv/Hbf/c24klzGleSyuEmZq/Yblby2YX8JoBiWUIIIafVNEOWNgH2+R8kb17PntqdZPurU9ioVn21tZnODnzkTSnFv287E2x6m5oD9qD1o/65mufyRPSc+d84dSqS56yaEQkrDzy3t3fU4DEqcJcVwI9CkasQly/exIrYZLsW2C+3RcCn6YD56kJLgJj6qPAfyuk/Q0rLY8uf7iTa2UHnlheh5SQRhJKjOnVvq1BJLA0uYSohYkb41ixs7j9XFB7BrwxtUNC+Pu0twRDlN06ZQ/r+3yF1jN8JLJ1pOs5NTe6qgbeYo052TOzSscITQcCKWE9STdhNLFEyl7UK60XCRIGz5kC+LF7C1IH6bg8aX3qZtyaeM+OZJeHcdl+RYIeUzyCt3BFBv0G1BFA32qUvrJyNPoN5dSVXj211eL/DYWosQVB97BMHyMioff6a9AWF7tFwqSaymR/3WQkmElhC2idDJHRpOOEJouCAltG1XZpAEYdCAMsGl2icoElKTSKpmOFDa1dE3sHLUiQl3KTxoDqMuOYPioxckPk40DEjIK0t+Pf2JFYVICN0KYlh+Jrn2ZlfXLAzLh2H5MSw/uhVEs8KZa7nQVwy3apcRDfR64o4Kk7fHfYf3xn6zy+v53o6FiXS52HT6SYhwhM3X3YOMRDvMcqnWljO9Suu2koRjGy67Gnxjr67FYfDh+ISGC8HmjuixRMio3dY7BQHUXhvOlVoktBWFVU/D7seA4UaK7ubASFMrmstE87opPmr/xMeKRtRElD8qeXRfNpEoQWhF7CcoE6fLS8DIB6FT7DoIAB8gkCAlmgyjWREMK4SQtgYQ8qvJM52Gf5nE5QVK1Pdp5PQqsl1VpQAz2sKo1s/4jN277RMqL6X6+CPZy2OBZp9EN5X5LJXacpqutG5/PeQmqRnnss1yrtzUUgscBjWOEBoORMN2TlAPZjhfg/qbymQYbFYTcKotGlY8Bisfg/wKGLdPt80yarH5hnuw2vxU/fFHiRNSLQusEOSNHJh23JGwmghBrcy9hUqb0Ix2jSxsqAz+Nkv5MHK1+AETgiiaFVEh5bHEYaEpYdbfAsmVpxYKgaaefydJ2L32v+xW/xpf5l3AWnNSt+3N06ZQNqsSABmNInS9o0yP4enZD2l4lB/JzEm8oBKa8hO2boPCMY6pdogzSOwcDr1GSnUzCiODZrhgei0atq1QAmjCgXEFEEDtIy/i+/gLio9ekFgASanMRjllqQdNZAIplbYSalO18/JGQPE4pYm585XQiPPZvtd0F+81JS4ML9GJam4lyAorVaO2nBIlDEK+1HwlmcRTCK78Pjn2V5UfTZNrJMe1Pkqulbjage+TNaz5zu8IVdd0iparTa0TayxaLqlZzq2iJh2z3JDHEUJDnWBzz32CrEgaZjjL7pSaohku0ARv3gT5o2HORXF3af1wFbUPv0jhwXMpPHR+gvOiqiF4ivsvITEmfMIBW1CMVYLHlZu4zUBf0EyVcFk0BgoqlCkp1KbO319+dm+RPYEHe/X2qObi/bFfx0WIE1ofQxDf92WUFRFtaWPTtXdjhcIdhXH9TXH370JMS0w1Wq6/hblDRnGE0FAmElJ9glzx22O346+3V6MpTKz+JjU5p1JNG1TF5WArLPhhXEEYrmtky4334x47ilHfOT1xf6CIHzwFPWfZZwKJ8p+FA8pUVjxW/e03859QC4K8kUrwuXIg3Kb68mT91MIOd9d7PXm3uEbzsvdYJkS/Ype6RXH3cY0sZfQV5xL8ahPbbn9cvWi4IdSiFk09YdpmuWCS2nKdzXJOtNyQxfEJDVViZrgEpqJ2gi3qZu5JUIGaHEIt6fkMZp0LVftDcVX87VELd1UFo759OpongYktElAReJ6i1M/bW8IBpe15i5TQy4bGkw66CbnlylTma7BNgu7UFwG9QWgq6rBlW6+TWZe7Z+OWfoJFcxPukz9nKqWnHkbdv18hZ8pECg+eq0KsffV2oEYP12h4lQnP9ChhE3cft/p9+xsHRTUNh/RxNKGhir9RreaT+U6ssLqJUwmxtiLKDKen6Afy2aaSvFEJ/UBIiTmihPG/+z7usaPi7xMJKX9WlrPghYyqCd7wKGe2t3jgBVBndJeqLpBfoYRkyJ/d1b1mKOFnhXuXzCrgA+8CQno+ggi65Y+7W/k5x5AzdRJtH39unzdWkLQHU1v7vrr6DSezV7pyVHpCL02MDgOLowkNRSJBOyk1mXZjt2gQemq1tnz1qe8bbIIXfgrj94HZ34i7S/4nn1Hw8Uqie1yEnpNAsLXnAmUxGVVKDBlEIpS/pw+RYTuye85hGTtWO6ZXBTEEWsBvV6rIVhiy4VI5RL4aEN5efgeS/dbdSlQzeWfcxezoSBS6zthffhvh6XQNhkstoAItPZtfY7Xl/M3KbxcPoanFU4sdLefUlhtSON/WUMOyoHmrMtkkM8P5m5XpyUwhyizQogRbKpOdtFQgQqAZJiyMu0tw41ZGPfkcRmsbwkygbVhRtQrPG5G1cGUhoxjST0jPodUzMqMCCKDSPZNK98yMHhNQ36u30A4/1tQknC2lyOVVwSARfy/PIagumM7otk+ZVPda3D00rxshBKEtNdT8026El1Y1BS/4a5NrOoarozmjw5Aia0JICHG3EGK7EGJlp9dKhBCvCCFW238dI266+OrsMjpJBEYkqG7aVDqlRoLpFSdd/hhsWw5zvwmlE7ttjvr8bPrdHUjTYMlRR6KZcWz5UqpSMrkjshYMoFtBNBnG5xpBwCxGknlB1xLdSkt0a8aP247ugsLRqnBruA0tlr+UaTz5yj8WiW9S64k1JQewOW8aU7c/TbF/XcL9Wt75mNpHXqTh+TeV1qW5VHRbT2HbQqhFV+t2lXCdCDNHafSpBD44DBqyqQndCxy5w2tXAv+VUu4K/Nd+7pAqwVZVRTjZil5G7YCFHqpog9JG2mpTD8fe/BGsfFTlA+3a3RQlpaT6pn8Sqq5l8+kn8VUkjrU3FortLUk9ETYdpIVh+YhqLlo9I4lo2cs3+qD5AT5ofiBrx1cIpRUVVCKQ6FaW/B6eIrUQ6ZVfRbCk8mz8RjFzN92jfG9xKDnpYHJnT2H7nU/iX71BLUCktHtf9YBuqt+2L8m+QqiFV0u104l1CJE1ISSlfAPYUTc+AbjP/v8+4MRsnX/YEQ1Dy1Z1kyUTLrE6XaloGL4GJXxSjcTSdRg5FeZeHHdzpKYB36q1jPj6CfgnJChMGgvFdvfQ66gXCCuCKf34jSJ8rlLkcHJ5Gm5aPSOIai5My5edoIXcEuUXjKYfuh3Rcvmg8nwimqlMtXEQmsboy89DL85n87V3EW1ps8O22xIKri6YXjsvLkmRU80AhFPkdAjR33fpSCllNYCUsloIkbBAlBDiIuAigHHjklRa3hlor4qgJw9rTSccO9ACEV96fpJR09VjBx7/cHP7/9vPPJP1fknzp8pMdeXjqvz/6GIv4woMoprJKQvGpn7OFFGdRC1aXSNVlYJhiETH5yrFjLSSE2kgggeZyQi/LqHbkbTDxOtzJvLaLj/j5IKKhPsYBbmMufIbrPvpTdQ89AKjLjpVaWC+OmV+7GnxZHhV4nWhO/G+plcJKidse0gwaAMTpJS3SylnSylnl5cnbxo27PHVq5DdZOaraAjaalITKulWx15yDyx/NOFms7GJskWLwbIYUVHKnIllVBarY8+ZUMqcCaWMKzCxEATMHopY9gLD8mMJjTb38BVAHQjCRj5trpHoMoyeaT9Re+h2qFftHyS60trf+wds+SjuPt7dqhh39XcYcd7x6gUhVBRgW03P/iFNU2Ns3Z58XzNXHS/cOz+XQ//R30JomxCiAsD+u72fzz/0CPmU3yaZdiOjavWquVL0A9Wkti/A+rfg82cTmkusYIgxDz9B6dvvYzbHz24XVgSwCLhKMhsgICWm5SOs5dDmLscSw8j81gMRzU2rewSSmBaYQQyXqt8X7WU5IWlBzWew+M9Kg49D7ozd0Twuoj4/bStW2z4fKzX/kOFSgi5ZrpGwq1I0O/6hwU5/C6H/AOfb/58PPN3P5x9aRMPQvEUl4yUUGNL2A1lgpOgHgtR8Ro3r4e2/Q+lusNe53c8sJVv/9jDe6m1sPuVYwkVdtZwCj4GQUTTC+F2lmRUSloUh/fiNQvyuYgZCqZ+SewxTco/p9/PGsIRJm6ecqGZgJEgW7TWuHDt0uxeRZoYbDvyp0qTeuCGpj2nbnU+w8epbCazZqIRGqC25zydGzD+UrKyPbgLSKeszyMna0lEI8RCwECgTQmwCfgVcCzwqhLgQ2ACclq3zD3ksS63ievID+ZvT9AP5UwvdDrbAomuVCfCAH8f1D9Q/+V+aXv+A7QftT+vuu3bbnu8x0WSIgFlMVGQwFNuKQsSP3ywhrKdw3VlilGvKgJ07RsxP5Ak34Yq2Eul10mkcPPlgRdhe/xXrm+Kbvj5YW9fttTElXuZVjYJ9vw9vXAcf3A7zL437/hHnHkfbR5+x8Xd3MOHPP8YoyOsooJssb03QUdbHcCdOWTBz1G/ZV5/1qhwOvSNrQkhKeVaCTYdk65zDirZatYJMJlwiAZVVn4ofKBxQiXxmipP29s9UMdNDf6naD+x4uPomah58nvz99+LTA/eLewi3DBA0SohoaXRm7YloWD0KKgnrWzJ33F7QENkAQLEx0IEzGgGzCEvoeCONhPFmrhutt4iRJYVUFEWIaqlVbphXZf9exs6FPU+Bz56Fqaepgq07YBQXMOb/vsX6n97E5uvuYdxvvovQ3MpknD8qeSKzpiltp2WrqjKRKEjDZQs2w52VqEyHvjFoAxN2agJNEGhMLoCssPIDpZIPFI3YK0ZP6l01x86Bk2+F8j3ibjZLChl/7Q8Yfdk5cc9vWkH8mpewnsGbPhJSWlBhZeq9jrLIRy2P8FHLIwM9DBtByCjAZ5ZiSn/m2osLgd8sBETvkmVnnAnH/imuAIrhnTSOikvPwrdiNTUPPtehdafSFrzd5JakvpwQaqHWUu20fRiEOEJosBH2q5Vd0kAEC1q2A6LnMFpp2QmpWmrlcda/BRvfU/+7u0eyRZrbaF68DADvruPQPN1Xx7oVIiJMfCKfXvWSjkc4oI5VODp5tYidnLCeS5trBLoMZswhLzHwuYoRWHaQSRoITRW5BfjiJWiNX2Gi8KA5jPz2aRQfvUC9YLiVwAim0H/I8Kimjcka3Gm6MvE1b+ldwVaHrOEIocFELBDBSGZOsQMRouHU6sL5G1T0XCoTd91XsPivsOrZ+GeORNl87V1sufE+wrXxo5hUyLCgRStEZsokFA4ov1jBqMQl/R3aiWge2lwjlDk3mpkQbkuY+M0SNMK907ICTbDsn8rPmCBsuuSYAzDLi5FRi+D6aru+XHOKiaw56reeLKjBcHVEkjqBCoOGnSemdbBjRTsCEZJFrvntiKA4mtKbX9ay+Mva9ueeSAsN9XUEhZsxxd39MnMmlHTY7/2NaoJwF8ABP+q2r5SSrbc9im/FakZffi5mWfckwFgots9VihQZqt8Vy4/KKx9crRcGOVHNrcyWzdXKSpVK5GSPx3QRMIvwhBuJ4kkvAMJTCAt+BK/9BhbfrKLnEry/5p/P0vDcGzRf9g3eDxkY1le0GKVENJP1dep3Nb60ux90v12KWTAeKDATtzgxc1Qgj68OcstSH79D1nCE0GBASpV8Z4WTBxmE2lTBR1d8P8uCSWUsmGTfWCEf+Gp5coXKBzppr9GJjxsJwuu/Vw3tDv+dmjB2oO6xl2l86W1KTztcNSfbgfZQbLMMK1ORcCE7ki+vPHOO9p0J3QUFo5V2HZEZaQkR0XIIGlE8kRYViZcOFTNU648ld8PHD8HMs+PuVnzsgTS9/gFFdz7MRTf+ELMwV2kw+SP51wfVAJw9L0HVjVh5q8LRibVmV66KltPNuL91h/7FubMHA746tTpLJoAiQSWozJye3SyRkBJWeoqr1XVvQcMa2P/yuJWxA2s3U/PAsxQcOJvyc4/t/n5pdYRipxhB1SMhH7hzIH/EoBVA0/JOYlreSQM9jOToppqQIWPtw8N6HiE9F90KpP/m3Y+BXQ6BTx5XwjEOZmkhY39xEdGWNjZdcztW2FK/49bank2BugkI5TNNtK8QShC1bHMqbg8CBufdvTPha1Crsh4j4baqG6ynhl1WpKMiQqrNvSYdAkf9UYXUxsEzoZIx//ctKi47G7GjUJOgyyABoyAzodgSWwDlq/IxmQpsyALl5iTKzUkDPYye0UwoqFDO+XAvBEc3BEEjn6jm6l21hjkXwSG/UlpaAjy7jKXyxxcQ+HIj1X/5l9LqZBR3pKlnf47pVvdB0og5TWnZzZudjqwDjCOEBpJgi9JuXLmJNRYZVas6ofVc5UBaqrijSCFqDmDtG9CwTv1fMqH78DZU4/9Cbc+fPz1+b6CIn5CeRyQTodgSCLep1gW5pQxmAQRQE/6SmvCXAz2M1NAM1T5cMzIkiOzcJET69et0o6MQbvXH0LQh7m7586Yx6rtnUHzcgeoFw40RDeCOxq/S3QXTjphrqyOhINIMJdyaNmcsgMMhfRwhNFCE/SpvwZWb2NwkLbtQY6Rne76023ljpRYJV70M3r4ZVjwWd3OkoZmNv76NzTfch4wkMGtEAmB6CRk9tGhOhXYBVATeYgalAJIWREPoVgg9GuCL5n/zRfNjGFE/RjSAHg2gWXb02GCMvmqPMDTUd9dHJDoBVwkge5dDFA3D23+D//4uYc244iP2I2cPtUAKrNlIVPdiRtpU9Y+eiFVL8CcRWrpL3X9NW5wacwOEI4QGgnAAmjapUOyEuTtS+XXCwdQSM/12q2Q9hbDt+rWqplfhONjnu902R1t9bPjV34k0tVL5068jjDhjjASVmScnAxpLuwAqHjwCSEp1jbFeN6FWFfIsdMKah5CZT5Nm0qi58ZvF+M1CQmae7ROTGDKIabVhRu3396JHT1YQuhJEIjMakSUM/K4SlUOUbui2bsLCnypB8frvko6n5f0VrP3B9RR98BFRzaO6AfcUui2wQ7frVJ25RBhuZXFwcogGBCc6rr+JhJQdWnclqQknlRkh2KYKSfZEoFH5UVKpCdeyVYXJGrlw0P91C4awAkE2/uY2ghu3MvYXF+OdFKckTTSkTH65ZX2vU9ZNAA0gVlRpCFKCJsDMU5+/4VYC1/ax+d3Kmd1ma6chM74mKGQEzYooQd25MKcw1DEzVeMtXWKCqHW7mvj72OE2aucQecO1WKQZmFK6iwrd/t/v4K0/woE/i7swy5s1hdzZUxj57EtE8nJhrwPVIk3oyccfE0RttUrjSRBZiulV5ruWrZBfkbo/1aHPOJ90fxIJKQ1IGMlNZv4GtXJLSQC1qIS+VHsDffKEMisd9qu4eRL1z/wP/+frqPzRBeTNmtz9/dGwmqTzRqRWgSEZg0EASakm4lCbCgDJKYWisVCyi5qoPQW2EEr/VpHCIKp7VGO1ojH2MSvV9xr2dWhXA4HQ1XdouDKiEakcomJVqSFdU+SYWapb75aPYO2b8Ydr6Iz56TcIVFYw+rGnaFv5pYr+bKvpObBACLsZ3rbk/YXMnI6KJb3opeTQOxxNqL+IhpUGJLTk/h1/g0ocTaXQaLBVmSWMFMK2Y8y5CKacoCbDOJSefAg5kyeSMzVO1JcVVVFH+aOSV/ZOhXYBVKL8QP1NTOtBqgRdT4FdWy+L2ommKQHkylGRf2GfveBoUZ9nts+/IzFBlCGNKKJ58RuFqlJ7Or9JgF0PV/XlKmYk3EXzuNn0tdMYd9c/2fTb29nl9l9i5HvV+PNHJQ/c0TQliFqqVVReIhO3K1ctSNq2q/EMlLa6E+EIof4gGlYaECJxJjeoqtX+BiWAevrth9qUH8jw9rivsELw/u2qmKS7oJsAklJS+/CLFB2+D2ZpUQIBZEE0qG7MHW72DfVtbKyPn2+x+Muabq+NL/WyYHzOwAggK6pWu5qmtB53ftptrGPMK/x678ehaaqisztPreT9TapOmrAny/6a/LoIohTbfCQhouepBNBAU/qCKCaAWrbCxvdhyvHddonm5rDxgrM4wGrGKIqZQe1k77yRyb9LTbcb3W3pWRAFbdOpI4iyjiOEsk06AihWDaFHAeSzS9P3PFkJGWbexrug7VOomAVjZ3fZLqVk251P0PCfRWheN6UnHtz9IFKqLps55XGvYVxJLuNKUmwRIaUtgEpVKHZ/IS0luIWmQpVdeX22+5ea3cPae4XhVkm5OcXKtOqrsyfMfhJGGRZEeArVoiXUklqbkR1Z/TKselp9Z3ue2G1zpCCfwr12B8C38ku0XC+ecSOUuS0VQaS7UxNEsQg8RxBlFccnlE1iPqCMC6AUqyFEI8zZeA8bmiMw96L4AuiOx2n4zyKKj19IyQkHdT+GlMq8klMKrj5OTnY77n4VQFJiWn5MK6jK/xRXKdNbBhzPm4Mfszn4cd/HGEM3VX5UyQRw5SufUX8lUgpdTbaGO7nfJFVyitVE3ht/017nwLh9YdkDSiDtwKZGpXXLaJTqvz7Ehl/8leCWOrXAaN3Wc6i1bnQIomSh6u48pRG1bnN8RFnEEULZIhKEpo2AlkQASdsHVJe+AEqlcsKbN1LZuoJXXIfCrkd0PbOUbPvHYzQ88z9KTjiIkd88OW41BCJ+1eY5lc6tybAFkN8s6jcBpFshTOkjYOTT5K1UK/QMRj193PI4H7c8nrHjtaObSjMqHq8mzGCL+j6zjdAyK4i8JUrTSDcnSWiw72Uwei94/x+wZlGXzZsb1NiErjPmFxeDEGy46i+EtjXY5X1SFESGRxV4TXatMdOcE6yQNRwhlA3CAWjcaIfiJgpCkF2DEFLxAaUqgEDdOI3rWTbyFFpdI7pttvwBfJ+soeSkQxhx4UkJBJBPCSBPH5NRpcSUfvxmESGjoG/HSgEho5jRNizNoMU9moBZjByKFbgNNxSOUdWwo2GlRWY7CTaTgkgIpdnprvQ1Ot2AA34CI6fCmlcTXre7cgTjr/keMmKx7sqbCWyqTV0Qabq6zpbq5DlHrlx1L7Q4eUTZwPEJZZqQz05E9SSO1pFW16KlPQqgVrvqb3IB9PiHm9FkBCkEmxqDVEcvommbBYS58vHlAIwu8lBZ4AFd58TrLkfzupMIoMLMCSCjsF8EkGEFQVq0ucoJ6znDw5bvyoXi8QSMGjzhJqLChZXNvkoxQRTzEfXpWHY+WVutCmxJJZna5vGPa9CLLgCibHj9S7Y0+GgOKCHQ/nsu9lJZlIPr3DPZ46FHaXzxbUZdfKryxbZuU76upFFzekf4dm65ClSJh5mrrAJNm1Udvp5KaDmkjCOEMkmgWantpjdxCLOMqsKKYX8aeUANKTmpNRlm9ub7iQoXsvIcKoty+GxrM5sb/MyZUApRi9H/eQEtEGDT6Seh58RxyrYLoIK+l7mXFqYM4DOKCGeitE8ShIxiWgGCeh5+VxFSDLOftqYTMIsJ6zl4g3WYUR9hLYuBC0KzgxVqMCx/34rTCk0JopbtaQuiWFX2cQUhzmh4kGWe3Xk6sJf6PXciNKKMqj/9uD1iTqKrj6ZlW0c+VCI0zU5orVEmt0TmYsPbYeUorEzu53VImQExxwkh1gkhVgghlgkhlgzEGDKKtCscxNpyJxJAVgSatyobeSpO/kCTHYadwoo+EmT+xjsZ2/IxjZ5KdlSvRCjM2Icfp2jZCgKjRqqKAN2ug04CqKjn8SWjHwWQEQ2gWyFaXCPwucuGnwDqRFRz0+oZRcAswGW19a5mW6oIDfLKiWgeDKuvGpGm/FzovQq2kEIQNnM5IvA8pbJ72D+AWVqE0HXCdU2s/cF1tH78pboXW7f2HCAhhNJ2/HXqnktU9NT0KO2pcYPTBiJDDKRP6CAp5Uwp5eyedx3EWJZS5f31dnBBgo80GuqoTdVTUqCUqsVDoCl1f9F/f82otlUsHXUGa0q7hlmXE2H8/Q+Tt3oN1cceQe1B+3cXahkWQEZ/CCBp4Yq2EdHdNHtHEzF6EQ7cB/Yrupj9ii7u13MCIDQCZjHN7goEFmZfBUQP5/K5SjIriISuNKI0sITJB5UXsNKcyrHWIiZvf5ZEgkLoGgjBpt/eTvO7K5Tm1bYtxVpzuSpaNVk/It2lzO1Nm9Q96tAnnMCEvhDLAQq12gIogbQI+5QtWWiq10kypKXs56G21PxFoLqi1q7mvcoLWFe83w7Hkxy66FU81VvZdPpJNMzZK845ybgA8hvF2RVAkRCm9NNmltLmKh8Q7afQqKTQiF95oj+I6h5aPBUE9Vxckdb0C4imjBJEYS0nsxpR2oLI4D+5p+P25jKl7mWmbH8u7n5GUT7jf/99PLuOZfN191D31CKk7lV+2J6qbwuUmTwSUFaLRJqmZnQ0xmutcSLn+sBA2S0k8LIQQgL/kFLevuMOQoiLgIsAxo2LU0RzoAn77c6QWpISO1KVkffXqZVTT7XWohElgGQK2lJnZpwOlsXmrd2j4BCCFfvsw4QiD76qOJ9jlyCEvvmAhIwqDcgsJaz3MaQ7GeE2ECYt7tGZ6+TaCzYEPgBgnGfOgI1BCh2/q5SI5iE3VIcl9Cx9Jhp+VzGEwLTaiGh90DpjgqjVrvuWhm9FCo3nck7BKBjN5sKZCffT83IY99tL2fLnf7L9nqcQhk7JcQuVf9WKqEodyUzcpsfO89uikpvjJbXGCqIGmuxqIqN6XX1jZ2agNKH9pJSzgKOA7wohDthxBynl7VLK2VLK2eXl5f0/wkTETGWNG5VankhYyKhydPrrlJDqSQBFgsqnhJXaTbltJXzxovp/1AyVU9GJvC/WUPrGOwA0jBqVQADJTmHYQ0AAWVEV/OEqgKKxAyqAAFa2PsPK1mcGdAwxwkYuzZ4KLIQK5c4KAr+rmJCeh9HXcHHb34Rupp/QKuDz8iNodY0EJLvW/hdXpLuGo7ldVP7kAkZ8/QQKFs7pqKgd8tlBCD1ojoZd6b55s/rdxTP/xVqFR0LQuD4z+VU7GQMihKSUW+y/24Engfh9pQcb0YhKbmurSR6AEA1BUzWE/HbTuh6OG2pVztNYp8ee+Op/8Oqv4fMXuneElJKSxe8x9l+PUfDpZ2iRBLkS0upIRO1jGLaQUXQZxGeWZU8ARUNqvAUVahXtlNrvhqWZtHpGEdTzcEVas5TTIgiYRQT1PEzp77sgyi1XWkYvK3nnhbYzueY5Dlr7R2hcF+cUGqUnH4pRkIsVDlP990cIN/rUfdOyVQmPZOiGWkT66tR9n8jkaXpVu4/GjWqROhibGg5S+v1OFkLkCiHyY/8DhwMr+3scaRP2q4iYiF/lEiQKQAi2KP8PpB6A4LMj4HpS5aUFH94H79wCIybDEb/rkq8gwhEqnn6BUS+/Tsseu7HuG1/DMrofM6a1kFPa9zygaEQJINcIlZeTDcI+de1F45TfyiExQsPvKqXVXa603Ky0ihAEzSICRgFGnwWRgLwyu72FP2FQWiJaXSN5s+oydBmFF38OGxMH24Y2bqN50RLWXnEDbas22GV+tqYQsCDU+EI+ZZ5LFN2nm8o856tV2pPTMjwlBmI5ORJ4SwjxMfA+8JyU8sUBGEdqWJYKv27cqDSVRAUZZVSVf2+rUSs7o4dktmhEOTXDsQCEHtQlKWHR72HVf1QJnkN+2aVBl7Qsxt//MMUfLafmgH3ZdPpJSFd3rUpYETQZwmeW9r0UTzQM0RBtrhFEtCzkTEiptETDqwSQk5eRMmEjD4rGd2i8GUcQNAoJGIW2RtRHx3xOiVoQRdI38zV4x/PahB+pgqRv/AFWxTeReiaOoerGH6Lnetlw1V+pf24xUvcoLScV7cW087KaN6sIuoTmuTylYTWsV4tSh6T0uxdNSvkVkLhpyGAiElSCIhrqIfrNr7LLIbWJPdRmV0o2UmvdDercFTNhzHzY9dDumzWNhr1nULfvXFom7xb3EJoVRmDhd5UTFX3M+I6EleAtGE1U29y3Y8XDiqjPNbe8ZyeyQ3wMtxLeLdvUZJjsN9xLQkYBEkFOpIEInr6VR/IUAboKHtDT620UMIvg8N/Be7dB/uiE+7nHjqLqTz9iy58eYNsdjxNpbGHEecepxWAkoJJqk1VD0E113/rrlHDPLVVmuB0xvXZeYDV4fOq4fW0COUxxQjniYVkq4qVtu7oZEgkWaanW2v5GlYuQijnN36CEUKo14L54SVWdHjsb9ji26+EiUWr+9Tzu8RUUHjibppnTEh5Gt0KAwOcqxeqrAAoHlZ+rYHR2ypdEQmCFVFZ6X7W1LHJA8fcHegg9o+nKj+arVwsfMyfjk2HYyMcndLzhOqL0UVv15Kvx+WoRVgSZTvNEww37Xdbx/MtXVTfbkq4tN/QcL2N+/k3qHn+VvDlT7fd6bD9RtSq86k7QBhw6ghHCAWjcpASMO49uzl/NUKb7UKt65I8a1L/ngcIRQjsS9quVoxVR5fSTaT9tNUqwpJJQGg6oSQBS67ES9sG7t8KGt2H8/t3aMIS21LD5xvsIfLGe4uMOpPDAxDm/phXE0vLwm4XIvn7lkYC6ufJHxF8B9pWwTxV+LRqfvNTKICBP794efVASKyRquNTK3HCnFgCTBmE9Byl0ckI1ytTcl1BlVw5oo9DYCpYk2pvfWSQIyx9Ti75Z56oFXKd7WWgaZacd3v68+u+P4K4cQfFxByL89ep3mFOSvIOw6bHN9TXKX5RbkkArylHzSdNmJZR60rZ2MhwhFCMaUavFQKO6SROtWKyIqo4QbFXaT09+Cmkp+3GoBTRXajfn9k9h8c3KwTnja7DnSR2Hk5Km/77Httv/DbpO5ZUXUrDfzITndlsB/JoXn1lMn12AIb9Kts2zs94ziZTKJOLKU8cfAqaLr/yLAZjo3a+HPQcJ7nwodqn8trCvdw3nkhDR3LS5R6h7RFp9W0QYLnyuMjyhenQrSDRdn6PhhmNugHf+Bh/eC1s+gn2/H7eTrwxHiNQ10vjCW7Qu+ZSK75+NWWwogZ1Tklx70TRbKwoqrchbAt4CEmpFYT80rIOcsoy3FhmqOEIoZnrz1amVUqIqukhlV2+r71DHe6Jd+5Gptzqu+wpe+QXkjlQ27vI9umz2r1pL9c0PkjN1EqOvOA+zvLh92+ZGH1vsXisaUVxEaCYPv+Zly5d13U41viyHBZNSWM1L1KTlzlU3T6LIwN5iRe2VZ5lqhjZE/D+ftb0EDCEhBHZ7iLEqKizYaqcQZO7zjgqXMtO22PXa0km63oF1DUE21UnyZRseagniJnYTfbC2++95TImXeVUlHS+4C2Dhz1Q+3ZJ74fkfwQl/67ZwFKbBmKsuovHFxWy7+0m+uvT3jPzWKRQeNAfhq1Pmc29J8gWk6VYLKX+9mifyyuL7e01vR1WUQJNacKVSyHgYs/MKISnVxNe6XWlBrpzEk2vYryLkrIi6qXq4ad/6opqPvliHywoS1UwsdDbZwmFMcffCpXMmlDBvlK5WRqUTYe63oWpB+w0sLYvgui14Jo4hZ8pExv7q2+TOmozYYRVVWZRDZVGO7f+R+M2SpAmdqQkg+3PyFturyAwLiGhIPQoqkiwAHDKKbkBBpR0VVm8LogwuLHRTfZ9tNXauXO8qcHe0jR+BEW3DG2kiKtwJgx+6CKDO7HYkjNwTaj7vEEA7tDEXQlB81P7k7rUHn193Lxtue4x3gjmEc3LQZZitDW20ilzKy0q7fVb7TSpT91IslDsaVtqmO0/l4e1oehOa2hYr++XKg9xSnl9Vz/Mrqrvs+sW2VnUJI7v7qI6eVsHR0yqSfIJDg51TCIX9dsmQgPohJopQi4bUTRr2KRt6TzeTtCDYxv4jIuw/sqtP48mPtgBw0l47RO6E2+Cjf8I7/4Pj/qRKf+x6WPvmwNrNbL31Ufyr17PL336Oa/QI8mbvmeD8El0GsDQ3AaMQq6/11KxoR8SQOwv5OWG/unGd8Ov+J9bnR3cpZ7zhzayfQjNUTyJfrd03KwW/aeLBEtHz8AkX3nA9yGj61TIKx6oHwOYP4e2bYa9zYZdDuiwqXaPKmHrjFYQ2bWP6uAqkZdG06ANeH1eBpkU4YUahWpAl0150Uz1CfnXt3mK1wNrRvxTbLxyAhvUcXVXE0VOmdfkebn51NQCXHbpretc7hNi5hFDYr1Z/IZ+a9BKtvKNhFSYabFWOxlTDrv2NSnMwvKndcOvehKX3qXPtfiy4izqG4PNT+68XqH/mf+h5OVRccgbmqCSaixXFkH6CRj4hPZ8+ayyRsNL88iu6rBgzQky7Mj1Ova2BxlOgJr1UK7ynQ6wigua26yd6++QDiWou2txleMNN6JafqOjZKhGX/JFKE3zvVvjyvzD3oi4RdELTcI9TGkbbss+p/vM/qRpRzvajDoW9xtuC1a0sA8kWT6ZHmbIDTWpu8JaAJ6+7P9X0gHQrv3GwUWlP3qKdJnhh+N/9UnYIn3BArfwSCp+Q+sEEm1WEVqp+H3+9uoF1d2o3mRWFV38FNaugeAIc+FMo61jpWMEQX13yeyL1TRQdsS8jzjsOPT/JWOwMbp9Znpl6auGAulEKKzN/I0hLhat6i1XFBscxO/CYdjJwc7WdPJ3JMGKhmsTpurI+xFb/vURi4DNLMKJteCLNSGmk32W2oFL5W9e8Dh/dDy/8CKaeCjPO6rZr3qzJjLnqIr7626OMu+8hNn62khFfPxH3KFO1cDFy7OtLMAaB7QeSKlLP36B++zsKIyFUoIiUav4JNICnGM0KZ7eL7iBg+Aohy1KrbV+9qnCruRLH/keCKiou1KZ+GGZeaiHX/kaVz6K5wOx58ndFmoHRKvJr5J4wcSHscjAIDSkl/s/WkjN5IprbRdkZR+CZNA7vrkkqiFsWRANKWHqK+p40GgtAMHOUYzXDEXDCimDIoFqJ9rVr6yDg4OIfDfQQModuQuEYlRsXaM68n8iVBwX2xN3HgIWYec6vufGGGtCtAFGRpjlXCJh0MIybC8v+pYJiQC0QrUgXDSd/3jTWikKK312C9tY7bLrmdib+7ecI01Q5bS3V6vNyFyQRRkKZ87sIo0I1J3UWMjsIo4LAJkJ6HoTHZlZLHUQMPyEUjajoFH+DyujX3V3K23Rga0j+RiWEhJ6i2c2vBJYVSV7GpzP+RqZVP0ptXQ1Unasi3uxVl7QsWt5ZRu0jLxH8ahPjr7ucnCkTKT5q/+THjATVDzWnLDPRNTH/T5YCEHQriJAWLe5Rw0IAAXj0YVbHTtOUH0d3q6CCZG3qe4PhVpFzsbyaVPtlJSAqTFrdZbgirbijLcqMnq6W5cpT5rgYX74KKx6DPU+ESYe1CyNpGiyfNp39zj2UcG0jQtOwAkFqHnqBkuMPwiwMQqhaXZO7IHF4emdhFGhW8487z36Pu+t+Zg5hPY+N2+2yYYZbhYybOcPKgjA8hJCUagL1NypTjxB2RYI4K3krDEGfEiTSsosO9jCJS1ur8jep/zUzNT9Jaw2sehpW/5cJMsS/tW8wP1f1/JFRi5a3l1H7yIsE11fjGl1OxQ/Owbvb+OTHtKJKs7O1n4zk00RC6rj5ozKeOwJgWj4iwoXPPQork5PaALPa9xoAu+Yc3MOeQwghVJi84VZ+Ii3Dlbg1Q/3O/I1qoZhKn63kByRkFBDVPKiFZUCNvbdh58Xj1fiW3gMrH4fJJ6joOmBzgx+juACjWC0+fJ+sof7p12l45n8UHbkfpScfgllkQHirEuSewsTaixAdPqOQXy2cddvPZHq7aKFf1EfsaDpb6xKauvfd+YM+oTsVhvaMEAkpE1q7ZmLaUTg7tq62lJAKNNsZ+Vpq/ptISNnIg62AUP6kVFcgVgRe/In6cU1cyKtifxq36EhvMQLl99l66yPoRQWM/uH5FCyYpdoSJ0JiVyvQIHdEZlRzidIGdROKRmW+AoJUraeDRj5+szjz+UUDzGrfImCYCaEYrhy77txWTMtPWMtkcIpQGrfhVm20U+k43ANRzaWCaIKtaj4Qeu8m6PI94PBrVML48kdh2QNQ/TGUfqvbrnl7T2GX235B7aMv0fDcmzQ89yYF+8+k4gfnoEm7koLQlJbjSlAuSWBfu1tpcq3b1YuefKWldS6qqrvUI1b+y1enBJan0G4tM/gTvOMx9IRQNNyhlUQC6oPX3XHCrKUyWYXalCCQMrVIt5hZKtCitCahpRbtFo3Axndg0xIezzkV0BhZfhYrQyNYt0VD317LzNWfsfg/T/DssScyuiSHCed/jVBpMafMHpv82JGg+uF5itSKKAPJhUJGlYD12g3tMiwgYv4fn1lCyBxmZqudBcMFhZWEtE24o62ENG9mfydmDhRV2u3s+26eQ9iTt+lVmlbYp+75NE10j3+4GSiEkm/RFFnHugaN2m1b0Ylw3YPP0OYqp6y0lMoiZTU45ftfo+z0I2h47g1Cm2vQTHU+/+pNeCZUIoJNSjCaXiVYEmhq//zAzhGSFtV1X7CloY16v0WIHM75x5tEhMGEslwmlCn3wkUHTFTaUet2QKpjewr6HIXY3wwJIaQLqaLWAs12WXotQZSbLXjCfrWvtGwh0kMop4wJrFb1w0Wz/T2pmNy2qiib1a9AsAlyR5BT0YjPLGG7a1cmrfmcvT9eSd5X64lqGi3TpjBvdB6Wx0OIHkxfUTtM2pVjl/jIzNelklmt7IRfA3o0gEDS4h5FNM1qyA6DDE3H5yolEnGRG64jLDzpFRXt8fimbZ5rBn+tXQi4jxq5bqjAmljwUNhvWzHS1xQKR1QxYwR8trWZxobtXM6/cIVC1NdX8hUL2FSgOhq7RpUx8sKTkbbmEmloZv3PbkHL8VBw4N4UHjwPT9UIRCw3zpWnhG48bU1oVJQVU1FWzKeb6qmub+GICTlEkYQMiMpwRxHidu3InsOatwB2RRdPAfoQKD4yJITQKG9YqbaaSxUV7YwVUauBUJt6SGkLnh7swtKyzXmtdldHqX6kqZbXARVi/fJV6v/Rs2C3o5Ajp+P/YCMAuWvXUfnkc4SKC/lk7lzeGT2Rffeq6vm40bDSwgyPCmPOkN031gE1rOUScBVkJf/HtPxEdA9trlJkX5NlHQYHQhAyC4hqLvKC25FWNP1abslPoCLFTI9a1Yf8mfltmh4wRyktK9AI4VDqaRRxaCOf53e/hnGN7zOx4U1mVz/MjOrHYfpd7RYWYc85emEeY/7vmzS++h6NLyym4Zn/4R5fwahLziRn8oQOC02seoJh+7R2mHwsoRMQLiJaDkJG8UaaIGKpQsT+EnWNuqtjzsPdIZAaN1GRP/hvwkE/QABNyI7cBWnZJrkAhFvtHBktNbNZNKJCmoNtyrmPpn6QhqdnweNvgE1LlcmtdDeYcQaU7gozz8MauTetX9TS+thKWj94jLK9ZlBz8AJaJ+3C2gvPwT+2ks+3tWD6e+i0GA0poaq7VZRSBqsI6FYQVcqnNCsdUIWMYlgBAmYhAaNoyNR/c0idqO6hxTOanHAdZtSn/ESZ/J4Nt8pNCzSBvyFzOTKuHCXUwj6lGUWtXmlGBR6DqOZhbckBrC1ZQIlvLcWBdcyMmfjfukk1s6yci6iYTt7sPcmbvSfRVh/Nb35I43/fwyhWlfnbVq7F/8V68udNwzVaImgmL7iNiOYloruJakZ7xfsir/orhU7ETpsQMtqRBIu0rzFHzWW6Yc8dwtGEMoUppPIBhds6tdYVdmRbEh+PFVGZ/xG/3TrYUjeN0FOPAvvkaVj/FjR8pZ57y1RzOdSPYvN/aml9/0ZkJIKW6yVv9p74qmwfj6HjHzem/VD53jg3VGzVIqNKiOaUZlT4CCuCTpiwlkPAVdD3Vg5x0K0QmozS6hpBxNh5ijEeXvJ/Az2EfsfSDFpdI/BEGvGEm4hofWxktyNCU35KMwfJVgzLR0R4+u6LihUdNnOUMAo0Kc1IM1I2/3W9fwX1OROpz5nIzNhLugHrFsOa10CYUDEVJh2OPnYuxUft3yXtom3Fl9Q9+hI1DzyLq3IEeXtPwZtbin/X8RiWajcuMfBaQQyPG0EUScfnLIXeEZwkUfNcyC6WLDS1TTMw9UwXe8w8Q0IIlbmjShPRjcTh1FJ2mLHCARVcIKO0azu6q2fzXP062P6J6iM//2L1ev2XSM0gXHkibQ2F+L6sJ7psA+N+o1Rvo6SQ4mMWkDdvGjmTJyIMnXc/7DlpVMgomgwDdnKaOy+z1QmkhSGDSHR8rvLstd8O+7CERqt7xLDP7N4RIxuf6VBACAJmMRHNQ16wBktoGTbPAYabNs8IzEgb3kgzUorMnCMmjGJN6QJNdsSs3vMc0RP7XApzL1aRdZuXwqYPoGkjjJ2rzG8f3gPlU2HkFEaceyzFR+9P63sraHl3OQ0vvMXIwgLW7KFylgpWfEo0x0O+4SZqusgNKiEUNjxYwoXVpdoCSgPDvv9iQincSJk30z1XMs+QEEII2bV4qLSUac2KKLNaOKiED0J9ISIFoRNj4wfwxfNQu9oOeoCINhJ9lg/hyqGuZgZ1/36VaMv7AOhF+eRO3w0ZjSJ0nVEXn5r6dchOQQHoBIwCorpb5WVkEN0KIpD4jULCRi597iMUj1j77ZxSWt2RYRd+nQqr2l4EYHLukQM8koEhontpbjfPtRHWcjJshtUIG/lEdS/ucDOm1YaFmbnFjulRj1jEbbAVpNU3M6BuQsUM9Zj9DXWfgAoY2PCeCmIC8JZhlk2ieN/TKD56AVYozHOLPlPbpGTU869g+PyMA5pLSrCqxtC6xyR8u1YBbYCEJrddrcXTYV7UDPUdGCZETbRBrwcNESHkEnZ2cSSo/CbSsrcINflpemJHZtgH9WtUw6nmzap0esM6OOxqKKoitLWG1iX1hPwTCTYZBLf7iDa3sctBPlyjcjDLS8ibPx3vbuPJ2XMXXGNGtjsfU0JKdBnGZQVVXppeRFTPISoMMl+VIIwgTEjPI2Tk972KdiLCPvW3cIzdAqN7b5edgbX+t4GdVwhBh3nOHWnGG25Q/YQyfQ5h4HeVELJy8YQbMSxf0pYOaaOboBeqfJ5IEKnVolsBACR634ReLJKwbFc47X5oXK+sLds+hbo1SmsBtI2LObDxIZr9FTS7K2i5YC7Bep2aL5sp2b6d8lVfEPV6adl9N0Q4zIQ77mfz7uNwjx+Fe8wIXKNKMUeWoLlsQaS7IRrB5ZjjMkOhO9IRSbKjhhONQMtmaKkFX42qfdVWC7sfCeWTia5fju/JWwi36YR9bsLBXMK+MkZWbSVnThWBUAXb3rLQcvy4x1WQN28S7nGj0DxK9S9YMIuCBbPSGm91QwvbGtvQsJAI/HgIiTzCGKxdFwJCXfZPublcIiIhsMJENYOgWZKViQDoaD4X637qVL92ABCCoFlIRPeQG6xRvlsj01oRRO3OraYVwBNqQshgZoWRXcXgsxYPm+oiuGUQD62YhBFIIugsWRtE7mBZ6NZML9nxi6vUY/djum7zFrJNG0lB6xZ2b13R7v357Z4/pm3qXsyKfMTE0Jfkr66jxZ9LqwnuFV8g3ljafoiR3zqFkuMXEt5WS+2jL+GpKqfI7Zjj4iKEOBK4GdCBO6WU1yZ9Q6ANa/mTWI11RJsaiTa3YE5ZgDnnGKLb1lF/yzVEghoRv0Y0oBMJmpSf8iGFp04mFCxi0xul6rymgTmiBHN0qao6AOTNmsKke36DUVqUnobTGctS5kApAYvKQg8jSouICg9RzaQnjafXAihsBzSYXsgbgS+bFTzCAZARldPhznei3xy6EdXctHgqwF3YUZIn4+0IBGHNS9jjwYz68YRtYUTmfvwdzfRiZ4yiWRGMaADTCgBR1D2tExV6agKoJ0bvxfKJak7SZZic0HbygzVMKVBBTrvXhJnYsAFvZAXCBPYHNDfR4+4ktGU7oXeexhv8Lyz+mHCNQcviNYTXFw+J27TfhZAQQgf+BhwGbAI+EEL8R0r5aaL3nLqukM9/tbjLa2VHb6R8DlhGIbWf5KPneVRdpzFFeIsLMXadDYB7l4lU3XAF5shS9MK8bt1INa8bzZuGw9Oy1GRsxX6IVkclBvuma8umv7pzNF2sflw2G8JZUeX7ceVC3pidpseJQ++QQoe8cvV7admqqsxnQSsCQVjPIax7Ma0ArlCTygfSjIzXU5PoRDWdqOYmSAGajKDJMEY0hGEFbfO07RoQWodfppdEhUmLu5IWd2X7a5+XH8Xn5UehyzCecD25oXoWjHOj5+Xg3a0Kr38UVFfDtk/J8Tew23FRZJEJn2fgA8gyA6EJzQW+lFJ+BSCEeBg4AUgohFYXWRxw5OHo+XloeV70XC+uyjIItWEUetnj0d8h9DhaZ6gNTQOvvcKIBR6khrQFjtX1tViZINPuRKkZXVseRMNo0UAa58FOlk1lLLbgi9WKigmEqMo/EjEnaKpEk+QtSWkHe6BaLzjaj0M6uHJUMVBfvdKKNMPOy0n9N5Tq7zmCScRVqjrFBptVxJvQ1H3Zi4CZVO5fC42Q5iGkeVSzulg6iBVUBZI7zxsx4SQEya5fj/hSGl9AyyfgyYfRo1XUHcDkY9UDAAkhHyLQCLd9P6VjDiRCdi6Q1x8nFOJU4Egp5Tft5+cC86SUl+6w30XARQALx2t7jMgVfWyWk5hQhAKXQXPn16RERiwIRaSMSIhKZNiivz8uACyJDEWRgYikLYxlJRiDJSnRBPWZOGfEQgYiyIaAjEasnvfPMGVAbb+fNbMMh2uADFyH10AUeoTm1hGayL6j3NQQ+W4hXDrC1NT5wlGKTJ3GbJ87hq6BIRC6QJi6wNAQhpJDfbr+eHNVMlbVRMuXb5cZsBdmj4HQhOJ9Cd2mVSnl7cDt2R8OCCGWSCln98e5solzHYOH4XAN4FzHYGO4XEdnBiK5YxPQuWz0GGDLAIzDwcHBwWGAGQgh9AGwqxBighDCBZwJ/GcAxuHg4ODgMMD0uzlOShkRQlwKvIQK0b5bSvlJf49jB/rF7NcPONcxeBgO1wDOdQw2hst1tNPvgQkODg4ODg4xdr6CXw4ODg4OgwZHCDk4ODg4DBjDWggJIY4UQnwuhPhSCHFlnO1CCHGLvX25EGJWp213CyG2CyFW9u+ou9Pb6xBCjBVCvC6EWCWE+EQIcVn/j77LOHt7HR4hxPtCiI/t6/h1/4++yzh7/buyt+tCiI+EEM/236i708f7Y50QYoUQYpkQYkn/jrzbOPtyHUVCiH8LIT6z75N9+nf07ePo7b2xu/0dxB7NQogf9PsF9AUp5bB8oIIe1gATARfwMTBlh32OBl5A5S7NB97rtO0AYBawcqheB1ABzLL/zwe+2PG9Q+Q6BJBn/28C7wHzh9p1dNp+BfAv4Nmh+Luyt60DygZq/Bm8jvuAb9r/u4CioXYNOxxnKzB+oL+XdB7DWRNqLw8kpQwBsfJAnTkBuF8q3gWKhBAVAFLKNyAz1Qf6SK+vQ0pZLaX8EEBK2QKsAioZGPpyHVJK2WrvE+veNVARNX36XQkhxgDHAHf256Dj0KfrGET0+jqEEAWoxeZdAFLKkJSysR/HHiNT38UhwBop5frsDzlzDGchVAls7PR8E90n4FT2GWgych1CiCpgL5QWMRD06TpsE9YyYDvwipRySF4HcBPwE1Rnw4Gkr9chgZeFEEuFKrE1UPTlOiYCNcA9tnn0TiFELv1PpuaqM4GHMj66LDOchVAq5YFSKiE0wPT5OoQQecDjwA+klCnXncowfboOKWVUSjkTVWFjrhBiamaHlzK9vg4hxLHAdinl0jjb+5u+/q72k1LOAo4CviuEOCCTg0uDvlyHgTK53yql3AvVsrSbP6YfyMQ97gKOBx7L4Lj6heEshFIpDzQUSgj16TqEECZKAD0opXwii+PsiYx8H7a5ZBEwUO1M+3Id+wHHCyHWoUwuBwsh/pm9oSalT9+HlDL2dzvwJMqkNBD05To2AZs6adX/Rgml/iYT98ZRwIdSym1ZGWE2GWinVLYeqFXOV8AEOpx9e+6wzzF0dfa9v8P2KgY+MKHX12E/vx+4aSh/H0A5tsMY8AJvAscOtevYYZ+FDGxgQl++j1wgv9P/b6Mq4w+p67C3vQnsbv9/NXDDULsGe/vDwNcH6vfUp+sf6AFk+cs9GhURtgb4P/u1bwPftv8XqAZ7a4AVwOxO730IqAbCqFXIhUPtOlD9FyWwHFhmP44egtcxHfjIvo6VwC+H6u+q0zEWMoBCqI/fx0R7ovwY+CT23qF2Hfa2mcAS+7f1FFA8BK8hB6gDCgfye+jtwynb4+Dg4OAwYAxnn5CDg4ODwyDHEUIODg4ODgOGI4QcHBwcHAYMRwg5ODg4OAwYjhBycHBwcBgwHCHk0GeEEFG7gu9KIcRjQoicJPsuFELs24dzFQkhLun0fLQQ4t89vKdKDIJq6Olg1zbrc5Vt+/N+Z4fXDCHEtk717PYRQtyR5BhVQoizOz2fJoS4t69jc3AARwg5ZAa/lHKmlHIqEELlNyRiIdBrIQQUAe1CSEq5RUp5ah+Ol3GEEEYGDnMFkFAwpMEbwBi7dmCMQ1FJ2NX28yOBF5McowpoF0JSyhX2McdlYHwOOzmOEHLING8Ck4QQJUKIp+zeJ+8KIabbE+G3gcttzWmBEKJcCPG4EOID+7EfgBDiaqF6Oi0SQnwlhPi+ffxrgV3s99/QWcux/39TCPGh/ehR2Akhfmyfd7mw+xTZx1klhLhDqP5FLwshvPa2XYQQL9qFO98UQuxhv36vEOJPQojXgevs/d61j/0bIUSrvd8DQogTOp3/QSHE8XGGdgq2YBBCXGB/ls8IIdYKIS4VQlxhF918VwhRkmhsUkoLVU/sjE7H3rHQ5SHAq0k+v2uBBfZnfrn92jP2cRwc+sZAZ8s6j6H/AFrtvwbwNPAd4C/Ar+zXDwaW2f9fDfyo03v/Bexv/z8OWNVpv7cBN1CGygg32aGUUufnqMxxj/3/rsCSHffZYdyHA7ejstE14FlUaf8qIALMtPd7FDjH/v+/wK72//OA1+z/77Xfr9vPnwXOsv//dqfP6EDgKfv/QmAtYOwwrgnA0k7PLwC+RPWEKgea6Mik/zOqMG2ysc0BPrL/d6MqkRfbz8uA13v4/BayQ3UHVB28Zwb6t+c8hv4jE2YDBwev3WYBlCZ0F6plxCkAUsrXhBClQojCOO89FJgiRHuR4AIhRL79/3NSyiAQFEJsB0b2MA4T+KsQYiYQBXbrYf/D7cdH9vM81OS7AVgrpYxd01Kgyq5Gvi/wWKfxujsd7zEpZdT+fx/gRPv/fwF/BJBS/k8I8TchxAjgZOBxKWVkh3FVoFoMdOZ1qXpCtQghmlCaCKgSLtOTjU1K+YEQIk8IsTswGXhXStnQ6TN42f4/nc9vOzA6yXYHh5RwhJBDJvBL1WahHdFpJuxEvBpRGrCPlNK/w/sBgp1eitLz7/VyYBswwz5uoIf9BfAHKeU/djh3VZxze+1jNu54rZ1o6+F8MR4AvoYyZ30jznY/4Nnhtc7jsTo9t1CfS09je9g+32S6muKOAv5k/5/O5+exx+ng0Cccn5BDtngDNdEihFgI1ErVy6gFZVaK8TJwaeyJvQpPxo7v70whUC2VH+RcVLvjZLwEfMPWIhBCVNoaSlzs8a8V4v/bu0OXyKIojuPfX/IPsJvEuuD4Dwh2MQgLKxi2y2AyGAw2YaNRbbttQbHYDMuKSUecIPoHGAQNukU4G857MIzzdEZ9vPL7xHn3zhxumMO95757tVi0l6QvFc1PKGaCvKyd7AHt4jsvB/S9IpcEhzZEbD+BJXJpdL9sQx4Oe1a0qRq/QWM+RR4ma/YhTkJWlw1gRlKHLGwvF58fAAvlxgRgpWwnqcvrO+uIiDvgj3I7+Fbf421gWdIJ+Sf56swkIo7IpbK/ki7I+2SqElzpG/BdUnmCdP81zKU2sCrplFxee+j53VvyqvXdirgegRtJk2/EMnRsEdEFnsg6UTkuLbJWVM5Qq8avAzxLOu/ZmDALHI4Yn9kLPkXbrAbKd6X+RURI+kpuUpjveXYBTEfEQ0X/BaAVEes1xrgOXEfErxH7jQHH5IaS/nqW2UhcEzKrR4ss8gu4p6j9SJoDdoAfVQkIICJ+SxqvM8CI2Hxn1wlgzQnIPoNnQmZm1hjXhMzMrDFOQmZm1hgnITMza4yTkJmZNcZJyMzMGvMfAwx2je22XjgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_stats(load_dfset(base_dir='example/VASP_3C-SiC_calculated/1x1x1/T_300K/', \n",
    "                      dfsetfn='DFSET.dat'), T=300, sqrN=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def monitor_stats(T=300, directory='phon', dfset='DFSET', plotchi2=False, sqrN=False, once=False):\n",
    "\n",
    "    prev_N = get_dfset_len(f'{directory}/{dfset}')-1\n",
    "\n",
    "    if get_dfset_len(f'{directory}/{dfset}') < 3:\n",
    "        print('Waiting for the first samples (>2).', end='')\n",
    "        sys.stdout.flush()\n",
    "        while get_dfset_len(f'{directory}/{dfset}') < 3:\n",
    "           sleep(15)\n",
    "           print('.', end='')\n",
    "           sys.stdout.flush()\n",
    "        print('done.', end='')\n",
    "    print('Calculating the plots.',)\n",
    "    sys.stdout.flush()\n",
    "    clear_output(wait=True)\n",
    "\n",
    "    while True :\n",
    "        N = get_dfset_len(f'{directory}/{dfset}')\n",
    "        if N > prev_N :\n",
    "            plot_stats(load_dfset(base_dir=directory, dfsetfn=dfset), T=T, plotchi2=plotchi2, sqrN=sqrN)\n",
    "            show()\n",
    "            if once:\n",
    "                break\n",
    "            clear_output(wait=True)\n",
    "            prev_N = N\n",
    "        else :\n",
    "            sleep(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEGCAYAAADMsSqUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACE1klEQVR4nO2dd5gcZf3AP++ULdfvcpdcKgktJKSRTg/Si/QuECyAIqhgA0WKoiJKVfkhTVAQBOmgSEvoCAkkEAgQAiE9uVzdu+0z7++Pd/buEu52Z/d2r2U+z7PPbZl55529nfm+3y6klHh4eHh4ePQFWl9PwMPDw8Nj+8UTQh4eHh4efYYnhDw8PDw8+gxPCHl4eHh49BmeEPLw8PDw6DOMvp6AG6qrq+XYsWP7ehoeHluxqmUVAGPLxvbpPDw8umPx4sVbpJQ1fT2PdAwIITR27FgWLVrU19Pw8NiKrz/zdQD+ethf+3gmHh5dI4T4oq/nkAnPHOfh4eHh0WcMCE3Iw6M/cu6Uc/t6Ch4eAx5PCHl45MieI/bs6yl4eAx4PCHk4ZEjHzV8BMBuVbv18Uz6L4lEgrVr1xKNRvt6Ktslzz333OSlS5eu6sMp2MCyZDL5rRkzZmzuagNPCHl45Mjv3vod4AUmpGPt2rWUlpYyduxYhBB9PZ3tDsuykpMmTdrSV8e3bVvU1dVN3Lhx4x3A0V1t4wUmeHh4FIxoNMqQIUM8AbSdommarKmpaQYmdbtNL87Hw8NjO8QTQNs3mqZJ0sgaTwh5eHh4ePQZnk/IwyNbpATbQrdtBBKSMUCApquHR7fc8NwneR3vooN3zbjNN77xDZ566imGDh3KsmXL2t9vaGjglFNOYdWqVYwdO5YHH3yQyspKAH77299y5513ous6N998M4ceemhe550tV155JSUlJfzoRz/q03kUAk8T8vBwQzIG0WZoWgsNK6Hxcy7e+UQu2ulEaFoDTV9Aw2fq0bIRYq1gJft61h7A2WefzTPPPPOl96+55hoOPPBAVqxYwYEHHsg111wDwIcffsgDDzzABx98wDPPPMP555+PZVm9Pe3tBk8IeXh0h21DtAUaV0PTamjdDNICowh8JUysncnE4bPAVwy+EvXQ/ZCMQmijI5A2QCKitCePPmG//fajqqrqS+8//vjjzJ8/H4D58+fz2GOPtb9/6qmn4vf7GTduHDvvvDNvvfXWl/a/5JJLmDhxIlOmTGnXUJ588knmzJnDHnvswUEHHcSmTZsApcnMnz+fQw45hLFjx/LII4/wk5/8hMmTJ3PYYYeRSCQAVaLspz/9KbNnz2b27Nl8+umnXzruypUrOeyww5gxYwb77rsvH32kUgUeeughJk2axNSpU9lvv/16/sX1Ep4Q8vDYFttWWk/j5xBSN5EOIeMDx9G+rP4DltV/sPW+QgPDrwSTv1RpUM1rlbYUD/fyiXikY9OmTQwfPhyA4cOHs3mzSmNZt24do0ePbt9u1KhRrFu3bqt9GxoaePTRR/nggw947733uOyyywDYZ599ePPNN3n33Xc59dRTufbaa9v3WblyJU8//TSPP/44Z5xxBgcccADvv/8+wWCQp59+un27srIy3nrrLS644AJ+8IMffGne5557Ln/84x9ZvHgxf/jDHzj//PMB+OUvf8l///tfli5dyhNPPJGfL6kX8HxCHh6diYeVxmMlwAyC2b2P58737wDghnk3dD+e4Qf8YMWVMPKVQEkN6GaeJ+6RL2QXWuu2EX5lZWUEAgG+9a1vceSRR3LUUUcBKi/qlFNOYcOGDcTjccaNG9e+z+GHH45pmkyePBnLsjjssMMAmDx5MqtWrWrf7rTTTmv/e9FFF2113NbWVl5//XVOOumk9vdisRgAe++9N2effTYnn3wyxx9/fA++gd7F04Q8PABsC0KblaAQGvhL8htkoPsczSgKjauUmc8z0fUpw4YNY8OGDQBs2LCBoUOHAkrzWbNmTft2a9euZcSIEVvtaxgGb731FieccAKPPfZYu0C58MILueCCC3j//ff5y1/+slWlCL/fD4CmaZim2S7YNE0jmezwH3YWeNsKP9u2qaioYMmSJe2P5cuXA3Drrbdy9dVXs2bNGqZNm0Z9fX3PvqBewhNCHh7JmPL5xFsck1sBtRQzCGaRMvO1blLCz6NPOProo7nnnnsAuOeeezjmmGPa33/ggQeIxWJ8/vnnrFixgtmzZ2+1b2trK83NzRxxxBHceOONLFmyBIDm5mZGjhzZPmYu/POf/2z/u+eeW9cnLCsrY9y4cTz00EOA0tqWLl0KKHPfnDlz+OUvf0l1dfVWgrQ/45njPLZvYiEVRKD7wAj0zjFTmla8DRJRKBsBhq93jt3HuAmpzjennXYaCxcuZMuWLYwaNYqrrrqKb37zm1xyySWcfPLJ3HnnnYwZM6b9xr777rtz8sknM3HiRAzD4M9//jO6vrVWHAqFOOaYY4hGo0gpueEGZZK98sorOemkkxg5ciRz587l888/z3q+sViMOXPmYNs2999//5c+v++++/jOd77D1VdfTSKR4NRTT2Xq1Kn8+Mc/ZsWKFUgpOfDAA5k6dSoffPBBF0foX4iu7J/9jZkzZ0qvqZ1HXpESIk3QVqc0kxxMbxctVPb6tD6hTCTjYMehbCT4inIfp5+yfPlyJkyY0NfTGDCkGnhWV1fnZbxly5aFJ02atDwvg/WApUuXVk+dOnVsV595mpDH9oeU0FoH0SZlfsuxrMx3p32353MxfGBr0LxGaUT+0p6P6eExgPCEkMf2hZTKFxMN9fiGv3PFzvmZk2aAWaxyikolBMryM67HgKNzlNz2gheY4LH9kBJAsVblk+khizYtZtGmxXmYGMoc6CtW/qloS37G9PAYAHiakMf2gZQq/ycayosAArhv+b0AzBw2Iy/jIbQOQSSEZ5rz2C7wNCGP7YNwPcRa8iaACobQVKBEy3oVPefhMcjxhJDH4CfcCOEG5XcZCGh6hyBKeG2xPQY3nhDyGNzEQioMuwdRcH2CZqhiqC3rVAkhj5yor69n2rRpTJs2jdraWkaOHNn+Oh6P5/VYTU1N3HLLLXkdc3vAE0Ieg5dEFEIblFYxkARQCt0EoauoOa+yQk4MGTKkvbzNt7/9bS666KL21z5f9wnCncvouKUvhFAu8+xveELIY3BiJZU5Sw8UrNHcRTMu4qIZF2XesCcYfrATSpsbAInlA4Hbb7+dWbNmMXXqVE444QTCYVXd/Oyzz+biiy/mgAMO4Kc//SkrV65k7ty5zJo1i8svv5ySkg5/4u9//3tmzZrFlClTuOKKKwDV2mHlypVMmzaNH//4x1867r333svs2bOZNm0a5513XnuPopKSEn7+858zdepU5s6d297+oa6ujhNOOIFZs2Yxa9YsXnvtNUBVZTj33HM55JBDOOuss6irq+Pggw9m+vTpnHfeeeywww5s2bKFX/ziF/ztb39rDz678MILR1599dVDC/S15kzBouOEEAHgZcDvHOdfUsorhBBVwD+BscAq4GQpZWOh5uGxHWLb0LoREAWtAzemdEzBxt4Ks0iZFXU/FFX2zjELxNef+fqX3jt07KGcutupRJIRzn/+/C99fszOx3DszsfSGG3k4oUXb/XZXw/7a9ZzOP744znnnHMAuOyyy7jzzju58MILAfjkk094/vnn0XWdo446iu9///ucdtpp3Hrrre37P/vss6xYsYK33noLKSVHH300L7/8Mtdccw3Lli1rryPXmeXLl/PPf/6T1157DdM0Of/887nvvvs466yzaGtrY+7cufz617/mJz/5CbfffjuXXXYZ3//+97nooovYZ599WL16NYceemh7sdLFixfz6quvEgwGueCCC/jKV77CpZdeyjPPPMNtt90GwDe/+U0OO+ww49prr8WyLB577LHKt99+u8+rJ2xLIUO0Y8BXpJStQggTeFUI8R/geOAFKeU1QohLgEuAnxZwHh7bG5EGZYrzFTYQ4fX1rwOw14i9CnocQAVVtNU5vYoGX3mf3mTZsmVcdtllNDU10draulXr7pNOOqm9Ttwbb7zR3uju9NNPb29c9+yzz/Lss8+yxx57AKqY6YoVKxgzpvtFyQsvvMDixYuZNWsWAJFIpL1qt8/na28FMWPGDJ577jkAnn/+eT788MP2MVpaWgiFQoAqshoMBgF49dVXefTRRwE47LDD2luUjx07loqKCl577bXghg0bzN133z1cW1vb7+y6BRNCUhWla3Vems5DAscA85z37wEW4gkhj3wRC6lIOF/hQ7Ef+kQVvOwVISSEqsDdsh4qdxiw/YjSaS5BI5j288pAZU6az7acffbZPPbYY0ydOpW7776bhQsXtn9WXJx54SKl5NJLL+W8887b6v101Q6klMyfP5/f/va3X/qsc1sHXdfb/Ty2bfPGG2+0C5vOdJ5nuvqfxx9/fOKOO+6o3rx5s/n1r3+9X/Z2KKhPSAihCyGWAJuB56SU/wOGSSk3ADh/u7RRCiHOFUIsEkIsqqurK+Q0PQYLybhK9ByogQiZ0AwVqNC6yfMP9YBQKMTw4cNJJBLcd9993W43d+5cHn74YQAeeOCB9vcPPfRQ7rrrLlpb1Rp73bp1bN68mdLS0nZNZVsOPPBA/vWvf7V3b21oaOCLL75IO89DDjmEP/3pT+2vuzLzgerm+uCDDwJKS2ts7PBuHHTQQdaCBQvKly5dWnzCCSc0pz1gH1FQISSltKSU04BRwGwhxKQs9r1NSjlTSjmzpqamYHP0GCTYthJAmlmwQIR+gRmAeETlPnnkxK9+9SvmzJnDwQcfzG677dbtdjfeeCPXX389s2fPZsOGDZSXlwNKOJx++unsueeeTJ48mRNPPJFQKMSQIUPYe++9mTRp0pcCEyZOnMjVV1/NIYccwpQpUzj44IPbG+p1x80338yiRYuYMmUKEydO3Mov1ZkrrriCZ599lunTp/Of//yH4cOHU1qqqm34fD722muvlqOPPrrBMPpngZxea+UghLgCaAPOAeZJKTcIIYYDC6WU49Pt67Vy8MhIWz1EGgvuB+pMXlo55IKUEG+FitHKRNePGcitHMLhMMFgECEEDzzwAPfffz+PP/54X0/rS8RiMXRdxzAM3njjDb7zne+0a03vvfde+OSTTxYPPfTQysmTJ8f6ao590spBCFEDJKSUTUKIIHAQ8DvgCWA+cI3zt//9Vz0GFomIKsvTC36gfoEQqgFfaCNUjBncml8fsnjxYi644AKklFRUVHDXXXf19ZS6ZPXq1Zx88snYto3P5+P2228H4MMPP+Twww8PHnHEEZv7UgBlomCakBBiCirwQEeZ/R6UUv5SCDEEeBAYA6wGTpJSNqQby9OEPLrFtlRrbqH3urN+c0TZ94cG+yj1ItEGvjIo7XepH+0MZE1oMLBdN7WTUr4H7NHF+/XAgYU6rsd2Rls9SKv3WnN3os+ETwqzWDXm8xf3qhnSwyOfeBUTPAYu8bC6CRt9kzezYM0CFqxZ0CfHbscMKrOcNfDLt3hsn/TPcAkPj0zYlhOOHezdcGxpg50E2+K/nzwCwAHDZgFC+WY0XYVS00tz0gyw4sonVjqsd47p4ZFHPCHkMTAJNwC2c8MvMFZCaV2JNkjGSAmYkkRUPWuto0PoyI7EUl+JqnBQ6DmaRRBtVr2SPLOcxwDDE0IeA49ExAnHLmQ0nFRN5SLNStMQQuUgdSqZk0gJl23L6EggmYB4HWCrDqn+MiWQCoUZhNAmVU3Bi5bzGEB4QshjYGHb6mZr+AtkhpMQa1VCzrZUxF22tdoEYDiVqiQquTQWUoIiWFUYYaQZYMXUvIur8z++h0eB8ISQx8Ai2qx8MoUwOyWjKtrOiitBkQ9hIQDTD/iVdtSyXkWzBavyb6Yzijrq5pm9Hy3oiqY16nvOF0ZAJe3miTVr1nDWWWexceNGNE3j3HPP5fvf/37W40SjUfbbbz9isRjJZJITTzyRq666Km/zHEx4Qshj4JCMQ7gOzDyb4aStouwiTaD7XGs+F8/8YXbHSWlH8QjE1kLxEGWqyxdCKMHZulndmPtj/bxkNL9m1Hhr5m2ywDAMrrvuOqZPn04oFGLGjBkcfPDBTJw4Matx/H4/L774IiUlJSQSCfbZZx8OP/xw5s6dm9f5Dga8EG2PgYGU0LYZNF9+b65WHJrXQaRF5d1kkfBa7iul3JeDEDEDSli01alzknmsrq/7VPBEtF/Wquwzli5dyn777cfEiRPRNA0hRHszus4MHz6c6dOnA1BaWsqECRNYt25dVmMACCHam+AlEgkSiUR7pexsxtke8DQhj4FBvFVpEP48rqJjISUIdB/4sq/BtmDNQgAOGD0v+2NrmjIpxiOQWAeltWoe+cBXpM7Ll51QHaxEo1FOOeUU/va3vzF79mx+8YtfEI1GufLKK9Put2rVKt59913mzJmT0xiWZTFjxgw+/fRTvvvd7+Y8zmDH04Q8+j+2pUxMeSvWKVXju7Y6Fd6c44164ZqFLFyzsGdTMQOABs1rVTRePhCaKmMU7pftY3qd559/nunTpzN79mwApkyZQkNDQ7tm0hWtra2ccMIJ3HjjjZSVleU0hq7rLFmyhLVr1/LWW2+xbNmynMYZ7HiakEf/J9IIyPyEHktbCZ94mzK/9Ydr3zDBdvoEBYdAsLznY5pBiLao0PDtvBPrsmXLmDx5cvvrd955p93k1hWJRIITTjiBr33taxx//PE5jdGZiooK5s2bxzPPPEMymcx5nMGKpwl59G+SMbWiz0dpHmmpG308okxV/UEApdA0pZVF6lWABHkoLJwKUrDtno81gBkyZAjvvfceAJ988gmPPPIIp556KqCazaV8PqC6lH7zm99kwoQJXHzxxa7G6Gqcuro6mpqaANXK+/nnn2e33XbLOM72iKcJefRfpHR8NnnICZIWtGx0wrv7aQ8eIZR2FmkEbAhW0iNJqftUzlOsBYIVeZpkDzEC+Y1oc1G49rTTTuOJJ55g0qRJVFdXc//99zNkyBBs2+bTTz+lqqqqfdvXXnuNv//970yePJlp06YB8Jvf/KbbMYAux9mwYQPz58/Hsixs2+bkk0/mqKOOorW1tdtxtlc8IeTRf0mEVbmcnoYxpwSQTPbf/JkUAkcjalKveyqIfEXQtkWFRev94HLPY06PW0pKSnjyySe/9P6HH37ICSecQDDYsSjZZ5996K69TVdjdDfOlClTePfdd13PZXumH/wqPTy6wLbzE4zQWQDlud3DJXMuzet47QgcjahJvQhW9mAsTWlY4YZ+3XeoL5g0aRLXX399vxlne8XzCXn0T6LNKiquJ1UFpCPI7PwLIICg7ieoF6geXLsgalT163qCEVTJuMl+21zTYzvGE0Ie/Q8rAeEtPdSCpBojESuYCe6ZVc/yzKpnCzI20EkQbVE5TTmPI5R/qK1O+dk8PPoRnhDy6H+EGx0zUg9+npFG5ZQvYBDCG+tf5431rxdsfKBDELXVqerhuWL4nXYU4bxNzcMjH3hCyKN/kYg63VJ7IDxiLcqfYg6S3jpCKHNiaKMqM5QrRsAL2fbod3hCyKP/IB0Tmt6D+nCJiIoGM4v6Vx5QT9F0Vdkh5ISZ54JuqjbgsZb8zs3Dowd4Qsij/5AKyc61hYKVcHoNBfpnBemeopsduVO5JrOmQratHAWZh0ee8YSQR/9AStUmO9coNunUl9P0wd1Z1AyoYItIY277C01piNtJle36+nqmTZvGtGnTqK2tZeTIke2v4/EemDa7oKmpiVtuuSWvY+aTkSNHTt6wYUPacNObb755yFlnnTUG4Nprr63505/+1G0m7VNPPVX63HPP9djmXbA8ISHEaOBvQC1gA7dJKW8SQlwJnAPUOZv+TEr570LNw2OAEAv1oFmdVKV97EQei5xm5qq9ruy1Y22FGezU+yiHquKGUx7IXwpGnip391OGDBnCkiVLALjyyispKSnhRz/6Ucb9kskkhpHd7TElhM4///xcppr3+fSUn/zkJ3XpPn/xxRdLS0pKrIMPPrhHlXcLqQklgR9KKScAc4HvCiFSnaFukFJOcx6eANresS1lYspVC4qGINbWqwKoTxGowI3WutwCFYQAYahK4tsht99+O7NmzWLq1KmccMIJhMMqYvDss8/m4osv5oADDuCnP/0pK1euZO7cucyaNYvLL7+8vT8QwO9//3tmzZrFlClT2nsBXXLJJaxcuZJp06bx4x//eKtjrlq1igkTJnDOOeew++67c8ghhxCJqGjHJUuWMHfuXKZMmcJxxx1HY6PScufNm8fPfvYz9t9/f2666SbmzZvHRRddxH777ceECRN4++23Of7449lll1247LLL2o917LHHMmPGDHbffXf++c9/ZpRcN91005CxY8dOmjVr1vjXX3+9/SQvvvjiEZdffvkwgKuvvnroTjvttPuuu+468aijjtrx448/9v3tb3+rufXWW4fttttuE5955pmce6wUTLRKKTcAG5znISHEcmBkoY7nMYCJtihzXC5mtGTM6bba+5Fwj69U5VeO2emrvX5sNA00E0KboXxE9uHsqSrbgYreLWX01yO//N7ux8Lsc5Q/8L6Tvvz5tNNhj6+p1usPnrX1Z19/OuspHH/88ZxzzjkAXHbZZdx5551ceOGFgCoq+vzzz6PrOkcddRTf//73Oe2007j11lvb93/22WdZsWIFb731FlJKjj76aF5++WWuueYali1b1q55bcuKFSu4//77uf322zn55JN5+OGHOeOMMzjrrLP44x//yP7778/ll1/OVVddxY033ggo7eqll14CVNkgn8/Hyy+/zE033cQxxxzD4sWLqaqqYqedduKiiy5iyJAh3HXXXVRVVRGJRJg0aZL53e9+V6+tre2yc+IXX3xhXnPNNSMWL168vKqqytprr73GT5o06Utx/DfffHPtF1988X4wGJRbtmzRq6urrbPOOquupKTE+uUvf7kp639CJ3rFJySEGAvsAfzPeesCIcR7Qoi7hBBd1iQRQpwrhFgkhFhUV5dWK/QYyFjJ3BNTU1Wx81HgNAfe2bSYdzYt7vXjtmOYqhxROEeNRvep7347Y9myZey7775MnjyZ++67jw8++KD9s5NOOgldV4uhN954g5NOUkLx9NNPb9/m2Wef5dlnn2WPPfZg+vTpfPTRR6xYsSLjcceNG9deFHXGjBmsWrWK5uZmmpqa2H///QGYP38+L7/8cvs+p5xyylZjHH300QBMnjyZ3XffneHDh+P3+9lxxx1Zs2YNADfffDNTp05l7ty5bNq0SXzwwQfdrjJefvnl4rlz54ZGjBiRDAQC8vjjj+/yxzR+/PjIcccdN+6WW26pMk0zrxnPBTcyCiFKgIeBH0gpW4QQ/wf8ChXe8yvgOuAb2+4npbwNuA1g5syZXpr3YCXSCIjcElPD9UqDMrbj7qFmUIVcm4Hs/UOG3+lYG+69nkPpNBdfUfrPi4fkpPlsy9lnn81jjz3G1KlTufvuu1m4cGHHIYoza9RSSi699FLOO++8rd5ftWpV2v38/o6oT13X281x6dh2PqkxNE3bajxN00gmkyxcuJDnn3+eN954g6KiImbOnGlHIpG0F5ebhnoLFixY8Z///Kf0scceq7j22mtHrFixYlnGnVxSUE1ICGGiBNB9UspHAKSUm6SUlpTSBm4HZhdyDh79mGQcoo0qpydbYiFVEaG/V8XuDdr9Q4ns99V9KmR7OyrnEwqFGD58OIlEgvvuu6/b7ebOncvDDz8MwAMPPND+/qGHHspdd91Fa6tqSbFu3To2b95MaWkpoVB25ZXKy8uprKzklVdeAeDvf/97u1aUC83NzVRWVlJUVMRHH33E+++/n/Yev99++7W9+eabpRs3btRjsZh49NFHv2SZsiyLlStX+r761a+GbrnllrWhUEhvbm7WS0tLrVAo1ONQ1EJGxwngTmC5lPL6Tu8Pd/xFAMcBeZOoHgOMSINykGdrSrMSTiDDdhKIkAlNU4VeWx3/UDZZurpPCfR4a89bZgwQfvWrXzFnzhx22GEHJk+e3K3guPHGGznjjDO47rrrOPLIIykvVx1vDznkEJYvX86ee+4JqPYM9957LzvttBN77703kyZN4vDDD+f3v/+9q/ncc889fPvb3yYcDrPjjjvy17/+NedzO+yww7j11luZMmUK48ePZ/LkyWnLY+ywww6Jn/70p+vnzp07oaamJjFlypSwZVlb/YCSyaQ4/fTTx4VCIV1KKc4777xN1dXV1gknnNB04okn7vSf//yn4sYbb1xdX1+vv/3228U33njj+mzmLLrrndFThBD7AK8A76NCtAF+BpwGTEOZ41YB53USSl0yc+ZMuWjRooLM06OPSEShaXUONz4JzetVhezeDi+WthKAtgV2ghsWXYdA8oM9fqAEgdBVsECqukFv+6niYQhWZd8e3E6qR8UO6jzyyPLly5kwYUJex+wtwuEwwWAQIQQPPPAA999/P48//nhfTysrli1bFp40adLyvp7H0qVLq6dOnTq2q88KGR33Kl0vybyQbI+O8jzZEmlRYcm95cOwEkpgxttUHlIKoXHRlG8DQs3Hksqk1b6okyrk3Fek/vZGAq1ZpCpum4Hsqk5ohnOOIQhkKcAGMYsXL+aCCy5ASklFRQV33XVXX09pUOI1tfPofRKR3DqmJmPOTTbnlAR3SOkUUm12BI9Qmk02EXwStW/YadXtKwZfgZNDBSpSsLUu+7BtM9jRgXUwV5zIgn333ZelS5f29TQGPZ4Q8uhd2svzZFkfTtrqJqn7C1uYNN62dUO9NILnyc9UpNZXd+wi90WgBJduKoGUiHXUxQuU514fLxO6CfGIOodsOrJqOiSlyh0q6kEnVw+PLPFqx3n0Lomw0miyNcVFm5VpTC9QOHYyptqAhxuUb8cMZjzWxw0f8XHDR5nHFigNyAwq4da6ySkzVKAiomZQhb4no9nvF673ipt69CqeEPLoPXLVgpJOwc5ClOWRtjKZtW4CpDpGIc1Ruql8N4kotGxQUWn5RqCEfGudOj/X+21fxU09+geeEPLoPeKtyk+SjTaTMsNpvvyb4VKCIBFWgqFQWlZXGH4VsBBudBrN5Vn70E2ldWUrUIwiRxvKIefIwyMHPCHk0TvYjjDJtkhpNKRuiPmuihANQdtmZXorlH8mE0IozctKKFNgIkvzWSbazXKx7Oak6apK9yBB1/X29g3Tpk1j1apV7LXXXoCqcvCPf/yjfdslS5bw739nH8A7b948+mMayaWXXur761//mtbJ9/HHH/t22WWX3QFefvnlorPPPnt0um1vvfXWqnzO0RNCHr1DPKR8DVoWsTBWXLUcyKcZTtrQusVpIV4Eeu6xOcVmCcX5iNQz/EpzadushGO+ECgNsm1LdmY5I6gqWSTz22+nrwgGgyxZsqT9MXbsWF5//XUgf0KoECSTve+b22+//cJ33333mu4+X7Fihf+f//xnXoWQFx3nUXhsW1VAzkqYSMcMZ+bPDGclVaUFaedFsH1n6nmZN3KLpoMoUjd/O6Ei2/KR7GqYKiovGnKfxNre6qERSof1fA79kJKSElpbW7nkkktYvnw506ZN47TTTuPPf/4zkUiEV199lUsvvZSjjjqKCy+8kPfff59kMsmVV17JMcccQyQS4etf/zoffvghEyZM6LYO3NixY5k/fz5PPvkkiUSChx56iN12242Ghga+8Y1v8Nlnn1FUVMRtt93GlClTuPLKK1m/fj2rVq2iurqaXXfdlc8//5wNGzbwySefcP311/Pmm2/yn//8h5EjR/Lkk09imia//OUvefLJJ4lEIuy111785S9/SVsT7pVXXin61re+NTYYDNpz5sxpd0w+9dRTpdddd92wBQsWfPr000+X/PCHPxwDqr7c66+//tHPf/7zkZ999llgt912m3jaaadtueKKKzb39H/hCSGPwhNtVjf+bBz+sZAyI+UrKTUZVwJIiL4zv2VCoHxT8bDyERVX51bYdVvMoCqR5MvC72UGnTDvirx+X1+cedaX3is9/DCqTj8dOxJhzblfFuzlxx1HxfHHkWxsZN33vr/VZzv8/W8ZjxmJRNqrV48bN45HH320/bNrrrmGP/zhDzz11FMADBs2jEWLFvGnP/0JgJ/97Gd85Stf4a677qKpqYnZs2dz0EEH8Ze//IWioiLee+893nvvPaZPn97t8aurq3nnnXe45ZZb+MMf/sAdd9zBFVdcwR577MFjjz3Giy++yFlnndXeAmLx4sW8+uqrBINBrrzySlauXMmCBQv48MMP2XPPPXn44Ye59tprOe6443j66ac59thjueCCC7j88ssBOPPMM3nqqaf46le7bzHyzW9+c+wNN9yw+sgjj2w977zzRnW1zXXXXVd78803f3HIIYe0NTc3a0VFRfavf/3rdSkhlfGLd4lnjvMoLLalboDZaB52IgfNKQ2JKLRuVP6fXKo0dMMjKx7lkRWPZt4wW8yA+g5Cm9X311OEUOce3oJKWnKJZubeJqIf0dkc11kAueHZZ5/lmmuuYdq0acybN49oNMrq1at5+eWXOeOMMwCYMmUKU6ZM6XaM448/Huho3wDw6quvcuaZZwLwla98hfr6epqbVRDJ0UcfTTDY8ds//PDDMU2TyZMnY1kWhx12GKDaOaTGW7BgAXPmzGHy5Mm8+OKLW7Wn2Jb6+no9FArpRx55ZCvAN77xjfqutps7d27rj370o9FXX3310C1btuimWZjAHU8T8igsqYZ12azoI41OuHAezFHxCIQ3g57/0jmfNX+W1/G2Qvcr7a11E5QMzc6X1hWmX2lYsSwKlZoBpZEmKvNWrTyd5qIFg2k/NyorXWk++URKycMPP8z48eO/9JmbFgjQ0X5B1/V2P09XNTtT46Vr32CaZvt2qfYN0WiU888/n0WLFjF69GiuvPJKotHug1yklK7m/pvf/Gbjscce2/z444+X77XXXhOeeeaZT1ycbtZ4mpBH4bCSKtw3G40mEc5fi4Z2AVTg3J9CkSrxk68QbiOgNJtsxtJ96n84SNm2/cK2rw899FD++Mc/tguNd999F4D99tuvvQ3EsmXLeO+997I6buf9Fy5cSHV1NWVlZTmdQ0rgVFdX09rayr/+9a+021dXV1slJSXWf//73xKAu+++u8tAgw8++MA/e/bsyK9//euNkydPblu2bFmgvLzcam1tzevF5Akhj8IRbVZ+DrdaUOfSPD1lKwE0gH/meh4FUep7yCb82vCrUkaJzA3YBiJTpkzBMAymTp3KDTfcwAEHHMCHH37ItGnT+Oc//8kvfvELEokEU6ZMYdKkSfziF78A4Dvf+Q6tra1MmTKFa6+9ltmzs2uLduWVV7Jo0SKmTJnCJZdcwj333JPzOVRUVHDOOecwefJkjj32WGbNmpVxnzvvvHPV9773vTHTpk3bLRgMdmmjvfbaa4fusssuu48fP35iMBi0TzzxxObZs2dHDMOQ48ePn3jVVVcNXbVqlbn//vvvnPPkKWArh3zitXIYgFhJaPwczGL3ZrVII0SawddDX1AiqsKd9UBBBdAfFl0HwI9m/rBgx2gnGVffY8nQnml1Eki0QdkI9zlbybj6Hiu6TR/ploHcymEwMBBaOQzgJaJHvybalJ1fx0rkpzRPMu4IIH/BNaCKQCUVgV4q9mn4nBynLEvxbEuqpE9bPa6DFAxfR+VzD4884wUmeOSfdoFSnHlbAKSK3OppTpCVVGYrzdcrPqBvTfpGwY+xFYZfha23bYHimtwDN/ROuUMBl34Iw6+Oa47u/WZ9HoMaTxPyyD+Rpuy0oJTPoSe9dmxLaUCa3qMqCP0eww/JhAp77wmpitlu/Uy6T1XlTmSvDQ0Ek79H4bBtW9DRXftLeELII78k4yrr33BpVpOWMg1lW1NuqzGc6txS9moR0gc+fpAHPn6w147XTircuifVroVQC4VshJnhV/+rLIRKIBCgvr7eE0TbKbZti7q6unJgWXfbDOIlo0efEGlUJV/cakHRFkD2zHwWaQCZ7Jkgy4G1oW5LbBUew6looBmqa2sumAEnb6jM3Xen+9T2ibDrY44aNYq1a9dSV1eX2xw9esTGjRsNy7Kq+3AKNrAsmUx+q7sNPCHkkT+SMXVj9Lks6mnFs/QddUE0pLSCTgEND7+zzvXuJ0wfmfux+xKBEkTheiWIci2tkwpSKB+BK4ec4bQPryxytdAwTZNx48blNjePHjNx4sT3pZQz+3oe6fDMcR75I9LoBBe40YKkcwPtQTBCPJyd6W+wIURHlexcc4h0Uy0GYi6rd+umKilUiGZ8HtslnhDyyA/JmDKtua10EA/3LBjBSqiIOj2Q/2Z3AwndAIRqT5Gr38UMQlsWlRSMgNMewvPzePScgpnjhBCjgb8BtSi74G1SypuEEFXAP4GxwCrgZCllY6Hm4dFLtDlajRtkjg3uttq/Th2vF6ohpGptWW0RNt32LxL1TditEb4S2oLUBQ3rF1J19DykbRN67V2Cu+2IWdNL+UPg5PFEVVRiUQ7HFUI9Ik2qcncmNENFysWyCPH28OgGV0JICHEtcDUQAZ4BpgI/kFLem2a3JPBDKeU7QohSYLEQ4jngbOAFKeU1QohLgEuAn/bgHDz6mkRUmWfcFsaMhtQqOtdghLYGtX9PQrozkNjcQPPLiwm9+R7BXcZQe95JaAE/kY9XoZcWY1SWMXRoFdK20cuVDyyxsZ51194NgDm8htI5kynbfybBnbOvNJA1ZkA1DjR8uQUqGAGItaj/oRv/Ukob8pUM7LJIHn2OW03oECnlT4QQxwFrgZOABUC3QkhKuQHY4DwPCSGWAyOBY4B5zmb3AAvxhNDAJrzFfYsEO+F0S82xT1A0BMlIfrutdiL8wUoaHl9A6H/vgS0J7DwG/+haAISusdOtv+h2X3NoFWNv+AmRD1fSuvhDGp56iYbHXmTUL86ldPbkgsx3K9oDFczsBXTnLqxughQ0w1l8hCDgslmeh0cXuBVCKTvLEcD9UsoGt2XMAYQQY4E9gP8BwxwBhZRygxBiaDf7nAucCzBmzBjXx/LoZRIRVSzU7zIiLtygetvkknWfjDmBCHlqdNcFTS/8j/CyFQw5/iAqDt0LX61jnoqFoGUt+CugbDiEG1j3/M8AycjS0WCWIIqrCY6eS/DoeVQdPQ+rNUzzwrcpnrYbAM0LF2HHYlQcNBehF6CigxBOD6AtUFqbfUM8w1SJw/E2dxGOZtDRhko9bcgjZ9wKoSeFEB+hzHHnCyFqgO4bVnRCCFECPIwy37W4FV5SytuA20AVMHU5T4/epi0LLShVfyyXbqm25bT79uU1ECG2ZiOb//o41acdTnCXMQyd/1VqzzsRze+DLZ/C0udh7dvQtErtMOEYmH6WCgiIhwGhiq42rYUvtkDJMKjeGVo3or91O1W7zIR4E/hqCL2xlNDrS2h84iWGfv0YimdMdN2TxjW6qYR1pBGKhmS/vxHsaCgoMghKTYekVGH5ufiiPDxwKYSklJcIIX4HtEgpLSFEGGVWS4sQwkQJoPuklI84b28SQgx3tKDhQI97lHv0EfGwclC7yguS6uaWa0WDsBO7kq+SPEmLmpde47NX30QL+Ehs3EJwlzEYZSVKo5A2LLhamZuqx8OUU6FqJ6jaUe3vL+P+sVOBTlW0rSTt1UkijRDaAIvuUI9hkxj5tSMI7T+dzXc/yZqrbqV0z6nUfvcUjHKXvjS3pNovGH73OVspNA0sqSIdgy4ES6p1eKBsYPZs8uhz3AYmFAHfBcagTGQjgPHAU2n2EcCdwHIp5fWdPnoCmA9c4/x9PKeZe/QtUmanBcVCKgQ4lzYNefYD+TfXMeKRpwhu2ET5gXMY+vVjMHwWvH07bFoOR16nTFnzLlEtD/wuI8A6C8iaCXDMLdCyHr54DT55DvHK7yk7/nZKZ/2M+scXsOW+fxP55AtKZ03Ky3lthRFUpk/dn73gN4JKiPpKMu8rNPVbiLZ42pBHTrhdVv4VWAzs5bxeCzxEGiEE7A2cCbwvhFjivPczlPB5UAjxTWA1KsjBY6CRCCuzjxtfkJ1UEW25ZPUn4xBtyKsfqPTDjzFbQqw+7XgOPWlP+PBx9bATsON+kAyrKg41u/X8YGUjYPJJsPvxUP8ZBCsRQPWozyi/6gjMKUoAtS35iODEndB8eap9l/IPtdVl7x8SqNJL4QYoHZZ5+1QxVH/p4C4e61EQ3P5idpJSniKEOA1AShkRGYzZUspX6d56f2AWc/Tob0ipkiPdCpVos3NTzNJ5LW2nxYOvx+0DRCKJ2dhEfGg1W/bdk8ZZ0wmaLfD4+RBrhtFzYfqZUFLresxRpVmEXms61OyinsdbYctyzOZnYf2LJHb9GquvuoPAuJGM+tm3MKvzpFHopjKXRnLw2ZidOqpm0kCFpq70aDMU5+CH8tiucSuE4kKIIE4XLCHETkCsYLPy6N/EW8GKucsLsuKqwZ2ZpW8C1M1T2j0uTGq0hBj9wCMYzS18+r3zkH4Dq7iINvwwag7sNA9qxmc97qnjTwZyqFXnK4Ejr4eVL8K792K+8UtGnfoV1j/8KZ9f9HtGXfJNinbfKev5dImeyh/yZx8QojtVsytGkjEaxChyfEPlnjbkkRVul6ZXoJJURwsh7gNeAH5SsFl59F9sp9qBW/9MuD63iLZ42Ll59kwA+TdsYtxf7sZXV8/Gow5leOxDDlz5W3xWCIkOc8/LSQD1GKHBzgfBV2+C0TMp1d9i7G+/g14c5Iuf30zzwjy1sxcoQRTOoixPCt1QJko3deVSrSGiTbnM0mM7xm103HNCiHeAuaif9fellFsKOjOP/km8VYVLuxEOCac+XLYZ/HbScar3TAAVrVrN6H88jO338cU3T2MXXmfXdS/R4B+FaUWJ6z2LSrtj2V0AVHJo7oMEymG/n0CkEX+wkrHX/oC1V/2R+IY8Bo1qGtgCwk1QkmVVfyOg/HlmkUpQTbttsJM21Ht9nTwGNml/VUKI6du8tcH5O0YIMUZK+U5hpuXRL0lpQW4EUKo+nJ5DMEK4KTcf0jZU/W8xydJiNp5+JLNaHmRodCUrKvdj2bBjsEXPb5JNURU2XpmPPE0nHFrf8Apjpi6GXUeCbZHY0oxRU9nzfCLDrxYF8dbswrY1p0NutDlz3pEQKrco3AilXeage3h8iUya0HVpPpPAV/I4F4/+TjykhIubfJBoSAktX5ZCKN7qRKf1IBrOtkHTWHfcUWiJBHs0/YvK6Gr+N+Is1pb369YqsPPBiKY18NETJNZ8yuf/jFM2bxbDzj0hD4IoqASEHsjOb2MElJnNX5o5JN8IqqoWwYqC1vbzGDyk/SVKKQ/orYl49HNSFQvcaEF2Upllss3rsZLqJtmD/kAtr77LuHv/wxdnnIws8mH5ing3eBIf1xxEyDci53F7Dd2EOefBkJ0x3vw/ynccQcNTLyE0wdBvHd8zQSSEWkC0bYEy91GA7XXlwvVQOtzFMXzq/1+axTE8tlvcJqsGgPOBfVAa0CvArVJKV6V7PAYB0Rb3la8jTU7YbpY3zLYtavwcb7Qtbyxl3e/vRo4awZjQIsZsWsRrO3wXSwsSyiVJti/Z+UBESTVD/Tchh+5FwxMLQdcZ+vVjeiaIdJ/y00VDEMjCJ2b4nJBtF1qqGVC/l0CF+/5SHtstbnXyvwEh4I/O69OAv+Mlmm4fWEmVr+PGRJaMqZYA2bbsjoZUJFaOVRFCby1j3bV/JbjLGKz9y5jU8BAbisbTXkanAOxYrkr4SJdNSbOmdiri+FsZppngf4iGR1/AP6aWioPm9mxcI+gkAPuzM5npfqfK9qjMya+6ozmVD9D26R69hlshNF5KObXT6wVCiKWFmJBHPyTa7Gg2mTzwObbsthI9atMd/vAz1l1zJ4GxIxh94lD0lY/xRdlM3hlxOrYoXM7K8bscB2SXJ5Q1ug8BDDt4KP7NrZTtlgfNot285lREcKtZ6QbEE2rBEMzQvsHwQ6zVXbKrx3aN27ied4UQ7csvIcQc4LXCTMmjX2EllGBxIyDibSpDP5vVtZQdXVlzNDMZVWUU77Ebo08bg77yMT4rn8PikV8rqADqbcSIaVTuUYr28u+wPv0fre8s79mAuql8d7GW7PYzA8rf4ybnSPdBa53XBtwjLW6v0jnAWUKI1c7rMcByIcT7gJRSTinI7Dz6nnCjOz+NtNXKOtvk0lhrzmY4K9SGVhzEV1vN6F+cB6GN4LN4194P9+ur3Pm/pX8BYChHFfxYBCvgkF/Bc1ex+U+307SqmNFXfIeSPXpQ384IKC3XCLpfOKSSUt20Ajd8KtE1Ec6t26vHdoHbK/UwYBywv/MYh2pwdxTw1cJMzaPPaW8i50JApEKysynnbyUc30QOAqg1zBeX3MTGWx+Cde+qN0trYdrX6A0BBNCWaKUt0dorxwLAXw4HX8XQ/crwl8ZZ++vbiK5ck/t4AqcJXn122kqqFXjSReUuI6DqDHrakEc3uLpapZRfAC1AOTAk9ZBSfuF85jEYCTe4Kx5qJ5yQ7Cy0oHYzXPYlfaRlse73dxNbv5myUW2w8Gr47KXsBhmo+EvRj7yC0Rcfg15azJqrbyPR0Jz7eLqpwu+zMct1Fl5kEC66CXbcXekfj+0StyHavwLOBlbS8avzklUHM4mounG4KVKaIST73jdXf+k90wrhT4awtK21oBOmZ46m2nzXY7S9s5zak6dQHHoaxuwF4/bNPM/Bgr8cc4/DGP2LSaz6yXVsvPluRl/5/dzHy8UslwrZdtMK3Ag4bcBLvDbgHl/CrU/oZFQ7h3ghJ+PRj2irc1dyJxlTfp0sKjTrMqEEkMg+0qvp2TdoeGIhlQfsSqX2DIyYAXt/P7t+OYOEwNjhjD4wga/4PWjdDCU5lspp12yyjJYzAk4r8KL0379mdITuBytym6PHoMXtlbsMqCjgPDz6E6m23RlXxU53VS2bOmwSf6IJiZFTNJw5rIqyvacwrPZN1XZ7vx9lLqxZIMZX7cb4qjw0vssVTaf49J9g+mPI539F69vv5j5We7RcFj4uTVf2kKgLU5sZVL8V28p5ih6DE7dX729RYdrL6NRHSEp5dEFm5dF3SKnCat1oQTGnr1AWkU+G1YYmE18yw2XCTiTQTJPiqeMpnjoeNu+tkibdthcvAF/d8UigwHlCmSgfA/MuoeFPf2DznXcx6uffpHTutNzGMgIqUMQMuK+CbfohUg/+ovSLkVTju0iT1/jOYyvcCqF7gN8B71PIFHSPvicWUo7kTHZ+aTkh2e6FiSYTBJItWCK7oqZ2JMaqn95A+X5TGTJnCIyeDUMnZjXGoGbYJCrPnE/LFw+w/rp7GHfTCHwjcjDNpZJY2+rd15ZLVc6ONEJxhmMaRSqYIVDmtXrwaMetOW6LlPJmKeUCKeVLqUdBZ+bR+2TTqiHaAsgsHM2SQLIFMLLy30gp2fCn+4l9sR5/4wJ45Qa1mu4H3PjuTdz47k19PQ0AtF0PYNRPzkSYftb+5g7saI6Nj3VTRTu6MbGlMANKK05mKCWZKqAabsxtbh6DErd3g8VCiN8KIfYUQkxPPQo6M4/eJ+q0087kY7ESauWbhRZkWGE0O4aVlf8IGv/9Ci0vL6Zm7zJKilbCPhf2G+d20kqStLLsVlpAzEn7MOLHZxNbvZGN1/1f7gMZAZUfls256Y4GlSlk2wiqthBucow8tgvcmuP2cP52rpzohWgPJqykMpW4qVwQbgBhuM7v0WSSQDKUtRku8vEqNt3xCCXjKxky4gOYdpYKx/bolpI9dqP24CoCvA11K6Bml+wHEUItRML1KlrODbqpAlpirenD+oVwBNYWr7ipB+A+WfWALh6eABpMRJs6SrKkIxGBRJtySLvEn2wBXIy9DbE1G/FVlzFi9w8ROx8Aux+T1f7bK5XnXkxwVAm89FuSG76co+UK3QdWXDUZdIsRUAsUmSECzvA7OUbh3ObmMahwfVcQQhwphPiJEOLy1CPD9ncJITY7EXWp964UQqwTQixxHkf0ZPIeeSIZVxUPMpnXpO0+fyhFvA3djmJp2UexVRw0l3G3/AL9wB/CrPOy3r9gSMC2ETgPO4mQlvp++gP+MjjgZ9S/B5//4Pckt2zJbRwjoPw3boqVguMflMqsm3Fspy2EV85nu8eVEBJC3AqcAlyIMsKcBOyQYbe7UTXntuUGKeU05/HvLObqUSjatoBwUcU6VR/ObWto24JwI7bITgDVP/oCoVf+B02r0EwTxuzZd9FUEuUDS0YdLTACVhSkxeTqKUyqngKaACS6TKLb0U6PGMLtDTzfVOxA8RGnYsUk66++AWnnICBTkW/hJvf7GEHlL7Qy5LXrPhXe75Xz2e5x6xPaS0o5RQjxnpTyKiHEdcAj6XaQUr4shBjb4xl6FJZERJlcMpXnaa8Pl0V+j2Pik8J9UdO2pR+z+a+PU76bSenGOjj2VnfN9PKJRN0gpVRLLiMAZpnKg9E7ovsO2f1rwLZliSSatBDSRsg4PkcggVQJur1IYO4hDDtlIxv/8Tb1Dz9P9UmHZD+I4XRijbe5ywfbqhV4LWkdh4aTwOorzq7wrcegwu1VEXH+hoUQI4AGVCXtXLhACHEWsAj4oZSyy3hNIcS5wLkAY8aMyfFQHmlJJaYaLsxrkcbsWnbHnRtXFgJEbwuz/o6/4xvio3b3VbDH+b0rgGzLWcEL8BerY+suCrhuhVB9jASAj6RegsBCtxP4LKc9ttBzGDc3Kk49k/Bai7p7n6ZoTBFFc/bJfhDdr8xyRsCdsHBbV07TIWmrcP+iyuzn5TEocOsTekoIUQFcCywGPgceyOF4/wfsBEwDNgDXdbehlPI2KeVMKeXMmpqaHA7lkZFYSJmZMlUdSEZV1JPbKtnSVlpTNr4jKRn+xH+wmkOMnLEGbfcjYOdein2xkmq1Ly0oqoLyERCsVMI5jaC44vUrueL1KzMOL9FJagHC5hAoG6F8NskYJGIF94kIIaj97qn4qw2Sr9ytFh3ZojmLj2iT+31SdeUyBSmYQeVntBLZz8tjUJBWCAkhZgkhaqWUv5JSNgElqKoJDwE3ZHswKeUmKaUlpbSB24HZOczZIx/YljKFZNQ0nJYL2ZTHiTiO6SxMLMWfraLsoxXUTA0R2GUXmD7f/fFyREhLaSYAxTVKQPiKC1sMVTMgUKoEXcARRsloxvSanqAXBRl33cWUjYrAy7/P7YbfHtEWybwtOHXlpJPUnAahOQmsDdnPyWNQkOlq+wsQBxBC7Adc47zXDNyW7cGEEMM7vTwOVRjVoy+ItjiJqRkERSykTFRuAwMSUYiH3Jn4OtG241jWnHECVUcfAPteXFgfgZTodhRNWlBUrUrUZNMLKR8IrUMY+UshGVZRioU6XMVI5F4X0rhoAw23/D63QXS/0nDdRgGaqSCFDELPLFIRdYkMFRc8BiWZhJAupUwtUU4BbpNSPiyl/AWwc7odhRD3A28A44UQa4UQ3wSuFUK8L4R4DzgAuKiH8/fIBSsBYRdakJ10SvW7DEZImeGyCMcWySS+LQ3oMkFol10QM85UJrECodkJdBklbhTT5q/OqgVFQRAaBMqhdIQS9ImwikAsBKPn0NY2jk3PryeyIG1cUdekFgYRFyHYoPxiwlC/tYzN7/yqHYUXsr3dkVEICSFSwQsHAi92+ixtUIOU8jQp5XAppSmlHCWlvFNKeaaUcrKUcoqU8mgp5YaeTN4jR8L1yjmeyTGebTBCzNGu3IZwA0Off4md/u8O5r3/a4rjOeazuMHRfiQQ8dUQ18uQ9KOILN2AkhooGqI0zwKUtRFCMPznP8YoNll//9Lc6ssZfqXputVaTH9HdF3acX3qnLNJjvUYFGQSQvcDLwkhHkdFyL0CIITYGWWS8xhoJCLKFJdJu0k6nVXdmqmSMTWu7t6sVbziM4a88TYV41qJlVUR9hUmQkrYSQwZIW4UE/ZXY4n85BztOWIv9hyR5zJCvmJlHtSd0Og8awZ6WSkjLvkO8Y2NbLrjkdz8Q5ovO7Oc6yCFgAqc8HoObVdk0mZ+LYR4ARgOPCtl+xWhoRJXPQYSbkOypQ2tW9wHI0iptCvNdF1PTm9tY+SjT2FUSMr3SLJw9NcLopnotvKzhM2anKo2pOOwsSrv5t6NOZbG6Q7N0YqiIdXfR/NvpV1m07+oq3bpxVN2ZcjxB1L/8PNU1q4kcOJl2c1PN5xFSgsEKjJvr+lKu4u2qKjDbrdzOrB6PYe2KzLaTaSUb3bx3ieFmY5HQYm1OE3oMvQKioZUcqpbf0msBWwJpnuhNeKxpzGiYUbtW8fbY88lapS529ctUqLLKJbmJ2qWFyRRNGIVuBJ0oFSZqdq2QNLKOtgjHTVfO5LiqnoCoafh43/D+CwraOkBJVSMoLt5GUGlPfmK0y9uzCKnSV6pi86+HoOBAsaievQrrKSjBWUww2VbGSEZV5FN2dwgbZtEVTlFM00+2/lA6orz3CJb2hgyStwoIWJWFaxSwTX/+y3X/O+3BRm7HcOvKg/oZl6jx4RpUPzVr8OI6USf/zty84osB0BpvuF6dybD9oZ5GYIUhFAlpNoK6B/06Fd4Qmh7IdzQ0VQsHW0N7oMRpFTjZmGGA0DX2XjEoSw+4gd8XNNVecEeYCXRZYyIWUFcLyO7ifVTNF3lMfmKVfRcvvxEQhAdcRKf/3cIDf93Q+bggW3RTRXJF8uQC5TC8HUkPqfDDKgAhWzn4zEg8YTQ9kAiqrLdM2lBibDTpsFlcEGsVYVxu8whEokEY+57kGnv3kUg2YJEz68fyEqAnSRiVpPU+jj0Ot8IoUrbBCrRZVQl2uYB/667UDp9ZzYv1ogsfSf7AYyA0oTd5jgZQRWkkKmwqxFQIduFClf36Dd4QmiwI6W6mPX0JWiQlmOucymArIRymrvdHhj27AJKPvmMoY3LKYlvcr2fK5JxQEJpbd4DEPoVgVKiZiWajOdFEAkhGP7DczCqKlh/58Lsw7azNctpmqo6nqlFu26qKDk3bSE8BjSeEBrspIIRMjl5oy2AdF+poK1e2fhdWrtKPv6UqrfeoWrXVr6YfCBbinLo+NkdVkyZEEuGZZWjNFBJakEiZjWaTOSlVYReWsyIH84nvr6OTddeD1uy9A+lBIZrs1xAbZvM4OMygyrRtYCVJDz6Hk8IDWbcBiMkYyoxNdN2KVLRcy7NcEaolZGPPomvIoE1t5aPa3JoKdAdyahaiZcO7fV2APNGz2Pe6Hm9eswUluYj4qtGw0Kze178s3jyLgw54QCM5Grky3/I3h9jBJTJ163AaK+QkMbcJjRVccELUhjUeEJoMOMqGEGqi9ytVmMlIJqFwAKq3lyEHo9RvU+cxWPnk7efXSKqQoWLawpbdLQbDhg9jwP6SAgBWMIk4huCwM6LIBo6/zhqLvwBIlIPr/8pu53b+wg1uDPL6YYK689U4DQVpOA1vxu0DH7bxfZKIqJWppma1UVaVCKhm5wg6VTUzjIabvOB+9E2YSxrawzieob5uCURVavv4iG90penK5rjfX9jtIRJ2DeEongDup3A0npYDaJmPK1FR9D80GuMqHkCsfvR7vdNhZHHWlQ9vEyYASd3qChD7lBQafRmkdf8bhDiaUKDEduG0KbMQQNWljlBsdaszHCxNRtJfroYXVi0jRpLi3+Uu+NkIhFVN7A+FEAA1y+6jusXddsSq9ewhUnYVwVI9DxoRMniCbSsLqLh/kezDwzIJloupT211pE2d0gzAFuZjD0GHZ4QGoxEm12ETkvl9HVTyBQ6RcO5E1h2NM7aX93C6qtuZ7eNT7ibtxtSAqiobwVQf6OzIOqpaa784LmUzp3E5mVlRNZmKYSyjZYzTKekTwat0ihSpj6v3cOgwxNCg41kDMJ1mds0xFqdG7qLSgftZjj30XCbbn+Q+IZGhs2WfDI0T4EIKROcJ4C6JCWIRA8FkRCC4d87A6OyjPV/uAd73fLsBsg2iTUVBZeumKoQqnqElzs06PCE0GAilROk+dLfpO2ECkZwm+MTa8kqKTX0xlKanv0fVbu1UnzqRSTy4QdK9r0PaCDQIYh6FqyglxYz4uKziK/fTNPtv4HlT2Y3gOF3zHIu8o6EcNd3SHfaPXi5Q4MKTwgNJlK5F2nruMmO0jyai39/e204dwIrUd/MhhvvIVAZZ+gZR0NNHurCJWPqBuQJIFfYnaLmsHLPIyqevAtj//BDKvefCO/cA3VZ1C3eyiznQnNJ9R3KVNLHV6Q0fS93aNDgCaHBgpVQWpArM5zL0jzS8RtlEQ2n+QyKp+3CiNP3QEw9zt1O6UjGlN+quLrfCaBDxh7CIWPzmPOUR1Lh29jxHvXnCe46FrH3hcQTQ0g8c112odK6qQSQW83FCCoNPZ0GJzSnEKrXhXWw4AmhwYCU0LpJmTTS5cvYCeXbcZ2U2uR0SnVnhpPJOHppMSMv/Q7+w77t7hjpsOJOJYS+yQPKxN4j9mLvfDe1yyOWMFUVCSvWIz+KjZ9Vz1ewfqFEvv7H7HbWA0pwuQko0DS14GirJ61ZzkhpTS59Th79mv53ZXtkTywE8UgG7SZlhhPuzHAJp7Oqy06pkU/X8MV3LyXx35vczTkTVkIJ15KafpsbsiVaz5ZofV9PIz2GH4pqwIrmrDloPpOaM48hvNlPw8dZ9n0SqOoI4Xp3GpnpV4V0M2lcZpFKQ/DMcgMeL1l1oGMllBaUKdcnZYbzFXPvm+k7gQqSFMe2YAsDKfQuu3N2xo7EWP/rP2G3RdCqR2d7Bl/GSqobVmmtkyPST5CoQq+2DdLm1kU3ALCTNt/ZQCARSKEhRT8SnL4gMEQJAqMoq0Tj9i6u1WMYNWFX5DMreHPUYqIjhrLtGrbb34mmg0wqzbrIRcdUs6hDY+9OCxeaMhO3bYaykf3OVOvhHk8TGsh0NsOl0xbao+HcmOEk/kQrIFzfSDfedAfxLW2M+Oow9BnHu9qnW2xLzbdkaN8XI7VtpRHGI6qWWjICSFUM1l9CWPcT1v1EjVLiRjFxI4itGQhpY9iR9oduJ9w55wuJrxgClZAMp7V0dYsQrD/6cJLFxezw0EPM++h3GFaGIILO6H71HbqpSSeEMstlipYzA05lEC9abiDTj5aZHlkTa4F4OENpHqc2nNBdmeHUjTOMpbnzGzU/9xLNr31E9TRJ8Sk/cjnx7qZqK/9F8dC+a+1sW8oXJW210g6UKZOWbjpaWceKO2KqVXpC6+r7t9FlEk1aGMkohu3c/BNO9GJfrNwDpYClEkPdVsnohF0UZN3xR1Gz5C3KEq8yc/0/eHP0ObhWrfSASjjV/ZkXGKZf/bYjLRBMUwLILII2Jy/Oawc+ICmYJiSEuEsIsVkIsazTe1VCiOeEECucv5WFOv6gJxlXNnFfcfrtIi2uk1I1mSCQaMYW7i5mKSUNT71CcKhF9fcuzjyX9IOhy5gy17htqpcvUsIhHlbaT1EVlI+GitEQrHRqlplk16VVwxI+ElqQiK+SUGAEYX+NEmrJmDpWuuTMQhGoUGHOOVYeCI/bgS+OO4kPRh7HyNZl7Fz/ovudNU1p7G6rYptBiNSrRUF3pMxyoY1etNwApZDmuLuBbXs3XwK8IKXcBXjBee2RLVKqi073pY8as+LqIna16rUJJJoA3bUZTgjBDr/7CaOu/RWiZidX+3SJlOgySswo7Zkgy/q4jrnN8ZVRNgIqRoG/zHVEoHsElvApoVYxGkqHqf9dPAyJLBvJ9ZRgldIarNyPuzYxkU8WjGH8qiepinzufkfdp8ytmapngxNE49SWS2fONPxOO5Im9/Pw6DcUTAhJKV8GGrZ5+xjgHuf5PcCxhTr+oCbSpC66dEmp0oaQi+oJDj6rFU0mXVdhbnnqMex3H0EL+DCGuHA2p0GXURJ6cX4qK7hBSnQ7hiHjSjOpGKPykNIm+X6Zk3Y9iZN2PSn74wtNaVdlI9TDdPwlvRXpJYQ6X/T0WkYapG6QaBCsf7OG2uZlmXfojBF0eg+5EIKGqQJVok3pt/MVqyAFr7bcgKO3AxOGSSk3ADh/h3a3oRDiXCHEIiHEorq6ul6bYL8nEVU28EwaQ7RJldoxMgsV3Y7jT4awhLubcOvrb7HuLy9Q/9iCnG9iHceOYWkBYkaWob89OJ4hI8SNYkKBYUozyTECb68Re7FXT/OEDL/K5SkbqbSveFuPqhy4RmhQUq206hzMgvHqKjZ89VBim3XqPnDRtmGrY+MknG5xGbYdUAuvRCTNmEJV9Qht6FFyrkfv02+j46SUt0kpZ0opZ9bU1PT1dPoHtqUuskyO7URYXbQuzHCCJMFEgxJALjSmZH0D62+8F395kiEXXpy+D0wGdDuOLXQiZgXZ+VuyR9hJTDtMUvPT6q8lZpQjexiXszq0mtWh9OHurjH8KiS9dLgSDPFw4X0cmqGiEO1ETjfu5qmTaNxjCtWvvE7NR2/Dpy+431k31L887KI9g0AJmNbNamHV7Zg+FULf1s9ztzy2orej4zYJIYZLKTcIIYYDm3v5+AObNsc2nq6Om51QZjgj4OK+LgkkWnAbji1tm/W/ug47bjHyoiPRqsdlM/utSBXXjPgqKehaSEoMGUWi0+YbSlLLzuSWjhsWqzyhCdqFeRsTMwjlI1QEW6Rehd+7qXSeK7qpohHbNoMIZh21t/GIgwiuXU/t66+B/EwlF9dOcXnsVGJqK/hL0m+r6Upja6tXrdy7+3GbxcoK4Atmbujo0S/obU3oCSCV2TcfeLyXjz9wibaoR9racKlwbM1VlQHTakW3o1iaO22m4cEnaVvZwrDDRuHf8yiXE/8yQloIbMK+qh5rI+nQ7ASmDBMzSgkFhuVVABUUoamw5PLRypwabyusVmQGVERgMpJ1DpH0+Vh91il8eOY5UDoSXrneffQbKP9QpN6dP8wMOFp+ppbgRSpwx6umMCAoZIj2/cAbwHghxFohxDeBa4CDhRArgIOd1x6ZSMacqggZipNGml2HY+t2nECyBUu4D4cuPWAfqo+eScW3fuh6n20R0kKTcSJmFbbIdwSag5QYdhQB7aa3fmx57h7dVCa64mp18y3kTdVXolpyJ9P4XbohWVaKZRaTnPE9WlYKeOla977CVARceIvLattBiGxR1eK7Q9M7wra93kP9noItQ6WUp3Xz0YGFOuagJNWqW/Ol124SEadVd+YQ53Y/EO4i56yWRrTwGny1U6g5Z37G7btFSjQZJ2pWuNa+ssa2MWWEmF5C1BygwmcrhAoZN4LKJxKPKK2oEMmugXLlc8lYh7Br6v+9lIbXS9B9ayje4WnY3WUVdd10mjE2OFF7aRBCJb2GNkH5yO6DSgwn4jDiYkyPPmWgX6GDn7Ytqhx/umxwO6EuSld+IAgmmgGBdBMVZidZe+mvWf+rP/bY4Ws4uUBJLYNGlyvJOCSjhM0hRM0C+5p6G92EsuEQLMOQEYQsUARYKofITfj0NtSccRS+EUNZ/85okrXzstvZcPxDmdp8Q0e1hUzVts0iVS8vm/YTHr3OILpKByHRZuVkTafdSFsl87n0AxFtRrNjrjWR3Z6+l/DqGMVzZqimcjmi21HiejEJPYMDOlfiEfUdlI8ioRdIyG3D1yacwdcmnNErxwIcX1EVYbMaXSZUTbq8H8PJIRJa1uH3WtDPyJ9+AyuSZP2N/0CGm6BxlfsBjCBEG90JQDf+ISHUteP5h/o1nhDqrySiyg/kJh8oGXcXQRVXxR7d+oFGvrcAbdFGyieXUnFW7v2BdDuOpflURYR8h2JLqcwu/hKlKeS90kH3zBw2g5nDZvTa8VIk9SCt/qHYgGEXIDlTaFBcA7aNSBcS3QWBcSMZdu6JtC35iPqbfgMvXu2+wGjKP9S2JX0odgqzSAU1pMsf0nQ1Zmijlz/UT/GEUH/ESkLLehXCmq4sTyzkOh8IK6mcv3rAlT+huH4dlU++gVEJtT//mfu5f+m4cWwEUbOCvP/cbEuthoNDOlbvvcinTZ/yadOnvXrMFLYwCQdqSGhBpzhqnqPndJVDpJHI2vRXeeheDDvnBCpOm6/Cr1/+gzuhkjouQpnaMp1TKn8otCl9wq3hUybttjqvvlw/xBNC/Q0poXWjep4uETQZUxeVEcysXEhbbasZ7hraAVZUkCwq4otTT0QrztGE5jSmi/qqkOS5v04yoTTA0uHpqywXkD8v+TN/XvLnPjk2gEQn4qskapQXxk9k+IialaqwbJatKKqOnocxdnfkrPNIrP4IFv81q+NiJdzVgtN09WjdnH6OZrFKcfDqy/U7PCHU3wjXO2HWabQbO9lRwNSNUGlrcN2mW2BRHN9CdOQIPvre92iu3TmLyXeeo9MXqLgGW+Q5CDMRBaSKjsqhJcHgQhA3yoiYQ9BlLGvzWSaSWpCIUe4Iouy1iHX/+ozVr4/F/uC/8GkWFbeNAMRDEHfRs8hwiqJm6j/kK1GLMTc9jTx6DU8I9SeizSpMNV0+UEqrAXf+j2izyv1wWZxz4rN3M/XRmwnGG0DP8echZeH6AsUjasxe9v/0dxJ6EW2+oegk8x6wkNRLiOsl6DJ7/1PlEfsQr4+x/uNJyFFZ+M8ESssP17sMVAgq01/GQIWgMnXnEP3nURg8IdRfiIeVduMrTuOzkUpIJeLu8jjiYSWEXHVUhXHv/xv5+maawzVEjBxbPUmU0Cuqym9fIInTwK9IFfzsT22/+wmW5qfVr2oC63Z+b7JxoxRLC2Y9bvHU8QydfzShZfU0/HuR8k26DfUXAjS/ClRwU9TVLFaBCum0J81QgQot63unUKxHRjwh1B9IxtVFYRald65HWlQ3VZ+Lm3sy3ikQIfPmNRuWUvLku2hlOh+d9g3QcoxiS4ZV0qMvj6HYUqqeP8FyFbXVywEIAwlbmLQGqrGFnufIOUHELMcWOrqdXbhz1XEHUrr3NDbf8zit9/4Gnr/cnZkNnEKnQgmiTH6plPbUujm9pmP4nNQGr6JCf8C7mvsaKwkt61SZkXSr+3ir06DORdM3O+kEIrjzGRW3bWLkP5/EtjVWnvk1rECOjeUSUfCVKiGUL6TsiIALVlHoatvZ8M3J3+Kbk7/V19P4EhKDNn+NqpUXz74MT/doTsHZjgK0bhBCMOL7Z1AyaxLG5AOUkHj591lEzKWqY2/bnqyrKWoqqjS0UfmJusMsUs0E27Z4EXN9jCeE+hLbVhqQlOl9J8mounDdRsK1bnHKm7gzWYlNbcRbTDYddwiRmpHu57/VHJ0me8GK3PbvCttWAqi4ps8i4NIxacjuTBqye19Poxs0wr4hKn8qHs66MGl3SAzCvioE2eUQaUE/oy87l8D0fWHOd7DXL4O3bnV/YMOvzLyZmtuBoz1pqpp8uohBX5EaL+xCuHkUDM+w3lekWnTb8fTajZWAlo1qdZdJq5HSya+wXAUi6HYcTcZo3XFHPvnB+STLcix9n4yD0FVFhXzVNLOSKrihtDZz4dY+Yln9B86z/hqhJ9T/RNOcfLLivCiStjCJ+IYQjG/Blu7agHRm00tNRJfsxhgWIIprYfKJ7nY0gk5bcB0CGX6rhq9j8VYyrHsTrq9EBT/oRn41eA/XeJpQXyClKrUTb0svgOykamKnGe60Gjdtvx0EFlMX/IVZC25Gl4ncBZCVUDe2kjwmi1oJVTKmbES/FUAAd75/B3e+f0dfTyMDQpkxg0OUXy1PpidLmERyzCEK7LID4VUtbPh0CnLUTPc7CsAoUqV94uHM2xsBFcSTLnRbCBUMFNrk1ZjrIzwh1BeEG5QZIF0jL2mpC0PiqkU30WaVV+EqIk0y6e17ib7SSsPqKqxcFWIrqfKBiofmL1otmVDjlo1I37zPIzuC5VDktITIlyDSAkSMCowsc4jK95tB9amH0fzuFhoWOBUnUgnamRCoYJtwnZMvlgFfAGJtEGmke0GkqcVOy3p3ws0jr3jmuN4m3KjU/3TRY9J2MsCT7oRKNNRhbnHBrh8/ifzvOqgM8MnpX8/NhGZbypRYMsy17ykTup1Qwrd8pJcDVAgCZU6kmVNpw2X1jHQk9WKi2PiTIXWPd/lTqj7tcGJrNrH57icwY59RxgI4+JdQvUvmnTUNCKjzKBmWORfNLFK9ttC691lqutqueS1UjMlveoFHWjxNqDcJN6oLx1fS/Y0/VRU7GXOnCcTblHnCLHJ1A0guW0DwsfeQhsEn87+BHcjhYtsqGTU/3UpV2K9UGpAngAqHvxRKhqpQ+jwV9EzopcT1kqwa4glNY8RFZ1A8bTxi1GQlIBf8WvlJ3aDpSvtuq8uc7yNwBFGDI4y6G9NQ11zzGi+ZtRfxhFBvEW12BFCaZNRUNYRExF05mnhEaVVG0LU2E1qZIBE3+fyM00hWVLiff/scpZOMOiRvq0XdjiERtPlr+q8AspPqxpQIK8Efb8NnJfFZSQwrjGm1YVphTCuMbkXR7OwLf/YavhJVcy8ZzZsgihulKtrMjYnMQfP7GH3V+ZTuuycceKWK2H7uCmUFcINuKlNa66bM4d4CJ5l1i8q1Szem7oMmTxD1Fp45rjeINiv/jq84jfNeKgdqPAI+FwIoEYXwZtBdCqD6lVC5A5XHHMIrxSNJVFVkcwbtUyQZgUBl5hYTLtHtGDYa4UA1sr/8HKVUgRHtNzapohPNoLpBOSHA39j7MgD++47mKKESIW10mUCzExh2DM2Og5BIKfJfQ68nmEVK60xVa++xSdUJgLDrlXBz6c8Tzm+3+d11bP7PcHbYfx2+F38FR97gbk66qf5XqSi4dD21UoKobYt67S/rZkzHvNe8FspH5U3b9+iafnRVDFLcCqC2OoiF3Qugts3KQevCri/XLWXD7/5I1aEzCBz5nR4IoDAEKjKHx7pEt2PYQifsr85/le1skbZa+UrLqTFWAv5iR+h03QZ9/NBpAPxb/2Sr9zunSAppockkuhXDZ4edSgHCiXh01169YBgBRxBtUK97KoiEExLeVqfMtbr7m3dgx1HIpM3q18cw9uKDMLKZi+5TaQKtm5WpMa0gEkoAuxVETWugYrQniAqIZ44rJOHGzAIoZYLLWgD5XXVSlZs+Yv01f6b58yAR24XTtzuSEWW3D3Rz0WaJbsewNbPvBVAyrgpfJiPqhlQ+Cqp2grJa5T8x/N0KijfWv8Eb699IO7wUOpbmJ26W0eqvhaodVfFVw69CpuNh95UDCkFKEFmJ9D153JLqzIqelTnLP7qWUZd/m2RTmDV3vI3VFoENS91XvG4vxVOX2cTYWRClNc35lKbVtCYrM6NHdnhCqFC0ByFkEECtBRRA9Z+x6ffX0bLKT82pB1L51YOyPInUcSPKjxCoyG3/bTDsqBJAhegz5AYpMawoptWmVv/lI6ByR5XrZLr3r9323m3c9t5t2R1b09Vvomw4VI5Tx7QtpSH1lQ/C8KvvwLZUiHxPERqUODX+smirXbTbOEb+9BtEP1/Lmsv/iP3CNfDir9J3Tu2M4VeabLaCKF2wQkoTbl7jfh4eWeEJoXwjpfpht0fBdSeALOVQTWbhA8pCAGFb1N14HY0fB6g6ci5DTj8mu/Nw0O2ocjgX5VhVexsMO0JS89PmG9L7AkjamHYE044SM4ppCYxU4eC+4ryEK2dNKku/cgelgemmSpjsi1W37lMakcyTINJ0RxChfDYuKZ09iZE/+TrFs6ag7f8DqF8BL/xSaY1uMPwdKQ4ZgxVER4vwtILI7Iiai7ksvOrhmj4RQkKIVUKI94UQS4QQi/piDgXBtpVgiTSmD8O2k6oUTzLuMgouDG2b3AsglCyMWuMoP2AqQ887vd0BnA26HSWpBVUkXB5QAihI2FdFr/70pMSwo5h2jKhRRnNwBFFfFbbWTyLxUn1uykcqgWQGlDDqbc1INzsJouwqZXeJplqEA1kJorK9plFz6mEweg6xXc7B3rQSnrvSfUUDw+dUJdnsInxbdLSAiDTQbUKrZqhqDS3rlZ/XI2/0ZWDCAVLKLX14/PxiW6rETiKWPhHVSqjtpHTZE6hVVVjQ3SUXyg1LkBs/Q9vjeEZd9X2ErucsgCzNT9TMjw/IsMMktGIivgp6UwDpVhRNWESNcmJGKbI/Rah1heFXprpEpQq/j4XUe+lavecT3VSmuRbnt9xTNEPlk7VuUr/9LELwrXCEL25+CX/tNEZPXor22UKY8FV3Oxs+J2pukxKE6Y7bHr7drBaIRd2UoEqZUkOb1IIzWNG3gSWDBM8clw+sBDSvUz96X5paZ8mY2g7NfSWEcIPr7Ha5ZhF1N97MFzcvwAq1oJkmIgczU0oARcwKevwTkUoAxfUSpw1A7/zkNDuBabWR1IO0+EcSNSv7vwDqjBlQmlH5KKehX1vWNdpyRjOVIBRafprj6YYKn5Z2VsEPelGQYd84jvAnG1n70WzscYdmeVyfU017U2atUqCu3Vg4fd6R0DrahLfWef2I8kBfCSEJPCuEWCyEOLeP5pAfEhFoWq0usHQdTOOtqm+QbrqsBdekKiEYRa5WW3LVq2y68RbqPywmMHUmWnFuTeXyLYBItBHXS4iaFfRKLyBpq4ADIOSvJeyvLpjZ7fI9L+fyPS8vyNjt+IqgcozyryQiKlm2N9BMKKvFphCCyL1prvwrsxn+vdNp++ALVl95C9bGz+CZS9zXmtOdPl2tm9zVhfMFnSaTG7oXmEKoyMlYCEJeh9ae0ldCaG8p5XTgcOC7Qoj9tt1ACHGuEGKREGJRXV1d78/QDdFmJYA0M00egVRqfqofUKb8B2mrwIZoyBFAmachP3yaDdf9lcZPiqk8cm9qv3dmDzSgQF4FEMGK3hNAiSjE2wibVYQCtVh6Yet/jSsfx7jycQU9BqBueoFyqBwLvmJ8VltWvXxyRjMIB2qcbqp5FESQlSCqOGguI388n8hHn7PlkZfVjf+ZS2HLp+6PqwdUMnjUhV/JDABCWS3SRcT5ipTAalrthXD3gD4RQlLK9c7fzcCjwOwutrlNSjlTSjmzpqamt6eYHttWDbNCm5Rq3p29WVpKbY80KJtzJsFgJ9W4yagTKuxuOpuf/oTmVUVUn3oIw847JUcBpIIG8iaA4q1OC4Fe6IZqW47vxITKscTNsl5pAb5wzUIWrllY8OO0oxtQWkvIX4tOEtOOFLwrqERv1ybz0i5c7xSskEXwQ9m+M9jhN9+n5ltfg0N/q7oGP3c5rHvH3QCa5vQjanQqamfAMJVfKbTBiZzr5ns2g8pX1Lza6XXkkS29LoSEEMVCiNLUc+AQYFlvzyNnrIQq5xELKZW8u5udlYDmDcoE4HOh0STjTusG213JEysBdSsAqJo/n9oLTqPma1/NPghBSnQ7QkIvJmqW02OBISUkWpVzN59dVrsjEVVh7qW1UDYyc0XlPHLPB/dwzwf39NrxUiT1AC2BEUT1Ykw7XPAadRKdNt8Qkpofw85DroxmdDSayyICsGjijmg+E0urZPU7E4nGh8HCX8OGJe4GEEJZF2JtTgh3hu8tVVk7Uu+0Ae/G/6P71LihjZ6fKAf6QhMaBrwqhFgKvAU8LaV8pg/mkT2xVmj8Qv140wUgxFuVoAKXIdhtysYtdHc30UgT8X9dxqbfX4Nsq8ccWkXloXu5O4fOSIkuo8T1EmJGGT0WQL3Zjlva6nt2tJ/2NgXbCVLoRH1DaPXXostEfrSUtGiEfVUktCIMOw9+qfY8ouwqKwAkm0LE1m3hi6c0Ws15MHSi+50FytxmJToi9tJu7zS9i7el9xNpurKKRFvUtZ+PEPfthF4XQlLKz6SUU53H7lLKX/f2HLLGdiobtKxXWkp3kW3SVmG1rZvVdpkCEKR0KitscecvAmj4nMg9P2XVQ600r64g0ZLjqkvaGDJKzCglnhcBZKnaciVDu6/HlS+SceVvKqpW2k9/rbzdC6S0ooQWxLTaCqwVaUR8lcT1Ekw7D83xNB1KhypNIotwcP+oYYz9ww8xh9ew5r5PaXz+bbX4eesvWZT58QNCaS9u9jGD6nyb13SfrySEWpxKG5q+8MxzLhlAMat9RCLq5AUk0iegWnElqKyE8v9kuqdbSSV87KT7itRfvErTvXew8a0SjOoKRv/ye/iGZ+8vE9JClzEiRgVJPQ/VsK2kOv/SEe40v56QaANhQsUOXlFJB+kUgU0kAxTHG7CFjqUVyiwpiJoVSKHhT7aQxH2Zo66H01StuUiDMl0b7nyh5pAKdrjmB6z73V1s/NMDyE2fU6U/BxuXwbyfQtmozIPoptLEwvVqYZMp78cwwdaVnzfVzkR0kTxu+ME2lIBLhNViKU+NHwcjXp5Qd9i20lKavnByCLrrA+REvzWvVSsln4uLKN6mHJ7YWfXkqXvoZTa8WUpw4k6Mvf4S/COHZnFCDlYCTSYIm9V5EkBO4cuy4YUVQKngA1+ZV9W4GxJGCS2B4dhCw7Ty18b7ywhiRjlRoxxDRnqevySEuqEHStUiw+W89aIAoy8/j+rTDqf0yKNUZ9ZYCJ7+Cax6xd2xNae1d7xNLTYzmec0zTHPRaBpbffRc5qufMbxNhU951ZD2w7xxHNXJKJOyY84+ErTaD8JZ1UUczqbZpA+0laRObFW9+2VWzerH3DVOIoPPQG7chlDzz4WoedQd82Kg5REfNVYIg8mrGRcnVP5iMJm9CdjShMtG64u7H7Cb/f9LQD3v95/zC62ZtLqH0Yg0UQw2ayEt8tST9kSN8qQQqcoUU+SPITDBypAGEorclmiSug6NacfAYC0yti4bh4VVcsJvnajMoftdqS7Y6f8RKENKqLTnyHPzgwoC0Bog9KgAuVda0VmkbJ2NK9T51dU5WlF2+B9G52xLYg0KfVc96cxk0mItKioGc1MH6SQIhFV44J789vq12l79HaioVKG/PwmiiaPp2jyeHf7bksyqi7w0mossT63Mboar6xWfQeFQEplztB9jvmt9yLf3FBbXOs86z9CCAChEfVVkdQCymyUNo+tZyT0YtrQKErUq5tyT2+w/pKOtt3SyMrfl6xvou39z2lqiDH04L2oGjUnO09nKrE10uCY0arU6263N9Tn0Ra1sCyu6doaoBmOVhRSj5Kh6U372xmeOQ7UzS4WUpFvUaf4aHc3vGQUmtc7uT9FmW+M0jHrtW12iiC6uBkkIsjXbmbzn29j9fMlNK+rxk7kmJwo1XgYAeUETndRuSUecaouF1AApdobpHr89DMBBPDM58/wzOf9N7AzaRQp4S00p+xPYcxzST1Im3+o0iTyUYHbDKiQe2RWkXPm0CrG3fRTSmZMZPN/VrH2xkewWlph4e/gk/+6GyRVWdtKqmi4eIaq2QInp09XPqC2NNW7zSK1uG3ZqKqneO3DAU8TUj+EtjrlFDWD3d+k7aTSkmIhtWJyo/3Ewx2JcS6rHxBuIPbQz9mwwCZSX0r5QXOoPe8kNF8ON3vbBiuq2nHnoxuqRK0Q/cXdF3nMB8k42HF1I8pTE71C8M+P/wnAFP3SPp5JGgyfEuJtW1QpKLOoIOY5S/hUrbvQRqX1Z+Hr7BLdVLlE4QZnEeUuYEEvKWLUz8+h8cmX2PTXx1jzy1vZ4bAI4u3bYN1imPNtpeFkIlWJO9yoruNgZXqtTDfUIx6B2Brl4/J30cpF09X7yRg0rnJ8YRXbtYlu+z1zK6F+YNEmtarvztcgbVXqI9LgFC90IXwsR2Alw1m1XwCwZJBVT6mOjiN+fBrl+81wve/WA8WVNlFU465fUSZSprFgpZOEWiBTQiKsVpUVY7zgg3yhaUoLNgMqH00vUFVu3Sl82lanbsY9/d2lcomizerhsp29EIKqo+dRtPtO2PEEYrdx2B88CYvvR9v8PZh+NuziosFjqsWGFVe+n0B5+gR1UN+xlMr0Hm1WkX9dmehSldGjTepeUVKj/M990deqj9n+hJCVVD+OcH1Hgll3UW/xNmird9ouuAhFlbbSlKLNyl9iuhBYAGveJvHaPzCO/gV6cRUjLv4GgV3HYg7JIeFTokyGmuHkYOThX2xbyrdQXFO4HCAp8dlh9Z2VDC2YM327JlCmbn4t65Wwd/v7zIZUNYRUGwo36QqZCJQrs2+4HqTuWoAGdhrd/nzLIknLyztRu6dFyXsPwth9wAzw8Dvr3M1BSk6YoinfT7Ay/WI0lS9kJZVmaAZUsMO2i6pUL6NUh2WtQQmt7cxftP2IXdtSmk/jqg6/T5cRbc6Kv2md+mHophN2neZHIaVTJXuDU8PMha8IILQR+7mrqfvzn1l5f5LQ66oOVumeU3MTQLattC9fMZQOy4sA0u2EMo+VjSyYABLSwme30WZWKROcJ4AKh+F3tMyg+s0Wwk+Uyv0JDlEh1/koY+MrgtLhgKbMfVlOu2TWJIQ/wJonWln70R4kQjGwEuy65Xl020V1AyGUMNEMVQi1ZWNmn45udAijlvVOhYYujiWcsG/NyS1q/ELdRwpcF7C/MPg1ISsJsRZlW253InYle6WyPYcb1Q9F97n3+0SblJDT/a7VabnkH7Q++282LSkjESqlbL/pFM34Uh1X9yRj6kdbVO1u3i7Q7Rg2mrL1F6gqgW7H0KRNyD+cpB7YrlaAfYamK7NZpEmZztL5QnNGqNJNuq7SDHRfz39DuqEWV7EWZW3QfK4XWkW778S4m39K/cMvUP/gf2k7/9cM/9psJrc8wU4NL/HesGNZV74HGdflmg5akVP2x6mMEihPv+g0fIBTFSK+VgmcQMWXNSPNUItjK6GEkWYon5GvZFCb6QavEErGVOhkpNEpRNiN8JG2MrtFmlTwgWvhE1EalZ1Uwsd0oflI2X6T3fDA+zR/UIlv+BBG//hUSvbYLbvz6zxm0nHcFlXm52YiJYaMkNSKVCfUAgkg0w6TFD5a/bXYeb8JFp7r510PwN2v9NNWI+kQQv1eUuY5zSqMD85XAuU+dVNNRt0V501Hqq2FEVDBFkn389ZMk5pTD6N8v+lsvO1hfFP2YuGaYUxZ+yBz1t1DY/0LLBt2DHXFLtIgdFM9rLjysxlFjrkzzX3A9AN+RxitV/ekYEVHCaFtx7aTKoFW26JMgIGyQWklGHhXfjqko81EGpVgSefzsZOqmm60yalc7XcRbu2Y6qJOG2DN586ubifh81eJv/EvjIN/gDZsZ0qOPJHA3q1UHr4Pwsjxh5WMq3YRRVXpW4pngZAWhowSMSqIG6UUJABB2ph2hJhRSsSs7JW2C4WgMlDpPBuAQiiFr0iZ51I11NzmsGWD7kTOtW5xolBdRoqmI9UGPdKscm+y0Ip8I4Yy5srvAFBfL/lk6W6UNJYzfOJGducpFu64K64nqPvUI9VK3Ah0+N66w3Q+S8bVAkD3K63RLNr6WtAMFUlnW8ofFt6iNKhM4w8wBocQshLKYRhtchLmuot2c/IOYi1KAAmhfjSZTEC2pS6eWIsSWLrPneZjxWHliyT+9yhbFido+ryYYdo7VJ2+M2X75hj1Bh2h10aRWknlKbxTt+MIbNp8w0hqhfmRa3YCnThtvmoSRn4EZ1/x2KePOc+yqOLcH9k2jNtXnP+FgXCKlaaSvI1Az1f1QlPanBl0EkwdTSsbAScl4R3GULRqDZufCxLZcShB3wbs2hJmrLuXFdUHUFe8GxkH1X2g4wijzSqQIlCWPqApZaazkkpAp7Q8X/HW1gdNV+9J6Zgim9R5BiudkPuBuYhLMXCFkG0rM1SkySlqqakQzq7UfTuhesfHnDImQne3GkvG1Oow3uqM73fvs5A20b9dRMM7UZpXFyG0AJVH7EPpoYdke6adxkQJH0RefT/K/BbF0nxEzBpsUZifhWFFkUIQ8g/HKpCQ600e//RxAKboA1wIQacwbr9a0RckjNvxExk+pyxW4kv5RPe+udr1aGfMHaOemAEwap3I1CYlANzOXQgaZ0+nadpkqt5+hyGvvsm42/9GyyHTqahZz36r/496/xg+rj4E5PDMwrldGCUdP7RQAT1mmir5qRwjKdX9LNLQ4WsyAx3lgFKJtKCEXcuGDsHlL1Xa0QD0qQ4sIWTbyq4cCznl1KVSw31daD3SUqa5aItTYsYRIpnUWDupVlSxlg6B5TbRtOFz5Oo3ENNOB6Gx8d1qoutbqDxyb4Yc+xXMoS6S5LrDiithGihXP+o8rVSFncQgRtQoz09Poa6QEtMOk9CLCPuqkAUSch55IFDuZPWvh3z0DeoKM6jqDbZt6UgS7+nNU2jOTbtImeMT2eXoSZ9J/d5zaJw5jcpFS2iePJENpQewwwfPMvzjxcxtuxMefxK+eqM7AddZsMRalP/YCDpRud1YX4ToyK2yEio6F1RyuK/EETLOdZ8yA6aqvUSblPkuUKEWp7pvwAikAXE38Gu2ctDFQ2BL9WV3FV5tJx1zW0j9CNGcMMkMdm7bUvbZWAtYMfWP1kx3JrdkDFa/SeLdf9O8eCPNq4rY4XezMEbvwvAffwe9vBSjrAd2diuhhI8RVHk6eQwS0G0VYlpI8xt2Ep/dRtisUkJugFwY2zVmQPmJWjdjWm0kNBfFebNFM1U4fqQFIlsczSsPv23dVHlmiahjootnFbVq+/3U7z2n/XVkpaDxHR8NvtFU7VFNxdQ6AuNGwgePQs14qJmQ/rsRomPhayU61Y8scsp+daO9pIITJE4VhlZHqyrttJ/WkVALnXxHTg5koByfXqis8vwxIITQ0EDS6TWy7cUgnXpVjnaUjAGaE+qYSfA4AiseVmY9tI52vi6xN3xM699/S/NnBq3rAyDLKJq8E5YswQD8o2szjtEtVlIJH82E4qE9L4PSCSEtDFt1VI36ypCF+hkkoiCtjvBrj4GDbkDZcCLmOooSDSS0ILKrKtE9wjHPmQG16o9Htoog7RFmAIzhylQfaQZL5mSu2nD04TRNm0zloiVo73xM4/+uoWzfqYwc+5JTQaQaxu4NY/eFqnHpB+ssWBJRZepPdW41gl3PT9ARyCBxgqmcgrm+oo46l5rZ4TsCJZDa6hlWnPd/Wt4ZEEIIITtuwqnkyURYCRAp1T9KMzMLnmS8w89jx8lK8Fhx2PQB9mdvYunVmHuehKVXse6VUvTSAEOO34eKQ/bENyKHHj9bHSfhRN6ZTlXePN+8E1F0mSDsqyahFyBjHpzk3TZlWigZRVL/vDDH8SgsQhAzy7E0P8WxzUgEViEWE4ZfmeeizZhyFZb0YeejMK4QHUnpsTYnsEgqU5VbZ74QRHYYTWSH0cza+UyaX/wfemkx7H8m9qevse7/nqZk2XOUjnoCY79vw84HOddvmluroCMSt13TaVMvzKDjP/J/2Ye0rUBKJiBeB9jqfpHyCxk+dV/TfWj9Xg8aIELIr0kVnhgPKwmPUKpoppWNlVDCIxFR2lIqAVk33Ws8n71EYtkrtC79lLb1Bm0b/QTHVTJmz5Mwh9Yw9rofEthpdG79fbaaq1PrrVDCx0oqU6O/hFBgWOG0Hzupvu+iISp0fBCb32456BYAbl2wpo9nUliSeoCW4AiK4/WFM88JDYKVtPprCcYaMOwISdHJB9LTsQOlyreSCDuakZ11BXijrJghx36l/XWieCKxyKu0vh1n4+JKij55j5JZUDY2grnqcRg5HUbOgNop3Td83FYgWQmlJSE78huNYEfbiK32MwHnHGxLpY5IqfZ1koP9hmeOywtV/qSKbtPT9EWRdkeXz0RE3XBTQkfTVeRcpn9H2xbY/CHWps/R584HYN3/PU3LRxGgFKOyiLIDplC2/8z2XYK7js39xKR0Gs1ZSigWleY//j/VykHTlA3eLEKSoTx9riQi6oDlo/IXudePCRoFbmXej5DCoNU3FH+yRZnnRABZgARjS/hoDQzFTLYRTDYjpZa/SEqhOZpRcXvQkm5HAU01ecxSsPpH17LTX35BdOVaQq+9S+viD9l89+MUXXYcZs2uRN55k/jzrxOsSWLuMA5x8C8dc1w3JkdBh8kOHIuC4w8CdR8znAjgVO+j1DiavnUQhm1BpJEhAc8clz/MTj9EO6m+ZDuhso+tmKMhAQj3QgeQa94l9vZ/ia74gvC6GJEtfuJtOuPvOwatpIKSQ48hsHeMkj0m4BtTi+jpClCiTIG2pS4Kf5kyWxWiYkAiqoRzsEqtBAuVFCptpyBmEEpqt5uy9A989IDzbHqfzqPXcMxzSc1PSbwObKtAAS0aCaMUSw8QSDRj2mEszPyY6KCjwKiviIivDTMZxrCVP8rGyEq4CiEI7jya4M6jGTr/aBL1zRiVpaB9heblD9L4P9VmXC+OUPTu31hRNoRdRr5DUbKJxsAONAXH0BgcQ8g/DHubbscnTB+5dQK9lMqiE09FLUongMrvBF84EXnCcTNo5oAwRAyIu0VAOG2xkzElgKSkXcJomlOxuptINisOjauRzWtJrv2M+BdriK3ZTNmZF2CMmUDDf99i8+PrAAO9KEhw19GUT5mIdC6u8oP27vkJSJTATDW7MoMQLMm/yS1FMqEEnb9E9RIqUNkdQH2/yagKnghWDGrz27b8d5VqlDZF306EkIOlB2gJjCAYb8RntZLUAgUIWgBbmIR91Rh2jEC8EcMOYwl/Xo9lCRPLLEdQgm7H8CVa0W0VqGQLI+tjdS48POxbJ1Jx6D5EPlxJePlnRD78jCGRT2n4+gTM0AcEFi5nWHg5O1QmiNRUs2TCfOKVFYwN/Y+oUQatOhQP67imhOjIQ0phW05Ah+NTAicdxQe2hd+LjssPpX7ni04lpHb+Wm0LwpuhbQtW3XoSG9aT2LiJwMwDMCfMJvLGQjb85RHiIR1ppTQBDXPvDZSOmUDpEcdh7DyDwE6j8I0cishX9rF0Iveko6GZQRXDb+bJzt0VKXOk4YeSET2v05WOVAkjzVTdOwslUD36JVLohH1DSFgBiuMN2AUsvZTU/LQGhmFaYQKJZoSMqyZ6eUSik9SKSPqL0GQCzY7jS7ahy0T753Ra+7pB6BqBsSMIjB1B5RH7AvDoayuxgwFWDDmQYR8+R9mmTzDXhIA4O79wO+FRw5m2zzvoSOqfK0YYGmZVKcaEuZgzv4peXoJY+yYUDYXSGvCVwbZlv6RUC95oMxWeOa5rhBCHATehZPodUspr0m1f1hKh7d9/x2psIdncgtUSpmj2bIoPOYH45x+x9oo/kmjTsZMdF0Kt+JDKCbPRR+2CMXIHikePwLfDWHyja/GNqsWoUm0JfCOG9jyiDTr+8bYFSNob4BlBpVIXsj5aKgBD9zt+nwAFSTrtfLxEBIqHKFPfAC8b4pEjQpAwSmjR/BQlGlRlkW3rn+XvYCT0YhJ6EDMZJpBsVgtTw30OkFtsYWLrJkm9GE0m0WQc04o6qRxqLilzV7aXmR3sWKxtOuJgNh1xMFokir++AV99A7Zp8NT4kyiLbmTUM0+QbIqqjZ97E3iTsr0nM3K00sDX/68CTZfoxT70cXug7zQd/4gKAuG3kYEypBgYftleF0JCCB34M3AwsBZ4WwjxhJTyw+72OWJ9OatvX7LVezVD6igGtKoRmGPGUDS0CrO2FnPESMyh1fhGDQPAN3YcY67Jc/tliWMWtJQ/BDqi9QJlXYdXFoJETM3B8KteK4UWPlKqfkVCh8ox3Uf8eGxX2JpJq2+oKiUV3qJuzgUrsKmRMEpIGEVQXKVarySdeo4FuOZsYWALg6RWpAJurIRyCyQiWwuldj9M9oqHHQwQGTWCyKgR7e81FO3IvL9di9UUIrGlieSWRhL1TZhDK2Dng5HNG4i88DhWawwrYsHSZcAyqg6fSaD8CeyEYM3rI/PyHRSavtCEZgOfSik/AxBCPAAcA3QrhF6rjXPaOd9Cr6zCKCtGLy1SlafjbRjFPkZf+f0u9nJyVXqKlI5JTXY0mWq3zQbBMDoSxVICwE52+H/yjqNxSVvlRfk7dWy03B1TZDu3VPi4nVRh18HKQVlS3qMHpFpD+IqhbbNKHtdMR0txvzDK6rdpBKHU70S6NTut4bWsgnw0K+r+eKkmdu1Jp3ZHgFTS8Y22N61L+XG0rb4DPem+FJJIhDGKdYziIbDDkK0/KxnKTrdOBUBaNlY4gtUSRgv4oOo4RHMTVbXL4M8fuz+/PkLIXu7eJ4Q4EThMSvkt5/WZwBwp5QXbbHcucC7AvB203YYWC5d9eLMnnqTMZ9DS1WeWjUzayHhSYoFM2pC0s+3rmD8kEE5gh2JSJraZhy2p0gQN+T6mLZGRBDIUl3Yk2SvnXg1s6YXjFJLBcA6Qw3kIoMyPVmQKETB6r3tzQEeU+rs+ZsKiwtRpKvQchABDIHQNYQowdIFPRwiRHxNFuntVVyyvs2re2yx7ULSy8PSFJtTVP+NLNzYp5W3AbYWfDgghFkkpZ2besn/jnUf/YTCcA3jn0d8YLOfRmb7wKK8FRnd6PQpY3wfz8PDw8PDoY/pCCL0N7CKEGCeE8AGnAk/0wTw8PDw8PPqYXjfHSSmTQogLgP+iQrTvklJ+0Nvz2IZeMfv1At559B8GwzmAdx79jcFyHu30emCCh4eHh4dHCi/L0MPDw8Ojz/CEkIeHh4dHnzGohZAQ4jAhxMdCiE+FEJd08bkQQtzsfP6eEGK68/5oIcQCIcRyIcQHQoiusmF7jVzPo9PnuhDiXSHEU7036y/Tk/MQQlQIIf4lhPjI+b/s2buz32qePTmPi5zf1DIhxP1CiD4ruufiPHYTQrwhhIgJIX6Uzb69Sa7n0Z+u8578L5zP+8U1nhNSykH5QAU9rAR2BHzAUmDiNtscAfwHlbs0F/if8/5wYLrzvBT4ZNt9B8J5dPr8YuAfwFMD8f/hfHYP8C3nuQ+oGGjnAYwEPgeCzusHgbP78XkMBWYBvwZ+lM2+A+Q8+sV13pNz6PR5n1/juT4GsybUXh5IShkHUuWBOnMM8DepeBOoEEIMl1JukFK+AyClDAHLUTeQviDn8wAQQowCjgTu6M1Jd0HO5yGEKAP2A+4EkFLGpZRNvTj3zvTo/4GKSA0KIQygiL7Lkct4HlLKzVLKt4FEtvv2IjmfRz+6znvyv+hP13hODGYhNBLo3Hd5LV/+gWXcRggxFtgD+F/+p+iKnp7HjcBPALtA83NLT85jR6AO+KtjcrhDCFFcyMmmIefzkFKuA/4ArAY2AM1SymcLONd0uDmPQuybb/Iylz6+znt6DjfSP67xnBjMQshNeaC02wghSoCHgR9IKV3Xa8ozOZ+HEOIoYLOUcnH+p5U1Pfl/GKj2pf8npdwDaAP6yg/Rk/9HJWqFOw4YARQLIc7I8/zc4qp8VgH2zTc9nks/uM5zPod+do3nxGAWQm7KA3W7jRDCRP0w75NSPlLAeWaiJ+exN3C0EGIVSsX/ihDi3sJNNS09OY+1wFopZWqV+i/6rqd2T87jIOBzKWWdlDIBPALsVcC5pqMn5bP6U+mtHs2ln1znPTmH/nSN58RgFkJuygM9AZzlRDPNRZlHNgghBMr/sFxKeX3vTvtL5HweUspLpZSjpJRjnf1elFL21cq7J+exEVgjhBjvbHcgaVp/FJiczwNlhpsrhChyfmMHovwQfUFPymf1p9JbOc+lH13nOZ9DP7vGc6OvIyMK+UBFKX2Cijz5ufPet4FvO88FqsHeSuB9YKbz/j4odfg9YInzOGKgncc2Y8yjjyNnenIewDRgkfM/eQyoHKDncRXwEbAM+Dvg78fnUYtapbcATc7zsu72HWjn0Z+u8578LzqN0efXeC4Pr2yPh4eHh0efMZjNcR4eHh4e/RxPCHl4eHh49BmeEPLw8PDw6DM8IeTh4eHh0Wd4QsjDw8PDo8/whJBHjxFCWEKIJU5l6IeEEEVptp0nhMg5QVOoatrnd3o9Qgjxrwz7jBVCLMv1mH2BUzOvxxWRne/7jW3eM4QQmzrVF9xTCHF7mjHGCiFO7/R6shDi7p7OzcMDPCHkkR8iUsppUspJQByV39Ad8+hZlYAKoF0ISSnXSylP7MF4eccpTtpTLga6FQxZ8DIwyqmNluIgYJlUCbQAhwHPpBljLNAuhKSU7ztjjsnD/Dy2czwh5JFvXgF2FkJUCSEeE6qfzptCiCnOjfDbwEWO5rSvEKJGCPGwEOJt57E3gBDiSiHEXUKIhUKIz4QQ33PGvwbYydn/9521HOf5K0KId5xHRmEnhPixc9z3hBBXdRpnuRDidqH6zDwrhAg6n+0khHhGCLHYOdZuzvt3CyGuF0IsAH7nbPemM/YvhRCtznZ/F0Ic0+n49wkhju5iaifgCAYhxNnOd/mkEOJzIcQFQoiLhSrm+qYQoqq7uUkpbeAh4JROY58K3N/p9YHA82m+v2uAfZ3v/CLnvSedcTw8ekZfZ8t6j4H/AFqdvwbwOPAd4I/AFc77XwGWOM+vZOueLv8A9nGej0GVUElt9zrgB6qBesBErcqXddq//TWqNULAeb4LsGjbbbaZ9yHAbagKBxrwFKplxFggCUxztnsQOMN5/gKwi/N8DqpMCsDdzv668/op4DTn+bc7fUf7A485z8tR/YWMbeY1Dljc6fXZwKeonjc1QDMdmfQ3oApvppvbLOBd57kf2IxTccL5bhdk+P7msU0mPqpm2ZN9/dvzHgP/kQ+zgYdHUAixxHn+Cqoe1/9Qq3mklC8KIYYIIcq72PcgYKIq4wVAmRCi1Hn+tJQyBsSEEJuBYRnmYQJ/EkJMAyxg1wzbH+I83nVel6BuvqtRhUZT57QYGCtUteW9gIc6zdffabyHpJSW83xP4Fjn+T9QLRyQUr4khPizEGIocDzwsJQyuc28hqNaV3RmgVQ9b0JCiGaUJgKqLNCUdHOTUr4thCgRqvbeBOBNKWVjp+8g1U4im+9vM6oSuIdHj/CEkEc+iEgpp3V+Q3S6E3aiqxpRGrCnlDKyzf4AsU5vWWT+vV4EbAKmOuNGM2wvgN9KKf+yzbHHdnHsoDNm07bn2om2DMdL8Xfgayhz1je6+DwCbNv2u/N87E6vbdT3kmluDzjHm8DWprjDgVTxzmy+v4AzTw+PHuH5hDwKxcuoGy1CiHnAFql6tYRQZqUUzwIXpF44q/B0bLt/Z8qBDVL5Qc5EtU1Ox3+BbzhaBEKIkY6G0iXO/D8XQpzkbC+EEFO72fxNHE2QL/tO7gZ+4Iz5QRf7foIyCbrGxdzuB85AmUafSG0DTEEV7oTuv7+uvvNdUUVYPTx6hCeEPArFlcBMIcR7KMf2fOf9J4HjUoEJwPdS2wkhPiR9ZB1SynrgNaHCwX+/zce3APOFEG+ibpJpNROpupr+A3hDCPE+qk9RdwIuxdeAbwohlgIf0H1b6x8AFwsh3kKZ15o7HXcTqoXDX7uZVxuwUgixc4a5uJ6blPJDIIzyE6W+lxkoX1FKQ+3u+3sPSAohlnYKTDgAeDrL+Xl4fAmviraHRwEQKlcqIqWUQohTUUEKx3T67H1gupSyuZv9jwNmSCkvK+AcLwM+lVI+kOV+fuAlVEDJtv4sD4+s8HxCHh6FYQbKyS9Q/V++ASCEOAi4C7i+OwEEIKV8VAgxpJATlFJeneOuY4BLPAHkkQ88TcjDw8PDo8/wfEIeHh4eHn2GJ4Q8PDw8PPoMTwh5eHh4ePQZnhDy8PDw8OgzPCHk4eHh4dFn/D9khLI6plFh2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "monitor_stats(T=600, directory='example/VASP_3C-SiC_calculated/1x1x1/T_600K/', \n",
    "              dfset='DFSET.dat', once=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "def moving_average(x, w):\n",
    "    return convolve(x, ones(w), 'valid') / w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "def ewma(data, window):\n",
    "\n",
    "    alpha = 2 /(window + 1.0)\n",
    "    alpha_rev = 1-alpha\n",
    "    n = data.shape[0]\n",
    "\n",
    "    pows = alpha_rev**(arange(n+1))\n",
    "\n",
    "    scale_arr = 1/pows[:-1]\n",
    "    offset = data[0]*pows[1:]\n",
    "    pw0 = alpha*alpha_rev**(n-1)\n",
    "\n",
    "    mult = data*pw0*scale_arr\n",
    "    cumsums = mult.cumsum()\n",
    "    out = offset + cumsums*scale_arr[::-1]\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "def plot_hist(v, el, n, l='', alpha=0.2, normal=True, df=3):\n",
    "    if normal:\n",
    "        rfun = stats.norm\n",
    "    else:\n",
    "        rfun = stats.chi2\n",
    "    hist(v, label=f'{el}{l}', bins='auto', alpha=alpha, density=True)\n",
    "    if normal or df is None:\n",
    "        fit = rfun.fit(v)\n",
    "    else:\n",
    "        fit = rfun.fit(v, f0=df)\n",
    "    m = v.mean()\n",
    "    s = v.std()\n",
    "    axvline(v.mean(), ls='--', color=f'C{n}', label=f'{el}{l} = {m:.2f}+/-{s:.2f}')\n",
    "    if normal:\n",
    "        x = linspace(m-3*s, m+3*s, 100)\n",
    "    else:\n",
    "        x = linspace(0, m+3*s, 100)\n",
    "    plot(x, rfun.pdf(x, *fit), ls='-', lw=1, color=f'C{n}', alpha=alpha+0.1)\n",
    "    if normal:\n",
    "        xlim(m-3*s, m+3*s)\n",
    "    else:\n",
    "        xlim(0, m+3*s)\n",
    "    return m, s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "def plot_virial_stat(cryst, smpl, T):\n",
    "    elems = cryst.get_chemical_symbols()\n",
    "    vir = array([abs(s[2]*s[3]) for s in smpl])/(un.kB*T)\n",
    "    nat = len(elems)\n",
    "    m, s = plot_hist(vir.mean(axis=(-1,-2)), 'Total', 0, normal=True)\n",
    "    for n, el in enumerate(set(elems)):\n",
    "        elmask = array(elems)==el\n",
    "        plot_hist(vir[:, elmask, :].mean(axis=(-1,-2)), \n",
    "                  el, n+1, normal=False, df=3*sum(elmask))\n",
    "    axvline(T/T, ls=':', color='C5', label=f'{T:.0f} K')\n",
    "    xlim(m - 5*s, m + 7*s)\n",
    "    legend()\n",
    "    title('Virial distribution in the sample')\n",
    "    ylabel('Probability density')\n",
    "    xlabel('Virial/Temperature');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "\n",
    "def plot_acceptance_history(smpl):\n",
    "    figure(figsize=(10,4))\n",
    "    na = array([n for i, n, x, f, e in smpl])\n",
    "    na = cumsum((na[1:]-na[:-1])%2)\n",
    "    plot(100*(na/arange(1,len(na)+1))[1:])\n",
    "    xlabel('Step')\n",
    "    ylabel('Acceptance ratio (approx., %)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "def plot_dofmu_stat(cryst, dofmu, skip=10, window=10):\n",
    "    symm = get_symmetry_dataset(cryst)\n",
    "    dofmap = symm['mapping_to_primitive']\n",
    "    dof = set(dofmap)\n",
    "    dofmul = Counter(symm['mapping_to_primitive'])\n",
    "    elems = dict(zip(symm['std_mapping_to_primitive'],symm['std_types']))\n",
    "    elmap = array(sorted(elems.items())).T\n",
    "    xdof = array(dofmu)\n",
    "    skip = min(skip, len(dofmu)//2)\n",
    "    window = min(window, len(dofmu)//2)\n",
    "\n",
    "    figure(figsize=(10,4))\n",
    "\n",
    "    for i, el in enumerate(set(elmap[1])):\n",
    "        n = len(xdof)\n",
    "        elmask = elmap[1]==el\n",
    "        semilogy()\n",
    "        plot(xdof[:,elmask,:].reshape((-1,3*sum(elmask))),\n",
    "                 '.', color=f'C{i}', ms=1, alpha=0.2)\n",
    "\n",
    "        asx = moving_average(xdof[:,elmask,:].mean((-1,-2)), window)\n",
    "        plot((n-len(asx))//2+arange(len(asx)), asx, '--',\n",
    "                 label=f'{chemical_symbols[el]} (ma, w={window})', color=f'C{i}');\n",
    "\n",
    "        asx = ewma(xdof[:,elmask,:].mean((-1,-2)), window)\n",
    "        plot((n-len(asx))//2+arange(len(asx)), asx, \n",
    "                 label=f'{chemical_symbols[el]} (ewma, w={window})', color=f'C{i}');\n",
    "\n",
    "    axvspan(0, skip, color='k', alpha=0.2)\n",
    "    title('Virial history')\n",
    "    xlabel('Steps')\n",
    "    ylabel('Virial/Temperature V$_{DOF}$/T')\n",
    "    legend();\n",
    "    show();\n",
    "\n",
    "    mi, ma = -1, -1\n",
    "\n",
    "    for i, el in enumerate(set(elmap[1])):\n",
    "        elmask = elmap[1]==el\n",
    "        m, s = plot_hist(xdof[skip:,elmask,:].mean((-2, -1)), \n",
    "                         chemical_symbols[el], i,)\n",
    "                         # normal=False, df=3*sum(elmask))\n",
    "        if mi < 0 or mi > m-3*s:\n",
    "            mi = m-3*s\n",
    "        if ma < 0 or ma < m+4*s:\n",
    "            ma = m+4*s\n",
    "    legend();\n",
    "    xlim(mi, ma)\n",
    "    title('Virial distribution')\n",
    "    ylabel('Probability density')\n",
    "    xlabel('Virial/Temperature')\n",
    "    show() ;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "\n",
    "def plot_xs_stat(cryst, xsl, skip=10, window=10):\n",
    "    symm = get_symmetry_dataset(cryst)\n",
    "    dofmap = symm['mapping_to_primitive']\n",
    "    dof = set(dofmap)\n",
    "    dofmul = Counter(symm['mapping_to_primitive'])\n",
    "    elems = dict(zip(symm['std_mapping_to_primitive'],symm['std_types']))\n",
    "    elmap = array(sorted(elems.items())).T\n",
    "    elmap = cryst.get_atomic_numbers()\n",
    "    xdof = array(xsl)\n",
    "    skip = min(skip, len(xsl)//2)\n",
    "    window = min(window, len(xsl)//2)\n",
    "\n",
    "    plt.figure(figsize=(10,4))\n",
    "\n",
    "    for i, el in enumerate(set(elmap)):\n",
    "        n = len(xdof)\n",
    "        plot(xdof[:,elmap==el,:].mean((-2,-1)),\n",
    "                 '.', color=f'C{i}', ms=2, alpha=0.25)\n",
    "\n",
    "        asx = moving_average(xdof[:,elmap==el,:].mean((-1,-2)), window)\n",
    "        plot((n-len(asx))//2 + arange(len(asx)), asx, '--',\n",
    "                 label=f'{chemical_symbols[el]} (ma, w={window})', color=f'C{i}');\n",
    "\n",
    "        asx = ewma(xdof[:,elmap==el,:].mean((-1,-2)), window)\n",
    "        plot((n-len(asx))//2 + arange(len(asx)), asx, \n",
    "                 label=f'{chemical_symbols[el]} (ewma, w={window})', color=f'C{i}');\n",
    "\n",
    "    axvspan(0, skip, color='k', alpha=0.2)\n",
    "\n",
    "    title('Amplitude correction history')\n",
    "    xlabel('Steps')\n",
    "    ylabel('Virial/Temperature V$_{DOF}$/T')\n",
    "    legend();\n",
    "    show();\n",
    "\n",
    "    mi, ma = -1, -1\n",
    "    for i, el in enumerate(set(elmap)):\n",
    "        m, s = plot_hist(xdof[skip:,elmap==el,:].mean((-2,-1)), \n",
    "                  chemical_symbols[el], i)\n",
    "        if mi < 0 or mi > m-3*s:\n",
    "            mi = m-3*s\n",
    "        if ma < 0 or ma < m+3*s:\n",
    "            ma = m+3*s\n",
    "    legend();\n",
    "    xlim(mi, ma)\n",
    "    title('Amplitude correction distribution')\n",
    "    ylabel('Probability density')\n",
    "    xlabel('Virial/Temperature')\n",
    "    show() ;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}