{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "30a3bbb2",
   "metadata": {},
   "source": [
    "# Variational circuit preparation I\n",
    "\n",
    "In this tutorial, we demonstrate two complementary classical pre-processing techniques to prepare a quantum circuit that approximates the ground state of a given Hamiltonian. Such a circuit can be used as an initial state in quantum phase estimation (QPE) or other quantum algorithms.\n",
    "\n",
    "We consider the 1D transverse-field Ising model and build a brick-wall circuit of general two-qubit gates, where each of them can be parametrized by at most 3 CNOT gates and 15 elementary one-qubit gates (for more details, read Vatan and Williams, [Phys. Rev. A 69, 032315 (2004)](https://doi.org/10.1103/PhysRevA.69.032315)).\n",
    "\n",
    "The two approaches are:\n",
    "\n",
    "1. **Global optimization** - all parameters of the circuit are optimized simultaneously using gradient-based methods (L-BFGS) via auto-differentiation. The parameters are the 15 real entries of the SU(4) gates, and the unitarity is enforced by construction.\n",
    "2. **Local optimization** - the circuit is represented as a tensor network, and the individual tensor entries are optimized locally. The optimization is performed by sweeping through the tensors, similar to DMRG, and often yields high-fidelity states.\n",
    "\n",
    "We compare the performance of both strategies and also illustrate a **sequential depth optimisation** that starts from a shallow circuit and progressively increases the depth, reusing parameters from shallower layers.\n",
    "\n",
    "The target state is obtained by DMRG (via `quimb.tensor.DMRG2`), which provides a reference energy and the numerically exact ground state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f5e0b8a2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T05:58:49.078406Z",
     "iopub.status.busy": "2026-09-15T05:58:49.078241Z",
     "iopub.status.idle": "2026-09-15T05:58:50.631839Z",
     "shell.execute_reply": "2026-09-15T05:58:50.631001Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"NUMBA_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"MKL_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"JAX_ENABLE_X64\"] = \"True\"\n",
    "\n",
    "import autoray\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import quimb.tensor as qtn\n",
    "\n",
    "# Local imports from qpe_toolbox\n",
    "from qpe_toolbox.circuit import ansatz_circuit_su4, su4swap_gate_param_gen, tn_fit\n",
    "from qpe_toolbox.hamiltonian import Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "018cdb05",
   "metadata": {},
   "source": [
    "## Hamiltonian: 1D Transverse-Field Ising (TFI) model\n",
    "\n",
    "We consider a chain of $n$ spins with open boundaries. The Hamiltonian reads\n",
    "\n",
    "$$ H = g_x \\sum_{i} X_i + g_{zz} \\sum_{i} Z_i Z_{i+1}, $$\n",
    "\n",
    "with $g_x = -1.1$ and $g_{zz} = -1.0$. We take $n = 8$ sites.\n",
    "\n",
    "First, we compute the ground state using DMRG (via `qtn.DMRG2`). This gives us the target state `GS` and the exact energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "820ab0b9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T05:58:50.633702Z",
     "iopub.status.busy": "2026-09-15T05:58:50.633339Z",
     "iopub.status.idle": "2026-09-15T05:58:50.712419Z",
     "shell.execute_reply": "2026-09-15T05:58:50.711716Z"
    },
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** DMRG reference energy: -10.4564576523\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# --- Hamiltonian definition ---\n",
    "n_qubits = 8\n",
    "gx, gzz = -1.1, -1.0\n",
    "terms = []\n",
    "for x in range(n_qubits):\n",
    "    terms.append((gx, \"x\", [x]))  # transverse field\n",
    "for x in range(n_qubits - 1):\n",
    "    terms.append((gzz, \"zz\", [x, x + 1]))  # nearest-neighbour coupling\n",
    "ham = Hamiltonian(terms, n_qubits)\n",
    "mpo = ham.to_mpo()\n",
    "\n",
    "# --- DMRG reference ground state ---\n",
    "dmrg = qtn.DMRG2(mpo)\n",
    "dmrg.solve(max_sweeps=16, tol=1e-8, bond_dims=64, verbosity=0)\n",
    "GS = dmrg.state\n",
    "dmrg_energy = np.real(dmrg.energy)\n",
    "print(f\"*** DMRG reference energy: {dmrg_energy:12.10f}\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88174f8c",
   "metadata": {},
   "source": [
    "## Global Optimization\n",
    "\n",
    "The circuit is built from SU(4) gates, each parameterized by 15 real numbers. We define a loss function that computes the expectation value of the Hamiltonian MPO with respect to the state produced by the circuit. The gradient is obtained via automatic differentiation (JAX).\n",
    "Ref: Haghshenas et al. [Phys. Rev. X 12, 011047 (2022)](https://doi.org/10.1103/PhysRevX.12.011047), [Tensor Network Training of Quantum Circuits](https://quimb.readthedocs.io/en/latest/examples/ex_tn_train_circuit.html)\n",
    "\n",
    "Three global optimization strategies are compared:\n",
    "1. **Standard L-BFGS** on a fixed-depth circuit.\n",
    "2. **Basin-hopping** - a global optimisation method that helps escape local minima.\n",
    "3. **Sequential layer-wise optimization**, where a shallower circuit is optimized first, and its parameters are used to initialize a deeper one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32f2e239",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T05:58:50.714088Z",
     "iopub.status.busy": "2026-09-15T05:58:50.713933Z",
     "iopub.status.idle": "2026-09-15T05:58:50.717123Z",
     "shell.execute_reply": "2026-09-15T05:58:50.716682Z"
    },
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def loss_circ(circ, mpo):\n",
    "    \"\"\"\n",
    "    Loss function: expectation value of the MPO Hamiltonian with respect to\n",
    "    the state produced by the circuit.\n",
    "    \"\"\"\n",
    "    psi = circ.psi\n",
    "    psiH = psi.H\n",
    "    norm_tn = psiH & psi\n",
    "    psi.align_(mpo, psiH)\n",
    "    energy_tn = psiH & mpo & psi\n",
    "    energy = autoray.do(\"real\", energy_tn.contract())\n",
    "    norm = autoray.do(\"real\", norm_tn.contract())\n",
    "    return (energy / norm).real\n",
    "\n",
    "\n",
    "def make_circuit_optimizer(circ, mpo):\n",
    "    \"\"\"\n",
    "    Factory to create a TNOptimizer for a circuit with the given MPO.\n",
    "\n",
    "    Uses L-BFGS-B (or basin-hopping) and JAX for gradients.\n",
    "    \"\"\"\n",
    "    return qtn.TNOptimizer(\n",
    "        circ,\n",
    "        loss_circ,\n",
    "        loss_constants={\"mpo\": mpo},\n",
    "        autodiff_backend=\"jax\",\n",
    "        optimizer=\"L-BFGS-B\",\n",
    "        progbar=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b4a935a",
   "metadata": {},
   "source": [
    "### 1. Direct L-BFGS optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "436c1ad0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T05:58:50.718814Z",
     "iopub.status.busy": "2026-09-15T05:58:50.718673Z",
     "iopub.status.idle": "2026-09-15T05:59:24.620248Z",
     "shell.execute_reply": "2026-09-15T05:59:24.619455Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Global L-BFGS optimization\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  630  Energy = -10.45638540  Error =  6.910e-06  1-F =  1.332e-05\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"*** Global L-BFGS optimization\")\n",
    "depth = 6\n",
    "circ = ansatz_circuit_su4(n_qubits, depth)\n",
    "circ_opt = make_circuit_optimizer(circ, mpo)\n",
    "optimal_circ = circ_opt.optimize(n=10000, tol=1e-8)\n",
    "ovlp = (dmrg.state.H & optimal_circ.psi).contract()\n",
    "print(\n",
    "    f\" # parameters = {circ_opt.d: 4d}\",\n",
    "    f\" Energy = {circ_opt.loss: >12.8f}\",\n",
    "    f\" Error = {np.abs(1 - circ_opt.loss / dmrg_energy): >10.3e}\",\n",
    "    f\" 1-F = {1 - np.abs(ovlp) ** 2: >10.3e}\",\n",
    ")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc4e6ed6",
   "metadata": {},
   "source": [
    "### 2. Basin-hopping optimization\n",
    "\n",
    "Basin-hopping performs random steps in parameter space followed by a local minimization. Here we use 5000 iterations with 10 hops per step and a temperature of 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "81c10074",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T05:59:24.621969Z",
     "iopub.status.busy": "2026-09-15T05:59:24.621595Z",
     "iopub.status.idle": "2026-09-15T06:00:21.596470Z",
     "shell.execute_reply": "2026-09-15T06:00:21.595849Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Global L-BFGS optimization + Basin-hopping \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  630  Energy = -10.45638978  Error =  6.491e-06  1-F =  4.192e-06\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"*** Global L-BFGS optimization + Basin-hopping \")\n",
    "depth = 6\n",
    "circ = ansatz_circuit_su4(n_qubits, depth)\n",
    "circ_opt = make_circuit_optimizer(circ, mpo)\n",
    "optimal_circ = circ_opt.optimize_basinhopping(n=5000, nhop=10, temperature=0.1)\n",
    "ovlp = (dmrg.state.H & optimal_circ.psi).contract()\n",
    "print(\n",
    "    f\" # parameters = {circ_opt.d: 4d}\",\n",
    "    f\" Energy = {circ_opt.loss: >12.8f}\",\n",
    "    f\" Error = {np.abs(1 - circ_opt.loss / dmrg_energy): >10.3e}\",\n",
    "    f\" 1-F = {1 - np.abs(ovlp) ** 2: >10.3e}\",\n",
    ")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bff974d",
   "metadata": {},
   "source": [
    "### 3. Sequential layer-wise optimization\n",
    "\n",
    "We start from depth 1 and increase the depth one layer at a time, always reusing the parameters from the shallower circuit as initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8348e481",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T06:00:21.598216Z",
     "iopub.status.busy": "2026-09-15T06:00:21.598070Z",
     "iopub.status.idle": "2026-09-15T06:02:12.732520Z",
     "shell.execute_reply": "2026-09-15T06:02:12.731891Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Global L-BFGS sequential optimization \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  105  Energy = -10.26060113  Error =  1.873e-02  1-F =  2.611e-01\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  210  Energy = -10.44128128  Error =  1.451e-03  1-F =  6.541e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  315  Energy = -10.45641137  Error =  4.427e-06  1-F =  9.179e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  420  Energy = -10.45642036  Error =  3.566e-06  1-F =  7.253e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  525  Energy = -10.45642342  Error =  3.273e-06  1-F =  6.531e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # parameters =  630  Energy = -10.45642066  Error =  3.538e-06  1-F =  6.245e-06\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"*** Global L-BFGS sequential optimization \")\n",
    "depths_global = []\n",
    "errors_global = []\n",
    "rng = np.random.default_rng(42)\n",
    "\n",
    "circ = ansatz_circuit_su4(n_qubits, 1, rng=rng)\n",
    "circ_opt = make_circuit_optimizer(circ, mpo)\n",
    "optimal_circ = circ_opt.optimize(n=10000, tol=1e-8)\n",
    "ovlp = (dmrg.state.H & optimal_circ.psi).contract()\n",
    "err = np.abs(1 - circ_opt.loss / dmrg_energy)\n",
    "depths_global.append(1)\n",
    "errors_global.append(err)\n",
    "print(\n",
    "    f\" # parameters = {circ_opt.d: 4d}\",\n",
    "    f\" Energy = {circ_opt.loss: >12.8f}\",\n",
    "    f\" Error = {err: >10.3e}\",\n",
    "    f\" 1-F = {1 - np.abs(ovlp) ** 2: >10.3e}\",\n",
    ")\n",
    "\n",
    "for ii in range(2, depth + 1):\n",
    "    # grow the optimized circuit by one SU4SWAP layer initialized close to identity\n",
    "    for start in range(2):\n",
    "        for q in range(start, n_qubits - 1, 2):\n",
    "            optimal_circ.apply_gate(\n",
    "                \"SU4SWAP\",\n",
    "                *(0.1 * rng.random(15)),\n",
    "                q,\n",
    "                q + 1,\n",
    "                gate_round=ii - 1,\n",
    "                parametrize=True,\n",
    "            )\n",
    "    circ_opt = make_circuit_optimizer(optimal_circ, mpo)\n",
    "    optimal_circ = circ_opt.optimize(n=10000, tol=1e-8)\n",
    "    ovlp = (dmrg.state.H & optimal_circ.psi).contract()\n",
    "    err = np.abs(1 - circ_opt.loss / dmrg_energy)\n",
    "    depths_global.append(ii)\n",
    "    errors_global.append(err)\n",
    "    print(\n",
    "        f\" # parameters = {circ_opt.d: 4d}\",\n",
    "        f\" Energy = {circ_opt.loss: >12.8f}\",\n",
    "        f\" Error = {err: >10.3e}\",\n",
    "        f\" 1-F = {1 - np.abs(ovlp) ** 2: >10.3e}\",\n",
    "    )\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79db0fd8",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "#### Plot: Error vs depth for global sequential optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "526ec0d1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T06:02:12.734197Z",
     "iopub.status.busy": "2026-09-15T06:02:12.734026Z",
     "iopub.status.idle": "2026-09-15T06:02:12.935337Z",
     "shell.execute_reply": "2026-09-15T06:02:12.934636Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAGGCAYAAABmPbWyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYDJJREFUeJzt3XdYU9f/B/D3TSAJW1GmIqC4EBGLirNura3W0WFtVZxtEa0Wbavtr44Ote33ayfO1tHWfrXO2qF17w1ii7jBCQKKbFnJ/f1BSY0MQ0hyIbxfz5PnSe4999xPDon5eM655wqiKIogIiIiogrJpA6AiIiIqCZg0kRERESkByZNRERERHpg0kRERESkByZNRERERHpg0kRERESkByZNRERERHpg0kRERESkByZNRERERHpg0kQWbe7cuRAEwaBjfXx8MHDgQKPFcu3aNQiCgNWrVxutTrJs+/fvhyAI2L9/f6WPrcpn31RMEdOYMWPg4+Nj1Dqr83lJWkyaqMZJSEjA5MmT0axZM9ja2sLW1hb+/v4IDw/HX3/9JXV4VVbyQ7lx48ZKHysIgs7Dzs4O/v7++Oijj5Cbm6tTdsyYMaXKlzx27NihUzYzMxMff/wx2rVrBycnJyiVSnh7e2P48OH4/fffS8Vx7do1jB07Fk2aNIFKpYK7uzuefPJJzJkzp9LvqTZYvHixxSTTubm5mDt3rkGJXnWTmJiIuXPnIiYmRupQqJqwkjoAosr47bffMHz4cFhZWeGVV15BmzZtIJPJcOHCBWzevBlLlixBQkICvL29pQ5VMn379sXo0aMBANnZ2Th06BDef/99nD17Fhs2bNApq1Qq8e2335aqo02bNtrnV65cQf/+/XH9+nUMHToUo0ePhr29PW7evIk//vgDAwcOxPfff49Ro0Zpy7dv3x42NjYYN24cfHx8kJSUhOjoaHzyySeYN2+eCd99zbR48WLUr18fY8aM0dn+5JNP4sGDB1AoFNIEZoDc3Fzt37hHjx46+/7v//4PM2fONOr5VqxYAY1GY9Q6SyQmJmLevHnw8fFBUFCQ2c5L1ReTJqoxrl69ipdeegne3t7Ys2cPPDw8dPZ/8sknWLx4MWSy2t2B2qxZM4wcOVL7+vXXX0dBQQE2b96MvLw8qFQq7T4rKyudso8qKirC0KFDkZycjAMHDqBLly46++fMmYOdO3dCrVZrt33++efIzs5GTExMqeQ1JSWlqm+vVpHJZDp/r5rOysoKVlbG/dmxtrY2an3V/bwkrdr960I1yqeffoqcnBysWrWqVMIEFP+D/MYbb8DLy6vCeoqKivDhhx+iSZMmUCqV8PHxwbvvvov8/Pwyy+/cuRNBQUFQqVTw9/fH5s2bdfanpaVhxowZaN26Nezt7eHo6IgBAwbg7Nmzhr9ZI3N3d4cgCJX+wdqwYQNiY2Px/vvvl0qYSvTr1w8DBgzQvr569SoaNmxYZm+fq6vrY895584djB07Fg0bNoRSqYSHhwcGDx6Ma9eu6ZTbvn07unXrBjs7Ozg4OOCZZ57BuXPnStW3detWBAQEQKVSISAgAFu2bCk1H6W8uUPlzUO7cOECnn/+eTg7O0OlUqFdu3bYtm2bTpnVq1dDEAQcOXIEERERcHFxgZ2dHYYOHYrU1FRtOR8fH5w7dw4HDhzQDo+W9NCUFdehQ4fwwgsvoFGjRlAqlfDy8sKbb76JBw8ePLZty7NhwwYEBwfDxsYG9evXx8iRI3H79m2dMmPGjIG9vT3i4+PRv39/2NnZwdPTEx988AFEUdS2l4uLCwBg3rx52vczd+5cAGXPaRIEAZMnT8aGDRvg7+8PGxsbdOrUCX///TcAYNmyZfDz84NKpUKPHj1KfQ7Kmlu0bt06BAcHw8HBAY6OjmjdujW+/PJL7X59vrP79+9H+/btAQBjx47VvpeSz0JZ583JycH06dPh5eUFpVKJ5s2b4z//+Y+2fR59zyWfTaVSiVatWpUaFqfqhz1NVGP89ttv8PPzQ0hISJXqmTBhAtasWYPnn38e06dPx4kTJ7BgwQKcP38eW7Zs0Sl7+fJlDB8+HK+//jpCQ0OxatUqvPDCC9ixYwf69u0LAIiPj8fWrVvxwgsvwNfXF8nJyVi2bBm6d++OuLg4eHp6VineysrLy8Pdu3cBFP8jfuTIEaxZswYvv/xymUlTSdkS1tbWcHJyAgD8+uuvAFBhb9SjvL29sXv3buzduxe9evWqdPzPPfcczp07hylTpsDHxwcpKSnYtWsXbty4of2R+uGHHxAaGor+/fvjk08+QW5uLpYsWYKuXbvizJkz2nI7d+7Ec889B39/fyxYsAD37t3TJmSGOnfuHLp06YIGDRpg5syZsLOzw88//4whQ4Zg06ZNGDp0qE75KVOmoG7dupgzZw6uXbuGL774ApMnT8b69esBAF988QWmTJkCe3t7vPfeewAANze3cs+/YcMG5ObmIiwsDPXq1cPJkyfx9ddf49atW6WGX/WxevVqjB07Fu3bt8eCBQuQnJyML7/8EkeOHMGZM2dQp04dbVm1Wo2nnnoKHTt2xKeffoodO3Zgzpw5KCoqwgcffAAXFxcsWbIEYWFhGDp0KIYNGwYACAwMrDCGQ4cOYdu2bQgPDwcALFiwAAMHDsTbb7+NxYsXY9KkSbh//z4+/fRTjBs3Dnv37i23rl27dmHEiBHo3bs3PvnkEwDA+fPnceTIEUydOhWAft/Zli1b4oMPPsDs2bPx6quvolu3bgCAzp07l3leURTx7LPPYt++fRg/fjyCgoLw559/4q233sLt27fx+eef65Q/fPgwNm/ejEmTJsHBwQFfffUVnnvuOdy4cQP16tWrsL1IQiJRDZCRkSECEIcMGVJq3/3798XU1FTtIzc3V7tvzpw54sMf85iYGBGAOGHCBJ06ZsyYIQIQ9+7dq93m7e0tAhA3bdqkE4eHh4fYtm1b7ba8vDxRrVbr1JeQkCAqlUrxgw8+0NkGQFy1alWF73Xfvn0iAHHDhg0VlisLgDIfQ4YMEfPy8nTKhoaGllm2e/fu2jJt27YV69SpU+o82dnZOm2ekZGh3RcbGyva2NiIAMSgoCBx6tSp4tatW8WcnJzHxn///n0RgPjZZ5+VWyYrK0usU6eOOHHiRJ3td+7cEZ2cnHS2BwUFiR4eHmJ6erp2286dO0UAore3t3ZbSZvv27dPp86y/ma9e/cWW7durdOeGo1G7Ny5s9i0aVPttlWrVokAxD59+ogajUa7/c033xTlcrlOTK1atdJp94rievjzXWLBggWiIAji9evXtdse/eyXpaCgQHR1dRUDAgLEBw8eaLf/9ttvIgBx9uzZ2m0ln5cpU6bovO9nnnlGVCgUYmpqqiiKopiamioCEOfMmVPqfGXFBEBUKpViQkKCdtuyZctEAKK7u7uYmZmp3T5r1iwRgE7Z0NBQnb/l1KlTRUdHR7GoqKjc963vd/bUqVPlfmcfPe/WrVtFAOJHH32kU+75558XBUEQr1y5ovOeFQqFzrazZ8+KAMSvv/663LhJehyeoxohMzMTAGBvb19qX48ePeDi4qJ9REZGllvPH3/8AQCIiIjQ2T59+nQAKHUlmKenp07PgaOjI0aPHo0zZ87gzp07AIonU5fMo1Kr1bh37x7s7e3RvHlzREdHV/atVtngwYOxa9cu7Nq1C7/88gtmzZqFHTt24OWXXy41TKBSqbRlSx7//e9/tfszMzPLbPP33ntPp81ffvll7b5WrVohJiYGI0eOxLVr1/Dll19iyJAhcHNzw4oVKyqM3cbGBgqFAvv378f9+/fLLLNr1y6kp6djxIgRuHv3rvYhl8sREhKCffv2AQCSkpIQExOD0NBQbc8ZUDxR3t/f//ENWYa0tDTs3bsXL774IrKysrTnvnfvHvr374/Lly+XGtZ69dVXdYakunXrBrVajevXrxsUg42NjfZ5Tk4O7t69i86dO0MURZw5c6ZSdZ0+fRopKSmYNGmSztypZ555Bi1atCjzysjJkydrn5cMMxUUFGD37t0GvJtivXv31hnqKulNfu655+Dg4FBqe3x8fLl11alTBzk5Odi1a1e5ZUzxnf3jjz8gl8vxxhtv6GyfPn06RFHE9u3bdbb36dMHTZo00b4ODAyEo6Njhe+NpMfhOaoRSv7hzM7OLrVv2bJlyMrKQnJy8mOHka5fvw6ZTAY/Pz+d7e7u7qhTp06pHzI/P79SczCaNWsGoHj+hru7OzQaDb788kssXrwYCQkJOpOiTdHNrlardebEAICzs7P2CquGDRuiT58+2n3PPvss6tWrhxkzZuC3337DoEGDtPvkcrlO2Uc5ODjg3r17pbZPmjRJu4ZVWW3erFkz/PDDD1Cr1YiLi8Nvv/2GTz/9FK+++ip8fX3LPadSqcQnn3yC6dOnw83NDR07dsTAgQMxevRouLu7AygeMgVQ7tCfo6MjAGj/lk2bNi1VxtAfxytXrkAURbz//vt4//33yyyTkpKCBg0aaF83atRIZ3/dunUBoNyk8HFu3LiB2bNnY9u2baXqyMjIqFRdJW3UvHnzUvtatGiBw4cP62yTyWRo3LixzraHvw+GerSNSpLcR+cnlmyvqO0mTZqEn3/+GQMGDECDBg3Qr18/vPjii3jqqae0ZUzxnb1+/To8PT11kjwAaNmypXb/wx59z0DxZ8PQzwWZB5MmqhGcnJzg4eGB2NjYUvtK/vdZmX+0jbnA3vz58/H+++9j3Lhx+PDDD+Hs7AyZTIZp06aZ5JLkmzdvwtfXV2fbvn37Sl3e/bDevXsDAA4ePKiTND1OixYtEBMTg9u3b+skAs2aNdP+WFZ0dZdcLkfr1q3RunVrdOrUCT179sTatWsrTNSmTZuGQYMGYevWrfjzzz/x/vvvY8GCBdi7dy/atm2rbdMffvhBm0g9zJCrs8r7PDz8YwpAe+4ZM2agf//+ZR7zaEIul8vLLPdor58+1Go1+vbti7S0NLzzzjto0aIF7OzscPv2bYwZM6bGXgJfXhsZ0naurq6IiYnBn3/+ie3bt2P79u1YtWoVRo8ejTVr1gAw/3e2LMb8XJD5MGmiGuOZZ57Bt99+i5MnT6JDhw4G1eHt7Q2NRoPLly9r/wcIAMnJyUhPTy91xVdJz8LDP6qXLl0CAO1wwsaNG9GzZ0989913Osemp6ejfv36BsVZEXd391JDDw+vq1SWoqIiAGX31FVk4MCBWLduHdauXYu33367coE+ol27dgCKh80ep0mTJpg+fTqmT5+Oy5cvIygoCP/973/x448/aoc0XF1dK0y+Sv6WJT1TD7t48aLO65Len/T0dJ3tj/YOlPSyWFtbV3juytI3if/7779x6dIlrFmzRrsWF4AKh6IqUtJGFy9eLNVzd/HixVLfB41Gg/j4eG3CDJT+PlSHVcgVCgUGDRqEQYMGQaPRYNKkSVi2bBnef/99+Pn56f2drcx7KbkAIisrS6e36cKFC9r9VPNxThPVGG+//TZsbW0xbtw4JCcnl9qvz//Qnn76aQDFVyw9bNGiRQCKE7OHJSYm6lxRl5mZie+//x5BQUHaXg65XF7q3Bs2bCg1t8VYVCoV+vTpo/Mo+dEvT8lVcI9Lrh714osvwt/fHx9++CGOHz9eZplH3/uhQ4dQWFhYqlzJfLKyhoJK5ObmIi8vT2dbkyZN4ODgoF0Son///nB0dMT8+fPLPE/J0KWHhweCgoKwZs0anWGrXbt2IS4uTucYb29vyOVyHDx4UGf74sWLdV67urqiR48eWLZsWZnJ36PDpvqys7MrlbCVpaR34uE2F0VR53L6ymjXrh1cXV2xdOlSnSU3tm/fjvPnz5f6PgDAN998o3Pub775BtbW1treTFtbWwClE1BzeXQ4WSaTaa/eK3mP+n5n7ezsAOj3Xp5++mmo1Wqd9gGK1y0TBEFnWQ6qudjTRDVG06ZN8dNPP2HEiBFo3ry5dkVwURSRkJCAn376CTKZrMLLydu0aYPQ0FAsX74c6enp6N69O06ePIk1a9ZgyJAh6Nmzp075Zs2aYfz48Th16hTc3NywcuVKJCcnY9WqVdoyAwcOxAcffICxY8eic+fO+Pvvv7F27dpScz8qa9OmTdr/pT4sNDS0wrWoLl26hB9//BFAcRJy/PhxrFmzBn5+ftpVu/VlbW2NLVu2oH///ujatSuGDRumXRvp9u3b2LZtG27cuKHz4/rJJ58gKioKw4YN0/5YRUdH4/vvv4ezszOmTZtWYey9e/fWJmtWVlbYsmULkpOT8dJLLwEonrO0ZMkSjBo1Ck888QReeukluLi44MaNG/j999/RpUsX7Q/XggUL8Mwzz6Br164YN24c0tLS8PXXX6NVq1Y6vW5OTk544YUX8PXXX0MQBDRp0gS//fZbmYtxRkZGomvXrmjdujUmTpyIxo0bIzk5GceOHcOtW7cMWp8rODgYS5YswUcffQQ/Pz+4urqWOWerRYsWaNKkCWbMmIHbt2/D0dERmzZtMngejLW1NT755BOMHTsW3bt3x4gRI7RLDvj4+ODNN9/UKa9SqbBjxw6EhoYiJCQE27dvx++//453331Xuz6TjY0N/P39sX79ejRr1gzOzs4ICAhAQECAQTFW1oQJE5CWloZevXqhYcOGuH79Or7++msEBQVpe5f1/c42adIEderUwdKlS+Hg4AA7OzuEhISUGh4HgEGDBqFnz5547733cO3aNbRp0wY7d+7EL7/8gmnTpulM+qYazPwX7BFVzZUrV8SwsDDRz89PVKlUoo2NjdiiRQvx9ddfF2NiYnTKlnWJc2FhoThv3jzR19dXtLa2Fr28vMRZs2aVuiTf29tbfOaZZ8Q///xTDAwMFJVKpdiiRYtSSwHk5eWJ06dPFz08PEQbGxuxS5cu4rFjx8Tu3bvrXEZe2SUHynscOnSo3GMfLSuXy8WGDRuKr776qpicnKxTNjQ0VLSzs6swlhLp6eniBx98ILZt21a0t7cXFQqF6OXlJT7//PPir7/+qlP2yJEjYnh4uBgQECA6OTmJ1tbWYqNGjcQxY8aIV69erfA8d+/eFcPDw8UWLVqIdnZ2opOTkxgSEiL+/PPPZbZT//79RScnJ1GlUolNmjQRx4wZI54+fVqn3KZNm8SWLVuKSqVS9Pf3Fzdv3lzqcnFRLL5U/rnnnhNtbW3FunXriq+99poYGxtb5t/s6tWr4ujRo0V3d3fR2tpabNCggThw4EBx48aN2jIlSw6cOnWqVNx4ZBmBO3fuiM8884zo4OCgs+xDWWXj4uLEPn36iPb29mL9+vXFiRMnai9XfzhOfZYcKLF+/Xqxbdu2olKpFJ2dncVXXnlFvHXrlk6Zks/L1atXxX79+om2traim5ubOGfOnFKX7x89elQMDg4WFQqFzvID5S05EB4errOt5Lvy6NITZS3H8ejfcuPGjWK/fv1EV1dXUaFQiI0aNRJfe+01MSkpSVtG3++sKIriL7/8Ivr7+4tWVlY6bVzWZygrK0t88803RU9PT9Ha2lps2rSp+Nlnn+ksOVHeexbF4n9zQkNDS22n6kMQRc46I6LaZcyYMdi/f3+VrviqbcaMGYONGzdWel4ckSXhnCYiIiIiPTBpIiIiItIDkyYiIiIiPXBOExEREZEe2NNEREREpAcmTURERER64OKWlaTRaJCYmAgHB4dqcbsAIiIiqhpRFJGVlQVPT0/IZOX3JzFpqqTExMQKV2MmIiKimunmzZsV3lWCSVMlldyI8ebNm3B0dDRq3RqNBqmpqXBxcakw06WqY1ubF9vbfNjW5sX2Nh9TtnVmZia8vLx0brZcFiZNlVQyJOfo6GiSpCkvLw+Ojo788pkY29q82N7mw7Y2L7a3+ZijrR837YZ/YSIiIiI9MGkiIiIi0gOTJiIiIiI9MGkiIiIi0gOTJiIiIiI9MGkiIiIi0gOTpmpCrRFxPP4edl5Iw/H4e1BreB9lIiKi6oTrNFUDO2KTMO/XOCRl5P2zJQEeTirMGeSPpwI8JI2NiIiIirGnSWI7YpMQ9mP0QwlTsTsZeQj7MRo7YpMkioyIiIgexqRJT5GRkfD390f79u2NVqdaI2Ler3EoayCuZNu8X+M4VEdERFQNMGnSU3h4OOLi4nDq1Cmj1XkyIa1UD9PDRABJGXk4mZBmtHMSERGRYZg0SSglq/yEyZByREREZDpMmiTk6qAyajkiIiIyHSZNEurg6wwPJxUquqeym6MSHXydzRYTERERlY1Jk4TkMgFzBvkDQLmJk0wQODxHRERUDTBpkthTAR5YMvIJuDvpDsG52CvhZGOFpIw8DFt8FBfvZEkUIREREQFc3LJaeCrAA3393XEi/i6u3EqFX0MXhDSujzuZeQhdeRJXUrLxwtKjWDG6HUIa15M6XCIiolqJPU3VhFwmoGPjeujXwhkdG9eDXCagQR0bbHy9E4K96yIzrwijVp7E9r+52CUREZEUmDRVc3VsFVg7IQT9/N1QUKTBpJ+iseboNanDIiIiqnWYNNUAKms5lowMxsiOjSCKwJxt5/DJjgsQRa4UTkREZC5MmmoIuUzAh4MDMKNfMwDAkv1XMX3DWRSqNRJHRkREVDswaapBBEHA5F5N8enzgZDLBGyOvo3xa04jO79I6tCIiIgsHpOmGujFdl74dnQ72FjLcfBSKkYsP47UrHypwyIiIrJoTJpqqJ4tXPG/VzvC2U6Bv29n4LklR3Htbo7UYREREVksJk01WJBXHWwK64xGzra4kZaL55YcRczNdKnDIiIiski1Lmm6efMmevToAX9/fwQGBmLDhg1Sh1QlvvXtsCmsMwIaOOJeTgFGLD+OfRdTpA6LiIjI4tS6pMnKygpffPEF4uLisHPnTkybNg05OTV7WMvFQYl1r3ZCt6b18aBQjQlrTmPD6ZtSh0VERGRRal3S5OHhgaCgIACAu7s76tevj7S0NGmDMgJ7pRW+C22PYW0bQK0R8dbGv/DN3stcy4mIiMhIql3SdPDgQQwaNAienp4QBAFbt24tVSYyMhI+Pj5QqVQICQnByZMnDTpXVFQU1Go1vLy8qhh19aCwkuG/L7ZBWI8mAID/7LyE93+JhVrDxImIiKiqql3SlJOTgzZt2iAyMrLM/evXr0dERATmzJmD6OhotGnTBv3790dKyr/zeIKCghAQEFDqkZiYqC2TlpaG0aNHY/ny5SZ/T+YkCALeeaoF5j3bCoIA/Hj8BiatjUJeoVrq0IiIiGo0QazG4zeCIGDLli0YMmSIdltISAjat2+Pb775BgCg0Wjg5eWFKVOmYObMmXrVm5+fj759+2LixIkYNWrUY8vm5/+7BlJmZia8vLxw//59ODo6Vv5NVUCj0SA1NRUuLi6Qyaqez/7xdxIifj6LArWIYO+6WDHqCdSxVRgh0prP2G1NFWN7mw/b2rzY3uZjyrbOzMxE3bp1kZGRUeFvu5VRz2piBQUFiIqKwqxZs7TbZDIZ+vTpg2PHjulVhyiKGDNmDHr16vXYhAkAFixYgHnz5pXanpqairy8PP2D14NGo0FGRgZEUTTKB6KdmxxfDG2Kt7ddRdT1+xi2+Ai+GNIU7o5MnIzd1lQxtrf5sK3Ni+1tPqZs66ysLL3K1aik6e7du1Cr1XBzc9PZ7ubmhgsXLuhVx5EjR7B+/XoEBgZq50v98MMPaN26dZnlZ82ahYiICO3rkp4mFxcXk/Q0CYJg1Cz6KVdX+Hq6YuzqU7iWlofXNl7CqjHt0cLdwSj111SmaGsqH9vbfNjW5sX2Nh9TtrVKpdKrXI1Kmoyha9eu0Gj0v8mtUqmEUqkstV0mk5nkCyIIgtHrbunphM2TuiB05UlcTsnG8GXHsXx0O3RqUs9o56iJTNHWVD62t/mwrc2L7W0+pmprfeurUX/h+vXrQy6XIzk5WWd7cnIy3N3dJYqqZvCsY4ONr3dGBx9nZOUXIXTlSfz2V+LjDyQiIiIANSxpUigUCA4Oxp49e7TbNBoN9uzZg06dOpn03JGRkfD390f79u1Neh5TcrK1xvfjO2BAgDsK1BpM+d8ZrDycIHVYRERENUK1S5qys7MRExODmJgYAEBCQgJiYmJw48YNAEBERARWrFiBNWvW4Pz58wgLC0NOTg7Gjh1r0rjCw8MRFxeHU6dOmfQ8pqayluObl5/A6E7eEEXgg9/iMP+P89BwLSciIqIKVbs5TadPn0bPnj21r0smYYeGhmL16tUYPnw4UlNTMXv2bNy5cwdBQUHYsWNHqcnhVD65TMC8Z1vB3UmFT3dcxPKD8UjJzMOnz7eBwqra5dFERETVQrVLmnr06PHYW39MnjwZkydPNlNExSIjIxEZGQm12jIWiRQEAZN6+MHNQYV3Nv2FrTGJuJtdgCUjn4CDylrq8IiIiKoddivoyVKG5x71XHBDfBvaDrYKOQ5fuYvhy44jJcu4608RERFZAiZNhB7NXbHu1Y6ob69AXFImhi0+ivjUbKnDIiIiqlaYNBEAILBhHWwK6wzvera4df8BnltyFNE37ksdFhERUbXBpIm0vOvZYVNYZwQ2dML93EK8vOI49pxPfvyBREREtQCTJj1ZwjpN+qhvr8T/JnZEj+YuyCvUYOL3p7Hu5A2pwyIiIpIckyY9WepE8LLYKa2wYnQ7PB/cEBoRmLn5b3y5+/Jjr2okIiKyZEyaqEzWchk+ez4Qk3v6AQA+330J7275G0Vq/e/bR0REZEmYNFG5BEHAjP7N8eHgVhAE4H8nb+L1H6PxoMAy1qoiIiKqDCZN9FijOvlgySvBUFjJsPt8Ml759jju5xRIHRYREZFZMWnSU22ZCF6epwLcsXZCCBxVVoi+kY7nlh7FzbRcqcMiIiIyGyZNeqpNE8HL097HGZvCOsPTSYX41BwMW3IU5xIzpA6LiIjILJg0UaU0dXPA5kld0MLdAalZ+Ri+7DiOXLkrdVhEREQmx6SJKs3dSYX1r3VCiK8zsvOLMGbVSfwSc1vqsIiIiEyKSRMZxMnGGmvGdcAzrT1QqBYxdV0MVhyMlzosIiIik2HSpKfaPhG8LCprOb4e0RZju/gAAD7+4zw+/C0OGg0XwSQiIsvDpElPnAheNplMwOyB/pg1oAUA4LvDCZi6Pgb5RVzLiYiILAuTJqoyQRDwWvcm+Hx4G1jJBPx6NhFjVp5CZl6h1KEREREZDZMmMpqhbRti1dj2sFPIcSz+Hl5cegzJmXlSh0VERGQUTJrIqLo1dcH61zqhvr0SF+5kYdjio7iSkiV1WERERFXGpImMLqCBEzaHdYZvfTvcTn+A55ceQ9T1NKnDIiIiqhImTWQSjerZYlNYZwR51UF6biFeXnECO8/dkTosIiIigzFpIpNxtlPgp4kh6NXCFflFGrz+YxTWnrgudVhEREQGsapMYV9fXwiCUOmTTJs2DW+88Ualj6tOIiMjERkZCbWal9JXhq3CCstHBeO9LbFYf/om3tsSi+SMPLzZt5lBnyUiIiKpVCppWr16tUEn8fHxMei46iQ8PBzh4eHIzMyEk5OT1OHUKFZyGRY+1xpuTip8tecyvtp7BcmZ+fh4aACs5OzsJCKimqFSSVP37t1NFQdZOEEQENG3GdwclXh/a3GvU2p2Pr55uS1sFZX6GBIREUmiyv/N12g0xoiDaolXQryxdGQwlFYy7L2QghErTuBedr7UYRERET1WlZOmt956CzNmzNC+HjduHObOnYtffvkF2dnZVa2eLFC/Vu74aWII6tha4+zNdDy/9Bhu3MuVOiwiIqIKVTlp2rVrFz7++GPt6+PHj6NevXr48ccfMX/+/KpWTxYq2NsZG1/vjAZ1bJBwNwfDlhxF7O0MqcMiIiIqV5WTJqVSCaVSqX0dGBiIKVOm4H//+x/2799f1erJgvm52mPzpM5o6eGIu9n5GL7sGA5dTpU6LCIiojJVOWlSKBRITk7Wvl63bh0AwMrKCgUFBVWtniycm6MK61/riM5N6iGnQI2xq05hy5lbUodFRERUilHmNA0dOhQ3b97U2X737l3k53OCLz2eo8oaq8a2x6A2nijSiHhz/VksPXAVoihKHRoREZFWla/1HjJkCDIzM9GuXTt07twZAQEBEAQBP//8M2bOnGmMGKsFLm5pWkorOb4cHgQ3ByW+PZyAhdsv4E5GHt4f6A+5jItgEhGR9IyysuDo0aNx6dIlvPDCCxBFEUqlEj/++CPGjBljjOqrhfDwcMTFxeHUqVNSh2KxZDIB/zfQH//3TEsAwOqj1zDlf9HIK2SiSkRE0jPaqoJOTk54+eWXjVUd1WITujWGi4MSMzacxR9/38G97JNYProdnGyspQ6NiIhqsSr3NO3evRsdO3ZE9+7dsXv3bgBAUlISVq1ahREjRlQ5QKqdBgc1wJqxHWCvtMKJhDS8uPQYkjIeSB0WERHVYlVOmiZPnozp06dj0aJF2LJlC8aOHYvmzZtj9+7dePbZZ40RI9VSnf3q4+fXOsHVQYmLyVkYtvgoLiVnSR0WERHVUlUenrOxscELL7wAAAgKCoKLiwvi4uLQsGHDKgdH5O/piE1hnRG66iTiU3Pw/JKj+G5Me7T3cZY6NCIiqmWq3NOUmpqK9evXIzo6Gg8ePICPjw8TJjIqL2dbbHq9M55oVAeZeUV45dsT2BGbJHVYRERUy1Q5aYqIiMDOnTsRFhYGT09PnDt3DoMHD8acOXOwZcsWY8RIhLp2Cqyd0BF9WrqhoEiDsLXR+P7YNanDIiKiWqTKw3Pjx4+Hk5OT9nVCQgJiY2MRGxuLjRs3YujQoVU9BREAwEYhx9KRT+D9X87hfydvYPYv55CcmYcZ/ZpDELiWExERmVaVk6a6devCy8sLrVq1QkBAAFq3bo2AgAD069dP5550RMZgJZdh/tAAeDipsGjXJUTuu4o7GflY+FxrWMuNsuwYERFRmar8KxMWFob69eujS5cuaN68OaKiovDWW2/B19cXLVq0MEaMRDoEQcAbvZti4bDWkMsEbIq+hQlrTiMnv0jq0IiIyIJVOWmKjIzE1q1bcf36daxbtw7Dhg3D7t27kZiYiEOHDhkjRqIyvdShEZaPCobKWoYDl1IxYsVx3M3m/Q6JiMg0jDKe4eXlheXLl2P58uVYs2YNBgwYgLi4OLi4uBijeqJy9W7phv9N7Ii6ttb461YGnltyFNfu5kgdFhERWaAqz2k6f/48Ll68iIsXL+L8+fO4evUqcnJycO7cOfj7+xsjxmqBN+ytvto2qotNYZ0xeuVJXL+Xi+eWHMWqse0R2LCO1KEREZEFEURRFKtSgUwmQ2BgIF588UUMHDgQLVu2hLW15d4jLDMzE05OTsjIyICjo6NR69ZoNEhJSYGrqytkMk5qrqyUrDyMXXUK5xIzYauQY/ErT6BHc9cyy7KtzYvtbT5sa/Nie5uPKdta39/2Kp/1P//5D4KDg/HLL7+gb9++6Ny5M8aOHYv//ve/+PPPP6taPZHeXB1UWPdqR3T1q4/cAjUmrDmNjVG3pA6LiIgsRJWH5yIiInReP7xO0w8//ID+/ftX9RREenNQWWPlmPZ4e+NZbI1JxIwNZ5GcmYdJPZpwLSciIqqSKidND9NoNPD19YWvry8GDRpkzKqJ9KawkmHRi0Fwc1Rh2cF4fPbnRdzJyMPcZ1tBLmPiREREhqny8Nz9+/fx8ssvw8nJCUqlEn5+fpg5cybS09ONEB6RYWQyAbOebonZA/0hCMAPx68jfG008grVUGtEHI+/h50X0nA8/h7UmipN6yMiolqiyj1NM2fORKNGjXDr1i0olUpcu3YNK1euRKdOnbBv3z64u7sbI04ig4zr6gtXRyUi1p/FjnN38MxXh5CdX4TkzJL1nBLg4aTCnEH+eCrAQ9JYiYioeqtyT9OJEyewcOFCODg4QKFQoFmzZli4cCFmz56NuXPnGiFEoqoZGOiJNeM6QGUlw9XUnIcSpmJ3MvIQ9mM0dsQmSRQhERHVBFVOmsq77G/EiBE4ffp0VasnMooOvs6wV5XdsVoyODfv1zgO1RERUbmqnDSlpKRg48aNOH/+fKmFH3m1ElUXJxPScDe7oNz9IoCkjDycTEgzX1BERFSjVHlO0/Tp07F9+3Z89tlnuHz5Mjw9PdGqVSv4+/sjJSXFGDESVVlKVp5RyxERUe1T5aTpzTff1Hn98DpNXbt2rWr1REbh6qAyajkiIqp9Kp00jRkzBosXL4atrW2Z+7lOE1VHHXyd4eGkwp2MPJQ1a0kA4O6kQgdfZ3OHRkRENUSl5zT98MMPyM7O1r4OCwsrtSZTUVFRlQMjMia5TMCcQcU3kC5rpp0IYM4gfy5+SURE5ap00vTo/X3Xrl2LtLR/J88mJycb/Ua2RMbwVIAHlox8Au5OpYfg+vq7cZ0mIiKqUJXnND2aRAFAXh4n01L19FSAB/r6u+NE/F1cuZWKTI01/rPzEg5fvot72fmoZ6+UOkQiIqqmqrzkQFmq81ID6enpaNeuHYKCghAQEIAVK1ZIHRKZmVwmoGPjeujXwhlh3RujdQMnPChUY8WhBKlDIyKiasygpOmnn35CdHQ0CgsLjR2PyTk4OODgwYOIiYnBiRMnMH/+fNy7d0/qsEgigiDgjd5NAQDfH7uGtJzy13IiIqLardJJU7du3TBnzhy0a9cO9vb2yM3NxZw5c7B06VIcP35cZ5J4dSSXy7VX/uXn50MUxTKHGKn26NPSFf4ejsgtUOO7w/FSh0NERNVUpZOmAwcOICMjAxcuXMDq1asxffp0JCUl4d1330Xnzp3RrFmzKgV08OBBDBo0CJ6enhAEAVu3bi1VJjIyEj4+PlCpVAgJCcHJkycrdY709HS0adMGDRs2xFtvvYX69etXKWaq2R7ubVpz9DrSc9nbREREpVV6Ivjs2bMxePBgBAcHo1mzZhgxYoR2X0JCAk6fPo0zZ84YHFBOTg7atGmDcePGYdiwYaX2r1+/HhEREVi6dClCQkLwxRdfoH///rh48SJcXV0BAEFBQWUue7Bz5054enqiTp06OHv2LJKTkzFs2DA8//zzcHNzMzhmqvn6+buhhbsDLtzJwsoj1xDRt2rJPxERWZ5KJ023bt3CgAEDoFAoMGjQIAwePBi9evWCQqHQLmz5wgsvGBzQgAEDMGDAgHL3L1q0CBMnTsTYsWMBAEuXLsXvv/+OlStXYubMmQCAmJgYvc7l5uaGNm3a4NChQ3j++efLLJOfn4/8/Hzt68zMTACARqOBRqPR6zz60mg0EEXR6PVSaWW19eSeTTD5fzFYdSQB4zp7w9HGWsIILQs/2+bDtjYvtrf5mLKt9a2z0knTypUrodFocOTIEfz666+YOnUqkpKS0LdvXwwePBgDBw6Es7NpVlUuKChAVFQUZs2apd0mk8nQp08fHDt2TK86kpOTYWtrCwcHB2RkZODgwYMICwsrt/yCBQswb968UttTU1ONvrSCRqNBRkYGRFGETGaSCxvpH2W19ROuMjSup0L8vTxE7jqH8R09JY7ScvCzbT5sa/Nie5uPKds6KytLr3IGrdMkk8nQrVs3dOvWDZ9++inOnz+PX3/9FcuWLcOrr76KDh064Nlnn8WIESPQoEEDQ05Rprt370KtVpcaSnNzc8OFCxf0quP69et49dVXtRPAp0yZgtatW5dbftasWYiIiNC+zszMhJeXF1xcXIy+iKdGo4EgCHBxceGXz8TKa+tpfTV4Y10M1sekIrxfKziq2NtkDPxsmw/b2rzY3uZjyrZWqfS772iVF7cEgJYtW6Jly5Z4++23kZqaim3btmHbtm0AgBkzZhjjFEbToUMHvYfvAECpVEKpLL3goUwmM8kXRBAEk9VNuspq62cCPfHlnsu4mpqDH4/fwOReTSWM0LLws20+bGvzYnubj6naWt/6Kn3WMWPGIDc3t9z9Li4uGD9+PH755RejJ0z169eHXC5HcnKyzvbk5GS4u7sb9VyPioyMhL+/P9q3b2/S85C05DIBU/5JlL49nIDsfN5HkYiIitWoG/YqFAoEBwdjz5492m0ajQZ79uxBp06dTHLOEuHh4YiLi8OpU6dMeh6S3qA2nmhc3w7puYX4/tg1qcMhIqJqotrdsDc7OxsxMTHaIbSEhATExMTgxo0bAICIiAisWLECa9aswfnz5xEWFoacnBzt1XREVSWXCQjv6QcA+PZQAnLY20RERKiGN+w9ffo0evbsqX1dMgk7NDQUq1evxvDhw5GamorZs2fjzp07CAoKwo4dO7jOEhnV4CBPfLX3Mq7fy8WPx6/jte5NpA6JiIgkVu1u2NujRw/tlW0PP1avXq0tM3nyZFy/fh35+fk4ceIEQkJCjBB1xTinqXaxksu0vU3LD8bjQYFa4oiIiEhqte6GvYbinKbaZ2jbBvBytsG9nAKsPXFd6nCIiEhite6GvUT6spbLEN6juLdp6YF45BWyt4mIqDar9JymAwcOAAAuX76MqKgoREdHIzo6Gu+++y7S09OrNDRHVN0Me6Ihvt57BbfTH+CnEzcwrquv1CEREZFEDJ4I3rRpUzRt2hQvvfSSdpsxbthbXUVGRiIyMhJqNXsbahOFlQyTejbBe1tisfTAVbwc0ggqa7nUYRERkQSMOhG85Ga98+fPN2a11QLnNNVezwc3hIeTCilZ+fj59E2pwyEiIolwzXeix1BayTGpR/GSA0v2X0V+EXsbiYhqIyZNRHp4oZ0X3ByVSMrIw4bTt6QOh4iIJMCkiUgPKms5wrr/29tUUKSROCIiIjI3oyRNly5dMtn95qoLLm5JL3VoBBcHJW6nP8CmaPY2ERHVNkZJmlq2bIn4+HhjVFVtcSI4qazleO3JxgCAyH1XUKhmbxMRUW1ilKSprPvPEVmiV0K8Ud9egVv3H2BL9G2pwyEiIjPinCaiSrBRyPHqP71N3+y7giL2NhER1RpMmogqaWRHbzjbKXAjLRdbYxKlDoeIiMyESRNRJdkqrDCx279zm9jbRERUOzBp0hOvnqOHje7kjbq21ki4m4Nf/2JvExFRbcCkSU+8eo4eZqe0woR/epu+3nsFag0vhiAisnRMmogMNLqTN5xsrBGfmoPf2NtERGTxjJI0vfPOO6hXr54xqiKqMRxU1hjf1RdAcW+Thr1NREQWzShJ04IFC5g0Ua0U2tkHDiorXEnJxvbYO1KHQ0REJsThOaIqcLKxxrguxb1NX+25zN4mIiILxqSJqIrGdfGFg9IKF5OzsDOOvU1ERJaKSZOeuOQAlcfJ1hpjuvgAAL7cw7lNRESWyuhJ082bNzFu3DhjVys5LjlAFRnXxRd2CjnOJ2Vi9/lkqcMhIiITMHrSlJaWhjVr1hi7WqJqra6dAqM7+wAAvtxzmTexJiKyQFaVPWDbtm0V7o+Pjzc4GKKabGK3xlhz9BrOJWZi74UU9G7pJnVIRERkRJVOmoYMGQJBECr8n7QgCFUKiqgmcrZTYFRHbyw7GI+v9lxGrxau/C4QEVmQSg/PeXh4YPPmzdBoNGU+oqOjTREnUY0w8cnGUFnLcPZWBvZfSpU6HCIiMqJKJ03BwcGIiooqd//jeqGILFl9eyVGhngDAL7czblNRESWpNJJ01tvvYXOnTuXu9/Pzw/79u2rUlBENdmr3RtDaSVDzM10HLp8V+pwiIjISCqdNHXr1g1PPfVUufvt7OzQvXv3KgVFVJO5OqjwckgjALySjojIknBxSz1xcUuqjNe7N4HCSoao6/dx9Oo9qcMhIiIjYNKkJy5uSZXh5qjCiPZeAIp7m4iIqOar1JIDvr6+Bl1CPW3aNLzxxhuVPo6oJnu9RxP87+RNnExIw/H4e+jYuJ7UIRERURVUKmlavXq1QSfx8fEx6DiimszDyQYvtm+IH4/fwJe7L6Pjq0yaiIhqskolTZzgTVQ5YT38sP7UTRyLv4eTCWno4OssdUhERGQgzmkiMqEGdWzwfHDx3KavOLeJiKhGY9JEZGKTejSBlUzA4St3EXU9TepwiIjIQEyaiEzMy9kWw55oAAD4cs8ViaMhIiJDMWkiMoPwnn6QywQcvJSKMzfuSx0OEREZwKCkqbCwEL1798bly5yjQaQP73p2GBJU3Nv09V72NhER1UQGJU3W1tb466+/jB0LkUWb3MsPMgHYeyEFf91KlzocIiKqJIOH50aOHInvvvvOmLEQWTTf+nYY/E9v01ec20REVONUap2mhxUVFWHlypXYvXs3goODYWdnp7N/0aJFVQ6OyNJM7uWHrTG3sft8MmJvZyCggZPUIRERkZ4MTppiY2PxxBNPAAAuXbqks8+QW61Ud5GRkYiMjIRarZY6FKrBmrjYY1CgJ7adTcTXey9j2ah2UodERER6Mjhp2rdvnzHjqPbCw8MRHh6OzMxMODmxd4AMN6WXH379KxF/nkvG+aRMtPRwlDokIiLSg8FJEwCkp6fju+++w/nz5wEArVq1wrhx45hUEFWgqZsDnm7tgd//SsLXey9j8SvBUodERER6MHgi+OnTp9GkSRN8/vnnSEtLQ1paGhYtWoQmTZogOjramDESWZwpvfwAAH/8fQcX72RJHA0REenD4KTpzTffxLPPPotr165h8+bN2Lx5MxISEjBw4EBMmzbNiCESWZ4W7o4YEOAOAPh6L9c7IyKqCarU0/TOO+/AyurfET4rKyu8/fbbOH36tFGCI7JkU3o1BQD8/ncSrqSwt4mIqLozOGlydHTEjRs3Sm2/efMmHBwcqhQUUW3g7+mIfv5uEEWuEk5EVBMYnDQNHz4c48ePx/r163Hz5k3cvHkT69atw4QJEzBixAhjxkhksd7oXdzb9OvZRFxNzZY4GiIiqojBV8/95z//gSAIGD16NIqKigAU314lLCwMCxcuNFqARJYsoIET+rR0xe7zKYjcewWLhgdJHRIREZXD4J4mhUKBL7/8Evfv30dMTAxiYmKQlpaGzz//HEql0pgxElm0krlNW2Nu49rdHImjISKi8hiUNBUWFqJ37964fPkybG1t0bp1a7Ru3Rq2trbGjo/I4rXxqoMezV2gEYHIfZzbRERUXRmUNFlbW+Ovv/4ydixEtVbJ3KbNZ27jxr1ciaMhIqKyGDw8N3LkSHz33XfGjIWo1nqiUV10a1ofao2IxfvZ20REVB0ZPBG8qKgIK1euxO7duxEcHAw7Ozud/YsWLapycES1ybQ+TXHo8l1sjLqF8J5+8HLmcDcRUXVicNIUGxuLJ554AgBw6dIlnX2CIFQtKqJaKNjbGV386uHIlXtYcuAq5g9tLXVIRET0EIOTpn379hkzDiICMLV3Mxy5cgwbTt/E5J5+8KxjI3VIRET0jypfPVdT5ebmwtvbGzNmzJA6FCKtDr7O6NjYGYVqEUv2X5U6HCIiekitvXru448/RseOHaUOg6iUkivp1p+6iTsZeRJHQ0REJWrl1XOXL1/GhQsXMGDAAKlDISqlU+N66ODjjAK1BksPsLeJiKi6MDhpKioqwpIlS9CuXTu89tpriIiI0HkY6uDBgxg0aBA8PT0hCAK2bt1aqkxkZCR8fHygUqkQEhKCkydPVuocM2bMwIIFCwyOkciUBEHQ9jb9dPIGUjLZ20REVB1Uu6vncnJy0KZNG4wbNw7Dhg0rtX/9+vWIiIjA0qVLERISgi+++AL9+/fHxYsX4erqCgAICgrS3g/vYTt37sSpU6fQrFkzNGvWDEePHjU4TiJT6uJXD8HedRF1/T6WHojH7EH+UodERFTrVbur5wYMGFDhsNmiRYswceJEjB07FgCwdOlS/P7771i5ciVmzpwJAIiJiSn3+OPHj2PdunXYsGEDsrOzUVhYCEdHR8yePduo74OoKkp6m0JXnsTaE9fxeo/GcHVQSR0WEVGtZnDSBACHDh3CsmXLEB8fjw0bNqBBgwb44Ycf4Ovri65duxorRq2CggJERUVh1qxZ2m0ymQx9+vTBsWPH9KpjwYIF2qG51atXIzY2tsKEKT8/H/n5+drXmZmZAACNRgONRmPI2yiXRqOBKIpGr5dKqwlt3bWJM4K8nBBzMwMrDsZj1oAWUodksJrQ3paCbW1ebG/zMWVb61unwUnTpk2bMGrUKLzyyiuIjo7WJhYZGRmYP38+/vjjD0OrLtfdu3ehVqvh5uams93NzQ0XLlww+vmA4iRr3rx5pbanpqYiL8+4c000Gg0yMjIgiiJkMoOnm5Eeakpbj37CBTE3M/DDsesY2tIBzrbWUodkkJrS3paAbW1ebG/zMWVbZ2Vl6VXO4KTpo48+wtKlSzF69GisW7dOu71Lly746KOPDK3WrMaMGfPYMrNmzdKZ2J6ZmQkvLy+4uLjA0dHRqPFoNBoIggAXFxd++UysprT1YBcXrDqdir9vZ+CXC9l456nmUodkkJrS3paAbW1ebG/zMWVbq1T6TX8wOGm6ePEinnzyyVLbnZyckJ6ebmi1Fapfvz7kcjmSk5N1ticnJ8Pd3d0k51QqlVAqlaW2y2Qyk3xBBEEwWd2kq6a09dTeTTHh+9P44fh1vNa9CZztFFKHZJCa0t6WgG1tXmxv8zFVW+tbn8FndXd3x5Urpe/GfvjwYTRu3NjQaiukUCgQHByMPXv2aLdpNBrs2bMHnTp1Msk5S0RGRsLf3x/t27c36XmIHtW7pStaeToit0CN7w7HSx0OEVGtZXDSNHHiREydOhUnTpyAIAhITEzE2rVrMWPGDISFhRkcUHZ2NmJiYrRXwCUkJCAmJgY3btwAAERERGDFihVYs2YNzp8/j7CwMOTk5GivpjOV8PBwxMXF4dSpUyY9D9GjHl63ac3R60jPLZA4IiKi2sng4bmZM2dCo9Ggd+/eyM3NxZNPPgmlUokZM2ZgypQpBgd0+vRp9OzZU/u6ZD5RaGgoVq9ejeHDhyM1NRWzZ8/GnTt3EBQUhB07dpSaHE5kSfq2dEMLdwdcuJOFlYcTENGvZs5tIiKqyQRRFMWqVFBQUIArV64gOzsb/v7+sLe3N1Zs1VJmZiacnJyQkZFhkongKSkpcHV15di4idXEtt7+dxLC1kbDQWmFwzN7wcmm5lxJVxPbu6ZiW5sX29t8TNnW+v62V/msCoUC/v7+6NChg0UnTJzTRFLr38odzdzskZVfhFVHEqQOh4io1mFarCfOaSKpyWQCpvQqntu08nACMvMKJY6IiKh2YdJEVIM83doDfq72yMwrwpoj16QOh4ioVmHSRFSDyGUCpvTyAwB8ezgB2fmlb0xNRESmwaRJT5zTRNXFwEBPNK5vh4wHhVhz9JrU4RAR1RoGJ02hoaE4ePCgMWOp1jiniaoLuUzA5JLepkPxyGFvExGRWRicNGVkZKBPnz5o2rQp5s+fj9u3bxszLiKqwLNtPOFTzxb3cwvx4/HrUodDRFQrGJw0bd26Fbdv30ZYWBjWr18PHx8fDBgwABs3bkRhIa/qITIlK7kM4T2Le5uWH4xHbgF7m4iITK1Kc5pcXFwQERGBs2fP4sSJE/Dz88OoUaPg6emJN998E5cvXzZWnET0iCFtG8DL2Qb3cgrw04kbUodDRGTxjDIRPCkpCbt27cKuXbsgl8vx9NNP4++//4a/vz8+//xzY5xCcpwITtWNtVyGyf/0Ni09EI8HBWqJIyIismwGJ02FhYXYtGkTBg4cCG9vb2zYsAHTpk1DYmIi1qxZg927d+Pnn3/GBx98YMx4JcOJ4FQdDW3bEA3q2OBudj7+d5K9TUREpmTwDXs9PDyg0WgwYsQInDx5EkFBQaXK9OzZE3Xq1KlCeERUEYVV8dymd7f8jaUHruLlkEZQWculDouIyCIZ3NP0+eefIzExEZGRkWUmTABQp04dJCTwHllEpvR8cEN4OqmQkpWP9aduSh0OEZHFMjhpGjVqFFQqlTFjISIDKKxkCPtnbtOS/VeRX8S5TUREpmDw8FxERESZ2wVBgEqlgp+fHwYPHgxnZ2eDgyMi/bzYriEi917Bncw8/Hz6FkZ19JY6JCIii2Nw0nTmzBlER0dDrVajefPmAIBLly5BLpejRYsWWLx4MaZPn47Dhw/D39/faAFLJTIyEpGRkVCr+b94qn6UVnK83r0x5v4ahyX7rmB4Oy8orHiXJCIiYzL4X9XBgwejT58+SExMRFRUFKKionDr1i307dsXI0aMwO3bt/Hkk0/izTffNGa8kuHVc1TdvdShEVwdlEjMyMPGqFtSh0NEZHEMTpo+++wzfPjhh3B0dNRuc3Jywty5c/Hpp5/C1tYWs2fPRlRUlFECJaKKqazleK17EwBA5L4rKFRrJI6IiMiyVOnecykpKaW2p6amIjMzE0Dx1XMFBQWGR0dElfJyh0aob6/E7fQH2BzN3iYiImOq0vDcuHHjsGXLFty6dQu3bt3Cli1bMH78eAwZMgQAcPLkSTRr1sxYsRLRY9go5HjtycYAgG/Y20REZFQGJ03Lli1D79698dJLL8Hb2xve3t546aWX0Lt3byxduhQA0KJFC3z77bdGC5aIHu+Vjo1Qz06Bm2kP8EtMotThEBFZDIOSpsLCQgwePBhvv/027t27hzNnzuDMmTO4d+8eli9fDjs7OwBAUFBQuQtfEpFp2CqsMLGkt2nvZRSxt4mIyCgMSpqsra3x119/AQDs7e0RGBiIwMBA2NvbGzW46oQ37KWaZFRHb9S1tca1e7n49S/2NhERGYPBw3MjR47Ed999Z8xYqjUuOUA1iZ3SChO6Ffc2fb33CtQaUeKIiIhqPoMXtywqKsLKlSuxe/duBAcHa4fkSixatKjKwRGR4UZ38sbyg/GIT83Bb38lYnBQA6lDIiKq0QxOmmJjY/HEE08AKF4J/GGCIFQtKiKqMgeVNSZ09cV/d13C13uvYGCgJ+QyfjeJiAxlcNK0b98+Y8ZBRCYQ2sUHKw7F40pKNrbHJmFgoKfUIRER1Vi8ORWRBXNUWWNcV18AwNd7rkDDuU1ERAarUtJ06NAhjBw5Ep06dcLt27cBAD/88AMOHz5slOCIqOrGdvaFg9IKF5Oz8Oe5O1KHQ0RUYxmcNG3atAn9+/eHjY0Nzpw5g/z8fADFt1eZP3++0QIkoqpxsrXGmC4+AIAv91xmbxMRkYEMTpo++ugjLF26FCtWrIC1tbV2e5cuXRAdHW2U4IjIOMZ39YWdQo4Ld7Kw63yy1OEQEdVIBidNFy9exJNPPllqu5OTE9LT06sSU7XExS2pJqtjq0BoZx8AwFd7LkMU2dtERFRZBidN7u7uuHLlSqnthw8fRuPGjasUVHXExS2pppvQrTFsFXKcS8zEnvMpUodDRFTjGJw0TZw4EVOnTsWJEycgCAISExOxdu1azJgxA2FhYcaMkYiMwNlOgVGdvAEAX+1lbxMRUWUZvE7TzJkzodFo0Lt3b+Tm5uLJJ5+EUqnEjBkzMGXKFGPGSERGMrFbY3x/9Dr+upWB/ZdS0bO5q9QhERHVGAb3NAmCgPfeew9paWmIjY3F8ePHkZqaig8//NCY8RGREdW3V2Jkx0YAgC93s7eJiKgyqry4pUKhgL+/Pzp06AB7e3tjxEREJjTxycZQWskQczMdhy7flTocIqIaw+DhOQDYs2cP9uzZg5SUFGg0Gp19K1eurFJgRGQarg4qvBLijZVHEvDlnsvo1rQ+7xdJRKQHg3ua5s2bh379+mHPnj24e/cu7t+/r/Mgourrte6NobCSIer6fRy9ek/qcIiIagSDe5qWLl2K1atXY9SoUcaMh4jMwM1RhZc7NMLqo9fw5e7L6NykHnubiIgew+CepoKCAnTu3NmYsRCRGb3WvTEUchlOXkvD8fg0qcMhIqr2DE6aJkyYgJ9++smYsRCRGXk42WB4ey8AxauEExFRxQwensvLy8Py5cuxe/duBAYG6tx/DgAWLVpU5eCIyLRe79EE607dwLH4eziZkIYOvs5Sh0REVG0ZnDT99ddfCAoKAgDExsYaKx4iMqMGdWzwQjsv/HTiBr7acxk/TgiROiQiomrL4KRp3759xoyj2ouMjERkZCTUarXUoRAZVVj3Jvj51E0cvnIXUdfTEOzN3iYiorJUek7T008/jYyMDO3rhQsXIj09Xfv63r178Pf3N0pw1Qlv2EuWysvZFs890RAA8OWe0jfhJiKiYpVOmv7880/k5+drX8+fPx9paf9eeVNUVISLFy8aJzoiMovwnn6QywQcvJSKMze4zhoRUVkqnTQ9eq8q3ruKqOZrVM8WQ9s2AMAr6YiIylPle88RkWWY3NMPMgHYdzEVf91KlzocIqJqp9JJkyAIpVYO5krCRDWfT307DAlibxMRUXkqffWcKIoYM2YMlEolgOL1ml5//XXY2dkBgM58JyKqWcJ7+WFrzG3sPp+C2NsZCGjgJHVIRETVRqV7mkJDQ+Hq6gonJyc4OTlh5MiR8PT01L52dXXF6NGjTRErEZlYExd7DGrjCYC9TUREj6p0T9OqVatMEQcRVRNTevlh29lE7IxLRlxiJvw9HaUOiYioWuBEcCLS4efqgGdaewAAvt7L3iYiohJMmoiolCm9mgIAtsfewcU7WRJHQ0RUPTBpIqJSmrs74OnW7gDY20REVIJJExGVqaS36fe/k3A5mb1NRERMmoioTC09HNG/lRtEEfhmH+9JR0TEpImIylXS2/Tr2URcTc2WOBoiImkxaSKicgU0cEKflq7QiEDkXvY2EVHtViuTJh8fHwQGBiIoKAg9e/aUOhyiau2N3sW9TVtjbiPhbo7E0RARSadWJk0AcPToUcTExGDfvn1Sh0JUrQU2rIOezV2Ke5s4t4mIarFamzQRkf5Kepu2nLmNG/dyJY6GiEga1S5pOnjwIAYNGgRPT08IgoCtW7eWKhMZGQkfHx+oVCqEhITg5MmTlTqHIAjo3r072rdvj7Vr1xopciLL1bZRXTzZzAVqjcjeJiKqtapd0pSTk4M2bdogMjKyzP3r169HREQE5syZg+joaLRp0wb9+/dHSkqKtkxQUBACAgJKPRITEwEAhw8fRlRUFLZt24b58+fjr7/+Mst7I6rJpv7T27Qp+hZuprG3iYhqn0rfsNfUBgwYgAEDBpS7f9GiRZg4cSLGjh0LAFi6dCl+//13rFy5EjNnzgQAxMTEVHiOBg0aAAA8PDzw9NNPIzo6GoGBgWWWzc/PR35+vvZ1ZmYmAECj0UCj0ej9vvSh0WggiqLR66XS2NaV19bLCV386uHIlXtYvO8KPh4aoPexbG/zYVubF9vbfEzZ1vrWWe2SpooUFBQgKioKs2bN0m6TyWTo06cPjh07plcdOTk50Gg0cHBwQHZ2Nvbu3YsXX3yx3PILFizAvHnzSm1PTU1FXl5e5d9EBTQaDTIyMiCKImSyatcJaFHY1oYZ1bY4adoQdQvDW9eBu6NCr+PY3ubDtjYvtrf5mLKts7L0u+tBjUqa7t69C7VaDTc3N53tbm5uuHDhgl51JCcnY+jQoQAAtVqNiRMnon379uWWnzVrFiIiIrSvMzMz4eXlBRcXFzg6OhrwLsqn0WggCAJcXFz45TMxtrVh+rm6olPUXRyLT8OGcxn4cHArvY5je5sP29q82N7mY8q2VqlUepWrUUmTMTRu3Bhnz57Vu7xSqYRSqSy1XSaTmeQLIgiCyeomXWxrw7zRuxmOxR/HhtO3MLmXHzycbPQ6ju1tPmxr82J7m4+p2lrf+mrUX7h+/fqQy+VITk7W2Z6cnAx3d3eJoiKqXTo1qYcOvs4oUGuw7EC81OEQEZlNjUqaFAoFgoODsWfPHu02jUaDPXv2oFOnTiY9d2RkJPz9/SscyiOqLUqupPvp5A0kZxp3bh8RUXVV7ZKm7OxsxMTEaK+AS0hIQExMDG7cuAEAiIiIwIoVK7BmzRqcP38eYWFhyMnJ0V5NZyrh4eGIi4vDqVOnTHoeopqgc5N6aOddFwVF7G0iotqj2iVNp0+fRtu2bdG2bVsAxUlS27ZtMXv2bADA8OHD8Z///AezZ89GUFAQYmJisGPHjlKTw4nIdARB0K4SvvbEdaRksbeJiCxftZsI3qNHD4iiWGGZyZMnY/LkyWaKqFhkZCQiIyOhVqvNel6i6qpb0/oI8qqDmJvpWHEwHu894y91SEREJlXtepqqKw7PEekSBAFT+xT3Nv1w/DruZuc/5ggiopqNSRMRGaxHMxcENnRCXqEGKw5xbhMRWTYmTURkMEEQtFfS/XDsOtJyCiSOiIjIdJg0EVGV9GrhioAGjsgtUONb9jYRkQVj0qQnrtNEVDZBEPBGr+LepjVHryE9l71NRGSZmDTpiRPBicrX198NLT0ckVOgxneHE6QOh4jIJJg0EVGVFfc2+QEAVh+5hozcQokjIiIyPiZNRGQU/Vu5o7mbA7Lyi7DyCHubiMjyMGkiIqOQyQRM6V3c27TySAIy89jbRESWhUmTnjgRnOjxng7wQFNXe2TlFWHNkWtSh0NEZFRMmvTEieBEjyeTCZj8z9ymbw8nIIu9TURkQZg0EZFRDQz0RGMXO2Q8KMT3x65LHQ4RkdEwaSIio5LLBEwp6W06FI+c/CKJIyIiMg4mTURkdIMCPeFTzxb3cwux5tg1HI+/h50X0nA8/h7UGlHq8IiIDGIldQA1RWRkJCIjI6FWq6UOhajas5LLMLlXU8zYcBb/+fMi/s2TEuDhpMKcQf54KsBDyhCJiCqNPU164kRwosqxsS7+5+XRjqU7GXkI+zEaO2KTJIiKiMhwTJqIyOjUGhEf/X6+zH0lOdS8X+M4VEdENQqH54jI6E4mpCEpI6/c/SKApIw8fPbnRbT3qYs6tgo42ylQ19YajipryGSC+YIlItITkyYiMrqUrPITpoctPXAVSw/obpMJQB1bBerYWqOureKfhzXq2v37/OEkq6SstZwd50RkWkyaiMjoXB1UepVr3cARgID7uQW4n1OAnAI1NCKQllOAtJwCADl6n9NBZQVnOwXq/JNYOdv++/zhhKvkeR1ba6is5Ya9QSKqlZg0EZHRdfB1hoeTCncy8lDWrCUBgLuTClvDu0L+0FBcfpEaGbmFSMstwP2cQqTnFiAttwDpuYW4n/PQ83+SrPu5hch4ULzqeFZeEbLyinD9Xq7ecdoq5NoE6uGEq3Tv1r9lbBVyCEL1HD5Ua0SciL+HK7fS4JctR0jj+jrtS1RTVZfPNpMmPXHJASL9yWUC5gzyR9iP0RAAncSp5J+5OYP8S/2jp7SSw9VRDldH/XqqAKBIrUHGg0Lcz/0nycr5N7FKyy1Aes4/SVZugbbM/dxCqDUicgvUyC14gNvpD/Q+n0IueyjJ+ifBsns42VKgrt3DQ4sKOKisTD5Pa0dsEub9GvfQXDIu70CWoTp9tgVRFHn5SiVkZmbCyckJGRkZcHR0NGrdGo0GKSkpcHV1hUzG+RmmxLY2j9L/2KFa/JCLoojMvKLSSdZDz+//09v1cMJVUKQx6HxymYA6Ntale7QeHjZ8JPlysrGGlZ7ztHbEJiHsx+hSvXoladqSkU8wcTIR/ltiWub6bOv7286eJiIymacCPNDX3x0n4u/iyq1U+DV0qRZDRoIgwMnGGk421vCuZ6fXMaIo4kGhuvwk65/hwocTrvTc4nlaao2IezkFuJdTgKup+s/TcnxknlZZc7McVVb4v62xZQ6Diij+cZn3axz6+rtL3u6WproMGVkqtUbEvF/jqtVnm0kTEZmUXCagY+N6aGyvhqtrvRq7nIAgCLBVWMFWYYWGdfU/Lr9I/VBi9W9SlZ5biLQc3eclQ4cl87Qy84qQmVcEVGKe1qNKlnd49pvDcHVQQmElg9JKDqWV7N/n1jIo5DIorYtfF29/+KFfeSuZUG3nexlbdRoyqk5EUUR+kQYFag0Kioofhf88f3R7wSOv8x/Zl3A3W6+lS04mpKFTk3pmeX9MmoiITEhpJYeboxxuBs7TKunFSi+ZIP/PPK3iSfEFuJGWi+TM/MfWeS4xE+eq8kb0IBNQRlL2SCJmLX8o4SorKXu0vAwKufyh58V1lFveSmbyxK28IaOS1e7NNRwqiiIK1WIZiYi6OEEpIzEpUGsq3FdqW7mJjbrMcoVq88/40XeJE2Ng0kREVM1YyWWoZ69EPXvlY8seu3oPI1Ycf2y5Kb380MjZVvuDmV+kQX6RuuznhSU/rup/nxf+8/qR4x/+kdSIQF6hBnmFhs39MhaFXFapJKsyvW/WMgHvbil/OBQA3t38NzQaEYUasVTSoe11qSAx0Sux+ed5dWclE6D4p10VcpnO85L21t0n1z5Py87Hn3HJjz2HvkucGAOTJiKiGkzf5R2m9WlmknkfGo34b1KlVv+TXJUkVqWTrH8TM/VDyVhJYla6fL72eRnHPLT/YQXq4vqyHt8BZxJpuYWY9NMZs59XLhP0TEx0nyutZLCWP7K/zOPlpY5XllGfdclxclmVhuPVGhFdP9n72M92B19ng89RWUyaiIhqMEOXdzAWmUyASib/Z6FQa5Oc43FKhqkqTtL+GbYqI6nLfzhJe6Rn7dG6UjLzcUuPJSp869vB3VGlk4Qo5f8kJ2UkMmUlH48mNjpJi5UM1nIByoe2WdokdKk/22Vh0kREVMM9FeCBJSOfKLW8g3stmZgsCAIUVsXDQA4mPpe+w6Hzh7Y22+RkS1bdPttMmoiILEB1Xd7B0ug7HGrOISNLV50+21yJS0+RkZHw9/dH+/btpQ6FiKhMJcs79GvhjI6N6zFhMoGSISPg3yGiElINGdUG1eWzzaRJT+Hh4YiLi8OpU6ekDoWIiCRUMmTk7qR71Za7k4qrr1s4Ds8RERFVUnUaMiLzYdJERERkAEtZ7Z70x+E5IiIiIj0waSIiIiLSA5MmIiIiIj0waSIiIiLSA5MmIiIiIj0waSIiIiLSA5ccqCRRLF44PzMz0+h1azQaZGVlQaVSQSZjPmtKbGvzYnubD9vavNje5mPKti75TS/5jS8Pk6ZKysrKAgB4eXlJHAkREREZU1ZWFpycnMrdL4iPS6tIh0ajQWJiIhwcHCAIxl3ILDMzE15eXrh58yYcHR2NWjfpYlubF9vbfNjW5sX2Nh9TtrUoisjKyoKnp2eFvVjsaaokmUyGhg0bmvQcjo6O/PKZCdvavNje5sO2Ni+2t/mYqq0r6mEqwQFYIiIiIj0waSIiIiLSA5OmakSpVGLOnDlQKpVSh2Lx2NbmxfY2H7a1ebG9zac6tDUnghMRERHpgT1NRERERHpg0kRERESkByZNRERERHpg0lQNHDx4EIMGDYKnpycEQcDWrVulDsliLViwAO3bt4eDgwNcXV0xZMgQXLx4UeqwLNKSJUsQGBioXVOlU6dO2L59u9Rh1QoLFy6EIAiYNm2a1KFYpLlz50IQBJ1HixYtpA7Lot2+fRsjR45EvXr1YGNjg9atW+P06dNmj4NJUzWQk5ODNm3aIDIyUupQLN6BAwcQHh6O48ePY9euXSgsLES/fv2Qk5MjdWgWp2HDhli4cCGioqJw+vRp9OrVC4MHD8a5c+ekDs2inTp1CsuWLUNgYKDUoVi0Vq1aISkpSfs4fPiw1CFZrPv376NLly6wtrbG9u3bERcXh//+97+oW7eu2WPhiuDVwIABAzBgwACpw6gVduzYofN69erVcHV1RVRUFJ588kmJorJMgwYN0nn98ccfY8mSJTh+/DhatWolUVSWLTs7G6+88gpWrFiBjz76SOpwLJqVlRXc3d2lDqNW+OSTT+Dl5YVVq1Zpt/n6+koSC3uaqFbLyMgAADg7O0sciWVTq9VYt24dcnJy0KlTJ6nDsVjh4eF45pln0KdPH6lDsXiXL1+Gp6cnGjdujFdeeQU3btyQOiSLtW3bNrRr1w4vvPACXF1d0bZtW6xYsUKSWNjTRLWWRqPBtGnT0KVLFwQEBEgdjkX6+++/0alTJ+Tl5cHe3h5btmyBv7+/1GFZpHXr1iE6OhqnTp2SOhSLFxISgtWrV6N58+ZISkrCvHnz0K1bN8TGxsLBwUHq8CxOfHw8lixZgoiICLz77rs4deoU3njjDSgUCoSGhpo1FiZNVGuFh4cjNjaWcxFMqHnz5oiJiUFGRgY2btyI0NBQHDhwgImTkd28eRNTp07Frl27oFKppA7H4j08nSIwMBAhISHw9vbGzz//jPHjx0sYmWXSaDRo164d5s+fDwBo27YtYmNjsXTpUrMnTRyeo1pp8uTJ+O2337Bv3z40bNhQ6nAslkKhgJ+fH4KDg7FgwQK0adMGX375pdRhWZyoqCikpKTgiSeegJWVFaysrHDgwAF89dVXsLKyglqtljpEi1anTh00a9YMV65ckToUi+Th4VHqP1otW7aUZEiUPU1Uq4iiiClTpmDLli3Yv3+/ZJMJayuNRoP8/Hypw7A4vXv3xt9//62zbezYsWjRogXeeecdyOVyiSKrHbKzs3H16lWMGjVK6lAsUpcuXUotDXPp0iV4e3ubPRYmTdVAdna2zv9QEhISEBMTA2dnZzRq1EjCyCxPeHg4fvrpJ/zyyy9wcHDAnTt3AABOTk6wsbGRODrLMmvWLAwYMACNGjVCVlYWfvrpJ+zfvx9//vmn1KFZHAcHh1Lz8uzs7FCvXj3O1zOBGTNmYNCgQfD29kZiYiLmzJkDuVyOESNGSB2aRXrzzTfRuXNnzJ8/Hy+++CJOnjyJ5cuXY/ny5eYPRiTJ7du3TwRQ6hEaGip1aBanrHYGIK5atUrq0CzOuHHjRG9vb1GhUIguLi5i7969xZ07d0odVq3RvXt3cerUqVKHYZGGDx8uenh4iAqFQmzQoIE4fPhw8cqVK1KHZdF+/fVXMSAgQFQqlWKLFi3E5cuXSxKHIIqiaP5UjYiIiKhm4URwIiIiIj0waSIiIiLSA5MmIiIiIj0waSIiIiLSA5MmIiIiIj0waSIiIiLSA5MmIiIiIj0waSIiIiLSA5MmIqqWBEHA1q1bTXqO/fv3QxAEpKenV7muuXPnIigoqMr16KNHjx6YNm2aWc5FRP9i0kREZnfnzh1MmTIFjRs3hlKphJeXFwYNGoQ9e/ZoyyQlJWHAgAEmjaNz585ISkqCk5MTAGD16tWoU6eOSc9ZGcZM6oio6njDXiIyq2vXrqFLly6oU6cOPvvsM7Ru3RqFhYX4888/ER4ejgsXLgAA3N3dK6ynsLAQ1tbWVYpFoVA89jxERCXY00REZjVp0iQIgoCTJ0/iueeeQ7NmzdCqVStERETg+PHj2nIPD89du3YNgiBg/fr16N69O1QqFdauXQsAWLlyJVq1agWlUgkPDw9MnjxZ55iYmBhtnenp6RAEAfv37weg25Ozf/9+jB07FhkZGRAEAYIgYO7cueW+j4ULF8LNzQ0ODg4YP3488vLySpX59ttv0bJlS6hUKrRo0QKLFy/W7iuJb926dejcuTNUKhUCAgJw4MAB7f6ePXsCAOrWrQtBEDBmzBjt8RqNBm+//TacnZ3h7u5eYaxEZCSS3CaYiGqle/fuiYIgiPPnz39sWQDili1bRFEUxYSEBBGA6OPjI27atEmMj48XExMTxcWLF4sqlUr84osvxIsXL4onT54UP//8c51jzpw5o63z/v37IgBx3759oiiK4r59+0QA4v3798X8/Hzxiy++EB0dHcWkpCQxKSlJzMrKKjO29evXi0qlUvz222/FCxcuiO+9957o4OAgtmnTRlvmxx9/FD08PLTxbtq0SXR2dhZXr16tE1/Dhg3FjRs3inFxceKECRNEBwcH8e7du2JRUZG4adMmEYB48eJFMSkpSUxPTxdFURS7d+8uOjo6inPnzhUvXbokrlmzRhQEQdy5c2fl/iBEVClMmojIbE6cOCECEDdv3vzYsmUlTV988YVOGU9PT/G9994r8/jKJk2iKIqrVq0SnZycHhtbp06dxEmTJulsCwkJ0UmamjRpIv700086ZT788EOxU6dOOvEtXLhQu7+wsFBs2LCh+Mknn5QZX4nu3buLXbt21dnWvn178Z133nls7ERkOA7PEZHZiKJYpePbtWunfZ6SkoLExET07t27qmFV2vnz5xESEqKzrVOnTtrnOTk5uHr1KsaPHw97e3vt46OPPsLVq1fLPc7Kygrt2rXD+fPnHxtDYGCgzmsPDw+kpKQY8naISE+cCE5EZtO0aVMIgqCd7F1ZdnZ22uc2NjYVlpXJiv9P+HCiVlhYaNB5Kys7OxsAsGLFilLJlVwuN8o5Hp0ELwgCNBqNUeomorKxp4mIzMbZ2Rn9+/dHZGQkcnJySu2vzKX1Dg4O8PHx0Vmm4GEuLi4AipcuKPHwpPCyKBQKqNXqx567ZcuWOHHihM62hyexu7m5wdPTE/Hx8fDz89N5+Pr6lntcUVERoqKi0LJlS208APSKiYhMjz1NRGRWkZGR6NKlCzp06IAPPvgAgYGBKCoqwq5du7BkyRK9hqZKzJ07F6+//jpcXV0xYMAAZGVl4ciRI5gyZQpsbGzQsWNHLFy4EL6+vkhJScH//d//VVifj48PsrOzsWfPHrRp0wa2trawtbUtVW7q1KkYM2YM2rVrhy5dumDt2rU4d+4cGjdurC0zb948vPHGG3BycsJTTz2F/Px8nD59Gvfv30dERIROezRt2hQtW7bE559/jvv372PcuHEAAG9vbwiCgN9++w1PP/00bGxsYG9vr3f7EJFxsaeJiMyqcePGiI6ORs+ePTF9+nQEBASgb9++2LNnD5YsWVKpukJDQ/HFF19g8eLFaNWqFQYOHIjLly9r969cuRJFRUUIDg7GtGnT8NFHH1VYX+fOnfH6669j+PDhcHFxwaefflpmueHDh+P999/H22+/jeDgYFy/fh1hYWE6ZSZMmIBvv/0Wq1atQuvWrdG9e3esXr26VE/TwoULsXDhQrRp0waHDx/Gtm3bUL9+fQBAgwYNMG/ePMycORNubm7a5RSISBqCWNWZmUREVGnXrl2Dr68vzpw5Y7bbrxBR1bCniYiIiEgPTJqIiIiI9MDhOSIiIiI9sKeJiIiISA9MmoiIiIj0wKSJiIiISA9MmoiIiIj0wKSJiIiISA9MmoiIiIj0wKSJiIiISA9MmoiIiIj0wKSJiIiISA//D1oqnrUkCiWDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6, 4))\n",
    "plt.plot(depths_global, errors_global, marker=\"o\", linestyle=\"-\")\n",
    "plt.xlabel(\"Circuit depth\")\n",
    "plt.ylabel(\"Energy error $|1 - E/E_{DMRG}|$\")\n",
    "plt.title(\"Global L-BFGS sequential optimisation\")\n",
    "plt.yscale(\"log\")\n",
    "plt.grid(visible=True, alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "870630bc",
   "metadata": {},
   "source": [
    "## Local Optimization\n",
    " Ref: Causer et al. [Phys. Rev. Research 6, 033062 (2024)](https://doi.org/10.1103/PhysRevResearch.6.033062), Gibbs and Cincio, [Quantum 9, 1789 (2025)](https://doi.org/10.22331/q-2025-07-09-1789).\n",
    "\n",
    "In this approach, the circuit is represented as a tensor network (each two-qubit gate is a 4-leg tensor). We optimize the tensor entries directly using a local (sweeping) optimization based on the polar decomoposition that maximize the fidelity.\n",
    "We use the function `tn_fit` from `qpe_toolbox.circuit`, which performs a local optimisation of the tensor network.\n",
    "\n",
    "We first optimise a fixed depth-6 circuit, and then demonstrate a **sequential depth optimisation** that starts from depth 1 and increases depth, reusing the tensor entries from the previous depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "16329a26",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T06:02:12.936817Z",
     "iopub.status.busy": "2026-09-15T06:02:12.936667Z",
     "iopub.status.idle": "2026-09-15T06:03:02.390218Z",
     "shell.execute_reply": "2026-09-15T06:03:02.389428Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Local optimization\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  6   Energy = -10.45626977   Error =  1.797e-05   1-F =  1.794e-05\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# --- Fixed depth 6 ---\n",
    "print(\"*** Local optimization\")\n",
    "depth = 6\n",
    "circ = ansatz_circuit_su4(\n",
    "    n_qubits=n_qubits, depth=depth, param_scaling=1.0, parametrize=False\n",
    ")\n",
    "\n",
    "tn = circ.psi\n",
    "tn_fit(tn, GS, tags=\"SU4SWAP\", steps=100000, tol=1e-8)\n",
    "\n",
    "ovlp = (dmrg.state.H & tn).contract()\n",
    "tnH = tn.H\n",
    "tn.align_(mpo, tnH)\n",
    "energy_tn = tnH & mpo & tn\n",
    "ene = autoray.do(\"real\", energy_tn.contract(all))\n",
    "\n",
    "print(\n",
    "    f\"Depth = {depth:2d}   Energy = {ene:12.8f}   Error = {np.abs(1 - ene / dmrg_energy):10.3e}   1-F = {1 - np.abs(ovlp) ** 2:10.3e}\"\n",
    ")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67c4717c",
   "metadata": {},
   "source": [
    "### Sequential depth optimization (local)\n",
    "\n",
    "Starting from a depth-1 circuit, we optimise it, then add a layer, reusing the tensor entries from the previous circuit as initialisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "389cec2d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T06:03:02.391753Z",
     "iopub.status.busy": "2026-09-15T06:03:02.391587Z",
     "iopub.status.idle": "2026-09-15T06:03:37.775890Z",
     "shell.execute_reply": "2026-09-15T06:03:37.775146Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Local sequential optimization\n",
      "Depth =  1   Energy = -10.18106100   Error =  2.634e-02   1-F =  8.268e-02\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  2   Energy = -10.44499676   Error =  1.096e-03   1-F =  3.466e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  3   Energy = -10.45641720   Error =  3.868e-06   1-F =  6.445e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  4   Energy = -10.45643009   Error =  2.636e-06   1-F =  4.525e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  5   Energy = -10.45643662   Error =  2.011e-06   1-F =  3.416e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth =  6   Energy = -10.45643808   Error =  1.872e-06   1-F =  3.065e-06\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"*** Local sequential optimization\")\n",
    "# --- Modified: store errors vs depth ---\n",
    "depths_local = []\n",
    "errors_local = []\n",
    "\n",
    "circ = ansatz_circuit_su4(\n",
    "    n_qubits=n_qubits, depth=1, param_scaling=1.0, parametrize=False\n",
    ")\n",
    "\n",
    "tn = circ.psi\n",
    "tn_fit(tn, GS, tags=\"SU4SWAP\", steps=10000, tol=1e-8)\n",
    "\n",
    "ovlp = (dmrg.state.H & tn).contract()\n",
    "tnH = tn.H\n",
    "tn.align_(mpo, tnH)\n",
    "energy_tn = tnH & mpo & tn\n",
    "ene = autoray.do(\"real\", energy_tn.contract(all))\n",
    "err = np.abs(1 - ene / dmrg_energy)\n",
    "depths_local.append(1)\n",
    "errors_local.append(err)\n",
    "\n",
    "print(\n",
    "    f\"Depth = {1:2d}   Energy = {ene:12.8f}   Error = {err:10.3e}   1-F = {1 - np.abs(ovlp) ** 2:10.3e}\"\n",
    ")\n",
    "\n",
    "rng = np.random.default_rng(42)\n",
    "\n",
    "for ii in range(2, depth + 1):\n",
    "    # grow the optimized network by one brick-wall layer of SU4SWAP gates,\n",
    "    # initialized close to the identity (small parameters)\n",
    "    tags = [\"SU4SWAP\", f\"ROUND_{ii - 1}\"]\n",
    "    for start in range(2):\n",
    "        for q in range(start, n_qubits - 1, 2):\n",
    "            gate = su4swap_gate_param_gen(1e-2 * rng.random(15))\n",
    "            tn.gate_(gate, (q, q + 1), tags=tags, contract=False)\n",
    "\n",
    "    tn_fit(tn, GS, tags=\"SU4SWAP\", steps=10000, tol=1e-8)\n",
    "\n",
    "    ovlp = (dmrg.state.H & tn).contract()\n",
    "    tnH = tn.H\n",
    "    tn.align_(mpo, tnH)\n",
    "    energy_tn = tnH & mpo & tn\n",
    "    ene = autoray.do(\"real\", energy_tn.contract(all))\n",
    "    err = np.abs(1 - ene / dmrg_energy)\n",
    "    depths_local.append(ii)\n",
    "    errors_local.append(err)\n",
    "\n",
    "    print(\n",
    "        f\"Depth = {ii:2d}   Energy = {ene:12.8f}   Error = {err:10.3e}   1-F = {1 - np.abs(ovlp) ** 2:10.3e}\"\n",
    "    )\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a595d259",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "#### Plot: Error vs depth for local sequential optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e31a430f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-15T06:03:37.777273Z",
     "iopub.status.busy": "2026-09-15T06:03:37.777119Z",
     "iopub.status.idle": "2026-09-15T06:03:37.938332Z",
     "shell.execute_reply": "2026-09-15T06:03:37.937891Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAGGCAYAAABmPbWyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZCxJREFUeJzt3XlYVGX/BvD7sO8gsiuCu6ICikruouSSmWblUibuZbiFmppvqLlmbr2FkqZppmmumZU77rmBaG64IS4ooAgICgLz/P7wx7xOgA7DzJwB7s91zVVz5sw59zwzg995znOeIwkhBIiIiIjopYzkDkBERERUFrBoIiIiIlIDiyYiIiIiNbBoIiIiIlIDiyYiIiIiNbBoIiIiIlIDiyYiIiIiNbBoIiIiIlIDiyYiIiIiNbBoIq2SJAnTpk2TO0aZ4u3tjYEDB5b4eTdv3oQkSVi1apXWM2lKF5kOHDgASZJw4MABrW3TkPdbQKFQoGHDhpg1a5Ys+yfda9++PRo2bKjx89esWYN69erB1NQUDg4Oym22b9++VLn09bdl0qRJCAwM1Ok+tI1FUxm0atUqSJKE06dPyx2F1HTs2DFMmzYNaWlpckfRinXr1mHx4sVyx9CKJUuWGFThWeCXX37B7du3MXLkSOWy8vY50idDfZ81dfnyZQwcOBA1a9bE8uXLsWzZsiLXS0xMxLRp0xAbG1voMbm/x2PHjsXZs2exfft22TKUmKAy58cffxQAxKlTp+SOUggAMXXqVLljGJyvv/5aABDx8fGFHsvOzhbPnj0r8Tbj4+MFAPHjjz+WPmAJdevWTXh5eRVarlAoxNOnT0VeXp7W9pWfny+ePn0q8vPztbbNFzVo0EC0a9dO7/t9FT8/PzF8+HCVZS/7HNHLFfc+y6ldu3aiQYMGGj136dKlAoC4evWqyvKcnByRk5OjvH/q1Kli/07o83tcnN69e4s2bdrofD/awp4mIpmZm5vD1NRU7hhaIUkSLCwsYGxsrLVtGhkZwcLCAkZG+v1zJdd+AeDMmTM4e/Ysevfurfd960tWVpbcEWSjjdeenJwMAMrDcgXMzMxgZmZWqm3r4ntcnN69e+PIkSO4ceOGzvelFXJXbVRy6vY0xcTEiC5dughbW1thbW0tOnToIP7+++9C6z169EiMHTtWeHl5CTMzM1GlShXx4YcfipSUFCHE818uX3zxhWjSpImws7MTVlZWonXr1mL//v2FtgU1e5r++9//Ch8fH2FpaSkcHBxEQECAWLt2rco6d+7cEYMGDRIuLi7CzMxM+Pj4iBUrVhTa1u3bt0WPHj2ElZWVcHZ2FmPHjhU7d+4UAERUVJRyPS8vLxESElLo+e3atSv0CzQ7O1uEh4eLmjVrCjMzM1G1alUxYcIEkZ2dXej1hoaGiq1bt4oGDRooc/7111/KdaZOnSoAFLoV9Bb8O9fDhw/FuHHjRMOGDYW1tbWwtbUVXbp0EbGxsSr7LklP0/Xr18W7774rKlWqJCwtLUVgYKDYsWOHyjpRUVECgFi/fr2YPHmycHV1FVZWVqJ79+7i1q1bKu3179dS8Gu1qEwhISHC2tpaJCQkiG7duglra2vh4eEhvvvuOyGEEOfOnRNBQUHCyspKVKtWrdDnoCDXi+/llStXRK9evYSrq6swNzcXVapUEX369BFpaWnKdVauXCmCgoKEs7OzMDMzE/Xr1xdLlixR2baXl1eh11LwWShqv0II8euvv4omTZoICwsLUblyZfHBBx+IO3fuqKxT8Jrv3LkjevToIaytrYWTk5MYN26cWr/ew8PDhZmZmUoP5Ks+R0IIsWbNGmW2SpUqiT59+qi8d0L8r3fjwoULon379sLS0lJ4eHiIr776qlAOdb6n6vydKfibdeDAATFixAjh7OwsHBwcin39JWm//Px8sWjRIuHj4yPMzc2Fi4uLGD58uEhNTVWuU9z7/OjRI2FkZCS++eYb5bopKSlCkiTh6OgoFAqFcvnHH38sXF1dVfZdks/CtWvXRNeuXYWNjY3o0aOHynvxol27dglLS0vRt29fkZubW2T7FPV6Cv7uvvj3rOAz/O/bjz/+qPH3WJ335MGDB6J///7C1tZW2NvbiwEDBojY2Ngi/16lpaUJSZLEwoULi3ythsZEB3UYGYALFy6gTZs2sLOzw2effQZTU1N8//33aN++PQ4ePKgcfJeZmYk2bdrg0qVLGDx4MJo0aYIHDx5g+/btuHPnDpycnJCRkYEffvgB/fr1w7Bhw/D48WOsWLECnTt3xsmTJ+Hv71+ibMuXL8fo0aPx7rvvYsyYMcjOzsa5c+dw4sQJvP/++wCApKQkvPbaa5AkCSNHjoSzszP++usvDBkyBBkZGRg7diwA4OnTp+jYsSNu3bqF0aNHw8PDA2vWrMH+/fs1bjuFQoG33noLR44cwfDhw1G/fn38888/WLRoEa5cuYJt27aprH/kyBFs2bIFn3zyCWxtbfHf//4X77zzDm7duoXKlSujV69euHLlCn755RcsWrQITk5OAABnZ+ci93/jxg1s27YN7733HqpXr46kpCR8//33aNeuHS5evAgPD48SvZ6kpCS0bNkST548wejRo1G5cmWsXr0ab731FjZt2oS3335bZf1Zs2ZBkiRMnDgRycnJWLx4MYKDgxEbGwtLS0tMmTIF6enpuHPnDhYtWgQAsLGxeWmG/Px8dO3aFW3btsW8efOwdu1ajBw5EtbW1pgyZQo++OAD9OrVC5GRkRgwYABatGiB6tWrF7mtZ8+eoXPnzsjJycGoUaPg5uaGu3fvYseOHUhLS4O9vT0AYOnSpWjQoAHeeustmJiY4Pfff8cnn3wChUKB0NBQAMDixYsxatQo2NjYYMqUKQAAV1fXYl/HqlWrMGjQIDRr1gxz5sxBUlISvvnmGxw9ehRnzpxR+dWfn5+Pzp07IzAwEPPnz8fevXuxYMEC1KxZEyNGjHhpex07dgwNGzZU6YF81edo1qxZ+OKLL9C7d28MHToUKSkp+Pbbb9G2bdtC2R49eoQuXbqgV69e6N27NzZt2oSJEyeiUaNG6Nq1KwD1vqfq/p0p8Mknn8DZ2Rnh4eGv7G1Rt/0++ugj5fsyevRoxMfH47vvvsOZM2dw9OhRmJqaFvs+Ozg4oGHDhjh06BBGjx4N4Pn3WZIkpKam4uLFi2jQoAEA4PDhw2jTpo1Gn4W8vDx07twZrVu3xvz582FlZVXka96xYwfeffdd9OnTBytXriy2p2fx4sX46aefsHXrVixduhQ2Njbw9fUttF79+vXx5ZdfIjw8HMOHD1fmb9myJapUqaLR9/hV74lCoUD37t1x8uRJjBgxAvXq1cNvv/2GkJCQIrdpb2+PmjVr4ujRo/j0009fun+DIHfVRiWnTk9Tz549hZmZmbh+/bpyWWJiorC1tRVt27ZVLgsPDxcAxJYtWwpto+BXVl5ensoxciGe9065urqKwYMHqyyHGj1NPXr0eOVx/CFDhgh3d3fx4MEDleV9+/YV9vb24smTJ0IIIRYvXiwAiF9//VW5TlZWlqhVq5bGPU1r1qwRRkZG4vDhwyrrRUZGCgDi6NGjKq/XzMxMXLt2Tbns7NmzAoD49ttvlcteNhbl37mys7MLjaOJj48X5ubm4ssvv1RZBjV6msaOHSsAqLyex48fi+rVqwtvb2/lvgp+lVapUkVkZGQo1/31118FAJVf48WNhSjuFyoAMXv2bOWyR48eCUtLSyFJkli/fr1y+eXLlwt9hv7d43PmzBkBQGzcuPGlr7vgM/Kizp07ixo1aqgsK26sy7/3++zZM+Hi4iIaNmwonj59qlxvx44dAoAIDw8v9JpffL+EEKJx48YiICDgpbmFEKJq1arinXfeKbS8uM/RzZs3hbGxsZg1a5bK8n/++UeYmJioLC/oYfjpp5+Uy3JycoSbm5vKPtX5nqr7d6bgb1br1q3V6mlTt/0OHz4sABTq/SroaX5xeXHvc2hoqEoPUlhYmGjbtq1wcXERS5cuFUI87/2VJEn5HdDkszBp0qRC+36xp2nz5s3C1NRUDBs2TK1xdAU9jwVHBF7c5ouvU5MxTS/7Hr/qPdm8ebMAIBYvXqxclp+fLzp06FBsjk6dOon69eu//AUbCI5pKofy8/Oxe/du9OzZEzVq1FAud3d3x/vvv48jR44gIyMDALB582b4+fkV6m0Anh/XBgBjY2PlMXKFQoHU1FTk5eWhadOmiImJKXE+BwcH3LlzB6dOnSrycSEENm/ejO7du0MIgQcPHihvnTt3Rnp6unK/f/75J9zd3fHuu+8qn29lZYXhw4eXOFeBjRs3on79+qhXr57Kvjt06AAAiIqKUlk/ODgYNWvWVN739fWFnZ2dxsfozc3NleNo8vPz8fDhQ9jY2KBu3boatfeff/6J5s2bo3Xr1splNjY2GD58OG7evImLFy+qrD9gwADY2toq77/77rtwd3fHn3/+qdHrKTB06FDl/zs4OKBu3bqwtrZWGbdTt25dODg4vLTtCnqSdu3ahSdPnhS7nqWlpfL/09PT8eDBA7Rr1w43btxAenp6ifOfPn0aycnJ+OSTT2BhYaFc3q1bN9SrVw9//PFHoed8/PHHKvfbtGmj1ufi4cOHqFSpktrZtmzZAoVCgd69e6t8Zt3c3FC7du1Cn1kbGxv0799fed/MzAzNmzdXyfaq72lJ/s4UGDZsWInGybyq/TZu3Ah7e3u8/vrrKq87ICAANjY2hV53Udq0aYOkpCTExcUBeN6j1LZtW7Rp0waHDx8G8Lz3SQih7KnR5LPwst7FX375BX369MFHH32E77//XpZxdOp61Xuyc+dOmJqaYtiwYcplRkZGyt7dolSqVAkPHjzQflgdMNx3hjSWkpKCJ0+eoG7duoUeq1+/PhQKBW7fvg0AuH79ulrzhKxevRq+vr6wsLBA5cqV4ezsjD/++EOjf3wmTpwIGxsbNG/eHLVr10ZoaCiOHj2qkj8tLQ3Lli2Ds7Ozym3QoEEA/jcIMiEhAbVq1VIWeAWKeu3qunr1Ki5cuFBo33Xq1FHZd4Fq1aoV2kalSpXw6NEjjfavUCiwaNEi1K5dG+bm5nBycoKzszPOnTunUXsnJCQU+1koePxFtWvXVrkvSRJq1aqFmzdvlnjfBSwsLAodjrS3t0fVqlULvXf29vYvbbvq1asjLCwMP/zwA5ycnNC5c2dEREQUapujR48iODgY1tbWcHBwgLOzMz7//HMA0LgdgaI/W/Xq1SvUjkW95pJ8LoQQame7evUqhBCoXbt2oc/tpUuXCn1mi2r3f2dT53uq7t+ZAi8ecn327Bnu37+vcsvPz1c+rk77Xb16Fenp6XBxcSn0ujMzMwu97qIUFEKHDx9GVlYWzpw5gzZt2qBt27bKounw4cOws7ODn58fgJJ/FkxMTFC1atUi9x8fH4/+/fvjnXfewbffflvofTEk6rwnCQkJcHd3L3QIslatWsVuVwhh0K/7RRzTRK/0888/Y+DAgejZsycmTJgAFxcXGBsbY86cObh+/XqJt1e/fn3ExcVhx44d2LlzJzZv3owlS5YgPDwc06dPh0KhAAD079+/2OPgRR2/f5XivpT5+fkqv34VCgUaNWqEhQsXFrm+p6enyv3ifjmX5B+9F82ePRtffPEFBg8ejBkzZsDR0RFGRkYYO3assm3KmuLaSNO2W7BgAQYOHIjffvsNu3fvxujRozFnzhwcP34cVatWxfXr19GxY0fUq1cPCxcuhKenJ8zMzPDnn39i0aJFemnH0px5VLly5RIV3QqFApIk4a+//ipyv/8eq6JOu7/qe6qJF3v/jh07hqCgIJXH4+Pj4e3t/dKML1IoFHBxccHatWuLfLy4cYMv8vDwQPXq1XHo0CF4e3tDCIEWLVrA2dkZY8aMQUJCAg4fPoyWLVtq3AP0Yu/xv7m7uyt7ck+fPo2mTZtqtA990NXZdI8ePVKO0TN0LJrKIWdnZ1hZWSm7m190+fJlGBkZKf/hr1mzJs6fP//S7W3atAk1atTAli1bVAqPqVOnapzR2toaffr0QZ8+ffDs2TP06tULs2bNwuTJk+Hs7AxbW1vk5+cjODj4pdvx8vLC+fPnC/1SKeq1V6pUqchJARMSElQOL9SsWRNnz55Fx44dtfbrpyTb2bRpE4KCgrBixQqV5WlpaRr9YfHy8ir2s1Dw+IuuXr2qcl8IgWvXrqkUqobwq7BRo0Zo1KgR/vOf/+DYsWNo1aoVIiMjMXPmTPz+++/IycnB9u3bVXoCizpco+5rKWinuLg45aHaAnFxcYXasTTq1auH+Pj4QsuLy1qzZk0IIVC9enVlj6g2vOp7qu7fmaL4+flhz549Ksvc3NxKlK9mzZrYu3cvWrVqpVKQFeVl73ObNm1w6NAhVK9eHf7+/rC1tYWfnx/s7e2xc+dOxMTEqBSK2vwsWFhYYMeOHejQoQO6dOmCgwcPKgefa8PLXrcuvsdeXl6IiorCkydPVHqbrl27Vuxz4uPjlb14ho6H58ohY2NjdOrUCb/99pvKIZWkpCSsW7cOrVu3hp2dHQDgnXfewdmzZ7F169ZC2yn41Vnw6+LFX6EnTpzA33//rVG+hw8fqtw3MzODj48PhBDIzc2FsbEx3nnnHWzevLnIgi4lJUX5/2+88QYSExOxadMm5bInT54UOTtuzZo1cfz4cTx79ky5bMeOHYUOIfTu3Rt3797F8uXLC23j6dOnGs2xYm1tDQBqzeRsbGxcqKdl48aNuHv3bon3Czxvo5MnT6q8X1lZWVi2bBm8vb3h4+Ojsv5PP/2Ex48fK+9v2rQJ9+7dU55VBTx/PZoc4tKGjIwM5OXlqSxr1KgRjIyMkJOTA6Doz2x6ejp+/PHHQtuztrZW631p2rQpXFxcEBkZqdwPAPz111+4dOkSunXrpsnLKVKLFi1w/vx5lf0UZAUKf4569eoFY2NjTJ8+vdBnRwhR6DunDnW+p+r+nSlKpUqVEBwcrHJ7cXyQOnr37o38/HzMmDGj0GN5eXkq7fSy97lNmza4efMmNmzYoDxcZ2RkhJYtW2LhwoXIzc1VOXNO258Fe3t77Nq1Cy4uLnj99dc16sEvzsv+9ujie9y5c2fk5uaq/P1UKBSIiIgocv309HRcv34dLVu21GoOXWFPUxm2cuVK7Ny5s9DyMWPGYObMmdizZw9at26NTz75BCYmJvj++++Rk5ODefPmKdedMGECNm3ahPfeew+DBw9GQEAAUlNTsX37dkRGRsLPzw9vvvkmtmzZgrfffhvdunVDfHw8IiMj4ePjg8zMzBLn7tSpE9zc3NCqVSu4urri0qVL+O6779CtWzflAOS5c+ciKioKgYGBGDZsGHx8fJCamoqYmBjs3bsXqampAJ4PLP3uu+8wYMAAREdHw93dHWvWrCnylN6hQ4di06ZN6NKlC3r37o3r16/j559/VhnEDQAffvghfv31V3z88ceIiopCq1atkJ+fj8uXL+PXX3/Frl27StyFHhAQAACYMmUK+vbtC1NTU3Tv3l35B+1Fb775Jr788ksMGjQILVu2xD///IO1a9eq9IaVxKRJk/DLL7+ga9euGD16NBwdHbF69WrEx8dj8+bNhQ4bODo6onXr1hg0aBCSkpKwePFi1KpVS2VgZ0BAADZs2ICwsDA0a9YMNjY26N69u0b5Smr//v0YOXIk3nvvPdSpUwd5eXlYs2aNstgGnn/GzMzM0L17d3z00UfIzMzE8uXL4eLignv37qlsLyAgAEuXLsXMmTNRq1YtuLi4FOo9AABTU1N89dVXGDRoENq1a4d+/fopTzP39vbW6unSPXr0wIwZM3Dw4EF06tRJJStQ+HNUs2ZNzJw5E5MnT8bNmzfRs2dP2NraIj4+Hlu3bsXw4cMxfvz4EmVQ53uq7t8ZXWnXrh0++ugjzJkzB7GxsejUqRNMTU1x9epVbNy4Ed98843yJJGXvc8FBVFcXBxmz56t3H7btm3x119/wdzcHM2aNVMu18VnwcnJSdmWwcHBOHLkCKpUqVKa5gHw/Meig4MDIiMjYWtrC2trawQGBqJ69eo6+R737NkTzZs3x7hx43Dt2jXUq1cP27dvV/7N/nfv1t69eyGEQI8ePUq1X73R56l6pB0Fp+8Wd7t9+7YQ4vmkc507dxY2NjbCyspKBAUFiWPHjhXa3sOHD8XIkSNFlSpVlBM5hoSEKE/3VygUYvbs2cLLy0uYm5uLxo0bix07doiQkJBCp6tCjSkHvv/+e9G2bVtRuXJlYW5uLmrWrCkmTJgg0tPTVdZLSkoSoaGhwtPTU5iamgo3NzfRsWNHsWzZMpX1EhISxFtvvSWsrKyEk5OTGDNmTJGTWwohxIIFC0SVKlWEubm5aNWqlTh9+nSRk1s+e/ZMfPXVV6JBgwbC3NxcVKpUSQQEBIjp06er5MT/T275b0VNbzBjxgxRpUoVYWRk9NLJLbOzs8W4ceOEu7u7sLS0FK1atRJ///13oZyaTG7p4OAgLCwsRPPmzYud3PKXX34RkydPFi4uLsLS0lJ069ZNJCQkqKybmZkp3n//feHg4KD2pHj/VtwlJLy8vES3bt0K5Sp4L2/cuCEGDx4satasKSwsLISjo6MICgoSe/fuVdnO9u3bha+vr7CwsBDe3t7iq6++EitXrix0yv79+/dFt27dhK2trVqTW27YsEE0btxYmJubC0dHx5dOaPhvBaeJq8PX11cMGTKk0PLiPkdCPD/du3Xr1sLa2lpYW1uLevXqidDQUBEXF6dcp7h2//f3Wd3vqTp/Z0p66aeStt+yZctEQECAsLS0FLa2tqJRo0bis88+E4mJicp1inufC7i4uAgAIikpSbnsyJEjAkCxl/kozWdBiKLfi2vXrgl3d3dRv379QtMJFNUWr5pyQAghfvvtN+Hj4yNMTExUvp/a+B4X9Z6kpKSI999/Xzm55cCBA8XRo0cF/n/y3Bf16dNHtG7dutjXaWgkITQcrUpkwA4cOICgoCBERUWV+orfFUVBm23cuFFlCgeSx5o1axAaGopbt24VulQGUVmzbds2vP322zhy5AhatWoFALh//z6qV6+O9evXl5meJo5pIiIyQB988AGqVatW7FgQIkP19OlTlfv5+fn49ttvYWdnhyZNmiiXL168GI0aNSozBRPAMU1ERAbJyMjolWe2EhmiUaNG4enTp2jRogVycnKwZcsWHDt2DLNnz1Y5y3Hu3LkyptQMiyYiIiLSmg4dOmDBggXYsWMHsrOzUatWLXz77bcYOXKk3NFKjWOaiIiIiNTAMU1EREREamDRRERERKQGjmkqIYVCgcTERNja2hrEpSSIiIiodIQQePz4MTw8PF56jUEWTSWUmJj40uspERERUdl0+/ZtVK1atdjHWTSVUMHlA27fvv3S6yppQqFQICUlBc7OzhpfTZvUw7bWL7a3/rCt9YvtrT+6bOuMjAx4enoq/40vDosmNUVERCAiIgL5+fkAADs7O50UTdnZ2bCzs+OXT8fY1vrF9tYftrV+sb31Rx9t/aphN3yH1RQaGoqLFy/i1KlTckchIiIiGbBoIiIiIlIDiyYiIiIiNbBoIiIiIlIDiyYiIiIiNbBoIiIiIlIDpxyQ2a30W3jw5AGA56dTpqamwjHfUXk6pZOVE6rZV5MzIhEREYFFk9r+PU+TNtxKv4W639VFdl52setYmFggbmQcCyciIiKZ8fCcmnQxT9ODJw9eWjABQHZetrInioiIiOTDoomIiIhIDSyaiIiIiNTAoomIiIhIDSyaiIiIiNTAoomIiIhIDSyaiIiIiNTAoklGTlZOsDCxeOV6O67s0EMaIiIiehlObimjavbVEDcyrvCM4I6OkCQJ3538DqvOrsLUA1NR2bIyQpuHypyYiIio4mLRJLNq9tWUs30rFAokGyfDxcUFRkZGWNljJVxtXPHV0a8w8q+RMJKMMKLZCJkTExERVUw8PGfAJEnCnI5zMKHlBADAJ39+gu9Pfy9zKiIiooqJRZOaIiIi4OPjg2bNmul1v5Ik4avgrzCuxTgAwMd/fIxl0cv0moGIiIhYNKlNF9eeU5ckSfj69a8R9loYAOCjHR/hh5gf9J6DiIioImPRVEZIkoT5neZjbOBYAMCw34dhRcwKeUMRERFVICyayhBJkrCw80KMbj4awPPCaeWZlTKnIiIiqhhYNJUxkiRhcZfFGNV8FAQEhm4filWxq+SORUREVO6xaCqDJEnCN12+QWizUAgIDP5tMFbHrpY7FhERUbnGoqmMkiQJ33b9FiOajoCAwKDfBmHN2TVyxyIiIiq3WDSVYZIk4bs3vsPHAR9DQCBkWwh+Pvez3LGIiIjKJRZNZZyRZISIbhH4KOAjZeG07p91csciIiIqd1g0lQNGkhGWdFuCYU2GQSEU+HDrh/jln1/kjkVERFSusGgqJ4wkI0S+GYmhjYdCIRTov7U/NpzfIHcsIiKicoNFUzliJBnh++7fY7D/YCiEAh9s+QC/XvhV7lhERETlgoncAUi7jCQjLH9rORRQYFXsKry/+X1IkPBeg/fkjkZERFSmsaepHDKSjPBD9x8Q4heCfJGPfpv7YfPFzXLHIiIiKtMqXNF0+/ZttG/fHj4+PvD19cXGjRvljqQTxkbGWPHWCgzwG4B8kY++m/tiy6UtcsciIiIqsypc0WRiYoLFixfj4sWL2L17N8aOHYusrCy5Y+mEsZExVr61Ev19+yNPkYc+m/pg66WtcsciIiIqkypc0eTu7g5/f38AgJubG5ycnJCamipvKB0yNjLGqh6r8H6j95GnyEPvTb3x2+Xf5I5FRERU5hhc0XTo0CF0794dHh4ekCQJ27ZtK7ROREQEvL29YWFhgcDAQJw8eVKjfUVHRyM/Px+enp6lTG3YjI2MsbrnavRr2A95ijy8t/E9bI/bLncsIiKiMsXgiqasrCz4+fkhIiKiyMc3bNiAsLAwTJ06FTExMfDz80Pnzp2RnJysXMff3x8NGzYsdEtMTFSuk5qaigEDBmDZsmU6f02GwMTIBD+9/RP6NuyLXEUu3v31Xfwe97vcsYiIiMoMg5tyoGvXrujatWuxjy9cuBDDhg3DoEGDAACRkZH4448/sHLlSkyaNAkAEBsb+9J95OTkoGfPnpg0aRJatmz5ynVzcnKU9zMyMgAACoUCCoVCnZekNoVCASGE1rdbwAhGWN1jNfIV+dh4cSPe+fUdbHpvE96s86ZO9mfIdN3WpIrtrT9sa/1ie+uPLtta3W0aXNH0Ms+ePUN0dDQmT56sXGZkZITg4GD8/fffam1DCIGBAweiQ4cO+PDDD1+5/pw5czB9+vRCy1NSUpCdna1+eDUoFAqkp6dDCAEjI911Ai5stRDZ2dn4/cbveG/je1jZeSU6Vuuos/0ZIn21NT3H9tYftrV+sb31R5dt/fjxY7XWK1NF04MHD5Cfnw9XV1eV5a6urrh8+bJa2zh69Cg2bNgAX19f5XipNWvWoFGjRkWuP3nyZISFhSnvZ2RkwNPTE87OzrCzs9PshRRDoVBAkiQ4Ozvr/Mu3se9GfLD1A2y+tBmDdw/Glt5b0LVW8T185Y0+25rY3vrEttYvtrf+6LKtLSws1FqvTBVN2tC6desSde2Zm5vD3Ny80HIjIyOdfEEkSdLZtl9kbmSOX975RTl/0zu/voNtfbehS60uOt2vIdFXW9NzbG/9YVvrF9tbf3TV1upur0y9w05OTjA2NkZSUpLK8qSkJLi5uel03xEREfDx8UGzZs10uh99MjU2xfp31uPtem8jJz8HPdf3xO7ru+WORUREZJDKVNFkZmaGgIAA7Nu3T7lMoVBg3759aNGihU73HRoaiosXL+LUqVM63Y++mRqbYv2769GzXk/k5Oegx/oe2HN9j9yxiIiIDI7BFU2ZmZmIjY1VngEXHx+P2NhY3Lp1CwAQFhaG5cuXY/Xq1bh06RJGjBiBrKws5dl0VHJmxmbY8O4GvFX3LWTnZeOt9W9h7429csciIiIyKAZXNJ0+fRqNGzdG48aNATwvkho3bozw8HAAQJ8+fTB//nyEh4fD398fsbGx2LlzZ6HB4VQyZsZm2PjeRnSv0x3Zedno/kt37Lux79VPJCIiqiAMrmhq3749hBCFbqtWrVKuM3LkSCQkJCAnJwcnTpxAYGCgznOVxzFN/1ZQOL1Z501l4RQVHyV3LCIiIoNgcEWToSqvY5r+zdzEHJve24Q3ar+Bp3lP0W1dNxy4eUDuWERERLJj0USFmJuYY3Pvzehaq6uycDp486DcsYiIiGTFoomKZGFigS19tqBLrS54kvsEb6x7A4cSDskdi4iISDYsmtRUEcY0/ZuFiQW29tmKTjU7PS+c1r6BwwmH5Y5FREQkCxZNaqooY5r+zcLEAtv6bMPrNV5HVm4Wuq7tiiO3jsgdi4iISO9YNNErWZpa4re+vyG4RrCycDp2+5jcsYiIiPSKRROppaBw6li9IzKfZaLLz13w9+2/5Y5FRESkNyya1FQRxzT9m5WpFbb3244g7yA8fvYYnX/ujON3jssdi4iISC9YNKmpoo5p+jcrUyv83u93tPduryycTtw5IXcsIiIinWPRRCVmbWaNHf12oJ1XO2TkZKDTz51w8u5JuWMRERHpFIsm0oi1mTX+eP8PtPVq+7xwWtMJp+5W7F44IiIq31g0kcYKCqc21dogPScdr695HacTT8sdi4iISCdYNKmJA8GLZmNmgz8/+BOtq7VWFk4x92LkjkVERKR1LJrUxIHgxbMxs8Gf7/+Jlp4tkZadhuCfglk4ERFRuWNSkpWrV68OSZJKvJOxY8di9OjRJX4elR225rb464O/ns/fdOdvBP8UjH0D9qGxe2O5oxEREWlFiYqmVatWabQTb29vjZ5HZYuduR129t+pnL8peE0w9g/YDz83P7mjERERlVqJiqZ27drpKgeVE3bmdtj5wfPC6cTdE+j4U0fsG7CPhRMREZV5pR7TpFAotJGDyhF7C3vs6r8Lzas0x8OnD9Hxp444l3RO7lhERESlUuqiacKECRg/frzy/uDBgzFt2jT89ttvyMzMLO3mqYwqKJyaeTRTFk7nk8/LHYuIiEhjpS6a9uzZg1mzZinvHz9+HJUrV8bPP/+M2bNnl3bzBoNTDpScg4UDdn+4G009muLBkwfosLoDCyciIiqzSl00mZubw9zcXHnf19cXo0aNwi+//IIDBw6UdvMGg1MOaMbBwgG7++9GE/cmSHmSgg6rO+BC8gW5YxEREZVYqYsmMzMzJCUlKe+vX78eAGBiYoJnz56VdvNUDlSyrIQ9H+5BY7fGzwunnzrgYspFuWMRERGViFbGNL399tu4ffu2yvIHDx4gJyentJuncsLR0hF7B+yFv5s/krOS0WF1B1xKuSR3LCIiIrWVaMqBovTs2RMZGRlo2rQpWrZsiYYNG0KSJPz666+YNGmSNjJSOeFo6Yi9H+5Fx5864mzSWQStDsKBgQdQz6me3NGIiIheSSuXURkwYACuXLmC9957D0IImJub4+eff8bAgQO1sXkqRypbVcbeAXvh6+qLpKwkBK0OQtyDOLljERERvVKpe5oK2Nvb4/3339fW5qgcc7Jywr4B+9BhdQf8k/yPssepTuU6ckcjIiIqVql7mvbu3YvXXnsN7dq1w969ewEA9+7dw48//oh+/fqVOiCVTwWFU0OXhriXeQ9Bq4Nw9eFVuWMREREVq9RF08iRIzFu3DgsXLgQW7duxaBBg1C3bl3s3bsXb731ljYyGgTO06R9ztbO2DdgHxo4N0Di40QErQ7CtdRrcsciIiIqkiSEEKXZQOPGjXHmzBkAQH5+PpydnXHu3DlUrVpVKwENTUZGBuzt7ZGeng47OzutbluhUCA5ORkuLi4wMtLKcLMyITkrGUGrg3Ax5SKq2lXFgZADqOlYU6f7rKhtLRe2t/6wrfWL7a0/umxrdf9tL/VeU1JSsGHDBsTExODp06fw9vYutwUT6YaLtQv2D9iP+k71cSfjDoJWB+HGoxtyxyIiIlJR6qIpLCwMu3fvxogRI+Dh4YELFy6gR48emDp1KrZu3aqNjFQBuNq4Yn/IftRzqofbGbfRflV7Fk5ERGRQSl00DRkyBCtWrMCJEyeQkZGBy5cvY+jQobCwsMCmTZu0kZEqCDcbN+wfsB91K9fF7YzbCFodhJtpN+WORUREBEALUw5UqlQJnp6eaNCgARo2bIhGjRqhYcOG6NSpk8o16YjU4W7rjqiQKLRf3R5XHl5B+1XtcWDgAXg7eMsdjYiIKrhS9zSNGDECTk5OaNWqFerWrYvo6GhMmDAB1atXR716nOmZSq6gcKrtWBsJ6QkIWh2EhLQEuWMREVEFV+qiKSIiAtu2bUNCQgLWr1+PXr16Ye/evUhMTMThw4e1kZEqIA9bD2XhdDPtJoJWB+FW+i25YxERUQWmlXP2PD09sWzZMixbtgyrV69G165dcfHiRTg7O2tj81RBVbGrgqiQKNRyrIX4tHgErQ7C7fTbr34iERGRDpR6TNOlS5cQFxeHuLg4XLp0CdevX0dWVhYuXLgAHx8fbWSkCqygcGq/qj2uP7qO9qvb4+DAg6hqx2ktiIhIv0pdNDVo0AC+vr7o3bs3wsLCUL9+fZiammojGxEAoKpdVeXg8BuPbqD9queFUxW7KnJHIyKiCqTUh+fmz5+PgIAA/Pbbb3j99dfRsmVLDBo0CAsWLMCuXbu0kZEInvaeiAqJQnWH6rj+6DqCVgfhbsZduWMREVEFUuqeprCwMJX78fHxOH/+PM6fP481a9agc+fOpd2FQYiIiEBERATy8/PljlJhVbOvpuxxupp6FUGrg3Bg4AF42HrIHY2IiCqAUl977kUKhaLcX3uH156T3820m2i/qj0S0hNQp3IdHAg5AHdb9xJtg22tX2xv/WFb6xfbW3/KxbXnHj16hPfffx/29vYwNzdHrVq1MGnSJKSlpZV200RF8nbwxoGBB1DNvhquPLyCoNVBuJ95X+5YRERUzpW6aJo0aRKqVauGO3fuICsrC3/++ScAoEWLFrh/n/+QkW54O3jjQMgBeNp5Iu5hHAsnIiLSuVIXTSdOnMDcuXNha2sLMzMz1KlTB3PnzkV4eDimTZumhYhERateqToODDyAqnZVcfnBZXRY3QFJmUlyxyIionKq1EVTcccV+/Xrh9OnT5d280QvVaNSDRwIeV44XXpwCR1+6oDkrGS5YxERUTlU6qIpOTkZmzZtwqVLlwqdWSZJUmk3T/RKNR1rIiokClVsq+BiykV0WM3CiYiItK/URdO4cePw119/YeDAgXB2dkbDhg3Rp08fTJ8+HcnJ/IeL9KOWYy1EhUTBw9YDF1IuoONPHZGSlSJ3LCIiKkdKPU/Tp59+qnL/xXmaWrduXdrNE6mtduXaykuunE8+j44/dcT+kP1wsnKSOxoREZUDJe5pGjhwIJ48eVLs49WrV0f37t0xefJkrF27tlThiEqqTuU6iAqJgruNO/5J/gcdf+qIB08eyB2LiIjKgRIXTWvWrEFmZqby/ogRIwrNyZSXl1fqYESaqutUF/tD9sPNxg3nks4h+KdgPHzyUO5YRERUxpW4aPr3BOJr165Famqq8n5SUpLWZ8omKql6TvWwf8B+uFq74mzSWQSvCUbq09RXP5GIiKgYpR7TVNRVWLKzs0u7WaJSq+9cH/tD9iNodRBi78ei7Y9t8W3Xb2FvYQ+FQoHU1FQ45jsqp81wsnJCNftqMqcmIiJDVeqiqSicaoAMhY+zD6JCotBmZRtcSLmADj91KHZdCxMLxI2MY+FERERF0mjKgXXr1iEmJga5ubnazqNzaWlpaNq0Kfz9/dGwYUMsX75c7kikYz7OPlj65tJXrpedl81B40REVKwS9zS1adMGU6dOxePHj2Fqaoq8vDxMnToVrVq1gr+/P5ydnXWRU2tsbW1x6NAhWFlZISsrCw0bNkSvXr1QuXJluaORDtVyrCV3BCIiKuNKXDQdPHgQAHDlyhVER0fjzJkziImJweeff460tDSDPzRnbGwMKysrAEBOTg6EEEWOyyIiIiJ6UYkPz4WHhyM6Ohp16tRBv379MG/ePOzduxepqam4fv061q9fj4kTJ2oc6NChQ+jevTs8PDwgSRK2bdtWaJ2IiAh4e3vDwsICgYGBOHnyZIn2kZaWBj8/P1StWhUTJkyAkxMnPyQiIqKXK3HRdOfOHXTt2hVVq1bFiBEjsHPnTjx79gzA84kt33vvPcyePVvjQFlZWfDz80NERESRj2/YsAFhYWGYOnUqYmJi4Ofnh86dO6tcsqVgvNK/b4mJiQAABwcHnD17FvHx8Vi3bh2SkpI0zktEREQVQ4kPz61cuRIKhQJHjx7F77//jjFjxuDevXt4/fXX0aNHD7z55ptwdHTUOFDXrl3RtWvXYh9fuHAhhg0bhkGDBgEAIiMj8ccff2DlypWYNGkSACA2Nlatfbm6usLPzw+HDx/Gu+++W+Q6OTk5yMnJUd7PyMgAACgUCigUCrX2oy6FQgEhhNa3S1C7TXXxvhI/2/rEttYvtrf+6LKt1d2mRlMOGBkZoU2bNmjTpg3mzZuHS5cu4ffff8f333+P4cOHo3nz5njrrbfQr18/VKlSRZNdFOnZs2eIjo7G5MmTVbIEBwfj77//VmsbSUlJsLKygq2tLdLT03Ho0CGMGDGi2PXnzJmD6dOnF1qekpKi9fmoFAoF0tPTIYRQzh1E2vHiBKyvWi/ZmBea1jZ+tvWHba1fbG/90WVbP378WK31tDJPU/369VG/fn189tlnSElJwfbt27F9+3YAwPjx47WxCwDAgwcPkJ+fD1dXV5Xlrq6uuHz5slrbSEhIwPDhw5UDwEeNGoVGjRoVu/7kyZMRFhamvJ+RkQFPT084OztrfeZzhUIBSZLg7OzML5+W1TGvAwtjC2TnF1/oWhhboE7VOnCxd9FjsoqBn239YVvrF9tbf3TZ1hYWFmqtV+KiaeDAgViyZInyDLR/c3Z2xpAhQzBkyJCSblovmjdvrvbhOwAwNzeHubl5oeVGRkY6+YJIkqSzbVdk3pW8ETcqTjkPU8GM4I+NH6PPpj7IF/lY2WMlvCt5yxu0HONnW3/Y1vrF9tYfXbW1utsrUxfsdXJygrGxcaGB20lJSXBzc9PJPgtERETAx8cHzZo10+l+SHeq2VdDE/cmypuvsy/ervc2Pmn2CQBg3rF5yFfky5ySiIgMVZm6YK+ZmRkCAgKwb98+5TKFQoF9+/ahRYsWOtlngdDQUFy8eBGnTp3S6X5I/8LbhcPe3B6x92Ox5twaueMQEZGBKnX/lrYv2JuZmYnY2FjlIbT4+HjExsbi1q1bAICwsDAsX74cq1evxqVLlzBixAhkZWUpz6YjKiknKyf8p+1/AACf7/scWc+yZE5ERESGSCcHYEszK/jp06fRuHFjNG7cGMDzIqlx48YIDw8HAPTp0wfz589HeHg4/P39ERsbi507dxYaHK5tPDxXvo1qPgrVHarjXuY9zD82X+44RERkgAzugr3t27dXntn24m3VqlXKdUaOHImEhATk5OTgxIkTCAwM1HqOf+PhufLN3MQcXwV/BeD52KbEx4kyJyIiIkNT4qKp4IK9TZs2hY2NDZ48eYKpU6ciMjISx48fVxkkTlSWvOvzLlp6tsST3Cf4Yv8XcschIiIDo/EFe69evYro6GjExMSUqQv2EhVHkiQs6LQALVa0wI+xP2J04Gj4ufnJHYuIiAyExpNb1q5dG7Vr10bfvn2Vy+Lj43H69GmcOXNGK+EMSUREBCIiIpCfz1PSy7PXqr6GPg36YMOFDRi3exz2fLiHPwSIiAgAIImiTn+jYmVkZMDe3h7p6ek6mRE8OTkZLi4unCRNx17W1jfTbqLed/WQk5+DHf12oFudbjKlLD/42dYftrV+sb31R5dtre6/7XyHif7F28EbYwLHAAAm7JmAPIVuJmslIqKyhUUTURE+b/M5nKyccOnBJSyPXi53HCIiMgAsmoiKYG9hj2ntpgEAph6YivTsdHkDERGR7LRSNF25ckVn15szFJzcsuIZHjAcdSvXRcqTFMw5MkfuOEREJDOtFE3169fHjRs3tLEpg8XJLSseU2NTfP361wCAxccX42baTXkDERGRrLRSNPEEPCqv3qzzJjpU74Cc/Bx8vu9zueMQEZGMOKaJ6CUKJryUIOGX87/g5N2TckciIiKZsGgiegV/N3+E+IcAAMJ2hbFnlYiogmLRpCYOBK/YZgbNhJWpFY7ePootl7bIHYeIiGTAoklNHAhesVWxq4LxLcYDAD7b+xly8nJkTkRERPrGoolITRNaTYC7jTtuPLqBiFMRcschIiI9Y9FEpCYbMxvM7DATADDj0Aw8fPJQ5kRERKRPWimaJk6ciMqVK2tjU0QGLcQvBL6uvkjLTsOMQzPkjkNERHqklaJpzpw5LJqoQjA2MsaCTgsAABGnInDl4RWZExERkb7w8JyaePYcFQiuEYw3ar+BPEUeJu6dKHccIiLSExZNauLZc/Sir1//GsaSMbZd3oaDNw/KHYeIiPSARRORBnycfTA8YDgAYNzucVAIhcyJiIhI17ReNN2+fRuDBw/W9maJDM609tNga2aL6HvRWPfPOrnjEBGRjmm9aEpNTcXq1au1vVkig+Ni7YLP2zy/iO/kfZPxJPeJzImIiEiXTEr6hO3bt7/08Rs3bmgchqisGfvaWCw9vRS30m9h0d+LMKXtFLkjERGRjpS4aOrZsyckSXrpRUslSSpVKKKywsLEAnM6zsEHWz7A3KNzMaTJELjZuMkdi4iIdKDEh+fc3d2xZcsWKBSKIm8xMTG6yElksPo27IvmVZoj81kmpkZNlTsOERHpSImLpoCAAERHRxf7+Kt6oYjKGyPJCAs7LQQA/HDmB5xPPi9zIiIi0oUSF00TJkxAy5Yti328Vq1aiIqKKlUoQ8TJLellWlVrhXfqvwOFUGDCnglyxyEiIh0ocdHUpk0bdOnSpdjHra2t0a5du1KFMkSc3JJe5avgr2BqZIqd13Zi17VdcschIiIt4+SWRFpS07EmRjYfCQAYv2c88hX5MiciIiJtYtFEpEVftP0CjpaOOJ98HivPrJQ7DhERaVGJphyoXr26RtMJjB07FqNHjy7x84jKmkqWlRDeNhxjd43FF1FfoG/DvrA1t5U7FhERaUGJiqZVq1ZptBNvb2+NnkdUFo1oNgLfnfoO11KvYd7ReZjRYYbckYiISAtKVDSVxwHeRNpmZmyGecHz0OvXXljw9wIMDxgOT3tPuWMREVEpcUwTkQ70rNcTbaq1wdO8p5iyn5dWISIqD1g0EemAJElY2Pn5hJdrzq1BdGLxE8ISEVHZwKKJSEeaejRFf9/+AIBxu8dxpnwiojKORRORDs3uMBsWJhY4mHAQ2+O2yx2HiIhKQaOiKTc3Fx07dsTVq1e1nYeoXPG090TYa2EAgAl7JuBZ/jOZExERkaY0KppMTU1x7tw5bWcxaLz2HGlqUutJcLF2wdXUq4g8HSl3HCIi0pDGh+f69++PFStWaDOLQeO150hTtua2+LL9lwCA6Qen49HTRzInIiIiTZRonqYX5eXlYeXKldi7dy8CAgJgbW2t8vjChQtLHY6ovBjSZAi+PfktLqRcwKzDszC/03y5IxERUQlpXDSdP38eTZo0AQBcuXJF5TFNLrVCVJ6ZGJlgfqf56Lq2K749+S0+afYJalSqIXcsIiIqAY2LpqioKG3mICr3utTqgk41O2H39d2YtHcSfn3vV7kjERFRCWhcNAFAWloaVqxYgUuXLgEAGjRogMGDB8Pe3l4r4YjKm/mvz4f/DX9svLgRR28dRatqreSOREREatJ4IPjp06dRs2ZNLFq0CKmpqUhNTcXChQtRs2ZNxMTEaDMjUbnRyLURBvsPBsAJL4mIyhqNi6ZPP/0Ub731Fm7evIktW7Zgy5YtiI+Px5tvvomxY8dqMSJR+TKjwwxYm1rjxN0T2HBhg9xxiIhITaXqaZo4cSJMTP53hM/ExASfffYZTp8+rZVwROWRm40bJrWeBACYtHcSsvOyZU5ERETq0LhosrOzw61btwotv337NmxtbUsViqi8C2sRhiq2VZCQnoD/nviv3HGIiEgNGhdNffr0wZAhQ7Bhwwbcvn0bt2/fxvr16zF06FD069dPmxmJyh0rUyvM7jgbADDr8CykZKXInIiIiF5F47Pn5s+fD0mSMGDAAOTl5QF4fnmVESNGYO7cuVoLSFRe9fftj29OfIOYezGYdmAaIrpFyB2JiIheQuOeJjMzM3zzzTd49OgRYmNjERsbi9TUVCxatAjm5ubazEhULhlJRljQaQEA4Pvo73Ep5ZLMiYiI6GU0Kppyc3PRsWNHXL16FVZWVmjUqBEaNWoEKysrbefTmSdPnsDLywvjx4+XOwpVYO2926NH3R7IF/n4bO9ncschIqKX0KhoMjU1xblz57SdRa9mzZqF1157Te4YRJj3+jyYGJlgx5Ud2Hdjn9xxiIioGBofnuvfvz9WrFihzSx6c/XqVVy+fBldu3aVOwoR6lSugxFNRwB4PuFlviJf5kRERFQUjYumvLw8LF26FE2bNsVHH32EsLAwlZumDh06hO7du8PDwwOSJGHbtm2F1omIiIC3tzcsLCwQGBiIkydPlmgf48ePx5w5czTOSKRt4e3CYW9uj7NJZ/HT2Z/kjkNEREXQ+Oy58+fPo0mTJgCAK1euqDwmSZLGgbKysuDn54fBgwejV69ehR7fsGEDwsLCEBkZicDAQCxevBidO3dGXFwcXFxcAAD+/v7KM/petHv3bpw6dQp16tRBnTp1cOzYMY1zEmmTk5UT/tP2P5iwZwKm7J+C3g16w9rMWu5YRET0Ao2LpqioKG3mUOratetLD5stXLgQw4YNw6BBgwAAkZGR+OOPP7By5UpMmvR8luXY2Nhin3/8+HGsX78eGzduRGZmJnJzc2FnZ4fw8HCtvg6ikhrVfBSWnFqC+LR4zD82H1PbT5U7EhERvUCjoik3NxddunRBZGQkateure1MxXr27Bmio6MxefJk5TIjIyMEBwfj77//Vmsbc+bMUR6aW7VqFc6fP//SgiknJwc5OTnK+xkZGQAAhUIBhUKhycsolkKhgBBC69ulwgyxrU2NTDGn4xz03dwX847Nw5DGQ+Bh6yF3LK0wxPYur9jW+sX21h9dtrW629SoaJLr7LkHDx4gPz8frq6uKstdXV1x+fJlnexzzpw5mD59eqHlKSkpyM7W7jXDFAoF0tPTIYSAkZHGw81IDYba1m0rt0Uz12Y4lXQKn/31GRa2Xyh3JK0w1PYuj9jW+sX21h9dtvXjx4/VWk/jw3MFZ8+V5dm/Bw4c+Mp1Jk+erDKwPSMjA56ennB2doadnZ1W8ygUCkiSBGdnZ375dMyQ23rxG4vR6sdWWB+3HuPbjoe/m7/ckUrNkNu7vGFb6xfbW3902dYWFhZqradx0ZSXl4eVK1di7969CAgIgLW16qDVhQu1/wvZyckJxsbGSEpKUlmelJQENzc3re8PAMzNzYuc4dzIyEgnXxBJknS2bVJlqG3dslpL9GnQBxsubMBnez/Dng/3lOrkCkNhqO1dHrGt9YvtrT+6amt1t6fxXgvOnrO1tcWVK1dw5swZ5e1lA7FLw8zMDAEBAdi3738TACoUCuzbtw8tWrTQyT4LREREwMfHB82aNdPpfogAYG7wXJgbm2Nf/D78efVPueMQEREM8Oy5zMxMXLt2TXk/Pj4esbGxcHR0RLVq1RAWFoaQkBA0bdoUzZs3x+LFi5GVlaU8m05XQkNDERoaioyMDNjb2+t0X0TeDt4YEzgG847Nw4Q9E9C5VmeYGGn8dSUiIi0oVf/W4cOH0b9/f7Rs2RJ3794FAKxZswZHjhzReJunT59G48aN0bhxYwBAWFgYGjdurDzDrU+fPpg/fz7Cw8Ph7++P2NhY7Ny5s9DgcKKy7vM2n8PJygmXHlzC8ujlcschIqrwNC6aNm/ejM6dO8PS0hIxMTHK0/LT09Mxe/ZsjQO1b98eQohCt1WrVinXGTlyJBISEpCTk4MTJ04gMDBQ4/2pi4fnSN/sLewxrd00AED4gXCkZ6fLG4iIqILTuGiaOXMmIiMjsXz5cpiamiqXt2rVCjExMVoJZ0hCQ0Nx8eJFnDp1Su4oVIEMDxiOupXr4sGTB5hzhJf+ISKSk8ZFU1xcHNq2bVtoub29PdLS0kqTiYj+n6mxKeZ3mg8AWHx8MW6m3ZQ3EBFRBaZx0eTm5qYyYLvAkSNHUKNGjVKFIqL/6Va7GzpU74Cc/Bx8vu9zueMQEVVYGhdNw4YNw5gxY3DixAlIkoTExESsXbsW48ePx4gRI7SZ0SBwTBPJRZIkLOi0ABIk/HL+F5y4c0LuSEREFZLGRdOkSZPw/vvvo2PHjsjMzETbtm0xdOhQfPTRRxg1apQ2MxoEjmkiOfm7+SPEPwQAELY7DEIImRMREVU8GhdNkiRhypQpSE1Nxfnz53H8+HGkpKRgxowZ2sxHRP9vZtBMWJla4djtY9h8abPccYiIKpxSz0NuZmYGHx8fNG/eHDY2NtrIRERFqGJXBeNbjAcATNw7ETl5OTInIiKqWHihHDVxTBMZggmtJsDdxh03Ht1AxKkIueMQEVUoLJrUxDFNZAhszGwws8NMAMCMQzPw8MlDmRMREVUcLJqIypgQvxD4uvoiLTsNXx78Uu44REQVBosmojLG2MgYCzotAAAsOb0EVx5ekTkREVHFoHHRFBISgkOHDmkzCxGpKbhGMN6o/QbyFHmYuHei3HGIiCoEjYum9PR0BAcHo3bt2pg9ezbu3r2rzVxE9Apfv/41jCVjbLu8DQdvHpQ7DhFRuadx0bRt2zbcvXsXI0aMwIYNG+Dt7Y2uXbti06ZNyM3N1WZGg8Cz58jQ+Dj7YHjAcADAuN3joBAKmRMREZVvpRrT5OzsjLCwMJw9exYnTpxArVq18OGHH8LDwwOffvoprl69qq2csuPZc2SIprWfBlszW0Tfi8a6f9bJHYeIqFzTykDwe/fuYc+ePdizZw+MjY3xxhtv4J9//oGPjw8WLVqkjV0QURFcrF3weZvnF/GdvG8ynuQ+kTkREVH5pXHRlJubi82bN+PNN9+El5cXNm7ciLFjxyIxMRGrV6/G3r178euvv+LLL3lKNJEujX1tLKrZV8OdjDtY9Dd/pBAR6YrGRZO7uzuGDRsGLy8vnDx5EqdPn8bHH38MOzs75TpBQUFwcHDQRk4iKoaFiQXmdpwLAJh7dC7uZ96XORERUfmkcdG0aNEiJCYmIiIiAv7+/kWu4+DggPj4eE13QURq6tuwL5pXaY7MZ5mYGjVV7jhEROWSxkXThx9+CAsLC21mMWg8e44MmSRJWNhpIQDghzM/4HzyeZkTERGVPyaaPjEsLKzI5ZIkwcLCArVq1UKPHj3g6OiocThDEhoaitDQUGRkZMDe3l7uOESFtKrWCu/UfwebL23G+N3jsbP/TrkjERGVKxoXTWfOnEFMTAzy8/NRt25dAMCVK1dgbGyMevXqYcmSJRg3bhyOHDkCHx8frQUmouJ9FfwVtsdtx67ru7Dr2i50rtVZ7khEROWGxofnevTogeDgYCQmJiI6OhrR0dG4c+cOXn/9dfTr1w93795F27Zt8emnn2ozLxG9RE3HmhjZfCQAYPye8chX5MuciIio/NC4aPr6668xY8YMlbPl7O3tMW3aNMybNw9WVlYIDw9HdHS0VoISkXq+aPsFHC0dcT75PFaeWSl3HCKicqNU155LTk4utDwlJQUZGRkAnp899+zZM83TEVGJVbKshPC24QCAL6K+wOOcxzInIiIqH0p1eG7w4MHYunUr7ty5gzt37mDr1q0YMmQIevbsCQA4efIk6tSpo62sRKSmEc1GoJZjLSRlJeGro1/JHYeIqFzQuGj6/vvv0bFjR/Tt2xdeXl7w8vJC37590bFjR0RGRgIA6tWrhx9++EFrYYlIPWbGZpgXPA8AsODvBbidflvmREREZZ9GRVNubi569OiBzz77DA8fPsSZM2dw5swZPHz4EMuWLYO1tTUAwN/fv9iJL4lIt3rW64k21dogOy8bU/ZPkTsOEVGZp1HRZGpqinPnzgEAbGxs4OvrC19fX9jY2Gg1nCHh5JZU1kiShIWdn094uebcGkQn8qQMIqLS0PjwXP/+/bFixQptZjFooaGhuHjxIk6dOiV3FCK1NfVoiv6+/QEA43aPgxBC5kRERGWXxpNb5uXlYeXKldi7dy8CAgKUh+QKLFy4sNThiKj0ZneYjU0XN+FgwkH8FvcbetbrKXckIqIySeOi6fz582jSpAmA5zOBv0iSpNKlIiKt8bT3RNhrYZh9ZDY+2/MZ3qj9BsyMzeSORURU5mhcNEVFRWkzBxHp0KTWk/DDmR9wNfUqIk9HYnTgaLkjERGVORqPaSKissPW3BYzgmYAAKYfnI5HTx/JnIiIqOwpVdF0+PBh9O/fHy1atMDdu3cBAGvWrMGRI0e0Eo6ItGdw48Fo4NwAqU9TMevwLLnjEBGVORoXTZs3b0bnzp1haWmJM2fOICcnB8Dzy6vMnj1bawGJSDtMjEwwv9N8AMB/T/wX11Ovy5yIiKhs0bhomjlzJiIjI7F8+XKYmpoql7dq1QoxMTFaCUdE2tWlVhd0qtkJuYpcTNo3Se44RERlisZFU1xcHNq2bVtoub29PdLS0kqTiYh0aP7r82EkGWHTxU04euuo3HGIiMoMjYsmNzc3XLt2rdDyI0eOoEaNGqUKRUS608i1EQb7DwbACS+JiEpC46Jp2LBhGDNmDE6cOAFJkpCYmIi1a9di/PjxGDFihDYzGgReRoXKkxkdZsDa1Bon7p7Ahgsb5I5DRFQmaFw0TZo0Ce+//z46duyIzMxMtG3bFkOHDsVHH32EUaNGaTOjQeBlVKg8cbNxw6TWz8c0Tdo7Cdl52TInIiIyfBoXTZIkYcqUKUhNTcX58+dx/PhxpKSkYMaMGdrMR0Q6EtYiDFVsqyAhPQHfHP9G7jhERAav1JNbmpmZwcfHB82bN4eNjY02MhGRHliZWmF2x+fTg8w+MhspWSkyJyIiMmwaX0YFAPbt24d9+/YhOTkZCoVC5bGVK1eWKhgR6V5/3/745sQ3iLkXg2kHpiGiW4TckYiIDJbGPU3Tp09Hp06dsG/fPjx48ACPHj1SuRGR4TOSjLCg0wIAwPfR3+NSyiWZExERGS6Ne5oiIyOxatUqfPjhh9rMQ0R61t67PXrU7YHf4n7DZ3s/w+/9fpc7EhGRQdK4p+nZs2do2bKlNrMQkUzmvT4PJkYm2HFlB/bd2Cd3HCIig6Rx0TR06FCsW7dOm1mISCZ1KtfBiKbP51cbt3sc8hX5MiciIjI8Gh+ey87OxrJly7B37174+vqqXH8OABYuXFjqcESkP+HtwvHT2Z9wNuksfjr7EwY1HiR3JCIig6Jx0XTu3Dn4+/sDAM6fP6+tPEQkEycrJ3zR9guM3zMeU/ZPQe8GvWFtZi13LCIig6Fx0RQVFaXNHERkAEY2H4klp5fgxqMbmH9sPqa2nyp3JCIig1HiMU1vvPEG0tPTlffnzp2LtLQ05f2HDx/Cx8dHK+GISL/MTcwxt+NcAMC8Y/OQ+DhR5kRERIajxEXTrl27kJOTo7w/e/ZspKamKu/n5eUhLi5OO+mISO/e9XkXLT1b4knuE/xn/3/kjkNEZDBKXDQJIV56vyzw9vaGr68v/P39ERQUJHccIoMiSZJywstVsasQez9W3kBERAai1NeeK6uOHTuG2NhYjs0iKsJrVV9DnwZ9ICAwfvf4MvnjiIhI20pcNEmSBEmSCi0jovJlbvBcmBubY1/8Pvx59U+54xARyU6jw3MDBw5Er1690KtXL2RnZ+Pjjz9W3h88eHCpAh06dAjdu3eHh4cHJEnCtm3bCq0TEREBb29vWFhYIDAwECdPnizRPiRJQrt27dCsWTOsXbu2VHmJyitvB2+MCRwDABi/Zzxy83NlTkREJK8STzkQEhKicr9///6F1hkwYIDGgbKysuDn54fBgwejV69ehR7fsGEDwsLCEBkZicDAQCxevBidO3dGXFwcXFxcAAD+/v7Iy8sr9Nzdu3fDw8MDR44cQZUqVXDv3j0EBwejUaNG8PX11TgzUXn1eZvPsTJ2JS4/uIzlMcvxSbNP5I5ERCQbSRjwYAVJkrB161b07NlTuSwwMBDNmjXDd999BwBQKBTw9PTEqFGjMGnSpBLvY8KECWjQoAEGDhxY5OM5OTkqZwtmZGTA09MTjx49gp2dXYn39zIKhQIpKSlwdnaGkVGFHW6mF2xr9S05tQSjdo6Ck5UTroRegb2FfYm3wfbWH7a1frG99UeXbZ2RkYFKlSohPT39pf+2azy5pRyePXuG6OhoTJ48WbnMyMgIwcHB+Pvvv9XaRlZWFhQKBWxtbZGZmYn9+/ejd+/exa4/Z84cTJ8+vdDylJQUZGdnl/xFvIRCoUB6ejqEEPzy6RjbWn09qvbAYofFuJ52HeF7wjElcEqJt8H21h+2tX6xvfVHl239+PFjtdYrU0XTgwcPkJ+fD1dXV5Xlrq6uuHz5slrbSEpKwttvvw0AyM/Px7Bhw9CsWbNi1588eTLCwsKU9wt6mpydnXXS0yRJEn+x6AHbumQWdl6IHht6YPk/y/Fpm0/h7eBdouezvfWHba1fbG/90WVbW1hYqLVemSqatKFGjRo4e/as2uubm5vD3Ny80HIjIyOdfEEkSdLZtkkV21p93et2R4fqHbA/fj/+E/UfrHtnXYm3wfbWH7a1frG99UdXba3u9srUO+zk5ARjY2MkJSWpLE9KSoKbm5tO9x0REQEfH5+X9koRlVcFE15KkPDL+V9w4s4JuSMREeldmSqazMzMEBAQgH379imXKRQK7Nu3Dy1atNDpvkNDQ3Hx4kWcOnVKp/shMlT+bv4I8X9+9mzY7jBOeElEFY7BFU2ZmZmIjY1FbGwsACA+Ph6xsbG4desWACAsLAzLly/H6tWrcenSJYwYMQJZWVkYNGiQjKmJKoaZQTNhZWqFY7ePYfOlzXLHISLSK4Mrmk6fPo3GjRujcePGAJ4XSY0bN0Z4eDgAoE+fPpg/fz7Cw8Ph7++P2NhY7Ny5s9DgcCLSvip2VTCh5QQAwMS9E5GTl/OKZxARlR8GVzS1b98eQohCt1WrVinXGTlyJBISEpCTk4MTJ04gMDBQ57k4ponouQktJ8Ddxh03Ht1AxKkIueMQEemNwRVNhopjmoieszazxswOMwEAMw7NwMMnD2VORESkHyyaiKjEQvxC4Ovqi7TsNHx58Eu54xAR6QWLJiIqMWMjYyzotAAAsOT0Elx5eEXmREREuseiSU0c00SkKrhGMN6o/QbyFHmYuHei3HGIiHSORZOaOKaJqLCvX/8axpIxtl3ehoM3D8odh4hIp1g0EZHGfJx9MDxgOIDnE14qhELmREREusOiiYhKZVr7abA1s0XMvRisPbdW7jhERDrDoklNHNNEVDQXaxd83uZzAMDn+z/Hk9wnMiciItINFk1q4pgmouKNfW0sqtlXw52MO1j09yK54xAR6QSLJiIqNQsTC8ztOBcAMPfoXNzPvC9zIiIi7WPRRERa0bdhXzSv0hyZzzIRHhUudxwiIq1j0UREWiFJEhZ2WggAWHFmBc4nn5c5ERGRdrFoUhMHghO9WqtqrfBO/XegEAqM3z1e7jhERFrFoklNHAhOpJ6vgr+CqZEpdl3fhV3Xdskdh4hIa1g0EZFW1XSsiVHNRwEAxu0ehzxFnsyJiIi0w0TuAERU/vyn7X+w4swKXEi5gGkHpqFn3Z5ITU2FY74jjIye/1ZzsnJCNftqMiclIlIfiyYi0rrHzx4jKzcLADDr8CzMOjyr0DoWJhaIGxnHwomIygweniMirXvw5MErD8tl52XjwZMHekpERFR6LJqIiIiI1MDDc2qKiIhAREQE8vPz5Y5CVG4sj16OFp4t4GXvBS8HL1S1qwoTI/5ZIiLDxL9OagoNDUVoaCgyMjJgb28vdxyiciEyOhKR0ZHK+8aSMarYVVEWUV72Xir/X82+GixNLWVMTEQVGYsmIpLNW3XeQlZuFhLSE3Ar/Rae5T/DrfRbuJV+C4dvHS7yOS7WLvCy94K3g7dqcfX//7W34I8aItINFk1EJJup7aeiiXsTAIBCKHA/8z4S0hKQkJ7wv/++8P+ZzzKRnJWM5KxknEoseqJZe3P7Qr1ULxZYzlbOkCRJny+TiMoJFk1EZBCMJCN42HrAw9YDLTxbFHpcCIFH2Y9UiqqbaTdVCquHTx8iPScd55LO4VzSuSL3Y2liiWr21Yo8/Oft4A0PWw8YGxnr+uUSURnEoomItM7JygkWJhbIzssudh0LEws4WTmpvU1JkuBo6QhHS0c0dm9c5DqZzzJxK/1Wkb1VN9Nu4t7je3ia9xRxD+MQ9zCuyG2YGJmgql3Vl46rMjcxVzs3EZUfLJqISOuq2VdD3Mg45TxMCoXi+YzgjrqdEdzGzAY+zj7wcfYp8vFn+c9wO/12sYf/bqffRq4iFzfTbuJm2k0goej9uNm4vXRcla25rVZfFxEZBhZNRKQT1eyrKYsihUKBZONkuLi4KIsmOZgZm6GmY03UdKxZ5OP5ivzn46r+v2eqqMLqSe4T3M+8j/uZ93Hi7okit1PJopJKL5W3g7dKYVXZsrJWx1XdSr9VuEDlJWuItI5Fk5o4TxNR+Wds9HzKgyp2VdDSs2Whx4UQePj0YaHDf8qxVWkJeJT96Pnt/iPE3o8tcj9WplbFHv7zdvCGu607jCT1istb6bdQ97u6rzwUykvWEJUeiyY1cZ4mIpIkCU5WTnCyckKAR0CR6zzOeax6+O9fvVX3Mu/hSe4TXHpwCZceXCpyG6ZGpvC09yy2sPK094SZsRmA55eseVnBBPzvkjUsmohKh0UTEZEW2ZrboqFLQzR0aVjk49l52WqNq7rx6AZuPLpR5DYkSHC3dYe3gzdszTh+ikhfWDQREemRhYkFaleujdqVaxf5eL4iH4mPE186rio7LxuJjxOR+DhR7f1OOzAN9Z3qw8XaReXmbO0MZytnnhFIpAYWTUREBsTYyBie9p7wtPdE62qtCz0uhEDKkxRlAXX01lEsPrH4ldv9/crv+P3K78U+bm9uX6igKu7maOmo9pgrovKERRMRURkiSZKyeGlWpRlqVKqhVtE0stlImBqbIjkrGSlPUpQzqydnJSNPkYf0nHSk56TjaurVV27LSDKCs5UznK2d/1dMWRVfZNmY2XAWdioXWDQREVUAgxoPUl6y5kVCCKRlp6kUUSq3J6r3U5+mQiEUSMpKQlJWklr7tjCxKFxMWbmoFl0Fhwt5qJAMGIsmIqIKTJIkVLKshEqWlVDXqe4r18/Nz8WDJw+URdS/e63+fcvKzUJ2XrbyQszqMPRDhZwXq+Ji0UREVIbp4pI1L2NqbAp3W3e427qrtX7Ws6xXFlYvFl+aHCp0snIq1Iulq0OFnBerYmPRRERUhsl1yRp1WZtZw9rMGt4O3q9cV9NDhQX31fHioUJnq8KHB191qJDzYumXofXqsWgiIirjDPGSNZoo7aHCVx0u1MahQnV7qZIzk5GSlQIzYzOYGZvB1NgUxpIxB8SXgCH26rFoIiKiMklXhwoLiq+SHip8Udd1XQstkyApi6gXi6lCy4wKL1Oua/Tq5xe3jZLuS+4izxB79Vg0qYnXniMiKtu0cagw9n4slsUse+XzjSVj5AvVfy8EBHLyc5CTn6PpS9Crfxd5GhVoahR5xT0/IS1B7iYohEWTmnjtOSKiiqO4Q4Ux92LUKppODjsJfzd/5Obn4ln+M+Qqnv/337eCx1WWaWNdDbaRp8hTeQ1lrcjTBxZNREREOmAkGcHcxLzMzDulEAqVokruQi89Ox3XH12Xu1lUsGgiIiIigyvyYu7FIGBZgNwxVJStUyuIiIhkVDAv1stoc14sMizsaSIiIlKToc+LRbrFoomIiKgEysu8WIZO37Pdq4NFExERERkcQ+zVY9FEREREBsnQevXYl0hERESkBhZNRERERGpg0URERESkBhZNRERERGpg0URERESkBhZNRERERGrglAMlJIQAAGRkZGh92wqFAo8fP4aFhQUnSdMxtrV+sb31h22tX2xv/dFlWxf8m17wb3xxWDSV0OPHjwEAnp6eMichIiIibXr8+DHs7e2LfVwSryqrSIVCoUBiYiJsbW0hSZJWt52RkQFPT0/cvn0bdnZ2Wt02qWJb6xfbW3/Y1vrF9tYfXba1EAKPHz+Gh4fHS3ux2NNUQkZGRqhatapO92FnZ8cvn56wrfWL7a0/bGv9Ynvrj67a+mU9TAV4AJaIiIhIDSyaiIiIiNTAosmAmJubY+rUqTA3N5c7SrnHttYvtrf+sK31i+2tP4bQ1hwITkRERKQG9jQRERERqYFFExEREZEaWDQRERERqYFFkwE4dOgQunfvDg8PD0iShG3btskdqdyaM2cOmjVrBltbW7i4uKBnz56Ii4uTO1a5tHTpUvj6+irnVGnRogX++usvuWNVCHPnzoUkSRg7dqzcUcqladOmQZIklVu9evXkjlWu3b17F/3790flypVhaWmJRo0a4fTp03rPwaLJAGRlZcHPzw8RERFyRyn3Dh48iNDQUBw/fhx79uxBbm4uOnXqhKysLLmjlTtVq1bF3LlzER0djdOnT6NDhw7o0aMHLly4IHe0cu3UqVP4/vvv4evrK3eUcq1Bgwa4d++e8nbkyBG5I5Vbjx49QqtWrWBqaoq//voLFy9exIIFC1CpUiW9Z+GM4Aaga9eu6Nq1q9wxKoSdO3eq3F+1ahVcXFwQHR2Ntm3bypSqfOrevbvK/VmzZmHp0qU4fvw4GjRoIFOq8i0zMxMffPABli9fjpkzZ8odp1wzMTGBm5ub3DEqhK+++gqenp748ccflcuqV68uSxb2NFGFlp6eDgBwdHSUOUn5lp+fj/Xr1yMrKwstWrSQO065FRoaim7duiE4OFjuKOXe1atX4eHhgRo1auCDDz7ArVu35I5Ubm3fvh1NmzbFe++9BxcXFzRu3BjLly+XJQt7mqjCUigUGDt2LFq1aoWGDRvKHadc+ueff9CiRQtkZ2fDxsYGW7duhY+Pj9yxyqX169cjJiYGp06dkjtKuRcYGIhVq1ahbt26uHfvHqZPn442bdrg/PnzsLW1lTteuXPjxg0sXboUYWFh+Pzzz3Hq1CmMHj0aZmZmCAkJ0WsWFk1UYYWGhuL8+fMci6BDdevWRWxsLNLT07Fp0yaEhITg4MGDLJy07Pbt2xgzZgz27NkDCwsLueOUey8Op/D19UVgYCC8vLzw66+/YsiQITImK58UCgWaNm2K2bNnAwAaN26M8+fPIzIyUu9FEw/PUYU0cuRI7NixA1FRUahatarcccotMzMz1KpVCwEBAZgzZw78/PzwzTffyB2r3ImOjkZycjKaNGkCExMTmJiY4ODBg/jvf/8LExMT5Ofnyx2xXHNwcECdOnVw7do1uaOUS+7u7oV+aNWvX1+WQ6LsaaIKRQiBUaNGYevWrThw4IBsgwkrKoVCgZycHLljlDsdO3bEP//8o7Js0KBBqFevHiZOnAhjY2OZklUMmZmZuH79Oj788EO5o5RLrVq1KjQ1zJUrV+Dl5aX3LCyaDEBmZqbKL5T4+HjExsbC0dER1apVkzFZ+RMaGop169bht99+g62tLe7fvw8AsLe3h6WlpczpypfJkyeja9euqFatGh4/fox169bhwIED2LVrl9zRyh1bW9tC4/Ksra1RuXJljtfTgfHjx6N79+7w8vJCYmIipk6dCmNjY/Tr10/uaOXSp59+ipYtW2L27Nno3bs3Tp48iWXLlmHZsmX6DyNIdlFRUQJAoVtISIjc0cqdotoZgPjxxx/ljlbuDB48WHh5eQkzMzPh7OwsOnbsKHbv3i13rAqjXbt2YsyYMXLHKJf69Okj3N3dhZmZmahSpYro06ePuHbtmtyxyrXff/9dNGzYUJibm4t69eqJZcuWyZJDEkII/ZdqRERERGULB4ITERERqYFFExEREZEaWDQRERERqYFFExEREZEaWDQRERERqYFFExEREZEaWDQRERERqYFFExEREZEaWDQRkUGSJAnbtm3T6T4OHDgASZKQlpZW6m1NmzYN/v7+pd6OOtq3b4+xY8fqZV9E9D8smohI7+7fv49Ro0ahRo0aMDc3h6enJ7p37459+/Yp17l37x66du2q0xwtW7bEvXv3YG9vDwBYtWoVHBwcdLrPktBmUUdEpccL9hKRXt28eROtWrWCg4MDvv76azRq1Ai5ubnYtWsXQkNDcfnyZQCAm5vbS7eTm5sLU1PTUmUxMzN75X6IiAqwp4mI9OqTTz6BJEk4efIk3nnnHdSpUwcNGjRAWFgYjh8/rlzvxcNzN2/ehCRJ2LBhA9q1awcLCwusXbsWALBy5Uo0aNAA5ubmcHd3x8iRI1WeExsbq9xmWloaJEnCgQMHAKj25Bw4cACDBg1Ceno6JEmCJEmYNm1asa9j7ty5cHV1ha2tLYYMGYLs7OxC6/zwww+oX78+LCwsUK9ePSxZskT5WEG+9evXo2XLlrCwsEDDhg1x8OBB5eNBQUEAgEqVKkGSJAwcOFD5fIVCgc8++wyOjo5wc3N7aVYi0hJZLhNMRBXSw4cPhSRJYvbs2a9cF4DYunWrEEKI+Ph4AUB4e3uLzZs3ixs3bojExESxZMkSYWFhIRYvXizi4uLEyZMnxaJFi1Sec+bMGeU2Hz16JACIqKgoIYQQUVFRAoB49OiRyMnJEYsXLxZ2dnbi3r174t69e+Lx48dFZtuwYYMwNzcXP/zwg7h8+bKYMmWKsLW1FX5+fsp1fv75Z+Hu7q7Mu3nzZuHo6ChWrVqlkq9q1api06ZN4uLFi2Lo0KHC1tZWPHjwQOTl5YnNmzcLACIuLk7cu3dPpKWlCSGEaNeunbCzsxPTpk0TV65cEatXrxaSJIndu3eX7A0hohJh0UREenPixAkBQGzZsuWV6xZVNC1evFhlHQ8PDzFlypQin1/SokkIIX788Udhb2//ymwtWrQQn3zyicqywMBAlaKpZs2aYt26dSrrzJgxQ7Ro0UIl39y5c5WP5+bmiqpVq4qvvvqqyHwF2rVrJ1q3bq2yrFmzZmLixImvzE5EmuPhOSLSGyFEqZ7ftGlT5f8nJycjMTERHTt2LG2sErt06RICAwNVlrVo0UL5/1lZWbh+/TqGDBkCGxsb5W3mzJm4fv16sc8zMTFB06ZNcenSpVdm8PX1Vbnv7u6O5ORkTV4OEamJA8GJSG9q164NSZKUg71LytraWvn/lpaWL13XyOj5b8IXC7Xc3FyN9ltSmZmZAIDly5cXKq6MjY21so9/D4KXJAkKhUIr2yaiorGniYj0xtHREZ07d0ZERASysrIKPV6SU+ttbW3h7e2tMk3Bi5ydnQE8n7qgwIuDwotiZmaG/Pz8V+67fv36OHHihMqyFwexu7q6wsPDAzdu3ECtWrVUbtWrVy/2eXl5eYiOjkb9+vWVeQColYmIdI89TUSkVxEREWjVqhWaN2+OL7/8Er6+vsjLy8OePXuwdOlStQ5NFZg2bRo+/vhjuLi4oGvXrnj8+DGOHj2KUaNGwdLSEq+99hrmzp2L6tWrIzk5Gf/5z39euj1vb29kZmZi37598PPzg5WVFaysrAqtN2bMGAwcOBBNmzZFq1atsHbtWly4cAE1atRQrjN9+nSMHj0a9vb26NKlC3JycnD69Gk8evQIYWFhKu1Ru3Zt1K9fH4sWLcKjR48wePBgAICXlxckScKOHTvwxhtvwNLSEjY2Nmq3DxFpF3uaiEivatSogZiYGAQFBWHcuHFo2LAhXn/9dezbtw9Lly4t0bZCQkKwePFiLFmyBA0aNMCbb76Jq1evKh9fuXIl8vLyEBAQgLFjx2LmzJkv3V7Lli3x8ccfo0+fPnB2dsa8efOKXK9Pnz744osv8NlnnyEgIAAJCQkYMWKEyjpDhw7FDz/8gB9//BGNGjVCu3btsGrVqkI9TXPnzsXcuXPh5+eHI0eOYPv27XBycgIAVKlSBdOnT8ekSZPg6uqqnE6BiOQhidKOzCQiohK7efMmqlevjjNnzujt8itEVDrsaSIiIiJSA4smIiIiIjXw8BwRERGRGtjTRERERKQGFk1EREREamDRRERERKQGFk1EREREamDRRERERKQGFk1EREREamDRRERERKQGFk1EREREamDRRERERKSG/wNKLePdbwolrQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6, 4))\n",
    "plt.plot(depths_local, errors_local, marker=\"s\", linestyle=\"-\", color=\"green\")\n",
    "plt.xlabel(\"Circuit depth\")\n",
    "plt.ylabel(\"Energy error $|1 - E/E_{DMRG}|$\")\n",
    "plt.title(\"Local sequential optimisation (tensor-network fitting)\")\n",
    "plt.yscale(\"log\")\n",
    "plt.grid(visible=True, alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b34c87e0",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "The results above show that:\n",
    "- **Global optimisation** with L-BFGS can already find a good approximation, but basin-hopping and sequential layer-wise initialization often improve the fidelity and energy.\n",
    "- **Local optimisation** (tensor-network fitting) yields very high fidelity (low infidelity) and energies very close to the DMRG reference. The sequential depth variant further improves convergence.\n",
    "\n",
    "Both approaches provide a classical pre-processing step to generate a high-quality initial state for quantum algorithms such as QPE. The choice between them depends on the available infrastructure (automatic differentiation, global vs local optimizers) and the desired accuracy"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}