diff --git a/docs/api/ffsim.tenpy.rst b/docs/api/ffsim.tenpy.rst new file mode 100644 index 000000000..7e6361615 --- /dev/null +++ b/docs/api/ffsim.tenpy.rst @@ -0,0 +1,7 @@ +ffsim.tenpy +=========== + +.. automodule:: ffsim.tenpy + :members: + :special-members: + :show-inheritance: diff --git a/docs/api/index.md b/docs/api/index.md index 99a780e17..21c172dee 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -9,5 +9,6 @@ ffsim.linalg ffsim.optimize ffsim.qiskit ffsim.random +ffsim.tenpy ffsim.testing ``` diff --git a/docs/how-to-guides/index.md b/docs/how-to-guides/index.md index 57c22a882..155d4c147 100644 --- a/docs/how-to-guides/index.md +++ b/docs/how-to-guides/index.md @@ -4,6 +4,7 @@ :maxdepth: 1 lucj +lucj_mps entanglement-forging fermion-operator qiskit-circuits diff --git a/docs/how-to-guides/lucj_mps.ipynb b/docs/how-to-guides/lucj_mps.ipynb new file mode 100644 index 000000000..8dea5d9aa --- /dev/null +++ b/docs/how-to-guides/lucj_mps.ipynb @@ -0,0 +1,448 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bd5ac3333ca6e15b", + "metadata": {}, + "source": [ + "# How to simulate the LUCJ ansatz using matrix product states" + ] + }, + { + "cell_type": "markdown", + "id": "bdf3ae858d82fccb", + "metadata": {}, + "source": [ + "Following from the previous guide, we now show how to use ffsim to simulate the [LUCJ ansatz](../explanations/lucj.ipynb) using matrix product states. In this way, we can calculate an approximation to the LUCJ energy, which is itself an approximation to the ground state energy. This is particularly useful in complicated cases, such as for large molecules, where even the LUCJ energy cannot be computed exactly. \n", + "\n", + "As before, let's start by building the ethene molecule." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7561238774dbb8b", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:20.643194Z", + "start_time": "2024-10-27T13:48:19.868560Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "converged SCF energy = -77.8266321248744\n", + "Parsing /tmp/tmp942ed7ir\n", + "converged SCF energy = -77.8266321248744\n", + "CASCI E = -77.8742165643862 E(CI) = -4.02122442107773 S^2 = 0.0000000\n", + "norb = 4\n", + "nelec = (2, 2)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Overwritten attributes get_ovlp get_hcore of \n", + "/home/bart/PycharmProjects/ffsim/.ffsim_dev/lib/python3.12/site-packages/pyscf/gto/mole.py:1294: UserWarning: Function mol.dumps drops attribute energy_nuc because it is not JSON-serializable\n", + " warnings.warn(msg)\n", + "/home/bart/PycharmProjects/ffsim/.ffsim_dev/lib/python3.12/site-packages/pyscf/gto/mole.py:1294: UserWarning: Function mol.dumps drops attribute intor_symmetric because it is not JSON-serializable\n", + " warnings.warn(msg)\n" + ] + } + ], + "source": [ + "import pyscf\n", + "import pyscf.mcscf\n", + "\n", + "import ffsim\n", + "\n", + "# Build an ethene molecule\n", + "bond_distance = 1.339\n", + "a = 0.5 * bond_distance\n", + "b = a + 0.5626\n", + "c = 0.9289\n", + "mol = pyscf.gto.Mole()\n", + "mol.build(\n", + " atom=[\n", + " [\"C\", (0, 0, a)],\n", + " [\"C\", (0, 0, -a)],\n", + " [\"H\", (0, c, b)],\n", + " [\"H\", (0, -c, b)],\n", + " [\"H\", (0, c, -b)],\n", + " [\"H\", (0, -c, -b)],\n", + " ],\n", + " basis=\"sto-6g\",\n", + " symmetry=\"d2h\",\n", + ")\n", + "\n", + "# Define active space\n", + "active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)\n", + "\n", + "# Get molecular data and molecular Hamiltonian (one- and two-body tensors)\n", + "scf = pyscf.scf.RHF(mol).run()\n", + "mol_data = ffsim.MolecularData.from_scf(scf, active_space=active_space)\n", + "norb = mol_data.norb\n", + "nelec = mol_data.nelec\n", + "n_alpha, n_beta = nelec\n", + "mol_hamiltonian = mol_data.hamiltonian\n", + "\n", + "# Compute FCI energy\n", + "mol_data.run_fci()\n", + "\n", + "print(f\"norb = {norb}\")\n", + "print(f\"nelec = {nelec}\")" + ] + }, + { + "cell_type": "markdown", + "id": "c0bd6bd083d51e00", + "metadata": {}, + "source": [ + "Since our molecule has a closed-shell Hartree-Fock state, we'll use the spin-balanced variant of the UCJ ansatz, [UCJOpSpinBalanced](../api/ffsim.rst#ffsim.UCJOpSpinBalanced). We'll initialize the ansatz from t2 amplitudes obtained from a CCSD calculation and we'll restrict same-spin interactions to a line topology, and opposite-spin interactions to those within the same spatial orbital, which allows the ansatz to be simulated directly on a square lattice.\n", + "\n", + "The following code cell initializes the LUCJ ansatz operator." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "435b6d06934db617", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:20.978075Z", + "start_time": "2024-10-27T13:48:20.654739Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " does not have attributes converged\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "E(CCSD) = -77.87421536374028 E_corr = -0.04758323886584166\n" + ] + } + ], + "source": [ + "from pyscf import cc\n", + "\n", + "# Get CCSD t2 amplitudes for initializing the ansatz\n", + "ccsd = cc.CCSD(\n", + " scf,\n", + " frozen=[i for i in range(mol.nao_nr()) if i not in active_space],\n", + ").run()\n", + "\n", + "# Construct LUCJ operator\n", + "n_reps = 1\n", + "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n", + "pairs_ab = [(p, p) for p in range(norb)]\n", + "interaction_pairs = (pairs_aa, pairs_ab)\n", + "\n", + "lucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n", + " ccsd.t2, n_reps=n_reps, interaction_pairs=interaction_pairs\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "e2a567f699df4868", + "metadata": {}, + "source": [ + "## Convert the Hamiltonian to a matrix product operator (MPO)" + ] + }, + { + "cell_type": "markdown", + "id": "2824dff2829fccbf", + "metadata": {}, + "source": [ + "Currently, our Hamiltonian is an instance of the `MolecularHamiltonian` class. Using the `from_molecular_hamiltonian` method from the `MolecularHamiltonianMPOModel` class, we can convert this to a TeNPy `MPOModel`, which respects the fermionic symmetries. We can then obtain the MPO using the `H_MPO` attribute and use this `MPO` object as outlined in the [TeNPy MPO documentation](https://tenpy.readthedocs.io/en/latest/reference/tenpy.networks.mpo.MPO.html#tenpy.networks.mpo.MPO). For example, the `MPO` class attribute `chi` tells us the MPO bond dimension, which is an important indicator of how complicated the Hamiltonian is in an MPO representation." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "7faac9a01ef5ba0a", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:21.712055Z", + "start_time": "2024-10-27T13:48:20.997730Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "original Hamiltonian type = \n", + "converted Hamiltonian type = \n", + "maximum MPO bond dimension = 54\n" + ] + } + ], + "source": [ + "from tenpy.models.molecular import MolecularModel\n", + "\n", + "print(\"original Hamiltonian type = \", type(mol_hamiltonian))\n", + "model_params = dict(\n", + " one_body_tensor=mol_hamiltonian.one_body_tensor,\n", + " two_body_tensor=mol_hamiltonian.two_body_tensor,\n", + " constant=mol_hamiltonian.constant,\n", + ")\n", + "hamiltonian_mpo_model = MolecularModel(model_params)\n", + "hamiltonian_mpo = hamiltonian_mpo_model.H_MPO\n", + "print(\"converted Hamiltonian type = \", type(hamiltonian_mpo))\n", + "print(\"maximum MPO bond dimension = \", max(hamiltonian_mpo.chi))" + ] + }, + { + "cell_type": "markdown", + "id": "ad645d3446decfa8", + "metadata": {}, + "source": [ + "## Construct the LUCJ circuit as a matrix product state (MPS)" + ] + }, + { + "cell_type": "markdown", + "id": "5f989277d7cbbca8", + "metadata": {}, + "source": [ + "Our wavefunction ansatz operator, on the other hand, is an instance of the `UCJOpSpinBalanced` class. In a future guide, we will show in detail how we can use such an ansatz to build and transpile Qiskit quantum circuits. In this guide, we will use this ansatz operator to construct our wavefunction as a TeNPy MPS, which respects the fermionic symmetries. Behind the scenes, this executes the ansatz as a fermionic circuit using the TEBD algorithm. \n", + "\n", + "We can pass the `options` dictionary to the TEBD engine to control the accuracy of our MPS approximation, which is detailed in the [TeNPy TEBDEngine documentation](https://tenpy.readthedocs.io/en/latest/reference/tenpy.algorithms.tebd.TEBDEngine.html#tenpy.algorithms.tebd.TEBDEngine). The most relevant key for us in the `options` dictionary is `trunc_params`, which defines the truncation parameters for our quantum circuit. In particular, `chi_max` sets the maximum bond dimension, and `svd_min` sets the minimum Schmidt value cutoff. We also introduce the `norm_tol` parameter in the `apply_ucj_op_spin_balanced` function, which sets the maximum norm error above which the wavefunction is recanonicalized.\n", + "\n", + "In the example below, we set the maximum allowed bond dimension to 15, and after running the circuit, we can see that the maximum bond dimension reaches 15. This indicates that we have most likely truncated the bond dimension with our choice of `chi_max`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e9d8e1b09ee778c2", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:22.585845Z", + "start_time": "2024-10-27T13:48:21.730120Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "wavefunction type = \n", + "MPS, L=4, bc='finite'.\n", + "chi: [4, 15, 4]\n", + "sites: SpinHalfFermionSite('N', 'Sz', 1.000000) SpinHalfFermionSite('N', 'Sz', 1.000000) SpinHalfFermionSite('N', 'Sz', 1.000000) SpinHalfFermionSite('N', 'Sz', 1.000000)\n", + "forms: (0.0, 1.0) (0.0, 1.0) (0.0, 1.0) (0.0, 1.0)\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "from tenpy.algorithms.tebd import TEBDEngine\n", + "\n", + "import ffsim\n", + "from ffsim.tenpy.gates.ucj import apply_ucj_op_spin_balanced\n", + "from ffsim.tenpy.util import statevector_to_mps\n", + "\n", + "# Construct Hartree-Fock state\n", + "psi_mps = statevector_to_mps(ffsim.hartree_fock_state(norb, nelec), norb, nelec)\n", + "\n", + "# Construct the TEBD engine\n", + "options = {\"trunc_params\": {\"chi_max\": 15, \"svd_min\": 1e-6}}\n", + "eng = TEBDEngine(psi_mps, None, options)\n", + "\n", + "# Apply the LUCJ operator\n", + "apply_ucj_op_spin_balanced(eng, lucj_op)\n", + "\n", + "# Print the wavefunction\n", + "psi_mps = eng.get_resume_data()[\"psi\"]\n", + "print(\"wavefunction type = \", type(psi_mps))\n", + "print(psi_mps)" + ] + }, + { + "cell_type": "markdown", + "id": "6c97e4db54214ecd", + "metadata": {}, + "source": [ + "## Compare the energies" + ] + }, + { + "cell_type": "markdown", + "id": "9c8924340fc05c75", + "metadata": {}, + "source": [ + "Now that we have converted our `MolecularHamilonian` to an MPO, and our LUCJ ansatz to an MPS, we can contract the tensors to compute an estimate of the ground-state energy. In order of increasing accuracy, we can compare the LUCJ (MPS) energy, the LUCJ energy, and the FCI energy." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a6a7d85060f3d8a2", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:22.711608Z", + "start_time": "2024-10-27T13:48:22.629846Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "LUCJ (MPS) energy = -77.78140377489262\n", + "LUCJ energy = -77.84651018653342\n", + "FCI energy = -77.87421656438624\n" + ] + } + ], + "source": [ + "# Compute the LUCJ (MPS) energy\n", + "lucj_mps_energy = hamiltonian_mpo.expectation_value_finite(psi_mps)\n", + "print(\"LUCJ (MPS) energy = \", lucj_mps_energy)\n", + "\n", + "# Compute the LUCJ energy\n", + "hf_state = ffsim.hartree_fock_state(norb, nelec)\n", + "lucj_state = ffsim.apply_unitary(hf_state, lucj_op, norb, nelec)\n", + "hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec)\n", + "lucj_energy = np.vdot(lucj_state, hamiltonian @ lucj_state).real\n", + "print(\"LUCJ energy = \", lucj_energy)\n", + "\n", + "# Print the FCI energy\n", + "fci_energy = mol_data.fci_energy\n", + "print(\"FCI energy = \", fci_energy)" + ] + }, + { + "cell_type": "markdown", + "id": "76da0123-c376-484e-9f78-231d049fc051", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-27T13:48:22.764638Z", + "start_time": "2024-10-27T13:48:22.762324Z" + } + }, + "source": [ + "To illustrate the effects of the truncation parameters more clearly, we can plot the energies at different values of `svd_min` and `chi_max`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "bf98d538-c182-4ede-917f-1eed31969c9a", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAFzCAYAAABcurqFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACXL0lEQVR4nOzdeViVdf7/8ec5hx0ERBZFUNwVF1BR1LK0nBwzWyyzmiZ12svKn9WkTWk1pTNt43zLxsYyayqzTTMrW8yyRUUhTBTNDNyQTRFkX879+4M8SYACAjcHXo/rui849/mc+34dwXN4n89yWwzDMBAREREREZFmYzU7gIiIiIiISFujQkxERERERKSZqRATERERERFpZirEREREREREmpkKMRERERERkWamQkxERERERKSZqRATERERERFpZirEREREREREmpmL2QFaA7vdTlpaGu3atcNisZgdR0TEqRiGwYkTJwgNDcVq1eeDLYHe10REGq6u72sqxBpBWloa4eHhZscQEXFqBw8eJCwszOwYgt7XREQaw5ne11SINYJ27doBlf/Yvr6+JqcREXEueXl5hIeHO15LxXx6XxMRabi6vq+pEGsEJ4dt+Pr66g1LRKSBNASu5dD7mojI2TvT+5oG44uIiIiIiDQzFWIiIiIiIiLNTIWYiIiIiIhIM9McMRFpVhUVFZSVlZkdQ5qRzWbDxcVFc8Ca2dq1a7n33nux2+088MAD3HTTTWZHEhGRU6gQE5Fmk5+fz6FDhzAMw+wo0sy8vLzo1KkTbm5uZkdpE8rLy5k9ezYbNmzAz8+PoUOHcsUVV9ChQwezo4mIyK9UiIlIs6ioqODQoUN4eXkRFBSk3pE2wjAMSktLycrKIiUlhV69eumizc0gLi6O/v3707lzZwAmTJjAZ599xrXXXmtyMhEROUmFmIg0i7KyMgzDICgoCE9PT7PjSDPy9PTE1dWV/fv3U1paioeHh9mRWryNGzfy1FNPER8fz5EjR1i1ahWXX355lTaLFy/mqaeeIj09naioKJ577jmGDx8OVF6Q+WQRBtC5c2cOHz7cnE9BRETOQB9LikizUk9Y26ResPopKCggKiqKxYsX13j/ypUrmT17NvPnzychIYGoqCjGjx9PZmZmMycVEZGGcop3xq+++gqLxVLjtnXrVgAeeeSRGu/39vY+7bEPHDjAxIkT8fLyIjg4mPvvv5/y8vLmeFoAZBzaR9J3H5JxaF+znVNERFq2CRMm8Pjjj3PFFVfUeP+zzz7LzTffzIwZM4iMjGTJkiV4eXmxbNkyAEJDQ6v0gB0+fJjQ0NBaz1dSUkJeXl6VraH0viYiUjdOUYiNGjWKI0eOVNluuukmunXrRkxMDAD33XdftTaRkZFMmTKl1uNWVFQwceJESktL+f7773n11VdZvnw58+bNa5bnteXdZwlaOpQBn19P4NKhxL23qFnOKyIizqu0tJT4+HjGjRvn2Ge1Whk3bhybNm0CYPjw4SQlJXH48GHy8/P55JNPGD9+fK3HXLhwIX5+fo4tPDy8Qdni3ltU5X1ty9tPN+g4IiJtgVPMEXNzc6Njx46O22VlZXzwwQfcddddjmFOPj4++Pj4ONps376dXbt2sWTJklqP+9lnn7Fr1y6++OILQkJCiI6O5u9//zsPPPAAjzzySJOu7pVxaB/DdjyG1VK5epzNYjDkx0fJiJ1ESFiPJjuviIg4t+zsbCoqKggJCamyPyQkhN27dwPg4uLCM888w9ixY7Hb7fz1r3897YqJc+fOZfbs2Y7beXl59S7GMg7tY+iPj1R5Xxu+8+9kPvI82a6dKfDqTLlvF1w6ROAd0oOAsF4EdYrA5uIUf4qIiDQ6p3z1W7NmDUePHmXGjBm1tnnppZfo3bs3o0ePrrXNpk2bGDhwYJU3s/Hjx3P77bezc+dOBg8eXOPjSkpKKCkpcdxuyBCOrP27CLFUXcLbxWIne/9uFWIiInLWLr30Ui699NI6tXV3d8fd3f2szlfT+5rFAsHkEFyWA7lJkAsc/O3+UsPGEWswOW6dKPQOw+4XjmtgN9qF9CAwvDcBQaFYNL9QRFoppyzEXn75ZcaPH09YWFiN9xcXF/PGG28wZ86c0x4nPT29xk8UT95Xm4ULF/Loo4/WM3VVQV0jqTAs2E550yo3rAR27XtWxxURaQ7Hjx9n3LhxlJeXU15ezj333MPNN99sdqw2ITAwEJvNRkZGRpX9GRkZVUaPNLea3tcqDCs/nvs8FSUFlGenYs3dj1fhIdqXHiHYno2bpYIw4whhJUegJAGOASm/HbPQcCfTFsJx91BKfMIw/LviHtQd3049Ce7Sm3Z+AbXmyTi0j6z9uwjqGqkPONG/h0hLZGohNmfOHP75z3+etk1ycjJ9+/5WnBw6dIhPP/2Ut99+u9bHrFq1ihMnTjBt2rRGy3qqxhjCERLWg7hBjzD0x0ewWQwMA7b1f5ARenEUESfQrl07Nm7ciJeXFwUFBQwYMIDJkyfrgsHNwM3NjaFDh7J+/XrHkvZ2u53169czc+ZM03KdfF8b8uOjuFjslBtWEgbNZ/gf/lRj+4ryctLTUjh6aC8FGfuoOJaKS94BfAoP06HsCIFGDl6WEiLsB6DoABQBWcDe346RizdZto7keXSmtF04lvZd8QjuTmFqPMNTlxBiMagwLMQNeoThV85qjn+GFinuvUUM/fER/XuItDCmFmL33nsv06dPP22b7t27V7n9yiuv0KFDh9MOt3jppZe45JJLqvV2/V7Hjh2Ji4ursu/kJ4yn+1SxMYZwAAy/chZHYi7G45WxtLfk4xnYsMnRItL05syZw7/+9S+uvPJK3nzzTbPjmM5ms+Hl5QVUDtc2DAPDMM7wKKmr/Px8fv75Z8ftlJQUEhMTCQgIoEuXLsyePZtp06YRExPD8OHDWbRoEQUFBacdst8chl85i4zYSWTv301g174MP82HizYXFzp26UXHLr1qvL+kuJDMg3vJOfwzRZm/YM85gNuJA7QrTiOoPJ325OFHAX4V+6BgHxQA6UDyrwf49UoZNovBsB/ns3PPuxS7B1Hu7o/d3Q+Lpz82r/a4+HTAzScAT98OePl1wLd9EJ5e7RptSGRT9kTZKyooyM+lIPcoRSeOUZx/nNKCHMryj1NRlIu9OBdL3mFis1djOeXfI+bH+Ww+FA/eQRhu3ljdvCo3d29s7t64eHjj6uGDm6c3bp6+uHv54OHlg6dXu7Oe06eeOZHfmFqIBQUFERQUVOf2hmHwyiuvcMMNN+Dq6lpjm5SUFDZs2MCaNWvOeLyRI0fyxBNPkJmZSXBwMACff/45vr6+REZG1jnX2ejUtTdbAscTm/0eJT9+ABdc0yznFZH6mTt3LmFhYdx111089thj9OzZ0+xItarLxYDh9BcErovjx49z/vnns3fvXp566ikCAwMb8Vm0bdu2bWPs2LGO2ydHYUybNo3ly5czdepUsrKymDdvHunp6URHR7Nu3bozfgDZHELCejTKH9juHl6E94oivFdUjffn5+WQdfAnctP2UZz1Cxzfj3v+IYIL99GZqsM2LRboX7oDSoETZz53qWHjhMWHAqsPhVZfSlzbUebqS4W7P3YPf0cR5+rTAfd2lUWct18g7doH4eH522VzTtcTZdjtFBcVUJB7jIK8oxTl51Caf5yyghzKC3MxinMxivOwluRiKz2BS3k+bmUn8KjIx9NegDcF+BhFtLMYtDvTE/rd5RutFhiRswZyzvxv8XvFhivFFneK8aDU6k6pxYNSqyflNg/KbZ5UuHhid/HEcPECVy8MV08svxZ7xuF4hh398Ld/j4HzGX7V/6t/CJFWwmI40UeY69evZ9y4cdWGK57q4YcfZtmyZRw4cACbzVblvlWrVjF37lzHqlIVFRVER0cTGhrKk08+SXp6On/+85+56aabWLBgQZ1z5eXl4efnR25uLr6+vvV+Xju/+4j+n19HLt54PZiCq9vZ97aJtDTFxcWkpKTQrVs3PDw8zI7TIEVFRfj4+PDuu+/Wen2nluCTTz7hu+++Y+jQoUyePLnGQmzlypXccMMNLFmyhNjYWBYtWsQ777zDnj17HB9MRUdH13hdxc8++6zKNakyMjKYPHky77//fq2FwOl+/mf7GiqNz5l/JhmH9hG4dOjv5qpZiOv9/7BYrFCYg6X4OLbSXFxLc/EoP4Gn/QTe9nx8jXxcLRVndf5iw5UTFh+K8CDcOMKp17A3DEi3BOFBCT5GwVmf66RSw0aBxZsCizdFVm9KXHwoc/GhzNUXAwvDcj7CekoOu2EhrsOlWKw2rGWF2CoKsVUU41JRhJu9GFd7Me72YtwpwcMowctSUvvJz4JhQKo1nByvbpT4dsHWoQdenXoR2KUfwaHdsP7u7zgRZ1HX11CnWqzj5ZdfZtSoUbUWYXa7neXLlzN9+vRqRRhAbm4ue/bscdy22WysXbuW22+/nZEjR+Lt7c20adN47LHHmuw51KRv7HiOfe5LAHns2PwxA89ruX/gibQER3KLSMkuoFugN538PJvtvOXl5Xh5eZGUlNSiC7EJEyYwYcKE07Y59YLAAEuWLOGjjz5i2bJljoWOEhMT63S+kJAQoqKi+Oabb7jqqqvOKrvI2aptrtrIOsyJMux2CgryOHE8m8LcoxTlZlOSf4yKgmNUFOZgFOVgLT6OS2kermV5eJTn4VVxAm+jsoizWQw8LGV4nOxq+l1PlMUCnciqcl+FYSHf4uUoooptPpS6tKPc1Qe7azvsHn5YPHyxevhh8/LDzbs9bt7+eLQLwNsvAB/fANw9vGhvtdK+lucV996iav8eI+oxR8xeUUFJcSFFBXkUF+ZTVpRPSVHl1/KSAsqLC7CXFmAvKcAoLcAoK8JSWoClvAhreRG28kK8i4/Qp/ynav8e3YyDdCs4WDm09AiQVHlfieFKuq0jOR5hFLfriiWgG54hvejQpS8h4b1wcW26SwyJNBenKsTONC/DarVy8ODBWu+fPn16tTlpXbt25eOPP26MeA1mc3Fhb8AYYo+toXD7KlAhJm2AYRgUldX/0+D34g8xf81O7Ebl8JpHL+3PlUNrXkG1Np6uNsc1COvjoYceIj8/n6SkpHo/tiEWLFhwxt75Xbt20aVLl3od9+QFgefOnevY9/sLAp9JRkYGXl5etGvXjtzcXDZu3Mjtt99erxwiTaU+c9VOZbFa8W7nj3c7fwiv3/Bje0UFeSeOk388m8LcbI4f2sWQuPsd11WDyqLrx5H/R/uukXi2a4+3bwDePn74Wa341ets9dPQf4+TrDYbnt7t8PQ+4yDIWmUc2kdFDT2V8YP/jr0oD8uxX/A4sZ+AksN0tGfgbimjq/0gXQsPQuEmyMAx/6/csHLIGswx9zCKfLpgBHTHI7gn7cP6ENK1T5XhoSItmVMVYq2ZV9QVsGENPY9+RUV5uS5wKa1eUVkFkfM+Patj2A14+IOdPPzBzno9btdj4/Fyq9//sfj4eJYsWcLEiRObrRC77bbbuPrqq0/b5tQhgnVVlwsCn8n+/fu55ZZbHIt03HXXXQwcOLDeWUSaSmPNVasrq82Gr38HfP07AH0g6hziiguqryL5xxuaLdOpmvvfo6bz17iq5uV3VWtbXlZK2qFfOHowmcL0vRhHf8H9xH7aFx+iY8URPCxlhBnphBWnQ/E2yAZ+7WyzGxbSLYEcdQulwKcLdv9uuIf0xDe0Dx0j+lYW2WjREGkZ9Nd+C9F35ERyN3jTgVx2xn1G/1EXmx1JRH5lt9u59dZbmTlzJrGxsVx//fWUlZXVumhQTdLS0rj//vt544036vyYgIAAAgJqv06SmYYPH17noYsibdXZ9kS1NnX993BxdSO0W19Cu1WfimKvqCAz/QBZ+5MpOPITFUf34Z63H9+iQ3QsT8PHUkRHsuhYmgXHtldem+6X3x6fjT9FFk/C7EcIsVQWbpu63MTgax/Bw8uniZ65SM1UiLUQrm7u/OR/HsOOf0L+D++BCjFp5Txdbex6bHy9HpOeW8y4Z7/GfsoSQ1YLfDH7fDr61X0BEE/X+k0Af+6558jOzuaxxx7jwIEDlJWVsXv37nr1AIWGhtarCIOmG5rYUi8ILNIamd0T1dKc7b+H1WYjuHM3gjt3A6r+rWTY7RzNSiPrwG5OpP1EefY+XHNTaVd4kODyNNpzgkCOg3HcMUfPajEYeXAp9n8u5YgliCyPLhS16waBvfEO7UtQtwEEh3ZrtMsZiJxKhVgL4jrwcvjmE7plfYm9okKrBUmrZrFY6j08sHuQDwsnD+TB95OoMAxsFgsLJg+ge1DTfYp5+PBhHn74YVasWIG3tze9evXC3d2dpKQkBg4cSGpqKpdddhkDBgwgLi6OcePGMX78eBYuXEhBQQGrVq2iV69epKamctVVV/Huu+9y2WWXER0dTVxcHIMGDeKtt96qcc5aUw1NbKkXBBYRORsWq5UOIWF0CAkDxlW7P/dYFslfvMKIXU9Uu8/660IqnYqzoDi+8uLhv85JKzTcSXMJI9erK6X+PXAN6Y1/eCSdug9wDHUUaQgVYi1Iv3MuJX/jLIItx9idsIG+w6q/iIi0dVOHdeG83kGkZhcSEejV5Ksm3n333UyYMIGJEycC4OLiQr9+/arME0tOTubtt9+mZ8+eDBgwAB8fH7Zs2cKLL77I888/z7///e8qx0xOTmbFihX069ePsWPH8u233zJ69Ohq527o0MQzXQwYaLEXBBYRaSp+AUF0G3UlFTsXVFk0pNywknrlR5QW5pF/OBl71l488n6hQ9F+OtnT8bKU0LNiH5zYBye+hIPAtsrHZhJApnsXCnwiMDr0wqtTXwIjBhAS3lPz/eWM9BvSgrh7eLHD7xxi8r7g+LZ3QYWYSI06+Xk2y7L1a9eu5csvvyQ5ObnK/oEDB1YpxPr06UOfPn0A6NevH+PGjXO0q2lV1j59+jguGj948GBSU1NrLMQa6kwXAwZa9AWBRUSaSq2Lhgwa9WuLP1ZpX1ZawoHUZI7t30Vx+m6sx36mXX4qIWUHCSCPYI4RXHIMShLhKJWLhnxdufz+AVsoOZ5dKfHvjktQb3zD+tGxxyD82gc6jq9FQ9o2FWItjK3/pbDpC7pmfIFht2tMsoiJLrnkEnJycqrtf+2116rcdnf/7SLsVqvVcdtqtVJRUX2J/lPb22y2GtucjTFjxmAYxhnbzZw5U0MRRaTNqc8iKq5u7nTpHU2X3tHV7ss9msGRlCROHEqmPPMn3HN/IaAoldCKI7hbyuhm30+3gv1QsBEOA4mVjzuKHxmu4QD0LU0ixPLrRccHPcLwelzfTZyfCrEWpu+5kyn8/n46WbLY++N39IpuvE/JRURERKRxFlHx6xCCX4cQiLmwyv6K8nIOH9hDdupOio7sxnLsZ3xOpBBcepAgcuhALh3Kcisb/zo92GYxGPLjo2TETlLPWBuiQqyF8fRuR0K7WIbkbyQ77m0VYiIiIiJOxObiQufu/encvX+1+07kHiP9lySyt77HyLTlVe5zsdjJ3r9bhVgbYjHqMn5FTisvLw8/Pz9yc3Px9fU96+Nt+2gpMVvv46AllLCHd2p4orQKxcXFpKSk0K1bNzw86r7UvLQOp/v5N/ZrqJw9/UxEmlbGoX0ELh1abdGQozdvUyHWCtT1NVR/4bdAfUdfRYnhSriRRmryVrPjiIiIiEgjCgnrQfygR7AblWMT7Qa8FTJbRVgbo0KsBfLxbc8u72EApG9+2+Q0IiIiItLYhl85i7wLFgLwo9GNZYWj67TQkrQeKsRaqPI+kwDodPgzk5OIiIiISFPw7zkcgE6W4/ySXcCWlGMmJ5LmpEKshep93tWUGjYi7AfYvyfR7DgiIiIi0tj8IwAIseTgTikr4g6Ym0ealQqxFsqvfSC7PYcAkLbpLZPTiIiIiEij8woANx8AOluy+WRHOjkFpSaHkuaiQqwFK+51CQBBBzU8UURERKTVsVjAvysA53YooLTCzvs/HDY5lDQXFWItWK/zplJuWOlZsY/DvySbHUdEREREGlv7ykJsYpfKnrAVcQe0aEcboUKsBWsf1IndHoMAOPjdCpPTiIiIiEij+7VHbHC7XDxdbfycmc+2/Tkmh5LmoEKshSvoMRGA9vvXmZxERERERBrdrz1ibicOMimqE4AW7WgjVIi1cD3Ouwa7YaFP+R7SD/5sdhwRERERaUy/9oiRs59rh3cB4KMfj5BbWGZiKGkOKsRauMCOXdjjFglA6jdaPVFEWo6UlBTGjh1LZGQkAwcOpKCgwOxIIiLO59ceMY7vJzrcn74d21FSbmfVD4fMzSVNToWYE8jtdjEAvimfmJxEROQ306dP57HHHmPXrl18/fXXuLu7mx1JRMT5nOwRK8rBUnLC0Sv21taDWrSjlVMh5gQiRl8DQN/SnWSna8ywiBnmzJmDu7s71113ndlRWoSdO3fi6urK6NGjAQgICMDFxcXkVCIiTsjdB7w6VH5/fD+XD+6Mu4uV3ekn+OHgcVOjSdNSIeYEOob35CeX3lgtBvs2aniiiBnmzp3LM888w4oVK/j555Y9X3Pjxo1MmjSJ0NBQLBYLq1evrrHd4sWLiYiIwMPDg9jYWOLi4up8jr179+Lj48OkSZMYMmQICxYsaKT0IiJt0CnzxPw8XblkUCgAK7boA/jWTIWYkzjWdQIA3vs+MjmJSNvk5+fHjTfeiNVqZceOHWbHOa2CggKioqJYvHhxrW1WrlzJ7NmzmT9/PgkJCURFRTF+/HgyMzMdbaKjoxkwYEC1LS0tjfLycr755hteeOEFNm3axOeff87nn3/eHE9PRKT1OWWeGMC1w8MBWPvjEfKKtWhHa6VCzEmEn3MtAH2LfyQn64jJaUTapvLycry8vEhKSjI7ymlNmDCBxx9/nCuuuKLWNs8++yw333wzM2bMIDIykiVLluDl5cWyZcscbRITE0lKSqq2hYaG0rlzZ2JiYggPD8fd3Z2LL76YxMTEZnh2IiKt0Ck9YgBDu7anV7APRWUVfJCYZmIwaUoqxJxE5+792GfrjovFzt6NK82OI2Ku3MOQsrHyazN66KGHyM/Pb7ZCbMGCBfj4+Jx2O3Cg/sNWSktLiY+PZ9y4cY59VquVcePGsWnTpjodY9iwYWRmZpKTk4Pdbmfjxo3069ev3llERIRqPWIWi8WxaMebWw5o0Y5WSjOrnUhm+Hh6pP4H971rgVlmxxE5O4YBZYX1f1zim/DJX8Gwg8UKE56E6HouoOHqBRZLvR4SHx/PkiVLmDhxYrMVYrfddhtXX331aduEhobW+7jZ2dlUVFQQEhJSZX9ISAi7d++u0zFcXFxYsGAB5513HoZhcNFFF3HJJZfUO4s0jYMHD/LnP/+ZzMxMXFxcePjhh5kyZYrZsUSkNr/rEQOYPKQz/1i3m+Qjefx4KJeocH9zskmTUSHmREJHXgOp/6FfUQK5Odn4tQ80O5JIw5UVwoL6FxFVGHb4+L7KrT4eTAM37zo3t9vt3HrrrcycOZPY2Fiuv/56ysrKcHV1rfMx0tLSuP/++3njjTfq/JiAgAACAgLq3L65TZgwgQkTJpgdQ2rg4uLCokWLiI6OJj09naFDh3LxxRfj7V3333sRaUbtIyq/Ht9f+UGlxYK/lxsXD+jI6sQ03tp6QIVYK6ShiU6ka59oUq3huFkq+Gnj22bHEWkznnvuObKzs3nssccYOHAgZWVlde45Oik0NLReRRg03dDEwMBAbDYbGRkZVfZnZGTQsWPHeh9PWp5OnToRHR0NQMeOHQkMDOTYsWPmhhKR2vmFAZbKDykLsh27Tw5P/CAxjfyScpPCSVNRj5iTOdJ5PBEHX8Jlz4fAHWbHEWk4V6/Knqn6yEuDxcMre8JOstjgzi3gW4/eNVevOjc9fPgwDz/8MCtWrMDb25tevXrh7u5OUlISAwcOJDU1lcsuu4wBAwYQFxfHuHHjGD9+PAsXLqSgoIBVq1bRq1cvUlNTueqqq3j33Xe57LLLiI6OJi4ujkGDBvHWW29hqWGoZFMNTXRzc2Po0KGsX7+eyy+/HKjs9Vu/fj0zZ86s9/Gk/jZu3MhTTz1FfHw8R44cYdWqVY6fxUmLFy/mqaeeIj09naioKJ577jmGDx9e73PFx8dTUVFBeHh4I6UXkUbn4l75PpZ3uLJXzCcIgOHdAuge5M0vWQWsSUzjutguJgeVxqRCzMmExF4NB18ismAr+Xk5+Pi2NzuSSMNYLPUaHghAYC+Y9G/4cBYYFZVF2KRFlfubyN13382ECROYOHEiUDnkq1+/flXmiSUnJ/P222/Ts2dPBgwYgI+PD1u2bOHFF1/k+eef59///neVYyYnJ7NixQr69evH2LFj+fbbbx0XRj5VQ4cm5ufnV7nWWUpKComJiQQEBNClS+Wb+OzZs5k2bRoxMTEMHz6cRYsWUVBQwIwZM+p9Pqm/k5cY+Mtf/sLkyZOr3X/y8gJLliwhNjaWRYsWMX78ePbs2UNwcDBQeXmB8vLqn5B/9tlnjgL92LFj3HDDDSxdurRpn5CInD3/rpWFWE4qhMUAvy7aMawLT3yczFtbD6gQa2VUiDmZbpHDOGTpRBhHSPrmPYZOvMnsSCLNa8gN0ONCOPYLBHQHv85Ndqq1a9fy5ZdfkpycXGX/wIEDqxRiffr0oU+fPgD069fPsRrhwIED+fjjj6sdt0+fPkRGRgIwePBgUlNTayzEGmrbtm2MHTvWcXv27NkATJs2jeXLlwMwdepUsrKymDdvHunp6URHR7Nu3bpqC3hI0zjT/LpTLy8AsGTJEj766COWLVvGnDlzAM54uYCSkhIuv/xy5syZw6hRo87YtqSkxHE7Ly+vjs9ERBpN+65w4HvHyoknXTk0jKc+3cOPh3JJOpzLgM5+JgWUxqY5Yk7GYrVysNMfKr9PXmNyGhGT+HWGbqObtAgDuOSSS8jJyak2b+q1117jgw8+cNx2d3d3fG+1Wh23rVYrFRUV1Y57anubzVZjm7MxZswYDMOotp0swk6aOXMm+/fvp6SkhC1bthAbG9uoOaRhGuPyAoZhMH36dC644AL+/Oc/n7H9woUL8fPzc2waxihighpWTgQI8HZj/IDK96EVcfWfFywtlwoxJxQ4rHIJ4r4nNlNUcMLkNCIi0phOd3mB9PT0Oh3ju+++Y+XKlaxevZro6Giio6PZsWNHre3nzp1Lbm6uYzt48OBZPQcRaYDfXUvsVNcOq/xw5IPENApLtWhHa6GhiU6oZ9S5HPkgiE6WLBK+Xc2Q8Wf+tFNERNqOc889F7vdfuaGv3J3d6/SUysiJqilRwxgRPcORHTwIvVoIWu3H+HqYeq1bg3UI+aELFYr+0Mqh6zYd642N4xIGxcREcG2bdsct999913GjBkDwIgRI1i7dm2Vdr9v//TTTzN9+vTmjCwtnC4vINJGnewRyz0E9qpD1q1WC9f8upT9mxqe2GqoEHNS/kOvBKBP7neUFBeanEZERBrLqZcXOOnk5QVGjhxpYjIRaVLtOoHVFexllZdr+Z0rh4ThYrWQePA4yUe0oE5roELMSfUeegGZBNDOUsTu77Voh4iIM8nPzycxMdGx8uHJSwycvED37NmzWbp0Ka+++irJycncfvvturyASGtntYH/r0MOa5gnFtTOnYv6V84dfUu9Yq2CCjEnZbXZSAmsXJ665McPztBaRERakm3btjF48GAGDx4MVBZegwcPZt68eUDl5QWefvpp5s2bR3R0NImJibq8gEhbcJp5YgDX/jo88f0fDlNU2rgr7krz02IdTsxn8JXw+Xv0Of41ZaUluLpporWIiDM4eYmB05k5cyYzZ85spkQi0iKcZuVEgHN6BBIe4MnBY0V8tOMIVw0Na8Zw0tjUI+bE+saO5xi++FHA7s3VLxorIiIiIk7kDD1iVquFa4ZV9oppeKLzUyHmxGwuLuwNGANA4fZV5oYRERERkbNzhh4xgClDw7BZLWzbn8NPGbqerDNTIebkvKKuAKDn0a+oKNcF/kREREScln9E5ddaesQAgn09GNcvGIAV6hVzairEnFzfkRPJxZsO5LI77jOz44iIiIhIQ53sETtxBMpLam128ppi7yccprhMi3Y4KxViTs7VzZ2f/M8DIP+H90xOIyIiIiIN5tUBXL0BA44frLXZeb2C6OzvSW5RGeuS0psvnzQqFWKtgOvAywHolvUl9gp9KiIiIiLilCyWU+aJpdbazGa1MHVY5TXH3tTwRKelQqwV6HfOpeQbngRzjJ8SNpgdR0REREQa6gwrJ540JSYMqwXiUo6xLyu/GYJJY1Mh1gq4e3ix2+8cAI5ve9fkNCKt0/nnn4/FYqm23XDDDWZHExGR1qQOKycCdPLz5IK+lYt2aCl756RCrJWw9b8UgK4ZX2DY7SanEWldDMPghx9+4Omnn+bIkSNVthdeeMHseCIi0prUsUcM4NpfF+14N/4QJeWanuJsVIi1En3PnUyh4U4nsvj5x+/MjiPSquzdu5cTJ05w3nnn0bFjxyqbj4+P2fFERKQ1qWOPGMD5vYPo6OtBTmEZn+3MaOJg0thUiLUSnt7t2N0uFoDsuHdMTiNSdwUFBbVuxcXFdW5bVFRUp7YNER8fj4uLC4MGDWrw8xQREamTevSIudisXP3roh26ppjzUSHWitj7VQ5PDDvyuYYnitPw8fGpdbvyyiurtA0ODq617YQJE6q0jYiIqLFdQyQkJFBRUUGHDh2qHOvWW28FIDAwsNpjUlNTiYmJqbJv+vTprF27FoBDhw4xefJkevToQUxMDFOmTCEjI6PW44mISBtxskes6BiUnDhj86nDwrFY4Pt9R0nJbtgHjmIOFWKtSN/RV1FiuBJupJGavNXsOCKtRkJCAtdeey2JiYlVtoULFzboeIZhcNlllzFx4kT27dvHtm3buPvuu8nKymrk5CIi4nTc24FnQOX3degV6+zvyfm9gwB4a6t6xZyJi9kBpPH4+LbnB+9hDC78nvTN79Ctf6zZkUTOKD+/9iV3bTZblduZmZm1trVaq36ulJqaela5TpWQkMATTzxBz549G+V469evx8fHhxtvvNGxb/To0Y1ybBERaQXad63sETu+HzoOOGPza4d34as9WbwXf4h7/9AHNxf1tTgDp/gpffXVVzUuG22xWNi6tbLn55FHHqnxfm9v71qPu337dq699lrCw8Px9PSkX79+/Pvf/26up9UkyntfAkCnw5+anESkbry9vWvdPDw86tzW09OzTm3r65dffuH48eNERUWd1fM81a5duxgyZEijHU9ERFqZeswTA7igbzDB7dzJzi/li2Qt2uEsnKIQGzVqVLUlo2+66Sa6devmmINx3333VWsTGRnJlClTaj1ufHw8wcHBvP766+zcuZO//e1vzJ07l+eff765nlqj633e1ZQaNiLsBzjwU6LZcUScXnx8PAAhISGkp6dX2eynmYtpsVjqtV9ERMShHisnArjarEyJCQO0aIczcYqhiW5ubnTs2NFxu6ysjA8++IC77rrL8UfN7yfib9++nV27drFkyZJaj/uXv/ylyu3u3buzadMm3n//fWbOnNnIz6J5+AUE8aPnEAYVbyXt+7fp0jva7EgiTi0hIQGAXr16Vdnv7u5OXl4ebm5uNT6uQ4cO5OTkVNl37NgxAgMDcXNz4/3332+awCIi4vzq2SMGcM2wLizesI9v9mZz8Fgh4QFeTRROGotT9Ij93po1azh69CgzZsyotc1LL71E79696z3vIjc3l4CAgNO2KSkpIS8vr8rWkhT3mghA4MF1JicRcX4LFy7EMIxqW3Fxca1FGFR+OOTv78/3338PVK6SuGPHDvr378+4cePIy8tj+fLljvbffvstSUlJTf10RETEGdSzRwwgPMCL0b0qV93Voh3OwSkLsZdffpnx48cTFhZW4/3FxcW88cYbVSbC18X333/PypUrueWWW07bbuHChfj5+Tm28PDwep2nqfU67xrKDSs9K/Zx+Jdks+OItHo5OTmEhYU5thUrVgDw6quvMmfOHKKjo7n88st58cUX8fHxwWKxsHr1alavXk2PHj3o378/zz33HEFBQZSXl+Pu7m7yMxIREVP5R1R+zdkPhlHnh103vAsAb287RFmFLmXU0pk6NHHOnDn885//PG2b5ORk+vbt67h96NAhPv30U95+++1aH7Nq1SpOnDjBtGnT6pwlKSmJyy67jPnz53PRRRedtu3cuXOZPXu243ZeXl6LKsbaB3UiyWMQA0oSOfjdCjp3f8zsSCKtWkVFRY37BwwYwMaNG2u8r0uXLqxevbra/u3bt9OtW7fGjCciIs7GPxywQFkBFB4F77pdX/LCfiEE+riRdaKEL3dnMr5/xzM/SExjaiF27733Mn369NO26d69e5Xbr7zyCh06dODSSy+t9TEvvfQSl1xyCSEhIXXKsWvXLi688EJuueUWHnrooTO2d3d3b/GfWBd0vxiSE2m/fx2gQkzEGbzyyis8/fTTTr96q4iInCUXd2jXCU6kVfaK1bEQc3OxctXQcJZ8vY8VcQdUiLVwphZiQUFBBAUF1bm9YRi88sor3HDDDbi6utbYJiUlhQ0bNrBmzZo6HXPnzp1ccMEFTJs2jSeeeKLOWVq6HqOvwb5rIX3K95B+8Gc6hjfO9Y9EpOnMmDHjtHNfRUSkDWnftbIQO54KYUPr/LBrhlUWYl//lMWhnELC2mvRjpbKqeaIffnll6SkpHDTTTfV2mbZsmV06tSJCRMmVLtv1apVVYY5JiUlMXbsWC666CJmz57tWJI6KyurSfI3p8DQruxxiwQg9Zu3TE4jIiIiIvXSgJUTASICvRnVowOGUTlXTFoupyrEXn75ZUaNGlWlmDqV3W5n+fLlTJ8+HZvNVu3+3Nxc9uzZ47j97rvvkpWVxeuvv06nTp0c27Bhw5rsOTSn3G6VxahvyicmJxERERGRemnAyoknXXty0Y6tBynXoh0tllMVYm+++SbfffddrfdbrVYOHjxY6xDD6dOnY5yy8swjjzxS47LUqampjR3dFBHnXgNA39KdZKdrGVMRERERp9HAHjGAi/qHEODtRnpeMV/tcf6RXq2VUxViUj8du/TiJ5feWC0G+zZqeKKIiIiI0ziLHjF3FxtXDukMwIo4fRjfUqkQa+WOdfkjAN77PjI5iYiIiIjU2ckeseMHwV7zZVJO55pfhydu2JPJkdyixkwmjUSFWCsXfnJ4YvGP5GQdMTmNiIiIiNSJbyhYXcFeBifq/zdcjyAfYrsFYDfg7a1atKMlUiHWynXu3p99tu64WOzs3bjS7DgiIiIiUhdWG/iFVX7fgHli8NuiHSu3HqDCbpyhtTQ3FWJtQGbYRQC4711rchIREWkuhYWFdO3alfvuu8/sKCLSUGcxTwzgjwM64ufpSlpuMRv3atGOlkaFWBsQOupaAPoVJZCbk21yGhERaQ5PPPEEI0aMMDuGiJyNs1g5EcDD1caVQyp71VZs0aIdLY0KsTaga59oUq3huFkq+Gnj22bHERGRJrZ37152797NhAkTzI4iImfjLHvEAK4dHg7A+t2ZZOYVN0YqaSQqxNqI9NDK4Ykuez40OYmIczr//POxWCzVthtuuMHsaOJkNm7cyKRJkwgNDcVisbB69epqbRYvXkxERAQeHh7ExsYSFxdXr3Pcd999LFy4sJESi4hpzrJHDKBXSDtiuranwm7wTrwW7WhJVIi1EcEjpgIQWbCV/Lwck9OIOBfDMPjhhx94+umnOXLkSJXthRdeMDueOJmCggKioqJYvHhxjfevXLmS2bNnM3/+fBISEoiKimL8+PFkZmY62kRHRzNgwIBqW1paGh988AG9e/emd+/ezfWURKSptI+o/HoWPWLw26IdK+IOYNeiHS2Gi9kBpHl0ixzGIUsnwjhC0jfvMXTiTWZHEnEae/fu5cSJE5x33nl07NjR7Dji5CZMmHDaIYPPPvssN998MzNmzABgyZIlfPTRRyxbtow5c+YAkJiYWOvjN2/ezFtvvcU777xDfn4+ZWVl+Pr6Mm/evFofU1JSQklJieN2Xl5ePZ+ViDSJkz1ieWlQXgIu7g06zMRBnXj0w50cyini25+zOa93UCOGlIZSj1gbYbFaOdjxD5XfJ68xOY3IbwoKCigoKMAwfvuErrS0lIKCgip/GJ7a1m63O/aVlZVRUFBAcXFxndo2RHx8PC4uLgwaNKhBjxepq9LSUuLj4xk3bpxjn9VqZdy4cWzatKlOx1i4cCEHDx4kNTWVp59+mptvvvm0RdjJx/j5+Tm28PDws3oeItJIvAPB1QswILfhwwo9XG1cMbgzAG9t1aIdLYUKsTYkcPgUAPqe2ExRwQmT04hU8vHxwcfHh+zs31b0fOqpp/Dx8WHmzJlV2gYHB+Pj48OBA7+9iSxevBgfHx9uvPHGKm0jIiLw8fEhOTnZsW/58uUNypiQkEBFRQUdOnRw5PXx8eHWW28FwMXFhejoaMdWVFQEwKFDh5g8eTI9evQgJiaGKVOmkJGRAUBgYGCDskjrlp2dTUVFBSEhIVX2h4SEkJ6e3mTnnTt3Lrm5uY7t4MGDTXYuEakHi+WUeWKpZ3Woa2Mrhyd+tjODrBMlZ2gtzUFDE9uQnlHncuSDIDpZskj4djVDxv/Z7EgiTiEhIYFrr72WRx99tMr+gIAAAPz9/asNFTMMg8suu4w77riD999/H4BvvvmGrKysan9kizSV6dOn16mdu7s77u4NG/IkIk2sfVfISj7rQqxvR18Gd/HnhwPHeTf+ELeP6dE4+aTBVIi1IRarlf0h4+iUsQL7ztWgQkxagPz8fAC8vLwc++6//35mzZqFi0vVl6iTixV4eno69t15553cfPPN2Gy2Km1TU1Orta3rH6W/l5CQwBNPPEHPnj3r/Jj169dX66kbPXp0g84vbUdgYCA2m83Rc3pSRkaG5ieKtFX+Z7+E/UnXDuvCDweO89bWA9x6XnesVstZH1MaTkMT2xj/oVcC0Cf3O0qKC01OIwLe3t54e3tjsfz2ZuDm5oa3t3e1T+hPtrVaf3vpcnV1xdvbGw8Pjzq1ra9ffvmF48ePExUVVWub48ePO4Yl3nRT5UI4u3btYsiQIfU+n7Rtbm5uDB06lPXr1zv22e121q9fz8iRI01MJiKmaX/2S9ifdElUJ3zcXdh/tJDNvxw96+PJ2VGPWBvTe+gFZH4cQLDlGNu/X0PUBdeYHUmkRYuPjwdqnqMTHByM1WqtcWiiSG3y8/P5+eefHbdTUlJITEwkICCALl26MHv2bKZNm0ZMTAzDhw9n0aJFFBQUOFZRFJE2phF7xLzcXLh8cCivbz7Am3EHGNVT85XNpEKsjbHabKQEjiU4+z1KfvwAVIiJnFZCQgIAvXr1qrLf3d2dvLw83Nzcanxcv379HHPDRE61bds2xo4d67g9e/ZsAKZNm8by5cuZOnUqWVlZzJs3j/T0dKKjo1m3bp3mFoq0VY3YIwZwzbAuvL75AJ/uTOdofgkdfDQ/1CwamtgG+QyuHJ7Y+/hGykq1ao7I6SxcuBDDMKptxcXFtRZhAOPGjSMvL6/KSo3ffvstSUlJzZBaWrIxY8bU+Dt16u/KzJkz2b9/PyUlJWzZsoXY2FjzAouIuU72iBUdg5KzX/V6QGc/BoX5UVZh8H7C4bM+njScCrE2qG/seI7hiz/57N78sdlxRFoli8XC6tWrWb16NT169KB///4899xzBAXpIpoiIlIPHr7g2b7y+0bqFbt2eOVS9iviDlS5jqc0LxVibZDNxYW9AWMAKNy+ytwwIq3AqddAO1WXLl1YvXo1+/btY+fOnaxcudIxvKy2x4iIiFTTiPPEACZFheLlZuOX7AK2pBxrlGNK/akQa6O8oq4AoOfRr6goLzc5jYiIiIjUqpHnifm4u3BZdCgAb8UdaJRjSv2pEGuj+o6cSC7edCCX3XGfmR1HRERERGrTyD1i8NvwxI+T0skpKG2040rdqRBro1zd3PnJ/zwA8n94z+Q0IiIiIlKrRu4RAxjY2Y/+ob6Ultt5/wct2mEGFWJtmOvAywHolvUl9ooKc8OIiIiISM38Iyq/NmKPmMVi4Zpfe8Xe0qIdplAh1ob1O+dS8g1PgjnGTwkbzI4jIiIiIjU5tUesEQumy6JD8XS1sTczn1e/T+VIblGjHbshjuQW8f2+7DaTQ4VYG+bu4cVuv3MAOL7tXZPTSFuhT9zaJrvdbnYEERHn5Rde+bWsAAqPNtphfT1c6R/qC8AjH+7inH98ycqt5izesXLrAc75x5dct3RLm8nh0mRHFqdg638pbPqCrhlfYNjtWKyqzaVpuLq6YrFYyMrKIigoCIvFYnYkaQaGYVBaWkpWVhZWq/W0F8EWEZFauHpAu05w4khlr5h3YKMc9khuEfEHchy37QbMeW8HG3Zn4elma5Rz1EVRaQXrdqa3yBwPvp/Eeb2D6OTn2ejnUyHWxvU9dzKF399PJ0sWW959lohRVxAS1sPsWNIK2Ww2wsLCOHToEKmpqWbHkWbm5eVFly5dsOrDHhGRhvHvWlmIHU+FsKGNcsiU7IJqIx0NqFKMmKWl5KgwDFKzC1WISePz9G7HXpdwelX8TOyuv1Ox83HiBj3C8CtnmR1NWiEfHx969epFWVmZ2VGkGdlsNlxcXNQLKiJyNtp3hYObG3XlxG6B3lgtlT0/J1kscNfYnvh6ujbaec4kr6iM5778mVNrwpaSw2axEBHo1STnUyHWxmUc2keP8n3w699HNovBkB8fJSN2knrGpEnYbDZstuYbZiAiItIqNMG1xDr5ebJw8kAefD+JCsPAZrGwYPIApg7r0mjnqKvO7T1bbI6m6A0DFWJtXtb+XYRYqvZJu1jsZO/frUJMREREpKVogmuJAUwd1oXzegeRml1IRKBXkxUdylGdCrE2LqhrJBWGBdspxZjdsBDYta+JqURERESkiiboETupk5+naYVPW86hWdNtXEhYD+IHPUK58duvgtVisH/LGhNTmS/j0D6SvvuQjEP7zI4iIiIi8luP2PGDYK8wN4s0ChViwvArZ3H05m3s/MObbA68CoAhPz5G4vq3TE5mjrj3FhG4dCgDPr+ewKVDiXtvkdmRREREpK3z7QxWF7CXVa6eKE5PhZgAlT1j/c+ZSOwdS9nqPwEXi50+G+9iz7YvzY7WrDIO7WPoj484hmo6Fi9Rz5iIiIiYyWoDv7DK7xt5npiYQ3PEpAqL1Ur0Ha/y47MTGVS8leC10zjot5bwXlFmR2sW6Xu21bh4SVbqLi1eItJKdOvWrUFL6c+aNYu77767CRKJiNSRf1fISf11ntg5ZqeRs6RCTKpxdXOnx53vsnfRhfSq+JmiN6eQfet6AjuGmx2tSZ3IPYbn90/VeJ/t6wVkRESqGBNpBZYvX96gx0VERDRqDhGRemvfFVJQj1groUJMauTdzp/2N6/m8IsX0tnI4Oell+Fxz3p8fNubHa1J5OZkk774YvpU7KXQcMOdMmwWA7thoQwb/cp2kfvSaOKHP8HQi2eYHVdEzsL5559vdgQRkYZpwpUTpflpjpjUKrBjOPY/vUcOvvSs2McvL1xJWWmJ2bEaXe7RDDKfH0+f8j0cx4e0yavIvjmenX94k6yb48m8/kv2uvTCjwKGxs1i67+u5kTuMbNji4iISFvTPqLyq3rEWgUVYnJa4T0HknXp/yg03BlUHE/i4j9j2O1mx2o0OVlHyH7hj/Sq+JkcfDl65Xv0jDrXsXhJSFgPwntFEfHX79gU9hcqDAvDcj/lxKJYkrd8anZ8EWkkQUFBBAcH17iFh4dz3nnnsWHDBrNjikhbpx6xVkWFmJxR7yFj2Hv+/1FuWBmW+ymbX/5/ZkdqFEczDnH8P3+kR8UvZONP7tRV9Bg4osa2rm7ujLzpX+y9+G3SLCGEGpn0/ngqm5bOapW9hCJtTVZWFpmZmTVuqampPPfcc1qoQ0TMd/JaYnlpUK6/P5ydCjGpk6gLriFh0HwARh5ezpaV/zQ50dnJTttP/ot/pJs9lSzaU3DtB0T0iznj4/rGXkS7WZvZ6j8Bm8Vg5OFXSH3yHA78lNj0oUWkyWRlZbFr165q+3ft2sWxY8eIiori3nvvNSGZiMgpvIPA1QswIPeQ2WnkLKkQkzobfuUsNnW5FYBhuxaS8On/TE7UMJmHUyha+ke62g+SQQeKr/+Qrn2i6/z4dn4BDJv1Fgmxi8jFm17lewl64w9sefupVjVsU6QtmTlzJjk5OdX25+TkOHrCpk+f3sypRER+x2IB/y6V3+ekmhpFzp4KMamXEdP/wZaAS7FaDPp///+cbp5U+oG9lL70R8KNNNIJovyGjwjvObBBxxoyYQYlN3/HDvfBeFpKid31ONufmsDRDH1CJeJsUlJSOOec6tfkOeecc0hKSjIhkYhILTRPrNVQISb1YrFaGXr7yyR6jcTdUkboJzPYnxxvdqw6SUvdg/2Viwkz0kmzhGCf/hGdu/c7q2MGd+5G/7+uZ3Pv+ygxXIku2gz/GUXi+rcaKbWINIeaesNOKioqasYkIiJncHKemFZOdHoqxKTeXFzd6HPnO+xx6YMfBbivvJrMwylmxzqtw7/sxLp8IqFGJocsnbD+5WNCI/o0yrGtNhsjrnuYtKs/IsUaQQdyif7mVrY8N42ighONcg4RaVqDBg2q8ULPr732GgMHNqzXXESkSahHrNWwGIZhmB3C2eXl5eHn50dubi6+vr5mx2k2OVlHyH/hAsKNNH6xRhB495f4+ncwO1Y1B/dux/2NywnmGAesnfG48SOCO3drknMVFxWQuPxeRmSsAOCAtTMlk5bQa/B5TXI+kdagJbyGHj58mMsuu4x27doxZMgQABISEjhx4gSrV68mLCzMlFxmaQk/ExGpRfKHsPJ6CB0Ct+iyGi1RXV9D1SMmDdY+qBO2G1aRjT/d7akcfOEKSooLzY5Vxf7dCXi8cRnBHCPVGo7XLeuarAgD8PD0ZsTtS0i68DUyCaCL/TARqy9n8/IHqSgvb7LzisjZ6dy5M9u2beOhhx4iIiKCiIgI/va3v7Ft27Y2V4SJSAunHrFWQ4WYnJXQbn05fsUbFBge9C/dTtLiP2GvqDA7FgApu7bi89blBJFDijWCdreuI7Bjl2Y594DRl+F+12YSfM7H1VLBiNTF/PTP0aSl7G6W84tI/Xz88cd8/PHHlJSU0KNHD3r27Imfnx+FhS3rw6X6SElJYezYsURGRjJw4EAKCgrMjiQijeHkHLHCo1CSb24WOSsuZgcQ59cz6lx2HP8Pfb+8iaEnvmTzf+9kxO1LTM20b8dmAt6bQnvy2GfrTsBtH9M+qFOzZvDrEMLg2avZ+uES+iU8Rr+yXeQvH8PWIfOImXQbFqs+BxFpKd55551q+44dO0ZSUhL//e9/ufDCC01IdXamT5/O448/zujRozl27Bju7u5mRxKRxuDhBx7+UHy8slcspL/ZiaSBVIhJoxh4/mS2Hk9n2A9zGZGxgs1vhDLiT/NMyfLz9m8JXDUVf/LZ69KL4Ns/wq9DiClZLFYrwy67g7SoCzn4xgz6le1k2A9zid/7KT1n/Ne0XCJS1SuvvFLj/kOHDjF58mTi4uKaOdHZ2blzJ66urowePRqAgIAAkxOJSKNq3xWOHK9cOVGFmNPSR/LSaIZddgebulde+HTE3meI/+ilZs/wU8JXBK+6Gn/y2ePSh+A717WIYic0og+9H9jI5og7KTNsDM3/ipLnRpD0zQdmRxOR0wgLC6OsrKzRj7tx40YmTZpEaGgoFouF1atXV2uzePFiIiIi8PDwIDY2tl7F4N69e/Hx8WHSpEkMGTKEBQsWNGJ6ETGd5om1CirEpFGNuP5RtgReCcDAuAdI+u7DZjv37q1f0OmDa/ClgGTXSELvWodf+8BmO/+Z2FxcGDF9AamXr+agJZRgjjFg/Q1s/s9tFBdp7oZIS7Rp06YmWTWwoKCAqKgoFi9eXOP9K1euZPbs2cyfP5+EhASioqIYP348mZmZjjbR0dEMGDCg2paWlkZ5eTnffPMNL7zwAps2beLzzz/n888/rzVPSUkJeXl5VTYRacF0LbFWQUMTpVFZrFZibvsvCf/KZEjBN3T97GZS/D+gW//YJj3vrs3r6PrJNLwtxex0G0jEXWvxbuffpOdsqF6Dz6Ow12a2vHI3sUdXMyJjBSlPbYIr/9vk/04iUrNhw4ZhsViq7Dt27Bjt27fn1VdfbfTzTZgwgQkTJtR6/7PPPsvNN9/MjBkzAFiyZAkfffQRy5YtY86cOQAkJibW+vjOnTsTExNDeHg4ABdffDGJiYn84Q9/qLH9woULefTRRxv4bESk2alHrFVwih6xr776CovFUuO2detWAB555JEa7/f29q7TOY4ePUpYWBgWi4Xjx4834bNp/WwuLkTe+RbJrv1pZynC+51rSD+wt8nOl/Tdh0R8cgPelmKS3KPpdvdHLbYIO8nLx4/Yu14lcfSLHMOXbvZUQt+eyOY3/95iVp0UaUveffdd3nnnHcf27rvvsn37drZu3cqOHTuaNUtpaSnx8fGMGzfOsc9qtTJu3Dg2bdpUp2MMGzaMzMxMcnJysNvtbNy4kX79+tXafu7cueTm5jq2gwcPnvXzEJEm1D6i8qt6xJyaUxRio0aN4siRI1W2m266iW7duhETEwPAfffdV61NZGQkU6ZMqdM5brzxRgYNGtSUT6NN8fDyIfT2D0i1hhPMMYqXX0Hu0YxGP8+Ojavo+dkMvCwl/OgRQ897PsLLx6/Rz9NUoi+8Bvtt35PoOQJ3SxkjfnqanU9eSObhFLOjibQpXbt2rbJ16dLF8UHe/fff36xZsrOzqaioICSk6vzWkJAQ0tPT63QMFxcXFixYwHnnncegQYPo1asXl1xySa3t3d3d8fX1rbKJSAt2ao+YYZibRRrMKQoxNzc3Onbs6Ng6dOjABx98wIwZMxxDSXx8fKq0ycjIYNeuXdx4441nPP5//vMfjh8/zn333dfUT6VN8QsIwmP6KjIJIMJ+kLQlVzTqXKjtG96h9/qb8bCUsd0zlt73rMHDy6fRjt9cAjuGE3X/J2zp/zBFhhsDS37Afek5JHzyChmH9pH03YdkHNpndkyRNstw0j9yJkyYwI4dO0hKSuLZZ581O46INCb/X6+LWpoPhcfMzSIN5hSF2O+tWbOGo0ePOsbO1+Sll16id+/ejqV7a7Nr1y4ee+wxXnvtNay6rlOj69ilF/lXvcUJw5N+ZTvZ9fw1VJSXn/VxE79YQb+vbsPdUsYPXufQb9YaPDzrNgy1JbJYrcROuY+sP33OXpde+FHAkC2zCF46hAGfX0/g0qHEvbfI7JgibdLv5441tcDAQGw2GxkZVUcRZGRk0LFjx2bNIiItlKsH+Pz6enA81dQo0nBOWXm8/PLLjB8/nrCwsBrvLy4u5o033jhjb1hJSQnXXnstTz31FF26dKnz+bW6VP10HxDL/ouWUmq4MKRgI9uW3IJhtzf4eD98+ir9v7kTN0s5CT7nMWDWKtzcPRoxsXm69I4m4q/fsTnwKgwDTv79Z7MYDPnxUfWMiTSRoKAggoODq21BQUEcOXKkWbO4ubkxdOhQ1q9f79hnt9tZv349I0eObNYsItKCaeVEp2dqITZnzpxaF+E4ue3evbvKYw4dOsSnn3562iJr1apVnDhxgmnTpp32/HPnzqVfv35cf/319cq9cOFC/Pz8HNvJVamkdgPOmcSO4f8EIDb7PTa/Pr9Bx4n/+GUGfj8LV0sF29pdyKB73sPVzb0xo5rO1c0dn8GX8/sP4V0sdrL37675QSJyVrKyssjMzKy2ZWVlNcl1xPLz80lMTHSsfJiSkkJiYiIHDhwAYPbs2SxdupRXX32V5ORkbr/9dgoKCk47EkRE2hitnOj0LIaJg9+zsrI4evToadt0794dNzc3x+2///3vPPfccxw+fBhXV9caH3PhhRfi6+vLqlWrTnvs6OhoduzY4Rh2YhgGdrsdm83G3/72t1qX8i0pKaGkpMRxOy8vj/DwcHJzczXB+Qw2v/EYI/Y+A8DWwQsZdtkddX7stjVLGBw/B5vFYKvfRQy5awU2l9Z5BYaMQ/sIXDoUm+W3/54VhoXsm+MJCethYjKRxpeXl4efn1+zvoZOnz6dF154AS8vr2Y53+999dVXjB07ttr+adOmsXz5cgCef/55nnrqKdLT04mOjub//u//iI1tnktcmPEzEZF6+vJx2PgUDJ0BkxaZnUZOUdfXUFMLsfoyDIMePXowefJknn766RrbpKSk0KNHD9asWXPaFaIA9u3bR1FRkeP21q1b+ctf/sL3339Pjx49CA4OrlMuvWHVz+b/3MaIjBWUGTZ2X/gyA8+74oyP2br6eYb+8BBWi0Gc/8UMnfm/VluEnRT33iKG/PgoLpbKYZyplnC6PvwjFs1llFbGjNdQm83GkSNHHK/zt99+OwsXLsTf39/Rpry8HJdW/jpTG72viTiBhP/BmpnQ4wL48+k7H6R51fU1tM5/0U2fPp3CwsJGCddQX375JSkpKdx00021tlm2bBmdOnWq8UKZq1atom/fvo7bPXr0YMCAAY6tW7duAPTr16/ORZjU3/BbFhPf7gJcLRV0X38bP2//9rTt4977l6MI29LhMmLuer3VF2EAw6+cxdGbt7FtyJOUGC5EGAdJ/OJNs2OJtAq//wzyjTfe4Nix31Yey8jIUAEiIi2b5og5vToXYv/73//Iz8933L799turXfi4vBFWwzudl19+mVGjRlUppk5lt9tZvnw506dPx2azVbs/NzeXPXv2NGlGOTOrzcaAO99gp1sU3pZi/Ff9ibSUmuc+bXn7SYbveKSyCAu8kuF3Lsdaw8+2tQoJ60HMpbeSEFY5jzF402ONegkAEalU0+CQ4uJiE5KIiNTRyTliuQfhLBZBE/PUuRBrCZ8evvnmm3z33Xe13m+1Wjl48CBPPPFEjfdPnz79tNeDGTNmDIZhVBmaIk3D3cOL8DtWkWKNIJDjVLx2BTlZVVcm27ziCWJ3Vf4sN4dcw/A7Xmqzw/IGXfMomQTQ2cjgh7dr/v0WkcbV3MvWi4jUi29nsNigohRONO/qrtI4GvxXrT49lLPl698B7xtXk04Q4UYamS9ezsGffiTpuw/ZtPT/MWLPkwBs6nQ9sbf+p80WYQDe7fw5MOQBAKJ+eYnMwykmJxJxfm+++SYJCQlNsiqiiEiTs7mA36+XctLKiU6pUf+y1aeHUl/BnbtRcs3b5OJNn/LdhL0xmgGfX8/Iw8sA2BT2F0bc/FybLsJOGnrJLex2jcTLUsKBlfebHUfEqY0ePZr58+cTExODj48PhYWFzJ8/nyVLlrB58+YqQ/FFRFoszRNzavX661afHkpT6Np3CPtGLKxyAWMAu2Gh+x9nqgj7lcVqxTbxSeyGhZi8z9kd97nZkUSc1tdff+2YN/zqq69y7733cuTIER588EFGjRpF7969zY4oInJmupaYU6vz0nMnPz08ceIErq6ulJeXM3/+fM455xyio6MJCgpqypzSyrm1C6h2AWOrxSB7/25dN+sUvaJHE/f1xQzP+Qjbp3OwD72gTS1eItLYevXqRa9evbjmmmsc+1JSUti2bRs//PCDiclEROpAPWJOrc6F2Ndffw3A3r17iY+PJyEhgYSEBB588EGOHz+uYYlyVoK6RlJhWKpcwLjcsBLYteYVMtuyHtc8yYkXvqRXxc/ErX6O4VfOMjuSSKvSrVs3unXrxpQpU8yOIiJyev4RlV/VI+aU6n0xJn16KE0hJKwHcYMecVzAuNywkjBoPsPVG1ZNh5AwNve+nRF7n6XHjmfJu/DP+Pp3MDuWiIiINDf1iDk1i3G69dylTup69Ww5s4xD+8jev5vArn01JPE0SkuKOfLPoXS1H2JzyLWMuH2J2ZFEGkyvoS2PfiYiTuJEBjzTG7DAQ5ng4mZ2IqHur6FaBUFalJCwHvQ/Z6KKsDNwc/cg97zHABia/jb79ySaG0ikFfjpp58oLy83O4aISN35BIOLJ2BUXthZnIoKMREnNWjMlSR6jcTVUkHuqvsw7HazI4k4tX79+vHLL7+YHUNEpO4sFvDvUvm95ok5HRViIk4scPLTlBouDCreyvYNb5sdR8SpaaS+iDglzRNzWirERJxYWM8BxIdeB0Dgt/MpKS40OZGIiIg0K11LzGmpEBNxcgOvfYws2hNmpJPw9kKz44iIiEhzUo+Y01IhJuLkfHzbkxp9PwCD9v2X7DS9EIuIiLQZ6hFzWirERFqBoZNuY49LH7wtxaSsvN/sOCIiItJc1CPmtFSIibQCVpsNy8VPAjAs91P2bPvS5EQiIiLSLE72iBVmQ0m+uVmkXlSIibQSvYeMIc7/YgAs6x7AXlFhciIRERFpcp7+4OFX+f3xA6ZGkfpRISbSinS/5p/kG570Lv+JbWteMDuOiFN54IEH6NChg9kxRETqT/PEnJIKMZFWJLBjF5J63QpA9+1PcyL3mMmJRJzHwoULVYiJiHPSPDGnpEJMpJUZMmUuBy2hBHKcnSseMjuOiIiINDX1iDklFWIirYybuwfHRj8KwJAjb3Fw73aTE4mIiEiTah9R+VU9Yk5FhZhIKxR1wdVs9xyOm6WCY+/fZ3YcERERaUrqEXNKKsREWqmAyU9TatiIKopj+5dvmx1HpEWZPn06hYWFZscQEWkcp84RMwxzs0idqRATaaXCe0WR0OkaAAK+mU9pSbHJiURajv/973/k5/92vZ3bb7+d48ePV2lTXl7ezKlERBrIv0vl19ITUJRjbhapMxViIq1Y/2sfJxt/wo00Et5ZaHYckRbD+N0nxm+88QbHjv22ymhGRga+vr7NHUtEpGFcPcEnpPL7nFRTo0jdqRATacXa+QXwS1TlHLEBe18kO10XehSpye8LM4DiYvUii4gT0Twxp6NCTKSVi7n0Dn5y6Y2PpYiUt/5qdhwRp2GxWMyOICJSd7qWmNNRISbSylltNow//hOAYcc/4aeEr8wNJNJCvPnmmyQkJFBWVmZ2lCbxr3/9i/79+xMZGcndd99dY6+fiLQi6hFzOirERNqAPjEXsNVvPADGx3/FXlFhciIRc40ePZr58+cTExODj48PhYWFzJ8/nyVLlrB58+YqC3k4o6ysLJ5//nni4+PZsWMH8fHxbN682exYItKU1CPmdFzMDiAizaPb1KcoePFr+pTvYeuHSxh2+Z1mRxIxzddffw3A3r17iY+PJyEhgYSEBB588EGOHz/eKoYllpeXO+a5lZWVERwcbHIiEWlS6hFzOuoRE2kjAkO78mOPWwDolvgk+Xla3lakV69eXHPNNTz55JN88cUXHDt2jJ9//pm33nqLBx54oMnOu3HjRiZNmkRoaCgWi4XVq1dXa7N48WIiIiLw8PAgNjaWuLi4Oh8/KCiI++67jy5duhAaGsq4cePo0aNHIz4DEWlxTvaIHT8Adru5WaROVIiJtCFDrp7LIUsnAjnOjhXzzI4j0iJ1796dKVOmsGDBgiY7R0FBAVFRUSxevLjG+1euXMns2bOZP38+CQkJREVFMX78eDIzMx1toqOjGTBgQLUtLS2NnJwc1q5dS2pqKocPH+b7779n48aNTfZ8RKQF8A0Diw0qSiE/3ew0UgcamijShrh7eJF9znzCvr2NoWlvcvDnWwjvOdDsWCLNqlu3bg0aejhr1izuvvvuRskwYcIEJkyYUOv9zz77LDfffDMzZswAYMmSJXz00UcsW7aMOXPmAJCYmFjr49955x169uxJQEAAABMnTmTz5s2cd955NbYvKSmhpKTEcTsvL6++T0lEzGZzAb/OlT1iOfvBN9TsRHIGKsRE2pioC6by47aXGFS8jaPv3Uf4A5+aHUmkWS1fvrxBj4uIiGjUHLUpLS0lPj6euXPnOvZZrVbGjRvHpk2b6nSM8PBwvv/+e4qLi3F1deWrr77illtuqbX9woULefTRR886u4iYzL9rZSF2fD90HWl2GjkDFWIibYzFasXvimcoe/MCoos28+OGdxk09iqzY4k0m/PPP9/sCKeVnZ1NRUUFISEhVfaHhISwe/fuOh1jxIgRXHzxxQwePBir1cqFF17IpZdeWmv7uXPnMnv2bMftvLw8wsPDG/YERMQ87btC6jdaOdFJqBATaYO69olmc8erGZGxAv9v5lE66hLc3D3MjiUijeiJJ57giSeeqFNbd3d33N3dmziRiDQ5/4jKr1o50SlosQ6RNiry2ic4ih9d7IdJePefZscRkV8FBgZis9nIyMiosj8jI4OOHTualEpEnIKuJeZUVIiJtFG+/h3YN7ByKFL/n/5DdvpBkxOJCICbmxtDhw5l/fr1jn12u53169czcqTmfIjIaehaYk5FhZhIGxZz+V3stfWknaWIX1Y23TWTRKSq/Px8EhMTHSsfpqSkkJiYyIEDBwCYPXs2S5cu5dVXXyU5OZnbb7+dgoICxyqKIiI1OtkjlncYKsrMzSJnpDliIm2Y1WajYvw/4OOriDn2MXsTv6FX9GizY4m0etu2bWPs2LGO2ycXypg2bRrLly9n6tSpZGVlMW/ePNLT04mOjmbdunXVFvAQEanCJwRcPKC8GHIPQkB3sxPJaagQE2nj+g7/A9u+HUdM3hdUrL0fY9D3WKzqLBdpSmPGjMEwjNO2mTlzJjNnzmymRCLSKlgs4N8Fsn+qnCemQqxF019bIkKXqU9TaLjTtzyZ+LUvmh1HREREGkrzxJyGCjERIbhzN37sdhMAXRKepODEcXMDiYiISMNo5USnoUJMRACInvo3DltCCOYYP66YZ3YcERERaQj1iDkNFWIiAoCHpzeZIysLsKGH3+DwLztNTiQiIiL1ph4xp6FCTEQcosddxw73IbhZysl89z6z44iIiEh9qUfMaagQExEHi9WK7+VPU25YGVz4PTs2rjI7koiIiNTHyR6xgiwoLTA3i5yWCjERqaJrv6FsC7kKAN+vHqastMTkRCIiIlJnnu3B3a/y++MHzM0ip6VCTESq6XftQnLwpav9IPHvPml2HBEREamP9l0qv2qeWIumQkxEqvFrH8hP/WcBEPnTCxzLPGxuIBEREak7zRNzCirERKRGMVfcw8+2HvhSyM9vPWB2HBEREamr9hGVX9Uj1qKpEBORGtlcXCj9wwIAYo5+yJZ3nibj0D6TU4mIiMgZqUfMKagQE5FaRY74I3tcemO1QOzOvxO4dChx7y0yO5aIiIicjq4l5hRUiIlIrTIO7aNn2V7HbZvFYMiPj6pnTEREpCU7tUfMMMzNIrVSISYitcravwubpeoLuIvFTvb+3SYlEhERkTPy/3XVxJI8KMoxN4vUSoWYiNQqqGskFYal2v68/QkmpBEREZE6cfMC7+DK7zVPrMVyikLsq6++wmKx1Lht3boVgEceeaTG+729vc94/OXLlzNo0CA8PDwIDg7mzjvvbOqnJOIUQsJ6ED/oEcqNypcK+6+dYzF7/sX2De+YmExEREROS/PEWjwXswPUxahRozhy5EiVfQ8//DDr168nJiYGgPvuu4/bbrutSpsLL7yQYcOGnfbYzz77LM888wxPPfUUsbGxFBQUkJqa2qj5RZzZ8CtnkRE7iez9uwkI68Xh9+YSk/cFvb66kz3tOtAn5gKzI4qIiMjv+XeFQ1vVI9aCOUUh5ubmRseOHR23y8rK+OCDD7jrrruwWCqHTfn4+ODj4+Nos337dnbt2sWSJUtqPW5OTg4PPfQQH374IRdeeKFj/6BBg5rgWYg4r5CwHoSE9QCgw51v8OO/JjKoeBvBa29gf7u1dO0TbW5AERERqUo9Yi2eUwxN/L01a9Zw9OhRZsyYUWubl156id69ezN69Oha23z++efY7XYOHz5Mv379CAsL4+qrr+bgwYNNEVukVXBz96DHne/xk0tv2nMC9xVXkXk4xexYIiIicipdS6zFc8pC7OWXX2b8+PGEhYXVeH9xcTFvvPEGN95442mP88svv2C321mwYAGLFi3i3Xff5dixY/zhD3+gtLS01seVlJSQl5dXZRNpS7zb+RN4y2oOWkLpSBYFL19G7rEss2OJiIjISeoRa/FMLcTmzJlT6yIcJ7fdu6suk33o0CE+/fTT0xZZq1at4sSJE0ybNu2057fb7ZSVlfF///d/jB8/nhEjRrBixQr27t3Lhg0ban3cwoUL8fPzc2zh4eH1e+IirUBAcGds01aTRXu62feT9p/LKC7MNzuWiIiIwCk9YgfAbjc3i9TI1ELs3nvvJTk5+bRb9+7dqzzmlVdeoUOHDlx66aW1Hvell17ikksuISQk5LTn79SpEwCRkZGOfUFBQQQGBnLgwIFaHzd37lxyc3Mdm4YySlsVGtGHE1etJA8v+pXtJPn5KZSX1d6bLCIiIs3ELwwsVqgogfwMs9NIDUxdrCMoKIigoKA6tzcMg1deeYUbbrgBV1fXGtukpKSwYcMG1qxZc8bjnXPOOQDs2bPHMczx2LFjZGdn07Vr11of5+7ujru7e51zi7Rm3QfEsuvEMnqs+zODC78n7oUZDLvrf1isTjnyWUREpHWwuYJvGOQeqJwn5tvJ7ETyO071l9KXX35JSkoKN910U61tli1bRqdOnZgwYUK1+1atWkXfvn0dt3v37s1ll13GPffcw/fff09SUhLTpk2jb9++jB07tkmeg0hrFDlyArvOWUSFYWF4zlo2L7vX7EgiIiKieWItmlMVYi+//DKjRo2qUkydym63s3z5cqZPn47NZqt2f25uLnv27Kmy77XXXiM2NpaJEydy/vnn4+rqyrp162rtcRORmg2+6Hq2DXgYgJGHlrHlrYUmJxIREWnjtHJii2YxDMMwO4Szy8vLw8/Pj9zcXHx9fc2OI2KqTa88wMj9S7AbFn6IfYahF59+9VIRvYa2PPqZiLQSXz8JG56A6Ovh8sVmp2kz6voa6lQ9YiLS8o2YtpAtgZOxWgwGbrmfpG8+MDuSiIhI26QesRZNhZiINCqL1UrMbUtJ8DkPN0sF3b64hZ+3f2t2LBERkbZHc8RaNBViItLobC4u9J+5kp1uUXhbimm/6joO/ZxkdiwREZG25WSPWN4hqCgzN4tUo0JMRJqEu4cX4XesYp+tOx3IxfLGZLLTa78+n4jU3xVXXEH79u256qqrqt23du1a+vTpQ69evXjppZdMSCcipvMJAZs7GHbIPWR2GvkdFWIi0mR8/Tvgd/MHHLaE0NnIIHfpZZzIPWZ2LJFW45577uG1116rtr+8vJzZs2fz5Zdf8sMPP/DUU09x9OhRExKKiKmsVvDvUvm95om1OCrERKRJBXbsgvGn9zmKHz0qfuHA4sspKS40O5ZIqzBmzBjatWtXbX9cXBz9+/enc+fO+Pj4MGHCBD777DMTEoqI6TRPrMVSISYiTS6s5wByrniTAsOD/qXb2fn8VCrKy82OJdKkNm7cyKRJkwgNDcVisbB69epqbRYvXkxERAQeHh7ExsYSFxfXKOdOS0ujc+fOjtudO3fm8OHDjXJsEXEyWjmxxVIhJiLNomfUuaSM+y+lho0h+RvZtuRmDLvd7FgiTaagoICoqCgWL6752j0rV65k9uzZzJ8/n4SEBKKiohg/fjyZmZmONtHR0QwYMKDalpaW1lxPQ0ScnXrEWiwXswOISNsxYPRlxJ94isFb7iU2+302vRrMyBn/NDuWSJOYMGECEyZMqPX+Z599lptvvpkZM2YAsGTJEj766COWLVvGnDlzAEhMTGzQuUNDQ6v0gB0+fJjhw4fX2r6kpISSkhLH7by8vAadV0RaIPWItVjqERORZjX04huJ6/tXAEbuX0Lcu8+anEik+ZWWlhIfH8+4ceMc+6xWK+PGjWPTpk1nffzhw4eTlJTE4cOHyc/P55NPPmH8+PG1tl+4cCF+fn6OLTw8/KwziEgLoR6xFkuFmIg0uxHXPsimzpW9AEN3PMYPn71uciKR5pWdnU1FRQUhISFV9oeEhJCenl7n44wbN44pU6bw8ccfExYW5ijiXFxceOaZZxg7dizR0dHce++9dOjQodbjzJ07l9zcXMd28ODBhj0xEWl5TvaIFWRCqRbLakk0NFFETDHixmeJey6T4TkfEfndLHa160DkyNqHcYlIdV988UWt91166aVceumldTqOu7s77u7ujRVLRFoSz/bg7gsleXD8AAT3NTuR/Eo9YiJiCovVypA7lvOD1yjcLWWEffoXUnZuMTuWSLMIDAzEZrORkZFRZX9GRgYdO3Y0KZWItEoWi+aJtVAqxETENC6ubvSb+Q7Jrv3xpRCfd6aSlrrH7FgiTc7NzY2hQ4eyfv16xz673c769esZOXKkiclEpFXSPLEWSYWYiJjKw8uH0Ns/INXahSByqHj1cnKyjpgdS+Ss5efnk5iY6Fj5MCUlhcTERA4cOADA7NmzWbp0Ka+++irJycncfvvtFBQUOFZRFBFpNOoRa5E0R0xETOcXEETJjWtIX/oHwo00fnrxUtzu+Rzvdv5mRxNpsG3btjF27FjH7dmzZwMwbdo0li9fztSpU8nKymLevHmkp6cTHR3NunXrqi3gISJy1hw9YqmmxpCqLIZhGGaHcHZ5eXn4+fmRm5uLr6+v2XFEnNb+PYn4rZiIP/n86BFDv9kf4+qmBQRaO72Gtjz6mYi0MnvWwYqp0HEg3Pat2Wlavbq+hmpoooi0GF37RJNxyf8oNNwZVLyN7c//CXtFhdmxREREnJujR+yAuTmkChViItKi9Im5gL1jFlNuWInJ+5y4/95pdiQRERHn5t+l8mtJLhTlmJtFHFSIiUiLEzV2Cj8MfhyAERkr2Pz6fJMTiYiIODE3b/AOqvxeKye2GCrERKRFGnb5nWzuOQuAET8vYuvqxeYGEhERcWZaObHF0aqJItJijbj+UTb/J4MRGSsY/MNDbAU8g7oQ1DWSkLAeZscTadPsdjulpaVmx5Bm5Orqis1mMzuGNFT7rnB4m3rEWhAVYiLSog2/ZTHb/p1NTN7nxPzwIBYLVBgW4gY9wvArZ5kdT6RNKi0tJSUlBbvdbnYUaWb+/v507NgRi8VidhSpL/WItTgqxESkRbPabHS64u8Yyz/n5Pu+zWIw5MdHyYidpJ4xkWZmGAZHjhzBZrMRHh6O1apZDm2BYRgUFhaSmZkJQKdOnUxOJPXmWDlRhVhLoUJMRFq8nLSf6fy7D19dLHay9+9WISbSzMrLyyksLCQ0NBQvLy+z40gz8vT0BCAzM5Pg4GANU3Q26hFrcfQxloi0eEFdI6kwqlZidsNCYNe+JiUSabsqfr22n5ubm8lJxAwni++ysjKTk0i9newRO34ADMPcLAKoEBMRJxAS1oP4QY9Qbvz2klWGjQr9ISBiGs0Rapv0c3difuFgsUJ5MeRnmJ1GUCEmIk5i+JWzOHrzNpLGvc4eW2/cLeXkrLgF+6+fzouIiMhp2FzBt3Pl95on1iKoEBMRpxES1oMB507C50/LKTTc6V+6na3vPmV2LBEREeegeWItigoxEXE6nbv3Z0fkbAAG7nqWw7/sNDmRiIiIE9DKiS2KCjERcUrDrrqfnW5ReFlKyNUQRRGRBjl+/DgxMTFER0czYMAAli5danYkaUqOHrFUU2NIJRViIuKUrDYb7a9bSoHhQWRZEnErF5odSUTE6bRr146NGzeSmJjIli1bWLBgAUePHjU7ljQV9Yi1KCrERMRphUb0IWnA/QBE7fk3B/duNzmRiLR0c+bMwd3dneuuu87sKC2CzWZzLElfUlKCYRgYWtq89dIcsRZFhZiIOLXhV85mh/tgPC2lFLx9GxXl5WZHEpEWbO7cuTzzzDOsWLGCn3/+2ew4p7Vx40YmTZpEaGgoFouF1atX19hu8eLFRERE4OHhQWxsLHFxcfU6z/Hjx4mKiiIsLIz777+fwMDARkgvLdLJHrHcw1Ch90uzqRATEadmsVoJ+tNS8g1P+pbtYutbj5sdSURaMD8/P2688UasVis7duwwO85pFRQUEBUVxeLFi2tts3LlSmbPns38+fNJSEggKiqK8ePHk5mZ6Whzcv7X77e0tDQA/P392b59OykpKbz55ptkZOgaU62WT0ewuYNRAXmHzE7T5qkQExGn17FLL3YNmgPA4L3Ps39PormBRKROjuQW8f2+bI7kFjXrecvLy/Hy8iIpKalZz1tfEyZM4PHHH+eKK66otc2zzz7LzTffzIwZM4iMjGTJkiV4eXmxbNkyR5vExESSkpKqbaGhoVWOFRISQlRUFN98802TPScxmdUK/uGV32uemOlUiIlIqzDsirv50WMY7pYySt65mfKyUrMjibQJhmFQWFpe7+1/m1I55x9fct3SLZzzjy/536bUeh+joXOZHnroIfLz85utEFuwYAE+Pj6n3Q4cOFDv45aWlhIfH8+4ceMc+6xWK+PGjWPTpk11OkZGRgYnTpwAIDc3l40bN9KnT596ZxEnonliLYaL2QFERBqDxWol5PoXyXvpXHqX/8SmFX9n5A1/NzuWSKtXVFZB5LxPz+oYdgMe/mAnD39Qv2sC7npsPF5u9ftTJj4+niVLljBx4sRmK8Ruu+02rr766tO2+X3vVF1kZ2dTUVFBSEhIlf0hISHs3r27TsfYv38/t9xyi2ORjrvuuouBAwfWO4s4Ea2c2GKoEBORViMkrAdbo//GsMS/MXTfC6TsuoxukTFmxxKRFsJut3Prrbcyc+ZMYmNjuf766ykrK8PV1bXOx0hLS+P+++/njTfeqPNjAgICCAgIaEjkJjd8+HASExPNjiHNST1iLYYKMRFpVWIuvYPEPR8SXbSZ8vdupazn97i6uZsdS6TV8nS1seux8fV6THpuMeOe/Rr7KSMLrRb4Yvb5dPTzqNe56+O5554jOzubxx57jAMHDlBWVsbu3bvr1QMUGhparyIMKocmLliw4LRtdu3aRZcuXep13MDAQGw2W7XFNTIyMujYsWO9jiVtiHrEWgzNERORVsVitdL5zy+Size9Kn5m2xvzzY4k0qpZLBa83FzqtXUP8mHh5IHYLBYAbBYLCycPpHuQT72OY/n18XVx+PBhHn74YRYvXoy3tze9evXC3d3dMTwxNTWVqKgo/vSnP9GrVy9uv/12Vq9eTWxsLAMGDGDv3r2OdjExMY7206ZNo1+/fkydOrXWOWu33XYbiYmJp90aMjTRzc2NoUOHsn79esc+u93O+vXrGTlyZL2PJ22EesRaDPWIiUirExQawbYh84hJeIChqf9l344r6DEw1uxYInKKqcO6cF7vIFKzC4kI9KKTn2eTnu/uu+9mwoQJTJw4EQAXFxf69etXZZ5YcnIyb7/9Nj179mTAgAH4+PiwZcsWXnzxRZ5//nn+/e9/VzlmcnIyK1asoF+/fowdO5Zvv/2W0aNHVzt3Q4cm5ufnV7nWWUpKComJiQQEBDh6z2bPns20adOIiYlh+PDhLFq0iIKCAmbMmFHv80kb0T6i8mt+BpQVgWvT/t+T2qkQE5FWaeglt/DD7jUMLvwOVt9Gae/NuLlriKJIS9LJz7PJCzCAtWvX8uWXX5KcnFxl/8CBA6sUYn369HGsGNivXz/HaoQDBw7k448/rnbcPn36EBkZCcDgwYNJTU2tsRBrqG3btjF27FjH7dmzZwMwbdo0li9fDsDUqVPJyspi3rx5pKenEx0dzbp166ot4CHi4Nke3NpB6Qk4fgCCtEqmWVSIiUirZLFaCb/hRXKWjKRHxS9sev0hRt74lNmxRMQEl1xyCTk5OdX2v/baa1Vuu5/yYY3VanXctlqtVFRUVHv8qe1tNluNbc7GmDFj6rRE/8yZM5k5c2ajnltaMYulcp5YRlLlPDEVYqbRHDERabUCO4bzy7BHAIg58DI/b//O3EAiIiItgeaJtQgqxESkVRsy4S8k+JyHq6UC2wd3UFJSZHYkERERczlWTkw1NUZbp0JMRFo1i9VKxA1LOIYv3eypJLz2oNmRRBrVFVdcQfv27bnqqquq7D948CBjxowhMjKSQYMG8c4775iU0DlERESwbds2x+13332XMWPGADBixAjWrl1bpd3v2z/99NNMnz69OSOLNJx6xFoEFWIi0uoFBHcmdcTfARh2aDk/JWw0OZFI47nnnnuqzXWCylUBFy1axK5du/jss8+YNWsWBQUFJiQUkRZH1xJrEVSIiUibMOSP04lvdwEuFjtua++guKjQ7EgijWLMmDG0a9eu2v5OnToRHR0NQMeOHQkMDOTYsWPNnE5EWiT1iLUIKsREpM3occMLHMWPCPtBfnjtAbPjSBuwceNGJk2aRGhoKBaLhdWrV1drs3jxYiIiIvDw8CA2Npa4uLhGzxEfH09FRQXh4eGNfmwRcUL+ldehozgXio6bGqUtUyEmIm2Gf1AnDo56AoDhaf9j97YvTU4krV1BQQFRUVEsXry4xvtXrlzJ7NmzmT9/PgkJCURFRTF+/HgyMzMdbaKjoxkwYEC1LS0trU4Zjh07xg033MB///vfRnlOItIKuPuAV2Dl9+oVM42uIyYibUr0RX9mW9IHxOR9jufHd1EcuQUPLx+zY0krNWHCBCZMmFDr/c8++yw333wzM2bMAGDJkiV89NFHLFu2jDlz5gCQmJjY4POXlJRw+eWXM2fOHEaNGnXadiUlJY7beXl5DT6niDiJ9l2hMLtynlinKLPTtEnqERORNqfXtMVk409X+yF+eO1+s+NIG1VaWkp8fDzjxo1z7LNarYwbN45Nmzad9fENw2D69OlccMEF/PnPfz5t24ULF+Ln5+fYNIRRpA3QPDHTqRATkTbHr0MIh8/9BwCxR1aQHPe5yYmkLcrOzqaiooKQkJAq+0NCQkhPT6/zccaNG8eUKVP4+OOPCQsLcxRx3333HStXrmT16tVER0cTHR3Njh07ajzG3Llzyc3NdWwHDx5s+BMTEeeglRNNp6GJItImRY27lq1Jqxl2fB0+n9xNUf84PL2rrzwn0tJ98cUXNe4/99xzsdvtdTqGu7s77u7ujRlLRFo69YiZzil6xL766issFkuN29atWwF45JFHarzf29v7tMfeunUrF154If7+/rRv357x48ezffv25nhaImKy3tMWk0kA4UYa21+91+w40sYEBgZis9nIyMiosj8jI4OOHTualEpE2gz1iJnOKQqxUaNGceTIkSrbTTfdRLdu3YiJiQHgvvvuq9YmMjKSKVOm1Hrc/Px8/vjHP9KlSxe2bNnCt99+S7t27Rg/fjxlZWXN9fRExCR+7QNJP/9JAIZnvM3OTZ+YnEjaEjc3N4YOHcr69esd++x2O+vXr2fkyJEmJhORNsHRI3YADMPcLG2UUwxNdHNzq/LpYFlZGR988AF33XUXFosFAB8fH3x8flv5bPv27ezatYslS5bUetzdu3dz7NgxHnvsMcfE5Pnz5zNo0CD2799Pz549m+gZiUhLMWjsFLb+uIphOR/h/9k9FAwYhXc7P7NjSSuRn5/Pzz//7LidkpJCYmIiAQEBdOnShdmzZzNt2jRiYmIYPnw4ixYtoqCgwLGKoohIk/ELByxQXgT5mdAu5IwPkcblFD1iv7dmzRqOHj162jeql156id69ezN69Oha2/Tp04cOHTrw8ssvU1paSlFRES+//DL9+vUjIiKi1seVlJSQl5dXZRMR59V32nOkE0hnI4Mdr/4/s+NIK7Jt2zYGDx7M4MGDAZg9ezaDBw9m3rx5AEydOpWnn36aefPmER0dTWJiIuvWrau2gIdIU0pJSWHs2LFERkYycOBACgoKzI4kzcHFDXw7V36veWKmcMpC7OWXX2b8+PGEhYXVeH9xcTFvvPEGN95442mP065dO7766itef/11PD098fHxYd26dXzyySe4uNTeWahlfkVal3b+Hci+4BkARmS/R9K3H5qcSFqLMWPGYBhGtW358uWONjNnzmT//v2UlJSwZcsWYmNjzQssbdL06dN57LHH2LVrF19//bUWbmlLNE/MVKYWYnPmzKl1EY6T2+7du6s85tChQ3z66aenLbJWrVrFiRMnmDZt2mnPX1RUxI033sg555zD5s2b+e677xgwYAATJ06kqKio1sdpmV+R1mfAeZcT1+EyAALWzyY/L8fkRCLSFObMmYO7uzvXXXed2VFahJ07d+Lq6uoYQRQQEHDaD6OllXHME0s1NUZbZWohdu+995KcnHzarXv37lUe88orr9ChQwcuvfTSWo/70ksvcckll5xxaMebb75Jamoqr7zyCsOGDWPEiBG8+eabpKSk8MEHH9T6OHd3d3x9fatsIuL8Iqf9myOWIEKNTHa+OsvsOCLSBObOncszzzzDihUrqszfa4k2btzIpEmTCA0NxWKxsHr16hrbLV68mIiICDw8PIiNjSUuLq7O59i7dy8+Pj5MmjSJIUOGsGDBgkZKL05BPWKmMvUjj6CgIIKCgurc3jAMXnnlFW644QZcXV1rbJOSksKGDRtYs2bNGY9XWFiI1Wp1LPgBOG7X9dorItJ6+Pi2J/XCZ+n0xZ+JPbqaH79exaDzrzA7log0Ij8/P2688UbuueceduzY0aIX5iooKCAqKoq//OUvTJ48ucY2K1euZPbs2SxZsoTY2FgWLVrE+PHj2bNnD8HBwQBER0dTXl5e7bGfffYZ5eXlfPPNNyQmJhIcHMwf//hHhg0bxh/+8IcmfW7SQuhaYqZyqjliX375JSkpKdx00021tlm2bBmdOnViwoQJ1e5btWoVffv2ddz+wx/+QE5ODnfeeSfJycns3LmTGTNm4OLiwtixY5vkOYhIyzbg3EvZElj5B0/whvvIO37U5EQirVjuYUjZWPm1GZWXl+Pl5UVSUlKznre+JkyYwOOPP84VV9T+gdCzzz7LzTffzIwZM4iMjGTJkiV4eXmxbNkyR5vExESSkpKqbaGhoXTu3JmYmBjCw8Nxd3fn4osvJjExsRmenbQI6hEzlVMVYi+//DKjRo2qUkydym63s3z5cqZPn47NZqt2f25uLnv27HHc7tu3Lx9++CE//vgjI0eOZPTo0aSlpbFu3To6derUZM9DRFq2gdP/xWFLCB3JZverd5kdR6RlMwwoLaj/FrcUFg2AVydVfo1bWv9jNPDaRw899BD5+fnNVogtWLDAcZmd2rYDBw7U+7ilpaXEx8czbtw4xz6r1cq4cePYtGlTnY4xbNgwMjMzycnJwW63s3HjRvr161fvLOKkTvaI5R6Ciuq9ptK0nGo25ptvvnna+61W62kXzpg+fTrTp0+vsu8Pf/iDut9FpAovH3/2X/RvOn96DcNzPmL7hneIGlv7xeFF2rSyQlgQenbHMOzw8X2VW308mAZu3vV6SHx8PEuWLGHixInNVojddtttXH311adtExpa/3/D7OxsKioqqs2JDwkJqbbYWW1cXFxYsGAB5513HoZhcNFFF3HJJZfUO4s4qXadwOYGFaWQd/i3HjJpFk5ViImINJd+Iyew+YerGZH5Np2+/iu5UWPwC6j7nFYRaXnsdju33norM2fOJDY2luuvv56ysrJa553XJC0tjfvvv5833nijzo8JCAggICCgIZGbxYQJE2qc0iFtgNVaeWHnY/sq54mpEGtWKsRERGoRNe1ZDj39DWHGEba+OpNh/2+l2ZFEWh5Xr8qeqfrIS4PFwyt7wk6y2ODOLeBbj54hV696nfa5554jOzubxx57jAMHDlBWVsbu3bsZOHBgnY8RGhparyIMKocmnmk1wl27dtGlS5d6HTcwMBCbzUZGRkaV/RkZGXTs2LFex5I2rH3XykIsZz90MztM2+JUc8RERJqTp3c78if8H3bDwrDcdSR+scLsSCItj8VSOTywPltgL5j078riCyq/TlpUub8+xzll1eMzOXz4MA8//DCLFy/G29ubXr164e7u7hiemJqaSlRUFH/605/o1asXt99+O6tXryY2NpYBAwawd+9eR7uYmBhH+2nTptGvXz+mTp2KUcuctdtuu43ExMTTbg0Zmujm5sbQoUNZv369Y5/dbmf9+vWMHDmy3seTNkorJ5pGPWIiIqfRd/hFbE64lhHpbxL27RxyB1+AX4fTX6NQROpgyA3Q40I49gsEdAe/zk16urvvvpsJEyYwceJEoHJuVL9+/arME0tOTubtt9+mZ8+eDBgwAB8fH7Zs2cKLL77I888/z7///e8qx0xOTmbFihX069ePsWPH8u233zoujHyqhg5NzM/Pr3Kts5SUFBITEwkICHD0ns2ePZtp06YRExPD8OHDWbRoEQUFBcyYMaPe55M2SisnmkY9YiIiZxB9w1McsHYmkOPsffUOs+OItB5+naHb6CYvwtauXcuXX35ZrZAaOHBglUKsT58+9OnTB5vNRr9+/RyrEQ4cOJDU1NRqx+3Tpw+RkZFYLBYGDx5cY5uzsW3bNgYPHszgwYOByqJr8ODBzJs3z9Fm6tSpPP3008ybN4/o6GgSExNZt25dtQU8RGqlHjHTqEdMROQMPLx8KLr4OSo+vJKYvC9I+PR/DBn/Z7NjiUgdXXLJJeTk5FTb/9prr1W57e7u7vjearU6blutVioqKqo9/tT2NputxjZnY8yYMbUOdzzVzJkzmTlzZqOeW9oQ9YiZRj1iIiJ10CfmQrZ2riy+um2aS8Jn/yPj0D5TM2Uc2kfSdx8qh4iINJx/ROXX/HQ4avLruEkXeTcrh3rEGlFBQUGNF5K22Wx4eHhUaVcbq9WKp6dng9oWFhbW+smZxWLBy8urQW2Lioqw2+01tgXw9vZuUNvi4uLTfnpYn7ZeXl5Yfp20XVJSQnl57RclrE9bT09PrNbKzytKS0spKytrlLYeHh6O35X6tC0rK6O0tLTWtu7u7ri4uNS7bXl5OSUlJbW2dXNzcyzvXJ+2FRUVFBcX19rW1dUVNze3ere12+0UFRU1SlsXFxfHp9qGYVBYWFhr28ipj5L1rw8JIgf/7+4kf4OFDQGX4Nbr/GptrVYL7r9mACgqrv3frD5tLRYLHu5ulP60geFHV+NTZsdu1JzjZNu6HBfA0+O3T/fr2rb0pw0MTF+FDwYVhoWvBzxIzOVVP5k/02vE6V7nRESkie1e+9v3zw2FmL9A9zHNn+OXr2DbMsAALC0jh8VaubDQkBua5FQWoy593nJaeXl5+Pn51Xr/xRdfzEcffeS47e3tXesfe+effz5fffWV43ZQUBDZ2dk1to2JiWHr1q2O2xEREezfX3O3cmRkJDt37nTc7t+/P7t27aqxbdeuXauMcx82bBjbtm2rsW1gYCBZWVmO22PGjOHrr7+usa2Xl1eVP7gmTpzIxx9/XGNboEqhOGXKFN59991a2+bn5zv+2Js+fTqvvvpqrW0zMzMJCqq8HtSdd97JCy+8UGvblJQUIiIiALj//vt5+umna22blJRE//79AXjkkUd49NFHa20bFxfHsGHDAHjqqaf461//WmvbDRs2MGbMGAAWL1582uEna9eudUxEX758+Wkna7/99ttMmVJ5keJ33nnntBcbfeWVVxwXQ//oo49Oe7HP559/njvvvBOAr776irFjx9ba9sknn+T+++8HYOvWrQwfPrzWtvPnz+eRRx4BYOfOnQwYMKDWtvfddx9PPfUUULnCWbduta/He8cdd7B48WIAsrKyCA4OrrXt1VdNZkXkeqwWg4JSA5+FJ2pte1WkC+9M+e0DDcujebW2vbiXCx9d91tb7wV5FNZSm5/f1cZX038rbIKeOkF2Yc0v4zGhVrbe7OO4HbHoBPtza24bGWRl5x2/te3/Qj67smr+UKWrn4XUWe0ct4ctzWdbWs1t6/MakZubi6+vb433SfM6+b5W08+kuLiYlJQUunXrVuVDRmkb9PNvZXIPw6IBVS8lIVVZbDBrR73msp7uNfRU6hETEamj4oLjWC11++zqhMWX3a69T9mzuda2+RYfdrv2ddw2iANqflMstHiRYu1CN/uBM+e1eLLbNdJxu8ySANTcS1pica/StsSyHai5J7HM4spu10jcK/LrlENERFqoY/tqLsKC+oFH7Z0Mja44F7KSW2YOo6JyddcmWFRIPWKN4GTVm5aWVmPVq6GJNbfV0EQNTXS2oYlH0/fT+bVR2CxGZdsyKDesHJ3xDSGdu1dpW5//9/V9jcg7mkbg0qHYfu2Zg5pzNPVrRMahfQQuHUppuR37r01rynGm14i8vDxCQ0PVI9aCqEdMaqOffytTU49YA3qAlKMq9YiZwNvbu8ofHKdrV59j1tWpf0Q1ZttT/5BrzLb1eQGvT1t3d/cqK1k1Vls3NzfHH/dmtXV1dXUUOY3Z1sXFxVGUNWZbm81W59/h+rS1Wq1N0tZisZy2rXePSOIGPcKQHx/FxWLH3dXKzkHzGd574BmP3Zj/7z3DejhyeLvZKTfqlqOxXyNCTsnhYqlbjppeIxp7pTkREakjv86Vc6A+nFXZ83Py4urNWfy04RzqEWsEda16RaR1yDi0j+z9uwns2peQsB7KcZY59Bra8qhHTGqjn38rlXu42S6u3hZyqEdMRKSJhIT1MLXwUQ5pCfQ5btukn3sr5dfZ3MKnjebQdcRERESkzk6d3yptz8l5tHUd+i4itVOPmIiIiNSZi4sLXl5eZGVl4erq6likSFq3k4sZZWZm4u/vX+N1U0WkflSINaKCggLatWvnWJHv5Gp4p67KdrIdVF1l7+QKd7WtnlaXtidXOTt1lb2TK9zVtnpaXdqeXOXs1FX2Tq5wV5+2v1+V7eRKiDWtsleXtqeuhnfq4gYnV0KsbeW8M7U9deW8U1dYPPnzrE/buvzsG+P3pKafZ2P8npz8eZ7t70ltK3E29Pektp/n2f6enPrzPNvfk9p+ng39PWntrxHiPCwWC506dSIlJaXWa1dK6+Xv70/Hjh3NjiHSOhhy1nJzcw0qLwNuZGZmOvY//vjjBmDcdNNNVdp7eXkZgJGSkuLY969//csAjOuuu65K28DAQAMwkpKSHPv++9//GoBx2WWXVWnbtWtXAzDi4uIc+15//XUDMMaNG1elbWRkpAEYGzZscOxbtWqVARijRo2q0jYmJsYAjLVr1zr2ffbZZwZgREVFVWl7/vnnG4Dx9ttvO/Z9++23BmD07NmzStuLL77YAIxXXnnFse+HH34wACM0NLRK26uuusoAjOeff96x76effjIAw8/Pr0rbadOmGYDx5JNPOvYdOnTIAAwXF5cqbe+44w4DMObPn+/Yl5OT4/h5lpaWOvbfd999BmDcd999jn2lpaWOtjk5OY798+fPNwDjjjvuqHI+FxcXAzAOHTrk2Pfkk08agDFt2rQqbf38/AzA+Omnnxz7nn/+eQMwrrrqqiptQ0NDDcD44YcfHPteeeUVAzAuvvjiKm179uxpAMa3337r2Pf2228bgHH++edXaRsVFWUAxmeffebYt3btWgMwYmJiqrQdNWqUARirVq1y7NuwYYMBGJGRkVXajhs3zgCM119/3bEvLi7OAIyuXbtWaXvZZZcZgPHf//7XsS8pKckAjMDAwCptr7vuOgMw/vWvfzn2paSkGIDh5eVVpe1NN91kAMbjjz/u2JeZmen4eZ7qnnvuMQDjwQcfdOzLz893tM3Pz3fsf/DBBw3AuOeee6ocQ68RlWp6jfj0008NwMjNzTWkZTj5vna6n0lFRYVRVFSkrQ1t5eXlzfhbKOK86vIaahiGoR4xERERqTer1apV80REzoKWr28Ep17QuWPHjhp2pKGJGpqooYkamliP14icnBwCAgK0fH0LoksKiIg0XF1fQ1WINQK9YYmINJxeQ1se/UxERBqurq+hWupIRERERESkmWmOWCM42amYl5dnchIREedz8rVTAzRaDr2viYg0XF3f11SINYITJ04AEB4ebnISERHndeLECfz8/MyOIeh9TUSkMZzpfU1zxBqB3W4nLS2tyjXE6iovL4/w8HAOHjxo6jh85VAO5VAOs3IYhsGJEycIDQ3VxYFbCL2vKYdyKIdyNP37mnrEGoHVaiUsLOysjuHr69siJkQrh3Ioh3KYkUM9YS2L3teUQzmUQzma/n1NHz2KiIiIiIg0MxViIiIiIiIizUyFmMnc3d2ZP39+lYu5KodyKIdyKIc4q5by+6AcyqEcytHSc2ixDhERERERkWamHjEREREREZFmpkJMRERERESkmakQExERERERaWYqxERERERERJqZCjGTLFy4kGHDhtGuXTuCg4O5/PLL2bNnj6mZ/vGPf2CxWJg1a5Yp5z98+DDXX389HTp0wNPTk4EDB7Jt27ZmzVBRUcHDDz9Mt27d8PT0pEePHvz973+nqde02bhxI5MmTSI0NBSLxcLq1aur3G8YBvPmzaNTp054enoybtw49u7d26w5ysrKeOCBBxg4cCDe3t6EhoZyww03kJaW1qw5fu+2227DYrGwaNEiU3IkJydz6aWX4ufnh7e3N8OGDePAgQPNmiM/P5+ZM2cSFhaGp6cnkZGRLFmypFEzQN1et4qLi7nzzjvp0KEDPj4+XHnllWRkZDR6Fml59L5Wnd7X9L5Wlxy/p/e1tvO+pkLMJF9//TV33nknmzdv5vPPP6esrIyLLrqIgoICU/Js3bqVF198kUGDBply/pycHM455xxcXV355JNP2LVrF8888wzt27dv1hz//Oc/+c9//sPzzz9PcnIy//znP3nyySd57rnnmvS8BQUFREVFsXjx4hrvf/LJJ/m///s/lixZwpYtW/D29mb8+PEUFxc3W47CwkISEhJ4+OGHSUhI4P3332fPnj1ceumljZrhTDlOtWrVKjZv3kxoaGijZ6hLjn379nHuuefSt29fvvrqK3788UcefvhhPDw8mjXH7NmzWbduHa+//jrJycnMmjWLmTNnsmbNmkbNUZfXrf/3//4fH374Ie+88w5ff/01aWlpTJ48uVFzSMuk97Wq9L6m97W65jiV3tcqtZn3NUNahMzMTAMwvv7662Y/94kTJ4xevXoZn3/+uXH++ecb99xzT7NneOCBB4xzzz232c/7exMnTjT+8pe/VNk3efJk409/+lOzZQCMVatWOW7b7XajY8eOxlNPPeXYd/z4ccPd3d1YsWJFs+WoSVxcnAEY+/fvb/Ychw4dMjp37mwkJSUZXbt2Nf71r381WYbackydOtW4/vrrm/S8dcnRv39/47HHHquyb8iQIcbf/va3Js3y+9et48ePG66ursY777zjaJOcnGwAxqZNm5o0i7Q8el/T+9pJel+rWw69r/2mrbyvqUeshcjNzQUgICCg2c995513MnHiRMaNG9fs5z5pzZo1xMTEMGXKFIKDgxk8eDBLly5t9hyjRo1i/fr1/PTTTwBs376db7/9lgkTJjR7lpNSUlJIT0+v8vPx8/MjNjaWTZs2mZYLKn9vLRYL/v7+zXpeu93On//8Z+6//3769+/frOc+NcNHH31E7969GT9+PMHBwcTGxp52uElTGTVqFGvWrOHw4cMYhsGGDRv46aefuOiii5r0vL9/3YqPj6esrKzK72rfvn3p0qWL6b+r0vz0vqb3tdrofa06va9V1Vbe11SItQB2u51Zs2ZxzjnnMGDAgGY991tvvUVCQgILFy5s1vP+3i+//MJ//vMfevXqxaeffsrtt9/O3XffzauvvtqsOebMmcM111xD3759cXV1ZfDgwcyaNYs//elPzZrjVOnp6QCEhIRU2R8SEuK4zwzFxcU88MADXHvttfj6+jbruf/5z3/i4uLC3Xff3aznPVVmZib5+fn84x//4I9//COfffYZV1xxBZMnT+brr79u1izPPfcckZGRhIWF4ebmxh//+EcWL17Meeed12TnrOl1Kz09HTc3t2p/wJj9uyrNT+9rel87Hb2vVaf3tarayvuay1kfQc7anXfeSVJSEt9++22znvfgwYPcc889fP75540+9re+7HY7MTExLFiwAID/396dB1VV/n8AfxOIXBZBQJarLNcNgVARHQQN3HKZhqBM0AhRmXIEVMzQRsWo3EJFcctcUdMIx0EDU8ACcnAXcUFCRBMaF9Ik1FDg8nz/8Mf5eQMEEQ86vl8zd6azPedzTo/n7XPuuUdXV1dcuHAB69evR3BwsGx1JCYmYufOndi1axecnZ2Rm5uLiIgIKJVKWet42VVVVcHf3x9CCHz77bey7vv06dOIi4tDTk4OtLS0ZN33k2pqagAAvr6+mDFjBgCgd+/eOHLkCNavXw9vb2/Zalm9ejWOHTuGn376CXZ2dvjtt98QFhYGpVL5wr4RaK3rFr0amGvMtVcNc4251hrXLX4j1srCw8ORkpKCjIwMdOrUSdZ9nz59GqWlpejTpw90dHSgo6ODrKwsrFq1Cjo6OlCr1bLVYm1tDScnJ415jo6OLf6WnsZERkZKdw9dXFwQFBSEGTNmtOqdVSsrKwCo84aeW7duScvkVBtW165dQ3p6uux3DQ8fPozS0lLY2tpK/fbatWuYOXMm7O3tZavD3NwcOjo6rd5vKyoqMGfOHMTGxsLHxwc9e/ZEeHg4AgICsGzZsheyz4auW1ZWVqisrERZWZnG+q3VV6l1MNceY641jLmmibmm6XXKNQ7EWokQAuHh4UhKSsKvv/4KlUolew1Dhw7F+fPnkZubK3369u2LwMBA5ObmQltbW7ZaBgwYUOd1oZcuXYKdnZ1sNQCP36D0xhuafyy0tbWlu0StQaVSwcrKCr/88os0r7y8HMePH4eHh4estdSGVWFhIQ4dOgQzMzNZ9w8AQUFBOHfunEa/VSqViIyMRGpqqmx16Orqol+/fq3eb6uqqlBVVSVLv23suuXm5oY2bdpo9NWCggIUFxfL3ldJfsw1Tcy1hjHXNDHXNL1OucZHE1tJWFgYdu3ahX379sHIyEh6ztTY2BgKhUKWGoyMjOo8u29gYAAzMzPZn+mfMWMGPD09sWjRIvj7++PEiRPYsGEDNmzYIGsdPj4+WLhwIWxtbeHs7IwzZ84gNjYWkyZNeqH7vX//Pi5fvixNX716Fbm5uTA1NYWtrS0iIiKwYMECdOvWDSqVClFRUVAqlfDz85OtDmtra3zwwQfIyclBSkoK1Gq11G9NTU2hq6srSx22trZ1grJNmzawsrKCg4NDi9XQlDoiIyMREBAALy8vDB48GAcPHkRycjIyMzNlrcPb2xuRkZFQKBSws7NDVlYWtm/fjtjY2Bato7HrlrGxMUJCQvDpp5/C1NQU7dq1w9SpU+Hh4YH+/fu3aC308mGuaWKuMdeaWgdz7TXOted+7yI1C4B6P1u3bm3VulrrNb9CCJGcnCzefPNN0bZtW9GjRw+xYcMG2WsoLy8X06dPF7a2tkJPT0907txZzJ07Vzx69OiF7jcjI6Pe/hAcHCyEePyq36ioKGFpaSnatm0rhg4dKgoKCmSt4+rVqw3224yMDNnqqM+Les1vU+rYvHmz6Nq1q9DT0xO9evUSe/fulb2OGzduiAkTJgilUin09PSEg4ODWL58uaipqWnROppy3aqoqBChoaGiffv2Ql9fX7z33nvixo0bLVoHvZyYa3Ux15hrTamjPsy11yPXtP6vCCIiIiIiIpIJfyNGREREREQkMw7EiIiIiIiIZMaBGBERERERkcw4ECMiIiIiIpIZB2JEREREREQy40CMiIiIiIhIZhyIERERERERyYwDMXqlxMfHw8TEpLXLeOXY29tj5cqVrbLvQYMGISIi4pm2iY6ORu/evaXpCRMmwM/Pr0XrehFa8zwT0auJudY8zDV5MNdeLA7E6JUSEBCAS5cutXYZTZaZmQktLS20b98eDx8+1Fh28uRJaGlpQUtLq876tR9LS0uMHj0aV65ckdY5e/Ys3n33XVhYWEBPTw/29vYICAhAaWmpbMclt7i4OMTHx7d2GY06efIkPvnkk9Yug4heIcw15trLjLn2YnEgRq8UhUIBCwuL1i7jmRkZGSEpKUlj3ubNm2Fra1vv+gUFBbh+/Tp2796NvLw8+Pj4QK1W46+//sLQoUNhamqK1NRU5OfnY+vWrVAqlXjw4IEch9IqjI2NX4k7xh06dIC+vn5rl0FErxDmGnPtZcZce7E4EKNmGTRoEKZOnYqIiAi0b98elpaW2LhxIx48eICJEyfCyMgIXbt2xYEDB6Rt1Go1QkJCoFKpoFAo4ODggLi4OGn5w4cP4ezsrHHnpaioCEZGRtiyZQuAuo9w1H7Vv2XLFtja2sLQ0BChoaFQq9WIiYmBlZUVLCwssHDhQmmbP/74A1paWsjNzZXmlZWVQUtLC5mZmQD+/w5eamoqXF1doVAoMGTIEJSWluLAgQNwdHREu3bt8OGHH+Lff/9t9HwFBwdLxwAAFRUVSEhIQHBwcL3rW1hYwNraGl5eXpg/fz4uXryIy5cvIzs7G//88w82bdoEV1dXqFQqDB48GCtWrIBKpXpqDffu3cO4ceNgYGCAjh07Yu3atRrLi4uL4evrC0NDQ7Rr1w7+/v64detWnXO9Y8cO2Nvbw9jYGGPHjsW9e/ekdR48eIDx48fD0NAQ1tbWWL58eaPnBgCWLFkCS0tLGBkZISQkpM5d1v8+wtGc/gcAFy5cwKhRo2BoaAhLS0sEBQXh9u3bGu1OmzYNs2bNgqmpKaysrBAdHS0tF0IgOjoatra2aNu2LZRKJaZNmyYt/+8jHC1xTolIHsw15hpzjbkmNw7EqNm2bdsGc3NznDhxAlOnTsWUKVMwZswYeHp6IicnB8OHD0dQUJB0Qa+pqUGnTp2we/duXLx4EfPnz8ecOXOQmJgIANDT08POnTuxbds27Nu3D2q1Gh999BHefvttTJo0qcE6ioqKcODAARw8eBA//PADNm/ejHfeeQd//vknsrKy8M0332DevHk4fvz4Mx9jdHQ01qxZgyNHjqCkpAT+/v5YuXIldu3ahf379yMtLQ2rV69utJ2goCAcPnwYxcXFAIA9e/bA3t4effr0aXRbhUIBAKisrISVlRWqq6uRlJQEIcQzHcvSpUvRq1cvnDlzBp9//jmmT5+O9PR0AI//3/j6+uLvv/9GVlYW0tPTceXKFQQEBGi0UVRUhL179yIlJQUpKSnIysrCkiVLpOWRkZHIysrCvn37kJaWhszMTOTk5Dy1rsTERERHR2PRokU4deoUrK2tsW7dukaP51n7X1lZGYYMGQJXV1ecOnUKBw8exK1bt+Dv71+nXQMDAxw/fhwxMTH46quvpPO0Z88erFixAt999x0KCwuxd+9euLi41FtfS51TIpIPc425xlxjrslKEDWDt7e3GDhwoDRdXV0tDAwMRFBQkDTvxo0bAoA4evRog+2EhYWJ0aNHa8yLiYkR5ubmIjw8XFhbW4vbt29Ly7Zu3SqMjY2l6S+++ELo6+uL8vJyad6IESOEvb29UKvV0jwHBwexePFiIYQQV69eFQDEmTNnpOV3794VAERGRoYQQoiMjAwBQBw6dEhaZ/HixQKAKCoqkuZNnjxZjBgxosHjq23n7t27ws/PT3z55ZdCCCEGDx4s4uLiRFJSknjyj+GT6wshxPXr14Wnp6fo2LGjePTokRBCiDlz5ggdHR1hamoqRo4cKWJiYsTNmzcbrEEIIezs7MTIkSM15gUEBIhRo0YJIYRIS0sT2traori4WFqel5cnAIgTJ04IIeo/15GRkcLd3V0IIcS9e/eErq6uSExMlJbfuXNHKBQKMX369AZr8/DwEKGhoRrz3N3dRa9evaTp4OBg4evrK003p/99/fXXYvjw4Rr7KSkpEQBEQUFBve0KIUS/fv3E7NmzhRBCLF++XHTv3l1UVlbWeyx2dnZixYoVQoiWOadEJB/m2mPMNebak5hrLxa/EaNm69mzp/Tf2traMDMz07iLYmlpCQAaP7Zdu3Yt3Nzc0KFDBxgaGmLDhg3S3bRaM2fORPfu3bFmzRps2bIFZmZmT63D3t4eRkZGGvt1cnLCG2+8oTGvOT/6ffIYLS0toa+vj86dOzer3UmTJiE+Ph5XrlzB0aNHERgY2OC6nTp1goGBgfSM/J49e6CrqwsAWLhwIW7evIn169fD2dkZ69evR48ePXD+/Pmn7t/Dw6POdH5+PgAgPz8fNjY2sLGxkZY7OTnBxMREWgeoe66tra2l4y8qKkJlZSXc3d2l5aampnBwcHhqXfn5+Rrb1FdrfZ61/509exYZGRkwNDSUPj169JBqr6/d/x7jmDFjUFFRgc6dO+Pjjz9GUlISqqurGzyu5z2nRCQv5hpzjbnGXJMTB2LUbG3atNGY1tLS0phX+9akmpoaAEBCQgI+++wzhISEIC0tDbm5uZg4cSIqKys12iktLcWlS5egra2NwsLC566jdl5tHbVBJp54BKKqqqrRthtrtzGjRo1CRUUFQkJC4OPj89QgPnz4MM6dO4fy8nLk5ubWuaCbmZlhzJgxWLZsGfLz86FUKrFs2bIm1fE8nuf45ajlaf3v/v378PHxQW5ursansLAQXl5eT223tg0bGxsUFBRg3bp1UCgUCA0NhZeXV4P9p7nH0VrnlOh1x1xjrjHXmGty4kCMZJOdnQ1PT0+EhobC1dUVXbt21bhjU2vSpElwcXHBtm3bMHv2bI27LC2hQ4cOAIAbN25I8578gfOLoqOjg/HjxyMzM/Opvw0AAJVKhS5dumjcUWqIrq4uunTp0ujbpY4dO1Zn2tHREQDg6OiIkpISlJSUSMsvXryIsrIyODk5NVoDAHTp0gVt2rTR+M3C3bt3G30ts6OjY53fOfy31pbQp08f5OXlwd7eHl27dtX4GBgYNLkdhUIBHx8frFq1CpmZmTh69Gi9d21b4pwS0cuNucZcqw9zjZpKp7ULoNdHt27dsH37dqSmpkKlUmHHjh04efKkxluR1q5di6NHj+LcuXOwsbHB/v37ERgYiGPHjkmPMDwvhUKB/v37Y8mSJVCpVCgtLcW8efNapO3GfP3114iMjGz0sZSGpKSkICEhAWPHjkX37t0hhEBycjJ+/vlnbN269anbZmdnIyYmBn5+fkhPT8fu3buxf/9+AMCwYcPg4uKCwMBArFy5EtXV1QgNDYW3tzf69u3bpNoMDQ0REhIiHZ+FhQXmzp2r8ShNfaZPn44JEyagb9++GDBgAHbu3Im8vDyNR2VaQlhYGDZu3Ihx48ZJb4+6fPkyEhISsGnTJmhrazfaRnx8PNRqNdzd3aGvr4/vv/8eCoUCdnZ2ddZtiXNKRC835hpzrT7MNWoqfiNGspk8eTLef/99BAQEwN3dHXfu3EFoaKi0/Pfff0dkZCTWrVsnPX+8bt063L59G1FRUS1ay5YtW1BdXQ03NzdERERgwYIFLdp+Q3R1dWFubq7xj10+CycnJ+jr62PmzJno3bs3+vfvj8TERGzatAlBQUFP3XbmzJk4deoUXF1dsWDBAsTGxmLEiBEAHj82sG/fPrRv3x5eXl4YNmwYOnfujB9//PGZ6lu6dCneeust+Pj4YNiwYRg4cCDc3Nyeuk1AQACioqIwa9YsuLm54dq1a5gyZcoz7bcplEolsrOzoVarMXz4cLi4uCAiIgImJiaNhmotExMTbNy4EQMGDEDPnj1x6NAhJCcn1/sXkJY6p0T08mKuMdfqw1yjptIS4hnfFUpERERERETPhd+IERERERERyYwDMSIiIiIiIplxIEZERERERCQzDsSIiIiIiIhkxoEYERERERGRzDgQIyIiIiIikhkHYkRERERERDLjQIyIiIiIiEhmHIgRERERERHJjAMxIiIiIiIimXEgRkREREREJDMOxIiIiIiIiGT2P4GXVWq2di3IAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from copy import deepcopy\n", + "\n", + "import matplotlib.gridspec as gridspec\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.ticker import MaxNLocator\n", + "\n", + "svd_min_list = [1e-3, 1e-6]\n", + "chi_max_list = np.arange(2, 21, 2)\n", + "lucj_mps_energy = np.zeros((2, len(chi_max_list)))\n", + "\n", + "# Construct Hartree-Fock state\n", + "initial_mps = statevector_to_mps(ffsim.hartree_fock_state(norb, nelec), norb, nelec)\n", + "\n", + "# Loop over cutoff and bond dimension\n", + "for i, svd_min in enumerate(svd_min_list):\n", + " for j, chi_max in enumerate(chi_max_list):\n", + " final_mps = deepcopy(initial_mps)\n", + " options = {\"trunc_params\": {\"chi_max\": int(chi_max), \"svd_min\": svd_min}}\n", + " eng = TEBDEngine(final_mps, None, options)\n", + " apply_ucj_op_spin_balanced(eng, lucj_op)\n", + " lucj_mps_energy[i, j] = hamiltonian_mpo.expectation_value_finite(final_mps)\n", + "\n", + "fig = plt.figure(figsize=(10, 4))\n", + "gs = gridspec.GridSpec(1, 2, wspace=0.3)\n", + "ax0 = plt.subplot(gs[0])\n", + "ax1 = plt.subplot(gs[1])\n", + "\n", + "for i in [0, 1]:\n", + " ax0.plot(\n", + " chi_max_list,\n", + " lucj_mps_energy[i, :],\n", + " \".-\",\n", + " label=f\"$\\\\lambda_\\\\text{{min}}=10^{{{np.log10(svd_min_list[i]):g}}}$\",\n", + " )\n", + "\n", + "ax0.set_xlabel(\"maximum MPS bond dimension\")\n", + "ax0.set_ylabel(\"$E$\")\n", + "ax0.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + "ax0.axhline(y=lucj_energy, color=\"k\", linestyle=\"dashed\", label=\"$E_\\\\text{LUCJ}$\")\n", + "ax0.axhline(y=fci_energy, color=\"k\", linestyle=\"dotted\", label=\"$E_\\\\text{FCI}$\")\n", + "ax0.legend(loc=\"best\")\n", + "\n", + "for i in [0, 1]:\n", + " ax1.plot(\n", + " chi_max_list,\n", + " np.abs(np.subtract(lucj_mps_energy[i, :], lucj_energy)),\n", + " \".-\",\n", + " label=f\"$\\\\lambda_\\\\text{{min}}=10^{{{np.log10(svd_min_list[i]):g}}}$\",\n", + " )\n", + "\n", + "ax1.set_xlabel(\"maximum MPS bond dimension\")\n", + "ax1.set_ylabel(\"$|E-E_\\\\text{LUCJ}|$\")\n", + "ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + "ax1.set_yscale(\"log\")\n", + "ax1.legend(loc=\"best\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b2f8fbd2-b019-4d38-a4f1-62afcf238e3c", + "metadata": {}, + "source": [ + "From the above plots, we can see that at an MPS bond dimension of 16 or above, the MPS representation of the LUCJ circuit is exact." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index c78e60552..ac71e80ca 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ dependencies = [ "numpy", "opt_einsum", "orjson", + "physics-tenpy >= 1.0.5", "pyscf >= 2.7", "qiskit >= 1.1", "scipy", diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py index 45fbbb409..806e7b2a4 100644 --- a/python/ffsim/__init__.py +++ b/python/ffsim/__init__.py @@ -10,7 +10,7 @@ """ffsim is a software library for fast simulation of fermionic quantum circuits.""" -from ffsim import contract, linalg, optimize, qiskit, random, testing +from ffsim import contract, linalg, optimize, qiskit, random, tenpy, testing from ffsim.cistring import init_cache from ffsim.gates import ( apply_diag_coulomb_evolution, @@ -175,6 +175,7 @@ "slater_determinant_rdms", "spin_square", "strings_to_addresses", + "tenpy", "testing", "trace", ] diff --git a/python/ffsim/tenpy/__init__.py b/python/ffsim/tenpy/__init__.py new file mode 100644 index 000000000..743b9fc87 --- /dev/null +++ b/python/ffsim/tenpy/__init__.py @@ -0,0 +1,43 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Code that uses TeNPy, e.g. for emulating quantum circuits.""" + +from ffsim.tenpy.gates.abstract_gates import ( + apply_single_site, + apply_two_site, +) +from ffsim.tenpy.gates.basic_gates import ( + givens_rotation, + num_interaction, + num_num_interaction, + on_site_interaction, +) +from ffsim.tenpy.gates.diag_coulomb import apply_diag_coulomb_evolution +from ffsim.tenpy.gates.orbital_rotation import apply_orbital_rotation +from ffsim.tenpy.gates.ucj import apply_ucj_op_spin_balanced +from ffsim.tenpy.random.random import random_mps, random_mps_product_state +from ffsim.tenpy.util import mps_to_statevector, statevector_to_mps + +__all__ = [ + "apply_ucj_op_spin_balanced", + "apply_diag_coulomb_evolution", + "apply_orbital_rotation", + "apply_single_site", + "apply_two_site", + "givens_rotation", + "mps_to_statevector", + "num_interaction", + "num_num_interaction", + "on_site_interaction", + "random_mps", + "random_mps_product_state", + "statevector_to_mps", +] diff --git a/python/ffsim/tenpy/gates/__init__.py b/python/ffsim/tenpy/gates/__init__.py new file mode 100644 index 000000000..7918d6a4b --- /dev/null +++ b/python/ffsim/tenpy/gates/__init__.py @@ -0,0 +1,9 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/python/ffsim/tenpy/gates/abstract_gates.py b/python/ffsim/tenpy/gates/abstract_gates.py new file mode 100644 index 000000000..dcb56d216 --- /dev/null +++ b/python/ffsim/tenpy/gates/abstract_gates.py @@ -0,0 +1,85 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy abstract gates.""" + +import numpy as np +import tenpy.linalg.np_conserved as npc +from numpy.typing import NDArray +from tenpy.algorithms.tebd import TEBDEngine +from tenpy.linalg.charges import LegPipe +from tenpy.networks.site import SpinHalfFermionSite + +# ignore lowercase argument and variable checks to maintain TeNPy naming conventions +# ruff: noqa: N803, N806 + +# define sites +shfs = SpinHalfFermionSite(cons_N="N", cons_Sz="Sz") +shfsc = LegPipe([shfs.leg, shfs.leg]) + + +def apply_single_site(eng: TEBDEngine, U1: NDArray[np.complex128], site: int) -> None: + r"""Apply a single-site gate to an MPS. + + Args: + eng: The TEBD engine. + U1: The single-site quantum gate. + site: The gate will be applied to `site` on the MPS. + + Returns: + None + """ + U1_npc = npc.Array.from_ndarray(U1, [shfs.leg, shfs.leg.conj()], labels=["p", "p*"]) + psi = eng.get_resume_data()["psi"] + psi.apply_local_op(site, U1_npc) + + +def apply_two_site( + eng: TEBDEngine, + U2: NDArray[np.complex128], + sites: tuple[int, int], + *, + norm_tol: float = 1e-8, +) -> None: + r"""Apply a two-site gate to an MPS. + + Args: + eng: The TEBD engine. + U2: The two-site quantum gate. + sites: The gate will be applied to adjacent sites `(site1, site2)` on the MPS. + norm_tol: The norm error above which we recanonicalize the MPS. In general, the + application of a two-site gate to an MPS with truncation may degrade its + canonical form. To mitigate this, we explicitly bring the MPS back into + canonical form, if the Frobenius norm of the `site-resolved norm errors array `_ + is greater than `norm_tol`. + + Returns: + None + """ + + # check that sites are adjacent + if abs(sites[0] - sites[1]) != 1: + raise ValueError("sites must be adjacent") + + # check whether to transpose gate + if sites[0] > sites[1]: + U2 = U2.T + + # apply NN gate between (site1, site2) + U2_npc = npc.Array.from_ndarray( + U2, [shfsc, shfsc.conj()], labels=["(p0.p1)", "(p0*.p1*)"] + ) + U2_npc_split = U2_npc.split_legs() + eng.update_bond(max(sites), U2_npc_split) + + # recanonicalize psi if below error threshold + psi = eng.get_resume_data()["psi"] + if np.linalg.norm(psi.norm_test()) > norm_tol: + psi.canonical_form_finite() diff --git a/python/ffsim/tenpy/gates/basic_gates.py b/python/ffsim/tenpy/gates/basic_gates.py new file mode 100644 index 000000000..d88e2224d --- /dev/null +++ b/python/ffsim/tenpy/gates/basic_gates.py @@ -0,0 +1,208 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy basic gates.""" + +import cmath +import math + +import numpy as np +import scipy.linalg +from numpy.typing import NDArray + +from ffsim.spin import Spin + + +def givens_rotation( + theta: float, spin: Spin = Spin.ALPHA_AND_BETA, *, phi: float = 0.0 +) -> NDArray[np.complex128]: + r"""The Givens rotation gate. + + The Givens rotation gate defined in :func:`~ffsim.apply_givens_rotation`, + returned in the TeNPy (N, Sz)-symmetry-conserved basis. + + The bitstring ordering of the TeNPy (N, Sz)-symmetry-conserved basis is: + + .. code:: + + 1010 # (2, -2) + 1000 # (1, -1) + 0010 # (1, -1) + 1011 # (3, -1) + 1110 # (3, -1) + 0000 # (0, 0) + 1001 # (2, 0) + 0011 # (2, 0) + 1100 # (2, 0) + 0110 # (2, 0) + 1111 # (4, 0) + 0001 # (1, 1) + 0100 # (1, 1) + 1101 # (3, 1) + 0111 # (3, 1) + 0101 # (2, 2) + + Args: + theta: The rotation angle. + spin: Choice of spin sector(s) to act on. + + - To act on only spin alpha, pass :const:`ffsim.Spin.ALPHA`. + - To act on only spin beta, pass :const:`ffsim.Spin.BETA`. + - To act on both spin alpha and spin beta, pass + :const:`ffsim.Spin.ALPHA_AND_BETA`. + phi: The phase angle. + + Returns: + The Givens rotation gate in the TeNPy (N, Sz)-symmetry-conserved basis. + """ + c = math.cos(theta) + s = cmath.rect(1, phi) * math.sin(theta) + mat = np.array([[c, s], [-s.conjugate(), c]]) + mat_a = mat if spin & Spin.ALPHA else np.eye(2) + mat_b = mat if spin & Spin.BETA else np.eye(2) + return scipy.linalg.block_diag( + 1, + mat_b, + mat_a.conj(), + 1, + scipy.linalg.block_diag(mat_a.conj(), mat_a.T)[[0, 2, 1, 3]][:, [0, 2, 1, 3]] + @ scipy.linalg.block_diag(mat_b.T.conj(), mat_b), + 1, + mat_a.T, + mat_b.T.conj(), + 1, + ) + + +def num_interaction( + theta: float, spin: Spin = Spin.ALPHA_AND_BETA +) -> NDArray[np.complex128]: + r"""The number interaction gate. + + The number interaction gate defined in :func:`~ffsim.apply_num_interaction`, + returned in the TeNPy (N, Sz)-symmetry-conserved basis. + + The bitstring ordering of the TeNPy (N, Sz)-symmetry-conserved basis is: + + .. code:: + + 10 # (1, -1) + 00 # (0, 0) + 11 # (2, 0) + 01 # (1, 1) + + Args: + theta: The rotation angle. + spin: Choice of spin sector(s) to act on. + + - To act on only spin alpha, pass :const:`ffsim.Spin.ALPHA`. + - To act on only spin beta, pass :const:`ffsim.Spin.BETA`. + - To act on both spin alpha and spin beta, pass + :const:`ffsim.Spin.ALPHA_AND_BETA` (this is the default value). + + Returns: + The number interaction gate in the TeNPy (N, Sz)-symmetry-conserved basis. + """ + phase = cmath.rect(1, theta) + alpha_phase = phase if spin & Spin.ALPHA else 1 + beta_phase = phase if spin & Spin.BETA else 1 + return np.diag([beta_phase, 1, alpha_phase * beta_phase, alpha_phase]) + + +def on_site_interaction(theta: float) -> NDArray[np.complex128]: + r"""The on-site interaction gate. + + The on-site interaction gate defined in :func:`~ffsim.apply_on_site_interaction`, + returned in the TeNPy (N, Sz)-symmetry-conserved basis. + + The bitstring ordering of the TeNPy (N, Sz)-symmetry-conserved basis is: + + .. code:: + + 10 # (1, -1) + 00 # (0, 0) + 11 # (2, 0) + 01 # (1, 1) + + Args: + theta: The rotation angle. + + Returns: + The on-site interaction gate in the TeNPy (N, Sz)-symmetry-conserved basis. + """ + return np.diag([1, 1, cmath.rect(1, theta), 1]) + + +def num_num_interaction( + theta: float, spin: Spin = Spin.ALPHA_AND_BETA +) -> NDArray[np.complex128]: + r"""The number-number interaction gate. + + The number-number interaction gate defined in + :func:`~ffsim.apply_num_num_interaction`, returned in the TeNPy + (N, Sz)-symmetry-conserved basis. + + The bitstring ordering of the TeNPy (N, Sz)-symmetry-conserved basis is: + + .. code:: + + 1010 # (2, -2) + 1000 # (1, -1) + 0010 # (1, -1) + 1011 # (3, -1) + 1110 # (3, -1) + 0000 # (0, 0) + 1001 # (2, 0) + 0011 # (2, 0) + 1100 # (2, 0) + 0110 # (2, 0) + 1111 # (4, 0) + 0001 # (1, 1) + 0100 # (1, 1) + 1101 # (3, 1) + 0111 # (3, 1) + 0101 # (2, 2) + + Args: + theta: The rotation angle. + spin: Choice of spin sector(s) to act on. + + - To act on only spin alpha, pass :const:`ffsim.Spin.ALPHA`. + - To act on only spin beta, pass :const:`ffsim.Spin.BETA`. + - To act on both spin alpha and spin beta, pass + :const:`ffsim.Spin.ALPHA_AND_BETA` (this is the default value). + + Returns: + The number-number interaction gate in the TeNPy (N, Sz)-symmetry-conserved + basis. + """ + phase = cmath.rect(1, theta) + alpha_phase = phase if spin & Spin.ALPHA else 1 + beta_phase = phase if spin & Spin.BETA else 1 + return np.diag( + [ + beta_phase, + 1, + 1, + beta_phase, + beta_phase, + 1, + 1, + 1, + 1, + 1, + alpha_phase * beta_phase, + 1, + 1, + alpha_phase, + alpha_phase, + alpha_phase, + ] + ) diff --git a/python/ffsim/tenpy/gates/diag_coulomb.py b/python/ffsim/tenpy/gates/diag_coulomb.py new file mode 100644 index 000000000..db085783e --- /dev/null +++ b/python/ffsim/tenpy/gates/diag_coulomb.py @@ -0,0 +1,67 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy diagonal Coulomb evolution gate.""" + +import itertools + +import numpy as np +from numpy.typing import NDArray +from tenpy.algorithms.tebd import TEBDEngine + +from ffsim.tenpy.gates.abstract_gates import apply_single_site, apply_two_site +from ffsim.tenpy.gates.basic_gates import num_num_interaction, on_site_interaction + + +def apply_diag_coulomb_evolution( + eng: TEBDEngine, + mat: NDArray[np.complex128], + time: float, + *, + norm_tol: float = 1e-8, +) -> None: + r"""Apply a diagonal Coulomb evolution gate to an MPS. + + The diagonal Coulomb evolution gate is defined in + :func:`~ffsim.apply_diag_coulomb_evolution`. + + Args: + eng: The TEBD engine. + mat: The diagonal Coulomb matrices of dimension `(2, norb, norb)`. + time: The evolution time. + norm_tol: The norm error above which we recanonicalize the MPS. In general, the + application of a two-site gate to an MPS with truncation may degrade its + canonical form. To mitigate this, we explicitly bring the MPS back into + canonical form, if the Frobenius norm of the `site-resolved norm errors array `_ + is greater than `norm_tol`. + + Returns: + None + """ + + # extract norb + norb = eng.get_resume_data()["psi"].L + + # unpack alpha-alpha and alpha-beta matrices + mat_aa, mat_ab = mat + + # apply alpha-alpha gates + for i, j in itertools.combinations(range(norb), 2): + if mat_aa[i, j]: + apply_two_site( + eng, + num_num_interaction(-time * mat_aa[i, j]), + (i, j), + norm_tol=norm_tol, + ) + + # apply alpha-beta gates + for i in range(norb): + apply_single_site(eng, on_site_interaction(-time * mat_ab[i, i]), i) diff --git a/python/ffsim/tenpy/gates/orbital_rotation.py b/python/ffsim/tenpy/gates/orbital_rotation.py new file mode 100644 index 000000000..d041c2f3d --- /dev/null +++ b/python/ffsim/tenpy/gates/orbital_rotation.py @@ -0,0 +1,65 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy orbital rotation gate.""" + +import cmath +import math + +import numpy as np +from numpy.typing import NDArray +from tenpy.algorithms.tebd import TEBDEngine + +from ffsim.linalg import givens_decomposition +from ffsim.tenpy.gates.abstract_gates import apply_single_site, apply_two_site +from ffsim.tenpy.gates.basic_gates import givens_rotation, num_interaction + + +def apply_orbital_rotation( + eng: TEBDEngine, + mat: NDArray[np.complex128], + *, + norm_tol: float = 1e-8, +) -> None: + r"""Apply an orbital rotation gate to an MPS. + + The orbital rotation gate is defined in :func:`~ffsim.apply_orbital_rotation`. + + Args: + eng: The TEBD engine. + mat: The orbital rotation matrix of dimension `(norb, norb)`. + norm_tol: The norm error above which we recanonicalize the MPS. In general, the + application of a two-site gate to an MPS with truncation may degrade its + canonical form. To mitigate this, we explicitly bring the MPS back into + canonical form, if the Frobenius norm of the `site-resolved norm errors array `_ + is greater than `norm_tol`. + + Returns: + None + """ + + # Givens decomposition + givens_list, diag_mat = givens_decomposition(mat) + + # apply the Givens rotation gates + for gate in givens_list: + theta = math.acos(gate.c) + phi = -cmath.phase(gate.s) + apply_two_site( + eng, + givens_rotation(theta, phi=phi), + (gate.i, gate.j), + norm_tol=norm_tol, + ) + + # apply the number interaction gates + for i, z in enumerate(diag_mat): + theta = cmath.phase(z) + apply_single_site(eng, num_interaction(theta), i) diff --git a/python/ffsim/tenpy/gates/ucj.py b/python/ffsim/tenpy/gates/ucj.py new file mode 100644 index 000000000..f7189e12b --- /dev/null +++ b/python/ffsim/tenpy/gates/ucj.py @@ -0,0 +1,71 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy unitary cluster Jastrow gate.""" + +from __future__ import annotations + +import numpy as np +from tenpy.algorithms.tebd import TEBDEngine + +from ffsim.tenpy.gates.diag_coulomb import apply_diag_coulomb_evolution +from ffsim.tenpy.gates.orbital_rotation import apply_orbital_rotation +from ffsim.variational.ucj_spin_balanced import UCJOpSpinBalanced + + +def apply_ucj_op_spin_balanced( + eng: TEBDEngine, + ucj_op: UCJOpSpinBalanced, + *, + norm_tol: float = 1e-8, +) -> None: + r"""Apply a spin-balanced unitary cluster Jastrow gate to an MPS. + + The spin-balanced unitary cluster Jastrow gate is defined in + :class:`~ffsim.variational.ucj_spin_balanced.UCJOpSpinBalanced`. + + Args: + eng: The TEBD engine. + ucj_op: The spin-balanced unitary cluster Jastrow operator. + norm_tol: The norm error above which we recanonicalize the MPS. In general, the + application of a two-site gate to an MPS with truncation may degrade its + canonical form. To mitigate this, we explicitly bring the MPS back into + canonical form, if the Frobenius norm of the `site-resolved norm errors array `_ + is greater than `norm_tol`. + + Returns: + None + """ + + # extract norb + norb = eng.get_resume_data()["psi"].L + + # construct the LUCJ MPS + current_basis = np.eye(norb, dtype=complex) + for orb_rot, diag_mats in zip(ucj_op.orbital_rotations, ucj_op.diag_coulomb_mats): + apply_orbital_rotation( + eng, + orb_rot.conjugate().T @ current_basis, + norm_tol=norm_tol, + ) + apply_diag_coulomb_evolution(eng, diag_mats, -1, norm_tol=norm_tol) + current_basis = orb_rot + if ucj_op.final_orbital_rotation is None: + apply_orbital_rotation( + eng, + current_basis, + norm_tol=norm_tol, + ) + else: + apply_orbital_rotation( + eng, + ucj_op.final_orbital_rotation @ current_basis, + norm_tol=norm_tol, + ) diff --git a/python/ffsim/tenpy/random/__init__.py b/python/ffsim/tenpy/random/__init__.py new file mode 100644 index 000000000..7918d6a4b --- /dev/null +++ b/python/ffsim/tenpy/random/__init__.py @@ -0,0 +1,9 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/python/ffsim/tenpy/random/random.py b/python/ffsim/tenpy/random/random.py new file mode 100644 index 000000000..eff665984 --- /dev/null +++ b/python/ffsim/tenpy/random/random.py @@ -0,0 +1,76 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +import random + +from tenpy.algorithms.tebd import RandomUnitaryEvolution +from tenpy.networks.mps import MPS +from tenpy.networks.site import SpinHalfFermionSite + +from ffsim.tenpy.util import _bitstring_to_product_state + + +def random_mps( + norb: int, nelec: tuple[int, int], n_steps: int = 10, chi_max: int = 100 +) -> MPS: + """Return a random MPS generated from a random unitary evolution. + + Args: + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + n_steps: The number of steps in the random unitary evolution. + chi_max: The maximum bond dimension in the random unitary evolution. + + Returns: + The random MPS. + """ + + # initialize Hartree-Fock state + n_alpha, n_beta = nelec + product_state = _bitstring_to_product_state( + ((1 << n_alpha) - 1, (1 << n_beta) - 1), norb + ) + shfs = SpinHalfFermionSite(cons_N="N", cons_Sz="Sz") + mps = MPS.from_product_state([shfs] * norb, product_state) + + # apply random unitary evolution + options = {"N_steps": n_steps, "trunc_params": {"chi_max": chi_max}} + eng = RandomUnitaryEvolution(mps, options) + eng.run() + mps.canonical_form() + + return mps + + +def random_mps_product_state(norb: int, nelec: tuple[int, int]) -> MPS: + """Return a random MPS product state. + + Args: + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + + Returns: + The random MPS product state. + """ + n_alpha, n_beta = nelec + + n_alpha_list = [1] * n_alpha + [0] * (norb - n_alpha) + random.shuffle(n_alpha_list) + n_beta_list = [1] * n_beta + [0] * (norb - n_beta) + random.shuffle(n_beta_list) + + s_alpha = sum(j << i for i, j in enumerate(n_alpha_list)) + s_beta = sum(j << i for i, j in enumerate(n_beta_list)) + product_state = _bitstring_to_product_state((s_alpha, s_beta), norb) + + shfs = SpinHalfFermionSite(cons_N="N", cons_Sz="Sz") + mps = MPS.from_product_state([shfs] * norb, product_state) + + return mps diff --git a/python/ffsim/tenpy/util.py b/python/ffsim/tenpy/util.py new file mode 100644 index 000000000..95cf4e18a --- /dev/null +++ b/python/ffsim/tenpy/util.py @@ -0,0 +1,253 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""TeNPy utility functions.""" + +from __future__ import annotations + +from typing import cast + +import numpy as np +import tenpy.linalg.np_conserved as npc +from numpy.typing import NDArray +from tenpy.algorithms.exact_diag import ExactDiag +from tenpy.models.hubbard import FermiHubbardChain +from tenpy.networks.mps import MPS +from tenpy.networks.site import FermionSite, SpinHalfFermionSite + +import ffsim + + +def mps_to_statevector(mps: MPS) -> NDArray[np.complex128]: + r"""Return the MPS as a state vector. + + Args: + mps: The MPS. + + Returns: + The state vector. + """ + + # generate the ffsim-ordered list of product states + norb = mps.L + n_alpha = round(np.sum(mps.expectation_value("Nu"))) + n_beta = round(np.sum(mps.expectation_value("Nd"))) + product_states = _generate_product_states(norb, (n_alpha, n_beta)) + + # initialize the TeNPy ExactDiag class instance + charge_sector = mps.get_total_charge(True) + exact_diag = ExactDiag(FermiHubbardChain({"L": norb}), charge_sector=charge_sector) + + # determine the mapping from TeNPy basis to ffsim basis + basis_ordering_ffsim, swap_factors_ffsim = _map_tenpy_to_ffsim_basis( + product_states, exact_diag + ) + + # convert TeNPy MPS to ffsim statevector + statevector = cast(NDArray[np.complex128], exact_diag.mps_to_full(mps).to_ndarray()) + statevector = np.multiply(swap_factors_ffsim, statevector[basis_ordering_ffsim]) + + return statevector + + +def statevector_to_mps( + statevector: NDArray[np.complex128], + norb: int, + nelec: tuple[int, int], +) -> MPS: + r"""Return the state vector as an MPS. + + Args: + statevector: The state vector. + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + + Returns: + The MPS. + """ + + # check if state vector is basis state + basis_state = np.count_nonzero(statevector) == 1 + + # generate the ffsim-ordered list of product states + if basis_state: + idx = int(np.flatnonzero(statevector)[0]) + string = ffsim.addresses_to_strings( + [idx], + norb, + nelec, + concatenate=False, + bitstring_type=ffsim.BitstringType.INT, + ) + bitstring = (string[0][0], string[1][0]) + product_states = [_bitstring_to_product_state(bitstring, norb)] + else: + product_states = _generate_product_states(norb, nelec) + + # construct the reference product state MPS + shfs = SpinHalfFermionSite(cons_N="N", cons_Sz="Sz") + mps_reference = MPS.from_product_state([shfs] * norb, product_states[0]) + + if basis_state: + # compute swap factor + swap_factor = _compute_swap_factor(mps_reference) + + # apply swap factor + if swap_factor == -1: + minus_identity_npc = npc.Array.from_ndarray( + -shfs.get_op("Id").to_ndarray(), + [shfs.leg, shfs.leg.conj()], + labels=["p", "p*"], + ) + mps_reference.apply_local_op(0, minus_identity_npc) + + mps = mps_reference + else: + # initialize the TeNPy ExactDiag class instance + charge_sector = mps_reference.get_total_charge(True) + exact_diag = ExactDiag( + FermiHubbardChain({"L": norb}), charge_sector=charge_sector + ) + statevector_reference = exact_diag.mps_to_full(mps_reference) + leg_charge = statevector_reference.legs[0] + + # determine the mapping from ffsim basis to TeNPy basis + basis_ordering_ffsim, swap_factors_ffsim = _map_tenpy_to_ffsim_basis( + product_states, exact_diag + ) + basis_ordering_tenpy = np.argsort(basis_ordering_ffsim) + swap_factors_tenpy = swap_factors_ffsim[np.argsort(basis_ordering_ffsim)] + + # convert ffsim statevector to TeNPy MPS + statevector = np.multiply(swap_factors_tenpy, statevector[basis_ordering_tenpy]) + statevector_npc = npc.Array.from_ndarray(statevector, [leg_charge]) + mps = exact_diag.full_to_mps(statevector_npc) + + return mps + + +def _bitstring_to_product_state(bitstring: tuple[int, int], norb: int) -> list[int]: + r"""Return the product state in TeNPy SpinHalfFermionSite notation. + + Args: + bitstring: The bitstring in the form `(int_a, int_b)`. + norb: The number of spatial orbitals. + + Returns: + The product state in TeNPy SpinHalfFermionSite notation. + """ + + # unpack bitstrings + int_a, int_b = bitstring + string_a = format(int_a, f"0{norb}b") + string_b = format(int_b, f"0{norb}b") + + # relabel using TeNPy SpinHalfFermionSite notation + product_state = [] + for i, site in enumerate(zip(reversed(string_b), reversed(string_a))): + site_occupation = int("".join(site), base=2) + product_state.append(site_occupation) + + return product_state + + +def _generate_product_states(norb: int, nelec: tuple[int, int]) -> list[list[int]]: + r"""Return the ffsim-ordered list of product states in TeNPy SpinHalfFermionSite + notation. + + Args: + norb: The number of spatial orbitals. + nelec: The number of alpha and beta electrons. + + Returns: + The ffsim-ordered list of product states in TeNPy SpinHalfFermionSite notation. + """ + + # generate the strings + dim = ffsim.dim(norb, nelec) + strings = ffsim.addresses_to_strings( + range(dim), norb=norb, nelec=nelec, bitstring_type=ffsim.BitstringType.STRING + ) + string_tuples = [ + ( + int(string[len(string) // 2 :], base=2), + int(string[: len(string) // 2], base=2), + ) + for string in strings + ] + + # convert strings to product states + product_states = [] + for bitstring in string_tuples: + product_states.append(_bitstring_to_product_state(bitstring, norb)) + + return product_states + + +def _compute_swap_factor(mps: MPS) -> int: + r"""Compute the swap factor due to the conversion from TeNPy to ffsim bases. + + Args: + mps: The MPS. + + Returns: + The swap factor (+1 or -1). + """ + + norb = mps.L + fs = FermionSite(conserve="N") + alpha_sector = mps.expectation_value("Nu") + beta_sector = mps.expectation_value("Nd") + product_state_fs_tenpy = [ + int(val) for pair in zip(alpha_sector, beta_sector) for val in pair + ] + mps_fs = MPS.from_product_state([fs] * 2 * norb, product_state_fs_tenpy) + + tenpy_ordering = list(range(2 * norb)) + midpoint = len(tenpy_ordering) // 2 + mask1 = tenpy_ordering[:midpoint][::-1] + mask2 = tenpy_ordering[midpoint:][::-1] + ffsim_ordering = [int(val) for pair in zip(mask1, mask2) for val in pair] + + mps_ref = mps_fs.copy() + mps_ref.permute_sites(ffsim_ordering, swap_op=None) + mps_fs.permute_sites(ffsim_ordering, swap_op="auto") + swap_factor = cast(int, round(mps_fs.overlap(mps_ref))) + + return swap_factor + + +def _map_tenpy_to_ffsim_basis( + product_states: list[list[int]], exact_diag: ExactDiag +) -> tuple[NDArray[np.int_], NDArray[np.int_]]: + r"""Return the mapping from the TeNPy basis to the ffsim basis. + + Args: + product_states: The ffsim-ordered list of product states in TeNPy + SpinHalfFermionSite notation. + exact_diag: The exact diagonalization class instance. + + Returns: + basis_ordering_ffsim: The permutation to map from the TeNPy to ffsim basis. + swap_factors: The minus signs that are introduced due to this mapping. + """ + + basis_ordering_ffsim, swap_factors = [], [] + for i, state in enumerate(product_states): + # basis_ordering_ffsim + prod_mps = MPS.from_product_state(exact_diag.model.lat.mps_sites(), state) + prod_statevector = list(exact_diag.mps_to_full(prod_mps).to_ndarray()) + idx = prod_statevector.index(1) + basis_ordering_ffsim.append(idx) + + # swap_factors + swap_factors.append(_compute_swap_factor(prod_mps)) + + return np.array(basis_ordering_ffsim), np.array(swap_factors) diff --git a/tests/python/tenpy/__init__.py b/tests/python/tenpy/__init__.py new file mode 100644 index 000000000..7918d6a4b --- /dev/null +++ b/tests/python/tenpy/__init__.py @@ -0,0 +1,9 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/tests/python/tenpy/gates/__init__.py b/tests/python/tenpy/gates/__init__.py new file mode 100644 index 000000000..7918d6a4b --- /dev/null +++ b/tests/python/tenpy/gates/__init__.py @@ -0,0 +1,9 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/tests/python/tenpy/gates/basic_gates_test.py b/tests/python/tenpy/gates/basic_gates_test.py new file mode 100644 index 000000000..36c764db0 --- /dev/null +++ b/tests/python/tenpy/gates/basic_gates_test.py @@ -0,0 +1,268 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the TeNPy basic gates.""" + +from copy import deepcopy + +import numpy as np +import pytest +from tenpy.algorithms.tebd import TEBDEngine +from tenpy.models.molecular import MolecularModel + +import ffsim +from ffsim.spin import Spin +from ffsim.tenpy.gates.basic_gates import ( + givens_rotation, + num_interaction, + num_num_interaction, + on_site_interaction, +) +from ffsim.tenpy.util import statevector_to_mps + + +@pytest.mark.parametrize( + "norb, nelec, spin", + [ + (3, (0, 0), Spin.ALPHA), + (3, (0, 1), Spin.ALPHA), + (3, (1, 2), Spin.ALPHA), + (3, (2, 2), Spin.ALPHA), + (3, (0, 0), Spin.BETA), + (3, (0, 1), Spin.BETA), + (3, (1, 2), Spin.BETA), + (3, (2, 2), Spin.BETA), + (3, (0, 0), Spin.ALPHA_AND_BETA), + (3, (0, 1), Spin.ALPHA_AND_BETA), + (3, (1, 2), Spin.ALPHA_AND_BETA), + (3, (2, 2), Spin.ALPHA_AND_BETA), + ], +) +def test_givens_rotation(norb: int, nelec: tuple[int, int], spin: Spin): + """Test applying a Givens rotation gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + linop = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mpo_model = MolecularModel(model_params) + mpo = mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate random Givens rotation parameters + theta = rng.uniform(0, 2 * np.pi) + phi = rng.uniform(0, 2 * np.pi) + p = rng.integers(norb - 1) + + # apply random Givens rotation to state vector + vec = ffsim.apply_givens_rotation( + original_vec, theta, (p, p + 1), norb, nelec, spin, phi=phi + ) + + # apply random orbital rotation to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_two_site(eng, givens_rotation(theta, spin, phi=phi), (p, p + 1)) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, linop @ vec) + mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec, spin", + [ + (3, (2, 2), Spin.ALPHA), + (3, (1, 2), Spin.ALPHA), + (3, (0, 2), Spin.ALPHA), + (3, (0, 0), Spin.ALPHA), + (3, (2, 2), Spin.BETA), + (3, (1, 2), Spin.BETA), + (3, (0, 2), Spin.BETA), + (3, (0, 0), Spin.BETA), + (3, (2, 2), Spin.ALPHA_AND_BETA), + (3, (1, 2), Spin.ALPHA_AND_BETA), + (3, (0, 2), Spin.ALPHA_AND_BETA), + (3, (0, 0), Spin.ALPHA_AND_BETA), + ], +) +def test_num_interaction(norb: int, nelec: tuple[int, int], spin: Spin): + """Test applying a number interaction gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate random number interaction parameters + theta = 2 * np.pi * rng.random() + p = rng.integers(0, norb) + + # apply random number interaction to state vector + vec = ffsim.apply_num_interaction(original_vec, theta, p, norb, nelec, spin) + + # apply random number interaction to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_single_site(eng, num_interaction(theta, spin), p) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, hamiltonian @ vec) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (3, (2, 2)), + (3, (1, 2)), + (3, (0, 2)), + (3, (0, 0)), + ], +) +def test_on_site_interaction( + norb: int, + nelec: tuple[int, int], +): + """Test applying an on-site interaction gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate random on-site interaction parameters + theta = 2 * np.pi * rng.random() + p = rng.integers(0, norb) + + # apply random on-site interaction to state vector + vec = ffsim.apply_on_site_interaction(original_vec, theta, p, norb, nelec) + + # apply random on-site interaction to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_single_site(eng, on_site_interaction(theta), p) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, hamiltonian @ vec) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec, spin", + [ + (3, (2, 2), Spin.ALPHA), + (3, (1, 2), Spin.ALPHA), + (3, (0, 2), Spin.ALPHA), + (3, (0, 0), Spin.ALPHA), + (3, (2, 2), Spin.BETA), + (3, (1, 2), Spin.BETA), + (3, (0, 2), Spin.BETA), + (3, (0, 0), Spin.BETA), + (3, (2, 2), Spin.ALPHA_AND_BETA), + (3, (1, 2), Spin.ALPHA_AND_BETA), + (3, (0, 2), Spin.ALPHA_AND_BETA), + (3, (0, 0), Spin.ALPHA_AND_BETA), + ], +) +def test_num_num_interaction(norb: int, nelec: tuple[int, int], spin: Spin): + """Test applying a number-number interaction gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate random number-number interaction parameters + theta = 2 * np.pi * rng.random() + p = rng.integers(0, norb - 1) + + # apply random number-number interaction to state vector + vec = ffsim.apply_num_num_interaction( + original_vec, theta, (p, p + 1), norb, nelec, spin + ) + + # apply random number-number interaction to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_two_site(eng, num_num_interaction(theta, spin), (p, p + 1)) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, hamiltonian @ vec) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) diff --git a/tests/python/tenpy/gates/diag_coulomb_test.py b/tests/python/tenpy/gates/diag_coulomb_test.py new file mode 100644 index 000000000..bac023726 --- /dev/null +++ b/tests/python/tenpy/gates/diag_coulomb_test.py @@ -0,0 +1,78 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the TeNPy diagonal Coulomb evolution gate.""" + +from copy import deepcopy + +import numpy as np +import pytest +from tenpy.algorithms.tebd import TEBDEngine +from tenpy.models.molecular import MolecularModel + +import ffsim +from ffsim.tenpy.util import statevector_to_mps + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (3, (2, 2)), + (3, (1, 2)), + (3, (0, 2)), + (3, (0, 0)), + ], +) +def test_apply_diag_coulomb_evolution(norb: int, nelec: tuple[int, int]): + """Test applying a diagonal Coulomb evolution gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate random diagonal Coulomb evolution parameters + mat_aa = np.diag(rng.standard_normal(norb - 1), k=-1) + mat_aa += mat_aa.T + mat_ab = np.diag(rng.standard_normal(norb)) + diag_coulomb_mats = np.array([mat_aa, mat_ab, mat_aa]) + time = rng.random() + + # apply random diagonal Coulomb evolution to state vector + vec = ffsim.apply_diag_coulomb_evolution( + original_vec, diag_coulomb_mats, time, norb, nelec + ) + + # apply random diagonal Coulomb evolution to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_diag_coulomb_evolution(eng, diag_coulomb_mats[:2], time) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, hamiltonian @ vec) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) diff --git a/tests/python/tenpy/gates/orbital_rotation_test.py b/tests/python/tenpy/gates/orbital_rotation_test.py new file mode 100644 index 000000000..fb6447a8c --- /dev/null +++ b/tests/python/tenpy/gates/orbital_rotation_test.py @@ -0,0 +1,75 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the TeNPy orbital rotation gate.""" + +from copy import deepcopy + +import numpy as np +import pytest +from tenpy.algorithms.tebd import TEBDEngine +from tenpy.models.molecular import MolecularModel + +import ffsim +from ffsim.tenpy.util import statevector_to_mps + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (3, (2, 2)), + (3, (1, 2)), + (3, (0, 2)), + (3, (0, 0)), + ], +) +def test_apply_orbital_rotation( + norb: int, + nelec: tuple[int, int], +): + """Test applying an orbital rotation gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + original_vec = ffsim.random.random_state_vector(dim, seed=rng) + + # convert random state vector to MPS + mps = statevector_to_mps(original_vec, norb, nelec) + original_mps = deepcopy(mps) + + # generate a random orbital rotation + mat = ffsim.random.random_unitary(norb, seed=rng) + + # apply random orbital rotation to state vector + vec = ffsim.apply_orbital_rotation(original_vec, mat, norb, nelec) + + # apply random orbital rotation to MPS + eng = TEBDEngine(mps, None, {}) + ffsim.tenpy.apply_orbital_rotation(eng, mat) + + # test expectation is preserved + original_expectation = np.vdot(original_vec, hamiltonian @ vec) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = original_mps.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) diff --git a/tests/python/tenpy/gates/ucj_test.py b/tests/python/tenpy/gates/ucj_test.py new file mode 100644 index 000000000..4d7634805 --- /dev/null +++ b/tests/python/tenpy/gates/ucj_test.py @@ -0,0 +1,97 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the TeNPy unitary cluster Jastrow gate.""" + +import numpy as np +import pytest +from tenpy.algorithms.tebd import TEBDEngine +from tenpy.models.molecular import MolecularModel + +import ffsim +from ffsim.tenpy.gates.ucj import apply_ucj_op_spin_balanced +from ffsim.tenpy.util import statevector_to_mps +from ffsim.variational.util import interaction_pairs_spin_balanced + + +@pytest.mark.parametrize( + "norb, nelec, n_reps, connectivity", + [ + (3, (2, 2), 1, "square"), + (3, (1, 2), 1, "square"), + (3, (0, 2), 1, "square"), + (3, (0, 0), 1, "square"), + (3, (2, 2), 1, "hex"), + (3, (1, 2), 1, "hex"), + (3, (0, 2), 1, "hex"), + (3, (0, 0), 1, "hex"), + (3, (2, 2), 1, "heavy-hex"), + (3, (1, 2), 1, "heavy-hex"), + (3, (0, 2), 1, "heavy-hex"), + (3, (0, 0), 1, "heavy-hex"), + (3, (2, 2), 2, "square"), + (3, (1, 2), 2, "square"), + (3, (0, 2), 2, "square"), + (3, (0, 0), 2, "square"), + (3, (2, 2), 2, "hex"), + (3, (1, 2), 2, "hex"), + (3, (0, 2), 2, "hex"), + (3, (0, 0), 2, "hex"), + (3, (2, 2), 2, "heavy-hex"), + (3, (1, 2), 2, "heavy-hex"), + (3, (0, 2), 2, "heavy-hex"), + (3, (0, 0), 2, "heavy-hex"), + ], +) +def test_apply_ucj_op_spin_balanced( + norb: int, nelec: tuple[int, int], n_reps: int, connectivity: str +): + """Test applying a spin-balanced unitary cluster Jastrow gate to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random LUCJ ansatz + lucj_op = ffsim.random.random_ucj_op_spin_balanced( + norb=norb, + n_reps=n_reps, + interaction_pairs=interaction_pairs_spin_balanced( + connectivity=connectivity, norb=norb + ), + with_final_orbital_rotation=True, + seed=rng, + ) + + # generate the corresponding LUCJ circuit statevector + hf_state = ffsim.hartree_fock_state(norb, nelec) + lucj_state = ffsim.apply_unitary(hf_state, lucj_op, norb, nelec) + + # generate the corresponding LUCJ circuit MPS + dim = ffsim.dim(norb, nelec) + wavefunction_mps = statevector_to_mps(np.array([1] + [0] * (dim - 1)), norb, nelec) + options = {"trunc_params": {"chi_max": 16, "svd_min": 1e-6}} + eng = TEBDEngine(wavefunction_mps, None, options) + apply_ucj_op_spin_balanced(eng, lucj_op) + + # test expectation is preserved + original_expectation = np.vdot(lucj_state, hamiltonian @ lucj_state).real + mpo_expectation = mol_hamiltonian_mpo.expectation_value_finite(wavefunction_mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) diff --git a/tests/python/tenpy/util_test.py b/tests/python/tenpy/util_test.py new file mode 100644 index 000000000..203a4c532 --- /dev/null +++ b/tests/python/tenpy/util_test.py @@ -0,0 +1,194 @@ +# (C) Copyright IBM 2025. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from copy import deepcopy + +import numpy as np +import pytest +from tenpy.models.molecular import MolecularModel + +import ffsim +from ffsim.tenpy.random.random import random_mps, random_mps_product_state +from ffsim.tenpy.util import mps_to_statevector, statevector_to_mps + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (2, (2, 2)), + (2, (2, 1)), + (2, (1, 2)), + (2, (1, 1)), + (2, (0, 2)), + (2, (0, 0)), + (3, (2, 2)), + ], +) +def test_mps_to_statevector_product_state(norb: int, nelec: tuple[int, int]): + """Test converting an MPS to a statevector using a product state.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random MPS (product state) + mps = random_mps_product_state(norb, nelec) + + # convert MPS to state vector (product state) + statevector = mps_to_statevector(mps) + + # test expectation is preserved + original_expectation = np.vdot(statevector, hamiltonian @ statevector) + mps_original = deepcopy(mps) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = mps_original.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (2, (2, 2)), + (2, (2, 1)), + (2, (1, 2)), + (2, (1, 1)), + (2, (0, 2)), + (2, (0, 0)), + (3, (2, 2)), + ], +) +def test_mps_to_statevector(norb: int, nelec: tuple[int, int]): + """Test converting an MPS to a state vector.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random MPS + mps = random_mps(norb, nelec) + + # convert MPS to state vector + statevector = mps_to_statevector(mps) + + # test expectation is preserved + original_expectation = np.vdot(statevector, hamiltonian @ statevector) + mps_original = deepcopy(mps) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = mps_original.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (2, (2, 2)), + (2, (2, 1)), + (2, (1, 2)), + (2, (1, 1)), + (2, (0, 2)), + (2, (0, 0)), + (3, (2, 2)), + ], +) +def test_statevector_to_mps_product_state(norb: int, nelec: tuple[int, int]): + """Test converting a state vector to an MPS using a product state.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector (product state) + dim = ffsim.dim(norb, nelec) + idx = rng.integers(0, high=dim) + statevector = ffsim.linalg.one_hot(dim, idx) + + # convert state vector to MPS (product state) + mps = statevector_to_mps(statevector, norb, nelec) + + # test expectation is preserved + original_expectation = np.vdot(statevector, hamiltonian @ statevector) + mps_original = deepcopy(mps) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = mps_original.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation) + + +@pytest.mark.parametrize( + "norb, nelec", + [ + (2, (2, 2)), + (2, (2, 1)), + (2, (1, 2)), + (2, (1, 1)), + (2, (0, 2)), + (2, (0, 0)), + (3, (2, 2)), + ], +) +def test_statevector_to_mps(norb: int, nelec: tuple[int, int]): + """Test converting a state vector to an MPS.""" + rng = np.random.default_rng() + + # generate a random molecular Hamiltonian + mol_hamiltonian = ffsim.random.random_molecular_hamiltonian(norb, seed=rng) + hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb, nelec) + + # convert molecular Hamiltonian to MPO + model_params = dict( + one_body_tensor=mol_hamiltonian.one_body_tensor, + two_body_tensor=mol_hamiltonian.two_body_tensor, + constant=mol_hamiltonian.constant, + ) + mol_hamiltonian_mpo_model = MolecularModel(model_params) + mol_hamiltonian_mpo = mol_hamiltonian_mpo_model.H_MPO + + # generate a random state vector + dim = ffsim.dim(norb, nelec) + statevector = ffsim.random.random_state_vector(dim, seed=rng) + + # convert state vector to MPS + mps = statevector_to_mps(statevector, norb, nelec) + + # test expectation is preserved + original_expectation = np.vdot(statevector, hamiltonian @ statevector) + mps_original = deepcopy(mps) + mol_hamiltonian_mpo.apply_naively(mps) + mpo_expectation = mps_original.overlap(mps) + np.testing.assert_allclose(original_expectation, mpo_expectation)