{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling and calculating observables\n",
    "## Generate new samples\n",
    "\n",
    "Firstly, to generate meaningful data, an RBM needs to be trained. Please refer to the tutorials 1 and 2 on training an RBM if how to train an RBM using qucumber is unclear. An RBM with a positive-real wavefunction describing a transverse-field Ising model (TFIM) with 10 sites has already been trained in the first tutorial, with the parameters of the machine saved here as *saved_params.pt*. The *autoload* function can be employed here to instantiate the corresponding *PositiveWaveFunction* object from the saved RBM parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from qucumber.nn_states import PositiveWaveFunction\n",
    "\n",
    "from qucumber.observables import ObservableBase\n",
    "\n",
    "import quantum_ising_chain\n",
    "from quantum_ising_chain import TFIMChainEnergy\n",
    "\n",
    "nn_state = PositiveWaveFunction.autoload(\"saved_params.pt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *PositiveWaveFunction* object has a property called *sample* that takes in the following arguments.\n",
    "\n",
    "1. **k**: the number of Gibbs steps to perform to generate the new samples\n",
    "2. **num_samples**: the number of new data points to be generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[1., 1., 1.,  ..., 1., 0., 0.],\n",
      "        [1., 1., 1.,  ..., 1., 0., 0.],\n",
      "        [0., 0., 0.,  ..., 0., 1., 1.],\n",
      "        ...,\n",
      "        [1., 1., 1.,  ..., 1., 1., 1.],\n",
      "        [0., 0., 1.,  ..., 0., 1., 1.],\n",
      "        [0., 0., 0.,  ..., 1., 1., 1.]], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "new_samples = nn_state.sample(k=100, num_samples=10000)\n",
    "print(new_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the newly generated samples, the user can now easliy calculate observables that do not require any details associated with the RBM. A great example of this is the magnetization. To calculate the magnetization, the newly-generated samples must be converted to $\\pm$ 1 from 1 and 0, respectively. The function below does the trick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_pm1(samples):\n",
    "    return samples.mul(2.).sub(1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the magnetization is calculated as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Magnetization = 0.54968\n"
     ]
    }
   ],
   "source": [
    "def Magnetization(samples):\n",
    "    return to_pm1(samples).mean(1).abs().mean()\n",
    "\n",
    "\n",
    "magnetization = Magnetization(new_samples).item()\n",
    "\n",
    "print(\"Magnetization = %.5f\" % magnetization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact value for the magnetization is 0.5610. \n",
    "\n",
    "The magnetization and the newly-generated samples can also be saved to a pickle file along with the RBM parameters in the *PositiveWaveFunction* object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn_state.save(\n",
    "    \"saved_params_and_new_data.pt\",\n",
    "    metadata={\"samples\": new_samples, \"magnetization\": magnetization},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *metadata* argument in the *save* function takes in a dictionary of data that you would like to save on top of the RBM parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate an observable using the *Observable* module\n",
    "### Custom observable\n",
    "Qucumber has a built-in module called *Observable* which makes it easy for the user to compute any arbitrary observable from the RBM. To see the the *Observable* module in action, an example observable called *PIQuIL*, which inherits properties from the *Observable* module, is shown below. \n",
    "\n",
    "The *PIQuIL* observable takes an $\\sigma^z$ measurement at a site and multiplies it by the measurement two sites from it. There is also a parameter, *P*, that determines the strength of each of these interactions. For example, for the dataset $(-1,1,1,-1), (1,1,1,1)$ and $(1,1,-1,1)$ with P = 2, the *PIQuIL* for each data point would be $\\left( 2(-1\\times1) + 2(1\\times-1) = -4 \\right), \\left( 2(1\\times1) + 2(1\\times1) = 4 \\right)$ and $\\left( 2(1\\times-1) + 2(1\\times1) = 0 \\right)$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PIQuIL(ObservableBase):\n",
    "    def __init__(self, P):\n",
    "        super(PIQuIL, self).__init__()\n",
    "        self.P = P\n",
    "\n",
    "    # Required : function that calculates the PIQuIL. Must be named \"apply\"\n",
    "    def apply(self, nn_state, samples):\n",
    "        to_pm1(samples)\n",
    "        interaction_ = 0\n",
    "        for i in range(samples.shape[-1]):\n",
    "            if (i + 3) > samples.shape[-1]:\n",
    "                continue\n",
    "            else:\n",
    "                interaction_ += self.P * samples[:, i] * samples[:, i + 2]\n",
    "\n",
    "        return interaction_\n",
    "\n",
    "\n",
    "P = 0.05\n",
    "piquil = PIQuIL(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *apply* function is contained in the *Observable* module, but is overwritten here. The *apply* function in *Observable* will compute the observable itself and must take in the RBM (*nn_state*) and a batch of samples as arguments. Thus, any new class inheriting from *Observables* that the user would like to define must contain a function called *apply* that calculates this new observable. \n",
    "\n",
    "Although the *PIQuIL* observable could technically be computed without the use of the *Observable* module since it does not ever use the RBM (*nn_state*), it is still nonetheless a constructive example.\n",
    "\n",
    "The real power in the *Observable* module is in the ability for the user to easily compute statistics of the observable from the generated sample. Since we have already generated new samples of data, the *PIQuIL* observable's mean, standard error and variance on the new data can be calculated with the *statistics_from_samples* function in the *Observable* module. The user must simply give the RBM and the samples as arguments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "piquil_stats1 = piquil.statistics_from_samples(nn_state, new_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *statistics_from_samples* function returns a dictionary containing the mean, standard error and the variance with the keys \"mean\", \"std_error\" and \"variance\", respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean PIQuIL: 0.1419 +/- 0.0014\n",
      "Variance: 0.0188\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    \"Mean PIQuIL: %.4f\" % piquil_stats1[\"mean\"], \"+/- %.4f\" % piquil_stats1[\"std_error\"]\n",
    ")\n",
    "print(\"Variance: %.4f\" % piquil_stats1[\"variance\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if the user did not have samples generated already, that is no problem. The *statistics* function in the *Observable* module will generate new samples internally and compute the mean, standard error and variance on those samples. Since the samples are not an argument in the *statistics* function, the user must now give the following additional arguments to the *statistics* function to generate the new samples. \n",
    "\n",
    "- **num_samples**: the number of samples to generate internally\n",
    "- **num_chains**: the number of Markov chains to run in parallel (default = 0)\n",
    "- **burn_in**: the number of Gibbs steps to perform before recording any samples (default = 1000)\n",
    "- **steps**: the number of Gibbs steps to perform between each sample (default = 1)\n",
    "\n",
    "The *statistics* function will also return a dictionary containing the mean, standard error and the variance with the keys \"mean\", \"std_error\" and \"variance\", respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean PIQuIL: 0.1409 +/- 0.0014\n",
      "Variance: 0.0186\n"
     ]
    }
   ],
   "source": [
    "num_samples = 10000\n",
    "burn_in = 100\n",
    "steps = 100\n",
    "\n",
    "piquil_stats2 = piquil.statistics(nn_state, num_samples, burn_in=burn_in, steps=steps)\n",
    "print(\n",
    "    \"Mean PIQuIL: %.4f\" % piquil_stats2[\"mean\"], \"+/- %.4f\" % piquil_stats2[\"std_error\"]\n",
    ")\n",
    "print(\"Variance: %.4f\" % piquil_stats2[\"variance\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TFIM Energy\n",
    "Some observables cannot be computed directly from samples, but instead depend on the RBM as previously mentioned. For example, the magnetization of the TFIM simply depends on the samples the user gives as input. Whereas the TFIM energy is much more complicated. An example for the computation of the energy is provided in the python file *quantum_ising_chain.py*, which takes advantage of qucumber's *Observable* module.\n",
    "\n",
    "*quantum_ising_chain.py* comprises of a class that computes the energy of a TFIM (*TFIMChainEnergy*) that inherits properties from the *Observable* module. To instantiate a *TFIMChainEnergy* object, the $\\frac{h}{J}$ value must be specified. The trained RBM parameters are from the first tutorial, where the example data was from the TFIM with 10 sites at its critical point ($\\frac{h}{J}=1$). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = 1\n",
    "\n",
    "tfim_energy = TFIMChainEnergy(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To go ahead and calculate the mean energy and its standard error from the previously generated samples from this tutorial (new_samples), the *statistics_from_samples* function in the *Observable* module is called upon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: -1.2362 +/- 0.0005\n",
      "Variance: 0.0022\n"
     ]
    }
   ],
   "source": [
    "energy_stats = tfim_energy.statistics_from_samples(nn_state, new_samples)\n",
    "print(\"Mean: %.4f\" % energy_stats[\"mean\"], \"+/- %.4f\" % energy_stats[\"std_error\"])\n",
    "print(\"Variance: %.4f\" % energy_stats[\"variance\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact value for the energy is -1.2381. \n",
    "\n",
    "To illustrate how quickly the energy converges as a function of the sampling step (i.e. the number of Gibbs steps to perform to generate a new batch of samples), *steps*, the *Convergence* function in *quantum_ising_chain.py* will do the trick. *Convergence* creates a batch of random samples initially, which is then used to generate a new batch of samples from the RBM. The TFIM energy will be calculated at every Gibbs step. Please note that the samples generated previously (new_samples) are not used here; different samples are generated.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'% Error in Energy')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXuYFNWZ/78vA8NlmGGYYQQREJCbaAR1RN1cNGqMmlU03pPsJtHEGC8xibkYk7iuG5Nf4m4uZt0kmpgYE29Zo5L1HqMxYlQQ8YKCwMBwExhgGGBAYJj398fbb87pmqqu6p6u7hl4P8/TT3dXV1edqjrnfM/7nnPeQ8wMwzAMw8hFn3InwDAMw+j5mFgYhmEYsZhYGIZhGLGYWBiGYRixmFgYhmEYsZhYGIZhGLGkKhZEdAoRLSKiJUR0TcQ+5xHRm0S0gIjuSjM9hmEYRmFQWvMsiKgCwNsAPgRgFYA5AC5k5je9fSYCuA/ACczcSkT7MfP6VBJkGIZhFEyalsUMAEuYuYmZdwG4B8DMwD6fBXALM7cCgAmFYRhGz6Rvisc+AMBK7/sqAEcH9pkEAEQ0G0AFgOuZ+bHggYjoEgCXAEBVVdWRU6ZMSSXBhmEYeysvv/zyBmZuKPT/aYpF0vNPBHA8gFEAniWi9zDzZn8nZr4VwK0A0NjYyHPnzi11Og3DMHo1RNTcnf+n6YZaDWC0931UZpvPKgCzmHk3My+D9HFMTDFNhmEYRgGkKRZzAEwkonFEVAngAgCzAvs8CLEqQETDIG6pphTTZBiGYRRAamLBzB0ArgDwOIC3ANzHzAuI6AYiOiOz2+MANhLRmwCeBvBVZt6YVpoMwzCMwkht6GxaWJ+FYRhG/hDRy8zcWOj/bQa3YRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRS6piQUSnENEiIlpCRNeE/P4pImohovmZ12fSTI9hGIZRGH3TOjARVQC4BcCHAKwCMIeIZjHzm4Fd72XmK9JKh2EYhtF90rQsZgBYwsxNzLwLwD0AZqZ4PsMwDCMl0hSLAwCs9L6vymwLcjYRvUZE/0tEo1NMj2EYhlEg5e7g/hOAscx8GIAnAdwRthMRXUJEc4lobktLS0kTaBiGYaQrFqsB+JbCqMy2f8DMG5l5Z+brLwEcGXYgZr6VmRuZubGhoSGVxBqGYRjRpCkWcwBMJKJxRFQJ4AIAs/wdiGh/7+sZAN5KMT2GYRhGgaQ2GoqZO4joCgCPA6gAcDszLyCiGwDMZeZZAL5ARGcA6ACwCcCn0kqPYRiGUTjEzOVOQ140Njby3Llzy50MwzCMXgURvczMjYX+v9wd3IZhGEYvwMTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIxcTCMAzDiMXEwjAMw4jFxMIwDMOIJVYsiOh0IjJRMQzD2IdJIgLnA1hMRD8goilpJ8gwDMPoecSKBTN/AsDhAJYC+A0R/Z2ILiGi6tRTZxiGYfQIErmXmHkLgP8FcA+A/QGcBWAeEV2Z639EdAoRLSKiJUR0TY79ziYiJqLGPNJuGIZhlIgkfRZnENEDAJ4B0A/ADGY+FcA0AFfn+F8FgFsAnApgKoALiWhqyH7VAK4C8GIhF2AYhmGkTxLL4mwAP2Lm9zDzTcy8HgCYeTuAi3P8bwaAJczcxMy7IFbJzJD9/gPA9wG8m1/SDcMwjFKRpM/ik8z8bMRvT+X46wEAVnrfV2W2/QMiOgLAaGZ+OFcaMn0kc4lobktLS1ySDcMwjCKTxA21lYi2BF4riegBIhpf6Ikzw3F/iByuLIWZb2XmRmZubGhoKPSUhmEYRoH0TbDPjyFWwV0ACMAFAA4CMA/A7QCOj/jfagCjve+jMtuUagCHAniGiABgBIBZRHQGM8/NmaLHHwdaWoCDDgKOPTbBJRiGYRjdgZg59w5ErzLztMC2+cw8Pew3b5++AN4GcCJEJOYA+BgzL4jY/xkAX4kTisYjjuC5r70G7NkjGxYtAiZNynkNhmEY+zpE9DIzFzziNEkH93YiOo+I+mRe58F1RkcqDTN3ALgCwOMA3gJwHzMvIKIbiOiMQhOMXbtEKC68UL43Nxd8KMMwDCMZSdxQHwfwEwD/AxGHFwB8gogGQsQgEmZ+BMAjgW3XRex7fIK0iFgAwEknAXffDaxbl+hvhmEYRuHkFIvMXImZzHx6xC7PFT9JMezcKe9HHy3vJhaGYRipk9MNxcx7AFxYorQkY9cuoF8/YMoUoLLSxMIwDKMEJHFDzSai/wZwL4B23cjM81JLVS527gQOPBCoqAD22w9Yv74syTAMw9iXSCIW0zPvN3jbGMAJxU9OAnbtAsaOlc/Dh5tlYRiGUQJixYKZP1iKhCRm585ssXjnnbImxzAMY18gyQzu4UT0KyJ6NPN9KhHligmVLh0dZlkYhmGUmCTzLH4DmSsxMvP9bQBfTCtBifDFYv16oLOzrMkxDMPY20kiFsOY+T4AncA/JtvtSTVVcfhi0dEBtLaWNTl5s2QJ8L3vlTsVhmEYiUkiFu1EVI/MbG0iOgZAW6qpisMXC6D3uaLuugu49lqgrby30TAMIylJxOLLAGYBOIiIZgP4LYCcK+SlChGw//7yubeKxYYN8r5jR3nT0dt56CFgxAjgXVsKxTDSJsloqHlEdByAyZCos4uYeXfqKYuishLok9E4FYtSzbV46ilg5Ejg4IO7d5yNG+V9+/bup2lfZt48aShs3iyiYRhGaiSZZwHIqndjM/sfQURg5t+mlqpc+JVCqS2Lz3wGOO444De/6d5xzLIoDtpIsPtoGKkTKxZEdCdk/Yr5cB3bDHFHlZ5hw9znujqZyV0qsdi8Gdi2rfvHUbEwy6J7qFiYG8owUieJZdEIYCrHLXxRDvr0kZAfpRALZhGKYrRizbIoDiYWhlEyknRwvwFZxa5nUqqJebt2yTDdYopFmGWxZ49b2MnIjYmFUUouvRS4uHzzkctNEstiGIA3ieglADt1IzMXvoBRMSmVWGzdKu/drZi2b3ciESYWH/+4RNW9887unWdfwMTCKCWvvLJPBy5NIhbXp52IbjFiBPDmm+mfR/squmtZ6EioqGMtWSIjvozc7NolfUiAiYVRGrZtA1askPh0/fuXOzUlJ1IsiGgKMy9k5r8SUX9m3un9dkxpkpeAkSMlmGBnpxtSmwZpiEWYZbF9u7i7jNy0tLjPJhZGKdi2TeqZ5cuByZPLnZqSk6t2vcv7/PfAb/+TQloKY+RIqVy1H6BQVq0Cjj1W3sMolhvKT2eY8PhuKiMa3x1gYmGUAm0wLl5c3nSUiVxiQRGfw76Xj5GZ+IZr1iTb/3OfAx58sOv2P/0JeOEFYP788P8Vy7LwxSLKsrBRUvGYWCTnyiuBb3yj3Kno/WgdsGSJ2/bUU8Ajj5QnPSUml1hwxOew7+UjH7HYswf45S+Bxx/v+tvs2fKuFkSQNMQi7Fjt7b3DsvjiF4E//KF85zexSM6jjwJ/+1txj8kMXHUV8NJLxT1uT2XXLnkB2WJxww3At77VvWM/+CCwcmX3jlECcnVwjyKimyFWhH5G5vsBqacsKfmIhYYzD6ukVSy2bAn/r4pIscSisrKrKDDLtnJMadm5UwpDdXWy/X/1Kyk0556bbrqi8MWiXJbY66+L//r008tz/iQwi2t14MDiHnf7duDmm4GaGmDGjOIeuyfS3u4++2KxYUN0AzMJHR3A2WcDX/taj49EnUssvup9nhv4Lfi9fGj4D18sli6VyuTYY7P31VX1gpX0mjVS6IF4y6KjQ159k0ZKCbBhAzB0qHTGBys5bSHv2JF+h32Qb35TTOpXXonft6ND7sdrr6WfrihaWmT2/p495bMsvv1t4LHHJF8NHVqeNMSxcaM0BKIaQYWi5cGvRHsazzwjUR4OO6z7x/IjN/h9Fhs3du/ebt4sZd0f+NJDiazxmPmOUiakYCorgYYGYPVqt+3ii4G333YCsnQpcNBBTiyClbRaFUC8WABSOQ0eXFh6N2yQkCU7dnQVLf/7u+8CgwYVdo5CaG6Wyj/JsEAtHCtXyloi5ago16+XhsKaNemKBTOwbBkwfnzX3+bPl/t1333SF9YT0QEb+6JYfPazwKGHAg880P1j6fVOnCj1ya5dMh9q40ZpPLW3A1VV+R9306bs9x5MCZuuKTJypBOGFSuAv/5VhGHLFuDpp4EJE4BXX3X7hInFwIHyihILf3t33B4bN4pYDBrU9Ti+WCQ5R2cn8IUvRHfK58P27XK8pqboc/3pT1J5+utwvPFG989dCOvXy4TMAQPSFYs//1nyz9Kl2dtbW0VgAeCOHtyu8sWimO5NLQ89WSzWrnUNxO6iYjF9upSF5ma5BzrMvdDRmCYWJcYXi7u8Eb9LlgBzMx6zN96IdkPNni1+19ra6BZY0LIoFLUsBg7MbVkk6eRuagJ++lPg4YcLS8s11zirSgt91LDAZ58FzjgDeO45NxkOiHZFzZsH/PjHhaUrCevXS1ywtMVi2TKpZIMdkHrdJ5wA/P3vYsn2RFQsOjuLO3BCy0NPHYyxfbuksVjRHXyxAKRu8V1HJha9BBULZgmToaHLlywB3npLPjc1RbuhFi8Wv2Z1dTI3VHcsiw0bgPr6cMvCb6UlKYRaYRXSutu2Dfj+952JrueLqvS00K1bl21ZRInFL34BfPnLbgRJsSmVWGiF4Ask4Ky5//xP6Vu699700tAd/HlDxXRF9XTLQidtrltXHIvKxCJeLIiogYiuJaJbieh2fSU5OBGdQkSLiGgJEV0T8vulRPQ6Ec0noueIaGohF4EDDpBMMXeuhP742tdk++LFTiyWLg23LDo7pRDV1cnIjjTdUMzJLYsk5+iOWKgLRf8bZ1noOucbN7qKs7ZWRgSFsXy5XG/S+S/5wFx6sQiu8/7qq3L+ww+X/Bd0U/UUfIuomMv49vQ+C3+tk2IsK6DHGDdO+ipWr84WiELFQvPV3iAWAB4CMATAnwE87L1yQkQVAG4BcCqAqQAuDBGDu5j5Pcw8HcAPAPwwj7Q7Ro6UCuS//kse5Kc+JdsWL3Zxo6IsC/Xl1tamb1ls3y4VW5I+i7QtCxWLoDshSiz8FpBWOu99r4hFZ2fX/XV0WRrjx9vb5d51RyyYk41AiRKL+fNdK3PECPGP90S6a1nceCNw9dVdt/cWsQCK44rS662pkbpl9eriWhbt7elZ4UUiiVgMYuavM/N9zHy/vhL8bwaAJczcxMy7ANwDYKa/AzP7ubcKhU7207kWf/gD8JGPiJUwcaL42bdskSGWvlj4FbHfSq6uzt1nQZmJ64W2ZDVz1df3PMsiTizCLIv3v1/uiwrDq6+6OF16/DTE4pln5H3SpMLF4rHHJN/EVSRhYrF7N7BgATBtmnzv6WKhw8sLEYuHHw7vE+vpbqigWGzf3jXg6NKlwA8Ttk9VLAYPFktyzZpssfBjleWDb1EEGyQ9jCRi8X9EdFoBxz4AgF9TrELIZD4iupyIlkIsiy+EHYiILiGiuUQ0tyXsoahYdHYCn/iEfJ44UTonAeB975OWQNhoKH1AtbXxbqj6+q7/zwc9dnV19y2Lbduc66MYlkV7u4jhqlXh5/bNZRWLYzLxJLWf44wzxAW4dq0MKQWiY211h//+b3nmp51WuFg0NUlLLi592mL0+ywWLpT/qmWx//7FG3VTTHRC3tSMQV+IWKxbF26B9XTLwq8n1q0DfvQjoLFRhF657TaxmpK45/R6q6qyLQsiqRe6a1kEP/dAkojFVRDB2EFEW4hoKxEVraeMmW9h5oMAfB1A6Lx5Zr6VmRuZubGhoaHrDioWtbViWQAiFso//7O879kjgrBjh+v0CloWudxQeu5CxUILVlVVuGWRTwe3P2Q1rMC+/rpU5lEVhFoD7e1u5rjeM3+GqqIZeeNGKVyDBwOjRsm2detEqFevBl5+2R0byLYsOju7dhTny9tvS7iWz31OXI4DBxb2PPQ5x6UnzLJYuFDeDzlE3keMkMopSbTgFSuAMWNkEmTaM883b5bn2h2xWLtWnn3Q1VhMy+KXvwQuvLD7x/EJWhavvy7323+O2shJMgN72zZpmPTtm21Z1NbKgBoTC4CZq5m5DzMPZOaazPeaBMdeDWC0931UZlsU9wA4M8Fxu7LfflJpnHeePFDAVXw1NeJbV8aPl8pRW775uKFULPyWbHu7O1YcKgBVVd2fZ6EuqEmTwgvsk08CL74YXvED2W6od9+Ve6It5TBXVNCyqK2V+w5Iwdy0ScR40SI3qGDAACcWu3cDZ50lz6WQ0SkPPijuxQ9/WETikkvcOQqxLLSlWIhY6DWNGSPvI0bINSVxRcybJ///7neBD30ofJ+33iqOW0utpkLFYts2N/8m2PrOZVl0dopLMilPPQXcf39431ehrF8vjRkiEQvN0/5zXLRI3pOKhU7EPeAA+c/y5WJVDBuWLRZPPgn88Y/J0rlpk5vUmoZYtLZGD0LJk0ixIKIpmfcjwl4Jjj0HwEQiGkdElQAuADArcA6v+Y+PACgs9m9FBfD888BNN7ltKhYHHyyztxX9rJWxLxY1NdkBw3yiLIuTT5YhoknQgjVokLMs/IozHzfUa6+JuB1ySHiB1cl1USa274bSc6kPPmz4rN9n0dYGDBkihWfAACmY2pLr7HQ+7qOOkoqRGbjoImDWLClUhbSo58yR89bUSAA79cMXKhZJLIvOTnfdQbGoqpI8A7i0JKng1RX6yU/KHJewhsZHPyqharrbMdtdsfDPH3RF+eFvfNcOIJXl9OnSr5OEjRvlGN1dZsBn/XrxONTXy3VontbnuGePa0jlKxbqyXj99XCxuOkm4N/+LVk6N21ydVIaYnHjjcBxxxXlULksC60B/yvk9Z9xB2bmDgBXAHgcwFsA7mPmBUR0AxHpkqxXENECIpqfOd8nC7sMSOas8QwefQAHHyyVvE7F1+1aQQYtC6Br5mGOFovFi5PFU/LPqZZFZ2d2QcvHsnj1VeA975EMHDY0UMUirDLcudP52NvbndgMHy6Zf8UKd4xnn5XPYZYFkVgXvlgA4iZqaACmTJEK629/A373O1dpFTKEc/16Oearr2Y3CqLE4t57gVNOkUr378HlWJBMLDRuT3C/lSuB0aPdgId8xaKiQgYH6PcgGzZIq3XmzO4NC1axGDdOGif53ndfLIIVuV9Ggo0VDb2jLfc4VIhW53I85IkOrR4+XPKMlhHNxytWuEZhIZaFHiNMLLZsSX6vN22SCAH6udi89ZZccxFcnpFiwcyXZN4/GPI6IcnBmfkRZp7EzAcx842Zbdcx86zM56uY+RBmnp45bsKmSAIGDpQJU5ddJoVaRUJj/PiWBZEITZRYaGC/oFgwywPWVnocwT4LoKs1oW60XJZFZ6cM3TziCDlWLssirDJUN8rQodmWRVVVdmftDTcA55wjn/2hs5s3i2UBuDXQ/Ypl+3apoEaPlkL76KMyce3KK+X3QsSipcW5vXyixOJ//kfWJ3nhBXFzBNFnnGsEilZiAwd2tSxGex7W/feX9ySd3GvWiLioCyusg33rVlmJ7cUXgSeeiD+mj2+pqniNGCH5Ow3LAuia/zTPNTdLeo46SjqTo0hLLBoaJH9qFAfAPUffei5ULIBssdCGxZYtyfrmmCU9Y8dK+UhDLLQeKEKgwr1jBncUV18tGRUQkaipkQcLZFsWNTXysNQyCRYqLRjBPostW8Scfeed6H6L5mbpwAOy3VAaJNBX/PZ2SUO/fpK+ZcuAr3+9qy93yRJJU5RY7NnjRoKFVczaAa0uLD9dvlisXCmVdFubnK+6Wiyh1audCyZoWaj7b+xYV6H+/vcyee3AA+V7IZ3cLS3u/vtEicX27cA//ZMU8LBCmMSy0AJ20EG5xUIjBgQti9WrJRyILyJr1ogbQwcHBMVi927JS9rPlm8FcuaZwKWXyucNGyQ/VVYWJhb+9QQrm1yWhea55mbJF3PnSlyxKLRV7ovF7t2Fhf5euzZ70ubw4dnWe5hYJLkvYW4owInFnj3uurdsyY4bFcXWrfK/YcOk4VZssejsdPVAEY69d4uFz+c+B1x7rWvR+5aFdjBFWRb6fcgQqcj1v/oAwmIHKd/4hkS/9JdLzWVZ6G87dkgn2Q9+0LVCmTdP3g8/XPbXTkhlzRpnYodVhmoJHXJI9uS0oGWhhVc7yNRcXrcuXCz69JHKEcgWi5UrgQ9+0Fkj3XFDBcklFoMGSWEOa1XlIxYTJkgFvmOH3Ne1a7PFYuBAubagWMyaJYEsn3vObYsTC02X/p5PBc8sbkPtXNZoAUA6loWG0I+yLJYvd5Xyyy+Hn2PnTvd/Xyy+9S3gyCPzS++yZXJvf/tbEQgVC8AtKeCLhaY/X8uiqsrl5fp6ly9V9PQ+R93vtjaJm6YDIurq5BVl5QZDzCcZdQdIXtNGbNqWBQmjc+3TazjlFGmlh4mFVny+WMye7Vrg/oScAQO6igUQ7opqa3Oxl7ZsCbcsgmKhv23f7h5wMDPPmyetxalTXV+M7xrzI8dGiUWfPm7RebUKBg1yE8yYXeHVkVf+cGQtLCoWa9dKoTn8cNk+dqyr8AARC73PhVoW+bihuiMWzNLi88VC99UYZKMDxSJsYt7zz8u7/zxULKqrpQKPEgt1deRTwW/YIGnU51kMsaivl7wSJhZ67FyWhYrFmjXhfTr+cX2x+NOfxILesyd5etXtpQEsfbGYMEHyg1bIixZJnxqQv1gA7vkMG+bug7qitL6IahTdcQfwpS+5EVN1dbkti29+063NM2+eXEfUKEcfP9+lLRbMzAD2rgVmg5W0Lxa+G+rMM4F//3f5rg+/ujp7XL//AMLE4r77XEW2dasUqgEDpPAFRUvTpCOlduxwmSdYyOfNk87tykonFu3tstTpaae5TFJREZ5hm5sls+t1+2Kx//7SKmtudtetYuGPKtP/qpn/9ttSOP/pn6QP6LDDXIWqHbqFWha7dsl/oiwLXZDKR+9lXV3+YnHllcCJJ2a7oQCpaNSCLEQsdu6UY6obY9SoaLEYOjT/TmntUA4TiyFDChOLESPC7+HWra5jP1efhT8MWy1iH/+4ei/WrpWOWe0TTIoKgQZ59MVi4kS5p75lMW2a5E295w89FN1XGBQLfYbqhgLkfuu8JSC6UaSDRh56SN7Vsoi6Vg1btGGDWKq7d2fPZYqilGKRYR4RHdXtM/UUklgWzc3yYLRi0Mw0eLD8XwXAf7hhD+83v3Gft251biYgP8vCL+TMMvrqiMzoZV8sXnlFwlg89pgUgokToy2LAw90/9XKRd1QQHanoLo1tIUNZFsWgEwSHD5cRGztWvG5a2Xd2Cj3tlDLQs37KMsC6NpnFGZZrFrlfLi5Orife04K9LJlIuxjx7p9o8QiOIt77VpXWPVdf88lFr4Vm681oK34bdskf2/Y4AS2piZeeB56KDvE/7p1boRcmGWhFXGUWLS2Sj5SCzPMFaXHHTLEWRYazgXIL4xG8FlGicWOHTKSafJkNxG3uVkaiNq/GCTKsvDFoqUl+3lt3iyWkR+Zmdmth66NiTix0Pz/yitOcJOMbmpqciP2SiQWRwP4OxEtJaLXMlFiy7ieZjfJZVmoWOgD0WGNfgH2LQt9uAMHdrUsVq6UzHDyyfJdLQutoHNZFioWenyt2J5+Whbj2bTJuXt8sdD9//AHGW3T0BBeMS9fLmKhmT9oWQDZBTvYZwFk91novdDPfqV+7bXy0uNHWTu50PSFWRZh91FnpKtY6H35/OeBf/kX+RxlWXR2SqXLLPNF6uulIAP5WRY6XHfsWCcWmp+SWBbqpspHLPyhqi0t+buhvvc9WWVSr2Pt2nCx6OiQ+62WRbA13tbmKqnnnpOGzaRJuS2LadOcWPzlL9nXkRR9zlrGo8RC1yiZMMFNxNX+mTlz5P3ii8VVBEiFv317vFhs2NBVLB58UIb163D0xYslP9fXOwtk6NDcYqH3aN48dw+TBBpdulTK+aBBJROLDwM4CMAJAE4H8M+Z995JsHJpbU0uFtXV4X0Whx3WVSy0lXfSSfKuYqEZOcyy0N9VkHzLorVVXCMqPmGWhbasmMV1MmRI14q5o0MqqKSWxZAhruXou6GClkXws3L11RIzCpAKpLY2XCzOO0/Ck7z4YtfftMKIckMB2f0Wu3ZJpT9woBTK1lY3MmTVquyJl0GxWLXKPd+mJvm/DoBQsait7bqs7ogRkk80rzz/vLgJzzlH8kZHR7hYvPNO9mid7oiFP8Jn+XLJW0GxyDV7fsUKuY86j0XdUEGx8Ofl+N+VzZtd/9bOnfL5yCOlAbJpU/acHD3uYYe58CRPP+2suXwti4oKmdTYp49c+/TpMkHuvPOcWGjFPXassyw0HS+/LHnj7ruBl16SbVpG/Wc+frycY8QIKTeVlXJt/vNqa5M85IfqV6vii190+6lYqCUSRC2L555zDYKklsX48dH9dnmSJNxHM4BaiECcDqA2s6134lfSHR2SUVQs+vaVCkYLXVubFIQoN9TGjbJtwoSuYqHfDz1U3qPcUEHLQifs+ZbFli2SYZjFVP7858PFYtMmV9mPHy/XFawM16yRDDl2bLhloa1FtSxmzHD/3W8/9x+/z8L/PY4hQ8KtnaeeEqE45hiZG+GjFUYuN5QvFlq41Q2mMal0DQJ9ng0Nsq8/Y18Lo46e8cVi82YRk6BVAbj7pi3U558X99vBB8v9XrkyXCyYsy2SXGKxaZNehC/tAAAgAElEQVRUfmEtdE27pkNDrvhisWdPdCWjEzX79QN+9jMR1vb2cMtCBTGXWGg0AECsiiOOkHtwwAFuxBzgjvue98i7hqg57zz5HiYWixeHb9eG33e/K/2FlZUiHtdf7zqRfbEYM8YFD9UKeeFCcYPt2OGehe9ZUP7lXySfDh8ujSDt1/E7y/3BBlqWn31W8t1nPiPfdTnnurquSxYD8sy0Efjoo07s/QEtt90mnoaglapikctqyYMkix9dBeD3APbLvH5HRFd2+8zlwrcstCBqxQdkzwIHpADlckPV10srfeXK7E5WHXHkj7gIc0Ml7bPQCvaii2TCWb9+8l2Pt3GjpOuTn5RKVZeJDVbMKmJBy6JfP3lVV8v21lap2DX91dXyu7pk1LLQygjIFo4owiyLLVvkXl52mXwPrimeyw0VJxYaKXjVKrkX7e3ueFrp++nRIIFqDdXXu/yhloU/yksJzrV4/XWpIHUSaFOTiEW/fi5N/vBZ7XPJJRbz5kn/0a9+1fX8HR1Syb7vffI9TCwAdzzmbCtDXUBXXy1pufxyd11BsdA0qnj7YqGLiU2aBPTvL9smTQJOPVUq55EjsxtWGzbIc1Kr9YYb5F3dhWGicMop4eE0WltFEEaPBs4+u+vvvlhUVEjDKmhZMAM/+Yl81nIfJhaVlW4OF+DuUdANpY0Hraz/9jd5RiNGyH3R8hQVH0qjCOhcDkXz+KWXSpy0+fOz43Ft2yb5/KCDSmdZALgYwNGZmdfXATgGwGe7feZy0b+/tAR27MgO9aGoK0oz+po1kpn69pVtQTdUXZ1UvHv2ZIduaG6WgqGZISgWUZaFuqFaW12G2LrVtS40Uyl6PPWljx0rBf/ii50byq8UtCPeF4uWFpcewFknI0e6ClXPqxWdb43ptnwsi7Y2GU67aJGrPI4+Wt6DGbulRQq3/5yUMLHQe+qLhR+lVzu59dqCweVqasSVodfbt69UFFrRhFkWep7WVjehrKGhq1iMHOn8+SoWv/61nHPOnGwrNigW2vfx4INdJ2o2N8t5VSx07QZ/NBTgXFHnnitWqh5H889JJ4loP/qofB8+XI6xY4fLj1p51tWJ+PlisW2bHHPoUDcJc+JEmdPT3Ax8+tOyj1pzGzfKvdM+gGeekXQdeqg876BYdHRIHg6LI9Xa6spbGBqxoKlJzte3b1exAIBHMgNAc1kWQZKIhU60bWyUbZ/+tAgf4CzC4Cx2TZe6s7Ucah6/+253PL/Ro3mllG4oAATAd6TtyWzrnRC5IH65xELXalizRiqI/fd3//WHztbVOf+q32LSEUeawbZscWIAdLUs/E7ZQYOyXRPaZwFEi4Wa1kOHOhdKbW32uG8/jWPGuLT56QKcWBxwQFexCFoWQHjHdhRqWbzyilQMjzziBGzyZLmeYOtKZ2/3CcmuSS0LP/Kmnk+vbdMmCZH95JMiFpMnuxnU+v+hQ2WU2YYN4ZPF/Jah5qu6OrmH/fpli4WiYnHbbVJ5Llokz6qiQq4rONxVRW7NmuzRaoBznx15pPw3KBZqWbS1SQiR+++XSYO3Z1ZI1vwzerS07vV/alkAXef9DB7cNYKAVli1ta5z1b9mf7CAHtMXC0DmFQDyzINisXat5Okwd5ofwTUM/e2111y4Fe3g1sEAvtWYy7IIEhSLQYPkXqhYtLa6a1Er9Jpr3OgrnfMUjKeloqgRimfMkPyhyyy0tzsXnp9XNI+PHVtSsfg1gBeJ6Hoiuh7ACwBC7OBehFb4YWKhheoDH5D3NWukVap9D8Ghs35GD1oWBx4oD7aqqqtl4Vs4gJj+nZ3OsvBbjr4bKti61uOpv9JvWYUNVW1udiHd9b/+cYBwy0KPW1/v1pFQ8hELtSy0JfvWW07Axo4N969Gzd4GkvVZANlioZWuVhivvALcc49ED164UFxvBx4ohfncc2WfoUPlt5oa4GMf65oOvxPcF/aKCjnW0qVyzX7FWVsraVRLQ/tTqqtdvDLfMmxqkoqmb1832VPRSmbyZLlX2kINisXmzTI5ddw4yeNf+Yq4Wn2xGDpUJrYNG+ZapoCrcPwBH0Gx8PPpeedJ9ALy2pb6PPz1Uerr5VjDhkkIem0pNzR0tSA0n0ct0JVELBYudM9e+yw0HdoQOPhgyVMdHcnEQvOtCuno0V0ti1wDNUaPlryg7kPF79N573vFPapuah3IoXnKtyx8d7M/IrAbJOng/iGATwPYlHl9mpl/3O0zlxO92Zqx/QymlkVjo1SIK1bIA1SxCHND6X80U2mHpprhaur6Hdy+hQNkV3B+Kx9IZllo5euLRdgkOBUxQPyu2vfhn1NN4jDLQsN4+BWAtpSSikVbm6uc3nxTWkEDBsj/w8QiKi4UUJhlEXRD6bj+N96QymjyZLm+733PdfDr9X/60+GVhgrzpk1dn9X48dKSf/tt52oD5Bwf+5gEvKyokOtUsQC6dko3NUnH8fHHiyvKZ8ECuUcNDe45ELk0qFj8+7+Lb/s73xGLpq1NrIsVK+S/2gj4+MdFpGtro8UizLLQMjVkiHTi/jhQVUSJBSAWz+9+5/YdNqyrZaEi2B2x2LMn27LQDu76evEoELlF1LZuzc+yaGuThuB++2ULRJxYaEQF7TPT+G96z4cNk9FQl13mGrt63xsaJP/4lkVzsytT9fVFWSskLtxHBREtZOZ5zHxz5pUwHncPJpdloQV14kRR7GefFQX3LQs1AVUsNBNppnrnHWmRBMXCHzoLZC+AFAxfrgwdKv/dvFkqd79FD0jG7NMnXCyiLAt1m+n5/Hcg2w01cqQcXwvadddlxzoCpOANG9ZV5MKorXULxwAiFsuWSZp0VIm/iP2770aH+gDixaK2VtKv4gQ4f66KxV//KvtoC01dAsF0A64TPkhFhVSQra0u/XrPpk+X33/8Y2nJ+9x2m1g0ujRnUCwAVwk0NYlFcOaZUqloxQKIGKo7Qu9Vfb2c1z/W7NkySOKCC6SD9ZBDZFswOCLgGgRxbii/4vbdUGGEiYVaP4cfnj1gIswNpZZF0A3FnB3nLQz/N18s9uwRERo2TGbvP/+8G9jhD4eOE4tdu6Ts19TI9Tc1uU7pOLEA5JwLF8q1HX64NFbUstJnAHQVi6qqrv1bzc1yjbrsaxGIC/exB8AiIhpTlLP1FIKWRdANRSStwZEj3VoVQTeURpwNEwu/XwBwflHfstBjhVkWviCMHessC11DwodIzq99HLnEorMz27IAuna4A9luqL59pYKcOdPtr78r117rQhjEMWSIFGxdGKe1Vcaza5p8sTj9dAkTonGnwogTC1/oVJCClsW6dSIQX/2qfNdK1+fjHxehnDQp+to07WpZ6LO4/nqpjK+6quvzUzTMdZRYtLXJscePd89CrYvOTrmfQbHwK94xY6Slf/fdMppK+3/e9z6pHJctc/k1SBI31Nq10hr2LYuoewS4pVpbW6MrM3VDhY3a0me8dKmsXaIRXOM6uBV99nqvddW7qiqxLrRMb92aLY5R6DUsX+7Ewh/KmlQsmpul72zbNrEAN250oxQVrb/8WHPBOVV+OS+FWGQYCmABET1FRLP0VZSzlwvfsujTJzsTfPSjwNe+Ji12bWkSiQ9T//vuu67g1NdLi7+y0mUq318IyIP2J74pYZZF0A01bpwTi6hWU1WVFKiKiuxMFXRDrV8vfSO+WPiRNJVp0+R6tPL56U+l4o5i6FB3f+JQAXvjDVcZrlrlrB1fLBYulI5cHVkUhopFcFQZ4O6jViCjR7vJT4C4z9QNN3068IUvSOMgzLI491wXKywKHZoZdEMNHBg/rFgrxq1b3TPxxUIFbvx46YQ96ijXb6FzIvR56b3yxaJvX7FiLrgg+7zvfa/kD9+PH6S+XsqANkg0n1dVObG4/HJZ6z6qb03xxUKHheYSi927syvBYJ/Fz38ugxM0lEq+loXe4127stPhBxXV5xklgEBXsQgOANEO7n79ug7PV6ZMkXL8s5/J97ffdu4xv5GR1LLQMlUkseibYJ9vF+VMPQlt0be2urUslJNPdrOkVSwmTHCtfa2ctDNbM7+/Yl2YWCxdKp99IfAti7CItP37SyWjHdxRBVAr+qFDszNV0LIIpsv/r5+u6dMlPX2TZI880UK0dasIs0be9DP2xo3ixlu3zlXAUW4ofS5RloUec/FieZ4tLa611r+/G545bZrkA12HvBA0zHTQDZWEYcOkb6xPHycsvljoM9ShuGedJRbd6tWuPyaXZRGFDrUFwocEA9JwGDPG5eFt21zolqoqqcDnzBHrSVv+URWrlrdNm8JdLD4qei0tLi8HLYutW6WC1TAdhbihFP9++d6C1lbZL1d50Hpg9Wo3IVaZMkWerfa9RVmX2uB67DF5b2oSSzb4HKPEQkV1xw5pGPrWehGI7bMAcD0z/zX4KsrZy4W26KPGzCsqFuqCAlzlpJk2SizUpAXkQWqrzG/BjxjhfPdhbqj6ejdaI86y8NOiaIHVFtwtt8h3v+UcZlkA6QgFkF2IjjnGFVbfsti9W1rLnZ3i+vne98RPH0acGwpwldHIka7g6Xn9foXuomGmW1vlvqrVkgTtzNVFpoBssfDHzQPufjz0kBMLXbY2H7EYO9bl8yjLApAGk4bF9oPqDRokZUH7zGbPlvxbWRl+HHULbtqUbZ2H4YuFEuyz0DKnYWJyiUX//pK26mpXNnyxyGVZxAm//pfZuaGUgw92YU6iLGRA+klVSMaPF7fanDld74/WX75Y+MOstX+ulG6oTJ9FJxHlsL96IarMS5e6whdGErHQBxEUC7/1Xl3tZnf7lfKMGeJr9pc29S2LujrJeLt2uVZ2GFFioZMI29rEH3/nneJK8f3uYZZFmvgtzjFjXAXniwXgJtGNHStDWKPcODp5MkwsfNEFpMM+KBZaqP3wFIWilkWSyiXIsGFuNI2mzZ9I19QkadX0HnywvH7xC/Ftjx/vKvB8xILIzSlJKhZ+v0pVVfYAipdeiraAFXU1agMqymrUivXttyWW0saNrty9+640JrTCTCIW+rt2/ALFFwsgWywqKmQWta5emUssBgwQtzPgVjtcvz7csvD7LIKWRdCDENbXWQBJ+iy2AXidiH5FRDfrq9tnLieDBsmN1tgpUegEHV8stCUbtCx0xBMgrR+/4AU7p5QZMyTDz5sXvoqeWhZ6vjg3VJi5WVsrHZg//KFkwG9/O/y/pRIL/xp8sQiazCoW/ryEMPr2lVdQLCornXWkxwyzLGprpbLS4cLdQVvMcZPDwmhocB2+QcuirU0qmmBevf56mWD2wAPZnfJhfRa5OPFEuVe5ysKECVJZt7aK+0iFzG/8VFZKn1gu3z7gxEKtaq0gg+h1XHuthOC46SY5vuYJv3WtIWLiXC7Dh2dfZ5RYBN1Qcc/TP68vFjp0FRBXaNzw8ilTpPz/67+GpwsId0P5lkVQLCoq8s+PISTxNfwx89p7GDjQBdTzI6kGee97ZYal7wLRilzNYX0IvmURHOHhZ8igZQFIq0gzW5hlAbgQCmHEicXs2dICv+GGri2MKDdUWgQti9NPl4pQLQe9Bh0tFRx5FUZwtbzgjPRcbqjPfS6/yKa5GDpULMiVKwuzLBR9JprGLVukVR90lZ17rnTwPv10tliMHSv5NOmgg898RkKv5KrINDz9kiUy6MAfHQfIczrwQAmul8SyaGmR515dHX2vVCy08/rWW+V94kQpv9u3uzKnkXvj7vudd2bndb+z2X8GQcsi1yg4QIRS6wC/g3v4cJend+3KbVkAspzsJz4h/9Ohw0HRVzeU38D0O7ibm0Ug/IbWlClu/YwCiRQLIqph5i3MfEfIb717KO2gQW78c67WVEWFxFjyUbF4+GEZo66+2cGDXaUT7IyOEouGBmlVvfSSc4P48yx8ywKIF4uw3zXTXnBBeEYtlxuqb19pzZ91lryUoGWRpMUfJxZ63WFuKI0BVQw07UuXSms9H8Iqqv795fXOO3JMDa6nEAE33yyNmuOPzz5WMN5XLioq4itDFYsnnhALQ8PhaP6ZNk36/5KKxaJFIhbjxkW7SDTKQJ8+YhVr6PSJE2VujN+6BiRPxTV6Djkk+3uUZaEz61UsknQS19e7Pic/KrP/3zixOPZYt4Tq5MnyHMMsi6AbasgQsbp27hSxGDUqu99x9uxuu6JyWRbPADgCAIjoKWb2c/+D+luvxJ/HkEsscv23vV1aAYq2Knbvlt+ixCJYgI8+WuYozJ4tcwq01QVkWxZA4W4oALjiivD/ltqy6NdP7oHOOg2i17BokbR0k3QSx4nF+eeLqB90kCus/jMpFirWW7Z0z7Lw01ZTIy1CZlnzIcihh0plFoybVeznqeVEZ1jrTHQ9z/TpriWb1A21bFn2glphXHSRRO894QQnFipsvmUBdB0NmAR/2LxffnT+UlI3FCCVenNzthsqX7HwmTRJJsDGjYYaNCh7MESwz7RI5BIL/64Ha6HeG0gQcBU+UfZs5iRon8XUqS5uEOD6LMJmsEZZFoC4ou65Rz7ffbe8+5aF/99C3FDHHCOtU423E/XfUlkWgFQmUZ2peg0dHclcUEC8WNTWSpgOoKtlUUz8+19In4USFAsNPR3VCR8WYLHYDBokltnChZJntIXuWxb6vJJYFps3i1tGo6lGcbPXPXrsseKyVeHS1vXw4bkHgOSiTx+5hr59uzZMqqulf2bHjuRiAXTts/D/m49Y6KjFXKOh+veXRpffv9XcDBx3XPLzJCSXWHDE57DvvQutSEaNcqNpkqLm3Y03ZreMtRWSK4QI0FUstIV24onuAQ8fLgHVjj+++5bF9dfnvp6wcB9pM3VqdMWni8Hs2BHfua3EiYVPmmLhVwr5jm3PZVl0dsq2FFqLeTFhggy0OOool/fHjROr7ZhjXF5NYlkA8pyiOrfDuO46WVZYz6OWxfvfL3MTCp1PUFPTNYwOIPdchwQnEQs9v4rFYYeJwBVqWaglGRzeP3CgCzek5Vbv+aZN8oxKbFnsR0RfhlgR+hmZ73lccQ9EM0a+LihAHkJbW9fKaPDg7KVQoxZUCv6vsVE6GP1lFisr3cQc7dwDCrMs4vDHy5eKJ5/M/XtdnWT4fCyLsHVBwuiploW/nK7vGtG8c9hhpbEgcjFhgvQV+MEQjztO3Kaazp/8JN5a8O9TPpb9KafIa/Zs+d7eLs962jTg8ccLH/Hjz7vwGTw4O/R/HGoB6AQ+fzEibdDkIxYf/rCsWBl0P2r9tWFD9lwuQCb/7dlTcrG4DUB1yGcA+GXRU1JKtCIpRCz8//toAdchtUn7LCorJQxDFPl0cBciFuVwQ8X5lVUskloWNTUySqejQwrp9u3R92q//WSfpMNK88E/ZyEVV0ODVE5hIVvC+itKjfYvBCPn+nn0C1+IP46fT/OxLBStLDdulBZ2XZ1Yq7nmieTCHyXnU10tQ5OB/N1QQerqZARXPmJB5JZP9tGy2tLS1bLQ9JZSLJg5JhBOPER0CoCfAKgA8Etm/n+B378M4DMAOgC0ALioJOt7a2bLNWw2X1Qs1GwNE4uBA/NvHWowvM7OaDdU2BKnScl3TH4p8OdFJOGyy4BzzpEYVl/6Um7LYvBgGVAQHBVTDAYPFiHq6ChMLIYN6yoWWvEUY9Jgd/ngB8WPrmu9FEqhloXiV5aA3Pc//7nwBs8994RHLKiuTj4kF8gtFkOHSr9KXH9OEnzLIjgnR2fzpyAWqdm1mVAhtwA4FcBUABcS0dTAbq8AaGTmwwD8L4AfpJWeLLprWYShD03nX4SJRSGZWVtuRNG+4LPPlsibcSNLwjjpJBlxkUblWShamSR1Q330o8Bpp8mEw1WrxJUT5oNWjj02Ophbd/DXjyjEygtzkfUksTj6aOng7m74CL03wdF+SQmKRVWVDLEu9Jn6qwH6+O7AJGJx7LFuCHGQujp5vsVwJYa5obRuULEo1MrKQZpO0BkAljBzEzPvAnAPgJn+Dsz8NDNrMPwXAIxCKZg0SVrUOimuGGjGyiUWhXYiq081KqNVVcmqZIXQp48L99BTyNeyIBJfeXu7iGYuyyJttFIp1LIAssWirk6ekR9FoLejz7cQFxQQblmkQZKRiD4zZshM8rD0HHhg4dcbRK+/tbVrn8XatSJ+OmqziKQULQ4AcACAld73VQCOjtgXAC4G8GjYD0R0CYBLAGBMMRRz4kQXMrxY+GKh0TiV/v2TTRiKoqYmfE7C3oq28pJaFoBYVbrgTDnFQivCQsRC55X4I/Quv1w6kdOqEMuBNqQKrTy1Ze1bFmng3/Puuo9uvlkmzBUDvX7m7AjVlZUyHDmlUXOJxYKIjgFwPYABAH7MzA/m/kdyiOgTABoBhA4OZuZbAdwKAI2NjT1z2K7fZxEM3EXkFoophJqa6CieeyPTpokZnW+8pnHjyi8WKhKFVC6XXSatUz/vjBhRnLhVPYm+fSUUydG52o45CIpF2pZFXHjyJBQhNtM/8F2sfp0yZIjck1KLBRGNYOa13qYvAzgLMnT2Rcgs7lysBuA770ZltgXPcxKAbwI4jpmLJL1lQDNW1BjnmprCK7CZM6XFsK/wsY/JK1/Gj3frGpTTshgypDBLcOJEee0LaOyvQujTR9wsaVsWwTD2PQU/bwfjXJVDLAD8nIjmAfgBM78LYDOAcwB0AtiS43/KHAATiWgcRCQuAJBVAxDR4QB+AeAUZi6yX6jEaOumoyO8VVlTU3im/vrXC0/XvsS4ccD998vnconFueeWf/Jcb6C7IbMHDXKu5LQsCz1uTxOLXJYFUHqxYOYzieh0AP9HRL8F8EVIZT8IQMRKNFn/7yCiKwA8Dhk6ezszLyCiGwDMZeZZAG4CMBjAH0gyzwpmPqO7F1UW4vyb3/lOOiNwDIfvAy+XWMyc6SKyGukxcKCb07SvWRZRYqH1Szn6LJj5T0T0CIDLADwA4EZmfjbpwZn5EQCPBLZd532OmerZi4gTizN6pwb2Kvyh0OUSC6M0+M837T6LniYWudxQQGpiETl0lojOIKKnATwG4A0A5wOYSUT3EFERZ7PtJfijWIox8cbIn55gWRilIarCLCbmhsoil2XxHchciYEAHmfmGQCuJqKJAG6E9EEYPoMHy/A4E4vy4BcSE4u9G60wBwxIb1h5T7Us/DkUfj4fNUrmJsUFciyQXJPy2gB8FMDZAP7R+czMi5nZhCIMbYmYWJSHAQMkjDZgYrG3o883zWjJPVUsiJxY+td/7bVuLfIUyCUWZwGoh1gfBYxj3AcxsSg/6ooysdi70eeb5mTFujqpmPOZHFoqwsRi8GCxLlIi12ioDQB+mtqZ90a0JWJiUT7Gj5dYVyYWezelsCwaGiQc+uGHp3eOQgkTi5RJM9zHvodZFuXHLIt9A60s0w6Douth9zRKIZYBTCyKiYlF+fnIR2TN6r0tRIaRTRkqyx6FWRa9HBOL8nPUUcATT5Q7FUbalKLPoidTBrEs8zqNexkmFoZRGsyykHcTi16KdXAbRmkoVZ9FT0Wvv4R9cyYWxeTII+W1r7Z2DKNU7OuWxaBBMqw3hUWOojCxKCbnnw/Mndv9iJqGYeRmX++zGDhQhLKEdY2JhWEYvY8y+Ox7FKNGha/1nSImFoZh9D72dcviuutk8mkJsaGzhmH0Pvb1PosBA0raXwGYZWEYRm9kX7csyoCJhWEYvQ8N7pdi4DwjG3NDGYbR+5g4EVixouSdvPsyZlkYhtE7MaEoKSYWhmEYRiwmFoZhGEYsJhaGYRhGLCYWhmEYRiwmFoZhGEYsJhaGYRhGLCYWhmEYRiwmFoZhGEYsqYoFEZ1CRIuIaAkRXRPy+weIaB4RdRDROWmmxTAMwyic1MSCiCoA3ALgVABTAVxIRFMDu60A8CkAd6WVDsMwDKP7pBkbagaAJczcBABEdA+AmQDe1B2YeXnmt84U02EYhmF0kzTdUAcAWOl9X5XZZhiGYfQyekUHNxFdQkRziWhuS0tLuZNjGIaxz5GmWKwG4IeFHJXZljfMfCszNzJzY0NDQ1ESZxiGYSQnTbGYA2AiEY0jokoAFwCYleL5DMMwjJRITSyYuQPAFQAeB/AWgPuYeQER3UBEZwAAER1FRKsAnAvgF0S0IK30GIZhGIWT6kp5zPwIgEcC267zPs+BuKcMwzCMHkyv6OA2DMMwyouJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEYuJhWEYhhGLiYVhGIYRi4mFYRiGEUuqYkFEpxDRIiJaQkTXhPzen4juzfz+IhGNTTM9hmEYRmGkJhZEVAHgFgCnApgK4EIimhrY7WIArcw8AcCPAHw/rfQYhmEYhZOmZTEDwBJmbmLmXQDuATAzsM9MAHdkPv8vgBOJiFJMk2EYhlEAfVM89gEAVnrfVwE4OmofZu4gojYA9QA2+DsR0SUALsl83UlEb6SS4uIyDIHr6KFYOotHb0gjYOksNr0lnZO78+c0xaJoMPOtAG4FACKay8yNZU5SLJbO4tIb0tkb0ghYOotNb0pnd/6fphtqNYDR3vdRmW2h+xBRXwBDAGxMMU2GYRhGAaQpFnMATCSicURUCeACALMC+8wC8MnM53MA/IWZOcU0GYZhGAWQmhsq0wdxBYDHAVQAuJ2ZFxDRDQDmMvMsAL8CcCcRLQGwCSIocdyaVpqLjKWzuPSGdPaGNAKWzmKzT6STrCFvGIZhxGEzuA3DMIxYTCwMwzCMWHqVWMSFDykHRDSaiJ4mojeJaAERXZXZfj0RrSai+ZnXaRtl9/4AAAcNSURBVD0grcuJ6PVMeuZmttUR0ZNEtDjzPrTMaZzs3bP5RLSFiL7YE+4nEd1OROv9eT5R94+EmzN59TUiOqLM6byJiBZm0vIAEdVmto8loh3eff15mdMZ+ZyJ6BuZ+7mIiD5c5nTe66VxORHNz2wvy/3MUQ8VL38yc694QTrJlwIYD6ASwKsApvaAdO0P4IjM52oAb0PCm1wP4CvlTl8grcsBDAts+wGAazKfrwHw/XKnM/DM1wI4sCfcTwAfAHAEgDfi7h+A0wA8CoAAHAPgxTKn82QAfTOfv++lc6y/Xw+4n6HPOVOmXgXQH8C4TF1QUa50Bn7/LwDXlfN+5qiHipY/e5NlkSR8SMlh5neYeV7m81YAb0FmpvcW/JArdwA4s4xpCXIigKXM3FzuhAAAMz8LGbXnE3X/ZgL4LQsvAKglov3LlU5mfoKZOzJfX4DMeyorEfczipkA7mHmncy8DMASSJ2QOrnSmQlPdB6Au0uRlihy1ENFy5+9SSzCwof0qEqZJGru4QBezGy6ImPi3V5u904GBvAEEb1MEkIFAIYz8zuZz2sBDC9P0kK5ANmFsKfdTyD6/vXk/HoRpFWpjCOiV4jor0T0/nIlyiPsOffU+/l+AOuYebG3raz3M1APFS1/9iax6NEQ0WAA9wP4IjNvAfAzAAcBmA7gHYipWm7ex8xHQCIBX05EH/B/ZLFPe8RYapKJnGcA+ENmU0+8n1n0pPsXBRF9E0AHgN9nNr0DYAwzHw7gywDuIqKacqUPveA5B7gQ2Q2ast7PkHroH3Q3f/YmsUgSPqQsEFE/yAP6PTP/EQCYeR0z72HmTgC3oUQmcy6YeXXmfT2AByBpWqfmZ+Z9fflSmMWpAOYx8zqgZ97PDFH3r8flVyL6FIB/BvDxTMWBjFtnY+bzy5C+gEnlSmOO59wT72dfAB8FcK9uK+f9DKuHUMT82ZvEIkn4kJKT8Vn+CsBbzPxDb7vv/zsLQFkj5RJRFRFV62dIh+cbyA658kkAD5UnhV3IarH1tPvpEXX/ZgH418yok2MAtHnugJJDRKcA+BqAM5h5u7e9gWTtGRDReAATATSVJ5U5n/MsABeQLJg2DpLOl0qdvgAnAVjIzKt0Q7nuZ1Q9hGLmz1L32nezx/80SC//UgDfLHd6Mml6H8S0ew3A/MzrNAB3Ang9s30WgP3LnM7xkNEkrwJYoPcPEhL+KQCLAfwZQF0PuKdVkICSQ7xtZb+fEPF6B8BuiI/34qj7Bxllcksmr74OoLHM6VwC8VFrHv15Zt+zM/lhPoB5AE4vczojnzOAb2bu5yIAp5YznZntvwFwaWDfstzPHPVQ0fKnhfswDMMwYulNbijDMAyjTJhYGIZhGLGYWBiGYRixmFgYhmEYsZhYGIZhGLGYWBh7PUQ0nIjuIqKmTKiTvxPRWZnfGono5szn64noKyH/P56I/q/Ac08momcyEUjfIqJbM9unUw+IRGwYSUltWVXD6AlkJis9COAOZv5YZtuBkFAiYOa5AOammISbAfyImR/KnPs9me3TATQCeCTFcxtG0TDLwtjbOQHALmb+x7oCzNzMzD8FQq2GaRnLYzERfdbbXkNED5OspfBzIupDRBVE9BsieoNknZAvhZx/f8hELj3365kIBDcAOD9jcZyfmWF/OxG9lAlCNzOTvk8R0UMZ62QxEf1bEe+NYSTGLAtjb+cQyEzapBwGie9fBeAVIno4s30GZH2AZgCPQWICLQNwADMfCgCUWVAowI8A/IWIngfwBIBfM/NmIroOMmv2isx/vwvgL8x8UeY4LxHRn71zHwpgO4A5RPRwxiIyjJJhloWxT0FEtxDRq0Q0J2KXh5h5BzNvAPA0XCC7l1jWUtkDCf/wPkjMn/FE9NNM7KUtwYMx868BHAyJnns8gBeIqH/IeU8GcA3JimvPABgAYEzmtyeZeSMz7wDwx8y5DaOkmFgYezsLIKucAQCY+XLIokoNEfsH499w1HZmbgUwDVK5Xwrgl6EHZF7DzLcz80xIePBDQ3YjAGcz8/TMawwzvxWTJsMoGSYWxt7OXwAMIKLPe9sG5dh/JhENIKJ6iCWgFsiMTMTjPgDOB/AcEQ0D0IeZ7wfwLXiipJCsG98v83kEJLDbagBbIctfKo8DuDLTIQ8iOtz77UMkaykPhKx0NjvhtRtG0TCxMPZqWCJlngngOCJaRkQvQZaX/HrEX16DuJ9eAPAfzLwms30OgP+GLFe5DLIeyAEAnsm4jn4H4BshxzsZwBtE9CpEEL7KzGsz55iqHdwA/gNAPwCvEdGCzHflJcg6Ba8BuN/6K4xyYFFnDaMHk1mw6B8d4YZRLsyyMAzDMGIxy8IwDMOIxSwLwzAMIxYTC8MwDCMWEwvDMAwjFhMLwzAMIxYTC8MwDCOW/w+MMxE5a/WeOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 200\n",
    "num_samples = 10000\n",
    "\n",
    "dict_observables = quantum_ising_chain.Convergence(\n",
    "    nn_state, tfim_energy, num_samples, steps\n",
    ")\n",
    "\n",
    "energy = dict_observables[\"energies\"]\n",
    "err_energy = dict_observables[\"error\"]\n",
    "\n",
    "step = np.arange(steps + 1)\n",
    "\n",
    "E0 = -1.2381\n",
    "\n",
    "ax = plt.axes()\n",
    "ax.plot(step, abs((energy - E0) / E0) * 100, color=\"red\")\n",
    "ax.set_xlim(0, steps)\n",
    "ax.set_ylim(0, 0.6)\n",
    "ax.set_xlabel(\"Gibbs Step\")\n",
    "ax.set_ylabel(\"% Error in Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can see a brief transient period in the magnetization observable, before the state of the machine \"warms up\" to equilibrium.  After that, the values fluctuate around the calculated mean.\n",
    "\n",
    "### Adding observables\n",
    "\n",
    "One may also add / subtract observables with the new observable also retaining the same properties in the *Observables* module. For instance, a new observable can be defined by adding the TFIM energy observable multiplied by an arbitrary constant to the PIQuIL observable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_obs = 0.01 * tfim_energy + piquil"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same statistics of this new observable can also be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 0.1295 +/- 0.0014\n",
      "Variance: 0.0188\n"
     ]
    }
   ],
   "source": [
    "new_obs_stats = new_obs.statistics_from_samples(nn_state, new_samples)\n",
    "print(\"Mean: %.4f\" % new_obs_stats[\"mean\"], \"+/- %.4f\" % new_obs_stats[\"std_error\"])\n",
    "print(\"Variance: %.4f\" % new_obs_stats[\"variance\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Template for your custom observable\n",
    "Here is a generic template for you to try using the *Observable* module yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "from qucumber.observables import ObservableBase\n",
    "\n",
    "\n",
    "class YourObservable(ObservableBase):\n",
    "    def __init__(self, your_constants):\n",
    "        super(YourObservable, self).__init__()\n",
    "        self.your_constants = your_constants\n",
    "\n",
    "    def apply(self, nn_state, samples):\n",
    "        # arguments of \"apply\" must be in this order\n",
    "\n",
    "        # calculate your observable for each data point\n",
    "        obs = torch.tensor([42] * len(samples))\n",
    "\n",
    "        # make sure the observables are on the same device and have the\n",
    "        # same dtype as the samples\n",
    "        obs = obs.to(samples)\n",
    "\n",
    "        # return a torch tensor containing the observable values\n",
    "        return obs"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}