diff --git a/jitterbug/model/linear.py b/jitterbug/model/linear.py index 294ce58..a39638c 100644 --- a/jitterbug/model/linear.py +++ b/jitterbug/model/linear.py @@ -12,34 +12,48 @@ from .base import EnergyModel -def get_displacement_matrix(atoms: Atoms, reference: Atoms) -> np.ndarray: - """Get the displacements of a structure from a reference - in the order, as used in a Hessian calculation. +def get_model_inputs(atoms: Atoms, reference: Atoms) -> np.ndarray: + """Get the inputs for the model, which are derived from the displacements + of the structure with respect to a reference. Args: atoms: Displaced structure reference: Reference structure Returns: - Vector of displacements + Vector of displacements in the same order as the """ - # Compute the displacements - disp_matrix = (reference.positions - atoms.positions).flatten() - disp_matrix = disp_matrix[:, None] * disp_matrix[None, :] + # Compute the displacements and the products of displacement + disp_matrix = (atoms.positions - reference.positions).flatten() + disp_prod_matrix = disp_matrix[:, None] * disp_matrix[None, :] # Multiply the off-axis terms by two, as they appear twice in the energy model n_terms = len(atoms) * 3 off_diag = np.triu_indices(n_terms, k=1) - disp_matrix[off_diag] *= 2 + disp_prod_matrix[off_diag] *= 2 - # Return the upper triangular matrix - return disp_matrix[np.triu_indices(n_terms)] + # Append the displacements and products of displacements + return np.concatenate([ + disp_matrix, + disp_prod_matrix[np.triu_indices(n_terms)] / 2 + ], axis=0) -class LinearHessianModel(EnergyModel): - """Fits a model for energy using linear regression +class HarmonicModel(EnergyModel): + """Expresses energy as a Harmonic model (i.e., 2nd degree Taylor series) - Implicitly treats all elements of the Hessian matrix as unrelated + Contains a total of :math:`3N + 3N(3N+1)/2` terms in total, where :math:`N` + is the number of atoms in the molecule. The first :math:`3N` correspond to the + linear terms of the model, which are known as the Jacobian matrix, and the + latter are from the quadratic terms, which are half of the symmetric Hessian matrix. + + Implicitly treats all terms of the model as unrelated, which is the worst case + for trying to fit the energy of a molecule. However, it is still possible to fit + the model with a reduced number of terms if we assume that most terms are near zero. + + The energy model is: + + :math:`E = E_0 + \\sum_i J_i \\delta_i + \\frac{1}{2}\\sum_{i,j} H_{i,j}\\delta_i\\delta_j` Args: reference: Fully-relaxed structure used as the reference @@ -50,10 +64,9 @@ def __init__(self, reference: Atoms, regressor: type[LinearModel] = ARDRegressio self.reference = reference self.regressor = regressor - def train(self, data: list[Atoms]) -> ARDRegression: + def train(self, data: list[Atoms]) -> LinearModel: # X: Displacement vectors for each - x = [get_displacement_matrix(atoms, self.reference) for atoms in data] - x = np.multiply(x, 0.5) + x = [get_model_inputs(atoms, self.reference) for atoms in data] # Y: Subtract off the reference energy ref_energy = self.reference.get_potential_energy() @@ -68,12 +81,28 @@ def train(self, data: list[Atoms]) -> ARDRegression: return model - def mean_hessian(self, model: ARDRegression) -> np.ndarray: + def mean_hessian(self, model: LinearModel) -> np.ndarray: return self._params_to_hessian(model.coef_) - def sample_hessians(self, model: ARDRegression, num_samples: int) -> list[np.ndarray]: + def sample_hessians(self, model: LinearModel, num_samples: int) -> list[np.ndarray]: + # Get the covariance matrix + if not hasattr(model, 'sigma_'): # pragma: no-coverage + raise ValueError(f'Sampling only possible with Bayesian regressors. You trained a {type(model)}') + if isinstance(model, ARDRegression): + # The sigma matrix may be zero for high-precision terms + n_terms = len(model.coef_) + nonzero_terms = model.lambda_ < model.threshold_lambda + + # Replace those terms (Thanks: https://stackoverflow.com/a/73176327/2593278) + sigma = np.zeros((n_terms, n_terms)) + sub_sigma = sigma[nonzero_terms, :] + sub_sigma[:, nonzero_terms] = model.sigma_ + sigma[nonzero_terms, :] = sub_sigma + else: + sigma = model.sigma_ + # Sample the model parameters - params = np.random.multivariate_normal(model.coef_, model.sigma_, size=num_samples) + params = np.random.multivariate_normal(model.coef_, sigma, size=num_samples) # Assemble them into Hessians output = [] @@ -83,15 +112,21 @@ def sample_hessians(self, model: ARDRegression, num_samples: int) -> list[np.nda return output def _params_to_hessian(self, param: np.ndarray) -> np.ndarray: - """Convert the parameters for the linear model into a Hessian""" + """Convert the parameters for the linear model into a Hessian + + Args: + param: Coefficients of the linear model + Returns: + The harmonic terms expressed as a Hessian matrix + """ # Get the parameters - n_terms = len(self.reference) * 3 - triu_inds = np.triu_indices(n_terms) - off_diag_triu_inds = np.triu_indices(n_terms, k=1) + n_coords = len(self.reference) * 3 + triu_inds = np.triu_indices(n_coords) + off_diag_triu_inds = np.triu_indices(n_coords, k=1) # Assemble the hessian - hessian = np.zeros((n_terms, n_terms)) - hessian[triu_inds] = param + hessian = np.zeros((n_coords, n_coords)) + hessian[triu_inds] = param[n_coords:] # The first n_coords terms are the linear part hessian[off_diag_triu_inds] /= 2 hessian.T[triu_inds] = hessian[triu_inds] return hessian diff --git a/notebooks/proof-of-concept/1_compute-random-offsets.ipynb b/notebooks/proof-of-concept/1_compute-random-offsets.ipynb index 77ead94..8bde676 100644 --- a/notebooks/proof-of-concept/1_compute-random-offsets.ipynb +++ b/notebooks/proof-of-concept/1_compute-random-offsets.ipynb @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "c6a28419-6831-4197-8973-00c5591e19cb", "metadata": { "tags": [] @@ -38,7 +38,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "c6be56c5-a460-4acd-9b89-8c3d9c812f5f", "metadata": { "tags": [ @@ -64,7 +64,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "0b6794cd-477f-45a1-b96f-2332804ddb20", "metadata": {}, "outputs": [], @@ -83,23 +83,12 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "ad9fd725-b1ba-4fec-ae41-959be0e540b3", "metadata": { "tags": [] }, - "outputs": [ - { - "data": { - "text/plain": [ - "Atoms(symbols='O2N4C8H10', pbc=False, forces=..., calculator=SinglePointCalculator(...))" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "atoms = read(Path('data') / 'exact' / f'{run_name}.xyz')\n", "atoms" @@ -111,7 +100,7 @@ "metadata": {}, "source": [ "## Compute many random energies\n", - "Compute $3N(3N+1)/2 + 1$ energies with displacements sampled [on the unit sphere](https://mathoverflow.net/questions/24688/efficiently-sampling-points-uniformly-from-the-surface-of-an-n-sphere). This is enough to fit the Hessian exactly plus a little more" + "Compute $3N + 3N(3N+1)/2 + 1$ energies with displacements sampled [on the unit sphere](https://mathoverflow.net/questions/24688/efficiently-sampling-points-uniformly-from-the-surface-of-an-n-sphere). This is enough to fit the Jacobian and Hessian exactly plus a little more" ] }, { @@ -124,7 +113,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "23502eea-0974-4248-8f19-e85447069c61", "metadata": { "tags": [] @@ -137,7 +126,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "bf1366fc-d9a7-4a98-b9c9-cb3a0209b406", "metadata": { "tags": [] @@ -157,7 +146,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "d4f21e81-5ec3-4877-a4d1-402077be2ee8", "metadata": { "tags": [] @@ -179,22 +168,12 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "0915595d-133a-43df-84fc-4ff6a3b538ea", "metadata": { "tags": [] }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - " Memory set to 3.815 GiB by Python driver.\n", - " Threads set to 12 by Python driver.\n" - ] - } - ], + "outputs": [], "source": [ "calc = Psi4(method=method, basis=basis, num_threads=threads, memory='4096MB')" ] @@ -209,42 +188,26 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "e2a28593-2634-4bb7-ae5b-8f557937bda1", "metadata": { "tags": [] }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Need to run 2629 calculations for full accuracy.\n" - ] - } - ], + "outputs": [], "source": [ "n_atoms = len(atoms)\n", - "to_compute = 3 * n_atoms * (3 * n_atoms + 1) // 2 + 1\n", + "to_compute = 3 * n_atoms + 3 * n_atoms * (3 * n_atoms + 1) // 2 + 1\n", "print(f'Need to run {to_compute} calculations for full accuracy.')" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "8bf40523-dcaa-4046-a9c6-74e35178e87f", "metadata": { "tags": [] }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Already done 1. 2628 left to do.\n" - ] - } - ], + "outputs": [], "source": [ "with connect(db_path) as db:\n", " done = len(db)\n", @@ -256,41 +219,19 @@ "execution_count": null, "id": "a6fa1b33-defc-4b35-895d-052eb64453fb", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - " 0%| | 0/2629 [00:00 1e-7).sum()} are nonzero')" + ] + }, + { + "cell_type": "markdown", + "id": "aa509659-701d-4001-8cc7-980c9d999976", + "metadata": {}, + "source": [ + "Compare the forces estimated at a zero displacement to the true value" + ] + }, + { + "cell_type": "code", + "execution_count": 120, + "id": "70d80f87-9983-4bd5-a6ae-b9c966b0d838", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "actual_forces = data[0].get_forces()" + ] + }, + { + "cell_type": "code", + "execution_count": 121, + "id": "f548b145-0aa8-47f7-802b-6b7232a74bd3", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pred_forces = -hess_model.coef_[:actual_forces.size].reshape((-1, 3))" + ] + }, + { + "cell_type": "code", + "execution_count": 122, + "id": "d7cd7762-6e12-4dcd-b564-67a33b18d9e0", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Maximum force: 5.57e+00 eV/Angstrom\n" + ] + } + ], + "source": [ + "print(f'Maximum force: {np.abs(pred_forces).max():.2e} eV/Angstrom')" + ] + }, + { + "cell_type": "code", + "execution_count": 123, + "id": "425b77a9-7fd7-40da-af6f-eaed197c9ab6", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAC+CAYAAADa6ROSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAASK0lEQVR4nO3df4wcdf3H8dfu3t3c3o/d0jvgaO+X0h98+X5LQZvUNjHQIHI5KNqKKaSGUolJqRFjoiSXGNqgeLFaetFENKE/hBC1KhDFbxtpK4hV/NHaFL+QUvxGeqWUa2mzu3fX7t3ufr5/4H2+rr2Gud2dzt7M85FMLjs7O/v+3H1mXzszn5uJGGOMAACQFPW7AABA9SAUAAAWoQAAsAgFAIBFKAAALEIBAGARCgAAi1AAAFiEAgDAIhQCLBKJ6Nlnn/W7DITYjh07NGPGDL/LmJLpWHMlEQoV8vvf/16xWEw9PT1Tel13d7cGBga8KQqYonvvvVeRSOSCyU2/nqwvr1q1Sq+//rpH1f6/sH+QV1KN3wUExbZt2/SFL3xBjz/+uI4dO6bOzk6/SwJK0tPTo+3btxfNcxynpHXF43HF4/FKlIVLhD2FChgZGdHOnTt1//336/bbb9eOHTuKnv/FL36hRYsWqb6+Xq2trVq5cqUk6aabbtKbb76pL33pS/YbmSRt3LhR119/fdE6BgYG1N3dbR//+c9/1i233KLW1lYlk0ndeOONOnjwoJfNREg4jqO2trai6bLLLpP0Xt/s7OyU4ziaNWuWHnjgAUkX78v//g1+om9v27ZNnZ2dampq0v333698Pq9Nmzapra1NV1xxhR555JGimh599FEtWLBAjY2N6ujo0Pr16zU8PCxJeuGFF7R27VqlUin73hs3bpQkjY2N6cEHH9Ts2bPV2NioxYsX64UXXiha944dO9TZ2amGhgatWLFC7777rge/1emDUKiAn/zkJ5o/f77mz5+vz3zmM9q+fbsmLj77q1/9SitXrtRtt92mv/71r9q7d68WLVokSXr66afV3t6uhx9+WG+//bbefvtt1++ZyWS0Zs0avfTSS3r55Zc1d+5c9fb2KpPJeNJG4Gc/+5m2bNmiH/zgBzp69KieffZZLViwQNLU+vLf//537dq1S7t379aPfvQjbdu2TbfddpuOHz+uF198Ud/85jf11a9+VS+//LJ9TTQa1Xe+8x397W9/0w9/+EPt27dPDz74oCRp6dKlGhgYUCKRsO/95S9/WZK0du1a7d+/Xz/+8Y91+PBhffrTn1ZPT4+OHj0qSfrjH/+oz372s1q/fr0OHTqkZcuW6etf/7pXv8LpwaBsS5cuNQMDA8YYY8bHx01ra6t5/vnnjTHGLFmyxKxevfqir+3q6jJbtmwpmrdhwwazcOHConlbtmwxXV1dF11PLpczzc3N5pe//KWdJ8k888wzU2oLwm3NmjUmFouZxsbGounhhx82mzdvNvPmzTNjY2OTvnayvrx9+3aTTCbt4w0bNpiGhgaTTqftvFtvvdV0d3ebfD5v582fP9/09/dftM6dO3ealpaWi76PMca88cYbJhKJmLfeeqto/s0332z6+vqMMcbcfffdpqenp+j5VatWXbCuMOGcQpmOHDmiP/3pT3r66aclSTU1NVq1apW2bdumj33sYzp06JA+97nPVfx9h4aG9NBDD2nfvn165513lM/nNTo6qmPHjlX8vRAuy5Yt02OPPVY0b+bMmRoZGdHAwIA++MEPqqenR729vVq+fLlqaqb2MdLd3a3m5mb7+Morr1QsFlM0Gi2aNzQ0ZB//5je/0Te+8Q29+uqrSqfTyuVyOn/+vEZGRtTY2Djp+xw8eFDGGM2bN69ofjabVUtLiyTptdde04oVK4qeX7JkiXbv3j2lNgUJoVCmrVu3KpfLafbs2XaeMUa1tbU6e/ZsSSfZotGoPfw0YXx8vOjxvffeq1OnTmlgYEBdXV1yHEdLlizR2NhYaQ0B/qmxsVFz5sy5YP7MmTN15MgRPf/889qzZ4/Wr1+vb33rW3rxxRdVW1vrev3/vmwkEpl0XqFQkCS9+eab6u3t1bp16/S1r31NM2fO1O9+9zvdd999F2wX/6pQKCgWi+nAgQOKxWJFzzU1NUnSBdsZCIWy5HI5PfHEE9q8ebM+/vGPFz33qU99Sk899ZSuu+467d27V2vXrp10HXV1dcrn80XzLr/8cp08eVLGGHvC7tChQ0XLvPTSS/re976n3t5eSdLg4KBOnz5doZYBk4vH47rjjjt0xx136POf/7yuueYavfLKK/rQhz40aV+uhL/85S/K5XLavHmz3ZvYuXNn0TKTvfcNN9ygfD6voaEhffSjH5103ddee23RuQtJFzwOG0KhDM8995zOnj2r++67T8lksui5O++8U1u3btWWLVt088036+qrr9Zdd92lXC6nXbt22ZNk3d3d+u1vf6u77rpLjuOotbVVN910k06dOqVNmzbpzjvv1O7du7Vr1y4lEgm7/jlz5ujJJ5/UokWLlE6n9ZWvfIWhf6iIbDarkydPFs2rqanRc889p3w+r8WLF6uhoUFPPvmk4vG4urq6JE3elyvh6quvVi6X03e/+10tX75c+/fv1/e///2iZbq7uzU8PKy9e/dq4cKFamho0Lx587R69Wrdc8892rx5s2644QadPn1a+/bt04IFC9Tb26sHHnhAS5cu1aZNm/TJT35Sv/71r0N96EgSJ5rLcfvtt5ve3t5Jnztw4ICRZA4cOGB+/vOfm+uvv97U1dWZ1tZWs3LlSrvcH/7wB3PdddcZx3HMv/45HnvsMdPR0WEaGxvNPffcYx555JGiE80HDx40ixYtMo7jmLlz55qf/vSnF5zoEyeaMUVr1qwxki6Y5s+fb5555hmzePFik0gkTGNjo/nIRz5i9uzZY187WV+e7ETzvw+iWLNmjfnEJz5RNO/GG280X/ziF+3jRx991Fx11VUmHo+bW2+91TzxxBNGkjl79qxdZt26daalpcVIMhs2bDDGGDM2NmYeeugh093dbWpra01bW5tZsWKFOXz4sH3d1q1bTXt7u4nH42b58uXm29/+dqhPNEeM4aAaAOA9/J8CAMAiFAAA1rQKhWw2q40bNyqbzfpdimfC0EYpPO2spLD8zsLQzmpu47Q6p5BOp5VMJpVKpYpG4gRJGNoohaedlRSW31kY2lnNbZxWewoAAG8RCgAAq+R/XisUCjpx4oSam5vtf916LZ1OF/0MojC0Ubr07TTGKJPJaNasWUXX2CkV/d87YWinH210uw2UfE7h+PHj6ujoKLlAwA+Dg4Nqb28vez30f0xX77cNlLynMHGVwyOvHy264mEQnTlf+eu5VJvWeOz9F5rGMpmM5s6dW7G+OrGe1ZqtuoAfhV33+n6/S/Bce3Od3yV4LpPJ6Jp5778NuA6FbDZbNHxq4mYuzc3NVXf2vNLGa4MfComGYIfChFIP9Vys/9cpGvhQaGoO9vYtSYlE8ENhwvttA657c39/v5LJpJ3YdUaY0P8RFq5Doa+vT6lUyk6Dg4Ne1gVUFfo/wsL14SPHceQ4jpe1AFWL/o+wKPt+Cs7ZY3JyTZWopWpd2VSZ68JXM6Ng/w0vzaDRYGqqDf75plQ2+OcNMy7bGOwzZACAKSEUAAAWoQAAsAgFAIBFKAAArLJHH51patd4wP/jsT4W/LErDLbExez53zN+l+C5u//rcr9L8FzNmLtRZOwpAAAsQgEAYBEKAACLUAAAWIQCAMAqe/RRc11UibpgZ0t0bNTvEjxnYg1+l4AqtfqaYN9ES5JqThz2uwTP1WVGXC0X7E9zAMCUEAoAAItQAABYhAIAwCIUAABW2aOPouczitZVopTqNe4E+9pOUgU6AgJrNBL8K2OZKxb4XYLn0vVpV8uxpwAAsAgFAIBFKAAALEIBAGARCgAAq+xBJxFjFDGmErVUrdrxEFz7qI5rH2Fy8Wiwt29JGjPB/35cE3V3B8ng/yYAAK4RCgAAy/Xho2w2q2w2ax+n0+7+EQIIAvo/wsL1nkJ/f7+SyaSdOjo6vKwLqCr0f4SF61Do6+tTKpWy0+DgoJd1AVWF/o+wcH34yHEcOc6F10Ap1DepUB/wOzNFgn/qJejjS8pt38X6fxhEssN+l+C5jBr9LsFzw2N5V8sF/9MOAOAaoQAAsAgFAIBFKAAALEIBAGARCgAAi7swumEKflfgvRAMu0VpwnA72haT87sEz9XWuRuYzScBAMAiFAAAFqEAALAIBQCARSgAAKyyRx+dPFfQSE2wR+e050/5XYLncomr/C4BVar+zD/8LsFz7zZ3+l2C5zI5d/sA7CkAACxCAQBgEQoAAItQAABYhAIAwCp79FFzbVTNdcHOlqH8FX6X4LkZAb8fZyHg7fNSrqXb7xI8d+hY2u8SPDcyPOJquWB/mgMApoRQAABYhAIAwCIUAAAWoQAAsMoefdRUOKfmQrBv4NasYF/bSZIKkWDfXSsa8buC6av25Kt+l+C5ZZe1+F2C59I151wtx54CAMBy/RU/m80qm83ax+l08Mf1AhPo/wgL13sK/f39SiaTduro6PCyLqCq0P8RFq5Doa+vT6lUyk6Dg4Ne1gVUFfo/wsL14SPHceQ4jpe1AFWL/o+wKHvYUCYal6INlailauWDP/hIM/wuwGMMPipdqvU//C7Bc03jIThHFBt3tRijjwAAFqEAALAIBQCARSgAACxCAQBglT36aCRXUHQ82MNz2ur9ruBS4PsBJnfmXM7vEjyXOPeO3yV4Ljo87G45j+sAAEwjhAIAwCIUAAAWoQAAsAgFAIBV9uij2aPHlYg1VaKWqjUU7fK7BM/NDMUIK5SiLhb87475mZ1+l+C5fI276zsF/68NAHCNUAAAWIQCAMAiFAAAFqEAALDKHn1k6upl6uKVqKVqtUTO+V2C54yCPYIMpWt77b/9LsFzZuEtfpdQNdhTAABYhAIAwCIUAAAWoQAAsAgFAIBV9uijd2ou12htohK1VK1WJ+Z3CZ6Ljbzrdwmeio5m/C5h2sp/eLnfJXguNnrW7xI8Fxk/72o59hQAAJbrPYVsNqtsNmsfp9PurrgHBAH9H2Hhek+hv79fyWTSTh0dHV7WBVQV+j/CwnUo9PX1KZVK2WlwcNDLuoCqQv9HWLg+fOQ4jhzH8bIWoGrR/xEWZY8+uiL/rhK5sUrUUr1CMHCl0Hyl3yV4qpCv9buEaSvoI9Mk6VRsht8leC4Tcbcco48AABahAACwCAUAgEUoAAAsQgEAYJU9+mgo1qJzNQG/9lE8+Nc+SmULfpfgqUzA2+elSD7gowslKfibuGvsKQAALEIBAGARCgAAi1AAAFiEAgDAKnv0UUs8pkTAR+ecOZ/3uwTPtdQH+28Ydfj+U6pMQ5vfJXiu7dT/+F2C5xoyw66WY0sBAFiEAgDAIhQAABahAACwCAUAgEUoAACssoekjhfem4KsxXF5H7tpLBrwWy5GR0NwT1WPNJ0/7XcJnjMjab9L8JwZHXG1HHsKAACLUAAAWIQCAMAiFAAAFqEAALDKHn2EYDhTc5nfJXgqEwv2Bf88FavzuwLPvXLZh/0uwXPDNe5GWLGnAACwXO8pZLNZZbNZ+zidDv64XmAC/R9h4XpPob+/X8lk0k4dHR1e1gVUFfo/wsJ1KPT19SmVStlpcHDQy7qAqkL/R1i4PnzkOI4cx/GyFqBq0f8RFmWPPqqLvjcFWj7ndwWem+EEe4QJt+MsXWR81O8SPHdtU/BHp6UL7j7H2FIAABahAACwCAUAgEUoAAAsQgEAYJV/7aP82HtTkIXg2i81p97wuwRP1WSG/S5h2krFr/S7BM/NHHrF7xI8V5vhzmsAgCkiFAAAFqEAALAIBQCARSgAACzuvAZJ0j/qu/wuwVOZce5/UKr6muB/d9yd+4DfJXhuNJ9xtVzw/9oAANcIBQCARSgAACxCAQBglXyi2RgjScpk3J28mNZCcJmLzPC43yV4avif/XSi35ZrYj1jKlRkfdUsnQ7+SfrR4eB/jk208f22gZJDYSIM5lzzn6WuArjkMpmMkslkRdYjSU/prbLXVe22X9XmdwmooPfbBiKmxK9OhUJBJ06cUHNzsyKRSMkFTkU6nVZHR4cGBweVSCQuyXteamFoo3Tp22mMUSaT0axZsxSNln/UlP7vnTC00482ut0GSt5TiEajam9vL/XlZUkkEoHtLBPC0Ebp0razEnsIE+j/3gtDOy91G91sA5xoBgBYhAIAwJpWoeA4jjZs2CDHcfwuxTNhaKMUnnZWUlh+Z2FoZzW3seQTzQCA4JlWewoAAG8RCgAAi1AAAFiEAgDAIhQAABahAACwCAUAgEUoAACs/wNUbjalucuDGwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(1, 2, figsize=(4, 2))\n", + "\n", + "for ax, l, h in zip(axs, ['Actual', 'Estimated'], [actual_forces, pred_forces]):\n", + " ax.matshow(h, vmin=-0.05, vmax=0.05, aspect='auto', cmap='RdBu')\n", + "\n", + " ax.set_xticklabels([])\n", + " ax.set_yticklabels([])\n", + " \n", + " ax.set_title(l, fontsize=10)\n", + "\n", + "fig.tight_layout()" + ] + }, { "cell_type": "markdown", "id": "46a0f2f8-f863-4de3-bd97-97ebd92676d4", @@ -173,7 +280,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 124, "id": "00a10907-667a-413c-851d-d47f0eff092b", "metadata": {}, "outputs": [], @@ -191,7 +298,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 125, "id": "d48893fd-df0d-4fa8-bfbe-0d04b71fbf1a", "metadata": { "tags": [] @@ -205,7 +312,7 @@ " [1.08009177e-03, 3.94902961e-03, 4.15881408e+00]])" ] }, - "execution_count": 9, + "execution_count": 125, "metadata": {}, "output_type": "execute_result" } @@ -216,7 +323,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 126, "id": "9b311dea-5744-4211-81cb-40aa1183301e", "metadata": { "tags": [] @@ -225,12 +332,12 @@ { "data": { "text/plain": [ - "array([[ 22.92078111, 0. , -0. ],\n", - " [ 0. , 104.76017451, -0. ],\n", - " [ -0. , -0. , 17.33479829]])" + "array([[3.26945590e+01, 4.95882374e+00, 8.59920621e-03],\n", + " [4.95882374e+00, 3.93490213e+01, 6.10182117e-01],\n", + " [8.59920621e-03, 6.10182117e-01, 2.73150623e+01]])" ] }, - "execution_count": 10, + "execution_count": 126, "metadata": {}, "output_type": "execute_result" } @@ -241,7 +348,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 127, "id": "addd7bef-854a-4b9f-96e9-2aa01b652495", "metadata": { "tags": [] @@ -249,7 +356,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAADJCAYAAAA3tRlxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzvElEQVR4nO2deZQc1X3vv9090z2LRjOa0S6EZFlCaBACCSFAsULYMWZL5AjbCiCIl0AgYMd+MScQtpxDXmzZ2M8cB4yQTeA940RADsQEgUGgREiyEEIGBo2E9mW0jmbRrN193x936arbVdW3uqtnpmt+n3PmzHTVrVu3aqpvfe/v/u7vF2GMMRAEQRChITrYDSAIgiCChTp2giCIkEEdO0EQRMigjp0gCCJkUMdOEAQRMqhjJwiCCBnUsRMEQYQM6tgJgiBCBnXsBEEQIYM6doIghgTLli3DjTfeONjNCAXUsftg2bJliEQiWT9XX331gJz/oYcewrnnnjsg5yJKl3Xr1iEWiw3YcxkUP/nJT/DLX/6y6OcZDi+QssFuQKlx9dVXY+XKlbZtiURikFpDENk888wzuPvuu/H0009j7969OP3004t6vv7+fpSXlxdcT21tbQCtIQAAjDDm1ltvZTfccIPjvrfffpuVl5ezd999V2374Q9/yBoaGtjBgwcZY4y99tpr7I/+6I9YbW0tq6+vZ1/60pfYjh07bPXs27eP3XTTTWzUqFGsqqqKnXfeeWz9+vVs5cqVDIDtZ+XKlcW6VKJE6ezsZDU1NezTTz9lN910E3v44YfVvrfffpsBYK+++iqbM2cOSyQSbMGCBWzr1q2qzMqVK1ltbS176aWX2IwZM1gikWCXX34527t3ryrz4IMPsnPOOYetWLGCfe5zn2ORSISl02m2Z88edv3117Pq6mpWU1PD/vzP/5y1tLQwxhhrampilZWV7Pnnn1f1rFq1iiUSCXV+/ft18cUXs7vuuovdc889rK6ujo0dO5Y9+eSTrLOzky1btoyNGDGCTZs2jf32t79VxySTSXb77bezqVOnsoqKCnbGGWewxx9/3NZ2/Xv09ttvM8YY279/P1uyZAmrq6tj9fX17Prrr2e7du0K5P8y0FDH7gOvjp0xxr73ve+xKVOmsJMnT7ItW7awRCLBXnzxRbX/3//939mqVatYc3Mz++CDD9h1113Hzj77bJZKpRhjjHV0dLBp06axRYsWsbVr17Lt27ezF154ga1bt451dXWxv/3bv2VnnXUWO3ToEDt06BDr6uoq9iUTJcaKFSvY/PnzGWOMvfLKK2zq1KksnU4zxjId+6xZs9jq1avZ1q1b2bXXXsumTp3K+vr6GGO8Yy8vL2fz589n69atY5s2bWILFixgCxcuVOd48MEHWXV1NbvqqqvY5s2b2YcffsjS6TSbO3cu+8IXvsA2bdrE1q9fz+bNm8cuvvhiddwTTzzBamtr2e7du9mBAwdYfX09+/GPf6z2O3XsNTU17NFHH2XNzc3s0UcfZdFolH3xi19kTz31FGtubmZ33HEHa2hoYKdOnWKMMdbX18f+4R/+gW3cuJHt3LmTPffcc6yqqoq98MILjDH+HVuyZAm7+uqr1feot7eXnTp1is2YMYPdfvvtbOvWreyTTz5hX/va19jMmTNZb29vMf5VRYU6dh/ceuutLBaLserqatvPI488whhjrLe3l82dO5ctWbKEnXXWWezrX/+6Z31HjhxhANgf/vAHxhhjTz75JKupqWHHjx93LC+VEkG4sXDhQqVQ+/v72ejRo9kbb7zBGMt07L/+9a9V+ePHj7PKykrV8cmR4fr161WZpqYmBoBt2LCBMcafw/LycnbkyBFVZvXq1SwWi9mU/ccff8wAsI0bN6ptX/rSl9iiRYvYZZddxq644gr10mHMuWP/whe+oD4nk0lWXV3Nbr75ZrXt0KFDDAB77733XO/JnXfeyRYvXux6Hsb4C3HmzJm29vT29rLKykr2+uuvu9Y9VCEbu08uueQS/PznP7dtq6+vBwDE43E899xzmDNnDqZMmYLHH3/cVu6zzz7DAw88gPXr1+PYsWNIp9MAgL1792L27NnYsmUL5s6dq+ojCD9s27YNGzduxIsvvggAKCsrw0033YRnnnkGl19+uSp30UUXqb/r6+sxc+ZMNDU1qW1lZWWYP3+++nzmmWeirq4OTU1NWLBgAQBgypQpGDNmjCrT1NSEyZMnY/LkyWpbY2OjOu78888HwO3/Z5xxBqLRKD766CNEIhHPa5ozZ476OxaLoaGhAWeffbbaNm7cOADAkSNH1LZ/+Zd/wdNPP409e/agu7sbfX19OZ0O3n//fezYsQM1NTW27T09Pfjss888jx2KUMfuk+rqakyfPt11/7p16wAAJ06cwIkTJ1BdXa32XXfddZg8eTJ+8YtfYOLEiUin05g9ezb6+voAAJWVlcVtPBFqVqxYgWQyiUmTJqltjDGUl5ejtbXV81i9g3XqcK3brM+1PI/TMfr2Dz/8EKdOnUI0GkVLSwsmTpzo2S59UjYSidi2ybqlSPrNb36Db3/721i+fDkuuugi1NTU4Ac/+AE2bNjgeZ50Oo3zzjsPzz//fNY+6wusVCB3xwD57LPP8O1vfxu/+MUvcOGFF+KWW25RD9zx48fR1NSE+++/H5dddhlmzZqV9WWbM2cOtmzZghMnTjjWH4/HkUqlin4dROmRTCbx7LPPYvny5diyZYv6+fDDDzFlyhRbh7V+/Xr1d2trK5qbm3HmmWfa6tq0aZP6vG3bNpw8edJWRqexsRF79+7Fvn371LZPPvkEbW1tmDVrFgAudpYtW4a///u/x2233YalS5eiu7s7kOuXrF27FgsXLsSdd96JuXPnYvr06VmK2+l7NG/ePGzfvh1jx47F9OnTbT+l6K1DHbtPent70dLSYvs5duwYUqkUbr75Zlx55ZW47bbbsHLlSnz00UdYvnw5AGDUqFFoaGjAU089hR07duCtt97Cd77zHVvdX/3qVzF+/HjceOON+J//+R/s3LkTq1atwnvvvQcAmDp1Knbt2oUtW7bg2LFj6O3tHfDrJ4Ymr776KlpbW/GXf/mXmD17tu3ny1/+MlasWKHKPvLII/jd736Hjz76CMuWLcPo0aNtft3l5eW4++67sWHDBmzevBm33XYbLrzwQmWGceLyyy/HnDlzsHTpUmzevBkbN27ELbfcgosvvliZdf7qr/4KkydPxv33348f/ehHYIzhu9/9bqD3Yfr06di0aRNef/11NDc344EHHsDvf/97W5mpU6di69at2LZtG44dO4b+/n4sXboUo0ePxg033IC1a9di165deOedd3DPPfdg//79gbZxQBhcE39pceutt2a5SgFgM2fOZA8//DCbMGECO3bsmCr/8ssvs3g8zj744APGGGNvvPEGmzVrFkskEmzOnDlszZo1DAB76aWX1DG7d+9mixcvZiNHjmRVVVVs/vz5atKqp6eHLV68mNXV1ZG7I2Hj2muvZddcc43jvvfff58BYMuXL2cA2CuvvMLOOussFo/H2fnnn8+2bNmiykp3x1WrVrFp06axeDzOLr30UrZ7925Vxm0S38vd8Ve/+hWrrq5mzc3NqvymTZtYPB5n//mf/8kYc548veeee2znmDJlis2ThjFm+w719PSwZcuWsdraWlZXV8fuuOMO9v3vf9/W3iNHjrArrriCjRgxwubueOjQIXbLLbew0aNHs0QiwaZNm8a+8Y1vsLa2Nsf7OpSJMEbJrAliOLBmzRpccsklaG1tRV1dnWOZX/7yl7j33ntx8uTJAW0bESxkiiEIgggZ1LETBEGEjII69t7eXjz00EOBT+KVUr2l1FaqN1hK7ZovuugiPPjgg55utcuWLfNthim1+zAc+oGCbOzt7e2ora1FW1sbRo4cmW81JV1vKbWV6g2WUrtmqrd49Q61tpIphiAIImRQx04QBBEy8g4pkE6nceDAAQB8uBAksr5SqLeU2lrq9TLG0NHRgYkTJyIaHVxNQs8/1VvsOp3qNf0OGNvYe3t7bQb8AwcOoLGxsZA2E0Re7Nu3D6eddtqAnpOef2Iokes7YKzYH3vsMTz88MNZ2/9j3VZ868drAABb/s8StT3C0uIP7a3itt1CSrxqoiJ2UKw78xZk5XEAQLqsgv8WZWOirHxLWcMRRTt55Lf0iLEAgKT2KotZCsePNAMA3lvGl/vP/+2rtja5HcdPzq+NeVybJOtQpzIF1JcP+cyiW8/rdnzEYX+u63eqt7OjA9NnzMiKwDcQuD3/27dvx6w/ewgAsP+N5QPbqCJS9um7AIDkmX88yC0hrHR0dGCGwXfAuGO/7777bLFN2tvbMXnyZFSPqEG0nLtPWWdtA+3YyzNfe1bO09D56tgjPNBQuoa3z7Nj7x4BAKiOldmuiTr23Octdscut+UK9VoM3J7/mpoaRGJcbAw1r51CKBtRBQBIhuiawkSu70De7o7SDaelpUU90D2W3q9C6/W8vsh+vqbyuJTWoZvU0SfeKeUe/aReT7TzKD/fiDFZ+/UOx+2zU/0mHZlbm0zundc/1at9hbQl1/n84FRvR3s7xo0fPyRcH52e/11tfWr/52rjg9U0IsS0t7djvMF3gLxiCIIgQgZ17ARBECGj4AxKEWRswRWxzHvitMU/AgB88K/3AADGVsUAANGOwwCAU1VjVVlptnEzaUS727JPXFnreYyTuSIumqfby61Wo6Rm4kkLE8y6C/gk0qI3/02VTdWMs9Xj1ha+0XtuwdZeUTYl3ru6Ld/JtOHHnmZignHb5nhtHttyHeO1T693qIchtZpfZt31HwCApp/dMFjNyUlnP3/ORnjZJgXRDasAAOkLFhe1TUQwkGInCIIIGYXnPGVpRxUqlXp3UihVcMXeXsEVcEU0o9FyKrJ0MvNnhd3Nx0R9ZnnZIG0ryyzvN+llUwa7R8qi33HFEu09lalXO6enotTukeeEoigbM6nXrQqHY3NNfJpMuOplC/WKCWJEMBQZykpdYqLUJUNBqZfKyG0oQIqdIAgiZBSu2C1K1KqopE1d6s5uIZurhUqwlpVuiHHd5V38TlU3ZJ22X3Nd9LJv6zZqL79w2QamvfOku6P8DQAxMV+QdrG1e2FiW3YrU6irqNsxXjZwk7bkqs/EVdSpXlJqBED/fz+QYicIgggZBSt2Bjh6fEjvF2lTl0r9zDtfBAB8/MSfqbJKqaf4Ao9UlHsXSKXda3FjiQlDuTwmkuwBACRjfCWqtKNbPV/K+joBAOkEX1UaWfMs/73wy/wa4lWqbKSvi58zxlfTxkUj4ge3AgA6x52dabdQ6sperK0ULT++W5VN1Y7nbRArZiUmNuWslZcOx0v0lbhqBTCASD+/V2nL9VqPkffOWr8ftezHKybXQi3r9jAptb/4v/w5eu5rcwa5JeHGj8dPIVi/XyYrxAeKodMSgiAIIhAC8WOXSt2qrKSfuvR+kUpMKvV5f/e6Krv1n6/ix4uYGzHYKbN40OjquzfCY8dENNUpy1nLSqILruPnOb4TANA3Ybbalyznajauqde2sVyp1zS9qcqmZl/B2w17WVVX3cRMu3dt4m2ZvtBW1lGNuvi8e9qzxTFlYtTD5MjAUgfTlLp+jD6asJV1aYOJV4wJJmEYSskrxg2p1G98ZrPa9vLt8walLX5iEUnKTuxRfyfrpwTepqAotlLvEyaBeGxoauOh2SqCIAgib6hjJwiCCBkFR3c8bIluZ0U3T7i5NALAkS6+1CfjIulch3WbPknodYxEX6gkifZmm21M3PpiXa38mKpRrmWyjhHhEVIuIRG88DRXSPONWMwlJ6D1a7XWI++HHinTek9zuUj6eXjyjYwpt7UP8eiO+dLWy/93tQnSWYQ3FN2RIAhimBKMu6MDuvLSXRrlRCmQUepnf/e/AABr/pFPStZX8O1OSS6UEhX1ISouRUwERUTQIgBgC/4UABCT+4SLJJKiLRWZN1/UxX1SjRCYJbyBUOplR3cAAE7WTQPgPHGjFlsJpS7dJ1Oj+ARUWmy3lpXo97LPckMSQl7LCbCIywS0E1FNoataU5m44lE5wSYmVv0E7dLxWnTkZ1FX2JBK/YKH3gIAbHjo0oLrlM8x4D0pToQTUuwEQRAhIxB3R91GC7irK7X4yGGfVOrn3vpzAMDeF+4CYLcT60ovKuzcUrHL8AP9CzJBixIn9/F9oybzY0XZSCSZ1dZu4T5ZoUlIpWoj2bfseK1Q6q/9lJe5/l5bW53aLZV67OhnAID06e4ub7qqTUSZZV/Esax+Xi/0+uX/yIpUABGLmgcyIwR+oH2fHJXlE/bAiTC4O3ohlboMeb1/1Xe8insyGCpdhtjQw1kXSkTMG8nvbbEoa+X9RFL0E6WM8Z3Ss7S3t7d7lCaIcEHPP1FKGHvFPPTQQ45Z2q1eAfkEtLLu09nX0Q8AmFxTrrbl8iLpd/C+Ueq+h38ZrTZ1wG7Dz1pWry0WMgnAFT11nNfrELzM7VhrMhGrvV0/p36sn6QZpnhdo54I3BZszTL/AAApMbpx8szRjzdR6m2GHgHFwOT5D5p24Uo20smVjBiWBO4Vc99996GtrU397Nu3L5CGEkQpQM8/UUoYm2ISiQQSiUTW9mhfN2LdXG9ZlaZSoNKvWqhXGdDLGiZA2q91P3Op1Ofdt1qV3fzYlQDc/cxVGF+Lb3oyLkIKCKUubXYqoYVlSXXsE+6ZkGy81Fam/EgzAKBnzBlZ7dZVZ1pcq58wuCnLvYu3fAIAaB8zizc75m5Hz6V4vZbom6h9vZ36NdvOJRS63CfvqvIosiY4kfdc/C/6xaOoxKkluUpfANGlC8Xt+S8mUqlP//rzatuOp5cW9ZwRB6+1gSbWeRSAPUQ24Q8a4xEEQYQM6tgJgiBCRuHujslusHJuRnAa2us5SmU8da/oi/owX5pfAPdFHPoxveWWOnV3TI8J0a6ZlwAAqoXpJTmWm16SY6YDABLb1qiyyTP/BE7k47Jn3d42mptgao82AQD6xzfa6vUylXjVa5oNyWQ23c9kuHKHtU6uiqG+dGHTB/5W1za5r9ugXWHEan659dd/AAD86itnuxUviFwmmJhwDAC8nQMKgUwwhUOKnSAIImQUHATM6u5ldRvUg3PpbojWk+oZfLIaaflbHne4i6u/8VVltu1eE4NugcO8kJPAugsikFkclRKhBUwCh0m8ltTr23Kdx+lc+bqe6nUVsrjI637IfdJ9siytZc9i2ZOnpzraMXHC0AoCdvjgftTU1Q9qW0qZXN99wg4FASMIghimBOJH5qSEdUWq3BBFcCKZ+QjIBLKSS9KVQvVYmiyV+pVPbAAAvPytBQCAyrLsV39rDw8LLIOKle1cDwA4MvF823Yg2x4slbpSuZYchzIIGN5cweudexkAINkw1bXdeshcE6WilHo6mbVPX2btZ8GPm3LPV6XnE9pX/rukbVcpDUvoBrk8zc9Ia6Ao9jL3XPxs4wEAwF0LJg1qO/KFlHpxIMVOEAQRMgKRG05Kys3eKsPhRpxknFQ/mgryWpAjlfriFTyn6GvfOj+r2gacAgCkwW1Sh4VSH7PtDQBA6pyrM9cCZ48ZdYkO+SHLz+WeNGg9KE441X6My7W4lXFT1CmHAGQxrWw+EyZ+PFz8qHuTeY9cHjq5jh/uSKV+zZO/BwD81uH5J4YfpNgJgiBCRiCK3Ul1uSkwaVOz2taUEhNqWIXeFeLZKeG4PFza1KVSl2n2xljS7MmgX7K+BmFTTwulLpcwA0Ba86HV5w+cVGNyNA/bC/E76rAkOqL91utLWyp2C1WQFaDMgpuq9aPCvcIPmOCnXrf6vVQ94Q4p9eIyUKGDg4IUO0EQRMgo+PXTn/ZWhbpSld4g1pWnTCaQFunsZJKMuENAL7mi1EnFAxmlfvZ3XlPbPv6ni/kx8Sr+W6SlS1dyb5OUJbB+/MCHAIAjo/nKPj3BsJOvfkSEA26N8LbVCaXueV+EAogIBeDl8RG0jdlNJZvYwnPV5XWMl6++n7JhRaazCyJJxjn/63X194f/fFXB9Q13SkWpS0ixEwRBhAzq2AmCIEJGweOLaASIilyHacuCIt30Ik0N0gRjDfylhuML/hRAJkepzD2o4qkDWWN2ufhIujTKiVJpfgGACTc9AQA48NL3eJtGTuD1vvsbfv5r/lqVPTqGm2DGtGwBAPSLXKRqoVUqk/2diSFzb5yfs75bZE5KGARHEgudnBYq+ZkI1TFxG3Q7plgmDz8mnnzMQWEhyDylVvPL+Q+8CQD4/aOX28oMZtzzmAjVkXII1UEUDil2giCIkBFIELBaEYymz5IgRw/2lU9wLjkpma7IDnbj5obo5SJ5VLhCjrW4QgIZ90Qgo15UAC6PIGBuN04fpTgRVHCtQiY1/QQBM1HW+Uzyuo0WnOptb2/HuEHKearjFASPIIqNaRAwY1MMZWknhjP0/BOlhHHH/thjjzlmaY+wtAq9alXJSrWKJfpMLD6KrHkWABBdcF2msAxyJdy9pGsR03OUAso2HRNBo2RALxkmQC4+ki6NQMamPlao8bHXPAIA2P7y/QCAkRYb4ykh+UeIi5FKXY4eOmIZe3+1KFPWyucEpNukcoO0BAxTyT2k2xSz3xcrJkkzcuEnWYZev8l5vMINu+FVtth2/kJxe/5LEekKSW6Q4cXYxk5Z2onhDD3/RClhrNjdsrQ7KU4vIgu/DACIHd+ptsmwtEjysL2RiFi+Kz1nrOfQzndEC+glwwTIxUdAxvslKrxfpFKffxffvv2pr6iy8ZgYWQiFLkcN0vMlms7WpR01pwEA+sUkg8wuH+nPeNCkxeIo5S0k3ql6EC/r335Uq5Hi1UYJ+jFeIQUKScphEqrAS+0PBW8Yt+ffjU5t5KdjHc35/Q75Rfd+kUr9yp/x0e7quy4s6vkLJeIxuiWcoTtFEAQRMgr2ijnc0oIah9lZP94UUREyQCp0PTCW1bskV5hXqU6swbz0fbrnS3Nrnyo7Y1Tcs34TFe2VRCOIRBiFKmA/3it+bN/5BBfzU28HecUUBTm6ANxHGMTQgFLjEQRBDFOoYycIgggZBYcUKD/SjHg3N6H0jztTbZcukNKcoiI19nXx/eVVmUqECUZGt+sW+VArhA0m9slbqmjXTJ6tSOZJ1d0qpQlGRmkEMmECpFujdGmUE6XS/AIAV/z0PQDAS3dcAACoKbdPNMbaD6my0o1SIdwyo8Kl0Y8roBXp9imzTTmkcc3Cz0KlXGX9xGU3mTx12p/rflAGpQyRVMZUKHPDBonV/HLZ4+sAAL+7d2Hg5yEGDlLsBEEQISOQIMP/veROAMAF72SUtQr6JRcoiXdIb6wSABB3CHolFWqFliko2XipKlt9pJlvG3sGPw8Ti5d0N0gRTx3IDuglFYru0ghklPqki+8GALSt4wHEomKk0VuTUenSVVG1U+QkdbqpbmrWaYK4X9yHtMHMYq5J03wmUQdzcZCfbEtDhTIxorQ+p0FSDJXuhlTqtQu5a7B8/omBIajRKSl2giCIkBGIu6N0u4k65A4txPbrtDDBLTiX7iLpFAbXK6CXW1s2t3Clft74qqw26u0LKsdnrrKFZhfyY2PP1V4/4Qe8IHfHockJERa7viKWoyQxEJC7I0EQxDAlEBu79ExhlkVB6y74YwDAot/xPKZyUVBcBOdqG5uxgVcKt48se7NQwuXCrg4AyTHTeVlNqSuvFQfpp9S8FtBLhgmwrslQ6r7PrtRHCpvj4bU/y7Q7Zn8v6srUsu5DtSsrNHF2c1X4Y1nWaaFTEJiMHvyMCArxthlKdv4giIr8vekLFg9ySwpDKnWyuZcWpNgJgiBCRsGKvS8NVIu/rSpr0Zv/BgCI9vKUdVKxd47jSr2m6U1VNjn7CgAZDxcWsTerZ8wZ6u/EtjW8vjP/hJ9TBlMS6l4qv5RFAsp0djKVnQy9KwN6xS2SWPqpS+8XaVmUSv3KH/63Krv27xbZ2qm3xaq0y7Q2eNmd5QhCjhr0AGK2c4rf+dix/ahkk/AGuXzd8/XrL0VyKfWyE3vU38n6Kcb1xk6J9IvVBukXA0QqdbnOAwDe+JuLBrQNpYpJ4p2gIcVOEAQRMqhjJwiCCBmBujt6VRTrOAwASNWM4ye27It2tQIA0iIuu0me1Jh2DDzK5jJXeLle6uYVK+v286iUC0+zR6V0qtftmvKZlDRxdwzK/BGEu6PTsblynTq129TVayAYLu6OXuSKOU8ED7k7EgRBDFMKnjw1DQKV1pS69TipusuO7gAAHK+dBgCojdsX/lhRWZfeXAEAKD+XBwdLjubHSpdGINutUeYolZmPqhyibKl2emRtkUpdukJ2vvmPtrZ5uQ8q1840D/BkXTaut0Z9Fvlij16wVO0bU+m8cMRksZFb27z2FbIQKqgRDDE0kEqdXCGHHsYdO2VpJ4Yz9PwTpYRxx+6WpT0CZOXSVNvhHibApt7EcSfruNoe8dpP+Y7r73Wsw3p82dzL+B+tB/lvodhbIyNU2fpuu4tYatRkAJkcpbbWiNC7KqBXxH4dTrkqpVI/59HfAwA+eOzKrPZmLaRCOqtMrmOOLOBKfUJfiyqTrJykX4Hvev2UNVHhOl5zHIW4UQ4kbs9/MdDno4Y6UqlP//rzAIAdTy/1Kj5kcQpHUqoY29gpSzsxnKHnnyglAs156qXe/ATIUjlKxWKMtGUxRi61GfXIeeq2UCBtaZR8W+ejUOW+7SKHqjWBR06bsmUkoC+28lKz+iImPwTtSWNSn1u95BUTHr7/X9sBAP909YxBbkn4IK8YgiCIYUqgXjEmdtHy47sBAMm6iZnjtEQC8lhpE/ej+GToAqdjpFJXdv5+vswfFrXr1m4Z0MsWJkCzv0ukUl/0v9eqbf9971wAQFqkAcxSqJb5CTeV7GijFm2XdtlTVWMBAL0iN+HId55WZdNXfJOXlcHFtGOqe07wcpYRUq4wAV77TIKMFRISgRiaSKU+777Vattmh3knoniQYicIgggZ1LETBEGEjEDisTvishQ/VTuen3jXJrWtfzrPsxjEIhWjY0XbnCYc3dqgcriKKI0AkBaRGvVjJdL8AgDL/mM3AOCZJbON2+vHTNFZyc0pNae4K6R09kxd+U1VJqbleNWPaSnnn0d7tMV0O1Ca+UuJ4LCaX2RUSIoIOTCQYicIgggZgSp2mwpzWYovVW5aqHQAiIlcpCktK5Ikn+BRnio3WpbdXhf0zEzMotJzLYeXE6VARqlf9cQGAMDqv77AtZ252uWYdUn4cqr8q/s+5p8bJ2QKCaUuwy2MSPUDAI7EuVIfv/MdAECy8dKcbTGZ0PbC7Z7lEwqBGPpIpX7GN38NAGh+6iuD2ZySJALz7wIpdoIgiJBRsGK3vkXyzWkplbrMh5oaxTPK6HlNnY5PaTuVS6MIDcArkElE47bPKWRnOnI7j6rKoYwK6AV7aAWnaqVS/+KTPPzA6zeLHK5a+GG3c7ltr0vwc6YTYhl6o/ty9HSFfUFZvfgtlXrs/VdU2YOzvggAGFfFHxU5uuqr4P+bMhHEDABifd22ffK+yOBrTouOgsipSnijj4gHE6nUb3xmMwDg5dvnDWZzjIk4hE0ZaBjMvwOk2AmCIEJG8bxiNHTl7mRDlUo9dvQzAED69Hm2/Vbk8W4hACLRzKWpUAKyjHjrxrRjrOgBgZzaoOzvFtUKABExMvCaG5BK/ap/5aGKX/vW+Q5ncD6fyeIgr3257OaHG69R2yaJgGMp8GBjPQmu+qqO7wQA7Kk4XZUdU2Xfx6L8Djvl9HS7FvKkMcNPwKqhoNR1pFJv/JtXAQCf/PTawWxOTgZTqedDabWWIAiCyEkgit1P8Cgv9S1t6lKpRx1sgxHtt46yd1sKuKkaL3WovGA8yqh69JAIMmSBx1te2tSlUr/ckv39TeFB4CdwmttnL9zqHW1J3pGSYYGFJ02ijF/rvkqu1E/f864q2zXzEl6vUOpo4SMvCMUedCKP4UwYQssCGaX+k/X71bZ7LjxtsJoTGkixEwRBhAzq2AmCIEJG8aI7uoQUyBR2jz8ukaaZeMsnalvb6FkAgEotT6m+kMhpn/7ZafgfSfKQAf2xClt9MtmSU0L2rFO6XI9jWcGblqXW0hXy+Vu5SWpURcy1Pj94uZzmql93kRwvnhzrYqaE+C0nSyPid9mJPeKEmZuXFFmsyo/tFMdw084BHl4eU3r2ZsqKfWUnMtuGEzGRYwDIRC8NG1bzy3de3QYA+NG1MwerOYES7e0EYF+w6EZbL+9oahOFaW5S7ARBECGj4GTWETg778vFP7pLoVKJPtyH2sfMUn/XHm0CAPSPb7SVMcnWY7J8PSmUupyElQMDFcNcZCwCMrHQ9fM4kUsdW/dLpf4X//oBAOC335hvfB4/ZXRMApLlVVYk6Urv2pLZeB5X7MmGqbyMUDXjqvnIIFk9NXO82GeN4T/QDGYy67CqdDekUpfhN14Xi/pKFROlLilUqUuMa3nsscdQW1urfiZPnhxIAwiiFKDnnygljHOeOimWyZMnu+Y8lfgJzuWnTLSrFQCQEm6DJjkzc9XpRaELZwpR0Ff8bD0A4I27LsyqK1fIBpMAal74yYaUq01WosJ98mCKj3oaKvngsfLTt3iB8Z9XZVO1PJDZqfdXY/SfLBmUnKeuz/+hg6iprQv0XDIcBosWd/3gUFgmb8LFP+CZyN753qJBbsngEk32oL29HWNPm5rzO2D85CQSCSQSidwFCSKE0PNPlBLGil3HKUu7H/u2n1C8Xkh1k4rwd5TTwo18lKrEyxau6lvzLADgyIKlAICxVeZeLH7u2b81HQMAfGViRjkmayc51mdyrfn8b3T8LEoz8RLyanebYYb2gUA+/3sPHMLoUe5L9v08y4Q3z27lOXpvmeMe4C7stBt+B4b2GIwgCILwTSBheyVeyRFMvFZksohElP+OOPiD55OEIte+fO3nct/RC7hSnyACZiUxyeUI72Boudollfqh8vFq21iXY/zY2HP5tztt8wrmlssLyakePVSEn//5YBJ3WjhhgZR6cAxnpe4XUuwEQRAhgzp2giCIkFG8nKcu+7yGpgkxrGVaTY4TlvKzcAmLITd+hsUmZgq5b4yIhpgUkRDlIiZ9AZPJ+byQE6VW84sMdRDXshT5mTz12l7IZGw+Zggvd9WhbJIhiKEEKXaCIIiQEYhi95oUcz3GGgRMHuOyUMLE/S4fVWfSbhNXQP1zWij1WMdhta+zkutsOUFc52PpsNe9lEr9RE8KANAQ5ROs0b5TvC2W5egyk5T8nZWL1GHBSq4MR07K2o9rpMmIgFwGibATFcEH02UVwdQXSC0EQRDEkGHAcp5KVHCtVCZPqFSITrlCc2HiPpjLzdHPIiE/nKrKWMNrTnFXSHmt6cQ43+fxcgWUSv1gP7+HsShfJWkNHyUXb+l5XFOqQnvgNqf2+Vl0VUiwMqdrJBs7kYsPD3cDAM4ZVznILfFHUEpdQoqdIAgiZASq2E3UssolanlDScWoe7YU6mWT65hCbbe51HZvMlOzDNwZ2fcx/6Mx92ILEzu03CZt6lKp7zjBlcvYqhFZZSWZ+y7mO1RylMxjkY9KNhmF5JorybWtVOns5/d4hFPGFqJgpFIPW7IOv9DTRRAEETICTY1nst8pZZ4euMuPJ4pOofbyIG3JI995Wv2duvKbvGzjBIOa7ZiECZDeL9KmLpX6M1synjm3n2sfJajRk3y/y7kOj7YU8r8x8SjKp95SgpT6wCCV+rz7Vqttmx+7crCaY4waRRf40NNTRhAEETKoYycIgggZRXd3zJr46+eO+H6W21sJaiESryzbLOSnrlzxzNNXfFNti4mMQagY6VjWa5tJlEQ5hIsq8wrHan7pEBN3NcIcoB/jJ6OSyT4/pi+T7WGaRCUGBqv5ZXMLD/Uxb3x+fc9AUKgJRkKKnSAIImQYK3a3LO3WF4xXNh19ub3TPje8lq3nE/ddfXaYLDRxv3NDL2t9+zKh1CNCuacr3LNOudXrNVmrq28npFLvEm6YVZF+0VD7wrB8A3AFObHtZ6HZQOD2/BOlg1Tqq3eeBABcOa1u8BpTZIwVO2VpJ4Yz9PwTpYRxzlPXLO2WnKe+FgtZgoClxPvFze3RJDNTVv0e5zZRwH4WM7mdK+oQBGxE7wkAQLq6wbFNXvV5KWr9WCe3KbUtzcM5nOjnS8JkOILjab64aVSFSRDk3G3xKuuG1zWa5nssBm7PvzXnL0EUG9PvgLEphrK0E8MZev6JUmJIBAFLRXl4gcyCGXNMlHUuBex3nylOQcCOxPm2+jzqMxk16GECGBwWggmbekOKewnIwGGnJXkbUxWZnK1+wuvqbdHrcEKvj0L0EoNFJJ0EkEneM9SIsLRjuHMnyCuGIAgiZAQSUsCPjVqqRmuYyhhzLut2PqcyJoGyTFRhPvZtN3VZ3XNC/d1SzpX6+J3vAACSjZe6ni+vkMHqTR61f7b46OvXfULY1KVSX9tRCwBYWGveJpORUaE+6WEKKWBCWes+AEByVGlN0MrRuJdn1lBkqCt1CYtEXZMR6ZBiJwiCCBkFv6Lc/NhNykuU54awu6ei3ObrtCIylx+7Xs5aJp8RQT62XuWzb/F8GS1+S6Uee/8VAMDhxmv4/srcniieo4aIliRDhN51ula5TXq/SJu6VOqv7mhVZW8Yw+3wqdqM3d16XpN2et13tzJ+vG3CRqkpdUmpKXWJrtT3tverv08fWT7QzQkEUuwEQRAhgzp2giCIkBFoPHaThUReE2lRzZXHzyRnPm54hcZnd2uDienh4KwvAgAm9QkXw8pJDqXN6zM1UVnLuH2W5hcAOFDOg4jJUGLSnbLM4QRJbV8+LqhelOhInygximV+kZO0QPEnakmxEwRBhIwBD9vrpChVGeECqb9t8lFq+WZAyuWyZzIpa6Lux1XxW58CV+oyOBiQHSDMxJXT7TxO5Gq3daJUKvVuLXCYPsENZNzdALurZT7/C1Ln+VEqrntyoY2p+14YGMj/yfC5qwRBEMOEQF8hfpSkp/IVbo9JoQqdgs/nYy93w4892mvhkx+3vFh3GwCgJ8F9DBNlcaP26LiNiJzqcmufl0qWNnWp1DvS/JGpiXLF1ZvKHBVX/yi7XnC6P4UEbyPcGepKXVKqSn1HK++bpo+K5yhZGNbQAfncq9K8uwRBEIQrRXu956NiTZakF6LaBjrAlJetva+CK/Wq4zsBAPsqT1dlxpfZyxaS/MMEr9GU9HCRNnWp1Pd0pAAAn+/bp8omR0/jZWUgMuZu7/UTKCyM9KUyVxwPKh8aUXSkUj/axZ//MVXmIa79UOiIhhQ7QRBEyAgkpICJN4gf229EhJWNOFRo6r9uElwrHy8TLw8aP14xZSLZxZ4KrtRP3/Ou2ifDDrjVkU/wMq/6TMpK7xdpU5dK/d2eTGjihTnqMVljYOKzX8pI22k8RpqqlJFK/VR/xhZeXT50/qdDpyUEQRBEIFDHThAEETKMTTFuWdpN47Hn2m5DuDvG5AKXiHsz/ZgTTBYU+YnHXkiZWF83AGBMFZ9E7Zp5iSrjloAtn5jwftwzndotwwTIxUdyok9OlC60lO0Qw9IaOfkrl1CLydN8Q07oZQYDt+ffL6Xq5kc4YzW//OEo/06fPaZysJqjMH7KKEs7MZyh558oJSKMMSMh5DdLez6KXak1mQ9VqhuLu5yfyc1cFFuxe7VFertVCHdHFs24TSXrpziey2RkZDIxnGui0vF+aBmZZPttrpGwu0JKL75JNe5BlUwmTeU+0wztxcDv808MX1789DgA4M/ObMhR0j+m3wFjUwxlaSeGM/T8E6XEgK0/9uOqxoS7o8qhaFBfIcv7C3Wjc1tA5OVGKa9NKfWWzzLtEYrd7ViveoNwF3SsX7MNy8VHtuOETT0W4de0t40r3NOEYs83gBpBlBJSqb/w8VEAwE1njRnwNtBMDkEQRMgIdIGSFRP7rVedQGZRTNpykLTtxqPO5/HjoRO0Pd6P8peT6dKeDotKLzuxh9cnpj/6G6b6rt8PJl5CWWWcgk2JbZNq+Eep1A90Zqv7SSO0nKzChh9x8hrRErAQg4clEoJjcL4giPZ2AgDSiRHFOcEAIZV6aw+fc5I5hgcCUuwEQRAhIxAbu2fyDO1zPn7KMWQUmwosBXvYTF8hCzzO67YvHy+TvJfzC9Wa3rWFfxaKPR/7udcIJp+ga/nY+yVbD3eqvyeNqLOXjTiH+rXuC0NIgVJnIOKVDSWlHkTiEqnUj3Wn1LbRlcVV76TYCYIgQgZ17ARBECEjEFNMoVH7ck60WofpMvJjjjZ5mUG8zEO5JlS9Jha96tDLet2P5CixqvE8/jsq8qHKXKi52mWt15ebqcM2P0v/3f7XcqLUan7pErEKEmJsLyfl1Apty4SpGg57tJ0gikGQGamKbX6xQoqdIAgiZBRtgZKbKs4rZns64y4n36BS4ckMPyYqPJ+Y5Vlt8ajXJKiWm3K3Un6MhxlIiknTg6kqAMAEj2NM6jWd5C10fsxklCaV+v4OnktVukbKa+9tmKbKRuUoLd0JgggzMt+zXKSZL6TYCYIgQkbR3R31Mk77c9nY+y3NlO8xmYGIudjcnRSraRAsJ/woYK9j9GOdyibreValiFio0VBpt617jR5yncdpm59gY3oZI7dSYS+3ujTKEZdU6t1JXiY26jQAmcVpNgK0dxLhIprsAQCkyyoGuSX50SlCXo8oL0ypS0ixEwRBhIyCJZAfLxLAbGm+fkzc4fWTivI3m77Lj9r0Ip8wvSb1myjfA13897hqrtQrm94CkMmF6me044XbsfmGiHC9NocwAdL7RdrUpVJfd4iPxBZMyigXtSgmnR2agCgOgxmUra2Xq9fahLnudFPqAxECoRCkTV0q9f/3EQ8c9tXZhQUOI8VOEAQRMgINApavt0ZOVeDgFaPCxkay065Zz5trn2lb/Njh86nfypSevQCAZPVUvmH85231+hnt5BMq10/IAj8iyFZW2N2l94u0qUulfkB4ywCWYGJd5Mk+UAzmnfaj1HMRtEqPOMwXFYLu/VKoUpeQYicIgggZ1LETBEGEDGNTjFeWdpMJy3xcDCV9Du6OcpvMpulngQ4KLONGPlEdndDdHVO1E2z7vSZP/dxnk4ntQHAYvsowAVFtKCqHzqdZ8qRKV8gxVe65U4uN1/NPDB+CMsEUG+NWUpZ2YjhDzz9RSkQYY0YCzS1L+2GXLO35TND5iWcusypJhZfPAignCl205NYGP+cpO74bAJAULoCxT9/ln4W7o9M58nHPzCeUgp9Jaq+6VFvkxLhcfCQ+7z+VqW1sNVfq7zQfwvXzPp8zQ3sxcHv+W1yef8KZfEaFxV58FPSEaLH4xzW70dvVieVfXpDzO2BsitGztMv3QUdHh2N56tgL6Ng7uAkmGePD/Vgnd2xPOgz/w9qxd1g8YCpSvGPv6uTPmqEWCRS/zz/hTGEde1/g7QFKp2Pv7epEbxfvG3J9B/J2d5QP9PQZM/KtgiDyoqOjA7W1tYPeBgCYQc8/MQjk+g4Ym2J00uk0tm3bhsbGRuzbty/Q4agc5pZCvaXU1lKvt6amBh0dHZg4cSKi0cFVV/T8U73FrtOpXsaY0Xcgb8UejUYxadIkAMDIkSOLYmcspXpLqa2lXO9gK3UJPf9U70DVqddr8h0Y2kYlgiAIwjfUsRMEQYSMgjr2RCKBBx980OYtEASlVG8ptZXqDZZSu2aqt3j1DrW25j15ShAEQQxNyBRDEAQRMqhjJwiCCBnUsRMEQYQM6tgJgiBCBnXsBEEQIYM6doIgiJBBHTtBEETIoI6dIAgiZPx/FhptCZiVez0AAAAASUVORK5CYII=", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAADJCAYAAAA3tRlxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZXklEQVR4nO29e7RdVXk3/Ftr7cu55OTkTgiEROQaIAgCSipSBBRv6Peh0JZXCGpr9dUPbe1obbWgvqN839tibd/XYVERtTpGrUV9B1bKRUGx4SKXkCKRcE8ICSEh5JxczmXvPb8/5nzmeuaz5lp77XNOknN25m8MRjhrzzXXXGuvtfZvPvP3/J5IKaUQEBAQENA1iA/2AAICAgICphbhxR4QEBDQZQgv9oCAgIAuQ3ixBwQEBHQZwos9ICAgoMsQXuwBAQEBXYbwYg8ICAjoMoQXe0BAQECXIbzYAwICAroM4cUeEBAwLbB69Wq8973vPdjD6AqEF3sHWL16NaIoyvx30UUXHZDjX3vttXjd6153QI4VMHOxZs0aJElywO7LqcI//MM/4Fvf+tZ+P86h8ANSOdgDmGm46KKLcNNNNznb6vX6QRpNQEAW3/zmN/GJT3wC3/jGN7Bx40YcddRR+/V44+PjqFark+5ncHBwCkYTAABQAaVx5ZVXqve85z3ez+666y5VrVbVL3/5S7vt7/7u79T8+fPViy++qJRS6tZbb1W/8zu/owYHB9W8efPUO9/5TvXUU085/WzatElddtllau7cuaqvr0+9/vWvV/fdd5+66aabFADnv5tuuml/nWrADMXu3bvVwMCA+u1vf6suu+wy9fnPf95+dtdddykA6ic/+YlauXKlqtfr6qyzzlLr1q2zbW666SY1ODiofvSjH6ljjz1W1et1dcEFF6iNGzfaNtdcc4069dRT1Y033qhe85rXqCiKVKvVUs8//7y6+OKLVX9/vxoYGFDvf//71datW5VSSq1fv1719vaq733ve7afm2++WdXrdXt8+Xyde+656uMf/7i6+uqr1Zw5c9SiRYvUDTfcoHbv3q1Wr16tZs2apY4++mj105/+1O7TaDTUBz/4QbV8+XLV09OjjjvuOPXlL3/ZGbt8ju666y6llFIvvPCCuvTSS9WcOXPUvHnz1MUXX6yeffbZKfleDjTCi70DFL3YlVLqz/7sz9SyZcvUq6++qtauXavq9br64Q9/aD//t3/7N3XzzTerDRs2qEceeUS9+93vVqeccopqNptKKaWGh4fV0Ucfrc455xx1zz33qCeffFJ9//vfV2vWrFF79+5Vf/qnf6pOOukktWXLFrVlyxa1d+/e/X3KATMMN954ozrjjDOUUkrdcsstavny5arVaiml0hf7iSeeqG6//Xa1bt069a53vUstX75cjY2NKaX0i71araozzjhDrVmzRj344IPqrLPOUqtWrbLHuOaaa1R/f79629veph5++GH16KOPqlarpU477TT1pje9ST344IPqvvvuU6effro699xz7X5f+cpX1ODgoHruuefU5s2b1bx589Tf//3f2899L/aBgQH1xS9+UW3YsEF98YtfVHEcq7e//e3qa1/7mtqwYYP66Ec/qubPn6/27NmjlFJqbGxM/fVf/7V64IEH1DPPPKO++93vqr6+PvX9739fKaWfsUsvvVRddNFF9jkaHR1Ve/bsUccee6z64Ac/qNatW6cef/xx9Qd/8Afq+OOPV6Ojo/vjq9qvCC/2DnDllVeqJElUf3+/898XvvAFpZRSo6Oj6rTTTlOXXnqpOumkk9SHP/zhwv62bdumAKj/+q//UkopdcMNN6iBgQG1Y8cOb3tiSgEBeVi1apVlqOPj42rBggXqjjvuUEqlL/Z/+Zd/se137Nihent77YuPZob33XefbbN+/XoFQN1///1KKX0fVqtVtW3bNtvm9ttvV0mSOMz+N7/5jQKgHnjgAbvtne98pzrnnHPU+eefry688EL7o6OU/8X+pje9yf7daDRUf3+/+sAHPmC3bdmyRQFQ9957b+41+djHPqYuueSS3OMopX8Qjz/+eGc8o6Ojqre3V9122225fU9XhBh7hzjvvPPw1a9+1dk2b948AECtVsN3v/tdrFy5EsuWLcOXv/xlp93TTz+Nz33uc7jvvvuwfft2tFotAMDGjRtx8sknY+3atTjttNNsfwEBneCJJ57AAw88gB/+8IcAgEqlgssuuwzf/OY3ccEFF9h2Z599tv3/efPm4fjjj8f69evttkqlgjPOOMP+fcIJJ2DOnDlYv349zjrrLADAsmXLsHDhQttm/fr1WLp0KZYuXWq3rVixwu535plnAtDx/+OOOw5xHOOxxx5DFEWF57Ry5Ur7/0mSYP78+TjllFPstsMOOwwAsG3bNrvtn/7pn/CNb3wDzz//PPbt24exsbG2ooOHHnoITz31FAYGBpztIyMjePrppwv3nY4IL/YO0d/fj2OOOSb38zVr1gAAXnnlFbzyyivo7++3n7373e/G0qVL8fWvfx1LlixBq9XCySefjLGxMQBAb2/v/h18QFfjxhtvRKPRwBFHHGG3KaVQrVaxc+fOwn3lC9b3wuXb+H1Nx/HtI7c/+uij2LNnD+I4xtatW7FkyZLCcclF2SiKnG3UN5Gkf/3Xf8WnPvUpXH/99Tj77LMxMDCAv/3bv8X9999feJxWq4XXv/71+N73vpf5jP+AzRQEueMU4umnn8anPvUpfP3rX8cb3/hGXHHFFfaG27FjB9avX4/PfvazOP/883HiiSdmHraVK1di7dq1eOWVV7z912o1NJvN/X4eATMPjUYD3/nOd3D99ddj7dq19r9HH30Uy5Ytc15Y9913n/3/nTt3YsOGDTjhhBOcvh588EH79xNPPIFXX33VaSOxYsUKbNy4EZs2bbLbHn/8cezatQsnnngiAE12Vq9ejb/6q7/CVVddhcsvvxz79u2bkvMn3HPPPVi1ahU+9rGP4bTTTsMxxxyTYdy+5+j000/Hk08+iUWLFuGYY45x/puJap3wYu8Qo6Oj2Lp1q/Pf9u3b0Ww28YEPfABvfetbcdVVV+Gmm27CY489huuvvx4AMHfuXMyfPx9f+9rX8NRTT+HnP/85/uRP/sTp+/d///exePFivPe978V//ud/4plnnsHNN9+Me++9FwCwfPlyPPvss1i7di22b9+O0dHRA37+AdMTP/nJT7Bz50586EMfwsknn+z89773vQ833nijbfuFL3wBP/vZz/DYY49h9erVWLBggaPrrlar+MQnPoH7778fDz/8MK666iq88Y1vtGEYHy644AKsXLkSl19+OR5++GE88MADuOKKK3DuuefasM4f//EfY+nSpfjsZz+LL33pS1BK4dOf/vSUXodjjjkGDz74IG677TZs2LABn/vc5/DrX//aabN8+XKsW7cOTzzxBLZv347x8XFcfvnlWLBgAd7znvfgnnvuwbPPPotf/OIXuPrqq/HCCy9M6RgPCA5uiH9m4corr8xIpQCo448/Xn3+859Xhx9+uNq+fbtt/+Mf/1jVajX1yCOPKKWUuuOOO9SJJ56o6vW6Wrlypbr77rsVAPWjH/3I7vPcc8+pSy65RM2ePVv19fWpM844wy5ajYyMqEsuuUTNmTMnyB0DHLzrXe9S73jHO7yfPfTQQwqAuv766xUAdcstt6iTTjpJ1Wo1deaZZ6q1a9fatiR3vPnmm9XRRx+tarWaestb3qKee+452yZvEb9I7vjtb39b9ff3qw0bNtj2Dz74oKrVaurf//3flVL+xdOrr77aOcayZcscJY1SynmGRkZG1OrVq9Xg4KCaM2eO+uhHP6r+4i/+whnvtm3b1IUXXqhmzZrlyB23bNmirrjiCrVgwQJVr9fV0Ucfrf7wD/9Q7dq1y3tdpzMipUIx64CAQwF33303zjvvPOzcuRNz5szxtvnWt76FT37yk3j11VcP6NgCphYhFBMQEBDQZQgv9oCAgIAuw6Re7KOjo7j22munfBFvJvU7k8Ya+p1azLRzPvvss3HNNdcUympXr17dcRhmpl2HQ+E9MKkY+9DQEAYHB7Fr1y7Mnj17ot3M6H5n0lhDv1OLmXbOod/91+90G2sIxQQEBAR0GcKLPSAgIKDLMGFLgVarhc2bNwPQ04WpBPU3E/qdSWOd6f0qpTA8PIwlS5Ygjg8uJwn3f+h3f/fp67fsM1A6xj46OuoE8Ddv3owVK1ZMZswBARPCpk2bcOSRRx7QY4b7P2A6od0zUJqxX3fddfj85z+f2f5/1qzDR/7+bgDA2v91qd0eqZb5H/GrkredoWl+amLjHZTsS38FVbUGAGhVevS/pm1i2tKvFLcjindr57fWrEUAgIb4KUtY49q2DQCAe1frdP8zfvoTZ0x5++mD63NTBedGyOzqazOJ/iaCiayi8+Pm7R95Pm93/r5+dw8P45hjj8048B0I5N3/69ZvwH9s0n4nV6w8zG6n8UeNEQCAiivOv/x+qqiG/h/zfUfGx0TVjHql1Uj7bYzrz6q6atdoSx8pMQ8LPQ+cr/WM7wYAjNf1das29ZjkMwQACfQY4j3b9WfmmaGxNT3RW9q9Yvbl47VtEv3c2veCAfUXsy+crg2NczTW46yaQ3vv9Zz3StNzbvLYcnvUHEu7NdfI/k1tzL/8XVJp6f3oHMcTva+9LvzczX1A+9t7wIy/wa4z7d9AjOHhYZxwXPtnoPSL/TOf+YzjbTI0NISlS5eif9YA4qq+Afmq7ZS+2Kvp1aMbuqMXe6QfvNaAHl/hi33fLABAf1Jxzim82Nsfd3+/2O3Lso3V6/5A3v0/MDCA3n73XgH4i12/0Dp7seu/Va1Pb3de7Prloar6/i/3Ytefjdf1+KpNPzkC2Is91rMTemYOzotd7xNe7K3MtnbPwITljiTD2bp1q72hR9hV7BFvvaIHuZPHlPZrihd6mT7GzHWtFrwnZT/x7pf18WYtzHwuv+S8v339l3mR5Y2pzLUr+lKLxjeZsbQ7Xifw9Ts8NITDFi+eFtJH3/3/q03D9vNzlmpGZRm7eUHI+xZI70tCDe5Dzl+mtJ/tB+4LzUdAGubNXY/0TKAZuXwuUemLmD6jbZF5SdOL2YF4mdKLpxJlz6sW0Y+WfvnRj0rRS1r+iFSa7uyHH9u+pMWPZCNOx23fFcp90dpzNi/mJtuHrmfNDMU+4+IHivc3Zviy/B6dH52W+Ix+1D0/iiORJrP1qImhoSEctuTIts9AUMUEBAQEdBnCiz0gICCgyzDpCkoR0mlJT5L+Thx5yZcAAI/889UAgEV9CQAgHn4JALCnb5FtS2GbvJBGvG9X9sC9g4X7+MIVNJ2S01U+LW6IqXLLhGDWvOHNAIBz7vyBbdscSBfLisaiNxavLTjjFVNRGcv3TV87iaeVCcHkbfOeW8G2dvsUfSb7nY42pFFzzI6Pwi8AcP0aXXDiY2/QyoW+Ub2AmVQpZpveBzUbkNfbxlr6sazSGhOf9rfo3qg4+yTUqyc0WVEmXm5ivgkEWGgmEfcpxYnrO58HADTmLUuHQqEXEyKJEjcebUMRYOsGMtbuCSElJhzRgK6UVI/dbz4e3WP/f6xHvwcoDu2EaQBEbNf0GLGzD12zUXO8HnM+AJCYkJG9T2ls5po56xPmXBRtjPPj/RUZ4qKYO4WF2BdYp2OPjdnF83YIjD0gICCgyzD5mqeq5WWhxNT3NYhtaJ4w1KMZcA9bBm/LyNiCQqvHlfmUYZ8ZlY1YBVfs941+bC0DMOd2zs9u1n0wtiCL1BUySnGNChcUBQubjFqlEyVKmQVX2XayqpipmBEcTDTiml0g5SCmfttTuvzhe4+bAyCdhVXGdtu2rfosZ9+aMqqMlod3GaZbi4WSw6hlmoY1OyoTYuqCJdO4G4xpx/QZKUPMQmLTMHXf900LocS0aTaRsHueFo+tCkSo2DibjQ2LrZhtTXNURSqTJLsgSh3Z8ZmxJL4ZcuSuhBILrxjW3Ip6srvQtTNtKnR9kuxCa91K9NwngrNwO7MQyrdYXBc+vlatD61adnHVh8DYAwICAroMk2fs7BeR/5pTTJ145z7zc9ZvtIa8LcmiauJnxv6a98/PHHZcSBeL4tsyRl2kC09lTSI+ZmLt9C8AJGa9oJUTay9CmdhyXpvJSkXz9imKgZcZS7v+ykhFff1O5xh7JWLsi8V3KaZOTP0nT+t1ore/di4AYLSasvQaMV0TX42FtJBL4Ei+Z9ef6H8qLisfY2zfyntbLju0TJ4zSXtQNw5N/fG2pNu2rNXONLP5F3aYVq9dcc6txq4djSFRfrlgUzGNN8lJxTWzEkY+O6d+DNuma0nnVBnbq7dX++wumTFQ/7Qv+24SkaNgZzCetrI/m4/jyQWw+6uWX2bpQWDsAQEBAV2GSTN2BXgVH6R+oZg6MfUTPvZDAMBvvvJ/27aWqTfdBAH6JR1lATjKsKN9ZJyQGIyzAm3imRTLjO7+jv531fv0OdTSX+jI/GqPJjqbtmYGUXtxHQBg92GnpOM2TN3Gi0W8rLrjOdu2ObhYj6FNJhtH3iykKK4tM3H5L3w0blLJ2fnyfXyKik7YcieqmHaJWnz7dGTqFq2GZeyOQsKoXyimTkz9prVbAQAfPnWBbWsTWgSbtfFodji6PpbZGbUJsU9SW9TZ964Mw91nVB+95jlLbPZjNgFKqkvqzWw8PrO2RizT7BubZwlInzFSfdCXWiElCT9H0Y/MZM2oepDe25SRm3gSq8ZtZMDEyRWpb8zsoZqNrUNcBxu7p+vL+pf3qb3eylXSAOkag42xe87JDsEmfMUYU+W4eGDsAQEBAV2GKdGxw/OrSzp1Ur/QDxQx9dP//Dbbdt3/fJve3/z6yV/kClvil+x71KTbRoJ1FqkO4rPerY+z4xkAwNjhJ9vPKL5WE+x11yLN1AfW32nbNk++UI8bblvb15wl6biffVCP5ZhVTlvvL3WO5r0wnm32odV663HBlQmCqct95GzCaZszhjKqmDIoY8MwHVUxKq7Y2WHFGaCrfqGYOjH1v/nVZtvyL9+sFSejTX3n95hZaIU01JzpUUxZxJ1t7No8Q8qTml6z/iKGoZp1gEpPNjWdvofYHI/S2ntH0pySlsklIaSxZNfPhh8rMc8isWeQdQFnvoIljyl9Xer0kBeobeyMhvTyzPelZvcTVgXCPsFR0tCsxuwT0TqIaRIzRRTNZmjtoWaux5gyMzJ2f+w1A+6tmNm9md2Mm/dPxK0Q0DkCYw8ICAjoMoQXe0BAQECXYdLuji8xdzsOGZ7IkzQCwLa9OtUnlUj6++Db5CJh0T4EmahEiEezYZsysr5kr04+afXNzW2T2cfYIzRzLBGKUBiuUK5MiqbF8lx5P3Q9pOOgT/7WbmG0DMrIHYuSmYamobvjS1texOwBfc+MscgmuRlSyISm+XahlD0Htz3zKgDgbct1P1b26FnMRo4zoW1rE4uyDoiQn4m+gPS+GTMPWK+5GVIJXxrasP3kWOX67qNYyBO91rnkFimFAKIvgJ0vjUGEcbitrg3TCBdHKWn0yTQlpGsnH5eUJFp5KRsLec2TZQO9z3zWxPzaDQ0NYdGRy4O7Y0BAQMChhqmRO3ogf12lpJHLkIipn/Lp/wAA3P0/9KLkvB693ecxbZmoZCHEkO6/OR3jWf8XgHRRxKaAU9ECtngU58gn7QyBeVcTU6+8/BQA4NU5RwMAZnkM322ylWHqJJ9sztULZ3whqp0EcIxdEEpftswwZwHah1gwdNsrW3CKSY4lzZAK+i3D7osSyvIwHWWPKorZYjf7QLLYyJU00kIpkDL1b67bAQC46nVGRmtY/t5W2tZOas39nsjniVLTmXHeSI++TxNiqMJSg7PcyrhexFO0iEfPSuJZWDfjI4basNYF5ImeNo2IDdM4zb42GauVleVWSCJJi9NKyCCRnZnGI0N6H1NUpCLPle1DC9sjlX4AQJ2EHlwqasYtIwR0rjzeIYtlFIGKh9T4DAhAU4lZBNJ7p5n0YDxx2+chMPaAgICALsOUyB2LYmoSNq7l+YyY+uuu/CoAYOP3Pw7AjRNLphebODf9ipP9wPhZl9h96q9qC9Xm3KV6X9M2ihqZse4zsq4eQSEtq42yl2zHoGHqt/6jbnPxJ52x+sZNTD15+WkAQOuo0zP9EiSr5TamSnDdIjli2f55BRkCMQAuHwNcmRrEZ9aitcQY5Fh8mI5yxwiwlrS1SsqTyHqXDL0sMzX3Tw+XyRn2Rkz9f/5qIwDgz1dpuWyNPVjWEEskL8VCBtk0s0kAqFHiX1Rz2libYD4DrOg2VcNmx2t6NmHtAxhbpnuB2Gt1j642Ntaf2m5IWAYt7w1uKZC460VkzkWywSrrz8bfjVyQJL2ZWQnfx5wLrafV6TwUmY+x2QPF+w1lT8wBM2UxwczPRCk86sOZPcQUdzfHpOpQ5uNRlb4h68yzzHdcH0q/2GWV9qGhobK7BgTMeIT7P2AmobQq5tprr/VWaec1H4t+TIpiqXkD2DSsTeWXDqS/0e1UJOMe9Y1l9yb+pkRChmMZKpUIIlmojAFXvEfHSn3mZXn78pioTPyYSFLQZNht0TnKQuD8cx4XBLKKDR/KqIJsktjQEBYfJFVM0f0/SGPxxGZtzVBhEMWTmeyzIZKKNuzS/R03N2sNS7DJQIIB82bEyKnmqVTFeFUgQjGTqeMJ2GdC1iaVY/LBXoeWJ2ZMhb+FhbBvBiiVXXLdyKvEyqmpKuul6o3mmDQjMElXlNrPGXTisQ7gY+CziLzYvU+pZCMMrcbU1zz9zGc+g127dtn/Nm3aVHbXgIAZj3D/B8wklA7F1Ot11Ov1zPZ4bB+SffpnhzNNy0BJV23YKxl6cZsAWXmdPiKmfvpnbrdtH77urQDydeZWB8q06Q0TJ4Rh6pYZkc0oYyzJ4z/X+6x4i9Omum0DAGBk4XGZcUuW1DLnWjQ7ycS12bWrbX0cADC08EQ9bBFYK1KX5B2Po90+RXpzHxOybcmK1fxt3WKJlfgsVM13MS713Yy9jk2Bu/RkkXf/A0jvcV+JOQNbzMH8zddq6HqR+oVi6seZMPm3171k2/7eSWlJSSBNnadrbFVoLPW/TvcWDUlov1uM31VYWTggjZ/XbNw4/WxvU4+zL3bVKpEwJgPYfWOOTSnzYybuz4Vk1jws0f1Rybq6sNsF0ri1VOJY5ctoGjKjmTqxZVl4JL3Z08HY9cOKUPOY/p1niJRK1E3LnQHwGQxdhwqVDCS7D5GfAKTPUTOqONuLEFQxAQEBAV2G8GIPCAgI6DJMXu7Y2AdV1VM939Re1iglP/Ui90U5zafwCwC84VodKrn/2rcU7sMr1CBvgcKzILr3+PMAAP0m9NJYpEMvjYXHAADqT9xt2zZO+F34MBHJHt++a4EOwQy+vB4AML54hdNvUaikqN+8MFDe50XoZDHcTmf54iq5EMa0KCf6YPI3+mxfiXEdaESArTea1Njj1HIXwWTlowqTO9IUnZKPSNJI142HX771qPZz/9Bph+v+TNgiSdwqTs4CPIUCYjdURtv54jY5fFIYoSnCaBHzWK+bJKBRCiFRtyas4PjIm9CoMs86MUq6J6IGqxJlkqNoEZbkvcqEbSInFml6EolPVj7IHE1jYQNAIRkruBB1Tfn4aEGUDDLJEsBxojTjalk/fQqtVZzjA2mIq0FhJnE+Pol3HBULETgCYw8ICAjoMkzaBIzLHbkUSwrppQyRH1RW8MkMkv0/7ffSXv1Luriv4mwvWhjMMw4rAi0CSwkikCZHUTJIHgP2oSilXm5rdxzfsSYqPZV9TSa5qOh60Gc2XTzPmAnp4ume4SEsOXx6mYC9uGUr5s425l3ILrqRGVjG7MpX/9JWNDKbzccNVpqJZrx3GOOwi147xzs+p3oRLfxRQpKi5Bjdby/JIIGMFzoxd2KsTpINWQqYGQslM1k2y+WOcpYs0u554puVNVKCj7RN4AvTsXj+ZU1Qdhyqu9wbiztVMPUx1oU1xhN2BlZCza4zeamTlNW2EbMVAOkCtjwnYYDGUYnKG+EFxh4QEBDQZZgSHZmPCUtGamWIJs5ElY+A1MiKfsUsQzU1RX0gpv7Wr9wPAPjxR84CAPRWshxz54hmJGQqVnnmPgDAtiVnOtuBbDyYmLqPEZAJGO68Ufd72vkAgMb85bnjtvLGNrMUDsvUPVVxMgyLthf0VzRb8P1dhDLsvmgslt0Y5mKZBpN1UXpaJzOtA4UkAkZbemD12CPnNEhN67JGUSRhI0Mv+VAmXCZnGCIx9ZvWaink5afoOLyVEbPYckt8AfVxtwpZQ7Fr3XDtZCPB1B2mLZk1Hc/EsH2zBpq5UJzAzta5hFFYU9jKRGMe5ivro4r1BEpKBIAeSkxs6sRHOzukckhm31qcnXlZGaK4lq0qq5fsfpTONCpZaWQkEp9oXaHpsze3z33szlYKEBh7QEBAQJdhShi7j0nlxVvp1zfy0bjYXdn39SXZJjH1S27UNUVv/ciZmW7nYw8AoAX9i/2SYeoLn7gDANA89aL0XCBigeK4PkvO6uu0kgY7XzQHXO7uk3MueW3yGLUvOSERbSeyYNKJwqUTdl9m3aOdQqfd/gcbCmnc20omwJJUTPwZVJuTEmfYGkIsZy2GNRPz5Yzbql8MiKl/9de6hurVbzwys09FGHpRwp6dITgGXK4JmK1bKtL7fZD1iBuMzVoTPWmpYK1us/e2jJfbcbNtVukjUvHp75annivBfm1kNWJmOXytxKpixPdGaw78OtN6iuXLNmGLTMbStk2zY1TRx6xG7ph8iU+Z9YMCBMYeEBAQ0GWYUlVMGfbWCUhJw9ON85gdbacyewv7ssbAPoMwAIh3v2z/vzXLtRyVlsRFKfqyvybrq1382WFlkb9Nu77Kts1rU+bcijAV/RbdJ0MH0QTMNxa6/3tn6bHw7y1ToEWqYXx2x5Eb66U2nOnJWLdUlP3sWa3iuuCo3rR7iiVnHh5XdQIwHXuOuoQrRqwhmFDzSIUTP3YkY/h0XszKgJeb4+P0GvEJkzUJfh5WmSRYtyzg4VxvYWxXFE3IHFsoX3zjtv0JW4o8w7ChoXLKsMDYAwICAroMk46xj7eKWWEk/rV2nSzzlFa5qZwdFckgZs0NvSij1FN9DkDK1E/5k1vttt/8v+fqfUwMrWrK0rV6tdqECnAAQG3zowCAbQtOAQAM1t0D+bT6FKPbGemxzTFMvfC6kAaYjIMKqPFUx5jbzXqKYuHt+irap0ir30nb6YRItWAnvSzGzrMFAWCspe+jOqk2uMKJFBxU6JwVyQBcQy9SaVnFibmnibkTU/9fD6Wz0E+codVlVAjClr8z8foGK3tXESoNOVPlKi4yMqOxRBRTJ4bqJKsYZm3YeHXMHYOjeRf6+JFEn1PieeZpfYDWJfYatV3NDFRxgzO4ihlps2VnHIwtWwZtJT8u03bWMnzZqEivHS9oYsdt/o7E306pQGtWWF4ZFhh7QEBAQJchvNgDAgICugyTXjzdvGUrFsTanqnlSSiSi48UVuHGX9KUKHlVS7caJkTCwx809aFQDCUfWUmjkTfx5IjDL/sKAGDzj/5M928WNxu//Fd9/Hf89/S8zOrQgq1rAQDjphapnYo2sos8tKDUs6995SQCLazIKutAufBEWXRiD9CJ7NHXZ7s2ZW60orbTcfF085atmG9y7UaqqeFdj9Jl9CjMYaV1cqEUyKSyU51U60PuFBN2ZX0UKqDQppUEskW/v/mVfp7+8s3L9KHNMxi19LOzp5ZeS5LhDbT080PPaSwkgXzcctGU/NgphAKk1ZtkaMMX/rDSSgrXVFiYhl8DpH7uTmUnBsfXnL4DKaNEGuoA3PeNrPBUqipUzgK3rya0FXTQ+MXCOW8bIVgKBAQEBByymBK5I9V8dKRQsigJ3L858sy5aFHSl2Qg95G/gL7F1ZeNFHKRkEJyuWNTLHwWmYDlXTg5S/Fhqsy1JrOo2YkJWBkp40QWefOSsXz9lmUrBwKcsdNYuDQuL3FmH1m58gpiHnkgkK0h6+vXpuibzyse1kn34a1PvwogtSPw1f50qlyxcfO6m4Ro3LXGbQrm67wPKHlHyCatyRafsTaEiZjYx1d3VY7Xl0gozb4gLDrk+QDMLoIZpQGsBqwn9T+v1q83wSgvEVJlZzDooOZpaVVMqNIecCgj3P8BMwmlX+zXXXedt0p7pFo2TsZZsjTotxXH7/4OACA+691pYzK5arhyISVrlAJpHN78mpOhF9kEzDeGXiRpBIDmbF2UYJFh44ve8QUAwJM//iwAYDZLJNpjKP8sczLE1Gn2MJykawP9pk1lpy5sTLLJxBfLk2zBUyGekMekJ8Pyy7SR0tQy+5aZPcj+i449XaWNefd/tTmCatMkEjHZoJ0TCpvaXpsMxB49WytT3BMUH2YJRHRdyHqXDL1sXV/qfTxdY6J7jZj6jY/oYh2Xr9RrYr1gcWR7bGGra55BXn+2RjJlWs8yckc/s3ZZcs30b2cL7FmxCUpmG/VHDLjJxiDvIzL08kYBqOgJxcKFBS+q7nEBoG7kmE1jlEaSxhppL5ksUQnJokxOa3ikl/TdRiKJib87aL9qq+E1AvShdIw9VGkPOJQR7v+AmYTSjD2vSruPcRYhWvU+AECy4xm7zSZkmISEKDLJDKSc4ccQx9smDL1axtCLko+AVP0SG/ULMfUzPq63P/m137Nt6ZfYGgOZWcOoUQ7E0gMVwPCANl4aN9Ritgm8UcwOAFombmcTtWzhBQ3e60RYaynGK2YJcp+i1P9O4vFFbfLaFrH96WD+lXf/j8Y9Ng3fN1CahTYNm08oVs2S7hylCZC1FuD3vLWP1VeQVCtk6GULQTCzsHhMq9aIBRJT/44ps/eR0zz22IIZUrw74UZWZGdgxp+I7Q5EsQwaCzHg0Tid7SjzjNUjKt1nZsZmLaLB1iIydh7mekfmOxllg6kblkxkO8OOPbYE1MbOpky/vnU0OUuQhoJVrqijmRxdDzkW9p3TebcqPWhV8k3YOIIqJiAgIKDLMGlVzEtbt2LAszrbiZpClo6Sxlg+/Wce4yONOjfzkp9J5cuGnemv4LFz3V/Qdsf1oaiIRhnlSDv2PVkGPBld+UQsBIrG20m/w9NQFbNx8xbMm6PXYRwVS0yxWf1nRrfNT1iwbVJi1MxOfJJotdg5JeusRSyL94+ZDkiJQyocOt6dzw3btue/xqwp0YacQhb8WErG5T3l7zJ6bblPMxvntwZfFI8nMzOe8i+4qVUHKXEcfi7Un7AK9urchWpFTtid/BPRvzR8c5RKECUT5Xg95f+aKpiABQQEBByyCC/2gICAgC7DpN0dq9s2oLZPh1DGDzvBbqcpJ01drFOjcI8DAFDasllc2Gcc2nrMVDR5/Oe26d7jz9O7mM+krJJCMOTSCAAvL9ROjSRrJEkjLZRS+AUALvzHewEAP/roGwAAA1V3oTEZ2mLbkozSwky9YlogYx91lLRjroOt9Vgi/tNJotJEUv8ns3jq+7zd9Zhs4tOBQk9zBAm0lQCf2ifKXfiz4UDPRaIQA0kU6zaF3pXYAalfOknrqiKZhz6vMEuNfREt3JpRCAEChV+AtIbqB081oUzybKdwKVvotS6GMkxD+7AF2BolGyk31GAX8nmqPlU/IssO85An4nOAXU8KrygxNoaWCOXYkJQIlSR84ZjGZ/6NhfDA8W6nUBqNl75XOVYATeUmcykRForYd8QtG5JWWDwNCAgIOCQxJTVPf3XpxwAAb/hFyqzpl8ga4pjfkFFjDFTzmF4RQ+2JxPYVb7Ft+7dt0NsWHaeP41skQeqnDgALhaEXJR9JSSOQMvUjzv0EAGDXGm0gRqZiowMpS7d1LWmcdiErizw261sgpuoydqGmA2sCyWonsoh6MJODyox3OkFV66klBTd/kzVCibnLpBiki5vKzGJtvVEpgwQyFYPoXqF9yJ+dz4gHqD/0OWOxvueMLRNT/5t7tE7/z8/RxmFJ3WN+lVNlaTett0apdUfdCgpcxht7jL6U8DqnmrJ2gbiRslaqybrPsOQeseDMF1oTMaOwMw1aBDYM22f0RQuq8plxvNspaczYDxAL94kp7PNOCUrUr53ZsEVqso1IaoUmZByBsQcEBAR0GaZE7kiyG1/t0MnEfiNP2n2eOZeUI/lscIsMvfLG8vBWzShev7gvM0Y5vk6kjPI4nbT1xaqnwm6gyIArb99O7AeKMFPlji+9+AJmz5kHIF0bAVImTZBWrj77Z1kPlCR7PHZP7I9kg9YGV85c2TNj7xGZjl6UAGXY6y0bXgEAvOf4ec7xnXMTNUOpibetlPkhO4bMeMlawFj01lvZRB85I/LVQiXJJdkwjAsbBruOJo3QGMjKgUzNeBIizTqUWGOzFgucRtO1l3JSOxj/9ShrAhYYe0BAQECXYUpi7HYlmiUFrXnDmwEA5/xM1zGlpKCaMefatSiNgfdWRFxMrHRXTVwdABoLj9FtBVO3v7Ye6mfZvDD0IpsAbl4mY3/E1Gev0nYEL93zv9NxiyKMkqCMsx9+nz0pPx6HtDT1JTpNBcrMHjqZEUxGbTOd4vydIGo1UoUHi3/Wdz4PAGjO0zFqW/PUKFxGotSeoGIlM4apG+VDWkSCxYnJUEokJmWSg9gYY2F/bdlskk3GI/ULxdSJqV/3S30+n37TsvQcG67CrSnMuhwrXsGgSbUzbvaNPDdHYs6FzpWuHZ8NWcWY6Z/eIRTfrjgKGnNtjAqvkrNGwI29ZCGTmk2iNH2xdRDL0M0LkZR7sRgjkF6HmIzJZKISj7HTuOKK1/bAh8DYAwICAroMk46xb3pxK+YOaibAf6HjYa2HJQOi8fnLAaS/Zn3r77RtGydfqAfjiYsBbqyu/sTdetsJv6s3iDihL75XperhJpa5e5w0xnrAfUwoTjp1Ur/QOY2YDt/6d7+ybe/583OcY0oDn0bBGMpcdJo1SAMxHyZiD5D3eVH/RSy/nda9aG2giLFP6xj7lhfRP3sOADc2K9ddZLy5OrLL/n+r1ut+6EnfJ9hyccQchfqGWD6f5VVEwQ6y3s2YVTkHcu9hYsI3PLjZNvnEmW4eR6YUXFGhEFLkENOupCw8z35AnivgmXUIszEOW7qOtON0qsKSuMF8A/ripjNuAq2RcPMyGi+pdkixI+8FPpYMQ6fj+MwVQ4w9ICAg4NBFeLEHBAQEdBmmVO5Y1FFiQjPNAe39zKcl8d6dAICW8WVvN00HgETsg4K2EwlB5IVXONa8oBdUVh3pulL6+s07p4ksSpaROxaFPzrBVMgdffvmXYeicQ8NDWHxdA7FsNADLbrZRJaWm6jkLIIJ6VskFmP3NtOrQQtylFpuQxE5iVBAei0j4aJqQxPMskCGHGjhsmrsDii8AAA/+K2WQl58nE7MqhvLAl/IhNL2pV2Gty6wDEGRBDPKhkrqGNf/Qw6IoNR/c86+7sXfNvRjvhoeQs1IIGWNUn6OFGaicBBZIlBeFHvVUhKTrI/qrUPLpKxDQ0M47PAlIRQTEBAQcKhh0nLHsiZQLcHU+X7EuisvPwUA2DF4NABgsJZddCDYqkt33ggAqL5Om4M1Fuh9SdIIZGWNVKOUKh/1eVy27DgLKkQRUycp5O47/4cztiL5oF1gamUXe+Ro7N+mXuzLb7jcfrawN4EPZZKN8sZW9NlkEqGmagYzrdBqpOfJq94YFmwXBYtSwYnhkvRNJDf1xelMYLSlv+8KVRsT7Jbuo5iZgJGksEKL8OYzqnzk3OOCsVfI971KVZLS4xFTv/4/NwIA/vLNWgqZLlJmk2x2Kz0+sjmgClBcGlyTNwFJJCmpi12fJoyxl7iRvDVPDaKcBVaaTVT4NbAL2XRAvW9iRCGtej/rQJ+vtQShmYDKVn5qigVsCKbuLGjzhdWSFetKv9hDlfaAQxnh/g+YSSgdY7/22mu9Vdpf2roVswdc6RVQIGPzSLhoP7LT7b31HwEAycWf9PbB+6/seE7/z84XAQCNY1YBAF4dTY8zr6Hj8U1j0kT77jLBr0Ge6ysr00hplMpK2ijef+oXfw0AeOS6t2bPUYxbxlV9jF3us22vjssdPr7Vtm0MHuH2m3Pcon59bSU6kSWW3TevH9nWfl8HMcZe5v7nVehtbFam+pvvu8nYG7E2YrrVPdqaY7x/ofM5kJXz2udgzE2T59eTZoX8mE6/HptayQwpTuwk84n1pxse1lLhPz7ZJA8yYy8rATSzBJlm79Qmjf13Cl3fqm9NgGY7sSt/5LJPkihSchiZi9n1g0r7WZWdXZFdMrtOtsKTuJmLZqO2HxGPJ5bP+42j8utMpWPsoUp7wKGMcP8HzCRMac3TIvbWCWuzGdZ7dgAAWswOtV0sNi6oeeqrLA7k1JTMGW+Zc3zS1FDlBTzaxpQL6iwWsVmZxNQJplpJU6a/vH5nqirmxS1bMXdWfhxdFoKxYCzZMj0lZouGaUfNtC19z8Q+eWKPbusWjwBYUpNgywQytNKN9P/vbup7sL/qFrBw1GxypmGY7x2btVKFF/DIU+3Y/lpZRZFXMQM32YvWMuS6hIxdOyiwVADce3NE2AETrJImyj63kXgmfe+dPANBp/YrjYeNc6hkkl5QxQQEBAR0GaZUFVNGgVE1MfHGnCXpfnJ12vwrY+J5x3X2MUzdtw/9YlrTfLLc9Bj5SNAvtJOqLeLvBGLq5/x/99htv/rkaQBSFUCGoXoUBBJelY0ZO+UJ7OlbBAAYNWLc2b/4hm3buvCPdFtz8Fjs0z+idcl8hlQmdp/3WZkYfpmZXLuxHEzw+4EzYUovT0Q6vC0swfTgUtFBoOeCqykq4rNYlpqkYhGsr6YpbkPk0qbqU1/MFpjS66lIhmSbjvWvkX2Q+oVi1Oe/Ro/lS2vScNWfnn2EOVjLGTddB6esoJ3dmGdCsNhmxNYnYte2mGwBaNw9+3amwyW1Gl07KrRDzNozU0iM7Yi0PE4oT4GvKxILN+dk1z/ofNjMK5MfQ5bHokwfkBb0USg/gw6MPSAgIKDLEF7sAQEBAV2GKfFj9yInFb85uFgf+NkH7bZxI1GciiSVUvtSFXTPgmPeGGwNV0/FeLkvgcIvALD6/zwHAPjmpSeXHm8nYYrdvTqcMrBHSyGpNkzzrX9k2ySixqvcZ2tV/72gYCxltwOdyRtnKiKk0+ca99A2025bHclsp2k/pfcDLDyXk3ziLLqRLQA5B1b8Eju+GGnz76xkz6316SQHmXOhGqXSAdEdWM1pY5OiEt2hDb8A+OfHtgMA/uDkRWb85JdurhM7930tPeCehIbthj34dW6IV1hdLKY2eR1aOhfzLyUuUkWlxFxT3ifJMZvKPU4iajkDLPFJeqbTtWMhZ0WyVwr/UD1X2oU9PFV2Tq2SsZjA2AMCAgK6DFPK2B3SkMM+iOW2DEsHgMTUIm2KqkiEiZhHFf6wSalVAWyiCS2UMpbeLuGHmBiQMvW3feV+AMDt//0NueNsNy7f5+Rzb+VTm36j/17BPLOpgo5h7rOaWpa2raZZ1OJnfgEAaKx4S9uxlFnQLkLeNZuIFcJBRXMsXYz3+ISToVRFSBkTdm+QdJH2qQg26PRL+4kUdFoslb7/gEduJ+R+Nc7GFZlouTYHGRbKxk02AbNoBixqlQIpU79prZ4dfvhUd15IUkmALfbaMel+xs3YqjGv0Wr2MVS2ajZnLB3YuVB/xPyVGHeFv7rIQoAOZxc3zXViTelY0vfdMnVPIhh9J1Rhq2bukxp7f3LDtDwJqERg7AEBAQFdhkkz9gjFLLlMog8xdaqH2pyrzYRkXVPf/rIyjVeWRb/U9Msp2E5RTdEykj1rOCQq0fu6Jab+9hu0/cBtHzA1XIX9cN6x8rbPqZvqNXVttoYVh+XszWpfmr/nmX+JqScP3WLbvnji2wEAh/UZpmlmV2M9+ruxdrRIjZHoM7oulIbuSzrqZF1lQpl0+xuMWfnYYRrfds2kfHFtGYu1SUFM7miPFrlMnY5DBlRVZgKWyDR+wcad2bVyY730jFCNUp52T0ySDL1GW3o2UaV6nmwMFFMnpv6DJ14FALz/BH33kVEZANTNwe1sx/RXRRZ03hWy7zUSTBpn5LNLoJg9naNdJND78Gtnq5fRdaF3h2HsfM0tHZTL3H0VqqS5oHU1abo1bfk5NlX2fZeHwNgDAgICugz7TxUjIJm7L4ZKTD15+WkAQOuo053POeyPbI4FQMRigukvp2lD8S3PWAjENosM+20IreUmUERCLeAbHzH1t/2ztiq+9SNneo7gP16Z5KCiz9rFzV9a8Q677YgxHRNtQiscRuqajffteAYA8HzPUbbtwj73MxXrK9yYtyz3WEXrKXn7TCeMocKsVhmkLQTFb4mpOUU5TMKMSHGn2dAYT8ghFYm5v61KwmyPKHEmSe9/Yr7WbAwu2+TxeJlIRUoRYtTcgCupmDGYuL80VOBJWLLuKjH1Gx7RSXIfOo2tBdnZjss7fXYJdhfaZvaxahP+bJprRteD4vqtqquOa7C/K2Y9qlF3iwmRSZdT8zRyk68i8x1YywgPc8/c02RRzJVQZgxJz+wQYw8ICAg4VDEljL0T86gi9k0xdWLqsVDL8P3zfrh8Bvt5MfQidmhVMAVtbD/SEoHicQWm+BRTJ6Z+wT/eaz+78/85W+9fYpx5qpIyyOt3ASve0ew1tsCGNdRN7HJTr2bqRz3/S9t27/Hn6X4NU8dWPfOCYeyTmU1MV9Ri//1Fs0T6tybiulHCWThZOItK9QaOVW7DVZzEFH8n4zChsAFYSrudLei/iW3WmOrGzjYrLgOOzMBbHlUYrReQn6As+wakMwFiyTQDIKb+vf96ybb9b6fo9aFMXN/04Vg3RCL2LZk6fzatSsfMphJXQ+4rckE5H0lTWB+T3p/F420xErIUMPvQ/cH18VTAJJPnQ+t+fCPlOTRGrM1vOwTGHhAQENBlCC/2gICAgC7D/nN3zLEUSBvn+48TKDRT2/q43bZrwYkAgF5Rp1QmEvk+k3/7pv801RkXVdR9FWTy+s87H29bAwq/AKkU8ntX6pDUXJNbPVm5X5HktF3/UiK52Nw5PJmpbv6lxdLI/Ft55XlzQJZ0MXcpAKC6/Rmzjw7tbDYz22UjG9O25rPKK+m2aQPVShNS2BQ+IZdEsZBOi6g84YdCduQ2qhI3JBN7pvs0va8I50O7EMplfpSQJ6r1KIpbeuSOvoQkQNYmNYk+Ij5q63my60E2AZR8VLcPoW5D4RcA+PajOixzxakmJCPa8pqoNlGIkpboQnuSguicKnSOtFht/pYhGSANk/EqRhzcliTh7zSkctCKsYGIuIUJfSfsuwXSOrQxq9lMz14j6cF4kvVr9yEw9oCAgIAuw6SLWUdAphoIgIwELMOOS1bbBoChhSfa/x98eT0AYHzxCqdN0eJbJ+nr9CtLZIbYgvUw50kXwkSsk0XOIvZMTP2//fMjAICf/uEZpY/TSRuJMoZkE2prFtVaz65NN75eM/bG/OW6jWE1h/UbdtK/PN3ffMY9/A808u5/FcWOSROhYdJpKrSYT3I/DxO2iTgitd0+T0y6aBcFRT/2GaS/G5zZGXMqkh+atnVaeGSs1lp+iHEnUXZhMY8XprPm9HNr6GX+tuccuQlQQMrU/+1xXQ3tspN0jYUWBCtnkOOlWUOTG3pRXVSSe5rtdJ2lxFl34M687OSkkf3OI3rn0TUzswiygeDGb1SbFkJGauWrPZOrEFb67XrddddhcHDQ/rd06dJJHTggYCYh3P8BMwmla576GMvSpUtza54SOjHn6qRNvHcnAKBpZINFcfN2JzjZdPb9lQ5P/V74v+8DANzx8Tdm+mpn2VDUtpNxdxKPL/M9UgzxxaZmLPN7NX/q/e3PdYPFr7Vtm4NaErfnodux4HcvPSg1T3Pv/80b0TOo0+Sdyjt5EjqSGqpUUlozO2bMo0z8fJQl09dj3YikdJlan+Y4o610MCQ/JJZP60U2Ts+SeOhYifmQPqMqTnwdwa6T0exBrqs1s8lBmX0InjUBun++s07H3FevMBYj3C5bsHD7fvDc3DTblslWMrnJ98xIYzNZp9aBtUfW/8oEMQCIaRZqmLu1QvFU2LL7jAxjaGgYC19zXNtnoHQopl6vo16vt28YENCFCPd/wEzClKpiyqATJunr2xr3EFMXhvX0Sz1ZpiqPJ/tyPrv7OwCAbWddDgBY1FdexVI0yyEQU//Bel2s4PeWpMyxMXiE07bMubZrU2ZGID/nbcrMHmil/3DRllQ2vnE3TnhzzigOHlRSs0qp2JM8oohRizWlepS9ypR2P2aKOlQpCYbb9kZugkxlzGV+1g6WH84cKo+pN5nJGM0jiJmTnUHdY4NLMe+KUJL5XioUx7bWu3TmdG4eZQ71T0z9jk36vr/gKJ6xZcYQu+zYTnr4ZTZrGDa2njPT4DbJdJ2plmrNXCCKc3D1ELFwWnuTz5djv2wYOX0XNE5lZhN8FmVnSz0DaI2Ve+MGVUxAQEBAl2FKbHsJ3tiU+KyIoVKxCBtHLIp5lRhP2c8mGj+nz15+g2bqhxvDrAaOyNkjy147MfQipr6luthuW5Szz2RmRkUsXx6nk5lRmTWYonFPR0SqlY7ZYW97AAAVUZjBRtYZQ7XnKuxprYKEMWoi+qSZJsVFkaFdqv92FSNFcXNS1dQNE7ZsnLFOUv6ArIJlWTpmXkaKFKdIBpAyYkcf754MxdSJqa/Zko5h1ZH6/OX6hC/vhLodNY1rRm0kC3BwOwbahzT7VPzEGogxZRzNmqS5mo2t8zwOWxIPThsbY2dKqEreg1qAwNgDAgICugzhxR4QEBDQZdh/NU9zPiuaTdTNfESJnnzhHjn9dbywc9DJQm+ZMAV9ttC4ITaMEyLJqmQCU5njFYEWSnn4xVbFEVWKyiwU512PiS48dyJtLdtH2WMfLHBZIZ9qUxUpmmJXxMKqUztXOAda33SPVUEqpTOPrrAosL7vLCxEC6BpOMF0Yf51UuhJfEByO+sjb7azMAUt6jaNE7tN2TehiISNwTobSkmnAY1Nn7fgm2Jxk8IvAHDPpmGzbcBpm9oOZBdla4lbQaopEou8FguyBiyrJCVBC8zSVoWH6qoNdzHaPiumXy6HJS2WiuJM5aU8BMYeEBAQ0GWYEsZeRrKX2UcY5gDI/TUqI7+bCKsrM+4yUkD5Ny2oJMOpx/TuXs2zaYGYapSWQdG1JKb+yohOQpkf6wXWeEwv3rVmLbRtpUd4phapxxoi7zoTysgdi9pORJ45nVCPFWQNXYCxbFoApJR0Y/TVYCdDafW2WpeZ8TUr5j7iSTCyOlfsLgASK+TSy5pliiXuOTOWfeYm6aV7wTB3XpmMjm0XB42fPC2e+pJ3miZnnmqUSiMuIDX0ovR6awFgtnMJIzH1NS9o5v6mpfpvKy3k4xU1Ze3x7OKmWExF+iyksyczJvK/9/STObaQbwJAi2Sw4pmkkfWoVNKskM4AfO9NHwJjDwgICOgyHLCapwRrruUkXZj4ladWaDuUkQ+2kzmWkeFNBHv60mj4wB4thbSGQ/XDOj5OkRSQmPqL4yZOG+vI3ELWlhiFrOOaMqBs7c6JxM07adtJnH86xtgBpElBHovYzIyPquvwjcSyjc0Crc3EvgtJzFlUQ7L/mmdob5Rmyfai6bQhc6pYJBYB6b3QY9e73NqqCXtubWKTTdZxjbLG2OuFEpysOULLjT/zakvUnzX0itxapU583pwTMfVbNrwCALj42Dn6Yzb7lBXN0ti6tCVgaw7mfBP7bjKzB9qHrQdIWaOdaZjticc2Qc6eQf21GDPnSVyBsQcEBAQcmphSxl6GLdtV90qWJUhly2RVNu32mWzsth3bHmWBVFrHjzb9Rv/PisOyO7Tpv0gxQjF1YupPvbIPALCob1amLSG97oL1MfY0EZZcZhbSbq2k3bZpA86gfGxKWF4Qa0uSrKqiUTf1Nc33YdUxLEGpIoyrVJyNCwNAjbtg2S/EjV1H5hmseGYa0gZXmfFzVQyxQtuf+VuqegBWdMIm75jELVmjFLCzDkqkkjNMXzITxbWJqa8xFVvecER6/xMzz7svK/I5AHLrGRPb56n/1DaNv7tJTePcQpjqocr7gOLxHkUdt1duh8DYAwICAroM+90ELKOi8JTMkxabnShRJCYbL5/KWPLsX3zD/n/zrX+k266QtlftUcYmgNQvFFMnpv7Ntaky54Ovc2cJdvZEv++01lEwlsl8N2UURRPp96AiilMDJ1/qvyhQYVUgTY/JE1zWTeddGU3LpNncCFLSmLi8LV2oXJYIMFYsZgKj1sIjy+/seGU5OdavZOj2OzNte/btTMfQP1+PkzTeNHMhxsuZa8u1PpBGWW55SnJgMzMKc12Iqf/wt9tty/efMM8Zp68UHuCqeSJh15tRAnH7AfOdRub7tDF7UrxweTyEGkho3fkMzK4FxBVvoRYfAmMPCAgI6DKEF3tAQEBAl2G/yx0zC39Uib2DdHuOqUpE0p1lw0Kd9JU3FtreuvCP7LaEqo6LWoZlFgvLuCTSNFX60fPwy7ApMT9Qjb37lLmmRSGqduP2Ia/NRJLeDgqiOA1b8Jq/mYthFtDIAiLKhjQkrISX3TMkUdxrUs57xP1kF/VYWMeOq+UuDtapDYsAtWxlIL2R/NMppMTrrUlbABobhZtapmYCRyYUKz3RkYYjyBEyGjPJVuad4RyW9hMyUBo/hV+AdEH1jCU6TENOk/Lp38cqVvW1tIx4NO5x2pLQw/meKaQmrQlaWauCzHUQ90fC2lK4SiW10s9AYOwBAQEBXYbSjD2vSjv/9fSxWfmZV8bT5thFaet5bLkM41OexcIy8rs8yLa8xiGxrkgudnXQb9FirWTfPhBT32tkmH3RuBmomxg2UQOuqVzY7iTR7EAg7/53FnkdyaHLxCjJxJfen4hKRvR3JBYEgVQm3EvHao4740zEIpwzrthN9KEqS1xGJ6Wv9YapClRPZYN2LNS/ObcRkxRVb7l2AQD77nLGwq9HRSQMgaoKmc9HGU0mQy9pE+B7Voip//RJncT03uM1m7d2AcZjvZfVGx036fw1iDHZsabXjmSO0pDNihLYdba2C3agLXEc9mo25xaP7UU8tg9lUJqxhyrtAYcywv0fMJMQKaVKkdLcKu1bt9pq2R0lC3nkPHmyxyLpWxnDKYkyDLiTZKa8Y8UeE7BZo5ottIz8S46pqL8iRi33tclHrLHdZhjhK+M6Ukh2BDtamnHN7Sljgtx+LEVt81B0jkNDQ1i8eHHbCu37A0X3f/+ASSzyVaw3sEkrPmMs+szUL20ZdmyThDzxcpkoIyWBvrR7yZbtWhD74qz8kBi1iNkX2cbKe46kmAAwWtPXqKbSeDHALK5ZxSArCTV/RxS7j11GDGQTszLSRd+1M3/f+ewuAMB5ywedffh7iPofM2sNNbjfH7dCoGPZtpH/ugOsxqn5uyKYO5fD8jqzZZ+B0qGYUKU94FBGuP8DZhKmhQlYM3bjTVNtE9COAXf6WVn4TMC21fS2ed49ilFm1iBtAnhCiWUihhHNb2q2RMZhRzb0GJs9ac3WMrYGeWORffgg+5uszcMBh2qliS6MvVmlESXimL8pdsqZOzH9kUo/gLSABaksuJrCqmoo2cgoO2jtJrJrWOk+GdZttqdsl7FcilXLMdA+LE48Jl4f9YZbYIYr3+rju51tdJ+2TPycM2ubZBWJ+Dmx2YLZflLA4On/6NjE1B/equ04Tlusr78zg4mIqdMMSSR78RmYaZsmUOn/Sdk+79fsIory2GQkxtKtfXFLWdvvdgiqmICAgIAuw5RYCnQSo6ZfVh43SpS/bd7xfG3KGGWVYYUTiW/n/Yb2j7xi/39rVTP1xc/8AgDQWPGW3ONNRKsPyb48Gn153q+YmDox9XuGNYNZ5YYcC8dUZmY0WU36tLYUYAyRG1mNGi00lairGGYWVd2UeiCNKdfNw9FQQg/NLhZpuyWLtew4ozbh7NDowUn1QXbB7HSkPYK99pHLPgFWjrHhGnvZvnj8WZTakzMBpzyfYcN2jUFovZ1rbjXj/jUBH0j90jCzBWLqP3v2VQDABa+ZkzY2421E7qvSFuDwlBWUCqCap0CI1aab6yLXvxrs+lC/PWoUY6wARxECYw8ICAjoMkyasefp2Mu0J1jlhijq68uIbKdjl+14m4nMCCYS67WafaZ8WWD+JaaePHQLAOClFe/Qn/e2V6IUzhqs2RQ18hd64NtI/UIxdWLqP3kqNW96z0JTpm0wjbvz45YZZ9F1z2vTidrmoCKuWAVDgxWs6DH2tK1IbyN2SOCZhVLpUpEqFq5EMSw7LWLuMkfSuHHGSsUgbPxdxKEjZ9ZhbHplyToqe9dK79M40t9SRRhTWZMt9iXatQWSpo+5rJnH2GXcXOrlHStdce0iMWPlMXDKKCWduo0emAMQU6diHQDwHqN1t2fYcmcIzsyg6V47aarnnCMZhZmD1yO3bF91dHfar1EMtSo9aFWYvXEBAmMPCAgI6DKEF3tAQEBAl2FK/djLJBIVLaTFcoro6Xci/eft00noqJMxlAk9vHji2wEAR4wZiWHvEZ7W5fsrG6LibfL+pvALAGyuahMxshLzJsAYyJqUE5GgFmE6LqI2VZrm32S3L5erASwkQGELnmRj/rUhGVmbky+6KZJPugtyMjmKL1xSlR9a0K0nboKSc7ycBCTqry9Oj0PSP3suJkwxYgaVROk3Zr3ITTiFrEUqVOeVmZlZyWLTDU9AJC7pHV1pZCLCNXyRmgy9yCaAQl42LGbaUvgFAP7dhCXfdcxc51y9Vd+ElzqFjqwHPXu/xSaEVo+azr7W7oFZOJC5WicIjD0gICCgy3DAbXt9jNK2qbjWmPLzTlBmwa0TJlnUdiK2A4f1mV9+aKYesfRraRBWRsqZdxwf2o2bL5QSU98njMPkAjfAEztEVZhJjGm6o8GyWdwKORq+imGAayZFDyFPRAGAnqY2fEqYjDBjTUCGWcIyNvFcQGKHTeUm2XCbgEgsRtoxCrsDgDFre7J63x5zU3C2nB4TTv+2zquTsGjuLcG6KSHKpuqzsVsGnVNDFEitd2sy9Z+Mt+yCaApi6mQ/cMFrtMLASlv5InhO/dl0lpu2td+1sV+uUBth6QCw92JjpDR7D4w9ICAgoMswpYy9EyZZyHybrkjfxz4mEi/PQyfx6KLEp05keck+zQBG6poB1Bkr6+Rc8mZEvr7yxlfEkoltEFMfbulbZiDWrGS0yeKo9oty+YLv+uSNpZPvYjqgjnH/DK0lYuoEWxAi++jR5UvM/6hYJPUAsGFrEeutyBi+JwGKxlRpUkKR7p+n0MdComdllR7bXmlJLK1zfZAySjtLZAk5GZkzJfoYuu/UXSWJKNzPpFUxkN6VNFsiQ6+mfA2ya0f7E1P/+fN6Zv2mpcbUjC9PiKIfFFu3cfgm75eS0pQZg/nOfesIbN1AlbjGQGDsAQEBAV2H/RZjnwiLLZOSPhnWNpmko4mgKNY+1qMZQN+OZwAAm3qPsm0WV9y27VQsk0XRbIoULsSoiKk/P6zjta8d22TbNhYcrdsKpYaSrLXNMWcUojhbRAIsXiuS7iguHXOmZ+KmsvgCxYBjxmarlPgkVDeZ+LDnmtsYdZJNcbdjoVMx/1I8m5gwj/FaRmqKP6iqiV1bU7+spW1Txc6/ZOPrY/vyulDyVXU8VW2RjUG69uDaGfPZfiYF0CpbxFoDG7dM9Sem/qtNmrm/5SjPTMaci7QxjnxWCOZ7SkQimC8pDaqVSS7LQ2DsAQEBAV2GKbEUKKMG6ST2G9k4XH5/ct8yMfcyqe4TWRuQf5dRxdDq9/M9mqkf9fwv7WdkO5DXR5n+CUVMuBPWTHFPiqkTU//lSGpNvKpNP2VyDMpo9qcVojjV77PNNm6dULk7V5PNmS/FumXJG5+NNc1+qEAF2d7S8ZRIkwfS787Gg1uk36Z9mMWvsCSQxedJWQIANdJp1/tNRyYWHgkGDMb4zd8U56f+6mMpCyeNOzH1ipj5+cpr0rUn+wSrLmEzAbkeQdp8Ys++QuCpRYeZYZgTIKZ+53PDtu35Jg5PzLxK187MKvjM1c6eaNwlCppEbBztEBh7QEBAQJchvNgDAgICugylQzFFVdonEoIpXCwzUzxbcd0jDZP9EcqEHiaS6t5Jqn6ZNrTgtLBPT9/2Hn+ebZNXgK0odFQm/NFOnukbN03NaWpPkkZaKF3F2g6b8j4DtPgrJH9F13kiVggHEnn3P6I4nfazaXIm5V9cA+4ESdN9GyKAGyJwFkJNB+MmBEPXSy5yUqgDAMbzQhp2YTCV4VWsjJIq+bjyO6fmJ41LnqtnsVCm1xNo3Nz9MpFOjXQ443jYqGUXLG0FJQovUTiXJT6RxLAhnRkjsdDazJ5jXvIRhV8A4GeUxLRMh6boO5IWG/oczaGtT7/rp8+fA5KcVuO4sOYsR2nGHqq0BxzKCPd/wExCpJRcsvEjr0r71q1bvdWyJ8LY7a9UUzAgnl6bt09Bv3nY34y9aCzEbnqM3FHFqRirMW+Z91hlZkZlFobbLVR6r4dIjZdV1oGUAZEUkhjQEQNVz8jLjYV/VrZC+/5A3v3/0uaNGBzQDJIvLNJjVTcXQUoLebUlnzQRSBfUuGSvJVixhL2OPMlGVPvJSwACgNGW3kr2A8RUqdYm2RzojnR/dkZA/VIdzypj52QhQAuWpm3sayvHJVgqZ+HWYEu8MxqeegRSfjtumTCcts79L+WFYhHckR+a/7/1WT2zOHeZZvO95sKMs6bWFkEkefnMy/gi79DQEA5bcmTbZ6B0KCZUaQ84lBHu/4CZhP1uAkboRKpGv1rETnzxonYx2qlIjCqLvASiIhmllacRU9/6dDoew9jz9i3qdyrkgt7+BWuSNrFAyjKSSJ/Txl2a4R5pGLtvJjCZ721aIIotY63xa0IZSIa5UxybqvjUeE1LWZ9TMFXO3sg0KjaU1zJVis22st+LlUiavy3rp+OyfRIzLiUYb13aHLDxUU4NNSY5ou9+swlvFGM2FsKchdtErdi95+x6EovTS+sAKRvkZmvyetbsPV1QJ5jH6JHOEGj8JGkE0pg6MfV7NuqY+0VHmfNJsjMYazucY3gGpPdVI6pkaq/mIahiAgICAroMU5qgxFEmflvUJ5CtSQikv5SUKFCG8U1Fan6ZeHwnzJ/iehRPB2PplVee1/0Ztjc+f3nH/XeCMiqhTBtfXNhsO2JA/0lMffPuLIs8YpaIgZISwrfqXzKN+qAgijPqCh/IUqDXsK8xxcypBDNNrXPNBnb+NcNwSSFCyS+2HiglO7H+7MyXlCNyvCyJh0ZFz1lT7svOkWLGZAug4KbSO/VARWzaWgx47J8rdG7GeIyUWdTEWV+IXLadKlDENeRjEPtY+97YNUADsgydVsKs+oaZ99GxE/M/xNQfeFmvV5x6mG2KWBblMNspnl7j147ZRfgK3PgQGHtAQEBAl2FKYuw+dp7H9CaiU+YxSLuyDdc0qCPLgoLj5n02EZXJhNP5KXb57Fr9t2HsE4mfF81g8vYpE8OfSNt1L6WV14+YNcdtK5UPns+mo6VA1BxLGTWbxUhNN8XYKUZaYyfTEBeKmC6xQ27VSlYExGatLl4wyjEmw6lbeunGnavCZIuPV+rCU7bLFGpmRmmLfIg+HD04nYNYA0hsLJyVCqQYvWHuVCZOKl8ApsknBZHpf5xS9vltJWLsNLuhq5uycqa/tzMBUq80nO3Ko9ijmQzF1Imp3/7MTtuWCnjI+z4RdgeAW+6vbJQhMPaAgICALkN4sQcEBAR0GaYkFFM0Re4kASU3PMGnK74KIwV9Fn3mC0m0m+oULSwW9SHbFl2PxlyT1fh6/W9s6qG2evITEvJCJB3JTD3bOkn9z/uuaaGUh1/2NvwJPLSozBcM7dS5YOwHCyqpeasWVVpuMgzdw7T4tZfFX/qhZaE2wSl2a8V6pYBmGzmhNs0KacUkFlU8IQL6PqSDIA89WCsB+2W6IQcuZLA1VMVrxC4Isu9QiYVLGYKoeKoWjRvrgCotppokJi75o3BSi8JJZl9ynuQhW1tnldpmwmWudz6QLmzb87a++llZItkE2OQj672v21D4BQBue+ZVAMD5r5njjMH6s4OBJQX6PPR9CIw9ICAgoMuw3xKU8lhxJ3JEywp5/UZZ69HKm/zH832WN9YilJlptFuk9bX1Hbu6XdsMNMyi6YtNvZh0eME+Zfotu8g72UXKMrM0YuovDOtaqiSNpHMfnX+0bZtWq9+NaQfux84X9Xwp5+zv3gqTDSqziEeMUbnCgApbhGyYBKeErgmZSFVcEynO2KwskWSVIn2ds/B04dNNzImJhfMFYmEPQBWPrNEXux6yLcH+yZOCyJvdJPyQ6Vcq10xhjyXlmaK2Kj8GjcHKKIWtQYNdEJqVUGIY1Si1MwN+LtQP+cWLmQufpRBT/4+n9YLqO4nN00yMjZsS4PKk5T4Exh4QEBDQZdjvckfZxvd5uxj7OE+vNf+SVaq3qrfneL7jtGPy7fpr189E5Y6NebqqEsm95ve6sfWi2UO74/i2lRl3mVh77nfsYSzErIip72sYid/cIwG4CSsWOWZZBxsVj72CtXuN3BgtXaMqqxg0Tmyb+mtS5R2TbMQqBlm3GooXG1ZbFbJEzqzJaleZ5ygSzwy3h5BWAjJVn5uXxWI9hMYr9wXSOP5Yy61EVISquA5SQqrPCc5nVdkvj+VnkplMvN9cX+qjHmXj/RWaEZizS4R0FIC13k2koZfn3OgciKnf9bxeR3vzUbMz+1SY9YNjHleAwNgDAgICugyTpkCdqEgAf1y3HWP0/brbVOQS4+mEScs2ZRKqpoKp8+2bDZk7rF//eveu/zmAtBZqJ7OdIuTtOyHlUsFnPpsAYlYUUyemvmaLZiRnHcFMsmSSzDRCUyFNfOHqFcN8M2n8BuOssASxbYqt03WrCKMofZDY+bdJuVHUR5SNscu0dSW2+6wQlIhR+4p+0GdkB0xjsRYDLLGKjm2N0uzAK04fABDX3DquFRsDN4lLbLYDk9Jvj6RcVs4TiOxag/gumnI2wsZt109EhICczxzbBJqlCUMvWVuVt6WYOjF1irm//bWpgoZmBVFccxQ7RQiMPSAgIKDLMKUmYBNVa7Rllx5VjP31i7Jl1/hx231WdiydxOEn0j/HspGNAIBG/3K9YfFrnX47me0U6e4nouYpu71tW8OoSP1CjJOY+majlgGYmdje6adkTyKWoi+2A+zaitlGFGfZLN3To0aBUTF6aJ5ebpk5aaaFGV5qxcuYpFDoZJg6zxmQ2/L+ZmOgY0orZ665tqUVRSEMUrVxhkl5GzB5G6RWSVlzej3oWvUoUwTFnLe1DmbnRgZmZF9AMfzExPDJGoF080BqZ5ApWefR41tFnrDepRmCYxNAbcVsipj6Dx5/2ba97KSFAPSMI+jYAwICAg5RhBd7QEBAQJehdCgmt0o7yi1YTkRiSBjzyB1pG01/O0nQwSTb5KEoJJXXvw9S7tgcPNz5vGjxtJPrXGZhe0rgkTvSdDVO3MUgCmEcyeqkkhRyYV9+7dT9jaL7nxbQRqK0dB79n11AE77jPMnGSvbMVL1uwyuehVA5sJa7qGmdMHmyFO1rQxHG3qBgsTCi+qhiXw4beqFF1IJzpLBHQ9ZHpbALhTyQb51hXSTZRaDrrGASlURiFV+UlYujlNyUqSPLQj3WTbMi+vdIb8fI1ZEcPEW4OOFji0S9WFHxicIvAHDvC/o9sGpJj5WutkNpxh6qtAccygj3f8BMQqTIVLkNcqu0b93qrZY9kQW6dvtwyF/ZiSRA+TDZpKW8MXRynMqO5wAADSMBTH77S/23kTv6jjEReeZErBQ6WaQu6suORTBOYqAv7El7W9SvmfovNmzBxae/tm2F9v2BvPv/xS1bMWDGwhOVooZmZ3IRj6r1cNRbmhUSe7Mp+kIuB6RsWCZ/Eax0z3cjUPKS+VcyViBltVKy5zM6k9LLovtLPq8y2YjsCID0OhBqI7t0Hz0DmVOyMxSfNJSPm40nNnJJ+m7Iu54WYFsV9/gAY+7SKoLPjITfu1wg5tfD1nw1f0sDMr5ISiz9x08OYe/uYXzonBPbPgOlQzGySjv9HgwPD3vbhxf7JF7sw3rq1Uj0NDXZrW/EBpv+y2N024t9mClgepr6xb53t77XSnKRKUXR/Z/qzvmLXat6WjW3MEPxi10/3AfuxW5+fKbdi90Np9RG9PfeGst+7xN7se/T/Znvhl7sY/bFns3unA4v9r27h7Fvj3G6bPMMTFjuSC/0Y449dqJdBARMCMPDwxgcHDzoYwCA448L93/AgUe7Z6B0KEai1WrhiSeewIoVK7Bp06YpnRrTNHcm9DuTxjrT+x0YGMDw8DCWLFmSKQB9oBHu/9Dv/u7T169SqtQzMGHGHscxjjjiCADA7Nmz90vMcyb1O5PGOpP7PdhMnRDu/9DvgepT9lvmGQg69oCAgIAuQ3ixBwQEBHQZJvVir9fruOaaaxy1wFRgJvU7k8Ya+p1azLRzDv3uv36n21gnvHgaEBAQEDA9EUIxAQEBAV2G8GIPCAgI6DKEF3tAQEBAlyG82AMCAgK6DOHFHhAQENBlCC/2gICAgC5DeLEHBAQEdBnCiz0gICCgy/D/A2RbA+9UQ3prAAAAAElFTkSuQmCC", "text/plain": [ "
" ] @@ -259,7 +366,6 @@ } ], "source": [ - "%matplotlib inline\n", "fig, axs = plt.subplots(1, 2, figsize=(4, 2))\n", "\n", "for ax, l, h in zip(axs, ['Exact', 'Approximate'], [exact_hess, approx_hessian]):\n", @@ -283,7 +389,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 128, "id": "abbbbfd6-7d17-4b93-880a-3352903b56c4", "metadata": { "tags": [] @@ -295,17 +401,17 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 129, "id": "fdd80af3-8c18-40d8-b971-4a473bc91498", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "5.451357701087528" + "7.268532902551082" ] }, - "execution_count": 13, + "execution_count": 129, "metadata": {}, "output_type": "execute_result" } @@ -316,7 +422,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 130, "id": "6b1af348-4bc9-4ced-9a12-44b3e49abe9c", "metadata": { "tags": [] @@ -328,7 +434,7 @@ "5.5067174465850215" ] }, - "execution_count": 14, + "execution_count": 130, "metadata": {}, "output_type": "execute_result" } @@ -356,7 +462,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 131, "id": "bce41a81-6c88-4b0c-9d8d-0891d1832fd6", "metadata": { "tags": [] @@ -366,7 +472,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Plotting at 16 steps: 5, 21, 37, 54, 70, ...\n" + "Plotting at 16 steps: 5, 54, 104, 154, 203, ...\n" ] } ], @@ -377,7 +483,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "id": "fe39ce86-1806-4367-8c86-e3ef58f81f84", "metadata": { "tags": [] @@ -387,19 +493,16 @@ "name": "stderr", "output_type": "stream", "text": [ - " 6%|███████████████▍ | 1/16 [00:00<00:06, 2.21it/s]/home/lward/miniconda3/envs/jitterbug/lib/python3.10/site-packages/sklearn/linear_model/_coordinate_descent.py:628: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations, check the scale of the features or consider increasing regularisation. Duality gap: 1.968e-09, tolerance: 5.772e-10\n", - " model = cd_fast.enet_coordinate_descent(\n", - " 25%|█████████████████████████████████████████████████████████████▌ | 4/16 [00:02<00:09, 1.32it/s]/home/lward/miniconda3/envs/jitterbug/lib/python3.10/site-packages/sklearn/linear_model/_coordinate_descent.py:628: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations, check the scale of the features or consider increasing regularisation. Duality gap: 1.086e-08, tolerance: 1.906e-09\n", - " model = cd_fast.enet_coordinate_descent(\n", - "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [00:27<00:00, 1.72s/it]\n" + " 62%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 10/16 [00:00<00:00, 16.03it/s]" ] } ], "source": [ "zpes = []\n", "for count in tqdm(steps):\n", - " model = LinearHessianModel(reference=data[0], regressor=LassoCV)\n", - " hess_model = model.train(data[:count])\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " hess_model = model.train(data[:count])\n", " \n", " approx_hessian = model.mean_hessian(hess_model)\n", " approx_vibs = VibrationsData.from_2d(data[0], approx_hessian)\n", @@ -416,23 +519,12 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "id": "1c6706a9-a27f-448f-81d4-957939bb2ca8", "metadata": { "tags": [] }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVQAAAC+CAYAAABqOvflAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlE0lEQVR4nO3de1yTdf8/8NeAbWxzDGQgTI4eEBVBETTNQjJQyfPXNLOk+/bbN1NMb+3OvLvLQxlqD+3w9c7MCi3rh96Ffk3LwsQTZhroDYkiIkcBkeNgY4Ntn98fyHRxcOBgbHs/H4893K7rc23vz654de06fC4OY4yBEELIQ7MzdwGEEGItKFAJIcREKFAJIcREKFAJIcREKFAJIcREKFAJIcREKFAJIcREKFAJIcREHMxdwMPQ6XQoKSmBWCwGh8MxdzmEEAvBGENdXR1kMhns7Ey3XWnRgVpSUgJvb29zl0EIsVBFRUXw8vIy2ftZdKCKxWIAzV+Kk5OTmashhFgKuVwOb29vfYaYikUHasvPfCcnJwpUQkinmXpXIR2UIoQQE7HoLVRCSO+nUGtQVK1EjbIJw2ROcHLkmrukbkOBSgh5KBqtDqW1KhRWKVFUpWz+t7oBhVVKFFcpUalo1LflcIChHk4Y498XY/z7ItyvL9zEfDNWb1oUqISQDjHGUKVo1AdlkUFwKlFSo4JW1/Gwys5CLkQ8B9yqaUBWqRxZpXLsOZcPABggFRkErJeLwGJPg7SKQFUoFLC3t2813d7eHo6Ojgbt2mNnZweBQNCltkqlEu2N083hcCAUCrvUtqGhATqdrt06RCJRl9qqVCpotVqTtBUKhfr/+NVqNTQajUnaCgQC/fmBjY2NaGpqMklbR0dH/X8rnWnb1NSExsbGdtvy+Xw4ODh0uq1Go4FarW63LY/HA5fL7XRbrVYLlUrVblsulwsej9dmW52O4UppLVJzKpB6oxJXyxVo0DZ/v4zpwJpa943nYAcvZ0d4uznB310Cbxch+js7wk0IeLkI9T/zy+UqpBVU4/eCKqQXyZFTqcbNCgVy79Tjm9QcAICnhI/Rvn0x2s8F4b4uGODWB1wuF3w+/24NDEqlst2+GfN339Hf90NhFqy2tpYBaPcRExNj0F4oFLbbNiIiwqCtVCptt21YWJhBW19f33bbDhs2zKDtsGHD2m3r6+tr0DYsLKzdtlKp1KBtREREu22FQqFB25iYmA6/t/vNnTu3w7b19fX6trGxsR22LS8v17ddunRph23z8vL0bV999dUO2/7xxx/6tuvWreuw7YULF/Rtt27d2mHblJQUfdsdO3Z02PbIkSP6tgkJCR22PXDggL7tgQMHOmybkJCgb3vkyJEO2+7YsUPfNiUlpcO2W7du1be9cOFCh20ljy5gfq8fYY+8e5xFv7mvw7avvvqq/n3z8vI6bLt06VJWo2hkx7PK2D++Odth24AJT7Hdp3PZpcJqlpFX1mHbydNmsWulcnatVM6ulnacEbW1tcyUrGILlRDSORodw/mblTh1/Q4OJ1/qsO2icb547+0p4DvY48qVKwh623R1SIRcTBraD8FSO7zbQbtbNSq8c/QqAEDX2P6WNwCcybmDyR+cNl2RncBhzHLvKSWXyyGRSFBSUtLmeaj0k7/ttvST3zZ/8hdVKnE2twJnb1TgQn4NlC0/43VaME0ThsnEmDBIigmDpBjp4wKuffP8+3cP6HQ6NDQ0tFtDZ9o6ODg88Gd8o0aHKyW1SC+S43KJAhnFNWjU6MCa1MDd3awte1tbdrvacexhx+XpX2vvBjDnvvZatRL/2fIMamtrTXoOu1UEqqm/FEKsgbJRg/M3K3H6egVOXb+DvArDjQRXEQ+PDZYiYogbJgxys6qj7Q/SXdlBP/kJsRJNWh0uF9Xg3I1KpOZW4FJhNZq097aXHOw4CPVxQcQQNzw+2A3DZU6ws7PMo+m9FQUqIRZKp2O4WibXB+iFvCooGw13z/R3FiBiiBsiAtwwfqArxFZ8Un1vQIFKiIVgjKGgUonU3Aqcu1GJX29WokphuK/WRcjF+IFSjB/kikcHSuHrKrTYczotEQUqIb1YuVyFc7mVSL1RgXO5lbhVY3iQR8izxxj/vnj0bogO9aCf8eZEgUpIL6Fq0qKoSokb5fX4La8KqTcqkFNeb9CGa8/BKB8XjB/oikcHSRHi5QyeA41x1FtQoBLSg1RNWhRXK5FXoUR+hQL5lXcfFUqU1Dbgz+fccDjAcJnT3S1QKcL9XCDk0Z9tb0Vrhli9Ro0OFfVq3Km7+7j7/P5p9WoN+vAdIOI7oI+jA8Qtz1sejobPRTwHiO9OE/EdDLYS1ZrmLc28CiUKKhXIq+g4NO/Xh+8AP6kQI72d8ehAKR4Z4AoXEa8HviViChSoxCIxxlDb0IQyuQrl8ntBWXFfYLZMq1G2f/K+qfAc7NCH7wCuPQfldeoOQ1PEs4efVAQ/qQj+riL4ugrhf/e1q4hHB5EsGAUq6XWatDqU16lRVqvCbbnq3r93n5fJm1+rmtq/MuzPHOw4cBPzmx99+Prn0rvPRXwHKNUa1N99KNQa1N39t151b3rzPC3qVBrUq5v0NTRqdKjS3Dvifn9o+rkK4efa8lwEaR8KTWtFgUp6XL1ag8ziWhRXK+8LSjXK5A0oq1WjUtHxFt79nIVc9BM7wt2pOSilfwrMlgCVCLjdcvRbo9VBodaivrE5eFVNWsicBRSaNooClXQrxhhuViiQXlCNS0U1SC+oxvXbdXjA8Jng2nPgLnaEh8QRHk6O6OfkCA8Jv/lfp+bp/Zwc4chtPWxjT3Kwt4NEaAeJkE6YJxSoxMTqVE24XFSDS4U1SC+sxqXCGtQ2tN6H2d9ZgIHufeDhxG8OTIPgdERfIY/OpyQWhwKVdJlOx5B7p94gPK+X17X6uc53sEOwlwShPi4Y5eOMUT4u6Ofk2PabEmLBKFBJp9yWq7D/YhF+L6jG5cJqyFWth+Dz7itoDk9vZ4T6uiDQw4lOPic2gQKVGOVOnRqfnMrFvvMFUGvuHV0XcO0R7CXBKB8XhPo4Y6SPM9zFtPVJbBMFKulQlaIRu07n4stzBWhoah7JaLSvC2aNlGGUjwsCPcRwsKetT0IAClTSjlplE3afuYmE1Dwo7g4JF+IlwaroIXh8sJROCSKkDRSoxIBc1YQvzubh8zN5qFM37x8dLnPCqqgAPBHoTkFKSAcoUAkAQKHWYM+5fHx6+qb+NKch/cT4W1QAJg/vR0FKiBEoUG1cQ6MWX/6aj12nb+oHKx7oJsLfogIQE+RJ54IS0gkUqDZK1aTFN78V4uOTuaiob76Tpp+rECueHIwZIf1hT0FKSKdRoNoYtUaLAxeLsCPlBm7Lm4PUy0WAVyYNxpxR/emIPSEPgQLVRjRpdfg2rRg7TtzQ30ZDJnFE3BODMXe0F514T4gJdCpQa2trcfDgQZw5cwb5+flQKpVwc3PDqFGjMHnyZIwfP7676iQPobS2AS/vS8flohoAgLuYj7gnBmF+uDf4DuYdXIQQa2LUZklpaSlefPFFeHp6YuPGjVAoFBg5ciQmTZoELy8vpKSkICoqCsOGDcP+/fu7VEh8fDw4HA5WrlzZpeVJ287lVmDaR2dxuagGTo4OeHPaMJx+LRKLxvlRmBJiYkZtoYaEhGDRokW4cOECgoKC2mzT0NCAQ4cOYfv27SgqKsKrr75qdBEXL17Ep59+iuDgYKOXIR1jjOGzM3nYfOwatDqGoZ5O2PXcaPi4Cs1dGiFWy6hAvXLlCtzc3DpsIxAIsGDBAixYsAB37twxuoD6+nosXLgQu3fvxjvvvGP0cqR9CrUGr32XgaMZpQCAOaP6Y9PsERDwaIuUkO5k1E9+Nzc3HDp0CFqt1qg3fVD43m/ZsmV46qmn8OSTTxq9DGnfzTv1mP1xKo5mlMLBjoMNM4Zj27wQClNCeoDRB6Xmzp0LqVSK2NhY/OUvf0FgYOBDf3hiYiLS09Nx8eJFo9qr1Wqo1Wr9a7lc/tA1WJOfr5Rh9YH/oE6tgbuYj48XhiLMr6+5yyLEZhh9rkxhYSGWL1+OgwcPYvjw4ZgwYQISEhKgUCi69MFFRUVYsWIF9u3bB0dH44Z7i4+Ph0Qi0T+8vb279NnWRqtj2PZzNv7nqzTUqTUI93PBkeUTKEwJ6WEcxoy9Hdo9p06dwhdffIGkpCRwOBzMmzcPixcvxrhx44x+j0OHDmH27Nmwt7/3U1Sr1YLD4cDOzg5qtdpgHtD2Fqq3tzdqa2vh5OTU2W5YhRplI15JvIzT15v3W78w3g9vPDUUXDpBn5B2yeVySCQSk2dHlwK1RX19PRITE5GQkIDz588jMDAQV65cMWrZuro6FBQUGExr2ZWwZs2ads8muF93fSmW4o9btViyLw3F1Q1w5Nohfs4IzB7lZe6yCOn1uis7HupKqT59+iAyMhL5+fm4du0arl+/bvSyYrG4VWiKRCK4uroaFaa27ru0YvzjYCbUGh18+grxyXOjMUxme/9TIaQ36dLvQqVSib179yIiIgIBAQHYv38/Vq1ahfz8fBOXR/6sUaPDm4f+wOp//wdqjQ4Th7jh+7gJFKaE9AKd2kJNTU3FF198gX//+9/QaDSYM2cOjh8/jsjISJMUc/LkSZO8j7W6LVfh5X1pSC+sAQC8MmkwVk4aTEPsEdJLGB2oAQEByM3NxahRo7BlyxY8++yzkEgk3Vkbuc+FvCos/TodFfVqiB0d8MH8kZg0tJ+5yyKE3MfoQJ0yZQoWL16MkJCQ7qyHtGHvuXy8fSQLGh3DkH5i7Hp+NPykInOXRQj5E6MD9aOPPtI/12g0OHnyJHJzc/Hss89CLBajpKQETk5O6NOnT7cUaou0Ooa3j2Rhz7l8AMCMEBk2/9cICHk06iIhvVGn/zILCgowZcoUFBYWQq1WIyoqCmKxGFu3boVKpcInn3zSHXXaHGWjBisSLyM56zYA4LUpQ/ByxEC6txMhvVinj/KvWLECYWFhqK6uhkAg0E+fPXs2fvnlF5MWZ6vu1Kmx4NPzSM66DZ6DHf53wSgsnTiIwpSQXq7TW6hnz55FamoqeDyewXRfX1/cunXLZIXZqhvldXgh4SKKqxvgIuRi96IwuoSUEAvR6UDV6XRtjjpVXFwMsVhskqJs1bncCiz5Kg1ylQZ+rkIk/GUM/OngEyEWo9M/+aOiovDBBx/oX3M4HNTX12PdunWIiYkxZW025bu0YsR+cQFylQajfV2QtPRRClNCLEynr+UvKSlBZGQk7O3tkZOTg7CwMOTk5EAqleL06dNwd3fvrlpbsYZr+Rlj+OiXG3j/ePNlu0+N8MS2eSFw5NL4pYR0l15zLb9MJsPly5eRmJiItLQ06HQ6LF68GAsXLjQ4SEUerFGjw9qkTHyXXgwAeCliANZMDqQrnwixUA812pS5WfIWam1DE17el4ZzuZWw4wAbZwbhuUd8zV0WITahu7LDqH2ov/76q9FvqFAojB7Cz1YVVysxd+c5nMuthJBnj89jwylMCbECRgXqokWLEBUVhQMHDqC+vr7NNllZWfjHP/6BQYMGIT093aRFWpOM4hrM/vgccsrr0c+JjwMvjUNkYM/tdyaEdB+j9qFmZWVh165deOutt7Bw4UIEBARAJpPB0dER1dXVuHbtGhQKBebMmYPk5GQaz7Qdx7NuY/n/u4SGJi0CPcRI+Es4PCW035kQa9Hpfajp6ek4c+YM8vPz0dDQAKlUilGjRiEyMhJ9+/bsCeiWtA/1y1/zsf7wFegY8NhgKT5eGAqxI9fcZRFik3rNUf7Q0FCEhoaarABrp9MxvPvDVXx2Ng8A8Ey4N96eFUT3fCLECtGwRd2ooVGLv+2/jGNXygAAf588BEsn0gAnhFgrCtRuomrSYuFn55FeWAOevR3eezoYM0f2N3dZhJBuRIHaTT4/m4f0whpIBM0DnIzxpwFOCLF2tCOvG5TUNGDHiRsAgI0zh1OYEmIjKFC7wbs/XEVDkxbhfi6YESIzdzmEkB5idKDGxMSgtrZW/3rTpk2oqanRv66srMSwYcNMWpwlOn+zEkcySmHHAdbPGE4HoAixIUYH6k8//QS1Wq1/vWXLFlRVVelfazQaZGdnm7Y6C6PR6rD+cPNlt8+O9cFwGd0VlhBbYnSg/vn8fwseU6XbfP1bIa6V1cFZyMXqqCHmLocQ0sNoH6qJVCkase3n5i301dFD4CLiPWAJQoi1MTpQORxOq/2BtH/wnvd+yoZcpcFQTyc8O8bH3OUQQszA6PNQGWN44YUXwOfzAQAqlQpLliyBSNR8m47796/amsziWiReLAQAbJgxHPY0QDQhNsnoQI2NjTV4/dxzz7Vqs2jRooevyMIwxrDu8B9gDJg5UkbnnBJiw4wO1ISEhO6sw2IdvHQL6YU1EPLssXbqUHOXQwgxo05delpQUICff/4ZTU1NmDhxos2fd1qnakL8j9cAAHFPDIKHxNHMFRFCzMnoQD19+jRiYmKgVCqbF3RwwN69e7FgwYJuK663+98TN3CnTg1/qQiLJ/ibuxxCiJkZfZT/zTffRGRkJIqLi1FZWYm//vWveO2117qztl7tRnk9vrg7xulb04aB70C3fSbE1hkdqJmZmYiPj4dMJoOLiwu2bduGkpISVFdXd/nD4+PjER4eDrFYDHd3d8yaNcsirrZijGHjkSxodAyTAt3pnlCEEACdCNSamhq4u98LDpFIBKFQaHA9f2edOnUKy5Ytw/nz55GcnAyNRoPo6GgoFIouv2dPSM66jdPX74Bnb4c3p9n2fmRCyD2dOiiVlZWFsrIy/WvGGK5evYq6ujr9tODgYKPf79ixYwavExIS4O7ujrS0NDz++OOdKa3HqJq0ePtoFgDgvx/zh59UZOaKCCG9RacCddKkSa2u4Z82bRo4HA4YY+BwONBqtV0upmU0q56+2V9n7D59E0VVDfBwcsSyyEHmLocQ0osYHah5eXndWQcYY1i1ahUmTJjQ7m2o1Wq1wRVZcrm8W2v6s1s1DfjXyeaBo9fGBELEpxseEELuMToRfH19u7MOxMXFISMjA2fPnm23TXx8PDZs2NCtdXTk3aNXoWrSYYx/Xxo4mhDSitEHpSZNmoSkpKR251dUVGDAgAFdKmL58uU4fPgwUlJS4OXl1W67tWvXora2Vv8oKirq0ud1xbncChzNvDtw9HQaOJoQ0prRgZqSkoJ58+Zh3bp1bc7XarUoKCjo1IczxhAXF4ekpCScOHEC/v4dnxzP5/Ph5ORk8OgJGq0OGw43H4haONYXw2Q987mEEMvSqfFQd+7ciQ8//BCzZ89GfX39Q3/4smXLsG/fPnzzzTcQi8UoKytDWVkZGhoaHvq9Temr8wXIvl0HFyEXq6MDzF0OIaSX6lSgzpw5E7/++iuysrIwbtw43Lx586E+fOfOnaitrcXEiRPh6empf+zfv/+h3teUKuvV2J58HQDw6uQhcBbSwNGEkLZ1esT+oUOH4sKFC/D29kZ4eDiOHz/e5Q9njLX5eOGFF7r8nqb23k/ZqFNpMFzmhGfCaeBoQkj7unQLFIlEgqNHj+LFF19ETEwM3n//fVPX1StkFNdg/+/NB75o4GhCyIMYfdpUW7c/2bx5M0aNGoXFixfjxIkTJi/OnHQ6hrf+7woYA2aP6o8wv957sQEhpHfo8l1PW8yfPx9nz55FZmamyYrqDZIu3cLlohqIePZ4fWqgucshhFgAo7dQU1JS2r0kdOTIkUhLS8PRo0dNVpg51amasPnuwNHLJw1GPycaOJoQ8mBGB2pERESH811dXa3mnlJbjl1DRX3zwNF/fZQGjiaEGKdLB6Ws2W83K7HvfPMdTDfNCgLPgb4iQohxKC3uo2rS4vWk5n3Bz4R7Y/wgqZkrIoRYEgrU+3z4Sw7yKhRwF/OxNobuYEoI6RwK1Lv+uFWLT083X/n1zqwgSARcM1dECLE0FKgAmrQ6vPZtBrQ6hqeCPRE93MPcJRFCLBAFKoBPT99EVqkcEgEX66cPN3c5hBALZfOBmnunHh/+kgOg+XbQbmK+mSsihFgqmw5UnY7h9e8y0KjR4fEAN8wJ7W/ukgghFsymA/Xr3wpwMb8aQp493p0dRKPwE0Ieis0G6q2aBv3lpWumBMLLRWjmigghls4mA5UxhjcOZkLRqMVoXxc8/0j33oCQEGIbbDJQ/+9yCU5m3wHP3g5b/msE7GicU0KICdhcoFbUq7Hh+ysAgFcmDcIgd7GZKyKEWAubC9QN32ehWtmEQA8xXooYaO5yCCFWxKYC9XjWbXz/nxLYcYCtc4PBtbep7hNCupnNJIpc1YR/HvoDAPDiYwMQ7OVs3oIIIVbHZgJ184/XUCZXwc9ViJVPBpi7HEKIFbKJQP01txLf/NY8aHT8nGAIePZmrogQYo2sPlBVTVqsTcoAADw71gfjBrqauSJCiLWy+kB9//h15Fcq0c+JT3cvJYR0K6sO1IziGuy+O2j0plkj4ORIg0YTQrqP1QZqy6DROgZMD5HhyWH9zF0SIcTKWW2g7jqVi2tldXARcrFu+jBzl0MIsQFWGag3yuvx0S83AABvTR8GaR8aNJoQ0v2sLlB1OoY132WgUavDxCFumDWSBo0mhPQMqwvUr84XIK2gGiKePTbNHkGDRhNCeozVBWqVohF2HOD1qYHo7ywwdzmEEBti9kD9+OOP4e/vD0dHR4wePRpnzpx5qPf7W1QAflzxOBaOpUGjCSE9y6yBun//fqxcuRJvvPEGLl26hMceewxTp05FYWHhQ73vEA8xDRpNCOlxHMYYM9eHjx07FqGhodi5c6d+2tChQzFr1izEx8c/cHm5XA6JRILa2lo4OTl1Z6mEECvSXdlhti3UxsZGpKWlITo62mB6dHQ0zp071+YyarUacrnc4EEIIb2Fg7k+uKKiAlqtFv36GV7B1K9fP5SVlbW5THx8PDZs2NBqOgUrIaQzWjLD1D/QzRaoLf58WhNjrN1TndauXYtVq1bpX+fl5WHkyJHw9vbu1hoJIdaprq4OEonEZO9ntkCVSqWwt7dvtTVaXl7eaqu1BZ/PB59/76onX9/mI/mFhYUm/VLMTS6Xw9vbG0VFRVa1b9ga+2WNfQKsv1+FhYXgcDiQyWQmfX+zBSqPx8Po0aORnJyM2bNn66cnJydj5syZRr2HnV3zLmCJRGJVK72Fk5MT9ctCWGOfAOvtV3dlhll/8q9atQrPP/88wsLCMG7cOHz66acoLCzEkiVLzFkWIYR0iVkDdf78+aisrMTGjRtRWlqKoKAg/PDDD/qf8oQQYknMflBq6dKlWLp0aZeW5fP5WLduncF+VWtA/bIc1tgngPrVVWY9sZ8QQqyJ2a/lJ4QQa0GBSgghJkKBSgghJmLRgWrqof960vr168HhcAweHh4e+vmMMaxfvx4ymQwCgQATJ07ElStXzFhx206fPo3p06dDJpOBw+Hg0KFDBvON6Ydarcby5cshlUohEokwY8YMFBcX92AvWntQv1544YVW6++RRx4xaNPb+hUfH4/w8HCIxWK4u7tj1qxZyM7ONmhjaevLmD715Lqy2EDtrqH/etLw4cNRWlqqf2RmZurnbd26Fdu3b8eOHTtw8eJFeHh4ICoqCnV1dWasuDWFQoGQkBDs2LGjzfnG9GPlypU4ePAgEhMTcfbsWdTX12PatGnQarU91Y1WHtQvAJgyZYrB+vvhhx8M5ve2fp06dQrLli3D+fPnkZycDI1Gg+joaCgUCn0bS1tfxvQJ6MF1xSzUmDFj2JIlSwymBQYGstdff91MFXXOunXrWEhISJvzdDod8/DwYJs3b9ZPU6lUTCKRsE8++aSHKuw8AOzgwYP618b0o6amhnG5XJaYmKhvc+vWLWZnZ8eOHTvWY7V35M/9Yoyx2NhYNnPmzHaXsYR+lZeXMwDs1KlTjDHrWF9/7hNjPbuuLHILtStD//VGOTk5kMlk8Pf3xzPPPIObN28CaB70payszKB/fD4fERERFtU/Y/qRlpaGpqYmgzYymQxBQUG9vq8nT56Eu7s7AgIC8OKLL6K8vFw/zxL6VVtbCwDo27cvAOtYX3/uU4ueWlcWGahdGfqvtxk7diy+/PJL/PTTT9i9ezfKysowfvx4VFZW6vtgyf0DYFQ/ysrKwOPx4OLi0m6b3mjq1Kn4+uuvceLECWzbtg0XL17EE088AbVaDaD394sxhlWrVmHChAkICgoCYPnrq60+AT27rsx+pdTD6MzQf73N1KlT9c9HjBiBcePGYeDAgdi7d69+h7kl9+9+XelHb+/r/Pnz9c+DgoIQFhYGX19fHD16FHPmzGl3ud7Sr7i4OGRkZODs2bOt5lnq+mqvTz25rixyC7UrQ//1diKRCCNGjEBOTo7+aL+l98+Yfnh4eKCxsRHV1dXttrEEnp6e8PX1RU5ODoDe3a/ly5fj8OHDSElJgZeXl366Ja+v9vrUlu5cVxYZqPcP/Xe/5ORkjB8/3kxVPRy1Wo2rV6/C09MT/v7+8PDwMOhfY2MjTp06ZVH9M6Yfo0ePBpfLNWhTWlqKP/74w6L6WllZiaKiInh6egLonf1ijCEuLg5JSUk4ceIE/P39DeZb4vp6UJ/a0q3rqlOHsHqRxMRExuVy2eeff86ysrLYypUrmUgkYvn5+eYuzSirV69mJ0+eZDdv3mTnz59n06ZNY2KxWF//5s2bmUQiYUlJSSwzM5MtWLCAeXp6MrlcbubKDdXV1bFLly6xS5cuMQBs+/bt7NKlS6ygoIAxZlw/lixZwry8vNjx48dZeno6e+KJJ1hISAjTaDTm6laH/aqrq2OrV69m586dY3l5eSwlJYWNGzeO9e/fv1f36+WXX2YSiYSdPHmSlZaW6h9KpVLfxtLW14P61NPrymIDlTHG/vWvfzFfX1/G4/FYaGiowakSvd38+fOZp6cn43K5TCaTsTlz5rArV67o5+t0OrZu3Trm4eHB+Hw+e/zxx1lmZqYZK25bSkoKA9DqERsbyxgzrh8NDQ0sLi6O9e3blwkEAjZt2jRWWFhoht7c01G/lEoli46OZm5ubozL5TIfHx8WGxvbqube1q+2+gOAJSQk6NtY2vp6UJ96el3RaFOEEGIiFrkPlRBCeiMKVEIIMREKVEIIMREKVEIIMREKVEIIMREKVEIIMREKVEIIMREKVEIIMREKVEIA7NmzB87OzuYug1g4ClTSY9q6tw+Hw8GUKVPMXRrmz5+P69evm7sMYuEsejxUYnmmTJmChIQEg2l8Pr/bPq+xsRE8Hu+B7QQCAQQCQbfVQWwDbaGSHsXn8+Hh4WHwaBkpncPh4LPPPsPs2bMhFAoxePBgHD582GD5rKwsxMTEoE+fPujXrx+ef/55VFRU6OdPnDgRcXFxWLVqFaRSKaKiogAAhw8fxuDBgyEQCBAZGYm9e/eCw+GgpqYGQNs/+b///nuMHj0ajo6OGDBgADZs2ACNRqOfv379evj4+IDP50Mmk+GVV17phm+MWBIKVNKrbNiwAfPmzUNGRgZiYmKwcOFCVFVVAWgeozIiIgIjR47E77//jmPHjuH27duYN2+ewXvs3bsXDg4OSE1Nxa5du5Cfn4+5c+di1qxZuHz5Ml566SW88cYbHdbx008/4bnnnsMrr7yCrKws7Nq1C3v27MGmTZsAAN9++y3ef/997Nq1Czk5OTh06BBGjBjRPV8KsRwPP4AWIcaJjY1l9vb2TCQSGTw2btzIGGseiu2f//ynvn19fT3jcDjsxx9/ZIwx9uabb7Lo6GiD9ywqKmIAWHZ2NmOMsYiICDZy5EiDNmvWrGFBQUEG09544w0GgFVXVzPGGEtISGASiUQ//7HHHmPvvvuuwTJfffUV8/T0ZIwxtm3bNhYQEMAaGxu7+G0Qa0T7UEmPioyMxM6dOw2m3X+HyuDgYP1zkUgEsVisv0NlWloaUlJS0KdPn1bvm5ubi4CAAABAWFiYwbzs7GyEh4cbTBszZkyHdaalpeHixYv6LVIA0Gq1UKlUUCqVePrpp/HBBx9gwIABmDJlCmJiYjB9+nQ4ONCflC2jtU96lEgkwqBBg9qdz+VyDV5zOBzodDoAgE6nw/Tp07Fly5ZWy7XczqLlM+7H2rjZGnvAMMA6nQ4bNmxo8yZujo6O8Pb2RnZ2NpKTk3H8+HEsXboU7733Hk6dOtWqD8R2UKASixEaGorvvvsOfn5+ndoSDAwMxA8//GAw7ffff3/gZ2VnZ3cY/gKBADNmzMCMGTOwbNkyBAYGIjMzE6GhoUbXRqwLBSrpUWq1utVdNR0cHCCVSh+47LJly7B7924sWLAAf//73yGVSnHjxg0kJiZi9+7dsLe3b3O5l156Cdu3b8eaNWuwePFiXL58GXv27AHQ+pbJLd566y1MmzYN3t7eePrpp2FnZ4eMjAxkZmbinXfewZ49e6DVajF27FgIhUJ89dVXEAgE8PX17dwXQqwKHeUnPerYsWPw9PQ0eEyYMMGoZWUyGVJTU6HVajF58mQEBQVhxYoVkEgksLNr/z9lf39/fPvtt0hKSkJwcDB27typP8rf3jmwkydPxpEjR5CcnIzw8HA88sgj2L59uz4wnZ2dsXv3bjz66KMIDg7GL7/8gu+//x6urq6d/EaINaF7ShGbtGnTJnzyyScoKioydynEitBPfmITPv74Y4SHh8PV1RWpqal47733EBcXZ+6yiJWhQCU2IScnB++88w6qqqrg4+OD1atXY+3ateYui1gZ+slPCCEmQgelCCHERChQCSHERChQCSHERChQCSHERChQCSHERChQCSHERChQCSHERChQCSHERChQCSHERP4/WxR5VrPlxN8AAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(3.5, 2))\n", "\n", diff --git a/tests/models/test_linear.py b/tests/models/test_linear.py index de2e89a..22a96bf 100644 --- a/tests/models/test_linear.py +++ b/tests/models/test_linear.py @@ -2,7 +2,7 @@ from ase.build import molecule from ase.vibrations import VibrationsData -from jitterbug.model.linear import get_displacement_matrix, LinearHessianModel +from jitterbug.model.linear import get_model_inputs, HarmonicModel def test_disp_matrix(): @@ -11,18 +11,20 @@ def test_disp_matrix(): # With a single displacement only the first term should be nonzero atoms.positions[0, 0] += 0.1 - disp_matrix = get_displacement_matrix(atoms, reference) - assert disp_matrix.size == 21 - assert (disp_matrix != 0).sum() == 1 - assert np.isclose(disp_matrix[0], 0.01) + disp_matrix = get_model_inputs(atoms, reference) + assert disp_matrix.shape == (27,) # 6 linear, 21 harmonic terms + assert (disp_matrix != 0).sum() == 2 # One linear and one harmonic term + assert np.isclose(disp_matrix[0], 0.1) # Linear terms + assert np.isclose(disp_matrix[6], 0.01 / 2) # Harmonic terms # With two displacements, there should be 3 nonzero terms atoms.positions[1, 0] += 0.05 - disp_matrix = get_displacement_matrix(atoms, reference) - assert (disp_matrix != 0).sum() == 3 - assert np.isclose(disp_matrix[0], 0.01) # (Atom 0, x) * (Atom 0, x) - assert np.isclose(disp_matrix[3], 0.1 * 0.05 * 2) # (Atom 0, x) * (Atom 1, x) * 2 (harmonic) - assert np.isclose(disp_matrix[6 + 5 + 4], 0.0025) # (Atom 1, x) * (Atom 1, x) + disp_matrix = get_model_inputs(atoms, reference) + assert (disp_matrix != 0).sum() == 2 + 3 + assert np.isclose(disp_matrix[[0, 3]], [0.1, 0.05]).all() # Linear terms + assert np.isclose(disp_matrix[6], 0.01 / 2) # (Atom 0, x) * (Atom 0, x) + assert np.isclose(disp_matrix[6 + 3], 0.1 * 0.05) # (Atom 0, x) * (Atom 1, x) * 2 (harmonic) + assert np.isclose(disp_matrix[6 + 6 + 5 + 4], 0.0025 / 2) # (Atom 1, x) * (Atom 1, x) def test_linear_model(train_set): @@ -31,12 +33,17 @@ def test_linear_model(train_set): assert reference.get_forces().max() < 0.01 # Fit the model - model = LinearHessianModel(reference) + model = HarmonicModel(reference) hessian_model = model.train(train_set) + assert hessian_model.coef_.shape == (54,) + + # Get the mean hessian + hessian = model.mean_hessian(hessian_model) + assert hessian.shape == (9, 9) # Sample the Hessians, at least make sure the results are near correct hessians = model.sample_hessians(hessian_model, num_samples=32) - assert len(hessians) + assert len(hessians) == 32 assert np.isclose(hessians[0], hessians[0].T).all() # Create a vibration data object