{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[](https://mybinder.org/v2/gh/PyLops/pylops_transform2022/HEAD?labpath=PyLops_v2.ipynb)\n",
"\n",
"[](https://colab.research.google.com/github/PyLops/pylops_transform2022/blob/main/PyLops_v2.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "fvzDiWgxrEbB"
},
"source": [
"# Building Operators in PyLops V2: What's New? \n",
"\n",
"**Author: C. Costa, PyLops Dev Team**"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "JXwvtkqirEbE"
},
"source": [
"The aim of this second tutorial is to introduce you to some features in our latest creation (still in the making): **PyLops V2**\n",
"\n",
"As any new Major version, there will be breaking changes! \n",
"\n",
"Nevetheless, this is the very first time since the initial creation of PyLops that we decided to sit back and revisit some the early design choices. We hope the benefit introduced by such changes outweights the pain of making a few changes to your codes.\n",
"\n",
"In the following we want to give you a taste of some new features and highlight along the way if (and how) these have required breaking changes.\n",
"\n",
"For more details on how to quickly migrate your PyLops v1.x codes into PyLops v2.x codes please refer to this [document](https://github.com/PyLops/pylops/blob/dev/MIGRATION_V1_V2.md)."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "hvtyDt5KrEbF"
},
"source": [
"## Useful links\n",
"\n",
"- Tutorial Github repository: https://github.com/PyLops/pylops_transform2022\n",
" \n",
"- PyLops Github repository: https://github.com/PyLops/pylops\n",
"\n",
"- PyLops reference documentation: https://pylops.readthedocs.io/en/latest/"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "d4BKiX3mrEbG"
},
"source": [
"## Let's Get Shifty\n",
"\n",
"To get started, let's suppose we want to implement a simple shift operator.\n",
"Stating the problem more descriptively, we want an operator that\n",
"\n",
"1. Depends on a (possibly fractional) parameter `shift`\n",
"2. Can be applied to a vector along a certain given `axis`.\n",
"3. A positive `shift` should shift the vector to the right; a negative `shift` to the left.\n",
"4. The shift is not circular, that is, when shifting to the right, the last samples will not appear in the beginning of the array.\n",
"\n",
"PyLops already features a much fancier [Shift operator](https://pylops.readthedocs.io/en/stable/api/generated/pylops.signalprocessing.Shift.html) that works in the Fourier domain, but it is circular.\n",
"Numpy offers the `roll` function, but on its own, it does not support fractional shifts and is also circular.\n",
"Note that both of these options could be made non-circular by padding, but for now we're gonna \"roll\" our own! 😉\n",
"\n",
"\n",
"---------------------\n",
"\n",
"Before we jump on creating the operator, let's code a function to do the shift.\n",
"After we have the workhorse, it will be easier to write the boilerplate code required for the `LinearOperator`."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uumt2jRqrI_L",
"outputId": "141a4b16-64cf-45f7-fad0-93260ceb87e4"
},
"outputs": [],
"source": [
"!pip install --upgrade git+https://github.com/PyLops/pylops.git@dev numba==0.55.1"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"id": "rVj---WSrEbG"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import logging\n",
"import matplotlib.pyplot as plt\n",
"import numba\n",
"import numpy as np\n",
"import numpy.typing as npt\n",
"import pylops\n",
"\n",
"from typing import Tuple\n",
"\n",
"logger = logging.getLogger()\n",
"logger.setLevel(logging.INFO)\n",
"try:\n",
" import cupy as cp # Remove if no GPU\n",
"except ImportError:\n",
" logger.warning(\"Could not import CuPy\")\n",
"\n",
"plt.style.use('tableau-colorblind10')\n",
"np.random.seed(0)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "x6gkOnozv0Pj",
"outputId": "4fd3bb10-4164-4dee-edc5-1f805d465085"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"numba : 0.55.1\n",
"google : 2.0.3\n",
"matplotlib: 3.5.1\n",
"numpy : 1.21.5\n",
"IPython : 5.5.0\n",
"cupy : 9.4.0\n",
"pylops : 1.13.1.dev405+g7eda999\n",
"\n",
"Watermark: 2.3.0\n",
"\n"
]
}
],
"source": [
"#!pip install -q --upgrade watermark\n",
"#%load_ext watermark\n",
"#%watermark --iversions -w"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "K8BNlzHcrEbI"
},
"source": [
"We're going to start with an integer-based shifting function, which we will exploit later to convert into a fractional shifter. We want this function to be pretty fast, so we're going to pad the original array and pad it with the correct number of zeros.\n",
"\n",
"\n",
"\n",
"![](\"figs/crop_pad_concatenate.png\")
\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "EKKC25SFrEbI",
"outputId": "85a6f5c6-da35-4d41-ec16-dc025a43ebf5"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Original array: [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
"Positive shift: [0. 0. 0. 1. 2. 3. 4. 5. 6.]\n",
"Negative shift: [4. 5. 6. 7. 8. 9. 0. 0. 0.]\n",
"No shift : [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
"Big pos. shift: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
"Big neg. shift: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
]
}
],
"source": [
"# Original array\n",
"x = 1.0 + np.arange(9)\n",
"\n",
"print(\"Original array:\", x)\n",
"\n",
"# Tests for our (integer) shifting method\n",
"shift = 3\n",
"shifted = np.concatenate((np.zeros(shift), x[:-shift]))\n",
"print(\"Positive shift:\", shifted)\n",
"\n",
"shift = -3\n",
"shifted = np.concatenate((x[-shift:], np.zeros(-shift)))\n",
"print(\"Negative shift:\", shifted)\n",
"\n",
"shift = 0\n",
"shifted = np.concatenate((x[-shift:], np.zeros(-shift)))\n",
"print(\"No shift :\", shifted)\n",
"\n",
"shift = 12\n",
"shifted = np.concatenate((np.zeros(shift), x[:-shift])) # Oops!!\n",
"print(\"Big pos. shift:\", shifted)\n",
"\n",
"shift = -12\n",
"shifted = np.concatenate((x[-shift:], np.zeros(-shift))) # Oops!!\n",
"print(\"Big neg. shift:\", shifted)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Np6MOMTurEbK"
},
"source": [
"Now that we know what to do, we can package it into a simple 1D function."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "lthU4sVprEbK",
"outputId": "485e2bf3-9b3f-4f18-c949-29165bde23ef"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0. 0. 0. 1. 2. 3. 4. 5. 6.]\n",
"[4. 5. 6. 7. 8. 9. 0. 0. 0.]\n",
"[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
"[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
"[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
]
}
],
"source": [
"def shift1d(x: npt.NDArray, shift: int):\n",
" shift = max(\n",
" min(shift, len(x)), -len(x)\n",
" ) # Fix the oops: ensure shift between -len(x) and len(x)]\n",
" if shift > 0:\n",
" return np.concatenate((np.zeros(shift), x[:-shift]))\n",
" else:\n",
" return np.concatenate((x[-shift:], np.zeros(-shift)))\n",
"\n",
"\n",
"for shift in [3, -3, 0, 12, -12]: # This should be probably be a unit test...\n",
" print(shift1d(x, shift))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "A08Ip7VIrEbL"
},
"source": [
"Looking good! Now we need to make sure this works for N-dimensional arrays.\n",
"The trick here is to rely on swapaxes (it returns a view, cheap!) and the ellipsis notation.\n",
"Let's have a go:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "XiNGPNUrrEbL",
"outputId": "2d511cf8-4c95-46aa-ad07-d7c368042a9e"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(2, 3, 9)\n"
]
},
{
"data": {
"text/plain": [
"array([[[1., 2., 3., 4., 5., 6., 7., 8., 9.],\n",
" [1., 2., 3., 4., 5., 6., 7., 8., 9.],\n",
" [1., 2., 3., 4., 5., 6., 7., 8., 9.]],\n",
"\n",
" [[1., 2., 3., 4., 5., 6., 7., 8., 9.],\n",
" [1., 2., 3., 4., 5., 6., 7., 8., 9.],\n",
" [1., 2., 3., 4., 5., 6., 7., 8., 9.]]])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_3d = np.tile(x, (2, 3, 1))\n",
"print(x_3d.shape)\n",
"x_3d"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_U3B8sNLrEbM",
"outputId": "d0b196a5-8ca6-498e-c47d-98951f39fef2"
},
"outputs": [
{
"data": {
"text/plain": [
"array([[[0., 0., 0., 0., 1., 2., 3., 4., 5.],\n",
" [0., 0., 0., 0., 1., 2., 3., 4., 5.],\n",
" [0., 0., 0., 0., 1., 2., 3., 4., 5.]],\n",
"\n",
" [[0., 0., 0., 0., 1., 2., 3., 4., 5.],\n",
" [0., 0., 0., 0., 1., 2., 3., 4., 5.],\n",
" [0., 0., 0., 0., 1., 2., 3., 4., 5.]]])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"shift = 4\n",
"axis = -1 # Last axis\n",
"\n",
"x_3d_tmp = np.swapaxes(x_3d, axis, -1) # It's view, chill\n",
"x_3d_crop = x_3d_tmp[..., :-shift] # Same as x_3d_tmp[:, :, :-shift] for 3D arrays\n",
"\n",
"pad_shape = list(x_3d.shape)\n",
"pad_shape[-1] = abs(shift)\n",
"pad = np.zeros(pad_shape)\n",
"\n",
"shifted = np.concatenate((pad, x_3d_crop), axis=-1) # Beware, default is axis=0!\n",
"shifted"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "64_c_JHQrEbM"
},
"source": [
"Packaging it:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "T44TpVOFrEbM",
"outputId": "30b15ebb-44ac-49d3-f97f-e370413502a6"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Shift = 1\n",
"Axis = 0 :\n",
"[[[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]]\n",
"\n",
" [[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]]]\n",
"Axis = -1:\n",
"[[[0. 1. 2. 3. 4. 5. 6. 7. 8.]\n",
" [0. 1. 2. 3. 4. 5. 6. 7. 8.]\n",
" [0. 1. 2. 3. 4. 5. 6. 7. 8.]]\n",
"\n",
" [[0. 1. 2. 3. 4. 5. 6. 7. 8.]\n",
" [0. 1. 2. 3. 4. 5. 6. 7. 8.]\n",
" [0. 1. 2. 3. 4. 5. 6. 7. 8.]]]\n",
"\n",
"Shift = 0\n",
"Axis = 0 :\n",
"[[[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]]\n",
"\n",
" [[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]]]\n",
"Axis = -1:\n",
"[[[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]]\n",
"\n",
" [[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
" [1. 2. 3. 4. 5. 6. 7. 8. 9.]]]\n",
"\n",
"Shift = -2\n",
"Axis = 0 :\n",
"[[[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]]\n",
"\n",
" [[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0. 0.]]]\n",
"Axis = -1:\n",
"[[[3. 4. 5. 6. 7. 8. 9. 0. 0.]\n",
" [3. 4. 5. 6. 7. 8. 9. 0. 0.]\n",
" [3. 4. 5. 6. 7. 8. 9. 0. 0.]]\n",
"\n",
" [[3. 4. 5. 6. 7. 8. 9. 0. 0.]\n",
" [3. 4. 5. 6. 7. 8. 9. 0. 0.]\n",
" [3. 4. 5. 6. 7. 8. 9. 0. 0.]]]\n",
"\n"
]
}
],
"source": [
"def shiftnd(x: npt.NDArray, shift: int, axis: int = -1):\n",
" if shift > x.shape[axis] or shift < -x.shape[axis]:\n",
" return np.zeros_like(x)\n",
" x = np.swapaxes(x, axis, -1) # Relax, original x is unchanged\n",
" pad = np.zeros((*x.shape[:-1], abs(shift)))\n",
" if shift > 0:\n",
" out = np.concatenate((pad, x[..., :-shift]), axis=-1)\n",
" else:\n",
" out = np.concatenate((x[..., -shift:], pad), axis=-1)\n",
" return np.swapaxes(out, axis, -1)\n",
"\n",
"\n",
"for shift in [1, 0, -2]: # This should be definitely be a unit test...\n",
" print(f\"Shift = {shift}\")\n",
" print(\"Axis = 0 :\")\n",
" print(shiftnd(x_3d, shift, axis=0))\n",
" print(\"Axis = -1:\")\n",
" print(shiftnd(x_3d, shift, axis=-1))\n",
" print()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "nnMWEc-7rEbN"
},
"source": [
"We're almost done. We figured out how to shift by an integer amount in either direction. To shift a fractional amount, we will linearly interpolate the result of shifting by the floor and ceil of the fractional shift! Let $L_\\delta$ represent the fractional shift operator by $\\delta$ shift, and let $S_\\delta$ represent the integer shift operator by $\\delta$. Let $w = \\delta - ⌊\\delta⌋$ represent the fractional amount of the shift.\n",
"\n",
"$$\n",
"L_\\delta x \\equiv (1 - w) S_{⌊\\delta⌋}x + w S_{⌈\\delta⌉}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "4TCmU2vDrEbN",
"outputId": "48fdcb86-0206-4505-8108-2ed6f552d148"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n",
"[0. 0. 0.9 1.9 2.9 3.9 4.9 5.9 6.9]\n",
"[0. 0. 0. 0. 0. 0. 0. 0. 1.]\n",
"[0. 0. 0. 0. 0. 0. 0. 0. 0.1]\n",
"[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
]
}
],
"source": [
"def shiftnd_linear(x: npt.NDArray, shift: float, axis: int = -1):\n",
" shift_floor = int(np.floor(shift))\n",
" w = shift - shift_floor\n",
" return (1 - w) * shiftnd(x, shift_floor, axis) + w * shiftnd(\n",
" x, shift_floor + 1, axis\n",
" )\n",
"\n",
"\n",
"print(x)\n",
"print(shiftnd_linear(x, 2.1))\n",
"print(shiftnd_linear(x, 8))\n",
"print(shiftnd_linear(x, 8.9))\n",
"print(shiftnd_linear(x, 9))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "PsEHsV7KrEbO"
},
"source": [
"## A Simple `LinearOperator`\n",
"\n",
"In this section we will take what we learned about shifting and package that in a `LinearOperator` class.\n",
"What follows is a fairly generic template that can serve as a basis for implementing a wide array of linear operators."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"id": "8uaZorqtrEbO"
},
"outputs": [],
"source": [
"class LinearShift(pylops.LinearOperator):\n",
" \"\"\"LinearShift shifts an N-dimensional array by a single real value along a given axis.\n",
"\n",
" Since we are responsible coders, this is the entirety of our very long\n",
" and detailed documentation for this class. At least we're annotating...\n",
" \"\"\"\n",
"\n",
" def __init__(\n",
" self,\n",
" shift: float,\n",
" dims: Tuple[int, ...],\n",
" axis: int = -1,\n",
" dtype: npt.DTypeLike = np.float32,\n",
" name: str = \"L\", # Not required but good practice since v2.0.0\n",
" ):\n",
" # We need to store these variables for computing the shifts\n",
" self.shift = shift\n",
" self.axis = axis\n",
"\n",
" # Every operator needs a dtype and a shape. Since v2.0.0, we can instead give\n",
" # dims (model shape) *and* dimsd (data shape). `self.shape` is automatically\n",
" # inferred as (np.prod(dimsd), np.prod(dims)).\n",
" # Calling `super()` will also populate `self.clinear` (defaults to True) and\n",
" # `self.explicit` (defaults to False). See the latest documentation:\n",
" # https://pylops.readthedocs.io/en/latest/api/generated/pylops.LinearOperator.html\n",
" super().__init__(dtype=dtype, dims=dims, dimsd=dims, name=name)\n",
"\n",
" def _matvec(self, x: npt.NDArray):\n",
" x = x.reshape(self.dims)\n",
" y = shiftnd_linear(x, shift=self.shift, axis=self.axis)\n",
" y = y.ravel()\n",
" return y"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "fzMFZWBbrEbO"
},
"source": [
" In the new version, we use the name of the operator to `describe` it: "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 57
},
"id": "TMIpJC6arEbO",
"outputId": "18680c36-fd59-4466-e89f-ab8c4cab2135"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle L$"
],
"text/plain": [
"L"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"where: {'L': 'LinearShift'}\n"
]
}
],
"source": [
"from pylops.utils.describe import describe\n",
"\n",
"\n",
"LOp = LinearShift(0.5, dims=x_3d.shape)\n",
"describe(LOp)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uAde_-PYrEbP",
"outputId": "1e0fd95b-4646-474f-cab9-46ee5085b7e8"
},
"outputs": [
{
"data": {
"text/plain": [
"array([[[0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],\n",
" [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],\n",
" [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5]],\n",
"\n",
" [[0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],\n",
" [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],\n",
" [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5]]])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Taking it for a spin\n",
"shifted = (LOp @ x_3d.ravel()).reshape(LOp.dimsd)\n",
"shifted"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "d5M_097ErEbP"
},
"source": [
"Great! But what if we try the following?"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "fCvpbfmtrEbP",
"outputId": "4fa02a2d-069f-4552-9343-68d00c2c68c2"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Oops!\n"
]
}
],
"source": [
"try:\n",
" LOp.T @ shifted.ravel()\n",
"except AttributeError:\n",
" print(\"Oops!\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "jz7P4udHrEbQ"
},
"source": [
"This is because we need to define both the forward and the adjoint. The good news for us is that the adjoint should just be a negative shift."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "-tk6FvWPrEbQ",
"outputId": "00788517-30bb-45c5-bea5-5d69cea857f7"
},
"outputs": [
{
"data": {
"text/plain": [
"array([[[1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25],\n",
" [1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25],\n",
" [1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25]],\n",
"\n",
" [[1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25],\n",
" [1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25],\n",
" [1. , 2. , 3. , 4. , 5. , 6. , 7. , 8. , 4.25]]])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"class LinearShift(pylops.LinearOperator):\n",
" \"\"\"LinearShift shifts an N-dimensional array by a single real value along a given axis.\n",
"\n",
" Since we are responsible coders, this is the entirety of our very long\n",
" and detailed documentation for this class. At least we're annotating...\n",
" \"\"\"\n",
"\n",
" def __init__(\n",
" self,\n",
" shift: float,\n",
" dims: Tuple[int, ...],\n",
" axis: int = -1,\n",
" dtype: npt.DTypeLike = np.float32,\n",
" name: str = \"L\", # Not required but good practice since v2.0.0\n",
" ):\n",
" # We need to store these variables for computing the shifts\n",
" self.shift = shift\n",
" self.axis = axis\n",
"\n",
" # Every operator needs a dtype and a shape. Since v2.0.0, we can instead give\n",
" # dims (model shape) *and* dimsd (data shape). `self.shape` is automatically\n",
" # calculated as (np.prod(dimsd), np.prod(dims)).\n",
" # Calling `super()` will also populate `self.clinear` (defaults to True) and\n",
" # `self.explicit` (defaults to False). See the latest documentation:\n",
" # https://pylops.readthedocs.io/en/latest/api/generated/pylops.LinearOperator.html\n",
" super().__init__(dtype=dtype, dims=dims, dimsd=dims, name=name)\n",
"\n",
" def _matvec(self, x: npt.NDArray):\n",
" x = x.reshape(self.dims)\n",
" y = shiftnd_linear(x, shift=self.shift, axis=self.axis)\n",
" y = y.ravel()\n",
" return y\n",
"\n",
" def _rmatvec(self, y: npt.NDArray):\n",
" y = y.reshape(self.dimsd)\n",
" # negative shift!\n",
" x = shiftnd_linear(y, shift=-self.shift, axis=self.axis)\n",
" x = x.ravel()\n",
" return x\n",
"\n",
"\n",
"# Taking it for a spin\n",
"LOp = LinearShift(0.5, dims=x_3d.shape)\n",
"shifted_back = (LOp.T @ shifted.ravel()).reshape(LOp.dimsd)\n",
"shifted_back"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "c5q8NwIxrEbQ"
},
"source": [
"Are we done now? Nope! We need to ensure that the operations in `_matvec` and `_rmatvec` are actually transposes of each other. To do so, we can use the [`dottest`](https://pylops.readthedocs.io/en/latest/api/generated/pylops.utils.dottest.html) utility function."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Osnv0zh0rEbQ",
"outputId": "9da28a64-eb74-463e-848a-6d947cba326e"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Dot test passed, v^H(Opu)=-1.8444506555541245 - u^H(Op^Hv)=-1.844450655554124\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from pylops.utils import dottest\n",
"\n",
"dottest(LOp, verb=True, raiseerror=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "9IAXekvmrEbQ"
},
"source": [
"Since it passes the dot test, we're done!"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ts0N0wwlrEbR"
},
"source": [
"Interestingly enough, we see saw that applying the adjoint in succession to the forward is not an inverse"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 57
},
"id": "biq89j4vrEbR",
"outputId": "c6dbc86b-15e2-4134-f49d-3322a039fc1b"
},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle L^{T} L$"
],
"text/plain": [
"L.T*L"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"where: {'L': 'LinearShift'}\n"
]
}
],
"source": [
"describe(LOp.T @ LOp)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "XGVXG1d8rEbR",
"outputId": "ed0acd09-1de1-4df3-fc89-ae2644b3a7cf"
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"shifted_back_adj = (LOp.T @ LOp @ x_3d.ravel()).reshape(x_3d.shape)\n",
"np.allclose(shifted_back_adj, x_3d)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 17,
"referenced_widgets": [
"17fe01e151c14d2cb7a7fb6c46d23a32"
]
},
"id": "6n93iNbBrEbR",
"outputId": "a1c11263-e7dd-46cf-f57a-bab48a0d6388"
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "478d4a7f41ec43ffb5fd450979dc4b26",
"version_major": 2,
"version_minor": 0
},
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAEsCAYAAAAfPc2WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9cElEQVR4nO3dd3xUdf798ddNDyGFDoEkEEoSAgFCivRQBGmii3VZF4VVQZSqK7pfF8siuiIWVAQLYsFeKFKkF2mB0ElC751AEpKQMnN/f2SX1Z+gQCa5k8x5Ph55aG7mzpw7CB7ufc/nGqZpmoiIiIiIw7hZHUBERESkolHBEhEREXEwFSwRERERB1PBEhEREXEwFSwRERERB1PBEhEREXEwFSwRERERB1PBEhEREXEwFSwRERERB1PBEhEREXEwFSwRERERB1PBEhEREXEwFSwRERERB/OwOgCA3W7n+PHj+Pv7YxiG1XFERETExZmmSXZ2NsHBwbi5Xf/5KKcoWMePHyckJMTqGCIiIiK/cuTIEerVq3fd+zlFwfL39weKDyIgIMDiNCIiIuLqsrKyCAkJudxRrpdTFKz/XhYMCAhQwRIRERGncaOjSxpyFxEREXEwFSwRERERB1PBEhEREXEwp5jBuhZ2u52CggKrY5Rbnp6euLu7Wx1DRETEJZSLglVQUMCBAwew2+1WRynXgoKCqF27ttYaExERKWVOX7BM0+TEiRO4u7sTEhJyQ4t9uTrTNMnNzeX06dMA1KlTx+JEIiIiFZvTF6yioiJyc3MJDg6mUqVKVscpt3x9fQE4ffo0NWvW1OVCEREp1y7mF1LZ29PqGFfl9AXLZrMB4OXlZXGS8u+/BbWwsFAFS0REyh273WTOjqNMWpqGzW6yenR3qyNdldMXrP/S3FDJ6T0UEZHyKCe/iBnr9/PasjT2nskGwMPNYP/ZbMKr39hK66Wt3BQsERERcS3HL+Ty9srdvLt6Dxm5xSsJBPl68XD7RjzWKYK6Qc47OqSCJSIiIk5l69HzTFqayuebDlFoK15BILx6ZUZ1juT+m8Kdevbqv1SwRERExHJ2u8mC1ONMWprGkvSTl7e3b1iD0V2iuLV5XdzL0UoCKlil6PPPP2fQoEHs37//8tIIDzzwAJs2bWLVqlUEBgZanFBERMRaeQVFfJp8kNeWpZJ6MgsAdzeDO1qGMrpLJAn1q1uc8MaoYJWie+65h5deeokXX3yRyZMnM27cOBYvXsy6detUrkRExKWdzr7EOyt3886q3Zy5mA+Av48HD7ZtxPCkCMKqVrY4YcmUu4Jlmia5BTZLXruSl/t1fRLPMAzGjx/PHXfcQe3atZk8eTKrVq2ibt26pZhSRETEee06kcmkpal8mnyA/KLi+aqwqn6MSIpgcJtGBPg6/3zVtSh3BSu3wEblMV9a8toXX70bP+/re8v69OlD06ZNef755/npp5+Ijo4upXQiIiLOyTRNFqedZNKyVBbsOnF5e0JYNcZ0jeJPLULwcC8/81XXotwVrPJmwYIFpKWlYbPZqFWrltVxREREykx+oY3PNx1k0tI0th+/AIBhwO0xIYzuEknb8BoVdo3GclewKnm5c/HVuy177euRkpLCXXfdxQcffMBHH33EM888w9dff11K6URERJzDuYv5vLt6D2+tTOdk1iUA/Lw8GNSmISOSImhYwzkXB3WkclewDMO47st0Vjh48CC9e/fm6aef5t577yU8PJw2bdqQkpJCbGys1fFEREQcbvepLF5fnsZH6/aTV1g8Lx0c6MvwpAgeateIKpW8LU5Ydpy/qZRDGRkZ3HLLLfTr14+xY8cCkJiYSM+ePXn66adZsGCBxQlFREQcwzRNVuw5zaSlqczZcezy9lb1qjCmaxR3tgrFy8P17n+rglUKqlatSlpa2m+2//jjjxakERERcbxCm52vUg4xaWkaKUcyLm/v26wuo7tE0alxzQo7X3UtVLBERETkmp3Pzee9n/fy5op0jl3IA8DX052BieGM7BxJRK0AixM6BxUsERER+UP7z2bzxvJ0Plizj5yCIgBq+fvwWKcIHm7fiOqVfSxO6FxUsEREROSKTNNk7YGzvLoklR+2HcVumgA0Dw5idJdI7m1dH29P15uvuhYqWCIiIvIrRTY73209wqSlqaw/eO7y9lua1mF05yi6RdZ26fmqa6GCJSIiIgBk5RXywdq9vLE8nUMZOQB4ebhxX3wDRnWJJLpOkLUByxEVLBERERd3OCOHN5an8d6avWRfKp6vql7Zm2EdmjC0Q2NqBfhanLD8UcESERFxUcmHzvHqklS+2XIYm714viqyVgCju0Txl/j6+HqpJtwovXMiIiIuxGa3M3v7MSYtTWX1vjOXt3eNqM3oLpHcEhWMm5vmq0pKBUtERMQFXMwv5KN1+3l9WRr7zl4EwNPdjXtbhzGqSyQt61W1OGHFooIlIiJSgR27kMvkFelMXb2XC3kFAFSp5MWQ9o15tGMTgoMqWZywYlLBEhERqYA2H8lg0tJUvth0iKL/zFc1quHPqM6RDEwMx89bFaA0uZX0CWw2G8888wwNGjTA19eXhg0b8sILL2D+ZzEyubJ3332Xli1b0rx5c7y8vGjZsiUtW7bk7bfftjqaiIiUU3a7ydztR+ny5mJiX57Pp8kHKbKbdGxUkx8e6kjaM314pGMTlasyUOJ3+OWXX2bKlCnMmDGD6OhoNm7cyAMPPEBgYCDDhw93RMYKaciQIQwZMoRt27bx4IMPsn79eqsjiYhIOZVbUMQnGw7w2rI00k9lAeDuZnBXq1BGd4kiLqyaxQldT4kL1po1a+jXrx+9e/cGoH79+nz++eds2LChxOFcwc6dO4mOjrY6hoiIlEMns/J4Z+Vu3lm1h3M5+QAE+HjyULtGDE+KIKSKn8UJXVeJC1bbtm2ZNm0au3fvpkmTJmzdupXVq1czadKkq+6Tn59Pfn7+5e+zsrJKGsNptWvXjujoaKZNm3bFn+/YseOKBeuP9hMREde1Z8d6TiyahO3cAT7M+ivn7FWoX82PkUmRDGrTEH8fT6sjurwSF6yxY8eSlZVFZGQk7u7u2Gw2xo8fz4ABA666z4QJE3juuedK+tJOz263s3XrVu67776rPmbnzp0MHTr0uvcTERHXYtrtbFr+Dax7m7jCLTQG8IRPa87ibO8PuC2mHh7uJR6tFgcp8a/EV199xWeffcbMmTNJSUlhxowZTJw4kRkzZlx1n6eeeorMzMzLX0eOHClpDKeUnp5OTk4OsbGxV33Mlc5g/dF+n3/+Ob6+vpw4ceLytgceeICYmBgyMzMdE15ERJzCpUIbM1dsZO8LzYlb9SBxhVuwmQYbfG7CNNxJ8jvOHRF+KldOpsRnsJ544gnGjh3LPffcA0Dz5s05dOgQEyZMYODAgVfcx9vbG29v75K+tNNLSUnBw8ODmJiYK/48Ly+P8+fPU69eveva75577uGll17ixRdfZPLkyYwbN47Fixezbt06AgMDHX4cIiJS9s6cz+TddYd5a+VuTmfnsSnIINvDm03V+tCw9+MkhDeFPT9BWDvw0qyVsylxwcrNzcXN7det2d3dHbvdXtKn/n0FOVf/mZs7ePhc22MNN/D0/ePH3sB/vCkpKTRt2hQfH58r/jw1NZXIyMjr3s8wDMaPH88dd9xB7dq1mTx5MqtWraJu3brXnVFERJzL/tRNHF04kcbnV/PvjGe5aPpQL8iPrQkv0bB9IklVa/zvwY27WxdUfleJC1bfvn0ZP348oaGhREdHs3nzZiZNmsSgQYMcke/qJgRf/WeNu8Ofv/7f9xMbQWHulR8b1h7u//F/37/RHHLP/fZx467/0ltKSsrvXh6MjY3l559/vu79APr06UPTpk15/vnn+emnn/RJRBGRcsy029m8ejZFP08moWAj4QBuMLrOXiK7P8wdrULx/L1LgKYdNn8CDbtAYEhZxZbfUeKCNXnyZJ555hkeeeQRTp8+TXBwMA8//DD//Oc/HZGvXNuyZQv9+/cvlf0WLFhAWloaNpuNWrVq3WhEERGxUEF+HhvmvUf17e8Tax4CwG4aJPsk4N1xOM/e1AvD7Rpmq+Y9ARvfh6b94M6PSzm1XAvDdIIl17OysggMDCQzM5OAgIBf/ezSpUscOHCABg0a/PqSmZNfIty3bx+NGjXi559/pm3btg7dLyUlhaSkJKZOncpHH31EQEAAX3/99RUf+0tXfS9FRKRMZeTkM+3nvXy/4mfWeP4dd8Mkx/RiY5VbCOv5OPWbtLi+Jzy5HaZ1LD6Tdd8PEN65VHK7kt/rJtei/K6Vfz2Fp7Qe+ztSUlKA4nm0HTt2/O/pvbxo0qTJDe938OBBevfuzdNPP829995LeHg4bdq0uabLiiIiYq1De7exZsm3/C09itwCG+DHjKpdCK/fiBa3jqZT9do39sS1m0P8g7BhKsz/Owz5Gdy9HJpdrk/5LVhO7r9F6aabbvrV9vbt27Nq1aob2m/WrFnccsst9OvXj7FjxwKQmJhIz549efrpp1mwYIEjD0FERBzAtNvZtm4BeSvfIOHSeu4Gnikah1/dJozpEsU9re/By8O95C/U+WnY+R2c3Q3rpkC7ESV/Trlh5fcSoVw3vZciImWnsCCf5AUfErjlPaLNfZe3J3u1xtb1eRLj22EYhmNfdMtnMOsR8PSDRzdCwO98IEx+l+teIhQREXFCmXkF/LBwHt22jKKtkQFAnulJcuDN1L1lDPFRcaX34i3uhU0fwdENsOj/oP+Hpfda8rtUsERERBzg4OkM3lh1gPfX7MWen8ORapc4QwA7691JdN/RdKxV74+fpKQMN+g1EX4cDTcNK/3Xk6tSwRIRESmBHRsWkb38dXyyD/L6+acAg6a1a7Kq9VR6JHUlybeMV1mv0wIGLwZHX36U66KCJSIicp1sRYVsWPgxfilTibGnF2/0gGGNLtKn+630iKrj+Pmq6/HL1y7I0a10LKCCJSIico2yMzPYNGcy4Xs/oY1xBoB804MN/l2o1X00bzVvY3HCX7AVwPKXIGVG8bIN/je4BITcEBUsERGRP3D0fC5vrkhj77of+M73DTDgnFmZ7XX+RFTfx+kQHGZ1xN9y84ADyyH3LCz+J9w+zepELqXcFCwnWE2i3NN7KCJyfVJTlvPT2nU8nh5Kkd0EmrDItzXeETcT13cYSX7X//H9MvPfgff3usC2LyH2fgi79juLSMk4fcFydy9efK2goABfX98/eLT8ntzc4htee3p6WpxERMR52W02khd/hnfyFFradlHT7sfT9n/RvnEoo7tE0jV6AG5u5WSAPDgWYgdCykcw73F4eGXxmS0pdU7/Lnt4eFCpUiXOnDmDp6cnbtdy00v5FdM0yc3N5fTp0wQFBV0urSIi8j85F7PYNOctQnbPIJGTABSabqT6J7D2z+2IiYiwOOEN6joOUmfB6Z2Q/D4kDrE6kUtw+pXcofjs1YEDB7Db7RakqziCgoKoXbu2tZ9sERFxMicy81g8+0N67xtPVSMHgAtmJbbU6kdE78epE9rI4oQOsPGD4rWxvAPg0U1QuabViZyeS6zk7uXlRePGjSkoKLA6Srnl6empM1ciIr+w7dApJq3cz8yNB6lPPvdVzeGQWZODje+j9a2PkeRfxeqIjhN7f/GnCc/shuMp0OQWqxNVeOWiYAG4ubnp/nkiIlIidpudTcu+xG39O+y56MmM7MEA1AxvxoqWH9O+Yy/CPCrgnKqbO9w2FbwrQ2CI1WlcQrkpWCIiIjfqUl4OyXPeoU7qdOI5BkC0tweDwoN46OZEEutXtzhhGagZZXUCl6KCJSIiFdaZU0fYMXsSzY59TQcjG4As04eUGn1p1GsMHzRw0dJx6Ge4eBqib7c6SYWlgiUiIhXOrhOZvLYslSpb3+Pflb4BA46a1dgbPoBWt44gKcgFzlhdzd7F8Fl/8AmE+h3Az4Xfi1KkgiUiIhWCabeTsuoHvtp8hH/vLS4NAUYb7qycRmHL+0i4ZRD1PL0sTukEwpOgVnM4tR2WPAu3vmV1ogpJBUtERMq1/LxckudNo8bOD2htHobCEF4xxnJ7TPHCoHHhg7U8zS+5eRSv8D69B2z+BFo/AHVbW52qwlHBEhGRcinjzAm2zZ5E1JEvaW9kAnDR9OZirTj2DutBeJ0aFid0YqE3QYt7YevnMG8MDF5S/ElDcRgVLBERKVd2n8piy3cv0efUNJKMQjDghBlEeti9tOw3ik5Va1kdsXzo9jyk/QjHN8Pmj4vPZInDqGCJiIjTM+12Vu05zsRle5m78xh9PG3cFVhImhHGueYPEd9rMEneul/tdalcE5KehoVjYclz0PQ28K1Ai6taTAVLREScVmFBPhvmvU/Vbe+x+GI0c3J7AWA07sHm5h1o2bYXhu5Re+MSHoS9i4ovF/oEWZ2mQlHBEhERp5OZcYbNsycRcfBz2hnnARjkk8XZVo8xoktTImpd/73h5ArcPOAv31mdokJSwRIREadxeN8ODsybSOtz80gy8sGA02YAu0Luonnf0bxTs67VESu2/CzwqgyGzgqWlAqWiIhYbs3+M0xamkrPfRMY7LMGDNhDPU5GDyah9xCSfCtZHbHi2/EtLBgLXf8Jre6zOk25p4IlIiKWKCosIHnhDN7c5cYXh30ASHXvQgv/PIw2w4jt1J/Gmq8qO1nHIec0LB4HkX008F5CKlgiIlKmsi5ksHnOG4Tv+5Q2xln2XErgO49B/CW+AaM696ZZ8BNWR3RNiUNgy6dwJg2WjS9ejFRumAqWiIiUiWOH0tkzdyKtzsyhk5EHBpw1K1M3PIbDd99GrQAts2Apd0/o+Qp83Bc2fgCxf4XaMVanKrdUsEREpFQlHzrHmW9H0T1rDnUNOxiwnzocjXyA+D6P0NXP3+qI8l8NOkL0n2DndzDvcXhggQbeb5AKloiIOJytqJA5O44xadluVu07w1hfk16V7aR4NMeW+AitO99NuLtuzeKUuv8Ldi+EI+th6xfQ8s9WJyqXVLBERMRhcrIz2ThnMmF7PubTzL6sKojF092NC03/Qnqrh4ht2cHqiPJHAupCxydgybNwaqfVacotFSwRESmxk0f3kzpnIi1PfU8nIxeAoX5raZJ0H8M6NqFukJZZKFfaDIP67aFevNVJyi0VLBERuWHpW1dzZvEkErKX09mwgQGHzJocbDKQm/o+Slf/IKsjyo1w91K5KiEVLBERuS52u8n8XceZtDSVf51+nPaeB8CAre5R5MUNJb7bnwnz8LQ6pjjKhUOQ8gl0/gcYhtVpyg0VLBERuSZ5Odkkz32HJ9LD2HDaBsCrPjczKjCdKl1G0qJ1F4sTisMV5sK0JMjLgOpNIOYuqxOVGypYIiLyu86cOMzOOa/S7Pg3dDQu0vliP9J8+vBg20YMT7qN0Kp+VkeU0uJZCW56BJb9Cxb9H0TcAt660fa1UMESEZEr2rtzPScWvkpC1hKSjCIw4IhZnS6xMTx92+0E+OoyoEtoOxy2zoSM/bDiZeg+3upE5YIKloiIXGaaJotTj+P/w1+5qXAjjQAM2O7WmOxWD5PQYyAhnl5Wx5Sy5OENt/wbZt4B698tvhF0jUirUzk9FSwRESE/P5+ZKUeZtDSVHScy+djfjXhvgw2+bamcNILmiT2sjihWanwzRPSG9B+LV3j/6xwNvP8BFSwRERd27vQxts+ZRNSRr5l4/jF22YLx8/LgQPPhHEtsQJuGzayOKM7ilgmwbwkcXAWps6FpP6sTOTUVLBERF3QgLYUjCyYSf+EnkoxCMGBkYDIXOj7Lg20bEVRJlwHl/xMUBp3GQmEeNO5udRqnp4IlIuIiTLudLT/PpnD1WyQUJNMAwIBdRjjnWzzI/T0H4+nlbXVMcWbtR1mdoNxQwRIRqeAKimx8lXKYt5ZsZfalYdR0u4jdNEj2ice7wwhatOmF4eZmdUwpb+w2uJQJlapancQpOaRgHTt2jCeffJL58+eTm5tLo0aNmD59OnFxcY54ehERuQEXMk6xYe5UBu9uytEL+QC8Ubk7PeoUEnLL4yRGtLQ2oJRfp3bCrKHgWxX+8r0G3q+gxAXr/PnztGvXjs6dOzN//nxq1KjBnj17qFKliiPyiYjIdTq0dzsH500kLmMe3Y0CWuYOoSgggcc6RfBwuzuoVlmXAaWEPH3hdBrY8jXwfhUlLlgvv/wyISEhTJ8+/fK2Bg0alPRpRUTkOph2O9vXLyR35Rsk5K0jzDDBgHQjlOFdmtOxx214e7pbHVMqiqrh0G4ErPw3LHwaGnUDL63o/0uGaZpmSZ6gadOm9OjRg6NHj7JixQrq1q3LI488woMPPnjVffLz88nPz7/8fVZWFiEhIWRmZhIQoCX4RUSuVZHNztzk7TRaNIhm9r2Xtyd7xeLe5lFadbxd81VSOgpz4e1EyDwM7cdA139ancihsrKyCAwMvOFuUuLfdfv372fKlCk0btyYhQsXMnToUIYPH86MGTOuus+ECRMIDAy8/BUSElLSGCIiLiUrJ49JS1Jp+Owsbv90OwX5+VwyPVgZcAt771xM/FPLiE3qr3IlpcezEvR4sfjf106Gc/uszeNkSnwGy8vLi7i4ONasWXN52/Dhw0lOTmbt2rVX3EdnsEREbszRg2nsm/sKDc8soWnG/5Ft+lKjsjfPtXbjjo7x1Kilv7BKGTJN+Kx/8QKkjbrBn7+pMAPvJT2DVeIZrDp16tC0adNfbYuKiuLbb7+96j7e3t54e2vIUkTkWu1MXkzmstdJzF1NPcMENxhRYxdh3R5hQFx9fL206o5YwDCg57/hnZsg7wLkZ4FPoNWpnEKJf0e2a9eO9PT0X23bvXs3YWFhJX1qERGXZisqJPmnj6m0aSox9v/8OWvAJo8Y7ImP8Fznu3Bz1+C6WKxaIxi8COq0AEOXpP+rxAVr1KhRtG3blhdffJG77rqLDRs2MG3aNKZNm+aIfCIiLudifiHT1+7nq+WrWc4Y3A2TAtOd9f5dqHXzaFrHtLU6osivBbeyOoHTKfEMFsDcuXN56qmn2LNnDw0aNGD06NG/+ynC/19Jr3OKiFQEJw7vZfXir3gotTEX8goA+DDoSxqEhBLV93FqBde3NqDIHym4CCsnQuxfi5dyKMdK2k0cUrBKSgVLRFxZ6uaVZCx5lYSLK3HHJOL8OIyqDRnVJZK/JoTj5635KiknvnsQtn8FTW6Be7+0Ok2JWD7kLiIi189us7FxyUy8NkyhpW1n8UYDNrtH817/ZnRs3w03t4rxaSxxIR0eh53fwe4FxV9NbrE6kWVUsEREylBuQRFzliwiYcNIEjgJQKHpxvrKnajWZSStYpOsDShSEjUi4KZHYM2bsOBJCE8CDx+rU1lCBUtEpAycvHCRt1fvY8qqPeTlZHGkWhaZhi+ba95Gkz5jaB/a2OqIIo7R8e+w/Ws4fxB+fhM6/d3qRJbQDJaISCnas20NpxZPwvd8OnHnnwQMGlSrzEux+fTsfDP+AVWsjijieDu+gW8HF5+9GrYBgsrf0k2awRIRcTKm3c6m5V/jtu4tYgu30RjAAx4Ky6B7t9u4rUU93HULG6nIovvDpo/g4CpY8hz0/9DqRGVOBUtExEEu5eWQPHcKdXZ9SBzHALCZBusrtScgaQRTE262OKFIGTEM6PkKrHoFbn7B6jSWUMESESmhM9mXmLJ6DztWf8dXXpMAyDJ9SKneh0a9H6dtgyiLE4pYoGaUS565+i8VLBGRG7QvdSMLV69kzK5gLhXagIYsrRaDW3gnWvUbSVJQdasjijiPrOMQEGx1ijKjgiUich1Mu53Nq2ZhWzOZ+IJN3GOvxBOF44kLrcOYLlF0aPVnPN01XyVyWUEOzB4GuxcWD7wHhlidqEyoYImIXIOC/Dw2/DiNGjs+INY8BIDdNNjtG8Oyh1oT3zwGw9DCoCK/4VkJsk9BYS789A+482OrE5UJFSwRkd+RkZPPT3M/JSn1edobFwDIMb3YWLUn9Xs9zk2NYqwNKOLsDAN6vQJTO8KuWbB/GYR3tjpVqVPBEhG5gj0nz/P6ij18tH4/dW2ZpFXJ5IQZRHrYPbToO4pO1WtbHVGk/KjVDOIfhA3vwrwnYOgacPeyOlWpUsESEfkP025n29p5XFr1JvuzbLyTNQgAv3qRLGnxLp269CPJ29filCLlVOenYOe3cG4PrHsH2o20OlGpUsESEZdXWJBP8vwPCdo6jRbmfgBivdyYG/Uof7u5DUmNa2m+SqSkfIKg23Mw6xFY8W9oficE1LU6ValRwRIRl5WZcYYts1+j8cHPaWtkAJBnepIc1J16PcbwWVRrixOKVDAt7i1e4T3zKFw4ooIlIlKRHDh7kTeWp+G9aSov+3wFBpwxA9hZ7y6a3zqajjUr7h/6IpYy3KD/B1CpKnhVtjpNqVLBEhGXsX39Qr7etJ/xqYHYTRN/I5F7fDZzMfrPxPd+mCRfP6sjilR8QaFWJygTKlgiUqEVFRaQ/NPH+Ke8S3P7HgoKQ3jBHEv3yGDGdI2iZeRgzVeJWMG0w5aZYNogdqDVaRxOBUtEKqTszAxS5rxJ+N5PaGOcBSDf9CC3alN2PNSV6NA6FicUcXFpc4tXePfyh8Y9wL9iLX2igiUiFcqR8zms+/ZVuh99h05GHhhwzqzM9uD+NO37OB3quMblCRGnF9kH6raGY5tg8T/h9mlWJ3IoFSwRqRA2HjjNpOW7+WrzYXp65HBnYB77qcPRyAeI7/MISX7+VkcUkV8y3KDXRHivC2z7EmLvh7C2VqdyGBUsESm37EVFJC/+DJ+NU5iTFc7nub0ByA3ryvpmrYlP6k+4u7vFKUXkqoJji+evUj6CeY/DwyvBrWJUk4pxFCLiUnKyM9k0ZzKhez4mkVMAVPM5wcHooYzsEk2rkKoWJxSRa9Z1HKTOgtM7Ifl9SBxidSKHUMESkXLj5NH9pM2dSIuT39PRyAXgglmJLbVuI6Lv48yo19DihCJy3SpVhS7/hB9HwfIXodVfKsQaWSpYIuL0th49z6SlqXRM/ReDfX4GAw6ZNTnYZCCt+w4jyb+K1RFFpCRiBxYPu8cOrBDlClSwRMRJ2W02Ni37kik7Cpm+t/iPqo3unYnzy+BS6yHE3TyAMA9Pi1OKiEO4uUO/t61O4VAqWCLiVPJyskn+cQp1U6cTz3HSLiXwsdsD3NkqlFGde9Ci/hNWRxSR0nbhEATUKy5e5ZQKlog4hTOnjrBz1qtEH/+GjkY2AJmmL8Ehjdh3d1/CqmmZBRGXsOZNWPov6DEe4h+0Os0NU8ESEUvtOpHJga//TteMb0gyisCAo2Z19jYcQGzfkXQN0icCRVyKhy/Y8mHpC9D0dvCrbnWiG6KCJSJlzrTbWZx2gknL01mw6wRP+ubRu3IRO9wakd1qCPE9BlLP08vqmCJihbhBsHkGnNwOS56FW9+yOtENUcESkTKTn5dL8ryp1NzxAVMze7CgIBY3w+BE43vYHjOA5ondQTdeFnFtbu7QcyJM7wGbPyle4b1enNWprpsKloiUuozTx9k2ZxJNj3xJeyMLgEf9VlO37b2MSIogvLrmq0TkF0Jvghb3wtbPYd4Y+NvScjfwroIlIqXmQPpmjsyfSNyFhSQZhWDACbMK6fXvpdWto0iqWtPqiCLirLo9D2k/woktkDKj+NJhOaKCJSIOZZomK/acZtLSVJ48OoqOnvvBgFSjARkxD5LQ62/U8fK2OqaIOLvKNSHp6eI5rKJLVqe5bipYIuIQhQWX2DDvff4vrRbLjxYA4OndFc/K1fFqP5wWbXtjuLlZnFJEypWEByGyNwSFWp3kuqlgiUiJXMg4xZZZrxFx6AvaGedJvNiP9Z69uP+mcEYm9aVJrQCrI4pIeeXmUS7LFahgicgNOrxvOwfmvUrrc/NIMvLBgFNmIJ2aNeaJ/rdTrbIuA4qIAx1eC5s/hVsng+H8Z8NVsETkmpmmydr9Z7B/M4i2easJNUwwYDchnG72N+J7PURP30pWxxSRiiY/C2beVfzPkESI/avVif6QCpaI/KGioiK+23aMSUtTWX/wHDP882nvY5LsFYt7m0dp1fF2mmi+SkRKi3cAdPw7LPq/4qH3qD7g69x3eVDBEpGryrqQwebZr9Fw/0yePz+EnbZgvD3cSIsYyr748cRHJ1gdUURcReIQ2PIpnEmDpeOh96tWJ/pdKlgi8hvHDqaz58dXiD0zh07GJTBgZMAajrV5jqEdmlDT38fqiCLiatw9oecr8HFf2PRh8WXCOi2sTnVVKlgictmujUu5sPQ1EnJXU9ewgwH7COZ41CAG9B6Cr59WXBcRCzXoCNF/gp3fwaaPoM9rVie6KhUsERdns9uZvf0Yby7ZxheZD9PULRsMSPGIwZb4CK0730VD9/J1iwoRqcC6j4fwJGh1n9VJfpcKloiLuph9nvU/vs+QtEbsPZsLwJt+XelZ8yI1uo0mtkU7ixOKiFxBQDDEDrQ6xR9SwRJxMSeO7iN9ziu0PDWLrkYuUVlDyKgUx9AOjRnW8U/UCfS1OqKISLnn8M9Vv/TSSxiGwciRIx391CJSAulbVvHzxNup/n4cSac/J8jI5SC1GNq+EYdfuJ1/9W2pciUi4iAOPYOVnJzM1KlTiYmJceTTisgNsttNftqcSq0FD9KqaAcRAAZscW9KQdxQ4rr9mfoeOpEtIuJoDvuT9eLFiwwYMID33nuPf/3rX456WhG5Abn5BXySfIjXlqWRfiqTjUGZFHm4scGvA1W6jqJlbGerI4qIVGgOK1jDhg2jd+/edOvW7Q8LVn5+Pvn5+Ze/z8rKclQMEZd2+sQhds15lbBj83giYyzZpi+Bvl783PT/qNM+lrahTayOKCLiEhxSsL744gtSUlJITk6+psdPmDCB5557zhEvLSLAnh3rOPnTqyRkLSXJKAI3GFF1GzU6P8IDNzXE38fT6ogiIi6lxAXryJEjjBgxgkWLFuHjc22rOz/11FOMHj368vdZWVmEhISUNIqISzHtdlJWfANr36J14VYaAxiwzS2CnNiHebbHfbh7eFkdU0TEJRmmaZoleYIffviB22+/HfdfLERos9kwDAM3Nzfy8/N/9bMrycrKIjAwkMzMTAICAkoSR6TCu1RoY+bGg8xcspqFRSNwN0xspsGGSu3wTxpJs4SbrY4oIlLulbSblPgMVteuXdm+ffuvtj3wwANERkby5JNP/mG5EpFrc/bUUVYu+oqhO8M4nX0J8OCTwHY0qFOLhr2foE2DKKsjiojIf5S4YPn7+9OsWbNfbfPz86NatWq/2S4i129/6kaOLnyV+AuLuI0ixub+E+8qDRiRFMntbe8k0FeXAUVEnI0WwBFxQqbdzpbVsyj6eTLxBZsIBzBglxHOW30b07lzbzzdHb5OsIiIOEipFKzly5eXxtOKVHgFRTbmr1xBxOqRtDIPAWA3DTb4JOLbcTgxN/WkqZuKlYiIs9MZLBEnkHExj2lr9jN5RToXMi9wpNoZcgwvNlbtSf2eY7ipcQurI4qIyHVQwRKx0KG92zg4byKVzmzlqfNPAAbBgUHMb/YKvTp3o1P12lZHFBGRG6CCJVLGTLudbevmk7fyTRIurSfMMMEDBtc5Saeb/8TdsWF4eejTtyIi5ZkKlkgZKSrIZ8OCDwnaOo0W9v3FGw1I9mqNR7vHeK99PwzNV4mIVAgqWCKlLDOvgPfX7GPTim+Z6f4KAHmmJ8lBN1O3xxjio+IsTigiIo6mgiVSSo4eSGX+iqWM3l6di/lFQH2GVI3CHtqWZn1H07FWPasjiohIKVHBEnGwHRsWkb38dRJyf6a/6cvI/PFE16nJ6C5RJMTdi4+n5qtERCo6FSwRB7AVFbJh4cf4pUwlxp5evNGAA16NmTewOR3j4jAMw9qQIiJSZlSwREog+1IhPy34ivgt42hjnAGgwHRnQ0AXat48htbN21icUERErKCCJXIDjmbk8ObKdKb9vJcaBadJr3KWDNOPbXX607TvGNoH17c6ooiIWEgFS+Q6pKYsJ2PJaxw+n8MrWQ8AUKtmQxbEvEbSzf1J8guwOKGIiDgDFSyRP2C32Uhe/BneyVNoadsFwE1eBt82epj7u7ahV3Rd3Nw0XyUiIv+jgiVyFbkXM9k45y1Cdn9MIicBKDTd2FC5E1W7juabVh0tTigiIs5KBUvk/3MiM4+3V6ZTtPZdXvL+HIALZiW21OpHRO/HaRfayOKEIiLi7FSwRP4jfevPfL0hled3+FFos+NvxPEXr7Wca9yf1rc+RpJ/FasjiohIOaGCJS7NbrORsuxL3NZPIbZoG70KQ3jGNpZ24TUZ0zWKqOaDcNf9AUVE5DqpYIlLupSXQ/KctwlOnU4cxwEoMt24FFCfDQM7EN84zOKEIiJSnqlgiUs5nX2J1d9PpsP+N+hgZAOQZfqQUqMvDXs9TtsGkRYnFBGRikAFS1zCruMXeG15Gp9sOEA3t3P8KTCbo2Y19oYPoNWtI0kKqmZ1RBERqUBUsKTCMu12Nq/8Hvuat5h9IZj3c3sBcCa0I6uiI2jTfQD1PL0sTikiIhWRCpZUOPl5uSTPm0aNnR8Qax4GoK7vXnY1+hujukbTNryGbrwsIiKlSgVLKoyMMyfYNnsSUUe+pL2RCcBF05tN1XpRv+cYvmnU3OKEIiLiKlSwpNzbfSqL15enEbv1Bf7mvRIMOGEGkR52Ly37jaJT1VpWRxQRERejgiXlkmm3s3XtPKZtyeTdNDBNiHLvRFufI5xv/jfiew0mydvX6pgiIuKiVLCkXCksyGfDvPepsu19Wpr7aXspninmA/RpVpcxXboR1XiM5qtERMRyKlhSLmRmnGHz7NdocnAm7YzzAOSZntStXZe0EX2IqB1ocUIREZH/UcESp7b/bDa7vv4nnU59RpKRDwacNgPYFXIXzfuOpnPNulZHFBER+Q0VLHFKa/adZtKyNL7fepTHfTLoUzmfPdTjZPRgEnoPIcm3ktURRURErkoFS5xGUWEByQtnUHnzVF4734lvC2IB2Bfan03N+hDbqT+NdeNlEREpB1SwxHJZFzJImfMGDfd9ShvjLACPVYKA1ncyqnMkzYKDrA0oIiJynVSwxDLHDqWzZ+5EWp2ZQ5KRBwacNSuzI/gOovuO4YM6oVZHFBERuSEqWFLmkg+dY9LSVB7ZN4Ikz71gwD6CORZ5P/F9HiHJz9/qiCIiIiWigiVlwlZUyMZFn/H8Ln/mHbgEQIFXZypV88We+AitO99NQ3d3i1OKiIg4hgqWlKqc7Ew2zplM2O6PSTROEXPxVha59+Le1mGM6nwLLUOqWR1RRETE4VSwpFScOLqP9DkTaXHqBzoZuWDAebMSHZrU49E7+1E3SMssiIhIxaWCJQ615cg5Lnw1lHbZi6hj2MGAg9TiUOOBxPV9lF7+WnFdREQqPhUsKTG7zc781BO8ujSVZbtP8ZH/BZJ87Gx1j+JS3FDiuv2Z+h6eVscUEREpMypYcsNyc7LZOOdt6qbP4MmMQey0BePuZrC1wWB2tR5Li9ZdrI4oIiJiCRUsuW6nTxxm15yJNDv+LR2NiwCMqryKtLhneaxTBKFV/SxOKCIiYi0VLLlme3as5+RPE0nIWkqSUQQGHDGrs6/RfdzVdzj+gVWtjigiIuIUVLDkd5mmyaK0k7y+eDsfnh1MY7csMGC7W2Muxg4hvvtfCfH0sjqmiIiIU1HBkivKz8tl7fzpPLYrhB0nswF4q1ISvaufo3KnETRP7GFxQhEREeelgiW/cvb0MbbPnkT00a9IMrJokDmEg96xDG7TkEGd+hJeI8DqiCIiIk5PBUsA2J+WwpEFE0m48BOdjUIw4LhZlb8l1OXjW28nqJIuA4qIiFwrFSwXZpomq3btw2f2wyQUbCQcwIBUowHnWzxEfM/B3OrlbXVMERGRcqfEBWvChAl89913pKWl4evrS9u2bXn55ZeJiIhwRD4pBQWFRXy5+TCTlqax5WgGyUGnsXsYJPvE49NhODFtemO4uVkdU0REpNwqccFasWIFw4YNIz4+nqKiIp5++mm6d+/Orl278PPTekjO5Py5U2yZ/Rp1D83i0XNjyDJ98fX0YGmjJ6jWJobEiJZWRxQREakQDNM0TUc+4ZkzZ6hZsyYrVqygY8eO17RPVlYWgYGBZGZmEhCgIWpHO7R3GwfnTSQuYz5+RgEA/yj8M34dhvFwu8ZUq6zLgCIiIr9U0m7i8BmszMxMAKpWvfqik/n5+eTn51/+Pisry9ExXJ5pt7N9/QJyV7xBwqX1hBkmGJBuhHKm2WD+2fMhvH0rWR1TRESkQnJowbLb7YwcOZJ27drRrFmzqz5uwoQJPPfcc458afmPIpudbzYf5uMla5idNwwPww4GJHu1xr3tY7Tq0I8IzVeJiIiUKodeIhw6dCjz589n9erV1KtX76qPu9IZrJCQEF0iLIHMC2dZvvArhu+ow+HzuQDMCPiE+jWCCO4xhkZN4y1OKCIiUn44zSXCRx99lLlz57Jy5crfLVcA3t7eeHtr7scRjh5IZe+PE4k9O5d+xiWeyBpHjcqhDOvYhFvaf0/NAF+rI4qIiLicEhcs0zR57LHH+P7771m+fDkNGjRwRC75AzuSF5O97DUScn+m3n/mq/ZSl1d7hHLzzbfj4+ludUQRERGXVeKCNWzYMGbOnMmsWbPw9/fn5MmTAAQGBuLrq7MnjmSz21m0di3BS0cRY08v3mjAJs8WmDcNo3XSnTTSfJWIiIjlSjyDZRjGFbdPnz6d+++//5qeQ8s0/L6Llwr4cN1+Xl+WzslzGRyp9g/8jUus9+9CrZtH0ySmrdURRUREKhTLZ7AcvIyW/MKJw3tJnzsR75PJjMgYBRhUrVSZ2U2ep1fnbnQIrm91RBEREbkC3YvQCaVuXknGkkkkXFxBHcMO7vDXGkdo06U/f00Mp5KXftlEREScmf5P7STsNhvJi2filTyFVradxRsN2OwRTWH8UKZ3/TNu7hpcFxERKQ9UsCyWW1DEjPX7+Xnp93zKiwAUmm5sqNyJql1G0io2ydqAIiIict1UsCxy6vgh5i/7iTFbg8jILQDq8kjVxhTUjiOizxjahTa2OqKIiIjcIBWsMrZ7+1pO/fQqCdnL6Gd68mjueBpUq87IzhHE3LSGyj5eVkcUERGRElLBKgOm3c6m5V9jrH2b1kVbaQLFN152b8i390TQrW1b3LV+lYiISIWhglWKLhXaWLToB6KSnyGOYwDYTIMNldrhnzSSmISbibE4o4iIiDieClYpOJOVxzur9/D2yt0E5B5hd9XjZONNSvU+NOz9BG0aRFkdUUREREqRCpYD7UvdyLEFEzlyJoNnswYC4FOlPnOiJpDU4046BVW3OKGIiIiUBRWsEjLtdjavmoVtzWTiCzbRELB5GXwe8lfu69ae/i1D8XDXfJWIiIgrUcG6QQWXctkwbxo1dnxArHkYALtpkOybiG/H4cxJ7ImhwXURERGXpIJ1nTJy8nl39R6yV73LBM9PAMgxvdhYtSf1ez1OYiONrYuIiLg6FaxrdHD3Vr5es5lnt/mQW2CjstGK+6su4UT9W2nRdxSdqte2OqKIiIg4CRWs32Ha7WxbO49Lq94g/lIy3Yrq8veCp2hZryqju0TSoNVAIjz1FoqIiMivqR1cQWFBPsnzPyBo63u0MPcXbzSg0K82K+5JpEN0IwzDsDakiIiIOC0VrF+4kFvAytlTiUubRFsjA4A805PkoO7U6zGGhKjWFicUERGR8kAFCzhwJps3VqTzwdp9dDKPcWtgBmfMAHbWu4vmt46mY826VkcUERGRcsSlC9b29Qu5uOIN5p+tyhu5PQE4VLsty5qG0qbnQJJ8/SxOKCIiIuWRyxUsW1EBGxZ+jH/KuzS37wGgoW9lNoXcx4huMdwcWVvzVSIiIlIiLlOwsjMzSJnzJuF7P6GNcRaAfNODDQFdqd19DD82S7Q4oYiIiFQULlGw3l+zF/PHMTzotRwMOGdWZntwf5r2HUOHOmFWxxMREZEKxiUKVm1/H/5+sQNdqu7mWMRA4vo+QpJfgNWxREREpIJyiYLVK7ouXkMG0KDJaBrqxssiIiJSylyiYLm5GXSPqmN1DBEREXEROp0jIiIi4mAqWCIiIiIOpoIlIiIi4mAqWCIiIiIOpoIlIiIi4mAqWCIiIiIOpoIlIiIi4mBOsQ6WaZoAZGVlWZxERERE5H+d5L8d5Xo5RcHKzs4GICQkxOIkIiIiIv+TnZ1NYGDgde9nmDdazRzIbrdz/Phx/P39MQyjVF4jKyuLkJAQjhw5QkCA69yH0FWPG3TsrnjsrnrcoGN3xWN31eOGsjl20zTJzs4mODgYN7frn6hyijNYbm5u1KtXr0xeKyAgwOX+QwTXPW7QsbvisbvqcYOO3RWP3VWPG0r/2G/kzNV/achdRERExMFUsEREREQczGUKlre3N+PGjcPb29vqKGXKVY8bdOyueOyuetygY3fFY3fV44bycexOMeQuIiIiUpG4zBksERERkbKigiUiIiLiYCpYIiIiIg6mgiUiIiLiYC5RsN5++23q16+Pj48PiYmJbNiwwepIpW7lypX07duX4OBgDMPghx9+sDpSmZkwYQLx8fH4+/tTs2ZNbrvtNtLT062OVeqmTJlCTEzM5YX32rRpw/z5862OZYmXXnoJwzAYOXKk1VFK3bPPPothGL/6ioyMtDpWmTh27Bh/+ctfqFatGr6+vjRv3pyNGzdaHavU1a9f/ze/5oZhMGzYMKujlSqbzcYzzzxDgwYN8PX1pWHDhrzwwgs3fK/A0lbhC9aXX37J6NGjGTduHCkpKbRo0YIePXpw+vRpq6OVqpycHFq0aMHbb79tdZQyt2LFCoYNG8a6detYtGgRhYWFdO/enZycHKujlap69erx0ksvsWnTJjZu3EiXLl3o168fO3futDpamUpOTmbq1KnExMRYHaXMREdHc+LEictfq1evtjpSqTt//jzt2rXD09OT+fPns2vXLl599VWqVKlidbRSl5yc/Ktf70WLFgFw5513WpysdL388stMmTKFt956i9TUVF5++WX+/e9/M3nyZKujXZlZwSUkJJjDhg27/L3NZjODg4PNCRMmWJiqbAHm999/b3UMy5w+fdoEzBUrVlgdpcxVqVLFfP/9962OUWays7PNxo0bm4sWLTI7depkjhgxwupIpW7cuHFmixYtrI5R5p588kmzffv2VsdwCiNGjDAbNmxo2u12q6OUqt69e5uDBg361bY//elP5oABAyxK9Psq9BmsgoICNm3aRLdu3S5vc3Nzo1u3bqxdu9bCZFKWMjMzAahatarFScqOzWbjiy++ICcnhzZt2lgdp8wMGzaM3r17/+r3vCvYs2cPwcHBhIeHM2DAAA4fPmx1pFI3e/Zs4uLiuPPOO6lZsyatWrXivffeszpWmSsoKODTTz9l0KBBGIZhdZxS1bZtW5YsWcLu3bsB2Lp1K6tXr6Znz54WJ7syp7jZc2k5e/YsNpuNWrVq/Wp7rVq1SEtLsyiVlCW73c7IkSNp164dzZo1szpOqdu+fTtt2rTh0qVLVK5cme+//56mTZtaHatMfPHFF6SkpJCcnGx1lDKVmJjIRx99REREBCdOnOC5556jQ4cO7NixA39/f6vjlZr9+/czZcoURo8ezdNPP01ycjLDhw/Hy8uLgQMHWh2vzPzwww9cuHCB+++/3+oopW7s2LFkZWURGRmJu7s7NpuN8ePHM2DAAKujXVGFLlgiw4YNY8eOHS4xkwIQERHBli1byMzM5JtvvmHgwIGsWLGiwpesI0eOMGLECBYtWoSPj4/VccrUL//2HhMTQ2JiImFhYXz11VcMHjzYwmSly263ExcXx4svvghAq1at2LFjB++++65LFawPPviAnj17EhwcbHWUUvfVV1/x2WefMXPmTKKjo9myZQsjR44kODjYKX/NK3TBql69Ou7u7pw6depX20+dOkXt2rUtSiVl5dFHH2Xu3LmsXLmSevXqWR2nTHh5edGoUSMAWrduTXJyMm+88QZTp061OFnp2rRpE6dPnyY2NvbyNpvNxsqVK3nrrbfIz8/H3d3dwoRlJygoiCZNmrB3716ro5SqOnXq/OYvDlFRUXz77bcWJSp7hw4dYvHixXz33XdWRykTTzzxBGPHjuWee+4BoHnz5hw6dIgJEyY4ZcGq0DNYXl5etG7dmiVLllzeZrfbWbJkiUvNpbga0zR59NFH+f7771m6dCkNGjSwOpJl7HY7+fn5VscodV27dmX79u1s2bLl8ldcXBwDBgxgy5YtLlOuAC5evMi+ffuoU6eO1VFKVbt27X6z/Mru3bsJCwuzKFHZmz59OjVr1qR3795WRykTubm5uLn9ura4u7tjt9stSvT7KvQZLIDRo0czcOBA4uLiSEhI4PXXXycnJ4cHHnjA6mil6uLFi7/6G+yBAwfYsmULVatWJTQ01MJkpW/YsGHMnDmTWbNm4e/vz8mTJwEIDAzE19fX4nSl56mnnqJnz56EhoaSnZ3NzJkzWb58OQsXLrQ6Wqnz9/f/zYydn58f1apVq/Czd48//jh9+/YlLCyM48ePM27cONzd3bn33nutjlaqRo0aRdu2bXnxxRe566672LBhA9OmTWPatGlWRysTdrud6dOnM3DgQDw8Kvz/ygHo27cv48ePJzQ0lOjoaDZv3sykSZMYNGiQ1dGuzOqPMZaFyZMnm6GhoaaXl5eZkJBgrlu3zupIpW7ZsmUm8JuvgQMHWh2t1F3puAFz+vTpVkcrVYMGDTLDwsJMLy8vs0aNGmbXrl3Nn376yepYlnGVZRruvvtus06dOqaXl5dZt25d8+677zb37t1rdawyMWfOHLNZs2amt7e3GRkZaU6bNs3qSGVm4cKFJmCmp6dbHaXMZGVlmSNGjDBDQ0NNHx8fMzw83PzHP/5h5ufnWx3tigzTdNIlUEVERETKqQo9gyUiIiJiBRUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxMBUsEREREQdTwRIRERFxsP8H2+5vQN5/nHYAAAAASUVORK5CYII=",
"text/html": [
"\n",
" \n",
"
\n",
" Figure\n",
"
\n",
"
\n",
"
\n",
"
\n",
" Figure\n",
"
\n",
"
\n",
"
\n",
"
\n",
" Figure\n",
"
\n",
"
\n",
"
\n",
"
\n",
" Figure\n",
"
\n",
"
\n",
"