diff --git a/magpy/core/function_product.py b/magpy/core/function_product.py index d692e55..9232427 100644 --- a/magpy/core/function_product.py +++ b/magpy/core/function_product.py @@ -1,4 +1,6 @@ from numbers import Number +from copy import deepcopy +import torch import magpy as mp @@ -29,7 +31,6 @@ def __init__(self, *funcs): self.funcs = {} self.scale = 1 - self.__name__ = "FP" for f in funcs: try: @@ -42,7 +43,7 @@ def __init__(self, *funcs): self.scale *= f except TypeError: # Other type of function. - self.funcs = self.__add_func(f) + self.funcs = FunctionProduct.__add_func(self, f) def __eq__(self, other): return self.funcs == other.funcs and self.scale == other.scale @@ -52,17 +53,17 @@ def __mul__(self, other): return other * self out = FunctionProduct() - out.funcs = self.funcs.copy() - out.scale = self.scale + out.funcs = deepcopy(self.funcs) + out.scale = deepcopy(self.scale) - if isinstance(other, Number): + if isinstance(other, Number | torch.Tensor): out.scale *= other else: try: out.scale *= other.scale out.funcs = self.__merge_funcs(other.funcs) except AttributeError: - out.funcs = out.__add_func(other) + out.funcs = FunctionProduct.__add_func(out, other) return out @@ -74,24 +75,37 @@ def __neg__(self): def __call__(self, arg): out = 1 for f in self.funcs: - out *= f(arg) + try: + out *= f(arg.clone().detach()) + except AttributeError: + out *= f(torch.tensor(arg)) return out * self.scale + def __repr__(self): + return f"{str(self.scale)}*{str(self.funcs)}" + + def __str__(self): + return (str(self.scale) + "*" if isinstance(self.scale, torch.Tensor) or self.scale != 1 else "") \ + + '*'.join([f.__name__ + (f"^{str(n)}" if n > 1 else "") for f, n in self.funcs.items()]) + def __hash__(self): return hash(tuple(self.funcs)) + hash(self.scale) - def __repr__(self): - return str(self.scale) + "*" + str(self.funcs) + def is_empty(self): + """Return true if function product contains no functions. + """ + return not self.funcs def __merge_funcs(self, funcs): # Combine funcs dict with own funcs dict, summing values with shared keys. return {f: self.funcs.get(f, 0) + funcs.get(f, 0) for f in set(self.funcs) | set(funcs)} - def __add_func(self, f): - # Add function to own funcs dict, adding new key or incremented existing value accordingly. + @staticmethod + def __add_func(out, f): + # Add function to funcs dict, adding new key or incrementing existing value accordingly. try: - self.funcs[f] += 1 + out.funcs[f] += 1 except KeyError: - self.funcs[f] = 1 - return self.funcs + out.funcs[f] = 1 + return out.funcs diff --git a/magpy/core/hamiltonian_operator.py b/magpy/core/hamiltonian_operator.py index 9f602de..c54eb9d 100644 --- a/magpy/core/hamiltonian_operator.py +++ b/magpy/core/hamiltonian_operator.py @@ -1,4 +1,3 @@ -from itertools import chain from copy import deepcopy from numbers import Number import torch @@ -46,7 +45,7 @@ def __init__(self, *pairs): except AttributeError: self.data[pair[0]] = [self.data[pair[0]], pair[1]] - HamiltonianOperator.__simplify(self.data) + self.data = HamiltonianOperator.__simplify_and_sort_data(self.data) def __eq__(self, other): return self.data == other.data @@ -55,14 +54,7 @@ def __mul__(self, other): out = HamiltonianOperator() try: - self_data = HamiltonianOperator.__expand(self.data) - other_data = HamiltonianOperator.__expand(other.data) - - for p in self_data: - for q in other_data: - out += mp.FunctionProduct() * p[0] * q[0] * p[1] * q[1] - - return out + out = sum((p[0]*p[1]*q[0]*q[1] for q in other.unpack_data() for p in self.unpack_data()), out) except AttributeError: out.data = deepcopy(self.data) @@ -80,7 +72,8 @@ def __mul__(self, other): for coeff in list(out.data): out.data[mp.FunctionProduct(coeff, other)] = out.data.pop(coeff) - return out + out.data = HamiltonianOperator.__simplify_and_sort_data(out.data) + return out __rmul__ = __mul__ @@ -91,7 +84,7 @@ def __add__(self, other): out.data = self.data | other.data except AttributeError: # other is PauliString; add it to constants. - out.data = self.data.copy() + out.data = deepcopy(self.data) try: out.data[1].append(other) @@ -114,7 +107,7 @@ def __add__(self, other): except TypeError: out.data[coeff].append(other.data[coeff]) - HamiltonianOperator.__simplify(out.data) + out.data = HamiltonianOperator.__simplify_and_sort_data(out.data) return out def __neg__(self): @@ -127,42 +120,60 @@ def __sub__(self, other): return out def __repr__(self): - return '{' + ', '.join((str(f) if isinstance(f, Number) - else f.__name__) + ': ' + str(q) for f, q in self.data.items()) + '}' + return str(self.data) + + def __str__(self): + out = "" + for f, p in self.data.items(): + try: + p_str = str(p) + scale_pos = p_str.find('*') + + if isinstance(p.scale, torch.Tensor) or p.scale != 1: + out += p_str[:scale_pos] + '*' + + out = HamiltonianOperator.__add_coeff_to_str(out, f) + out += p_str[scale_pos + 1:] if scale_pos > 0 else p_str - def __call__(self, t=None, n=None): - if n is None: - pauli_strings = chain.from_iterable(p if isinstance(p, list) else [p] for p in self.data.values()) - n = max(max(p.qubits) for p in pauli_strings) + except AttributeError: + out = HamiltonianOperator.__add_coeff_to_str(out, f) + + out += '(' if f != 1 else "" + out += " + ".join([str(q) for q in p]) + out += ')' if f != 1 else "" + + out += " + " + + return out[:-3] + + def __call__(self, t=None, n_qubits=None): + if n_qubits is None: + n_qubits = max(max(p.qubits) if p.qubits else 0 for p in self.pauli_operators()) if self.is_constant(): - try: - return sum(p(n).type(torch.complex128) for p in self.data[1]) - except TypeError: - return self.data[1](n).type(torch.complex128) + return self.__call_time_independent(n_qubits) if t is None: raise ValueError( "Hamiltonian is not constant. A value of t is required.") - out = 0 - for coeff, ps in self.data.items(): - try: - out += coeff(torch.tensor(t)) * ps(n).type(torch.complex128) - except TypeError: - if isinstance(coeff, Number): - out += coeff * ps(n).type(torch.complex128) - else: - for p in ps: - out += coeff(torch.tensor(t)) * p(n).type(torch.complex128) + # Convert input to tensor. + try: + t = t.clone().detach() + except AttributeError: + t = torch.tensor(t) - return out.to(_DEVICE_CONTEXT.device) + return self.__call_time_dependent(t, n_qubits) def is_constant(self): "Return true if the Hamiltonian is time-independent." for coeff in self.data: - if not isinstance(coeff, Number): - return False + if not isinstance(coeff, Number | torch.Tensor): + try: + if not coeff.is_empty(): + return False + except AttributeError: + return False return True def is_interacting(self): @@ -190,20 +201,108 @@ def pauli_operators(self): """All Pauli operators in H.""" return [u[1] for u in self.unpack_data()] + @property + def n_qubits(self): + """Number of qubits in operator.""" + return max(p.n_qubits for p in self.pauli_operators()) + + def qutip(self): + """Convert to QuTiP form.""" + if self.is_constant(): + return self.__qutip_constant(self.unpack_data()) + + return [self.__qutip_constant(self.constant_components())] \ + + [[c[1].qutip(self.n_qubits), c[0]] for c in self.time_dependent_components()] + + def __call_time_independent(self, n_qubits): + # Return matrix representation when constant. + out = 0 + for p in self.pauli_operators(): + try: + p_val = p(n_qubits) + + # If p is a batch, repeat the current value to agree with its shape. + if out.dim() == 2 and p_val.dim() == 3: + out = out.repeat(len(p_val), 1, 1) + + except AttributeError: + pass + + out += p(n_qubits) + + return out + + def __call_time_dependent(self, t, n_qubits): + out = 0 + for coeff, p in self.unpack_data(): + p_val = p(n_qubits).to(_DEVICE_CONTEXT.device) + + # Evaluate coefficient if it's a function. + try: + coeff = coeff(t).to(_DEVICE_CONTEXT.device) + except TypeError: + pass + + # Evaluate next term in data. + next_term = 0 + try: + next_term = coeff.reshape(-1,1,1) * p_val + except AttributeError: + next_term = coeff * p_val + + # If p is a batch, repeat the current value to agree with its shape. + try: + if out.dim() == 2 and next_term.dim() == 3: + out = out.repeat(len(next_term), 1, 1) + except AttributeError: + pass + + out += next_term + + return out + + def constant_components(self): + return [state for state in self.unpack_data() if isinstance(state[0], Number | torch.Tensor)] + + def time_dependent_components(self): + return [state for state in self.unpack_data() if state not in self.constant_components()] + @staticmethod - def __simplify(arrs): + def __simplify_data(arrs): # Collect all PauliStrings in all lists in arrs. for coeff in arrs: arrs[coeff] = mp.PauliString.collect(arrs[coeff]) @staticmethod - def __expand(data): - # Expand all functions and lists of qubits into pairs of functions with single qubits. - expanded_data = [] - for pair in data.items(): + def __add_coeff_to_str(out, f): + try: + out += f.__name__ + '*' + except AttributeError: try: - for qubit in pair[1]: - expanded_data.append((pair[0], qubit)) - except TypeError: - expanded_data.append(pair) - return expanded_data + if f != 1: + out += str(f) + '*' + except (RuntimeError, AttributeError): + out += str(f) + '*' + + return out + + @staticmethod + def __sort_data(data): + # Move all constant keys to the start of the dictionary. + const_keys = [] + other_keys = [] + + for key in data: + (const_keys if isinstance(key, Number | torch.Tensor) else other_keys).append(key) + + return dict((key, data[key]) for key in const_keys + other_keys) + + @staticmethod + def __simplify_and_sort_data(data): + HamiltonianOperator.__simplify_data(data) + return HamiltonianOperator.__sort_data(data) + + @staticmethod + def __qutip_constant(states): + n = max(max(states[i][1].qubits) for i in range(len(states))) + return sum(state[0] * state[1].qutip(n) for state in states) diff --git a/magpy/core/pauli_string.py b/magpy/core/pauli_string.py index 3bea981..89b3d44 100644 --- a/magpy/core/pauli_string.py +++ b/magpy/core/pauli_string.py @@ -1,6 +1,8 @@ from numbers import Number +from copy import deepcopy import torch import magpy as mp +import qutip as qt from .._device import _DEVICE_CONTEXT @@ -45,6 +47,13 @@ class PauliString: 'Z': torch.tensor([[1, 0], [0, -1]]) } + qutip_states = { + 'X' : qt.sigmax(), + 'Y' : qt.sigmay(), + 'Z' : qt.sigmaz(), + 'Id' : qt.qeye(2) + } + def __init__(self, x=None, y=None, z=None, scale=1): """A multi-qubit Pauli operator. @@ -70,61 +79,42 @@ def __eq__(self, other): return self.qubits == other.qubits and self.scale == other.scale def __mul__(self, other): + if isinstance(other, mp.PauliString): + return self.__mul_ps(other) + if isinstance(other, mp.HamiltonianOperator): - return -other * self + return self.__mul_hop(other) - s = PauliString(scale=0) + if isinstance(other, Number | torch.Tensor): + return self.__mul_num(other) try: - if other.scale == 0: - return s - - s.qubits = self.qubits | other.qubits + # other is FunctionProduct. + self *= other.scale + other.scale = 1 except AttributeError: - if isinstance(other, Number): - if other == 0: - return s + # other is another type of function + pass - s.scale = self.scale * other - s.qubits = self.qubits - else: - try: - # other is FunctionProduct. - self *= other.scale - other.scale = 1 - except AttributeError: - # other is other type of function. - pass + return mp.HamiltonianOperator([other, self]) - return mp.HamiltonianOperator([other, self]) - - else: - # other is PauliString. - s.scale = self.scale * other.scale - - for n in list(set(self.qubits.keys() & other.qubits.keys())): - if self.qubits[n] == other.qubits[n]: - del s.qubits[n] - else: - scale, spin = PauliString.__pauli_mul(self.qubits[n], other.qubits[n]) - s.scale *= 1j * scale - s.qubits[n] = spin - - return s - - __rmul__ = __mul__ + __rmul__ = __mul__ # This may have to be changed def __add__(self, other): try: - if self == other: - s = PauliString(scale=self.scale + other.scale) - s.qubits = self.qubits - return s - if self == -other: - return 0*PauliString.Id() + if self.qubits == other.qubits: + scale = self.scale + other.scale + if torch.is_tensor(scale) and torch.all(scale.eq(0)) or isinstance(scale, Number) and scale == 0: + return 0 + + out = PauliString(scale=scale) + out.qubits = self.qubits + return out + return mp.HamiltonianOperator([1, self], [1, other]) + except AttributeError: - # other is HamiltonianOperator + # other is HOp return other + self def __neg__(self): @@ -133,10 +123,20 @@ def __neg__(self): return s def __sub__(self, other): - return self + -1*other + return self + -other def __repr__(self): - return str(self.scale) + "*" + "*".join(q[1] + str(q[0]) for q in sorted(self.qubits.items())) + out = "" + + if isinstance(self.scale, torch.Tensor) or self.scale != 1: + out += str(self.scale) + "*" + + if self.qubits.items(): + out += "*".join(q[1] + str(q[0]) for q in sorted(self.qubits.items())) + else: + out += "Id" + + return out def __call__(self, n=None): if n is None: @@ -146,7 +146,49 @@ def __call__(self, n=None): for index, qubit in self.qubits.items(): qubits[index - 1] = PauliString.matrices[qubit] - return self.scale * mp.kron(*qubits).type(torch.complex128).to(_DEVICE_CONTEXT.device) + scale = self.scale.view(len(self.scale), 1, 1).to(_DEVICE_CONTEXT.device) \ + if isinstance(self.scale, torch.Tensor) else self.scale + return scale * mp.kron(*qubits).type(torch.complex128).to(_DEVICE_CONTEXT.device) + + def __mul_ps(self, other): + # Right multiply by PauliString + out = PauliString() + out.scale = self.scale * other.scale + out.qubits = self.qubits | other.qubits + + for n in list(set(self.qubits.keys() & other.qubits.keys())): + if self.qubits[n] == other.qubits[n]: + del out.qubits[n] + else: + scale, spin = PauliString.__pauli_mul(self.qubits[n], other.qubits[n]) + out.scale *= 1j * scale + out.qubits[n] = spin + + return out + + def __mul_hop(self, other): + # Right multiply by Hamiltonian + return mp.HamiltonianOperator([1, self]) * other + + def __mul_num(self, other): + # Multiply by scalar + out = deepcopy(self) + out.scale *= other + + return out + + @property + def n_qubits(self): + return max(self.qubits) if self.qubits else 1 + + def qutip(self, n=None): + """Convert to QuTiP form.""" + if not self.qubits: + return qt.qeye(2 ** n) + + n = (self.n_qubits if n is None else n) + 1 + + return self.scale * qt.tensor(self.qutip_states.get(self.qubits.get(i, 'Id')) for i in range(1, n)) @staticmethod def collect(arr): diff --git a/tutorials/quickstart.ipynb b/tutorials/quickstart.ipynb new file mode 100644 index 0000000..8f1df3e --- /dev/null +++ b/tutorials/quickstart.ipynb @@ -0,0 +1,175 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# MagPy quickstart guide\n", + "\n", + "MagPy provides functionality for evolving quantum systems under a magnetic field, utilising the Magnus expansion and the the highly parallelised nature of CUDA computation via PyTorch.\n", + "\n", + "Specifically, we may simulate the evolution of density matrices with respect to the Liouville-von Neumann equation,\n", + "\n", + "$$\\frac{\\partial \\rho(t)}{\\partial t} = -i\\,\\big[H(t),\\,\\rho(t)\\big],$$\n", + "\n", + "with initial condition $\\rho_0.$\n", + "\n", + "## Pauli Strings\n", + "\n", + "These objects are the fundamental building blocks of quantum objects in MagPy. The Pauli spin operators are accessed with the functions `X`, `Y`, and `Z`. These provide a symbolic representation which we may manipulate algebraically. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(-3+0j)*Z1*Y2\n", + "{1: [5*X1, 5*Y1]}\n" + ] + } + ], + "source": [ + "from magpy import X, Y, Z\n", + "\n", + "print(3j * X(1) * Y(1, 2))\n", + "print(5 * (X() + Y()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The matrix representation of these matrices may also be accessed:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{1: [1*X1, 3*Y2]}\n", + "tensor([[0.+0.j, 0.-3.j, 1.+0.j, 0.+0.j],\n", + " [0.+3.j, 0.+0.j, 0.+0.j, 1.+0.j],\n", + " [1.+0.j, 0.+0.j, 0.+0.j, 0.-3.j],\n", + " [0.+0.j, 1.+0.j, 0.+3.j, 0.+0.j]], dtype=torch.complex128)\n" + ] + } + ], + "source": [ + "a = X() + 3*Y(2)\n", + "\n", + "print(a)\n", + "print(a())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Hamiltonians\n", + "\n", + "The functions `X`, `Y`, and `Z`, enable us to write Hamiltonians and density matrices mathematically:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "from magpy import X, Y\n", + "\n", + "H = torch.sin*X() + 4*Y()\n", + "rho0 = Y()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Simulating a quantum system\n", + "\n", + "Given a Hamiltonian and initial condition, we can simulate its evolution over a specified discretisation of time using the `System` class:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import magpy as mp\n", + "\n", + "tlist = mp.timegrid(0, 10, 0.5**6)\n", + "qsys = mp.System(H, rho0, tlist)\n", + "qsys.evolve()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The results of which are stored in the `state` attribute of the `System` instance." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "tensor([[[ 0.0000e+00+0.0000e+00j, 0.0000e+00-1.0000e+00j],\n", + " [ 0.0000e+00+1.0000e+00j, 0.0000e+00+0.0000e+00j]],\n", + "\n", + " [[ 2.4350e-04+0.0000e+00j, 1.5239e-05-1.0000e+00j],\n", + " [ 1.5239e-05+1.0000e+00j, -2.4350e-04+0.0000e+00j]],\n", + "\n", + " [[ 9.7013e-04-6.9389e-18j, 9.1190e-05-1.0000e+00j],\n", + " [ 9.1190e-05+1.0000e+00j, -9.7015e-04+0.0000e+00j]]],\n", + " dtype=torch.complex128)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsys.states[0:3]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/solver/simulating.ipynb b/tutorials/solver/simulating.ipynb index 2597950..f8fcb79 100644 --- a/tutorials/solver/simulating.ipynb +++ b/tutorials/solver/simulating.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Simulating a quantum system (WIP)\n", + "# Simulating a quantum system\n", "\n", "MagPy uses a two-term Magnus expansion to evolve solutions to the Liouville-von Neumann equation,\n", "\n", @@ -43,7 +43,7 @@ "metadata": {}, "outputs": [], "source": [ - "qsys = mp.System.create(H, rho0, tlist)\n", + "qsys = mp.System(H, rho0, tlist)\n", "qsys.evolve()" ] }, @@ -57,25 +57,35 @@ "\n", "$$\\langle A, B\\rangle_\\text{F} := \\text{Tr}\\big(A^\\dagger B\\big),$$\n", "\n", - "we can measure the magnitude of a single spin component, giving us a scalar value to plot over time.\n", + "we can measure the magnitude of a single spin component, giving us a scalar value to plot over time. MagPy provides a batch function which will calculate the inner product of all state matrices with the given spin component matrix.\n", "\n", - "MagPy provides a batch function which will calculate the inner product of all state matrices with the given spin component matrix.\n", + "Note: the double parentheses indicate that we are calling the `Y` function and then calling it's matrix representation.\n", "\n", "For the $Y$ component of the system's spin," ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 3, "metadata": {}, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/user/Documents/code/magpy/.venv/lib/python3.10/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return math.isfinite(val)\n", + "/home/user/Documents/code/magpy/.venv/lib/python3.10/site-packages/matplotlib/cbook.py:1398: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return np.asarray(x, float)\n" + ] + }, { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 6, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, @@ -97,13 +107,6 @@ "plt.plot(tlist, y_component)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: the double parentheses indicate that we are calling the `Y` function and then calling it's matrix representation." - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -115,14 +118,15 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import torch\n", "\n", "mp.set_device('cuda') # Default CUDA device\n", - "mp.set_device('cuda:0') # First available GPU" + "mp.set_device('cuda:0') # First available GPU\n", + "mp.set_device('cpu') # Return to CPU" ] }, { @@ -134,12 +138,69 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "H = torch.sin*X() + 4*Y()" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Multiple qubits\n", + "\n", + "MagPy provides support for multi-qubit systems. This is done by specifying the position of the qubit in each call of `X`, `Y`, or `Z`. The default value is 1, so this parameter can be omitted in single-qubit systems, as above." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/user/Documents/code/magpy/.venv/lib/python3.10/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return math.isfinite(val)\n", + "/home/user/Documents/code/magpy/.venv/lib/python3.10/site-packages/matplotlib/cbook.py:1398: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return np.asarray(x, float)\n" + ] + }, + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABS2UlEQVR4nO3deXhU5d0+8PvMnsm+L2QPEEBW2QSqguC+vNZWbWsVl9oNWxXbCm3V+mst2tq+vlaLy6vV12rV2qq1rSgFBTf2RdZAICEhIfsyk20mmTm/P2bOBCRAljnznHPm/lzXXG0g5HwN4Zn72SVZlmUQERERCWASXQARERFFLwYRIiIiEoZBhIiIiIRhECEiIiJhGESIiIhIGAYRIiIiEoZBhIiIiIRhECEiIiJhLKILOB2/34/a2lrEx8dDkiTR5RAREdEgyLIMt9uNnJwcmEynH/PQdBCpra1FXl6e6DKIiIhoGKqrq5Gbm3vaz9F0EImPjwcQ+A9JSEgQXA0RERENhsvlQl5eXuh9/HQ0HUSU6ZiEhAQGESIiIp0ZzLIKLlYlIiIiYRhEiIiISBgGESIiIhKGQYSIiIiEYRAhIiIiYRhEiIiISBgGESIiIhKGQYSIiIiEYRAhIiIiYSIWRB5++GFIkoS77rorUo8kIiIijYtIENm8eTOefvppTJ48ORKPIyIiIp1QPYh0dHTghhtuwLPPPovk5GS1H0dEREQ6ovqld0uWLMHll1+ORYsW4Ve/+tVpP9fj8cDj8YQ+drlcqtRU3tCBlz6rhNlkgtkEmEwSzJIEs0mCSZJgMUkwmSTYLSYkO21IibUhK9GBorRYOKxmVWoiIn1o7+pFdWsXGtw9aHR74O7pg6fPH3z54PPJA/65L979JUmBNsZhNcNhNSPGakZKrA05SQ7kJMUgNdY2qAvDiPRO1SDy6quvYtu2bdi8efOgPn/FihV48MEH1SwJAFDT1o0XPzsy5D9nkoCC1FhMyU3EuWPSMb80HalxdhUqJCKtqGnrxgf7G7D+QCN217Sjtr0nIs9NjbVhUm4izilOxYUTMlGSHheR5xJFmiTL8sDxfYSqq6sxY8YMrF69OrQ2ZP78+Zg6dSoee+yxAf/MQCMieXl5aG9vR0JCQthqq2jqxN+2HoVPluHz97/88on/29PrR2uXFy2dXlS3dMHV03fC17GaJVw0IQs3zSnA7OLUsNVHRGL5/TJW76vHnzccwUcHm076/fR4OzIT7MiIdyDeYQmNbNgtJphNJ894yzi5mZVlwNPrQ3evD929fnR7+9DY4cWxtm40dnjwxZb5rJwE3PalIlwxOQc2Czc8kra5XC4kJiYO6v1btSDy1ltv4ctf/jLM5v6pDJ/PB0mSYDKZ4PF4Tvi9gQzlP0Rtsiyj0e1BWb0bnx1qxodljdh7rH/q6IJxGfjpZeMwOiNeYJVENFJbj7Tg/72zFzuPtgMITKnMKEjG/NIMzCpKQWlWPBIcVlVr6On1YX+dG9urWrF2fwM2HG5Gb3DKZ1RSDH562XhcNimLUzekWZoIIm63G0eOnDj9ccstt2DcuHG49957MXHixDN+DS0FkYHsrXXhpQ1H8Nct1ejzy7CYJPzo4lJ8+9ximExsIIj0pKfXh4ff3Y8XPq0EAMTazLhxTiFumJ2PvBSn0NpaO714ZVMVXvy0Eg3uwKjxRRMy8fBXJiMl1ia0NqKBaCKIDORMUzNfpPUgojjc2IGH/rUPa/Y3AAAWjc/E41+fCqdN9bXARBQGDe4efOvFLfg8OApy/Yw83HPxWGTEOwRXdqJurw9PrTuEP35Yjl6fjFFJMXj+5pkozeJILGnLUN6/OdEYBsXpcfjfxTPw8DWTYLOY8J999fj6MxvQ3t0rujQiOoOq5i58+clP8fnRdiQ7rfjTLTPxyFcnay6EAECMzYy7LxyLN78/D4WpTtS0deOrKz/F1iMtoksjGraIjogMlV5GRI639UgLvvXiFrR29WJKXhL+fNssxKs8n0xEw1Pb1o1rn/oMNW3dKEx14oVbZqEwLVZ0WYPS2unFd17aik2VLYizW/Dnb83G1Lwk0WURAeCIiFDTC1Lwyu3nIMlpxc7qNvzwL9vh82s26xFFrQ5PH27502bUtHWjKC0Wr393jm5CCAAkx9rw4q2zMLsoBR2ePnzrxc042toluiyiIWMQUcH47AT8362z4LCa8EFZI36zar/okojoOH6/jLtf24Gyejcy4u14+VuzNTkVcyYxNjOev3kmJmQnoKnDi2+9uAXdXp/osoiGhEFEJZNzk/DotVMAAE+vP4wPyxoEV0REihc+rcTqvfWwWUx4+sbpyEmKEV3SsMXaLfjfxTOQFmfH/jo3fvWvvaJLIhoSBhEVXTE5BzfPLQQA/OivO9Hc4Tn9HyAi1ZXVufFwcJTyvsvHY1q+/u/AykmKwX9fH+j4vLyxCqv31guuiGjwGERUtuzScSjNjEdThxcP/Xuf6HKIoprfL2PZ3z+Ht8+PBaXp+OY5BaJLCptzx6Tj2+cVAwDue2s33D3ctUf6wCCiMofVjIe/MgmSBPx9Ww0+PXTycdFEFBmvbq7G9qo2xNrMWHHNZMOdTLr0wrEoSHWiztWD371/QHQ5RIPCIBIB0/KTccPsfADAg//Yy100RAK0d/fiN+8FpmTuuagUWYn6W5x6Jg6rGb+6OnBq9UsbjqC8wS24IqIzYxCJkB9fNA4JDgvK6t34+7ajosshijpPrTuEtq5ejMmIw01zjDMl80XnjknHhRMy4fPLePjdMtHlEJ0Rg0iEJDqtuOOC0QCA371/AD293GJHFCl17T14/uMKAMC9l4yDxWzspu/eS8bBbJLwn3312FzJU1dJ24z9r1FjbppTiJxEB+pcPfjrlmrR5RBFjWc/OgxPnx8zCpKxcHyG6HJUNzojDtfNyAUAPL7moOBqiE6PQSSCHFYzvnN+CQDgqXWH0evzC66IyPjaurz4y6YqAMAPFo4x3ALVU/n+/NEwmyR8dLAJO6rbRJdDdEoMIhF2/cw8pMXZUdPWjbd31Iouh8jwXvrsCLq8PozPTsB5Y9JElxMxeSlOXD11FADgyQ/KBVdDdGoMIhHmsJpx65cKAQDPf1wBDd85SKR7Pb0+vPBpJQDgu+cXR81oiOJ78wMjsP/ZV4+qZt5DQ9rEICLAN2blw2E1Ye8xFzZXtoouh8iw/r6tBs2dXuQmx+DySdmiy4m40RlxOG9sOmQZ+L/PKkWXQzQgBhEBkpw2fHlaYMj0hU8rBFdDZFyvbDoCALh5bqHhd8qcyi3zCgEAr22pRqenT2wxRAOIzn+ZGrA4eAfN+3vq0cQ7aIjCbndNO3bXuGAzm3DN2bmiyxHm/DHpKEx1wt3Th3/vOia6HKKTMIgIMi4rAVNyE9Hnl7lolUgFr24O7JS5eGIWUmJtgqsRx2SScO2MPADAG1t5mCJpD4OIQF+dHuilsXEgCq8ubx/e3h4I+F+bmSe4GvG+PG0UJAnYWNHCRaukOQwiAl01ZRRsZhP2HXNhd0276HKIDOPfu+rg9vQhP8WJOcWpossRLicpBl8aHdi6/MZWHqZI2sIgIlCi04oLz8oEwFERonB6a3sNAOC6GbkwmaJry+6pKCOwf9tWAz8v3iQNYRARTGkc3t5Rw5NWicKgucODTw81AQCunJIjuBrtuPisLMQ7LKhp68aGimbR5RCFMIgIdu7oNKTG2tDa1YsNh9k4EI3Uqj118MvApFGJKEiNFV2OZjisZlxyVhYAYNXuOsHVEPVjEBHMYjbhouD0zLtsHIhG7J87A1tUr5gcfQeYncllwUPdVu2u4/QMaQaDiAZcMjHQOLy/pw4+Ng5Ew9bg7sHG4LTDZVF4kuqZzB2dini7BQ1uD7ZV8VRn0gYGEQ2YU5yKBIcFTR1ebD3CxoFouN7bHZiWmZKXhLwUp+hyNMduMWPh+AwAHIEl7WAQ0QCbxYRFE5TpGZ58SDRcypvr5ZOyBFeiXZceNz3DSzdJCxhENOLS4PTMe2wciIbF3dOLTRUtAIALJzCInMr5Y9PhtJlR09aNXTy/iDSAQUQjzh2ThlibGbXtPdhT6xJdDpHufHywCX1+GcVpsShK426ZU3FYzZhfmg4A+M++BsHVEDGIaIbDasbc4MmH6w40Cq6GSH/W7g+8qS4YlyG4Eu2bXxr4Hq0rYxAh8RhENETppXzIxoFoSPx+GR8E/91cwCByRvPHBtqaz2va0czbv0kwBhENOW9MoHHYVtWG9u5ewdUQ6ceumnY0dXgRZ7dgZmGK6HI0LyPBgQnZCZBlYP1BjsCSWKoGkZUrV2Ly5MlISEhAQkIC5syZg3fffVfNR+paXooTJemx8PllfFLeJLocIt1QpmW+NDoNNgv7V4NxfmgElkGExFL1X2xubi4efvhhbN26FVu2bMEFF1yA//qv/8KePXvUfKyu9c/dsnEgGqwPg+uqOC0zeMr0zPoDjTxIkYRSNYhceeWVuOyyyzBmzBiMHTsWDz30EOLi4rBhwwY1H6tr5wcbh3UHGrmNl2gQ2rt7setoGwDg3LFpYovRkbMLkhFvt6C1qxefB79/RCJEbAzT5/Ph1VdfRWdnJ+bMmTPg53g8HrhcrhNe0WZWUQocVhPqXD0oq3eLLodI8zZVtMAvA8VpschOjBFdjm5YzSbMC+7UW3+AU8EkjupBZNeuXYiLi4Pdbsd3v/tdvPnmm5gwYcKAn7tixQokJiaGXnl5eWqXpzkOqxmzilIBAJ+W8zZeojNR1lPNHZ0quBL9mTcmEER48zeJpHoQKS0txY4dO7Bx40Z873vfw+LFi7F3794BP3f58uVob28Pvaqrq9UuT5POKQ6s+mfjQHRmnx0K/DuZW8JpmaGaUxwIb1urWtHT6xNcDUUri9oPsNlsGD16NABg+vTp2Lx5M/7nf/4HTz/99Emfa7fbYbfb1S5J85TGYWNFC/x+GSaTJLgiIm1qdHtCU5jnFHNEZKhK0mORHm9Ho9uDHdVt/B6SEBHf5+b3++Hx8ACd05k4KhGxNjPau3uxv47rRIhO5bPgqOGE7ASkxNoEV6M/kiSFwocyskQUaaoGkeXLl2P9+vWorKzErl27sHz5cnz44Ye44YYb1Hys7lnNJswo5PQM0Zl8dii4PqSEPfnh4lQwiaZqEGloaMBNN92E0tJSLFy4EJs3b8Z7772HCy+8UM3HGoLSS2HjQHRqnyrrQ7hQddiUqeDtVW1cJ0JCqLpG5LnnnlPzyxua0kvhOhGigTW4enCkuQuShNAIIg1dUVosMuLtaHB7sK2qlYt+KeJ4FrJGTeI6EaLT2nKkFQAwLisBCQ6r4Gr0S5IkzAlObW3gOhESgEFEoyzHrRPZVMHGgeiLtlQGgsiMgmTBlejf7ODZRUq4I4okBhENUxrYbVVtYgsh0qCtR1oAADMKGURGanqwrdlR3YY+n19wNRRtGEQ07Oxg47CVvRSiE3R5+7C7NnAFBNeHjNyYjDjE2y3o8vp4tQRFHIOIhk3JS4JJAmraulHv6hFdDpFm7Khug88vIzvRgVFJvF9mpEwmCVPzkwAA29jxoQhjENGwOLsFpVkJANg4EB1va3B9yHSuDwmb6RyBJUEYRDRuekESAGBbFRsHIsXm4JvlTE7LhM10rkkjQRhENO7sfPZSiI7n98vYfoQjIuE2NS8JkgRUtXShwc2pYIocBhGNU4LI7hoXPH089ZDoUGMH3J4+xFjNGJcVL7ocw4h3WFGaGfh+bjvSJrYYiioMIhpXkOpEaqwNXp8fu2tcosshEm7n0XYAgUP/LGY2YeGk7NTbzqlgiiD+K9Y4SZIwLTgqwgWrRMDO6jYAwJS8RLGFGJAyAss1aRRJDCI6MC24re7zmnaxhRBpwOdH2wAAk3OThNZhRFOD4W5PrQs+vyy4GooWDCI6MGlUoHHYFWyAiaKVp8+HvccCU5RT85LEFmNARWlxiLWZ0eX14VBjh+hyKEowiOiAEkQqm7vQ3tUruBoicfYdc6PXJyPZaUVuMg8yCzezScJZwfbm86McgaXIYBDRgeRYG/JTnACA3bVsHCh6KdMyU/KSIEmS2GIMakquEkTaxBZCUYNBRCcm5bKXQrQjuFCV60PUMyn4vWVbQ5HCIKITk0exl0KkvDlO5Y4Z1Shtzd5jLvTyJl6KAAYRneCICEU7d09vaAElR0TUU5DqRILDAm+fH2V1vImX1McgohMTg72UmrZuNHd4BFdDFHl7al2QZSAn0YG0OLvocgxLkqRQ0NvFIwMoAhhEdCLBYUVxWiwANg4UnfbUBrbtKrs6SD0cgaVIYhDRETYOFM32BHeMnZWTILgS4+OaNIokBhEdUc4T2cMtvBSF9iojIjkcEVGbMhV8sL4D3j4uWCV1MYjoyIRgT1A5WZIoWvT0+nCwIbBQlSMi6stNjkG83QKvz88TVkl1DCI6MiE70ABXt3SjvZsnrFL0OFDvhs8fOFE1O9EhuhzDkyQJ45WOTy07PqQuBhEdSXLakBNshPdzVISiiLJQdUJOAk9UjRCl48MRWFIbg4jOKNMz+9g4UBTpX6jK9SGRMoEjIhQhDCI6w14KRaPQ1l2uD4mY49saWZYFV0NGxiCiM1ywStHG55ex/1jghE8GkcgZkxkHi0lCe3cvatt7RJdDBsYgojMTsgND0wfqOngPBEWFiqYOdPf6EGM1oygtTnQ5UcNuMWN0RuD7zekZUhODiM5wWx1Fm33B0ZDSrHiYTVyoGklcJ0KRwCCiMyaThPHZbBwoeigXr43LihdcSfTpXyfCQxRJPaoGkRUrVmDmzJmIj49HRkYGrr76apSVlan5yKjAnTMUTfbX9Y+IUGRxTRpFgqpBZN26dViyZAk2bNiA1atXo7e3FxdddBE6OzvVfKzhjc8ONMhsHCgaHKhnEBFlfFb/IYodnj7B1ZBRWdT84qtWrTrh4xdeeAEZGRnYunUrzjvvPDUfbWilwcbhQD3XiJCxdXr6UNXSBQAozWQQibTkWBsy4u1ocHtwsN6NafnJoksiA4roGpH29sA8Y0pKyoC/7/F44HK5TnjRycYEV7I3uj1o6fQKroZIPcpoSFqcHalxdsHVRKexwQCo/F0QhVvEgojf78ddd92FefPmYeLEiQN+zooVK5CYmBh65eXlRao8XYm1W5CXEgOAjQMZm/LzzYWq4ihBpKyOI7CkjogFkSVLlmD37t149dVXT/k5y5cvR3t7e+hVXV0dqfJ0p5S9FIoCXKgqXmlWYASWbQ2pRdU1Ioo77rgD//znP7F+/Xrk5uae8vPsdjvsdg6/DsaYzHj8Z19DaGsjkRGVMYgIFxoRYRAhlag6IiLLMu644w68+eabWLt2LYqKitR8XFRRRkQOcsEqGVgoiHChqjBjgt/7RrcHrVyTRipQNYgsWbIEf/7zn/HKK68gPj4edXV1qKurQ3d3t5qPjQrH91J4IRUZUVOHB82dXkhS/887RV6c3YLcZK5JI/WoGkRWrlyJ9vZ2zJ8/H9nZ2aHXa6+9puZjo0JxeizMwQupGtwe0eUQhZ0yGlKQ4kSMzSy4mujGNWmkJlXXiLCnrh6H1YzCVCcONXairM6NzASH6JKIwooLVbVjbFY81uxv4DoRUgXvmtExpYFmL4WM6ADXh2hGaESEW3hJBQwiOjYmg0GEjGt/6Gj3BMGVENekkZoYRHRMGREp484ZMhi/X8ZB3jGjGVyTRmpiENGxsaEtvG74/eylkHHUtHWjy+uDzWxCYapTdDlRT1mTBoBnF1HYMYjoWGGqEzazCV1eH2rauCWajKO8ITDKV5QWC4uZzZQWcE0aqYX/wnXMYjahJHgBHnspZCRKEBkd/Pkm8frvnGFbQ+HFIKJzpZnBIMJeChmIEkRKGEQ0g2eJkFoYRHROOX75UAMXrJJxlDdyRERrlL+LQ42d3DlDYcUgonMl6UrjwCBCxiDLcv/UTDqDiFYUpAZ2znR4+rhzhsKKQUTnRmfEAmAvhYyjqcOL9u5eSFJg2yhpg81iQkFKYOdMOUdgKYwYRHQuP6W/l1LvYi+F9E95k8tLdsJh5R0zWlLMEVhSAYOIztksJhSkspdCxsH1Idql/J2wraFwYhAxAK4TISMpD+7KYBDRnpJ0ZSqYbQ2FD4OIATCIkJGERkS4UFVzSjgiQipgEDEADpeSkYR2zGQyiGiN0umpd3ng7ukVXA0ZBYOIAXC4lIzC1dMbWnTNqRntSYyxIj3eDiCwU48oHBhEDEAZLq13eeBiL4V0TDmYLyPejgSHVXA1NBBlyoyHKFK4MIgYQILDioxgL+UweymkY7xjRvtKgmcXlXMElsKEQcQgSthLIQPg1l3t44gIhRuDiEGEFqyyl0I6dogjIppXwraGwoxBxCBCC1bZSyEd4x0z2qeExKrmLvT6/IKrISNgEDEI9lJI77x9flS3dgPo/3km7clKcCDWZkafX8aR5i7R5ZABMIgYBHsppHfVrV3w+WU4bebQ4mvSHkmSeLAZhRWDiEGwl0J6VxHc8VWUFgtJkgRXQ6fD05wpnBhEDIK9FNK7iqb+IELaxkMUKZwYRAykONiAKw06kZ4cDv7cFjOIaF5RWqDTU8m2hsKAQcRACkNBhL0U0h/l57YonUFE64rY6aEwYhAxEKVxqGziGhHSn8OhNSLcMaN1hWlOAEBrVy9aO72CqyG9YxAxkOJgA36YvRTSmQ5PHxrcgcvuilI5IqJ1TpsFWQkOAEBFM9sbGhkGEQNReilNHbyim/RFWWuQGmtDopOX3elB/wgsgwiNDIOIgcQ7rEiLC5y/wOkZ0pPQQlWuD9ENZS0P14nQSKkaRNavX48rr7wSOTk5kCQJb731lpqPIwBFwVGRw1ywSjpy/BkipA/KFBqngmmkVA0inZ2dmDJlCp588kk1H0PH4YJV0qPQjhkuVNUNTs1QuFjU/OKXXnopLr30UjUfQV+gNOTcwkt6wsPM9Of4qRlZlnkaLg2bqkFkqDweDzweT+hjl8slsBp9UqZmOG9LeiHLMteI6FBeshNmk4Qurw8Nbg8yg7toiIZKU4tVV6xYgcTExNArLy9PdEm60z8iEuilEGldc6cX7p4+SBKQn+IUXQ4Nks1iQm5yDID+M2CIhkNTQWT58uVob28Pvaqrq0WXpDsFqU5IEuDq6UMLDxoiHVBG70YlxcBhNQuuhoYitE6EZ4nQCGgqiNjtdiQkJJzwoqFxWM3ISQz0Ujg9Q3pwuFFZqMppGb3hUe8UDpoKIhQebBxIT3jZnX4pbQ2nZmgkVF2s2tHRgfLy8tDHFRUV2LFjB1JSUpCfn6/mo6NaUVosPi5vYhAhXeAZIvrFqRkKB1WDyJYtW7BgwYLQx0uXLgUALF68GC+88IKaj45qhRwRIR0Jbd1N5xkieqMEkSPNnfD5ZZhN3MJLQ6dqEJk/fz53bghQzCBCOuHzyzjSHDh8j1Mz+pOTGAObxQRvnx81rd3IT+WuJxo6rhExoOOHS/1+BkHSrtq2bnh9ftgsJuQkxYguh4bIZJJQmMprJWhkGEQMKDc5BhaThJ5eP+pcPaLLITolZaFqYaqTw/o6xaPeaaQYRAzIYjaFDoZi40BaVhHculuYymkZvTr+EEWi4WAQMajQtjo2DqRhlcH1Idwxo1/9N36zraHhYRAxKO6cIT2oagkEES5y1C9lRIRbeGm4GEQMivO2pAfKmxenZvRLaWuOtnbD0+cTXA3pEYOIQfF0VdI6n1/G0ZZuALzsTs/S4myIs1sgy0BVcKqNaCgYRAxKmZqpbu2Cj1t4SYOOtQe27lrNErfu6pgkSSgITq0dYRChYWAQMajsBAdsFhN6fTJq27pFl0N0EqX3nJfMrbt6FwoiLQwiNHQMIgZlMknISw70MtlLIS1SdswUcKGq7hWk9h/1TjRUDCIGpiwA5Gp20qIjLYGfywIuVNU95XTVSnZ6aBgYRAxM2RJZxeFS0qAjTRwRMYr8lECYrGKnh4aBQcTAQiMi3DlDGqSsJ2AQ0b/C4KFmR1u70evzC66G9IZBxMC4kp20Spbl0HoCTs3oX2a8A3aLCX1+Lo6noWMQMbDQArKWTsgyt/CSdjR2eNDl9UGSApc0kr6ZTFLoLBh2fGioGEQMbFRSDMzBW3gb3B7R5RCFKFt3cxJjYLeYBVdD4cCdMzRcDCIGZrOYkJPkAMB1IqQt3LprPNw5Q8PFIGJwhaHpGTYOpB1VXB9iOFyTRsPFIGJw/Y0DR0RIO5RecyFHRAyDUzM0XAwiBleQojQO7KWQdnDrrvEoo69VLV3w834rGgIGEYPjcClpEbfuGk9OkgMWkwRPnx/17h7R5ZCOMIgYXMFxx7xzCy9pQXtXL9q6egEgtOWT9M9iNoW2Ylc2seNDg8cgYnBKQ+/u6Qs1/kQiKXfMpMfbEWu3CK6GwqkgND3DdSI0eAwiBhdjMyMrIbiFl4vISAOUacICjoYYTgG38NIwMIhEgXyuEyEN4foQ4+LOGRoOBpEoUMggQhrCw8yMi20NDQeDSBRgL4W0pIpBxLCO36XHxfE0WAwiUSDUOPB0VdKASk7NGFZushOSBHR4+tDc6RVdDukEg0gUKOSICGlEl7cvdAEjT1U1HofVjJzEwBZeTs/QYDGIRAFlsWpThxcdnj7B1VA0qwqOyiU4LEhy2gRXQ2pQjgxgx4cGi0EkCiQ4rEiJDTT6bBxIJKWXXJjGaRmjKkzjFl4amogEkSeffBKFhYVwOByYPXs2Nm3aFInH0nF41DtpAbfuGl/oUDN2emiQVA8ir732GpYuXYoHHngA27Ztw5QpU3DxxRejoaFB7UfTcZTDo3ioGYnEw8yMr7+tYaeHBkf1IPL73/8et99+O2655RZMmDABTz31FJxOJ55//nm1H03H6e+lsHEgcY5w667h8bgAGipVg4jX68XWrVuxaNGi/geaTFi0aBE+++yzkz7f4/HA5XKd8KLw6D96mY0DicOtu8antDWtXb1o7+b9VnRmqgaRpqYm+Hw+ZGZmnvDrmZmZqKurO+nzV6xYgcTExNArLy9PzfKiCkdESDRvnx+1bd0AuHXXyGLtFqTF2QGwvaHB0dSumeXLl6O9vT30qq6uFl2SYSgNf217D3p6fYKroWh0tLULfhmIsZqRHm8XXQ6pqJAjsDQEqgaRtLQ0mM1m1NfXn/Dr9fX1yMrKOunz7XY7EhISTnhReKTE2hAXvHK9mieskgDKyb4FqU5IkiS4GlJTaASWbQ0NgqpBxGazYfr06VizZk3o1/x+P9asWYM5c+ao+Wj6AkmSuIWXhDrSFOgd53PHjOGF1qQ1cUSEzsyi9gOWLl2KxYsXY8aMGZg1axYee+wxdHZ24pZbblH70fQFBalO7Kl1cbiUhFBGRHiYmfGx00NDoXoQuf7669HY2Ij7778fdXV1mDp1KlatWnXSAlZSX/+2OjYOFHnKzx1HRIwvdL9VCzs9dGaqBxEAuOOOO3DHHXdE4lF0GoW8hZcEUs6VKOTWXcNTRkTqXR50e32IsZkFV0RapqldM6Su/BQeNERi+PwyqlsCW3d5mJnxJTltSIyxAuCCVTozBpEoolxGVdPajV6fX3A1FE3qXD3w+vywmiVkJzpEl0MRwEMUabAYRKJIZrwDNosJfX45dLAUUSQoO2bykp2wmNnsRAMe9U6DxRYhiphMUuhCKi5YpUhSLkDL57RM1GBbQ4PFIBJl+rfVsZdCkaPsnuBC1ejBLbw0WAwiUaZ/wSobB4qcI03cuhttlPNiuIWXzoRBJMooC1a5hZciqf8wMwaRaKFMzdS0dsPbx8XxdGoMIlEmP4VTMxRZsiyHft6UETkyvvR4O2KsZvhloIaL4+k0GESiTOFxl1H5/bLgaigaNHV40eX1QZKAvJQY0eVQhBx/vxW38NLpMIhEmVHJMTCbJPT0+tHg9oguh6KAMhqSkxgDu4UnbEYTJYhUcU0anQaDSJSxmk0YlRTolXJ6hiJBWRjNE1Wjj3KWCEdE6HQYRKIQt9VRJCmBl0Ek+nBEhAaDQSQKhYIIt9VRBCg7Zgp4hkjUKUjhiAidGYNIFCrgWSIUQcqpqgU8QyTqKJ2e6pZu+Lg4nk6BQSQK5XNqhiKof2qGIyLRJicpBlazBK/PjzpXj+hySKMYRKJQ4XELyGSZvRRST3tXL9q6egFwjUg0Mpsk5CUHOz5NnJ6hgTGIRCHlUDN3T1/oTYJIDco6pLQ4O2LtFsHVkAihEVie5kynwCAShWJsZmQm2AGwcSB1KdN/hRwNiVqF3MJLZ8AgEqX6F6yycSD1hI52ZxCJWtzCS2fCIBKleJYIRUJlaESEC1WjVf8x72xraGAMIlGKd0BQJFTxVNWop+yWquLieDoFBpEo1d84sJdC6qnk1t2ol5scA0kCOr0+NHV4RZdDGsQgEqUKuJKdVNbl7QtdrMjFqtHLbjEjJzFwv1UVT3OmATCIRCllsWqj24NOT5/gasiIqoIhNzHGiiSnTXA1JFJoKriJHR86GYNIlEp0WpHktALof8MgCiflTYfrQ0iZmuMuPRoIg0gUU+7+YONAalCG4bk+hDgVTKfDIBLF+nspbBwo/Cp5mBkFFXILL50Gg0gU4/5+UlPoMDPeuhv18lP6t/ASfRGDSBQLbeHlSnZSQeh49zROzUQ7pdPT2tWL9m7eb0UnYhCJYjxdldTi6fOhtq0bABerEhBrtyAtLnC/Fc8uoi9iEIliymLV2rZuePv8gqshIzna2g2/DDhtZqQH34AouhXyNGc6BdWCyEMPPYS5c+fC6XQiKSlJrcfQCKTH2xFjNcMvA0db2Uuh8FF6vfkpTkiSJLga0gLl4kMeF0BfpFoQ8Xq9uPbaa/G9731PrUfQCEmSxOkZUoXS6+Vld6RQfhYqmzgiQieyqPWFH3zwQQDACy+8oNYjKAwKUp3YX+fmWSIUVkd42R19Ac8SoVNRLYgMh8fjgcfjCX3scrkEVhMdlJ0z3MJL4XSEl93RF/B0VToVTS1WXbFiBRITE0OvvLw80SUZXgHnbUkFR3iYGX2B8rNQ7/Kg2+sTXA1pyZCCyLJlyyBJ0mlf+/fvH3Yxy5cvR3t7e+hVXV097K9Fg6NcfsdeCoWLzy+jOrj4OZ9BhIKSnDYkOAKD8Oz40PGGNDVzzz334Oabbz7t5xQXFw+7GLvdDrudW/0iSRkRqW7phs8vw2ziDgcamdq2bvT6ZNjMJmQHr38nAgKH231+tB2VzZ0ozYoXXQ5pxJCCSHp6OtLT09WqhQTISYqB1SzB6/OjztWDUUl846CRUaZl8lJiGGzpBPkpTnx+tJ2HmtEJVFusWlVVhZaWFlRVVcHn82HHjh0AgNGjRyMuLk6tx9IQmU0ScpOdqGjqxJGmTgYRGrEjvHWXTiG0hZdTwXQc1YLI/fffjxdffDH08bRp0wAAH3zwAebPn6/WY2kYClKDQaSlC3NFF0O6x627dCo81IwGotqumRdeeAGyLJ/0YgjRHuWod/ZSKByUA6t4mBl9EUdEaCCa2r5LYoRu4eW8LYWB0tvljhn6ImWUrKaV91tRPwYR4jHvFDayLPN4dzqljHg7HFYT/DJQE7ydmYhBhE448VCWZcHVkJ41uD3o6fXDbJK48JlOIkkSzy6ikzCIEPJSYiBJQKfXh+ZOr+hySMeUUbWcJAdsFjYvdDKOwNIXsaUg2C1m5AQPnmIvhUaC0zJ0JoVpyogIgwgFMIgQgMBBQwAbBxqZ/svuuFCVBtbf1rDTQwEMIgSg/42Dt/DSSITOEEnhiAgNjFt46YsYRAjA8Vt42TjQ8PEwMzqTL95vRcQgQgA4IkIjd8LW3TSOiNDAshMdJ9xvRcQgQgD6h0s5b0vD1dbVC3dPHwAgL5kjIjQwi9mE3GSuE6F+DCIEAChMCzQMrV29aOUWXhqGw8Gj3bMTHYixmQVXQ1rGLbx0PAYRAgA4bRZkJTgAABXspdAwVASDSBGnZegMCrhLj47DIEIhyhuIcmkZ0VBUMojQIBVwKpiOwyBCIUXpgcahgkGEhoEjIjRYylQw2xoCGEToOMXBN5DDbBxoGJSfm+J0BhE6vaK0OACBqRk/t/BGPQYRClF2zlQ0MojQ0Pj98nFTM3GCqyGty02OgcUkobvXxy28xCBC/Y6fmuEtvDQU9e4edPf6YDFJyE3mrbt0elazKXTUO6dniEGEQvKSnTAHeyn1Lo/ockhHlFG0/BQnrGY2K3RmRZwKpiC2GBRis5iQF+zNHm7qEFwN6YnyZsITVWmwlLVEnAomBhE6Qf8WXu7vp8HjjhkaKmUtETs9xCBCJ1B6tBVsHGgIeIYIDVVRGo8LoAAGETpBMRsHGgbl56WYQYQGqSQ4NVPd0gVvn19wNSQSgwidoH+4lEGEBqfX50dVS2Aqr4hniNAgpcfbEWszwy8DVS1sb6IZgwidQHkjqWruQp+PvRQ6s6Ot3ejzy4ixmpEZ7xBdDumEJEmh9uYwF6xGNQYROkF2ggN2iwl9fhlHW7tFl0M6oKwnKkyLhckkCa6G9KQ4OALLqeDoxiBCJzCZpP4TVnkhFQ2C0pstCt4fQjRYXLBKAIMIDSDUOHC4lAaBW3dpuIo5NUNgEKEB8BZeGorKZt4xQ8PD01UJYBChAXC4lIaiopEjIjQ8ys9MU4cHrp5ewdWQKAwidBKeJUKD1e31obY9cHsqzxChoYp3WJEebwfAqeBoxiBCJ1FOV61p60ZPr09wNaRlyrRMktOK5Fib4GpIjzgCS6oFkcrKStx2220oKipCTEwMSkpK8MADD8Dr9ar1SAqT1Fgb4h0WAMCRZt45Q6fGhao0UsoJq1wnEr0san3h/fv3w+/34+mnn8bo0aOxe/du3H777ejs7MSjjz6q1mMpDCRJQnFaLHYebcfhxg6UZsWLLok0ikGERiq0YLWR91tFK9WCyCWXXIJLLrkk9HFxcTHKysqwcuVKBhEdKE6Pw86j7TjExoFOI3SGSCqDCA1PEQ81i3qqBZGBtLe3IyUl5ZS/7/F44PF4Qh+7XK5IlEUDGJ0RaBwOcQEZnUZ5MKiWZHDrLg1P8XHHBciyDEni6bzRJmKLVcvLy/GHP/wB3/nOd075OStWrEBiYmLolZeXF6ny6AtK0gNvLOUNHBGhgcmyjEPBn4/RDCI0THnJTphNErq8PjS4PWf+A2Q4Qw4iy5YtgyRJp33t37//hD9TU1ODSy65BNdeey1uv/32U37t5cuXo729PfSqrq4e+n8RhcXojEAv5VBjB/x+WXA1pEX1Lg86PH0wH3ctANFQ2Swm5CXHAACngqPUkKdm7rnnHtx8882n/Zzi4uLQ/6+trcWCBQswd+5cPPPMM6f9c3a7HXa7faglkQoKUmNhCfZSjrl6MCopRnRJpDEHG9wAgIJUJ2wWngRAw1eSHofK5i4cauzE3JI00eVQhA05iKSnpyM9PX1Qn1tTU4MFCxZg+vTp+NOf/gSTiY2VXljNJhSkOnGosRPlDR0MInQSZdpudDqnZWhkRmfEYc3+BpTXu0WXQgKolgxqamowf/585Ofn49FHH0VjYyPq6upQV1en1iMpzJR5f64ToYGUc30IhYnyM3SQbU1UUm3XzOrVq1FeXo7y8nLk5uae8HuyzDUHejA6Iw7v7alnEKEBMYhQuIzJDJxVxCASnVQbEbn55pshy/KAL9KH/i28bBzoZMrPBYMIjZTyM9To9qC9i5ffRRsu2qBTGp0e6KUcYi+FvqCty4umjsB1DSVcI0IjFGe3IDvRAQAob+Q6kWjDIEKnVBLcwtvc6UVrJ+8Ion7KtExOogOx9oiei0gGFVonUs+OT7RhEKFTctosod0y5ZyeoeMoQYQnqlK4jMkIjMByTVr0YRCh0yrhzhkaABeqUrhx50z0YhCh01Ku6GYQoeOVc6EqhdmYTHZ6ohWDCJ0WzxKhgfAwMwo35Weppq0bnZ4+wdVQJDGI0GkpjQO38JKi2+tDTVs3AI6IUPgkx9qQFmcDwPYm2jCI0GkpbzQ1bd3o9voEV0NacKixA7IMJDutSI3j3VAUPtw5E50YROi0UuPsSHZaIcvspVAADzIjtSg7Z7hgNbowiNAZjQ0ev7y/jgcNEXAgeDEZgwiFG9ekRScGETqjcVmBIFJW5xJcCWlBWV3gTaI0GFCJwmVMKIiw0xNNGETojEqzEgBwRIQCyuoDgVT5uSAKl9HBLbxVLV3o6eWatGjBIEJnVBoaEWEQiXYdnj5UtwR2zCg/F0Thkh5nR2KMFX6Z0zPRhEGEzkh5w2lwe3jnTJQ7GFwfkh5vR0qsTXA1ZDSSJLHjE4UYROiM4uwW5KUE7pzh9Ex0U94cxnE0hFQyXgki9WxrogWDCA1KaWZgPQAXrEY35c1hLBeqkkqUtUf7jrGtiRYMIjQo49hLIfSPiHB9CKllXDanZqINgwgNivLGw6mZ6MapGVKbsi28we1BC9ekRQWL6AJIH8Ydt4DM75dhMkmCK4oMv1/GhopmfFLehIP1Heju9SHBYcXYzHgsHJ+Bs3ISIEnR8b1o6vCgudMLSeo/AZMo3GLtFuSnOFHV0oX9dS7MLUkTXVLEVDR1YvXeOuypdaGtqxcWk4TCtFjMKU7FuWPTYLeYRZeoCgYRGpTCtFjYzCZ0eX042tqN/FSn6JJU5fPL+Nu2o3h8zUEcbe0+6ff/tesY/vs/BzAhOwH3XDQWF4zLMHwgUUZDClKciLEZs0EkbRiXFR8IIsfcURFEdlS34bfv7ccn5c0D/v5zH1cg2WnFt88rwc1zCw33749BhAbFajahJCMO+465sL/OZeggUtXchbtf34GtR1oBAPEOCy6akIVJoxIQ77CipdOLbVWtWLu/AXuPuXDbi1tw2aQs/OrqSYbe0rqf60MoQsZlxeP9vfWGXyfS0+vDI6v240+fVAIATBIwb3QazilORWaCA54+H/Ydc2H13nrUuzx4ZNV+/HVLNX533RRMy08WW3wYMYjQoI3Pise+Yy6U1blx0VlZostRxaflTfjOn7fC3dOHOLsFP1w4GjfNKYTDenIPpLXTi6fWH8JzH1Xg37vq8PnRdvzp5pkYY9AdJfuDuxh4tDupbVy2cpqzcXfONLo9+NaLm7HzaDsA4Jppo7D0orHITT65k/eLK8/CP3bW4pFV+3G4qRPXPf0ZHvryJFw3Iy/SZauCi1Vp0EILVg26c2bV7jrc9PwmuHv6MC0/Ce/eeS6+fV7JgCEEAJJjbVh+6Xi8tWQeClOdONrajWtWfood1W2RLTxC9tQG3hQm5CQKroSMTmlrDtR3wOeXBVcTfrVt3bhm5SfYebQdyU4rnr95Bn5//dQBQwgAWMwmXHN2Lt6/63xcNikLvT4ZP3njc6z88FCEK1cHgwgNmpFPPFx/oBE//Mt29PllXD45G3+5/RzkpQxu+mniqET8/fvzMKMgGe6ePtz43EbsrmlXueLI8vb5cTB4EdlZObxjhtRVmBoLu8WE7l4fqlq6RJcTVk0dHnzzuY2obulGfooTf//+PFwwLnNQfzbRacUTXz8bSxaUAAAeWbUf//vRYTXLjQgGERq0ccGDhiqaOg11IdX2qlZ8+6Ut8Pr8uGxSFh7/2rRTjoKcSkqsDS/eOgszCwNh5LYXN6Pe1aNSxZF3oN6NXp+MBIcFuckxosshgzObpNBOvb21xpme6fD0YfHzm3C4sROjkmLw6rfPQVFa7JC+hskk4ccXj8Pdi8YCAH71r31YtbtOjXIjhkGEBi0zwY7UWBt8ftkw54k0uj347p+3oqfXj/PHpuOx66fBPMytybF2C56/eSbGZMSh3uXBt/9vi2EC295jyrRM9GxXJrHOGhWYAtxda4zRRVmWce8bn2NPrQtpcTa8dNss5CQNP9TfuWgMbppTAAC4+7Uduh6FZRChQZMkqb9x0PEPvaLP58cP/rIN9S4PStJj8eQNZ8NmGdk/iXiHFc8tnolkpxU7j7bjZ2/uDlO1Yim90rO4PoQiZGKOcdoaAHj+k0r8a9cxWM0SnrlpBorT40b8Ne+/YgLOHZOG7l4fvvPSVrR394ah0shjEKEhmTQqMD2zxwC9lMf+cxAbDrcg1mbG0zdOR5w9PJvI8lOd+OMN02GSgL9tO4p/fl4blq8rkvL3zfUhFCkTQ22NC7Ks7wWr26paseLf+wAAP798As4O09Zbi9mEJ75xNvJTnKhp68bP39qty+8VgwgNSX8vRd/ztjur2/DHD8sBAA9/ZTJGh/mk0DklqViyYDQA4Kd/34XatpMPRdMLv1/miAhF3NjMeFhMElo6vTjWrt/1Vj29Pvzo9Z3o88u4YnJ2aDolXBJjrHjsa1NhNkl4Z2ct3tpRE9avHwkMIjQkE4NTM2V1bnj7/IKrGZ6eXh/u+etO+GXgqik5uHJKjirP+eHCMZiSlwRXTx9+8sbnuuypAEBVSxc6vT7YLCaUpA9tYR3RcDms5tCZPHqennn0vTIcbupEZoIdD109SZU1VmfnJ+POhWMAAPe/tQd1OgtuDCI0JLnJMUiMscLr8+OATs8TeXzNQZQ3dCAtzo4HrzpLtedYzSY8dv1U2C0mfFzepMueCtB/fsi4rHhYzGwyKHImBqcCd+t058zWI6147pMKAMCKayYh0WlV7Vnfn1+CqXlJcHv68P/+uUe156hB1VblqquuQn5+PhwOB7Kzs3HjjTeitlb/8+XRTJKk4+Zu9ddLOdTYgWeD++5//eWJSFb5SPaitFj8MNhT+dU/96GtS3+3iXJ9CImijMDu0eGIiM8vB9dsANecPWrQZ4UMl8Vswq+/PAlmk4R/76rDmn31qj4vnFQNIgsWLMDrr7+OsrIy/O1vf8OhQ4fw1a9+Vc1HUgTodZ2ILMt48J296PXJWDguI2LH1N9+bjHGZMShudOLR1btj8gzw2lX8E2AJ6pSpCmdHj1u4f3LpirsO+ZCgsOCn102PiLPnJCTgG99qQgAcP/be9Dt1cfxAaoGkbvvvhvnnHMOCgoKMHfuXCxbtgwbNmxAb68+txhRgLKFd5fOeimr99Zj/YFG2Mwm3HfFhIg912Yx4dfXTAIAvLq5WlcjSbIsY2fwyPqpuUlCa6HoMz47AZIE1Ls8aHDrZ91Da6cXj75fBgBYeuFYpMbZI/bsOxeNwaikGNS0dYdGf7UuYhO+LS0tePnllzF37lxYrQPPk3k8HrhcrhNepD3KvO2+Yy70+fSxYLWn14df/msvAOBb5xahcIinGY7UzMIUXDklB7IMPPSvfbpZuFrZ3AVXTx9sFhPGZfOyO4osp82CkuB5G3t0NAL7+9UH0NbVi9LMeHzznPDukjkTp82Cey8dBwB4at0hNOjghGfVg8i9996L2NhYpKamoqqqCm+//fYpP3fFihVITEwMvfLyjHGzoNEUpsYizm6Bp8+PA/UdossZlOc+rkB1SzeyEhyhbbWR9pOLS2GzmPDpoWas3d8gpIah2lHdCiAQPq1cqEoCTA6OwOrlMsn9dS68vPEIAOCBqyYIWeB95eRsTMtPQpfXh9+9fyDizx+qIX+Hli1bBkmSTvvav79/HvzHP/4xtm/fjvfffx9msxk33XTTKXuDy5cvR3t7e+hVXV09/P8yUo3JJGFKXqBx2B58o9Ky1k4vngreUrns0nGIDdPBZUOVl+LELfMKAQC//vc+XYwm7awOTCNNyUsSWwhFrWn5SQCA7ToJIr9ZVQa/DFw6MQtzS9KE1CBJEn5+eWBdyutbqzV/X8+QW+R77rkHN99882k/p7i4OPT/09LSkJaWhrFjx2L8+PHIy8vDhg0bMGfOnJP+nN1uh90eubk0Gr6z85PxSXkztle14YbZkR16HKqV6w7B7enD+OwEXKXSmSGDtWTBaLy+uRqHGjvxl83VuDHCw7ZDpTT+UxlESJBpwVNId1S1wu+XYRrmXVCRsKWyBWv3N8BskvDji0uF1jK9IAWXT8rGv3Ydw6//vQ9//tZsofWczpCDSHp6OtLT04f1ML8/0AP0eDzD+vOkHaFeSpW2R0SOtXfjhU8rAQSmRkQ3YgkOK+5aNBYP/GMPHl9zEF89OxcxtqHd9Bspnj4f9gV7UgwiJEppVjwcVhNcPX043NSJ0Rkjv6NFDbIs4zerAgtUr5uRG5a7ZEbq3kvG4f29dfi4vAmfHWrGnJJU0SUNSLXJq40bN+KJJ57Ajh07cOTIEaxduxZf//rXUVJSMuBoCOnL1LxAL+VQY6emz8Z4fM1BePv8mFWYgvmlwwvQ4fb1WfkYlRSDRrcnNJesRfuOueH1+ZHktCI/xSm6HIpSVrMJk0clAdB2x+fDA43YVNkCu8UUOjtItPxUJ742Mx8A8N//OaDZRfKqBRGn04m///3vWLhwIUpLS3Hbbbdh8uTJWLduHadfDCAl1oai4M4TrS4iO9zYgde3HAUA/OSSUs1cX2+zmPDDhYEFs0+tO4Qub5/gigambNudkpukme8dRSetrxPx+2X8NjgasnhuIbITYwRX1O/7C0pgs5iwqaIFn5Q3iy5nQKoFkUmTJmHt2rVobm5GT08PKioqsHLlSowaNUqtR1KETQsO12+vahNax6k8sbYcPn/g8LIZhSmiyznBNWfnIj/FiaYOL176TJujIju5PoQ0on8quE1oHafy/t467D3mQrzdgu+dXyK6nBNkJ8bgG7MCoyK/X12myVER7sejYZtWEJie2abB4dLKpk68vTNwncBdi8YKruZkVrMJP7igf1Skw6O9UZGtwb9XBhESTVmwWlbnQqfG/q3IsozH1wRu8r55XqHq10YMx/fnl8BuMWFbVRvWHWgUXc5JGERo2JQRkR3VbfD7tZWy//hhYDRkQWk6JuVq82jyL08bhcJUJ1q7evFicEGtVtS19+BIcxckCZhemCy6HIpymQkO5CQ64JeBz49q62TiNfsasPeYC7E2M26dVyS6nAFlJDhCO/T+e7X21oowiNCwjcuKR4zVDHdPHw41audgs+qWLvx9W+Cm2x9oZNHYQCxmE+5cFKjv2Y8Ow92jnasPNlW2AAAmZCcgwaHejaFEg6WMimhpBFaWZfxh7UEAwI1ztDkaovju/BLEWM3YebRdcwcqMojQsFnMptCw/caKFrHFHGflukPo88s4d0wazs7Xdm/+qimjUJIei7auXvyfhtaKbA7+fc4q0tbaGopeys/ihsPaWXC57kAjdh5th8NqwrfO1eZoiCItzo6b5gZGRR5fW66pUREGERqRc4oD+9I/00jjUNvWjb9uCZzI+4MLtDsaojCbpFCdz350WDNrRTYpQURji3wpeiltzZbKVnj7xJ9KHBgNCawN+ebsAqRF8GK74br93GI4rCbsrG7Dx+VNossJYRChEVEOyNl4uFkTCfvpdYfQ65NxTnGKbnrzV0zORlFaYFRECztoWju9KKt3AwBm6uR7SMY3JiMOKbE2dPf6sKumTXQ5+OxQM7YeaYXNYsK3zys+8x/QgLQ4O74e3EHzRDBEaQGDCI3IlLxE2C0mNHV4Ud4gdp1Ig6sHf9kcGA35oQ5GQxQWswl3BC/ie/ajw8LPFdkcXB9Skh6ri14eRQeTScLs0PSM+Kng/1kTWBvy9Zl5yEhwCK5m8L59XjGsZgkbK1pCI5+iMYjQiNgtZswI7qoQPXf79PrD8Pb5MaMgWbNHGZ/Kf03NQUGqEy2dXry8oUpoLaFpmSJ9fQ/J+JTpGdFtzaaKFmysaIHVLOE7Gjs35EyyE2Pw1emBm+2f+EAboyIMIjRi5xSJXyfS1NF/XPoPFo7R3UmgFrMJS4KjIk+vP4Rur09YLcqOmdmcliGN0co6EWWnzLUz8pCTpJ1TVAfre+eXwGySsP5AIz4/2ia6HAYRGjll9GHD4RZh60Se/egwenr9mJKXhPPGiLl6e6S+PG0UcpNj0NThxSubxIyKtHV5sbsmcE7D7GIGEdIWLawT2VbVio8ONsFikjR3iupg5ac68V9TAzeRa2GtCIMIjdjk3CTEWM1o6fTiQH3k14m0dvYfk/7DC0brbjREYT1uVOSpdYfQ0xv5UZFPDzXDLwcafC3dl0EEBNaJnBMMyJ8dEjMC+4fg2pAvTxuFPB1fBvn9+aMhScD7e+uxv84ltBYGERoxm8UUWify0cHIHx/8/CcV6PL6cFZOAi4YlxHx54fTV87ODd3M+6qAURHl7+/cMdq4qZjoi+aUBEY81x+I/PbTXUfb8UFZI0wSQp0GvRqdEYfLJmYDAP74wSGhtTCIUFicPzbwxhXpewzau3vxwieVAIAf6Hg0RGGzmPC9+YHh3pURHhWRZTnUuJ87Vp/TW2R884NtzdaqVrR3R/Y0YmVtyFVTclAYvH1cz5YsGA2LSYLNYhJ6TQeDCIXF/NLASMTGwy0RvZTqxU8r4fb0YWxmHC6akBWx56rp2hm5yE50oN7lwV+3Ho3Ycw83daKmrRs2s4kLVUmz8lKcKEmPhc8v45MIHsq175gL7++thyQBd1yg79EQxYScBHy2fCEevXYKTCZxnTgGEQqLkvRY5KXEwOvzR2zutsPTh+c/qQAQSPYi/yGFk91i7h8V+aAcnr7IjIqs3Re4f2JmUTKcNktEnkk0HErH58OyyN2Zomx1vWxiNkZnxEfsuWpLjxd/VhCDCIWFJElYEGwc1kaocXjpsyNo6+pFcVosrpicE5FnRsp1M/KQEW9HbXsP/ra1JiLPXL2vHgBw4fjMiDyPaLjmlwamZz4sa4zIlEJ5gxv/3nUMgHFGQ7SEQYTCRgkia/bVq944dHn78L8fHQYAfH/BaJgNMhqicFjN+G5wa+CTH5Sj16fumQktnV5sCZ4fsmgCgwhp26yiFMTZLWhwe7AjAudgPLG2HLIMXDQhE+OzE1R/XrRhEKGwmTs6FfF2C+pdHmyvVveq7lc2VqG504u8lJjQfnij+cbsfKTF2VHT1o2/b1N3rcja/Q3wy8CE7ATkJut3SyJFB7vFjAXBHXKrdtep+qyKpk78Y2ctAH1cpKlHDCIUNnaLGQvHBxqHd3ep1zh0ePrwxw8D282WzB8Nq9mYP8aBUZHAZVpPqDwq8t6ewN/XhRwNIZ24dGJgcfq7u4+pepDif68+AL8MLChNx6TcRNWeE82M2YKTMJcE96W/u7tOtcbhuY8q0NLpRVFaLL46PVeVZ2jFDbMLkBZnQ3VLN97eUavKM9q7erGuLLDt+tJJxth5RMY3vzQdDqsJ1S3d2FOrzoFce2rbQ6Mh91xUqsoziEGEwmx+aTqcNjNq2rqxvbot7F+/pdOLZ4NrQ5ZeOBYWg46GKGJsZtx+bnBUZO1B9KkwKrJqzzF4fX6My4rHuCzOf5M+OG0WzB8bGIF9Z6c6If3R98oAAFdMzsbEURwNUYuxW3GKOIfVjEvOCvSq/6bCGRhPrTuEDk8fJmQn4PJJ2WH/+lr0zXMKkBJrQ2VzF97cHv4dNMpIy1UGXWtDxnX1tFEAgDe314Q9pG+ubMEHZY0wmySOhqiMQYTCTpku+cfO2rCeDFrX3oMXP60EAPz4klLDnBtyJrF2C75zXmBU5HfvH0CXN3wHxtW2dYduTb7SYFugyfguGJeBZKcVDW4PPg7j4WayLOM3q/YDCGylLzLAKapaxiBCYXdOcSpGJcXA3dOH9/fWh+3rPrJqPzx9fswqTAkd8xwtFs8tRG5yDOpcPXh2fUXYvu5fNlVBloFzilN0fYEXRSebxYSrpgQC9BthHIF9b08dNle2wm4x4c6F3CmjNgYRCjuTScJXzg4Mmf55w5GwfM0tlS14c3sNJAn4+RXjdX+nzFA5rGYsu3QcgMD0VL2rZ8Rfs9fnx6ubqwEEpn+I9OjaGXkAAuEhHP8uur0+/PKf+wAA3z6vGFmJjhF/TTo9BhFSxTdmF8BikrCpogW7jraP6Gv5/DJ+8c4eAMD1M/IwOTcpDBXqz+WTsnF2fhK6e3146F/7Rvz13ttTh0a3B2lxdsPc00PRZ+KoRMwsTEavTw5Lx+fp9YdQ09aNnEQHvj+fp6hGAoMIqSIr0YErJgcWkz738eERfa1XNh7B7hoX4h0W/Oji6F00JkkSHrxqIkxSYP3NByM4Sl+WZTy1LnAWyzdm5cFmYVNA+nXrvCIAwMsbq0a0Lu1IcydWBs8o+tnlExBjM4elPjo9tj6kmtu+FFhg+c7nx1DR1Dmsr1HV3IUV7wYWjd1z4VikxYm/oEmkSbmJoUb352/uHvZNxx+UNWB3jQtOmxk3B78ekV5dOCETuckxaOn0DntUxOeX8eO/fg5Pnx9zS1JxGc/UiRgGEVLNpNxELChNh88v49H3y4b85/1+GT96Yye6vD7MLkrBTXMKw1+kDi29aCxyk2NQ09aNh4MhbSj8fhn/85+DAIAbg1uDifTMYjbhB8HL6J78oBzunt4hf43nP67ApsoWxNrMeOQrk6NuHZpIDCKkqp9cMg6SBPzr82PYXjW0+2ee/egwNlW0wGkz47dfnRI123XPxGmzYMU1kwAAL204EjqefbDe2HoUO4+2w2kz41vBw9KI9O4rZ+eiOD0WrV29oemVwdpf58Jvg52l+66YwB1kEcYgQqoan52Aa6YFzhX5yRufD3r+9uODTXgkuI//55dPQH4qG4bjnTsmHd8Oni3y47/uRHlDx6D+XGunFw8Hv693LxqL9Pjonuoi47CYTbj3ksDOsqfXH8bumsEtkm/t9OL2/9sCb58fC0rTcf3MPDXLpAFEJIh4PB5MnToVkiRhx44dkXgkacjPLx+PtDg7DjZ04DerzjxFs6e2Hd9/eSv8cuBwtK/PYsMwkB9dVIrpBclw9fThlhc2ocF9+q2LPr+MH766HS2dXozNjMPN8wojUyhRhFx8VhYun5QNn1/G0td3oOMMa6i6vH34zktbUd3SjbyUGPz+uqmckhEgIkHkJz/5CXJyeGpjtEqOtYWmEp7/pAIvfHLqA7l2Vrfhxuc2wdXTh+kFyfjV1RPZMJyCzWLCMzdOR36KE9Ut3bj+6Q2obuka8HP9fhm//OdefHSwCQ6rCf/ztWmGvbWYotv/+6+zkBZnw4H6Dnzvz1vh7Rv46Pf2rl7c+sJmbKpsQbzdgmdvmoFkrpcSQvWW6N1338X777+PRx99VO1HkYZdOCETP7poLADgF+/sxa//vQ+evv5pGm+fH899XIFrn/4MLZ1eTM5NxJ9umQmHldvnTic1zo6XbpuFUUkxqGjqxJVPfIy3d9TA7++/+bil04u7XtuBF4LH4z/ylckYn83L7ciYUuPs+N/FMxFjNeOjg034xrMbUNPWfcLnfHqoCVc+8TE2HG5BnN2CF2+bxQsfBZJkte5qB1BfX4/p06fjrbfeQlpaGoqKirB9+3ZMnTp1wM/3eDzweDyhj10uF/Ly8tDe3o6EBP6Q6J0sy3hkVVno/IqUWBu+NDoNALDhcDMa3IG/+0XjM/H766cgwWEVVqve1LX34LYXN4euQy9IdWJ6fjJ8soy1+xrg9vTBbJLwyFcmh+4CIjKy9QcaseTlbXB7+mAxSfjSmDSkxdmxt9aFvccC/05yk2Pw1Den82ZdFbhcLiQmJg7q/Vu1ICLLMi677DLMmzcPP//5z1FZWXnGIPKLX/wCDz744Em/ziBiLO/tqcMv/rEHx9pPXNOQEW/HnYvG4Osz87lDZhi8fX788cNyPPdRBdxfmBs/KycB918xAbOLUwVVRxR5R5o78ZM3PsfGipYTft1qlvD1WflYeuFYJDk5HaMGVYPIsmXL8Mgjj5z2c/bt24f3338fr7/+OtatWwez2TyoIMIRkejR5/Pjo/ImlNd3wOvzY3JuImYWpnAqJgy6vH34sKwRR5q74JdlTM1LwpziVIY7ilr7jrmwubIFbV29KEyLxTnFKciI5x0yalI1iDQ2NqK5ufm0n1NcXIzrrrsO77zzzgkLDX0+H8xmM2644Qa8+OKLZ3zWUP5DiIiISBs0MTVTVVUFl8sV+ri2thYXX3wx3njjDcyePRu5uWeep2YQISIi0p+hvH9b1CoiPz//hI/j4uIAACUlJYMKIURERGR8PEiAiIiIhFFtROSLCgsLoeJOYSIiItIhjogQERGRMAwiREREJAyDCBEREQnDIEJERETCMIgQERGRMAwiREREJAyDCBEREQnDIEJERETCMIgQERGRMBE7WXU4lJNYj788j4iIiLRNed8ezInqmg4ibrcbAJCXlye4EiIiIhoqt9uNxMTE036OJGv4Ahi/34/a2lrEx8dDkqSwfm2Xy4W8vDxUV1ef8YpiGj5+nyOD3+fI4Pc5Mvh9jhy1vteyLMPtdiMnJwcm0+lXgWh6RMRkMiE3N1fVZyQkJPAHPQL4fY4Mfp8jg9/nyOD3OXLU+F6faSREwcWqREREJAyDCBEREQkTtUHEbrfjgQcegN1uF12KofH7HBn8PkcGv8+Rwe9z5Gjhe63pxapERERkbFE7IkJERETiMYgQERGRMAwiREREJAyDCBEREQkTlUHkySefRGFhIRwOB2bPno1NmzaJLslwVqxYgZkzZyI+Ph4ZGRm4+uqrUVZWJrosQ3v44YchSRLuuusu0aUYUk1NDb75zW8iNTUVMTExmDRpErZs2SK6LEPx+Xy47777UFRUhJiYGJSUlOCXv/zloO4roVNbv349rrzySuTk5ECSJLz11lsn/L4sy7j//vuRnZ2NmJgYLFq0CAcPHoxYfVEXRF577TUsXboUDzzwALZt24YpU6bg4osvRkNDg+jSDGXdunVYsmQJNmzYgNWrV6O3txcXXXQROjs7RZdmSJs3b8bTTz+NyZMniy7FkFpbWzFv3jxYrVa8++672Lt3L373u98hOTlZdGmG8sgjj2DlypV44oknsG/fPjzyyCP4zW9+gz/84Q+iS9O1zs5OTJkyBU8++eSAv/+b3/wGjz/+OJ566ils3LgRsbGxuPjii9HT0xOZAuUoM2vWLHnJkiWhj30+n5yTkyOvWLFCYFXG19DQIAOQ161bJ7oUw3G73fKYMWPk1atXy+eff7585513ii7JcO699175S1/6kugyDO/yyy+Xb7311hN+7ZprrpFvuOEGQRUZDwD5zTffDH3s9/vlrKws+be//W3o19ra2mS73S7/5S9/iUhNUTUi4vV6sXXrVixatCj0ayaTCYsWLcJnn30msDLja29vBwCkpKQIrsR4lixZgssvv/yEn2sKr3/84x+YMWMGrr32WmRkZGDatGl49tlnRZdlOHPnzsWaNWtw4MABAMDOnTvx8ccf49JLLxVcmXFVVFSgrq7uhPYjMTERs2fPjtj7oqYvvQu3pqYm+Hw+ZGZmnvDrmZmZ2L9/v6CqjM/v9+Ouu+7CvHnzMHHiRNHlGMqrr76Kbdu2YfPmzaJLMbTDhw9j5cqVWLp0KX76059i8+bN+OEPfwibzYbFixeLLs8wli1bBpfLhXHjxsFsNsPn8+Ghhx7CDTfcILo0w6qrqwOAAd8Xld9TW1QFERJjyZIl2L17Nz7++GPRpRhKdXU17rzzTqxevRoOh0N0OYbm9/sxY8YM/PrXvwYATJs2Dbt378ZTTz3FIBJGr7/+Ol5++WW88sorOOuss7Bjxw7cddddyMnJ4ffZwKJqaiYtLQ1msxn19fUn/Hp9fT2ysrIEVWVsd9xxB/75z3/igw8+QG5uruhyDGXr1q1oaGjA2WefDYvFAovFgnXr1uHxxx+HxWKBz+cTXaJhZGdnY8KECSf82vjx41FVVSWoImP68Y9/jGXLluFrX/saJk2ahBtvvBF33303VqxYIbo0w1Le+0S+L0ZVELHZbJg+fTrWrFkT+jW/3481a9Zgzpw5AiszHlmWcccdd+DNN9/E2rVrUVRUJLokw1m4cCF27dqFHTt2hF4zZszADTfcgB07dsBsNosu0TDmzZt30vbzAwcOoKCgQFBFxtTV1QWT6cS3JbPZDL/fL6gi4ysqKkJWVtYJ74sulwsbN26M2Pti1E3NLF26FIsXL8aMGTMwa9YsPPbYY+js7MQtt9wiujRDWbJkCV555RW8/fbbiI+PD801JiYmIiYmRnB1xhAfH3/SmpvY2FikpqZyLU6Y3X333Zg7dy5+/etf47rrrsOmTZvwzDPP4JlnnhFdmqFceeWVeOihh5Cfn4+zzjoL27dvx+9//3vceuutokvTtY6ODpSXl4c+rqiowI4dO5CSkoL8/Hzcdddd+NWvfoUxY8agqKgI9913H3JycnD11VdHpsCI7M3RmD/84Q9yfn6+bLPZ5FmzZskbNmwQXZLhABjw9ac//Ul0aYbG7bvqeeedd+SJEyfKdrtdHjdunPzMM8+ILslwXC6XfOedd8r5+fmyw+GQi4uL5Z/97Geyx+MRXZquffDBBwO2x4sXL5ZlObCF97777pMzMzNlu90uL1y4UC4rK4tYfZIs88g6IiIiEiOq1ogQERGRtjCIEBERkTAMIkRERCQMgwgREREJwyBCREREwjCIEBERkTAMIkRERCQMgwgREREJwyBCREREwjCIEBERkTAMIkRERCQMgwgREREJ8/8BLEDUeKyUV3wAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "H = np.sin*X(1) + Y(2)\n", + "rho0 = Y(1) + X(2)\n", + "tlist = mp.timegrid(0, 10, 0.5**6)\n", + "\n", + "qsys = mp.System(H, rho0, tlist)\n", + "qsys.evolve()\n", + "\n", + "y_component = mp.frobenius(qsys.states, Y(1)(2)) # Y in first position of two-qubit system.\n", + "plt.plot(tlist, y_component)" + ] } ], "metadata": {