01_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bcc1a161",
   "metadata": {},
   "source": [
    "# Ground-State Variational Search with NetKet\n",
    "\n",
    "## Filippo Vicentini (EPFL-CQSL)\n",
    "\n",
    "### Lemanic Quantum Science School, Geneva\n",
    "\n",
    "17 November, 2021\n",
    "\n",
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PhilipVinc/NetKet-lectures/blob/master/01_intro.ipynb)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae721781-2989-46c5-a92b-ca528b145e6c",
   "metadata": {},
   "source": [
    "In this Tutorial we will introduce the open-source package [NetKet](https://www.netket.org/), and show some of its functionalities. We will guide you through a relatively simple quantum problem, that however will be a good guide also to address more complex situations. \n",
    "\n",
    "Specifically, we will study the transverse-field Ising model in one dimension: \n",
    "\n",
    "$$ \n",
    "\\mathcal{H}=\\Gamma\\sum_{i}\\sigma_{i}^{(x)}+V\\sum_{i}\\sigma_{i}^{(z)}\\sigma_{i+1}^{(z)}. \n",
    "$$\n",
    "\n",
    "In the following we assume periodic boundary conditions and we will count lattice sites starting from $ 0 $, such that $ i=0,1\\dots L-1 $ and $i=L=0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39bddf19",
   "metadata": {},
   "source": [
    "## 0. Installing Netket \n",
    "\n",
    "Installing NetKet is relatively straightforward. For this Tutorial, if you are running it locally on your machine, we recommend that you create a clean virtual environment and install NetKet within: \n",
    "\n",
    "```bash\n",
    "python3 -m venv netket\n",
    "source netket/bin/activate\n",
    "pip install --upgrade \"netket>=3.2\"\n",
    "```\n",
    "\n",
    "**If you are on Google Colab**, uncomment the following cell and run it to install the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "43e97750",
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install --upgrade netket"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f15de72d-9454-4873-a88f-a02ec10459ed",
   "metadata": {},
   "source": [
    "We also want to sure that this notebook is running on the cpu. \n",
    "You can edit the field by changing \"cpu\" to \"gpu\" to make it runa on the GPU if you want. \n",
    "But you'll need to use much larger systems to see a benefit in the runtime.\n",
    "For systems with less than 40 spins GPUs slow you down remarkably."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f752288d-caf4-4f53-a104-d4b2efe6e0fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"JAX_PLATFORM_NAME\"] = \"cpu\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f424667",
   "metadata": {},
   "source": [
    "You can check that the installation was succesfull doing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5a4c83c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import netket as nk"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7deb0afe-29e4-40e0-b650-b6480a46a8ae",
   "metadata": {},
   "source": [
    "You should also check that your version of netket is at least 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "53d59c56-3a27-4f8b-a4fb-a9dde50e9bc5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NetKet version: 3.3.dev8+gd1eb9e981f.d20211202\n"
     ]
    }
   ],
   "source": [
    "print(f\"NetKet version: {nk.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8b7effe",
   "metadata": {},
   "source": [
    "## 1. Defining The Hamiltonian\n",
    "\n",
    "The first step in our journey consists in defining the Hamiltonian we are interested in. \n",
    "For this purpose, we first need to define the kind of degrees of freedom we are dealing with (i.e. if we have spins, bosons, fermions etc). This is done specifying the Hilbert space of the problem. For example, let us concentrate on a problem with 20 spins. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "65e5273d",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20\n",
    "\n",
    "hi = nk.hilbert.Spin(s=1 / 2, N=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4b2bf39c-ae7f-4cb0-b141-093c448138e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([[ 1., -1., -1.,  1., -1.,  1.,  1., -1., -1., -1., -1., -1.,\n",
       "              -1., -1., -1., -1.,  1.,  1., -1., -1.],\n",
       "             [ 1., -1.,  1.,  1., -1.,  1.,  1., -1.,  1., -1., -1., -1.,\n",
       "               1.,  1.,  1.,  1., -1.,  1., -1.,  1.],\n",
       "             [-1.,  1.,  1.,  1.,  1., -1., -1.,  1.,  1.,  1., -1.,  1.,\n",
       "               1.,  1.,  1.,  1., -1.,  1., -1.,  1.]], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import jax\n",
    "hi.random_state(jax.random.PRNGKey(0), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc2e67ac",
   "metadata": {},
   "source": [
    "We now need to specify the Hamiltonian. For this purpose, we will use NetKet's ```LocalOperator``` (see details [here](https://www.netket.org/docs/_generated/operator/netket.operator.LocalOperator.html#netket.operator.LocalOperator)) which is the sum of arbitrary k-local operators. \n",
    "\n",
    "In this specifc case, we have a 1-local operator, $ \\sigma^{(x)}_i $ and a 2-local operator, $ \\sigma^{(z)}_i \\sigma^{(z)}_j $. We then start importing the pauli operators. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b018e4d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netket.operator.spin import sigmax,sigmaz "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c35665",
   "metadata": {},
   "source": [
    "We now take $ \\Gamma=-1 $ and start defining the 1-local parts of the Hamiltonian "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b35bde0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gamma = -1\n",
    "H = sum([Gamma*sigmax(hi,i) for i in range(N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b41d740",
   "metadata": {},
   "source": [
    "Here we have used a list comprehension to (mildly) show off our ability to write one-liners, however you could have just added the terms one by one in an explicit loop instead (though you'd end up with a whopping 3 lines of code). \n",
    "\n",
    "We now also add the interaction terms, using the fact that NetKet automatically recognizes products of local operators as tensor products. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eb378d9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "V=-1\n",
    "H += sum([V*sigmaz(hi,i)*sigmaz(hi,(i+1)%N) for i in range(N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59025d8c-dede-496f-b2f0-6ad63966be1d",
   "metadata": {},
   "source": [
    "In general, when manipulating NetKet objects, you should always assume that you can safely operate on them like \n",
    "you would in mathematical equations, therefore you can sum and multiply them with ease."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65581d8b",
   "metadata": {},
   "source": [
    "## 2. Exact Diagonalization\n",
    "\n",
    "Now that we have defined the Hamiltonian, we can already start playing with it. For example, since the number of spins is large but still manageable for exact diagonalization, we can give it a try. \n",
    "\n",
    "In NetKet this is easily done converting our Hamiltonian operator into a sparse matrix of size $ 2^N \\times 2^ N $. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d36dd96c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1048576, 1048576)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp_h=H.to_sparse()\n",
    "sp_h.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0b786e6",
   "metadata": {},
   "source": [
    "Since this is just a regular scipy sparse matrix, we can just use any sparse diagonalization routine in there to find the eigenstates. For example, this will find the two lowest eigenstates  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1123bc8c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalues with scipy sparse: [-25.49098969 -25.41240947]\n"
     ]
    }
   ],
   "source": [
    "from scipy.sparse.linalg import eigsh\n",
    "\n",
    "eig_vals, eig_vecs = eigsh(sp_h, k=2, which=\"SA\")\n",
    "\n",
    "print(\"eigenvalues with scipy sparse:\", eig_vals)\n",
    "\n",
    "E_gs = eig_vals[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e51cc79",
   "metadata": {},
   "source": [
    "## 3. Mean-Field Ansatz\n",
    "\n",
    "We now would like to find a variational approximation of the ground state of this Hamiltonian. As a first step, we can try to use a very simple mean field ansatz: \n",
    "\n",
    "$$ \\langle \\sigma^{z}_1,\\dots \\sigma^{z}_N| \\Psi_{\\mathrm{mf}} \\rangle = \\Pi_{i=1}^{N} \\Phi(\\sigma^{z}_i), $$\n",
    "\n",
    "where the variational parameters are the single-spin wave functions, which we can further take to be normalized: \n",
    "\n",
    "$$ |\\Phi(\\uparrow)|^2 + |\\Phi(\\downarrow)|^2 =1, $$\n",
    "\n",
    "and we can further write $ \\Phi(\\sigma^z) = \\sqrt{P(\\sigma^z)}e^{i \\phi(\\sigma^z)}$. In order to simplify the presentation, we take here and in the following examples the phase $ \\phi=0 $. In this specific model this is without loss of generality, since it is known that the ground state is real and positive. \n",
    "\n",
    "For the normalized single-spin probability we will take a sigmoid form: \n",
    "\n",
    "$$ P(\\sigma_z; \\lambda) = 1/(1+\\exp(-\\lambda \\sigma_z)), $$ thus depending on the real-valued variational parameter $\\lambda$. \n",
    "\n",
    "In NetKet one has to define a variational function approximating the **logarithm** of the wave-function amplitudes (or density-matrix values).\n",
    "We call this variational function _the Model_ (yes, caps on the M).\n",
    "\n",
    "$$ \\langle \\sigma^{z}_1,\\dots \\sigma^{z}_N| \\Psi_{\\mathrm{mf}} \\rangle = \\exp\\left[\\mathrm{Model}(\\sigma^{z}_1,\\dots \\sigma^{z}_N ; \\theta ) \\right], $$ where $\\theta$ is a set of parameters. \n",
    "\n",
    "In this case, the parameter of the model will be just one: $\\gamma$.  \n",
    "\n",
    "The Model can be defined using one of the several *functional* jax frameworks such as Jax/Stax, Flax or Haiku. \n",
    "NetKet includes several pre-built models and layers built with [Flax](https://github.com/google/flax), so we will be using it for the rest of the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "967e1611",
   "metadata": {},
   "outputs": [],
   "source": [
    "# numerical operations in the model should always use jax.numpy \n",
    "# instead of numpy because jax supports computing derivatives. \n",
    "# If you want to better understand the difference between the two, check\n",
    "# https://flax.readthedocs.io/en/latest/notebooks/jax_for_the_impatient.html\n",
    "import jax.numpy as jnp\n",
    "\n",
    "# Flax is a framework to define models using jax\n",
    "import flax\n",
    "# we refer to `flax.linen` as `nn`. It's a repository of \n",
    "# layers, initializers and nonlinear functions.\n",
    "import flax.linen as nn\n",
    "\n",
    "# A Flax model must be a class subclassing `nn.Module`\n",
    "class MF(nn.Module):\n",
    "    \n",
    "    # The most compact way to define the model is this.\n",
    "    # The __call__(self, x) function should take as \n",
    "    # input a batch of states x.shape = (n_samples, L)\n",
    "    # and should return a vector of n_samples log-amplitudes\n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "        \n",
    "        # A tensor of variational parameters is defined by calling\n",
    "        # the method `self.param` where the arguments will be:\n",
    "        # - arbitrary name used to refer to this set of parameters\n",
    "        # - an initializer used to provide the initial values. \n",
    "        # - The shape of the tensor\n",
    "        # - The dtype of the tensor.\n",
    "        lam = self.param(\n",
    "            \"lambda\", nn.initializers.normal(), (1,), float\n",
    "        )\n",
    "        \n",
    "        # compute the probabilities\n",
    "        p = nn.log_sigmoid(lam*x)\n",
    "\n",
    "        # sum the output\n",
    "        return 0.5 * jnp.sum(p, axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ccb969-487e-43e6-ac87-50cde68c68e1",
   "metadata": {},
   "source": [
    "The model itself is only a set of instructions on how to initialise the parameters and how to compute the result. \n",
    "\n",
    "To actually create a variational state with its parameters, the easiest way is to construct a Monte-Carlo-sampled Variational State. \n",
    "To do this, we first need to define a sampler.\n",
    "\n",
    "In `netket.sampler` several samplers are defined, each with its own peculiarities. \n",
    "In the following example, we will be using a simple sampler that flips the spins in the configurations one by one.\n",
    "\n",
    "You can read more about how the sampler works by checking the documentation with `?nk.sampler.MetropolisLocal`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c2cb0ac3-bc4c-4743-adec-03bb4e2d7c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an instance of the model. \n",
    "# Notice that this does not create the parameters.\n",
    "mf_model=MF()\n",
    "\n",
    "# Create the local sampler on the hilbert space\n",
    "sampler = nk.sampler.MetropolisLocal(hi)\n",
    "\n",
    "# Construct the variational state using the model and the sampler above.\n",
    "# n_samples specifies how many samples should be used to compute expectation\n",
    "# values.\n",
    "vstate = nk.vqs.MCState(sampler, mf_model, n_samples=512)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "57cfd849-40f7-4d44-b563-4f0a749f0687",
   "metadata": {},
   "source": [
    "You can play around with the variational state: for example, you can compute expectation values yourself or inspect it's parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "578fae25-6b39-41c1-a25e-51de791b7aa3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FrozenDict({\n",
      "    lambda: DeviceArray([0.01556348], dtype=float64),\n",
      "})\n"
     ]
    }
   ],
   "source": [
    "# you can inspect the parameters which contain the single\n",
    "# variational parameter `lambda`\n",
    "print(vstate.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0ea20c03-dcfa-4cda-8f2f-b1d86a95872b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-19.77 ± 0.19 [σ²=19.40, R̂=0.9995]\n"
     ]
    }
   ],
   "source": [
    "# Expectation value: notice that it also provides an error estimate.\n",
    "E = vstate.expect(H)\n",
    "print(E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "df25e787-deb1-4c3c-a6a8-5150948a6f3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean                  : -19.774195046540363\n",
      "Error                 : 0.1946632634109754\n",
      "Variance              : 19.40161849436712\n",
      "Convergence indicator : 0.9994552457867388\n",
      "Correlation time      : 1.1102230246251565e-16\n"
     ]
    }
   ],
   "source": [
    "# the energy (expectation value) is a structure with a lot of fields:\n",
    "print(\"Mean                  :\", E.mean)\n",
    "print(\"Error                 :\", E.error_of_mean)\n",
    "print(\"Variance              :\", E.variance)\n",
    "print(\"Convergence indicator :\", E.R_hat)\n",
    "print(\"Correlation time      :\", E.tau_corr)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c161f0a6-3f95-4f6e-9dd2-d2a71527efce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-19.77 ± 0.19 [σ²=19.40, R̂=0.9995],\n",
       " FrozenDict({\n",
       "     lambda: DeviceArray([-0.23529235], dtype=float64),\n",
       " }))"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vstate.expect_and_grad(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b5b6302",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 4. Variational Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ab77267-892e-48d2-97e2-eee6c8eae796",
   "metadata": {},
   "source": [
    "We will now try to optimise $ \\lambda $ in order to best approximate the ground state of the hamiltonian.\n",
    "\n",
    "At first, we'll try to do this by ourself by writing the training loop, but then we'll switch to using a pre-made\n",
    "solution provided by netket for simplicity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d103e357-f079-4816-8323-d9a606748a04",
   "metadata": {},
   "source": [
    "### 4a. DIY Optimisation loop\n",
    "\n",
    "The optimisation (or training) loop must do a very simple thing: at every iteration it must compute the energy and it's gradient, then multiply the gradient by a certain learning rate $\\lambda = 0.05$ and lastly it must update the parameters with this rescaled gradient.\n",
    "\n",
    "You can do so as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ef4b0687-1f78-4c2e-9f17-16c3fac2f064",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [00:01<00:00, 75.10it/s]\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "energy_history = []\n",
    "n_steps = 100\n",
    "\n",
    "# For every iteration (tqdm is just a progress bar)\n",
    "for i in tqdm(range(n_steps)):\n",
    "    # compute energy and gradient of the energy\n",
    "    E, E_grad = vstate.expect_and_grad(H)\n",
    "    # log the energy to a list\n",
    "    energy_history.append(E.mean.real)\n",
    "    # equivalent to vstate.parameters - 0.05*E_grad , but it performs this\n",
    "    # function on every leaf of the dictionaries containing the set of parameters\n",
    "    new_pars = jax.tree_multimap(lambda x,y: x-0.05*y, vstate.parameters, E_grad)\n",
    "    # actually update the paramters\n",
    "    vstate.parameters = new_pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "747f2a0e-3408-443d-8f9b-13f70db15231",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f85f8164910>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsPklEQVR4nO3deXxU5dn/8c+VfZnsK9mBhCUgiwRBRFABRYs72ta6UG2pT221WmtL7fpr+7Rqq9bHfbdqXaviXgEF2THse0ggCRDIvu/L/ftjJiGBLMBMmGTmer9eeZG5z5mZ63DgO2fuc5/7iDEGpZRSrsnD2QUopZTqPxrySinlwjTklVLKhWnIK6WUC9OQV0opF+bl7AI6i4yMNCkpKc4uQymlBpWNGzeWGGOiuls2oEI+JSWFzMxMZ5ehlFKDiojk9bRMu2uUUsqFacgrpZQL05BXSikXpiGvlFIuTENeKaVcmIa8Ukq5MA15pZRyYS4R8ocr6vnbZ3s4Wtng7FKUUmpAcYmQr21s4ekVOXy5p8jZpSil1IDiEiGfFm0hPtRfQ14ppY7jEiEvIlw0KprV2SU0NLc6uxyllBowXCLkAS4aFU19cyvrD5R1u7yltY173trCqn0lJyzbc7SKQ+V1/V2iUkqdcS4T8ucOj8DP24OveuiyeXfjId7bfJhPthecsOx/XtvEDc+tp7axpUv7q2tzeWp5Tr/Uq5RSZ4JdIS8iD4nIHhHZJiLvi0hop2WLRCRbRPaKyCV2V9oHP29Ppg2P5Ms9RRx/c/L6plYeWZoFQG5J1yP2ppY28kpryS+r438/3d3R/vmOo/x28U7+8cVeSmoa+7t8pZTqF/YeyS8BxhpjxgFZwCIAEUkHvgOMAeYCT4qIp53v1acLR0WTX1ZHTnFtl/aX1hygsKqR1GgLuaVdl+WX1dFmYFhkIK+vz2dFVjHZRTXc+85WhkUG0tJmWLzlxKN/pZQaDOwKeWPMF8aY9j6OdUCC7fcrgTeNMY3GmANANnCOPe91Mi4aFQ3QpcumvLaJp5bnMHt0NFeMj+NIZUOXk7O5JdbQ/8vVZ5EabeGX727j9tc24uvlwWs/mMK4hBDe3Xiov0tXSql+4cg++VuBz2y/xwMHOy07ZGvrV/Gh/oyMCeoylPKJr7KpbWzhF5eMIjkiAIC80mNdNu1H9qOHBPHw9eMprmlkf3EN//fdicSF+jN/UgK7j1Sxs6Cyv8tXSimH6/POUCKyFIjtZtH9xpjFtnXuB1qA10+1ABFZCCwESEpKOtWnn+DCUdE8v3I/r6/P45NtR1i7v5T5ZycwMjaIxhbrEXxuaS0jY4MA2F9SS2iAN6EBPoQG+PDP70zAQ4RpqZEAXD4ujj9/vJt3Nx5iTFyI3fUppdSZ1GfIG2Nm97ZcRBYA84BZ5tgZz8NAYqfVEmxt3b3+s8CzABkZGaa7dU7FrNHRPL0ih/vf30FyRAA/vSiNhTOGAZAcEQgc66Jp/z3F1g4wb1xcl9cLC/Rhdno0i7cUsOjS0fh4df3ys7+4hpfX5CKAt6cHsSF+3HreUDw8xN5NUUopu9l1j1cRmQvcB8w0xnQetvIh8G8ReRiIA9KADfa818nKSA7jofnjSI22MCExFJFjYRvi7014oA+5nbtrSmqZMiyi19ecPymBT7cf5au9RVwypuuXmudXHeDNDfkE+3vT0NxKQ3MbE5NCmZQc7tgNU0qp02Bvn/zjQBCwRES2iMjTAMaYncDbwC7gc+AOY8wZuRRVRLguI5GJSWFdAr5dckRAx5F8Q3MrBZUNXY7kuzMjLYpIi+8JJ2CNMazYW8zs0TFs+d3FbLh/Nl4ewtLdPU+vUFTVwKNLs3RYplLqjLB3dE2qMSbRGDPB9nN7p2V/McYMN8aMNMZ81tvrnElDIwLJs51sbT8BOzSq95D38vTg6olxLN9bREVdU0d7dlENhyvquWCkdVRPsJ835wwNZ9nuwhNewxjD+5sPMeeRr3l06T4e+nxvl+VVDc1c9s+VfLhVh2sqpRzHZa54PVnJEYEU2IZRHiipAazB35crJ8TT3Gr4bMfRjrYVWcUAzBwZ1dE2a3QMWYU15HfqEqpvamXhqxu5+62tDI8K5PLxcby76RA5xTUd6/xz6T52HaninczOg5KUUso+bhfyKZHWYZT5ZXUcsF392t7WmzFxwQyLDGTxlmPnj5fvLe6YAbPd7NHWo/qlnY7mX1x9gCW7Cll06SjeuX0av5uXjq+XB48ssV6Fu6+wmlfW5BLo48m6/aXUHDe9glJKnS73C/lOI2xyS2qJtPgQ5Ofd5/NEhCsmxLH+QBlHKxuobWxhw4EyLuh0FA/Wbwqp0RaW7bGGfG1jC8+v3M+FI6P40czheHoIUUG+3HreUD7edoRdBVX88aNdBPh48vfrxtPcalhp+4aglFL2ct+QL63lQGltnyddO7tifBzGwMfbCli3v5Sm1raO/vjOZo2OZv3+MqoamvnX2jzK65q5c1Zal3V+OGMYwX5eLHw1k1XZJdwzZwRz0mMI9vNimc6Lr5RyELcL+ZAAb8ICvMktrbOOkY88+ZAfFmVhbHwwH24tYPneYgJ8PMlICTthvTmjY2hpM3y+4yjPrdzPjBFRTEzqul6Ivze3XzCcQ+X1jIixcOPUZLw8PbhgZDRf7Smirc3uSwaUUsr9Qh6sXSo7C6ooqm5k6CmEPMCV4+PZdqiSD7cWMG14BL5eJ867NjEpjPBAH/700S7Kapu4a1Zqt6+1YFoK8ycl8ND88Xh5WnfFrNHRlNY2sfVQxSlvl1JKHc8tQ35oZCDbbCF6qiE/b/wQRKCyvpmZ3XTVAHh6CBeMjKK6sYXpqZE9XhgV4OPF368bz/jE0I62mSOi8PQQvZWhUsoh3DLkkyMCaJ+A4VT65AGGhPhzToo1tC8YEdXjevPGDcFD4Gez03pcpzuhAT5MSg7r9YIqpZQ6WW4Z8p2P3k9m+OTx7p4zgh9fMJzE8J6fe9GoGDbcP5uMlFOf3mDWqGh2H6mioKL+lJ+rlFKduWXIt09UFhPsS4DPqU/fM3VYBPfNHdXnepEW31N+bbBeUAXoKBullN3cMuTbr3A91f74M2V4VCDJEQHdTo+glFKnwi1DPiTAm9hgP0bFBju7lG6JCBenx7AmW69+VUrZxy1DHuA/P57Gzy8e4ewyejQnPZam1jZW7NWrX5VSp89tQz4+1P+kpjNwlknJ1rH2X+w62vfKwPMr97PovW39XJVSarBx25Af6Dw9hFmjovlyTxHNrW29rmuM4cVVB1i8pYBjN+dSSikN+QFtTnoM1Q0trN9f1ut6OwuqKKhsoK6plZKapl7XVUq5Fw35Aez8tCj8vD1Y0keXTedpjdtviKKUUqAhP6D5+3hyfloUS3YV9toNs3R3ITHB1jH5eaV1Pa6nlHI/GvID3Jz0GAoqG9hZUNXt8iOV9ew4XMWNU5LxED2SV0p1pSE/wM0aFY2HwOc7uu+yWbrL2lVz6VlDiAv1J69Mj+SVUsdoyA9wERZfZo6I4vX1eVQ3NJ+wfMnuIoZGBjI8KpCUiEBytbtGKdWJhvwgcPecEZTXNfPS6twu7dUNzazNKWH26GhEhKSIAPK1u0Yp1YmG/CAwLiGUi9NjeO7r/VTUHRsiuXJfCc2thjnpsQCkRARQXtdMZf2JR/xKKfekIT9I3HPxCGqaWnhu5X4AymubeHl1LmEB3pydFApAUrh1wrV87bJRStmc+jy7yilGxQYzb1wcL63OZWikhQc+30NFXRN/vGJsx60D2+fGzy2t5ayEEGeWq5QaIPRIfhD52ew0GppbufedrUQE+rD4juncMCWpY3mS7SYm+TrCRillo0fyg8jwKAu/+VY69c2t/PD8Yfh4df2MDvDxIjrIl9wSPfmqlLLSkB9kbp0+tNflyREBOlZeKdXBru4aEXlIRPaIyDYReV9EQm3tESLylYjUiMjjDqlUnZTkiMAuV71+su0INzy3jrY2nZ1SKXdkb5/8EmCsMWYckAUssrU3AL8F7rXz9dUpSokIoLCqkfqmVtraDP9Yspc1OaUUVOpNwZVyR3aFvDHmC2NM+/3p1gEJtvZaY8wqrGGvzqAk2/1r88vqWJldwv5i61H9vqIaZ5allHISR46uuRX47FSfJCILRSRTRDKLi/VWd/ZKibCOsMkrreXl1QcIDbDe/Sq7UENeKXfUZ8iLyFIR2dHNz5Wd1rkfaAFeP9UCjDHPGmMyjDEZUVFRp/p0dZxk2wVRy7OK+WpvMQumpRBp8WVfUbWTK1NKOUOfo2uMMbN7Wy4iC4B5wCyj955zupAAb0IDvHlzQz7ensINU5JYv79Mu2uUclP2jq6ZC9wHXGGM0XF7A0RyeABtBuaNiyM6yI+0GAvZhTV6/1el3JC9ffKPA0HAEhHZIiJPty8QkVzgYWCBiBwSkXQ730udpGTbydcF01IASIu2UN3YwtEqPQ+ulLux62IoY0xqL8tS7HltdfrmT0pgSKgf4xNDAUiLCQJgX2ENQ0L8nViZUupM07lrXNCMEVEsunR0x+O0aAugwyiVckca8m4gwuJLeKAP2TrCRim3oyHvJlKjLezTsfJKuR0NeTeRFm1hX5GOsFHK3WjIu4m0aAuV9c0U1zQ6uxSl1BmkIe8mRnQaYaOUch8a8m4iNcY2wqZQT74q5U405N1ElMWXEH9vHUaplJvRkHcTItJx8lUp5T405N1IWoyFbA15pdyKhrwbSYsOoqy2iaOVOoeNUu5CQ96NnDM0HID1B0qdXIlS6kzRkHcjo4cEE+LvzZpsDXml3IWGvBvx9BCmDgtnzf4SZ5eilDpDNOTdzLThkRwsq+dgmd7jRSl3oCHvZqYNjwBgbY522SjlDjTk3UxqtIVIiy9rcrTLRil3oCHvZkSEc4dHsCanVGekVMoNaMi7oWnDIyiqbiSnuNbZpSil+pmGvBs61i+vXTZKuToNeTeUFB5AfKg/a/Tkq1IuT0PeDbX3y6/dX0pbm/bLK+XKNOTd1PTUSCrqmsnMK3d2KUqpfqQh76YuHhNDkJ8Xr67Lc3YpSql+pCHvpgJ8vLhuUiKf7zhCUbXOSqmUq9KQd2M3Tk2iudXw5oaDzi5FKdVPNOTd2LAoC+enRfLv9fm0tLY5uxylVD/QkHdzN01N5mhVA0t3Fzq7FKVUP7Ar5EXkIRHZIyLbROR9EQm1tc8RkY0ist3250UOqVY53KzRMcSH+vOvtXoCVilXZO+R/BJgrDFmHJAFLLK1lwCXG2POAm4BXrXzfVQ/8fQQbpiSxJqcUg6U6DQHSrkau0LeGPOFMabF9nAdkGBr32yMKbC17wT8RcTXnvdS/efSsbEAbNDbAirlchzZJ38r8Fk37dcCm4wxjd09SUQWikimiGQWFxc7sBx1slIiAgn282LroUpnl6KUcjCvvlYQkaVAbDeL7jfGLLatcz/QArx+3HPHAA8AF/f0+saYZ4FnATIyMvQaeyfw8BDGJ4ay9WCFs0tRSjlYnyFvjJnd23IRWQDMA2aZThOUi0gC8D5wszEmx846VT8bnxDKUytyaGhuxc/b09nlKKUcxN7RNXOB+4ArjDF1ndpDgU+AXxljVttVoTojxieG0tpm2FmgXTZKuRJ7++QfB4KAJSKyRUSetrX/BEgFfmdr3yIi0Xa+l+pH4xNCANhyUENeKVfSZ3dNb4wxqT20/xn4sz2vrc6s6GA/hoT4ab+8Ui5Gr3hVHcYnhLL1UIWzy1BKOZCGvOowPjGUvNI6ymubnF2KUspBNORVh/GJ1n55PZpXynVoyKsOZ8WHIALb9KIopVyGhrzqEOTnTWqURU++KuVCNORVF+NsJ187XdemlBrENORVFxMSQyipaSKnuMbZpSilHMCucfLK9UwdFoGXh3DJoyuZMjScS88awnWTEnSqA6UGKT2SV12kxQTx0U+nc/vMYRytauC3H+xg8ZbDzi5LKXWaNOTVCUYPCeYXl4xi2T0zsfh6saugytklKaVOk4a86pGIMCLGwp6j1c4uRSl1mjTkVa9Gxgaxt7BaR9soNUhpyKtejYwJoqKumeLqbm/spZQa4DTkVa9GxgYDaJeNUoOUhrzq1cjYIAD2asgrNShpyKtehQf6EBXkq0fySg1SGvKqT6Nig8gq1JBXajDSkFd9GhljDfnWNh1ho9RgoyGv+jQyNojGljbySmudXYpS6hRpyKs+6clXpQYvDXnVp7ToIERgr/bLKzXoaMirPvn7eJISEahH8koNQhry6qSMjAnSkFdqENKQVydlRGwQuaW1NDS3AtDaZnQ+G6UGAQ15dVJGxQbRZuDT7UdY9N420n/3OS+uznV2WUqpPmjIq5PSPsLmnre38t6mw/h4efDlnkInV6WU6ove/k+dlKERgdwwJYn4UH++e04Sjy7N4t2Nh2hpbcPLU48VlBqo9H+nOikeHsL/Xn0Wd1yYSnigD5OSw6hratU5bZQa4OwKeRF5SET2iMg2EXlfREJt7eeIyBbbz1YRudoh1aoBIyMlHIDM3DInV6KU6o29R/JLgLHGmHFAFrDI1r4DyDDGTADmAs+IiHYNuZD4UH+GhPiRmVfu7FKUUr2wK+SNMV8YY1psD9cBCbb2uk7tfoCOtXNBk5LD2Kghr9SA5sg++VuBz9ofiMgUEdkJbAdu7xT6XYjIQhHJFJHM4uJiB5aj+ltGchhHKhs4XFHv7FKUUj3oM+RFZKmI7Ojm58pO69wPtACvt7cZY9YbY8YAk4FFIuLX3esbY541xmQYYzKioqLs3yJ1xmi/vFIDX5/95MaY2b0tF5EFwDxglunmEkhjzG4RqQHGApmnWacagEbFBhHg48mmvHKunBDv7HKUUt2wd3TNXOA+4ApjTF2n9qHtJ1pFJBkYBeTa815q4PHy9GBiUqiefFVqALO3T/5xIAhYYhsu+bStfTqwVUS2AO8DPzbGlNj5XmoAmpQUxu4jVdQ0dnvKRSnlZHYNazTGpPbQ/irwqj2vrQaHSSnhtBnYkl/B9LRIZ5ejlDqOXvGq7DIxKRQRyMzTk69KDUQa8souwX7eDIsMZFdBlbNLUUp1Q0Ne2S012kJ2cY2zy1BKdUNDXtktNdpCXmkdTS1tzi5FKXUcDXllt9RoC61thrzSWmeXopQ6joa8sltqlPWGIjnaZaPUgKMhr+w2LCoQgOwiDXmlBhoNeWW3QF8v4kL8Tink9SbgSp0ZGvLKIYaf5Aib+qZWfvBKJlc/uUaDXqkzQENeOURqtIWcolra2o4F98GyOnYcruwI88q6Zm58YT1Ldxey5WAF+0v0RK1S/U3v1qQcIjXaQn1zKwWV9SSEBQDwo1c3sutIFanRFq6eGM9HWwvIKa7h/stG85dPd/N1VjHDoyxOrlwp16ZH8sohUm1h3d4vn19ax64jVcwdE0uovzcP/Xcv+WV1vLhgMj+cMYxhUYGsyNKbxCjV3/RIXjlEavSxkL9gZDRf7DoKwKLLRpEcEUh+aR0GQ3KEdSTOzBFRvLEhn4bmVvy8PZ1Wt1KuTo/klUNEWHwJC/Amp9jaz75kVyGjYoM6Qj0pIqDjd7CGfENzGxsO6MRmSvUnDXnlMMOjLOQU1VBW28Q3uWXMSY/pcd2pwyLw9fLQLhul+pmGvHKY9onKlu0upM3AxemxPa7r5+3JlGERGvJK9TMNeeUwqdEWymqbeDvzIENC/BgbH9zr+jPSIskuquFQeV2v6ymlTp+GvHKY4baTr9/kljMnPQYR6XX9C0ZGAfB1lt4ZUqn+oiGvHCa105j33rpq2g2PshAf6s+KrKL+LEspt6YhrxwmPtQff29Pgvy8mDIsvM/1RYQZIyJZk13a5UpZpZTjaMgrh/HwEM5LjeCaifF4e57cP62z4kOpbmzhcEV9P1enlHvSi6GUQz1/y+RTWn9EjLWLZ19RNYnhAf1RklJuTY/klVOlxVhvOJJVqHPRK9UfNOSVU4X4exMT7EtWYbWzS1HKJWnIK6dLiw5inx7JK9UvNOSV06XFWMguqtERNkr1Aw155XQjYoKob27VETZK9QMNeeV07SNstF9eKcezK+RF5CER2SMi20TkfREJPW55kojUiMi9dlWpXFpqtI6wUaq/2HskvwQYa4wZB2QBi45b/jDwmZ3voVxc+wibfXokr5TD2RXyxpgvjDEttofrgIT2ZSJyFXAA2GnPeyj3MCImiH1FeiSvlKM5sk/+VmxH7SJiAX4J/LGvJ4nIQhHJFJHM4mKdW9xdpUUH6QgbpfpBnyEvIktFZEc3P1d2Wud+oAV43db0B+ARY0yfh2bGmGeNMRnGmIyoqKjT3Aw12I2IsVDf3Mqhch1ho5Qj9Tl3jTFmdm/LRWQBMA+YZYxpPwybAswXkQeBUKBNRBqMMY/bV65yVcemN6gmKULnsFHKUewdXTMXuA+4whjTcXsfY8z5xpgUY0wK8CjwvxrwqjepthuOZBVZT76+ujaXGQ9+RU6x9tMrZQ97++QfB4KAJSKyRUSedkBNyg2F+HsTG+xHdmEN72Qe5LeLd5JfVsedb2ymsaXV2eUpNWjZNdWwMSb1JNb5gz3vodxHWoyFL/cW8cGWw5yfFsl3z0nix69v4qHP9/KbeenOLk+pQUmveFUDxoiYICrqmpmYFMYzN03isrOGcPO5yTy/6gDL9+otApU6HRryasC4fHwc10yM58UFkwnwsX7J/PVloxkZE8S972yltKbRyRUqNfhoyKsBY0JiKA9/ewIh/t4dbX7envzzuxOoqGvmr5/tcWJ1Sg1OGvJqwBsVG8wPZwzj3Y2HWL+/1NnlKDWoaMirQeHOi9KID/XnNx/soKmlzdnlKDVoaMirQcHfx5M/XjGGfUU1vLDqgLPLUWrQ0JBXg8bs9BguTo/hn8uyOFhW1/cTlFIa8mpw+f0VY/AU4dfvb+fYLBpKqZ5oyKtBJT7Un19eOoqV+0r4z6bDzi5HqQFPQ14NOjdOSSYjOYw/fbyLouoGZ5ej1ICmIa8GHQ8P4YH546hvbuUPH+o9aZTqjYa8GpSGR1m4a1Yan24/ynNf7+9xvZbWNjbll1NZ33wGq1Nq4LBrgjKlnGnhjGHsOFzJXz7dTUFlPb/5VjqeHkJzaxtrckr5bPsRvthVSFltEzdNTeZPV411dslKnXEa8mrQ8vb04IkbzuYvn+7mhVUHOFhWR0SgL//ddZSKumYsvl5cNCqafUU1bDhQ5uxylXIKDXk1qHl4CL+dl05cqD9//mQXgT5ezEmP4bKzhnB+WiR+3p48tmwfjyzNorK+ucu8OEq5Aw155RJumz6UK8bHEeTnhZ+3Z5dlk5LDMAa2HKxg5ohj9xHefqiSvLJaJiaFERfih4ic6bKV6nca8splRAX5dts+PjEUD4FNeeUdIW+M4fbXNnK4wnrj8NhgP26elsyPL+jzPjhKDSoa8srlWXy9GBkbzKb88o62HYerOFxRz50XpRIe6MMXuwp58PO9nJ0UxtRhEU6sVinH0iGUyi1MSg5lc34FrW3WqRD+u/MoHgILzhvKgvOG8sItk0kM9+fX72/Xe8oql6Ihr9zCpOQwahpbyCqsBuCLXUc5Z2g44YE+gHWWy79cdRb7i2t58qucbl+jsr6Z5XuLOj4olBoMNOSVW5iUFA7AxrxyDpTUklVYwyVjYrusM2NEFFdNiOPJ5dlkF1Wf8Bp/+2w3C176hsv/bxWr9pVgjGHd/lJ+/PpGzvvblxwq15kx1cCjffLKLSSG+xNp8WFTfjk1jS0AXHxcyAP8Zl46y7OKuf/9Hby5cGrHiJvaxhY+3FLA2UmhFFU3cuML64kJ9qWwqpEQf29qGlt4dW0eiy4bfUa3S6m+6JG8cgsiwtlJYWzKK+e/O48yNj6Y+FD/E9aLtPhy9+wRrD9QxvpOF1B9sv0ItU2tLLpsNMt+PpNfXzaK0UOCefDacaxbNIu5Y2J5Y0M+dU0tZ3KzlOqThrxyG5OSw8gtrWNzfgWXpJ94FN/u25MTibT48sRX2R1tb39zkGFRgWQkh+Hr5cnCGcN5+fvncP3kRPx9PFlwXgpVDS18sLngpGrZlF9OVcPpzaezbHchH249ufdRSkNeuY1JyWEdv18ytueQ9/P25AfnD2XlvhK2Hqwgu6iGzLxyrs9I7PGCqYzkMMbEBfPymgN93swkp7iGa55cw/df+qbLSJ7y2iauemI1i97b3u3J3eqGZu59Zyu3vZLJnW9s5tW1uX1ssVIa8sqNjI0PwdtTGBoZSFq0pdd1b5yaTLCfF08uz+adzIN4egjXnB3f4/oiwoJpKWQV1rA2pxSAb3LLuOPfm064VeH7mw4jYj0J/LsPdmKMobaxhQUvf8P2w5W8sSGfX7+3nbZOQb82p5TLHlvJe5sO8ZMLU5k9OprfLt7Je5sOndLfQXltE29syNdhom5ET7wqt2E9Qh/G8ChLn1MYWHy9WHDeUB5bto8gP+tEZ9FBfr0+5/Lxcfz1sz08tSKHj7YV8MaGgwAEeHvy0HXjAWhrM7y/+TAz0qIYGx/ME1/lkBZjYUVWMdsPVfD0jZPYcbiSx77MxstTuHZSAo8syWLlvhISw/15+0fnkpESTkNzK7e+/A2/eHcbAT5ezO3lm0m7fYXV3PZKJvlldew5UsUfr3TcrJwP/XcPPp6e3DU7zWGv2ZPK+mY25pWRkRJOsN/Am4vIGMOKrGImJYcRNADq05BXbuWXc0ed9Lrfn5bC8yv3U93QwncmJ/a5vp+3Jzeck8TjX2Xj6SEsnDGMyrpm3tt8iJ9fPJLYED/WHSjlcEU9980dyeXj4th7tJo/f7IbgAevHcfFY2KZkx5DU6vh6RU5vL4+n/BAH3592ShumpqCv49nx3s9d3MGN72wnrvf2sLY+BkkhAX0WNvyvUX89N+b8fX25FtnDeGVtXmcOzyCuWOHdLv+wbI64kP98fDo+mH4duZBRsQEMSExtKPti51HecJ2bUFMsC/fOSepY9kbG/IprWnkJxedevi3tRm2Ha6krrGFljZDRX0zn20/wrLdRTS1tjE0MpBnb5pEWkzQKb+2vcprm8gqrKayvpmLRkXj5XmsU+SZr/fzt8/2MCzKWl9qtLW+5tY2lu0uYmJSKDHBvR8wOJJdIS8iDwGXA01ADvB9Y0yFiKQAu4G9tlXXGWNut+e9lDrTwgJ9+OH5w/hk+5EuE5v15rbpQ6lpbOG6jATGxIVwsKyOdzYe5KU1B1h06Wj+s/EwQb5eXDImFg8P4ZFvT+Cnb2zmghFRXG/7IBERfjl3JBGBPrQaw01Tkwn0PfG/aqCvF/93w9nM/scK/vjRLp67OaPbmr7YeZTbX9vIyNhgnr8lgyiLL4fK6/jFu9sYExdCYrj1w6G1zbBkVyEvrNrPN7nl/PiC4dzX6UPxm9wy7nt3GxZfL9760VTGxIVQWdfMbz7YwajYIKKCfPnd4p2MjA1ifEIoD3y+h2dsN3SZNTqG0UOCT+rv0BjDf3cW8ujSLPYc7Xq9QqTFh+9NTWJCYih/+ng3Vz2xmoe/PYE5o2OoqG+mrLaRISH+3f59OcLXWcXc9+42jlYdu+3kJWNieOy7E/H18mR1dgkPfr6H81Ij2Hu0misfX81frx1HaU0jz329n4LKBmKD/Xjl1nMYGXtmPpzEnjvei8jFwJfGmBYReQDAGPNLW8h/bIw5pe+DGRkZJjMz87TrUWog+sm/N7FibzFLfz6TC/++nMvHxfHA/HEOe/1nVuTw18/28NzNGcxJj+myrKi6gUse+Zr4MH/eWnhuR/gdLKvjssdWkhwRwPTUKPYX17DjcCUFlQ0khPkTG+zH1kMVfHbX+aRGB9HWZrj6ydUcrWrAU4TmNsN7/zONx5bt473Nh1l8x3kkhPlz+eOraG4xTB4azkdbC/h2RiIfbStg7phYHv72hD63Jbuohrve3MzOgiqGRQZy+wXDSQoPwMtD8PXyZPSQoI6j5iOV9fzo1Y1sO1SJp4d0nKz29hQmJoUxPTWS+ZMSiDtuqOzq7BLW7y8lwNeLQB9Pgv29ibT4EhXkS3xozx8QOcU1XPXEamKD/bg+I5ERsUFkHa3mL5/u5vy0SH5/+Riuf2YtEYE+fHDHeVQ1NHP7a5vYerACgMkpYcyflMDDS7Kob2rluZszmOKgeZJEZKMxpttPebtC/rg3uRqYb4z5noa8UsdsP1TJ5Y+vYnxiKFsPVvD2j87lnKHhDnv95tY25j22iprGFpbcM4MAH2tIGWO47ZVMVmeX8Mmd0zu6Ddp9vuMI//P6Jrw8hOSIQFKjLFwxIY6L061HxRf9fTlj4kL49w+n8MGWw9z91lb+cd14xiWEMP/ptfh5e1BY1chPLkzl3ktGArCroIprnlpNQ3MbP58zgp9clMqfPt7Nv9bmsuK+C7u9NqFdXmkt1z+zltY2w6JLR3PlhLgu3SDdaWhu5YVVB6hraiHS4ktYgA97C6tZua+YnQVVRFp8ee22KR1HzR9uLeBnb26mp5kpfDw9uGBkFFdOiGfW6OiOaaurGpq56onVVNY18+FPp3fZjrczD/Kr/2zDQwR/b08++Ml5DI+ynthvbGnl1bV5jE8MZXKKdZ8fKq/jlhc3cLC8nmdumsSFI6N73caTcaZC/iPgLWPMa7aQ3wlkAVXAb4wxK3t43kJgIUBSUtKkvLw8h9Sj1EDyvefXsTq7lMRwf77+xYUOn7v+m9wyrnt6Lbecm8yiy0bj5+3Jmxvy+dV72/ndvHRunT602+eV1TYR7OfVbZi+ui6P336wgwevHccjS7OItPiy+I7z8PAQvskt43vPryc5PICP75yOr9exOfzX7y+lsr6544riwxX1zHzwK26ZlsJv56V3W8fhinquf3otdU0tvLnwXId0ZWQVVnPTC+tpbGnj5e+fQ35ZHT97czOTU8J5ccFkPESobWqhoq6ZkppGiqsb2ZRfzsfbjlBc3Yi/tyczRkQye3QMn24/wsp9Jbz+gyndHn1/su0Iv/9wJ3+95qwTvk11p7y2ie8+t46SmkaW3jOT0AAfu7bVrpAXkaVAd6fu7zfGLLatcz+QAVxjjDEi4gtYjDGlIjIJ+AAYY4yp6u299EheuaoVWcXc8uIG7pqVxt1zRvTLeyx6bxtvbDiIj5cHk1PC2JJfwfjEUF67bcoJJ1BPRquti2b74UqM4YRvINlF1YQG+BBp6X4e/87ueWsLn+88ytpfzSIkwJuSmkYyc8uoamihtrGFV9bkUlrbxBs/nMrY+JBTrrUnB8vq+N7z6ympaaShuZXJKeG89P3JHd92utPaZli/v5RPdxxh6a6ijv73P101lpumJvf4PGPMKX147yqo4vLHV3HNxPiO0Venq1+P5EVkAfAjYJYxptsZmkRkOXCvMabXBNeQV67KGMMXuwo5Py2y14CxR2ubYeW+YlbtK2HlvhLK65p4/47zeu0i6cvWgxVc9eRqLkmP5embJp326+w5WsXcR1cyf1IC9U2tfLHrKM2tx7InxN+bFxdM7nLBmqMUVjWw4KVviLT48MxNk07p798Yw86CKo5UNjB7dLTDv4E98Pkenlqew79/MIVpqZGn/Tr9FvIiMhd4GJhpjCnu1B4FlBljWkVkGLASOMsY0+vdlDXklRp4dhyuZFhUoN0fTgte2sDyvcWEBnhz7dkJXD4+johAHyy+Xlj8vPDuo//dHu05N9Bu8djQ3MrcR7/GAP/92YwTbl15svoz5LMBX6DU1rTOGHO7iFwL/D+gGWgDfm+M+aiv19OQV8p1FVY1dNxn93TDzBWtySnhhufWc/vM4fzq0pO/jqOz3kLero9mY0y3N8Q0xvwH+I89r62Uci0xwX4nzOGvYNrwSBZMSyE+tH8ukNIrXpVSysn+cMWYfnttnaBMKaVcmIa8Ukq5MA15pZRyYRrySinlwjTklVLKhWnIK6WUC9OQV0opF6Yhr5RSLsxhUw07gogUA/bMNRwJlDionMHCHbcZ3HO7dZvdx6lud7Ixptvblw2okLeXiGT2NH+Dq3LHbQb33G7dZvfhyO3W7hqllHJhGvJKKeXCXC3kn3V2AU7gjtsM7rndus3uw2Hb7VJ98koppbpytSN5pZRSnWjIK6WUC3OJkBeRuSKyV0SyReRXzq6nP4hIooh8JSK7RGSniNxlaw8XkSUiss/2p+PvhDwAiIiniGwWkY9tj4eKyHrbPn9LRHycXaMjiUioiLwrIntEZLeInOsO+1pE7rb9+94hIm+IiJ8r7msReVFEikRkR6e2bvevWD1m2/5tInL2qbzXoA95EfEEngAuBdKB74pIunOr6hctwM+NMenAVOAO23b+ClhmjEkDltkeu6K7gN2dHj8APGK7BWU5cJtTquo//wQ+N8aMAsZj3XaX3tciEg/cCWQYY8YCnsB3cM19/TIw97i2nvbvpUCa7Wch8NSpvNGgD3ngHCDbGLPfGNMEvAlc6eSaHM4Yc8QYs8n2ezXW//TxWLf1FdtqrwBXOaXAfiQiCcC3gOdtjwW4CHjXtopLbbeIhAAzgBcAjDFNxpgK3GBfY70lqb+IeAEBwBFccF8bY74Gyo5r7mn/Xgn8y1itA0JFZMjJvpcrhHw8cLDT40O2NpclIinARGA9EGOMOWJbdBSIcVZd/ehR4D6gzfY4AqgwxrTYHrvaPh8KFAMv2bqonheRQFx8XxtjDgN/B/KxhnslsBHX3ted9bR/7co4Vwh5tyIiFuA/wM+MMVWdlxnreFiXGhMrIvOAImPMRmfXcgZ5AWcDTxljJgK1HNc146L7OgzrUetQIA4I5MQuDbfgyP3rCiF/GEjs9DjB1uZyRMQba8C/box5z9Zc2P7VzfZnkbPq6yfnAVeISC7WrriLsPZXh9q+0oPr7fNDwCFjzHrb43exhr6r7+vZwAFjTLExphl4D+v+d+V93VlP+9eujHOFkP8GSLOdgffBeqLmQyfX5HC2fugXgN3GmIc7LfoQuMX2+y3A4jNdW38yxiwyxiQYY1Kw7tsvjTHfA74C5ttWc6ntNsYcBQ6KyEhb0yxgFy6+r7F200wVkQDbv/f27XbZfX2cnvbvh8DNtlE2U4HKTt06fTPGDPof4DIgC8gB7nd2Pf20jdOxfn3bBmyx/VyGtX96GbAPWAqEO7vWfvw7uAD42Pb7MGADkA28A/g6uz4Hb+sEINO2vz8AwtxhXwN/BPYAO4BXAV9X3NfAG1jPOzRj/eZ2W0/7FxCsIwhzgO1YRx+d9HvptAZKKeXCXKG7RimlVA805JVSyoVpyCullAvTkFdKKRemIa+UUi5MQ14ppVyYhrxSSrmw/w+SgQztsxvrqwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We now can plot the energy during those optimisation steps:\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(energy_history)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a10ade91-7ad7-4a8c-997d-e4ee9098b2d7",
   "metadata": {},
   "source": [
    "### 4b. Use NetKet's optimisation driver \n",
    "\n",
    "As writing the whole optimisation loop by yourself every time is.. boring, we can make use of a coupled of NetKet's built-in utilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4572b5c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No output specified (out=[apath|nk.logging.JsonLogger(...)]).Running the optimization but not saving the output.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:04<00:00, 72.51it/s, Energy=-25.056 ± 0.053 [σ²=1.435, R̂=1.0001]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized energy and relative error:  -25.018 ± 0.048 [σ²=1.171, R̂=1.0089] 0.018562865344048777\n"
     ]
    }
   ],
   "source": [
    "# First we reset the parameters to run the optimisation again\n",
    "vstate.init_parameters()\n",
    "\n",
    "# Then we create an optimiser from the standard library.\n",
    "# You can also use optax.\n",
    "optimizer = nk.optimizer.Sgd(learning_rate=0.05)\n",
    "\n",
    "# build the optimisation driver\n",
    "gs = nk.driver.VMC(H, optimizer, variational_state=vstate)\n",
    "\n",
    "# run the driver for 300 iterations. This will display a progress bar\n",
    "# by default.\n",
    "gs.run(n_iter=300)\n",
    "\n",
    "mf_energy=vstate.expect(H)\n",
    "error=abs((mf_energy.mean-eig_vals[0])/eig_vals[0])\n",
    "print(\"Optimized energy and relative error: \",mf_energy,error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4f4b1a2d-fcbf-4248-ae18-d0f527ac67e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final optimized parameter:  [2.63226906]\n"
     ]
    }
   ],
   "source": [
    "# we can also inspect the parameter:\n",
    "print(\"Final optimized parameter: \",vstate.parameters[\"lambda\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5bb5976b-da47-45af-88e8-690afb3980b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-25.49098968636478"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eig_vals[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9f0e656",
   "metadata": {},
   "source": [
    "## 5. Jastrow Ansatz\n",
    "\n",
    "We have seen that the mean field ansatz yields about 2% error on the ground-state energy. Let's now try to do better, using a more correlated ansatz. \n",
    "\n",
    "We will now take a short-range Jastrow ansatz, entangling nearest and next-to nearest neighbors, of the form \n",
    "\n",
    "$$ \\langle \\sigma^{z}_1,\\dots \\sigma^{z}_N| \\Psi_{\\mathrm{jas}} \\rangle = \\exp \\left( \\sum_i J_1 \\sigma^{z}_i\\sigma^{z}_{i+1} + J_2 \\sigma^{z}_i\\sigma^{z}_{i+2} \\right),$$\n",
    "\n",
    "where the parameters $J_1$ and $J_2$ are to be learned. \n",
    "\n",
    "Again we can write the model using flax. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4b390e0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "class JasShort(nn.Module):\n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "        \n",
    "        # Define the two variational parameters J1 and J2\n",
    "        j1 = self.param(\n",
    "            \"j1\", nn.initializers.normal(), (1,), float\n",
    "        )\n",
    "        j2 =self.param(\n",
    "            \"j2\", nn.initializers.normal(), (1,), float\n",
    "        )\n",
    "\n",
    "        # compute the nearest-neighbor correlations\n",
    "        corr1=x*jnp.roll(x,-1,axis=-1)\n",
    "        corr2=x*jnp.roll(x,-2,axis=-1)\n",
    "\n",
    "        # sum the output\n",
    "        return jnp.sum(j1*corr1+j2*corr2,axis=-1)\n",
    "    \n",
    "model=JasShort()\n",
    "\n",
    "vstate = nk.vqs.MCState(sampler, model, n_samples=1008)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1616acc8",
   "metadata": {},
   "source": [
    "We then optimize it, however this time we also introduce a stochastic reconfiguration (natural gradient) preconditioner. Also, we now log the intermediate results of the optimization, so that we can visualize them at a later stage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "11db148b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:04<00:00, 70.02it/s, Energy=-25.298 ± 0.019 [σ²=0.374, R̂=1.0126]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final optimized parameters: j1=[0.23400646], j2=[0.08184639]\n",
      "Optimized energy : -25.337 ± 0.020 [σ²=0.418, R̂=1.0026]\n",
      "relative error   : 0.006027448868855777\n"
     ]
    }
   ],
   "source": [
    "optimizer = nk.optimizer.Sgd(learning_rate=0.05)\n",
    "\n",
    "gs = nk.driver.VMC(H, optimizer, variational_state=vstate, preconditioner=nk.optimizer.SR(diag_shift=0.1))\n",
    "\n",
    "# Loggers that work together with optimisation drivers are defined in nk.logging.\n",
    "# RuntimeLog keeps the metrics in memory, JsonLog stores them to a json file which can be read\n",
    "# as a dict, TensorBoardLog can be used to log to TensorBoard.\n",
    "log=nk.logging.RuntimeLog()\n",
    "\n",
    "# One or more logger objects must be passed to the keyword argument `out`.\n",
    "gs.run(n_iter=300, out=log)\n",
    "\n",
    "print(f\"Final optimized parameters: j1={vstate.parameters['j1']}, j2={vstate.parameters['j2']}\")\n",
    "\n",
    "jas_energy=vstate.expect(H)\n",
    "error=abs((jas_energy.mean-eig_vals[0])/eig_vals[0])\n",
    "print(f\"Optimized energy : {jas_energy}\")\n",
    "print(f\"relative error   : {error}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9b18a6a",
   "metadata": {},
   "source": [
    "You can now see that this ansatz is almost one order of magnitude more accurate than the mean field! \n",
    "\n",
    "In order to visualize what happened during the optimization, we can use the data that has been stored by the logger. There are several available loggers in NetKet, here we have just used a simple one that stores the intermediate results as values in a dictionary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4b34166d-92eb-4ad8-a43f-cffc3bb717be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Energy': History(\n",
      "   keys  = ['Mean', 'Variance', 'Sigma', 'R_hat', 'TauCorr'], \n",
      "   iters = [0, 1, ... 298, 299] (300 steps),\n",
      ")}\n"
     ]
    }
   ],
   "source": [
    "data_jastrow = log.data\n",
    "print(data_jastrow)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00670bd1",
   "metadata": {},
   "source": [
    "These report several intermediate quantities, that can be easily plotted. For example we can plot the value of the energy (with its error bar) at each optimization step.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "6d5c22fb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmDklEQVR4nO3deZxcZZn28d9de3d6SUISEkhCwCABQkBsEFkCYVdRmFEHXkYFYURUxm2QF4yjuKBiQHQGFeMyOrworggzsiSBhAAaIEAIIUAMhNUsna33pZb7/eOc6lR3dXc6Id3Vnbq+n08lVefUcp861XXV8zxnMXdHRESkUKTUBYiIyPCjcBARkSIKBxERKaJwEBGRIgoHEREpEit1AXvCuHHjfNq0aaUuQ0RkRHniiSc2u/v43ubtFeEwbdo0li9fXuoyRERGFDN7pa956lYSEZEiCgcRESmicBARkSIKBxERKaJwEBGRIgoHEREponAQEZEiJQkHM5tnZs+b2Uozu8PMRhfMu8bM1prZC2Z2VinqExEpd6VqOSwEZrr7LGANcA2AmR0GXAAcDpwN/NDMooNZyPk//ivn//ivg/kSIiIjTknCwd0XuHsmvLkMmBxePxe43d073H0dsBY4thQ1ioiUs+Ew5nAJcE94fX/gtYJ5r4fTipjZZWa23MyW19fXD3KJIiLlZdCOrWRmi4CJvcya6+53hveZC2SA23b1+d19PjAfoK6uTuc6FRHZgwYtHNz99P7mm9nFwDnAab7jRNZvAFMK7jY5nDYoXt/WystbWphQnRyslxARGZFKtbXS2cBVwPvcvbVg1l3ABWaWNLMDgYOBxwarjq0tnWxs7KAjkxuslxARGZFKdcjum4EksNDMAJa5++Xu/qyZ/RZYTdDd9Cl3zw5WEdGIAeDqlBIR6aYk4eDu0/uZdx1w3VDUEYsEDSdlg4hId8Nha6WSybcc1HQQEemurMMhlu9WKnEdIiLDTVmHQ77lcPnJbylxJSIiw0tZh0MsGoRDNqe2g4hIobIOh3zLIa1wEBHppqzDIb+1Ujar/RxERAqVdziE3UoZtRxERLop73CIaMxBRKQ3ZR0O+TEHtRxERLor63DoGnNQOIiIdFPW4ZDfQVotBxGR7so6HMyMWMTI5rS1kohIobIOBwjGHdRyEBHpruzDIRYxslmFg4hIobIPB7UcRESKlX04xKIRMhpzEBHppuzDIRoxbcoqItJD2YdDPGJkNOYgItJN2YdDNKqWg4hIT2UfDrFIRAPSIiI9lH04aMxBRKRY2YdDLGLaWklEpIeyDwe1HEREipV9OMS0E5yISJGyDwe1HEREipV9OMQiEe3nICLSQ9mHQ1QD0iIiRUoSDmY2z8yeN7OVZnaHmY0Op+9jZovNrNnMbh6KWmJRjTmIiPRUqpbDQmCmu88C1gDXhNPbgX8HrhyqQmIacxARKVKScHD3Be6eCW8uAyaH01vc/WGCkBgSUY05iIgUGQ5jDpcA9+zqg8zsMjNbbmbL6+vrd/vF1XIQESkWG6wnNrNFwMReZs119zvD+8wFMsBtu/r87j4fmA9QV1e329/u0agGpEVEehq0cHD30/ubb2YXA+cAp7l7yX66q+UgIlJs0MKhP2Z2NnAVcLK7t5aihjydJlREpFhJwgG4GUgCC80MYJm7Xw5gZi8DNUDCzM4DznT31YNViFoOIiLFShIO7j69n3nThrCUYGslhYOISDfDYWulklLLQUSkWNmHQzRipLPaWklEpFDZh4NaDiIixco+HKI6tpKISJGyD4d4JKKWg4hID2UfDvmT/ZRwPzwRkWGn7MMhFjEAtR5ERAqUfThEo0E4aNxBRGSHsg8HtRxERIqVfThEI8FboJaDiMgOZR8OajmIiBQr+3CIRvJjDtpLWkQkr+zDQS0HEZFiZR8OXS0HnUdaRKRL2YdDTJuyiogUUTiEWytlNeYgItJF4RBRy0FEpKeyDweNOYiIFCv7cMiPOWhrJRGRHco+HLSHtIhIsbIPB+3nICJSrOzDIT/m8OU7V5W4EhGR4aPswyHfchARkR3KPhzyLQedCU5EZIeyD4f8TnDKBhGRHRQO4aasygYRkR3KPhwSseAtyKnpICLSReEQVbeSiEhPJQkHM5tnZs+b2Uozu8PMRofTzzCzJ8zsmfD/Uwe7FrUcRESKlarlsBCY6e6zgDXANeH0zcB73f0I4CLg1sEuJN9yuOj4aYP9UiIiI0ZJwsHdF7h7Jry5DJgcTn/K3f8eTn8WqDCz5GDWEg9bDumsDtktIpI3HMYcLgHu6WX6+4En3b2jtweZ2WVmttzMltfX1+/2i+dbDp0ZhYOISF5ssJ7YzBYBE3uZNdfd7wzvMxfIALf1eOzhwPXAmX09v7vPB+YD1NXV7faAQTzclFXhICKyw6CFg7uf3t98M7sYOAc4zQt2TzazycAdwEfc/cXBqq/g9UhEI3TqfA4iIl0GLRz6Y2ZnA1cBJ7t7a8H00cCfgavd/ZGhqicRi6jlICJSoFRjDjcD1cBCM1thZreE068ApgNfDqevMLMJg11MIhbRgLSISIGStBzcfXof078BfGOIyyEeNbUcREQKDIetlUouEYvQqZaDiEgXhQOEA9IKBxGRPIUDEI9qQFpEpNCAwsHMbgz3PdgrJbW1kohINwNtOTwHzDezR83scjOrHcyihlo8qq2VREQKDSgc3P2n7n4C8BFgGrDSzH5lZnMGs7ihov0cRES6G/CYg5lFgRnhZTPwNPB5M7t9kGobMtpaSUSkuwHt52BmNwHvBe4Hvunuj4WzrjezFwaruKGiAWkRke4GuhPcSuBL7t7Sy7xj92A9JaGWg4hIdwMNh6eBQ8yscFoD8Iq7N+zxqoZYUgPSIiLdDDQcfggcTdCCMGAmwcl4as3sE+6+YJDqGxLqVhIR6W6gA9J/B97m7nXu/nbgbcBLwBnAdwaruKGirZVERLobaDi81d2fzd9w99XADHd/aXDKGlrBfg46n4OISN5Au5VWm9mPgPxmq+eH05JAelAqG0JqOYiIdDfQlsNFwFrgs+HlJeBigmAY8TvC5bdWKjghnYhIWdtpyyHc+e1ud58D3NjLXZr3eFVDLBGeRzqddRIx28m9RUT2fjttObh7FsjtbcdTKpSIBW+D9nUQEQkMdMyhGXjGzBYCXTvCufunB6WqIZaIBuGQzuQgWeJiRESGgYGGwx/Dy14prpaDiEg3AwoHd/+lmVUAU919xB9Lqad8y0FbLImIBAZ6sp/3AiuAe8PbR5nZXYNY15DSmIOISHcD3ZT1WoID7G0HcPcVwEGDUlEJqOUgItLdQMMh3csB9vaab9KuloPCQUQEGPiA9LNmdiEQNbODgU8Dfxm8soZWPL+1krqVRESAgbcc/hU4HOgAfg00EuwpvVfItxy+9KdVJa5ERGR4GOjWSq3A3PCy18mHgw6fISISGOhpQt8KXAlMK3yMu586OGUNrfyAdE7ZICICDHzM4XfALcBPgeybfVEzm0dwTupO4EXgo+6+3cyOBebn7wZc6+53vNnX2xm1HEREuhtoOGTc/Ud78HUXAte4e8bMrgeuAf4vsAqoC6dPAp42s/9x98wefO0iqVgUUMtBRCRvoAPS/2NmnzSzSWY2Nn/Z3Rd19wUFX/jLgMnh9NaC6SlgSL6uU/F8t5LSQUQEBt5yuCj8/wsF05w9syPcJcBv8jfM7B3Az4EDgA/31Wows8uAywCmTp36pgpIJYKWw4eOO+BNPY+IyN5ioFsrHbirT2xmi4CJvcya6+53hveZC2SA2wpe61HgcDM7FPilmd3j7u291DSfcHyirq7uTf3kz3crtaff9HCKiMheod9uJTO7quD6B3vM+2Z/j3X30919Zi+XfDBcDJwD/LP3MhLs7s8RHCp85oCXZjfFo0bEoD2tneBERGDnYw4XFFy/pse8s3f3Rc3sbOAq4H3hPhT56QeaWSy8fgAwA3h5d19nF+ohFY+q5SAiEtpZt5L1cb2327viZoLT6iw0M4Bl7n45cCJwtZmlCY7d9El33/wmXmfAUvEo7RmFg4gI7DwcvI/rvd0eMHef3sf0W4Fbd/d534xULKJuJRGR0M7C4UgzayRoJVSE1wlvpwa1siGWikdpU7eSiAiwk3Bw9+hQFVJqyXiUDoWDiAgw8J3g9nqpuLqVRETyFA6hVExbK4mI5CkcQhUJba0kIpKncAipW0lEZAeFQ0jdSiIiOygcQsl4VC0HEZGQwiGUike0KauISEjhENJOcCIiOygcQqlYlEzOyWTVtSQionAI5c8G155ROIiIKBxCqbhO+CMikqdwCHW1HBQOIiIKh7wdLQd1K4mIKBxC6lYSEdlB4RDKh0OHjq8kIqJwyEvFgrdi7h2rSlyJiEjpKRxC+ZZDznf77KciInsNhUOoKxxyCgcREYVDqCoVnDE1q2wQEVE45NXkwyGnTVlFRBQOoVGJIBwy6lYSEVE45EUiRm1FnDMO3bfUpYiIlJzCoUBNRYzG9kypyxARKTmFQ4GaVJzGtnSpyxARKTmFQ4GaVJzGdoWDiEhJwsHM5pnZ82a20szuMLPRPeZPNbNmM7tyKOuqqYjR2KZuJRGRUrUcFgIz3X0WsAa4psf87wL3DHVRajmIiARKEg7uvsDd8z/RlwGT8/PM7DxgHfDsUNdVUxGnQWMOIiLDYszhEsJWgplVAf8X+OrOHmRml5nZcjNbXl9fv0cKqUnFae3MktZ5pEWkzA1aOJjZIjNb1cvl3IL7zAUywG3hpGuBm9y9eWfP7+7z3b3O3evGjx+/R2quqQh2hGvS5qwiUuZig/XE7n56f/PN7GLgHOA0965Dob4D+ICZfQcYDeTMrN3dbx6sOgvVpOIANLalGTsqMRQvKSIyLA1aOPTHzM4GrgJOdvfW/HR3P6ngPtcCzUMVDBCMOQAalBaRsleqMYebgWpgoZmtMLNbSlRHN7VhOFz1+5UlrkREpLRK0nJw9+kDuM+1Q1BKN/kxBx18T0TK3XDYWmnYyI8zpDPaWklEypvCocD4qiQRg45MttSliIiUVEm6lYYrM2P6hCoO2GdUqUsRESkptRx6mDq2kte2tu78jiIiezGFQw+TxwThsGPXCxGR8qNw6GHq2EpaOrNsbeksdSkiIiWjcOhh6thKAD7ys8dKXImISOkoHHqYEoZDu7ZYEpEypnDoYb/RKQA6ta+DiJQxhUMPVckYlYkoZxw2sdSliIiUjMKhBzNjYk2KjU3tpS5FRKRkFA69mFCTZGODwkFEypfCoRf7quUgImVO4dCLiTUpNjZ2aEc4ESlbCodeTKhJ0ZnJsb1VJ/0RkfKkcOjFvjVJAHUtiUjZUjj0YmJNsK/DZ379VIkrEREpDYVDL/YNw6Ezqx3hRKQ8KRx6sd/oCsygrVOH0BCR8qST/fQiGjGO2L+W6pTeHhEpT2o59GHGxGqeW9+kzVlFpCwpHPpw6KQatrZ0Ut/UUepSRESGnMKhDzMm1gDw4Z89WuJKRESGnsKhD4dOqgagVYPSIlKGFA59GF2ZIBGNKBxEpCwpHPpxwvR9GF+dLHUZIiJDTuHQjxmTali7qVlnhRORsqNw6Mehk2rI5Jy1m5pLXYqIyJAqSTiY2Twze97MVprZHWY2Opw+zczazGxFeLmlFPXlHToxGJT+1K+eLGUZIiJDrlQth4XATHefBawBrimY96K7HxVeLi9NeYGDxleRikVY39CmriURKSslCQd3X+DumfDmMmByKerYmWjE+NGH3k57OscZ332w1OWIiAyZ4TDmcAlwT8HtA83sKTN70MxO6utBZnaZmS03s+X19fWDVtycGROoTsVY39Cu1oOIlI1BCwczW2Rmq3q5nFtwn7lABrgtnLQemOrubwM+D/zKzGp6e353n+/ude5eN378+MFaDABuvvBoOrM5zvre0kF9HRGR4WLQDjvq7qf3N9/MLgbOAU7z8Oh27t4BdITXnzCzF4G3AssHq86BmH3wOCoTUdZvbyObc6IRK2U5IiKDrlRbK50NXAW8z91bC6aPN7NoeP0g4GDgpVLUWMjMuOGDR9KeyTH9i3dz8rzFpS5JRGRQlWrM4WagGljYY5PV2cBKM1sB/B643N23lqjGbs46fCLjq5KYwWtbW1nf0FbqkkREBo3tDecrqKur8+XLh6bn6bWtrcz+zmL2rUmy7Iv99pyJiAxrZvaEu9f1Nm84bK00okwZW8mcGROIxyI8tm4rx163iLM1UC0iexmFw2447dAJvLa1jX/68V/Z1NTBus0tbG/t7HdT14a2NA1t6SGscvdsae4gmxuZrcnOTI458xbz/h8+UupSRpQ5NyweUfvx1Dd1kBuhn9GRROGwG06bsS8A8WgwUN2RyXHU1xbyjm8uwt35xSPr+PnD62juCPbzc3eO//b9nPDt+/t93nWbW/jLi5t3uZ6efyjZnLN0TX1XWDW1p6n7xkLO+8HDXfdpKzgUeT4MGtrSHHPdIk68/oFdrqFnPX0FZVN7mjUbm7pN21nXZiZb/FxvbG8retxPHnqJdVtaeX17/+NBT726jYbWgQX1mTc9yFk37f4XZ0Nbuuu9aE9n2drS2TUvm3Ne2dJS9Jjmjgw3LVzDD5es7bae8pra0wM6fe3iFzbxvUVrOHju3Rx0zZ8566alRY/b0NDOy5tbeWVrK2//+kKOuW5R0fOc/t0HOanHZ+I9//EQc25Y0m2au9PUvuN9TWdzvLAhWNfv+t5SZvz7PWxsbO+a39aZ5clXt/HnleuL1nFf6yffWj/mukU7fQ9OmReEXi7nvO1rCzjyqwvY0NBOU3uaOTcsKVmLv6Ujw68efZVtBZ+Fts7ssDsl8aBtyro3m1ib4ivvPYyjp45h1uRarvvzajoyOba1ppl57X20dAR/0PPue4HRlXEy2VzXtLWbmmjuyNLcnqEyGeXJV7YxeUwl3/jf1V1fapNHV7ChsR0z2GdUgmgkwpbmDm758NtZu6mZiBl3rniDzc0dZHNQ39zBHz5xPLP2r+W5DY188rYneWVLK7PfOp6/bWwiFjE2N3cSjwa/BZa8sImP/tfj7FuT5E+fOpETr3+AqmSMy04+iJxDY1uabS2dLHtpC7UVcd75ln1o7sgQMaMyEWXl6w1c/v+eoCoZ44YPHkk6m6Nu2lia2tO8VN/C//nJMgx4+itncsoNS9jemuaGDx7JWYfvy0nfWUxDa5rx1Una01lyHnwZTqhOUpWMsaGxnQWfm83kMZUAnHrDEtZtbqGmIk5rZ4ZfXnIsV/1uJa9vb+O0GRNYvb6R8dVJ5n+4jhsXvIAZbGtN8/DfNnPl71aQiEW477Mn09KZ4czvPkgkfC9S8QiLrzyFNRub+eIfV1JbkeDuzwT7XJ7+3QdJRI1ffew41m5qJufwli/ezUHjRnHrpe8g686lv3icTM5Z+LnZtKWztKdzPP7yVpKxCKccMgGAE759P29sbycWMWor4mwJvwwm1aZ46Ko5HPvN+9na0smk2hSNbWl+etExvP2AMbz7+0t5dWvwWdjS3MmHjzuAdZtbOPqAMbywoYnzf/xXIhFjypgKfvHRY6lv7mDW5FoS0Qgnz1vCG9vbGF+VZEP4RVyVjJHJ5nhpczOHffk+KhNRHvi3U6ipiDF/6Us4QatrS6aTRCxCS0eGUckY7/zW/XRmcjS2p0lnnV8/9iqxiFHf3MFz6xvJOfx2+Wv8U90UGlrTnDxvMQ1twbqNRYy2dJZtrWluOv9IqivitKdz3LbsFS498SA+eMtfWFNwQMupYyu7Pu8fO+kgPnHbk0yqTbH0qjnc/vhrJGMRjpw8mgt/sgwz2NLSyak3LuGez8wmGYtgFmxe/r8r/87UsZVB6G0JNoR8+zcWsq01jRGEvRk0tGUwg1nX3kdTR4aJNSk+NWc6qXiU847aj2jEePzlbTS0pTn90An8bVMzOXeu/sNKtrWkScQirNu8I9gn1aa473OzWf7yNsaOSnD4fjWYGe7OA89v4po/PkN1Ksb5x0zhxgVr6Mjk+Pr/rub2y46jMhHlrO8txYCZ+9fym4+/kzO++yC1FXE2N3cwujLBvZ+d3fVaje1pzrppKVXJGDdfeDQ1FTEm1Vbske+2QhqQ3kM6MzmO+toC2jqzTBlbwdfOnclHf/E4+bc3v2/EnEMm8MDzG3GHSMS6deHUpGK0pbOks05+T4q+1k7EIP9QM0hEI1QlY2xp6SQWMcaOSrCpx/mvJ4+p4DcffycnXf8AETMyOScRi5DLOTl3cg5jKuNs6/GrzQx6+5hEIwYOWffu9YR1f+KUt3DLkhe7lqEmFaOxPdM1v2fdeRXxKLdeeiyfvX0Ff29oIxqx4P0yozP8hVmVjNLamcU9eK7aijgNbWlmTKzmbxubyRYUPK4qQXNHhvZ08NjqZIyWzgw9eyamT6ji1kuP5fhvPUAsalxywoH8eOlLjKtKsKW5EzOIRSLhe+W4Q3W4TIViEWN8dZL6pg5OOngcj67bSmtnlgnVSXLubG7upCoZo7kjQywSrIfC9Tp2VIID9hnFIROr+dWjr3ZfF0AsGrwfhY8bX5XgC2fN4Ko/rOx6f1PxCMuuOY3RlQlWvr6d9938SK/raUxlnO2t6a719IMLj2bymArO/cGO7rloj89qRTxKLGo0t2e6lhWK34941Ehng8eNSkRpKWgJHbF/LVecOp0fLl7L0683dKsp/3q9ffZm7V/Lus0ttKWzeFjL9y84inWbW/jGn5/rul9lIkoyFmFba5qaVIyLjp/GDxavxYHz66bwuydeJ5tzqlMxmnqsw0KVieCzZgSfs+1h97B1/RPUWLisEQtCORWPsqmpo2u58vOmjKnk9XC/qVjEqEhEaevMdi1z4WczYsGOuCcdPI53fe+hoNVM8PcTMaMqGePpr5zZZ/396W9AWuEwiFa90UDEjJc2N1OdivPZ259iW2uaRDSC42Rzzq8/dhyViRhTxlYwujLBu76/lMa2NIuvnEMiFuGsm5bS2J7m1kvfwZIXNvHuIybxjz/8Cxsb2/nRh45m7Kgk1961itXrg+b7qGSU1o4sv77sOI7Yv5av/s+z/Hb56yRjEToyua75R06u5egDxvLzR9YxvjrJvA/M4uL/epz9R1fwpfccyr/fuYqaVJwrTp3OtXc9SzIW/KFtbGpn8pgKrvuHI7jwJ49iBL+aNjS2U52K8/XzZnLE/rWce/PDXUHw4Bfm8KGfPcqrW1uJGCz43Gy2taY5ZtrYrvfq+G/dz6amDqaMrez2iywRjbDo8yczdZ9K/uEHj7CpqZ1PzpnOB98+hUQswqamdmZfv5j2TI7TD53ATy86hs3NHSxcvZETp4/j5gfW8pvlrwHBUXZ/dvExTKhO8vKWFs77wV+oTsW464oTmXPDEpo7MlSnYrSGf6QAoyvirPjKmazb3NLVjVK4D2TO4eLjpzF5TAWHTarhyt89zfqG9q4vuEe/eBo1qeAX4H6jK8jmnFnX3kdLZ5YL3zGVZ17fzoaGdn53+fE89Ld6vnXP87R2Zpn3gVmcNXMiZ9+0lGQswtfOm8m//fZp6ps7OGTfau75zEm8vKWVxc9vYvX6Rn7/xOtA8OX8xJfO6Prlede/nthV66LVG5k8toLNTZ2s+nsDtzz4IhXxKEuvmsPx336AiAUtlXzwRSNGVSpGNuv86YoTWPVGA0dNGc0+VUlGJaI0dWQ46fqgtVCTivGrjx3HzP1rcXde2tzCU69u54xD9+Vd/7GUjnSOcVUJ1m1uZeyoOOmsc/dnTmLfmhSLVm/kX/57OeOqErSls3RkchyxXy2nzpjAzYvXMqk2xY3/dCSX/mI5+9YEv9CfeGUr7//RX4lHjZzv6BqtTsWIGMSjEe797GzGVSVp6cgQj0ZYs7GJc/4z6Fq9+9MnkcnlaG7PcPz0cWxv7aSxLcO//PJx1m1pYUxlgivPOoTr/vwcje1pPj77Ldzy4ItAELoTqlMs+NxsUvEoAMd98342NLbzpfccym3LXuHVbW3kco4D5x61H9e/fxbP/r2RSbUp9hsd/Mpv7sjw84fX8b1Fazhgn1H89KI6zr354fAHZiUbGtsZU5lga0snHQXdtKMr49x0/lFc9t/LiUaM+//tFPYfvXstB4XDMJHJ5nhjexsTqlNsbu6gsT3N4fvV7tZzNbanqUnFu243tadZ8dp2TnjLODoyOSoSwYd2fUMbZ920lKljK3mxPvi1NaE6yWNzT8fdWfJCPcceOJZRyRgbGtoZX50c8B7gJ17/APGIsfgLc4rmtXZmeM/3HyIaMRb92yldtTS2ZTgkPBR6T5lsjlg0wnv/8yG2tnRyxydPYEJNaqd1PPHKNi6Y/1cOnlDF3Z+Z3W3ea1tbOWXeEibWpnjk6lP7fI7G9jQnfOsBmjoyjK9KUpmMkss59352NqOSQe/rKfMW48DPLqqjPZ3j879dQWcmx5Jelr+lI0NnJseYUYmiee/7z4epb+5g8ZWndH255G1p7uDuZ9Zz/jFTScQGNiTYmckx+ztBa/DOK058U2cvPO3GoBtvVCLG7R8/jsMm1QS/bqN919KezhYtx656bWsrk8dUdHUPDcS9q9Zz+H61jKtK8swbQcvjqCmj+3zf3J1DvnQvZvD8188e0Gtlc8E4yujKBId9+V5aO7P8/OI6Tg3HHQufu76po9vntaEtzaLVG3nPrEn9vj87O+pCOpvj0Ze28oXfP01VMsbCz58MBJ/ZVCw64M9JbxQOAgQB8u7vP8T46iR//OQJpS5nj3L3Pv/Yd+XLq7UzQyoWJaJDpOyV3vW9pWBwT48fEQNx2o1L2NDQzlNfPvNNfSEPJwoHEREpop3gRERklygcRESkiMJBRESKKBxERKSIwkFERIooHEREpIjCQUREiigcRESkiMJBRESK7BV7SJtZPfDKm3iKccCun0hh+NlblgO0LMOVlmV42t1lOcDdx/c2Y68IhzfLzJb3tQv5SLK3LAdoWYYrLcvwNBjLom4lEREponAQEZEiCofA/FIXsIfsLcsBWpbhSssyPO3xZdGYg4iIFFHLQUREiigcRESkSFmHg5mdbWYvmNlaM7u61PXsKjN72cyeMbMVZrY8nDbWzBaa2d/C/8eUus7emNnPzWyTma0qmNZr7Rb4j3A9rTSzo0tXebE+luVaM3sjXDcrzOzdBfOuCZflBTM7qzRVFzOzKWa22MxWm9mzZvaZcPqIWy/9LMtIXC8pM3vMzJ4Ol+Wr4fQDzezRsObfmFkinJ4Mb68N50/brRd297K8AFHgReAgIAE8DRxW6rp2cRleBsb1mPYd4Orw+tXA9aWus4/aZwNHA6t2VjvwbuAewIDjgEdLXf8AluVa4Mpe7ntY+FlLAgeGn8FoqZchrG0ScHR4vRpYE9Y74tZLP8syEteLAVXh9TjwaPh+/xa4IJx+C/CJ8PongVvC6xcAv9md1y3nlsOxwFp3f8ndO4HbgXNLXNOecC7wy/D6L4HzSldK39x9KbC1x+S+aj8X+G8PLANGm9mkISl0APpYlr6cC9zu7h3uvg5YS/BZLDl3X+/uT4bXm4DngP0Zgeuln2Xpy3BeL+7uzeHNeHhx4FTg9+H0nuslv75+D5xmZrarr1vO4bA/8FrB7dfp/8MzHDmwwMyeMLPLwmn7uvv68PoGYN/SlLZb+qp9pK6rK8Lulp8XdO+NiGUJuyLeRvArdUSvlx7LAiNwvZhZ1MxWAJuAhQQtm+3ungnvUlhv17KE8xuAfXb1Ncs5HPYGJ7r70cC7gE+Z2ezCmR60K0fktsojufbQj4C3AEcB64EbS1rNLjCzKuAPwGfdvbFw3khbL70sy4hcL+6edfejgMkELZoZg/2a5RwObwBTCm5PDqeNGO7+Rvj/JuAOgg/NxnzTPvx/U+kq3GV91T7i1pW7bwz/oHPAT9jRRTGsl8XM4gRfpre5+x/DySNyvfS2LCN1veS5+3ZgMfBOgm68WDirsN6uZQnn1wJbdvW1yjkcHgcODkf8EwQDN3eVuKYBM7NRZladvw6cCawiWIaLwrtdBNxZmgp3S1+13wV8JNw65jigoaCbY1jq0ff+DwTrBoJluSDcouRA4GDgsaGurzdhv/TPgOfc/bsFs0bceulrWUboehlvZqPD6xXAGQRjKIuBD4R367le8uvrA8ADYYtv15R6JL6UF4KtLdYQ9N/NLXU9u1j7QQRbVzwNPJuvn6Bv8X7gb8AiYGypa+2j/l8TNOvTBP2ll/ZVO8HWGj8I19MzQF2p6x/Astwa1roy/GOdVHD/ueGyvAC8q9T1F9R1IkGX0UpgRXh590hcL/0sy0hcL7OAp8KaVwFfDqcfRBBga4HfAclweiq8vTacf9DuvK4OnyEiIkXKuVtJRET6oHAQEZEiCgcRESmicBARkSIKBxERKaJwEAHMrDn8f5qZXbiHn/uLPW7/ZU8+v8hgUDiIdDcN2KVwKNhLtS/dwsHdj9/FmkSGnMJBpLtvAyeFx/r/XHjAs3lm9nh4sLaPA5jZKWb2kJndBawOp/0pPAjis/kDIZrZt4GK8PluC6flWykWPvcqC87LcX7Bcy8xs9+b2fNmdlv+qJpm9u3wHAUrzeyGIX93pGzs7BePSLm5muB4/+cAhF/yDe5+jJklgUfMbEF436OBmR4c4hngEnffGh7i4HEz+4O7X21mV3hw0LSe/pHgAHBHAuPCxywN570NOBz4O/AIcIKZPUdwyIcZ7u75QyqIDAa1HET6dybB8YNWEBzyeR+C4+4APFYQDACfNrOngWUEBz47mP6dCPzagwPBbQQeBI4peO7XPThA3AqC7q4GoB34mZn9I9D6JpdNpE8KB5H+GfCv7n5UeDnQ3fMth5auO5mdApwOvNPdjyQ4Fk7qTbxuR8H1LBDz4Nj8xxKcwOUc4N438fwi/VI4iHTXRHBaybz7gE+Eh3/GzN4aHgW3p1pgm7u3mtkMgtM45qXzj+/hIeD8cFxjPMHpRvs8Emh4boJad78b+BxBd5TIoNCYg0h3K4Fs2D30C+D7BF06T4aDwvX0furVe4HLw3GBFwi6lvLmAyvN7El3/+eC6XcQHJf/aYIjiF7l7hvCcOlNNXCnmaUIWjSf360lFBkAHZVVRESKqFtJRESKKBxERKSIwkFERIooHEREpIjCQUREiigcRESkiMJBRESK/H/Y/IATd6CwqAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.errorbar(data_jastrow[\"Energy\"].iters, data_jastrow[\"Energy\"].Mean, yerr=data_jastrow[\"Energy\"].Sigma)\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d216a4d6",
   "metadata": {},
   "source": [
    "## 6. Neural-Network Quantum State \n",
    "\n",
    "We now want to use a more sofisticated ansatz, based on a neural network representation of the wave function. At this point, this is quite straightforward, since we can again take advantage of automatic differentiation. \n",
    "\n",
    "Let us define a simple fully-connected feed-forward network with a ReLu activation function and a sum layer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "2b5c30e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "class FFN(nn.Module):\n",
    "    \n",
    "    # You can define attributes at the module-level\n",
    "    # with a default. This allows you to easily change\n",
    "    # some hyper-parameter without redefining the whole \n",
    "    # flax module.\n",
    "    alpha : int = 1\n",
    "            \n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "\n",
    "        # here we construct the first dense layer using a\n",
    "        # pre-built implementation in flax.\n",
    "        # features is the number of output nodes\n",
    "        # WARNING: Won't work with complex hamiltonians because\n",
    "        # of a bug in flax. Use nk.nn.Dense otherwise. \n",
    "        dense = nn.Dense(features=self.alpha * x.shape[-1])\n",
    "        \n",
    "        # we apply the dense layer to the input\n",
    "        y = dense(x)\n",
    "\n",
    "        # the non-linearity is a simple ReLu\n",
    "        y = nn.relu(y)\n",
    "                \n",
    "        # sum the output\n",
    "        return jnp.sum(y, axis=-1)\n",
    "    \n",
    "model = FFN(alpha=1)\n",
    "\n",
    "vstate = nk.vqs.MCState(sampler, model, n_samples=1008)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a72882a6-57e8-4462-a333-24a500e71020",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "76a17a81-31cd-44c7-8e5a-d7d8743df8c7",
   "metadata": {},
   "source": [
    "*Warning*: Flax has a bug with its layers, where they drop the imaginary part\n",
    "of complex numbers if the layer has real weights.\n",
    "This is not a problem in the above example, but if you try to work on more complex\n",
    "problems where you work with complex numbers you should rather use the equivalent \n",
    "`nk.nn.Dense` which contains a fix for this bug."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f29016f4",
   "metadata": {},
   "source": [
    "We then proceed to the optimization as before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "279c052e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:05<00:00, 55.24it/s, Energy=-25.434 ± 0.012 [σ²=0.140, R̂=0.9989]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized energy and relative error:  -25.420 ± 0.012 [σ²=0.145, R̂=0.9973] 0.0027729382276211805\n"
     ]
    }
   ],
   "source": [
    "optimizer = nk.optimizer.Sgd(learning_rate=0.1)\n",
    "\n",
    "# Notice the use, again of Stochastic Reconfiguration, which considerably improves the optimisation\n",
    "gs = nk.driver.VMC(H, optimizer, variational_state=vstate,preconditioner=nk.optimizer.SR(diag_shift=0.1))\n",
    "\n",
    "log=nk.logging.RuntimeLog()\n",
    "gs.run(n_iter=300,out=log)\n",
    "\n",
    "ffn_energy=vstate.expect(H)\n",
    "error=abs((ffn_energy.mean-eig_vals[0])/eig_vals[0])\n",
    "print(\"Optimized energy and relative error: \",ffn_energy,error)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26014bdc-f72b-466f-b7d7-761a7499ce25",
   "metadata": {},
   "source": [
    "And we can compare the results between the two different ansatze:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "10580934-c128-47ef-beb4-315420d63b9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy/UlEQVR4nO3deXxddZn48c9z99zsW9c0bVpoC5RupEDBQtmUcUNABsQZRX7awVFxREGwDgOOzKBIcRgURGd0ZkQtLnVjLxapYKEL3RfoQku6pmmaPXd9fn+cc9M0S5OG3Nyked6v133l3u8595zn3Jv2yXc536+oKsYYY0x7nkwHYIwxZvCx5GCMMaYTSw7GGGM6seRgjDGmE0sOxhhjOvFlOoD+UFJSohMmTMh0GMYYM6SsXr36sKqWdrXtlEgOEyZMYNWqVZkOwxhjhhQR2d3dNmtWMsYY00lGkoOIXCcim0QkKSKVHbbdJSLbRWSbiLwvE/EZY8xwl6lmpY3ANcAP2heKyJnADcBZwBhgqYhMVtXEwIdojDHDV0aSg6puARCRjpuuAn6hqhFgl4hsB84F/jqwERpjhppYLEZVVRWtra2ZDmXQCYVClJWV4ff7e/2ewdYhPRZY0e51lVvWiYgsABYAlJeXpz8yY8ygVlVVRW5uLhMmTOjqD89hS1WpqamhqqqKioqKXr8vbX0OIrJURDZ28biqP46vqo+raqWqVpaWdjkSyxgzjLS2tlJcXGyJoQMRobi4+KRrVGmrOajq5X14215gXLvXZW6ZMcb0yBJD1/ryuQy2oay/B24QkaCIVACnA69nOCZjjBl2MjWU9WoRqQLmAk+JyHMAqroJeBLYDDwLfC7tI5V+/AHnYYwx75LX62XmzJltj/vvv7/fjr127VqefvrpfjteTzI1WmkJsKSbbfcB9w1IIK310HoU/NkDcjpjzKktKyuLtWvXpuXYa9euZdWqVbz//e9Py/E7GmzNSgPr8FtwcCNEGjIdiTHmFFVXV8eUKVPYtm0bAB/72Mf44Q9/CMBnP/tZKisrOeuss/iXf/mXtvesXLmSCy64gBkzZnDuuedSV1fH3XffzeLFi5k5cyaLFy9Oe9yDbSjrwPJ43Se2VKoxp5r58+f36/FeeumlHvdpaWlh5syZba/vuusurr/+eh555BFuuukmvvjFL1JbW8tnPvMZAO677z6KiopIJBJcdtllrF+/nqlTp3L99dezePFi5syZQ319PeFwmG984xusWrWKRx55pF+vqzvDPDm4l2/raBtj+kF3zUpXXHEFv/zlL/nc5z7HunXr2sqffPJJHn/8ceLxOPv372fz5s2ICKNHj2bOnDkA5OXlDVT4x7HkADD/jszGYYzpd735S3+gJJNJtmzZQjgcpra2lrKyMnbt2sV3vvMdVq5cSWFhITfddNOgurt7ePc5eN1byZM2dZMxJn0eeughzjjjDH72s5/xqU99ilgsRn19PdnZ2eTn53Pw4EGeeeYZAKZMmcL+/ftZuXIlAA0NDcTjcXJzc2loGLj+0WFec3D7HJLxzMZhjDkldOxzuPLKK/nUpz7Fj370I15//XVyc3O56KKL+OY3v8m9997LrFmzmDp1KuPGjePCCy8EIBAIsHjxYr7whS/Q0tJCVlYWS5cu5ZJLLuH+++9n5syZbX0Z6TTMk4N7+ZYcjDH9IJHouhViy5Ytbc8XLVrU9vwnP/lJl/vPmTOHFStWdCpP1SYGwvBuVrLkYIwxXbLkANbnYIwxHQzz5GB9DsYY05VhnhysWckYY7piyQEsORhjTAeWHMCSgzHD1PU/+CvX/8BWIe7K8E4OkupzsA5pY8y7l5OTc9Lv+e53v0tzc3Maonl3hndy8HhAPFZzMMZkzImSQ3f3TQyE4Z0cwGlasuRgjOknjY2NXHbZZcyePZuzzz6b3/3udwA0NTXxgQ98gBkzZjBt2jQWL17Mww8/zL59+7jkkku45JJLAKf28eUvf5kZM2bw17/+lUWLFjFt2jSmTZvGd7/7XQAeeOABHn74YQC+9KUvcemllwLwpz/9iY9//OP9ch3D+w5psORgzCno3j9sYvO++h7327zf2ac3/Q5njsnjXz50Vo/7hUIhlixZQl5eHocPH+b888/nwx/+MM8++yxjxozhqaeeApx1HvLz81m0aBHLli2jpKQEcJLIeeedx4MPPsjq1av58Y9/zGuvvYaqct5553HxxRczb948HnzwQW699VZWrVpFJBIhFouxfPlyLrrooh5j7I1MLRN6nYhsEpGkiFS2K79CRFaLyAb356VpD8bjsz4HY0y/UVW+9rWvMX36dC6//HL27t3LwYMHOfvss3nhhRf46le/yvLly8nPz+/y/V6vl2uvvRaAv/zlL1x99dVkZ2eTk5PDNddcw/LlyznnnHNYvXo19fX1BINB5s6dy6pVq1i+fDnz5s3rl+vIVM1hI3AN8IMO5YeBD6nqPhGZBjwHjE1rJB6v1RyMOcX05i98OFZjWPwPc/vt3E888QTV1dWsXr0av9/PhAkTaG1tZfLkyaxZs4ann36ar3/961x22WXcfffdnd4fCoXwer1dHPkYv99PRUUFP/nJT7jggguYPn06y5YtY/v27Zxxxhn9ch0ZqTmo6hZV3dZF+Ruqus99uQnIEpFgWoOxZiVjTD+qq6tjxIgR+P1+li1bxu7duwHYt28f4XCYv/u7v+P2229nzZo1ACecinvevHn89re/pbm5maamJpYsWdJWM5g3bx7f+c53uOiii5g3bx6PPfYYs2bNQkT65ToGc5/DtcAaVY10tVFEFgALAMrLy/t+FksOxph+EI/HCQaDfPzjH+dDH/oQZ599NpWVlUydOhWADRs2cPvtt+PxePD7/Tz66KMALFiwgCuvvJIxY8awbNmy4445e/ZsbrrpJs4991wAPv3pTzNr1izASQ733Xcfc+fOJTs7m1Ao1G9NSgCiaVoiU0SWAqO62LRQVX/n7vMS8BVVXdXhvWcBvwfeq6o7ejpXZWWlrlq1qqfduvbQNKi4CD7y/b693xgzKGzZsuWkm1T6s1lp3bp1fOYzn+H1119/18dKh64+HxFZraqVXe2ftpqDql7el/eJSBmwBPhEbxLDu2Z9DsYMW/3V1/DYY4/x8MMPtw01PRUMqmYlESkAngLuVNVXBuSk1qxkjHmXbrnlFm655ZZMh9GvMjWU9WoRqQLmAk+JyHPups8DpwF3i8ha9zEircFYcjDGmE4yUnNQ1SU4TUcdy78JfHNAg7H7HIwxphObPsP6HIwxphNLDtasZMzw9eMPOA/TiSUHSw7GmH7i9XqZOXNm2+Ptt9/mpZdeIj8/v63s8sudgZz33HMP4XCYQ4cOtb2/L1N+p8ugGq2UEdbnYIzpJ1lZWaxdu/a4srfffpt58+bxxz/+sdP+JSUlPPjgg3zrW98aoAh7z2oO1udgjMmQm2++mcWLF3PkyJFMh9KJ1Rw8Poi1ZDoKY0x/euZOOLCh5/0OrHd+9qbfYdTZ8Df3n3CXlpYWZs6cCUBFRQVLljiDMpcvX95Wft1117Fw4ULAaUa6+eab+Y//+A/uvffenmMYQJYcrM/BGNNPumpWArptVgK49dZbmTlzJl/5ylfSHN3JseRgycGYU08Pf+G3SdUYPvVU+mLpQUFBATfeeCPf+973MhZDVyw5eLzWIW2MyajbbruNOXPmEI8Pnj9UrUPaag7GmAwrKSnh6quvJhLpcoWCjLCagyUHY0w/aWxs7FQ2f/585s+f36n8nnvuOe71okWLWLRoUZoiO3mWHCw5GDN8ZbCvYbCzZiXrczDGmE4sOVjNwZhTRrpWthzq+vK5WHKw5GDMKSEUClFTU2MJogNVpaamhlAodFLvsz4HSw7GnBLKysqoqqqiuro606EMOqFQiLKyspN6jyUHm3jPmFOC3++noqIi02GcMjK1TOh1IrJJRJIiUtnF9nIRaRSR9N9PbhPvGWNMJ5nqc9gIXAO83M32RcAzAxKJNSsZY0wnmVpDeguAiHTaJiIfAXYBTQMSTCo5qEIX8RhjzHA0qEYriUgO8FWgx7lrRWSBiKwSkVXvqgPK4+ZHTfb9GMYYc4pJW3IQkaUisrGLx1UneNs9wEOq2vke9A5U9XFVrVTVytLS0r4H6vE6P61pyRhj2qStWUlVL+/D284DPioi3wYKgKSItKrqI/0aXHupmkMyDgTTdhpjjBlKBtVQVlWdl3ouIvcAjWlNDNAhORhjjIHMDWW9WkSqgLnAUyLyXCbiANolB7vXwRhjUjI1WmkJsKSHfe4ZkGCsz8EYYzoZVKOVMsKalYwxphNLDpYcjDGmE0sOlhyMMaYTSw7WIW2MMZ1YcrAOaWOM6cSSgzUrGWNMJ5YcLDkYY0wnlhysz8EYYzqx5GB9DsYY04klB2tWMsaYTiw5WHIwxphOLDmkksOzX8tsHMYYM4hYckj1OahmNg5jjBlELDl4/O4TWybUGGNSLDl4A85PqzkYY0wbSw5et+Yw77bMxmGMMYOIJQefu250IprZOIwxZhDJ1DKh14nIJhFJikhlh23TReSv7vYNIhJKazCpZqV4JK2nMcaYoSQjy4QCG4FrgB+0LxQRH/BT4O9VdZ2IFAOxtEaSSg6J9J7GGGOGkkytIb0FQEQ6bnovsF5V17n71aQ9mLbkYDUHY4xJGWx9DpMBFZHnRGSNiNzR3Y4iskBEVonIqurq6r6fsa3PwWoOxhiTkraag4gsBUZ1sWmhqv7uBPG8B5gDNAMvishqVX2x446q+jjwOEBlZWXfx6F6vCAe63Mwxph20pYcVPXyPrytCnhZVQ8DiMjTwGygU3LoV96gjVYyxph2Bluz0nPA2SISdjunLwY2p/2s3oA1KxljTDuZGsp6tYhUAXOBp0TkOQBVrQUWASuBtcAaVX0q7QF5/dYhbYwx7WRqtNISYEk3236KM5x14PisWckYY9obbM1KmeH1Q9ySgzHGpFhyAOuQNsaYDiw5gNshbcnBGGNSLDkA+Cw5GGNMe71KDiLyoIicle5gMsZqDsYYc5ze1hy2AI+LyGsicouI5KczqAHnDViHtDHGtNOr5KCqP1LVC4FPABOA9SLyMxG5JJ3BDRirORhjzHF63ecgIl5gqvs4DKwDbhORX6QptoFjycEYY47Tq5vgROQh4EM4cxz9m6q+7m76lohsS1dwA8Y6pI0x5ji9vUN6PfB1VW3qYtu5/RhPZngDNiurMca009vksA6Y0mFxnjpgt6rW9XtUA80m3jPGmOP0Njl8H2fq7PWAANOATUC+iHxWVZ9PU3wDw/ocjDHmOL3tkN4HzFLVSlU9B5gF7ASuAL6druAGjE28Z4wxx+ltcpisqptSL1R1MzBVVXemJ6wB5vVbcjDGmHZ626y0WUQeBVLDVq93y4LA0G+s9wadDmlVOL5fxRhjhqXe1hw+CWwH/sl97ARuwkkMQ/9GOG8AUEgmMh2JMcYMCj3WHNyb355W1UuAB7vYpbHfoxpoXr/zMxEFb0bWPzLGmEGlx5qDqiaAZH/OpyQi14nIJhFJikhlu3K/iPyPiGwQkS0icld/nfOEfEHnpy0VaowxQO/7HBqBDSLyAtB2I5yq3trH824ErgF+0KH8OiCoqmeLSBinX+Pnqvp2H8/TO201h6HffWKMMf2ht8nhN+6jX6jqFgDp3PmrQLaI+IAsIArU99d5u+V1aw52l7QxxgC9TA6q+j8ikgWUq2o651L6FXAVsB8IA19S1SNd7SgiC4AFAOXl5e/urN6A89OGsxpjDND7xX4+BKwFnnVfzxSR3/fwnqUisrGLx1UneNu5QAIYA1QAXxaRiV3tqKqPuzflVZaWlvbmMrrns+RgjDHt9bZZ6R6c/7hfAlDVtd39p52iqpf3IZ4bgWdVNQYcEpFXgEqcobPpYzUHY4w5Tm/vc4h1McFesr+DAfYAlwKISDZwPrA1Dec5XqrPwTqkjTEG6H1y2CQiNwJeETldRP4TeLWvJxWRq0WkCpgLPCUiz7mbvgfkiMgmYCXwY1Vd39fz9FpqtJJ1SBtjDND7ZqUvAAuBCPBz4DngX/t6UlVdAizporwRZzjrwLJmJWOMOU5v15BuVtWFqjrH7QReqKqt6Q5uwKQ6pJ/9WmbjMMaYQaK3y4ROBr4CTGj/HlW9ND1hDbBUzSEt3SjGGDP09LZZ6ZfAY8CPcIaanlr8Yedn0pKDMcZA75NDXFUfTWskmZRKDnP/MbNxGGPMINHb0Up/EJF/FJHRIlKUeqQ1soHkz3J+xpozG4cxxgwSva05fNL9eXu7MgVOeCPckBHIdn5acjDGGKD3cytVpDuQjPL6weOHqCUHY4yBHpqVROSOds+v67Dt39IVVEb4w1ZzMMYYV099Dje0e95x4Z0r+zmWzApYcjDGmJSekoN087yr10ObP2zNSsYY4+opOWg3z7t6PbRZs5IxxrTpqUN6hojU49QSstznuK9DaY1soFmzkjHGtDlhclBV70AFknH+LGtWMsYYV29vgjv1+bOt5mCMMS5LDinWrGSMMW0sOaRYs5IxxrSx5JBizUrGGNMmI8lBRB4Qka0isl5ElohIQbttd4nIdhHZJiLvG7Cg/FmWHIwxxpWpmsMLwDRVnQ68iXv3tYiciXNX9lk4d2B/X0QGZsRUIAzJOMRtqVBjjMlIclDV51U17r5cAZS5z68CfqGqEVXdBWwHzh2QoPypmVmbBuR0xhgzmPV2yu50uhlY7D4fi5MsUqrcsk5EZAGwAKC8vLxPJ25ojfG3j/2VrICX35ybWtOhBbIK+3Q8Y4w5VaSt5iAiS0VkYxePq9rtsxCIA0+c7PFV9XFVrVTVytLS0j7FuKO6iS0HGmiMxI+t6WAjlowxJn01B1W9/ETbReQm4IPAZaqamqdpLzCu3W5lbllaZAec7oxkUtutBmfNSsYYk6nRSlcCdwAfVtX2f6r/HrhBRIIiUgGcDryerjiy3OSw4OJJx9aRjrWk63TGGDNkZGq00iNALvCCiKwVkccAVHUT8CSwGXgW+JyqJtIVRDjgVJya2zcr/fG2dJ3OGGOGjIx0SKvqaSfYdh9w30DEEXZrDs2xxLFmpfTlImOMGTKG9R3SQZ8Hj0BLNAGBHKcwmcxsUMYYMwgM6+QgIoQDPpoiCQjmOYXJ+InfZIwxw8CwTg7gNC21xOIQcpPDrBszG5AxxgwClhwCXpqjCfAFwRuE1vqe32SMMae4YZ8cslLNSuDUHiKWHIwxZtgnh+xUsxI4/Q6RhswGZIwxg8CwTw5ZqWYlcGoO1qxkjDGWHMIBL82pZqWgNSsZYwxYciA74KO5rVkp12oOxhiDJQeyAl7nJjiAUL7VHIwxBksOhAPeY6OVrEPaGGMASw6EAz5aYgln2u6QmxxsCg1jzDBnycGdfK81nppCQyFqtQdjzPBmycFNDk2RxLEpNKxT2hgzzFlycNd0aIm2m3zPOqWNMcOcJYdUzSEad4aygnVKG2OGvWGfHFJLhTZHE85QVoA/fimDERljTOZlag3pB0Rkq4isF5ElIlLgll8hIqtFZIP789J0x5IddJqVFi7ZAKECpzBhazoYY4a3TNUcXgCmqep04E3gLrf8MPAhVT0b+CTwf+kOJNWslEgq5I1xChORdJ/WGGMGtYwkB1V9XlVTf56vAMrc8jdUdZ9bvgnIEpFgOmMpCAcAiCcVAmHIGQkVF6XzlMYYM+gNhj6Hm4Fnuii/Flijql3+GS8iC0RklYisqq6u7vPJC7L8AFw/Z5xTUDgBat/u8/GMMeZUkLbkICJLRWRjF4+r2u2zEIgDT3R471nAt4B/6O74qvq4qlaqamVpaWmf4wwHvAS8Ho42x5yCgvFwdHefj2eMMacCX7oOrKqXn2i7iNwEfBC4TFW1XXkZsAT4hKruSFd87c5HftjP0eaoU1A4ATb+ChIx8PrTfXpjjBmUMjVa6UrgDuDDqtrcrrwAeAq4U1VfGah4CsP+YzWHwvGgSah7Z6BOb4wxg06m+hweAXKBF0RkrYg85pZ/HjgNuNstXysiI9IdTEFWgNr2NQeAX/x9uk9rjDGDVtqalU5EVU/rpvybwDcHOBwKwn72HHErMAXjnZ/x1oEOwxhjBo3BMFop4wrC/mM1h7wxgFhyMMYMaxmpOQw2heHAsT4HjxeKKmD0jMwGZYwxGWQ1ByA/7CcSTx5bLrRwAtTacFZjzPBlyQGnQxrgaIvbtFQw3m6EM8YMa5YccIayAtQ2pYazToCWI7bojzFm2LLkgNOsBHDbk2udgtRwVrtT2hgzTFlyAIqy202+B86NcAC/+nSGIjLGmMyy5ACU5jgTv8biSaegZAp4fNCwH47N7GGMMcOGJQecoaw+j/D+6aOdgkAYrrwfInXw6AWZDc4YYzLAkgPg8QiluUEO1bebHfycm8AbhLq9GYvLGGMyxZKDa0ReiEMN7e6K9vrhsn92ag+PnGfNS8aYYcWSg2tEx5oDwOxPgi8LDm+FR+ZkJjBjjMkASw6uEbnB42sOAKE8uHMP+MNO57QxxgwTlhxcI/NC1DbHiMQTx2/wBeDSf4ZoI/xgfkZiM8aYgWbJwTUi1xnOWt3QxZLVZ18H4oGatyDaNMCRGWPMwLPk4BqR5ySHQ10lh5xSuP6nTu3h4dkDHJkxxgw8Sw6uEbkhAL76q/Vd7zD1A3DeZ6HxADw6bwAjM8aYgZepNaQfEJGtIrJeRJa4a0e3314uIo0i8pWBimlCSTYALbFE9ztdcpdz78PB9fCNEnj80gGKzhhjBlamag4vANNUdTrwJnBXh+2LgGcGMqCcoI+ywixmjivofqdQPnxxLYQKIBmDmu2QdJNJYzUk4gMQqTHGpF+m1pB+vt3LFcBHUy9E5CPALmDAe36njsrlzYMNJ94pbwzcuRtefQSeXwjfKIJwiTPFd6gA7tgJIgMSrzHGpMtg6HO4GbeWICI5wFeBezMRyJRRueysbiKamoDvRCo/Bb6QM4qp+bBzB3XLEfjeeekP1Bhj0ixtNQcRWQqM6mLTQlX9nbvPQiAOPOFuuwd4SFUbpYe/vkVkAbAAoLy8vF9injwyl3hS2VHdyBmj8068cyAbvn4Qkkl4eCYE8+HIDji6xynzDIa8a4wxfSOaoTmDROQm4B+Ay1S12S1bDoxzdykAksDdqvrIiY5VWVmpq1atetcxvXWwgSseepmK4jDLbr/k5A+w4Vfw6/8H2SMgXASxFmf6jc/+BcRrCcMYM6iIyGpVrexqW0b6HETkSuAO4OJUYgBQ1Xnt9rkHaOwpMfSn00bkUFGSzUj3noeTduZH4PdfgKZDziPlX0sAgaJJ8I+vgq+PxzfGmAGSkeQAPAIEgRfc5qMVqnpLhmJpIyJcM2ssD77wJlW1zZQVhk/uAF4ffG0fJOPw6IXgDUCkwRnZ1FoHR7bDN0dCQTnc8AQkok7tYvQMCOam56KMMaYPMtas1J/6q1kJ4J0jzcz79jLKi7J4+Y5+vI9BFXYugydvcqYBb88XgsvudhLF3M+BP6v/zmuMMd0YdM1Kg9m4ojCTR+a03THdb0Rg0qVw1x6ofRv2rnHWjIi1wm8+A899zdlv+SK49Q1nwr9QgQ2LNcZkhCWHLsyfMoKfvPI2TZE42cE0fESFE5xHSiAbWo9CViH8/GPw4GSnPJgHt22GTb+FIzth5o1Qcnr/x2OMMR1Ys1IXXt1+mBt/9Bqnj8jhhdsu7rfj9squl+HtV5waw0v/3nl7IBfOvwUu/qpT8wCn9lG/10k4Hu+AhmuMGbqsWekkzakooqwwi3DQh6rS0z0X/ariIucBsP6XTkf2hx6CsnPhR5c7SeDlB+Dl7wAdErvHD+FiiDXDFd9w1sGONsHeVTBx/sBdgzFmyLOaQzd+8foe7vzNBqaMzOG5Lw1w7aEn/znHuSt75o1OR3ciAiOnwXMLIdbkdHDHW+GaH8Kzd0JzDdz0NEy4MNORG2MGEas59MG155TxyLLtBP3ega899OQLK7sun/lxZznT/DK4f7zT0Z3y6n8603uMv9C5Qe+xi5yO8Uu/DnM+bTfoGWOOY8mhG36vhy9cehpf/fUGHv3zDm65aBIezyBKEF3xBaBwvPP886/Daz9wEsGKR+HNZ5wHOM1PgWyINsAztzt9F5WfylzcxphBx5qVTiCWSPLZn65h6ZaD5IZ8PPPFeSd/Y9xg0HQY1vwvjJkF+96AFd937qn49Ivwx3+C6q1QcTFMvx4a9sH0GyCYk+mojTFpdqJmJUsOPVBVfrm6ijt/vR6PCOdPLOZ9Z43kxvPG4x3sNYne+MHFTnLwh51mJ3CG1N66Ftb+zLmLe8YNkNvVHIrGmKHMkkM/2HW4iXv/sIm9tS28daiRkpwAt79vCjuqm7j8jJHMmVB4wn6JZFIRYXD1XbTXWg9vPU+ybi+epXcDQttoKI8XPvYkjJgKG38D066F/LHONtVBc6NeNJ6koTVGcY7NXdVbsUSS6oYIYwrsrvzhyJJDP1JVHlr6Fg+/+FZbmUfA5/VwTnkhs8cXsHJXLQBzJxWTE/QR9Hv4t6e2EPB5+PPtl/DGO7Vk+X1E4gneOdKMz+uhviXGpn31ZAe9zCgrYF9dK36PkJflxyNQ3RDhrLH5VDdEyPJ7GVuYxZb99UTjSY40R7lhTjk+j/DWoQZqm2K8uPUg08sKyAv5yQ35+PEru/j2R2cwoTjMUxv2s2JnDReeVsL4omx+ufodKkqyCfg8/PNvN3Jtzkb+uehFVpd9An/pJM5Z+WVCNZtIpYBGTy7Lpj/AOeUFFP/hJvYGTyP2yWfYV9cCBzawqibImaefxuVnjuB/X91NLJlkZlkBE0tz+NHynXi9Qm7QRzSh1DVHAThvYjGHGyMkksrVs8aSHfRR2xTl9+v2UVXbQsjvJcvvdT/fIxyob6WiJIeckI8LJhVT3RDhk//9OtFEkq+8dwpjCkLsO9rK5JG5XDy5lD9tPciW/Q2oKlsPNHDBpGLOGpvP1gMNJJPKpVNHMK4ozNHmKM9tOsD7zx7Nunfq+P5L25lUmsPoghAfmTmWrQfqmViSw6rdtVwwqZg9R5ppaI0TSySpbY4ydVQe5UVhIvEE//WXXdQ2RTl9ZC41jVEONrQytiCLa2eX0RyN88vVVZQXhUkklXDAy4SSbCqKs3n0pR2sqzrK7PGFXHHmSOqaYwR9HsYVhfnzm9UcqGtlfHGYgM9DwOvh9JE5nDO+iNW7j7D8rcOE/F7qWmK8uqOGrfvriSWSzC4v5L1njeRwY5TK8YXkhJzuxs/87yqaIgk+MnMMb7xzlB9+opLJI3PZuLeOZVsPkR/28+Sqd/jbynHkZ/mpb42z/2gLq3bX8t4zR/K+s0ZR0xRl07463thzlLLCLII+L8XZAV7ZcZjZ5YUUZQf4zvPbuGZWGbPKC5xbeLZVU1XbTGMkzgenjyE76KMoHCDg8/D0hv2U5ga5eHIp2w810tAao7w4m7V7jvLIsrfwiPDv15zN9LJ8WqJJJo/Kwe/xsHz7YRpb4xRlB/jFyj3UNEa58LQS/rqzhh2HGvngjNFUji/iJ6/uIsvv47yKIvYcaWZ0QYj8LD+qMKHYWS64MRKjJZZgRlkBbx1qJOT3MrOsgOc2HyAc8LLjkLMWmc8r5AR9nD+xmC376/F5hXGFYUJ+L2MKQrz85mHermmiODtAYyTO3qMtbNlfz7Qx+Zw3sZjR+SF+tbqK7KCXi04vJTvoY13VUcqLwrx9uIny4mwumFSM3+uhNZZg07461r5Tx8TSbIJeD0G/h3PGF/X5/zNLDv1MVfm/FbuZVJrDzHEFfP+l7dS3xPn563uIJ5VpY/Pwez28sedo23t8HiGePPFn7fc6+/THV+L1CIkO5yvJCdASTdAUTSBC23na1RHwe4VY4vj35dLMAt8fadQsdupovuV/nCJpBKBRQ+RIK9WaTwQ/ZXKYuHrYoWNYmjyHczxv0qxB/pycQZAoMfEzjkNs03E8mbgYEWeUlE9jxPGSdNefOlN2UyaHWJOcTINkgyoR/N1en9cjlBeFqWmMUN96/HKtHfcNeD1EE50XdDp7bD5bD9QTSygegaTCyLwgDa1xmqMnWFu8GwKMKchi79EWPALlRWF21zS3fdbtv4Pj3ifwntNK+MtbhzveyQLQFlt7I3KDHGqIHFd2XkURJblBRuQG+emK3Z2+19Sxgj4vLbEEAZ+HcMDLmPwsth6obztHV79LE0uz2Vl9/GKNHX/H28fZ8RgiMHlELrFkstNxuvtcAC6ZUkpNU5T1VXVd79Du3AGvh9Z4Ep9HeM/pJby0rRqA8cVh9ta2EE8quUEfDZHeLe/bm3/DPRGBM0fnsXlf/bHfAzrdsdRJSU6QmsZIp/0Kwn7W3v3edxGPJYcBEXf/w/F5nf/gDjW04hHhaHOUEXkhtu5vYNXuI5xWmsPRlhg+j3DBpBIUJS/kJxzwcqghQm1zlNNH5JJUpaq2hdZYghG5QapqWxiRF+SdIy3sO9rC+ROLEYGkKi9sPkhzNMEZo/Oob4nxN9NG0dAaZ82eWl7YfJBrzynju0vfZFd1E19+7xSumjmGFTuPsH7vUa47Zxx1LTGWv1XNdZXj2F3TRHVDhIklOdQ2R9l7tIWReUFEhGg8yYhQkpx3/swr67eQM+MjTKl/Bd2zghyJEh83lyKtpXbbXxhx+DUiEsLnAW+ite1zUgRB0dzRkF0CjdVI4wHiEiCZM4qWaIy8yP5j/2jEAwixkTNokSwCOUWEat8k6c+ieWQlGw/F8DXtZ3ZoP55YE5Hc8WhzDf6SSTTtXs3O5Giyp15CRYGfRMtRQl6ooYCjNQcZ5W8mUPUqu/0VrG8ZSXl8J5N9h9gnIymkgeLyM/CJ0tDYwBsN+RSNmUTiyC6Ks/0cPlJLXn4Ro+o3EM8di2SX0HpoB/uzz6DkyBryC4vJLp1AqzcHf04h3mA2kd0r2XqoFU+4kCljCqlrTRBu3ofvyDYaw+PYIFOY4D/KhKIQdZIHu/9CMFZPc+FU9ocmUapHKA0liAYKSUZbiAbyWV/jZXuDn0nZrZw3sYjmYCn+N58mp2kPhAshdwyJ2j00B4rwT76cAzW1eGrfptWXy/isFqTxIHW1hwlk5bKiOkBu5CBFYS8Vid3UBMsYceY8GvdtRXNGk1W/A3/9brylU6hvamZLcDq5WQHKYrvIzc4hHiok6c+m4cghCqa8h30tPmpra5g8cQJvbt/JyH1L8UbryBl1GlmFo0lGGtl8OEFxdC8abSKaVcrIkmI8W/9AU8JHS8UV5IycRF3VVnKy/BSVjiHZWM2et9aj8RiaM4KD8Ryi+BmdlSArnE2jN5/xJTmE80poibTiO7oL34G17Np3iNcSU7hmahYekkSb6sguGUfEGyZRs5N4qJijMT8+jeAXJXh4A7XxIKOattDU2Mhv4hcyb9pERjVvIycUwFM8kVjFfKqbEtSue5pR2ULrmLm07NuEr+kA+2PZTPbuJ7+wmPoRc8hu3kcwegTJG0OspYF9WkRT9W7Gjh5DKBjkyLZXSOIh3xejllwCIyZTV1dLTVOc5tZWPOECCsZOYUJhgL17dhCOVDN28iyClX/f5/+zLDn0YP78+f0XjGlzWk4zB1sDNMS9FAXiRJNCwKM0xr3MKznK+cX1hL0J6mI+DrQGGJMVJeBxEuz+liCv1uQxp6iBbG+CsC/JyGCUkmCMbF+CzfVhigJxTstpIeBJUh3xs7clSFKFUaEorZEIY3JgzdF8JmXVMcYdfBVNOn+WBryQUEiosKY2lwnZrYwKRdnTHKSqOcioUJS6mI+ycIRIwkPMPW6WN0lLDBLiIZLwUBCIs63BiSXHl6A26mNcOMKBJohKkEJ/nGxfgtTYhYYoIF6yvAl87q0lNREfu5tDnJ7TQq4/QTQpJBVCXqW6GfZEcjgtp4V8f4J4EmJJyPJBLCn4PV3/+22IednVFGJEKEahP0ZVS5CRgVZyAsfvl1CojfpoiPvI9cUpDMSpjvjxCmxvzGJiqJ6R2dAY95LjS7C/JcD+1gDjwhFUYUQoBkBjFLxeD1neEy+x2xSDmliQ0aHocbE3xaAx4ac4EMPngR2NIWJJD1PznOVekgrtx38cbPUTTXooDMTI8fW8rO+RqI9YLO5ei/PBN8a9lLjnOxL1keeL4/NAqlK5uzlE2Jdgc302AY9yQXEdHoG3m0IoMDbUSsCdrSaaFJJJJeRzYm2Me8nzJ9jbEiDHlyDfnyCp0JzwHBdvIgnu35JUR/zEk0Jr0kOB3/kuUtfd8fpT53z2QBEffnxnj9ffHbsJzmTE9sZjw36PRP3HbXvxUBEvHuq5rXRTfU9DatWtYZyoU1wp8MdpTniJJgUPkONLUB/34hUnQQAEPEmiye5vBhSUwkCc2qiv3fkUOpx7VChCTcRPTD1t78v2Jcj2Jjkc9bedL7UtdSyfJAl7kzQlvCQUsrxJIgkPSfcKSwIxWpMeWhJevKJEk0KWN0meL0GeP05rwkOuP0FRIMaKmjzi6jnu8wl4kkwItxJXYU9zkLA3SUPce9xn5xXtFF9JMEZ1xE/Aox0+H2VsVpRYUjgU8YN7jrA3QSTpYUZ+I0mgJeGlwB+nMe5lc302kaQHryhZ3gStCSehNMW9JBE8KEFvkpaEU1ssCsTI9cWpjgRIqJDjS9AY9xJpF0fAk8QrSmvCQ8ibJNeXQIDCgJO4djZlEU0KXoGQx/l82783y5ugLuZHUPwedRNzksb48f89+iSJ36O0uO/P8iY4M6+JkCfJmqO5JFU4K6+JtxqzaIj72n6fPCgTc1rY1xKkOeH8YRBNeigNRqmJ+Al5k/hEqY35jvtdyvE536ni/LGQ40swxv28D0f81Me9gPDhbn9j3x1LDsBLL72U6RCMMWZQsTkTjDHGdJKR5CAiD4jIVhFZLyJLRKSg3bbpIvJXEdkkIhtEpJ9X3THGGNOTTNUcXgCmqep04E3gLgAR8QE/BW5R1bOA+UAsQzEaY8ywlZHkoKrPq2pqcPEKoMx9/l5gvaquc/erUdWTH2BujDHmXRkMfQ43A+50oUwGVESeE5E1InJHBuMyxphhK22jlURkKdDVbG0LVfV37j4LgTjwRLt43gPMAZqBF91xuC92cfwFwAKA8vLy/r8AY4wZxtKWHFT18hNtF5GbgA8Cl+mxO/GqgJdV9bC7z9PAbKBTclDVx4HHwbkJrv8iN8YYk6nRSlcCdwAfVtXmdpueA84WkbDbOX0xsDkTMRpjzHCWkekzRGQ7EARq3KIVqnqLu+3vcEYvKfC0qvbY7yAi1cDudxFSCXD4Xbx/sDhVrgPsWgajU+U6wK4lZbyqlna14ZSYW+ndEpFV3c0vMpScKtcBdi2D0alyHWDX0huDYbSSMcaYQcaSgzHGmE4sOTgez3QA/eRUuQ6waxmMTpXrALuWHlmfgzHGmE6s5mCMMaYTSw7GGGM6GdbJQUSuFJFtIrJdRO7MdDwnS0Tedqc1Xysiq9yyIhF5QUTecn8WZjrOrojIf4vIIRHZ2K6sy9jF8bD7Pa0XkdmZi/x43VzHPSKy1/1e1orI+9ttu8u9jm0i8r7MRN01ERknIstEZLM7Zf4X3fIh9b2c4DqG3PciIiEReV1E1rnXcq9bXiEir7kxLxaRgFsedF9vd7dP6PPJVXVYPgAvsAOYCASAdcCZmY7rJK/hbaCkQ9m3gTvd53cC38p0nN3EfhHO1Cgbe4odeD/O5IwCnA+8lun4e7iOe4CvdLHvme7vWRCocH//vJm+hnbxjQZmu89zcabTP3OofS8nuI4h9724n22O+9wPvOZ+1k8CN7jljwGfdZ//I/CY+/wGYHFfzz2caw7nAttVdaeqRoFfAFdlOKb+cBXwP+7z/wE+krlQuqeqLwNHOhR3F/tVwP+qYwVQICKjByTQHnRzHd25CviFqkZUdRewHef3cFBQ1f2qusZ93gBsAcYyxL6XE1xHdwbt9+J+to3uS7/7UOBS4FduecfvJPVd/Qq4TEROtMB6t4ZzchgLvNPudRUn/gUajBR4XkRWu7PUAoxU1f3u8wPAyMyE1ifdxT4Uv6vPu00t/92uaW/IXIfbHDEL5y/VIfu9dLgOGILfi4h4RWQtcAhnobQdwFE9tiZO+3jbrsXdXgcU9+W8wzk5nAreo6qzgb8BPiciF7XfqE7dckiOVR7KsQOPApOAmcB+4MGMRnOSRCQH+DXwT6pa337bUPpeuriOIfm9qGpCVWfiLIp2LjB1IM47nJPDXmBcu9dlbtmQoap73Z+HgCU4vzgHU1V79+ehzEV40rqLfUh9V6p60P0HnQR+yLEmikF/HSLix/kP9QlV/Y1bPOS+l66uYyh/LwCqehRYBszFacJLLbnQPt62a3G353NsgtOTMpyTw0rgdLfXP4DTefP7DMfUayKSLSK5qec4S6xuxLmGT7q7fRL4XWYi7JPuYv898Al3dMz5QF27Zo5Bp0O7+9U43ws413GDO6KkAjgdeH2g4+uO2zb9X8AWVV3UbtOQ+l66u46h+L2ISKmIFLjPs4ArcPpQlgEfdXfr+J2kvquPAn9ya3snL9O98Zl84Iy2eBOnDW9hpuM5ydgn4oywWAdsSsWP0774IvAWsBQoynSs3cT/c5yqfQynzfT/dRc7zoiN77nf0wagMtPx93Ad/+fGud79xzq63f4L3evYBvxNpuPvcC3vwWkyWg+sdR/vH2rfywmuY8h9L8B04A035o3A3W75RJwEth34JRB0y0Pu6+3u9ol9PbdNn2GMMaaT4dysZIwxphuWHIwxxnRiycEYY0wnlhyMMcZ0YsnBGGNMJ5YcjAFEpNH9OUFEbuznY3+tw+tX+/P4xqSDJQdjjjcBOKnk0O5O1e4clxxU9YKTjMmYAWfJwZjj3Q/Mc+f7/5I76dkDIrLSnbDtHwBEZL6ILBeR3wOb3bLfupMgbkpNhCgi9wNZ7vGecMtStRRxj71RnHU5rm937JdE5FcislVEnkjNrCki97vrFKwXke8M+Kdjho2e/uIxZri5E2fO/w8CuP/J16nqHBEJAq+IyPPuvrOBaepM8wxws6oecac5WCkiv1bVO0Xk8+pMnNbRNTiTwM0AStz3vOxumwWcBewDXgEuFJEtONM+TFVVTU2rYEw6WM3BmBN7L878QWtxpn0uxpl7B+D1dokB4FYRWQeswJn87HRO7D3Az9WZDO4g8GdgTrtjV6kzSdxanOauOqAV+C8RuQZofpfXZky3LDkYc2ICfEFVZ7qPClVN1Rya2nYSmQ9cDsxV1Rk48+GE3sV5I+2eJwCfOvPzn4uziMsHgWffxfGNOSFLDsYcrwFnacmU54DPulNAIyKT3VlwO8oHalW1WUSm4izlmBJLvb+D5cD1br9GKc6So93OBuquT5Cvqk8DX8JpjjImLazPwZjjrQcSbvPQT4D/wGnSWeN2ClfT9dKrzwK3uP0C23CallIeB9aLyBpV/Xi78iU4c/Ovw5lF9A5VPeAml67kAr8TkRBOjea2Pl2hMb1gs7IaY4zpxJqVjDHGdGLJwRhjTCeWHIwxxnRiycEYY0wnlhyMMcZ0YsnBGGNMJ5YcjDHGdPL/AY7RbIX/zFc0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_FFN = log.data\n",
    "\n",
    "plt.errorbar(data_jastrow[\"Energy\"].iters, data_jastrow[\"Energy\"].Mean, yerr=data_jastrow[\"Energy\"].Sigma, label=\"Jastrow\")\n",
    "plt.errorbar(data_FFN[\"Energy\"].iters, data_FFN[\"Energy\"].Mean, yerr=data_FFN[\"Energy\"].Sigma, label=\"FFN\")\n",
    "plt.hlines([E_gs], xmin=0, xmax=300, color='black', label=\"Exact\")\n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ff286a6",
   "metadata": {},
   "source": [
    "## 7. Translation Symmetry \n",
    "\n",
    "In order to enforce spatial symmetries we can use some built-in functionalities of NetKet, in conjunction with equivariant layers. \n",
    "\n",
    "The first step is to construct explicitly a graph that contains the edges of our interactions, in this case this is a simple chain with periodic boundaries. NetKet has builtin several symmetry groups that can be used to target specific spatial symmetries. In this case for example after constructing the graph we can also print its translation group. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "72f147d6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PermutationGroup(elems=[Id(), Translation([1]), Translation([2]), Translation([3]), Translation([4]), Translation([5]), Translation([6]), Translation([7]), Translation([8]), Translation([9]), Translation([10]), Translation([11]), Translation([12]), Translation([13]), Translation([14]), Translation([15]), Translation([16]), Translation([17]), Translation([18]), Translation([19])], degree=20)\n"
     ]
    }
   ],
   "source": [
    "graph=nk.graph.Chain(length=N, pbc=True)\n",
    "\n",
    "print(graph.translation_group())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64cabe68",
   "metadata": {},
   "source": [
    "Graphs are in general quite handy when defining hamiltonian terms on their edges. For example we can define our Hamiltonian as "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "ae71872a",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gamma=-1\n",
    "H = sum([Gamma*sigmax(hi,i) for i in range(N)])\n",
    "\n",
    "V=-1\n",
    "H += sum([V*sigmaz(hi,i)*sigmaz(hi,j) for (i,j) in graph.edges()])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0196f9f",
   "metadata": {},
   "source": [
    "We now write a model with an invariant transformation given by the translation group. Notice that we will now use NetKet's own ```nn``` module, instead of Flax, since it contains several additions and also an extended and seamless support for complex layers/parameters.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a0f65b95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "84"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import netket.nn as nknn\n",
    "\n",
    "class SymmModel(nn.Module):\n",
    "    alpha: int\n",
    "\n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "        \n",
    "        x = nknn.DenseSymm(symmetries=graph.translation_group(),\n",
    "                           features=self.alpha,\n",
    "                           kernel_init=nn.initializers.normal(stddev=0.01))(x)\n",
    "        x = nn.relu(x)\n",
    "\n",
    "        # sum the output\n",
    "        return jnp.sum(x,axis=(-1,-2))\n",
    "\n",
    "\n",
    "sampler = nk.sampler.MetropolisLocal(hi)\n",
    "\n",
    "#Let us define a model with 4 features per symmetry\n",
    "model=SymmModel(alpha=4)\n",
    "\n",
    "vstate = nk.vqs.MCState(sampler, model, n_samples=1008)\n",
    "\n",
    "vstate.n_parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e0ca35d",
   "metadata": {},
   "source": [
    "As it can be seen, the number of parameters of this model is greatly reduced, because of the symmetries that impose constraints on the weights of the dense layers. We can now optimize the model, using a few more optimization steps than before.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "16a45cf8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 600/600 [00:31<00:00, 19.19it/s, Energy=-25.4934 ± 0.0014 [σ²=0.0020, R̂=1.0011]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized energy and relative error:  -25.4906 ± 0.0016 [σ²=0.0027, R̂=1.0003] 1.5732015642812893e-05\n"
     ]
    }
   ],
   "source": [
    "optimizer = nk.optimizer.Sgd(learning_rate=0.1)\n",
    "\n",
    "gs = nk.driver.VMC(H, optimizer, variational_state=vstate,preconditioner=nk.optimizer.SR(diag_shift=0.1))\n",
    "\n",
    "log=nk.logging.RuntimeLog()\n",
    "gs.run(n_iter=600,out=log)\n",
    "\n",
    "symm_energy=vstate.expect(H)\n",
    "error=abs((symm_energy.mean-eig_vals[0])/eig_vals[0])\n",
    "print(\"Optimized energy and relative error: \",symm_energy,error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "0e2ed895",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f85b8188c70>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEGCAYAAAC+fkgiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABZrUlEQVR4nO29e7wcVZXo/13d55Vz8iYJBEJywiNAwiOBCKKCKIgZfygDv/GGUUcBGcDXoOOMiszPUUbm4kVF5jJXwjDj44oalQHUESWML9TwSEIeEIgSwis5QN4nOe/uXr8/qqpPdXVVdVV3VXefk/39fE5SXY9de1ft2muvtddeW1QVg8FgMBjqQabRGTAYDAbDoYMROgaDwWCoG0boGAwGg6FuGKFjMBgMhrphhI7BYDAY6kZLozPQ7MyYMUO7u7sbnQ2DwWAYM6xdu3aXqs70O2aETgW6u7tZs2ZNo7NhMBgMYwYReSHomDGvGQwGg6FuGKFjMBgMhrphhI7BYDAY6oYROgaDwWCoG0boGAwGg6FuGKFjMBgMhrphhI7BYDAY6oYROgaDwWCoG0bo1IHlK1azfMXqRmfDYDCkQDabZfHixcW/m2++ObG0169fz89+9rPE0msGTEQCg8FgqIEJEyawfv36VNJev349a9as4R3veEcq6TcCo+kYDAZDwuzfv58TTjiBLVu2APCXf/mX/Nu//RsAH/rQh1i6dCmLFi3iH//xH4vXPP7447zhDW/gtNNO48wzz2T//v187nOfY+XKlSxevJiVK1c2pCxJYzQdg8EwLhCRVNJV1dDjAwMDLF68uPj7+uuvZ/ny5dx+++1cfvnlXHfddezdu5e//uu/BuCmm25i+vTp5PN5zj//fDZu3MiJJ57I8uXLWblyJa973evo7e2ls7OTG2+8kTVr1nD77benUrZGYISOwWAw1ECQee1tb3sbP/zhD/nIRz7Chg0bivt/8IMfcOedd5LL5ejp6WHz5s2ICLNnz+Z1r3sdAJMnT65X9uuOEToGg2FcUEkjqTeFQoGnn36azs5O9u7dy5w5c9i2bRtf/vKXefzxx5k2bRqXX345g4ODjc5qXTFjOgaDwZACt956KyeddBLf/e53ueKKKxgZGaG3t5euri6mTJnCq6++ygMPPADACSecQE9PD48//jgABw4cIJfLMWnSJA4cONDIYiSO0XQMBoOhBrxjOsuWLeOKK67grrvu4rHHHmPSpEmce+65fPGLX+QLX/gCS5Ys4cQTT+Too4/mjW98IwBtbW2sXLmSj33sYwwMDDBhwgQeeugh3vKWt3DzzTezePHi4ljRWEeaTSVtNpYuXaq1LuLmzNFZec3ZSWTJYDAYmhoRWauqS/2ONcy8JiK3iMgzIrJRRO4Vkan2/m4RGRCR9fbfHRXS+aSIqIjMsH9PEZGfiMgGEXlKRK7wnD9ZRF4WkfHjDmIwGMYd43VSeSPHdFYBJ6vqqcAfgetdx7aq6mL779qgBETkaOBC4EXX7o8Am1X1NOA84Csi0uY6/k/AbxMqg8FgMBhi0DCho6oPqmrO/vkIMKeKZG4FPgW4bYQKTBLLaX8isAfIAYjIGcDhwIPV5ttgMBgM1dMs3mtXAg+4fs8XkSdE5Dcico7fBSJyMbBdVTd4Dt0OnATsADYB16lqQUQywFeAv6uUGRG5WkTWiMianTt3VlMeg8FgMPiQqveaiDwEHOFz6AZVvd8+5wYsTeRu+1gPMFdVd9uayX0iskhVe13pdgKfxTKteXk7sB54K3AssEpEHgbeD/xMVV+uNHNZVe8E7gTLkSBicQ0Gg8FQgVSFjqpeEHZcRC4HLgLOV9uNTlWHgCF7e62IbAUWAG4XsmOB+cAGW4DMAdaJyJnAFcDNdnrPisg24ETgbOAcEfkwltmtTUQOqupnkiqvwWBofow3aWNppPfaMqzxmHepar9r/0wRydrbxwDHA8+5r1XVTao6S1W7VbUbeBk4XVVfwXIqON++/nDgBOA5VX2vqs61z/874NtG4BgM6TFeva/8uOmmm1i0aBGnnnoqixcv5tFHH63LfS+//HI6OztLJpB+/OMfR0TYtWtX5HQ+//nP8+Uvf7nmc6LQyMmhtwPtWOYvgEdsT7VzgRtFZAQoANeq6h4AEbkLuENVwybO/BPwTRHZBAjwaVWN/vQNBoMhBqtXr+anP/0p69ato729nV27djE8PFy3+x933HHcf//9vO9976NQKPDLX/6So446qm73j0sjvdeOU9Wjva7RqnqPqi6y952uqj9xXXOVn8CxNZ5d9vYOVb1QVU9R1ZNV9Ts+539TVT+aZvkMBsOhQU9PDzNmzKC9vR2AGTNm8Mwzz/Dnf/7nxXNWrVrFJZdcAsDEiRP5+7//exYtWsQFF1zAY489xnnnnccxxxzDj3/8YwC++c1v8ruvf5rffO06uru7uf322/nqV7/KkiVLeP3rX8+ePXuKaV922WXFZQ9+/etf88Y3vpGWllF94qtf/Sonn3wyJ598Ml/72teK+2+66SYWLFjAm970puISDABbt25l2bJlnHHGGZxzzjk888wziT4vEwbHYDCMC77wk6fYvKO34nmbe6xzopj+Fh45mX9856LQcy688EJuvPFGFixYwAUXXMDy5ct5y1vewoc//GF27tzJzJkz+cY3vsGVV14JQF9fH29961u55ZZbuOSSS/iHf/gHVq1axebNm/nABz7Au971LgB6dzzH2274JivecxrHHXccX/rSl3jiiSf4xCc+wbe//W0+/vGPA7BgwQJ+/OMfs3fvXr73ve/xvve9rxjTbe3atXzjG9/g0UcfRVU566yzePOb30yhUOD73/8+69evJ5fLcfrpp3PGGWcAcPXVV3PHHXdw/PHH8+ijj/LhD3+YX/7ylxWfVVSM0DEYxhFmkLz+TJw4kbVr1/Lwww/zq1/9iuXLl3PzzTfzV3/1V3znO9/hiiuuYPXq1Xz7298GrDhry5YtA+CUU06hvb2d1tZWTjnlFJ5//vliujMXnE5rRxczZ85kypQpvPOd7yxes3HjxpI8XHrppXz/+9/n0UcfZcWKFcX9v/vd77jkkkvo6uoqnvfwww9TKBS45JJL6OzsBCgKuoMHD/KHP/yBd7/73cU0hoaGEn1eRugYANNYGcY+lTQShzTqejab5bzzzuO8887jlFNO4Vvf+hYrVqzgne98Jx0dHbz73e8umrxaW1uLC85lMpmiWS6TyZDL5UbTbBkNpBJ2HsDy5cs544wz+MAHPkAmU/2oSaFQYOrUqaktvw3NMznUYDD4cCh5gI1VtmzZwp/+9Kfi7/Xr1zNv3jyOPPJIjjzySL74xS9yxRVXhKRQO/PmzeOmm27iwx/+cMn+c845h/vuu4/+/n76+vq49957Oeecczj33HO57777GBgY4MCBA/zkJ9bQ+eTJk5k/fz4//OEPAWuNIvcCdElgNB2DwWCogYMHD/Kxj32Mffv20dLSwnHHHcedd94JwHvf+1527tzJSSedlHo+rrnmmrJ9p59+OpdffjlnnnkmAFdddRVLliwBLO3otNNOY9asWcUVSwHuvvtuPvShD/HFL36RkZERLrvsMk477bTE8mmWNqjAobK0wVjI46FI3PfSTO+xmfLipp75+uhHP8qSJUv44Ac/GPvaZn1+UQhb2sBoOnVic08vy1esHpMVyGAYT9TrGzzjjDPo6uriK1/5Sl3uB2NDUBmhkzLP7TzIo9v2MKHVDJ8ZDIcSa9eubXQWmhLTEqbMEy/uA2Akb8yYBoPBYISOwWAwGOqGMa+lyPIVq9l5INmJVQaDwTCWMZpOylRYusdgMBgOKYzQMRgMBkPdMELHMC4wM/cNhrGBEToGg8FgqBvGkSBlzJCOwWCohh37Bjg4lKt84hjDCB3DuGIszMg2GKLw0t6BRmchFYx5zWA4hDBjX4ZG0xChIyK3iMgzIrJRRO4Vkan2/m4RGRCR9fbfHRXS+aSIqIjMsH9PEZGfiMgGEXlKRK5wnTtXRB4UkadFZLOIdKdZRkMymEZy7NI3lOP53X2YoMIGN43SdFYBJ6vqqcAfgetdx7aq6mL779qgBETkaOBC4EXX7o8Am1X1NOA84Csi4qyE9G3gFlU9CTgTeC2x0oRhJuoYKjBeBevmnl5e7R2ibzjf6KwYmoiGjOmo6oOun48Af1FFMrcCnwLudycNTBJrWb6JwB4gJyILgRZVXWXf/2BVGR+HjMfG7lCmoMrgSKHR2SihkqZjxuEOLZphTOdK4AHX7/ki8oSI/EZEzvG7QEQuBrarqndJu9uBk4AdwCbgOlUtAAuAfSLyn3bat4hINvmiBJMrGBODIX1e2tPPpu37eXlvf6OzUsTU/OakURp2apqOiDwEHOFz6AZVvd8+5wYgB9xtH+sB5qrqbhE5A7hPRBapaq8r3U7gs1imNS9vB9YDbwWOBVaJyMNY5TwHWIJljlsJXA78e0DerwauBpg7d270Qvth7NmGOnJg0HKx3dM3zJxpnbGuTUvj0IQVr2bXjJo9f40mNaGjqheEHReRy4GLgPPV1r9VdQgYsrfXishWLC3FvXTnscB8YINlRWMOsE5EzgSuAG6203tWRLYBJwIvA+tV9Tn73vcBrydA6KjqncCdYK0cGrfsJWnVcnHCpPExmA+sOWkmxboQseNl6tKhQUPGdERkGdZ4zJtVtd+1fyawR1XzInIMcDzwnPtaVd0EzHJd8zywVFV3iciLwPnAwyJyOHCCff1eYKqIzFTVnViaUG1rUEfEKDqGqCTS6Np+K83gMSYIoOSbIC+G5qFRk0NvB9qxzF8Aj9ieaucCN4rICFAArlXVPQAichdwh6qGCYt/Ar4pIpuwPr9Pq+ou+/q/A/7bdjJYC/xbOkVLFq/N1fQCozPWe865fIEXdvdz5NSO2Nc2haZjyRwKTZCZsV4XxhON8l47LmD/PcA9AceuCtjf7dregf9YD7bn2qlx81orjf/corFtVx+7+4brft9GNgbN3hD9estOXukdZDgX3eVYioGXmqfmjSVNp9nrxHjAhMFJmyb84Pw+rNcSXGzOfLjJkCtYI/DV1KAmUC6KNFNexjv5grLr4BCqijTpHMFmcJke15jvrTkYixMwa2msm6mv0wzmtUOFF/f08/zufn61pT5z36vBaDop4/7chnPNNWlv+YrVbO7pZeHsyYHnnPL5XwCw6fNvr1e2UmUsaWGO11ec/qpzblSPsXqQN0InFrV0jnJ5q43pb+IoEEbTSRvX9/an10wgBEN0inIjjpmk6L2WeHaqZiyN6RjSxwidlHF/bo3s8TlajaE6GmGeq0Vb0SYy7DaD+/YhQ5OO47gxQidFrICHg8XfY6A+GJqQasxrYe385p7eogCthzC1LT5jclytGahGaDeznDdCJ2Xc7z5jhE5D2dzTO6a0PbemE9RgB+1vhkbHqe7VavhGSFk0w7tMEiN0UmbE5TwwMJwf9x9Ro8x4I/kCO/YNxOoVNrsQKtTgd9JM4yjN5NSQNPUQjOPt+RmhU0+MfS01tu3q46W942tN+aL3WhXVppnclMdbo1lvCjq+tD7jMn2I0sw9/Gpwlo7IF8ZPI1dNMRwBleZSGnHdzq+/ZxMTO1rGXZ2rF9XU52b+AozQaTBjad7IWMjrllcP1JxGUDmdAfh6lb+aeTrO2X7jKK/1DvL483voaLWWkooyTysJ4nrSjZcefVKMkz5UEWNeqyNBjYfbm8hgcBhta+Lb1/yEzq+37KSgMJyv7yTlcdZmlpH29zteNHcHo+mkiLOgVlqMBc2jVgqq9A6MFH+PxzIXCgoCGc/gTTVjOkWPMZ+GKuO4T8ZowxJ53g1sM6sVBu5yJ1Xnqk2nmZxCksBoOikxVjSXA4M5Nu9oXlv7S3v62fLqQda9uDfWdWPl+QM8/sJeNry0v2x/LcMyfo4E2QZ97ZWKMZIvsN/VsTCUUs3Kq808IdcInXpSB+e1arxcDjSxx9fgiPXF7e0bTsUdu38oV1NvOMq1Uc7zNXlV5Ulg/efnSODVpJIgrGxTO1uL25t7eukPqGfPvHKAZ1450NQNZSNJy7y2r3+YA4P1F/ZG6NSRseAwXYtrZppzdBrhAdxoN1WnzHHqzeu6pwPwr796tuxYGkInCpWESTMHp6yWJKM+pCV0trx6kM09tTvexMWM6YxjqrEhN8M6HCX5rnL5ZUf41eqZtbmntxhpOw5JjAMUG5sqxnT8tKSsPabj1YKcBvK1A4Ns29VP/3COzrb6Nw1K/TpmzT426O68xelwjYWOrRE6TU41QqAWb5rNPb2x7rdpe/lYRBqMF8OL3wC1Q76gJQKhGu1OQnwFKoVh2rHPihO488AQ8w5LommQwLz4UVDllf2DzJ7SwdOvWD1wv07DSL7AS3sGGBzJF92/m4nXegcTDe473syORug0MU9u389jz+/lhMMnMrWzrdHZ8SUp00jQnJEoASzHC0/u2F8cw4LRxiZewE+7ofd5Xo0zr/nv9wrdHfsG6dk/SEsF6fjy3gF2HhzinnUv896z5gWm2ygtZtvu/kTTM95rCSEit4jIMyKyUUTuFZGp9v5uERkQkfX23x0V0vmkiKiIzLB/TxGRn4jIBhF5SkSucJ37v+x9T4vIv0iD7UhhYyDLV6zmI3evA2BfBc+eqDbjgiojdZ6jkRzxPjynwU5yDkX/UC6RMaugPLkFDrht+dGrabimU7/qvnzFavYPDJftD3PcKEaVGGeNbK00UUSjRGikprMKuF5VcyLyJeB64NP2sa2qurhSAiJyNHAh8KJr90eAzar6ThGZCWwRkbuBpcAbgVPt834HvBn4dQJliUTV31JCla5vKM+6F/dx1vzp0W+dcANQbS90LLRDSc/wLw7pxJAVjmBxL6nh4B2QdjzKOttLm4Gkn3Xk5BK8r1uwP7VjP0O5Qtl78Qr/ldecTd9Qjr6hHLMmdySXmRqpJo6e867jfG/10hAbJnRU9UHXz0eAv6gimVuBTwH3u5MGJtlazERgD5Cz93cAbVhdx1bg1SrumRpBL12BvqEcrx0YolDQ0Ul+VfLC7j4KqhV7voMjeQ4O5WlrSU8hPjA4wsT2lorjSEnLnJ79AxwYzJWo+ml7qsX17quqhxui6aQ9sz2o/r68t5/WSJOE4mt2le4N1nMcySvPRli590l7zlo9hU6lxv6j313Hc7v6YnVmaolQnjbN4jJ9JfCA6/d8EXlCRH4jIuf4XSAiFwPbVXWD59DtwEnADmATcJ2qFlR1NfAroMf++4WqPp10QXzzSpUBpl3XbHn1AK8dGOK1A0M15+eV3iGGRgoVTW0D9njNcC5eDS4UlMGR8LGe5StW87av/obNPQcilama9nIkX+DAYI6RfIGX95Yue/DkdqtxqeXbTHtphLsffQGAixcfFfmasI5Es1tW1WfLy/IVq9nbX262i5d+cNr1Jq1pBn4myihm+DXP76nKWzMOqWo6IvIQcITPoRtU9X77nBuwNJG77WM9wFxV3S0iZwD3icgiVS2+GRHpBD6LZVrz8nZgPfBW4FhglYg8DMzCEkZz7PNWicg5qvqwT76vBq4GmDt3brxC+yGj8qPeQSMdvGYypXwMwUu1tvW+4TwbXq7s1TaUswRTkIBavmJ1MZRQ3KCRyqiw7Nk/OkDd0ZphcCTPlAmtodc3U0TkOIpt2Klh73NzT2/JBFV349SssQElYQfhZnrnbuLU/DcdP4P71+9oao+3VIWOql4QdlxELgcuAs5X+ymp6hAwZG+vFZGtwAJgjevSY4H5wAbbLDMHWCciZwJXADfb6T0rItuAE7HGbx5R1YP2vR8AzgbKhI6q3gncCbB06dKGvb0kP6lqXDjjmmPCZp37E72EtXxD7muHcwU2vLyfnTE1xoJq4mMdQY1c0YQaYUwnsKH0yWzSa+xEHsOKeNsmbifryv6BkdK4jVU8l2bWahvpvbYMazzmXara79o/U0Sy9vYxwPHAc+5rVXWTqs5S1W5V7QZeBk5X1VewnArOt68/HDjBvv5F4M0i0iIirVhCqC7mtSSp1gGpKqFTcG+XXp/GbP0w77DYuXdd4H5mjofUKz4D7aF5G84HmuJqeRZ+QtrJo1MEx2QWReMIe07uOpAveoqlHxGgWllSUGV7yGqwcbXftFm+YjX/444/xDrfqe/ud7unr9R8WE05ax2/SzMaRyPHdG4HJmGZudyu0ecCG0VkPfAj4FpV3QMgIneJyNIK6f4T8AYR2QT8N/BpVd1lp7UVa5xnA7BBVX+SdKF80ZL/EiPOmEI1i3q5K26ai4JFSTmuucB9tvjsi+s67ty+HmaLUeFg/R/H1TmssXGb117Y3edzP/ddqydonEI1XJvxBv0cHLHG4p4KCEj7b799ruFmv4JqohNBIRnX9mZeDqGR3mvHBey/B7gn4NhVAfu7Xds78BnrUdU8cE01ea03zkd7WJc1ITSs+jjaQZiJY/mK1Wzd1Rd43I07FIlfzzhJ4nxaYd/Q8hWrI/XW3WUYiekcUU9yhQJQ3Ux7Ryg6JXWbwP58yahDwsBIPtBjMC3BGveJO/kIqns9+wfZ3RfdqSCNYm3fN0D/cJ6+oRxd7S2JCOw9faWm3zj5dt5oHFNqvceYm8V7bdwzoS3LxPZ4Mn40jpbPPg8VTS+emhulwS/VdAqR7pN33SbJXmhcE8OgSwj5ta1JTZJNasKoG6eRLWpXFc49MJijZ/8AEO4q6268w3rTcfsXSTsavNpbu4dmvRiynXGKnZ4Kzy7Ss0pA03F/h3739BNKm3t6S65LCyN0UsLdW7S0B6F3YIS+oVxNA7r1nF/qzmYamk6RCEnX0kv183IaqbI8QyEaUnxHCn+8YzqoovafF0d4rtpsTTnzM6s4gtF9LBviEpeUpuN+7u40n9ox6tnYtB5jEZ+BM//IeQ9Rn1yYq3QSDkSV8t/IqA9G6NQJwaqQBYWnon5oPqpOVFtt2XhPVR4wtY3pVGxQYixkWanYnW2l5qgCLqHpp+lUaV4byVuD21GoVgPwG2N5tXeIg0P5opu5l5aM9SmHvSZ3utmMlPx2N1K19i+8wldVa+pBJxW9Z2Akn/gclLYWK3OO0HnmleSFqGI9w3390de+qdRJDDuehvbuxgT8rBcBH04uXwiJMFC+32kblq9YHatX7a1iQUJkYChPV0eLfa/SMZ1Cwb+3XU1+4lJLO+j3dH+ycUfV6X307nW0tWRSK2/OY14rqLLHnhA55Jlb5ZzTmpXiuUF4zWvucbD+oXzxQVXq2OzrH+bFPQOcclS4q7TTEA/nFfKj96r0LodzBatTkICweXL7/qq88wZG8pGiKDhmypEU7VKq1jPc8uoBfvvHnZy7YGbguc44XV41VJtK0zGoEkbopIT3lUrAsQG7Edk/MBI4YdF9vrtBCGobnMbQG1MrLvmSMR1l3Uv7QntIcTX2eI4EVuIFW/jVSi3WhT++diDV4Jkf/s5a7v/om3De/L3rthcbzqBsOw1kmMCIal6rJHSe29XHSF5Ta2jV/nOa/FomgfZV6Q6+aXsvLRlh+YrVbNq+n6OmTig7Z/mK1fTst1zvfVd+9TA4kq/Sq2y0zgdF73A06iPtfAbdxok4sGDWxCrykQxG6KSF66W3tWRKK79PjXjmlQOcNX96Sc/E71Nzu+5GMRDVoia7B6U//J21RYGT1AJpcbBMk0rfcJ6tu/rKZszHTq9EeJeaL5/a0Yuq0hIpXlg41TwrDfgfghuTUaETfJ67XQyLcuC+dmgkX3nlzwgelEHpl+wPOr9B83FyBbVMhcN5/lQhbpsjTPxymisoWaFilI4o9TjoXTjXOkKnknmtcXqOETqp4R5OaG/J+Kr41QzYRvFoSooSrcq13y8ysR+VPhDHJTwU1yQb5zvaH8O2HYT7m3zs+b3F7cGRAv3DeTJCsNBRkhntDUBDBIf3zTu/WrLOOjqlx90eSVF72cUGVJX1L+8P1IoGRvJlQVOTxMltmpagxDpQAXk8OJRj7Qt7mT2l+gCi1azPU+ldN3Iaj3EkaAB+PdiojDYI/sf39g+T12QqVelAc4S81X7LMqbZJsc9/cOsfWFvhbP98cu6d99oYx39wakGN4i5fKFip0IDQ+sE95qDUnTMfWGNzfceHV0BxHtWAfcYUumxoF7z7oPDoXkKImpnyzktVc/JKvGOlwTlcJ89Fuc8q6qo4mP+4ZqXwpOsNi8JYIROSkQyCXhOCQwB4zovTI0HijHFFDgwmKNvMP5gd3H8pKSyl98xyEQV+36eNPf0DbvSsxrTWpbF9suaN79+7VpB/Z0nij1wgt/D2hf3FbXboZE8jz+/p3isfyjH49v2cHAo7y9Y1Lvhc3MPTsPsmET9BZa6fwQS1rHZ3NNLzladqnW7HRop+Hv1BaQXJkjDFoVrNI9u28OfXvU3y4W519ciECp1Skdv0jixY4ROWnjeqfgcq25I0dnwv7qrraXkvLj9dwUODuXZ3z9S8rFv3Vka0aB/KMcal1kqSh1evmI1a57fU/zY/Ca/9g1Z9nPvgGltQs1PYJb+vuxOq+Fyl7NvKE/fsLWmUByczoMjyIbz6hprKdVuwrQZt+xxntWLez2mFk8Hwfu/X7pB9x1NMrxj4xDk0KFqlTfo+rj709B0hnL5RNMNqp6fvmcjkNxk5EoUJxZXOM9oOuOQ0I+6+H/wWbl8ged8Qtf4NQgDw/ni8gC1Dro6V+8bGI5lS6/UyMWl+JE6bryubzavGst5wK9BKMtvDZlVOz9RXKgPDgUHDnWnZ/1fOXPOEa+btaqPoK5QRim6TJfmA/zj/AVpOvlCskP/cUO6RKkbm7ZbDgKJh/zxPvJShb2EvPoHfK1EWDDOqMVp+jEdEfmKiCxKOzOHGmHr2XjdUb32dnelyRVG3VeTqkwj+cqBDN25r6bTGH6JuP5Nfga1N7/NNGpQTV6chtktTA96GjSvGTOIz933ZKR8BNWPqPnfeWAo0gKBeY/wy4W4ahdUGRiurMV4Qw3VSpCYrcpFuoZMRTWvefOb9jw7N1G9154G7hSRFuAbwPdUtXoj+yFApZf+au9gYPRcL284bgYPbOopSTeuiSJuRf7Eyie4/s9O8k1XtXyWuG8jFvmWWnSAcJKNlL6LsAmAfteWjdOoAlISsy0uiTVgnnccJdm8j9B51bN8QyWzXlA+yrDDa5RFTlANDRPkxavJt7VkfDticSYyDgznyRWUQoWVa9tbMgzlCuRVyUR0RXTWOXJrU0masZxGP2wOlRu3tuN4lDrvf3ffEH3DucDJ543UdCIJHVW9C7hLRE7AWiRto4j8Hvg3Vf1VmhkcjyhUDGnh/tC880jsnWlkzXUfjzZQ4X5R24XR7Ct99ljJ3oGR4kzqIGtELfZ3vyu9yW159UDNZiF1/lzPyt17jL+cgn9unMZm5TVnF/PraAPucvWWOZFEK11UYeftxffsH+SlvQPF6AhxCapiXvNaS1ZKLAGbe3rpG8yRzQrHz5oU6V5Ow6628Nzw8j5aKjT2bgFVFBC2W31QZzAJ893QSKFoVfj6r7ey8+BQoIv3o9tGHVZUlYGRQqyVZ+tB5Hk69sJqJ9p/u7DWpPlbEblGVS9LKX9jlmdeORB8UDVWPCl3tQ3zXtvc0+u79HP/UC62+cuKlxX9orAxHXePzGmYh/PKwSHLWSCfV/Z51lLxPp9aJr/7RS325lcZjQ5RibKGhdHYehA8C37di/uipV/hufuNWYyai0av9Zqh3ALdWSjMnW83Q7l8Ranj7Qg4Qq7adjbIFBWlw6FYJmFnWYCCEsl0B9Z5QdqUu9Pwnjsf4YcfekOkNB08SyOV4Qg9t3bjPTVqV6V/KOcbqM6vaP4encHHkiSS0BGRW4F3Yi2K9s+q+ph96EsisiWtzI1lBuzGv73F6gm564ISPCPcz666euvu0WvrVDGsOShxhE689CsNDntDn6x/cW/ZOTVFW6jxAYYJhihJh51SbKfc79pvJTrXb6dhftY1cz7n0azcXnBh9A6M8Mwrg7S1jA75qiq5gpZoA95X6IwhRTE5BS3y5kfJEhv5Qtl4Z/9QzveefoF1a3Gx9t4jr1CwBVtQ3r3jalDqLl3AMg13tY8GrA0bs6oW7/N26/SqiojUNH8wDlE1nY3AP6iq30pgZyaYn/FH4FceR9XRss20K4YSz2uoZjOC93qBdS/ujS7MYt6/puyqoinGXvN7x0Eyx8FPA467fINzX2d8zK1hFNRyfJnQmp7Da9A7ccvOOEWKqulEJbyjEKLKVLgWKJqaq7RMArY1IGLFDjvNayJOmqhCZwNwgmelwf3AC8ahID6q/prO4EievJZXPC3ZtqLHbt0ZHguqVrzL8Ma170e5xo239/rK/sGSfUl/ApU0nYqNm+u4OzJ3EoxqOuWZGMoVaMlI2eCwX+/4pv/aXNzuH8qVmF5Di1eUcG6hk3Y3BwoBuXKbeSvlwruSaJD5sKq2PeQZ5PLKC7v7AycSJ03c8D2Ok4dfh8atSDukGQQi6pfyf4DTsTQeAU4GngKmiMiHVPXBlPI3bvn4yvVl+8LcnqdMaGV33zDX/+emwBVIqxm7CcIyrwUf9x7zPTfBhqpiSJnY6dWWgPv0AlQV+SEw7Qq941xBaXMGwp08+LjKhk7TsX+02oPylbQqx6vMaaDCHle1Th/Bmo6/xpIrqK8npRuRyr36yPmrcOyV3sFIY7X9Q7nUV+j0Jt83nC8xr5ZMucgrbS1SWqdT7GRE1ZV3AEtUdamqngEsAZ4D3gb8r7g3FZFbROQZEdkoIveKyFR7f7eIDIjIevvvjgrpfFJEVERm2L+n2eltFJHHRORk17nLRGSLiDwrIp+Jm+daWDh7MtNdwS1zBS2xvXvx+8S8dn4//CqyBuyvxEt7+/nyg8HDdWVtdlS1PuL9ved527H+oVxN8wpq/aTKwujUmF5p2s5GyX+h+DX0ZauNOpvehjGgoYwtOyrUgWotR//7v5/1vUe+UB5p3avx1bIsgh8VJ58m0Fa3VLCxRar73g5HWV0a3XripX3ll6coFKMKnQWq+pTzQ1U3Ayeq6nNV3ncVcLKqngr8EbjedWyrqi62/64NSkBEjgYuBF507f4ssN5O9/3Abfa5WeBfgT8DFgJ/KSILq8x7VSRX9dM3c4zeKvq9/BqowaAYW1XcO7S3msAXErfHGzaxt9b77bBXJi3EkDr5onY0enKuoOUhfATashkmO2s3STLNcpphXrymoDgkOfTmV80quYdXI7fDNMUd+wbIa2k0A3ck8aj3rvTJ5EMWbKyVqEJns4h8XUTebP/9H3tfOxA7zryqPqiqjqh+BJgTNw3gVuBTlD7bhcAv7Xs8A3SLyOFYzg7PqupzqjoMfB+4uIp7Vo0kPPAcuzLHPV/9zTFBBKnjm3t6WfP8ntieZmURGUIyMJQrlAi9tlpGY5uAoZwtrItl8vbey72w8nlrrM89ObN0rpfrekl2oNgRwrEEd4z774i4PHhVxMhHNU8szDU56B71CKpd6Ra5grKzlsjYIUQVOh8AngU+bv89B1yOJXDeUmMergQecP2eLyJPiMhvROQcvwtE5GJgu6pu8BzaAFxqn3MmMA9LoB0FuGN9v2zv80VErhaRNSKyZufOnbELlDR+Pb2hGmbOR2W762Ov1NtMe6A57ENMc6lgP7xLRqdFYKlsmeoW5N9/vDyUvfuVOBGxBcvk5DeGE/S72QmLtBDU16uqxiRUx8NqT5I1OUhj6h/K8ZI3cKwPXpf7pKjoSGCbpn6mqm8BvuJziu/ghIg8BBzhc+gGVb3fPucGIAfcbR/rAeaq6m4ROQO4T0QWqWrx6xKRTiwz2oU+ad8M3CYi64FNwBNA7NZZVe8E7gRYunRpIvXAW/erTlStXmJ9mr1SwgaQ/cZ6leCAhl4hVqmhC+sZVxrYjkKU8DeVBqX9qKXXGs89vvyssgmwPq73QHnBwrTyGsozY2IbvYM5y4sq5rWOVhf3+SfZF/JLqmxsM7nbRc6DHwVVBkdqy01rAivn+lFR6KhqXkQKIjIljnu0ql4QdlxELgcuAs5X+2tQ1SFgyN5eKyJbgQXAGtelxwLzgQ22yWoOsE5EzlTVV7DC9CDWwW1YWtkE4GhXGnOA7VHLUgtStmERy7zhaZXD1n1PtNLHSMwZU6hWAFS6pg4euxWZ0Jq1IhPHuKaWbMcJyuP3fPyiLgil5jVf8RI28TVyjnxIcW5TEJM6WtjrCjlVUMfBwvrdBNUqFcI6O36v1601ZwhfK6pWorpMHwQ2icgqoDhBVFX/ppqbisgyrPGYN6tqv2v/TGCPLeiOAY7HEhpFVHUTMMt1zfPAUlXdZXvB9dvjNlcBv1XVXhF5HDheROZjCZvLgPdUk/dqKdN0YvfYnIHiaAIrbuPvPd97bUFHm0DfHl/CUsGbml+YkuItk1B1KiDU3+zkjkhQUSiXbfhon67Oi3Pa4ZPb2bFvsOy0NHAEXi03iVun/aKJ5woaWaCHuaA3C3G/vQIRTcQpFTiq0PlP+y8pbgfagVW2tvKI7al2LnCjiIxgPZtrVXUPgIjcBdyhqmsC0gQ4CfiWiCjWPKIPAqhqTkQ+CvwCyAL/4fbGqwdJum5GcUrISG3xyrwTICqNm9R7RWFVq4Lk8oVymSMpSaE6Sx3nm/dOeCw6tQWYy5zSB2o6JeeLT0Tv8Hx12BpfXNy3GQi4PpuRRBdXcz+D4nMJW2HOe71rO8pigI0giSz4RmJPIF0/okaZ/paITMAab6k51pqqHhew/x7gnoBjVwXs73Ztr8Yyx/md9zPgZ3Hzmhg1mNcODOVKLq9X2xcnjEjc3lYtDUuLq2Hyeq5VSxO0HSVs7ukNdEP2y6sTyt5NoHOHx7yW8Qjpap1CKl4lletuRqCzo4WDg6Px1IZzhaL5J1a9Ceh75LX2ReacRxRn2YXUSCkLaZUs6iJu7wTWAz+3fy8WkR+nlKdxSVKOBKnFRAoJveMmaOKa8+050XLrJRj9vvlZE9tj378Jmo7IVOoMOO/I79l0trVw4hGTS73XIr774vlBBypcmFadCLqtZV4r1XSAkqW0R/LRo6mrjgbrbKb6klpeUko4qnvC57HmuuwDUNX1wDGp5GicUq8x1IzAxPYWOgNC5dRKUDE05FjSVOxdjjWf3wBq/eb9NBYB/vjKgdE4bALHzpxYet+QGw/lCqGOLGGISOLz1cIo6Gi0dzfP7TzIiEtwR+3HFcc0o1vnxjQN1XSAER/PtUZ47Y5hSj+2WhSWsEvdjX+lRanc1NoUlGlgzoC1a/xhS9gaQzZJOWmm0bSNNVnmK3TsQR23D0ZXW7bknKFcdUIlTpUOajySfsZJNpzb7JVOGylwknLiiajbxU43ClG/8adE5D1AVkSOF5H/DfwhlRyNM7pndAE1mtc8UXbDQ8JAv927C+pURlkON251C5A5RUZ8Fmrzw8lzxXk7YWkA7SmG4K+VpBrWfm+IGzweiD4PSTxhb8RnX1qTbS89/aj6C26fohwKWkoSpGXJj/plfgxYhDWH5ntAL1ZkAkMF/unik1l5zdkVW5qww3n1mpSi1YagNCt9+NU4fxXcIWxrJBHfs4TNOJWCMCZJNWX3E/p+6yEVS6HeHekjUTwJQnodtQR4dVNtYzo4nG+q8Zy8wpPb95eteFvNK/XVoKrKVWUiCR1V7VfVG1T1dXak6RtUdbDylQbHhl0+T6dKDyGN8dFU2fC2t0TXElo8IfYTocaGUALatqQXIKs2m20Rnm8SDiN+Q18iXpf7+kkdETjh8EmRzm1LaDZ80FP0rkc0VvA6kfiNr1Xz2ff5aM0NdSQQkQUicqeIPCgiv3T+0snS+MB575ngjpv/BRGILHOC9le4VyZGrXVMdX0et+5aqDaduNdVihBcKx0xhHdc8ho8LnJYVzutLRmGfVyuHVOaU4c+fN6xaWXRF7+65ScC47yZ2Kbgai7CnqXfLGpOHajdsdyfqF/FD7HimP0D8Peuv3HP2rVri143cf5yfXsBeP1ZZyIiPP5//ycAhcGD5AcOMLy/NJDo8K7ygI0OQ688S37AGoT//Yrr+c1tHx+9bucLZecf2PIIT35hGdvvu8U3vV2//2HZvv6XRleZ3Pv80xz40+Psffg7AOT797uu/UFxe+f9N7P9x18FYKDnWQZf3cbuP1hpH9xixcp6acVV5JWS6MdB9G35Pc/dfBG5wcrBCL2M7H6J4d0vA7D+h7fxx9veX3bO89/5bMnv/p3b2X5X4OoZJexZ93OeuvHPSvb1bvhFye+hV54t+b3tro+x/w/fL0tr+49vLW4XhkrLOrK3h6Edf2TzP/952XWat8xLCvQ+di/bV4xOXdu/+bc8+YVlaG6ELQ/+X/pee5nCiGWM2PEfH0Nz1njaS489yI5Nvy9ed8Xll/ODa9/Aa/ffHFT0mnDyDPBPX/wi93ykPIZv7uDe4nbPA/+HJ7+wjAOvjtbr/Wt+zJNfWMZzN1+EFqI7OeQHSiObb/uSdf3O336XkaH4dQzg+Vv/B3lXfmvhhS9dxIEnqp86mOv1D0bc/1Llee8Dz2/w3b/rp1+1npMWeOpn36yq7avkoRjVrzanql+PeK4Bl3mk+L/joiroyCDSNqH0/IDVEQEkky2qJyIZWg4bXQlC8z6D82qlVRj2t4AOv+azDJK6XEgLI5BtAcnYv10fekno4kIx1opkW9H8CPt+9x1ye7dTGB6kc8HZSMvo4nVBFAYPkumYCM59tArHSMmMXqegOZ+w7N5nrPmyff1/XE3ngrOj3dNzbcHTyGmh4N9IuvZ5jxeG+pDWNshkvVeVqKiTz7yEAxt+7jpmvStpaUXzObSQI9PaYd8kj+aGrGOFHJms67N3nllI/asFzQ0jfvdzn+Oue45gdc+vyWScDetbiIz1vAZfepKDmx6y0h0ZQlo7is/Ly/CuF2ibMS84xUzGqqtNQLZzqu/+/ME9Fa/NtLX77i/Wx0IBCXhGtSKR4niJfB54DbgXOyAngBOiZjyzdOlSXbMmLPKOP/M/818o8F9/8yYWHTmF87/ya7bu7As8PyPRQsm0ZKVsdcSyczLChLYsuXyhbJCxsy1LNiMc8Cyv7A4/4mw7Zhi3OaYtKwy77t/RmikuaJYR6GpvYeHsyWx4eR+DI4Vi8MAoOPkeGM5XNdPbeYbzpneyu2+obAEzJ+3i79YsCw6fyIaXRzW5rvasr327NSucPncaj24brfJeh4eWjJTku9O+n7ck7mcWWBbizUloyQgFeyXNOVMnsKd/2DdUzcyJ7Qzm8sX3f9tli/nuoy/y+LY9qc+B+Ju3Hsej2/aUPEMofY4dLRlaWzLk8lqcYyPAxI4W+gZzVeVxzrQJ7B8YKS4THeaoctIRk3g6xLXfeS9JBVoKSmdSRwtzp3fy1I5461CBNSYbxbIQlB8nhNbsKR2svv786tIRWauqS/2ORdV0PmD/7zapKWaCaEWijI90tGTIFZSCasXKXEngVCKKu3QSVHUX+6KaY28F3Nxvt3dfLU/H+6qD0pozdQLPhnRAqsEbYSBsPK90FYN064NffQ6t484YaMxsOW7yta7oWu9A2E6+hzyL4J10xKSq300t33h7a4bWbKasU5okUWOvzU8tB+OcKGFhslmpOoZTWIDEes7+Bn9NLUqp3BpVEjhzTyqeJ9TUyrRkhY7WbPEDlai+5iH3TKIHbdnVg2/ttm6k3QfpbMsyUtASrytvGdtq6Jm7ac36C504RaxUb5LWBkWsfI/kCiUBeqv9dh1NxUtWygMAi1gasntelvvShixXLSKfcm2/23Psn1PJ0XjBfntF77UKDY37cJyGIFSYZaLHtp7kCpvjXKPeHRDaYGYzwsLZk+3T0hd4QXe46pxjivmodH2ZplNDvsuujKBxpdHoW5qOf8LeVUOzdnlTm0wbse6EYU2ITq4BrDSxOAq1vjbBEgQAC2dPTixsVWB9CvAaDNOKGjVP5zLX9vWeY8sSzsu4xNuIZcTq/SWFpQ4LbQHuv5EbtSq+orDKc4wdiSEOSVVyx6TZ1V7hOUv1io5QOt8mqIcZdG2aZAgul4gnCGaEBxBlXlEtRHGZjmpy9ctpnHdcTX2opuOQZh0Ieqd+z+aoqaUOTXOnd4122Bo0T8evPgT9NoQQ9iDdxJ35nhHLxNPeOtrAVvPhJP0ynTzEqbe1rubotEtOI1BpPM1P04lKRkbTz4o10D1/RkSvpgg3jW3G8QiSEs3ZfWvPmE6UBrOmcS5P3qq9QcrzoUevT/BLCFQ6KM/nwtmTE4s7GERXR0vZPco1/dH9jdJ0wmKvHELTpOJTNFE1INxINS8mAStICXEmmDpndrRmIpnFKqYXwaxl/Rbf91L9x1/65MManSRxGghHUFnlDzGvlQidov034Vz53RtWXnM2nR4NNEpDH8e6lsErTGPUxQQfQ0tWQjs2nSlGg/filwelXCv6wrsWlRxPg0olPk1EerHyPMHexv7dkVKexhXetejV/cOHXF5rr/hVSZ3wMScnycMntfPS3oGSY5W2o1KLZ547j1HHZYLGPhzvL3dDl6TXn7uUE1qzVS8V4JDxOJN4NZ3O9mzRfVykdKZ5Ju3udVzsjLv7akrlKu3tNGUyQiGvZcfqiWNy9VoG08hPmEk1EJ+HWlLPG+FIoKpZVZ2sqpNUtcXedn63ppKjcUJ7i9WbcxaIqtQQVlsRF86eXFE7iGKxO+f4GVXmwJ/WKsYBanUHdwjUsjy7vdGVg+g+rJPT5ky1Ard6CBoEFvtYxQH6RLRKz++QZDMe81pQbMCkiOxFaDP/sC7f+jycK1TvORZnTKf6S30SK7/6jLlT6eqoj3ZTiWym3MsxmxFWXnM2rdlMw9fTMVRJ2URvTe4DDxqs72rPsunzb49lqnL3cMLyF7U3lYnYoLs5ZmY854NazVeZENuHe7f72XS1ZwM/Gqej4fRwwxrcDNYAfRJjCF4h652ns3D2ZJfrvpQ4EkQxgyZlcoqSTlIKpVfrjuwh5rpwYogjirsscRrnloQCmcbG57n6BYZvsVXfw6e0c87xM1PJSkOegIjcIiLPiMhGEblXRKba+7tFZEBE1tt/d1RI55MioiIyw/49zU5vo4g8JiIn2/uPFpFfichmEXlKRK5Lu4yOx09bS+WeZLVuutM6/ZXNaiKaZCPnIfk+sYiljU3vKg2ZU+2dMhlr7KCS0PVqOn7ebjMntnGYK18nHzklsKfa1pLhdd3TIuVbxJo5nkTPulzolMa/Kt0ubSQbEWzZW9/9xhMbFVjT2wmI6rnnV3caZdbzw8mLezytoFr2LpzOSVaC5//VSqM0nVXAyap6KvBHSt2xt6rqYvsvMBqjiBwNXAi86Nr9WWC9ne77gdvs/Tngk6q6EHg98BERWZhccco5dmYX82d0cdysSXZ+g89dOHsyMydVMUQWkKjfipEZCf+ASmy5FcZ3gggSgkniLkNHgOt5VCcG71l+1x0xZYJvp6DTDvfjd2/36XOndfreu3i9fW5Ha4YM5fNloqz+6h2Xcd/fGbx3//Z1JAghK0I2I7EE1JFTOmjLRnQMKclv8mqVs7Xp82+P5R0a1fyKxnOcSYMJ7dlImtyiI6fQYdcxvwnZzvPJZKQ4NJA0DRE6qvqgqjpxFh4B5oSdH8CtwKco7bgtBH5p3+MZoFtEDlfVHlVdZ+8/ADwNHFVt/qPQms0wa1J5UL1MRjjJ9SG6XRTj4jWhOC/Tr4PS1d4Suk6O3/oiQnkP7oo3dhe3yxp914e38pqzI/USw0wYlQgygUU5H8JNYM7zqNmBIOLlAr4REqI8Q6+WmmFU8HnrV/m8McuGf2yIaTOTkchzy7JimX2Pnt7JkrlTS46FTVh1+N/vWcLKa87myKkJ+imJ76b/qZ4TjvTMY4lzuzA55NdpmVDDtxA1T1BaX/zGdJxz/BYBTIJmGNO5EnjA9Xu+iDwhIr8RkfI46ICIXAxsV9UNnkMbgEvtc84E5uERaCLSDSwBHg3KkIhcLSJrRGTNzp3+4cPj8s+XngpQDKTp4F2fPg5BldrRdFZeczZLu6eX9ICCHA+CzGvOPKDi75BGuBpTkbcR9BusD0rTGbvIQLH3BjE0nZDzOloydLZlfQV1FOeNqFR6ZnEW1CumISENvGf36OuML1z98tbWkmGmT2crDPednTGPZDvZ1XccolSlMTF3xKv++v8oatYZEV+LSRKkJnRE5CERedLn72LXOTdgmb7utnf1AHNVdQnwt8B3RWSyJ91OLDPa53xuezMwVUTWYy2x/QSQd107EbgH+LiqBoZvVdU77RVSl86cmcxgWmnvwqcm27tas8KE1miCKOh7qFRZ/Br2sB59NiPFyujO+/GzJtJ92GgPuaxYCVschNIxCBFhYkcLXR0tJfmq1rxWcsw2KdVKpRSKA8sidLa3lGkcjnOCg58jg3dwWsSl4fiMn7g7Opmig0FlvNUqsUgFrps7db/WBq9kaYTqslLx2qDSV5vzWmqb1+oRdtxvDM3BqfMv7e0viwaeFKkJHVW9QFVP9vm7H0BELgcuAt6rdg1R1SFV3W1vrwW2Ags8SR8LzAc2iMjzWJrMOhE5QlV7VfUKVV2MNaYzE3jOvl8rlsC5W1X/M61yB1HP+RDVaMXF2fUZYebE4DVwytrhGKaLJEjU28o+L8wjLQ0+9fYTiw3D/MM6mTe9s+jW7i2ft2H3jhlNaM2UaRwljgOee193wQJOPGJ0yWjnff7LXy6pmG+/anXiEZNCTXNRcOfRMeP17B8sOxaHINNQEnVUsMxjldzhxfN/PfGuWhuUh7IxHXdD1aAwOKkgIsuwxmPepar9rv0zRSRrbx8DHI8tNBxUdZOqzlLVblXtBl4GTlfVV0Rkqog4LeZVwG9VtVesr/DfgadV9aupF9CHNJYUCDIPHTmlsj3ce6kzYVAIN21kSjQ2b68pWhmjmoz8Uotyi6i9ZCcpb0NeDQtnTy5qkKUzzV0aGNCWzfC67unFfS3ZDEdM6Sg+u0pCx0vGZ4DfPVHQW66CaolmJK7Ohh/uvX6nTJnQWuLd5/cgnbTPOmZ62TEv1ZgTfdMJCgtVab6c97jP+V3t2YgOEpb22pZWMNUQyubJeSwEDt4vxR0Vf7zN07kdmASs8rhGnwtstM1jPwKudRaKE5G7RMR3USAXJwFPisgW4M8AxzX6jcBfAW91uWO/I9kihVPJJfk9Z80tbletnos1kFvN4KfTTntdar2UmLcq9OGCjrbGjC9XkqYI7fZ4S9D9gjS9k47wNhTl+Ugq4q8jhEpM6QJL5k5lSoiXn7dhr9QIC+XvwfK6chqP0mNegewIOb+xOvck04WzJ3P63Gn+zz1EswJLsC6dN43XH3NYMS1nnlJQWkF1MKjmTGgrFQRJdfJCUwn4pjtaM7Rlpai9Jt3dnNTRUjbh+0RP3fYKRT8DxcLZk8o0wqLQkQYtbZAWqnqcqh7tdY1W1XtUdZG973RV/YnrmqtUtWwJT1vj2WVvr1bVBap6gqpeqqp77f2/U1VR1VNd96x+cfIqCBqAH50RHvyxpaGee+vTYV1tdLZlfSecLpw9mckTrIbS3cCseP8Z3HTJKaP5jJBRv/kMcWNQtbVkyhoV9wRI98fibiS/e/XrS64RiTafpxb+v4vieeZ7NZ0Xdpcu9lbWmPhoaXd9YGlx34h3SW3PHC7nMfp1itxjQ+A/gz0q3vfldJCsbZ9E3VETfPJbhr2/Yl2q0JAGJV8yedrO7/zD/N3hMyK0t2ZZdOSUWHXLSTeKmzx4ArlWTtyTR5jU0VrmQj5azvSMgs3gvXZIUEnTqYebv9sE5KW1JcMpR01hameb73fZOzAClDaKWZFSc5u3ovoUytuoJukJNsUWjO5Jbe7GohFzKeLes9y6U/l67zmtrkmn3quDNZ2I+Yt2Ws1oTH3/hMMnBR5z59m7eFzYGGUlap1T5O3wLJw9maXzppV4YsbITFWHO9tKI0+XaDrxcxEJI3TqRCV1P4mPuZaotX928hHFbUcb8R1TcW17e75lDaZnO2jNn6Rw7h80k9r7CurRgM6ze8NR7+U3j6YS3nK1ZjKBnoReoeOcFz0iRTjVpJLEnasxg0GNno4pVCDrmwrLr//+WrLivl1x0UkYX44EhyJh81tg9MUf1tVe4oYM6TeOC2dPYp7rnod1tdHVlvX9VksqaEY42j3jXq3em582NbGjhfZWy+6+6MgpCeS5XENyshZkQSlzH3b9DIowUCtHTp3A0nnRQuP44fcOOttLbfrec6wxM09Z7f87PO74RU0nKS0w5coa1KkSe5JrJXfhkqrhE3zMeXLehrGS2Tu07lT5bP2+o66Q8oeNRwblwFpqwtJ2jpnRRWfb6PVG0xnjBMkcx+Q1+tHHNSxUJq4JS0RCxqBGt7MiTOtqKw52D+eTWUE+6jylcqzMRQ/fUXsLGWVMKO5YSCxbPeVmHr/7tbdYrtUXnXpk2blQuVMUlTAzVxzijmF/969fH3jskxeeUGzA3c+2JVse++4H176BiR0tgd5/DiuvOZtb3n1avEzaVPoe/czg86Z3MsE1PmlZNcq/k0ATeiXzm1AyqVdEUmiJLIzQqRNxzBdpeY3Exa/377fGjNNeDeWSETqnHDWZie3ZkjA75X33ciqZ18rO99mXtmNBEHE6Bt4erVteWBpq+WiDiJQ4YIjnWj/zb1CvuTUrgQ1HlFru9K6jEBRfzy9PQbiL1tZixbib1NFS0eS9cPZkPr3sxMDjUQf8w84Ksgx4mTWpPfL9vGQzwqK4HU8w5rWxjrsn6VfJKrmKVkvUSh0Vf6Fj/T/BMwBavUmp3K5txUmLmseIT7H+fgWpUBKNwTUQ7CZIqBXn6XgumDvdcrv3djxErLBI3vSLYwGeA2kLcb9IGV7cpsNFIVHC/SgZs/SmGzK3Kcxpx49K32mUuh+UxoTWTOQ5dG6MeW2MU9F7zWdfe0uGOVMnROrt1auH7q6ITpkyGavne/R0jwtpjY16pcu9ZXbOr0XTiUrSwtxN3Pahlrn3ft5rE9qyzJ4Sb67XhLYsHT6u7JHwuSRqg7foyMkcP2ti6Dm1RANx64zeokW1XgStexUrHwmMucWps2k6ehqhUyeiOhK4OXp6J0dNiz/RMwmimHvcZRKpfjA6isCMswBXQguQNgVBT3Th7MmhK8IWxyQqpO83TyeKGcerAWVEYq8WG/qanMnKFdLoaM2WrMPkV5dqabC9l2Zl1GkgqoD1phGn8Q8LOJuEIAojLSu/ETp1IrLLtLqiA0RIN80etx+qytSI6+bU/EnETGDUe200yra3ATp97tSi/T+pbzbpd1AiYKNk0qd1iFo0d8w9L0nOoapEW0umZKDcGcQ+btbEWGvg+FGLZ16o2S6yVlf9/RcdOZmFs8udM2qxbHjHBP283owjwTigkiruuCz7mdLKBoUD0qiHiU2xoksvOXpqqvdZOHsyv/rkeTGvsp5MWOy1VnfcsYAn2ShngiQplizqPJRiZIzG0N6S8dWw/vW9p5fGdouA9/2FyYZK5XUicvudF/hNu3avvOZsvnbZ4gp3CaY1m2FSR+VOXtIdT+NIMA6oZFd+y4mzWDh7Eoe73RbLNiy62rM1LX5WFUU1wuo5Rglrf1JIw93VlvUNieMm7viAE/Jm7vRwG7oT8PKuD1QK5WdRT23S22BGcYr4mwusQOzOADbAh95yXOx7HzuzK/JibUkS9patqMelZ8Rdg6oWTeeMedOYP6PLN0pA1LGiar3O/IjroADVmxeNI8EYJ8qg46SO1lK1Nqgj5ePdFZd69Gh/EPJxWNGRw3PREnMEeMbENhYdOZllrugK/ukKZ82fzvEJzSlxs/Kas9n0+bdXbBiSFGRH+ayyWU1DN2Nie2LzdSrhmO6yIr5OC46sdc85cnJ28lHxJhfH/VTc5iYRYdakdt/vzV0/w96lu57X0xReC2mGwak9nK4hEtV49Yjn/7SIU7ni2HlFhK62LEdM6WDGxHY29wSum+f7Mca15YsIExOIEN0IghqjsKftDOafMW86i46cXBIctGJ9EyvxtFaHjEJrNsPS7mmh5ySpJfgi9oPwEEWjiOq9VuuYVK1UpR1BalJnbH6hYxBvT2nGxDZ2Hxz2PfelPQOJ37/WHtboIH2868J6pZU+BnfD6TY5hQmvJPHLX717qlHdvye2l66eGrWT0yTzkMtwstWSFf7f0+dw+6+eLSnThNZMaQimELyRtZMiqiKeVFy7+mIcCcYdx86cyMQKk9TCTGjV9F7i4mcCqmcblXov10Xaz7PaYKy5Kv2/Kz0752iU8Q6/ehD0vJIe/2rJZIomvykTWoudj5ZshmkRHQxyIVKnlugfQYLdGwqoXibLqERxlMlk0ovKboROHWnJSHGmNwR/uN4ehoiEzslIm5XXnM3UCdYHXs+ecZTeelqeZvV2RQ8iqqbjpVKDMaEty7zDOpkbsCZMs5B1rYwatw106k9SJkRvxyFqo5xW5ynN+tl9WBeLU/JQNea1OnLGvHD7tRenqi6cPZnNPb0cGMwlnykiqv/2KWmp3L639MnXymvOZvmK1WX7gLL9Y5nOtizHz5rIpu37A88J086ctxT0ZjMiHDG58rLm1RK3QQzqOLRkpGbzVFgcWstVWxkYiW+Di7LENzRO02nJCLkqOy1pYoROM1OHutrekgkMl+6XlWYdA2h24mpj2YwV48zbZkQVsPUKGuvk55TP/yLRdI+aOoHt+wYsL8cqG23nKq+26LyLzT29iIg10F+N0Amcp1O6v55mYjdVLQZXB4zQaWLCQqAkRVh0Xj+MzGkccTzzwjq4ftpi0Hl+v/2uTdrEOWfaBObYIaCqGVtYec3ZLLnxQfb2j1RtoqxEWMDPkvMa5EiQdpicammI0BGRW4B3AsPAVuAKVd0nIt3A08AW+9RHVPXakHQ+CXwZmKmqu0RkGvAfwLHAIHClqj7pOj8LrAG2q+pFyZcsaaxK4zQSj27b0+D8jE8aOXYT9d5d7VlOOiL6vKJGukJHIY5JNFtlh91pdKOsr9TVnuXEhOZt/d8PnlUS0qeqIKjjmEZpOquA61U1JyJfAq4HPm0f26qqiyslICJHAxcCL7p2fxZYr6qXiMiJwL8C57uOX4cl1FKPcZJEQ1ZNRyX1BjSBxixqT3u8Uc27OXxSO68eGCLjWVivUlpO537yhGhx8poZR1M49/iZfHX54sh1x3lahQBNx21my4gUI1WEEcXLsTxIqBE6bhoidFT1QdfPR4C/qCKZW4FPAfe79i0Ebrbv8YyIdIvI4ar6qojMAf4f4Cbgb6vLeX1pxqpa7/7zSUdMKgu5E3edku7P/FfS2aob3TO66B/Jx77OGdNphjpUa0fIETpx697MSe3s7hvm9LnxHHj8iGM+9JrTms1lGkrfSb07gc0w0nQl8IDr93wReUJEfiMi5/hdICIXY5nINngObQAutc85E5gHzLGPfQ1LSFUcMRSRq0VkjYis2blzZ6zCjBWOnTUxNBCilxkTLZfpNx43o+xY2NrsDlHDw3iZPKGVjqqXrz50aXbzWhyqNU9NmdDKWfOn190tvAllTFORmqYjIg8BfkGwblDV++1zbgBywN32sR5grqruFpEzgPtEZJGqFqegi0gnlhntQp+0bwZuE5H1wCbgCSAvIhcBr6nqWhE5r1LeVfVO4E6ApUuX1v/rdUKvpVh5D+tq49Ve/9fvJxgmdVgf8LEzwxfMMjQHTgTt8dAA1rsMSWlmSd87zgRmR3upV/SOOKQmdFT1grDjInI5cBFwvtq2AFUdAobs7bUishVYgDX473AsMB/YYA8UzgHWiciZqvoKcIWdvgDbgOeA5cC7ROQdQAcwWUS+o6rvS6i4daElI0WPnqRwV+RmmAxpSIZ3nDKb//nA0xw1tTGLACaJY5767R+bz+pw/KyJTPBo4kEypxbhOZ6+zUZ5ry3DMnW9WVX7XftnAntUNS8ixwDHYwmNIqq6CZjluuZ5YKntvTYV6FfVYeAq4Le2lnS9/Yet6fxdMwuc0Yl9pbV0QluWw1Oc0NcsjKcPDOKXJwkbe1tLhu7Dal8m2Y96v59GuRxHYbpPKB4/V+UTDi8XTvUgqnZUz3faKO+124F2YJX9ghzX6HOBG0VkBGvs5VpV3QMgIncBd6jqmoA0AU4CviUiCjwFfDDFMhjqSC0fxZKjp9Y1kkIzMpYFeZLeX97nUK8B9Kmd8Raii0pGxp5LdqO813xXmFLVe4B7Ao5dFbC/27W9GsscF3bvXwO/jpbTxtLEHbwSmt0FOsqCc+OBegSBbQTN4v3VjM92aczQWs2AiUjQxHztsiXMn5GOiSRpmvGDNFRHNcIrzfcfJnIO9XrXrFEHwjBCp4kZe9Wp+Wl2razRNEpbCrvnGGxXDSEYoWMwGAxNynjU5IzQaSKcCnbCP1hzZYN6eI2uiI2+v8FQT0x9TxYjdAyGccCh1DAaE+nYxgidJsY7T2c8MRYiOx9qmOdiqAdG6DQI84EbopLGctzNgvkOmpM034sROk3IEZM7eGFPP4dNTGdCWT2pZ6NiGjBDNYxnod6MGKHThBwxpYMjpnSULCNt7NiGQxU/hxrTwRi7GKFTR8yHYohKnJU1Dc3BSUdMom84/tpHadGs7Y0ROgbDIUizNkhjmckTWsfFSq1pY4TOIcx4jdVlGPvUq14ajbL+HBqREA0Gg8HQFBhNx9DUGE3MYBhfGKFjMDQxxvxTH0znpn4YoWMwjAFMoxgf88yaEyN0DJEwH7DBYEgCI3QOUYwQMRhqx3xH8TFCZwxhKrjhUGTudGv13K72bINzYkiChrhMi8gtIvKMiGwUkXtFZKq9v1tEBkRkvf13R4V0PikiKiIz7N/T7PQ2ishjInKy69ypIvIj+75Pi4hpwQ8RVl5zthHYY5gz5k3jlKMmc8TkjkZnxZAAjZqnswo4WVVPBf4IXO86tlVVF9t/1wYlICJHAxcCL7p2fxZYb6f7fuA217HbgJ+r6onAacDTyRTFYDCkTWdbC2LWrR4XNMS8pqoPun4+AvxFFcncCnwKuN+1byFws32PZ2zN6XBgEDgXuNw+NgwMV3FPg2HMYLS75iGtdzEW33EzRCS4EnjA9Xu+iDwhIr8RkXP8LhCRi4HtqrrBc2gDcKl9zpnAPGAOMB/YCXzDTvsuEekKypCIXC0ia0Rkzc6dO6svmcFgMBhKSE3oiMhDIvKkz9/FrnNuAHLA3fauHmCuqi4B/hb4rohM9qTbiWVG+5zPbW8GporIeuBjwBNAHkujOx34up12H/CZoLyr6p2qulRVl86cObOq8hsMBoOhnNTMa6p6QdhxEbkcuAg4X1XVvmYIGLK314rIVmABsMZ16bFYmssG28Y7B1gnImeq6ivAFXb6AmwDngM6gZdV9VE7jR8RInQMBoPBkA4NGdMRkWVY4zFvVtV+1/6ZwB5VzYvIMcDxWEKjiKpuAma5rnkeWKqqu2wvuH57zOYq4Leq2gv0ishLInKCqm4Bzgc2p1pIg8FgMJTRqHk6twPtwCpbW3nE9lQ7F7hRREaAAnCtqu4BEJG7gDtUdU1AmgAnAd8SEQWeAj7oOvYx4G4RacMSZFckXCaDwWAwVKBR3mvHBey/B7gn4NhVAfu7XdurscxxfuetB5bGzKrBYDAYEqQZvNcMBoPBcIhghI7BYDAY6oaJvdaEjMUJXwaDwRAFo+kYDAaDoW4YTcdgMBgSwFgoomE0HYPBYDDUDSN0DAaDwVA3jNAxGAwGQ90wQsdgMBgMdcMIHYPBYDDUDSN0DAaDwVA3jNAxGAwGQ90wQsdgMBgMdcNMDjUYDE2PmXg5fjCajsFgMBjqhhE6BoPBYKgbRugYDAaDoW4YoWMwGAyGumGEjsFgMBjqRkOEjojcIiLPiMhGEblXRKba+7tFZEBE1tt/d1RI55MioiIyw/49zU5vo4g8JiInu879hIg8JSJPisj3RKQj1UIaDAaDoYxGaTqrgJNV9VTgj8D1rmNbVXWx/XdtUAIicjRwIfCia/dngfV2uu8HbrPPPQr4G2Cpqp4MZIHLkiyQwWAwGCrTEKGjqg+qas7++Qgwp4pkbgU+Bahr30Lgl/Y9ngG6ReRw+1gLMEFEWoBOYEc1eTcYDAZD9TTDmM6VwAOu3/NF5AkR+Y2InON3gYhcDGxX1Q2eQxuAS+1zzgTmAXNUdTvwZSytqAfYr6oPBmVIRK4WkTUismbnzp1VF8xgMBgMpYiqVj6rmoRFHgKO8Dl0g6reb59zA7AUuFRVVUTagYmqultEzgDuAxapaq8r3U7gV8CFqrpfRJ7HMpvtEpHJWCa1JcAm4ETgr4EXgHuA5cA+4IfAj1T1OxHKsdO+fiwwA9jV6EykhCnb2GU8l288lw2qL988VZ3pdyC1MDiqekHYcRG5HLgIOF9tyaeqQ8CQvb1WRLYCC4A1rkuPBeYDG0QELNPcOhE5U1VfAa6w0xdgG/Ac8HZgm6rutI/9J/AGoKLQCXpwzYiIrFHVpY3ORxqYso1dxnP5xnPZIJ3yNST2mogswxqPebOq9rv2zwT2qGpeRI4BjscSGkVUdRMwy3XN84xqOlOBflUdBq4CfquqvSLyIvB6W0saAM6nVJAZDAaDoQ40akzndmASsMrjGn0usFFE1gM/Aq5V1T0AInKXiFSSuCcBT4rIFuDPgOsAVPVRO711WGa3DHBnskUyGAwGQyVSG9Mx1B8RuVpVx6UwNWUbu4zn8o3nskE65TNCx2AwGAx1oxlcpg0Gg8FwiGCEjsFgMBjqhhE6YwQR+Q8ReU1EnnTtmy4iq0TkT/b/0+z9IiL/IiLP2nHoTm9czisjIkeLyK9EZLMdH+86e/94KV+HHQtwg12+L9j754vIo3Y5VopIm72/3f79rH28u6EFiICIZO1J3T+1f4+nsj0vIptsp6c19r7xUjenisiPxIqF+bSInJ122YzQGTt8E1jm2fcZ4L9V9Xjgv+3fYHnuHW//XQ18vU55rJYc8ElVXQi8HviIiCxk/JRvCHirqp4GLAaWicjrgS8Bt6rqccBe4IP2+R8E9tr7b7XPa3auA552/R5PZQN4ix0P0vGgHS918zbg56p6InAa1jtMt2yqav7GyB/QDTzp+r0FmG1vzwa22NsrgL/0O28s/AH3A28bj+XDivu3DjgLa6Z3i73/bOAX9vYvgLPt7Rb7PGl03kPKNMdunN4K/BSQ8VI2O5/PAzM8+8Z83QSmYE2gF8/+VMtmNJ2xzeGq2mNvvwI4wU2PAl5ynfeyva/psc0tS4BHGUfls81P64HXsKKsbwX26WjgW3cZiuWzj+8HDqtrhuPxNazJ3gX792GMn7KBFVT4QRFZKyJX2/vGQ92cD+wEvmGbRu8SkS5SLpsROuMEtboeY9r/XUQmYsXI+7i64u3B2C+fquZVdTGWVnAmVlzAMY+IXAS8pqprG52XFHmTqp6OZV76iIic6z44hutmC3A68HVVXQL0MWpKA9IpmxE6Y5tXRWQ2gP3/a/b+7cDRrvPm2PuaFhFpxRI4d6vqf9q7x035HFR1H1bA2rOBqWIttQGlZSiWzz4+Bdhd35xG5o3Au8QKR/V9LBPbbYyPsgGgVpR6VPU14F6sTsN4qJsvAy+rFbEFrKgtp5Ny2YzQGdv8GPiAvf0BrLEQZ//7bW+T12Mt5dDjl0AzICIC/DvwtKp+1XVovJRvpoyujjsBa7zqaSzh8xf2ad7yOeX+C+CXdo+z6VDV61V1jqp2Yy2M+EtVfS/joGwAItIlIpOcbayFI59kHNRNtQIkvyQiJ9i7zgc2k3bZGj2YZf4iD/p9D2stoBGsHsoHsWzh/w38CXgImG6fK8C/Yo0bbMIKiNrwMoSU7U1YKvxGYL39945xVL5TgSfs8j0JfM7efwzwGPAs1nIb7fb+Dvv3s/bxYxpdhojlPA/46Xgqm12ODfbfU1hLszCO6uZirODHG7GWkpmWdtlMGByDwWAw1A1jXjMYDAZD3TBCx2AwGAx1wwgdg8FgMNQNI3QMBoPBUDeM0DEYDAZD3TBCx2BIERE5aP/fLSLvSTjtz3p+/yHJ9A2GNDBCx2CoD91ALKHjmtEfRInQUdU3xMyTwVB3jNAxGOrDzcA59posn7ADgN4iIo/ba5NcAyAi54nIwyLyY6zZ4YjIfXawyaecgJMicjMwwU7vbnufo1WJnfaT9jowy11p/9q1fsrddjQIRORmsdYz2igiX6770zEcMlTqSRkMhmT4DPB3qnoRgC089qvq60SkHfi9iDxon3s6cLKqbrN/X6mqe+wQOo+LyD2q+hkR+ahaQUS9XIo10/w0YIZ9zW/tY0uARcAO4PfAG0XkaeAS4ERVVSdkj8GQBkbTMRgaw4VYcazWYy3jcBjW4lgAj7kEDsDfiMgG4BGsgIvHE86bgO+pFdn6VeA3wOtcab+sqgWscEPdWMsLDAL/LiKXAv01ls1gCMQIHYOhMQjwMbVWo1ysqvNV1dF0+ooniZwHXIC18NlpWDHcOmq475BrO4+10FoOK3Lyj4CLgJ/XkL7BEIoROgZDfTgATHL9/gXwIXtJB0RkgR3F2MsUrOWd+0XkRKzlvB1GnOs9PAwst8eNZgLnYgXX9MVex2iKqv4M+ASWWc5gSAUzpmMw1IeNQN42k30Ta82ZbmCdPZi/E/hzn+t+Dlxrj7tswTKxOdwJbBSRdWotJ+BwL9Z6PRuwond/SlVfsYWWH5OA+0WkA0sD+9uqSmgwRMBEmTYYDAZD3TDmNYPBYDDUDSN0DAaDwVA3jNAxGAwGQ90wQsdgMBgMdcMIHYPBYDDUDSN0DAaDwVA3jNAxGAwGQ934/wElDAsSMafn0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(log.data[\"Energy\"].iters[50:],log.data[\"Energy\"].Mean[50:],yerr=log.data[\"Energy\"].Sigma[50:],label=\"SymmModel\")\n",
    "\n",
    "plt.axhline(y=eig_vals[0], xmin=0, xmax=log.data[\"Energy\"].iters[-1], linewidth=2, color=\"k\", label=\"Exact\")\n",
    "plt.xlabel('Iterations')\n",
    "plt.ylabel('Energy')\n",
    "plt.legend(frameon=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68f501c7",
   "metadata": {},
   "source": [
    "## 8. Measuring Other Properties\n",
    "\n",
    "Once the model has been optimized, we can of course measure also other observables that are not the energy. For example, we could decide to measure the value of the nearest-neighbor $X-X$ correlator. \n",
    "Notice that since correlators do not enjoy the zero-variance principle as the Hamiltonian instead does, it is important to use a larger number of samples to have a sufficiently low error bar on their measurement. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "6cf7a569",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = sum([sigmax(hi,i)*sigmax(hi,j) for (i,j) in graph.edges()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "10fb783e",
   "metadata": {},
   "outputs": [],
   "source": [
    "vstate.n_samples=400000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "e3b3b4ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.851 ± 0.017 [σ²=31.585, R̂=1.0000]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vstate.expect(corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1daae42d",
   "metadata": {},
   "source": [
    "And we can further compare this to the exact ED result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "c89a3597",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.852248713127903\n"
     ]
    }
   ],
   "source": [
    "psi = eig_vecs[:, 0]\n",
    "exact_corr=psi@(corr.to_sparse()@psi)\n",
    "print(exact_corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a02b475-eff7-49f3-8406-8913d62194f7",
   "metadata": {},
   "source": [
    "# Extra: Assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92a0939c-2f00-4aec-aeb9-da69fd1403d3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Netket development)",
   "language": "python",
   "name": "dev-netket"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}