meeting_04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 1 solution\n",
    "\n",
    "Let's use the following convention for numbering legs:\n",
    "```\n",
    " 1--A--3\n",
    "    |\n",
    "    2\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using TensorOperations\n",
    "using LinearMaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rand_UMPS"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    rand_UMPS(d, D; keep_it_real=true)\n",
    "\n",
    "Return a random three-valent tensor A, that defines a uniform MPS (UMPS).\n",
    "The bond dimension of the physical leg should be d, and the bond dimension\n",
    "of the two \"virtual\" legs (the horizontal ones) should be D.\n",
    "keep_it_real is keyword argument, for whether the matrix should be real or\n",
    "complex.\n",
    "\n",
    "This means you can call\n",
    "`rand_UMPS(2, 9)`\n",
    "or\n",
    "`rand_UMPS(2, 9; keep_it_real=true)`\n",
    "and they both give a you real A, but you can also call\n",
    "`rand_UMPS(2, 9; keep_it_real=false)`\n",
    "to get a complex A.\n",
    "\"\"\"\n",
    "function rand_UMPS(d, D; keep_it_real=true)\n",
    "    shp = (D, d, D)\n",
    "    if keep_it_real\n",
    "        A = randn(shp)\n",
    "    else\n",
    "        A_real = randn(shp)\n",
    "        A_imag = randn(shp)\n",
    "        A = complex.(A_real, A_imag) / sqrt(2)\n",
    "    end\n",
    "    return A\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm(A)\n",
    "\n",
    "Return the transfer matrix of A:\n",
    " --A---\n",
    "   |  \n",
    " --A*--\n",
    "\"\"\"\n",
    "function tm(A)\n",
    "    @tensor T[i1,i2,j1,j2] := A[i1,p,j1]*conj(A)[i2,p,j2]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_eigs_dense"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function eig_and_trunc(T, nev; by=identity, rev=false)\n",
    "    S, U = eig(T)\n",
    "    perm = sortperm(S; by=by, rev=rev)\n",
    "    S = S[perm]\n",
    "    U = U[:, perm]\n",
    "    S = S[1:nev]\n",
    "    U = U[:, 1:nev]\n",
    "    return S, U\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "    tm_eigs(A, dirn, nev)\n",
    "\n",
    "Return some of the eigenvalues and vectors of the transfer matrix of A.\n",
    "dirn should be \"L\", \"R\" or \"BOTH\", and determines which eigenvectors to return.\n",
    "nev is the number of eigenpairs to return (starting with the eigenvalues with\n",
    "largest magnitude).\n",
    "\"\"\"\n",
    "function tm_eigs_dense(A, dirn, nev)\n",
    "    T = tm(A)\n",
    "    D = size(T, 1)\n",
    "    T = reshape(T, (D^2, D^2))\n",
    "    nev = min(nev, D^2)\n",
    "    \n",
    "    result = ()\n",
    "    if dirn == \"R\" || dirn == \"BOTH\"\n",
    "        SR, UR = eig_and_trunc(T, nev; by=abs, rev=true)\n",
    "        UR = [reshape(UR[:,i], (D, D)) for i in 1:nev]\n",
    "        result = tuple(result..., SR, UR)\n",
    "    end\n",
    "    if dirn == \"L\" || dirn == \"BOTH\"\n",
    "        SL, UL = eig_and_trunc(T', nev; by=abs, rev=true)\n",
    "        UL = [reshape(UL[:,i], (D, D)) for i in 1:nev]\n",
    "        result = tuple(result..., SL, UL)\n",
    "    end\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_r"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm_l(A, x)\n",
    "\n",
    "Return y, where\n",
    "/------   /------A--\n",
    "|       = |      |  \n",
    "\\- y* -   \\- x* -A*-\n",
    "\"\"\"\n",
    "function tm_l(A, x)\n",
    "    @tensor y[i, j] := (x[a, b] * A[b, p, j]) * conj(A[a, p, i])\n",
    "    return y\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "    tm_r(A, x)\n",
    "\n",
    "Return y, where\n",
    "-- y -\\   --A-- x -\\\n",
    "      | =   |      |\n",
    "------/   --A*-----/\n",
    "\"\"\"\n",
    "function tm_r(A, x)\n",
    "    @tensor y[i, j] := A[i, p, a] * (conj(A[j, p, b]) * x[a, b])\n",
    "    return y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_eigs_sparse (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tm_eigs_sparse(A, dirn, nev)\n",
    "    if dirn == \"BOTH\"\n",
    "        SR, UR = tm_eigs_sparse(A, \"R\", nev)\n",
    "        SL, UL = tm_eigs_sparse(A, \"L\", nev)\n",
    "        return SR, UR, SL, UL\n",
    "    else\n",
    "        D = size(A, 1)\n",
    "        x = zeros(eltype(A), (D, D))\n",
    "        if dirn == \"L\"\n",
    "            f = v -> vec(tm_l(A, copy!(x, v)))\n",
    "        else\n",
    "            f = v -> vec(tm_r(A, copy!(x, v)))\n",
    "        end\n",
    "\n",
    "        fmap = LinearMap{eltype(A)}(f, D^2)\n",
    "        S, U, nconv, niter, nmult, resid = eigs(fmap, nev=nev, which=:LM, ritzvec=true)\n",
    "        U = [reshape(U[:,i], (D, D)) for i in 1:size(U, 2)]\n",
    "\n",
    "        return S, U\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_eigs (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tm_eigs(A, dirn, nev; max_dense_D=10)\n",
    "    D = size(A, 1)\n",
    "    if D <= max_dense_D || nev >= D^2\n",
    "        return tm_eigs_dense(A, dirn, nev)\n",
    "    else\n",
    "        return tm_eigs_sparse(A, dirn, nev)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "normalize!"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    normalize!(A)\n",
    "\n",
    "Normalize the UMPS defined by A, and return the dominant left and right\n",
    "eigenvectors l and r of its transfer matrix, normalized so that l'*r = 1.\n",
    "\"\"\"\n",
    "function normalize!(A)\n",
    "    SR, UR, SL, UL = tm_eigs(A, \"BOTH\", 1)\n",
    "    S1 = SR[1]\n",
    "    A ./= sqrt(S1)\n",
    "    \n",
    "    l = UL[1]\n",
    "    r = UR[1]  \n",
    "    #We need this to be 1\n",
    "    n = vec(l)'*vec(r)\n",
    "    abs_n = abs(n)\n",
    "    phase_n = abs_n/n\n",
    "    sfac = 1.0/sqrt(abs_n)\n",
    "    l .*= sfac/phase_n\n",
    "    r .*= sfac\n",
    "    return l, r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 2 solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_l_op :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n",
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_r_op :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tm_r_op"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm_l_op(A, O, x)\n",
    "\n",
    "Return y, where\n",
    "/------   /------A--\n",
    "|         |      |  \n",
    "|       = |      O  \n",
    "|         |      |  \n",
    "\\- y* -   \\- x* -A*-\n",
    "\"\"\"\n",
    "function tm_l_op(A, O, x)\n",
    "    @tensor y[i, j] := (x[a, b] * A[b, p2, j]) * (conj(A[a, p1, i]) * conj(O[p1, p2]))\n",
    "    return y\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "    tm_r_op(A, O, x)\n",
    "\n",
    "Return y, where\n",
    "-- y -\\   --A-- x -\\\n",
    "      |     |      |\n",
    "      | =   O      |\n",
    "      |     |      |\n",
    "------/   --A*-----/\n",
    "\"\"\"\n",
    "function tm_r_op(A, O, x)\n",
    "    @tensor y[i, j] := (A[i, p1, a] * O[p1, p2]) * (conj(A[j, p2, b]) * x[a, b])\n",
    "    return y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'expect_local :: NTuple{4,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "expect_local"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    expect_local(A, O, l, r)\n",
    "\n",
    "Return the expectation value of the one-site operator O for the UMPS state\n",
    "defined by the tensor A.\n",
    "\"\"\"\n",
    "function expect_local(A, O, l, r)\n",
    "    l = tm_l_op(A, O, l)\n",
    "    expectation = vec(l)'*vec(r)\n",
    "    return expectation\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'correlator_twopoint :: NTuple{6,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "correlator_twopoint"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    correlator_twopoint(A, O1, O2, m, l, r)\n",
    "\n",
    "Return the (connected) two-point correlator of operators O1 and O2 for the\n",
    "state UMPS(A), when O1 and O2 are i sites apart, where i ranges from 1 to m. In\n",
    "other words, return <O1_0 O2_i> - <O1> <O2>, for all i = 1,...,m, where the\n",
    "expectation values are with respect to the state |UMPS(A)>.\n",
    "\"\"\"\n",
    "function correlator_twopoint(A, O1, O2, m, l, r)\n",
    "    local_O1 = expect_local(A, O1, l, r)\n",
    "    local_O2 = expect_local(A, O2, l, r)\n",
    "    disconnected = local_O1 * local_O2\n",
    "    \n",
    "    l = tm_l_op(A, O1, l)\n",
    "    r = tm_r_op(A, O2, r)\n",
    "    \n",
    "    result = zeros(eltype(A), m)\n",
    "    result[1] = vec(l)'*vec(r) - disconnected\n",
    "    for i in 1:m\n",
    "        r = tm_r(A, r)\n",
    "        result[i] = vec(l)'*vec(r) - disconnected\n",
    "    end\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'correlation_length :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "correlation_length"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    correlation_length(A)\n",
    "\n",
    "Return the correlation length ξ of the UMPS defined by A. ξ = - 1/ln(|lambda[2]|),\n",
    "where lambda[2] is the eigenvalue of the MPS transfer matrix with second largest\n",
    "magnitude. (We assume here that UMPS(A) is normalized.)\n",
    "\"\"\"\n",
    "function correlation_length(A)\n",
    "    S, U = tm_eigs(A, \"L\", 2)\n",
    "    s2 = S[2]\n",
    "    ξ = -1/log(abs(s2))\n",
    "    return ξ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "expect_O1 = -0.8891572886397 - 3.885780586188048e-16im\n",
      "expect_O2 = 1.4205060807800933 + 1.1102230246251565e-16im\n",
      "  0.300390 seconds (40.52 k allocations: 95.788 MiB, 2.62% gc time)\n"
     ]
    }
   ],
   "source": [
    "d = 2\n",
    "D = 100\n",
    "A = rand_UMPS(d, D; keep_it_real=false)\n",
    "l, r = normalize!(A)\n",
    "\n",
    "O1 = randn(d, d)\n",
    "O2 = randn(d, d)\n",
    "O1 = (O1 + O1')/2\n",
    "O2 = (O2 + O2')/2\n",
    "expect_O1 = expect_local(A, O1, l, r)\n",
    "expect_O2 = expect_local(A, O2, l, r)\n",
    "@show expect_O1\n",
    "@show expect_O2\n",
    "\n",
    "m = 80\n",
    "@time corrs = correlator_twopoint(A, O1, O2, m, l, r);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.03828417702439"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ξ = correlation_length(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Pkg.add(\"PyPlot\")\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGgCAYAAACnqB1FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4lfX9//HnfU7GyQ4hEBIIAspGEmREZCphKSjg3oILvlXb0mq1/altv1bbr1VbixVxgeJGwY0KggyRHRmyZ0hIQoDsnXN+f9w5ISE7ORkneT2u61xX73Pu8bm1ysvPeH8Mh8PhQERERKQVsTR3A0RERERcTQFHREREWh0FHBEREWl1FHBERESk1VHAERERkVZHAUdERERaHQUcERERaXUUcERERKTVUcARERGRVkcBR0RERFodBRwRERFpdTyauwFNxW63k5iYSEBAAIZhNHdzREREpBYcDgeZmZlERERgsdS+X6bNBJzExEQiIyObuxkiIiJSD/Hx8XTp0qXW57eZgBMQEACYf4ECAwObuTUiIiJSGxkZGURGRpb+OV5bbSbgOIelAgMDFXBERETcTF2nl2iSsYiIiLQ6CjgiIiLS6ijgiIiISKujgCMiIiKtjgKOiIiItDoKOCIiItLqKOCIiIhIq6OAIyIiIq2OAo6IiIi0Ogo4IiIi0uoo4IiIiEiro4AjIiIirY4CTgMdSEzl4L8m89XSd5q7KSIiIlKizewm3ljsG16md9qPXJT2I/jthSseBw+v5m6WiIhIm6YenAbqMO4h3ioabx78+CK8OQnOHGneRomIiLRxCjgNFBIcxFshD3J/wW8p9AyChK3wymjYuaS5myYiItJmKeC4wNBuIXxjH8r8vguh63DIz4CP74ZPfwUF2c3dPBERkTZHAccFhnVvB8CKk95w5xcw5g+AAdsXw4KxkLSzWdsnIiLS1ijguMDQbiEA7EpIJ7sIuPyPcOfnEBABqfvh1XGwcQE4HM3bUBERkTZCAccFurTzJSLIRrHdwfbjaeaX3UfB7HXQazIU58PXD8P7t0LOmeZtrIiISBuggOMiQ7ubvTibjpYJMH7t4eb3YNI/wOoF+76E+SPh2I/N1EoREZG2QQHHRZzDVJuPnNdDYxhw6Wy4ZwW0vwgyEmDhVbD672AvboaWioiItH4KOC4SU9KDsz3+LAVF9oonhEfBfT9A9K3gsMPqZ2DRVEhPaOKWioiItH4KOC5yUUd/2vl6kldoZ1dieuUnefvDtP/CjNfAKwCOrYf5I2DvV03bWBERkVZOAcdFDMNgSFXDVOcbeD3MXgMRgyD3LLx/M3z1CBTmNUFLRUREWj8FHBca5gw4R2uxUiqkB8z6FoY/YB5vegVei4VT+xuxhSIiIm2DAo4LOVdSbT56Fru9FjVvPLxg4t/g1iXgGwrJO2HBGLNAoGrmiIiI1JsCjgv1jwjEx9NKem4hB1Kyan9hz/EwZz30GAuFOeYWDx/fA3kZtQtKIiIiUo4Cjgt5Wi1cckEwAJuOnK7bxQGd4LalMO5JMKywawm58y7jhidfYs7irZzOym+EFouIiLROCjgu5qyHs+no2bpfbLHAqLkw6xsI7opPVjzvWZ4gcs+rTHrhB1btTXFxa0VERFonBRwXG1ZmJZWjvvNoIody+Nqv+aI4Bk+jmD96vsc/C/6Xhxeu4E9Ld5JTUOTCFouIiLQ+CjguNqhrOzwsBkkZeZw4m1vut2K7gxe+28+v399OVn71IeWdnzN4oPAhFrWfi8PDhzHWHXzt/SjHNn/JVS+uY/vxevQQiYiItBEKOC7m42Xl4i5BAGwqUw8nr7CYB97dxr9XHuDTuEReX3ukynvkFRbz8bYTgEHk+NkY962Gjv3oYKSz2OsZbkh7jZvnr2V/cmbjvoyIiIibUsBpBOfXw8nIK+TONzbx9a4kDMM857W1h0nLKaj0+m92J5GWU0h4kI0xvTpCxz5w7/cw5G4A5nh8znsef2HL9u2N/zIiIiJuSAGnEZybaHyGlIw8bnzlJzYeOYO/tweL746hT6cAMvOLWLDmcKXXv7vxOAA3Do3EailJRJ4+MOV5uOFt8jwCGGQ5yPRNN8Guj5vknURERNyJAk4jGNKtHQCHT2Uz7aX17DmZQai/N+/fdykjLgrldxN6A/Dm+qOknrf8+9CpLDYeOYPFMANOBf2uZuvkz9ls74WPPRuWzIJPH4CC7EZ/LxEREXehgNMIgn296B0WAEBieh4XtPflkzmXMaCzOTcntm9HoroEkVtYzH9XHSp37fubzN6bK/p0JDzIp9L797ioDzcVPM5/imfgwIDtb8OCsZC0q/FeSkRExI0o4DSSkT1DARjQOZAlsy+ja3vf0t8MwyjtxVm88Rgn083VVvlFxSzZegKAm4d1rfLenQJtBPraeK7wOo5e9R4EhEPqfnj1Ctj0qrZ5EBGRNk8Bp5H8OrYn824ZxAf3DadDgHeF30f1DGVY9xAKiuzM+/4gAN/sTuZs6eTiDlXe2zAM+oYHArCZATB7PfSaBMX58NXv4YPbIKcWG36KiIi0Ugo4jSTQ5smUgRH4eXtU+rthGPxufC8APtgcT/yZHN4rmVx8w5BIPKzV/63pVxJwfjmZAX7t4eb3YdLfweoFe7+A+aPg2I8ufCMRERH34TYBJy0tjSFDhhAdHc2AAQN49dVXm7tJDRbToz2jeoZSZHfwyJIdbDh8GosBN1Q2ufg8fcsGHADDgEvnwD0rIORCyDgBC6+C1f8Ae3FjvoaIiEiL4zYBJyAggDVr1hAXF8fGjRt5+umnOX26jhtatkDOuTgbDpvvMrZ3RzoHVz65uKx+EWbA2XMyo/yWEOFRcP8aiLoFHHZY/TQsuhrSE1zfeBERkRbKbQKO1WrF19ecqJufn4/D4aj/Xk8tSHRkMOP7hZUeVze5uKwLO/jjaTXIzCuqsCUE3v4w/WWYvgC8/OHYOpg/AvZ+5cqmi4iItFguCzhr1qxh6tSpREREYBgGy5Ytq3DOSy+9RLdu3bDZbMTExLBp06Y6PSMtLY2oqCi6dOnCww8/TGhoqKua36zmju+Fl9XCBe19ubx31ZOLy/LysNCzo7kUfY9zmOp8UTeavTnh0ZB7Ft6/Gb56BArzXNV0ERGRFsllASc7O5uoqCjmzZtX6e8ffPABc+fO5cknn2Tbtm1ERUUxceJEUlJSSs9xzq85/5OYmAhAcHAwP//8M0eOHOHdd98lOTm5yvbk5+eTkZFR7tNS9Q0P5JvfjmbJ7MtqnFx8/nVQZh5OZdpfiOPub3Fc+ivzeNMr8HospB5oSJNFRERaNJcFnMmTJ/PUU08xY8aMSn9//vnnuffee5k5cyb9+vVj/vz5+Pr68sYbb5SeExcXx65duyp8IiIiyt0rLCyMgQMHsnbt2irb88wzzxAUFFT6iYyseeJuc+oe6lfpcvLq9A2voQcH2J+cSf+/ruJvRbfBrUvANxSSdsIrY2D7O6qZIyIirVKTzMEpKChg69atxMbGnnuwxUJsbCwbNmyo1T2SkpLIzDR3z05PT2ft2rX07t27yvMfe+wx0tPTSz/x8fENe4kWyDnRuLoenPc3xZNTUMyiDUc51Wk0zFkP3cdAYTZ8+j/wyb2Q13J7t0REROqjSQJOamoqxcXFhIWFlfs+LCyMpKSkWt3j+PHjjBo1iqioKEaNGsWDDz7IxRdfXOX53t7eBAYGlvu0Ns5aOPFncsnIK6zwu8Ph4Jvd5l/fwmKHuQ1EQCe4fSmMewIMK+z8CF4ZDQlbm7TtIiIijalZV1E5HA4Mw6jVucOGDSMuLo6ff/6ZHTt2cP/99zdy61q+YF8vIoJsAOw9mVnh950J6SSknVth9c7G4xQV28FihVG/g1nLIagrnD0Cr0+A9S+C3d5k7RcREWksTRJwQkNDsVqtFSYFp6SkVOjVkbpxTjSubB7O8l1m701s3zBC/b1Iysjju1/K/D2IHAaz10K/aWAvgu8eh3evh6xTTdJ2ERGRxtIkAcfLy4vBgwezcuXK0u/sdjsrV65k+PDhTdGEVqtswb+yHA5HacC5OjqCm4aa9XUWbTha/gY+wXD9QpjyL/CwwcEVZs2cQ6saueUiIiKNx2UBJysri7i4OOLi4gA4cuQIcXFxHD9u7q80d+5cFixYwKJFi9izZw9z5swhOzubmTNnuqoJbVJVS8UPpmRxODUbL6uFy3t34JaYrlgtBj8dPsO+pPOGswwDhsyE+1ZDh76QlQxvT4cVf4biinN7REREWjqXBZwtW7YwaNAgBg0aBJiBZtCgQTzxxBMA3HjjjTz33HM88cQTREdHExcXx/LlyzVE1UDOgLMvKdOcX1Pi65Lem1E9QwmweRIR7MOEkorJb/90tPKbdewL962CIbMAB6x7Ad6cDGerOF9ERKSFclnAGTt2bOn2CWU/CxcuLD3ngQce4NixY+Tn57Nx40ZiYmJc9fg264IQX3y9rOQX2TmSml36vXN4auKATqXf3T78AgA+2ZZQ6aorADx9YMoLcMNbYAuCE5vNncl3fdJ4LyEiIuJibrMXlVTOYjHo08ks+Occpjp+OodfTmZgtRjE9j3XQza8R3t6dvQnp6CYj7eeqP7G/a6B2esgMgbyM2DJTPjsISjIabR3ERERcRUFnFbg/IJ/zto3Md1DCPHzKj3PMAzuuKwbAG9vOIbdXkMV4+CucNdXMOr3gAHbFsGCsZC829WvICIi4lIKOK3AuaXi5uThr3edBGBymeEpp+mDOuPv7cHh1GzWH0qt+eZWDxj3ONzxKfh3gtR9sOBy2PyatnkQEZEWSwGnFXBWNP4lMYPkjDy2HU8DYEL/igHH39uD6wZ3AWDRj8dq/5AeY8xtHnpOhOJ8+PJ38MFtkHOm4S8gIiLiYgo4rUDvTgEYBqRm5fPOT2ZouaRrMGGBtkrPv+1Sc7Lx93uTiT9Thzk1fqFwywcw8RmweMLeL8wJyMdqt5+YiIhIU1HAaQV8vTzo3t4PgDfWHwVgUiXDU04XdfRn5EWh2B3w+rojdXuYYcDw/4F7VkDIhZBxAhZeCT/8H9iL6/sKIiIiLqWA00r0LZlonJVfBMCk/uHVnj9n7IUAvLvxOIll9quqtYhojlz7Fet9Y8Fhh1V/g7eugYzEut9LRETExRRwWgnnPBzn/+7a3rfa8y+7sD0x3UMoKLYzb9XBej3zne1nuPXMLBa0fwS8/OHoWnh5BOxbXq/7iYiIuIoCTitRNuBUNzzlZBgGv5vQG4APN8dz/HTd69scLbnmn8mXkDvrewiPgtwz8N6N8PWjUJRf53uKiIi4ggJOK+GshQO1CzgAw7qHMKpnKEV2By9+f6DOzzx+xqycXFBkZ92ZYLj7O7j0V+aPG1+G12IhtX69QyIiIg2hgNNKhAXa+N34Xvwmtic9O/rX+jpnL84n205w+FRWra+z2x0cK9Pr8/3eZPDwhklPwy0fgW97SNoBr4yGuHdVM0dERJqUAk4r8uC4nvwmtheGYdT6mujIYMb16YjdAf9aUftenJTMfPKLzm3uuXJPCg5niOk1AWavh+6joTAbls2BT+6D/Mwq7iYiIuJaCjjCb8f3AuDzHYnsS6pdCDl22hyeCg+y4etlJSUzn10JGedOCAyH25fBFY+DYYWdH5o1cxK2ubz9IiIi51PAEQZ0DmLygE44HPCvFftrdc2xkgKBF3X0Z1TPUABW7k0uf5LFCqN/DzO/hqBIOHsEXp8AP/4H7PbzbykiIuIyCjgCmL04hgFf70piV0J6jec7e3AuaO/LuD7mjuXf702p/OSuMTB7LfS9GuyF8O3/g3evh6xTLmu/iIhIWQo4AkCvsACujooA4IXvau7FcU4wviDEj7F9OgCw40Q6KRl5lV/g0w5ueAumvAAeNji4AuaPgMOrXdJ+ERGRshRwpNSvx/XEajFYuTeltIemKsdLhqi6tvelY4CNqMhgAFbtq6IXB8xtHobMgntXQYe+kJWM461pLH32Pj7bdtRVryEiIqKAI+f06ODPgJJ6OrsTM6o992iqGYC6leyBNa5PR8BcTVWjsH5w7/cweCYGDqZnf8AFn11HZtKhBrReRETkHAUcKadPJzPg7D1ZdcBJyykgI8/c86priLklxBUlAWftgVTyCmux6aaXL3uH/pU5Bb8mw+FLFAfwenU07F7awDcQERFRwJHz9AkPAGBPNcvFnfNvOgZ44+NlBaB/RCBhgd7kFhbz0+HTtXrWv1cc4Gt7DDO9X2CrvSfexVnw0V3w2UNQUPetI0RERJwUcKSc0h6cpKp7cJxLxC8os6GnYRhcUdNqqjJ+Sczg611JGAY8NfNKHgv8O/OKrsGBAdsWwauXQ/IvDXkVERFpwxRwpJw+ncwenPgzuWTlF1V6zrFU5xJxv3Lfx/Y9Nw/HUcPWDP9eaa7UuuricPqGB/Kr2L78s+hG7uNx7P5hcGqvGXI2v65tHkREpM4UcKScdn5ehAV6A1RZ1bi0ByfEt9z3l10YireHhYS0XPYnV72v1e7EdL7ZnYxhmCu3AKYMjODCDn58l9eH1/svhp4ToCgPvpwLH94BuWdd8XoiItJGKOBIBTUNUx0/fW6JeFk+XlZGXFRFVeMy/l2y59WUgRH0DDN7jKwWg4dKws68TWlkzFgME58Giyfs+czc5uH4Tw14KxERaUsUcKQC5zBVVT04R0+XXyJe1hU1LBfflZDOt784e28uKvebsxcnPbeQRT8eh+G/gnu+g5AekB4Pb14JPzwL9lqs0hIRkTZNAUcqcK6k2nuyYsDJLSgmJTMfKD/J2GlcyTycbcfPsu5AKnZ7+fkzzh3Lr46K4KKOAeV+K9uL89q6I2TkFULEILh/DQy8CRzFsOopeOsayDjZwLcUEZHWTAFHKnAOUe1JyqgwWdhZwTjQ5kGwr1eFa8ODfIiODMbhgNte38joZ1fx/Hf7OXY6m50n0lmxJxmLAQ9e0bPSZ08ZGEEPZy/O+qPml94BMOMVmDYfPP3g6Fp4+TLYt9x1Ly0iIq2KAo5UcGEHfzwsBpl5RZxML7+31LlNNisOTzn999ZLuHlYJAHeHpw4m8uLKw8w5tnV3Pb6RsDZe+Nf6bVWi1E68bi0F8cp+mazN6fTQMg9A+/dCF8/CkX5DXldERFphRRwpAIvDwsXdjADyPkTjUs32axkeMopItiHZ2YMZPP/i+XfN0UzqmcohgHpuYVm7824yntvnMr24ny05UT5H0MvgntWwKX/Yx5vfBlei4XUg3V8SxERac0UcKRSpRWNz5uHc+yMswen6oDjZPO0ck10Z96+O4YfH72Cx6f045Xbh5SGp6pYLQZ3XdYNgKXbT1Q8wcMbJj0Dt3wIvu0haQe8Mhri3qvFm4mISFuggCOVOrdU/LyA4+zBCal6iKoy4UE+3D2yO+P7hdXq/KsuDsfDYrArIYMDyVVsG9FrIsxeD91GQWE2LJsNn9wH+VVvMyEiIm2DAo5U6txS8boPUblCe39vxvTqAMCyuISqTwwMhzs+hSv+HxhW2PGB2ZuTuL1R2yciIi2bAo5UyjlEdehUNvlFZt2ZwmI7CWm5QPWTjF1l+iWdAVi2PbHCcvNyLFYY/TDM/AqCIuHMYXhtPPw4D+z2Rm+niIi0PAo4UqlOgTaCfDwptjs4mGJuu5CYlkux3YG3h4WOAd6N3obYvmH4e3uQkJbL5qNnar6g66Uwey30vRrshfDtn8yVVtmpjd5WERFpWRRwpFKGYVSoaOwcnuoa4ovFYjR6G2yeViYP6ATUMExVlk87uOEtuOp58LDBgW/h5RFw+IdGbKmIiLQ0CjhSJWfA2VsacGqugeNqzmGqL3acJK+wlls0GAYMvRvuXQUd+kBWkln9eOVfobiw5utFRMTtKeBIlfqEl1Q0PmlONG6qCcZlXdq9PeFBNjLzili1t/L9raoU1s8MOYPvAhyw9jlzP6uzxxqjqSIi0oIo4EiVKgxRnWn6gGOxGFwdHQHA0u21HKYqy8sXpv4brl8I3kFwYpO5M/nuZa5tqIiItCgKOFKlXmEBGAakZOZzOiuf42Xm4DSl6YPMYapV+1JIyymo8PvuxHQm/WsNz327r+qb9J9uTkDuMgzy0+GjO+Hz30BhbmM1W0REmpECjlTJz9ujNMzsTcosrWLcrQnn4IBZdLBveCCFxQ6+3Fl+F/EtR89w04Kf2JuUyfwfDnE2u2IAKtXuAnMp+ci5gAFb34QFl0PyL437AiIi0uQUcKRazmGqNQdOkVdox2ox6NzOp8nbMX1QyTDVtnPDVD/sP8Vtr28kM68IgMJiB5/WtNrK6gmxT8Idy8A/DE7tgVcvhy1vgKOaWjsiIuJWFHCkWs4tG77dnQxARLANT2vT/9/m6qjOGAZsOXaW46dz+GrnSe5ZtJm8Qjtje3fgD5P6ALBkWyV7V1Wmx1hzm4eLxkNRHnzxW/jwDsg922jvICIiTUcBR6rl7ME5kto8w1NOnYJsjLgwFIDff/QzD7y7jcJiB1MGhrPg9iHcODQST6u5d5Vz1VeN/DuYG3ZO+BtYPGHPZ+YE5OMbG/FNRESkKbhVwOnWrRsDBw4kOjqayy+/vLmb0yY4l4o7NfUE47KmlUw23nT0DHYH3Dwskn/fNAgvDwshfl6M62Nu5Llkay17cQAsFrjsAbj7W2jXHdLj4c3JsOZZsNey7o6IiLQ4bhVwAH788Ufi4uJYtWpVczelTega4ouPp7X0uCmXiJ9v0oBO+HqZbbl/TA+enn4x1jIVla8f0gWAZdsTKCyu4x5UnS+B+9fAxTeAoxi+f8osDphxsuZrRUSkxXG7gCNNy2ox6FUyTAXQNaR5hqgA/L09ePvuYbxy+2Aem9wXwyi/XcSYXh0I9ffmdHZB3YsCAtgCYcYCmPYyePrB0bUwfwTs/8ZFbyAiIk3FZQFnzZo1TJ06lYiICAzDYNmyioXUXnrpJbp164bNZiMmJoZNmzbV6RmGYTBmzBiGDh3KO++846qmSw36hJ0LON1Cm68HB2DwBSFM7N+p0t88rBZmlGzt8FFdhqnKMgyIvsXszek0EHJOw7s3wPI/QlF+fZstIiJNzGUBJzs7m6ioKObNm1fp7x988AFz587lySefZNu2bURFRTFx4kRSUs79l3Z0dDQDBgyo8ElMTARg/fr1bN26lc8++4ynn36anTt3uqr5Uo0+4WV7cJo34NTkusHmMNWqvSmkZjUgkIReBPesgJg55vFPL8Hr4+H0IRe0UkREGpvhcLi++IdhGCxdupRp06aVfhcTE8PQoUNLA5DdbicyMpIHH3yQRx99tM7PePjhh+nfvz933XVXpb/n5+eTn3/uD7iMjAwiIyNJT08nMDCw0mukcj8dPs1NC34iLNCbjX+Mbe7m1Oiaeev4+UQ6j0/px90juzf8hvuWw7I5kHsGvPzhqucg6qaG31dERGqUkZFBUFBQnf/8bpI5OAUFBWzdupXY2HN/OFosFmJjY9mwYUOt7pGdnU1mprknUlZWFt9//z39+/ev8vxnnnmGoKCg0k9kZGTDXqING9YthNljLuQvV1f917slcfbifLQlHpfk996TYM56uGAkFGTB0vvhk/shP7Ph9xYRkUbRJAEnNTWV4uJiwsLCyn0fFhZGUlJSre6RnJzMyJEjiYqK4tJLL+WOO+5g6NChVZ7/2GOPkZ6eXvqJj49v0Du0ZRaLwaOT+zBpQHhzN6VWro7qjJfVwt6kTHYn1rImTk0CI+DOz+DyP4FhgR3vwytjIDHONfcXERGX8mjOhzscjgorYarSo0cPfv7551rf29vbG29v7/o2TdxYkK8n4/uH8eWOkyzZeoIBnYNcc2OLFcY8At1Gwcf3wJlD8FosjP8rXDrHnKAsIiItQpP04ISGhmK1WklOTi73fUpKSoVeHRFXuL5kmGpZXAL5RS4u2HfBcHNn8j5TwF4I3zxmrrTKTnXtc0REpN6aJOB4eXkxePBgVq5cWfqd3W5n5cqVDB8+vCmaIG3MqJ4dCAv0Ji2nkA82N8LwpG8I3LjYnHBs9YYD38LLI+DIGtc/S0RE6sxlAScrK4u4uDji4sw5CUeOHCEuLo7jx48DMHfuXBYsWMCiRYvYs2cPc+bMITs7m5kzZ7qqCSKlrBaDm4d1BeCJT3cz98M4MvIKXfsQw4Ch98B9qyC0N2QlwaKrYeX/QnGRa58lIiJ14rJl4qtXr650f6g777yThQsXAjBv3jyeffZZkpKSiI6O5sUXXyQmJsYVj69RfZeZifsqLLbzwnf7mf/DIewO6Bzsw7PXD+Sykk07XaogB5Y/CtsWmceRMXDtaxDc1fXPEhFpQ+r753ej1MFpiRRw2q4tR88w98OfOX4mB4C7R3bn4Ym98fawkJFXRGpWPqcy8zmdVUBUZBBd2jWgmOGuT+DzX0N+BtiC4Or/QL9rXPQmIiJtjwJODRRw2rbs/CKe+vIX3ttkzscJsHmQX2SnoKj8ppw2TwsLZw7j0h7t6/+ws0dhyd2QsMU8HjwTJj0Dnj71v6eISBvVogv9iTQ3P28PnpkxkNfvHEKovzeZeUWl4SbA5kGPDn50a+9LXqGdWQs3s/nomfo/rF03mLUcRv4WMGDrm/DqFZCyxyXvIiIiNVMPjrQ52flFHDqVRTtfLzoEeGPztAKQV1jMvW9tYe2BVPy8rLx9TwyXdG3XsIcdWmVWPs5KBg8fsydn8F2qmSMiUkvqwRGpJT9vDwZ2CSYyxLc03ADYPK0suH0Iw3u0J7ugmDtf38TP8WkNe9iFl8Ps9XBRLBTlwhe/gY/ugtwG3ldERKqlgCNSho+XldfvGsKw7iFk5hdx++sb2ZWQ3rCb+neAWz6CCU+BxQN+WQbzR0H8Jtc0WkREKtAQlUglsvKLuPONTWw9dpZgX09i+4aRllPA2ZxCzuYUkJZTiKfV4K1ZMfTuFFD7GydsNScgnz0ChhUu/6M5V8dirflaEZE2SKuoaqCAI3WVmVfI7a9vIq6aYaorL+7Ef28dXLcb52XAl3Nh50fmcfcxMGMBBHRqQGtFRFonBZwaKOBIfWTkFfLRlhPkFRbTzteLdr6eBPt6kVdYzMyFmzEM+P53Y+ke6le3GzscEPcufPV7KMzWlbGpAAAgAElEQVQB31CYPh96jm+cFxERcVMKODVQwBFXm7VwM9/vTeGWmK48Pf3i+t0k9QAsmQlJO83j4Q/AuCfBw8t1DRURcWNaRSXSxGaPuRCAJVtPcCozv343Ce0Jd6+AmNnm8YZ58Pp4OH3IRa0UEWmbFHBE6mlot3YM6hpMQZGdRT8erf+NPG0w+R9w8/s4fELgZBz5L43EHve+y9oqItLWKOCI1JNhGNw/2uzFeWvDUbLyG7aDeHpkLI90+C8/2fvibc/Bsux+WDob8rNc0FoRkbZFAUekAcb3C6NHqB8ZeUW8v+l4ve+z80Q6U+at5aP9dm4p+BPPF16HHQv8/B68MhoS41zYahGR1k8BR6QBrBaD+0b3AOD1dUcoLLZXel5Vc/kdDgdvbzjKtS//SPyZXCJDfHh86gBeLJ7BLJ7EERABZw6Z83J+etlcfSUiIjXyaO4GiLi7aYM689x3+zmZnsdncYlcO7hL6W+Fxeb8nP+uPkSx3UF4kI3OwT5ElHx2JaTz5c6TAEzoF8az10cR4O3By6sPsTqzJ+tmfMqoX/4Ce7+A5Y/C4dVwzX/BrwG7nYuItAEKOCINZPO0MnNEN/5v+T5eWXOIGZd0xjAM1h1I5c+f7+Zgyrk5NOm5hexNyix3vYfF4NHJfbh7ZHeMkk04Jw3oxFsbjvH5/lxG3bgYNr8G3/wJ9i+H+SNgxqvQfVSTvqeIiDtRHRwRF0jPLWTE378nK7+Iv00fwPqDqXy1MwmA9n5ePDKpN5d0bUdCWi6JaXkkpuWSmJZLfpGdWSO7M/iC8ruWbzh0mptf/YlgX082/ykWT6sFknbBklmQug8wYPTvYcyjYNV/p4hI66VCfzVQwJHG9vRXe1iw5nDpscWAO4Z347exvQjy9azTvYrtDmKeXkFqVgFvzRrG6F4dzB8KsuHrP8D2t83jyEvh2tcgONJVryEi0qKo0J9IM5s1ojteHuY/UsO6hfDlQ6P489X96xxuwJy8PLG/uTfVVyVzdADw8oNr5sF1b4B3IMT/ZA5Z/fKZS95BRKS1UA+OiAvFxaeRllPAmF4dSufT1Nf6g6nc+tpG2pUMU3lYz/vvkbNHzZ3JE7aYx0Puhol/A0+fBj1XRKQlUQ+OSAsQHRnM2N4dGxxuAGK6hxDi58XZnEI2HjlT8YR23WDWchjxG/N4y+vw6jhI2dvgZ4uIuDsFHJEWysNqYWL/MIDSpeQVWD1h/F/g9qXg1xFSdsOCsbB1oWrmiEibpoAj0oJdeXE4AN/sSqLYXk1gufAKmLMeLhwHRbnw+a/ho7sgN61pGioi0sIo4Ii0YJf2aE+wryenswvYeOR09Sf7d4Rbl8D4/wWLB/yyDF4ZBfGbm6axIiItiAKOSAvmabUwoZ85TPV1SV2dalksMOIhmPWtOUcn7Ti8MRHWPg/2yreREBFpjRRwRFo45zDV1zUNU5XVZTDcvxYGXAeOYlj5F3h7GmTWIiSJiLQCCjgiLdxlF4YSaPMgNSufLUcrWU1VFVugWQTwmpfA0xeO/AAvj4ADKxqvsSIiLYQCjkgL5+VhYXw/s+jf17vq2ANjGDDoNrjvBwi7GHJS4Z1rzX2tigoaobUiIi2DAo6IG7jyYjPgrNiTXL8bdOgF96yAYfebxxvmwRsT4PQhF7VQRKRlUcARcQP9IszqnSfT87DXdh7O+TxtcOX/wU3vgk87SNwOr4yGHR+6sKUiIi2DAo6IG2jv5w2Ym3CezWng0FKfq2D2erhgBBRkwSf3wtI5kJ/lgpaKiLQMCjgibsDLw0K7kk07T2XlN/yGQZ3hzs9h7GNgWODnd2HBGDi5o+H3FhFpARRwRNxExwAbAKcyXRBwACxWGPso3PkFBHaG0wfhtXHw03xt8yAibk8BR8RNdAgwh6lcFnCcuo2A2eug91VQXADL/wDv3QzZNVROFhFpwRRwRNxEowUcAN8QuOkduPKfYPWG/V/D/BFwZK3rnyUi0gQUcETcRGMFnOz8Ip78dBfrD52GYffCvSshtBdknoRFU+H7v0FxkUufKSLS2BRwRNxEB/+SgOOKScZlLP7pGIs2HGPWws3ExadBp4vhvtUw6HbAAWv+DxZNgfQTLn2uiEhjUsARcRPOHpyUDNcGnC92nAQgv8jOPYu2kJCWC15+cM08uPZ18AqA4xvMbR72fO7SZ4uINBYFHBE3UTpE5cIenKOp2exMSMdiQM+O/qRm5XP3ws1k5ZcMSV18HcxeC50HQ14afHAbfPk7KMx1WRtERBqDAo6Im2iMOThf7jR7by67MJSFs4YR6u/N3qRMHnpv+7mdy0O6w8zlMOLX5vHm1+DVcXBqn8vaISLiago4Im7COQcnPbeQ/KJil9zz858TAZgyMJzOwT68ducQvD0sfL83hb99uefciR5eMP6vcNsn4NcBUnbDK2Ng6yLVzBGRFkkBR8RNBPt64mk1AEjNavhO4AdTstiblImHxWDSAHMzz+jIYJ6/IRqAN9YfYfFPx8pfdNE4mPMjXHgFFOXC5w/BklmQl97g9oiIuJICjoibMAzj3EoqFwxTfbHD7L0Z2TOUYF+v0u+vGhjO7yf0AuDJz3bz0+HzCv75d4RbP4bYv4DFA3Z/AvNHQvzmBrdJRMRVFHBE3Igr5+F8WbJ6asrAiAq//eryi5g+qDPFdgePLNlBTsF5dXAsFhj5G5j1DQRfAGnH4c1JsO4FsNsb3DYRkYZym4Czb98+oqOjSz8+Pj4sW7asuZsl0qRcFXD2JWVyICULL6uFCf3DKvxuGAZ/vaY/EUE2jp/J4blv91d+oy5DzFVWA64FexGs+DMsngGZyQ1qn4hIQ7lNwOnduzdxcXHExcWxbt06/Pz8GD9+fHM3S6RJuSrgOIenRvfqQKDNs9JzAmyePD3jYsCcj7Pt+NnKb2YLMuvlXD0PPH3h8Cpzm4eDKxrURhGRhnCbgFPWZ599xrhx4/Dz82vupog0qXPVjPPqfQ+Hw1Fa3G9qVHi1547t3ZEZl3TG4YBHluyoevWWYcAlt5sVkMMGQPYpWHwtfPv/oKjhE6JFROrKZQFnzZo1TJ06lYiICAzDqHT46KWXXqJbt27YbDZiYmLYtGlTvZ714YcfcuONNza0ySJuxxU9OLsTMziSmo23h4VxfSsOT53viSn9CPX35mBKFv9ZebCGBvaGe1bC0HvN4x//A29MhDOH691eEZH6cFnAyc7OJioqinnz5lX6+wcffMDcuXN58skn2bZtG1FRUUycOJGUlJTSc6KjoxkwYECFT2JiYuk5GRkZrF+/niuvvLLa9uTn55ORkVHuI+LuXBFwnL03V/TpiL+3R43nB/t68b/X9Afg5R8OsTuxhiXhnja46p9w4ztgC4bEbTB/NOxcUu82i4jUlcsCzuTJk3nqqaeYMWNGpb8///zz3HvvvcycOZN+/foxf/58fH19eeONN0rPiYuLY9euXRU+ERHnVnl8+umnTJw4EZvNVm17nnnmGYKCgko/kZGRrnlRkWZUuh9VPQOOOTzlLO5XcfVUVSZfHM7kAZ1KV1UVFtdipVTfKTBnPXS9DAoy4eO7YdmvID+rXm0XEamLJpmDU1BQwNatW4mNjT33YIuF2NhYNmzYUKd71XZ46rHHHiM9Pb30Ex8fX+d2i7Q0HfzNYH8qMx9HPSoI/3winRNnc/H1snJFn451uvYv1/Qn2NeT3YkZLFhTyyGnoC5w5+cw5lEwLBC3GBaMhZM76tx2EZG6aJKAk5qaSnFxMWFh5cf7w8LCSEpKqvV90tPT2bRpExMnTqzxXG9vbwIDA8t9RNydswcnv8hOZn5RDWdX9EXJ1gzj+obh42Wt07UdA2w8MaUfAM99u4/fvL+dA8mZNV9o9YDLHzODTkAEnD4Ar42Dja9omwcRaTTNuorK4XBgGEatzw8KCiI5ORkvL6+aTxZphXy8rASUzJup6zwch8NRurnmlIHVr56qyvRBnblpaCR2ByyLS2TCv9YwZ/FWdiXUYquGbiPNIaveV0JxAXz9CLx/C+ScqVdbRESq0yQBJzQ0FKvVSnJy+eJfKSkpFXp1RKR69Z1onJyRz8n0PKwWgzG9OtTr2YZh8PdrB/L5AyOZ1L8TDgd8vSuJKf9Zx6yFmzl0qob5Nb4hcNO7MPlZsHrBvq/g5RFwdF292iMiUpUmCTheXl4MHjyYlStXln5nt9tZuXIlw4cPb4omiLQaofUMOCfO5gAQHmTD5lm34anzXdwliPm3D+ab34zmmugILAZ8vzeF+9/eWvPFhgEx95nLydv3hMxEWDQVVj0DxXUfdhMRqYzLAk5WVlZppWGAI0eOEBcXx/HjxwGYO3cuCxYsYNGiRezZs4c5c+aQnZ3NzJkzXdUEkTahvj04J87mAtClnY/L2tK7UwD/vmkQK+aOwdNqcDAli6Op2bW7OHwg3P8DRN8GDjv88Hcz6KSfcFn7RKTtclnA2bJlC4MGDWLQoEGAGWgGDRrEE088AcCNN97Ic889xxNPPEF0dDRxcXEsX75cQ1QidXSumnH9enC6tPN1eZt6dPBn8AXtAFh74FTtL/Tyg2kvmVs9eAXA8R/NIau9X7q8jSLStrgs4IwdOxaHw1Hhs3DhwtJzHnjgAY4dO0Z+fj4bN24kJibGVY8XaTPq24OTkOb6HpyyRvU05/X8sD+17hdffB3MXgMRl0Bemjn5+MvfQ2H9t6QQkbbNLfeiEmnLGjpE1Tm4cQKOc+LyhkOpNRYC/GhLPL96Zxu5BWX2tgrpAbO+gcseMo83v2ouJz+1r1HaKyKtmwKOiJtp+Bwc1w9RAfQLDyTEz4vsgmK2H0+r8ry8wmL+8vkvfLnzJN/+cl4dLA8vmPC/cNvH4NcBkneZhQG3vaWaOSJSJwo4Im6mY0Dd5+DY7Q4SGmGScVkWi8HIi0IBWLO/6nk43+9NIaukSOGek1UUCrwoFmavhx5joTAHPnsQlsyCvFrU2xERQQFHxO04e3BOZ+VTbK9dr0ZqVj4FxXasFoPwoOr3cWuI0SXDVNVNNP40LqH0f/9ysppNcAPC4LalEPtnsHjA7k9g/ig4UYul6CLS5ingiLiZ9n7eWAywO+B0du16ceJLem86BdrwsDbeP/ajepo9ODsS0jmbXVDh9/ScQlbtPRd+9lQXcAAsFhj5W3NuTnBXSDsGb0yAdf8Cey02/BSRNksBR8TNWC0GIX51m4fjXCLeuZGGp5zCAm30DgvA4YB1Byuuplq++yQFxXa6h/phGGb7UzJrsVKqyxCYvQ76zwB7Eax4EhbPgMzkmq8VkTZJAUfEDdV1onFjFPmryuheZi9OZcNUy7abm31eN7gL3UP9gGrm4ZzPFgTXvQFX/wc8fODwKpg/Ag6urPlaEWlzFHBE3FBdA865GjiNs4KqLGc9nLUHUnGUWfmUlJ7HT0dOA3B1VAT9wgMB+CWxhmGqsgwDLrnDrIAcNgCyT5k9Od8+DkUVh8REpO1SwBFxQ3WtZtyUPTjDuofg7WHhZHoeB1PObb75xY5EHA4YckE7IkN86RdREnBqmodTmQ69zb2sht5rHv/4Irw5Cc4cccUriEgroIAj4obqPkRVsk1DIxX5K8vmaWVY9xAA1hw4Nw9nWcnqqWuiIwDK9ODUc+m3pw2u+ifcuBhswZCw1VxltXNJA1ovIq2FAo6IG6pLwHE4ytbAafwhKoDRJcNUzno4B1Oy2JWQgYfF4KqBJQGnpAfnSGp2+YrGddV3qjkBuetwKMiEj++GT38FBbXc9FNEWiUFHBE31LEOASc1q4D8IjsWAzo1Yg2cspz1cDYeOU1eYTGflfTejOoZSoifFwAdA2yE+ntjd8C+5FpONK5KcCTc+QWM+QNgwPbFZgXkpJ0Nu6+IuC0FHBE31KEO1Yydw1OdAm14eTTNP/K9wvzpGOBNXqGdLUfP8unP5uqpaYM6lzuvb3gAUMeJxlWxesDlf4Q7P4eACEjdD6+Og40LtM2DSBukgCPihuoyRFW6yWYTTDB2MgyjdDXVvFUHOHY6Bx9PK7F9w8qdd26isQu3YOg+yhyy6jUZivPh64fh/Vsh54zrniEiLZ4CjogbcgaczLwi8gqrn7/S2JtsVsVZD+enw2awmNA/DD9vj3Ln1GapeHpOIY99soOtx+oQUPzaw83vwaR/gNUL9n0J80fC0fV1fAsRcVcKOCJuKMDbA++S4aaaenES0kpWUDVhDw5QuvGmk3P1VFn9S3pw9iZlYq9iX61X1x7mvU3xPLJkR7m6OjUyDLh0NtyzAtpfBBkJsGgKrP472BswqVlE3IICjogbMgyjtBcnpYaA05Q1cMpq7+/NgM5mgGnn61k6ZFVW91B/bJ4WcgqKOXYmp8LvDoeDz0rm7xw6lc2Gw6fr3pDwKLjvB4i+FRx2WP0MLJoK6Qk1XysibksBR8RN1XYeTukcnOCmHaICGN+3EwDXRHfGs5JNPq0Wg96dqh6m+vlEOsfLBJ93fjpev4Z4+8O0/8KM18ArAI6tN7d52Ptl/e4nIi2eAo6Im6pNNWOHw3GuyF8T9+AAzBl7IfNvu4RHJ/ep8pzSeTiVTDT+LM7svXH2BH2zO4mUjFpszlmVgdfD7DUQMQhyz8L7t8BXD0NhA+4pIi2SAo6Im6pND86Z7ALyCu0YBoQHN00NnLK8PCxMGhCOzdNa5Tn9qlgqXmx38MUOM+D8elwvBl/QjiK7gw+3xDesUSE9YNa3MPwB83jTAngtFk7tb9h9RaRFUcARcVMdA8zAUl3AcQ5PhQXY8PaoOmQ0p6r2pNp45DQpmfkE2jwY3SuU2y7tCsB7m+IprmJCcq15eMHEv8GtS8A3FJJ3woIxZoFA1cwRaRUUcETcVG16cJprgnFd9O4UiGFAckY+p8sMt31eMrl48oBwvD2sTB4QTjtfTxLSclm1N8U1D+85Huashx5joTDH3OLh43sgz4V1eUSkWSjgiLip2lQzds6/acoif3Xl7+1Bt/Z+AOw5aW7ZUFBk56udSQBcXbK83OZp5fohkQAs3njMdQ0I6AS3LYVxT4JhhV1LzE07T2x13TNEpMkp4Ii4KWfASa2mBychreX34EDFicZrD5wiPbeQDgHeXNqjfel5twwzh6l+2H+K+EqWldebxQKj5sKsbyC4K6QdgzcmwPp/g93uuueISJNRwBFxU2WHqKoqgNdcVYzrqnQeTslEY2ftm6suDsdqMUrP6xbqx6ieoTgc8M7Gei4Zr07kULh/LfSfDvYi+O4JeOc6yHLRkJiINBkFHBE3Fepv7spdUGwnI7eo0nOac4l4XZzrwckgt6CY735JBs4NT5V126UXAPDhlnjyi2pXkfhsdgH3LNrCu7UJRT7BcN2bMPVF8PCBQyvh5RFwcGUt30ZEWgIFHBE35e1hJcjHE4CUzIp1XMwaOM4ify074PQtCTiHTmXz5c6T5BQUExniw6DI4ArnjuvTkU6BNs5kF7B8V1Kt7v+vFftZsSeZZ77eQ2FxLYacDAMG3wn3rYaO/SA7BRbPMHt0igvr8GYi0lwUcETcWHUrqdJyCskpMHs4Ilp4wAkL9CbEz4tiu4MXVx4AYOrACAzDqHCuh9XCTcPMyca1qWx8JDW7dDgrM6+In+qy3UPHPnDv9zDkbvN4/b/hjYlw5kjt7yEizUIBR8SNVVfN2Nl70zHAu9pCey2BYRilw1TOrRkqG55yumloV6wWg01Hz1S7EznAs9/spcjuwJmVvt2dXLfGefrAlOfhhrfBFgQJW+GV0bDr47rdR0SalAKOiBurrgfHXebfODknGgP0CvOnT6fAKs/tFGRjUn9zn6u5H8aRW1D5XJxtx8/y1c4kLAY8OsncLuLbX5Kq3LncafFPx/jN+9vJzCszHNXvapi9HiIvhfwMWDILPn0ACrJr+4oi0oQUcETcWMdqdhQvnX/TwldQOTl7cACujqq698bpian9CPX3Zm9SJn9atrPCSjKHw8EzX+0B4PrBkdw1oht+XlaSM/LZkVB1Ib+z2QX89YtfWBaXyN+/3lv+x+BIuOtLGP0IYMD2t2HBWEjaVev3FJGmoYAj4sa6tjfDyze7kyqsKHKXGjhOZXtwptYi4IQF2ph3yyCsFoNPtiXw3qbye1R990sym4+exeZp4bfje+HtYWVsn46A+derKh9vO0FBkTkR+Z2Nx9lw6Lw5O1YPuOJPcOdnEBAOqfvh1Stg06va5kGkBVHAEXFjMy7pQscAb46dzuHN9UfL/eZuQ1Q9O/pz98ju/G58Ly4oqWxck0t7tOfhib0B+PNnu9lxIg2AomI7f19u9r7cM7IHnYLMfbsmlgxrfVtFwHE4HKVLySNDzL9uj36yo/IhsO6jzSGrXpOgOB+++j18cBvknKnlG4tIY1LAEXFj/t4e/KFkbsl/Vh4gJePccnF3KfLnZBgGj0/px4PjetbpuvtH92BCvzAKiu3MWbyNs9kFfLAlnsOnsgnx8+L+MT1Kzx3buwOeVoNDp7I5mJJV4V4bDp/mcGo2/t4efHT/ZXQKtHHsdA7Pf7ev8of7tYeb34dJfwerF+z9wtzm4diPdXoHEXE9BRwRNzd9UGeiIoPJLijm/74x/yB2pxo4DWUYBv+8IYpu7X1JSMvlofe388J35lLzX4/rSYDNs/TcQJsnl10YCpiTjc/nXE5+TXQEnYJsPD1jAACvrztCXHxaVQ2AS+fAPSsg5ELIOAELr4LV/wB77QoRiojrKeCIuDmLxeDPU/sBsGTrCX6OTyMjt4isfLO6sbsMUTVEoM2T/946GG8PC2sPpJKalU+39r7cXLJ3VVkT+ocB8M15y8VPZeaXDl3dGmNWS76iTxjToiOwO+CRJT+Xzs0py+FwkJiWS3HYQLh/DUTdAg47rH4aFl0N6Qmufl0RqQUFHJFWYFDXdsy4pDMAf/58N/El829C/Vt+DRxX6RcRyN+mX1x6/MikPnh5VPxX3Ph+YRgG/ByfRlL6uSG9j7bGU1jsIDoyuNyE5yem9qe9nxf7k7N4adXB0u9Ppufy0qqDXPHcD1z29+/NYSxvf5j+MkxfAF7+cGwdzB8Be79qpLcWkap4NHcDRMQ1/jCpD8t3JbH9eBrzvjf/IG4LvTdlXTe4C5l5hWTkFjF5QKdKz+kYYGNQZDDbjqfx3S9J3D68G3a7g/c2mcNTt8aU7/UJ8fPiL9f054F3t/PSqoME+3qyat8p1h04RdlyOu9sPM5D43ri7WGFqBuhyxCzVs7JOHj/Zhh2P4z/K3jaGu39ReQc9eCItBJhgTZ+dflFACwvGWppawEHYOaI7vw6tmel2zw4la6mKtnUc+3BVOLP5BJg82DKwIpL1K+6OJwJ/cIosjv4y+e/sGa/GW6GdQ/h2esGEhboTVpOIav2njp3UfsL4e7vYPgD5vGmV+D1WEg94LqXFZEqKeCItCJ3j+xO15Bzq6Y6t8GAUxsTSgLOhkOnSc8p5J2fjgFw7SVd8PGqOKRnGAZPTRtA52AfIoJsPHD5Raz+/Vg+vH841w+JZNogc3jwk20nyl/o4QUT/wa3LgHfUEjaaW7zsH2xauaINDIFHJFWxOZp5U9X9S09dpcl4k2te6gfPTv6U2R38N7m46zcmwJUHJ4qq2OgjfWPXsGPj43j9xN70y30XK2eGYO6ALBqXwpnsgsqXtxzPMxZD93HQGEOfPor+OReyKt+Hy0RqT8FHJFWZkK/MGL7dsQwYMgF7Zq7OS2Wc5jq+W/3U2x3MKxbCD3DAup1r96dAhjQOZDCYgdf7Eis/KSATnD7Uhj3BBhW2PmR2ZuTsLW+ryAi1VDAEWllDMNg/m2D2fjYOPqGV71hZVvnXC5eUGwu/b6lmt6b2phe0ovzybZqloVbrDDqdzBrOQR1hbNH4PUJsP5FsFdcgi4i9edWAeef//wn/fv3Z8CAASxevLi5myPSYnlYLXQM1Gqd6lzcOYjwki0c2vl6MqmKVVe1dXVUBFaLQVx8GodOVaySXE7kMJi9FvpNA3sRfPc4vHMdZKU0qA0ico7bBJydO3fy7rvvsnXrVjZv3sx//vMf0tKqqCwqIlIDwzBKN/W8aVjXBtcL6hDgzZheHQBYWl0vjpNPMFy/EKb8CzxscGglvDwCDn3foHaIiMltAs6ePXu47LLLsNls+Pj4EB0dzfLly5u7WSLixuaO78WC2wczd3wvl9zPWWxx6fYE7PZarJIyDBgyE+5bDR36QnYKvD0DVvwZigtd0iaRtsplAWfNmjVMnTqViIgIDMNg2bJlFc556aWX6NatGzabjZiYGDZt2lTr+w8YMIBVq1aRlpbG2bNnWb16NQkJKoEuIvVn87QyoX8nPK2u+VdhbN8wArw9SEjLZeOROuwq3rEv3LcKhswCHLDuBXhjEpw96pJ2ibRFLgs42dnZREVFMW/evEp//+CDD5g7dy5PPvkk27ZtIyoqiokTJ5KScm7MOTo6mgEDBlT4JCYm0q9fPx566CGuuOIKZsyYQUxMDB4eKsQsIi2HzdPKVQPDAVi6/UQNZ5/H0wemvAA3vAW2IEjYYu5MvuvjRmipSOtnOByurzZlGAZLly5l2rRppd/FxMQwdOjQ0gBkt9uJjIzkwQcf5NFHH63zM+655x6mT5/OVVddVenv+fn55Ofnlx5nZGQQGRlJeno6gYFaWSIijWPTkTPc8MoG/L092Pyn2EoLB9Yo7Th8fA/EbzSPL7kDJv0DvFTXSNqejIwMgoKC6vznd5PMwSkoKGDr1q3Exsaee7DFQmxsLBs2bKj1fZy9Pfv27WPTpk1MnDixynOfeeYZgoKCSj+RkZH1fwERkVoackE7IkN8yMov4ttfkup3k+CucNdXMOr3gAHb3oIFYyF5tyubKtKqNUnASU1Npbi4mLCwsHLfh4WFkZRU+38BTJs2jX79+nHbbbfx5ptvVjtE9dhjj5Genl76iY+Pr3f7RURqy2IxmB7t3LqhAdcm1aEAACAASURBVPMErR4w7nG441Pw7wSp+2DB5bD5NW3zIFILzTqJxeFwVLsh3vl+/PHHWp/r7e2Nt7d3fZolItIg0y/pwovfH2TtgVOkZOQ1rCZRjzHmNg/L/gcOfAP/v737Do+qSh84/r0z6aQnEBJCQjX0EEAwgBSJi8oqoLioqJQVBVFRdl1RVyyr4i6rP11xYdGluHZUUCyIhCqCSEmQ3gIJkEIC6X3m/P64mSG9h0wm7+d58kjmnnvnHAeT13Pe855v/wSnNsNtb4Obb+N1Wgg7c1VmcPz9/TEajSQnJ5d5PSUlpcKsjhBCtHSd/dswMNQHs4LXNxxv+APb+MM9n8LYhWBwhKPf6AnIZ2u/xC9Ea3NVAhwnJycGDhxIdHS09TWz2Ux0dDSRkZFXowtCCHFVPXVTDzQNPt2TwNbjFxv+QE2DyIfhgY3g2xUyz8HKW2DrP8BsavjzhbAzjRbgZGdnExMTQ0xMDABxcXHExMQQHx8PwLx581i2bBmrVq3iyJEjzJ49m5ycHKZPn95YXRBCCJsxuLMv04Z2AmD+FwfIzK++cJ+pNoUBAYL6w0NbIfxuUGbY/Aq8Px4yqzjkU4hWqtECnD179hAREUFERASgBzQREREsWLAAgMmTJ/P666+zYMEC+vfvT0xMDOvXr5clKiGE3XpybBihfm4kZuTz6rdHqmz3Vcx5+r+0gRe+ruUuKWcPmLgUJi4DJ3c4s10/5uHY943UcyFaviapg2OL6ruPXgghGuKX02lMXrYLgPdnDGZEyXlVoG+0WLr1NH9ff9T62pqHhxIR4lP7N0g7BZ9Ph8RY/fshs+DGl8BBNlkI+2DTdXCEEKK1GtLFr8xSVVbJUpXJrFjw1SFrcBPs4wrAq98doU7/3+nXFf74I1w3R//+l6Xw3hhIPdFoYxCiJZIARwghmthfbgojxNeNCxn5vPrdUfIKTcz6YC//23UWTYPnft+L1bMicXE08OuZy2w4nFzzQ0tzcIabXoV7VoObHyT9Bv8ZCTEfSc0c0WpJgCOEEE3MzcmBv9/RD4CPd8dz6+Kf+PFwMk4OBhbfPYA/Du9MoJcrDwzvAsBr3x+lyGSu+xtd8zuYtQM6j4CiHFg7G758EPIzG3M4QrQIEuAIIcRVENnVj6mRoQCcTMnGy9WRD/44xHo4J8BDI7vg18aJuNQcPt4dX7838gyE+9bCDc+BZoTfPoP/jIDz+xpjGEK0GBLgCCHEVfKXm3oQ3tGbawLc+WJ2JIM7l61E7OHiyOM3XgPAWxtPWPN16sxghBF/hunfg1dHuBwH//0d/Pw2mOsxMyRECyS7qIQQ4ioymxWaRpXH1BSZzIx9cxunL+YwZ3RXnhzbo2FvmHcZvn4Mjnytf98tCiYsBfe21d8nhI2QXVRCCNECGAxatWfwORoNzL9JD2re2x7HhfS8Wj1XKcWKHXEVqya7+sAf3off/x84uMDJjbB0mH6elRB2TAIcIYSwMTf2CmBwZ18Kis21Psvq0IVMXlx3mJmr9nAiOavsRU2DQTNg5mZo2xOyk+F/E2HjC2Cq5zKYEDZOAhwhhLAxmqbx7C09Afhy/zkOX6h5F9TZtFwACk1mnvz8QOVHPwT0gpmbYOB0QMFP/wcrbobLZxuz+0LYBAlwhBDCBoV39GZs7wCUgvUHE2tsX3opKyYhnRU74ipv6OQGt74Jd64CZy8496t+MvmhNY3VdSFsggQ4Qghho8I7egNw7nLNeTjnSwKcEF83AP654Rhn03KqvqH3BJj9E3QcAgUZsHqanoxcmNvgfgthCyTAEUIIG9XBWz++4VwtEo0tAc4D13dmaFc/8ovMPPXFAczVnVLuHQLTvoPr/wxosG8VvDsakmt56KcQNkwCHCGEsFHBPvpszPlazOBYlqg6eLvy2u39cHE0sOv0JT7+tYaCgUYHGPMc3P8VuLeHi0fh3Rvg1/fkmAfRokmAI4QQNspyAGdSZj7FNRzdYAlwgrxdCfFz48+/CwNg4XdHa7fVvMtImL0Duv8OivPh2z/BZ/fpdXSEaIEkwBFCCBvV1t0ZJ6MBk1mRmJFfZbvcwmIu5+rbvTuUBEXTh3UmIsSb7IJinl3zW+1OKG/jD3d/CmNfBYMjHFmnJyDH72qU8QhxNUmAI4QQNspg0KwBy/lqZmEsMzQezg54ujgCYDRo/OOOfjgZDWw+dpFvf6t5J1bJm0LkHHjgR/DtAhkJsOIW2LoIzKaGDUiIq0gCHCGEsGHWRONq8nDOp+uzO0ElbS26B3gwfVgnAH48nFy3Nw6KgIe2Qb+7QJlg88vw/njIvFC35wjRTCTAEUIIG2bJwzl3uert25YkZMtsT2l9OngB1PrIhzKcPeD2/+hnVzm2gTPbYckwOLa+7s8S4iqTAEcIIWyYZQanup1UVxKMXSreXxL0XEivOoenRv3v1mdz2veDvEvw8WT4fj4UF9T/mUI0MQlwhBDChgX71maJyrJF3K3CNUuAVJudWNXy7wYPbITrHta//2UJvBcFqSfr/0whmpAEOEIIYcMsQUt1Scbnq5nBKb0TKymzAbM4AA7OcNNCuOczcPODpAPwnxEQ83HDnitEE5AARwghbFiwdYkpr/IDNClb5K88g0EjsCTwadAyVWnXjIVZO6DT9VCUA2tnwZcPQkFWzfcKcZVIgCOEEDYswNMFB4NGsVmRXMkMjMmsSCqpkVNZkjFAkJdlq3kjnjPlGahXP77hr6AZ4cCn+mzOhf2N9x5CNIAEOEIIYcOMpWZgKlumSsnKp9isMBo02nlUXKKCRko0rozBCCOehOnfgVdHuHQa3rsRfl4M5gbk+wjRCCTAEUIIGxdckodT2VZxy/JUe08XjAat0vuDalFLp0FCroNZ26HnbWAugg3Pwkd/gOyLTfN+QtSCBDhCCGHjrNWMKwlQzlVTA8ci2PtKHk+TcfWBP7wP494ABxc4+SMsHQantzTdewpRDQlwhBDCxl0p9lcxQLEsO1WWYGxhmcGpbidWo9A0uPaPMHMztO0B2cnw/gSIfglMRU373kKUIwGOEELYuA7VBCjVFfmz3l9qJ1atDt1sqIBeepAzcBqgYPvr+nlWl882/XsLUUICHCGEsHHBPpYcnIoBTnVF/iwCvfTgJ7fQRHruVZpJcXKDW9+CO1eCsxec262fTH5ozdV5f9HqSYAjhBA2LrhUDo65XC2c2szguDga8Xd31p/R1MtU5fWeqCcgBw+GggxYPQ3WzYXCRtyyLkQlJMARQggb197LBYMGhSYzqdllz3+yBCzB1SQZQ6lE5asd4AD4hOpbyYfPAzTYuxLevQGSD1/9vohWQwIcIYSwcY5GA4ElxfoSSi1TZeYXkZVfDGC9XpUOllo6TbVVvCZGR4h6Hu5fC+4BcPEIvDsa9iyHq5EXJFodCXCEEKIFqCzR2LI85e3mSBtnh1rd36RbxWujyyj9mIduN0JxPnzzBHx2P+Rdbt5+CbsjAY4QQrQAV7aKX8ldqe4MqvKu2lbx2nBvqx/Y+btXwOAIR77WE5DjdzV3z4QdkQBHCCFagA6V1MKxLDcF1SLAsZkZHAuDAYY+An/cAD6dISNB30q+bRGYTc3dO2EHJMARQogWILiSasbna1Hkz8KmZnBK6zAAHtoGff8AygSbXob3x0NmYnP3TLRwEuAIIUQLcKUWTv2WqCwBUmp2IflFNjZD4uIJty+DCUvAsQ2c2a4f83D8h+bumWjBJMARQogWoHSSsaUa8fn02i9Rebk64uZkBKpeplJK8caGY7yz+WRjdLluNA3636PP5rTvB7lp+oGd65+G4oKa7xeiHAlwhBCiBQj0dkHTIL/ITFpOIVBqBqeGGjgAmqZVe+QDwOnUHP616SSLfjjGyZTsRup5Hfl3gwc2wpDZ+ve7/g3/vRHSTjVPf0SLJQGOEEK0AM4ORtp5lFQjvpxHkclMcqaeg1NdFePSSp9JVZl9Z69s1f469kJDutswDs5w82tw96fg6guJsfCfERD7SfP1SbQ4EuAIIUQLUfpMqqSMfMwKnIwG/Ns41+p+a6JxFcX+9iekW//8Vcz5q3MwZ3XCboLZO6DT9VCYDWsegi8fgoKs5u2XaBEkwBFCiBbCssR07nJumTOoDAatTvdbdl+VV3oG52xaLrHnMhrS3cbhGQT3fwWj/wqaAQ58os/mXNjf3D0TNk4CHCGEaCGCS50nVZcEY4srAU7Fgy6zC4o5nqzPjFzXxRfQZ3Ga2oX0PF5ad5gzqTlVNzIYYeSTMO078AyGS6fhvRth5ztyzIOokgQ4QgjRQpReorpQnwCnmgM3DySkY1Z6EPTgiC4ArItNxGRu2gDi7+uPsnxHHNNW7CYjr6j6xqGR+snkPX4P5iL44Rl9p1VOapP2UbRMNhngTJw4ER8fHyZNmlSna0IIYc86lCr2V5cifxaWYCgpI79C4LIvXl+eigjx5vrubfFxcyQ1u4Cdp9Iao+uVyswvYv3BJADOpOXyp89iMNcUULn5wuQPYNzrYHSGExtgyTA4vbXJ+ilaJpsMcB577DHef//9Ol8TQgh7Vvo8qvN1KPJnEeDhjNGgUWRSXMwqW1tmX7yeYDwgxAdHo4Fb+gYCTbtM9e2BRAqKzbT3dMHJwcDGIym1q8GjaXDtA/DgZvAPg+wkvfpx9EtgKm6y/oqWxSYDnNGjR+Ph4VHna0IIYc8swUxOoYnDFzL112pRA8fCwWigvae+pbz0MpVSiv2lZnAAxvfvAMD6g0lNVvn4873nAJg+rBMvj+8DwBsbj7PlWErtHhDQGx7cAgOmAgq2vw4rb4H0+Cbpr2hZ6hzgbNu2jVtvvZWgoCA0TWPt2rUV2rzzzjt06tQJFxcXhgwZwu7duxuls0II0Zq5OBrxd9e3hKdm6zMwdcnBASot9ncmLZfLuUU4ORjoHeQFwKBQH4K8XMgqKGbz0VoGHHVw+mI2e89exqDBxIgO/OHajtw9OASlYO4nMSRcqpgIXSknN7jtXzBpBTh7QsIvsHQ4HP6q0fssWpY6Bzg5OTmEh4ezePHiSq9/+umnzJs3j+eff559+/YRHh7O2LFjSUm58h9I//796dOnT4WvCxcar7BUQUEBmZmZZb6EEKKlKz9jE+hVuyJ/5e8vXezPsj28bwcvnBz0XwsGg8at/YMA+Cqm7j+ba8ql+WKfPnsz8pq2tCuZVXrhtl6Ed/QmI6+Ih/63t24zR31u1xOQOwyC/Az47H5Y9zgU2djhouKqqXOAc/PNN/Pyyy9z++23V3r9jTfeYObMmUyfPp1evXqxdOlS3NzcWL58ubVNTEwMBw8erPAVFBRU/5GUs3DhQry8vKxfHTt2bLRnCyFEcwkuFeD4uzvj4mis0/2Wqseli/1ZEowHlCxPWYwP15epNh1LqXmHU4ncwmKmvLeLsW9uIyWr8no7JrPiy316bs+kgVd+Njs7GFkyZQB+bZw4nJjJk58f4NCFjNoHOj6dYMZ6GP4EoMHeFbBsNKQcqd39wq40ag5OYWEhe/fuJSoq6sobGAxERUWxc+fOxnyrGj399NNkZGRYvxISEq7q+wshRFMoHeB0qOURDaV18Na3mpdeotpfKsG4tJ6BHnRv505hsZkfDiXV+GyzWTHv01h2nEzjREo2z609WGk15J9PpZKYkY+XqyNjerYrcy3I25W3747AoMG62AuM+9dP9Fywnuv/sYnpK3bzyreH2XPmUtWdMDpC1Atw3xpwD4CLR2DZKNizXGrmtDKNGuCkpqZiMpkICAgo83pAQABJSTX/x2ERFRXFnXfeyXfffUdwcHCZ4Ki6a6U5Ozvj6elZ5ksIIVq64FI5N3XNv9Hv0YMiyxJVTkExR5P0JfyIcgGOpmlMiNBncb6uxTLVmxuPs/5QEk5GAw4GjR8OJfPNgcQK7SzJxbeFB1U6AzW0mz9v3z2AQaE+eLs5ohQkXMpj87GLvLs9jqnLd5NdUMNuqa6jYdYO6BYFxfnwzROweirkpVd/n7AbDlfjTZRSaFrtSokDbNy4sV7XhBDC3lmK/UHdtohfub/seVSx5/QCf0FeLrSvJJ/ntvAgFv1wjJ9PpZKSmW/NlylvXewF/rVJ3+L96u19SbiUy1vRJ1jw1UEiu/pZk6NL176ZNDC4yn6O6xfIuH6BKKVIyynkZEo2J1OyWbLlFOfT8/juQCJ/uLaG1AP3tnDPatj1Dmx8QU88Pr8fJv0XOg6u/l7R4jXqDI6/vz9Go5Hk5OQyr6ekpFSY1RFCCFF3pZOM6zeDo9+TVVBMZn6RdXkqItSn0vYdfd0YEOKNWcGqnWcqrWx84Fw6f14dC8CDI7owaWAwc0Z3o0d7Dy7nFvH8V4esbS21b7q3c6dfsFeN/dU0DX93Z67r4se914Vy73WhAKzeW8u0A4MBhj4Kf9wAPp0hIx6W3wTb/gnmptn+LmxDowY4Tk5ODBw4kOjoaOtrZrOZ6OhoIiMjG/OthBCiVerQwCUqNycHfNwcAX0WZ781wbjyAAf0bdwA72w+xYh/bGbp1lNczikEIDkzn5nv76Gg2MzosLY8dVMPAJwcDPzzznCMBo1vf0vku9/0pSrL8tSkgcF1mtm3uH1ABwwa/HrmMnHVnV9VXoeB8NA26HsnKBNs+hv8bwJkVlxCE/ahzgFOdnY2MTExxMTEABAXF0dMTAzx8XphpXnz5rFs2TJWrVrFkSNHmD17Njk5OUyfPr1xey6EEK1QG2cH2nnoyz2d/N1qaF05S2B0/nKetYJxRLkdVKXdMySUR0Z3w9vNkfPpebz2/VGGLIzmz6tjmfn+HpIzC+jezp1/3R2BsdTJ5n06eDF7ZFcAnlt7kL1nL5WpfVMfAZ4ujLimLQBflARLtebiCbe/C+P/DY5uELcNlg6D4xvq1Rdh2+oc4OzZs4eIiAgiIiIAPaCJiIhgwYIFAEyePJnXX3+dBQsW0L9/f2JiYli/fr0sUQkhRCP5553hLPh9L8IC6lfV3TILtPN0GpdyCnEyGugdVPVGDKNB489jw9j19Bj+MakfvYM8KSw28/necxw4l4G3myPvTR2Eh4tjhXsfHdONawLcScspZNryX4GytW/q486SreVf7DtX98NANQ0ipuizOe37Qm4afHQn/PAsFBfWu0/C9miqsj18digzMxMvLy8yMjJkR5UQolV7cd0hVuw4Q3tPF5Iy8xkQ4s2XDw+r9f1KKfbFp/P+zjP8dj6DhRP7MqSLX5XtYxPSmfjvHVhikXfuGcC4foH17n9+kYkhr0aTkVfE+zMGW2d0yisymTmTmkP3qgLBonzY+Dz8slT/PrA/TFoOfl3r3TfR+Or7+9smz6ISQgjRdCwzOEmZeiG+6vJvKqNpGgNDfXjrrgg2/WlUtcENQHhHbx4coQcNni4OFWrf1JWLo5HxJVWWP69imcpkVtz/393c+H/b+Mf6o5U/yNEFbv473P0JuPpCYgz8ZwTEftqg/gnbIAGOEEK0MuW3l5evf9MUHo/qzqM3dOOtuyLqXH25MpZlqh8OJVVaZXnJlpPsPJ0GwL+3nOK97aerfljYzTB7B4QOh8JsWPMgrJkFBdkN7qdoPhLgCCFEK1N+99WA0KoTjBuLi6ORP/0ujNE9GjZ7Y9GngydhAR4UFJtZF1u2COG++Mv838YTAIwK05evXv72CGv2V5OU7BkEU7+GUc+AZoDYj/XZnAsxjdJfcfVJgCOEEK1M6Vo6gV4uBHrVfbt5c9M0jTsH6YUCSy9TZeUX8fgnMZjMitvCg1gx7VpmDOsMwJOrD1R/MrrBCKOegmnfgmcHuHQK/nsj7Py3HPPQAkmAI4QQrYxfGyecS04Nr2v+jS0Z378DRoNGTEI6J1OyAHj+q0PEX8ol2MeVlyf2QdM0/jquJxP6B1FsVsz+cC97S05Pr1LoUJj1E/T4PZgK4Yen4aPJkJN6FUYlGosEOEII0cpommbNw6mu/o2ta+vhzOgwfclr9Z5zrN1/ni/3n8egwVt39cezZNu6waCx6M5wRoW1Jb/IzIyVv3I8Oav6h7v5wuQP4JZ/gtEZTvwAS4bptXNEiyABjhBCtEI39g7Ay9WRG3u17BplpZep/rr2IACPjenOwFDfMu0cjQb+PWUAA0K8ycgrYvqKXykoruGoBk2DwTNh5ibwD4PsJFh1G2x6GUw1HPYpmp0EOEII0Qo9fXNP9j93I6F+bZq7Kw0yOqwdvm2cSMspJLugmEGhPjwyululbd2cHFg+7Vq8XPWKzEcTa5jFsWjfBx7cDBH3AQq2LYKV4yA9vvEGIhqdBDhCCNFKGQx1PwvK1jg5GJjQXz/2wcPFgTfv6o+Dsepfbd5uTtZDPg9dyKzDG7WB8Yv1QoDOnpCwC5YOh8NfN6j/oulIgCOEEKJFmzWqC7eGB7H03oEE+9R8PlfvID3AOXgho+5v1ucOmLUdOgyC/Az47D745gkoyqv7s0STkgBHCCFEi9bOw4W3745gWDf/WrW3nLtVpxmc0nw6wYz1MOxx/fs9y+HdGyDlSP2eJ5qEBDhCCCFalT4d9Bmco4mZFJvM9XuI0RFufBHuWwNt2kHKYVg2GvaulJo5NkICHCGEEK1KqK8b7s4OFBSbOXUxp2EP63qDfsxD1zFQnAfr5sLqaZCX3ih9FfUnAY4QQohWxWDQ6BVoWaaqRx5Oee7tYMrncOPfwOAAh9fC0ushYXfDny3qTQIcIYQQrU6vhubhlGcwwLDHYMYGPUcnIx6W3wTbXwdzPZfBRINIgCOEEKLVsSQaHzzfCDM4JUxmxUv7Xfl04EfQZxIoE0S/BP+bAFlJjfY+onYkwBFCCNHqWBKND1/IxGxunKTgX+LSWL4jjufWx1M4fhmMfwcc3SBuq37Mw4kfG+V9RO1IgCOEEKLV6dbOHScHA1kFxSRczq2y3c8nU7nlre0cOFdz0vAvpy8BUFhs5nhKNkTcCw9uhYC+kJsKH06CH56F4sJGG4eomgQ4QgghWh1Ho4GwAA+g+jyctzed5HBiJst/iqvxmbvjLln/fOBcydJX22vggY0w+CH9+52LYfnvIO1U/TsvakUCHCGEEK1Snw7V76RKzy1k9xk9aPnpZBqqmvo2BcUm9sVftn5fZsbH0QVu+Qfc9RG4+sCF/fCfEXDgs0YYhaiKBDhCCCFapV6WIxvOVz6Ds/lYCqaS/JzU7AKOJ2dX+awD5zIoKDaX+b6CHuNg1g4IHQaF2fDlTFgzGwqqfq6oPwlwhBBCtEp9gq7M4FQ2O7PhUHKZ7386mVrls345nQZA/47eABxLziK/yFSxoVcHmLoORj0NmgFiP4JlIyExtr7DEFWQAEcIIUSr1KO9JwYNUrMLSckqKHMtv8jE1uMXARjXLxCAHdUFOCX5NxP6B+Hv7oTJrKrO7TEYYdR8Ysd8QLpDW0g7Ce9Fwa4lcsxDI5IARwghRKvk6mSka1t3oGIezs5TaeQWmmjv6cLskV0B2HU6jcLiikX7ikxm9p7V82+GdPGjX7A+i/NbDTuvntjVhlHZL3M+4AYwFcL6+fDx3ZCT1uCxCQlwhBBCtGKWejiHyuXhbDisL09F9WpHr0BP/No4kVtoIiahYtBy8HwGuYUmvN0cCQvwoF+w/sxK83BKJFzK5XRqDul48EmXhXDLP8HoDMe/h6XDIG57Yw2x1ZIARwghRKtlrWhcagbHbFZsPKIHODf2ao/BoDG0mz9QeR6OZXnq2k6+GAzalQCnmirJ209ceU785TwYPBNmRoP/NZCVCKtuhU2vgKm4gSNsvSTAEUII0WpVdiZV7Ll0LmYV4O7swHVdfAEY3s0PqDwPx1L/ZkhnvW3fDvoS1amL2WQXVB6gbD9x0frns2klhQbb94UHt0DEfYCCbf+AVb+H9IR6j681kwBHCCFEq9W7ZKv4uct5pOfqFYYty1Ojwtri7GAEYFjJDE5MQjpZ+UXW+01mxa/WAEcPgtp6OBPk5YJS8Fsly1TFJnOZQCnhUqlKyk5tYPxiuOO/4OQB8Tth6XA4sq6xhtxqSIAjhBCi1fJydaSjryugn0sF8ONhy/JUgLVdsI8bnfzcMJmV9UgGgCOJmWQVFOPh7GCdDQKuJBqfr5izc+B8Bpn5xbRx0oOntJzCijM9fSfBrO3QYSDkp8On98K3f4KivEYYdesgAY4QQohWrU/JLM6hC5nEpeZwMiUbB4PGqLB2ZdoNqyQPZ1dJ/ZtBnXwwGjTr631L8nBiK5nB2Vay/XxkWFt82zgBcDYtp2LHfDvD9PUwbK7+/a/vwbtjIOVofYbZ6kiAI4QQolUrnWj84+EkAK7r4oeXq2OZdtd3rxjgWPJvBpcsT1mEW7eKVwxwLAnG13dvS4ivG1Bumao0Bye48SW490to0xZSDsGyUbB3ldTMqYEEOEIIIVq13qVmcCzLU7/rHVChXWQXfzQNTqZkk5SRj9msrGdVDSlJRrboW7L9PP5SLpdzrpwenpFXZN1qfn13f0L99ADHmmhclW5jYPbP0PUGKM6DdY/B59Mhv+qdWq2dBDhCCCFatd4lh26evphtLdgX1bNigOPl5ki/ksBlx8lUjqdkkZ5bhJuT0RrQlG7bqSR4Kb1dfOepVExmRZe2bQj2cbPO4JytaganNPd2MOULiHoRDA5waI2egJzwa90H3QpIgCOEEKJVa+fhQlsPZ8wKzEo/ZTzI27XStpY8nB0nU63JxgNDfXA0Vvx1WllF420ly1MjurcFqHmJqjyDAYY/DjN+AO9QSI+HFTfB9jfAXLHKcmsmAY4QQohWr3epHVA39mxfZbvhpRKNf4nTE4wHd/KttG2/conGSilrgrEln8c6g1PTElV5wYP0XVZ97gBzMUS/CB9MhKzkmu9tJSTAEUII0epZdlJB5fk3FgNCfXBxNJCSVcDGIymAfv5UZfqVSzQ+k5bLuct5OBo1riu5J9SvDQDn0/MoMtVxBsbFS6+Xc9ticHSD01tgyVA4wcSNMAAAFY5JREFUsbFuz7FTEuAIIYRo9Szbujv6utKjvUeV7VwcjVxbMmNTWGzG2cFAeEevStv2DtJPK0/KzCclM99avXhgqA9tnB0AaOfhjLODAZNZkZieX/eOaxoMuE+vgBzQB3JT4cM7YMNfobiwprvtmgQ4QgghWr2ongE8OTaMNydHoGlatW0ty1QAESHe1mrH5bVxdqBbO/208gPnMth2/Mr2cAuDQaOjNdG4klo4tdU2DB6Ihmtn6t///DYsHwuXTtf/mS2cBDhCCCFaPaNBY87obgwM9amx7bBSAU75+jflWZap9sZfZucpPcAZeU3bMm1CSwKc+NomGlfF0QXG/RMmfwgu3nBhHywdAQdWN+y5LZQEOEIIIUQd9Ar0xN9dr0AcWUX+jYUl0fizXxPIKTTh28aJXoGeZdpYZnDi65poXJWev4fZOyBkKBRmwZcPwNqHoSC7cZ7fQkiAI4QQQtSBwaDxzj0D+NuEPtbTxqtimcFJKyn2N7ybPwZD2SWwWhf7qwuvYJi6DkbOB80AMR/qFZATDzTee9g4CXCEEEKIOhrSxY/7rgutMV+nR3sPHEoFNCPKLU/BlQCnwUtU5RkdYPTTeqDjEQRpJ+C9MbBraas45kECHCGEEKKJuDga6RF4ZVeWpf5NaSGlcnBUUwQenYbrS1Zht4CpENY/BR/fDTlpjf9eNkQCHCGEEKIJWZapwgI8CPB0qXA92EcPcLILirmU00Rbu9184a6P4OZFYHSC49/rxzyc+alp3s8GODR3ByozceJEtmzZwpgxY/j888+tr6enpxMVFUVxcTHFxcXMnTuXmTNnNmNPhRBCiOr9vl8gH++O554hIZVed3E00t7ThaTMfOIv5eLn7lyr5xabzLwVfYJvDiTy1E09uKlP1RWYAdA01OCZaCHXwecz9CWrVbfCiCdhxF/0Ja0qKKW4lFPIuct5nE/PIzu/GJNSmMwKc8k/TWZFR183xvauoR9XiaaaZD6sYTZv3kx2djarVq0qE+CYTCYKCgpwc3MjJyeHPn36sGfPHvz8qs9iB8jMzMTLy4uMjAw8PT1rbC+EEEI0Fsuv2qpydv7wn53sjrvEW3f1Z3z/DjU+70J6Ho99vJ89JYeDGg0a/ze5P7eFB1V5z88nU5n3WSzBPq68eHNnese8DDEf6BdDhsId7+rJyejB02d7zrHhcJIe1FzOI6/IVGO/Roe1ZcX0wTW2q4v6/v62yRmc0aNHs2XLlgqvG41G3Nz0qbyCggKUUk2zXimEEEI0opqSkUN93dgdd6lWW8V/PJzMk5/Hkp5bhIezAwNCfdh6/CKPf7KfomIzdwwMLtNeKcXKn8/w8rdHMJkVSZn53Pqfy9wf+TB/ue163Nb/GeJ/hiXDUOMX86P5Wl5bf5TTF8sWHtQ0vfJysI8bXq6OGDQNo0EPrvQ/axW2wDenOgc427ZtY9GiRezdu5fExETWrFnDhAkTyrR55513WLRoEUlJSYSHh/P2228zeHDjRHTp6emMHDmSEydOsGjRIvz9KyZsCSGEEC2J9dDNanZSFRSbeO37o6zYcQbQa+wsvnsAwT6uPLv2Nz7encCfP4+lyGTmrsH6clh+kYm/rj3I53vPATChfxAmBetiL7Dy5zN85+HPq6M+Z8yhp9Eu7EP79F6Sim/kfPEUfNu489CILvQO8iLYx5VAb5cqqzbbojoHODk5OYSHhzN9+nTuuOOOCtc//fRT5s2bx9KlSxkyZAhvvvkmY8eO5dixY7Rr1w6A/v37U1xcXOHeDRs2EBRU9fQagLe3N7GxsSQnJ3P77bczadIkAgIqHoxWUFBAQUGB9fvMzMy6DlUIIYS4KkL8qi/2l5lfxL3v/cKBkoM7Hxjemb/c1AMnB32v0CsT+uJoNPD+zrPM//I3Ck1mxvZuz0P/20tMQjoGDZ65pSd/HN4ZTdP4w6BgFnx1iLjUHB5YV0CY/zNMLF7JLIdvuN/hR27yjKPNPe/TJrjr1fkX0AQalIOjaVqFGZwhQ4Zw7bXXsnjxYgDMZjMdO3bk0UcfZf78+bV+9pYtW1i8eHGZHJzyZs+ezZgxY5g0aVKFay+88AIvvvhihdclB0cIIYSt2R9/mYn//pn2ni7semZMhetLt57ite+P4u3myD8nhRPVq+L/2CuleOXbI7z3UxwAXq6OZOQV4eXqyOJ7IsqcgQX67M6ybadZvPkkhcVmNA2e7n6BP6a+hjE3VT+h/Oa/Q8R9+vpUM6lvDk6jbhMvLCxk7969REVFXXkDg4GoqCh27tzZ4OcnJSWRlZUF6IHK9u3bCQsLq7Tt008/TUZGhvUrISGhwe8vhBBCNIVQvzaAfvJ4fiXJvF/FXADgL2N7VBrcgD7p8Oy4njw8Sp91ycgr4poAd75+ZFiF4Ab03VuPjenOj0+M4Imoa/jm0eE8OGMmxtk/Q5dRUJQLXz+q77jKz2icgV5FjZpknJqaislkqrBkFBAQwNGjR2v9nKioKGJjY8nJySE4OJjVq1cTGRlJfHw8Dz74oDW5+NFHH6Vv376VPsPZ2Rln59pttRNCCCGak4+bIx7ODmQVFJNwKZfuAVeKAx5PzuJIYiaORo1b+la/BVvTNJ4cG0awjxtn03J4dEx33J2r/1Uf6teGuVHdr7zgEQD3roGf34JNL8OhL+H8Xpi0HIIHNWicV9NV2UWllKoxg7y0jRs3Vvr64MGDiYmJaaxuCSGEEDZB0zQ6+rpxODGT+HIBzlcx5wEYeU07vN2cavWsqmru1JrBAMOfgE7Xw+fTIf0sLB8LN/wVhs7Vr9u4Ru2hv78/RqOR5OTkMq+npKRUmggshBBCCF1lh24qpazLU+P7V78Jp0kED4JZP0Hv28FcDBtfgA9uh6zkGm9tbo0a4Dg5OTFw4ECio6Otr5nNZqKjo4mMjGzMtxJCCCHsSukzqSz2xV/m3OU82jgZierZTBMFLl768tRtb4ODK5zeDEuHwcnKV1tsRZ0DnOzsbGJiYqxLRXFxccTExBAfHw/AvHnzWLZsGatWreLIkSPMnj2bnJwcpk+f3rg9F0IIIexISCWniltmb8b2bo+rUzPWoNE0GHA/PLQVAvpAzkX44A7Y8BwUN9H5WQ1U5xycPXv2MHr0aOv38+bNA2Dq1KmsXLmSyZMnc/HiRRYsWEBSUhL9+/dn/fr1skQlhBBCVCPUV99JdTZNryBcZDLzzYFEAMZH1Hx8w1XRNgweiIYNf4Vf34Wf/6Uf2DlpOfh2bu7elWGTZ1E1BTmLSgghhC2LT8tlxKLNODkYOPrSTWw9cZHpK37F392JXU+PwcFoY4m9R9bBV49Afjo4ecCtb0LfinXpGsom6uAIIYQQon6CvF1wMGgUFptJzsrnq/367qnf9wuyveAGoOetegJySCQUZsEXf4S1c8BU8aSC5mCD/8aEEEKI1sfBaKCDjysAR5Oy2HBY36l0W3Psnqot744w9RsY+RRoBjAXgdE2zvG2jV4IIYQQghBfN86m5bL8pzhyC02E+LoR0dG7ubtVPaMDjH4GuoyG9n2auzdWEuAIIYQQNsKyVXz7iVRAr31Tl0K5zSrUtsrByBKVEEIIYSMsAY5FsxT3sxMS4AghhBA2wlLNGKB3kCfd2nlU01pURwIcIYQQwkaElNTCAZm9aSgJcIQQQggbEernhrODAaNB49ZwCXAaQpKMhRBCCBvRxtmBFdOuRQGBXq7N3Z0WTQIcIYQQwoYM7ebf3F2wC7JEJYQQQgi7IwGOEEIIIeyOBDhCCCGEsDsS4AghhBDC7kiAI4QQQgi7IwGOEEIIIeyOBDhCCCGEsDsS4AghhBDC7kiAI4QQQgi7IwGOEEIIIeyOBDhCCCGEsDsS4AghhBDC7kiAI4QQQgi702pOE1dKAZCZmdnMPRFCCCFEbVl+b1t+j9dWqwlwsrKyAOjYsWMz90QIIYQQdZWVlYWXl1et22uqriFRC2U2m7lw4QIeHh5omlavZ2RmZtKxY0cSEhLw9PRs5B7aBhmjfZAx2ofWMEZoHeOUMdafUoqsrCyCgoIwGGqfWdNqZnAMBgPBwcGN8ixPT0+7/QtqIWO0DzJG+9AaxgitY5wyxvqpy8yNhSQZCyGEEMLuSIAjhBBCCLtjfOGFF15o7k60JEajkVGjRuHgYL+rezJG+yBjtA+tYYzQOsYpY7y6Wk2SsRBCCCFaD1miEkIIIYTdkQBHCCGEEHZHAhwhhBBC2B0JcIQQQghhdyTAEUIIIYTdkQCnDt555x06deqEi4sLQ4YMYffu3c3dpXrbtm0bt956K0FBQWiaxtq1a8tcV0qxYMECAgMDcXV1JSoqihMnTjRTb+tu4cKFXHvttXh4eNCuXTsmTJjAsWPHyrTJz89nzpw5+Pn54e7uzh133EFycnIz9bh+lixZQr9+/ayVQyMjI/n++++t1+1hjKUtXLgQTdN4/PHHra/ZwxhfeOEFNE0r89WjRw/rdXsYI8D58+e599578fPzw9XVlb59+7Jnzx7r9Zb+c6dTp04VPkdN05gzZw5gH5+jyWTiueeeo3Pnzri6utK1a1f+9re/lTkI02Y+RyVq5ZNPPlFOTk5q+fLl6tChQ2rmzJnK29tbJScnN3fX6uW7775Tzz77rPriiy8UoNasWVPm+muvvaa8vLzUmjVrVGxsrLrttttU586dVV5eXjP1uG7Gjh2rVqxYoQ4ePKhiYmLULbfcokJCQlR2dra1zaxZs1THjh1VdHS02rNnj7ruuuvU0KFDm7HXdff111+rb7/9Vh07dkwdO3ZMPfPMM8rR0VEdPHhQKWUfY7TYvXu36tSpk+rXr5+aO3eu9XV7GOPzzz+vevfurRITE61fFy9etF63hzFeunRJhYaGqmnTpqlffvlFnT59Wv3www/q5MmT1jYt/edOSkpKmc/wxx9/VIDavHmzUso+PsdXXnlF+fn5qW+++UbFxcWp1atXK3d3d/XWW29Z29jK5ygBTi0NHjxYzZkzx/q9yWRSQUFBauHChc3Yq8ZRPsAxm82qffv2atGiRdbX0tPTlbOzs/r444+bo4sNlpKSogC1detWpZQ+HkdHR7V69WprmyNHjihA7dy5s7m62Sh8fHzUe++9Z1djzMrKUt27d1c//vijGjlypDXAsZcxPv/88yo8PLzSa/YyxqeeekoNHz68yuv2+HNn7ty5qmvXrspsNtvN5zhu3Dg1Y8aMMq/dfvvtasqUKUop2/ocZYmqFgoLC9m7dy9RUVHW1wwGA1FRUezcubMZe9Y04uLiSEpKKjNeLy8vhgwZ0mLHm5GRAYCvry8Ae/fupaioqMwYe/ToQUhISIsdo8lk4pNPPiEnJ4fIyEi7GuOcOXMYN25cmbGAfX2OJ06cICgoiC5dujBlyhTi4+MB+xnj119/zaBBg7jzzjtp164dERERvPvuu9br9vZzp7CwkA8++IAZM2agaZrdfI5Dhw4lOjqa48ePAxAbG8tPP/3EzTffDNjW59j8tZRbgNTUVEwmEwEBAWVeDwgI4OjRo83Uq6aTlJQEUOl4LddaErPZzOOPP86wYcPo06cPoI/RyckJb2/vMm1b4hh/++03IiMjyc/Px93dnTVr1tCrVy9iYmLsYoyffPIJ+/bt49dff61wzV4+xyFDhrBy5UrCwsJITEzkxRdf5Prrr+fgwYN2M8bTp0+zZMkS5s2bxzPPPMPu3bt57LHHcHZ25v7777e7nztr164lPT2dadOmAfbzd3X+/PlkZmbSo0cPjEYjJpOJV155hSlTpgC29ftDApwGUEqhaVpzd+OqUUphMLS8Sb85c+Zw8OBBfvrppxrbtsTPNCwsjJiYGNLT0/niiy+YOnUqW7durbJ9SxpjQkICc+fOZcOGDbi4uNT6vpY0RsD6f78A/fr1Y8iQIYSGhvLZZ5/h6upa6T0tbYxms5lBgwbx6quvAhAREcGhQ4dYsmQJ999/f5X3tdSfO//973+5+eabCQoKqrZdS/scP/vsMz788EM++ugjevfuTUxMDI8//jhBQUFMnTq1yvua43NseX9rmoG/vz9Go7FCtntKSkqFKNUetG/fHsAuxvvII4/wzTffsHnzZoKDg62vt2/fnsLCQtLT08u0b4ljdHJyolu3bgwaNIiFCxcSHh7OW2+9ZRdj3Lt3LykpKQwcOBAHBwccHBzYunUr//rXv3BwcCAgIKDFj7Ey3t7eXHPNNZw8edIuPkeAwMBAevXqVea1nj17Wpfi7OnnztmzZ9m4cSMPPPCA9TV7+RyffPJJ5s+fz1133UXfvn257777eOKJJ1i4cCFgW5+jBDi14OTkxMCBA4mOjra+ZjabiY6OJjIyshl71jQ6d+5M+/bty4w3MzOTX375pcWMVynFI488wpo1a9i0aROdO3cuc33gwIE4OjqWGePx48eJj49vMWOsitlspqCgwC7GOGbMGH777TdiYmKsX4MGDWLKlCnWP7f0MVYmOzubU6dOERgYaBefI8CwYcMqlGo4fvw4oaGhgH383LFYsWIF7dq1Y9y4cdbX7OVzzM3NrTATYzQaMZvNgI19jlc1pbkFs2wTX7lypTp8+LB68MEHlbe3t0pKSmrurtVLVlaW2r9/v9q/f78C1BtvvKH279+vzp49q5TSt/l5e3urr776Sh04cECNHz++RW3XnD17tvLy8lJbtmwps20zNzfX2mbWrFkqJCREbdq0Se3Zs0dFRkaqyMjIZux13T399NNq27ZtKi4uTh04cEDNnz9faZqmNmzYoJSyjzGWV3oXlVL2McY//elPasuWLSouLk7t2LFDRUVFKX9/f5WSkqKUso8x7t69Wzk4OKhXXnlFnThxQn344YfKzc1NffDBB9Y2Lf3njlL6DtuQkBD11FNPVbhmD5/j1KlTVYcOHazbxL/88kvl7++v/vKXv1jb2MrnKAFOHbz99tsqJCREOTk5qcGDB6tdu3Y1d5fqbfPmzQqo8DV16lSllL7V77nnnlMBAQHK2dlZjRkzRh07dqx5O10HlY0NUCtWrLC2ycvLUw8//LDy8fFRbm5uauLEiSoxMbH5Ol0PM2bMUKGhocrJyUm1bdtWjRkzxhrcKGUfYyyvfIBjD2OcPHmyCgwMVE5OTqpDhw5q8uTJZerD2MMYlVJq3bp1qk+fPsrZ2Vn16NFDLVu2rMz1lv5zRymlfvjhBwVU2m97+BwzMzPV3LlzVUhIiHJxcVFdunRRzz77rCooKLC2sZXPUVOqVPlBIYQQQgg7IDk4QgghhLA7EuAIIYQQwu5IgCOEEEIIuyMBjhBCCCHsjgQ4QgghhLA7EuAIIYQQwu5IgCOEEEIIuyMBjhBCCCHsjgQ4QgghhLA7EuAIIYQQwu5IgCOEEEIIu/P/Wu/gHQ2o6TgAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7fda51524710>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x7fda4f4f3690>\n",
       " PyObject <matplotlib.lines.Line2D object at 0x7fda4f4cf4d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = abs(corrs[1])\n",
    "semilogy(1:m, abs.(corrs), [1,m], [C, C*exp(-m/ξ)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# But I don't care about random MPSes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you know how to compute correlation functions for MPS states, we should start thinking about how do we get a useful MPS state in the first place. So far all our states have been defined by tensors $A$ with random elements.\n",
    "\n",
    "Some fun little analytical examples exist, such as how to represent GHZ or W states as MPSes. You can google them up. The real question is though, given an arbitrary local Hamiltonian, can we find an MPS for the ground state. We know that for gapped Hamiltonians (i.e. most Hamiltonians) correlations in the ground state decay exponentially and entanglement is limited (area law), so there's hope that the answer would be \"yes we can\".\n",
    "\n",
    "You can think about this problem in a few different ways. One is to ask, what's the lowest energy state we can find on the manifold of MPS states. The question is then a matter of optimizing for minimal expectation value $\\langle \\psi | H | \\psi \\rangle$ under the constraint that the state $|\\psi \\rangle$ is an MPS. Another way to approach this would be writing down some expression for the exact ground state (an expression we obviously can't just go an evaluate, otherwise there would be no problem), and then try to massage that expression into an MPS.\n",
    "\n",
    "Various algorithms for implementing these different approaches in different situations (finite-size, infinite-size, translation invariant or not, what kind of boundary conditions) exist. Some examples are the Density Matrix Renormalization Group (DMRG), Time-Dependent Variational Principle (TDVP) and Time-Evolving Block Decimation (TEBD). We'll implement an infinite-lattice version of TEBD. Why TEBD? Because it's simple, and because it you can also use it to do time-evolution of MPS states. (I maintain that the choice to pick TEBD is independent of the fact that Guifre invented it.)\n",
    "\n",
    "Implementing iTEBD will be the rest of this minicourse. However, we'll only really get to that next time. First, we need to talk a bit about gauge freedom in MPS and how to fix this freedom, since TEBD will be making use of a specific gauge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauge freedom and canonical forms\n",
    "In any tensor network, we can always take any contracted leg, and insert on that leg a pair of invertible matrices $g g^{-1} = \\mathbb{1}$. This doesn't affect the value the network contracts to in anyway: we've just changed the basis in which we perform one of the sums (traces).\n",
    "\n",
    "This means that an MPS representation of a state is not unique. We can insert any $g g^{-1}$ on any of the virtual legs of the MPS, and the state stays the same. Obviously, we would like to stay within the framework of translation invariant MPSes, so we allow only for gauge transformations (as these changes of basis are called) of the form\n",
    "```\n",
    "--A-- -> --gi--A--g--\n",
    "  |            |\n",
    "```\n",
    "where by `gi` I mean the inverse of $g$.\n",
    "\n",
    "Since we have some freedom in picking $A$, we should probably try to make use of it somehow, while fixing the freedom. How? Well, we could try to impose some nice properties on $A$, if such properties can be achieved by a gauge transformation. And it turns out there's a really useful property that we can impose, which goes as follows.\n",
    "\n",
    "First, we have to complicate our setup slightly. Previously our how MPS was defined by $A$. Now, let's think of an MPS that is instead written as\n",
    "```\n",
    "... --λ--Γ--λ--Γ--λ--Γ--λ-- ...\n",
    "         |     |     |\n",
    "```\n",
    "Here $\\Gamma$ is some tensor, and $\\lambda$ is a diagonal matrix, with real, non-negative elements. Obviously we could just absorb $\\lambda$ into $\\Gamma$ and a get single tensor $A$. The point is, by using the gauge freedom, we can choose $\\Gamma$ and $\\lambda$ in such a way, that *$\\lambda$ always has on its diagonal the Schmidt coefficients of partitioning the state into two parts at $\\lambda$*. What are Schmidt coefficients? https://en.wikipedia.org/wiki/Schmidt_decomposition\n",
    "\n",
    "Another way of stating the same condition is to say that the left dominant eigenvector of\n",
    "```\n",
    "--λ--Γ--\n",
    "     |\n",
    "--λ--Γ--\n",
    "```\n",
    "and the right dominant eigenvector of\n",
    "```\n",
    "--Γ--λ--\n",
    "  |\n",
    "--Γ--λ--\n",
    "```\n",
    "are both the identity matrix.\n",
    "\n",
    "Choosing our gauge so that this property holds turns out to be quite handy. You'll see next time how it's used in the TEBD algorithm, but you can probably already imagine that it's kinda nice that on each virtual legs you have these Schmidt coefficients $\\lambda$, that immediately tell you everything you would want to know about the entanglement between the two sides of the state separated by this leg. The bond dimension is also obviously the Schmidt rank, which gives us a very concrete, physical characterization of what high bond dimension and low bond dimension mean, and in what sense are MPSes \"low entanglement\" states.\n",
    "\n",
    "So, given $A$, how do we find a gauge transformation such that the above property holds? I'm getting tired of typing, so just go read Section II of https://arxiv.org/pdf/0711.3960.pdf, that will explain it all. Then, it's time for your homework."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 3\n",
    "#### -1)\n",
    "Make sure that, for the expectation values and correlation functions that you implemented last time,<br>\n",
    "a) Your expectation values for Hermitian operators are real.<br>\n",
    "b) Your correlation functions decay like they should (make a plot).<br>\n",
    "c) You can generate two-point correlators with $D=100$ for distances $1, \\dots, 100$ or so in a few seconds.<br>\n",
    "Feel free to make use of my code above, as long as you understand what's going on. You don't necessarily need to send this part of the homework to me once your done, just see yourself that it works.\n",
    "\n",
    "#### 0)\n",
    "Read Section II of https://arxiv.org/pdf/0711.3960.pdf. Note that there they start from some $\\Gamma$ and $\\lambda$, that are not the canonical form, and try to turn them into $\\Gamma'$ and $\\lambda'$ that would be in the canonical form. We, on the other hand, start from just $A$. You can easily translate to their language by setting their $\\Gamma$ to your $A$, and their $\\lambda$ to just the identity matrix.\n",
    "\n",
    "#### 1)\n",
    "As explained in above reference, the form of MPS transfer matrix $T$, i.e. the symmetry it has because of how it consists of $A$ and $A^\\ast$, guarantees that the dominant eigenvectors $l$ and $r$ are Hermitian and positive semi-definite. However, as we discussed in the last meeting, when we numerically ask for the eigenvectors, they come with an arbitrary phase. Last time we modified `normalize!` so that it partially fixes this phase (and also the norm of $l$ and $r$), imposing `vec(l)'*vec(r) == 1`. Your first job is to further modify `normalize!` so that it fully fixes the phases, so that the `l` and `r` it returns are Hermitian and positive semi-definite, in addition to having the property that `vec(l)'*vec(r) == 1`.\n",
    "\n",
    "How do you do that? Well, you should definitely start by getting the `normalize!` function from above, and only make small modifications to that. Since we know that $r$ is actually $a \\cdot \\tilde{r}$, where $\\tilde{r}$ is Hermitian and pos. semi-def. and $a$ is a scalar, we can just take the trace $\\mathrm{Tr} r = a \\mathrm{Tr} \\tilde{r}$, and we know that $\\mathrm{Tr} \\tilde{r}$ is something real and positive. Any part of $\\mathrm{Tr} r$ that is not real and positive, is $a$, and you divide $r$ by $a$ to get $\\tilde{r}$. And similarly for $l$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'normalize! :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "normalize!"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    normalize!(A)\n",
    "\n",
    "Normalize the UMPS defined by A, and return the dominant left and right\n",
    "eigenvectors l and r of its transfer matrix, normalized so that they are\n",
    "both Hermitian and positive semi-definite (when thought of as matrices),\n",
    "and l'*r = 1.\n",
    "\"\"\"\n",
    "function normalize!(A)\n",
    "    SR, UR, SL, UL = tm_eigs(A, \"BOTH\", 1)\n",
    "    S1 = SR[1]\n",
    "    A ./= sqrt(S1)\n",
    "    \n",
    "    l = UL[1]\n",
    "    r = UR[1]  \n",
    "    # ???\n",
    "    return l, r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2)\n",
    "Now that we can trust our $l$ and $r$ to be Hermitian and $\\geq 0$, we can go ahead and implement a function that takes in $A$, $l$ and $r$, and outputs $\\Gamma$ and $\\lambda$, that define the same MPS, but now in canonical form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "canonical_form (generic function with 1 method)"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    canonical_form(A, l, r)\n",
    "\n",
    "Return a three-valent tensor Γ and a vector λ, that define the canonical\n",
    "of the UMPS defined by A. l and r should be the normalized dominant\n",
    "left and right eigenvectors of A.\n",
    "\"\"\"\n",
    "function canonical_form(A, l, r)\n",
    "    # ???\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3)\n",
    "Prove to yourself that<br>\n",
    "a)<br>\n",
    "the way you've chosen to transform the MPS, going from $A$ to $(\\Gamma, \\lambda)$, is actually a gauge transformation. In other words, prove that they both define the same state.\n",
    "\n",
    "b)<br>\n",
    "that the MPS defined by your new $\\Gamma$ and $\\lambda$ is actually in canonical form, i.e., that the left dominant eigenvector of\n",
    "```\n",
    "--λ--Γ--\n",
    "     |\n",
    "--λ--Γ--\n",
    "```\n",
    "and the right dominant eigenvector of\n",
    "```\n",
    "--Γ--λ--\n",
    "  |\n",
    "--Γ--λ--\n",
    "```\n",
    "are both the identity matrix. Do this both numerically and analytically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}