HMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f67b692e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "### Util for the KMC code with numpy vectorisation :>\n",
    "import numpy as np\n",
    "import sys\n",
    "import math\n",
    "import random\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "# Init File for the KMC Python package\n",
    "Xup = 0\n",
    "Tup = 1\n",
    "Tdn = 2\n",
    "Xdn = 3\n",
    "Ndirs = 4\n",
    "\n",
    "COLDLAT = 0\n",
    "HOTLAT = 1\n",
    "PLUS = 1\n",
    "MINUS = -1\n",
    "\n",
    "FRESH = 0\n",
    "FORGET = 1\n",
    "FRESHHOT = 2\n",
    "MAXFILENAME = 256\n",
    "\n",
    "lattice = None\n",
    "\n",
    "Nx = None\n",
    "Nt = None\n",
    "N = None\n",
    "Volume = None\n",
    "Nf = None\n",
    "G = None\n",
    "Step = None\n",
    "Warmup = None\n",
    "NMdtraj = None\n",
    "Mdstep = None\n",
    "max_hmc_iter = None\n",
    "Nmeas = None\n",
    "Measlen = None\n",
    "CgiterH = None\n",
    "CgiterP = None\n",
    "ResidueH = None\n",
    "ResidueP = None\n",
    "\n",
    "sites_on_node = None\n",
    "\n",
    "readflag = None\n",
    "saveflag = None\n",
    "read_binary_sigma = None # Although defined as char array we take string\n",
    "save_binary_sigma = None # Although defined as char array we take string\n",
    "\n",
    "av_sigma = None #Defined as pointer check later !!!!\n",
    "\n",
    "iseed = None\n",
    "junk_id = None\n",
    "\n",
    "N_POINTERS = 4\n",
    "\n",
    "gen_pt = None\n",
    "neighbor = None # neighbour stores the index of the neighbouring site in place of the pointer\n",
    "\n",
    "#-----------------------\n",
    "IM1 = 2147483563\n",
    "IM2 = 2147483399\n",
    "AM = 1.0 / IM1\n",
    "IMM1 = IM1 - 1\n",
    "IA1 = 40014\n",
    "IA2 = 40692\n",
    "IQ1 = 53668\n",
    "IQ2 = 52774\n",
    "IR1 = 12211\n",
    "IR2 = 3791\n",
    "NTAB = 32\n",
    "NDIV = 1 + IMM1 / NTAB\n",
    "EPS = 1.2e-7\n",
    "RNMX = 1.0 - EPS\n",
    "#-----------------------\n",
    "# These are the static variable defined inside the codebase\n",
    "iset = 0\n",
    "gset = None\n",
    "iseed = -772\n",
    "iseed2 = 123456789\n",
    "iy = 0\n",
    "iv = [0] * NTAB\n",
    "\n",
    "\n",
    "class LATTICE:\n",
    "    def __init__(self,volume_):\n",
    "        self.x = np.array(list(0 for i in range(volume_)))\n",
    "        self.t = np.array(list(0 for i in range(volume_)))\n",
    "        self.index = np.array(list(0 for i in range(volume_)))\n",
    "        self.sign = np.array(list(0 for i in range(volume_)))\n",
    "\n",
    "        self.mom = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.pi = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.sigma = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.phi = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.rho = np.array(list(0.0 for i in range(volume_)))\n",
    "\n",
    "        self.eta = np.array(list(list(0.0 for i in range(4)) for j in range(volume_)))\n",
    "        self.zeta = np.array(list(list(0.0 for i in range(4)) for j in range(volume_)))\n",
    "\n",
    "        self.omega = np.array(list(list(0.0 for i in range(4)) for j in range(volume_)))\n",
    "        self.xi = np.array(list(list(0.0 for i in range(4)) for j in range(volume_)))\n",
    "\n",
    "        self.p = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.r = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.mp = np.array(list(0.0 for i in range(volume_)))\n",
    "        self.mmp = np.array(list(0.0 for i in range(volume_)))\n",
    "\n",
    "    def __getitem__(self,key):\n",
    "        return self.__dict__[key]\n",
    "    def __setitem__(self,key,value):\n",
    "        self.__dict__[key] = value\n",
    "\n",
    "class params:\n",
    "    def __init__(self):\n",
    "        self.stopflag = None\n",
    "        self.iseed = None\n",
    "        self.nx = None\n",
    "        self.nt = None\n",
    "        self.nflv = None\n",
    "        self.g = None\n",
    "        self.step = None\n",
    "        self.eps = None\n",
    "        self.krmax = None\n",
    "        self.krtraj = None\n",
    "        self.cgiter1 = None\n",
    "        self.cgiter2 = None\n",
    "\n",
    "        self.residue1 = None\n",
    "        self.residue2 = None\n",
    "\n",
    "        self.read_sigma = None\n",
    "        self.save_sigma = None\n",
    "\n",
    "        self.read_bin_double = None\n",
    "        self.save_binary_double = None\n",
    "\n",
    "par_buf = params()\n",
    "\n",
    "def setup():\n",
    "    prompt = initial_set()\n",
    "    setup_layout()\n",
    "    make_lattice()\n",
    "    make_nn_gathers()\n",
    "\n",
    "    return prompt\n",
    "\n",
    "\n",
    "def initial_set():\n",
    "    global Nx\n",
    "    global Nt\n",
    "    global N\n",
    "    global Volume\n",
    "\n",
    "    prompt = getprompt()\n",
    "\n",
    "    print(\"GN model with HMC algorithm\\n\")\n",
    "    print(\"Alpha machine, Version 1\\n\")\n",
    "    print(\"type 0 for no prompts, 1 for prompts or 2 for list of prompts\\n\")\n",
    "\n",
    "    prompt = getprompt() #Using the direct value inplace of the pointer in the original code base\n",
    "\n",
    "#     Nx = get_i(prompt, \"nx\")\n",
    "    Nx = 16\n",
    "#     Nt = get_i(prompt, \"nt\")\n",
    "    Nt = 16\n",
    "    if Nx == Nt:\n",
    "        N = Nx\n",
    "    else:\n",
    "        print(\"Nx and Nt must be same, exiting\")\n",
    "        sys.exit()\n",
    "\n",
    "    Volume = Nx*Nt\n",
    "    print(f\"Lattice dimensions = {Nx} {Nt}\\n\")\n",
    "\n",
    "    return prompt\n",
    "\n",
    "def readin(prompt):\n",
    "    global Nf, G, Step, Eps, NMdtraj, KRmax, Measlen, CgiterH, ResidueH, CgiterP, ResidueP, max_hmc_iter, Mdstep\n",
    "#     Nf = get_i(prompt, \"num_flavors\")\n",
    "#     G = get_lf(prompt, \"coupling\")\n",
    "#     Step = get_lf(prompt, \"step_size\")\n",
    "#     Eps = get_lf(prompt, \"noise_parameter\")\n",
    "#     KRtraj = get_i(prompt, \"kr_trajectory_length\")\n",
    "#     KRmax = get_i(prompt, \"kr_refreshing_steps\")\n",
    "#     Measlen = get_i(prompt, \"measurement_interval\")\n",
    "#     CgiterH = get_i(prompt, \"max_cg_iter_for_hamil\")\n",
    "#     ResidueH = get_lf(prompt, \"residue_cg_hamil\")\n",
    "#     CgiterP = get_i(prompt, \"max_cg_iter_for_piup\")\n",
    "#     ResidueP = get_lf(prompt, \"residue_cg_piup\")\n",
    "\n",
    "    Nf = 4\n",
    "    G = 0.54\n",
    "    Step = 0.02\n",
    "    Eps = 1.5\n",
    "    NMdtraj = 300\n",
    "    Mdstep = 50\n",
    "    max_hmc_iter = 1000\n",
    "    KRmax = 2\n",
    "    Measlen = 1\n",
    "    CgiterH = 20000\n",
    "    ResidueH = 0.000001\n",
    "    CgiterP = 10000\n",
    "    ResidueP = 0.0001\n",
    "\n",
    "def make_lattice():\n",
    "\n",
    "    global lattice, Volume, Nx, neighbor\n",
    "\n",
    "    lattice = LATTICE(Volume)\n",
    "    neighbor = np.array(list(list(0.0 for i in range(Volume)) for j in range(4)))\n",
    "\n",
    "    for t_ in range(Nt):\n",
    "        for x_ in range(Nx):\n",
    "            i = site_index(x_,t_)     # Function not defined yet !!!!\n",
    "            lattice.x[i] = x_\n",
    "            lattice.t[i] = t_\n",
    "            lattice.index[i] = x_+Nx*t_\n",
    "            if t_%2 ==0:\n",
    "                lattice.sign[i] = 1\n",
    "            else:\n",
    "                lattice.sign[i] = -1\n",
    "\n",
    "def setup_layout():\n",
    "    global Volume, sites_on_node\n",
    "    sites_on_node = Volume\n",
    "\n",
    "def get_lf(prompt,variable_name_string):\n",
    "\n",
    "    if prompt == 1:\n",
    "        x = float(input(f\"enter {variable_name_string}\"))\n",
    "        return(x)\n",
    "\n",
    "    else:\n",
    "        print(\"Change Prompt type currently not supported\")\n",
    "        system.exit()\n",
    "\n",
    "\n",
    "def get_i(prompt, variable_name_string):\n",
    "\n",
    "    if prompt == 1:\n",
    "        x = int(input(f\"enter {variable_name_string}\"))\n",
    "        return(x)\n",
    "\n",
    "    else:\n",
    "        print(\"Change Prompt type currently not supported\")\n",
    "        system.exit()\n",
    "\n",
    "\n",
    "def getprompt():\n",
    "#     prompt = int(input(\"Enter the prompt type (only 1 is supported for now):\"))\n",
    "    prompt = 1\n",
    "    return prompt\n",
    "\n",
    "def average_sigma():\n",
    "    global Volume, lattice\n",
    "\n",
    "    av_sigma = None\n",
    "    t_sigma = None\n",
    "    t_sigma = 0\n",
    "\n",
    "    t_sigma = np.sum(lattice.sigma[:])\n",
    "    av_sigma = t_sigma/Volume\n",
    "    return av_sigma\n",
    "\n",
    "\n",
    "def average_pbp():\n",
    "\n",
    "    print(\"check this function with sir\")\n",
    "\n",
    "\n",
    "def congrad(src,dest,cgiter,residue,cgflag,flv):\n",
    "    niter = None\n",
    "    i = None\n",
    "    size_src = None\n",
    "    cp = None\n",
    "    d = None\n",
    "    dsize_r = None\n",
    "    dsize_src = None\n",
    "    size_r = None\n",
    "    a = None\n",
    "    b = None\n",
    "    c = None\n",
    "\n",
    "    global Volume\n",
    "    global lattice\n",
    "\n",
    "    dsize_src = 0\n",
    "\n",
    "    dsize_src = np.sum((lattice[src][:, flv])**2)\n",
    "    size_src = math.sqrt(dsize_src)\n",
    "\n",
    "    if cgflag == 0:\n",
    "\n",
    "        lattice[dest][:,flv] = 0\n",
    "        lattice.r[:] = lattice[src][:,flv]\n",
    "        lattice.p[:] = lattice.r[:]\n",
    "\n",
    "        dsize_r = 1\n",
    "        size_r = dsize_r\n",
    "\n",
    "    if cgflag != 0:\n",
    "        matp2d(dest, \"mp\", PLUS, flv)\n",
    "        matd2d(\"mp\", \"mmp\", MINUS)\n",
    "\n",
    "        dsize_r = 0\n",
    "\n",
    "        lattice.r[:] = lattice[src][:,flv] - lattice.mmp[:]\n",
    "        dsize_r = np.sum(lattice.r[:]**2)\n",
    "        lattice.p[:] = lattice.r[:]\n",
    "\n",
    "        size_r = math.sqrt(dsize_r)/size_src\n",
    "\n",
    "    cp = dsize_r\n",
    "\n",
    "    niter = 0\n",
    "\n",
    "    while niter < cgiter and size_r > residue:\n",
    "        c = cp\n",
    "\n",
    "        matd2d(\"p\",\"mp\",PLUS)\n",
    "        matd2d(\"mp\",\"mmp\",MINUS)\n",
    "\n",
    "        d=0\n",
    "        d = np.sum((lattice.p[:]*lattice.mmp[:]))\n",
    "        a = c/d\n",
    "\n",
    "        cp = 0\n",
    "        lattice[dest][:,flv] = lattice[dest][:,flv] + a*lattice.p[:]\n",
    "        lattice.r[:] = lattice.r[:] - a*lattice.mmp[:]\n",
    "        cp = np.sum(lattice.r[:]**2)\n",
    "\n",
    "        b = cp/c\n",
    "        dsize_r = 0\n",
    "\n",
    "        lattice.p[:] = lattice.r[:] + b*lattice.p[:]\n",
    "        dsize_r = np.sum(lattice.r[i]**2)\n",
    "        size_r = math.sqrt(dsize_r)/size_src\n",
    "\n",
    "        niter = niter + 1\n",
    "\n",
    "    if size_r > residue:\n",
    "        # print(\"CG_MD Not Converged\")\n",
    "        system.exit(1)\n",
    "\n",
    "def matd2d(src,dest,isign):\n",
    "    i = None\n",
    "    n = None\n",
    "\n",
    "    global Xup,Xdn,Tup,Tdn, volume, lattice\n",
    "\n",
    "    gather(src, Xup, gen_pt)\n",
    "    lattice[dest][:] = 0.5*lattice.sign[:]*gen_pt[:]\n",
    "\n",
    "    gather(src, Xdn, gen_pt)\n",
    "    lattice[dest][:] = lattice[dest][:] - 0.5*lattice.sign[:] * gen_pt[:]\n",
    "\n",
    "    gather(src, Tup, gen_pt)\n",
    "    lattice[dest][:] = lattice[dest][:] + 0.5 * gen_pt[:]\n",
    "\n",
    "    gather(src, Tdn, gen_pt)\n",
    "    lattice[dest][:] = lattice[dest][:] - 0.5 * gen_pt[:]\n",
    "\n",
    "    lattice[dest][:] = isign*(lattice[dest][:]) + (lattice.phi[:])*(lattice[src][:])\n",
    "\n",
    "\n",
    "def matd2p(src, dest, isign, flavor):\n",
    "    i = None\n",
    "    n = None\n",
    "\n",
    "    global Xup, Xdn, gen_pt, Tdn, Tup, lattice\n",
    "\n",
    "    gather(src, Xup, gen_pt)\n",
    "    lattice[dest][i,flavor] = lattice.sign[:] * 0.5 * gen_pt[:]\n",
    "\n",
    "    gather(src, Xdn, gen_pt)\n",
    "    lattice[dest][:,flavor] = lattice[dest][:,flavor] - 0.5*lattice.sign[:] * gen_pt[:]\n",
    "\n",
    "    gather(src, Tup, gen_pt)\n",
    "    lattice[dest][:,flavor] = lattice[dest][:,flavor] + 0.5 * gen_pt[:]\n",
    "\n",
    "    gather(src, Tdn, gen_pt)\n",
    "    lattice[dest][:,flavor] = lattice[dest][:,flavor] - 0.5 * gen_pt[:]\n",
    "\n",
    "    lattice[dest][:,flavor] = isign*lattice[dest][:,flavor] + lattice.phi[:]*lattice[src][:]\n",
    "\n",
    "\n",
    "def matp2d(src,dest,isign,flavor):\n",
    "\n",
    "    global Xup,gen_pt, Xdn, Tdn, Tup, lattice\n",
    "\n",
    "    gather(src, Xup, gen_pt)\n",
    "    lattice[dest][:] = 0.5*lattice.sign[:] * gen_pt[:,flavor]\n",
    "\n",
    "    gather(src, Xdn, gen_pt)\n",
    "    lattice[dest][:] = lattice[dest][:] - 0.5 * lattice.sign[:] * gen_pt[:,flavor]\n",
    "\n",
    "    gather(src, Tup, gen_pt)\n",
    "    lattice[dest][:] =  lattice[dest][:] + 0.5 * gen_pt[:,flavor]\n",
    "\n",
    "    gather(src, Tdn, gen_pt)\n",
    "    lattice[dest][:] = lattice[dest][:] - 0.5 * gen_pt[:,flavor]\n",
    "\n",
    "    lattice[dest][:] = isign*lattice[dest][:] +lattice.phi[:]*lattice[src][:,flavor]\n",
    "\n",
    "\n",
    "\n",
    "def matp2p(src,dest,isign,flav):\n",
    "\n",
    "    global Xup, Xdn, Tup, Tdn, gen_pt, lattice\n",
    "\n",
    "    gather(src, Xup, gen_pt)\n",
    "    lattice[dest][:,flav] = lattice.sign[:] * 0.5 * gen_pt[:,flav]\n",
    "\n",
    "    gather(src, Xdn, gen_pt)\n",
    "    lattice[dest][:,flav] -= lattice.sign[:] * 0.5 * gen_pt[:,flav]\n",
    "\n",
    "    gather(src, Tup, gen_pt)\n",
    "    lattice[dest][:,flav] +=  0.5 * gen_pt[:,flav]\n",
    "\n",
    "    gather(src, Tdn, gen_pt)\n",
    "    lattice[dest][:,flav] -=  0.5 * gen_pt[:,flav]\n",
    "\n",
    "\n",
    "    lattice[dest][:,flav] = isign*(lattice[dest][:,flav]) + (lattice.phi[:])*(lattice[src][:,flav])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def nn_coords(x,t,dir):\n",
    "    global neighbor\n",
    "    xp = x\n",
    "    tp = t\n",
    "\n",
    "    global Nx, Nt\n",
    "\n",
    "    if dir == 0:\n",
    "        xp = (x+1)%Nx\n",
    "    if dir == 3:\n",
    "        xp = (x+Nx-1)%Nx\n",
    "    if dir == 1:\n",
    "        tp = (t+1)%Nt\n",
    "    if dir == 2:\n",
    "        tp = (t+Nt-1)%Nt\n",
    "    return xp, tp\n",
    "\n",
    "def make_nn_gathers():\n",
    "\n",
    "    #This will run only once hence not vectorised\n",
    "\n",
    "    global neighbor, lattice, Volume\n",
    "    i = None\n",
    "    j = None\n",
    "    dir = None\n",
    "\n",
    "    xpt = None\n",
    "    tpt = None\n",
    "    mul = None\n",
    "\n",
    "    for dir in range(0,4):\n",
    "        for i in range(Volume):\n",
    "            xpt, tpt = nn_coords(lattice.x[i],lattice.t[i], dir)\n",
    "            j = site_index(xpt,tpt)\n",
    "            neighbor[dir,i] = j\n",
    "    neighbor = neighbor.astype(np.uint16)\n",
    "\n",
    "\n",
    "def gather(field, index, dest):\n",
    "\n",
    "    global neighbor, lattice\n",
    "    global gen_pt\n",
    "    gen_pt = lattice[field][neighbor[index,:]].copy()\n",
    "\n",
    "\n",
    "def hmc():\n",
    "    global lattice, Nf, MINUS, max_hmc_iter, Mdstep, Step\n",
    "\n",
    "    i, j, nf, ctr = 0, 0, 0, 0\n",
    "    hold, hnew, deltah = 0.0, 0.0, 0.0\n",
    "    xx = 1.0\n",
    "    z = 0\n",
    "    max_count = 100\n",
    "    count = 0\n",
    "\n",
    "    while (ctr < max_hmc_iter and xx>z):\n",
    "\n",
    "        lattice.phi[:] = lattice.sigma[:]\n",
    "        lattice.pi[:] = np.random.normal(loc = 0, scale = 1, size = lattice.pi.shape)\n",
    "        lattice.zeta = np.random.normal(loc = 0, scale = (1/math.sqrt(2)), size = lattice.zeta.shape)\n",
    "\n",
    "        try:\n",
    "            hold = hamil(1,1)\n",
    "            piup(Step/2)\n",
    "\n",
    "            for j in range(Mdstep-1):\n",
    "                lattice.phi[:] = lattice.phi[:] + lattice.pi[:]*Step\n",
    "                piup(Step)\n",
    "\n",
    "            piup(Step/2)\n",
    "\n",
    "            hnew = hamil(1,1)\n",
    "\n",
    "            xx = random.random()\n",
    "            deltah = hnew-hold\n",
    "            z = math.exp(-1*deltah)\n",
    "            z = min(1,z)\n",
    "        except:\n",
    "            z = 0\n",
    "        # print(f\"Old Hamiltonian = {hold}\")\n",
    "        # print(f\"New Hamiltonian = {hnew}\")\n",
    "        # print(f\"DeltaH = {deltah}\")\n",
    "        ctr = ctr + 1\n",
    "\n",
    "    if xx<z:\n",
    "        # print(\"Accepted\")\n",
    "        lattice.sigma[:] = lattice.phi[:]\n",
    "        \n",
    "    # else:\n",
    "        # print(\"Rejected\")\n",
    "        \n",
    "\n",
    "    return ctr\n",
    "\n",
    "def layout():\n",
    "    global no_odd_sites, no_even_sites, volume\n",
    "\n",
    "    no_even_sites = volume/2\n",
    "    no_odd_sites = volume/2\n",
    "\n",
    "\n",
    "def site_index(x,t):\n",
    "    global Nx, Nt\n",
    "    i = None\n",
    "    xr = None\n",
    "    tr = None\n",
    "\n",
    "    xr = x%Nx\n",
    "    tr = t%Nt\n",
    "\n",
    "    i = xr + Nx*tr\n",
    "\n",
    "    return i\n",
    "\n",
    "\n",
    "def hamil(hflag, flag):\n",
    "    global lattice, Volume, G, Nf, CgiterH, ResidueH\n",
    "    i,n = 0,0\n",
    "    h = 0\n",
    "\n",
    "    hml = np.sum(((0.5/(G*G)) * lattice.phi* lattice.phi)) + np.sum((0.5 * lattice.pi * lattice.pi))\n",
    "\n",
    "    if hflag != 0:\n",
    "        for n in range(Nf):\n",
    "            congrad(\"zeta\",\"eta\",CgiterH, ResidueH, flag, n)\n",
    "            hml = hml + np.sum(lattice.zeta[:,n] * lattice.eta[:,n])\n",
    "    else:\n",
    "        for n in range(Nf):\n",
    "            hml = hml + np.sum(lattice.zeta[:,n] * lattice.eta[:,n])\n",
    "\n",
    "    return hml\n",
    "\n",
    "\n",
    "def zerolat():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    print(\"ZEROLAT: All zero initial config. of `sigma' field.\")\n",
    "    lattice.sigma[:] = 0\n",
    "\n",
    "\n",
    "def coldlat():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    print(\" COLDLAT: Cold initial config. of `sigma' field\")\n",
    "\n",
    "    lattice.sigma[:] = 1\n",
    "\n",
    "\n",
    "def coldlat2():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    print(\" COLDLAT.4: Cold initial config. of `sigma' field\")\n",
    "\n",
    "    lattice.sigma[:] = 0.4\n",
    "\n",
    "\n",
    "def hotlat():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    print(\" HOTLAT: Hot initial config. of `sigma' field.\")\n",
    "\n",
    "    for i in range(volume):\n",
    "        lattice.sigma[i] = 2 * ran2() - 1\n",
    "\n",
    "\n",
    "def filelat():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    print(\" Configuration to be read from the file sigma.in\")\n",
    "    print(\"currently not supported was feeling tired\")\n",
    "    # for i in range(volume):\n",
    "        # lattice[i].sigma = 0.4\n",
    "\n",
    "def funnylat():\n",
    "\n",
    "    global lattice, volume\n",
    "\n",
    "    for i in range(volume):\n",
    "        lattice[i].sigma = i\n",
    "\n",
    "def piup(t):\n",
    "    global lattice, Volume, G, Nf, CgiterP,ResidueP, PLUS\n",
    "    # print(\"Piup Calculation. \\n\")\n",
    "\n",
    "    lattice.pi[:] = lattice.pi[:] - ((1/(G*G)) * lattice.phi[:] * t)\n",
    "\n",
    "    for n in range(Nf):\n",
    "        congrad(\"zeta\",\"xi\", CgiterP, ResidueP, 1, n)\n",
    "        matp2p(\"xi\",\"omega\",PLUS,n)\n",
    "        lattice.pi[:] = lattice.pi[:] + (2 * lattice.omega[:, n] * lattice.xi[:, n] * t)\n",
    "\n",
    "\n",
    "\n",
    "def main():\n",
    "    global lattice,NMdtraj, Measlen\n",
    "    console = sys.stdout\n",
    "    print(\" ===================== HMC ALGORITHM =====================\")\n",
    "#     print(\" =============Shifting from console to out.dat file============\")\n",
    "#     f = open(\"./OUTPUT.dat\",\"w\")\n",
    "#     sys.stdout = f\n",
    "#     print(\" ===================== HMC ALGORITHM =====================\")\n",
    "    AVERAGE_SIGMA = []\n",
    "    prompt = setup()\n",
    "    readin(prompt)\n",
    "\n",
    "    coldlat()\n",
    "    l=0\n",
    "    numtrial=0\n",
    "    numaccp=0\n",
    "\n",
    "    for todo in tqdm(range(0,NMdtraj)):\n",
    "\n",
    "        numtrial += hmc()\n",
    "        numaccp += 1\n",
    "\n",
    "        # if todo%Measlen == 0:\n",
    "            # print(f\"trajectory = {todo} , av_sigma = {average_sigma()}\")\n",
    "\n",
    "        AVERAGE_SIGMA.append(average_sigma())\n",
    "    acceptance_rate = numaccp/numtrial\n",
    "    print(f\"Acceptance rate = {acceptance_rate}\")\n",
    "    print(\" =======================Success============================\")\n",
    "#     sys.stdout = console\n",
    "#     f.close()\n",
    "#     print(\" =============Shifting from out.dat to console ============\")\n",
    "    print(f\"Acceptance rate = {acceptance_rate}\")\n",
    "    return AVERAGE_SIGMA, acceptance_rate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f90c0d19",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GN model with HMC algorithm\n",
      "\n",
      "Alpha machine, Version 1\n",
      "\n",
      "type 0 for no prompts, 1 for prompts or 2 for list of prompts\n",
      "\n",
      "Lattice dimensions = 16 16\n",
      "\n",
      " COLDLAT: Cold initial config. of `sigma' field\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87581425c6ce4f9aaa22242ff4fcdd06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/10 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "2\n",
      "9\n",
      "5\n",
      "3\n",
      "1\n",
      "Acceptance rate = 0.4\n",
      " =======================Success============================\n",
      "Acceptance rate = 0.4\n"
     ]
    },
    {
     "ename": "SyntaxError",
     "evalue": "'return' outside function (4007444480.py, line 27)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[12], line 27\u001b[0;36m\u001b[0m\n\u001b[0;31m    return AVERAGE_SIGMA, acceptance_rate\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m 'return' outside function\n"
     ]
    }
   ],
   "source": [
    "AVERAGE_SIGMA = []\n",
    "prompt = setup()\n",
    "readin(prompt)\n",
    "\n",
    "coldlat()\n",
    "l=0\n",
    "numtrial=0\n",
    "numaccp=0\n",
    "NMdtraj = 20\n",
    "for todo in tqdm(range(0,NMdtraj)):\n",
    "    trial = hmc()\n",
    "    numtrial += trial\n",
    "    numaccp += 1\n",
    "    print(trial)\n",
    "\n",
    "    # if todo%Measlen == 0:\n",
    "        # print(f\"trajectory = {todo} , av_sigma = {average_sigma()}\")\n",
    "\n",
    "    AVERAGE_SIGMA.append(average_sigma())\n",
    "acceptance_rate = numaccp/numtrial\n",
    "print(f\"Acceptance rate = {acceptance_rate}\")\n",
    "print(\" =======================Success============================\")\n",
    "#     sys.stdout = console\n",
    "#     f.close()\n",
    "#     print(\" =============Shifting from out.dat to console ============\")\n",
    "print(f\"Acceptance rate = {acceptance_rate}\")\n",
    "return AVERAGE_SIGMA, acceptance_rate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40a7b100",
   "metadata": {},
   "outputs": [],
   "source": [
    "hmc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8d253d4f-411e-446a-8e8a-938a5ebb3dbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc6b70e9c90>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABeyklEQVR4nO3deXiU5b0+8Hv2yTYTsq8kYQ8EQgjKjrZqLCrVahVrK2qxpxywiLT9tZSeLp7TcrrIwQ2sVeRYl3Ja0GpFa7rIqgIhgbAl7AnZF5LJOuv7+2PmnUzIQmYyM+8s9+e65rrMZJYnJszc8zzf5/vIBEEQQERERCQRudQDICIiovDGMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUlKKfUARsJms6G2thYxMTGQyWRSD4eIiIhGQBAEdHR0IC0tDXL50PMfQRFGamtrkZmZKfUwiIiIyAPV1dXIyMgY8vtBEUZiYmIA2H8YnU4n8WiIiIhoJAwGAzIzM53v40MJijAiLs3odDqGESIioiBzvRILFrASERGRpBhGiIiISFIMI0RERCQphhEiIiKSFMMIERERScrtMLJ3714sXboUaWlpkMlkePfdd4e9/a5du3DbbbchMTEROp0O8+bNw9/+9jdPx0tEREQhxu0w0tXVhfz8fLzwwgsjuv3evXtx2223Yffu3SgpKcEXvvAFLF26FKWlpW4PloiIiEKPTBAEweM7y2R45513cM8997h1v2nTpmHZsmX4yU9+MqLbGwwG6PV6tLe3s88IERFRkBjp+7ffm57ZbDZ0dHQgLi5uyNsYjUYYjUbn1waDwR9DIyIiIgn4vYD1mWeeQVdXFx544IEhb7Nx40bo9XrnhefSEBERhS6/hpG3334bP/vZz7Bjxw4kJSUNebv169ejvb3deamurvbjKImIiMif/LZMs2PHDqxYsQJ/+tOfcOuttw57W41GA41G46eRERERkZT8MjPy9ttv49FHH8Vbb72FO++80x9PSSQZq03A/x68hPNNnVIPhYgoKLgdRjo7O1FWVoaysjIAwMWLF1FWVoaqqioA9iWW5cuXO2//9ttvY/ny5XjmmWcwd+5c1NfXo76+Hu3t7d75CYgCzAfldfjpeyfx43dOSD0UIqKg4HYYOXLkCAoKClBQUAAAWLduHQoKCpzbdOvq6pzBBAB+97vfwWKxYPXq1UhNTXVennzySS/9CESB5USNPWiXVbfBavN45zwRUdhwu2bk5ptvxnCtSbZv397v608++cTdpyAKahX1HQCAHrMVF5o6MTE5RuIREREFNp5NQ+RllQ0dzv8+foXLkURE18MwQuRF7d1m1LX3Or8ur2EYISK6HoYRIi+qbOzo9zXDCBHR9TGMEHmRWC8yNi4SAHCq1gCL1SblkIiIAh7DCJEXifUit09LRpRagR6zFeebuiQeFRFRYGMYIfKiM46ZkdxUHaal6QFwqYaI6HoYRoi8RBAE58zIpOQYTM9whJErbRKOiogo8DGMEHlJU4cRbd1myGXAhKRoTE/nzAgR0UgwjBB5SYVjViQ7PgpalcI5M3KqjkWsRETDYRgh8hJxJ83kFHvH1Zz4KERrlOg123COh+YREQ2JYYTIS8QwMsnR/l0ul2Famg4AO7ESEQ2HYYTIS8TiVXFmBICzbuQE60aIiIbEMELkBTabgMoG+1LMJJeD8cS6Ec6MEBENjWGEyAuqr3ajx2yFWilHdnyk83pxZuQ0i1iJiIbEMELkBWK9yITEaCgVff+ssuOjEKNRwmix4Wwji1iJiAbDMELkBYPViwCOItZ0exFrOZdqiIgGxTBC5AUVjnqRa8MIAMzIiAXA5mdERENhGCHygop6AwBgcvLAMJLnqBs5zjBCRDQohhGiUTJZbLjgOJl30mAzIy5FrGYWsRIRDcAwQjRKF5u7YLEJiNEokabXDvh+VnwkYrRKmCw2Z20JERH1YRghGqUzjiWaSSkxkMlkA74vk8nY/IyIaBgMI0SjJM52TBqkXkQkhhE2PyMiGohhhGiUKuodO2mSo4e8jdiJlTMjREQDMYwQjZJzZmSQ4lVRXyfWDpgsLGIlInLFMEI0Cl1GC6pauwEMvq1XNDYuEjqtEiYri1iJiK7FMEI0CmKL94RoDeKjNUPeTiaTOZdq2PyMiKg/hhGiUaisF9vAD10vIhKbnzGMEBH1xzBCNAoVI9hJI5qRHguAZ9QQEV2LYYRoFMTTeqcMU7wqEotYz9QbYLRYfTouIqJgwjBCNAruzIxkxkVAH6GC2Sqg0rEdmIiIGEaIPNbaZUJThxEAMHEEYcS1EyvrRoiI+jCMEHlI3KKbGReBaI1yRPfhjhoiooHcDiN79+7F0qVLkZaWBplMhnfffXfY29fV1eGhhx7C5MmTIZfLsXbtWg+HShRYxHqR4fqLXKtvZqTNF0MiIgpKboeRrq4u5Ofn44UXXhjR7Y1GIxITE7Fhwwbk5+e7PUCiQOVOvYhIDCMV9R0sYiUichjZ3LKLJUuWYMmSJSO+fXZ2Np599lkAwLZt29x9OqKA1ddjZORhJGNMBGIjVWjrNqOivgMzMmJ9NDoiouARkDUjRqMRBoOh34UokAiC4JwZcSeMsIiViGiggAwjGzduhF6vd14yMzOlHhJRP3XtvejotUApl2FcwvW7r7pyhhE2PyMiAhCgYWT9+vVob293Xqqrq6UeElE/4qxITkIU1Er3/hnN4I4aIqJ+3K4Z8QeNRgONZuhDx4ikJtaLTHJjiUaU51LE2mu2QqtSeHVsRETBJiBnRogCnTgzMsWNnTSi9NgIxEWpYbEJzu3BREThzO0w0tnZibKyMpSVlQEALl68iLKyMlRVVQGwL7EsX768333E23d2dqKpqQllZWU4derU6EdPJJGKUcyMyGQy5+zIcS7VEBG5v0xz5MgRfOELX3B+vW7dOgDAI488gu3bt6Ours4ZTEQFBQXO/y4pKcFbb72FrKwsXLp0ycNhE0nHahNwttF+tow7Dc9czUjXY29lE06wiJWIyP0wcvPNN0MQhCG/v3379gHXDXd7omBzuaULJosNWpUcmXGRHj0GZ0aIiPqwZoTITc4lmuQYKOQyjx5D3FFztsFexEpEFM4YRojc5Ekb+Gul6rWIdxSxnq5jUz8iCm8MI0RuEk/r9bReBHB0YnXMjpzgUg0RhTmGESI3VXhwJs1gxE6sx1nESkRhjmGEyA29ZisutXQD8F4YYSdWIgp3DCNEbjjf1AmrTYA+QoWkmNF1CRaXac42drKIlYjCGsMIkRtc60VkMs920ohSdFokRKthtQk4xSJWIgpjDCNEbqiodzQ7G+USDeAoYk1nESsREcMIkRsq6u0zGJ60gR8Mi1iJiBhGiNxS2TC6NvDXmp4RC4AzI0QU3hhGiEaoo9eMmrYeAF4MI46ZkcqGDvSYWMRKROGJYYRohMRZkRSdFvpIlVceM1mnQWKMBjYBLGIlorDFMEI0Qs4zabxULwL0L2Itv9LmtcclIgomDCNEI9S3rTfaq4/b1/yMMyNEFJ4YRohGyPW0Xm/qCyNtXn1cIqJgwTBCNELizMiUFJ1XH1fsxHqusRPdJotXH5uIKBgwjBCNQFOHES1dJshkwIQk7y7TJOu0SBKLWGu5VENE4YdhhGgExFmRrLhIRKgVXn/8GRk8NM+XWjqN6Og1Sz0MIhoCwwjRCIj1It5oAz+YPOeOGoYRb2vqMOILv/0E97/0KWw2QerhENEgGEaIRsAZRrxcvCrizIjv7KlsgqHXgjP1HfjsYovUwyGiQTCMEI1ARYP3e4y4EmdGzjV1osvIIlZv2ne2yfnfu47WSDgSIhoKwwjRddhsAs42+HZmJClGixSdFgI7sXqVzSbgwLlm59cfltex7T5RAGIYIbqOmrYedJmsUCvkyE6I8tnz5PEEX687U9+B5k4TIlQKZMZFoMtkxcen6qUeFhFdg2GE6DrEepFxiVFQKXz3T0asG+EJvt4jLtHMHReHewsyAAA7uVRDFHAYRoiuQ6wX8dVOGlFfJ1aGEW/Z71iiWTgxEffOSrdfd7YJjYZeKYdFRNdgGCG6DrHHiLfbwF9LXKY539SJThaxjlqv2YrPL7YCABZPTEBWfBQKs8bAJgB/KauVeHRE5IphhOg6xGWaKT6eGUmM0SBV7yhiZSfWUTt8qRUmiw3JOo2za644O7KrlEs1RIGEYYRoGGarDeebOgH4fmYEcC1ibfP5c4W6/WftSzSLJiZCJpMBAO6anga1Qo7TdQYGPqIAwjBCNIxLzV0wWwVEqRVIj43w+fPNSGcRq7fsdYaRBOd1+kgVbslNAgC8U3pFknER0UAMI0TDEItXJybHQC6X+fz58hw7ao4zjIxKU4cRpx39WhZMSOj3vXtn2XfVvFtWC4vV5vexEdFADCNEw6j0U72ISNxRc7G5iwe7jYLY6Gxqqg4J0Zp+37tpUiLGRKrQ1GHEgfNsD08UCNwOI3v37sXSpUuRlpYGmUyGd99997r32bNnDwoLC6HVajFu3Di89NJLnoyVyO/O1PtnJ40oIVqDNEcR60nWNHhs3yBLNCK1Uo4v56cBAN45yqUaokDgdhjp6upCfn4+XnjhhRHd/uLFi7jjjjuwaNEilJaW4kc/+hHWrFmDnTt3uj1YIn+r9FOPEVfT2fxsVARBwP5z9mZniyYmDnobcanmo5P13EZNFACU7t5hyZIlWLJkyYhv/9JLL2Hs2LHYvHkzACA3NxdHjhzBb3/7W9x3333uPj2R3/SYrLjc2g3AfzMjgH2p5m8nG9gW3kNnGzvRYDBCo5RjdvaYQW8zI0OPcYlRuNDUhQ/L63D/7Ew/j5KIXPm8ZuTTTz9FUVFRv+tuv/12HDlyBGYz18QpcJ1r7IQgAPFRaiTGaK5/By+ZnhELgDMjnhKXaG7MiYNWpRj0NjKZDPc5ZkfeYc8RIsn5PIzU19cjOTm533XJycmwWCxobm4e9D5GoxEGg6HfhcjfztTb/+78OSsC9BWxXmjugoFFrG4Tz6MZrF7E1d0z7XUjn15oQW1bj8/HRURD88tuGrHhkEgQhEGvF23cuBF6vd55yczkFCr5nxT1IgAQF6V29jQ5WcMg7g6jxYrPL9hbwC+cMHi9iChjTCTmjouDIADvlnF2hEhKPg8jKSkpqK/vf2R3Y2MjlEol4uPjB73P+vXr0d7e7rxUV1f7ephEA1Q02Duv+juMAK6H5rX5/bmD2dHLbegxW5EQrRnRdmzxJN9dR2ucH5KIyP98HkbmzZuH4uLiftd9/PHHmD17NlQq1aD30Wg00Ol0/S5E/lYh0TIN0LejppwzI24Rl2gWTogfUZO6JdNToFHKca6xk6clE0nI7TDS2dmJsrIylJWVAbBv3S0rK0NVVRUA+6zG8uXLnbdfuXIlLl++jHXr1uH06dPYtm0bXn31VXzve9/zzk9A5ANt3SY0GIwAgEnJ0X5/fufMCM+occt+R7OzhUNs6b1WjFaF26elALDPjhCRNNwOI0eOHEFBQQEKCgoAAOvWrUNBQQF+8pOfAADq6uqcwQQAcnJysHv3bnzyySeYOXMm/vM//xPPPfcct/VSQKt0LNGkx0YgRjv4DJ4viWHkUks3i1hH6GqXyTm7cb3iVVdfcZzk+/6xWpjZHp5IEm73Gbn55puHXVvdvn37gOtuuukmHD161N2nIpJMhUTFq6IxUWpkjInAlas9OFHTjvnjR/7mGq4OnG+GINhnspJ12hHfb9GEBCREa9DcacTeyibckpt8/TsRkVfxbBqiQUhZLyKaIdaNsPnZiOx39Be53i6aaykVcuc2Xy7VEEmDYYRoEJX14k4a/9eLiPKcO2oYRq5HEIS+82gmuT+LdK9jqab4dAPae7gsRuRvDCMhYP/ZZiz81T9x8PzgTeTIPYIgOJdpJJ0ZSY8FwDAyEhebu1DT1gO1Qo45OXFu339qqg5TUmJgstiwu7zOByMkouEwjISAtw9X4crVHrx+8LLUQwkJjR1GtPeYoZDLMD5RypkR+5b2yy3daO/mp/XhiLtoCrPGIFLtdikcZDIZvlJgnx3ZxZN8ifyOYSQEHHds//z8YgtsNjZuGq0z9fZZkez4yCHPNvGH2Eg1xsZFAgBO1HJ2ZDh7K8UtvZ4X+t5TkA65DDh86SqqWrq9NTQiGgGGkSB3tcuE6lb7uRpXu83O5QXyXGW9tDtpXE1n3ch1ma02fHahBYB7W3qvlazTYsEE+/15eB6RfzGMBLlr36TEF2XyXCDUi4icRazcUTOksuo2dBotGBOpwrQ0/ageSyxk3VV6he3hifyIYSTIiWFE7Hz96XmGkdESD8gbydkmvubc3suZkSGJu2jmT0iAYgQt4Idz+7QURKoVuNzSjaNVbV4YHRGNBMNIkDtW3QYAWDI9FQDw+cVW1o2MgtUmOMNIQMyMOD7pV7V2o63bJPFoAtN+x3k0i0exRCOKVCvxpTyxPTwLWYn8hWEkyImfmL9+41hEqRVo7zE7CzDJfdWt3eg126BRypEVHyX1cKCPVCEr3lHEykPzBmjvMaPMEchHeh7N9dw3y36S71+P18FosXrlMYloeAwjQayxoxd17b2QyYD8zFjMzrb3V/iUdSMeE+tFJiZHj3rK31vEupHjNW3SDiQAfXq+BTYBGJcQhfTYCK885txx8UjRadHeY8a/zjR65TGJaHgMI0HshGNWZEJiNKI0SswbHw+ARayjIe6kCYQlGtEMRxg5wbqRAfY5lmhGs4vmWgq5DPc4eo7sZHt4Ir9gGAlix6rtb07THUWOc8fZw8jnF1pgZd2IR86IB+QFUBgRt/ce546aAcRmZ95aohGJu2o+qWhEaxdrdYh8jWEkiIn1IvkZsQCAvDQdojVKGHotOF3H+gJPOGdGAmAnjWiaI4xcudqDq3xjdKpq6cbllm4o5TLMHed+C/jhTEqOQV66DmargL8er/XqYxPRQAwjQUoQBOcnZXFmRKmQ44bsMQC4VOMJo8WKi81dAAJjW69IH6FCtqOIlVt8++w7Z1+iKRgbixityuuPf2+BvZCVJ/kS+R7DSJCqN/SiudMIhVyGqak65/XiUg3DiPsuNHXBYhMQo1UiRaeVejj9THfMfjGM9Nnv6C+ycIJ3l2hEX56ZBoVchrLqNpxv6vTJcxCRHcNIkBLrRSYlx/Q7P0UsYv38YivrRtxU6VIvIpMFxk4a0XTHoXnsxGpntQk44KgXWTTJe8WrrhKiNbhpkj3ovMPZESKfYhgJUuWObZ75Gf3bX09N1SFGo0RHrwWnalk34o6KAKwXEU1PjwXAmRHR8SttMPRaEKNVOncb+YJ4ku87pTVsJkjkQwwjQeraehGRUiHHDTn2Yj4u1bgnkNrAX2uaY2akpq2HuzvQt0SzYHwClArfvYzdNjUZMRolatp6cOhSq8+ehyjcMYwEIUEQnJ+QZzg+Mbuax7oRj5wJwB4jIp1WhXEJ9o6wnB3pO49moRf7iwxGq1Lgzhn2oxa4VEPkOwwjQai6tQdt3WaoFfJBj7kXi1gPXWyFxWrz9/CCUqfRgitXewAEZhgB+jqxhnvzs06jBUerrgLwbrOzoYhLNR+U16HXzPbwRL7AMBKExLbguakxUCsH/gqnpukQo1Wiw2jBKfYbGZGzjiWaxBgN4qLUEo9mcOIJvsevtEk7EIl9dr4FFpuAsXGRfjk/6IbsOGSMiUCn0YKPTzX4/PmIwhHDSBAqH6JeRKSQyzDHUTfy6Xku1YxEINeLiPpmRsI7YPZ1XfX9rAgAyOWyvkJWnuRL5BMMI0HomOOT8WD1IiL2G3FPINeLiKal6SCT2YtYWzqNUg9HMuJ5NIv9FEaAvqWavWeb0dQRvv/viXyFYSTI2GyC85PxjMyhtzSKYeTwpausGxmBygA8k+ZaMVoVcsK8iLW2rQfnm7oglwHzxvsvjIxLjEbB2FhYbQLeO8b28ETexjASZC62dKHTaIFWJceExOghb5ebqoNOq0Sn0YIT7DdyXRX19g6bgxUEBxKxp0a4Nj8Tt/TOyIiFPsL7LeCHc69jdmQXl2qIvI5hJMiIb0LT0vTD9ldQyGWYw6WaEWnpNKLZsewxMXnogBcIxLqRcJ0Z2SvBEo3orhlpUClkOFlrcDbIIyLvYBgJMmK9yPQRdJ0Ul2pYxDq8CscSzdi4SESqlRKPZnjTwziM2GwCDjr+lhdO9M15NMMZE6XGFyYnAQB2lXJ2hMibGEaCjDgzkj9MvYhIPFb9yKVWmFk3MqTKICheFU1L10MmA+rae8OukPJUnQGtXSZEqRUoGBsryRjunWU/yffd0hqe/UTkRQwjQcRiteGko/5j+jA7aUS5KTroI1ToMlnDvlHWcCoa7PUigbytVxStUTo7sYbb71Rcopk3Ph4qH7aAH84XpiRCH6FCg8HIGUciL2IYCSLnm7rQY7YiSq1wviENR+7ab4R1I0MSd9IE4gF5g5mREQsg/JZqxOLVhRP8Xy8i0igVWJpvbw/PQlYi72EYCSJivUheuh5y+ciOuJ83Xixi5SFfgxEEwblME8jbel2JRazHw2hHTY/JiiOXHC3gJ/m/XsSVuFTz4Yl6dBktko6FKFR4FEa2bNmCnJwcaLVaFBYWYt++fcPe/sUXX0Rubi4iIiIwefJkvP766x4NNtz11YvEjvg+YhEr60YGV9veiw6jBUq5zNnDI9CJbeHDaZnm84stMFltSNNrRzQr6EsFmbHISYhCj9mKv52sl3QsRKHC7TCyY8cOrF27Fhs2bEBpaSkWLVqEJUuWoKqqatDbb926FevXr8fPfvYznDx5Ej//+c+xevVqvP/++6MefLg57njzGclOGtHk5BiMiVSh22QNq0/SIyXOioxPjB70nJ9ANDXV3om13tCLxo5eqYfjF+ISzaKJiZDJRjYr6CsyWV97+F08yZfIK9x+9d20aRNWrFiBxx9/HLm5udi8eTMyMzOxdevWQW//hz/8Ad/+9rexbNkyjBs3Dg8++CBWrFiBX/3qV6MefDgxWWw47Tj0bsYQZ9IMxl43wn4jQ3G2gQ+SehEAiNIonQ3vwmV2ZN9Z/55Hcz1iGDlwvhn17eERCIl8ya0wYjKZUFJSgqKion7XFxUV4eDBg4Pex2g0QqvV9rsuIiIChw4dgtlsdnO44auyoQMmiw36CBXGxkW6dV9xiy/DyEB9beADu9nZtZz9Rq6EfnfdRkMvKho6IJMBCyQsXnWVGReJG7PjIAjAu2WcHSEaLbfCSHNzM6xWK5KTk/tdn5ycjPr6wddOb7/9drzyyisoKSmBIAg4cuQItm3bBrPZjObm5kHvYzQaYTAY+l3CnbjEMiND7/Y0tXiGx5FLV2GysG7EVUUQ9RhxJZ7YXF7TJu1A/EA8pTcvTY+4KLXEo+lz76y+9vCCwJ4jRKPh0SL5tW+GgiAM+Qb5H//xH1iyZAnmzp0LlUqFu+++G48++igAQKFQDHqfjRs3Qq/XOy+ZmZmeDDOkiG867tSLiCYmRSMuSo0esxXHHTtyyN635VyT2GNEJ/Fo3BNOnVgDbYlGtGR6KtRKOSobOp39f4jIM26FkYSEBCgUigGzII2NjQNmS0QRERHYtm0buru7cenSJVRVVSE7OxsxMTFISBj8xWX9+vVob293Xqqrq90ZZkhynRlxl2u/ES7V9LnU0g2TxYYIlQIZYyKkHo5bpqbpIJcBDQYjGg2hW7MgCIIzjCwKsDCij1Dhtqn21z0WshKNjlthRK1Wo7CwEMXFxf2uLy4uxvz584e9r0qlQkZGBhQKBf74xz/irrvuglw++NNrNBrodLp+l3DWa7Y6lxOmOxpeuYv9RgZyNjtLjh5x35ZAEalWYkKSvc4llGdHztR3oLnTiAiVAoVZY6QezgDiSb7vHauBhVvniTzm9jLNunXr8Morr2Dbtm04ffo0nnrqKVRVVWHlypUA7LMay5cvd96+srISb7zxBs6ePYtDhw7hwQcfxIkTJ/DLX/7Sez9FiDtdZ4DFJiAhWo00vfb6dxiEs9/I5VYYLVZvDi9oBWu9iEg8EiCUt2yLW3pvzImDRjn4sq6UFk9KRHyUGs2dJucMDhG5z+0wsmzZMmzevBlPP/00Zs6cib1792L37t3IysoCANTV1fXrOWK1WvHMM88gPz8ft912G3p7e3Hw4EFkZ2d77YcIdeUu/UU87bEwMSka8VFq9JptIf3m5Q7nTpog2tbranq6fcYwlLf37jsXmEs0IpVCjqX5aQCAnWwPT+Qxj85LX7VqFVatWjXo97Zv397v69zcXJSWlnryNORwrNoRRjxcogHsRcdzx8Xjg/I6fHq+BTdkx3lpdMFLnBkJ2jDi+Hs4HqJhpNdsxeeOGqdFE6VtAT+c+2ZlYPvBSyg+1QBDrxk6rUrqIREFneBoORnmxJ00+R4Ur7piv5E+vWYrLrV0AQieM2muNTXVXsTa1GFEQwgWsZZcvgqjxYakGA0mBXAfmLx0HSYmRcNoseHD8jqph0MUlBhGAlyX0YJzjfbtp55s63UlFrHaX+TDu27kXGMnbAIwJlKFxBiN1MPxSIRa4ax3CcWlt71nmwDYt/RK3QJ+ODKZDF+ZxfbwRKPBMBLgTtUZYBOAFJ0WSTrPildF4xOjkRCtgdFiQ1lVm3cGGKT6dtLEBPQb3fXkhXC/kf0BuqV3MPfMTIdMBnx+sRXVrd1SD4co6DCMBLhj1W0A+jpujoa9bkRcqgnvLb7BXi8iEvvOlIdYM7uWTqOzkVigtIAfTlpsBOY5dqz9he3hidzGMBLgxE+8o60XEYlbfD+9EN7bECsagntbr6hvZsQQUi3JxRbwU1JikBQzuhlBf7l3VgYA+1JNKP0uiPyBYSTAlV8Z/U4aV2IYOVrVhl5z+NaNVDpmRqYE+czI1FQdFHIZmjuNqA+hIlZxiWbxpMDdRXOtL+WlIEKlwIXmLpQ5ZjSJaGQYRgKYodeMC832HR+jLV4VjU+MQmKMBiaLLWxfMA29ZtQ6jn2fGOQzI1qVAhPFTqwhUsTq2gJ+YRAs0YiiNUrcPs3eHv6dUi7VELmDYSSAnXC8uWSMifDaaaVivxEA+PR8eG7xFWdFUvVa6COCvydEqB2ad76pE/WGXqiVctyYE1z9cMSlmveO1fKEbCI3MIwEsOPOepFYrz6uWGgXrv1GQqVeROQsYg2RMCLOityYHQetKvBawA9nwYQEJMVo0NZtxicVjVIPhyhoMIwEsL56Ee8s0YjEHTWlYVo3Eir1IiJnEeuV9pAonHQu0QTBlt5rKeQy3FPAniNE7mIYCWDHHZ1XZ3ipXkSUkxCFpBgNTFYbjlZd9epjB4MzQX5A3rVyU3VQymVo6TKhrj24i1hNFptzxi4Y+osM5l5HA7R/nGlAW7dJ4tEQBQeGkQDV2mVCdWsPAGCal8OITCZzdmMNt34jgiAE/QF519KqFM5C3GBfqjladRXdJivio9TITdFJPRyPTEnRITdVB7NVwF+Psz080UgwjAQo8U1lXEKUT4osxSLWz8KsiLWp04ir3WbIZcCEpMA978RdM1yWaoKZuKV3wYQEyOXB2xn3Pmd7eJ7kSzQSDCMBSuyo6e16EZFYxFpW3YYeU/jUjVTW28/5yY6PCrriyOHkhUgR675zwdMCfjhfzk+DXGbv53PJsT2fiIbGMBKgxIPPvNVf5FpZ8ZFI0WnDrm7kTL29xXio1IuIZrhs7w3WIta2bhOOO0L4oonB0+xsMEk6rfNn2MWeI0TXxTASoMQwMsPL23pF/c+pCZ+lGucBeSFSLyKanBIDpVyG1i6Ts6FbsDl4vgWCAExMikaKPjhawA9HLGR9p/RK0AZEIn9hGAlAjYZe1Bt6IZcB09J8V8TXV8QaPmGkosG+TBMq23pFWpXCWZAbrIfm7TvbBCA4t/QOpmhqCqLUClS39uDI5fCZfSTyBMNIABLX/SckRSNKo/TZ88wNs7oRm03A2RBreOYqmDuxuraAD/Z6EVGEWoE7pqcCYCEr0fUwjASgvnqRWJ8+z9i4SKTptTBbBZSEwSe3K1d70G2yQq2QIzs+UurheJ1Y7Hw8CHfUXG7pxpWrPVApZJiTEy/1cLzmK46lmr8erwvLBoNEI8UwEoDEIr78TN8Ur4pcz6kJh6UasQ38+KRoKBWh96cvzoycCMIiVnGJZtbYMT6dDfS3uTnxSNNr0dFrwT9Osz08YJ8F27b/In723knuNCKn0HtFDnKCIDin2X21k8aV89C8MAgjYvFqqNWLiCanxEClkOFqtxlXrvZIPRy3iEs0iycF9y6aa8n7tYfnUo0gCPjVRxV4+q+nsP3gJdyyaQ/W/V8ZLjKUhD2GkQBT196L5k4TlHIZclN934FSLGI9Vt2GbpPF588npYoQawN/LY2yr4j1RBDVjVisNucJ0gsnhEa9iCtxV82eyiY0dxolHo10BEHAb/5WgZf2nAcAFIyNhdUmYNfRGtzyzCdYt6MMF5o6JR4lSYVhJMCISzSTkmP80pQrY0wE0mMjYLEJOHIptOtGxDAyOSV0Oq9eS6wzOh5EYeTYlTZ0GC3QR6ich/6FkglJMcjP0MNiE/D+sVqphyMJQRCwqbgSWz6xB5GfLZ2Kd1YtwF9WL8AtU5JgE+z9WG7dtAdP7SjDeYaSsMMwEmDE4kNf14uIwqVuxGSxOV/gQnVmBOhfNxIsnKf0TkiAIohbwA/nKwViz5HwbIC2+e9n8fw/zwEAfnLXVDy6IAcAkJ8Zi1cfvQHvPbEAt+baQ8k7pTW4bdMerP1jKUNJGGEYCTB99SKxfnvOcGh+dqmlCxabgGiNEumxEVIPx2dmuOyoCZYiVmcYCZEtvYNZmp8GpVyG41fandvLw8Xmv1fi2X+cBQD8+M5cfHNhzoDbzMiIxSuP3ID3n1iIW3OTYROAd8tqcdumPXjyj6U418hQEuoYRgKIIAgunVf9N10tzowcv9KOLmNo1o301YtEQyYLzU/fgH3WR62Qo70nOIpYDb1mlFW3AQjNehFRfLQGN08Ov/bwz/3jLDb/3R5ENtyRi8cXjRv29tMz9Hjlkdn463cW4rap9lDyl7Ja3PY/e7Dm7VKcawyvIBdOGEYCSFVrN9p7zFAr5H5dSsiMi0TGGEfdSIj2G+mrFwndJRoAUCvlzp8xGPqNfHa+BVabgJyEKGTGhV7vF1f3zsoAALxbWgObLThmrUbjhX+exabiSgDA+iVT8K3FwwcRV3npevx+uT2UFE1NhiAA7x2rxW3/sxffebs07GaXwgHDSAAR3zxy03RQK/37q3Fu8T0fmks1FSHcefVa04PoBF/XepFQ98UpSdBplahr7w3pJVEAePFf5/Dbj+1B5AdfmoJv3zTeo8fJS9fj5eWz8cGahbh9mj2UvH+sFkWb9+KJt446t+tT8GMYCSDim8cMCXYUzAvxIlbxRSvUZ0aA4Cpi3X8utFrAD0erUuDOGWkAQnupZusn5/Gbv1UAAL5/+2T8+82eBRFX09L0+N3Ds7F7zSJ8aVoKBMHe1fb2zXux+q2jzplPCl4MIwFE3NY73Y/1IqI5jiLW8pp2dIZY3Ui3yYKq1m4AwORwmBlxOaMmkItYq1u7cbG5Cwq5DHPHh04L+OHc5+g58mF5XUj29fndnvP41UdnAADfK5qE1V+Y4NXHn5qmw0sPF+LDJxdhSZ49lHwghpI3GUqCGcNIgLDZBJyoMQDwb/GqKGNMJDLjImC1CTh8qdXvz+9LZxs6IQhAQrQa8dEaqYfjc65FrNWtgVvEKs6KzMyMhU6rkng0/lGYNQZj4yLRZbLi45MNUg/Hq36/9wI2fmgPIutum4QnvjjRZ8+Vm6rD1m8U4qO1i3DH9BQAwAfl9lCy6s0SnKk3+Oy5yTcYRgLEheYudBotiFApMCFRmqZcobpUUxFGSzSAvYg1N9VRxFrTJu1ghrE/xE7pHQmZTObsORJKSzWv7LuAX+w+DQBYe+tErLnFd0HE1ZQUHbZ83R5K7pyeCpkM2F1ejy9t3oeVfyjBqVqGkmDhURjZsmULcnJyoNVqUVhYiH379g17+zfffBP5+fmIjIxEamoqHnvsMbS0hNYb3miVO940pqXpJDvEzdn8LMSKWCtDvA38YPLSA7uI1WoTwqpexJUYRvafbUKjoVfi0Yzetv0X8V8f2IPImlsmYu2tk/w+hikpOrz49Vn46MnFuHOGPZR8dLIedzy3D9/+wxGcrA3MfwfUx+13vR07dmDt2rXYsGEDSktLsWjRIixZsgRVVVWD3n7//v1Yvnw5VqxYgZMnT+JPf/oTDh8+jMcff3zUgw8l4k4aKepFRGIYKa9pR0evWbJxeJtzZiSMwoi41FceoNt7T9S0o73HjBiNEvkZsVIPx6+yE6JQmDXG2UMjmG0/cBFP//UUAOA7X5yAp271z4zIUCanxODFh2bhb2sX4y5HKPnbyQbc+dx+/NvrDCWBzO0wsmnTJqxYsQKPP/44cnNzsXnzZmRmZmLr1q2D3v6zzz5DdnY21qxZg5ycHCxcuBDf/va3ceTIkVEPPpRI0ezsWmmxEciKj4RNQEidU+NseBYmyzRA/5mRQCxi3Xe2CYD9oEapZgKlJB6etzOIT/J9/dNL+Nn79iCy+gvjse62SQHTUHBScgxeeGgWPl67GEvz0yCTAR+fsoeSb71+JCh2moUbt14FTCYTSkpKUFRU1O/6oqIiHDx4cND7zJ8/H1euXMHu3bshCAIaGhrw5z//GXfeeeeQz2M0GmEwGPpdQpnFanMm9hkSf0qcm+PoNxIidSNXu0xo7LCflBpOyzSTkmOgVsrR0WvB5ZZuqYczwL4wrBdxddf0NKgVcpyp7wjKuoY/fHoJP/nLSQDAv988Ht8rmhwwQcTVxOQYPP+1AhQ/tRh3z7SHkuJTDbjr+f14/H8ZSgKJW2GkubkZVqsVycnJ/a5PTk5GfX39oPeZP38+3nzzTSxbtgxqtRopKSmIjY3F888/P+TzbNy4EXq93nnJzMx0Z5hB51xTJ3rNNkRrlMiJj5J0LPPGh1YRq9hfJGNMBKI1SolH4z8qhRy5qToAgVc30mW04GiVfeZt0cREiUcjDX2kCrfkJgEA3ikNrtmRNz67jP9wBJFvLx6H/3d7YAYRVxOSYvDsgwUofuom3DMzDXIZ8PfTYig5HLDLmeHEo/nRa//wBEEY8o/x1KlTWLNmDX7yk5+gpKQEH330ES5evIiVK1cO+fjr169He3u781JdXe3JMIOGuESTl66DXOJTS8W6kRM17TCEQN1IONaLiGYEaBHr5xdbYLYKyBhjXxYMV2Ih67tltbBYbRKPZmTe+rwKP373BADgW4ty8MMlUwI+iLiakBSNzQ8WoHjdTfhKQbojlDRi6Qv7sWL7YWevJ/I/t8JIQkICFArFgFmQxsbGAbMloo0bN2LBggX4/ve/jxkzZuD222/Hli1bsG3bNtTV1Q16H41GA51O1+8SysR/AIFQyJei1yInIQo2ATh8Mfj7jYRjvYjI2fwswD719S3RJAbVG5m33Tw5CWMiVWjqMOJAEOxg23G4Cj96pxwAsGJhDn50R27Q/v7GJ0bjf5bNxN/X3YR7HaHkH2ca8eUXDuCb2w/jmOPwRvIft8KIWq1GYWEhiouL+11fXFyM+fPnD3qf7u5uyOX9n0ahUABAQBbWSaE8AHbSuJrr6MYaCks14jLNlHAMIxl9beED6WC2cK8XEamVcnw539EePsALWf/vcDV+uMseRB5bkI0f3xm8QcTVuMRobBJDySx7KPnnmUbc/eIBPPbaIeeJ0uR7bi/TrFu3Dq+88gq2bduG06dP46mnnkJVVZVz2WX9+vVYvny58/ZLly7Frl27sHXrVly4cAEHDhzAmjVrcOONNyItLc17P0mQMllsOF1nf8OckR4r7WAcnIfmBXkYEQShb2YkDJdpJiZFQ6OUo8NoweXWwChirWvvwbnGTshkwPwwaQE/nK84TvL928n6gD2G4U9HqvGDXcchCMCj87Pxk7umhkQQcTUuMRqbHpiJf3z3Ztw3KwMKuQz/qmjCPS8ewKOvHUJpVejsLgxUboeRZcuWYfPmzXj66acxc+ZM7N27F7t370ZWVhYAoK6url/PkUcffRSbNm3CCy+8gLy8PNx///2YPHkydu3a5b2fIohVNnTAZLUhNlKFzLgIqYcDoC+MnKw1oL0neOtG6g29MPRaoJDLMC5R2sJgKSgVckxNsy9xBspauNh1dUZGLGIj1RKPRnr5GXqMS4xCr9mGD8sHX7aW0s6SK/h/O+1BZPm8LPx0aegFEVc5CVF45oF8/GPdTfhqoT2UfFLRhK9sOYhHth1yFl6T93lUwLpq1SpcunQJRqMRJSUlWLx4sfN727dvxyeffNLv9t/5zndw8uRJdHd3o7a2Fm+88QbS09NHNfBQcUw8HC9dHzD/yJN1WoxLiIIgAIeCuG5EnBXJSYiCRqmQeDTSCLQTfJ1LNBPCe4lGJJPJcK/YHv5oYLWHf6f0Cr7352MQBOAbc8fi51+eFjCvUb6WnRCF396fj39+9ybc7wgleyqbcO+Wg1i+7RBKLjOUeFv4dRsKMOUB0OxsMHNDYItvZZidSTMYMYwcD4AiVptNwIEwbQE/nHscYeSziy2oaQuMgw3fLa3Bd//PHkS+Pmcsnv5yXtgEEVdZ8VH4zf35+Nd3b8YDs+2hZG9lE+7behAPv/o5Si4H74e1QMMwIjFnG/gAqRcRzQ2BQ/Mq6jsBhOe2XpFYxHqy1iB5EeupOgNaukyIVCtQMHaMpGMJJBljIjEnJw6CYA8BUvtLWQ3W/V8ZbALwtRvH4j/vzpO85YDUxsZH4tdftYeSZbMzoZTLsO9sM+7b+ikefvXzgAmRwYxhREK9ZquzD0bAzYzk2HfUnKozoK3bJPFoPFPRYO9sGY7Fq6IJidHQquToNFpwsaVL0rGIB+PNHRcPtZIvPa7ucxSy7jp6RdJdhu8fq8VTO+xB5MEbMvGLexhEXI2Nj8SvvjoD//rezXjwhr5Qcudz+/DPMw1SDy+o8RVBQqfqDLDaBCREa5Cq10o9nH6SdFqMTwzeuhGrTcDZBvvMSDhu6xUpFXJMdXRilbpuZD+39A5pyfQUaJRynG/qkqxJ3V+P12KtI4g8MDsDv/zKdAaRIWTGReK/75uBv6+7CdPT9WjrNuOb24/gvz88EzQN7AINw4iEXOtFAnE9Npi3+Fa1dsNosUGrkiMzLny7fAJ95x1J2fys12zFoUv2UMswMlCMVoWiaSkApClk3V1ehyf/WAarTcBXCzPw3/fOYBAZgeyEKPz53+fhkXn23aQv7TmPr/3+M9S390o8suDDMCKhvnqRwFqiEfWdUxN8MyMV9fYlmolJMVCE+YuqeILvcQlnRg5dbIXJYkOqXovxidGSjSOQiSf5vnesFmY/frr+sLwO33m7FFabgPtmZeBX9zGIuEOjVODnd+fhxYdmIVqjxOFLV3HHc/uwt7JJ6qEFFYYRCYm9HwKtXkQ0x3GC75n64KsbEYtXw7leRCSG3ZMSdmLdd9b+wrxwQkJAzgIGgkUTEpAQrUFrlwl7KvzzRvbRiXpnELm3IB2//uqMsA/vnrpzRire/85CTE3VobXLhEdeO4RnPq6ANYC6HwcyhhGJdBktONdkf8MMlDbw10qM0WBiUjQEIfhmR8K5Dfy1xidGIUKlQJfJigvN0hSxiv1FFnKJZkhKhRx3z3S0h/fDSb4fn6zHE28dhcUm4O6ZafjN/fkMIqOUkxCFXavm46E5YyEIwPP/PIdvvPI5Gju4bHM9DCMSOVlrgCAAqXotkmICq3jVVbBu8RV3KYXjAXnXcu3EKkURa2NHL844GtAtZLOzYYlLNX8/3Yj2bt91Py4+1YDVjiDy5fw0PMMg4jValQK//Mp0PPvgTESqFfj0QgvueHY/Djp2k9HgGEYkctyl82ogC8YwYrRYcdExAxDOPUZcSdn8TGx0Ni1Nh/hojd+fP5hMTdVhcnIMTBYbPvBRe/h/nG7AqjdLYLYKuGtGKjY9kA+lgm8F3nb3zHS898RCTE6OQXOnEd949XM894+zXLYZAv8CJXI8QDuvXmuO4wTfM/UdaO0KjrqR841dsNoE6LRKJOv45gdI2xa+75TeRL8/d7CRyWTO2ZF3fLBU868zjfj3N47CbBVw5/RUbF42k0HEhyYkRePd1Qtwf2EGbAKwqbgSj752CM2dRqmHFnD4VygRsZeAuO0yUCVEazAp2b774dDF4Jgd6asX0bFY0kEMvSdq2/36yUwQBPYXcdPdM9MhkwGHL13FZS82qvukohHf/kMJTFYb7piegs0PMoj4Q4Ragd/cn4/f3p8PrUrubJL2eRDNNvsD/xIl0N5jdi4jBPoyDQDME/uNnA+Ofzx99SLcQioalxiNSLUC3SYrLjZ3+u15Kxs60dhhhEYpR2EWW8CPRIpe66ytecdL7eH3VDbh3xxB5EvTUvDsgwVQMYj41VcLM/CX1QsxPjEKDQYjHnrlc2z55JzkxzQECv41SuCkY1YkMy4CY6IC/xj1vrqR4NhRI57Wy3qRPgq5DNMcRaz+rBsRt/TOGRcPrSo8T072RN9STc2o28PvrWzCt14/ApPFhtunJeP5hxhEpDI5JQbvPbEQXylIh9Um4NcfVWDF/x7G1SBZAvcl/kVK4JhYLxJgh+MNZY4jjFQ0dKAlCNY6xTDCHiP9ic3P/Nlu3Fkvwl00brl9Wgoi1QpcbunG0SrPj6vff7bZGURum5qM5782i0FEYlEaJTY9kI//vnc6NEo5/lXRhDue2xf2JwDzr1IC5TVtAAK/eFUUF6V29uv4PMDPqenoNTtP0JzMbb39iH9v/moLb7RY8bmjzmjRJIYRd0SqlfhS3ujawx8414wV/3sYRosNt+Ym4cWHZvGAwgAhk8nw4I1j8e7qBchJiEJdey+W/e4z/H7vBUkPSpQS/zIl4GwDHyRhBAieLb5nG+31EMk6DWIjA38JzJ+cnVhrDX4pYi25fBW9ZhsSYzRcMvPAvQX2k3zfP1YLo8Xq1n0Pnu8LIrdMScKLX2cQCUS5qTq8/52FWJqfBotNwC92n8a3Xi/xaY+ZQMW/Tj9r7TLhylX7J/e8ICheFc0NkiJWLtEMLSchGlFqBXrMVpxv8n0Rq7PrKlvAe2Te+Hik6LQw9Frwz9ONI77fp+db8M3th9FrtuELkxOx5RuzoFGyXidQRWuUeO7BmfjPe/KgVsjx99MNuPP5fThW3Sb10PyKYcTPxGZn4xKjoNOqpB2MG+bk2PuNnG3sDOg98ixeHZq9iNV/SzXc0js6CrkMdxeI7eFHtlTz+YW+IHLTpERs/UYhg0gQkMlkeHhuFnatmo+xcZG4crUHX33pIF47cDFslm0YRvys3Fm8GjyzIgAwxqVuJJCXasQeI6wXGZy4NOjrItbWLhNO1Nqfgy3gPScu1fzrTON1mw4eutiKx7YfRo/ZikUTE/C7hwu5gynI5KXr8dc1C/GlaSkwWwX8/P1TWPXmURh6Q3/ZhmHEz8Rj3KcHeLOzwcwbH/h1Iwwjw5vupx01B841QxDsBxUm6QL37KVANzklBnnpOlhsAv56vHbI2x2+1IpHXzuEbpM9iPx++WwGkSCl06qw9Ruz8NOlU6FSyPDhiXosfX6/JN2T/YlhxM/EZZpg2UnjKtD7jTR3GtHcaYJMBkxMYhgZjDgzcqrWAIvV5rPn2e9SL0Kj8xXH7MjOIXbVlFxuxaPb7EFk4QQGkVAgk8nw2IIc/GnlfKTHRuBySzfu3XIQb3x2OWSXbRhG/KjB0IsGgxFyGZwNqILJnJw4yGTAucbOgDwSu9JRL5IVF4kINV+MB5MTH4VojdJRxOq9VuOuBEFwNjtbyHqRUftyfhoUchmOVbcNKDwuuXwVj2w7jC6TFfPHxzOIhJiZmbH4YM1C3JqbBJPVhh+/ewJr/liGTqNF6qF5HcOIH4n1IhOTYhCpVko8GvfFRqqRm2IPUZ8H4OyIsw08i1eHJHfpxOqrpZoLzV2obe+FWiHHnJx4nzxHOEmM0WCxI9S94zI7crTqKh7ZdgidRgvmjYvHq4/cwBAegmIj1fj98tnYcEcuFHIZ3j9Wiy8/vx+n6wxSD82rGEb8qK9eJPiWaETOLb4BWDfi3EnDepFhOetGHEuG3rav0j4rMjt7DN8cveTeWfalmndKa2CzCSirbsMjr9qDyJycOLz66Gz+vw5hMpkM31o8Dv/37blI1WtxobkL97x4ADsOV4XMsg3DiB8Fc72IKJCLWDkzMjK+3lGz/5y4pTfRJ48fjm6bmowYjRI1bT14df9FPPzq5+gwWnBjThxee+yGoJxpJfcVZsXhgzWLcPPkRBgtNvxgZzm++3/H0G0K/mUbhhE/EQShb1tvEO6kEd2Yba8budDUhUZD4NSNCILgrBmZwpmRYYkzI6fqvF/EarbanI3x2F/Ee7QqBe6YngoA+MXu0+joteCG7DF47VEGkXATF6XGtkduwPdvnwy5zN6D5u4XDuCs48NYsGIY8ZPa9l60dJmglMuC+s1SH6ly1hwE0lJNTVsPukxWqBQyZCdEST2cgJbtKGLtNdtwzsudWEur2tBlsiIuSo2pqcFXpB3IxJN8AWB21hi89tiNiNIwiIQjuVyG1V+YgLe/NRdJMRqcbezEl184gJ0lV6QemscYRvxEXJ+fnBIT9NXuc3MCb4uvWC8yPjGap5Jeh1wuQ166PSgc93In1v2OXTQLJiRALmcLeG+6ITsOX85Pw+3TkrH9mzcimkEk7M0ZF4/dTy7CwgkJ6DFb8d0/HcMP/nwcvWb3zjIKBHzV9pNjziWa4K0XEQXioXmsF3GPuFTj7UZKe8UW8Owv4nVyuQzPfa0Av3t4NoMIOSVEa/C/37wRT906CTIZsONINe558YBfzp/yJoYRPwmFehHRDTlxkMuAi81dqG8PjLqRSu6kcYvYAdibMyPt3WZnkTb7ixD5j0Iuw5O3TsQbK+YgIVqNM/Ud+PLz+/GXspGdaRQIGEb8QBAE54v09CA7k2Yw+giV88C1zy8GxuxIRYP9UwAPyBsZ8e/wdJ0BZi8VsX56oRk2ARifGIW02AivPCYRjdyCCQnYvWYR5o6LQ5fJiif/WIYN75QHxbKNR2Fky5YtyMnJgVarRWFhIfbt2zfkbR999FHIZLIBl2nTpnk86GBT1doNQ68FaqU8ZJYRxC2+4s4JKZmtNpxvdIQRzoyMSFZcJGK0ShgtNpxt8M50rnOJhlt6iSSTpNPijRVz8J0vToBMBrz5eRXu23oQl5p903HZW9wOIzt27MDatWuxYcMGlJaWYtGiRViyZAmqqqoGvf2zzz6Luro656W6uhpxcXG4//77Rz34YCHWi+Sm6qBWhsZk1NxxcQACo27kcksXTFYbItUKpPMT+YjI5TLkpXm3bmS/M4xwiYZISkqFHN8tmoztj92IuCg1TtYasPT5/dhdXif10Ibk9jvjpk2bsGLFCjz++OPIzc3F5s2bkZmZia1btw56e71ej5SUFOflyJEjuHr1Kh577LFRDz5YiDtp8kOgeFU0O9teN3KppRt17T2SjqWi3v7JflJyDHdwuEEspj5e0zbqx7rc0oWq1m4o5TLMGccW8ESB4KZJifhgzULckD0GHUYLVr15FD/9ywkYLYG3bONWGDGZTCgpKUFRUVG/64uKinDw4MERPcarr76KW2+9FVlZWUPexmg0wmAw9LsEM7FIMBTqRUQ6rcr580g9OyLupGG9iHvyxLbwNaP/97XPMSsyK2sMd3oQBZBUfQTe/tZcrLxpPADgfz+9jPtf+hTVrd0Sj6w/t8JIc3MzrFYrkpOT+12fnJyM+vr6696/rq4OH374IR5//PFhb7dx40bo9XrnJTMz051hBhSbTXBOg4fCThpXzi2+56XtN1JRb38zncR6EbeIMyPeKGLdzy29RAFLqZDjh0umYNujsxEbqcLxK+2487l9+Pjk9d+3/cWjAgaZrP9UuCAIA64bzPbt2xEbG4t77rln2NutX78e7e3tzkt1dbUnwwwIF5o70WWyIkKlwISkaKmH41VzxwfGoXmVjgLMYO5sK4WxcZHQaZUwWWyoHEUraYvVhgPn7WGEW3qJAtcXpyTjgzWLUDA2FoZeC/7tDyX4r7+e8tqOutFwK4wkJCRAoVAMmAVpbGwcMFtyLUEQsG3bNjz88MNQq9XD3laj0UCn0/W7BCtxiSYvXQdFiNUzzM4aA4VchqrWbtS0SVM30mu24lKLvUo8VHYq+YtMJnMemjeaItbjNe3o6LVAp1WG3OwfUahJj43Ajn+bh8cX5gAAXtl/EQ/87lPJXsNFboURtVqNwsJCFBcX97u+uLgY8+fPH/a+e/bswblz57BixQr3RxnE+upFYqUdiA/EaFXOuoPPJNrie66xE4JgPzwqIXr4kEsDib+/0TQ/21dpnxVZMCEh5AI3UShSK+X48V1T8buHCxGjVaK0qg13PrcP/6polGxMbi/TrFu3Dq+88gq2bduG06dP46mnnkJVVRVWrlwJwL7Esnz58gH3e/XVVzFnzhzk5eWNftRBpLwmdNrAD2aexK3hz9SLbeCjR7RUSP3NcITk0cyM7D9nP4+G/UWIgsvt01Kwe80izMjQo63bDKOEzdHcLntftmwZWlpa8PTTT6Ourg55eXnYvXu3c3dMXV3dgJ4j7e3t2LlzJ5599lnvjDpIWKw2nKwN7TAyd1wcXtpzHp9J1IlVrHWYkhK8S3lS6uvE2gGTxeZ2H5yOXjNKq9oAsL8IUTDKjIvEn1bOw7/ONOJLeamSjcOjPXirVq3CqlWrBv3e9u3bB1yn1+vR3R1Y24j84WxjJ3rNNsRolMiOD81j7W/IjoNCLkN1aw+uXO1GxphIvz5/RT0PyBuNzLgI6CNUaO8xo7Khw7lsM1KfXWiFxSYgKz4SmXH+/d0TkXdolApJgwjAs2l8qtxZvKoP2WZcURqlc9bnswv+3+Jb4TwgL7R2KvmLTCZzzo6Ue7BUs/+suETDWREi8hzDiA8dc3ReDdUlGpHYb8Tf59S0d5tRb7CfGjyRMyMeE3fUeBJGxGZnCyewXoSIPMcw4kPlIdrs7FpSFbFWNtpnRdJjI6DTqvz63KHEOTPi5o6amrYeXGjugkIucx6cSETkCYYRHzFarDhdZ+8MGuozI4VZY6CUy1DT1uPXFsMVLjtpyHNiGDlTb3DrzApxiSY/Qw99BMMgEXmOYcRHKus7YbYKiI1UIWNMaJ8kG6VRIj8zFoB/u7E6wwg7r45KxpgIxEaqYLYKqHQcOjgSe8UlGm7pJaJRYhjxEbFeZHq6Piz6X8wdFwfAv0s1PCDPOzwpYrXaBBw4Zw8ji1m8SkSjxDDiI+L6e36I14uI+g7Na4EgCD5/PkEQnD1GJnNmZNT6wkjbiG5/srYdbd1mRLvMihEReYphxEeOOz5hTg/xehFRYdYYqBQy1Lb3orrV92ccNHUY0dZthlwGjE9kzchouTszIu6imTc+HioFX0aIaHT4KuIDvWar81N7qBeviiLVSucskD+WasQ28NkJUdCqFD5/vlAnhuaK+o4RFbHud4QR9hchIm9gGPGBk7UGWG0CEqI1SNFppR6O34jbO/1RxNrXBp5LNN6QHhuBMY4iVrEweCjdJguOXLY3uFs4gWGEiEaPYcQHyh3Fq/kZ4VG8Kprr0m/E13UjbAPvXTKZDNMdM1vXO8H384utMFsFpMdGICchNI85ICL/YhjxgXCrFxHNGmuvG6lr78XlFt/2G6nkThqvm55uP2zweif4ui7RhFPYJiLfYRjxAXEnTbjUi4gi1AoUZI4B4Nu6EZtNQGWDvR8Ge4x4z/T0WADXnxnZ5zyPhv1FiMg7GEa8rNNowbkm+xul+OIeTvzRb6T6ajd6zFaolfKQPQ1ZCuJMXmVDB3rNgxexNhh6UdnQCZkMmM8W8ETkJQwjXnayph2CAKTptUiM0Ug9HL+b61LE6qu6EbFeZGJSNBQhehqyFNL0WsRHqWGxDV3EKm7pnZ6ux5gotT+HR0QhjGHEy8rDtF5ENGvsGKgVcjQYjLjko7oR1ov4hkwmQ56j38jxIepG9juXaLiLhoi8h2HEy45fCY+TeoeiVSkwc2wsAODT875ZqjnDM2l8RqxzOjFI3YjNJmD/OfvvdOEE1osQkfcwjHjZcce23nArXnU1z2WLry+wDbzvDDczcqa+A82dRkSqFZiVFevnkRFRKGMY8aL2brNzaUJsrx2OfNlvxGSx4UJTFwAu0/iCGKLPDlLEuv+cfYlmTk4cNEp2vSUi72EY8aITtfZPk2PjIhEbGb7FfQVjY6FWytHYYcSF5i6vPvaF5k5YbAJiNEqk6sOnu62/pOi0SIi2F7GerjP0+55YvLqQW3qJyMsYRrzomGOJJlyLV0ValQKzHHUj3l6qqXCpF2HDLe+TyWTOWT3X5me9ZisOXbS3gF/M4lUi8jKGES8Sm53lh3kYAfqWarxdxMp6Ed8Tw4hr87Mjl67CaLEhWafBhCSekkxE3sUw4kXii3c4Nju7Vl8Ra6tX60Yq6u0N5Vgv4jviGTXlLjMjYtfVhRMSOSNFRF7HMOIlLZ1G1LT1AADyHGd8hLOZY2OhUcrR3GnE+Sbv1Y1UNNjrGHhAnu+IMyNnGzudRaxivcjiSVyiISLvYxjxEnEr5LjEKMRoVRKPRnoapQKFWfZzaj71Ut1Il9GC6lZ74OMyje8k6zRIjNHAahNwqs6Apg4jTjmKWRdMYBghIu9jGPGSvnqRWGkHEkDmernfyNlG+xJNYowGcWxF7jOuRazlV9px8Lx9VmRqqg4J0eF3xAER+R7DiJf01YuweFUkhpHPvdRvpLKebeD9RWx+Vl7T7lyiYQt4IvIVpdQDCBXlNW0Awrvz6rXyM/XQquRo7jThXGMnJo4yRDjbwDOM+NwMl5mRth4TAGAhwwgR+QhnRrygwdCLBoMRchkwNY3FqyLXuhFvLNWI23qnsF7E58ReORUNHWgwGKFRynFDdpzEoyKiUMUw4gXiEs2k5BhEqjnZ5Erc4uuNItaKBh6Q5y/JOi2SYvrqQ27MiYNWxRbwROQbDCNeUC52XmW9yABzvdRvpLXLhKYOIwBgIptu+YXr3zPrRYjIlxhGvEDc1st6kYFmZMQiQqVAa5cJlQ2dHj+O2AY+My4CURrOPvmD67EGCyfwPBoi8h2PwsiWLVuQk5MDrVaLwsJC7Nu3b9jbG41GbNiwAVlZWdBoNBg/fjy2bdvm0YADjSAIfTtpuK13ALVSjtnZo68bcbaBT2ZNjr+I29QTojWs0yEin3L7I+aOHTuwdu1abNmyBQsWLMDvfvc7LFmyBKdOncLYsWMHvc8DDzyAhoYGvPrqq5gwYQIaGxthsVhGPfhAUNPWg9YuE1QKGXJT+YI9mLnj4rHvbDM+u9CCR+Zne/QYFc4zabhE4y+LJyVizRcnYFbWGMjlbAFPRL7jdhjZtGkTVqxYgccffxwAsHnzZvztb3/D1q1bsXHjxgG3/+ijj7Bnzx5cuHABcXH2avzs7OzRjTqAiM3OJqfEQKNkgd9gXJuf2WyCR29sldzW63cKuQzriiZLPQwiCgNuLdOYTCaUlJSgqKio3/VFRUU4ePDgoPd57733MHv2bPz6179Geno6Jk2ahO9973vo6ekZ8nmMRiMMBkO/S6AS60V4ON7QZmToEaFS4Gq3GZWNHW7fXxAEl5kRhhEiolDjVhhpbm6G1WpFcnJyv+uTk5NRX18/6H0uXLiA/fv348SJE3jnnXewefNm/PnPf8bq1auHfJ6NGzdCr9c7L5mZme4M06+OO3bS5LN4dUgqRV/dyKfn3a8bqWvvRUevBUq5DOMSuExDRBRqPCpgvfYIcUEQhjxW3GazQSaT4c0338SNN96IO+64A5s2bcL27duHnB1Zv3492tvbnZfq6mpPhulz/YtXGUaGM2+85+fUiLMi4xKjoFZyAxgRUahxq2YkISEBCoViwCxIY2PjgNkSUWpqKtLT06HX971Z5+bmQhAEXLlyBRMnThxwH41GA40m8A/kutzSjY5eC9RKOWsZrsN5Ts3FVrfrRipYL0JEFNLc+pipVqtRWFiI4uLiftcXFxdj/vz5g95nwYIFqK2tRWdnX4+JyspKyOVyZGRkeDDkwCHWi0xN1UGl4Cf24UxP1yNKrUBbt9l5xsxI8YA8IqLQ5vY76Lp16/DKK69g27ZtOH36NJ566ilUVVVh5cqVAOxLLMuXL3fe/qGHHkJ8fDwee+wxnDp1Cnv37sX3v/99fPOb30RERIT3fhIJHK9uA8B6kZGw143Yd1O5u1TD4lUiotDmdhhZtmwZNm/ejKeffhozZ87E3r17sXv3bmRlZQEA6urqUFVV5bx9dHQ0iouL0dbWhtmzZ+PrX/86li5diueee857P4VEnDtp2OxsROZ6cE6N1SbgbKN9Vo1hhIgoNHnUV3vVqlVYtWrVoN/bvn37gOumTJkyYGkn2FltAk6yDbxbxCLWQ27UjVxq6YLJYkOESoHMMZG+HiIREUmAhQ4eutDUiS6TFZFqBcYncrvpSOSl6RCtUaK9x4zT9SPrHdPX7CyaXUCJiEIUw4iHxC29eWl6KPgmOSJKhRw3uNlvRKwX4U4aIqLQxTDiofIa9hfxRF9r+NYR3b6SxatERCGPYcRDYudV1ou4p6/fSAusNuG6tz/DHiNERCGPYcQDZqsNJ2vtNQ/T0xlG3DEtTYcYjRIdvRacrhu+bqTXbMWl5i4A4BH2REQhjGHEA2cbOmG02BCjVSI7Pkrq4QQVpUKOG3JG1m/kfFMnbAIQG6lCYkzgd+QlIiLPMIx4oLymDYB9VoQ7PNw3T+w3cp0i1kqX4tWhzj4iIqLgxzDiAR6ONzpi3cihi63D1o2cYRt4IqKwwDDiATGMzEiPlXYgQWpqmg4xWiU6jBacrG0f8nbOM2lYL0JEFNIYRtxktFhxxtGwiztpPKOQyzBnBHUjlQ1sA09EFA4YRtxUUd8Bs1XAmEgVMsYE90F/Urpev5GOXjNq2noAAJOSGEaIiEIZw4ib+upFYllUOQqudSMWq23A98Xi1RSdFvpIlV/HRkRE/sUw4iax2Vk+l2hGJTdVB51WiU6jxdmzxVVFPZdoiIjCBcOIm5wzI2x2NioKuQw35ji2+A5SN8I28ERE4YNhxA09JivONto/sc/IiJV2MCFg3nixbmRgGBGLhNkGnogo9DGMuOFUnQFWm4DEGA2SdewIOlpzx9l31By+pm5EEARUOLb1sg08EVHoYxhxg2u9CItXRy83RQd9hApdJqvzFGQAaO404Wq3GTIZMCEpWsIREhGRPzCMuKHcWS8SK+1AQoS8X7+Rvi2+Yr1IdnwUtCqFJGMjIiL/YRhxw3HHp3c2O/MecYuvaxGr2AZ+UjJnRYiIwgHDyAh1Gi0432QvXuWZNN4jFrEeudQKs6NupK8NvE6ycRERkf8wjIzQiZp2CAKQHhuBhGgWr3rL5OQYjIlUodulbqSigQfkERGFE4aRESpnfxGfsNeNOJZqzrfAZhNw1tljhMs0REThgGFkhMR6ES7ReJ+4xfezCy2oaetBl8kKtUKOrPgoiUdGRET+wDAyQuK2Xhavet9cZ93IVZxwhL7xSdFQKfjnSUQUDvhqPwLt3WZcbukGAMzgtl6vm5QUg7goNXrMVvy55AoAYDJ30hARhQ2GkREQCyuz4iN5gqwPuPYb+WdFIwBgEjuvEhGFDYaREThe0waAxau+JG7xFQT719xJQ0QUPhhGRuB4NZud+ZrY/EzE03qJiMIHw8gIlDs7r8ZKO5AQNjEpGvFRagBAlFqB9NgIiUdERET+wjByHc2dRtS09UAmA6alsSOor8hkMufsyKSUGB5ESEQURhhGrkOcFRmXEIUYLYtXfWnJ9BQAwKIJCRKPhIiI/Ekp9QACXV+9SKy0AwkDd81IQ26qDpljIqUeChER+ZFHMyNbtmxBTk4OtFotCgsLsW/fviFv+8knn0Amkw24nDlzxuNB+1O5YycNi1f9Y3xiNNRKTtgREYUTt1/1d+zYgbVr12LDhg0oLS3FokWLsGTJElRVVQ17v4qKCtTV1TkvEydO9HjQ/nT8CnfSEBER+ZLbYWTTpk1YsWIFHn/8ceTm5mLz5s3IzMzE1q1bh71fUlISUlJSnBeFQuHxoP2lwdCLxg4jFHIZpqYyjBAREfmCW2HEZDKhpKQERUVF/a4vKirCwYMHh71vQUEBUlNTccstt+Bf//rXsLc1Go0wGAz9LlI4Vt0GwL7tNEId+OGJiIgoGLkVRpqbm2G1WpGcnNzv+uTkZNTX1w96n9TUVLz88svYuXMndu3ahcmTJ+OWW27B3r17h3yejRs3Qq/XOy+ZmZnuDNNr+vqLcFaEiIjIVzzaTXNtDwhBEIbsCzF58mRMnjzZ+fW8efNQXV2N3/72t1i8ePGg91m/fj3WrVvn/NpgMEgSSMR6kencSUNEROQzbs2MJCQkQKFQDJgFaWxsHDBbMpy5c+fi7NmzQ35fo9FAp9P1u/ibIAh9MyM8k4aIiMhn3AojarUahYWFKC4u7nd9cXEx5s+fP+LHKS0tRWpqqjtP7XdXrvagtcsElUKGKak8J4WIiMhX3F6mWbduHR5++GHMnj0b8+bNw8svv4yqqiqsXLkSgH2JpaamBq+//joAYPPmzcjOzsa0adNgMpnwxhtvYOfOndi5c6d3fxIvE2dFpqTooFGyeJWIiMhX3A4jy5YtQ0tLC55++mnU1dUhLy8Pu3fvRlZWFgCgrq6uX88Rk8mE733ve6ipqUFERASmTZuGDz74AHfccYf3fgof6KsX4RINERGRL8kEQRCkHsT1GAwG6PV6tLe3+61+5KHff4aD51vw3/dOx4M3jvXLcxIREYWSkb5/s+/2IGw2l+JV7qQhIiLyKYaRQVxu7UZHrwUapRwTk6OlHg4REVFIYxgZxPErbQCAqWk6qBT8X0RERORLfKcdhPNwPPYXISIi8jmGkUGUX2G9CBERkb8wjFzDahNwopZn0hAREfkLw8g1LjR1ottkRaRagXGJLF4lIiLyNYaRaxxzLNHkpemhkA9++B8RERF5D8PINcodO2m4RENEROQfDCPXOF7DNvBERET+xDDiwmy14VStAQB30hAREfkLw4iLyoYOGC02xGiVyI6PlHo4REREYYFhxEVffxE9ZDIWrxIREfkDw4gLZ71Ieqy0AyEiIgojDCMuXGdGiIiIyD8YRhyMFivO1IvFqwwjRERE/sIw4nCmrgNmq4C4KDXSYyOkHg4REVHYYBhx6KsXYfEqERGRPzGMOLDzKhERkTQYRhyOO4tXY6UdCBERUZhhGAHQY7KisqEDAGdGiIiI/I1hBMCpunbYBCApRoNknVbq4RAREYUVhhEAx6rZX4SIiEgqDCMAymtYL0JERCQVhhEAxx07aaZzZoSIiMjvwj6MdPSacaG5C4C9xwgRERH5V9iHkRM1BggCkB4bgYRojdTDISIiCjthH0bKa9oAsHiViIhIKmEfRsRmZ6wXISIikkbYhxHnTpr0WGkHQkREFKbCOoy0dZtwuaUbAItXiYiIpBLWYUScFcmOj4Q+UiXxaIiIiMKTR2Fky5YtyMnJgVarRWFhIfbt2zei+x04cABKpRIzZ8705Gm9rq9eJFbagRAREYUxt8PIjh07sHbtWmzYsAGlpaVYtGgRlixZgqqqqmHv197ejuXLl+OWW27xeLDeVi6e1MslGiIiIsnIBEEQ3LnDnDlzMGvWLGzdutV5XW5uLu655x5s3LhxyPs9+OCDmDhxIhQKBd59912UlZWN+DkNBgP0ej3a29uh0+ncGe6wik814OD5Ztw3KwN5DCREREReNdL3b7dmRkwmE0pKSlBUVNTv+qKiIhw8eHDI+7322ms4f/48fvrTn47oeYxGIwwGQ7+LL9w2NRk/XTqNQYSIiEhCboWR5uZmWK1WJCcn97s+OTkZ9fX1g97n7Nmz+OEPf4g333wTSqVyRM+zceNG6PV65yUzM9OdYRIREVEQ8aiAVSaT9ftaEIQB1wGA1WrFQw89hJ///OeYNGnSiB9//fr1aG9vd16qq6s9GSYREREFgZFNVTgkJCRAoVAMmAVpbGwcMFsCAB0dHThy5AhKS0vxxBNPAABsNhsEQYBSqcTHH3+ML37xiwPup9FooNHwnBgiIqJw4NbMiFqtRmFhIYqLi/tdX1xcjPnz5w+4vU6nQ3l5OcrKypyXlStXYvLkySgrK8OcOXNGN3oiIiIKem7NjADAunXr8PDDD2P27NmYN28eXn75ZVRVVWHlypUA7EssNTU1eP311yGXy5GXl9fv/klJSdBqtQOuJyIiovDkdhhZtmwZWlpa8PTTT6Ourg55eXnYvXs3srKyAAB1dXXX7TlCREREJHK7z4gUfNVnhIiIiHzHJ31GiIiIiLyNYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUm73GZGCuPvYV6f3EhERkfeJ79vX6yISFGGko6MDAHh6LxERURDq6OiAXq8f8vtB0fTMZrOhtrYWMTExg54O7CmDwYDMzExUV1ezmVqA4O8ksPD3EVj4+wgs/H1cnyAI6OjoQFpaGuTyoStDgmJmRC6XIyMjw2ePr9Pp+IcUYPg7CSz8fQQW/j4CC38fwxtuRkTEAlYiIiKSFMMIERERSSqsw4hGo8FPf/pTaDQaqYdCDvydBBb+PgILfx+Bhb8P7wmKAlYiIiIKXWE9M0JERETSYxghIiIiSTGMEBERkaQYRoiIiEhSYR1GtmzZgpycHGi1WhQWFmLfvn1SDyksbdy4ETfccANiYmKQlJSEe+65BxUVFVIPixw2btwImUyGtWvXSj2UsFZTU4NvfOMbiI+PR2RkJGbOnImSkhKphxWWLBYLfvzjHyMnJwcREREYN24cnn76adhsNqmHFrTCNozs2LEDa9euxYYNG1BaWopFixZhyZIlqKqqknpoYWfPnj1YvXo1PvvsMxQXF8NisaCoqAhdXV1SDy3sHT58GC+//DJmzJgh9VDC2tWrV7FgwQKoVCp8+OGHOHXqFJ555hnExsZKPbSw9Ktf/QovvfQSXnjhBZw+fRq//vWv8Zvf/AbPP/+81EMLWmG7tXfOnDmYNWsWtm7d6rwuNzcX99xzDzZu3CjhyKipqQlJSUnYs2cPFi9eLPVwwlZnZydmzZqFLVu24L/+678wc+ZMbN68WephhaUf/vCHOHDgAGdvA8Rdd92F5ORkvPrqq87r7rvvPkRGRuIPf/iDhCMLXmE5M2IymVBSUoKioqJ+1xcVFeHgwYMSjYpE7e3tAIC4uDiJRxLeVq9ejTvvvBO33nqr1EMJe++99x5mz56N+++/H0lJSSgoKMDvf/97qYcVthYuXIh//OMfqKysBAAcO3YM+/fvxx133CHxyIJXUByU523Nzc2wWq1ITk7ud31ycjLq6+slGhUB9hMe161bh4ULFyIvL0/q4YStP/7xjzh69CgOHz4s9VAIwIULF7B161asW7cOP/rRj3Do0CGsWbMGGo0Gy5cvl3p4YecHP/gB2tvbMWXKFCgUClitVvziF7/A1772NamHFrTCMoyIZDJZv68FQRhwHfnXE088gePHj2P//v1SDyVsVVdX48knn8THH38MrVYr9XAIgM1mw+zZs/HLX/4SAFBQUICTJ09i69atDCMS2LFjB9544w289dZbmDZtGsrKyrB27VqkpaXhkUcekXp4QSksw0hCQgIUCsWAWZDGxsYBsyXkP9/5znfw3nvvYe/evcjIyJB6OGGrpKQEjY2NKCwsdF5ntVqxd+9evPDCCzAajVAoFBKOMPykpqZi6tSp/a7Lzc3Fzp07JRpRePv+97+PH/7wh3jwwQcBANOnT8fly5exceNGhhEPhWXNiFqtRmFhIYqLi/tdX1xcjPnz50s0qvAlCAKeeOIJ7Nq1C//85z+Rk5Mj9ZDC2i233ILy8nKUlZU5L7Nnz8bXv/51lJWVMYhIYMGCBQO2u1dWViIrK0uiEYW37u5uyOX93z4VCgW39o5CWM6MAMC6devw8MMPY/bs2Zg3bx5efvllVFVVYeXKlVIPLeysXr0ab731Fv7yl78gJibGOWOl1+sREREh8ejCT0xMzIB6naioKMTHx7OORyJPPfUU5s+fj1/+8pd44IEHcOjQIbz88st4+eWXpR5aWFq6dCl+8YtfYOzYsZg2bRpKS0uxadMmfPOb35R6aMFLCGMvvviikJWVJajVamHWrFnCnj17pB5SWAIw6OW1116TemjkcNNNNwlPPvmk1MMIa++//76Ql5cnaDQaYcqUKcLLL78s9ZDClsFgEJ588klh7NixglarFcaNGyds2LBBMBqNUg8taIVtnxEiIiIKDGFZM0JERESBg2GEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQphhEiIiKSFMMIERERSYphhIiIiCTFMEJERESSYhghIiIiSf1/PfGL2beB0DwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(AVERAGE_SIGMA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c9f269e-4637-4109-ba6b-96ea157c12a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74338237-aefb-4ac4-b3a6-11708ffe6baf",
   "metadata": {},
   "outputs": [],
   "source": [
    "python3 autocorrel.py 0& python3 autocorrel.py 1& python3 autocorrel.py 2& python3 autocorrel.py 3& python3 autocorrel.py 4& python3 autocorrel.py 5&"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "NERSC Python",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}