MCRG with Cluster MC - Working Version.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MCRG Code following Swendsen Description circa 1982, Using Wolff's Algorithm to combat critical slowing down"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division #safeguard against evil floor division\n",
    "\n",
    "import numpy as np\n",
    "from scipy import linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit\n",
    "%load_ext line_profiler\n",
    "\n",
    "np.set_printoptions(precision = 3,suppress=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D Ising Model, periodic B.C., initial Hamiltonian contains only nearest neighbor interactions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\mathcal{H} = -\\beta H_\\Omega = K\\sum_{\\langle i,j\\rangle} \\sigma_i \\sigma_j + h\\sum_{i = 0}^{N_s} \\sigma_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Block spin transform, scale factor = b "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assignBlockSpin(total):\n",
    "    '''Rule for assigning block spin value. Random tiebreaker'''\n",
    "    if total > 0:\n",
    "        s = 1;\n",
    "    elif total < 0:\n",
    "        s = -1;\n",
    "    else:\n",
    "        s = np.random.choice([-1,1])\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RGTransform(S,b):\n",
    "    '''Take a spin config S and produce renormalized block spin config that \n",
    "    groups b*b spins together into one block. Does not modify input S'''\n",
    "    L = S.shape[0];\n",
    "    assert L//b >= 2, \"Renormalized lattice will have linear dimension <=1\"\n",
    "    newS = np.empty([L//b, L//b],dtype=int)\n",
    "    for j in np.arange(L//b):\n",
    "        for i in np.arange(L//b):\n",
    "            block = S[(b*i):(b*i+b),(b*j):(b*j+b)]\n",
    "            total = np.sum(block)\n",
    "            newS[i,j] = assignBlockSpin(total);    \n",
    "    return newS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First 7 short range even couplings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N1Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are nearest neighbors, dir = (1,0)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small for NN interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N2Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are next-nearest neighbors, dir = (1,1)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small for next nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,(j+1)%L] + S[(i-1)%L,(j+1)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N3Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are third-nearest neighbors,dir = (2,0)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for third nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[i,(j+2)%L] + S[(i+2)%L,j])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N4Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are 4th-nearest neighbors, dir = (2,1)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for 4th nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,(j+2)%L] + S[(i+2)%L,(j+1)%L] \n",
    "                           +S[(i-1)%L,(j+2)%L] + S[(i-2)%L,(j+1)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N5Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are 5th-nearest neighbors,dir = (2,2)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for third nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+2)%L,(j+2)%L] + S[(i-2)%L,(j+2)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PlaqCoupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSlSk for i,j,k,l that form 1 plaquette,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small to form plaquettes\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i+1)%L,(j+1)%L]*S[(i+1)%L,j]\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SubPlaqCoupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSlSk for i,j,k,l that form 1 sublattice plaquette,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to form sublattice plaquettes\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[(i+1)%L,j]*S[i,(j+1)%L]*S[(i-1)%L,j]*S[i,(j-1)%L]\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AllEvenCoupling(S):\n",
    "    '''for spin field config S, \n",
    "    Integrate measurement of first 7 even correlation functions in one vector'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to fit first 7 even couplings\"\n",
    "    val = np.zeros(7,dtype = float);\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += [S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j]), #nearest neighbor (1,0)\n",
    "                    S[i,j]*(S[(i+1)%L,(j+1)%L] + S[(i-1)%L,(j+1)%L]), #next nearest neighbor (1,1)\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i+1)%L,(j+1)%L]*S[(i+1)%L,j], # plaquette\n",
    "                    S[i,j]*(S[i,(j+2)%L] + S[(i+2)%L,j]), #3rd nearest neighbor (2,0)\n",
    "                    S[i,j]*(S[(i+1)%L,(j+2)%L] + S[(i+2)%L,(j+1)%L] \n",
    "                            +S[(i-1)%L,(j+2)%L] + S[(i-2)%L,(j+1)%L]),#4th nearest neighbor (2,1)\n",
    "                    S[(i+1)%L,j]*S[i,(j+1)%L]*S[(i-1)%L,j]*S[i,(j-1)%L], # sublattice plaquette\n",
    "                    S[i,j]*(S[(i+2)%L,(j+2)%L] + S[(i-2)%L,(j+2)%L])]#5th nearest neighbor (2,2)\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First 4 short range odd couplings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Magnetization(S):\n",
    "    L = S.shape[0];\n",
    "    assert L >=1, \"Lattice linear dim < 1\"\n",
    "    val = np.sum(S);\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletPlaq(S):\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small to have 3 spins on a plaquette\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[(i+1)%L,j]*S[(i+1)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[(i-1)%L,j]*S[(i-1)%L,(j-1)%L];\n",
    "            val += S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-1)%L];\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletAngle(S):\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to have 3 spins at an angle\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[(i+1)%L,j]*S[(i+2)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+2)%L];\n",
    "            val += S[i,j]*S[(i-1)%L,j]*S[(i-2)%L,(j-1)%L];\n",
    "            val += S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-2)%L];\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletRow(S):\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small have 3 spins in a row\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,j]*S[(i+2)%L,j] \\\n",
    "                           + S[i,(j+1)%L]*S[i,(j+2)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AllOddCoupling(S):\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to fit first 4 odd couplings\"\n",
    "    val = np.zeros(4,dtype = float);\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += [0, #magnetization\n",
    "                    S[i,j]*S[(i+1)%L,j]*S[(i+1)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[(i-1)%L,j]*S[(i-1)%L,(j-1)%L]+\n",
    "                    S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-1)%L], #3 spin plaquette\n",
    "                    S[i,j]*S[(i+1)%L,j]*S[(i+2)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+2)%L]+\n",
    "                    S[i,j]*S[(i-1)%L,j]*S[(i-2)%L,(j-1)%L]+\n",
    "                    S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-2)%L], # 3 spin angle\n",
    "                    S[i,j]*(S[(i+1)%L,j]*S[(i+2)%L,j] + \n",
    "                            S[i,(j+1)%L]*S[i,(j+2)%L])] #3 spin row\n",
    "    val[0] = np.sum(S);\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering for Wolff Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NNBonds(p):\n",
    "    '''returns set of bonds that connect site p to its 4 nearest neighbors'''\n",
    "    i = p[0]; j = p[1];\n",
    "    nbrs = [(i,(j+1)%L),(i,(j-1)%L),((i+1)%L,j),((i-1)%L,j)]\n",
    "    bonds = set();\n",
    "    for n in nbrs:\n",
    "        bonds.add(frozenset({p,n}))\n",
    "    return bonds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def buildCluster(S):\n",
    "    ''''Build Wolff cluster starting from random site for spin configuration S'''\n",
    "    #random seed location\n",
    "    L = S.shape[0];\n",
    "    init = (np.random.choice(L),np.random.choice(L))\n",
    "    Si = S[init[0],init[1]]\n",
    "    \n",
    "    #cluster starts with 1 element\n",
    "    cluster = {init}\n",
    "    \n",
    "    #nearest neighbors make up the frontier\n",
    "    bonds = NNBonds(init)\n",
    "    \n",
    "    #set of points already considered for adding to cluster\n",
    "    checked = set();\n",
    "    \n",
    "    #while the set of fresh bonds is nonempty, do...\n",
    "    while (len(bonds) > 0):\n",
    "        if len(cluster) == Ns:\n",
    "            break;\n",
    "        #take out one bond in fresh bond set\n",
    "        #frozenset ijbond represent unordered edge (i,j)\n",
    "        ijbond = bonds.pop() \n",
    "        #add to list of bonds that have been checked\n",
    "        checked.add(ijbond)\n",
    "        #pick out element j from edge (i,j)\n",
    "        jwrap = ijbond - cluster \n",
    "        #both i and j may already be in cluster, in this case skip to next iteration\n",
    "        if len(jwrap) == 0:\n",
    "            continue;\n",
    "        #otherwise, only i in cluster already, we are left with j\n",
    "        j = set(jwrap).pop() #convert to usable form\n",
    "        Sj = S[j[0],j[1]]\n",
    "        #if parallel to seed spin, activate bond with probability Pij\n",
    "        if Sj == Si:\n",
    "            r = np.random.random()\n",
    "            if r < Pij:\n",
    "                #add j to cluster, add nearest bonds of j to the fresh bond list\n",
    "                #also remove bonds already considered from the fresh bond list\n",
    "                cluster.add(j) \n",
    "                bonds |= NNBonds(j) \n",
    "                bonds -= checked \n",
    "    return cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrated MC + RG simulation function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Energy(S):\n",
    "    '''Brute force Find energy of spin configuration S for sanity check'''\n",
    "    L = S.shape[0];\n",
    "    E = 0;\n",
    "    for i in np.arange(L):\n",
    "        for j in np.arange(L):\n",
    "            E += K*S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j])\n",
    "    E += h*np.sum(S)\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RunMCRG(K,h):\n",
    "    '''Run MCRG simulation to find y_t, y_h exponent, \n",
    "    keeping Nc_even and Nc_odd coupling terms'''\n",
    "    print('running MCRG for linear size',L,'lattice.')\n",
    "    print('Setting K =', K, \" and h = \",h)\n",
    "\n",
    "    #measurement accumulators for y_t\n",
    "    evenK = np.zeros(7,dtype = float)\n",
    "    evenK_new = np.zeros(7,dtype = float)\n",
    "    mix_11 = np.zeros((7,7),dtype=float)\n",
    "    mix_01 = np.zeros((7,7),dtype=float)\n",
    "    #measurement accumulators for y_h\n",
    "    oddK = np.zeros(4,dtype = float)\n",
    "    oddK_new = np.zeros(4,dtype = float)\n",
    "    mix_11_odd = np.zeros((4,4),dtype=float)\n",
    "    mix_01_odd = np.zeros((4,4),dtype=float)\n",
    "    \n",
    "    #Run simulation\n",
    "    k = 0;\n",
    "    for n in np.arange(nmeas+nwarm):\n",
    "        # Every MC n-loop, build a Wolff cluster and flip it\n",
    "        # Result: A S-field config drawn with probability propto Boltzmann weight\n",
    "        cluster = buildCluster(S)\n",
    "        for p in cluster:\n",
    "            S[p[0],p[1]] = -S[p[0],p[1]]\n",
    "            \n",
    "        #Sanity checks\n",
    "        if n % interval == 0:\n",
    "            #energy[k] = Energy(S);\n",
    "            clustersize[k] = len(cluster);\n",
    "            k = k+1   \n",
    "            \n",
    "        # take measurements every (interval) steps if finished warmup\n",
    "        if n % interval == 0 and n >= nwarm:\n",
    "            if n % 100 == 0:\n",
    "                #monitor progress\n",
    "                print(\"iteration\",n)\n",
    "            if layer == 1:\n",
    "                S1 = RGTransform(S,b);\n",
    "                evenK += AllEvenCoupling(S)\n",
    "                evenK_new += AllEvenCoupling(S1)\n",
    "                oddK += AllOddCoupling(S)\n",
    "                oddK_new += AllOddCoupling(S1)\n",
    "                #A*B = C, B is unknown, A is symmetric\n",
    "                #Problem: C is not symmetric!?\n",
    "                mix_11 += np.outer(AllEvenCoupling(S1),AllEvenCoupling(S1))\n",
    "                mix_01 += np.outer(AllEvenCoupling(S1),AllEvenCoupling(S))\n",
    "                mix_11_odd += np.outer(AllOddCoupling(S1),AllOddCoupling(S1))\n",
    "                mix_01_odd += np.outer(AllOddCoupling(S1),AllOddCoupling(S))\n",
    "\n",
    "            elif layer == 2:\n",
    "                assert L > 8, \"Lattice too small\"\n",
    "                S1 = RGTransform(S,b);\n",
    "                S2 = RGTransform(S1,b);\n",
    "                evenK += AllEvenCoupling(S1)\n",
    "                evenK_new += AllEvenCoupling(S2)\n",
    "                oddK += AllOddCoupling(S1)\n",
    "                oddK_new += AllOddCoupling(S2)\n",
    "                #A*B = C, B is unknown, A is symmetric\n",
    "                #Problem: C is not symmetric!?\n",
    "                mix_11 += np.outer(AllEvenCoupling(S2),AllEvenCoupling(S2))\n",
    "                mix_01 += np.outer(AllEvenCoupling(S2),AllEvenCoupling(S1))\n",
    "                mix_11_odd += np.outer(AllOddCoupling(S2),AllOddCoupling(S2))\n",
    "                mix_01_odd += np.outer(AllOddCoupling(S2),AllOddCoupling(S1))\n",
    "\n",
    "            elif layer == 3:\n",
    "                assert L > 16, \"Lattice too small\"\n",
    "                S1 = RGTransform(S,b);\n",
    "                S2 = RGTransform(S1,b);\n",
    "                S3 = RGTransform(S2,b);\n",
    "                evenK += AllEvenCoupling(S2)\n",
    "                evenK_new += AllEvenCoupling(S3)\n",
    "                oddK += AllOddCoupling(S2)\n",
    "                oddK_new += AllOddCoupling(S3)\n",
    "                #A*B = C, B is unknown, A is symmetric\n",
    "                #Problem: C is not symmetric!?\n",
    "                mix_11 += np.outer(AllEvenCoupling(S3),AllEvenCoupling(S3))\n",
    "                mix_01 += np.outer(AllEvenCoupling(S3),AllEvenCoupling(S2))\n",
    "                mix_11_odd += np.outer(AllOddCoupling(S3),AllOddCoupling(S3))\n",
    "                mix_01_odd += np.outer(AllOddCoupling(S3),AllOddCoupling(S2))\n",
    "\n",
    "            elif layer == 4:\n",
    "                assert L > 32, \"Lattice too small\"\n",
    "                S1 = RGTransform(S,b);\n",
    "                S2 = RGTransform(S1,b);\n",
    "                S3 = RGTransform(S2,b);\n",
    "                S4 = RGTransform(S3,b);\n",
    "                evenK += AllEvenCoupling(S3)\n",
    "                evenK_new += AllEvenCoupling(S4)\n",
    "                oddK += AllOddCoupling(S3)\n",
    "                oddK_new += AllOddCoupling(S4)\n",
    "                #A*B = C, B is unknown, A is symmetric\n",
    "                #Problem: C is not symmetric!?\n",
    "                mix_11 += np.outer(AllEvenCoupling(S4),AllEvenCoupling(S4))\n",
    "                mix_01 += np.outer(AllEvenCoupling(S4),AllEvenCoupling(S3))\n",
    "                mix_11_odd += np.outer(AllOddCoupling(S4),AllOddCoupling(S4))\n",
    "                mix_01_odd += np.outer(AllOddCoupling(S4),AllOddCoupling(S3))\n",
    "            else:\n",
    "                raise ValueError(\"too many RG iterations\")\n",
    "                       \n",
    "    #Results\n",
    "    evenK /= ndata; evenK_new /= ndata; mix_11 /= ndata; mix_01 /= ndata;\n",
    "    oddK /= ndata; oddK_new /= ndata; mix_11_odd /= ndata; mix_01_odd /= ndata;\n",
    "\n",
    "    print('evenK = ',evenK)\n",
    "    print('evenK_new = ', evenK_new)\n",
    "    print('mix_11 = ', mix_11)\n",
    "    print('subtract ',np.outer(evenK_new,evenK_new))\n",
    "    print('mix_01 = ', mix_01)\n",
    "    print('subtract ',np.outer(evenK_new,evenK))\n",
    "    MatA_even = mix_11-np.outer(evenK_new,evenK_new)\n",
    "    MatC_even = mix_01-np.outer(evenK_new,evenK)\n",
    "    print('MatA_even (lhs) = ',MatA_even)\n",
    "    print('MatC_even (rhs) = ',MatC_even)\n",
    "    print('\\n')\n",
    "    \n",
    "    print('oddK = ',oddK)\n",
    "    print('oddK_new = ', oddK_new)\n",
    "    print('mix_11 = ', mix_11_odd)\n",
    "    print('subtract ', np.outer(oddK_new,oddK_new))\n",
    "    print('mix_01 = ', mix_01_odd)\n",
    "    print('subtract ', np.outer(oddK_new,oddK))\n",
    "    #TODO: cancellation bad, how to avoid?\n",
    "    MatA_odd = mix_11_odd-np.outer(oddK_new,oddK_new)\n",
    "    MatC_odd = mix_01_odd-np.outer(oddK_new,oddK)\n",
    "    print('MatA_odd (lhs) = ',MatA_odd)\n",
    "    print('MatC_odd (rhs) = ',MatC_odd)\n",
    "    print('\\n')\n",
    "\n",
    "    \n",
    "    return MatA_even, MatC_even, MatA_odd, MatC_odd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getExponent(MatA,MatC,Nc):\n",
    "    '''Get thermal or magnetic exponent based on output of RunMCRG\n",
    "    and desired number of coupling constants to consider Nc'''\n",
    "    LinRGMat = la.solve(MatA[0:Nc,0:Nc],MatC[0:Nc,0:Nc])\n",
    "    lmbd = la.eigvals(LinRGMat);\n",
    "    amplitude = np.absolute(lmbd)\n",
    "    imax = np.argmax(amplitude)\n",
    "    y = np.log(lmbd[imax])/np.log(b)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test for specific set of input parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['64', '1000', '3000', '10', '4', 'testmultipleRG']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x239b3d24128>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAD8CAYAAADqmhgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+0XlV55z/fhh8pCA0hoDEJElZThY4CmgZdzPJHABNtB5wZsGDrBBesLGZJa8exBXQWtnTsip01orOGqrcQwWoJiFUyGI2RH3V1WSA3CoFAgZgy5Bo0koDaIoHc+8wf51zm3Pe+73n3u999zvvr+bDOuu85Z5+9n/cHT/Z+9vNDZobjOM4o8Cu9FsBxHKcuXOE5jjMyuMJzHGdkcIXnOM7I4ArPcZyRwRWe4zgjgys8x3EqQ9J6SXslPdziviT9L0k7JW2X9MbCvTWSnsiPNSnkcYXnOE6V3AisLrn/LmBZfqwFPgsgaT7wceAMYAXwcUnHdCtMVwpP0mpJj+Xa+cpuhXEcZ7gws+8C+0uanAd80TLuBeZJWgisAraY2X4zexbYQrniDOKQ2AclzQGuA84BJoCtkjaa2SOtnjlMh9tcjmzb92+84fmW9x7ffkSQfI19xD5XJ0UZi3I0yl52r1W70HFjnwv9zlJ8vlXIG3svpo8UtOr/Bf6VF+2Auul71TuOtH37J4Pabtt+YAfwQuHSmJmNdTDcImB34Xwiv9bqeldEKzyyaeZOM9sFIGkDmbZuqfDmciRn6Ky2HW/e/EDLe6tefVqQcI19xD5XJ0UZi3I0yl52r1W70HFjnwv9zlJ8vlXIG3svpo8UtOr/Pruz67737Z/k/s0nBLWds/CJF8xseRfDNVPOVnK9K7pZ0gZpYElrJY1LGn+JA10M5zhOHRgwFfhfAiaAJYXzxcCekutd0c0ML0gD59PbMYCjNf/l+5v3RM7A9nQ/+yvrL0aOWTOOSBnL+gx5powyGWP7CO0/lNj3HPO9p3ovqX+3MZ9b8bkVq7o3GRjGSxa2pE3ARuDyfIV4BvAzM3ta0mbgLwobFe8Erup2sG4UXiUa2HGc3pNo9oakm4G3AwskTZDtvB4KYGafAzYB7wZ2As8DH8jv7Zf058DWvKtrzKxs8yOIbhTeVmCZpKXAj4ALgfd1K5DjOL3FMCYTpY0zs4va3Dfggy3urQfWJxEkJ1rhmdlBSZcDm4E5wHoz25FMMsdxesZU9/sDfUk3MzzMbBPZlDSI33jD81G7dK3sJJ3YVlq17cT+08rG1q3dpdN2MWNXYYdK8RmktmnGEjpWqLyxxH4eKcaexoBJV3iO44wKPsNzHGckMOClIS39oDprWhyt+TbteJxyCj5N2RKr1b3YpVKsm0Gd7iapqWI5l2KpGvNd1Pm5dUK3v8cVq3Yz/uALXUVavOHUQ+0bmxYEtT1h8Y+3del4XCs+w3McZyYGk8M5wXOF5zjOTLJIi+HEFZ7jOA2IyaaBVINPrQqvzC0ltX2ozEaTwp2ibKzQe6H2vRSfTdVuHWXjpXadSRFaVkaddsV2bTvt83HbF9RXGdmmhSs8x3FGgMwPzxWe4zgjwtSQzvB65pZSRuiSJUVWkjKq8HpPkdGlrL/QiJJQ+VMv/8tI7bJT1n/sZ1pFBE8rOlm6p3RLOeUNh9mX7nhVUNs3vWa3u6U4jjO4GGJySMvduMJzHGcWw7qk7RuFF7PESkHVS4/Ge6mXR7HL/9Cd5dQRFCkSLVQdwN9ITARPmVyx8teVTMEQL9qcyvrvJX2j8BzH6Q8yx2Nf0jqOMyK4W4rjOCOBmZg0n+ElpRNbSKi9ouqknDGREJ2QOtFp8V4Ku2gntqcYOcr6rNo9pl/6r9s+2YqphDM8SauBz5BlRr/ezNY13L8WeEd+egRwvJnNy+9NAg/l954ys3O7kcVneI7jzCDbtEijGiTNAa4DziEr/LVV0kYze7l+tZn9l0L7PwBOL3TxSzNLtkMznPNWx3Gimd60CDkCWAHsNLNdZvYisAE4r6T9RcDN3b+L5vQseUBsLYkisUHYsUu4Vv2nSvKZYmmcwi2lVbuq6+2mSNpadzRIisD/0GdCXIlSJA8AmEznh7cI2F04nyCrPzsLSa8BlgJ3FS7PlTQOHATWmdnXuxHGl7SO48ygw0iLBblCmmbMzMYK5800Z6t41guB28xmVAE/wcz2SDoJuEvSQ2b2w1DhGnGF5zjOLKbCd2mfaRNLOwEsKZwvBva0aHshDTVqzWxP/neXpHvI7HvRCs9teI7jzCBLD/UrQUcAW4FlkpZKOoxMqW1sbCTptcAxwD8Wrh0j6fD89QLgTOCRxmc7oS9neFWElsW4fKSow5rC5lOFzarOUL4Ubi8p+ihrl9rOWLUtscokqIZ4KVFomZkdlHQ5sJnMLWW9me2QdA0wbmbTyu8iYIPNTN90MvB5SVNkk7N1xd3dGPpS4TmO0zvMSOp4bGabgE0N165uOP/TJs99D3h9MkEIWNJKWi9pr6SHC9fmS9oi6Yn87zEphXIcp5eIqcBj0AiZ4d0I/G/gi4VrVwJ3mtk6SVfm51e06+jx7Ue8PPVOkSQzNktJGandNWKXxVVnhalzaR0bxVB1/eB+Ifa3VBVG2hleP9H2XZnZd4H9DZfPA27KX98EvCexXI7j9JCEmxZ9RawN75Vm9jSAmT0t6fhWDSWtBdYCzOWIyOEcx6kLQ54ANJbcCXEMspoWMX2krnWQor+ygPiysWKWiLH9pyY2iUGdS/cUlJkhUpXkDOkjZvm/YtXzwc+0IivTOJz7mbHv6ieSFuazu4XA3pRCOY7TS4a3EHfsInwjsCZ/vQa4PY04juP0GiOLtAg5Bo22MzxJNwNvJ4uZmwA+DqwDbpV0CfAUcEGVQjqOUy/DOsNrq/DM7KIWt9oXmG2gmC2ljBS2s7I+Y21KMba/WK/9qm1bdUZQlBHjetLJ91dlMaHG/sveS9Wfd0rbrZkGcvYWwnBaJh3HiSbbtPCqZY7jjARe0yI5nbhrtFqOVp1MsozQ5UsnQd4pEp3GvJ/Y2rap3YBiqbKObmOfob/bKupWhPzmUiQAzTYtRtSG5zjO6DGIURQhuMJzHGcGHmnhOM5IEVigZ+CoVeGFZkspI3UiyBTb+XUX8anaDpg6DKpqN6CyZ6rOCpPaNtyq73YkDS0zeGnKFZ7jOCNAtqR1hec4zogwspEWKQmNtCgjdtmQukZECjeG0LGrqL+buo5s1VEAVS85UyQwrdq1pa7MMsPsljKc81bHcbpASZMHSFot6TFJO/MM6Y33L5b0U0kP5MelhXtr8lIST0ha0/hsp/iS1nGcWaSqVyFpDnAdcA5ZjdqtkjY2qT52i5ld3vDsfLJkJcvJJp7b8mefjZWnZ7u0jdSZJDJ1cs0UtR7atY0ZL0UpxipLZDb2X/Zcahlja2t00083z7QjaaSFwUtTyWJpVwA7zWwXgKQNZCUiQsotrgK2mNn+/NktwGrg5lhhfEnrOM4Mph2PQw6ytHHjhWNtQ3eLgN2F84n8WiP/UdJ2SbdJWtLhs8H4ktZxnFl0sKR9xsyWl9xv1lFjqYf/A9xsZgckXUZWGGxl4LMd4TM8x3FmML1LGzjDa8cEsKRwvhjYM2M8s31mdiA//WvgTaHPdspAz/BSJIJs91xIu9iMKFXUcq1S/hQ1ZWNtmv1Y7AfSJIUN/U3U+XkkdDzeCiyTtBT4EXAh8L5ig+n6OPnpucCj+evNwF9IOiY/fydwVTfCDLTCcxwnPWbiYCKFZ2YHJV1OprzmAOvNbIeka4BxM9sI/KGkc4GDZDWwL86f3S/pz8mUJsA10xsYsbjCcxxnFikdj81sE7Cp4drVhddX0WLmZmbrgfWpZBkIhVdlrdVYqpYp1kWlyuQBVUcnpEimEFonpJOxUiRVLSPGhFBlFMYwR1oMhMJzHKdeXOE5jjMSeAJQx3FGilShZf1G32RLCbWvpAhBi7WPpa6bWmdmjxQ2n1ibUpWhgY2kCAeMlaXOpLOt3meqBKAHPQGo4zijwrAuaduqcUlLJN0t6VFJOyR9KL8+X9KWPG3LloJzoOM4A0yHsbQDhczKQ9MkLQQWmtn3JR0FbAPeQ+YcuN/M1uU5ro4xsyvK+jpa8+0MnQVUkx0lps8U9SiqyJaSoi5GimwpoWP3MhtLaheh2JrJ/SD/ilW7GX/wha400VGvfZW96a9+P6jt35/9P7e1iaXtK9rO8MzsaTP7fv76F2RhH4vIUrzclDe7iUwJOo4zBEyhoGPQ6MiGJ+lE4HTgPuCV0/FvZva0pONbPLMWWAswlyO6kdVxnBowG14bXrDCk/QK4KvAH5nZz6WwD8TMxoAxyJa0MUI6jlMnYnKUd2klHUqm7L5sZn+XX/7JdJaD3M63tyohY+wkVdiNUmf/TZGlJDQkKtbVp4yq+09R3CZ1YaQUdWlTuF216iNFxmPIEggMIyG7tAJuAB41s08Vbm0EpotqrAFuTy+e4zh1kzgfXl8RMsM7E3g/8JCk6X9WPgqsA26VdAnwFHBBNSI6jlMrltnxhpG2Cs/M/oHmqZYBzupksGKkRSdJEKtMdlhFEtEUS7gy6qyTWtZf6lqrsckv63RHSmFCSFEvuGoGcQc2BI+0cBxnBjbqmxaO44wWI7ukrYoUdUHLlpyhS5YUu7QpdlE7IcVytIzUiTdTLFX7cVlc1meKhA+d/HZSJg+AEd6ldRxntDDLFF7IEYKk1ZIek7QzD0NtvP9hSY/kdWnvlPSawr1JSQ/kx8Zu35svaR3HmUUqlxNJc4DrgHPIyi5ulbTRzB4pNPsBsNzMnpf0n4G/BH43v/dLM0u2a+kzPMdxZmEWdgSwAthpZrvM7EVgA1kcfmEsu9vMptfi95LVn62EgZvhpahxWjV1JixNbR9r12eRFC4UVSdEjbWJhfZfZQ3fsrGq/D0bYip8l3aBpPHC+VgeTjrNImB34XwCOKOkv0uAbxbO5+b9HwTWmdnXQwVrxsApPMdxqqeDTdpn2qSHarY2btq9pN8HlgNvK1w+wcz2SDoJuEvSQ2b2w3DxZuIKz3GcmVjSXdoJYEnhfDGwp7GRpLOBjwFvM7MDL4titif/u0vSPWTZmgZP4cUGV4f2Eds2Noi81b2qk0JWGYXSbqwUSVBD+w99roplcYyMVSQPCOk/VfKATqZ4bdgKLJO0FPgRcCHwvmIDSacDnwdWm9newvVjgOfN7ICkBWRhrn/ZjTA+w3McZxapZnhmdlDS5cBmYA6w3sx2SLoGGDezjcD/AF4BfCVPO/eUmZ0LnAx8XtIU2Qbruobd3Y5xhec4zgwMmJpK53hsZpuATQ3Xri68PrvFc98DXp9MEFzhOY7TiAFDGmnRN6FlodkkUidZjLUDpghFCn0uRWaZWNeWMup0v0k9Viypk5l202eVeCyt4zijgys8x3FGg/A42UGjb9xS6qyT2olcIX1UvURJsYRL4TZSdWLT0D46WeIXic24UnVESSghcqTKluIzPMdxRgMDS7hL20+4wnMcpwmu8Lrm8e1HvDwtT5EANIXHehUlEENJUeqxkdT1KFL0UYWJIjVVJAcNJeUudh9GWvQVPsNzHGc2rvAcxxkJ3PHYcZxRwh2PE1CsSxtLFVEYof2nGCuF7S/muU5cLWJquaaI+Ejx3aZwEao7sWxMsaLK8V1ax3FGBQ3pDK9tHmdJcyXdL+lBSTsk/Vl+famk+yQ9IekWSYdVL67jOJVjHRwDRsgM7wCw0sz+RdKhwD9I+ibwYeBaM9sg6XNkueg/W9ZR0S2lkaq3/nsVoF1FfdLUdWlT1HkNvVenqaHxuRTJCUJJ8b03EiJXmkgLDe2mRdsZnmX8S356aH4YsBK4Lb9+E/CeSiR0HKd+hnSGF1SaSNIcSQ8Ae4EtZDnlnzOzg3mTCbLqRM2eXStpXNL4Sxxo1sRxnH5jKvAYMIIUnplN5sVwF5PVmTy5WbMWz46Z2XIzW34oh8dL6jhOPUz74YUcA0ZHu7Rm9lxeOejNwDxJh+SzvKaViBpJ4ZZSZr+KyRzSic2krowX7drFjFV1xpVQOvncUha36aS/2GwsoYTaZKvI2hJKyl1aSauBz5DVtLjezNY13D8c+CLwJmAf8Ltm9mR+7yqy/YFJ4A/NbHM3soTs0h4naV7++leBs4FHgbuB8/Nma4DbuxHEcZw+IpENT9Ic4DrgXcApwEWSTmlodgnwrJn9OnAt8Mn82VPIqpz9JrAa+Ku8v2hClrQLgbslbScrubbFzO4ArgA+LGkncCxwQzeCOI4zlKwAdprZLjN7EdgAnNfQ5jyyjU/INkLPUla+7Dxgg5kdMLN/Bnbm/UXTdklrZtvJit82Xt/V7eCtqDuJY4xcVWQASb0kr5rULjahVOHaUmWC2MbnUtQXqfp772BJu0DSeOF8zMzGCueLgN2F8wngjIY+Xm6Tl3X8GdkkahFwb8OzTTdHQ/FIC8dxZmJ0Elr2jJktL7nfrKNGddqqTcizHRG0S+s4zoiRzg9vAlhSOG+2wflyG0mHAL8G7A98tiMGYoaXosRijOd/qByxpIgKSJFENLZkZord3TprbdSZyDOFSaXWZAENJNyl3Qosk7QU+BHZJsT7GtpsJNv4/EeyjdC7zMwkbQT+VtKngFcDy4D7uxFmIBSe4zg1k0jh5Ta5y4HNZG4p681sh6RrgHEz20i24fk3+QbofjKlSN7uVuAR4CDwQTOb7EYeV3iO48wmoR+emW0CNjVcu7rw+gXgghbPfgL4RCpZXOE5jjMD2fCmh+pZEZ9YYm1PdRaEibVLpbZFxbqDhHr+p6bMTSem9my7tqHtYqJj6rYNJy/i4wlAHccZFXyG5zjO6OAKLy0pliWxy4YUlC1Dqq6ZkTpyIVSORlIHs9cZNVJnIH4nY9eZyKElbsNzHGekcIXnOM6ooAFM7hmCh5Y5jjMy9M0ML0WyyqozXrRql8rGFmOXia0Hm9oGVGeBnLL+Q8dO9RursoZvWZ+Vh535ktZxnJHANy0cxxkpXOF1T1lNi9S1AqrwsG/Vtu5EpGWk8PxPUWMhdY3d1G5AKZat7cYOJYULT/LflSs8x3FGATG8u7Su8BzHmYnb8NIQmzyg1e5aisSSsbuosTvCMTttsUusfql3kbqeQyOhy93YKIaql92h1FrLxBWe4zgjgys8x3FGBV/SOo4zOrjCS0snNojU3vKxtpDQrC0p6urG2qVSJAANdeWo2vO/ygiEWBebFFRR0zgpVs8uraT5wC3AicCTwHvN7NmGNqcBnwWOBiaBT5jZLfm9G4G3AT/Lm19sZqX/Q3ssreM4s0lXprGMK4E7zWwZcGd+3sjzwH8ys98EVgOfljSvcP+Pzey0/Gg7ewlWeJLmSPqBpDvy86WS7pP0hKRbJB0W2pfjOP3NdF2LdkeXnAfclL++CXhPYwMze9zMnshf7wH2AsfFDtjJkvZDwKNkU0uATwLXmtkGSZ8DLiGberakLNIi1o0k5Jmy51J53MeQ+j2XkTopabvnYtyHqna1iKnj0e65mKQXsb+rWuvUhiuzBZLGC+djZjYW+OwrzexpADN7WtLxZY0lrQAOA35YuPwJSVeTzxDN7EBZH0EzPEmLgd8Grs/PBawEbsubNNXOjuMMIKHL2UwpPmNmywvHDGUn6TuSHm5ynNeJSJIWAn8DfMDMpi2MVwGvA34LmA9c0a6f0Bnep4E/AY7Kz48FnjOzg/n5BLCohaBrgbUAJyzyTWHH6XdEOrcUMzu75TjSTyQtzGd3C8mWq83aHQ18A/hvZnZvoe+n85cHJH0B+Eg7edrO8CT9DrDXzLYVLzdp2vQjMrOxae1/3LFz2g3nOE4fUJMNbyOwJn+9Brh9lhzZ3sDXgC+a2Vca7i3M/4pshflwuwFDplxnAudKejcwl8yG92lgnqRD8lneYmBPu45CQ8tSJFmMtUuFjNtNH6kTgIaSol5rioJBKcK26q6xm9qGmiLscUgSgK4DbpV0CfAUcAGApOXAZWZ2KfBe4K3AsZIuzp+bdj/5sqTjyCZgDwCXtRuwrcIzs6vI1spIejvwETP7PUlfAc4HNtBCOzuOM6DUoPDMbB9wVpPr48Cl+esvAV9q8fzKTsfsxg/vCuDDknaS2fRu6KIvx3H6hcDl7CCGn3W0i2Bm9wD35K93ASvSizSbFG4pRVJkWSmjiiwuoXKELutjIj7KiF3u9kuN1hS/q7L+Yl1iWvXZqr8Vq54Pkq8tA6jMQvBtU8dxZuEJQB3HGRkGcbkaQt8kDwjdhYtNAJo68L/smdilaZWlE2NrWvSS1Ev8FDvOqZe+ZX3G/BYft33B47YkTZxsX+IzPMdxZuMKz3GcUSBlpEW/4QrPcZxZaGo4NV7P6tLG2rZiE4DG2MeqcCFJ7UIRaj9MIUeKSItuxuv2mZjMJo2kcEtJUayoUtyG5zjOKOFLWsdxRgdXeP1P1VEYrZ6LjVRI4WpRRmpXi07GCnUHSTFWTNt+Wp6HylgnPsNzHGd0cIXnOM5IUFPVsl7gCs9xnBm4H14iQhOAVpElIybhZejYnYSWVd1/KFWHWcXYxPrRzSUVKb731G5Gpdhwajyf4TmOMwuf4TmOMxq443H1xGZBadVH2XN1unxUXcMihTtIbA2RMlL30cuar6mX2lXUI55+LlUC0Do2LSTNB24BTgSeBN5rZs82aTcJPJSfPmVm5+bXl5KVmJgPfB94v5m9WDZmNyneHccZUjQVdnTJlcCdZraMvJB2i3a/NLPT8uPcwvVPAtfmzz8LXNJuQFd4juPMxMg2LUKO7jgPuCl/fRNZqcUg8tKMK4HbOnm+Z8kDOiHUEz10ZzNmrDJSlEBs7Cf1MipVKclW7VKUHix7rkgV9SLK+qvazBFThrRqOti0WCBpvHA+ZmZjgc++crqYdl6M+/gW7ebmYxwE1pnZ18kKhz2Xl4kFmAAWtRuwb2x4juP0EeEK7xkzW97qpqTvAK9qcutjHUhzgpntkXQScJekh4CfN2nXVmpXeI7jzCCl47GZnd1yHOknkhbms7uFwN4WfezJ/+6SdA9wOvBVYJ6kQ/JZ3mJgTzt53IbnOM5MzNBU2NElG4E1+es1wO2NDSQdI+nw/PUC4EzgETMz4G7g/LLnG+nZDK+TBIkxNpSy/qssltPYfyf3YpJ3piog1IpYl48q3YxS2LZiM66E2g+rdkeqPNtLPX5464BbJV0CPAVcACBpOXCZmV0KnAx8XtIU2QRtnZk9kj9/BbBB0n8HfgDc0G5AX9I6jjOLOiItzGwfcFaT6+PApfnr7wGvb/H8LmBFJ2MGKTxJTwK/ACaBg2a2PNRp0HGcAcMAr2nBO8zsmcL5tNPgOklX5udXxApS9bIk9LkUiRpjl88xURKxS6wUtR/KiOkjtA5JFckUUjyXuo5ujBxJ6tLC0IaWdbNpEe006DhOfyMLOwaNUIVnwLclbZO0Nr82w2kQaOo0KGmtpHFJ4z/dN9m9xI7jVE5Nu7S1E7qkPTN3/Dse2CLpn0IHyL2uxwCWnzp38D4hxxk1Rj1bSsHxb6+kr5HtjAQ5DRYpSwAa41JSdYaLKjJ0pAiDqrMoTpXhbp3QL0lPG6kzK0yMHDFkjsfDqfHaLmklHSnpqOnXwDuBhwlwGnQcZ0CZCjwGjJAZ3iuBr2XJCTgE+Fsz+5akrTRxGnQcZ/AZ1hmerMY3tvzUuXb/5iVd9ZEi40WdERkpojA6WX6lcAdpNVYV9TmqJkV93CqjdBr775YVq3Yz/uAL6qaPo49abL+1/INBbe+656PbypIH9BseaeE4TgODuQMbgis8x3FmM6RLWld4juPMxAtxpyG0Lm0sKTJ0pCgcUzZWjP2wCleFFO4xKbLYhMoU+3mkyN6cugBUGbHvM31omc/wHMcZFYZT37nCcxxnNpoazjVt3xTxiU142YpYl4y6Eze2IoULRQo5QvuONRNUbV5IEQkRQ9UFeCot6GMMpFNxCD7DcxxnBsKG1vHYFZ7jOLNxhdc9xV3a2LqxVScM6JeogFD6NaKkFZXXYogkdbKGKnbIU9R9CaYGhReSNV3SO4BrC5deB1xoZl+XdCPwNuBn+b2Lzaz0B+VVyxzHmcm0Da/65AHTWdOXAXfm5zNFMbvbzE4zs9OAlcDzwLcLTf54+n47ZQeu8BzHaYKmpoKOLuk0a/r5wDfN7PnYAV3hOY7TgGVL2pCjO4Kyphe4ELi54donJG2XdO10/doy+nLTItbelMK+lzoRaehYZf2UtYvNHlNGisI0qW2hKTK6pC7UE9suBa3e54pV0ZOf/4/RiTJbIGm8cD6WZzkHQNJ3gFc1ee5jnYiUJxl+PbC5cPkq4MfAYWRZ1a8Arinrpy8VnuM4PSZ8tfpMWXooMzu71T1JnWRNfy/wNTN7qdD30/nLA5K+AHyknbC+pHUcZxYyCzq6pJOs6RfRsJzNlSTKshO/hywTeyk9m+Gl2LbvpP9Wy42ql2Kx9S7Klq0xclXtDhJbHzemj27kakWKWsKh7fplWVxKPX5462iSNV3ScuAyM7s0Pz8RWAL8fcPzX5Z0HFkZjgeAy9oN6Etax3FmYgaT1ceWmdk+4Kwm18eBSwvnTwKLmrRb2emYrvAcx5mNR1o4jjMyuMJLS4rQmE5sPinsXilInQwzdKwUtXNjCxKF9tFIiveZYqzUoWVlbauondsxBnhNC8dxRgMDG878UK7wHMeZiVHLpkUv6Bu3lJit+hQuJXUnAA0dK3TcWNNAzPssy3CTOsolVKbGeykykZS5NJX1WfZMilomIfe8pkU5PsNzHGc2Q6rwgiItJM2TdJukf5L0qKS3SJovaYukJ/K/x1QtrOM4dVBb8oDaCZ3hfQb4lpmdL+kw4Ajgo2S5rNZJupIsl9UVZZ0Ua1p0smMWs2woa5uiXYrdtBRJB2J3DVPv+KVYWleRCCGU2LFT193o1hSTLHnAkBbxaTvDk3Q08FbgBgAze9HMnqPzXFaO4wwKQzrDC1nSngT8FPiCpB9Iul7SkQTmspK0VtK4pPGf7ptMJrjjOFWRh5aFHANGiMI7BHgj8FkzOx34V5qkYm6FmY2Z2XIzW37csXMixXQcpzYMzKaCjkEjxIY3AUyY2X1fqtIIAAAFgUlEQVT5+W1kCq+TXFZtSW1jq7RuZxs5+sUeFGs3isnakiKKIdY2V/V3XTZWCpeblDbqZG4pQxpp0XaGZ2Y/BnZLem1+6SzgETrLZeU4ziAxpDa80F3aPyDLPXUYsAv4AJmynJXLynGcAcdsaHdpZTVq6eWnzrX7Ny9p267qJUrVNRBSB/53IkcKN48USShjkgf0MjFmL+sdp3VL2c34gy8oVk6AX5uzwN5y5L8Larv5FzduK0vx3m94pIXjOA0YNjmcHhWu8BzHmYmnh3IcZ6QYQJeTEGpVeI9vP+Jl20MndpGU2SQaX4eO267/1FQdZhUbNhfaLtTVoowUyTZj5Ii1waZOQFtGlb8/A6yGGZ6kC4A/BU4GVuS1LJq1W00W3joHuN7M1uXXlwIbgPnA94H3m9mLZWN6mUbHcWZieQLQkKM7Hgb+A/DdVg0kzQGuA94FnAJcJOmU/PYngWvNbBnwLHBJuwFd4TmOMwubnAw6uhrD7FEze6xNsxXATjPblc/eNgDn5bVoV5IFQkBgPH+tbimSfgr8X2AB8ExtAzenH2QAl6MRl2MmncrxGjM7rpsBJX0rHzeEucALhfMxMxvrcLx7gI80W9JKOh9YXahR+37gDLKl8L1m9uv59SXAN83s35SNVasNb/qLkDTea9+dfpDB5XA5+lEOM1udqi9J3wFe1eTWx8wsJDqrmU+hlVwvxXdpHcepDDM7u8suJoBitMJiYA/ZrHeepEPM7GDheiluw3Mcp5/ZCiyTtDQPbb0Q2GiZLe5u4Py8XVA8f68UXkdr/IroBxnA5WjE5ZhJv8iRHEn/XtIE8BbgG5I259dfLWkTQD57uxzYDDwK3GpmO/IurgA+LGkncCx5kuLSMevctHAcx+klvqR1HGdkcIXnOM7IUKvCk7Ra0mOSduaVzuoad72kvZIeLlyrvcykpCWS7s5LXe6Q9KFeyCJprqT7JT2Yy/Fn+fWlku7L5bglNxJXjqQ5eb2UO3olh6QnJT0k6QFJ4/m1XvxGvCRqhdSm8NqEiFTNjUCjb9GVZGUmlwF30kGdji44CPxXMzsZeDPwwfwzqFuWA8BKMzsVOA1YLenNRITqJOJDZAbpaXolxzvM7LSC31svfiPTJVFfB5xK9rn0Qo7hxMxqOch2YjYXzq8Crqpx/BOBhwvnjwEL89cLgcfqkqUgw+3AOb2UhazG8PfJvNefAQ5p9n1VOP5isv+JVwJ3kDmU9kKOJ4EFDddq/V6Ao4F/Jt9M7JUcw3zUuaRdBOwunE/k13pFUJnJqpB0InA6cF8vZMmXkQ+QFV/aAvwQeM4yNwCo7/v5NPAnwHQk+rE9ksOAb0vaJmltfq3u76WrkqhOe+pUeFGhIMOIpFcAXwX+yMx+3gsZzGzSzE4jm2GtIEvRM6tZlTJI+h1gr5ltK16uW46cM83sjWQmlw9KemsNYzbSVUlUpz11KrxWISK94id5eUlSlJkMRdKhZMruy2b2d72UBcDMngPuIbMpzpM0HW5Yx/dzJnCupCfJsmCsJJvx1S0HZrYn/7sX+BrZPwJ1fy/NSqK+sQdyDC11KrymISI1jt9I7WUm85Q2NwCPmtmneiWLpOMkzctf/ypwNplxvONQnW4ws6vMbLGZnUj2e7jLzH6vbjkkHSnpqOnXwDvJcrXV+r2Yl0StnjoNhsC7gcfJ7EUfq3Hcm4GngZfI/hW9hMxWdCfwRP53fg1y/Fuy5dl24IH8eHfdsgBvAH6Qy/EwcHV+/STgfmAn8BXg8Bq/o7cDd/RCjny8B/Njx/Rvs0e/kdOA8fy7+TpwTC/kGNbDQ8scxxkZPNLCcZyRwRWe4zgjgys8x3FGBld4juOMDK7wHMcZGVzhOY4zMrjCcxxnZPh/+jCfU0auIVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# READ SETTINGS FROM FILE\n",
    "lines = [line.rstrip('\\n') for line in open('setting.txt')]\n",
    "# Lattice and MC Parameters, Output File name\n",
    "L = int(lines[0])\n",
    "nwarm = int(lines[1]); \n",
    "nmeas = int(lines[2]);\n",
    "interval = int(lines[3]);\n",
    "filename = lines[5]\n",
    "print(lines)\n",
    "# Start on Critical manifold\n",
    "Kc = np.arccosh(3)/4; # Critical temperature Kc assumed to be known\n",
    "K = Kc; h = 0; \n",
    "\n",
    "# RG analysis setting\n",
    "b = 2; \n",
    "layer = int(lines[4]); #How many layers of MCRG to do\n",
    "\n",
    "# Derived constants\n",
    "Ns = L*L; #total number of grid points\n",
    "ndata = nmeas//interval\n",
    "Pij = 1-np.exp(-2*K) # Wolff add probability\n",
    "energy = np.zeros((nmeas+nwarm)//interval,dtype=float)\n",
    "clustersize = np.zeros((nmeas+nwarm)//interval,dtype=int)\n",
    "\n",
    "#Initialize 2d spin field\n",
    "S = np.random.choice([-1,1],(L,L))\n",
    "plt.imshow(S);plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=========== SETTINGS FROM FILE ==================\n",
      "Linear Dimension of Lattice: 64\n",
      "How many layers of RG Transformations?  4\n",
      "Number of warm up Monte Carlo sweeps: 1000\n",
      "Number of measurement Monte Carlo sweeps: 3000\n",
      "Interval between data measurements: 10\n",
      "Saving results to: data/testmultipleRG.txt\n",
      "================ RUNNING ======================\n",
      "running MCRG for linear size 64 lattice.\n",
      "Setting K = 0.4406867935097715  and h =  0\n",
      "iteration 1000\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-26-73e01a4845ce>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"================ RUNNING ======================\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m \u001b[0mMatA_even\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mMatC_even\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mMatA_odd\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mMatC_odd\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mRunMCRG\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mK\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[0myt_arr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mempty\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdtype\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcomplex\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-23-db86a7d81065>\u001b[0m in \u001b[0;36mRunMCRG\u001b[1;34m(K, h)\u001b[0m\n\u001b[0;32m     21\u001b[0m         \u001b[1;31m# Every MC n-loop, build a Wolff cluster and flip it\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     22\u001b[0m         \u001b[1;31m# Result: A S-field config drawn with probability propto Boltzmann weight\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 23\u001b[1;33m         \u001b[0mcluster\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbuildCluster\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mS\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     24\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcluster\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     25\u001b[0m             \u001b[0mS\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mS\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-18-d2cd609f76da>\u001b[0m in \u001b[0;36mbuildCluster\u001b[1;34m(S)\u001b[0m\n\u001b[0;32m     34\u001b[0m         \u001b[1;31m#if parallel to seed spin, activate bond with probability Pij\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     35\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mSj\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mSi\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 36\u001b[1;33m             \u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     37\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mr\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mPij\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     38\u001b[0m                 \u001b[1;31m#add j to cluster, add nearest bonds of j to the fresh bond list\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "print(\"=========== SETTINGS FROM FILE ==================\")\n",
    "print(\"Linear Dimension of Lattice:\", L )\n",
    "print(\"How many layers of RG Transformations? \", layer )\n",
    "print(\"Number of warm up Monte Carlo sweeps:\", nwarm )\n",
    "print(\"Number of measurement Monte Carlo sweeps:\", nmeas )\n",
    "print(\"Interval between data measurements:\", interval)\n",
    "print(\"Saving results to:\",\"data/\"+filename+'.txt')\n",
    "print(\"================ RUNNING ======================\")\n",
    "\n",
    "MatA_even, MatC_even, MatA_odd, MatC_odd = RunMCRG(K,h)\n",
    "\n",
    "yt_arr = np.empty(7,dtype = complex);\n",
    "yh_arr = np.empty(4, dtype = complex);\n",
    "\n",
    "for i in np.arange(7):\n",
    "    yt_arr[i] = getExponent(MatA_even,MatC_even,i+1)\n",
    "\n",
    "for i in np.arange(4):\n",
    "    yh_arr[i] = getExponent(MatA_odd,MatC_odd,i+1)\n",
    "    \n",
    "    \n",
    "print(\"==================== RESULTS ======================\")\n",
    "print(\"mean cluster size as fraction of lattice size: \",np.mean(clustersize)/Ns)\n",
    "print(\"y_t array = \", yt_arr)\n",
    "print(\"y_h array = \", yh_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x18620a770f0>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXe8H0W5/z/z/Z6T3pOTkJAeUkiABIihRFqoITSVagERRcGOguD16tWLGu9VRO/1J1LkIh3pAipFUGoqIQkESDvpvfdTvvP7Y3d2Z2en7X732w7zfr0g37M7O/Ps7uyzzz7zzDOEUgqHw+Fw1D65SgvgcDgcjmxwCt3hcDjaCE6hOxwORxvBKXSHw+FoIziF7nA4HG0Ep9AdDoejjeAUusPhcLQRnEJ3OByONoJT6A6Hw9FGqCtnY3369KFDhw4tZ5MOh8NR88yZM2czpbTBVK6sCn3o0KGYPXt2OZt0OByOmocQssKmnHO5OBwORxvBKXSHw+FoIziF7nA4HG0Ep9AdDoejjeAUusPhcLQRnEJ3OByONoJT6A6Hw9FGcArdklc+2IjV2/Yay324YRdmLt8a204pxZ9nr8KBltbI9hnLtmDJxl3aOv/yzlrs2NesrV/Hu2t3YO7KbZFtuw+04Mm31yiP2X2gBU/NU+8HgGWbduONpZsj25pbC/jtS4vx1Lw1wfH7m1vx6JzVsFnusLm1gEdmr0KhUJ6lETfu2o+/LVwPAHh9yWY8NHOltvxfF6zD1j1Nwd9vCffvuQXrsGX3gcgxT81bE9w/xq79zXhq3hq8uXQLlmzcXexpFMXCNTvwtt8/nhPOj8H6b1NLIbbvqXlrsHDNDry6eFNsn835sfYf8etfu30f/vH+hsTnsW1PE56dv067/7kF6v3FQCkN5Gc8/+56PDVvDZZuKt/9LevEolrm83fPQsf6PBb951nacmf8+l8AgMbp0yLb//7uBlz/6Hws3bQHN04dE2y/5Pa3pOUZyzfvwdcffBtTxvTFHz//MWX9Oqb99rXYMTc+Nh/PzF+HQ/p2wWEHd48dc9PjC/CXd9ZiRIN8PwBM+dU/Y/X+/pWluOWFD4O/RzR0waNzVuP/3mhEQ9f2OGmUfrLbba8sxa9e+BB1OYJPHjXQ+hzTcvldM/H++l1Y9JOz8Jk7ZwAAph7eH9071sfKbt59ANfcPxdHD+mJx645HgBwKXf/tu5pwrX3z8WRg3vgiWsnAwAWb9iFbz40D2eO64c/fG5iUNdNjy/AM5zySXI/s+ac//H6x+wfnIZr75+Ljw3tiT9/5fhImecWrMf1j87Hii178d0zRwfbl2zcjW8+NC/4WzyPy+7Q92++fQBYvXUvHpi5Cpt3H0h8Ta6+dzZmNW7Dx4adir5dO8T2f+W+OZixfCtmfP9U9OsW318Mzy1YjxsenY+V/vVpbi3g6nvnBPvLdX+dhZ6Afc2t5kIKdvoW2mbBejO22eS1uXb7vtRty1jj13dAYnEBwDp//96mZOcsnt++5lZs2LkfALB7f4vx+C2+dbh9b7OhZDas2up9dbVyXw+tiq8DZn2t2Sa/Fy2FQqROwPvSAYD1O6PXZf2O/SklLh3s/FZLzm/7Pu++bNkTPQ/VtUrL1r1NQR9KWvfyzXu0+9l5yb4yimXXfq+/btqVTvascAq9yiGkNPUyl0auRPXLSHIu5XociC9UgVPoBYVriG1VnUc+qCvcxuoSr3Op7msxBOcn2xecU3RvfT7bE8lzFyap4mXl84qLW8prnhP6kYV3sTRyVKbZjyBV9gC3+FqnLifvAiblpULsyElPm7Vn42/PAtYe77NXKvRAOcvPim3nrbOWVnadq6wDSGDXgEjOjypeTPX5bFVIjmtAHG8y0exfa1XPCe51CfpWWLf3L/taKzdOodcIWfdB1vEU+jxTkshOyvzmkylhlbym82C7+ZdDq+IlUO7ztKElUOjxfWlf8EnhrWuVO1BFc6tXXqWwS3nN2f1lL74K6XOn0MtG1gq5SB9d6HIpvWKh/snbtBRa6KWTR9Zeq4WFLh4jwo6LuG/8Bzsf87kkk7McMIUo6xPslEr9IuKv0/6EY1bshWR88ZagbzHDiN371gr5XJxCLzNpXRAiTa3FmQDsk1BVP7M0kspLhTcXX7/Nu4MVEespFYGFHvGhy8uyh1V9zRCri/0WFXoV6vNAoUstdMW5Z+2+yOVIcK2SWugMpYXOjIVUteoJfeje325QtIoplz83CQeai1Poth2veAOeJHO5lNlCZ3qW+boB9dePyUplL6HIoGhBodCrUKOzayATLXDRCYKL96nYZyVPSOB2SdvHVV1adItkiTi47hR6FVMN+ly0WJMOGIm0Wn6eZotZi7EHo3xiJfChsyMMFnrEh84Ueg340LUuF8Ux4vZi+1PUQk/Xx1Uv5PDrL3tygiHiXC5VTBa3Jq0LQfXg7y/SQmeKxhSil5RYlEtCvRU8dGW20KMuF3njYQiiwkKn0XJ8vbkasNBZlIisy6lcLqK1W6wiyxMSRASldbmYfejZdy4xbLFcM51jclSk1RojjZ9Q1WkShwEqVGuxFrr9p2HxmifR1Qv8nOUNW7QZFA3C+hR1hYOi4TalhV6VCl1toTNEA0PsPsW6GvK58OWX2kJXPa/BvU5VrRZ2VZzLpQhaCxT3z1gRdETAm1xw/4wVKBQolmzcHeSXePG9DVi5ZS/+8s7aYDYXD6UUD85ciXmrtuONpZvxyOxVwSw/voO8v34n3ly6BY/MCvcvWL0Dsxuj+VXeWLoFd/xrWUw+G1heiD0HWpThTz986l3sORDOvFyxZQ9+8bf38Yu/vY/H5qzGMkn+iJ/85T3cP8NbmjC00PWyPD53NW5+5j384Z9LUShQLFrnnb9YZvveJjw4c2XMqnpo5kq0+Oc/u3ErFq7ZAQDYub8ZD89aiftnrMAjs1dh+15vJiJTGHNXbMO8VduVcs0x7Ae8vCvrdnizAxs374nlB1mycRc27IzP7PvH+xsx/a/vY+f+Ztw/Y0UgPyuybPOe2PV9YMZKPDAjzAPz4MyVuOu15UHkRT5HsHVPE56atwaUUry+JHoNS8GfZ6/C//5jMWY1bsWcFeb8P8FAuWQfewaemb8Wm3cfwBtLN+OPry2H+Lo+0FzAAzNWolCg+GB9mOPm9n8tjVmtf12wLjZjNpfjLHTFVyilFA/PWhnp/zJZnxPqZ+d131srMLtxK+au3Ba5vzw79jXjsTmr8drizXhszurI9Xv5g42Yv3o7HpuzOqxbGBRtqZBCr+lcLg/NWol/e2Ihdu9vwZdPGgHA6zi/fP5D1OdyuOGx+QC8PApf/FO4OPWEQT3w5FcnR+r6+7vrcdPjCyLb3lq6BbdcMiHyCXfWra8Gv2cs34pfXTwe5/5vPFcKywsCADv3taBX53huEBWzGrfhhkfnY8ayrbhy8tDIvp6d6rFtbzPeXLYFP3r6XfzyovEAgJN/+UrsU1PMH/HH15d7ZUf35Wa06QcA7+eU1JQxfTH1N6/G6r7ukXeU5/LI7LDT3/nactz52nI0Tp+G7wu5TJ4/tC/uvOJjgeX64qKNeHHRRmUOjE/9/g3pOTIKBYpr7p+LQb064tUbpuDkX74SK3/aLf8KfvMK/eZnFwEAbvvnUgDAngMtuPrEEZEyU371z0hd338i2ndYX/ruGaMAeIrqGj+XSCkmtogUChTXPzo/ss2UT6SpRR3Fw0TeuOsAvnzvHMxbtR2tBYr7rjomUu7Wlz7E3a83onP7fCTHy8+eex+nHtoPIxq6AABaWgu45v65GNK7U+T4PCGchS5X6G8s3YLvPbYA81Ztx88/eURsf4F6ht2198/F8IbO+Md3TvbPy6v33rdW4N63wjWXm1oKuHLysEgd1//5HTz/XtQAYNfvyrtnBdvGDuiGQ/t3i/nQy3GPZdS0hb51t2fV7dwf5v3Y5ucA4beJyCz0nfvib/tNfk4J1c3ZuMsuH8e2vfHsdTpYzpiNu/bHlPRw/4EAonlTkvSfQoEGyinJp2GxoZI8G4XcJuyeZOWJYL5cVd4VEZ1FxfpUmod0l5+/Jk+Adb61uG1P6fPUpJGVWeimQdH1O/YH/aZZ+IRk93XPgbi7hI8kYr/4vDeA9yXD3FOqOHSWF0mWFRLwzp25a/h+pupbshxD6yxz7YgvHepcLukJY3zjp6Hrz7L8E7LBnPBtm06+oO6EFfADg2JuiCze/HV5Yu1y4cm0kwq3gEU2ZOVbZrLaTpzSnVsx136X7xbgB0XL8ainaUPnGuRPnX/cRDdK6IeP18FfP9WzlSPmOHRT2ooCpUFeF/5ZV+bgkegD1RiO+EUb+s7DtgGn0FMhG3BiHUn38MnyT+jKFzsqnlihcwODolx8VWlneeYJ4cIW7aNcSukXZA9nVuF87PxsL5Feoad/SJmfN09IWfPUpGmCRbnIcrnw/ZDvd+IladG8SG0SoOVzBHV5/aAouw91isRghQKw31fo7erCZ13Vt2R5dlRjV+IzEJvyT6Mylpu2odC5s5B1RhH+JjNkYUayiSJpSGrZhZ0k3javDNKqPgqu4yWQLctOKsqeuYVO1cpJhi6ZEjvrNBY6U+i88imHezWVy6WVKWN9Od6AEvtEYKEblKRKvFzOPLGItRFLp8DaoRQHfHeN7FkXkX3hq66f+NUgJuVyU/8tady8B4/MXoXH54aDbKFCl7hcNHXJLfR4uZZWinvfWlF0/uTWAg068BOKVYK2c352PgSKKXCVUtmwcz/+8s5aZdssFzkPpWGHY+f9/LvrsXKLfkUm3gdaLKKeZQpP9ojubWrBAzNWBtfi5fc3Sut8a9mWIIpGFWL4xNurg4gmHl0yJUoplm/egxfei7b7trAKlAzmQ3941qrUivydVdsxqzHZKlUyXnwvvgoQbyDYTP0X94v9UfdSiLhcFE+oTXKuwEJXKPR731wRjEO145511bv9ibdXR/qXjrkr5Pecndv+5lbc+9aKTJ+VJFhFuRBCGgHsAtAKoIVSOpEQ0gvAwwCGAmgEcDGl1NzDU3LKr8IojolDemFw705yC93/V3dv2kkUusz6nLF8K2Ys3xobuElKIfLQeMphWJ/OkTLfeGge/vSFSd4f3Ii5KBaveAgBPnvnDCzWLPF1+V0zpfKIE4uuvncO2uVz+PCnU8PGBbK10KNPV/BwSp66nz/3Pu59awX69+iAU0b3xZX/NytWBoiuHqTyoX/7YS8iR4xo0lroFDjFj5Lh+cT/e0N5DGNPEwt9BVb6/ShpjP35v3sdQLJVb2T9/4t/mh2rg7+lOpdaxIcecbkICl0zsBp1ucjbyefCgVaVy6U5cLnI7dGHZ69Cx3Z5AEC7unywXfW1tnDNTnz/iQUY1KsjThjZEJOV5/I/Rp8nMSXvrMZtmNW4DRcdXfrVtmQksdBPoZROoJSydbRuBPASpXQkgJf8v0sGf33ZjW6RWegWE1Pq6/SdTWSbYjTdFlERykbnN3KWdOBykcgV/ZsYR+NZHLZYR6DQOdlMUSylzPGcD3zocdgqOXslkRMqQoUe3ydbDUl2/9v7n+vFvMb2SVZ8qqTLReyL/D0Nxwpk9YW/dT700A+vl0k98Y6g1a9DaaH7AupyzLMIsKgPXQ8f125ruzDDRDydtZLnrhwU43I5H8A9/u97AFxQvDh2iO4HfmxEdYF5kg6KFu9Dj/69SxNSCfBfGeGgaDjoIv/0FWGj+zI/I3+qqnOTbc7UQhddLr6cWaXz1fnQZS3IPpGZQi8mskiWoqEcH+OqNsS+J8thY5q2zvvHxbLsBSHrdwWLfgcaWuCqsEV+spYKdt3bJ1iAozXi47e7S6pFM8RFwcuF7dlSAM8TQuYQQq72t/WjlK4DAP/fvqUQUIYYGpQ3+MnEmyNzuej6cLHT0AuFaA27DGtrEm5QVOxXvJwEaqXOztGs0O3PrTkjvyClNCY3Cx0zDcjJZvXJyCLKhX2uF2NRy5RSOSadqBSSON+iORIbrh7Q43swf0ljg6ItzMiSDYqaLfQCpcE9TutD946ND4qa+oIsrNKEmDaXUY65BjJsZ4pOppSuJYT0BfACIeR92wb8F8DVADB48OAUIsZhX4mysEXZPRNvTr1k5FtrfRb5/Im+SZNCD2adQe9yIUT9Gck6stGXWaDW1khWFrrMV1tnGeViO7mJ9RHZ+cusdpkSa28RIWFCtrB4ORI3qZoQJ9zZWujKOHThujUHefb111jV5SjC/qGKcgni0DXWNzuWV+imrz+bsEqbY4Eqt9AppWv9fzcCeALAJAAbCCH9AcD/Vxp6QCm9nVI6kVI6saGhIROhmbXw1DwvuoN/U7N7xk8rFy/2pl0H8NyCdWhpLeCeNxrR1FLQKjXbm3vvm43S7aKy+OvCddixrzmS++P99btiGe0KlErCFsPff383HrXAYB1ZZqG/umRz8LsgfAUsWL0j1g7jbj91AIAgJ0waXlq0ITYoGoYtRrf/5sXFWOBHrgDA3BX6/C0M9um/Y19zLM+O7IXypzfi52MT8mZCZqG/IIk2AYA3lmzGgtU7sHzzHtzy/AfKsZtVW/fi2flh/21uLeCvC9ZhxRZu1XtFl92xrxmUUtz31grc+eoyPDN/bewQvr+u27EvyD/D4BWj2E+WbfJkeP7d9bG2X/fPD1A/U4/PXR24v/hB0Q/W78KvX/gQuw+0xNZp3XOgBfe9Fb1/M/17bhPlIjsXmXR/erMxfgyiLlHGbiHPzE+ffU86Qz1rjBY6IaQzgByldJf/+wwAPwHwNIArAEz3/32qlILysGvHrB/ep8cUxaJ1O8PywvEzl2/FzOVb8ZPzx+FHT7+LAy2t2gxsKmtHnFn370+9Kz9eqGDJxt34/hML8Cz30gG8pFCnHtoP4MYBAh96IEu0rp0Ka59N1JFZJf/+5MJQNkojD/C5//uaMppiNhey9W9PLJSWseEr983FCSP7RLapVmr/9YsfRv7+7F0zpOVE+Ot04W1vRvbJoidmasICi3GRyPrOO6t3xDcC+LSf/+eSiYPw8OxVGNK7Mz4liZY4739fw7a9zfjhOWPxk2feQ3NrATc/uwjt63L44OapWpl3H2jBG0u34AdPxu9fYKFzx1559yy8v34XvsDlOuH7lCre+nFJeO7vXl6K3728FI3TpymfqVmNYR/jv8Z++twi/OvDTThiYHe0FtigqNfHf/bcokjOIR7V5CMZ/NeK7Gv0h5LnO7xm+rrveHU5Lps0GA1d21vLkwYbE6QfgNcIIe8AmAngWUrp3+Ap8tMJIYsBnO7/XRZ0D5hML6hcBeyNua+poJ8pqpLDMuhDbJ9Ank+GvdV5l4v45rdVLszi1Q0csfoqNatNRnZT/9X7ks4rMF3y8ycMSFSfCXaPtys+21lumTXbvUgK9lLnfc4qkVtaqTJLYWht8m15Xwm8C4HvUmn7js24FC/HZv95aSlQLmzRE0Rl1IiYupatD70vp5Rl68iqkAVjZI3RQqeULgMwXrJ9C4BTSyGUCc9NYO/vUnU69gC0r89hr6KTA+oBHNvZYDG3iaF8OChKYy8N25dIEoUuO41y5SNnBFEpGU/9l5F04QTTtVDlFEkLi6HeafDDsvBL2UCxqs+2FAqoy+el+9gh/KFd2tdhAw5EJr5FXS4pFbrFYXwZ9jXMJ5ZjfVuWm0mKwVrgRdKdF6+YbS108bhSUTMzRXlaC4WovysyUGg34AWE/s32dTmtclbtss1too8lj2/nJ0cFZYOHza5NdhlMUSOFQmWmKYu+5aRRKSa0Cj3hak+my2OtUBK2p8sYCoSzi2VJtVSn740XqdqN7+jSwUv7vE2h0NN+3Nl0ueiEPO8cWymN+dBlUWsMYgiYiLTHnYzuvGRjKzbPZRL3T1pqUqG3tEZtJjGUT6RVEW7HJn3U5XOpwhZtoxVE5UIppGY66xO6QVHbB4j5pE0j+60VcrmIkT7s5aiTNomy172k9idcCcd0ebJ+UNm1MEVKMJeMzIWk6rNe3fJ9oYUe7u/WwfuI5ydj8R8kafuOjYuCL8HCK1sLNPChMwvddvDaHLZoJx//Akmy5Fx9xl9yMmpSobcKoXYFSvHonNWRFUR4VJY0G1S9+7XlkQgOEdW9XbBmB255/gOjvKJbY+OuA9JBOPGTt0CBVz7wgoeWbd6Dvy1cn7kP/aVFG4xhaqVAdCcs3rALP3pqYSSixYbV2+JpGVZu2Yun56nz2+yXzN6UEfax8rpcWJ6ZnftasHa7esYhs5qbZAaLQuTmVo2F7v/Ld4euvkLfqrTQk3eUphb9mBUv6x9fW47m1kIQtURp+Dx/sH4X7p+xAs8tWKerJsD02v3rQj4yTl2On2nOTuOZ+WYZZDPUs6YmVyxqESbqFAoU3/2zl6PjW6eNjJVXWRHs023Z5j3S/UH9ms73238sMUhrb8WIES2UUjw0a1Ww/yv3zUGfLnaj5MTSQv/7uxvws0+UX6GLYV0fbtiNDzeoc9Iw+nZtj43cgPKFv38Tb30/OpRzzv+8qh0ok8WG6yi3y4Wx+0AzLrvjLfV+/xyTuFyaW9UjAqF3LyzRod7zt/NfVLYKvUv7uth9Brxkazb96+2V2/H2yu040FIQLHTv90OzVkWeDxOmzJuvLt7M/WXnQ2fnP1uRtIsn6xe/jBq20MO/+c4rU56qHCQm65VRrHKzVei8IlfLYmuhs3+9czx+RG/ccflEadlK+NDTrn40+qCukb/XS7JJmqIebBU6u22my6Ob4FIMTS0FabZMBrNUZQpdpbb1FrpvUERcmF7/4ftwNNuiUjwM7NlRul2c+2Biw879ER96Obqr1kLnB0UT1FmqFz9PTSp0z0IPLyXf2WSzMFUK1fZzsdj+Y92OIeY8SV3Mhx4uqqBZrb50ObeUpE1JnMXDrMoRIiJmpFRRquiF/c0FbdRPs38NeYUe9CGFyC2tBc1qPP6/3DaZzROJQ9doPtVLu7UQX7hFx54DLZEolyQRWLz4lvab146lD93WwKrLEevc/MVQkwpdjHLhO7QsMkDpcrG1nIvUIq0JO6EsfIxhOwbFJluxLlSgFKovPnnujtKSOjoiA8lkCbNkiC4wFfVJNEUCVOljGWyaPZ+PhYVkqhRSU6vawg1fBmEBmcsur0nOFZFPodALNNld3NvUGpxjUuueJ0lIrK6NemGxEhv9UI6QRaBGFbqoiHlrTxa7q7TQrRV6AuEk2Oa0CsVRK5KkFjqCmHa1D7EcuUWyIouvCVsLnfWzSrlcDrQUtJEZbAyIV5xMoessdFUfYltV7hVGztLlwhJ1iSS10HcfaAkXNU/4MoiQkYXOp+suUDvjpBwhi0CNDore+epyzOEGIfiJIi8uiqeUUVniaePIk2KrMO94dRnOPqw/7nptOQD5m99WFNGyKlCqdrkIlX7tgbmR1AnVwl8VUT5zVmxNFFtuq9DZfXtsrjx6ipGFb1R2rw+0FCL3bOOu/ZHVmmQ+9NmNW1GXz2FYb28BFUKifcZmEYtIGLDk1IjloKhsfIMd8+bSLcrjRPY2hW7Uf39yYUmnzxcKFLkc0T5n/AuNT3GtQ5cZMktqUqHPW7Ud81aFSZpMA2wqCz3pYGVabF8IyzbtwfifPB+2K3W52NXFjIiIy0Vh7onXwSYEqxL85Z21mDS0V2z7p37/pqS0GtuZorYv/CxyuMuaOiC8eK7+05xIv2fwLper7pkNAPjn9ScD8D71+S9Y7cSiYA1d3kKPn1vUQk/+dBQKkOaSUSHeryRJrnjxbXRqK6XIgWjdKOLEKhs9Ug7/OVCjLhcRseOLqB7M1gLFwT3kI/E8xbpc0h5fjMslJwyKFqha8djorSG9O1m1m4TDD+6e+JgsconbDsjaflllYXzJIrE8l0tYuUqRyaf+e/+KCzy0FDQuF2ah8y4XSTneh55mYpFNVBXfRlaLntiEDYYD4eoyfDUUdhZ6mQz0tqHQzRa6erSd5c3QUbTLJeXxMish6fMTpBGAeqacjXyqbIjFIFsZ3kQW3n7bkElrCz2Dp1WmGG2/JGQLj7B7Ks6ibG5Rzwxmm5O5XKxElMqmg3djZaUMbe5TS6DQ1TLy50+dhZ49Jv+p6vltKRS0eSAYFQjTBiB/WGwjbsI0AmGiL1WfsuuQVs0mIo3rOQsLXRV9IWIbn5/FpVG9PPiXj+rc5XHoHmJ0RXOhoLzfsmUOZU0Wm5zLapo8J3dWFrpNLSxNiN6HHlXoNi+1MhnobUShp1zcuLVArSYXlTvzoK5de7+/Vy6ND12G7SSsJKSZOZdFQI6ty8XanZCBwlHlG+JlVSl9qUL3i8Ys9NYw/WxMBol1KitabPpcmxclb2gV8wXEH2ljDDBdoSsrjiHYvKCyeikZ2ylLKyXGbKGro1xs+kqxRuGu/S2469XliY9btTWex8P2+SlQ4JFZq4IpybqJRT97bpGxvlJ0yDQzod+RDAomxfYa2iqrLN51i9abo4pUPnSZy4UqXC5/eWetcl1WmctFZoHz7rc0s4zZqkY6ohZ64iak2IhqM5lMTH1w9xuNxnrLpM/bhkI3+UR1g6I2vq1iP/PXbN9nzBeThH7dzGFbTS0F3PDY/ODv688crbR03lCEkPXp0i74XQofYI6Q2MpFtUgWOdw/fYfdSkwypIOi/r+ycLkNO+UvBtliDbK+XydMrEnKtffPNZbhE1lllSPfzkL3XS6aMht3heGY81fvwG9fWmys17lcEtBkmFGn+py1tdCrYd7Nl08cHvz+9KQhxvLiBKuTR/eNdapxA7pp63j9ximRv7P2uuQIwb1XHYMLMl7xp1Rc9fFhOKhbh9h203UZ0ruTNGlcVrBsi3zuFKa7eOU7tr93v1VrlQYzYynvT4+X69SuLnZM1nTz87AXCy9dEgtdLMvnQeI9uFsV11LEDYomwBQNoLPQbVwJ1TCTMhrGZS4vW75M7FSmyQ68j5sgez96sAhHuWK6ikR1vUxdKEcIunfMRkHJYD70TlzEVrCqDydcj06eDHwqXB6Z71ymBHk3TKly6fMKvRhdyMufxEIXiUTdcFpTlk1ShnO5JMA0yKW6ka3+rDATqkGkcsIrExuZt0seWvEw05T1yIw42PvRbfUze8GUIiSyFKheaCZPyJ6HAAAgAElEQVTriyA7i1MGc7l0rA8VOnND8jIHCl1hVfKem9Bal/jnES+XNZ3bh18BRSl0/reFqKrxBX4sgn8OdluuZ+oGRRNQnIVurl8Vx15O+PwRNn1Ddsqi4jFZ3GJ5205pa8mzUuXq7DxppmIrFbrhOEKAbiWy0PM5EgyKduAVektcobOvBKXLpRC3vGUKO2L1lujR4LtEMe8M03iAiNpCl6tKZ6GXgCUb9Qsj6CYW2SiTFtvsWiUkOoU5Xe8Q9VGSHCRJXC628lXS5ZImWZLq/E3nW0qXS/u6XGCN8y6X37+yJCYbe6moBuj5iJVXF2/GW8u2yOdC8KmryzBJoyi3DgXufn05Xl28yWoRij+92Yhlm+L6hDcA+BfDa0s2x8rKcBZ6hqgUcouv0E8d01d7fJoO9fUphyQ+Rgf/OZ1W/4mdKkkcOE3QblLFLzN+hvXpbClZOtKkM03rQyfEW70nCTdNHWNVjr/W/GDlyx9sAhB9cXVpp5eBV1Rf+tNsXHr7W2YLPYVC72QxO5uvthiFPm/Vdvz4L+/hc3fN1Jbr6bujHpy5ClN+9c/Yfv5ZKRTCAWZbXJRLhqiTc3npSf/jvHGpjldxzckjcO3JeoV+5OAewe/fffoobdnrzxydKLfF4F52eVeSuh1sLWlbn7jO5XL9maNtxUoFm7iS5BrkcznpZC9R/POFqJ0cIWhfb/+oTT6kN7580ghLmcLGeZcL33bw23CuNmvL9unSLrFfWqST4cUiovsKOG+8PkLKNgrlN5ceqd2fFyz0Wy4ZHyujfY6dyyU7dOlzc4QYO7ptTg9GjlhYbUJ5HfkciZQxKSFb10AStwOBvaLOW9YbuFwk9ZYyKgQIzz3JnVUZ9aK+Eb98CCFob7kyPRBPrKajLmKhxxU6v9/05WRyr3iQoq1nu9QLnO9b00avzu2U+wD75QaNEV9C7L0sNl5XhbPQM0QX5ZLPEaOiUo18A+olukwPD6/EbJR/LsGDqdoteliSLMxAUQILPXC5VEChp5immlccI+obUTkQAO3r1G6Gbh3kFquN35UvI0s0x8tsUloyS1gcfvKaK86HbptLh6EzqFjkTrGYnil+P4U8mEL3fFRdHDohJE8IeZsQ8oz/9zBCyAxCyGJCyMOEEP2rsoLofejmKei6DiVTDAR2A2VBeUNZKpQxKVZVfaJVkXTpNFtFbav4WalKKHQWhpYkuVRMUft/ihak+IWSy0HrcunRSf7o2FxFk8uFf2eb+qTMEhaVL0H0iyRNci6bXDq2XwE9FdcuKaav1eigqNwI0z0f1Zg+95sA+KQfvwDwa0rpSADbAFyVpWBZMqtxq3R7k59v2qSo9MtRxY8lhJhvIJH+lCLmYTFZWqq94mEqi1NVp/WgaMIoF1nx7hlZXqa2k3gMxHtd718/sX/ELXS9y0V8ebHqklroMpdLPsGXnUxximv0iisgpQlbTOrC1Cp0g8vFFtOzIPrQZUaT7vpmlb7AhNUTTQgZCGAagDv9vwmAKQAe9YvcA+CCUgiYBbrl1HLE/ODIkh+N7d8NB/foiLMOO0hSp3mF72i+Z5OFHg2vzBGCU0Y3KMur6hM324QtfvwQL9fKTWePsba8bQcAgygXibydEw6cJWXZpj0Y2bcLbr1kAo4ZFl8FSYb4wHbwz1N0O4jXP0egTdPcVXC5BJkybXzo3D2UtcErKuOgqERvbt8rKHSQsoQt8rXq2uiaMHpIhckIqRPWEZWV1lVRbXHotwK4AQB7H/cGsJ1SyqLqVwM4OGPZUvNrYQRa1yFsBkVlPr8Jg3vg9Run4OzD+8f22dw8PozNxlDmRcznCH7xqSOCv7t3rEfj9GkY6q8qpGo/NvXfQqHf98Vj0Dh9Gk4Y2WAdjtjD0l3CxJF+5VjVUBwvXHcSLjjyYPzoXH2UE0OUk1mHoiIUyxHDC17MiMisXtUhvOJmisgzTCQyc9tMt0/2JSqmkIhb6KWPQ9e1IV47HaP6dUHj9GlonD4NPzp3bGSf6Znln5WCYv6K1kKvFh86IeQcABsppXP4zZKi0qtOCLmaEDKbEDJ706ZNKcVMRpf2UYWiWn0csBvAlPng2YMkM7xs9F5karNBfVEatfryuehLSOzwqs5TTBy67HgVXRSDfCLsvGX1lnPyqG1boiuFvbhiPnSFr12FaFkzC9jK5ZJj/VBumERmGEPvCpQp9B0xCz36AivV1H/eN68fw7LvKOJXbhLEiUVJfejl6s42T/RkAOcRQhoBPATP1XIrgB6EEPbkDgSwVnYwpfR2SulESunEhga1myBLxJusWuAC8B4I0+eWLD0va0KujMy3r0uCXBXiakNeGGO4QfwCUVUnbk8ch25ZPPFM0ZTXMCtsmxIVNRvMFJVa3OWib0C0Mo0+dL4vBCGO8q+A6KCoXhaZr1rs+4REF1AulYEeyRejU+gJIrXqNMvame5R1IcuL6/70k+T+z8NxmYopTdRSgdSSocCuBTAPyilnwHwMoAL/WJXAHiqZFImRLzWulF1m5hxeRvFWZcRl4vFAfzDmhcGXUULRlVdPA49WS/Lfuq/+iunnPBfSLpTjLlc/IFbURGK52O6bDGFHshlhsmkdLlwG4lBodta20nzo6SBr1ZnodssIcngDTfx5WfqspGwRcU5665tVQ2KKvgegOsIIUvg+dTvykak4hEvrK5D2LhcZIQPkuRNbaHQkmST81wuXP3C57WoULIcFLWp19SOuj7/3wrkcpHJYSLmclFY6OLD22TIBRRTSn5xCwM96IfeS17v1yUg2reETcQKIcBiLnfSq4vtcpkkQTyNjYqVmoDoQhgmcsLLTddm7FiuQIHKs7Tq3i3VGLYISukrlNJz/N/LKKWTKKWHUEovopSqr3qG2AyC5AjBiIYwF4huIgMh6dK3hi4XSZ2GYwd074BPHhWOIZujXAQfuvDwMoX+4/MPw9DenXBwj45iFZ5cQjOl8qETACeNMrvXWG2VyLYYkYNrXudCEEPbWCiiaZ4MWzZv2hHhADofqigOTodRLurr8pljBuPKyUMDxZJThMpG+rZEn7evy2G4nzfHFLHyX586AoQA765VR411tsjTYoLAfgZvkpw8/LUY279bJEWG2AdFPcO/dFVRLrp+fNPZh1rLWQw1N1O0pxCf/NNPHBYrkyPAf54fbmdhh/dddYykLNG+nX/+ycOl24M3tORYk4J646ZT0b97qHSNCl2YyJDPyx/ek0Y14JXrT1GGDRYz9R9IYnkT3HxB/L7E6wuty2L44TljzYX0kliV4i30Oy6fGPSBmIWuqO43l0wIfvP9SlRKoQ9dLctPP3E4fnTuuEAmQuQvgIhViriivP3yibj2lEP8duVqdNKwXmicPg0Xf2yQ0XXw92+fqN1fDP914RG4lbuGABKlVOC/Vo4c3BP/uuGU4G/+rHp0qseHN08N/n7ky8dFgxAUC66rnuMfnTsWxw7vbS1nMdScQu8qLBQgG9gjilBE2fXOKR4EhsoCYDcviwgN0+eYGIeu+rwWZTNtTzooantennKxr69Yl0vPzsVNQrJ+UQnKkb2IYlFGiuNVq07FFHpQxixYGG0lH9yvi7gZJH2Vk0U1gYeX1SRSFoPZOWHgld8u9hVdSgUR27BCsVTMI6aaKaqov5zfnzWo0KMhcfI3pb2iNfnPVX7mYqNcksgQs9BzeoWu2hNzuSQcjbQf2DFPrOLrS5GaPIJq6rwtti6f6KpRnCK0HBiMKA3ut9jHmDJT+tC57cwL5Cm7eNmc4iUSbgu/UFXuppxG2cVkU2zvkCDbpHpQP/41lyQOXa/QdfJEd6rCFtVjV+VT6TWn0MW80rKblMvJXRIyhWSe0Sm/RGG4WHxfUoPTxkKNxaHrDrFQBEApwxbtrJKsLHTbiUxKOSzLiQOMoctFrDDZ/RTHMoIoFxsLPXC5KAZFIy+R+LnykS+qiJUkeYdUu5OkzFUZDoTEn/e0Lpd4m5rjZApdcoR6iUIr8TKh5hS6aKFLFTqRK4k0ylel0FnHlh2edJDP5EP2BkW58jlDZI7C0hI7YWI9miDKxaao7isnCcUm8tI1z7/06nIkvIa8IrR0uajKiBZ6wcKHzmAybN1zQKps64SXUEwObpPK5bK3KVxmLe2d6pBA8aoakUWkpQ1blNUdNC+UE9ssFJLpEudykcAs82+cOhKnHapfYUhlreQIwXWnj4ptA4DjFIMWKiuWt4zi7WvFk8qlRUgGxHLFjDmoq7T4qm17pdtjzUTqNMtpe1o5g48/rM93uRRpodfncxjdz7sWNtE1KjlkRFaK4r6McoTgwqMHYkD3DrjkY4OCQeBD+nZRRhnx8EZ5zGhgLheLK77UDyEsUPk91IXqsTbY9W9plaeFnedH6fgHaFHJ3CFB9AuBeuEM0SJP8qzp09tq5PH3sYVLfvqJwyL9+zunj8KYg7qiv+q+O5dLnO4d63Hh0QMx5qBuuPOKjwU3VjqyT+QDRIR4LwSxLAA8ePWx0nZV94Kf0BE/JqErw3AXYmGLfqPPfeMEaXlZMjEg+ix+94xRgewTh/TETz/hRV2Iq+1Ejrc8Ldmnva6+YqNcCPGiKxqnT8PIvl20ZWX5s3XN84qoLkdQz/odgAE9OuKNm07FoF6d8Nljh6Bx+jS8eN1JVn5dvo+IYxnhoKixmgiyl6jJrZYj4QulubWAYX06Y1Q/7xqylyNvuJte1KrdHSWpfdUyySvJERJbbDvJs6afmq/ex56331x6JBqnT8PJo/tGSn/91JH427dOTLzMYCmoGYXu+a1CglVdJGXzRB6KqBoU0mGapCMdFNXWGMdqUJT7mz2kKotD9ekcX3uGvZTC7pwmekbENIElKKcZh0hCdPEBPbJPdF37vCLK50hwvHZQOoGrxJMpegCzTm3y2tdzLw+ThS6ti5DA5dPU6qWTZvV3bh9XwqZTU+1PotAJka2U5Ct0yzxBMnQrafGXWiwlexEkcRM6l4sEMUFVEAYutZBVF1x2Y/TtKi10EirDeJ0JLXSjDz06qm4qb5Nvms/ZHunMuk9PY60erQX5oFFcBu/fYl0u/ANnCjiRjYnorDw+OqMulwus72JfQvwpxy10+/S50an98QNMXz+8hd7UUvBj1b32O9bHladRJsV+2eIbSSGkuPESrYWuOS/pSzHB/XeDohLEUKFwUFLhclEMlsa2GT9J1Z9/nhyyfdoqrdtgUBpVVOaFCuxWHZB9ZegUsW3HbFWEdcXq8/8tVqEniZKRuUN0R0ctdFgp9KQLU4huEXb7rFI4GNwhpmtLSBhl41nooYtFtmCG6UWt2p8kvNCLQ5dtj89DSULa9Laq0GhbypXHBagphR696OynVKHm5Bc8zQCmar+u/aRvZKPLBdGQMtPAvspCFx+SIMGYkJFPhW3HbGktWJUMX4pFKvQEx8vmFeiO7xBR6DmrqeZJXmaAJjmXRT18DLypz8v7f+hyaW4pgCCc1CNV6Cmfl0Rx6Jo9SV4MInp3ovo4+XhcApeLs9BlREfgTb5e1aw4WVkdqv3njfcGD0c0xAfhEg+KGop/etLgiEeRzylywYQB+DchT4RqDVUxbQJrNid5Ucqwt9Dl1+DKyUOl9RU7KKo6XhajLFtqjD/8tEP7YlCvMFphVL8wkqguFy4lp7rGgN0LJpI9U+gATKHy9ciUKxB90aexJAk4l4tvobMaZYtOq5gwqAduuXh87Bkb1qczjhzcI7YC1cQhPbUD2CoLHQBOH9tPegyfm0WGboCYN1bEyyidsKWo6junj8KJQqRVGfV57Sh0cWHWwEKWlM0RudUrs1xNHV72kBwxsDv6dusAwLPgPnHkwcZjkrbBmDKmLwb16hR1uXDlb730SHzpxOGRY1SDonX5HL5y0ggA0XGGrLMoUmEAGwBG9u2CC48eGK3P/7dIj0vkgWP+3x9MOxQf3Dw1lixKZmDzzQ/s2Qmv3jAl+JtPqJXPkYi/WSmPX+GZ4+SKhy8DqBUNu97HDu+lXCaP967J7o/p6y9HSJDTp7nVC4+lOpeLohM88KVj8MmjBsb2n3ZoXzxx7eSYHJOG9cIL150krUvVz1g/vePyidL9fG4W6fGaa6G10GVzWhRq+uunjsQvLzoiWtZZ6HGokBCH/ZK7POQzKWVZF03WoWy3+DCLeSeS3j+bGWx8G7IIAB7dorr8saEP3dZNYHdmqlwX4vFhlEuRFrp0Epm8TunLi9sk3ssOiigXXQZPVqFubDo6JV+00KNiiXH9JFJWb6HbzOzkI38IV6d0AFlRjy7qjLXDY/Jny/p4sYpR6y3T1J00X1PSLKZZUjMKXVwlJOyocp+o7CbIrCpTh5f1O9kKRpFjEl5Vm8GaJCvD6FZoCuoFkVrouigRawsd8igX8XjdOEQSEoWQSd0S4bb4gs/h7zwXh67rA+wY3bUU6+UpCC4XlQvRK8vXKbMk9XgWOqfQOZeL/EWpr0/crwrBTJNXpdgZxTpFG+2v0XaSKnTRSHSDohJUeSakg6JE/nklU+jmiRJ29USOSXgDbfzWvMViCs3Thi1yu0KFSqxktn2gKIVUk8QUOhcHXwy2LyRAngiM3yTqaTEahVmzBzR9gPUZVTpasVWxr4rZFj33WLx+wDwoamMH8APFhMRfKBGpFbfKqOiFv03rb8ouXSlf/EldLrq6YvHuzuUiQbTQhX95cgqXi25tUBU2ln6K3EwR9Dkm/Da4RkwPaathhRzAkzFcRs/sxmHH2KByuYjXMpxGb1eviiRhj6pID0Z8we1oO+38FXJsfOi6Kxqx0GMul2gcuiqVBV+WlZMUCPdL5SAR14oX5eLLlcB3zLaL+1URO/pJPoo2inbNqffpwxYl5TXtiGMiZdTnqPxcVUu8pPLh33xnF+EH/AAvuqNPl/aYfEgfAECfLu2weXcTAJs49Pi2X1403nBM9KCjh/TEJRMH4YbH5kvL53ME/33hEbj+0fh+9oDwXyimXCE6C53fw6Qs1kKOtxFfAICQeMcOffjm9gd074CO7fLYvLsJZ407CA/PXhXsS5KnLEeAqYcdFIlE4A+/+iRvgPnak0f4Mod78zmCL584ArMat+Gsww5StjlpWC+M7tcV150+Cj061eOwAd1jZXKE4MfnjcOcFdtiCkDMtuh9ccrbYrf69LH9UlnohER95TkSvgPyhOCTRx6MIwZ2j5SXESh/YT8bcBXvsS7ixDMw5NtlfH3KIdjb1Kqsj6F7bnQ9UL7cHMGxw3vhC5OHxcvH+n75VHoNKXSxMxHu/1F4H3q7uhze/uEZkf3/c9lRuOyOt7zjjZ+K0QLnjh8QC0uKHSPU+dg1xwOAUqETQnDRxEH4f68sxfLNe6Ln4j9r7CE7b/wAYyyuLj93YP1xcuZydm4i246pWqJLPFw3OUukR6d2eO6bYe6a5Vv2YObyrUq5VFXmCMHvP3u0VA4gDEO94awxAID314fLreVzBEP7dMaLiugMRtcO9cHKPbdcPEFapi5PcMXxQ3HF8UPx1rIt0Z1CtkXvXsnPiL3o/+O8cVi7fR8AYHhDZ7Svy2PRup1GFxQhwspVJIxDz+UIbhFWCFJd1zBZXXR7Vz+/SfzrTHfT5Rpd1Ue/c8ZoTV0hZ4xTv4S1LhTFV91DVx8nLV9JC71mXC6i1acbUIu4XCQdg3cvJI1ykZUWH5rE6XMt4mOTTAdvMS1yCeYrZb/DSm3dBFoULhdV105lwVgoKlmLpigXkagPPbvHha9LNSjKf8GoolxYRFN9LtrnQ1ed/kIRkHiUSyCX7ADD8yL8zWZ2iodpY8IVu4p1zWlzp2vqTjqTuZKLnteMQi9YDrQB/icqiSpCHr6Pp51YpD8mWXn2UpE9fMGgaGC1mSvXjonyvnj/j7zMHyLBXp/bRbmEg3+WFacg5nKRxaFrFbr8d7FEl6OTu1xUg6KRsn7h+nwu8mIM+78eCiosU2cYFDXUJ76cWQZCsT/kNQ5tAvlzW6zror0mnwxfdcyIK/K+l9HjUjsKXT0oKv8cYm9JmYFSiAwk6ZsVFYDdzUl2B3V1imGLWfUNAsI9uLAKh0jkclF8OfHIZkSqsI3aNFmkSZWUaep8WsSFM3jEOHTdoCi7h3V5ElHiNqGTgDcxKboSE/SDooZLIO7u4mdHjA2KaipS7Sr28uss9KQulyQ4hS5BNSgqexpznEUj689ZW+iygbckaBfLYG0En+HF9Q5eVlMSKBHb81KGmIqyBC4fu3ojx1qreEEGSWPaVLipWjHDK0v10mXe9tjEL+43u9b8soSUWxBFdS8YXibP0O3iTerx5UploUf/Dn3o0e2mafjyqf/F3Q19m2qKTR7n4tAliOv4hWFScdrV5UJrRdIzjh7SM/htk+viu2eMQp8u7ZTtxY5J6kPXyMq2nD62H8Yc1BVfPWWEsb7bPnuU1cBtgXtJnDHOq//ak0fgvPEDMP2Th8uOMrYNyMMWCeI56sOl1vT1ThjUI1gRyMQXTxiOMQd1xbnj5Qt1yBNYqevLOgKIwQ9Exl0u3JcTEJkEJvL7zxyN40f0Roe6fMRC/4/zxuGIgd0xdkC3oKysCtblwmgURAZFRYwzT7k+MmFQDxw5uKf0HFndl35sEPp1ax/Z96uL5VFkxbq8ZLJfePRA/Pi8ccoxCtVxydot6vBE1IxC5z8jAX5QNH61urSvC96qMn9y5/Z1QY4OmwUAvjZlJL4vJMCKyKaZXWiDTgZWd49O7fC3b52I4ZJkYCJnHdYff/rCJEV94W/e5cLX/9vLjsSlkwbHjrXtmN5wh8y6kysuU7VPfnVy5CUMqF0Jg3p1wt++dSL6dIkqiVv9aA25y0UtQakUus5CZ19OfDZMVRc5ZUxfPPClY/3l8aLK9Omvfdx6YQkWushfC5176uN+CDCAmA+e8eRXJ4cJvhQW+vRPHYHPHTsksm/yIX2k31+luBW/vGg8rjh+aFmVbikxKnRCSAdCyExCyDuEkHcJIT/2tw8jhMwghCwmhDxMCGlXKiEppb7VZ/cJSIh8xaJopd4/5gdWnGhivvNplYCs7nSOBTtsLWSGqpSobKgiH3psUJS1X8JR0XgEUryMfgwjW3kYuigXho2FHikvhLgC5v7KXurBbFHC+9Dj5Vl1vMx1mklCgWyCHPzxSdJa1CLljEO3sdAPAJhCKR0PYAKAswghxwL4BYBfU0pHAtgG4KpSCSnztbKf6hAn/UUMP2v15ZpaaEQGK7K8fxl3dt73XNB8WstQXSsxnE9cMo8RV+hR10Ip0d1vrculRMLpLHRxvIQQYXBWUWcaQ4KPkmF166NciF8+3FdvEc4pXkb+nGXJ5GTux1Iq/lJ9iZUb452gHrv9P+v9/yiAKQAe9bffA+CCkkiIeBgXwHV2Rfe2HZk2Pa9iVj1Z8figaHadI+3gnwlCCGeh2x4j3y4+z2ygTdYmjywO3pakV0U3AKxdpSlhO7bwborY1H+hbVUqC5EwEICL4jIcw0fJAN49DuPQZZ8zoUyMvJALRkYsbJGr2xSZxNBlES0WXdhi0XVnW50WKx86ISRPCJkHYCOAFwAsBbCdUtriF1kN4GDV8cUSDN5JZZMfY7optnHdLP9Lkq6U1qiz7djFEPGhF+zDBgH1NRUVUqEgG1iSRLlYfiVlQWvwNRLfV/FBUfGFKPRNz0IP9yfJdWJ8Dvx/eR86P/U/Vp//L6+QO9Tluf0q2aJ/83XL9LTsSTBF7BRDKXtg1Q2KUkpbKaUTAAwEMAmAbIRQerUJIVcTQmYTQmZv2rQplZAyX2swKMqVu/bkEbhx6hjWLk4e3aBMhi/WI2Ngz444/GAvj8Xph/bDYQd3w9emHBIr9+3TRuGwg7sFK6aoOvVVHx8WyGdLMX3455JIlS+eMAxj+3fD+RMGROPQLVCdlyxbICHAKaMb8I1TR4blREvUEFv/808ejosnDpTu+49zxwEIV45S8fvPHoWPH9IniEFOGraonV1YBLybiv993PDeuOUSL8qDbc4TYjVLNYzsCreNPihccYm/fzdNHYOjBvfA2P5eFAxzm+Ryelccu1T1+Vxw7HWnj4rtjx8nWOh5XqHHO/mPzxuHIwf3wG1cmgbLpXIBeLpgYM+OuGzSIFx/5mhcc7I+OqytuFwS5XKhlG4nhLwC4FgAPQghdb6VPhDAWsUxtwO4HQAmTpyYSj2JU6EjcNtY/g3G/10pj/Tw5PL+1d3Ix685PljgoHunejzz9ROk5Q7p2wXPfP0EnP2bV/065fX9+zljlW0BikHRIhT6ZZMG46bHF0S2DezZKciJktTlwRf71FED8djc1QDkPmBCCO6+chLmr96O3760OHY8377qHlw2aTAuk0TbAMDhA7ujcfo0o8wnj+6Lk0f3xVPz1ijb0p191w6lSXcUmYHK6eoHrz42+B36q3NWLxZZv2tfl8fZhx+E5xasj2wfP6gHvnxSqOQClwuXR0WXbTGfI5HcOuF+O9l4l5Ms99ARA3vgiWsnAwBOGNkHry7enMhCv+GsMTF9oCM6Ppetci/Dh3eATZRLAyGkh/+7I4DTACwC8DKAC/1iVwB4qlRChpaczEJPd/HZ57520kDa+5qhy6VUPnS+PXuXCzcQxllY+vS/usE8zYs6Y8IXeHyfrv06i0Wh0xBZU1TpQmEyEEuFHrfQTWUZgcuF2OVDt7XEg+2iD50rZ04g5pXVJZ0rllJGopQziMfG/OgP4B5CSB7eC+ARSukzhJD3ADxECLkZwNsA7iqVkOJECyDsIGnvg80sxbQvi0wHRUvYG5K7XEIiuUgUyaVMJMlPUyytmvGCcoaVyTDNFK3P5yJ5SLIQV2ySf0HrBkWDsEWl4rZr7wAXbCDmoFcdW44xplJgOr8sMSp0Sul8AEdKti+D508vObJPc5kPPQ3aad8JK5dF4yShUnHotlObefHEHNo80ThoXfvqgcqs0VmdlcYUh15vaaEnQexr9fzUf+2gaOhykVu+qq4AABUnSURBVNcrb0986e850BL8Nuk7JkcSH3oxZN1FSjmYK1ITM0VlPnT2M611ZaN8U3tcMnS5lBJ+6r8NkVA1i2yBgP4rR+ZKKxU6l0ulUSt0uQ9ddQqq28h3K6a4VS/hXfubtS9a1oYqPt+2L0UVur7fl8PlUkrKOXGqJha4CN0j6k/AtHXy/fLLJw1Hn87t8dPnFinb09cZdWH814VH4O2V25TlTfuB0sbeXjxxEP7+7gZ8/vihVuX5q8GH3YkLbthaJEnj4IvhzHEH4YGZK6VRSiZunDoG+5vNK+KkRWVUsK11+RzOHT8Av35xMbbuaTLWpzIMCAFuOnsMHp2zGkN7d47se9NfZOPtlduD+6n7ckuagZAVn3ZEf2zYsR9njesf7GP9pV1dDt+TDGTeOHU0tu45EKw4BgDfmHKINB3uj84dizXb9iWSTYSliZj+ycOxYM2OVHV8Y8oh+O0/lgAor6FWIwo9HoceTixKXWukHgC4aaoXjckUelpYnRdPHISLJw5SljPtB+ITm7KkT5f2eOqrk+0P4C42H6UQiwSx7L9Ucg9KRfdO9XgyyblyfOUkc0K0YlApTmaRtssT9OjUDo98+Ticdss/lfUEi1Nr2vrMMUPwmWOG6AXSuVy4r4YksJfWwB4d8btPHxXZx4yWfz9nbCyvCwAc0rcrHr82eu+uU6xSdKVkSbikHDO8NwDg0kmDcWnKOq47YzTW7diPP89ZXV1RLtWAzJo2jbbbUkVBLlJ0q8uXG941wsdFd/NXpWFY998qdoOUE5W1yxRdHReBUg60cej+v/UW+Vt4WFWyr7dyfqmVE6I551JREwpd1sGKvfc6N07QRspGshwU1a0uX25ylhY634F1YZfFDiK3FVT+6JZWX6GzfPn+dnVooEex+iOIctG4OJNa6OFi55L2qnjAuhhsV43KtM0ytpUa2Wo9YacublBUd3TawbosB0WrSaHz58XP9ItZ6JIeLJ80Vb449FqEWejtNLNcebK6jvyiGSJi/nRbdNZqGFKaqMqqJxjMLeOoaE0odJ2vtdhOrD2+zBa6jKYS+tCTwr/g+Ax73TqKLpekg6Jt7EnOCOZDZ+4t26tU7GQ0XVrjYFHqlD502cu+rfYDdoncoKiAbAIKb5//YNqh2LGvOWGd5oucOA7dYrJSUko5KJoU/rxyOYJjh/fCW8u24ksnDMfWPU0YN6Abbn52kTJeuHeX9jhldAMG9OiI+2eslIajnjqmL0Zx+UfKxQUTBuCsw/qbCxbJf15wGFZv3Stt/+zDo+23BMrTu0Amhde3awecMroB15wcjeSx0Se3ffYofOW+ubj5gsPwgycXApC7XFr8m5vUh86qks6GbqMuF52bqVTUhEKXxqEHg6IEXzxheOI6A5dLSQZFM7TQq8rlEn2hPnT1ccHfv75kAlZv24ubn1VHCOVzXn6Xp99Zi/tnrJT60O/6/MeyFtuKWy+NzZ0rCbIoDlX7rYHytBsUZddXhe7wsw7rH+TGYQpdFofO/PrJfegeMt1Wzglm5aQSM1xr4hJKo1zY1P8i69Tmwk49uJnqMCnVpdBN+5lFwg2KSvpyTni625pllhUq5VmuqyW30KORN7Yw942sP7S20X4QPg/la7MmFLpsRmNWYYtZPh2hrz+7OqtKofO/Jeeoi7KIHIuo4m9jz3FmhGGLvvFS5uskGxQNXC4JRzBlL3tGNadlKIac5pxL1mbZWioCXd7s9NkWzVSFy6VKfegywjAt01Ru799SjDm0JcKJRdEol3JdL9mgaFqXSxiHHt/XVn3oOaGfl6XN8jWVHl1WvqKjXHT7Ug6KZkk1KXQe6ctV89DyiEultbUHOWvCXOXlRe9ySWihg7lcPjphi+yFWE4LvYYHRYu7++JCvDLKHYcucu74AbjEkBqgnPDX44Ij4ysO9unSHmeM7YerTwwHqUf164pTRjfgO5Gp2lHfolPocr556kjc++YKjO7nRf2kvU43Th2DnfubceKohkTHyVwubFO7hBb6OeP74y/vrJXm0rlx6qHYfaAlsXzVjq2BkyU1odBl0RBZqYAsLfSs+Z/LyhN5YQu7Ht8+bRR6dGoX25/PEdwuLPnXri4Xi7wQXS5tzTLLivMnHIzzJ4QvzrT9cWifznjgS8eaCwrIXiCd2nkqI+mgaLcO9ZHVmHiG9emM+7+YXL5qRxwrKgc14XLRhS2mpRTXuDaTeyan+Pw50c/vcqTPbQtkb87okb1oO7XzMhwmdbl8FKmEoVITCl3nHilWMWvj0FPekLbe1bMKLArmAtREL/zoIXO5dPRT1jZXUfRVtRJEubip/1F02RbTTnMu5VqdbRV2+Yu20IUp0c6HnoxyXS6ZAdXRt9D3lTA/fFtBF9lTsjbL11R6pEvQBaPmxdWtnViU0Bat1TUPk1LsgLQ4Jdr50O2oht7FXC77mpxCN6GLvS8VNTEoGvjQuW2hhZ4Oq0Wi07pcilRQN00dg8UbdxdXic81J4/AIQ1dMqkrK44b0RunjG7AD84ZC8BZ6NXGPV+YhGfnr5Xuu+700di8uwnTjih93pta5/LjhmDuym24/DjDgiIZUhMKXZa7vFgVoJuslLaNrN7DJ45qwJczWiVHtqRXWrKa2NKhPq/NOeKoLCeNasBJihDCg7p3wB8rlG+n1ujdpT3uveqYsrZZIy4X9RTxtG4OGx96etdCcRqv2i3WrKNSqv18qw13tRwqakKhM3KSuMViXS66p6NSD061+5Sz1r/Vfr4OR61QEwo9TN4TbtMlgkqCPttiwsoy8rlUu8GatXjOQrfjIzLm7iiCGlHo3r/yiUVpXS6lo/hVlKpbwWUtXpWfbtXhrpdDhVGhE0IGEUJeJoQsIoS8Swj5pr+9FyHkBULIYv/fnqUSUjaxqNg+/bNPHI4zx/XDUUN6KMskVay/vexITDuiP4b06lSUbNX6vJYqbrzaX2DVQt+u7XH+hAH4w+cmmgtXkC9+fBhu++zRlRbjI4lNlEsLgO9QSucSQroCmEMIeQHA5wG8RCmdTgi5EcCNAL5XCiF1SZzSfoYe0rdL5g/GYQd3x+8+fVSmdVYTLSnXk3RkQy5H8JsyraxUDCwc1VF+jE8mpXQdpXSu/3sXgEUADgZwPoB7/GL3ALigVEJSaRx6cYOi1Uy1WqziggsOh6O6SGRqEUKGAjgSwAwA/Sil6wBP6QPom7VwDF22xbY4UFStM04DC72tLf7ocLQRrJ9MQkgXAI8B+BaldGeC464mhMwmhMzetGlTGhmD5DayQdFqVX5tEWehOxzVjZVCJ4TUw1Pm91NKH/c3byCE9Pf39wewUXYspfR2SulESunEhoZ0CeyDkPFI2KJzuZSbZn/1JFkWPofDUXmMg6LE0y53AVhEKb2F2/U0gCsATPf/faokEkK+iOz0Tx2OXz3/IY4anH1wzYNfOhZ/Xbgu83prncBCL4HL5crJQ3HK6JJ57RwOI987awz6dm1faTGKwibKZTKAzwFYQAiZ52/7PjxF/ggh5CoAKwFcVBoR5XlXhjd0we8+U5qIkuNG9MZxI3qXpO5ahvnQS2Gh/+jccZnX6XAk4ZqTs8mfVEmMCp1S+hrUodGnZiuOSgbvX9kq5I7yEVro7j44HNVITYQryKb+O8pP2hXfHQ5HeagphV69cyg/GrQWvEHRUvjQHQ5H8dTEkxnGoVdUjI88La2l86E7HI7iqQ2F/hFZe/K2zx2Ni44eiMFF5oIpFc7l4nBUNzWxYpH/pd/ms8yNOagb/vui8ZUWQ0kpo1wcDkfx1IaF7v/b1i30aif0obv74HBUIzWh0Mu5arZDjfOhOxzVTU0o9CAO3VnoFaXVpc91OKqamngyg0HRmpC27dLqfOgOR1VTEyoyWILOxaFXlD987mh8+pjBGNa7c6VFcTgcEmoiyoXCzRStBkb264qffeLwSovhcDgU1JaF7nzoDofDoaQmFHq4SHSFBXE4HI4qpkYUuvevi3JxOBwONTWh0AuSRaIdDofDEaUmFLqz0B0Oh8NMTSj0gvOhOxwOh5GaUOjBEnROoTscDoeS2lDo+Gikz3U4HI5iqAmFXnA+dIfD4TBSIwrd+dAdDofDRE0odOdDdzgcDjM1otBZHLrT6A6Hw6GiNhS6/69LzuVwOBxqakKhFwouysXhcDhMGBU6IeSPhJCNhJCF3LZehJAXCCGL/X97llLIgvOhOxwOhxEbC/3/AJwlbLsRwEuU0pEAXvL/LhnM5eLS5zocDocao0KnlP4LwFZh8/kA7vF/3wPggozlEmUA4HzoDofDoSOtD70fpXQdAPj/9s1OpDhhHLrT6A6Hw6Gi5IOihJCrCSGzCSGzN23alKqOMNtihoI5HA5HGyOtQt9ACOkPAP6/G1UFKaW3U0onUkonNjQ0pGrMLRLtcDgcZtIq9KcBXOH/vgLAU9mII8dN/Xc4HA4zNmGLDwJ4E8BoQshqQshVAKYDOJ0QshjA6f7fJcfFoTscDoeaOlMBSullil2nZiyLEjaxyOlzh8PhUFMTM0XDqf9OozscDoeKmlDoBReH7nA4HEZqRKF7/7o4dIfD4VBTEwodlDr/ucPhcBioCYVeoM5/7nA4HCZqRKFTN6XI4XA4DNSEQqdwFrrD4XCYqAmFXqAUzkR3OBwOPTWh0EFdyKLD4XCYqAmFXqDUuVwcDofDQI0odOdxcTgcDhM1odCpC1t0OBwOIzWh0AtuYpHD4XAYqQmFTil10/4dDofDQG0odLgoF4fD4TBREwq94Cx0h8PhMFITCp26OHSHw+EwUhMKvUBd6lyHw+EwURMKnbrkXA6Hw2GkRhS6i0N3OBwOEzWh0L2p/5WWwuFwOKqbGlHozofucDgcJmpCoVO4maIOh8NhojYUOoVT6A6Hw2GgJhS6S5/rcDgcZopS6ISQswghHxBClhBCbsxKKBEX5eJwOBxmUit0QkgewO8ATAUwFsBlhJCxWQnG4xaJdjgcDjPFWOiTACyhlC6jlDYBeAjA+dmIFYXC+dAdDofDRDEK/WAAq7i/V/vbMoc6H7rD4XAYKUahyzQsjRUi5GpCyGxCyOxNmzalaqhQcBa6w+FwmChGoa8GMIj7eyCAtWIhSuntlNKJlNKJDQ0NqRqicBa6w+FwmChGoc8CMJIQMowQ0g7ApQCezkasKG6mqMPhcJipS3sgpbSFEPI1AH8HkAfwR0rpu5lJFm3LRbk4HA6HgdQKHQAopc8BeC4jWTTtALmamALlcDgclaMm1KQXh+5sdIfD4dBRlIVeLiYO7YVd+1sqLYbD4XBUNTWh0L96yiGVFsHhcDiqnppwuTgcDofDjFPoDofD0UZwCt3hcDjaCE6hOxwORxvBKXSHw+FoIziF7nA4HG0Ep9AdDoejjeAUusPhcLQRCKWxFOala4yQTQBWpDy8D4DNGYqTFU6uZDi5klOtsjm5klGMXEMopcb842VV6MVACJlNKZ1YaTlEnFzJcHIlp1plc3IloxxyOZeLw+FwtBGcQnc4HI42Qi0p9NsrLYACJ1cynFzJqVbZnFzJKLlcNeNDdzgcDoeeWrLQHQ6Hw6GhJhQ6IeQsQsgHhJAlhJAby9z2HwkhGwkhC7ltvQghLxBCFvv/9vS3E0LIb3055xNCjiqhXIMIIS8TQhYRQt4lhHyzGmQjhHQghMwkhLzjy/Vjf/swQsgMX66H/YXFQQhp7/+9xN8/tBRy+W3lCSFvE0KeqRaZ/PYaCSELCCHzCCGz/W3V0Md6EEIeJYS87/ez4yotFyFktH+d2H87CSHfqrRcflvf9vv8QkLIg/6zUN4+Rimt6v/gLUC9FMBwAO0AvANgbBnbPxHAUQAWctv+C8CN/u8bAfzC/302gL8CIACOBTCjhHL1B3CU/7srgA8BjK20bH79Xfzf9QBm+O09AuBSf/ttAK7xf18L4Db/96UAHi7hNbsOwAMAnvH/rrhMfhuNAPoI26qhj90D4Iv+73YAelSDXJx8eQDrAQyptFwADgawHEBHrm99vtx9rKQXPKMLdRyAv3N/3wTgpjLLMBRRhf4BgP7+7/4APvB//wHAZbJyZZDxKQCnV5NsADoBmAvgGHgTKurEewrg7wCO83/X+eVICWQZCOAlAFMAPOM/4BWViZOtEXGFXtH7CKCbr6BINcklyHIGgNerQS54Cn0VgF5+n3kGwJnl7mO14HJhF4qx2t9WSfpRStcBgP9vX397RWT1P9eOhGcNV1w237UxD8BGAC/A+8LaTillC8PybQdy+ft3AOhdArFuBXADgIL/d+8qkIlBATxPCJlDCLna31bp+zgcwCYAd/tuqjsJIZ2rQC6eSwE86P+uqFyU0jUAfglgJYB18PrMHJS5j9WCQieSbdUamlN2WQkhXQA8BuBblNKduqKSbSWRjVLaSimdAM8qngTgUE3bJZeLEHIOgI2U0jn85krKJDCZUnoUgKkAvkoIOVFTtlyy1cFzNf6eUnokgD3wXBmVlstrzPNFnwfgz6aikm2Zy+X77M8HMAzAAACd4d1PVdslkasWFPpqAIO4vwcCWFshWRgbCCH9AcD/d6O/vayyEkLq4Snz+ymlj1eTbABAKd0O4BV4vssehBC2KDnfdiCXv787gK0ZizIZwHmEkEYAD8Fzu9xaYZkCKKVr/X83AngC3kuw0vdxNYDVlNIZ/t+PwlPwlZaLMRXAXErpBv/vSst1GoDllNJNlNJmAI8DOB5l7mO1oNBnARjpjxa3g/eZ9XSFZXoawBX+7yvg+a/Z9sv9kfVjAexgn4FZQwghAO4CsIhSeku1yEYIaSCE9PB/d4TX0RcBeBnAhQq5mLwXAvgH9R2LWUEpvYlSOpBSOhRe//kHpfQzlZSJQQjpTAjpyn7D8wsvRIXvI6V0PYBVhJDR/qZTAbxXabk4LkPobmHtV1KulQCOJYR08p9Ndr3K28dKOWiR4YDD2fCiOJYC+Lcyt/0gPJ9YM7y36lXwfF0vAVjs/9vLL0sA/M6XcwGAiSWU6+PwPtHmA5jn/3d2pWUDcASAt325FgL4ob99OICZAJbA+0xu72/v4P+9xN8/vMT382SEUS4Vl8mX4R3/v3dZ/670ffTbmgBgtn8vnwTQs0rk6gRgC4Du3LZqkOvHAN73+/29ANqXu4+5maIOh8PRRqgFl4vD4XA4LHAK3eFwONoITqE7HA5HG8EpdIfD4WgjOIXucDgcbQSn0B0Oh6ON4BS6w+FwtBGcQnc4HI42wv8HwgYgzrqYIXAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x18620ade160>]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztfXncHUWV9lP3fbOTkARCCAkQkKhE1hAhgKCiuKHCOOg4jsrMMOKoM26zxfGbcQeEcZ1RERdEhwFUVJBEMIR9CwTCnkBCCCEsyUv2/V1ufX90173V1bWcqq6+S+zn90vue/tWV52urjp16qlTpxjnHBUqVKhQoftRa7cAFSpUqFAhDiqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CCqFXqFChQp7CHpbWdi+++7Lp0+f3soiK1SoUKHr8cADD7zMOZ/kStdShT59+nQsXry4lUVWqFChQteDMfYsJV1FuVSoUKHCHoJKoVeoUKHCHoJKoVeoUKHCHoJKoVeoUKHCHoJKoVeoUKHCHoJKoVeoUKHCHoJKoVeoUKHCHoJKoVeoUOFPAhu292PeIy+2W4xSUSn0ChUq/Eng41c8gE/834N4afOudotSGiqFXqFChT8JvLApUeS7BobaLEl5qBS6grVbdmFl37Z2ixEdG7f348mXtkbNs2/rbqxYl82zXudYtHJ9cJ6cF7s/Fhav2oCBoXru+iNrNmFH/2Dj+9N927Bua3dbfC9s2onV63c40y1auR71Os9df3D1Rry0eZexfQ0O1XH/qg3O/O97ZoM2fyp29g/hoec2Nb4/9vxmbN010PjeU2MAgCEeXkano1LoCk44fyFO+8Zt7RYjOs747h1467dvj5rnKRfdjDd/M5vnFYuexV9cei9ufPyloDx/fk9y/4In1sYQMQiPrNmEsy+5B99c8FTm+uYdA3j3/9yFT1/1UOPam75xG+acv7DVIkbFSRfejFMvvsWa5ran+vAXl96LH9+5MnN99fodeM/378acCxYa29e3b1qO915yD5as3mjM/5Zl6/C+H96Dy+5e5S2/wKevXoKzvncXNu3ox8BQHe/87zvx0V880Pi9odALDBqdjkqh/4nghRJ4w10DeQv2mZcTS49i8enwzMvbAQBrNobdHwNrt+wGACxfm7U4d6ZTddkKBIA9WD808MKmnQCAlX3bM9e3SBawCY+/sBlAsihpwtPprPj5jTtDRcTdK5KZHefAjv7kXT26ZnPj9x6WKPT+wXy73VPwJ6fQ123dpZ1KA81G2wl4afMu8HRquHbLroZVsX33IDbvGMj83kno7Uk6zWBBLRfyaJt3DmD77kF3wlxZHC9u3omhOsfaLbtQTwtnqQIQEF9dom3Y3t/gaQeG6lZKZuP2/gaFs3lHmPwA8PK23dg92HpuWLxvG0Rb6O3Rq5ttuwfRtzUZREcOC1dJW9O6e2nLLuxO63+ElF8ttdB3eyj0TTv6MxSbL+T32wr8SSn0bbsHcfzXFuJLv38899uDqzfipAtvboNUeTzdtw1zLliIS29fiXVbd+GE8xfi4hufBAC8/uJbcfSX/9j4vdMgprX1wMGGqjR1OPpLf8QJAfTHz+95FidecDP+6sf34oTzF2LTjsSSrCm6Snx1PdqsryzAey+5BwDwxesex/FfW4id/Xple+xXFuCd370zkf/Lfwymb2Z/9Sb8vUQvtAq9aiVpIAwoU9ojvnAjfpi25VHDegrL9Pbv3IHbnuoDAIyU8uttKHT6wHfMlxfgHd+5I1iWY7+yAG/7dvj9vviTUugb0ynfLcv6cr89okyj24k16bTzzhUvY+P2ZEp787KEU3552+5GujtXvNx64RwYlnYa0yzIBZaqzdDZx7YAC/eep5Op+r0rk4U7MV2vMYNGJww3jz6fTPXvSt/R85bZ38qXmzTG1gD5RV3f8mS+XceC6XX01LIqRLeoOTiUWugE5T9qeHGFDgBL0v4sK3RhoftSLqsC6UOB1RtaRx/+SSl0MQ2WG41QHNsNFlQroCovwfUN1TnELFW3kGNa3GknFSM6uE42ilyqDm3HswjZcwrdAtNAMmXvUQD0dF6RZ5Pv3borKZugL3N5+MqgVkmPcmGgnleWA2l9Dut1q5uRESx0ANiRvg+ZwunVUC6dSFsWAUmhM8bGM8Z+zRhbxhhbyhg7kTE2kTG2gDG2PP2cULawvjj3Z/dj+tx5je9iUUu85L//xQM45HPzATQ7hQvvveTuTJ4q7lz+MqbPndf4t3ZLwp+uWLcV0+fOw02p98avFj/XSPOXP7o3k0dNUuKCx9Xp7rufXo91W/L87CGfm49/+L8HG9/P+t5dpGeT8ctUvrWa/G0QnOoVi1Zj+tx5GfmO/OIfnbII9fDVeUsxfe48HPK5+bi7xJnItxY8hRsUj5yvzlsKoGnRCYi+r9MBR3zhxgYPLGPK+JEAgBc35xX6jkAj4ovXPd5otwCwZWcyixszvHkAGee80b7ue2YDps+dh7nXPJL5/ZDPzcf585cGydDIR5mtCGscAM756X2YPnceBtMZhKr8dRg1rAfPbdiB6XPn4dqHnnemv/APyzB97rycASEMtJG9zQFClP/RXzzQqJtDPjcfD3fQ7LwoqBb6dwDcwDl/NYCjASwFMBfAQs75DAAL0+8dhYXL1mW+C68MwdPJHZm6GHX/KrPrFQAseCKrHMTq/UPPJVPw+Y8mW49/u6TZWMVUX6CHNXlo2VrXwTSdu17a4qx6ZVBwzQNrMvJTIfqs8GiQ5du2e9Api67PX/Ogu2OH4hf3mk/2Utf7xLqAyaYT3LuMsSMSJbttd155U40IFT9LXfsEFyzyGT2iqbzk9rIwpeuuuv+5xjXx84/ueMarbHUwU7/LCl3w2OIaxRYeOawHT6XeRdc+9IIz/Q9vfxpAvn+Ihcgs5aLP45HnN+t/6EI4FTpjbByAUwH8BAA45/2c800AzgRweZrscgBnlSVkLDQt9Oy0rl7nQdyrDqreFQ1eUCdCKQwzrPgnaZtWuct3VvXECOWujTJ4Zjc0lJVTtXJdUJ8H0Fu3sTBmhHmKr1Iu9YaFbqC6NL+JtqbbnUhx+dPLlXwKRS7ykS10efPMwGBe3tBFaxVqLnrKpU4us6fWHNQp6RuzJkWS7buzs3EA6DVo9KnpLGpPAMVCPxRAH4DLGGNLGGM/ZoyNATCZc/4iAKSf+5UopxOnf/M2fO+WFXjihS2YPncenl2f9Zc9+wd345yf3gcgv5K+dfdgsLuYwNotuzB97jw8qGyeUPnY3z30As763l0N60VA0BNPvLAFZ6ceEkN1jj//wd0AzI1b1jn9g3XM+PwfSPKe9/PF+NRVSxrfT7noZvzwtqcb34VC/+BPFuH+VcmUXbjfbdjen6Gd5N/VAUXMMHRK8Om+bZg+d15jt+mHfrJI67njG3vjsruewQnn30RKKytBFavWb8f0ufMaFqNY8DOpmbd863Y88Gz2/QsKamf/EK5/5AXM+Px87BoYwm1P9eEt39JvxHlf+v5NEO1XKHSxG3L0iB78+28fxfsuuSdjOfcPNQeTR9Zswvptu7Xt5Lan+nDUF2+0buq66v7nMlSI2i6XrN6E6XPnZRaBGxY6cRARg/rytc328cPbnsbp3zRv+FOzFhb6CM2iqIqRvT34zk3Lcdp/3WrM/7B/n99wqnjTN27NbToDgE9dtQTT587Da/7zBvzsruzM5/D/uAEr1pW/A52i0HsBzALwA875sQC2w4NeYYydxxhbzBhb3NdX3ir88nXbcPGNT+JXDyTTypuWZumWxVJHUy30rbsGCvtN35tuV3/8hS2Z66LBy9aeiXa4f9VG/DqlOsS961Je1mShy/nqpvwm/PGJtZkp7XMbduKCPyxrfO+RGv+P70iU7INpHT6oKK0f3d78vV+10FP5dHzxbx5MnvUPjyY01R3L9Vy5j98wAHzp909g7ZbdJAUyZoRZoT+4OnlPv0vpMRuHLvDtm5Znvou0O/qHcMH8ZRgY4ujbuhvfv2WFMY/7HNvkR6WDkODORZ331Gr4v0Wrcd+qDZn2IlvoV9y7uvFcKpav3YotuwYbHlUmXHTDk42/1boQbUFe9xAcOqWLcd5cRxGDwvWPvIgL/rAMy9dtM3qoqAOLGERk3t40V+QAvnXTUxlvIxWDdY4VKf34dN92fHfh8lwa0Z+29w/hN0uyNOHOgSEMI/jsFwVFoa8BsIZzvij9/mskCn4tY2wKAKSf63Q3c84v5ZzP5pzPnjRpUgyZrRDv1VZ1wsIZk3q7bNk5WHgKalqd1yl0E3pqLLNRQ+6UJvlko0PuMD0eVIdO8ckdYVBSGEB+05CQs7dWy1noYpar44u37EyujR1pVqqAnlengLLFezTBTU6lAGwDhbqmIUTYOTDUrKcehuEEjw8ThMyiTrlGLrm9qO/E1daf35SfEcmURlYxZfPakBoVY0cOa5afPjc1TotKu8lt0TRbUx9J94ym0qldf8LoYe5ElvJH9Mbx4LHB2ao45y8BeI4x9qr00psAPAHgOgDnpNfOAXBtKRIquPr+1fj6DYklOf/RF3Or9KJRf/n6J3LXBBY9sx7/eOWSRqPbsmsgaGeiDKNCT/sSRb/21lhG8csKiWKhy7ypzwYN3exEHhDE4rJw+xpUeFLRYXt6WM6CumD+Mlx0wzKc8d3m5ooL05mAoArkzq+DeMYlqzfib392f8Pic6Fv226c89P7tJ4637tlBX5+zyrsZbHQVcj1u+ylLTj3Z/fn0sgK/UM/WdRYuNw1MNSo554awwiHQv/s1Q/lrj24eiP+5rL7GoPBX192H+6UZjWyEpFfUb9SX66Zy/NS2IV7nl6Pt337dnz+t481rsltQ81K0BLjRjXr1ctCBzBUV40C1hjE1mzagTuW9+Hffv1IJo2qQMW7kq+anvuZ9WbLXIVtET1TvqaJut55DFBL+EcAVzDGHgFwDIDzAVwI4HTG2HIAp6ffS8e/XfMofnBrwvV+/IoHcentKxXLJH+Pem3V+h34/cMvNBZMdg/Wg3YmUtCw0AkavaeHZTZfyO3PvCgqlSWlsVnoasPWLaTq7jct0IoOOExjod+54mV8/9ansV6K43FJytULC5NqoX/yqiW4edk66yYdGb9avAa3PdWH/7k5T29cfOOT+M9rHzduR9dBrrd/u+bRnBeVijuWv4xH0lgiO/uHGvVUY24LXZ2yA8Anr1yCW57sa7hHDtY5PviTRY12Ilf9kMVCdxkvMj32d5ffj2VKFEV5cVFtluI9ywaFL4feryzi1lizX+zYPYQP/eQ+XL34ucyWelUOMSbIZZqK/7pENdrAOfAfv3vMnRD6Zx1RIKwBFSTzhHP+EIDZmp/eFFecMMhBf9TVbsDs+VFreHLkN1hwzrUeFyYMGLm95JPig9tbYxnFP+QYqIDmzkqAZtED+U1UOi8InUIXdNDAkDogNC1P1Rq0QSh012YSVRJmJdSaEFPkF5VpuomacKHh5QKQZwki/x0DQ5JiA4Z7DCQCjbAIJi8bI+WitG1HOXL2w3trufYi04K6/qai6eXiTArOgX7FQu+pNV2OORKqdHv/UCaQl67/qqDQlrGgKyvknfuia3aKPvzcptzKsYA8zdU1GpNyE0pWRzmo9yxZvdFY/sKla/F3P1+s/U28WAqn3VNjGcUvW9zbdg/iMk35V963GkCy0Cr8k0X6C/+wLKd4bnlyXSZGNJCfkpvkFdTHz+7OyiHK+PoNy/DyNtrC7PxHX2ws/i1cuha/f9jsc1xreMok3xkD/vDoi84QvcL6vknJXzYAfOJvNzopz/pb2yBu2dXfpFw4eBCHLuphi7IeIZSq3Gavl55XpcFMik13WSenaRYpQ+y3AJp1Jcp97PnN+NHtK7GzfwhfkahRIKHhvnhdNtZSjTUNhW/f9BQm7jUcAPDrB2UHgmz5upjnJscHqmutevdFNyxrej6plI9SVo2Zg5PFBJ1AbDPOVHYYyhW4UfLu0I3Mps4nUwg6Dk6unD/7/t1G2c69XK/MgWYjphj7vbVadlFUkelLv39CvQW/uPdZfOWsI7Q7MC+57WmccMjEzLW/uex+3PLPb8hc0yp0jcCi8T72fNaTR3QUnw1MH7+iuZP18nuexeX3WLjJhmXavPSx9P5VF55hvE0OiPWPVy7Bu44+AEASlVHAy0IX03jQ/f1Fe9w12FwUBQ9bIDN6aaTZys/yRamtqArdNRjJVrdWoffIlIs+L3nDkmgfIu07/zsJRrZ7cAg/uTNrHHxzwVO5MLuycfH4C1sazgw/vK3p4qrKIc+mGnIY3hll9qzD9299Gm8/YgqOnLa3Zp0im9bkAx8bXWOhq5ArUK5LdXoJ5BfxBITFM1Tn+VXySCGTG5QL1UKv6SmUUOjiiquNX0cX6eQ1WThUazUUoTPiXYaoerK0PvuwZC8X3QYaW1n9g/VGO+TQK0oXTJ5SLnfK3XJf4dzo+icUuZyPbuChWOg2OQV0u2cp6zm6uEs5he4R+4hsoWtuF3pada1V5aGEGY6BrlXou/plhd78W/bjFjApIvlIKvVl/e+9z+Kq+1aTF+BMuPzuVRgYqpPcFnsVyiXGZr41GvnvknyEX9q8S1tnOoVuojhMA2YsmOKSA8lmmKUvboEOuwixUnyCM8lJqYOY6NgDQ3Vpp2leoW+0HP7QgKEJCUlM1rIcq2TnwJAzfKyci473DVVOfVt341eLn2tw1r9c/Jz9hhSkIGnKozcW4aXrOmNPyCXw2PObcdeKl3HFomdzmw1vWpr3zxf1ow6SQ5zjOon28nEjLoKuoVxU7JS2Urs2Bbm4s3qd5xZ3vpa6Q44b2YtHvvjWYDkfeHYjLr19JWYd5I5dVoaFvkHDaf/ntU2O8m9+dr9WIeoa4BWLVuNrf3Zk7nrRTVkuCEmE8pUVvNj9q6NeTMGvZL3ntygqOHA6RPayMuHgjTDDAl+4Lh+jX4WZchGUhlue/sG6caOWztLXzSTkkLk+Rse/XpN1NbSdYJQtz60MKc9OockEHQQgF7RLnEkgQ4TwUOt090Adn7yyuRObEjo4BrrWQpddllyKz8ydpb/XubFBqAtQIdi4vZ+8kt5j8HIJhSsPOb66SQ4XyqZchIXmO27sJJzuTslTeNU0KRd6mFeRf38mZKtmEU+ziJaTw0S5NMqiVZB7562dQ5cHo1gxYWygcNwUOSjtQcZ6wgK/6Ce7lbxVY6IVC6JAVyv0ZoW5XqbRQmdNC91mdi02bMVeu2VX47xEV/lUnk5OpwvH6gtXOzftlDTztfkMyw7gT4mRrguBend6cIXA+sbg1byf8v4EmgttnKzQhTEhW4cc+TY5eVw2QFSdA6uUrejGJuTg0DNJuVmhy7ev7NuGZ9dv126GyWwscheZg81a1Q1aFMbF6LkjSWg6NcqEe1audydKcd8zWR2hBvtTZ2RloWsVuhy9zmUhuhZDBjWUi4yzDcGSTrxgIc747p3a32QM1utkCz321Mw1ezHtKDUpdB8/81gQnVy3aAckVqTqBQU0D5wW0HkqUdwsm+U2y6fu+hOKW643zvNeVZPGjsjd+wYlWJTJ/17IR7OWzYuiTfmA075xG15/8a1ODj3kgAhfDp4yW6SI4Wuh+8Srn/ubR62/91SLonbInKSrIRs3FqV1rFsUpYBKAfj4LIe6UJngqhtXDBoVlEBZsQelJodeLB8xk/DNp+lHncoBenyZAZ2FzvMDLWVGbiqzuVOU2s4M1qz0fAK63Y22naIU2EJH62BT6CLevKudc869FXpMUDfDFUXXKnS58boW5VxxUFa9vL0Rj6QMDAx5DBgB713nmijgKtdkoRsV+oBboRcJPKWDenC0Kpuvlej7pgeGOF7avAvbpXUbapFCkee8IHLhE9x5uTl0dx4Jf0+vAbeXi3+/8VXoNi8XES3T9EicA6vX78D2/qEoXmO6/Cmb0yLbaUZ0rZeL7CrnskxM7kpi5Fc3N1Dgo0TUYEPGPOF/xiMAvO7rt1jKtudn4tBN91FOTB85rCf4eDUd1J2iahXF8AayYbBex5wLFjYveBSnWveAXqlSlKxrYxG17ZgGD5VaAvSHNof6oQv4uvDZ0ov2a6q/Z17ejlMvvgV/fdJ0rzKp4KDtSWiRPt8zLHRXh6ae9uMDn4Y8UKdEvEjyjG1FuLxcRhoUehHKhRKS1gfN16Tnin31uW8d52Oh0N4noF9z4OC5Nkmh5YyUiweHzi3pdJTLxDHDc+myO0WdRebgE24BsFvo8uleOojAYktW24+OLALRPmyRO4voGh90rUJ/SQqJGuq2WITq9XEpVI9lsyG2senbeQSMFjqBcomt0LfuGszMDMq1x/NQ20+dm4Ox5e/NS6uz0Cmbs1wcOvVVuyx5+Wcd95u10P3fhu8xiTYLvTl7s8sxQTMwxQDnzfax/97mo+wqC90BOT6zU6E7gnOFwIeHHKzXaS5lgIftR4NLznnSYdIy922qUoqXy2jLsW4heObl7fjAjxZ5Uwsm+Nax7kCPJwy7U1XolBdHvs3+tybErwqjK2n6STuD07zngoqibou2/qp7QpvhVXNY6ALjR9EPp/ABR9NC910bKAPtlyACXNayyfppFeUyqAnPe+ikMdq0sS10rz0/UlqTZa877FhF7EVRINlxKyRqwV6WDIoMIDpjgnMetLvW2Fold0oKXP0lE0Nco7KLcugxHRCEKK7BrKek4FjJu0zPA2iRa6INe4ZCd2gtEz9ZhNbystA15U/fZ0xO8XEetihqg1esEqnzmjo9RaGXtYfCZ4u7PR/P9MWK0+YXQoXF2CnKYXNb5Jn8BFQ3VHmmENJefRexbUWoC+btgDhPwGqhV14udITGcikCnyzvXPEypu87OnOtp6b3TI2+KOoV7zv5/N4tKzIHSMu47K5VznzK9rktQktNnzsP5516qF95kd/J925egd8Z6leHT165BC9u3hnFDz1xs7PnI1ev7tl/dMdK/DA9DPqSD85ylqnCV6HbBqom5eKYdZS08nLaN27Dwn96PQC7hV5x6B5wvUxjuNAC79i3Uf7vvasz33X8vc0DIRSU/PYa0YuTD9unYW3pghD5oKzQzw1LtOBm1UtvX+lOJCFGTB0ZuuPlbLju4Rdw/6qNTqVAbTumdLoY4kB+Jis3/ec2FItGSoHtqcTkwfnoJVrwgwQOvfJy8YBLuZq3Ooe/ZdtUcyTh7MCeHtboKNMmjGpcNz3K7IPd0Rp1oIw7R03bG7MOmkCqDZq/dDmNt+lWV6x3+vpB+3gp6RCLWzWJ4UNFcZgHKDneu5zehtEjynJRbcLW14Rh1IogYSaIhW/bEXOVhe4Bl0JXw3YK3LTUfsivDbYihxFMVNlCZ5KVYWq8pJjQGlC42mE9NTDQZiwkhV4yh16079r8hXUoqiym7D3KnYiAvi27tNeFdJRZ4+1P9eGO5S9rf9O1vZueWGsdoEPbpQ9sj1UjKvQy1b04G9cWo6ZVO0X3CIVedjxuHWwNiBJ8KOHQk3RyhzFlG0pjUJTRsB5GbnEUuqPGGD5wwkGk/HwQy8vFt3MVVejDe2sYPzrMbS571GLz6Ly75p4GANh3r+He9WGK4qmjXJav2+aQz6/sEFgXRdN+4VIBZVrwIq67lXKpYrnQUfbWbx1sli8l9rE87Zd3QpoaXuiJJxS2oLdWyx0iYc6PZqGffvhkgnSeSIsu2jnVk2hcKNq8agw4Mz3P1BfyLlURXGrC6GGYOn5UY9CM1fqHTDMgS9OLvb6gg3VRtNF5HBZ6iWJu2J4MkFbKpbLQ6WiLQrcUaXuxAj2syaFTQh2VSbn0Sny+q+FT8qsxVipp+OOA2DsyTLF9TCjavrbvHgpeFLPtIKUOwlQ0T2Si5xfbzVbXcOyLorSNRWVqiA3bk5lT12wsYoytYow9yhh7iDG2OL02kTG2gDG2PP0MW7WLgFZYCSpsVgNlEaxWazZdJvnSmvINVegUZSTTP67URYJIFYUo+fcP013+QiGHLyg6I3CdS2t7tQNDPLfILg8OMZu+abC2vU9R/uyDJ+CYA8dHkCIvg23QkI+RtOaa5vGeWVPx6TfPyP3+ydMO8xEyg227E4Vu59A7j3J5I+f8GM757PT7XAALOeczACxMv7cFofFKCpVp5dDN1XrQROGPLlEu0u+mbEMpF6oCbrp/6dPP2G8vHDVtb9qRbawcxjDEGjSFB3Zh7Mjmwqmvhe7r1WKLHz84VM+FUmgaAslnLKUuH2KdKc/yOKJ9feFdr8Fpr94vjiAK7BuLhBy0PP7s2Kn49JtfieMUr7EZk8cGy7enbP0/E8Dl6d+XAziruDhhaMeiqK2R2Y/Yan4yhXP50R0rjfxuqIVO6exMon9sVclA93IpwyIJecuhu1bHjAhX6DpKx2WFmzBY5xjZq1rozb/Xb+/HYx7H6NkQcgi2qBrG6Kc42ZGvqJ/dvcqYuqcxu7VLLZ6tMRNV0hfx1hHxjcoIeeELqgQcwB8ZYw8wxs5Lr03mnL8IAOmndnhmjJ3HGFvMGFvc19dXXGIN2sGh28q0r3abr934+Fp8/9antfeFKibSIibyx7zpEzKSQq+VS6F7wXdgOfag8Rg9vAfvlhYx27qtfKhuHGSFcjLt6vWFbqeoXI7+nmbisizUh5TzYmW3U0bl0KWBB8gPWkXCVQi3xTGWoHSd5od+Mud8FoC3A/gEY+xUagGc80s557M557MnTZoUJKQL7VkUDePQRQNkyHPoZYDE/0rFu5JTIp/WJIs/JkyyyfSILw6cmPURf9XksXjiy2/DAZLvuGlQfMvMyblQwdf/4+u0M7RQEopz3Qw0bUOR61j0I59F0Ybly8JpQV88+sW3NKi0BlXo8nJRfs/TSuGyr9m4Az01hv33zp8L28w/OHsvkBQ65/yF9HMdgN8COB7AWsbYFABIP8N36RREOxT6QsumJIqFLr/gMt816Xgs0BQwo+ZX1sYiQ6ctMl02Klrpsk/7iv3sdSmaX9kQyll3oIf5nuSzxlhpQdl0EGW5DrgQaFBDjvxCsGT1Juw/bmTjrFVdXh2j0BljYxhjY8XfAN7Gpyt3AAAgAElEQVQC4DEA1wE4J012DoBryxLShXYo9K/NX2r8jRx1Lf27zN12tEVMmVu0p6Nx6CzYIrUh9oIxkO9ovlvPc/kZnjv0FQ/VuTH0QOwaNr3bXZZDTWQLvRZBo5MMC8YafYa8U1T8nOb/wTkHZ34v2gflNqhziuikjUWTAdzJGHsYwH0A5nHObwBwIYDTGWPLAZyefm8LKEpm6vg4268poJywIr9gWiMOk4Uy2DEpf9fUleo1UwrlYrhepDPm72XS/wl83GIZ08sZKmFdQ7mUNf6HTAQa3DRYqYaJClGUGETcg252UfTs46ZhwWeazHHRgHJ///pXNGTSBd7rmEOiOecrARytub4ewJvKEMoXpIW6Fi5A2xS6xkAvFRTrUt4HZLXQ0dzmbEOtJLdFk0YvshZnklPmVH1OTIvdcbmFcom99hLiby8ouBordgKYL4QiF11t6YtbrenF2aIyMrRnYQvd/nunLYp2NEi0Qgv9LmyWiuy22LxWnmw0Lxdpp6gpTeqKKMcTMebXYheXQoqEQLmY1g30lnhkyoWHnW4UWpYvZLfFVr53UZR4966Qz2s2Jpu7sjI2v8SYXdiyOH1mCaEwNNgjFDqJViC8r6+cdQRu+5c3FJbH6ocuTellj5eyQLEuGaMdtkuV07b1/4yjphBzycO4KEqcEfn8XoRy0V8Pe8uDQ7xlbpMh5TTfCTPOTMtQ9A3q0jPzMhZFk3zNFOr7Zk/Dx98QvhPVB3uEQo+1HX3M8J4oBxxbKReW/5vSJkM7NZVyEYhhDDKYLdUih/WaNuDYrCvX45gUQrZOfBZF42K3MZZ/fIQ4F8gWulGhe+RHTSveWxFFLL/jwhY6M7f5CaOHR1kwpmCPUOjU3ZDuNHF8aSl5yPKUuShKm66zpjyG5JzTZZB3nrYCRd4Z5dYYoSVCJTQfzlKOi6QvhMGQ7D1oJa2ZfPq+e8oAHh0t7At7hEKnenJQEGNhx26hS4q88VneGx8gcC6ZRVGLTUuV07ZTtJzpd/i9FM7bvFMz/5yM+cec+fKZrzH+pnt/ZbWbkFlgc0u9ue8wxnDmMWHhg00IplzkGbImv1AwlDwoELFHKPRYp+gwsCjeMHa3RU25JTYEm4XXKF+SwVSVPjJS08aahRbaWGTivKXubjIYvLzDLSLaTjSivL9YCPJykSgX0/tkSIK7xUTDRdCzv5o59PI6YUsdMlpWUomgbnZxp4lEuRC9XFqxKEqhXLIWui0hrUzb9Du7eBTnycvecu7FoTP9c9k6tU0p9VtmWLF1kK2tmGL81yXKxbo4HVlYdWNRCGSZioahaSXdZENXKvRf/f2JmchuNLdFGmKM1Havixab6AQwiUOP4eVCdWHrDAvdTbl4bf2HP+Vik7+VkURtC7C3GLy/5Ec1Uy7xLeDGxqJIlEtRhcwMf6tllo2uVOivnT4RBzbiisdzW2TM7HrlA1L4XEm1tzIGhg6yAo6hPmzx0LMdqhwL/YC9R5LvpUjgZ6H7L7rZ2pxtQZZaeycftg8p3c5+89F8pp3WYvBizLx5L9lFShKBjHALXZ8+jh96ahQVzikcXafQdfVOcs0jNP/hPbXyF0XFp8FSMN9XJscnUS6OWC46qM/b6gUiH5cwVS6KnD5Gcshj29pLK+MU7egf8r6nyaHbt/5TFabv+ovvQKFzGw7JR5ev0YgplrUXuk6h60ANQGXD4VPG4S0zJ0fxF1Ubb8ZC0rgrtpt/YxLnHeLlog6C1p2yxi95TB5nDkeaLZ+UDAAwTDEjKRuBYkRbtInYa1mJ121q8tm/kJRNS7jTodB//OHZuWtNDt3y3okUnA/E+4m1flLUYLIZaBXlYgFTPgEq5WKv1U+88RXRnP/Voj584vTmb8qn+rcPjj0oxhmOqQyCcgkwCFV9lFSjiXqQBjRHvkdO3RsTxwx3lu/TqdVzH41ui9LffhuLQigX82+6ts0lqzgmXBb6m2dOxt7KxjA5OBclKF1sFNkpGnvW2wnrol2n0HWIsVM05stV23WP1ipvShXaEGJJnEwXk9w27jAH3zLJqXZW6sYiSien0Gk+yiLEohMn0lAQ8i5t8segXKgy7bBw6CbI4XNts5Oy1om83RYD1jhI+Vp6Y+W2SID8Amg7RYv97gO10chKJGMhCOVOeOHvPuYAvPe4aThL2qARy0KTg3N9849PmdORFTpxXSBSnRdR6GbKhZYfVd3a3RbNv9mMldhqwhb3vFGmUmjmxCKbl0tJGt3by0X+O2qfz77ji84+ChNGh4e5CEXXKXSdEqMET2pprGa1bKkxa3e4EUQbOayGi997NL79/mN9biNBVsAh9qDaV2U3SGu5hLwp+ag0ig2udiBTCCHoFAv9PbOmet9j83kXUCWtEykXqvHhW+8x+0AsMMbwvtkH4sxjpkbP24WuU+g6xDgWLWad5yx0DeUilxladqyGwqS87PHQ9QWq1hfVGHMPsiw65ZIffOJDG1bXUpBtQLIqdEuepo1AscEzlItpPaLE4+kK+KHHRogHVWx0rUKXlUusaIuxoL5IedFQDp9rSq/NU/MEsbi57HTRMsU3FKdOtamUS6yX4tdx7Il9PUjy95uUmhm2fQsh2/GB8nbPqlEOmzMaS7RFBrxl5v6YOn5UtLjgoiRvt0XpTWSNq2L1ZZuBtFL3dJ1C11UOLaAgrSO7sP8496YVNaseg1nuw6GTCiqCAnnlLfTWRlv0KSqnABRBZQVVtiwCvjtFKe2mrPoX2TYPZ5a2/lvKnDR2BO6aexpevf/YKHKIWimyU5Ry3Svv3PfWm+hdp9AF5BcgpqX2gyWcOXqXS01jWhT1yVN3Yxn6PMhtMd+SSY3ZnYImjGod2e4qO+6L8V1aXnLoTlGrHBlLNP4zNw9nFmWYlat8NXb9F8nNN4S1rxyNPFto3XSdQtfVTSMmM/Fgidh45eQkkpyw3tWGLcslW1cN+iV4eh92nwrOaaenmxRD7nmJFjpF0dCiudOhGXuy34tSLgS/dhX2naL5a//wxuT0G5uMLvlVf3IqRL6qhW47JFp+zzajy0uO9LNVB0dQ0A7OXEXXKXQdhJWgvtuvnHVE42/X1Iz6MnTJTpkxCasuPAMHjB+pTZNZFNUo8VBrNuaUTsgTYhCqUlD7WLRFXY981EFEvbc5noUJF/JMVoWuDLCjh/fgQ+lGNWpRunTvOnoKVl14Bob3+qqAJDfRphsUVY1mffe08rR2DUzOZXF2iiptS1NO2dgjFPqQxOOZEKtSdVZlzqi1+KFnFkg9rEFduTH7Bm1h1nSv2pCJ9FWkND5vt2x9wqCnrWz1a2u3KuVSRIlnfy/WI2o5C928ySdLuTjk8hQrViz84pRLB5jn8FDojLEextgSxtj16fdDGGOLGGPLGWNXM8bce7QjQFdxotHnXq7Us1rptqgaKjWdhS6XHdiaolroaV4mA13efORCjYUPUiooEwYfV8TcJihDuuAOrrnvyKl722+xlGULn2unXOz8cFGar0m5NMszvk/pcmwL3fc5QiixorJ0qh/6pwAslb5/HcC3OOczAGwEcG5MwVyQGw/NbbG8Ws35nypl9Wg49CRd9tNaBqHcUDBJAZv8vuWgXaqnQk4Oxkj1HYv+9KkHctS/UFk0d/7m4ydZ68N6yHX6PkQbansgt/RTyMwlC50iWSwOXRRWKNoi8S1/7wOzPPPNX2sVSAqdMTYNwBkAfpx+ZwBOA/DrNMnlAM4qQ0DAHc+jzvXnWGbuclno1I5Osjyz33XTzESJMm36ToUQc+zI3uz13IyE+kxxHtxn2p0bfAtUPpVaGdZTC6ZchAeXjp8u4rZYtOZFm5bdFil1Gc3LpbFu5pdfiBU9ZkSPO19YZnstpGOoFvq3AfwrALHmvg+ATZxzEc1nDQD/vcZEuAzwep07X6yrHZG5SUJCK+UilLjsVhZYbkxrTeRF8XJxceZkDj2S+MccmI06aWsu9LjcoTSYP2wx+AcVl1yf/Pfdy8yCBj9fwzJO20u9eb2lFnoEZKsgbAbVaXAqdMbYOwGs45w/IF/WJNX2I8bYeYyxxYyxxX19fYFi2lHnnGCRlP9SRAVQg3M1rrVQeRBoTmd5Oc462EKPg4+ccij++qTppLRUvj2YcjG57tnusfTCukK5yBlZOXQwLPynN+DuuadZSvaH6EcNygX0dSrAbaGT610ZWOi3Gd6PJRvKrEKecfvkHRu97iQ4GcC7GWPvADASwDgkFvt4xlhvaqVPA/CC7mbO+aUALgWA2bNnB+2ScFEnda6vTNnYdCp8YqWTXAxzSoPlfjO5T/mUG9JQGOyxRkwGenZGoVrkWdRqNA49Vjuv1RheOZm2A1EXGbLdoOwU1Vm21kGCJb7me48aVoopo7q51izrJvLV2CcweS+KevY7gDpodEBDAsFC55x/jnM+jXM+HcD7AdzMOf8rALcAODtNdg6Aa8sSUsehZxoJ525KJVJ909z7LNym8knN05aXD0yNs7nIZSnPyD+GDTYxp7L7jaWdbkRFyFoJYLH4LfnZ2m5uUdQtFildUS8XdRFdXli3YbvjEI0XNu/ykqdI+NzMbMdyD5X3jz3bC0ERH6J/A/BZxtgKJJz6T+KI5A/OUw6d5a8LxNpYpL1X+W6jJPTTMn+r35yXfz6yBHaPIWbNQyD2+ZEUvOnw/fDns6Y50+XdFsOnyT7+5rZ3RVkU1d5fYD2gsB96Y81F5EfDjt3+h2jY5YiTj+39UAJXZm5vo7HupdA557dyzt+Z/r2Sc3485/wwzvl7Oee7yxExSxHo6qrOi9dhkc0wjn1FGTQaYGahlFS0URaf+108X0icMzVH6gJZzHbPGMPpM/fLXVdnd6XHcgl4Kgrl0nzXcQbLwm0u/Wxa6LR9Ci4L3V8gTws9Y5WHD4i5NBZRWknrdcVOUZclNETwcol2ug/p5VrZTcIVyl3S9Jdwv7sst0ZvlsesFomtcWfzi93SA2Y6xhlL63gw2xhj3DTnKMrn7FYfNMPnphSd9BslBvsI71ADBjnST//x2WDQWO6weSGRSmyhRu8KhS7DZxR0WfbZDIjlE67lLFjdbFm6Hv6+BQVCz0CrGCTrihs0Ogc3zgjylBPTXG0trvzIHPzyoyfmXSzpq99BkLP/f2ccjh99eHbuugor5aLEKYolfuG3k2YgUy77jRuJj77+UOttH3vDK/Avb31V0dIlMQpY6Ia/bfeY07Q2ZLQJXaHQdUom5ypXc9jFHVDZgF6JkxqlhUL1sVJcHDo1fG7GQNd4jpAWRSO3PrnME1+xD44/ZGK+TMdgVFgG6e+ZB4wjHehgq6u6hUO31nFJlItAbqdomt/7X3uQ9b6Rw3rwkVPsSt9Pjjj5xDjkuRPiuXSHQqccAo08W5BxWyTcTwJFUXlahcU59GIWupwHJZaLGrcjtBlTOkDggT1GlL7131S/NoVh5dDraRq7XL5tqOjGItUAaG6Yay28n9vwd2E5dNfaYEV2hUKXYTzX0lF5sVzkbLk0G7f5noYSlq6SpnQ67t3RyQHkLES3hW47Zd40GJjzKxtyXBlKma5NUc3rgQov83fxWpD9vPNlmQdV3/blC5EnKY6SZgYXTw5fyoUZ/rbdQ8kXxgZYLYoSkNvcwvLdJ8OhO6egRMuNkI4yeISstufyILgRfnDOwco9dnmoB0rY6S1aXA9yR7Hgw2lscFGuOz+l3UQeekI4Whuai6IiH70S8t0wFT4DyZbHuf9zxqzxQn7oeyC6QqGTKBfnm3LQHmRp3LAtimrlDFRsaufSIRdXxkA6+mwsys9AVCVJq09XGs7d79713tUZB/nwDVoyzX0myiUMQ9zMoWfyt8wKtaUHCqQ+n2kRvVUo0m8zs6mCCoAxi7HUwmGkKxS6DJNScXXUyFE7g9No46GHykJQ6K5t+uoPVspFonhsgxR1UTQGx+hbj1RLNnhdIzPzKg51447peb2ph4LSZSx0z3yjvPc0D9+F9bLoD99ge2WhKxQ6xQrQKTVZObkpF5ostnSm4FymRs4t/GiuXO01N+WS54yLW5BOa5Gcj0ehkeCu63gWp88hxP/xzpna600vF/v99t3J9nt90DAiUs3BDXGUWgVvt0XD2oLtGei6wdC3Wlg93aHQCd4qrnjMsaY9QcG5LHSJ6XeiMPKHQRaahe7ycknubc4ubKfKJ+sZYYOUL7z5W8K7AVo7TU7K02NI2Vhk5OhzMzE/ioaK5kwh+avOeaakVnPosZRl4WzaaZZL6AqFroNv/bkXiWg5hnqkqPf7uv3p43EksJ18nldg+nzEVUIoF9Jsp3WUi6eFRl78DpGGrnCpkA+QEDm5ys2XrUkfJE0+f64p25TW9L2VML0Tm0xU46QTdHpXKHTSuZKOJykzOFe+LCVvRznhPsEpj2hdFFWtZ1NeySfFDQ1gmYdSsyRTLsR0VFCqsSf3buJKYeJSQym/OjdTLnK7iXVyz6SxI6xnoKr+5q7TxMpGIaOgJA3cLp98oFsUeoYL15uJzgbdwtr15eNpFnr+GmWh1+79IF/3o0isyorRtkFHGUQdVrBr67+ZcgkUJ3I7kw9htuUfOqi68jH9LlN0Zc2SykAIxUadbTZmLW0c5LpCoVOgq/NSdoqSZHFrUWbSjgD2GkE5d0SiXCzC0zfSJJ92t0W9UqEOGrn8CCldncOVQ6jbYmiD8LHKKcUNKX7oxvs1g6q9PPPvNCeE5LOeaHRqsdFRpLiYC8cJ5dJ+0qUrFLouyJZadTWHVVhmtEW3gmPav03Xzn/Pkfj53x6vpDHLYl+hVyxSEwebfpI3Flks44RDJ1j80TkXQpKSNQ6zjdQBUCkX04Dhs45iuma73kwg0qUWem5RtDVKLXwGZW//+nuKldFKg707FHrgxqLMeYfODLxE0iO3gGUux7Yg08MYelWy15yVvTGq382cCwDX1n99HroBLFa3dm+o8SuJvrEo7AmM1Ruo6FQvF99yfdM7vWPSz1qjvbTeKpfhrWwNf9vvoRkn7awHga5Q6BS0atEziLd2TZc19+usXtON1o1FOQtdj+YJNJSNRfY8Y3q5+FAujVmG5RZq4LT4Xi6u+/QpbBy6/B5sbeC10/NRJ4tClJ2EVZZnoJ2N7PvxeEGufDvkybtDoWs4F4rSpDZ4NW1R+C4I6jho0sCBfCfPp1HLzSfmmnTWch1mPwOtPqMzLoSKj3mOqY6g8nFHPe7gCbj/82+2lyAoF0detnWMc193CG7/lzfi0H3HmG+wX5Z+z3I/ZdEJh8iy2uSJ1IpM+dz/+TfTnRyiSFIMXaHQqTtF1Qr1oVxi9nObBaujSWgcdF5Aisz08z1Ti8u2KGool0zrBKaLmQfd6ycudHKOGdGLSY7DrdVF0bCdjgwH7TPa78AXU17pp7woGnNxUWDUsB6aPN6Ui37ANeXjej+hcpSFrlDoFDitM48Kf9XkscbffCMpJt/19zS2/ufeAtHtz50Ee48alvlu8tdveLnYypOEyngQ5ZQkTf4Y1rLvrNl3gCuCGDHjhyjrMvCPVWRcFAX1/AFhAHC3seQlWXGcMmNf428hr5V8i9sGKx1dodAp7oe66/J9PotKV390DnoNPaQoN6x3/ctb6BR3Q53ng4qD9hmNC95zpFSSxtKX8jBx1lavFl16i0xlw/aOqG6cNrxn1lT/mwJhO1NURqtDLZhiuUSjQAJpjlNm7Isffug4UvqIFLo+3za4o3eHQiekce0U9Xlh40cPx1HTzLvlTGgE5yKUbbNwqWdyUjvPCZqj2HJ5OSz0TFQ9psaa1/AYpIHPnUbGEVPH5fPQUBAUX/rm/f6yNS1AupEQ7OXSODXFns53oDIfVkKcwWQWRePDzeXrr792+kSMHk7bx5HNz2aEUTKIN5gVgVOhM8ZGMsbuY4w9zBh7nDH2pfT6IYyxRYyx5Yyxqxljw8sXV5Ir993RuZwNhNqQCWksiSjWLNPkob0voP2Yp+YUzsUsiyaZE77i205tMv1uS28viyKHwwMnQv9OT6AjrAH5FWZL7uMmrKZtNZdMXRPRJSjdO6dDKZfdAE7jnB8N4BgAb2OMzQHwdQDf4pzPALARwLllCUnZSqtTVNlFoGIK3yehzeVQ97fOYyeEqzenY9q/dXmRoi2yrImeG5AY0RM9ohKiwve811hyhMrOiRy6bh0jl1fm92IQ9VZXA6JHAt3a9WxDhkXlomVQ143KhlOh8wTb0q/D0n8cwGkAfp1evxzAWaVICL2SyU/zIynsEuDerJEfAHRKUnenvyz265TBMy9b/nfaoqg7jQ9ilhnFAyeCBdhYFPUMPlf2SULN9qJQS7HyD8wolKop+r51fbYdIHHojLEexthDANYBWADgaQCbOOeDaZI1AFqyUmTbuqv+Ji9sxnKCKUq56NKEeokE0QeGe5obi8xlyTOKrEuoe0ByyqUBh9uq9N8pqshqTFm8e8YYFIZSyqVxQpBUI9lFf798iyq2xk5RKO+gzZRLcD62llCwjFYe00dS6JzzIc75MQCmATgewOG6ZLp7GWPnMcYWM8YW9/X1BQkZuvVfPiC5yOngvsjd6WjvuvQUxSNfmzklv2BIkk1cb3RQwtZ/R6aJ8qdQRpG0kEmuGGXGyiOw3GYsF9cMVG0vful9IVMuZcx+Q7MMPk2LWOCBE0cZb29nFEkBLy8XzvkmALcCmANgPGNMLCdPA/CC4Z5LOeezOeezJ02aFCRkxhoUVqKSRleVI4f1NCIXuqqazkfrrtHpHn1wL6rVaM73r0+eHlx+UmZqAQYYE/l3QXuCGFY8ZfYhwzf6ZLtRV3aKmryjolmqzL6O0jw2sfk9S7m423dM+ObscyygroxTZuh1WCcoc4Dm5TKJMTY+/XsUgDcDWArgFgBnp8nOAXBtWUJS4KNUSfl5XBe8c8PDrOC7rTFGOpgixE3KaaFbXf7EJ7NvLKJSLoREMRfyAN3iopnCMyFkwA2VfYh4pmir1UkzHnpJdILT0o77xFHWTDpAp1McNqcAuJwx1oNkAPgl5/x6xtgTAK5ijH0VwBIAPylNSk2b0S7EWbJwNwCiVUl6a34Wt1YhEoqRFSz12Yusrcpnipryt5ZDuM8mln52RM+PIldZsUlCF+uoB4nHclv0nV3V62FWLzV/ZzqC4WPKl7po7RsKup163anQOeePADhWc30lEj69dFCstNAOQ/3dfi/9Zjllc89Inv+kuKGFiGzeUOKViTNEQKs2WfiWQz+Sz5wvfapevA5sB1zYFkXdis2cgOTp1CEKzJ9yMf1d7Cn0dlLra6YrdopS4GvthZdDSJPrXPbWr7PQSZRLxjKiKSCXZWaNh86yaQVUlzrW+M+BSDRY43dCflR/7mKzvbw8oZ17SFkUNb0d/2Pg/K6bylMXRdvNOvi6CNPy7B7475FtA3Q6hnLMW+bnCB4SST6ENMS8RV66Q6VDBo4iaHq5uNMwppzzqs4wiHK5XO1UWSinRbmgo+oyZbYw/gZF9NDwuaGgZiPHcmmnyis0szb8bSvD1D46gT8HusRCz3q5GCgDRx7xGjyFT/O1lvJaxpeDDuUcm9eTT9sBF6aydPRQDA7dVa56jZJbjI1FRTyiZFDGDvXEItPr8Y1cGdodGvGKmMFC1w26gWXZEDzDKDgzCUErjYSusNC10NAUek8Q8emaaodZmRS4pt7qleR8VMrAof/bXr4eLoVhvztfXqtnGEl+7gzJHLrlCciDZwRVJjZ6OYPPFZypNK8T33H6mZtF+YkRjFCPMtPgY3/fJfCHJaE7LHRN3JD3Hjctk8ZU6eJW56npRFlsDUjMJLxfLYFysQ1WPjDGQ08/yZSLLR1oiuE9s6Y50+jKd12z59EZHc8XvrGI3IZqsXqQD0RpZ42aQjnvZziYoixZYy6wFkF3KHTNtffOPhCrLjwDw3uTR2hnP6VQEKbf1N+TNLEtdLcl0uDQPQ430OXfTGPHF941E2cdc4C7MEu5akm0RVFKnvEpl6Lt03V/D6EA3YJ3ru1S5ZHyjPmcVLjKue/zb8bV582x50HMz/d9N9Y92qCTukKhU2CeOhe7nwKXDnS6kGnoI4riiWsJuG30Bn3F7IeOdLIRTKVcYiCQztWnZfb3E+qHzpXvVGRiubTTRtf0nebf9j4TlVbNi9IWdIVCz3hUGKaWzFCloQ3WBNJin220J1xLnoRgoWf+pnG+OupJXsSk7RTVX8/kSKonQiLHJgTd+oBt92KcCI9xGpPPWplL7qLPxXJ/uMpLEubdFjtBrSXQt/VioKw9VJSLAyFKJp/OwUFSF/xIitaiXKWfTAs7ugVerfgh7cbhJURTMo66ZO56ihEzRr0WMoiEeEuEUC5FPR18d4r69gfTebEqGsaVD0XXwlmQk940WPAxZOyE9ZmuUOg2tGpDUSO/MkZ9XbOkKCcif0zxcmkuctkoF6YtSzvDIMlfDmxtgnxIdAssTT/KJfk0uy1S8pDai0EWqkw1qb1k7mm/TmvAbwc3zQgz398Z6EKFrq86U0c1nfOZy7XAGxFlcENhroVQ9fcao1lcITK71hpIXi6562ZrLyb0dZAlnlTkNycp3wO6YqOunNapfcD1o1x0z+ben5FJb1kU9fVjz0RbbKNlausn+hmdYYZalLJinaHUu0Khx6BcYiFktM5uAEotm8wBEWoZjGZxEWWibEDy8XIxyeEqw1RmEfhmEWdjEdHK9xDOaXA4fg894EJdY2KMkSIoylv/dfnGxjUfO9Eghxkxdhab8rGm9y8iGrpCocswL4o66BDfBm8ayXXeJup3SwPQunipVoamHF25tM1A+nvMedG8XDLXNd99PXtioKmkzINlTRdnwZJXp8A5y1Pbi3PNyJ0/RR4ecG8Ijjt4ovc9PoNc0UfomjNFOwEki8FJuTgaeATLzecem9VMUYhqvkWDcwnYZ0PJzTrOXy2vjBAJ+jxMUhjS5wbK8iDnHWtRVM5Gd/ALFWrd12hg33kAACAASURBVBrvlvruUhm4/zuIibwxYR7k7PlY+k+AHO1Cdyh0kpLpHPjKouN1KV4uVMpFzdsmA0XvOC30FlorvuWQKRcPpWb+3c9KLpKXLwduk8VmQKmxXNS0neDpIeCz7hRD6k5w2ewOhS79bbINXS/Kh2cuCpsXiMtaF2lI8li8FkzlmykGN+VCKUskaEfTDovlYh/gbPCxuosviooyTRuLPDLTpG8YRsR8hDx13l5FZptxUfoaqYz262kyukKh0+CiXBx3kymXkAaRaWbOsnWUi65UV+PVyuK4TovlYqctWtn+Y8UBLwMxi4ott47m87u/aQC4KZfyKt26eO1FucSXpR3jQFcodNtOUTDD9TbAdAKRC/lHyp8pqr2PWozc4YwWOjEvwr3UWDRR4FAm6gCVs9DV9Fx/3QSKB5bxd2IZgJ5Dz5ZVzERvzGTBvDYL0YJztS5+rPxYukB05vZv49Dzv/Uo3F0n6B+gWxS69LdpiuWyPGN5XVCS+Y7UWgvdkSbJl0npaQ3S5b1jVVBSeab8k+/tsU6a5ZtLj+K26ClPDLjKDN36L24T3j++/SDZ+i+1rxZXjp1q1MyGQ8rQ3PTpN81wltUOdIVCl2FquDolCNC9XKgo2mDl+5tb//Ojva+XS0j5uuu2Ay6o1j1Vfgpsg3lyLS+HNZaLw7IicefuJKls8Tp5zAVWOX0zvLQvdWWfMehLiwfKQTfaWC6RGqbPgmsr0RUKPRPZz7CoVTi+s3L/188+qkBeynfP2UHiJWLglgxXrBalPA11dARXnHMhnyl/+fcPn3gwDttvL0uOdqgLgDpF7eum5urQzSKL986YHdw3losv5LutSlocuNFYRFfaYastVWK7L6kIbVntVOxdodCtIV1Z9jP3e2CJh+w7Bv/vjMNp+TUsQ30iWwyN5FonWOjqQ9DzMIX2/fKZR2D2wRM8JVTzjgvivqKOsLZkyJy19ndDet/8g9KXUFeh9e+KehhLVO2s2uSY0cIz6LpEoTdhPHHH8Kao4XMLebko76uo50VNQx/pp3j2gYLym5q/ja5oDJ6OAmJaKyGUiw2uHZWdpsgFnCduFRacNfKh5CTL4/sOKKDqQN+jAoPk69A2oYNToTPGDmSM3cIYW8oYe5wx9qn0+kTG2ALG2PL0s5gpZkH2MAWlA0rXiyxkUd9ZLMUp/61TU0U3ipjLd1AulE1cDivXJnlsa8XXk8SlGCleLuS2UgLlUsQPXXenuE2uF9sbUjcWqX+3GrY1ENcBF1HLdeifVoFioQ8C+CfO+eEA5gD4BGNsJoC5ABZyzmcAWJh+LwVyA8tNmc1a0Xa5FFBmA3p6Iv+d4ilDtYQplrzPIpczuiIzfmk5cm6L6qKo4b6WK6kAjtx2ahQV6qIog//msna8YeN7k/6Oc5gJbcbdKYukToXOOX+Rc/5g+vdWAEsBTAVwJoDL02SXAzirLCFlmDa1mOquyWu7Ogy1fN1F61d3nppZB8kPXUyTPUp0rTWQgnMZrqty6cuP28qpMeGb6VuHmM8aizIUsIUVJi2Mm2Z6kRQbuT9afyvvbVMMs3bAi0NnjE0HcCyARQAmc85fBBKlD2A/wz3nMcYWM8YW9/X1BQmZ9XIxyhaUtz/iT9nyFjptqpi9ROMSTZY+aU3UWPfm72W/Ft/8cwNliTO7mI/uyiu/Ycqv9KaihrUR6GahZbzjYGYuQwXFkYUySNHmfeWDrNAZY3sBuAbApznnW6j3cc4v5ZzP5pzPnjRpUoiMmYU60yG/rikYxd6lgOZ9YrNQk0/7Bh5/a7PQoi6aCoDUkVzWIk2U6KAoMdOiei6vDrC29hrR2/jb7bbozs82sGeOoCPIll0UZdq/WwFqlNHo5WrKMnu5lCeHClLzZowNQ6LMr+Cc/ya9vJYxNiX9fQqAdeWIqMhiuGKrUPG7Nd8IUzxBV/guzuo8LVxb1F3lmNIZbyEpBD3FY1MOZXdv3/yLWrI+KKpQZKXZCM5lKovwHLZF0YZhxKgcenxLOJN/QNsGsu0hlmNBB4ztZFC8XBiAnwBYyjn/pvTTdQDOSf8+B8C18cVLQKJcyio8ABRZbEq2RnQdo5xEZCtX4LXTJxBnHvo8rMeARXwxrvWLmFxtJ2zlluOFOClFT3G58ikrP6uXS9oZbTu2OwWlWui5vDvDy6XXnQQnA/gQgEcZYw+l1/4dwIUAfskYOxfAagDvLUfErEI3Ui6uKamjxRdRiKEQz0XZolxksSn77NmbvveBWTjt1ZOxbfcgLbNcDu7vseAyHCn1oQZVKpKXM4+CNSG39cbfJW8sogbnynLV8d843Q9d+W6ggopAv0HJxBW0F06Fzjm/E2ZZ3xRXHH8w5dOYLlJtFz34gLLgmVAuHjIxeqdSk00YPSy5TrnXkIetjKKdKrbfeln+/ToUplwyFrojbcHCsha6bXOZoN0MvxeSgg7K48ZyW6SUn1BVtLRloit2itoWRRtgxSou7sG/fnnpKIy8e2Y+T/KZohZ6R5XFmo1kxZnyt5dSTEG7WQf3Q7RSoYdCLIaadmPq4G2RG77bFBMgUy7+dJ8Pwjl0j0qjykJMY5K54xZF2w0bh06dKoa821a9CFPnuumzr7ffF9B2Y3DG1k6kylWwU2W8LyK8jzKtNhUhRb12+gRcnAaG65EpF4fgvmXlDn6RMvCJh57/u3OoiHYN3u20GbpCocswKZPCiiOwfFkGWxqfTMXXw/bbC6OH91jKpZVHoUBobm/6RNSAVyHwsehp0/B4MzGnZAEVMWPyWIwcJt55nnKhHnBBLVq3U5SCzOyhrZMeizERqwTS7JVFLDEcXafQ1Q750pZdAIB6nVsr3kfpWdNpXlruoFzPctRLdKVD41jln3LpPNpgczakv66Vi5499hkznFR+EdAXRdvfOWVZY/ihW++XMrINotpYLgHeVrHgu14VVIZuAZQwgLbDU6orFDolZsVAvYVEVWTkZx1h91FgjIfukxeB3/fFlR+Zg/mfOqVQHkKKIsfC2RYE/eUppugybovppzE4F6EE3a3iLqrbYuM+E3XX5nFQLr7URVHdtfbbACS3xY4C5aQSn/ua95NNdOe9vucTuhcZ3aKQLZUCHDpTPk1lZ+WiPcyJr9iHlC4vE22WIkDZsGW77gO5KI91a+39pVvo0uyLxKF71ntZsLXFaBayrs/nZqVo5dGpRnSHhZ7xcjGl8bseCkoT8W1GrkVGXRr5IjPcU0QGHUzRJF3fWw1b+S7KxWdNpuzHlBdFnRw6IT/qe6HMUkw2QquoKtL7aaF2Y5Ye2Eo93x0K3XIEXTMNL9TBirTDHIdO5PJNsbfpfL5/OpOl5+MR4N6kFadT+3cEd7mu5/ShXEKs7szvjgQ+O0VNG+5cyC2KMpqFXstqcWvaMnS88TzeEvj8oi6hrURXKHQZvhZ6gyaINfsqoXVSOHRntMWCg0DMp7K5LbZ7sVFtP0Vio7hAP79Uf02+X8htfH+RKBeAGpwrvuIMga3sVm8i8zs4uxx0hUJ3HUOWS6S5HG0bsPZa+iJ59js9TzcHbyvXmb9Fwbqu01Bex/HpHA1awnKT25+bTrn4ICS7Hql3ivvNbovu/CiLoi4OXZtHCa8/eHZDaOsxZKFca4ft0h0KneDm4poqh7gt6k+Zt+fjSqPtLEr6HmJLCPG5zTW6hgLzoFw8OPTYrltU90wTYlAusSkxEzKxXJwuG8VKyz6Tuw7Ktn6/8/5j8aE5BzvT2dtDeTL29NTwsTe8wlBu+9AVCl2GkXIhBi0qE95cW5o+x39q3grFf91YTsYjoQTKyFKeCq+NQpq8i+4Wpbqy2ZKRA0cVrOoMrSGm84ay85u7PGeJDQ69PGPCBwdOHI2vnHWE930e1D49T8OM+d/e9upMWe12BgC6RKGTKBfoKz63xdnjXt+IaqHKRs0zyEKn3kPMTwdq8CEKxUOFH+XiLozqthhbS2knZrpBWrqW3Vhkzz+0nnX9g+a2KJfdAZpMg1aK1Sl10B0KXQ6fa5C4jLgresrF/eJCLXUBnWsdyX+dkL/dR96Wh+0+dQ2AJFYQYlq9ZcP1flyLotkDLlIL3TDEedMLSuGyH3qd0Jlqmk1PJpRBfZhDWMiDYJxyQ2bC7UJXKHQZxmOejA09/XRZOAVkKgpVNnJDDOHQbZxjoJVv86MvUq++OxYpZdGrNkKLKEhL+HiShOoucZuv8iuD2oiNMsXKHbJtS1tFW1TRrBFT4zFVGtXLpchCl02h0fLMptda6BYOvch0LwYX6lLwMVGEMqKmcaWLx83af89y6I68Agtv+qE3L1tPLII4ZlEebDpHo5dBBcVoM61CVyh00saiFsly9uxp+bK5kMF9pqgOeQudeh+9wwv8+ays/LbToLLp3IceNL5nfqPJ1Sp0wwEXoh319sjv12WQFKO9ZKVct8RFKsutswhed9i+2ustDZUccd2oCLpCoctQ6+ob7z0aQHEvF4qVserCM3DSK/SNJ5NXwRdK90On5tf82xozpYSG2MrFopgHXMSQumgeOj5YbufyIOtdzeJAc+V+Vx3WdOnaqMA4B/73707AzCnjElEyijUWh+7Op1OOpOsKhS7ranXUbS6S6jW6Kf5ITARbYobtyyHlhk57qZSLdVG0RS03pBx1oA+Z/ZQF08EtDStYSuv0conoh26b7YrFUJ/TlFqJTpo5tAPdodAttAD5GDYHYjYEfw49LJHsmRBjg4WrDoT3g8+aQcz+FeOQ6G5yZdNN44vsFLWWBakx2XbaNvzV45VdBK0om9quOmEw6QqFLiPv85xd3Mmlp+YbLFFxhHsoxFHU1PwG64LbzTabTmjIVMTw1y+C98yaiks+eFwqC/0+F1XkzReri6LpK3Xoc2nxtDMWRdVB3ibLsq+8rTQ5EqOq/R2hKxR6litULfR8msy9hvs6CcEcXeRpr6uKBobqAIBeH+0Rc+bDir/HmEfQOfPQXBtWq2FEr63b8aBDNrzrRfVDl6S1LYD7HlVXFvKGnfsecbTfngynQmeM/ZQxto4x9ph0bSJjbAFjbHn6OaFMIbM7RbOoOSx0032UBFQa59X7j6UlVItk2U//++Nam650QqEPVxSSeh9pikqUyQda184cPWTHkVP3LlWe0BORnAdcEPLQbpQTn1JbtEko2pyPS2UrUYYspDxZ+wc5gGah/wyAOleZC2Ah53wGgIXp95bA5OLHuaFDlyzPHz9zKs48ZmqhPEjtxfVswbSNXIY9k4HBpKsP86BczBvBwpAbPDzvtynGBZ85FR9MA0J1km81QJhZFOXQpfyt0SrTZPKO7U5S6O1Cp7QXp0LnnN8OYINy+UwAl6d/Xw7grMhyKTI0/84H8ida6M7+EPZCXjk5zDrPlB1soVPTxbHkB+p6yqVMOktVLhTr0QZbkhmTxzbziEK5aGiyot5IxuBcYbRd87dmGtssokm50Dn0MhV+K/afxF6rKhOhHPpkzvmLAJB+7hdPpDzkBpa30FOF7uBHWlHX4Z42lM6ouxZ3YcrVIAeHkgd0US4+eZaNnNtiC3ebhD27/ia326JaNq1wdacoQIsn3yleLp0ig82FOObB4y6UvijKGDuPMbaYMba4r68vRn7Z7573v3b6BPz7O16dux7VbdGSV2iMdVs57kMbqPk5KJfGoqjSbHK3xalMHY1mHTwIedLrwp3GJxRw5j7H5nq9PGImavo9SJRc/i4OvacDlLcNrThRTJsmeqlhCFXoaxljUwAg/VxnSsg5v5RzPptzPnvSpElhpVnOt/Dl8qZNGI3D011lMtr5Qmgcum76nqDGInllOPIQCn2Y0qvL5A99dCalDlq5U1QH0gEamtLdi6JuiW2LohmbgOCH3soQChS0isPurKfOI1ShXwfgnPTvcwBcG0ccPTJeLgYOvZtRdKeos7NTuXbH7wMp5TJMpVxsM5ICs03dBqair7td8T0y1x01bVO82a3/zb+LP5dsobs5dN3hG62EqcQyJKEaXLm23oZ6obgtXgngHgCvYoytYYydC+BCAKczxpYDOD393hLktv43OPRi+doOx6DC9zAGITOlM+qTJFd7anHsEyrlMqxG59BjwhWfnlQLkReSfWGM4110YbEgFUnl0HWutt1vUsVFO23MXlcCzvlfGn56U2RZLDI0/85v/U/TENVpJza+ohx6T41ZFRD1N5cYYlG0V6VcVEtaVg4lLwj58tgxKRdnBMSIrU0ndwylqsY6Srxc3HLI5ZXsUemFcvzQ3Zl2il7pjp2iUhPLucwhkoVOvNYu2PzQi/CZPt4K/Q0OPdtsqEfmhUDm66lWLEUhxYDTs8pQVNhOUF357rLI+RN3iopY/bbNfu1Eu/zBy6IdfdEVCl3giKnj8O/vODxzrdGJOf1l6ipYv6svDBQpMuEMCP7B2t+YRLkEtOOPnHIIjjtI3uRLpFwkJXv89Ik47/X6088BFHYU/vXHTsKBE0elWYUdCSiDHm3RK1syiu4UNd0dOlA1jAIRy4W5fP2TT8oxda1AqKeRD4jrxR2BrlDo4p198V2vwYQxwzO/NRaLOr6q8/Dpg7Z4y6EW8ufPmJk9G9KRjaBcZAv9S2e+BnuNyDJ3JLqCkIbzZOPWp9/0SkJqWt7kTSIx/PoLy9CE79Z/17vMrd9lLHTzfbo1q05yTGiXKF2zU7QTYI1pLlnoRdDqF5LpBAU5dMbiTHtd1qvYKSordPdiXUGhMlmxwu+JRWjxRXfo2owP40EtjXZelh86LV1zzaqJXQNDxQqPgDIVOSnvztDn3aHQm9BZqfap6B4DHYcuLYrGKcKez8VnH4U5h07EwfuMtt7DjV8C5YrYWWJEWwwxHkKf4fK/PR5vP2J/Z4RLXytZTU2tF93O7LufXu9Vtg6nvtJzj4pBXt9q/svjD8JXzzrC8y49OkEHdYVCt59nKRIVLCSC0gjl84oWHcqh5+Rw5HHcwRNx1Xkn+lnobYCVA26ZFAr3GtA0OIDXv3ISfvDB40qnNZqzPdrA4cOhU2S/5IOzyPnFxAXvObIRkM2E0BhB7ega3aHQ009bpXUjh+4D3bOnDEg0z42QXHSGY+xFpDIt9LZP1SPl5dsG8hx69tMEMZaLthcLsSjP9nHonYGuUOgCWgrdsbHoc6lXzOjhdpf7shqCoCc+edph2Hev4do0odaXsJISyoWex9Txo/DRUw+NJEeY7F6bsGxeQJ7Fq+ldnHUZCHke51KFksBXfGrQMkHvdYqXi4BL+pNsh6NHRjuVe3codMLONQ59Z/jACQdh1YVnOHnmmC9BluO2f3kjAOCzb3kVFv+/08Pz1FwTncrXOrtr7mmNga4onN4UsTs+USGbQN9YVHwQKWVR1FQWTSQjRvT2ZMoxgRqu2hehA6g5SEIW//eROWEFEGEzhlo59HWFQhcdwBagyqU4OsueSCBEDu2M6RGf6KnFsShteZh+aoU10krqoiyElMss30Lyz7gaKr+NHEZTBdRw1Xsq2tkPKOgKhS6gp1zsv1Oha55ym7WfA1kMlM6oG8yG6jLlUhw7+/3dz3Ry7Tt2RAxxiOX7pddZ6AdNHI1RynmTsRV/8VhD9nyKrqOI8zZFmzKh5pCj3dBVw/GHTGxJ2UaX0paUnsAZy6UTQGk8Zbav8049FHPflo+hHoKY27VlyiVGo9kZ4E+slvv0+e8gDTDtsmh0ot36z2/IXWvJzEP+20V1aK7Zd3TaM1TvFQbL7oHmaueyr7wNr/6PGzLpxMDh0Pve8F4LIWaw8vx3BMnjC6uba0skSNAVCv3JtVsBmCrNvihKhc36Z6AtGoWHCghcWGxEa7QH56IiSKErxarKPIYlF/MwbF1erTzFSEbMjl70/QsOffdQU6HrBuYmlx9XTUXzclG+t/aEqvYTL11BuVx0w5MA9C+dEpCJAhflUgYaHLr0WMN76K8kNuXi87xiu79PRzwhwtTXVlqrFOS0CYnn0gmHxvOcmLHfXgCAYw8a7y0PQB3MzL81LfTmoK6jcXpKstBjoV1KlYF1xLpCVyh0EQzK6rxPrUzWGQukDMDuwaTzyAtSi//jzeQ8hiJTLgLvPvoAZ5qxI/0mdx855RBc/rfHB8nTji5qK/OVk8fijn99o9b105ifre0CmD19Iu741zfifbMP9JYHoAUdsy+KJhZ6v2Sh6/Isa1FUrp+bPvv6qHm3Ah1gnAPoEoU+UlmwklGWG1Ux0OgZsQg5alhTOY4bOYxcSr1hoftJ58LUCaOcaYSFvmNg0JpOvJf9xo60vsd2IbQjHjhxdHRr0Jana0DwlUXtL8Ko6B9sKnQ9PZV81ks00SeO0e/X0EEdVzpEr7YN3a/Q089SZzsltRLBWY8aHqbohqSNRTF1i+rxAQCvnjI28/24g5Owu2psdBWN9xMw5B6W0hDN0LG88AYaKtppce0zJvESEs8P6OmPTExy5ecpe4/0KlNw6LJC1+HQfROZ9nfkL3uFkby43EmsmDY+MULU4xFbCdFvXjN1bwDJAA0A0wgGUix0xaKoqCitlSIt0pCslA4y5YWRMzpQoQsrSe3sN332VAzrqQUvNKkKffbBE3D+nx2ZufalM1+DM4+ZildM2gs2jE4t+e279QuuCz5zKsZqZiVXnTenMWiMHtFjzCPkoIX5nzwFNy1di28ueMqYhjGGaz9xMs783l3EXO3QGRympnj0geNx5UfmYPb0Zqx6nzZy6YeOw+kzJ9PkSj8bHLpDoZ/7ukMw84BxVlfAP3zqFIwbNQwnX3gzAGCMZpf2rf/8BvTUGE656BYAQK8cH4gkeRZfeNdMvPWI/XFigXWNmz77eowe3oOTUrl9wBiwz14jcM3HTsKr90+Mnz+fNRVT9h7Z0l2qXaXQdYh1YlEUBMoQSkXIi6Ky8j5sv7GmW2jyKMpj9vSJORlH9PbgREJDHZsq9K279NTMjMl6WedIHVPQUNt22+kdKmYeMA4Pr9nkTHf0geMxalhPkPdPUah1O26UnYqTleBbXrO/d3nDe/OUiw61GsPJh+2LwSFzusOnjGusDwH69Zbp+44x3k/bl5H9vt+4kaS1HxvkGVEohBECJEbByYftWzhPH3QH5TLcPB10bbgQEJzvPoZ4Krop7d5pJxo/is7pUSGXZhuwbBiS3BZjYpJSRxNG03l9FaIzb901EJzHOEseoXVXBoQs4zX1VTQan3MROrANiLuEhb5rcKjQ+27m25THdwGdMrOcPDahfKg7XMtGpxxw0SUWevLSqJbSF941Mxdf+a2vmYzz/+xIvGfWVCx6ZkPj+nX/cDJWb9iRO3UHAP7itQeCgxs9D2JhuIX3s41T9Xo5HPrpM/fHRX9+FN525P74/cMvFHr+vdLOXMS6FtapauUzBhyR8pW+8K2ueZ98Hc747p3G37/z/mMwKz3O7/p/fB0ef2FLkFwmCI7bhH08FhJ1ELt7OQd+T5DfZxF2L8tC/y/OPR6Tx2X5eHkA+N9zT8B+4xLZrjpvTuO3r599FN48czJec0DY+6fg2k+cTE7bKV4uXaLQk8ZsOxmFo9lJT33lpBy3yxjDB044KHPtlBn74qhp43HUNL3vb0+N4a9OsMdK1qFVLzcbbTEeemoM73ttosRDnl+G4MdNlAsFgnLZ0T+UscjPOHJKIdl8cIiFIgCA02dObkT0nDZhdMNXvVU4YHyxhTdZqcaWf5zFQj9lRv5gC3kz0OtmNCkLmYbbe9QwnH3ctEgS6nH0gXq90MkoNF9hjL2NMfYkY2wFY2xuLKFUCOpjYChvr/amPuojemsNntel3sTpL7Hjs1DKF2X29tCUsC3sb4/mOfZ2cK2txvhUHjH4iNlIj4evpTytlp9PtRLFdzm9aTFR3NpTs8vhWoykLlb2puUM66k1yrTNzCiQPYwOGO/2apEHQ9F2hPyuE5FUuFLLr8bHDbHTId6ZWl8dYqAnGwRC/gHoAfA0gEMBDAfwMICZtnuOO+44HoKN23fz8+c/wQcGh3K/1et1/q0FT/LnN+7gz23Yzr+94Cler9et+Q0N1flFNyzlfVt3Bclz5aJn+QXzl/IFj7+Uuf7chu38Ozcl5d/w2Iu53znnfOuuAX7+vCf4roFBfv8z6/nV963mnHP+s7ue4Y+u2ZRL//S6rfz7t6zIXFv64mb+4ztW8p39g/xr857g23cP8J39g3zuNQ9ry7TlL+Pah57nNy9ba01jwk/vXMkfeS6ff71e59/4Y/J+OG8+/+6B/LvknPO7lvfx3y1Zk7t+2Z0r+aNrNvGbl63l//zLh/iXrntc+/7+5+blfNXL2xrfV6zbyn9w6wr+0OqN/Bf3rGpcf27Ddv6Zq5fwx5/fbH0ucT/nSf3c/tQ64+829A8O8fPnP8E37+zng0N1/vU/LOUve7a/3y1Zw297sln+zv5Bfu7P7uPfuHEZr9frWvlkrF6/nZ8//4lM/S9fu5Vfksr/2wez+XPO+a8WP8fvffplbX7/deMy/q7/vsP4+3dueor/zWX38e27B8jPd+uTZvlbhSvufZY/+OyGzLXL7lzJH3t+E9+ys5+fP79Zf5ff/Qy/+IZlpcsEYDEn6GXGQ49NY+xEAF/knL81/f65dIC4wHTP7Nmz+eLFi4PKq1ChQoU/VTDGHuCcz3alKzLnmwrgOen7mvRahQoVKlRoA4oodFJET8bYeYyxxYyxxX19fQWKq1ChQoUKNhRR6GsAyP5s0wC8oCbinF/KOZ/NOZ89aVJ+RbtChQoVKsRBEYV+P4AZjLFDGGPDAbwfwHVxxKpQoUKFCr4I9kPnnA8yxv4BwI1IPF5+yjl/PJpkFSpUqFDBC4U2FnHO5wOYH0mWChUqVKhQAJ0RCKFChQoVKhRGpdArVKhQYQ9B8MaioMIY6wPwbODt+wJ4OaI4sVDJ5YdKLn90qmyVXH4oItfBnHOnm2BLFXoRMMYWU3ZKtRqVXH6o5PJHp8pWyeWHVshVUS4VKlSosIegfYYk6gAABP5JREFUUugVKlSosIegmxT6pe0WwIBKLj9UcvmjU2Wr5PJD6XJ1DYdeoUKFChXs6CYLvUKFChUqWNAVCr1VJyMZyv4pY2wdY+wx6dpExtgCxtjy9HNCep0xxr6byvkIY2xWiXIdyBi7hTG2lDH2OGPsU50gG2NsJGPsPsbYw6lcX0qvH8IYW5TKdXUa/weMsRHp9xXp79PLkCstq4cxtoQxdn2nyJSWt4ox9ihj7CHG2OL0Wie0sfGMsV8zxpal7ezEdsvFGHtVWk/i3xbG2KfbLVda1mfSNv8YY+zKtC+0to1RTsFo5z8EnIwUufxTAcwC8Jh07SIAc9O/5wL4evr3OwD8AUlo4TkAFpUo1xQAs9K/xwJ4CsDMdsuW5r9X+vcwAIvS8n4J4P3p9UsAfCz9++MALkn/fj+Aq0uss88C+D8A16ff2y5TWsYqAPsq1zqhjV0O4O/Sv4cDGN8Jckny9QB4CcDB7ZYLyVkQzwAYJbWtv251Gyu1wiNV1IkAbpS+fw7A51osw3RkFfqTAKakf08B8GT69w8B/KUuXQtkvBbA6Z0kG4DRAB4EcAKSDRW96jtFEtztxPTv3jQdK0GWaQAWAjgNwPVpB2+rTJJsq5BX6G19jwDGpQqKdZJciixvAXBXJ8iF5oE/E9M2cz2At7a6jXUD5dKJJyNN5py/CADp537p9bbImk7XjkViDbddtpTaeAjAOgALkMywNnHOBzVlN+RKf98MYB/Ex7cB/CuAevp9nw6QSYAD+CNj7AHG2HnptXa/x0MB9AG4LKWpfswYG9MBcsl4P4Ar07/bKhfn/HkA/wVgNYAXkbSZB9DiNtYNCp10MlKHoOWyMsb2AnANgE9zzrfYkmqulSIb53yIc34MEqv4eACHW8ouXS7G2DsBrOOcPyBfbqdMCk7mnM8C8HYAn2CMnWpJ2yrZepFQjT/gnB8LYDsSKqPdciWFJVz0uwH8ypVUcy26XClnfyaAQwAcAGAMkvdpKrsUubpBoZNORmox1jLGpgBA+rkuvd5SWRljw5Ao8ys457/pJNkAgHO+CcCtSLjL8YwxEa5ZLrshV/r73gA2RBblZADvZoytAnAVEtrl222WqQHO+Qvp5zoAv0UyCLb7Pa4BsIZzvij9/mskCr7dcgm8HcCDnPO16fd2y/VmAM9wzvs45wMAfgPgJLS4jXWDQu/Ek5GuA3BO+vc5SPhrcf3D6cr6HACbxTQwNhhjDMBPACzlnH+zU2RjjE1ijI1P/x6FpKEvBXALgLMNcgl5zwZwM0+JxVjgnH+Ocz6Ncz4dSfu5mXP+V+2USYAxNoYxNlb8jYQXfgxtfo+c85cAPMcYe1V66U0Anmi3XBL+Ek26RZTfTrlWA5jDGBud9k1RX61tY2UuWkRccHgHEi+OpwF8vsVlX4mEExtAMqqei4TrWghgefo5MU3LAHwvlfNRALNLlOt1SKZojwB4KP33jnbLBuAoAEtSuR4D8J/p9UMB3AdgBZJp8oj0+sj0+4r090NLfp9vQNPLpe0ypTI8nP57XLTvdr/HtKxjACxO3+XvAEzoELlGA1gPYG/pWifI9SUAy9J2/wsAI1rdxqqdohUqVKiwh6AbKJcKFSpUqEBApdArVKhQYQ9BpdArVKhQYQ9BpdArVKhQYQ9BpdArVKhQYQ9BpdArVKhQYQ9BpdArVKhQYQ9BpdArVKhQYQ/B/wdz7Q7/pN5naAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(clustersize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "save data to file:test.txt\n"
     ]
    }
   ],
   "source": [
    "#WRITE DATA TO TEXT FILE\n",
    "f = open(\"data/\"+filename+'.txt','a')\n",
    "print(\"======================= SETTINGS =======================\",file= f)\n",
    "print(\"L = \", L, file = f)\n",
    "print(\"K = \", K, file = f)\n",
    "print(\"h = \", h, file = f)\n",
    "print(\"b = \", b, file = f)\n",
    "print(\"layer = \", layer, file = f)\n",
    "print(\"nwarm = \", nwarm, file = f)\n",
    "print(\"nmeas = \", nmeas, file = f)\n",
    "print(\"interval = \", interval, file = f)\n",
    "print(\"ndata = \", ndata, '\\n', file = f)\n",
    "print(\"======================== RESULTS ========================\", file = f)\n",
    "print(\"Mean Wolff cluster size = \",np.mean(clustersize)/Ns, '*', Ns, \"\\n\",file = f)\n",
    "print(\"y_t array = \", yt_arr,\"\\n\",file = f)\n",
    "print(\"y_h array = \", yh_arr,\"\\n\",file = f)\n",
    "print(\"======================= SANITY CHECK ====================\", file = f)\n",
    "print('MatA_even (lhs) = ',MatA_even, \"\\n\",file = f)\n",
    "print('MatC_even (rhs) = ',MatC_even, \"\\n\",file = f)\n",
    "print('MatA_odd (lhs) = ',MatA_odd, \"\\n\",file = f)\n",
    "print('MatC_odd (rhs) = ',MatC_odd, \"\\n\",file = f)\n",
    "print(\"=========================== END =========================\", file = f)\n",
    "print('\\n\\n\\n\\n\\n',file = f)\n",
    "f.close()\n",
    "\n",
    "print(\"wrote run data to: data/\"+filename+'.txt')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}