diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..eb3dd6f --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,47 @@ +name: ci + +on: + push: + branches: + - main + pull_request: + branches: + - main + workflow_dispatch: + +concurrency: + # cancel previous runs of this workflow (only for the same pull request, not for main) + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.run_id }} + cancel-in-progress: true + +jobs: + build: + runs-on: ${{ matrix.os }} + + strategy: + fail-fast: false + matrix: + python-version: [3.9] + # os: [ubuntu-latest, macOS-latest, windows-latest] + os: [ubuntu-latest] + + steps: + - uses: actions/checkout@v3 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + cache: "pip" + cache-dependency-path: "pyproject.toml" + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e ".[dev]" + + - name: List installed dependencies + run: python -m pip list + + - name: Run CI checks + run: task ci diff --git a/.github/workflows/publish-doc.yml b/.github/workflows/publish-doc.yml new file mode 100644 index 0000000..ec0dd2f --- /dev/null +++ b/.github/workflows/publish-doc.yml @@ -0,0 +1,43 @@ +name: publish-doc + +on: + push: + branches: + - main + workflow_dispatch: + +permissions: + contents: write + +concurrency: + # cancel previous runs of this workflow + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + deploy: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + + - name: Set up Python 3.10 + uses: actions/setup-python@v4 + with: + python-version: "3.10" + cache: "pip" + cache-dependency-path: "pyproject.toml" + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e ".[dev]" + + - name: List installed dependencies + run: python -m pip list + + - name: Build Docs + run: mkdocs build + + - name: Publish website + run: mkdocs gh-deploy --force diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..9e7130c --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,79 @@ +# Contributing to floquet + +We welcome your contribution, whether it be additional cost functions, new functionality, better documentation, etc. + +## Requirements + +The project was written using Python 3.10+, you must have a compatible version of Python (i.e. >= 3.10) installed on your computer. + +## Setup + +Clone the repository: + +```shell +git clone git@github.com:dkweiss31/floquet.git +cd floquet +``` + +It is good practice to use a virtual environment to install the dependencies, such as conda. Once this environment has been activated, you can run + +```shell +pip install -e . +``` + +to install the package and its dependencies. As a developer you also need to install the developer dependencies: + +```shell +pip install -e ".[dev]" +``` + +## Code style + +This project follows PEP8 and uses automatic formatting and linting tools to ensure that the code is compliant. + +## Workflow + +### Before submitting a pull request (run all tasks) + +Run all tasks before each commit: + +```shell +task all +``` + +### Build the documentation + +The documentation is built using [MkDocs](https://www.mkdocs.org/) and the [Material for MkDocs](https://squidfunk.github.io/mkdocs-material/) theme. MkDocs generates a static website based on the markdown files in the `docs/` directory. + +To preview the changes to the documentation as you edit the docstrings or the markdown files in `docs/`, we recommend starting a live preview server, which will automatically rebuild the website upon modifications: + +```shell +task docserve +``` + +Open in your web browser to preview the documentation website. + +You can build the static documentation website locally with: + +```shell +task docbuild +``` + +This will create a `site/` directory with the contents of the documentation website. You can then simply open `site/index.html` in your web browser to view the documentation website. + +### Run specific tasks + +You can also execute tasks individually: + +```shell +> task --list +lint lint the code (ruff) +format auto-format the code (ruff) +codespell check for misspellings (codespell) +clean clean the code (ruff + codespell) +test run the unit tests suite (pytest) +docbuild build the documentation website +docserve preview documentation website with hot-reloading +all run all tasks before a commit (ruff + codespell + pytest) +ci run all the CI checks +``` diff --git a/docs/_static/custom_css.css b/docs/_static/custom_css.css new file mode 100644 index 0000000..cd62b6e --- /dev/null +++ b/docs/_static/custom_css.css @@ -0,0 +1,167 @@ +/* Fix /page#foo going to the top of the viewport and being hidden by the navbar */ +html { + scroll-padding-top: 50px; +} + +/* Fit the Twitter handle alongside the GitHub one in the top right. */ + +div.md-header__source { + width: revert; + max-width: revert; +} + +a.md-source { + display: inline-block; +} + +.md-source__repository { + max-width: 100%; +} + +/* Emphasise sections of nav on left hand side */ + +nav.md-nav { + padding-left: 5px; +} + +nav.md-nav--secondary { + border-left: revert !important; +} + +.md-nav__title { + font-size: 0.9rem; +} + +.md-nav__item--section > .md-nav__link { + font-size: 0.9rem; +} + +/* Indent autogenerated documentation */ + +div.doc-contents { + padding-left: 25px; + border-left: 4px solid rgba(230, 230, 230); +} + +/* Increase visibility of splitters "---" */ + +[data-md-color-scheme="default"] .md-typeset hr { + border-bottom-color: rgb(0, 0, 0); + border-bottom-width: 1pt; +} + +[data-md-color-scheme="slate"] .md-typeset hr { + border-bottom-color: rgb(230, 230, 230); +} + +/* More space at the bottom of the page */ + +.md-main__inner { + margin-bottom: 1.5rem; +} + +/* Remove prev/next footer buttons */ + +.md-footer__inner { + display: none; +} + +/* Change font sizes */ + +html { + /* Decrease font size for overall webpage + Down from 137.5% which is the Material default */ + font-size: 110%; +} + +.md-typeset .admonition { + /* Increase font size in admonitions */ + font-size: 100% !important; +} + +.md-typeset details { + /* Increase font size in details */ + font-size: 100% !important; +} + +.md-typeset h1 { + font-size: 1.6rem; +} + +.md-typeset h2 { + font-size: 1.5rem; +} + +.md-typeset h3 { + font-size: 1.3rem; +} + +.md-typeset h4 { + font-size: 1.1rem; +} + +.md-typeset h5 { + font-size: 0.9rem; +} + +.md-typeset h6 { + font-size: 0.8rem; +} + +/* Bugfix: remove the superfluous parts generated when doing: + +??? Blah + + ::: library.something +*/ + +.md-typeset details .mkdocstrings > h4 { + display: none; +} + +.md-typeset details .mkdocstrings > h5 { + display: none; +} + +/* Change default colours for tags */ + +[data-md-color-scheme="default"] { + --md-typeset-a-color: rgb(0, 189, 164) !important; +} +[data-md-color-scheme="slate"] { + --md-typeset-a-color: rgb(0, 189, 164) !important; +} + +/* Highlight functions, classes etc. type signatures. Really helps to make clear where + one item ends and another begins. */ + +[data-md-color-scheme="default"] { + --doc-heading-color: #DDD; + --doc-heading-border-color: #CCC; + --doc-heading-color-alt: #F0F0F0; +} +[data-md-color-scheme="slate"] { + --doc-heading-color: rgb(25,25,33); + --doc-heading-border-color: rgb(25,25,33); + --doc-heading-color-alt: rgb(33,33,44); + --md-code-bg-color: rgb(38,38,50); +} + +h4.doc-heading { + /* NOT var(--md-code-bg-color) as that's not visually distinct from other code blocks.*/ + background-color: var(--doc-heading-color); + border: solid var(--doc-heading-border-color); + border-width: 1.5pt; + border-radius: 2pt; + padding: 0pt 5pt 2pt 5pt; +} +h5.doc-heading, h6.heading { + background-color: var(--doc-heading-color-alt); + border-radius: 2pt; + padding: 0pt 5pt 2pt 5pt; +} + +/* Make errors in notebooks have scrolling */ +.output_error > pre { + overflow: auto; +} diff --git a/docs/examples/transmon.ipynb b/docs/examples/transmon.ipynb new file mode 100644 index 0000000..96bbcef --- /dev/null +++ b/docs/examples/transmon.ipynb @@ -0,0 +1,346 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "# Transmon Floquet analysis\n", + "\n", + "Here we perform a Floquet simulation of a transmon subject to off-resonant drives. We extract as a function of drive strength and drive frequency the probability of ionization for the two qubit states. We then compare the 2D plots (with the x and y axes drive frequency and induced ac stark shift on the qubit, respectively and the z axis ionization probability) with Blais-style branch crossing plots, which identify the states the qubit states leak to." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "68ebc2fd", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import qutip as qt\n", + "import scqubits as scq\n", + "\n", + "from cycler import cycler\n", + "\n", + "import floquet as ft\n", + "\n", + "color_cycler = plt.rcParams['axes.prop_cycle']\n", + "ls_cycler = cycler(ls=['-', '--', '-.', ':'])\n", + "alpha_cycler = cycler(alpha=[1.0, 0.6, 0.2])\n", + "color_ls_alpha_cycler = alpha_cycler * ls_cycler * color_cycler" + ] + }, + { + "cell_type": "markdown", + "id": "f5c02996", + "metadata": {}, + "source": [ + "## Define simulation parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "fde5cfeb", + "metadata": {}, + "outputs": [], + "source": [ + "filepath = ft.generate_file_path('h5py', 'transmon_floquet', 'out')\n", + "\n", + "# drive frequencies to scan over\n", + "omega_d_linspace = 2.0 * np.pi * np.linspace(7.5, 10.0, 120)\n", + "# induced ac stark shifts to scan over\n", + "chi_ac_linspace = np.linspace(0.0, 0.1, 59)\n", + "\n", + "num_states = 20\n", + "qubit_params = {'EJ': 20.0, 'EC': 0.2, 'ng': 0.25, 'ncut': 41}\n", + "tmon = scq.Transmon(**qubit_params, truncated_dim=num_states)\n", + "state_indices = [0, 1] # get data for ground and first excited states\n", + "\n", + "# express operators in eigenbasis of transmon\n", + "hilbert_space = scq.HilbertSpace([tmon])\n", + "hilbert_space.generate_lookup()\n", + "evals = hilbert_space['evals'][0][0:num_states]\n", + "H0 = 2.0 * np.pi * qt.Qobj(np.diag(evals - evals[0]))\n", + "H1 = hilbert_space.op_in_dressed_eigenbasis(tmon.n_operator)\n", + "\n", + "# to achieve same range of chi_ac for the various drive frequencies,\n", + "# need to drive at different strengths\n", + "chi_to_amp = ft.ChiacToAmp(H0, H1, state_indices, omega_d_linspace)\n", + "# resulting array has shape (a,w), where a is amplitude and w is frequency\n", + "amp_linspace = chi_to_amp.amplitudes_for_omega_d(chi_ac_linspace)\n", + "\n", + "options = ft.Options(\n", + " fit_range_fraction=0.5, # split the fit into segments based on this fraction\n", + " floquet_sampling_time_fraction=0.0, # what part of the period to look at\n", + " fit_cutoff=4, # polynomial cutoff\n", + " overlap_cutoff=0.8, # cutoff for the excluding from the fit\n", + " nsteps=30_000, # qutip integration parameter\n", + " num_cpus=6,\n", + " save_floquet_mode_data=False, # don't save floquet modes\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "7fb27b941602401d91542211134fc71a", + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running floquet simulation with parameters: \n", + "fit_range_fraction: 0.5\n", + "floquet_sampling_time_fraction: 0.0\n", + "fit_cutoff: 4\n", + "overlap_cutoff: 0.8\n", + "nsteps: 30000\n", + "num_cpus: 6\n", + "save_floquet_mode_data: False\n", + "state_indices: [0, 1]\n", + "omega_d_linspace: [47.1238898 47.2558895 47.38788919 47.51988888 47.65188857 47.78388826\n", + " 47.91588795 48.04788764 48.17988733 48.31188703 48.44388672 48.57588641\n", + " 48.7078861 48.83988579 48.97188548 49.10388517 49.23588487 49.36788456\n", + " 49.49988425 49.63188394 49.76388363 49.89588332 50.02788301 50.1598827\n", + " 50.2918824 50.42388209 50.55588178 50.68788147 50.81988116 50.95188085\n", + " 51.08388054 51.21588023 51.34787993 51.47987962 51.61187931 51.743879\n", + " 51.87587869 52.00787838 52.13987807 52.27187777 52.40387746 52.53587715\n", + " 52.66787684 52.79987653 52.93187622 53.06387591 53.1958756 53.3278753\n", + " 53.45987499 53.59187468 53.72387437 53.85587406 53.98787375 54.11987344\n", + " 54.25187314 54.38387283 54.51587252 54.64787221 54.7798719 54.91187159\n", + " 55.04387128 55.17587097 55.30787067 55.43987036 55.57187005 55.70386974\n", + " 55.83586943 55.96786912 56.09986881 56.23186851 56.3638682 56.49586789\n", + " 56.62786758 56.75986727 56.89186696 57.02386665 57.15586634 57.28786604\n", + " 57.41986573 57.55186542 57.68386511 57.8158648 57.94786449 58.07986418\n", + " 58.21186388 58.34386357 58.47586326 58.60786295 58.73986264 58.87186233\n", + " 59.00386202 59.13586171 59.26786141 59.3998611 59.53186079 59.66386048\n", + " 59.79586017 59.92785986 60.05985955 60.19185925 60.32385894 60.45585863\n", + " 60.58785832 60.71985801 60.8518577 60.98385739 61.11585708 61.24785678\n", + " 61.37985647 61.51185616 61.64385585 61.77585554 61.90785523 62.03985492\n", + " 62.17185462 62.30385431 62.435854 62.56785369 62.69985338 62.83185307]\n", + "amp_linspace: [[ 0. 0. 0. ... 0. 0.\n", + " 0. ]\n", + " [ 1.12426412 1.1338659 1.14344842 ... 2.14500354 2.15300056\n", + " 2.16098688]\n", + " [ 1.58994957 1.60352853 1.61708027 ... 3.0334931 3.04480259\n", + " 3.05609696]\n", + " ...\n", + " [ 8.4132223 8.48507542 8.55678449 ... 16.05173669 16.1115809\n", + " 16.17134508]\n", + " [ 8.48800797 8.5604998 8.63284629 ... 16.1944216 16.25479777\n", + " 16.3150932 ]\n", + " [ 8.56214045 8.63526541 8.70824376 ... 16.33586028 16.39676376\n", + " 16.4575858 ]]\n", + "num_states: 20\n", + "num_fit_ranges: 2\n", + "num_amp_pts_per_range: 29\n", + "exp_pair_map: {0: (0, 1), 1: (1, 0), 2: (0, 2), 3: (1, 1), 4: (2, 0), 5: (0, 3), 6: (1, 2), 7: (2, 1), 8: (3, 0), 9: (1, 3), 10: (2, 2), 11: (3, 1)}\n", + "calculating for amp_range_idx=0\n", + "calculating for amp_range_idx=1\n", + "finished floquet main in 0.5436457713445028 minutes\n" + ] + } + ], + "source": [ + "floquet_analysis = ft.floquet_analysis(\n", + " H0,\n", + " H1,\n", + " state_indices=state_indices,\n", + " omega_d_linspace=omega_d_linspace,\n", + " amp_linspace=amp_linspace,\n", + " options=options,\n", + ")\n", + "floquet_analysis.run(filepath=filepath)" + ] + }, + { + "cell_type": "markdown", + "id": "3ddf2651", + "metadata": {}, + "source": [ + "## Plot the probability of nonlinear transitions" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "acae54e37e7d407bbb7b55eff062a284", + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "data_dict, param_dict = ft.extract_info_from_h5(filepath)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "c8d0ce12", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n \n \n \n \n 2024-09-27T14:56:53.984268\n image/svg+xml\n \n \n Matplotlib v3.9.2, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgNTcyLjAwMDYyNSAyMzQuMDI1NjI1IF0gL0NvbnRlbnRzIDkgMCBSIC9Bbm5vdHMgMTAgMCBSID4+CmVuZG9iago5IDAgb2JqCjw8IC9MZW5ndGggMTIgMCBSIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nNVXS3PbNhC+41fg6By0wi7exziu3eTSuta0B08OqqLI9shxJTf1TH99FyApgg9R6jEHzpALYL9vsU/Or9b/PK7Wv91cyg93Yt5+rV4Fyid+NlLJJ37eJMobfjZC8dezsJ5AKeXI8ue2/CRtQJFNr1ve3Pl8EOKrmL9nNa987EYIpyBqo72VxkDgPaxZedAd2baUkYpgdBK2Z1tZBtjJgVpjHBiJzoBDuV/LP+Q3OX9PiQfbyc9b4iO7N7BjebIWh/pWz7XKvF6pTa8sn39EefUib8Wt3DUYiu8x4SgINRJLhCPwsWt8IULQ2XZxyW54EzuREGaJjdVAJpkdwDOiuFyI+TVKJLn4mt2z+CLu5YUH+05+lotP4qeFuBUZW6AJoG0XtJQdR0V2ZAwewznIMQyhSXtwPWeXsuPQRAbInQUdwMQhtCYHEbvQpew4tEYN1pwJHUegDTrAnpdL2XFooxy4M/wcx/xslQXj+5nVyiZgIwGqsyxGBaqELl1mFGDlWN6lc8okNUhgHSptkUuCnCGE/GaTxvk1pc02gVzs38nFkzDgyStO/8DYBMapmFTmzcyn2jqvt1qrPVYbMsMLygsYwHpD2aBKLhvdLijn+J6pXrjPC5ZJUQNaLdxUmhAcaVuTr1Z+blbIR4quBfk3L3gIvnfic17gC9txPVBoc+E4vEFtYqok1zS7Wj8tf/9+t/z2Ovvlz+3j7vt69vK83iyrEjPUoNmxzN7qmK06re2vx7pa1YFTlDoHRLnUc3EtJKNpwjFmcvnnyHLW24mo6QVNHxbJc+R0cGvRaGnwwEFDKapOgY5kSdlCAijTQa1F0+aS0pxTJ6CxmyRJwyzpwsjkOd2cro8HjhYXrUbW3M8NbHNjmSPIZaLk0GOsQ2tVL2DQqKxxIRxNmh87B1YPh7BNbfZ/93u+yDd5Ly3XjHRD7B4lv3BQWIXRsdFGD6YBvuxBuy5k7Qwi7srWT/JT0fa7Kkemn7GhhhUOh6Ln8aGI9545U3V2NgqOa1XJjnqMqcalTScPnWsaOaHNx6qQzm4tC3wVfr/mAIicYNagqprBMFSbgMfj4cI5pGzQht0WTYOVAwfLwFnu9y9v+8fNw99N3NRjaDniqnrEFZa7vqp9rVMI8b1wYQulaFuKiDu6ybL25EFWT70np03utqEaKHmaDjKp0i6hRWMbqz5Sxb8YLJvo2pSEfAB0McSKugXdF48OljwINucNaK/5glNTLYubmixu9fDR8kBu+J5ixNAlUshPMEEuw1FrZSii1ZNUaJqK4V4RIl9mj0orP0WFfWGiVl4Rl+lJKnqaSuBXE6PXPSqt/BQV9qSyWkXvdDCTVMwklTpKOzQOkTtJgdiFDnOvmoLvtN1cDQ8/Qt0EGM067u/maNZNJmJ7cjI5B6A5X7mEi/8AE0YerwplbmRzdHJlYW0KZW5kb2JqCjEyIDAgb2JqCjEwNjcKZW5kb2JqCjEwIDAgb2JqClsgXQplbmRvYmoKMTkgMCBvYmoKPDwgL0xlbmd0aCAxNzggL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPY87jkQxCARzn4IjGGg+Ps9brSbw3j9dYDSTuEtgQ9m30iaJOnwzhRz64cV+Bv+Wnhi6C65DxqBgJcMhT6FnWSq58Qzw7ZOGqE4TPMjCCBtkAlLf1UEIKXe1th8l3TH5LK5pTXexeGndNXb3K3XXa4kJObIE+RzSQFVV9EPxrkFqSSahboNbCCUM1Nr6FtInDSPU1G+yP7BLu5JPd1CkORN6pyZPPl+LNvr9BzTbPqYKZW5kc3RyZWFtCmVuZG9iagoyMCAwIG9iago8PCAvVHlwZSAvWE9iamVjdCAvU3VidHlwZSAvRm9ybSAvQkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0KL0xlbmd0aCAxNTggL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPdC5DcMwDAXQnlP8BQTwpjWPgCCFs38byo5T6UH4PCQyZgxRwYc0HRk4KRVDuVrD7aE2pY6WzAPOtxjuE5ICL4bIfS466kJHWBG+G4hlK3cr5ZZBsxAS0Op0t1hk0rflnQnJX13q08EsrlGtmhi5V7CeNIQnvAuGWC+UPyzyyTdP8srnJe6yuTexqE2HyU5W4P8fi970+gJl+jY3CmVuZHN0cmVhbQplbmRvYmoKMjEgMCBvYmoKPDwgL1R5cGUgL1hPYmplY3QgL1N1YnR5cGUgL0Zvcm0gL0JCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdCi9MZW5ndGggMTgzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDVQOw5DMQjbOYUvUCn8yXmeVHVo77+W8NohwokxsSFhw4MVH8oaoAuliIJk4KISqAtYHGYK3gG3bEK2DnqT/Lmy7mVwrhZr97C1UvuEIQSiecpFuuaubmDb0EpoD7nIWnrQm0xvznwhGM51S32vA8JWT26fY/T8FeOUEZnjJhrfTtPq5zQ74+GyfFJl8i9m9qw6Y3OW4CxT25HovZZeirZq8yTqDL7nYa3h/3u86EXPL2UCP5AKZW5kc3RyZWFtCmVuZG9iagoyMiAwIG9iago8PCAvVHlwZSAvWE9iamVjdCAvU3VidHlwZSAvRm9ybSAvQkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0KL0xlbmd0aCAxNjYgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNZDBDQMgDAP/TOEFkBKSEJinUtVHu/+3DqWvMw6yA20rwhOfNjWOerfYCx5HjX2VT4MOLRUbO+FTsISYSMGjeQ7Q8mUIBokWHi1knrOOgyFgHO2RJcJoL8ILZSv6rGLZ6Mo4NnU1xgYZt6frr78rF/V1yIn8lHGFnrDkxkbf5gALCrQsnbbWfcv7uGFxlfqEVH8cqPsd5F/c/3q3V3t+Acx1O5QKZW5kc3RyZWFtCmVuZG9iagoyMyAwIG9iago8PCAvTGVuZ3RoIDE2NiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9kEsOAzEIQ/c5hY8QMORznqmqLqb339ZkpG5iB+TkQZihg6kjBpGj42UtYuga+LZguYVw6XSEhTRxtei6zQluRecG5UuvRt/H+VJFD3qa8qGOW4JrwbaD3LCl7smYF8Ld+IiFI2PKOf9O4Rh+agLahI+BNCH0jmR9QGoGfUch1yzRx1HhKn0qAskUIEtXdaKcULOy+exB/vqv4W6f9v4BkwU4JgplbmRzdHJlYW0KZW5kb2JqCjE3IDAgb2JqCjw8IC9UeXBlIC9Gb250IC9CYXNlRm9udCAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRmlyc3RDaGFyIDAKL0xhc3RDaGFyIDI1NSAvRm9udERlc2NyaXB0b3IgMTYgMCBSIC9TdWJ0eXBlIC9UeXBlMwovTmFtZSAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdCi9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdIC9DaGFyUHJvY3MgMTggMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcgL0RpZmZlcmVuY2VzIFsgODAgL1AgMTE0IC9yIF0gPj4KL1dpZHRocyAxNSAwIFIgPj4KZW5kb2JqCjE2IDAgb2JqCjw8IC9UeXBlIC9Gb250RGVzY3JpcHRvciAvRm9udE5hbWUgL0dDV1hEVitEZWphVnVTYW5zLU9ibGlxdWUgL0ZsYWdzIDk2Ci9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM1MCA+PgplbmRvYmoKMTUgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTggMCBvYmoKPDwgL1AgMTkgMCBSIC9yIDIzIDAgUiA+PgplbmRvYmoKMjggMCBvYmoKPDwgL0xlbmd0aCAyNDcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicTVFJbsQwDLv7FfzAAJasxXlPikEP7f+vJR0U7cEQI0tc4u7ERBZetlDXQofjw0ZeCZuB74PWnPgaseI/2kaklT9UWyATMVEkdFE3GvdIN7wK0X6kgleq91jzEXcrzVs6drG/98G05pEqq0I85Ngc2Uha10TR8T203nNDdMoggT43IQdEaY5ehaS/9sN1bTS7tTazJ6qDR6aE8kmzGprTKWbIbKjHbSpWMgo3qoyK+1RGWg/yNs4ygJPjhDJaT3asJqL81CeXkBcTccIuOzsWYhMLG4e0H5U+sfx86834m2mtpZBxQSI0xaXfZ7zH53j/AJVPXCYKZW5kc3RyZWFtCmVuZG9iagoyOSAwIG9iago8PCAvTGVuZ3RoIDc5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM3NVIwULC0ABJmpiYK5kaWCimGXEA+iJXLZWhpDmblgFkmxgZAlqmpKRILIgvTC2HB5GC0sYk51AQECyQHtjYHZlsOVwZXGgDWlBwMCmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0xlbmd0aCAzMDcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPZJLbgMxDEP3PoUuEMD62Z7zpCi6mN5/2ycl6Yoc2RZFapa6TFlTHpA0k4R/6fBwsZ3yO2zPZmbgWqKXieWU59AVYu6ifNnMRl1ZJ8XqhGY6t+hRORcHNk2qn6sspd0ueA7XJp5b9hE/vNCgHtQ1Lgk3dFejZSk0Y6r7f9J7/Iwy4GpMXWxSq3sfPF5EVejoB0eJImOXF+fjQQnpSsJoWoiVd0UDQe7ytMp7Ce7b3mrIsgepmM47KWaw63RSLm4XhyEeyPKo8OWj2GtCz/iwKyX0SNiGM3In7mjG5tTI4pD+3o0ES4+uaCHz4K9u1i5gvFM6RWJkTnKsaYtVTvdQFNO5w70MEPVsRUMpc5HV6l/DzgtrlmwWeEr6BR6j3SZLDlbZ26hO76082dD3H1rXdB8KZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvVHlwZSAvWE9iamVjdCAvU3VidHlwZSAvRm9ybSAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0xlbmd0aCA3NgovRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxFzcENwCAMA8B/pvAIkNSELFT1Aft/S1qq/k6yZYv3gBkxJaWhGMJw1FJSFqjsS40K9Sf1H0bu7FX2jxZ7g+pL38OQS84bxQ0V/wplbmRzdHJlYW0KZW5kb2JqCjMyIDAgb2JqCjw8IC9MZW5ndGggNzMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicM7Y0UDBQsDBT0DU0NlQwsjRWMDczUEgx5AIKgVi5XDCxHDDLzBLEMjQ3Q2LpmhlCZZFYIONyuGAG58DMy+HK4EoDAB6JFpUKZW5kc3RyZWFtCmVuZG9iagozMyAwIG9iago8PCAvTGVuZ3RoIDY5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDO2NFAwULA0V9A1NDZUMDYwUTA3M1BIMeSCMXPBLLBsDhdMHYRlBmIYGZogscyAxoEl4QyQGTlw03K4MrjSAPqpFkUKZW5kc3RyZWFtCmVuZG9iagozNCAwIG9iago8PCAvTGVuZ3RoIDIzMiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UUluxDAMu/sV/MAA1u68J8Wgh/b/11LKFAhAJba4JWJjIwIvMfg5iNz4kjWjJn5nclf8LE+FR8Kt4EkUgZfhXnaCyxvGZT8OMx+8l1bOpMaTDMhFNj08ETLYJRA6MLsGddhm2om+IeGzI1LNRpbT1xL00ioEylO23+mCEm2r+nP7rAtt+9oTTnZ76knlE4jnlqzAZeMVk8VYBj1RuUsxfZDqbKEnobwon4NsPmqIRJcoZ+CJwcEo0A7sue1n4lUhaF3dp21jqEZKx9O/DU1Nkgj5RAlntjTuFv5/z72+1/sPTiFUEQplbmRzdHJlYW0KZW5kb2JqCjM1IDAgb2JqCjw8IC9MZW5ndGggMzk1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1SS27FQAjb5xRcoNLwm895UlXdvPtva0NSqSq8iTHGMH3KkLnlS10ScYXJt16uWzymfC5bWpl5iLuLjSU+ttyX7iG2XXQusTgdR/ILMp0qRKjNqtGh+EKWhQeQTvChC8J9Of7jL4DB17ANuOE9MkGwJOYpQsZuURmaEkERYeeRFaikUJ9Zwt9R7uv3MgVqb4ylC2Mc9Am0BUJtSMQC6kAAROyUVK2QjmckE78V3WdiHGDn0bIBrhlURJZ77MeIqc6ojLxExD5PTfoolkwtVsZuUxlf/JSM1Hx0BSqpNPKU8tBVs9ALWIl5EvY5/Ej459ZsIYY6btbyieUfM8UyEs5gSzlgoZfjR+DbWXURrh25uM50gR+V1nBMtOt+yPVP/nTbWs11vHIIokDlTUHwuw6uRrHExDI+nY0peqIssBqavEYzwWEQEdb3w8gDGv1yvBA0p2sitFgim7ViRI2KbHM9vQTWTO/FOdbDE8Js753WobIzMyohgtq6hmrrQHazvvNwtp8/M+iibQplbmRzdHJlYW0KZW5kb2JqCjM2IDAgb2JqCjw8IC9MZW5ndGggMjQ5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE1RSYoDMAy75xX6QCFek7ynQ5lD5//Xyg6FOQQJr5KTlphYCw8xhB8sPfiRIXM3/Rt+otm7WXqSydn/mOciU1H4UqguYkJdiBvPoRHwPaFrElmxvfE5LKOZc74HH4W4BDOhAWN9STK5qOaVIRNODHUcDlqkwrhrYsPiWtE8jdxu+0ZmZSaEDY9kQtwYgIgg6wKyGCyUNjYTMlnOA+0NyQ1aYNepG1GLgiuU1gl0olbEqszgs+bWdjdDLfLgqH3x+mhWl2CF0Uv1WHhfhT6YqZl27pJCeuFNOyLMHgqkMjstK7V7xOpugfo/y1Lw/cn3+B2vD838XJwKZW5kc3RyZWFtCmVuZG9iagozNyAwIG9iago8PCAvTGVuZ3RoIDk0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWNwRHAIAgE/1RBCQoK2k8mk4f2/40QMnxg5w7uhAULtnlGHwWVJl4VWAdKY9xQj0C94XItydwFD3Anf9rQVJyW03dpkUlVKdykEnn/DmcmkKh50WOd9wtj+yM8CmVuZHN0cmVhbQplbmRvYmoKMzggMCBvYmoKPDwgL0xlbmd0aCAzMjIgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVG7bcUwDOw1BRcwIH4lzeMgSJG3f5s72qlI07wfVV4ypVwudckqWWHypUN1iqZ8nmam/A71kOOYHtkhulPWlnsYFpaJeUodsZos93ALNr4AmhJzC/H3CPArgFHARKBu8fcPulkSQBoU/BTomquWWGICDYuFrdkV4lbdKVi4q/h2JLkHCXIxWehTDkWKKbfAfBks2ZFanOtyWQr/bn0CGmGFOOyzi0TgecADTCT+ZIBszz5b7OrqRTZ2hjjp0ICLgJvNJAFBUzirPrhh+2q75ueZKCc4OdavojG+DU7mS1LeV7nHz6BB3vgzPGd3jlAOmlAI9N0CIIfdwEaEPrXPwC4Dtkm7d2NK+ZxkKb4ENgr2qFMdyvBi7MxWb9j8x+jKZlFskJX10ekOytygE2Ieb2ShW7K2+zcPs33/AV8Ze2QKZW5kc3RyZWFtCmVuZG9iagozOSAwIG9iago8PCAvTGVuZ3RoIDgzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWMuw3AMAhEe6ZgBH4m9j5RlMLevw0QJW64J909XB0JmSluM8NDBp4MLIZdcYH0ljALXEdQjp3so2HVvuoEjfWmUvPvD5Se7KzihusBAkIaZgplbmRzdHJlYW0KZW5kb2JqCjQwIDAgb2JqCjw8IC9MZW5ndGggNTEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicMza0UDBQMDQwB5JGhkCWkYlCiiEXSADEzOWCCeaAWQZAGqI4B64mhyuDKw0A4bQNmAplbmRzdHJlYW0KZW5kb2JqCjQxIDAgb2JqCjw8IC9MZW5ndGggNzAgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicMzM2UzBQsDACEqamhgrmRpYKKYZcQD6IlcsFE8sBs8wszIEsIwuQlhwuQwtjMG1ibKRgZmIGZFkgMSC6MrjSAJiaEwMKZW5kc3RyZWFtCmVuZG9iago0MiAwIG9iago8PCAvTGVuZ3RoIDU0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM2NlcwAEJdSyMFYyDb3MhSIcWQy8jUBMzM5YIJ5nBZGINV5XAZQGmYohyuDK40APuEDh8KZW5kc3RyZWFtCmVuZG9iago0MyAwIG9iago8PCAvTGVuZ3RoIDE4IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM2tFAwgMMUQ640AB3mA1IKZW5kc3RyZWFtCmVuZG9iago0NCAwIG9iago8PCAvTGVuZ3RoIDM0MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UjluBDEM6/0KfSCAbtvv2SBIkfy/DanZFANxdFKUO1pUdsuHhVS17HT5tJXaEjfkd2WFxAnJqxLtUoZIqLxWIdXvmTKvtzVnBMhSpcLkpORxyYI/w6WnC8f5trGv5cgdjx5YFSOhRMAyxcToGpbO7rBmW36WacCPeIScK9Ytx1gFUhvdOO2K96F5LbIGiL2ZlooKHVaJFn5B8aBHjX32GFRYINHtHElwjIlQkYB2gdpIDDl7LHZRH/QzKDET6NobRdxBgSWSmDnFunT03/jQsaD+2Iw3vzoq6VtaWWPSPhvtlMYsMul6WPR089bHgws076L859UMEjRljZLGB63aOYaimVFWeLdDkw3NMcch8w6ewxkJSvo8FL+PJRMdlMjfDg2hf18eo4ycNt4C5qI/bRUHDuKzw165gRVKF2uS9wGpTOiB6f+v8bW+19cfHe2AxgplbmRzdHJlYW0KZW5kb2JqCjQ1IDAgb2JqCjw8IC9MZW5ndGggMjUxIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nC1RSXIDQQi7zyv0hGan32OXK4fk/9cIygcGDYtAdFrioIyfICxXvOWRq2jD3zMxgt8Fh34r121Y5EBUIEljUDWhdvF69B7YcZgJzJPWsAxmrA/8jCnc6MXhMRlnt9dl1BDsXa89mUHJrFzEJRMXTNVhI2cOP5kyLrRzPTcg50ZYl2GQblYaMxKONIVIIYWqm6TOBEESjK5GjTZyFPulL490hlWNqDHscy1tX89NOGvQ7Fis8uSUHl1xLicXL6wc9PU2AxdRaazyQEjA/W4P9XOyk994S+fOFtPje83J8sJUYMWb125ANtXi37yI4/uMr+fn+fwDX2BbiAplbmRzdHJlYW0KZW5kb2JqCjQ2IDAgb2JqCjw8IC9MZW5ndGggNzYgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPYw7DoAwDEP3nMJHaH4kB0KIgd5/pSm0i/30JNvF0WBakQK3wMnkPqnTcs8kO3wQmyHkVxtata7K0poMi5qMvw3f3U3XC6Y4F8AKZW5kc3RyZWFtCmVuZG9iago0NyAwIG9iago8PCAvTGVuZ3RoIDIxNSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UTkOAyEM7PcV/kAkjC94T6Iozf6/zYzRVh7BXIa0lCGZ8lKTqCHlUz56mS6cutzXzGo055a0LXOAuLa8L62SwIlmiIPBaZi4AZo8AUPX0ahRQxce0NSlUyiw3AQ+irduD91jtYGXtiHniSBiKBksQc2pRRMWbc8npDW/Xosb3pft3chTpcaWGIEGAVY4HNfo1/CVPU8m0XQVMtSrNcsYCRNFIjz5jqbVE+taNNIyEtTGEaxqA7w7/TBOAAATccsCZJ9KlLPkxG+x9LMGV/r+AZ9HVJYKZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvVHlwZSAvRm9udCAvQmFzZUZvbnQgL0JNUVFEVitEZWphVnVTYW5zIC9GaXJzdENoYXIgMCAvTGFzdENoYXIgMjU1Ci9Gb250RGVzY3JpcHRvciAyNSAwIFIgL1N1YnR5cGUgL1R5cGUzIC9OYW1lIC9CTVFRRFYrRGVqYVZ1U2FucwovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdCi9DaGFyUHJvY3MgMjcgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcKL0RpZmZlcmVuY2VzIFsgMzIgL3NwYWNlIDQ2IC9wZXJpb2QgL3NsYXNoIC96ZXJvIC9vbmUgL3R3byAvdGhyZWUgL2ZvdXIgL2ZpdmUgNTUgL3NldmVuCi9laWdodCAvbmluZSA3MSAvRyAvSCA5MSAvYnJhY2tldGxlZnQgOTMgL2JyYWNrZXRyaWdodCA5NyAvYSA5OSAvYyAxMjIgL3ogXQo+PgovV2lkdGhzIDI0IDAgUiA+PgplbmRvYmoKMjUgMCBvYmoKPDwgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9Gb250TmFtZSAvQk1RUURWK0RlamFWdVNhbnMgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM0MiA+PgplbmRvYmoKMjQgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMjcgMCBvYmoKPDwgL0cgMjggMCBSIC9IIDI5IDAgUiAvYSAzMCAwIFIgL2JyYWNrZXRsZWZ0IDMyIDAgUiAvYnJhY2tldHJpZ2h0IDMzIDAgUgovYyAzNCAwIFIgL2VpZ2h0IDM1IDAgUiAvZml2ZSAzNiAwIFIgL2ZvdXIgMzcgMCBSIC9uaW5lIDM4IDAgUiAvb25lIDM5IDAgUgovcGVyaW9kIDQwIDAgUiAvc2V2ZW4gNDEgMCBSIC9zbGFzaCA0MiAwIFIgL3NwYWNlIDQzIDAgUiAvdGhyZWUgNDQgMCBSCi90d28gNDUgMCBSIC96IDQ2IDAgUiAvemVybyA0NyAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YyIDE3IDAgUiAvRjEgMjYgMCBSID4+CmVuZG9iago0IDAgb2JqCjw8IC9BMSA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAwIC9jYSAxID4+Ci9BMiA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAxIC9jYSAxID4+ID4+CmVuZG9iago1IDAgb2JqCjw8ID4+CmVuZG9iago2IDAgb2JqCjw8ID4+CmVuZG9iago3IDAgb2JqCjw8IC9JMSAxMyAwIFIgL0kyIDE0IDAgUiAvRjItRGVqYVZ1U2Fucy1PYmxpcXVlLWNoaSAyMCAwIFIKL0YyLURlamFWdVNhbnMtT2JsaXF1ZS1vbWVnYSAyMSAwIFIgL0YyLURlamFWdVNhbnMtT2JsaXF1ZS1waSAyMiAwIFIKL0YxLURlamFWdVNhbnMtYXJyb3dyaWdodCAzMSAwIFIgPj4KZW5kb2JqCjEzIDAgb2JqCjw8IC9UeXBlIC9YT2JqZWN0IC9TdWJ0eXBlIC9JbWFnZSAvV2lkdGggMTIwIC9IZWlnaHQgNTkKL0NvbG9yU3BhY2UgWyAvSW5kZXhlZCAvRGV2aWNlUkdCIDExNQooCDBr9/v/9fn+9Pn+8vf98ff97/X87vX87PP76/P76PH65e/55O/54u344e34CEGF3uv33On22+n22ef12Of11eX00+Pzz+Hyzd/xzN/xyt3wyd3wx9vvxtvvw9nuvtfsvNfrtdPpq8/lqs/lps3jocviCDNvncngmsfggLnae7fZeLXYcrHXba/WZ6vUY6nTXqXRWqPPU53MUJvLR5XIP4/EPY3DOYnBN4fANYW/M4O+J3e4JXW3H2+zGWetFGGoC1WfCEWLCDh39vr+8/j98Pb87/b87fT77PT76vL66fL65/D55vD55O744+744ez34Oz33er22uj11OTz0eLyzuDxyNzvxdruxNruwNjtv9jst9TqttTpstLordDmrNDmp87ko8zjoMrhl8bfCEyWb7DWXaTQSJbIQJDFMoK+LHy7Knq5G2qvFWKpE2CnEl6mEFxcpFxyWKEMVqAIUpwpCl0KL0JpdHNQZXJDb21wb25lbnQgOCAvRmlsdGVyIC9GbGF0ZURlY29kZQovRGVjb2RlUGFybXMgPDwgL1ByZWRpY3RvciAxMCAvQ29sb3JzIDEgL0NvbHVtbnMgMTIwIC9CaXRzUGVyQ29tcG9uZW50IDggPj4KL0xlbmd0aCA0OCAwIFIgPj4Kc3RyZWFtCnicxdhnW9swEADgWPJITEoZpaS0dNJJ995LJ8cZ0JgA3ZPR3dL2/z9P5dhJrLPJcBznPhCME72RdDoJZ5RW3PNemBI3uP/KVPcHI2BM/5vNmZU6YVKjE+5Vpn294t2F2DBtwo0mBKzXbh4uGytHEWwQGS56d+diw63WmzDPn1uo6GcvUhkGKsPg3dVjw62gHkw18935JWvrKYIpyDAl8vceFGZUzR68U9PHtxFMNBlmVEko/B4QKM48WTGKryiRYJZDsJo4rDu/Vm1eAAQfYBKs2EnBzIcpX7r6PKd9RrDyAMH1yFaOd3VCU+TDYpLNjQ+m9lpF8ASCf3cleowmTCD38sIYX+YIthD8LBmwfSmya/X7OM9mEQxEhnnkx/sINOYCtiYzk7o+Q+U2CZVh//ZWbNhAMKPa8rf9hnUNwW5NDcL+gMTv8U4IVs3FeZuXAMG6DDMSm/SCo2sx1tljH0tc1eTSpUzJcHhd9Bn4i7s74+zlCocyyq4/CNZignt3+bsYa2vq5zKHrwjOMxmu9YFN9/CeRnb9dSwwEawi+EcfcE8hxtq8UjcAUHZRBJeHANsnXhSBcpTWRIb9U481KNf+jYFe2MgBeYNgKsPRy6nvdR1ohlE+9cUEUkKNaDKcwOHDjcBxkRE1fyqvEor6tMeDW2AyshmEoXTD0ShBZ9fbCB60dHkROKeKSS6erFuUzcl9sj24VbHin6mDESyAooQUjuiUobT2k6t15jGVJCJoiEl2Fm1gqIIw1oA3m9fbicCSQKD6qCwOXXivbsCthYu3l0D0NP0RuwyD7P0xAd9CzWWkRgdNrrcRMNXPOGqox5AJv3WAiFiOjGhrs5r4tw3tT8nCUd+FQG3BErklr5i5ocNiks27htiO5ANZwkMdCVPjYU7AaD0h2M+uZAqJTxDtU9Xdj2UZwYPuiJEy1B0XXusEF4YC59chlO/RQz1QrGOY2Je0UMPDTy7RG/7YcmGphqQBM3hvEzzWu8DXE4WpU6U9wskGKe1rwEE8FZjx06G0TgeGza5wYo+c5FbntW5w7NLV8YDKKsZIhlpMcnU0sAKHvNLVPlulBJMdXDNTgpmDszYlWOF4i08O7rweKD45p9Vjhh9zpAUri6OC8cPw1OARFZBwpAejLo+sx/8B72SMyAplbmRzdHJlYW0KZW5kb2JqCjQ4IDAgb2JqCjg2NwplbmRvYmoKMTQgMCBvYmoKPDwgL1R5cGUgL1hPYmplY3QgL1N1YnR5cGUgL0ltYWdlIC9XaWR0aCA0MCAvSGVpZ2h0IDI0MAovQ29sb3JTcGFjZSBbIC9JbmRleGVkIC9EZXZpY2VSR0IgMjM4CigIQIL3+//1+f70+f7f6/fy9/3x9/3u9fwIUJvr8/vq8/t8t9ro8frn8frl7/nk7/ni7fjh7fgIQYTe6/fc6fbb6fbZ5/XY5/XW5fTV5fTT4/PS4/PQ4fLP4fIIS5PM3/FfptHK3fDJ3fDH2+/G2+/C2e7B2e0IMm691+y81+u41eq00+mz0+iPwt6v0eeu0ecIRouqz+Wpz+WlzePj7vmhy+Kgy+IIM3AIUZzf7PeZx+CVxd+Rw96Nwd2Kv92HvdwINHEITJV5tdl3tdl0s9hxsddtr9dmq9RkqdNhp9Jgp9IINXNdpdEIR41bo9BYoc/c6vZSncxPm8tNmcpJl8lGlchEk8dBkcYINnQ9jcQ7i8I5icE3h8A1hb8zg74xgb0vf7wIQoUre7pcKXm5J3e4JXW3I3O2XHJXoQg3dh5tshxrsBtprxlnrRhlrL/Y7RZjqhVhqRJdphBbpA5ZoghNlgtVn1xuU54IMW0IOHcISI71+v7z+P7y+P3w9v3v9vzt9Pzs9Pvp8vrn8Prm8Pnj7vjg7PgIOXkIQ4fa6PbZ6PXX5vXW5vTU5PSfyuHQ4vLO4PLN4PHL3vHK3vDI3PDH3O/D2u4IOnq+2OwcarC61uu51uq31Oq11Omy0uiButut0Oar0OaozuSmzuQtfbukzOOjzOMIO3ydyuGayOCXxt+UxN+SxN5/udoISZCMwN3T5POJvtyFvNyEvNs8jMMIPH19uNp6ttl1tNhysthvsNcIRIhsrtZqrtZorNU/j8VlqtRjqNPN3/EIPX9cXKTQWqLPV6DOVqDOU57NUZzMTprLS5jKSpjJSJbIRZTHQpLGCE6YPo7ECEqROorCOIjBNobANIS/MoK+Ln68LHy6Knq5Jna4JHS3InK2IHC0H26zHWyxCEWKGmiuGWatF2SrFWKpFGCoEl6mEVxcpQ9aow5YogxWoFxuVJ4JUp0IMGshcbUpCl0KL0JpdHNQZXJDb21wb25lbnQgOCAvRmlsdGVyIC9GbGF0ZURlY29kZQovRGVjb2RlUGFybXMgPDwgL1ByZWRpY3RvciAxMCAvQ29sb3JzIDEgL0NvbHVtbnMgNDAgL0JpdHNQZXJDb21wb25lbnQgOCA+PgovTGVuZ3RoIDQ5IDAgUiA+PgpzdHJlYW0KeJyNyuMCGAQAgMEvLNvLttdqcdm2jS0s27Zt2za31WqtWlq29Qy9wt3v61/U8ajFUKug1kZthtoZNQx1Auo81HWo+1DPoN5CsZlQ+6DOR72E+gG1PGoL1ImoR1HjUfOi1kEdg/oQNQVqCOof1HGov1HHov5CHYT6E3U06g/UUajfUUeifkP9ijoC9QvqcNTPqMNQP6EORf2IOgR1A+pg1Peo4ajvUN+i/kN9gzoQ9TXqANRXqP1R+6G+RO2L+gJ1N2oCam/UXqjPUXuiPkPtgfoUtTvqE9RuqI9Ru6KeRu2C+gj1Gmon1AeoHVHjUDug3kdtjxqLeg+1Hepd1LaoMahtUO+g3kaNRm2FGoXaEjUStTlqftSmqE1Qb6A2Rr2O2gj1KuoV1MuoDVEvojZAvYBaH/U8aj3UuqjnUFOjnkU9groN9RTqSdRQ1BOotVCPodZELYNaA/Uw6iHU6qgHUauhHkDdj7oUtTJqMOpe1D2oQai7UHeiVkStgLoDdTtqOdSyqFtRS6OWQt2Cuhm1JOom1I2oJVCLo65HjUAtiloEdS1qYdRCqGtQV6MWRC2Augp1JWo+1JuoK1CXo+ZBzY26DDUXak7U46hLUANRc6AuRl2Emh01G+pC1AWoWVGzoLZGzYyaFLUq6lzUjKgZUOegVkJNj5oOdTbqLNS0qGlQZ6KmQk2JOgN1Ompy5PE01KmoyVADUKegTkZNgpoYdRJqIvQ/5V5eAAplbmRzdHJlYW0KZW5kb2JqCjQ5IDAgb2JqCjU1MwplbmRvYmoKMiAwIG9iago8PCAvVHlwZSAvUGFnZXMgL0tpZHMgWyAxMSAwIFIgXSAvQ291bnQgMSA+PgplbmRvYmoKNTAgMCBvYmoKPDwgL0NyZWF0b3IgKE1hdHBsb3RsaWIgdjMuOS4yLCBodHRwczovL21hdHBsb3RsaWIub3JnKQovUHJvZHVjZXIgKE1hdHBsb3RsaWIgcGRmIGJhY2tlbmQgdjMuOS4yKQovQ3JlYXRpb25EYXRlIChEOjIwMjQwOTI3MTQ1NjU0LTA1JzAwJykgPj4KZW5kb2JqCnhyZWYKMCA1MQowMDAwMDAwMDAwIDY1NTM1IGYgCjAwMDAwMDAwMTYgMDAwMDAgbiAKMDAwMDAxNDgwMiAwMDAwMCBuIAowMDAwMDExMzcxIDAwMDAwIG4gCjAwMDAwMTE0MTQgMDAwMDAgbiAKMDAwMDAxMTUxMyAwMDAwMCBuIAowMDAwMDExNTM0IDAwMDAwIG4gCjAwMDAwMTE1NTUgMDAwMDAgbiAKMDAwMDAwMDA2NSAwMDAwMCBuIAowMDAwMDAwMzQ0IDAwMDAwIG4gCjAwMDAwMDE1MDcgMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDAxNDg2IDAwMDAwIG4gCjAwMDAwMTE3MzQgMDAwMDAgbiAKMDAwMDAxMzIzOSAwMDAwMCBuIAowMDAwMDAzNDc3IDAwMDAwIG4gCjAwMDAwMDMyNjIgMDAwMDAgbiAKMDAwMDAwMjkyNiAwMDAwMCBuIAowMDAwMDA0NTMwIDAwMDAwIG4gCjAwMDAwMDE1MjcgMDAwMDAgbiAKMDAwMDAwMTc3OCAwMDAwMCBuIAowMDAwMDAyMDcwIDAwMDAwIG4gCjAwMDAwMDIzODcgMDAwMDAgbiAKMDAwMDAwMjY4NyAwMDAwMCBuIAowMDAwMDEwMDQ0IDAwMDAwIG4gCjAwMDAwMDk4MzcgMDAwMDAgbiAKMDAwMDAwOTM4MyAwMDAwMCBuIAowMDAwMDExMDk3IDAwMDAwIG4gCjAwMDAwMDQ1NzIgMDAwMDAgbiAKMDAwMDAwNDg5MiAwMDAwMCBuIAowMDAwMDA1MDQzIDAwMDAwIG4gCjAwMDAwMDU0MjMgMDAwMDAgbiAKMDAwMDAwNTYzMiAwMDAwMCBuIAowMDAwMDA1Nzc3IDAwMDAwIG4gCjAwMDAwMDU5MTggMDAwMDAgbiAKMDAwMDAwNjIyMyAwMDAwMCBuIAowMDAwMDA2NjkxIDAwMDAwIG4gCjAwMDAwMDcwMTMgMDAwMDAgbiAKMDAwMDAwNzE3OSAwMDAwMCBuIAowMDAwMDA3NTc0IDAwMDAwIG4gCjAwMDAwMDc3MjkgMDAwMDAgbiAKMDAwMDAwNzg1MiAwMDAwMCBuIAowMDAwMDA3OTk0IDAwMDAwIG4gCjAwMDAwMDgxMjAgMDAwMDAgbiAKMDAwMDAwODIxMCAwMDAwMCBuIAowMDAwMDA4NjIzIDAwMDAwIG4gCjAwMDAwMDg5NDcgMDAwMDAgbiAKMDAwMDAwOTA5NSAwMDAwMCBuIAowMDAwMDEzMjE5IDAwMDAwIG4gCjAwMDAwMTQ3ODIgMDAwMDAgbiAKMDAwMDAxNDg2MiAwMDAwMCBuIAp0cmFpbGVyCjw8IC9TaXplIDUxIC9Sb290IDEgMCBSIC9JbmZvIDUwIDAgUiA+PgpzdGFydHhyZWYKMTUwMTkKJSVFT0YK" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "tmon_idx = 1 # Plot leakage probability for the first excited state. Try 0 for the ground state\n", + "omega_d_idx = 100 # Drive frequency where we take a linecut for comparison with Blais plots\n", + "\n", + "plot_data = np.clip(\n", + " 1 - data_dict['displaced_state_overlaps'][:, :, tmon_idx].T ** 2, 0.0, 0.5\n", + ")\n", + "fig, ax = plt.subplots(figsize=(8, 8))\n", + "xticks = param_dict['omega_d_linspace'] / (2.0 * np.pi)\n", + "yticks = chi_ac_linspace\n", + "num_x_pts = len(xticks)\n", + "num_y_pts = len(yticks)\n", + "im = plt.imshow(\n", + " plot_data,\n", + " origin='lower',\n", + " cmap='Blues',\n", + " aspect=0.75, # 2.3, #0.15,\n", + " interpolation='none',\n", + ")\n", + "plt.axvline(omega_d_idx, color=\"grey\", ls=\"--\")\n", + "ax.set_title(f'$P_{tmon_idx}$' + r'$\\rightarrow$', fontsize=15)\n", + "xticklabel_locations = np.linspace(0, num_x_pts - 1, 6, dtype=int)\n", + "ax.set_xticks(xticklabel_locations)\n", + "ax.set_xticklabels(\n", + " np.array(np.around(xticks[xticklabel_locations], decimals=2), dtype=str),\n", + " fontsize=12,\n", + ")\n", + "yticklabel_locations = np.linspace(0, num_y_pts - 1, 3, dtype=int)\n", + "ax.set_yticks(yticklabel_locations)\n", + "ax.set_yticklabels(\n", + " np.array(np.around(yticks[yticklabel_locations], decimals=2), dtype=str),\n", + " fontsize=12,\n", + ")\n", + "ax.set_ylabel(r'$\\chi_{\\rm ac}$ [GHz]', fontsize=15)\n", + "ax.set_xlabel(r'$\\omega_r/2\\pi$ [GHz]', fontsize=15)\n", + "cax = plt.axes([0.91, 0.35, 0.05, 0.3])\n", + "cbar = plt.colorbar(im, cax=cax)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "source": [ + "## Now lets compare with the Blais branch plots" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 29, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n \n \n \n \n 2024-09-27T14:57:31.455992\n image/svg+xml\n \n \n Matplotlib v3.9.2, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgNDI0Ljg5NzUgMjc5Ljk4Mzc1IF0gL0NvbnRlbnRzIDkgMCBSIC9Bbm5vdHMgMTAgMCBSID4+CmVuZG9iago5IDAgb2JqCjw8IC9MZW5ndGggMTIgMCBSIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nK2dS49kSXKd9/krYkkuGO3vx5IjiS3MbsgGtBjMguh5cKgeAuRQIqBfr3PMPOOaWVRmZVb3AE12nY6615+fm5ub+f3uv//h//75xz/84/e/uv23f3r57vrTj399ybd/xT9/uqXbv+Kf/7rl2/f4508vCX/6y0sr7b727PjDT9cfytz3vSr+7Sf8zv7pX15e/vjy3d/jAX/FX/n+5aWOext9lH6r894bf/WXl7LmHc9x6k9WLbPc53nk9QSrnhcVfdGfUGQU/75QAbyWykubLOzay77ZiO2ezotffoUW+K+Xf8f/Tbe/S3gW3jlHT6uXMuutlPvutx//8vKrH16++4d8y+n2wx+lgX74/ctvb3+T7in97e13tx9+/fI/fnj5zYuU42XX+yhjr23fb8T33r/Kve2W58h1tI+8vzy/P9d9r30kV38jvvf+XNa9lFrGTK33jxSgfaEAq9/RhHm6AlziuwWYFa8stc/d24cKMJ4LUCr+f25lZDf4jPpeEUre94EftrVGHR8pwvpCEebgi2svrghGfbcIo99bTr3uOcv8QBGyG4fmSbvghQuv2fcyVy7vPebP//b7//PjH35/++cfb3/9z3/+j/99++u//PmPf4suuWc+4W/+8/bb7//n//vdc2WvmdowfDBy50RlK0od1S92eb33vtkUrd1b3SNh8qHW7zb5e4VAI+NppVZXiEt9vxBgTQds+iodvf92IfCz9wqRU7t3zDrMeFsKI79fjJzSfeaZMWBra++Uo7/fGBnvHr3M1n05Lvkr5QARVu9rppJ7facc8yvtMVCfnXvKvhyX/KVyzHu55T7uezfgM+9R3ilB/sq4yGvdd02l+9Fp5DeLAHjlWusqeb43MPNXBgX+OoC1U/ZD08hvlaCgRzArcl0ZaHyvCF8ZD6X1e0lzDj8ujfxmEVomlBbWp7rWe0UIQ+Hfb18wBkrPmAg3kP5e2u0//nD7X7d/u5Xbr2+QuZjfM8bbWKPMjjq0cf438V9m6mvhL5Z9+8fvb96qMQu9ZRFetlpKQD9VTEBhUe8wK0rB+6kOlI8d0Pd9llbXFHXlnAbUgQV8jK4w25gbIs47ltSJfoNac5qND5ioAmyIOUQtHdOMKpCOYVa3qC0VtCHUDeT1koqo+Gvy21XBqlmHvK2Orc+FmVQzFkepWl0N7Qh1ZxQ9jy1Vq3ujh6iyO+vqUjV02JTRvfc9t5GyVK2VhYdAzaliDVtZnttaUTihNveE8japXOt4QqGMWbZn6Ulq1zgkq8jtvvHsIdVrCz2lMukxVpX6YZZneTQWtIm3bKlfZz9PkWF9tZy7FKTXlAepkLFu4wm1SA17ayN3kRtMKixuUkXQKrUlMv4wJvpd5Flbl1c2TM20lz5jzSWPaBiXPe8pVeS6o8XrCQ2JRpMqom37ULndax1lSxVHxUyQUvfFpRV8ErmhMzk08iiwpNJZ/kbfsOtEHnf8p7mkjmikOaRVJ9fwtruUb6wlPMyTRgLXAKogcKlSkAl7ee0ypYozj6Y1Xxkcg90iak2jqtrvC0MpSR0nfppU3neiZEgdZ19bCg0zFf28i1QRIzBtKd2eeHviHyivqVUsKaO+peikQJEBUZE7Wrg1HU94SZ1b5H3vMGy0C1bZTSpDsDX0vxZ7tSp1IRHBmp2kiguDL7H10EGYJPWoM/ecRe2cklzaKWO1SirD5qkTBoXIWEJ0t1HrPcO67FLHjWFUpdQgE5Azi9QRjVtlABN6bMkldcSgwOwQGZhNg0OOMqZ4Upm/wSohddyYlFrsXvD+VKfUcaPuMvwKFjbMDLSPyHss4U/hkjg5fvc9YRbqk0djT8MiFRnv71LHIWMrlSFyPTOaGyZMc3SyyNjNpCnyw/ykPDDiZcX5p5ff3L4F0lgewLQFpoHVE5hOvWEV29xmvEdmELbtjH+xYIaIfoPV5cBMFfVY2YGZKv6z9P4DzBRrWTs7MFNtoG1zYKbaO2DjwEx1oPc8mKkuNP12YKa6QZnswAx1YIjIGLnATLXChvZgpkoLvzswU11ox+rAzKpjpPbqwUwZ9RjNg7lzbueWtgcz5cWJ6cEMeQPFpXgyYxSi0lkfcpEZMpoz1UBmyBu/ysuTeWRuuXIvnszoM8yvrcx5kHkQpAnT1KN5LkxMmCnNo3nj3Y3c9WjeS3AyPJlzKtjEJKxWHs00x/HuoWvExWbsvYAUQGB4OMMWQPfRfvB4lj4ZQFz3fMZijL4aKJEDNFr1zoVJl9+L0DmDbxhruXlESyeuIcPKMJq9xX7OntFoHYysJvg3jMZjMZrylIY0kM5roOx9SZcaSvOvLgBlLY/pgoEGY5k0c6AuKCQNyD4cqfmLumE8bI9qbrMxLYei8GI1egatXba0gYE1dSwEQiwDa74IA6wLIQytqa+OPtse18AGSo++7Z7X1Aemq7baBWysuRgURJQnNn+P4ZezQzbbJlXQtnhmF0432ESleGizUkAausBTmzqeV0b12Ka+sRkXnH0zt9HAAIWiGtQe3AZiFmCtsv/lXYRPKR+azjGcaq8KxYvhosL+GY7hVIFtZ1tTQ4uJSXQhXNS5kkc41YWN9HAIp4qWTR7hooIWxSEc6oK9LVuBC+GigmXNIZwqRpkMygvhVDGOhUEXwkVdTYH6QDjVSueeRzhljN9UPcJFxujKHuGUe5sj2NaUAQ01ly+Ei7yOPXohnDIMqlQ8wimjoLrGXAgXeenG0CCcMmb16A7hUNH6B78G4Vwd0DTdI1xkdGLzCKeMkV89wqmCENUDnCo4qg1y8RvDBIvOVAv4wjfkiuGp6+VFb8hojbGWhzefjU7QrnmwW1Q0djCuKe+uvXuBm08myroH9+SySGvNgxsyFjoU1YObD8F+twRuU+5Y06rHtshtKIEuaos8dYEy0KaMVVtKcjFbVOyRlkc2ZaCieGBL8dDu1fNa5J0CrqluTPZgW7O/UlNwGFhTxkxUal6splzQTMWjWmRQNtjWlLGH1vXnFdQUURMZqIbTMm6qLrQG05RHyWICGkpTnmm2YFuLvOsYP4/RS8xqFguMxmjGnMmcb/yDOEYwud5l9Coc1jQyLKNFnTC0HaOpiq/RMVpUjBoHaYrlWIAXpEXFlm47SIuK3ex2kKYK5MgguCAtKgyI7CBNFcVK3UFaVJjk20GaagdIhoO0qGj64iBNdWR1t1yQFrWv6AChjH1YtLNF7nsHO5vyysqei9Gicph5RlPeNbg/IBKXrXlCU874t+0JTbnq0mQATZVu0+kADRVTYe3qAY2WxXzbKdjYPI5Kc8pW0wAaMib61gXuIrS8se7aPaKXmt5K1wvRlIHInT2iVe4K9AvRIgPGwf8hz85DdrcXokXFMhkQLTLZ7Rktxcag2J7R2n5DF6eL0aw7mmoGBwhlWq2B0SJ3faVhNGXo4pMzjBZ5KjgMo7W5i3pRHowWFV05PaMpY5+oboAL0pTpaAseEJGBgOkpTXm0LGPHUJoyna/VU1rkeRaii9KUVx0ytg2lKQOq+uuL0iIDvMVSGiImoJo7htKUsRPLw1Na5KlrsKE0ZTbJ9pSmXLN2+8+gdD/uj7Hppk6vzKabWnjdSM/3KI0RVzBalndTUwVum3dTU6UrbDpKU0VDL+cNoYhNzvCUhloBFjWPH5SmClNTLL+L0lQxZsU7dlGaKhp/ezc11ZaDM4Qi/n/3zhCq2BALUy5Ib3mU+jcvSFPdsH28lxpqA9FXsKQp04M7PKQpNxj03UOa8ihD+X9RmjKG+jGZH5SGjNrXGpwhlDGCAqZ51p+xhwqYpoyNVame05BXSkuXiwenoW7MgBEM6d04A2fwUkOlT1HNSeMKGYR61+JdmEbzVry/FY9pNhksQbXYrCukoUmAMO0E4woZ9OLpGmAcIWQzLEdduowjZHFbhddnR+osi347DnbjCOE6NcbKwVedMeRgUOj+znpCoJeGXdv2tKZeMdFa9IUs9nPqB8zGF4IexZxSMl/AFn22Gn0hosPUa8Ujm/pIqanb9cFs/fnUEWOgLTpHXaC26HQyBmxrMbtuIAy3pVr0cw8PbtErB60ntzRPTzMHdEtzJszGwG52C/pTTycMvKVDwW7nC2Hvt3LOuQy9Ycrix+Usihe+Ofg5LNXfffE7N+6nMYOXBzh1bBlK/5m+kFef9VwgeK3q/yj4D8fi5q7kPYC7k/6L4CKD0MGhTXlUNL/3hojMgJNqIU51okLFHzZSxmRSi+DCuMiYwr07jlPeGRvt7UAuMgONHMjp+Ezo1FwcykVeuhW+UE4VvaYxFRfLRd6tVe/ZFp8qRLX4HjQXGY+eAefUq57IGpqL2nNawTGiOs0Nz3PRseWfwbstb8WMVk5dRNeKYhvZPdOlJ0BmPai8oM567DyOtX9RXRzO9LYGrJNMnKxP5jf+a9FCGv9IJ5HLOWZ8cH0QsdgThdNHzDBQqexw+jhY53Gs3gvrAzM2HdfLBfVB//M6p34X1ek5HtxLeqgP8V01dcE/mD7YaNh5BKQP2kqYAdsTnTLerf1wAZ3HDB0byuV5ThnzPgWcDzYZHjM8zVUexyh/wJwPKe14FS6Wa0mKbI0MyokzDMngIUGLFEYTdQ9yHmystpWbF8c353w7J3kG46NLiEIL5jeX1t5608NDg/HBg/iResT44KkE5n3waechlmFSQ9lgXM6NdtOTAIPxjtlXSvbHkCwlrOpSnzheZRMeKY7Owspaa6Q4liyMIG0xQ/EqR1tZVsqf4y15hXh7ChcRtfLM5l2MZ1QTY2V7lzZMCiy13N54jFdxiw7ZZxmM18x1X90KF8Yr9/5bd4cG45AbZ2bAeJVo11K8Zxsb/3viYZ93bVOG7VcDxqECLTXEjVCGkT2HN8kpw0DO4YCSmEsbpqwPHaGMDu/HSfLgeGX0HkZ7CxyHDhOpHVf2RfLKtY2RdoHkVYKNz3MMydH1CYuK0tCQHHpGu5Xg5aZeYaqMYJ1Tx4tSD0eV1FG/nXogOfTJcEDvSKGM7k96JnehnPrmOrE8y3koh2XlOG8umFPHFOx6cmhOLKHDvpvqXDJmOnSsf0mPVI2ZzkO/2pvuDYyZLiHN6ZwkGEMd+sIY1wXWGOqYcgl7+eldKpTRiDlSnTrI0mo8sYQO1OiQt4Y6dJRranXtmSWmBYqtDgpjqEOfFdbBDIY6dA7BFdhOHYPneHuNoY4OzQkrcvB/U8eEnjkY6pUxWr3k4F2h3jBlZvCBy3gB0yPgIU9MsnOWeQG+8jiwrhwBXxny8XomagAPnRbcE+ChwwwcKRxaUl8EeTi0pMnMxhnOy5LxosbdTQ6AxzK0R1l6ZmAQv+To/QTfGMSLx6eN6A+nPiuq8rOCTWyoye37l7c5XriIYQ1ugeMYSGCmFMJwvCzZRMisNxxH16IR1ca6OF6kWevy55OEC8o2xNK3HOepAeAdOY4yoVnk0ZbjDP0famJ7kIMWwf1NFWRTX7nlOBAMw0EKYjneyfG5vQdcwVZ28K4IeEYpIdSE8mZbBh84Jxy2tL2HKECZWDQmIsUxI/aJq7AQ36AdjbIIcQxBWHDqxTYQx0BGs/UUIc65NDGdR4A4dBjpGgZlII5/w3/auvcwEG+Mjl7noNBAvDEYADvmEiAOHaM592CSg5w8e5o9+lqgYwWqvQeIc8qPY/UahkPuPGfageHQ98O3ZRjepAv3Ds4WyBhVykbLcEyeDH0Gxzh1zJequmE49A17cEdniwTBPqIAL4ZLpGPSwyXLcIlkw98N3nHqLFiMD6RO9qbIcHQQqNmSDzyRfivrOE8Mw6HDGt8KTcPwzhPmWnSpMRAHRIAKdXlaiJNXCfDtAeLQ0f0pusnp90gMhhwB4p1+66GRDM5KZyZCPfi1Vjp+UrClHQ7itMbHCSz0RnpNu64QeEId2O16mGsYDn1yySqB4dA3pmH/mWb6fPWvvDgHeZLwev3f+1Y6sbISg0kc3Tl3ufR4uGOaYXRNHzyoDk09XLvQ3uiG32qNGrRzwtObmT3ahQN7yRbSoF1+PdUyNmjns8uhgEE7Ze61ffwJ5YqalOBrYakxSWTpNmyXbVVVl7BhO2W6fIKvhZxaeHa00Rv3raMrHg3coS/OuBbgfhpWgWToTt4xXjoEolDHmNgaNWjwDh3ySk94b3TD5xoCvYWb2Omu4EOnjqF5DmgN3pvE/cpgdnjnVqIfB4qhO/1WXa0sR3f6lpLGCFq4d3nPk4UOnYHdugZ5uK9W1Bh0dB88ZddTHUd3/NVdurqcAt0x5bXxHd07g9VzpPtgUFbRGBRLdyzoKMCThc7VoOtpnoP7ZALhcf1auE/pBaWghfuUXeyTgQ7wJYxd9f5auNMQ37WFqMIsyQd7hqhCWVPQehpfYuBOrtZX14mBO9aOzM1giFKh3iqP6wPcudasrsaRhTv0jfLUEFbIYVSwvajRQofOg4YVPemdQXEMVohwzzwQ0Vw5A3fUmxwKYYXUB+ydFV3pnd69pL4CR3ceaq+qS5ClO6cMWqf/PLon64Q5QYWwN+ieOajv7+fssHlhFjGkwuGdm20ebPjAQnEtrLzF4LBOGOaQjuTDVsQl0s+st8Y7T1phYs5gvNO51UYAPIz0NTR21JvusDlrDXxnih0w92S7o62Z2OeTd8T+BxXV2WJs9yIRwdF0L3JSMaLpzvPyWXN0pVcergNxIYFH9IWqRWe6PCcfD7M13nmUPk9GjbXeoY+lk85Z74We1lGjM/3sdWK0uDTN6vMsB8YFI1ErqUcXDHoDE+HJepe0tnwCHK0LBtMes2yFkHHqPNIeIZ5FnWjrRJ5bFwxgiRHRg0td/GUJI3YGwNNhtmpWj5N1wXD4Tg3KcC4Y9sQ8ZxbGBYNRgZbVM2HrguncfGjSgXPBMK57nNFgXTCDmRY7uNZlgzdaUpeH98CgYfNT2Dg9LaOo+eI8MIu74ROBYj0wW+KwNIjOemCYHQF4R8AnRr2MFjJ85DTwddlygM8AbdbNigM8kFfGAbkFvHhZNXjLAZ7WYTouFQN4mjblBD8avtPESLXHKBfxqazeNYTGemA4XFafwQODH8CM29EDw8UVM3fswHculpjaGsRuneySFlVy5LtgNesR4Cf5rmz/7Q0dj1/c2+13t3T7/YtN0Lx9W4KmSwk2sOdBDEBegsedxzZYd5cPUaQ8Xkl9wZ5HUyBkCx73wbDPrKaWgT1k7knEMjO0Z/YKANA87QmCdLZxhvaQMXDHCgeng5u1MqM1j5IylNaHk5/zmi1HUIb2g0MylxyM+cGYBLROCFbM2gon/tzQnodT2G7rxDe0p45fp+hwlyydkWWQWdpTTz2HmEV5bT/ne4b1LOXK6clTw0rB3Cshrpw65sHKPjdImgaDJcTDUJYA1ehtZ/uW15XQoJ7dwXjs6KjhKSkGg9raBvWDc2usJ1t+8OgmndhwGxTDvXw5YfEG9dAfWSQW9UMibE/XXqjnyEPfzGjLQ5+o+I6eGujY0Z70I4N6Duyaju1vWM/jS1RW/c8+RSgzRyGGxfAAGH39ZMxzmjEfdATWD0kPGTWynpMVEyOyXqY25nA05nmAi3muR7KG9YMZSCdO0h6n0o0w9wz5nNQTfr+fHDWbx15bn2Nt+c0UkKFnCNaWB3Xwd9sXjlMxSHRiOVueKY1rnYPThy2/eYat6cDOlt8wrqoe7DhbnvU+kdXOlicgVq8t2vJM7YScfynW3/KLc75/OM8Trc8JQJ+xAzwKiyGWd4iMwbQASaWehu9QFyO1l+P7XHSU6gQxfJ/M+ceeM5yoTubWN3UoGr5P9f1XH4dOGcNpdR+ILjI4Gbw1UzLESzxRpdxgdARrfkrOdw9hjiKvV6/4BXjKJ6rK8l3krfdVWL6LjjaJzhp56VbTxfKdutifge9TstuPV9kAfsoZ10laNISX9mqpRGcNdeZKRmeNNDvaJnvCs5NovATCU8Z4j7Y8e7oJPTzgmS2KJlZCGsBzwODpKeQPiQ6Olxj1SEZXRrcEwHM4YvikEJ8uv0f5UnTWsDxNcOMAj+FWeOAeE0DXYujUjAmg2EQRAur6NoAvneMsJoCSNZgE6tO46C4IwhKtnLroTn33dYIbTQIoqd+P39vkf1JG16yY/zkkLadrlv2V/0kZa1ANhrzqOceYR33M2D3EPIrOgRATQKknWIvBU6PF7PME0VwJoCcN9ikBVIKcwJ/q6X5a4Ti4TAIo9Vo18uJKAJWVElvGECwjbbzXiH540RezVjzdRR9ZrV2bAEq91fwtiftvW/JXzue3ZYNKHGkSCjnQQ8aUX8I1A/qFaTMBvBA7A7kwS9OTfvFMEyzYnvSL0ZWYaMExD3lgaoSsUMqzZI32NqSnPJbGtRnSL+7pksakGdRTxgIt5y4G9ZTrqCugnvJr4pJBvcj5pK1fqJeH9JOCaFgv72zl6dRVioiOGoH1rBFWrhZteeodK1L03EjDJA27sKxnO85aFXOG9dTx3BrDINkdsM9nDINk76Gc0THPrsacX9GcXxKIpLUysGfftMNiw3rJh8IYCAHuZO7i6Wdk/WDsCSZVuIhFj81nj7Z8ojHRZ422fEpM+zrr3YV68Bez89xZ8yA9g/AL9+oe9BgowF2P0ZCbMbm5Fk95JgOg2GoeXJjHjysjikIwJOTW6gjuGqq8XiN4a/gMXjYYvPGUz+UvhvBU0X96X8oD8FunrZLt4jtlTPESTllFnsdfctFdf32iAy+4y7P7jvHsUpCegg+eKqZkC8HsUhcsTNtznXLjcuixLg0yTpzQK9Up8jLAkC0qbbqezlYp16LXERmkbxllZYc8JMo8Lf6Gg9U3gX6ljt6+LXUUcwTPZF5KALoY9DUHy53ZDq1qtoMFeiUw9SIdA3SomKbhAi2R61zqcTdA5zNACn+Flsgg8RPQMWe2Fs/ynGfatFwDz7k6Db3xwfKc846mX+A5vSRpruCJp8ykOJ/sLzI426Nvhi21h0bkOp5Dn6/+DctzOlVmzzGsnXptGsjleM68k3l8P5bn9BX1ucK1LZpI8gVPvOTH1LJjMKTk05Sszzc8L+zTWkJcu5r64FlI/Nc0nnEeb4lemcpaSiQ6dWwj4+0ty3p5DNBRHIyC17iYi+gsJuNHYiwks41WSfF2LeodgDlWuk1awva6xPxSWY8ws6dnOv/imq9B9TZnKfO0Q6nicpYwFlt6vr+FKVTcl3mwU8bU06Xa+mYwLSbDZatHO/tk4385+ma4E8JfKtnRnZv1jpWzh+tbuLffTAnYnu8SYshbmMJB6+DVbVPjPg3g6S9jeFsw3umfwgKrW4+L8HQ3pZxycMyIV6mpo9IQnn5kcD8HwmtqQdP9wivhJbEgpxIOWOmAWufCCEN4ypieI9zaQhnTZ4cQd8pofCHHL0X4K+309m1ppxxBDLDI3vlemBUE8ob7EQvTXMZxPF6ERx+BYLBbnfOdYxbzp4Vwd5iZdKdrGKIhPGDKW4PkdMwQfvPGu6wuREN4zJMxmqYdGMTDzGKH1OB9B79Z1B3iJJkwWWOUJK2srUGthu/Mo81Zb1szfN+8/Qc/n4HveMnkpb078H0zFjQvjao2fN/oSh5OjcD3zZC2cfYDlu9yS1puMZBmMXAYPRRPWhfRvYfy2vKdbol1Qvgs36cEiyR/rQtlYH+NGOzOewNq04ng+I6SJWzOYrA7dLzzpNdavg/u4Ve4OUDYWtu5yMzyvdMVt0q8nWsx5J8hM5HvEgDb4zkrN1GA8fQXvMgeZ4z0RPcmp1yaX2Xxzpzcc4+BpbskT6QZ7g+QLU6p2pQhIxWQzTsGScqKnjVEKmSk8lDoCe+yFZ8nd99npGLwp/WckZrAshmDJLmE8qaWiHfGcnM1Kx7vsoQysTX63rmEolo5XPhCvWGKxtsUZW3trZybYExGatH7b4MVf/Q9ou8dHVTWuX7x4XqXOxOGek6t6136/jXa3wa6856ec0mMjXPnKQzMsuh5h06b5pd0zVzJqbdvSk51Nw4bzOfBqy/L9gE1lCe9jwHzkDsm83KeGaq8HS4HzGeGrS3NJrgwTzmlrI6uC/OFWVeMnfe3DFBGH+gtNRfmKVdsZLK/aIAyJuNSX8sD87yQjolHIT1V7qmDxdG9E54yqn7yiB6kp1wqRk0gPaccdiU1emaoY0U8wSQX6WWKYlmNaU1ysXM+Dp4L9JDofz5+gwv01Ec/wcsG9NSxi1czzYCeeuonzdeAHnOKj0zbX45LudG/sz3oqaexTiTlBXred4dZrzHWBvTUK4yeHrKaCDAwIvdgyQvY6AwOl8RQ59nxCCE1JXGSnkPlC/SUYcfmePtAQb+h9C3cE0M5swHDMSvKJBM5BNRQTgPLe7gql7TGWqqx1NYNz3xsmjYhooZ6aXMFDw1tLFi89dyta+5hLLwPrYxwyEq9vCamWT883sS7z2rwwye57q6mkNJEXXIORiA974ngChBJv2nJdzXjLOlZcab5lED6zevb6U0KpN+StLrOJTMX6TFTMUJ0S2pJr0m0JV7vxY7gVXnTn7JuhsOLSexRvxnytk58u0H9JpXXSYk1rN+M0zwH5Qb1tM9S10TfX8xr84iefAqo+VgKaylLrrrJ3g1PGeMhxNOIiumXA+sLp1AKFzOKCqOx+PNWyglFbP4+GZHzSf4zrC+MGMt6LmVYTznxjNezvjDD77i+DeuZx7X3DFd/idzayXG6WF848bcuaIb1lJnp4eMnKXd+iKYE1lPHYptCRA31xoDFYNVTr1wCg1UvetV8dMt6yLyrUO/8NaynzpKFkBrqbMUenPDUeU1ZCclPpXAGjhTuahQZnC0hfJL6ghV9Ymcu1rNkjFAIyU/UB5bJHIx66lioaglueOqSdhZiaqiTJTlkP5FVPDLpIfuJemMaUgmwL3ogk318PC9yXf31unoD+yIBgurVtrRnv2Ehb+GqGVlRW4mw578u7JCfElgnX9vDxeiS4tRaO1HzJnhSjsbOjenGqmebcY3JwapH4fvGi0JATZG7PMeKATVFjEhdli3roWMDXltMfZInbs0ztKzXy3LbjPfMFHF9j2jUc4TAOKnhjnTmU6I4q8XgScYHg+47XjPDMHD6JrJDPQNIc+ol3k/ASM75WBjcDQVopjJi+mrnQXLXDG0XUMMZtqqMy1+K9TaX9R2gV949DrPZu+FFRv9EJ00VI2p074YX+dzxcwGdmcM8WQzGO+WTAWN4ThVdUILtTjnN2ryLpmguZwpXEhT9AlUPAZIiwzjY/r5dkVPX1d7wvPI7Riejx/Cccqdx4HlOmf7hEEEjesonatDwnGnWs2kugOV5lVCrkz1qeE6dfRiuJKCOrbCGBligU4cBGT9rIXqjhyIAXdLBXxcGA/TKrwgMvYHbAJ1yXyft0wCdes16QmeBTh3tGq9Qp77mOQOzQKcOSyveMCM62nMFN3yRXF+sDOELF6LzwDxcwys6rzAIVxKIjk1s8zE0lAfvLgzprKK33VaMoaEOBus5rwG65O5nNYCt+S4NXs5NJsZ8p15f7+w29jv1VJ+CJKm3SYd+sN+pt5MVa8x3ymiz5q97pIyuOjmul/VOuZ/zWxNEQxlmWEx3ol4Y8xhyWUXvu7fgphG9nJgYE0PD2xT4KYHmgS46Fq8VL1GvjLKa5zZgE0PDW1RmPVuhRwwN73wgHLYHepG7Hfpa4T4C0deJSLYxNJKh8Hq4bWJoqKPR6i9pvE/jm/mWxFYA6K57SEd6qrmsSPou15yvkPgkMj860h3quwTjzmi7Uy7Hz2lY3yV+fzaf+UR5rJ4l1MywnvJoYwbbnTJGRwt+GsqFX+PwrO8Swq93BhrWdx5ztNX97eoij67ZVob1lJkoGHKfRC/nqgTLeup0UYd7C4ok7K1zXa1hPXWQQa1Zw3rqPdf4BSPRYUukyPouV1j1Gk5ci3B5lOiRF32qk82gnupIx49iUE8dmDpXQV6ol1z41uNV66KLfzCgnumM66SNW9RTZwR6tN2p95OTblFPHZMpXigmejlhbxb11FM910deqIec+f2EkPok+tgxHF5kpryG1FbR6ewOqU+i87KKSHrqzFWNnhpeOMA872C9i44RtaOnhjrWnRQjJqnXfS5MuFBPufBMN7CeeqI3IsCeOak7nw94GNhTn7S7A+ypYzuY4yczNLc1txAPL3rZzxGT1NE9OZy7UqcVMsPBq+gAW3a5T6KOFVOfRO713NJgYE+dlyJE2FMvLa2Q+iQ6Orl/Q2rrm7BP1lPzLXmuhXkfAHUJbvkhN4qr1WtwTxnVH8FVM3Qn6uNrRC1dsWlwT5mXkAXc8x5AdFX2AZMiTzSvz30SGbuz6QNsRMa8Wj42XmRGVfjcJ5HTIaPBPXNeVk/hnvZyUmFGMO0pw2RTo83gnnpNMfdJ5McxrqE9E2o2lofoladOWzTE14jeW44XAYuO98/oqpGEnaanepb20MeuLX4YSfSZ9Xt/BveU2349lr1wT7084mUu3DMfiAGQ4UZg0ZlSEF011JnEGA5giyTGFU2Ms7iXvKKH6+XCPe+PhI3aIu4Zow0T+An3kgFYz92fBvdEaj+t43FPv9cIrho6dfM+aUj2DHbSM9DTjBE2c+oNLeFDG9QxlvSsx3prgBPU4yRz2hAb9ETaen7gQmwGh2/u8aIaTDDeRxHukqTMuPMWIiipt5JPcIzx1kzeePd69Gu8NZPnXeucSxpvDfTBG2Cjt2YyUjSdTyoZdw10cYXHuwymHF2/cv1y10y5/ehE1Bh3zWSwWKvbH8JOxnKN83ER466ZzJSsGuVq3TX8IC/7PJ7CTt50X/Q7FNZdA7jtdq5PPMDnHQXp8aFwT+Uvfqj8y98ex9O+9AXzv7z1BXP8/hOfQXe/vh7z3tO/+/uq30H/NT/djn/kJszzIfeXsuWTx/INhsGrjfRL6JkDbBrtJ2hDP/D+0BpG7RCF9+95hQ4bKj++PDT2Y12vD3uI/Xwf/vF87DdO+R/FuKQfbYkvGb27+WWxIBe5Iv2Iry9y4muR+NyH/Cj9T059raV51aMtvtSOP/JD9L96/RC93IGB8fX87dqPpka/YCQWhjydwvMmVuZmcSzyjDXKP71UrkrPMkZFktF+hgMGA/7dfeOetzW8/sXCWxSn+yp4CR/btp/yfevbj1/JCbQ1E19ifarZQ/Y1M/LnakZDTnvqzZrlL9fsW/NhbCX58bn03H0P2VfSyJ+rZFbz7r1Klrcq+W0x4raSicvEc08+ZF9JI3+qkjyH1sn9diXrW5X8tjBJU0lGpvXniXjJrpJW/lwlJezp/Z5sb1byW2KEbB3HeqwOto4P2dfRyJ+rI9f6r0zJ/uZo/aazcVtJXqH0PCUv2VfSyJ+rJG/P+sqUHG9V8p1DIVuTykvcnrvrIfuaGPlzNSk8IXq/u+ZbNfkmb6itI6PFvjDtHrKvo5E/V8ecv7r+rbfq+G1OAFNJfqr2edo9VFdFo36qhkyg+sqc27aCb9wA8GGT5iUa141xzNGimeNJRW3X+oIqtdWa/va1MF+oMcPC80dNmpw+UOVPGTimnta+ueppzZurntG6+Uw9P2Lg5I917cfMnbe71lo7V5WtsXNVOdo6n6nyB8ydXD5U5Q8aP29X2do+V5Wt6XNVOVo+n6jyR4yfXD9U5Q+aQm9W2VlCjyo7Q+hR5Sc76DNV/oAplNvHqvxZw+jF19jaRVeNrVl01ThaRZ+p8QcMo9w/Nq4/Zia93cnWSrqqbI2kq8rRRvpMlT9gJuXxoSq/YzSZelmb6aqXNZmuekWL6TP1+oDRlOeH6vVpEyoMXmtBXTW2BtRV42g/fabGHzCh8vpQjT9oUL05eK099aixNaceFY7W1Cfq+wGDKjuL6jcv/x+HgBfVCmVuZHN0cmVhbQplbmRvYmoKMTIgMCBvYmoKMTAzNTQKZW5kb2JqCjEwIDAgb2JqClsgXQplbmRvYmoKMTcgMCBvYmoKPDwgL0xlbmd0aCAyNDcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicTVFJbsQwDLv7FfzAAJasxXlPikEP7f+vJR0U7cEQI0tc4u7ERBZetlDXQofjw0ZeCZuB74PWnPgaseI/2kaklT9UWyATMVEkdFE3GvdIN7wK0X6kgleq91jzEXcrzVs6drG/98G05pEqq0I85Ngc2Uha10TR8T203nNDdMoggT43IQdEaY5ehaS/9sN1bTS7tTazJ6qDR6aE8kmzGprTKWbIbKjHbSpWMgo3qoyK+1RGWg/yNs4ygJPjhDJaT3asJqL81CeXkBcTccIuOzsWYhMLG4e0H5U+sfx86834m2mtpZBxQSI0xaXfZ7zH53j/AJVPXCYKZW5kc3RyZWFtCmVuZG9iagoxOCAwIG9iago8PCAvTGVuZ3RoIDc5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM3NVIwULC0ABJmpiYK5kaWCimGXEA+iJXLZWhpDmblgFkmxgZAlqmpKRILIgvTC2HB5GC0sYk51AQECyQHtjYHZlsOVwZXGgDWlBwMCmVuZHN0cmVhbQplbmRvYmoKMTkgMCBvYmoKPDwgL0xlbmd0aCAzMDcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPZJLbgMxDEP3PoUuEMD62Z7zpCi6mN5/2ycl6Yoc2RZFapa6TFlTHpA0k4R/6fBwsZ3yO2zPZmbgWqKXieWU59AVYu6ifNnMRl1ZJ8XqhGY6t+hRORcHNk2qn6sspd0ueA7XJp5b9hE/vNCgHtQ1Lgk3dFejZSk0Y6r7f9J7/Iwy4GpMXWxSq3sfPF5EVejoB0eJImOXF+fjQQnpSsJoWoiVd0UDQe7ytMp7Ce7b3mrIsgepmM47KWaw63RSLm4XhyEeyPKo8OWj2GtCz/iwKyX0SNiGM3In7mjG5tTI4pD+3o0ES4+uaCHz4K9u1i5gvFM6RWJkTnKsaYtVTvdQFNO5w70MEPVsRUMpc5HV6l/DzgtrlmwWeEr6BR6j3SZLDlbZ26hO76082dD3H1rXdB8KZW5kc3RyZWFtCmVuZG9iagoyMCAwIG9iago8PCAvTGVuZ3RoIDczIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDO2NFAwULAwU9A1NDZUMLI0VjA3M1BIMeQCCoFYuVwwsRwwy8wSxDI0N0Ni6ZoZQmWRWCDjcrhgBufAzMvhyuBKAwAeiRaVCmVuZHN0cmVhbQplbmRvYmoKMjEgMCBvYmoKPDwgL0xlbmd0aCA2OSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwztjRQMFCwNFfQNTQ2VDA2MFEwNzNQSDHkgjFzwSywbA4XTB2EZQZiGBmaILHMgMaBJeEMkBk5cNNyuDK40gD6qRZFCmVuZHN0cmVhbQplbmRvYmoKMjIgMCBvYmoKPDwgL0xlbmd0aCAyMzIgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVFJbsQwDLv7FfzAANbuvCfFoIf2/9dSyhQIQCW2uCViYyMCLzH4OYjc+JI1oyZ+Z3JX/CxPhUfCreBJFIGX4V52gssbxmU/DjMfvJdWzqTGkwzIRTY9PBEy2CUQOjC7BnXYZtqJviHhsyNSzUaW09cS9NIqBMpTtt/pghJtq/pz+6wLbfvaE052e+pJ5ROI55aswGXjFZPFWAY9UblLMX2Q6myhJ6G8KJ+DbD5qiESXKGfgicHBKNAO7LntZ+JVIWhd3adtY6hGSsfTvw1NTZII+UQJZ7Y07hb+f8+9vtf7D04hVBEKZW5kc3RyZWFtCmVuZG9iagoyMyAwIG9iago8PCAvTGVuZ3RoIDIzMSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1TzmSBCEMy3mFPjBVGNtAv6entjbY+X+6kplOkPAhydMTHZl4mSMjsGbH21pkIGbgU0zFv/a0DxOq9+AeIpSLC2GGkXDWrONuno4X/3aVz1gH7zb4illeENjCTNZXFmcu2wVjaZzEOclujF0TsY11radTWEcwoQyEdLbDlCBzVKT0yY4y5ug4kSeei+/22yx2OX4O6ws2jSEV5/gqeoI2g6Lsee8CGnJB/13d+B5Fu+glIBsJFtZRYu6c5YRfvXZ0HrUoEnNCmkEuEyHN6SqmEJpQrLOjoFJRcKk+p+isn3/lX1wtCmVuZHN0cmVhbQplbmRvYmoKMjQgMCBvYmoKPDwgL0xlbmd0aCAyNDkgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPVA7jkQhDOs5hS/wJPIjcB5Gqy1m79+uA5opUEx+tjMk0BGBRwwxlK/jJa2groG/i0LxbuLrg8Igq0NSIM56D4h07KY2kRM6HZwzP2E3Y47ARTEGnOl0pj0HJjn7wgqEcxtl7FZIJ4mqIo7qM44pnip7n3gWLO3INlsnkj3kIOFSUonJpZ+Uyj9typQKOmbRBCwSueBkE004y7tJUowZlDLqHqZ2In2sPMijOuhkTc6sI5nZ00/bmfgccLdf2mROlcd0Hsz4nLTOgzkVuvfjiTYHTY3a6Oz3E2kqL1K7HVqdfnUSld0Y5xgSl2d/Gd9k//kH/odaIgplbmRzdHJlYW0KZW5kb2JqCjI1IDAgb2JqCjw8IC9MZW5ndGggMzk1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1SS27FQAjb5xRcoNLwm895UlXdvPtva0NSqSq8iTHGMH3KkLnlS10ScYXJt16uWzymfC5bWpl5iLuLjSU+ttyX7iG2XXQusTgdR/ILMp0qRKjNqtGh+EKWhQeQTvChC8J9Of7jL4DB17ANuOE9MkGwJOYpQsZuURmaEkERYeeRFaikUJ9Zwt9R7uv3MgVqb4ylC2Mc9Am0BUJtSMQC6kAAROyUVK2QjmckE78V3WdiHGDn0bIBrhlURJZ77MeIqc6ojLxExD5PTfoolkwtVsZuUxlf/JSM1Hx0BSqpNPKU8tBVs9ALWIl5EvY5/Ej459ZsIYY6btbyieUfM8UyEs5gSzlgoZfjR+DbWXURrh25uM50gR+V1nBMtOt+yPVP/nTbWs11vHIIokDlTUHwuw6uRrHExDI+nY0peqIssBqavEYzwWEQEdb3w8gDGv1yvBA0p2sitFgim7ViRI2KbHM9vQTWTO/FOdbDE8Js753WobIzMyohgtq6hmrrQHazvvNwtp8/M+iibQplbmRzdHJlYW0KZW5kb2JqCjI2IDAgb2JqCjw8IC9MZW5ndGggMTM2IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE2PQQ4DMQgD73mFn0AgQHjPVlUP2/9fS9h20wseyYBsUQaBJYd4hxvh0dsP30U2FWfjnF9SKWIhmE9wnzBTHI0pd/Jjj4BxlGosp2h4XkvOTcMXLXcTLaWtl5MZb7jul/dHlW2RDUXPLQtC12yS+TKBB3wYmEd142mlx932bK/2/ADObDRJCmVuZHN0cmVhbQplbmRvYmoKMjcgMCBvYmoKPDwgL0xlbmd0aCAyNDkgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicTVFJigMwDLvnFfpAIV6TvKdDmUPn/9fKDoU5BAmvkpOWmFgLDzGEHyw9+JEhczf9G36i2btZepLJ2f+Y5yJTUfhSqC5iQl2IG8+hEfA9oWsSWbG98Tkso5lzvgcfhbgEM6EBY31JMrmo5pUhE04MdRwOWqTCuGtiw+Ja0TyN3G77RmZlJoQNj2RC3BiAiCDrArIYLJQ2NhMyWc4D7Q3JDVpg16kbUYuCK5TWCXSiVsSqzOCz5tZ2N0Mt8uCoffH6aFaXYIXRS/VYeF+FPpipmXbukkJ64U07IsweCqQyOy0rtXvE6m6B+j/LUvD9yff4Ha8PzfxcnAplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9MZW5ndGggOTQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRY3BEcAgCAT/VEEJCgraTyaTh/b/jRAyfGDnDu6EBQu2eUYfBZUmXhVYB0pj3FCPQL3hci3J3AUPcCd/2tBUnJbTd2mRSVUp3KQSef8OZyaQqHnRY533C2P7IzwKZW5kc3RyZWFtCmVuZG9iagoyOSAwIG9iago8PCAvTGVuZ3RoIDE2NCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxFkMdxBTEMQ++qAiUwgAr1rMfzD+v+r4b000F6GEIMYk/CsFxXcWF0w4+3LTMNf0cZ7sb6MmO81VggJ+gDDJGJq9Gk+nbFGar05NVirqOiXC86IhLMkuOrQCN8OrLHk7a2M/10Xh/sIe8T/yoq525hAS6q7kD5Uh/x1I/ZUeqaoY8qK2seatpXhF0RSts+LqcyTt29A1rhvZWrPdrvPx52OvIKZW5kc3RyZWFtCmVuZG9iagozMCAwIG9iago8PCAvTGVuZ3RoIDcyIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiGeAmCBtEMUgFkSxmYkZRB2cAZHL4EoDACXbFskKZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvTGVuZ3RoIDgzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD3MORKAMAgF0J5T/COEyCL3cRyLeP9WMNEGHqt6oCE4g7rBreFgyrp0E+9T49XGnBIJqHhKTZa6C3rUtL7Uvmjgu+vmS9WJP83PF50Pux0Z3QplbmRzdHJlYW0KZW5kb2JqCjMyIDAgb2JqCjw8IC9MZW5ndGggMTYzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWQOxIDIQxDe06hI/gjAz7PZjIpNvdvY9hsUsDTWCCDuxOC1NqCieiCh7Yl3QXvrQRnY/zpNm41EuQEdYBWpONolFJ9ucVplXTxaDZzKwutEx1mDnqUoxmgEDoV3u2i5HKm7s75Q3D1X/W/Yt05m4mBycodCM3qU9z5NjuiurrJ/qTH3KzXfivsVWFpWUvLCbedu2ZACdxTOdqrPT8fCjr2CmVuZHN0cmVhbQplbmRvYmoKMzMgMCBvYmoKPDwgL0xlbmd0aCAzMjIgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVG7bcUwDOw1BRcwIH4lzeMgSJG3f5s72qlI07wfVV4ypVwudckqWWHypUN1iqZ8nmam/A71kOOYHtkhulPWlnsYFpaJeUodsZos93ALNr4AmhJzC/H3CPArgFHARKBu8fcPulkSQBoU/BTomquWWGICDYuFrdkV4lbdKVi4q/h2JLkHCXIxWehTDkWKKbfAfBks2ZFanOtyWQr/bn0CGmGFOOyzi0TgecADTCT+ZIBszz5b7OrqRTZ2hjjp0ICLgJvNJAFBUzirPrhh+2q75ueZKCc4OdavojG+DU7mS1LeV7nHz6BB3vgzPGd3jlAOmlAI9N0CIIfdwEaEPrXPwC4Dtkm7d2NK+ZxkKb4ENgr2qFMdyvBi7MxWb9j8x+jKZlFskJX10ekOytygE2Ieb2ShW7K2+zcPs33/AV8Ze2QKZW5kc3RyZWFtCmVuZG9iagozNCAwIG9iago8PCAvTGVuZ3RoIDgzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWMuw3AMAhEe6ZgBH4m9j5RlMLevw0QJW64J909XB0JmSluM8NDBp4MLIZdcYH0ljALXEdQjp3so2HVvuoEjfWmUvPvD5Se7KzihusBAkIaZgplbmRzdHJlYW0KZW5kb2JqCjM1IDAgb2JqCjw8IC9MZW5ndGggNTEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicMza0UDBQMDQwB5JGhkCWkYlCiiEXSADEzOWCCeaAWQZAGqI4B64mhyuDKw0A4bQNmAplbmRzdHJlYW0KZW5kb2JqCjM2IDAgb2JqCjw8IC9MZW5ndGggMTYwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWQORIDMQgEc72CJ0hcgvesy7XB+v+pB9ZHoukCNBy6Fk3KehRoPumxRqG60GvoLEqSRMEWkh1Qp2OIOyhITEhjkki2HoMjmlizXZiZVCqzUuG0acXCv9la1chEjXCN/InpBlT8T+pclPBNg6+SMfoYVLw7g4xJ+F5F3Fox7f5EMLEZ9glvRSYFhImxqdm+z2CGzPcK1zjH8w1MgjfrCmVuZHN0cmVhbQplbmRvYmoKMzcgMCBvYmoKPDwgL0xlbmd0aCAzMzQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicLVJLcsUgDNtzCl2gM/gH5DzpdLp4vf+2kpNFRg5g9DHlholKfFkgt6PWxLeNzECF4a+rzIXPSNvIOojLkIu4ki2Fe0Qs5DHEPMSC76vxHh75rMzJswfGL9l3Dyv21IRlIePFGdphFcdhFeRYsHUhqnt4U6TDqSTY44v/PsVzLQQtfEbQgF/kn6+O4PmSFmn3mG3TrnqwTDuqpLAcbE9zXiZfWme5Oh7PB8n2rtgRUrsCFIW5M85z4SjTVka0FnY2SGpcbG+O/VhK0IVuXEaKI5CfqSI8oKTJzCYK4o+cHnIqA2Hqmq50chtVcaeezDWbi7czSWbrvkixmcJ5XTiz/gxTZrV5J89yotSpCO+xZ0vQ0Dmunr2WWWh0mxO8pITPxk5PTr5XM+shORUJqWJaV8FpFJliCdsSX1NRU5p6Gf778u7xO37+ASxzfHMKZW5kc3RyZWFtCmVuZG9iagozOCAwIG9iago8PCAvTGVuZ3RoIDcwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMzNlMwULAwAhKmpoYK5kaWCimGXEA+iJXLBRPLAbPMLMyBLCMLkJYcLkMLYzBtYmykYGZiBmRZIDEgujK40gCYmhMDCmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0xlbmd0aCAzMjAgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVJLbgUxCNvPKbhApfBPzvOqqou++29rE70VTDBg4ykvWdJLvtQl26XD5Fsf9yWxQt6P7ZrMUsX3FrMUzy2vR88Rty0KBFETPViZLxUi1M/06DqocEqfgVcItxQbvINJAINq+AcepTMgUOdAxrtiMlIDgiTYc2lxCIlyJol/pLye3yetpKH0PVmZy9+TS6XQHU1O6AHFysVJoF1J+aCZmEpEkpfrfbFC9IbAkjw+RzHJgOw2iW2iBSbnHqUlzMQUOrDHArxmmtVV6GDCHocpjFcLs6gebPJbE5WkHa3jGdkw3sswU2Kh4bAF1OZiZYLu5eM1r8KI7VGTXcNw7pbNdwjRaP4bFsrgYxWSgEensRINaTjAiMCeXjjFXvMTOQ7AiGOdmiwMY2gmp3qOicDQnrOlYcbHHlr18w9U6XyHCmVuZHN0cmVhbQplbmRvYmoKNDAgMCBvYmoKPDwgL0xlbmd0aCAxOCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNrRQMIDDFEOuNAAd5gNSCmVuZHN0cmVhbQplbmRvYmoKNDEgMCBvYmoKPDwgL0xlbmd0aCAxMzMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRY9LDgQhCET3nKKOwMcf53Ey6YVz/+2AnW4TYz2FVIG5gqE9LmsDnRUfIRm28beplo5FWT5UelJWD8ngh6zGyyHcoCzwgkkqhiFQi5gakS1lbreA2zYNsrKVU6WOsIujMI/2tGwVHl+iWyJ1kj+DxCov3OO6Hcil1rveoou+f6QBMQkKZW5kc3RyZWFtCmVuZG9iago0MiAwIG9iago8PCAvTGVuZ3RoIDM0MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UjluBDEM6/0KfSCAbtvv2SBIkfy/DanZFANxdFKUO1pUdsuHhVS17HT5tJXaEjfkd2WFxAnJqxLtUoZIqLxWIdXvmTKvtzVnBMhSpcLkpORxyYI/w6WnC8f5trGv5cgdjx5YFSOhRMAyxcToGpbO7rBmW36WacCPeIScK9Ytx1gFUhvdOO2K96F5LbIGiL2ZlooKHVaJFn5B8aBHjX32GFRYINHtHElwjIlQkYB2gdpIDDl7LHZRH/QzKDET6NobRdxBgSWSmDnFunT03/jQsaD+2Iw3vzoq6VtaWWPSPhvtlMYsMul6WPR089bHgws076L859UMEjRljZLGB63aOYaimVFWeLdDkw3NMcch8w6ewxkJSvo8FL+PJRMdlMjfDg2hf18eo4ycNt4C5qI/bRUHDuKzw165gRVKF2uS9wGpTOiB6f+v8bW+19cfHe2AxgplbmRzdHJlYW0KZW5kb2JqCjQzIDAgb2JqCjw8IC9MZW5ndGggMjUxIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nC1RSXIDQQi7zyv0hGan32OXK4fk/9cIygcGDYtAdFrioIyfICxXvOWRq2jD3zMxgt8Fh34r121Y5EBUIEljUDWhdvF69B7YcZgJzJPWsAxmrA/8jCnc6MXhMRlnt9dl1BDsXa89mUHJrFzEJRMXTNVhI2cOP5kyLrRzPTcg50ZYl2GQblYaMxKONIVIIYWqm6TOBEESjK5GjTZyFPulL490hlWNqDHscy1tX89NOGvQ7Fis8uSUHl1xLicXL6wc9PU2AxdRaazyQEjA/W4P9XOyk994S+fOFtPje83J8sJUYMWb125ANtXi37yI4/uMr+fn+fwDX2BbiAplbmRzdHJlYW0KZW5kb2JqCjQ0IDAgb2JqCjw8IC9MZW5ndGggMTc0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE2QSQ5DIQxD95zCF6iEM8DnPL+qumjvv61DB3WB/OQgcDw80HEkLnRk6IyOK5sc48CzIGPi0Tj/ybg+xDFB3aItWJd2x9nMEnPCMjECtkbJ2TyiwA/HXAgSZJcfvsAgIl2P+VbzWZP0z7c73Y+6tGZfPaLAiewIxbABV4D9useBS8L5XtPklyolYxOH8oHqIlI2O6EQtVTscqqKs92bK3AV9PzRQ+7tBbUjPN8KZW5kc3RyZWFtCmVuZG9iago0NSAwIG9iago8PCAvTGVuZ3RoIDc2IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2MOw6AMAxD95zCR2h+JAdCiIHef6UptIv99CTbxdFgWpECt8DJ5D6p03LPJDt8EJsh5FcbWrWuytKaDIuajL8N391N1wumOBfACmVuZHN0cmVhbQplbmRvYmoKNDYgMCBvYmoKPDwgL0xlbmd0aCAyMTUgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVE5DgMhDOz3Ff5AJIwveE+iKM3+v82M0VYewVyGtJQhmfJSk6gh5VM+epkunLrc18xqNOeWtC1zgLi2vC+tksCJZoiDwWmYuAGaPAFD19GoUUMXHtDUpVMosNwEPoq3bg/dY7WBl7Yh54kgYigZLEHNqUUTFm3PJ6Q1v16LG96X7d3IU6XGlhiBBgFWOBzX6NfwlT1PJtF0FTLUqzXLGAkTRSI8+Y6m1RPrWjTSMhLUxhGsagO8O/0wTgAAE3HLAmSfSpSz5MRvsfSzBlf6/gGfR1SWCmVuZHN0cmVhbQplbmRvYmoKMTUgMCBvYmoKPDwgL1R5cGUgL0ZvbnQgL0Jhc2VGb250IC9CTVFRRFYrRGVqYVZ1U2FucyAvRmlyc3RDaGFyIDAgL0xhc3RDaGFyIDI1NQovRm9udERlc2NyaXB0b3IgMTQgMCBSIC9TdWJ0eXBlIC9UeXBlMyAvTmFtZSAvQk1RUURWK0RlamFWdVNhbnMKL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXQovQ2hhclByb2NzIDE2IDAgUgovRW5jb2RpbmcgPDwgL1R5cGUgL0VuY29kaW5nCi9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0NiAvcGVyaW9kIDQ4IC96ZXJvIC9vbmUgL3R3byAvdGhyZWUgL2ZvdXIgL2ZpdmUgL3NpeCAvc2V2ZW4KL2VpZ2h0IC9uaW5lIDcxIC9HIC9IIDkxIC9icmFja2V0bGVmdCA5MyAvYnJhY2tldHJpZ2h0IDk3IC9hIDk5IC9jIC9kIC9lIC9mCjEwNCAvaCAvaSAxMDcgL2sgMTEwIC9uIDExNCAvciAvcyAvdCAvdSAxMjIgL3ogXQo+PgovV2lkdGhzIDEzIDAgUiA+PgplbmRvYmoKMTQgMCBvYmoKPDwgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9Gb250TmFtZSAvQk1RUURWK0RlamFWdVNhbnMgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM0MiA+PgplbmRvYmoKMTMgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTYgMCBvYmoKPDwgL0cgMTcgMCBSIC9IIDE4IDAgUiAvYSAxOSAwIFIgL2JyYWNrZXRsZWZ0IDIwIDAgUiAvYnJhY2tldHJpZ2h0IDIxIDAgUgovYyAyMiAwIFIgL2QgMjMgMCBSIC9lIDI0IDAgUiAvZWlnaHQgMjUgMCBSIC9mIDI2IDAgUiAvZml2ZSAyNyAwIFIKL2ZvdXIgMjggMCBSIC9oIDI5IDAgUiAvaSAzMCAwIFIgL2sgMzEgMCBSIC9uIDMyIDAgUiAvbmluZSAzMyAwIFIKL29uZSAzNCAwIFIgL3BlcmlvZCAzNSAwIFIgL3IgMzYgMCBSIC9zIDM3IDAgUiAvc2V2ZW4gMzggMCBSIC9zaXggMzkgMCBSCi9zcGFjZSA0MCAwIFIgL3QgNDEgMCBSIC90aHJlZSA0MiAwIFIgL3R3byA0MyAwIFIgL3UgNDQgMCBSIC96IDQ1IDAgUgovemVybyA0NiAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YxIDE1IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMCAvY2EgMSA+PgovQTIgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMSAvY2EgMSA+PgovQTMgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMC44IC9jYSAwLjggPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgPj4KZW5kb2JqCjIgMCBvYmoKPDwgL1R5cGUgL1BhZ2VzIC9LaWRzIFsgMTEgMCBSIF0gL0NvdW50IDEgPj4KZW5kb2JqCjQ3IDAgb2JqCjw8IC9DcmVhdG9yIChNYXRwbG90bGliIHYzLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZykKL1Byb2R1Y2VyIChNYXRwbG90bGliIHBkZiBiYWNrZW5kIHYzLjkuMikKL0NyZWF0aW9uRGF0ZSAoRDoyMDI0MDkyNzE0NTczMS0wNScwMCcpID4+CmVuZG9iagp4cmVmCjAgNDgKMDAwMDAwMDAwMCA2NTUzNSBmIAowMDAwMDAwMDE2IDAwMDAwIG4gCjAwMDAwMjA3NjEgMDAwMDAgbiAKMDAwMDAyMDUyNCAwMDAwMCBuIAowMDAwMDIwNTU2IDAwMDAwIG4gCjAwMDAwMjA2OTggMDAwMDAgbiAKMDAwMDAyMDcxOSAwMDAwMCBuIAowMDAwMDIwNzQwIDAwMDAwIG4gCjAwMDAwMDAwNjUgMDAwMDAgbiAKMDAwMDAwMDM0MSAwMDAwMCBuIAowMDAwMDEwNzkyIDAwMDAwIG4gCjAwMDAwMDAyMDggMDAwMDAgbiAKMDAwMDAxMDc3MCAwMDAwMCBuIAowMDAwMDE5MDg5IDAwMDAwIG4gCjAwMDAwMTg4ODIgMDAwMDAgbiAKMDAwMDAxODM4MSAwMDAwMCBuIAowMDAwMDIwMTQyIDAwMDAwIG4gCjAwMDAwMTA4MTIgMDAwMDAgbiAKMDAwMDAxMTEzMiAwMDAwMCBuIAowMDAwMDExMjgzIDAwMDAwIG4gCjAwMDAwMTE2NjMgMDAwMDAgbiAKMDAwMDAxMTgwOCAwMDAwMCBuIAowMDAwMDExOTQ5IDAwMDAwIG4gCjAwMDAwMTIyNTQgMDAwMDAgbiAKMDAwMDAxMjU1OCAwMDAwMCBuIAowMDAwMDEyODgwIDAwMDAwIG4gCjAwMDAwMTMzNDggMDAwMDAgbiAKMDAwMDAxMzU1NyAwMDAwMCBuIAowMDAwMDEzODc5IDAwMDAwIG4gCjAwMDAwMTQwNDUgMDAwMDAgbiAKMDAwMDAxNDI4MiAwMDAwMCBuIAowMDAwMDE0NDI2IDAwMDAwIG4gCjAwMDAwMTQ1ODEgMDAwMDAgbiAKMDAwMDAxNDgxNyAwMDAwMCBuIAowMDAwMDE1MjEyIDAwMDAwIG4gCjAwMDAwMTUzNjcgMDAwMDAgbiAKMDAwMDAxNTQ5MCAwMDAwMCBuIAowMDAwMDE1NzIzIDAwMDAwIG4gCjAwMDAwMTYxMzAgMDAwMDAgbiAKMDAwMDAxNjI3MiAwMDAwMCBuIAowMDAwMDE2NjY1IDAwMDAwIG4gCjAwMDAwMTY3NTUgMDAwMDAgbiAKMDAwMDAxNjk2MSAwMDAwMCBuIAowMDAwMDE3Mzc0IDAwMDAwIG4gCjAwMDAwMTc2OTggMDAwMDAgbiAKMDAwMDAxNzk0NSAwMDAwMCBuIAowMDAwMDE4MDkzIDAwMDAwIG4gCjAwMDAwMjA4MjEgMDAwMDAgbiAKdHJhaWxlcgo8PCAvU2l6ZSA0OCAvUm9vdCAxIDAgUiAvSW5mbyA0NyAwIFIgPj4Kc3RhcnR4cmVmCjIwOTc4CiUlRU9GCg==" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "for curve_idx, sty in zip(range(param_dict[\"num_states\"]), color_ls_alpha_cycler):\n", + " plt.plot(chi_ac_linspace, data_dict[\"avg_excitation\"][omega_d_idx, :, curve_idx],\n", + " label=curve_idx, **sty)\n", + "ax.set_xlabel(\"induced ac stark shift [GHz]\")\n", + "ax.legend(fontsize=12, ncol=2, loc='upper left', bbox_to_anchor=(1, 1))\n", + "plt.tight_layout()\n", + "plt.show()" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "We see that the scar plots and Blais plots both predict two locations (in drive power) where the first excited state can ionize. Lets look at the quasienergies to confirm that the quasienergies have an avoided crossing for the states experiencing branch crossings" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 35, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n \n \n \n \n 2024-09-27T15:04:06.729351\n image/svg+xml\n \n \n Matplotlib v3.9.2, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgNDI0LjQzNSAyNzkuNjI1MTY2MjkxMSBdIC9Db250ZW50cyA5IDAgUiAvQW5ub3RzIDEwIDAgUiA+PgplbmRvYmoKOSAwIG9iago8PCAvTGVuZ3RoIDEyIDAgUiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJy1W02TFDcSvdevqKN9oJBSX6mjWe+y4Rs2sT4QPhDDGOMFNvCYdcT++n1P6m5J1R9UDzNEjN3kVEv5pMx8L1Xi6fe3/313c/vj82fz336anra/3dxNdv4dP29nM/+On79mOz/Hz9vJ4G8fJi9+8S7g8/vDZ0l5iRJsjDCa8a+/TdOv09PvMMAdvvN8moIsNjgrYXZpCd6lgFFF8dmO1ve9VRK+pil7D3sbYjDvppI61Vs4DQCLAgImpmWKsigGNKmfuzP6xeymnp5hDf6aPuG/Zn5iMFawS7Lqs3dWZRZZcphvPkzPXk5P/2Fna+aXv5YlevlmejV/YxZjvp1/mV/+MP395fRiKn5M1oQlem+c7x3orZc8yB6rYUVsyCZv8UBOeOAVWxZTtIMHnfWSBxbfccZlrxpd3OKCP+FCtosIxsiDC531oguKUXOkry76LS7EYxeEUwhWMg5B2FkvuYBHFtXMDTPZbnFBT7iQ4qKW4TS40FkvuhARMiqI2qhWN7hgh2js19Pgl5G5lBdJyg/nx3n38c3nm9s38+ub+e7P13/8e7777d2v3842MQnD/M2f86vn//zfL8doW8YGXYxFCEfWEl3i2noSLQZJZYbgFzwlaoJL8chRXVzSbFB7Dh7j+ebNJ5QDw1FN+3DzYcYAT76//f31vz7/9Prj3ZMP7z5+vpu//8/84oT7ZblccOoG/zvzZQA5Ls5pMCkqnt0EwDwsAC+LBEWmD/4frJfdtw6FItrkvLq8xX236MO6r35xCHmXRv+b+QsAUD4Cqni0WXPahiA+KAKRuAQTDIpXj6AzX0bA6pMMsGbrjd+GwD8sggivRSWPMdSZv4Ag+kWd9TbiWbsNwddl8aeJHjzhs9YvKbOsWuORstbDp2COU7EVvU+fX9+9u/14+0dX6N6+u707rnWf5hPiRlxYyJWo5yB5i7mNifMft/PP88dZ5h9muwRKlKXUP42SAlbcx92fhN8kE1QxjuT5x+fzqNamn6YX89dMbDlbVrA65udsJvhsbI5GT8zWxBJ2G7xhIgt5ROyCDiEFY14wkHOCGEiGjO6im/EpiqhnZCS3ZCQqoiRZjCWhyMSUlmiSlxm/RnD7RKNCEmT4rnNCepscXLH6JWvgEyksMMWstCL4kg/ICFBoqpH1fsqCJUBdhgsJ31JILVojQGJiuADSw7rnTGtekrcBAgQ1ArLWKIuJNUgf0KxNM9yJ8LY4YQ1CAYvkw6xYkoSNYum3FlFuFJ7McNNkn8sYCDp8SkZmxZ5YzljMWELlFxXrJz4XN/BLjBdS0BnyN7mQbLECa/RJ4wykGco7lZEFAEA7FsyCCPIuJ668dcJnFDICsgqEY7RY4xKhXjEg1iUAuONqoIwsToxi9oxip1jPAsUDgMlZlMJTRQo5Y2cWPJDxlxwWfCfbMiHUsYMOV0+GKxFSxghhMcCNZzLCEvPZ4jXoHiuGvYIVa54L76AcLD6gAubycIrOlfVAcMEXg5WEWVDATTVnFJES8ZgxgZHKGAifQIXsaBXFbhXkid/Eelt6nVLMviDHRmdjpFglINAKFuwcokBcAZ4k0iOalQJJEK4wgz4q2VO35pigXOvqmVxCzGJEBIsTodlmxESxIj18wDJwB4AQgcw+BzEGMelSptknbGgsZgBD+iLZaA5MTJoRY94KG68ySPRSBkGQWWN9LI64IueLGUucjc9Ct0EvPnChsEaYKIdi1ZBDKsaAnEZuKRckMt98MSseQUhGLmoIAT0IzQiygDoFvNgY5EiqfiDKsNHYvTmjlGNHynwuL1nwNJDnsheFJoBioZBH0hUJhYyx1c5GQxLCBZiKNpRiD9ixbJOPxZ5tWUvaPdc+Ie+sQT5aoyVKoGlQAFDdZGYiR++0IEKosaa6wgLAlkNxHao6IccCvaHANtV3hBqWU5XOYAEUk5alRbAh/uDEzHKA2mNLxHYiHn4sXhBLpMSvrdRYHyxxLc4TKyrrHbD54TePQBGoaYUdUBscJoYA5zalGDhx5SxG1SW2YOSjeY15pAsx7BFRDnXgC5oDQj/HgTAEW47tldwTBo2lqo2MIQgD8JMtG9sog+aYTY2mxhnMQJQWxk9PGjQjICqVHEhDEDNgaStxYA2aEWjsxwbaEISSc8bW4tJ4g3aEHvJ2JA6mPWpLNivqoD0DfkHUuEMMky8iY0byoB2FVSshNPoQhqRWTB19iMlIM4SAG/mDdtQB7MvIIKxCJhpbvWwUQnuAuqtc1jgEKY+9C7a62ViEdjgfKgU0GkHCQ54hI/PII7QnRM/OvicStMILedzLyCS0g0OlFuVGJayWyBhXybxxCe2IAldESccm+A8iy4e6iY1OaKdckYFOxLLNQ7ULI5/QjjhEbI2EwiLtkkulvneMQju+uSOaA6MgUsEAAD9SCs1ZMmTOyClYlwUpYcohWUcqtOcYUrU3VhFIdfR3RTN0rELKMPizI6EDrdDOw746zIFX2KcYBwp2I7PQjv8zogdqEXGkJ1uZqHEL7SjQIn4kFxggsJAmMtIL7ZH9URz5pbgAaCXFO36hHdkeooz8wobMeo1eR36hHfM4WfGLCDI7OaTEQDA0E3i0I8EQOepAzfKOYWhPBi1KGBlGmKpoVeoqN4ahPQmyXB+GYsKuE4iZfYHZ1312IaXm+2Qeh2L2XUdSTOxc5TPBL3ZkgzJ8mWGQ6TYi3sPAMBDbCAeI6LEjoRnRi8AeGAZJjXKL5NCeYlicYSyB3xiGVoctLrnfGAYaBsUk1Qkbw9AsilZTB4ahQEiYr9SbRjE0W5SbEkyNYsgeEXRQ60SjmMIqBtUhjRRT5ErE9GmkGNpBJDWGO4qh6kFqxsJsjWIokpDSElb9Ce3Y0lS7mUYx9AwdfIoycgztqKwhrzimSDZsRpGEHcfQ7jiBjhxD6Ydd33ckB46h3WmOFVbjGEpLcErWVatCezJed/YDx1S7uCIpGsdweHxbK7U1jilKVF20ceQYwkWwm1KPO44pJw9WNI8Uw8VHs2Fr29kohpsFDWmiDhxTYgFDlXrZcUyJqJyDW3EMAzDnKpc7imHEoz2JwY8cwwTBWlYZ3ZEM0ywgU8tJW0cythys8GR+JBnQFCtqbY87krE8S8n1xLcjGQYdtzCGkWQQ0Dy3c6vmxYKqEiRULfaNZBikqI65Ng2NZGwhGSl6vOMYmrHHNaA6jimxLuiUVy0MvVQeg+ZVD0OwstOifQ+DRTM5qF33MGU93K7N6HsY8j6A+XUPg/2x2KG8amJQ+iBl0opjGJigQ1/kRscx2E2U1FwfbxTDZjyQOkeGYSMdXG2CH6KV2Jd5f3TMVayOofUYBNMOuaZHGD3tGWsamJKvBfZ/HgWV6Zdz1xACJRd650R4uFPDOumrGZUeTy5+/mU285upP7ucH//s8owXdr7mIPOhF6D134/fmZ/1orXp8zVtOqmB50LsDtXs391/GJrSRDGADn7dlCZbztOzH5tSFKEASq6ngF1TGjIFiei6KQ2ooRm1e92VoobmFGR1qCm+kNNOB3VNKT5mLGSpm31T6rBkJKGVYKCWd9Hvzg67phQeo266lV5g+4YdScaue1I8jjjPo16QQtyyI/quJ+WJRQrGrfUCCjB22ZaOqdcLPP3DdGHVkwIytICtHUovGNTxJV2oqq8JhpTZMKmJK8GQyNyxHlB2egFiGUTj40ovwF8PfjrSC8CH1Ra7OufkcggUnHUrvQC2c1jioil7vRBYQIxxutILHrvsY66U2AkGhAKkpjq3Egxga5C1jnIBWRWTWZ11MrCQYVKn7MVC+Wjt6rTTImoMmF3dWi2A/ZFFdZxOLQj3VXedaqcW2ALCnepPpxaEq6daD187teCwIw7CPq/UAqI7YbnriWqnFpgwztad7dQCagqmqDvbqwWWHGGxHtUC8ytIynlUC7wLkgS7Fke5wFfEWaIv0uuhC1zrTufH707Pe3FoVecrWtWHrvUH4j+i3Es6qjXMobxkQWkY+uUSPtGW9xKtXUbZZqSmPHTLsObyfqRvliMcQPSV3q81y5EFuLbVrVVmSTSuls3WKcdyklf7mtYoRxRtl+pBceuTYzm9r81La5PL5RgUxjB2ybG+Yc+r93cRiecQ03nskcv7TanJ1hEezEidHEa+Y7EW9IB2pDu+94xRd6/qDmwHzgQRpHp+18gusdvVo+6YKZclrc5fE9cji1m1xjAjZ6k8BqJLyHKne/478ByyNrBIjTSHrfQ+1SapYzkoAHQqsc54ILnEyIBfK47jLS3wXO3cG8UpLwHkekbWMZyygGVegBsITgWl8tAoH/it0BtJaqA33gABv4QVuynLfNC4YjdlVwroq25YWcbg6+rAlfcQ1dcXsY3aOBqqpY0js8GMnjHWA7xGbEr2NfWkuOM1vs3l5KtXeEwHhFY5C+hYDWaDEhfTQGtwz4Pd653JxmsAw1tGlWEbrSnf24Iz8shqMBsKRBlJrbzGZhM5chr2AGGh1dwoDRuZeZy1YjRlUcqybn8RO8j4etrU8RkCjccm9eC40VmimgrGpoHNKF7MLtc7MkMuqWgqCq7jssS2qEqFxmS8txJtPZTviCzyNQwI5Sve3Z0t3l2f8tg8kTp+euzG9KwTpiera7tUtrvmcJd3fOTkbeIz94Mx3Kl7xh/O3jPGF664rTw83Y1zcfyn37l6XfkH3rHGTzki2d24ngQlKNcvKt8o7cbzlheoett72GK5i91sKAlMYlYHu7Ioo4KWm2lvE9Bu0sNge2MopFC+W8cXFvO9qbrRmW56lzs7EiujPh/ZoS4OK9PmGqwHtzj2wd4gvB/MB6zdfG1JTq7nDS+OP9tfHC9nKwi241tZW082JoSllKsPFUHAwkWpO+f4zuDI/n4CyZy0I0JMif5daOzusXfX6kCgh/VG+Am/WXOyu9Emq+vA/aW17ibYNQcoPUIoft3HZo+w2UeEvf06hKw6/OaXENrTCO97TtOD5dGAPbGdzT6C7e3XgUUSGrthO+Uc2PsdB/VgoWqCntjZZh/B9varwBZZpRt21p0De7/WsANrqQlPJGpnH8AO9uvAUu9uSVR/Fux9OtAeKzSPO5GynX3E2tuvwxp9LclfwhrORvG9Gt0eLBRnOpGynX0E29uvA+vRUG1J2XgO7AVh2CNyvCV7YvuafUTU269DxL5hS16mc4jupT97rDac5M/OPmLt7ddh5THulrTUc1jvJ3M7sDmdZNJmHqD25quQ8qbqlpzMPdDT4n67RFr9M0MHZRdCaSkGhcRL3WszUFNCH5sL6or41d6fE8ixUHb35Y0Safz3UmegXyWYOryDXmp4B7nU8B6ppWvwbhZMdttWb5NP57d6UE8N+iCeGvQj7XQN9K3yycom6BvF1Hnog5Zq0Acp1aAfKakroG8WU9Ztgr5RWp2FPiqrA/RRWB2gH+uqa6BvlVbWb4N+rdCaRuSDzmrIB5nVkB+prGuQbxVaNmyL922y6/ymD6qrQR9EV4N+pLmugb5Vdtm4CfoFEdbhGzRYwzdIsIbvSIFdg2+rCLNpE76rJdkqqAdF1pAPgqwhP9Jj1yDfKsmsbkK+UaCdDepBnx2QD/LsAPxInV2Be6tAs4NCezH9H9HObjEKZW5kc3RyZWFtCmVuZG9iagoxMiAwIG9iago0MjA1CmVuZG9iagoxMCAwIG9iagpbIF0KZW5kb2JqCjE3IDAgb2JqCjw8IC9MZW5ndGggMjQ3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE1RSW7EMAy7+xX8wACWrMV5T4pBD+3/ryUdFO3BECNLXOLuxEQWXrZQ10KH48NGXgmbge+D1pz4GrHiP9pGpJU/VFsgEzFRJHRRNxr3SDe8CtF+pIJXqvdY8xF3K81bOnaxv/fBtOaRKqtCPOTYHNlIWtdE0fE9tN5zQ3TKIIE+NyEHRGmOXoWkv/bDdW00u7U2syeqg0emhPJJsxqa0ylmyGyox20qVjIKN6qMivtURloP8jbOMoCT44QyWk92rCai/NQnl5AXE3HCLjs7FmITCxuHtB+VPrH8fOvN+JtpraWQcUEiNMWl32e8x+d4/wCVT1wmCmVuZHN0cmVhbQplbmRvYmoKMTggMCBvYmoKPDwgL0xlbmd0aCA3OSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNzVSMFCwtAASZqYmCuZGlgophlxAPoiVy2VoaQ5m5YBZJsYGQJapqSkSCyIL0wthweRgtLGJOdQEBAskB7Y2B2ZbDlcGVxoA1pQcDAplbmRzdHJlYW0KZW5kb2JqCjE5IDAgb2JqCjw8IC9MZW5ndGggMzA3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2SS24DMQxD9z6FLhDA+tme86Qoupjef9snJemKHNkWRWqWukxZUx6QNJOEf+nwcLGd8jtsz2Zm4Fqil4nllOfQFWLuonzZzEZdWSfF6oRmOrfoUTkXBzZNqp+rLKXdLngO1yaeW/YRP7zQoB7UNS4JN3RXo2UpNGOq+3/Se/yMMuBqTF1sUqt7HzxeRFXo6AdHiSJjlxfn40EJ6UrCaFqIlXdFA0Hu8rTKewnu295qyLIHqZjOOylmsOt0Ui5uF4chHsjyqPDlo9hrQs/4sCsl9EjYhjNyJ+5oxubUyOKQ/t6NBEuPrmgh8+CvbtYuYLxTOkViZE5yrGmLVU73UBTTucO9DBD1bEVDKXOR1epfw84La5ZsFnhK+gUeo90mSw5W2duoTu+tPNnQ9x9a13QfCmVuZHN0cmVhbQplbmRvYmoKMjAgMCBvYmoKPDwgL0xlbmd0aCA3MyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwztjRQMFCwMFPQNTQ2VDCyNFYwNzNQSDHkAgqBWLlcMLEcMMvMEsQyNDdDYumaGUJlkVgg43K4YAbnwMzL4crgSgMAHokWlQplbmRzdHJlYW0KZW5kb2JqCjIxIDAgb2JqCjw8IC9MZW5ndGggNjkgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicM7Y0UDBQsDRX0DU0NlQwNjBRMDczUEgx5IIxc8EssGwOF0wdhGUGYhgZmiCxzIDGgSXhDJAZOXDTcrgyuNIA+qkWRQplbmRzdHJlYW0KZW5kb2JqCjIyIDAgb2JqCjw8IC9MZW5ndGggMjMyIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDVRSW7EMAy7+xX8wADW7rwnxaCH9v/XUsoUCEAltrglYmMjAi8x+DmI3PiSNaMmfmdyV/wsT4VHwq3gSRSBl+FedoLLG8ZlPw4zH7yXVs6kxpMMyEU2PTwRMtglEDowuwZ12Gbaib4h4bMjUs1GltPXEvTSKgTKU7bf6YISbav6c/usC2372hNOdnvqSeUTiOeWrMBl4xWTxVgGPVG5SzF9kOpsoSehvCifg2w+aohElyhn4InBwSjQDuy57WfiVSFoXd2nbWOoRkrH078NTU2SCPlECWe2NO4W/n/Pvb7X+w9OIVQRCmVuZHN0cmVhbQplbmRvYmoKMjMgMCBvYmoKPDwgL0xlbmd0aCAyMzEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNU85kgQhDMt5hT4wVRjbQL+np7Y22Pl/upKZTpDwIcnTEx2ZeJkjI7Bmx9taZCBm4FNMxb/2tA8TqvfgHiKUiwthhpFw1qzjbp6OF/92lc9YB+82+IpZXhDYwkzWVxZnLtsFY2mcxDnJboxdE7GNda2nU1hHMKEMhHS2w5Qgc1Sk9MmOMuboOJEnnovv9tssdjl+DusLNo0hFef4KnqCNoOi7HnvAhpyQf9d3fgeRbvoJSAbCRbWUWLunOWEX712dB61KBJzQppBLhMhzekqphCaUKyzo6BSUXCpPqforJ9/5V9cLQplbmRzdHJlYW0KZW5kb2JqCjI0IDAgb2JqCjw8IC9MZW5ndGggMjQ5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1QO45EIQzrOYUv8CTyI3AeRqstZu/frgOaKVBMfrYzJNARgUcMMZSv4yWtoK6Bv4tC8W7i64PCIKtDUiDOeg+IdOymNpETOh2cMz9hN2OOwEUxBpzpdKY9ByY5+8IKhHMbZexWSCeJqiKO6jOOKZ4qe594FiztyDZbJ5I95CDhUlKJyaWflMo/bcqUCjpm0QQsErngZBNNOMu7SVKMGZQy6h6mdiJ9rDzIozroZE3OrCOZ2dNP25n4HHC3X9pkTpXHdB7M+Jy0zoM5Fbr344k2B02N2ujs9xNpKi9Sux1anX51EpXdGOcYEpdnfxnfZP/5B/6HWiIKZW5kc3RyZWFtCmVuZG9iagoyNSAwIG9iago8PCAvTGVuZ3RoIDM5NSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9UktuxUAI2+cUXKDS8JvPeVJV3bz7b2tDUqkqvIkxxjB9ypC55UtdEnGFybderls8pnwuW1qZeYi7i40lPrbcl+4htl10LrE4HUfyCzKdKkSozarRofhCloUHkE7woQvCfTn+4y+AwdewDbjhPTJBsCTmKULGblEZmhJBEWHnkRWopFCfWcLfUe7r9zIFam+MpQtjHPQJtAVCbUjEAupAAETslFStkI5nJBO/Fd1nYhxg59GyAa4ZVESWe+zHiKnOqIy8RMQ+T036KJZMLVbGblMZX/yUjNR8dAUqqTTylPLQVbPQC1iJeRL2OfxI+OfWbCGGOm7W8onlHzPFMhLOYEs5YKGX40fg21l1Ea4dubjOdIEfldZwTLTrfsj1T/5021rNdbxyCKJA5U1B8LsOrkaxxMQyPp2NKXqiLLAamrxGM8FhEBHW98PIAxr9crwQNKdrIrRYIpu1YkSNimxzPb0E1kzvxTnWwxPCbO+d1qGyMzMqIYLauoZq60B2s77zcLafPzPoom0KZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvTGVuZ3RoIDEzNiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxNj0EOAzEIA+95hZ9AIEB4z1ZVD9v/X0vYdtMLHsmAbFEGgSWHeIcb4dHbD99FNhVn45xfUiliIZhPcJ8wUxyNKXfyY4+AcZRqLKdoeF5Lzk3DFy13Ey2lrZeTGW+47pf3R5VtkQ1Fzy0LQtdskvkygQd8GJhHdeNppcfd9myv9vwAzmw0SQplbmRzdHJlYW0KZW5kb2JqCjI3IDAgb2JqCjw8IC9MZW5ndGggMjQ5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE1RSYoDMAy75xX6QCFek7ynQ5lD5//Xyg6FOQQJr5KTlphYCw8xhB8sPfiRIXM3/Rt+otm7WXqSydn/mOciU1H4UqguYkJdiBvPoRHwPaFrElmxvfE5LKOZc74HH4W4BDOhAWN9STK5qOaVIRNODHUcDlqkwrhrYsPiWtE8jdxu+0ZmZSaEDY9kQtwYgIgg6wKyGCyUNjYTMlnOA+0NyQ1aYNepG1GLgiuU1gl0olbEqszgs+bWdjdDLfLgqH3x+mhWl2CF0Uv1WHhfhT6YqZl27pJCeuFNOyLMHgqkMjstK7V7xOpugfo/y1Lw/cn3+B2vD838XJwKZW5kc3RyZWFtCmVuZG9iagoyOCAwIG9iago8PCAvTGVuZ3RoIDk0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWNwRHAIAgE/1RBCQoK2k8mk4f2/40QMnxg5w7uhAULtnlGHwWVJl4VWAdKY9xQj0C94XItydwFD3Anf9rQVJyW03dpkUlVKdykEnn/DmcmkKh50WOd9wtj+yM8CmVuZHN0cmVhbQplbmRvYmoKMjkgMCBvYmoKPDwgL0xlbmd0aCAzNDEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRVJLbkQxCNu/U3CBSOGXkPO0qrqY3n9bm0zVzeAJYGx4y1OmZMqwuSUjJNeUT30iQ6ym/DRyJCKm+EkJBXaVj8drS6yN7JGoFJ/a8eOx9Eam2RVa9e7Rpc2iUc3KyDnIEKGeFbqye9QO2fB6XEi675TNIRzL/1CBLGXdcgolQVvQd+wR3w8droIrgmGway6D7WUy1P/6hxZc7333YscugBas577BDgCopxO0BcgZ2u42KWgAVbqLScKj8npudqJso1Xp+RwAMw4wcsCIJVsdvtHeAJZ9XehFjYr9K0BRWUD8yNV2wd4xyUhwFuYGjr1wPMWZcEs4xgJAir3iGHrwJdjmL1euiJrwCXW6ZC+8wp7a5udCkwh3rQAOXmTDraujqJbt6TyC9mdFckaM1Is4OiGSWtI5guLSoB5a41w3seJtI7G5V9/uH+GcL1z26xdL7ITECmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0xlbmd0aCAxNjQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRZDHcQUxDEPvqgIlMIAK9azH8w/r/q+G9NNBehhCDGJPwrBcV3FhdMOPty0zDX9HGe7G+jJjvNVYICfoAwyRiavRpPp2xRmq9OTVYq6jolwvOiISzJLjq0AjfDqyx5O2tjP9dF4f7CHvE/8qKuduYQEuqu5A+VIf8dSP2VHqmqGPKitrHmraV4RdEUrbPi6nMk7dvQNa4b2Vqz3a7z8edjryCmVuZHN0cmVhbQplbmRvYmoKMzEgMCBvYmoKPDwgL0xlbmd0aCA3MiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzMrdQMFCwNAEShhYmCuZmBgophlxAvqmJuUIuF0gMxMoBswyAtCWcgohngJggbRDFIBZEsZmJGUQdnAGRy+BKAwAl2xbJCmVuZHN0cmVhbQplbmRvYmoKMzIgMCBvYmoKPDwgL0xlbmd0aCA4MyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9zDkSgDAIBdCeU/wjhMgi93Eci3j/VjDRBh6reqAhOIO6wa3hYMq6dBPvU+PVxpwSCah4Sk2Wugt61LS+1L5o4Lvr5kvViT/NzxedD7sdGd0KZW5kc3RyZWFtCmVuZG9iagozMyAwIG9iago8PCAvVHlwZSAvWE9iamVjdCAvU3VidHlwZSAvRm9ybSAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0xlbmd0aCAzOQovRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJzjMjQwUzA2NVXI5TI3NgKzcsAsI3MjIAski2BBZDO40gAV8wp8CmVuZHN0cmVhbQplbmRvYmoKMzQgMCBvYmoKPDwgL0xlbmd0aCAxNjMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRZA7EgMhDEN7TqEj+CMDPs9mMik2929j2GxSwNNYIIO7E4LU2oKJ6IKHtiXdBe+tBGdj/Ok2bjUS5AR1gFak42iUUn25xWmVdPFoNnMrC60THWYOepSjGaAQOhXe7aLkcqbuzvlDcPVf9b9i3TmbiYHJyh0IzepT3Pk2O6K6usn+pMfcrNd+K+xVYWlZS8sJt527ZkAJ3FM52qs9Px8KOvYKZW5kc3RyZWFtCmVuZG9iagozNSAwIG9iago8PCAvTGVuZ3RoIDMyMiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UbttxTAM7DUFFzAgfiXN4yBIkbd/mzvaqUjTvB9VXjKlXC51ySpZYfKlQ3WKpnyeZqb8DvWQ45ge2SG6U9aWexgWlol5Sh2xmiz3cAs2vgCaEnML8fcI8CuAUcBEoG7x9w+6WRJAGhT8FOiaq5ZYYgINi4Wt2RXiVt0pWLir+HYkuQcJcjFZ6FMORYopt8B8GSzZkVqc63JZCv9ufQIaYYU47LOLROB5wANMJP5kgGzPPlvs6upFNnaGOOnQgIuAm80kAUFTOKs+uGH7arvm55koJzg51q+iMb4NTuZLUt5XucfPoEHe+DM8Z3eOUA6aUAj03QIgh93ARoQ+tc/ALgO2Sbt3Y0r5nGQpvgQ2CvaoUx3K8GLszFZv2PzH6MpmUWyQlfXR6Q7K3KATYh5vZKFbsrb7Nw+zff8BXxl7ZAplbmRzdHJlYW0KZW5kb2JqCjM2IDAgb2JqCjw8IC9MZW5ndGggODMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRYy7DcAwCER7pmAEfib2PlGUwt6/DRAlbrgn3T1cHQmZKW4zw0MGngwshl1xgfSWMAtcR1COneyjYdW+6gSN9aZS8+8PlJ7srOKG6wECQhpmCmVuZHN0cmVhbQplbmRvYmoKMzcgMCBvYmoKPDwgL0xlbmd0aCA1MSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNrRQMFAwNDAHkkaGQJaRiUKKIRdIAMTM5YIJ5oBZBkAaojgHriaHK4MrDQDhtA2YCmVuZHN0cmVhbQplbmRvYmoKMzggMCBvYmoKPDwgL0xlbmd0aCAyNDMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicTVG7rQMxDOs9hRY4wPrZvnkueHjFZf82pJwEqURDFEnJw1O6ZMphfUpGSI4uD20aS2y6PDdCU4eKgqlrieqUq5mmzFMsTdDz3lmu5hjge1U31N/0iF4CkVGCVWGBDpA7uGD42WsmbFELIjGGUDOAacIKc7gSMQQZjLVnGJQqDE7VzypX+y+nZdgqsHgwnSI/sppop1+6HHjrKQdC2NyVu3ohTQjujQZjzCxcd6mynQAcTHSZiYxYvA3H0yEMDV6aBqxw1o2YILEbI6UPXgcZ07B3RR51txjxvlvGlLvVz31RfeZd7R8IwRsn+HsByhtdXgplbmRzdHJlYW0KZW5kb2JqCjM5IDAgb2JqCjw8IC9MZW5ndGggMTYwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWQORIDMQgEc72CJ0hcgvesy7XB+v+pB9ZHoukCNBy6Fk3KehRoPumxRqG60GvoLEqSRMEWkh1Qp2OIOyhITEhjkki2HoMjmlizXZiZVCqzUuG0acXCv9la1chEjXCN/InpBlT8T+pclPBNg6+SMfoYVLw7g4xJ+F5F3Fox7f5EMLEZ9glvRSYFhImxqdm+z2CGzPcK1zjH8w1MgjfrCmVuZHN0cmVhbQplbmRvYmoKNDAgMCBvYmoKPDwgL0xlbmd0aCAzMzQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicLVJLcsUgDNtzCl2gM/gH5DzpdLp4vf+2kpNFRg5g9DHlholKfFkgt6PWxLeNzECF4a+rzIXPSNvIOojLkIu4ki2Fe0Qs5DHEPMSC76vxHh75rMzJswfGL9l3Dyv21IRlIePFGdphFcdhFeRYsHUhqnt4U6TDqSTY44v/PsVzLQQtfEbQgF/kn6+O4PmSFmn3mG3TrnqwTDuqpLAcbE9zXiZfWme5Oh7PB8n2rtgRUrsCFIW5M85z4SjTVka0FnY2SGpcbG+O/VhK0IVuXEaKI5CfqSI8oKTJzCYK4o+cHnIqA2Hqmq50chtVcaeezDWbi7czSWbrvkixmcJ5XTiz/gxTZrV5J89yotSpCO+xZ0vQ0Dmunr2WWWh0mxO8pITPxk5PTr5XM+shORUJqWJaV8FpFJliCdsSX1NRU5p6Gf778u7xO37+ASxzfHMKZW5kc3RyZWFtCmVuZG9iago0MSAwIG9iago8PCAvTGVuZ3RoIDcwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMzNlMwULAwAhKmpoYK5kaWCimGXEA+iJXLBRPLAbPMLMyBLCMLkJYcLkMLYzBtYmykYGZiBmRZIDEgujK40gCYmhMDCmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0xlbmd0aCAzMjAgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVJLbgUxCNvPKbhApfBPzvOqqou++29rE70VTDBg4ykvWdJLvtQl26XD5Fsf9yWxQt6P7ZrMUsX3FrMUzy2vR88Rty0KBFETPViZLxUi1M/06DqocEqfgVcItxQbvINJAINq+AcepTMgUOdAxrtiMlIDgiTYc2lxCIlyJol/pLye3yetpKH0PVmZy9+TS6XQHU1O6AHFysVJoF1J+aCZmEpEkpfrfbFC9IbAkjw+RzHJgOw2iW2iBSbnHqUlzMQUOrDHArxmmtVV6GDCHocpjFcLs6gebPJbE5WkHa3jGdkw3sswU2Kh4bAF1OZiZYLu5eM1r8KI7VGTXcNw7pbNdwjRaP4bFsrgYxWSgEensRINaTjAiMCeXjjFXvMTOQ7AiGOdmiwMY2gmp3qOicDQnrOlYcbHHlr18w9U6XyHCmVuZHN0cmVhbQplbmRvYmoKNDMgMCBvYmoKPDwgL0xlbmd0aCAxOCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNrRQMIDDFEOuNAAd5gNSCmVuZHN0cmVhbQplbmRvYmoKNDQgMCBvYmoKPDwgL0xlbmd0aCAxMzMgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRY9LDgQhCET3nKKOwMcf53Ey6YVz/+2AnW4TYz2FVIG5gqE9LmsDnRUfIRm28beplo5FWT5UelJWD8ngh6zGyyHcoCzwgkkqhiFQi5gakS1lbreA2zYNsrKVU6WOsIujMI/2tGwVHl+iWyJ1kj+DxCov3OO6Hcil1rveoou+f6QBMQkKZW5kc3RyZWFtCmVuZG9iago0NSAwIG9iago8PCAvTGVuZ3RoIDM0MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1UjluBDEM6/0KfSCAbtvv2SBIkfy/DanZFANxdFKUO1pUdsuHhVS17HT5tJXaEjfkd2WFxAnJqxLtUoZIqLxWIdXvmTKvtzVnBMhSpcLkpORxyYI/w6WnC8f5trGv5cgdjx5YFSOhRMAyxcToGpbO7rBmW36WacCPeIScK9Ytx1gFUhvdOO2K96F5LbIGiL2ZlooKHVaJFn5B8aBHjX32GFRYINHtHElwjIlQkYB2gdpIDDl7LHZRH/QzKDET6NobRdxBgSWSmDnFunT03/jQsaD+2Iw3vzoq6VtaWWPSPhvtlMYsMul6WPR089bHgws076L859UMEjRljZLGB63aOYaimVFWeLdDkw3NMcch8w6ewxkJSvo8FL+PJRMdlMjfDg2hf18eo4ycNt4C5qI/bRUHDuKzw165gRVKF2uS9wGpTOiB6f+v8bW+19cfHe2AxgplbmRzdHJlYW0KZW5kb2JqCjQ2IDAgb2JqCjw8IC9MZW5ndGggMjUxIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nC1RSXIDQQi7zyv0hGan32OXK4fk/9cIygcGDYtAdFrioIyfICxXvOWRq2jD3zMxgt8Fh34r121Y5EBUIEljUDWhdvF69B7YcZgJzJPWsAxmrA/8jCnc6MXhMRlnt9dl1BDsXa89mUHJrFzEJRMXTNVhI2cOP5kyLrRzPTcg50ZYl2GQblYaMxKONIVIIYWqm6TOBEESjK5GjTZyFPulL490hlWNqDHscy1tX89NOGvQ7Fis8uSUHl1xLicXL6wc9PU2AxdRaazyQEjA/W4P9XOyk994S+fOFtPje83J8sJUYMWb125ANtXi37yI4/uMr+fn+fwDX2BbiAplbmRzdHJlYW0KZW5kb2JqCjQ3IDAgb2JqCjw8IC9MZW5ndGggMTc0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE2QSQ5DIQxD95zCF6iEM8DnPL+qumjvv61DB3WB/OQgcDw80HEkLnRk6IyOK5sc48CzIGPi0Tj/ybg+xDFB3aItWJd2x9nMEnPCMjECtkbJ2TyiwA/HXAgSZJcfvsAgIl2P+VbzWZP0z7c73Y+6tGZfPaLAiewIxbABV4D9useBS8L5XtPklyolYxOH8oHqIlI2O6EQtVTscqqKs92bK3AV9PzRQ+7tBbUjPN8KZW5kc3RyZWFtCmVuZG9iago0OCAwIG9iago8PCAvTGVuZ3RoIDc2IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2MOw6AMAxD95zCR2h+JAdCiIHef6UptIv99CTbxdFgWpECt8DJ5D6p03LPJDt8EJsh5FcbWrWuytKaDIuajL8N391N1wumOBfACmVuZHN0cmVhbQplbmRvYmoKNDkgMCBvYmoKPDwgL0xlbmd0aCAyMTUgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNVE5DgMhDOz3Ff5AJIwveE+iKM3+v82M0VYewVyGtJQhmfJSk6gh5VM+epkunLrc18xqNOeWtC1zgLi2vC+tksCJZoiDwWmYuAGaPAFD19GoUUMXHtDUpVMosNwEPoq3bg/dY7WBl7Yh54kgYigZLEHNqUUTFm3PJ6Q1v16LG96X7d3IU6XGlhiBBgFWOBzX6NfwlT1PJtF0FTLUqzXLGAkTRSI8+Y6m1RPrWjTSMhLUxhGsagO8O/0wTgAAE3HLAmSfSpSz5MRvsfSzBlf6/gGfR1SWCmVuZHN0cmVhbQplbmRvYmoKMTUgMCBvYmoKPDwgL1R5cGUgL0ZvbnQgL0Jhc2VGb250IC9CTVFRRFYrRGVqYVZ1U2FucyAvRmlyc3RDaGFyIDAgL0xhc3RDaGFyIDI1NQovRm9udERlc2NyaXB0b3IgMTQgMCBSIC9TdWJ0eXBlIC9UeXBlMyAvTmFtZSAvQk1RUURWK0RlamFWdVNhbnMKL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXQovQ2hhclByb2NzIDE2IDAgUgovRW5jb2RpbmcgPDwgL1R5cGUgL0VuY29kaW5nCi9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0NiAvcGVyaW9kIDQ4IC96ZXJvIC9vbmUgL3R3byAvdGhyZWUgL2ZvdXIgL2ZpdmUgL3NpeCAvc2V2ZW4KL2VpZ2h0IC9uaW5lIDcxIC9HIC9IIDkxIC9icmFja2V0bGVmdCA5MyAvYnJhY2tldHJpZ2h0IDk3IC9hIDk5IC9jIC9kIC9lIC9mCi9nIC9oIC9pIDEwNyAvayAxMTAgL24gMTEzIC9xIC9yIC9zIC90IC91IDEyMiAveiBdCj4+Ci9XaWR0aHMgMTMgMCBSID4+CmVuZG9iagoxNCAwIG9iago8PCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL0ZvbnROYW1lIC9CTVFRRFYrRGVqYVZ1U2FucyAvRmxhZ3MgMzIKL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvQXNjZW50IDkyOSAvRGVzY2VudCAtMjM2IC9DYXBIZWlnaHQgMAovWEhlaWdodCAwIC9JdGFsaWNBbmdsZSAwIC9TdGVtViAwIC9NYXhXaWR0aCAxMzQyID4+CmVuZG9iagoxMyAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzQyIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjMgNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjEyIDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTIgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwNQo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTgyIDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoxNiAwIG9iago8PCAvRyAxNyAwIFIgL0ggMTggMCBSIC9hIDE5IDAgUiAvYnJhY2tldGxlZnQgMjAgMCBSIC9icmFja2V0cmlnaHQgMjEgMCBSCi9jIDIyIDAgUiAvZCAyMyAwIFIgL2UgMjQgMCBSIC9laWdodCAyNSAwIFIgL2YgMjYgMCBSIC9maXZlIDI3IDAgUgovZm91ciAyOCAwIFIgL2cgMjkgMCBSIC9oIDMwIDAgUiAvaSAzMSAwIFIgL2sgMzIgMCBSIC9uIDM0IDAgUiAvbmluZSAzNSAwIFIKL29uZSAzNiAwIFIgL3BlcmlvZCAzNyAwIFIgL3EgMzggMCBSIC9yIDM5IDAgUiAvcyA0MCAwIFIgL3NldmVuIDQxIDAgUgovc2l4IDQyIDAgUiAvc3BhY2UgNDMgMCBSIC90IDQ0IDAgUiAvdGhyZWUgNDUgMCBSIC90d28gNDYgMCBSIC91IDQ3IDAgUgoveiA0OCAwIFIgL3plcm8gNDkgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAxNSAwIFIgPj4KZW5kb2JqCjQgMCBvYmoKPDwgL0ExIDw8IC9UeXBlIC9FeHRHU3RhdGUgL0NBIDAgL2NhIDEgPj4KL0EyIDw8IC9UeXBlIC9FeHRHU3RhdGUgL0NBIDEgL2NhIDEgPj4KL0EzIDw8IC9UeXBlIC9FeHRHU3RhdGUgL0NBIDAuOCAvY2EgMC44ID4+ID4+CmVuZG9iago1IDAgb2JqCjw8ID4+CmVuZG9iago2IDAgb2JqCjw8ID4+CmVuZG9iago3IDAgb2JqCjw8IC9GMS1EZWphVnVTYW5zLW1pbnVzIDMzIDAgUiA+PgplbmRvYmoKMiAwIG9iago8PCAvVHlwZSAvUGFnZXMgL0tpZHMgWyAxMSAwIFIgXSAvQ291bnQgMSA+PgplbmRvYmoKNTAgMCBvYmoKPDwgL0NyZWF0b3IgKE1hdHBsb3RsaWIgdjMuOS4yLCBodHRwczovL21hdHBsb3RsaWIub3JnKQovUHJvZHVjZXIgKE1hdHBsb3RsaWIgcGRmIGJhY2tlbmQgdjMuOS4yKQovQ3JlYXRpb25EYXRlIChEOjIwMjQwOTI3MTUwNDA2LTA1JzAwJykgPj4KZW5kb2JqCnhyZWYKMCA1MQowMDAwMDAwMDAwIDY1NTM1IGYgCjAwMDAwMDAwMTYgMDAwMDAgbiAKMDAwMDAxNTU2NyAwMDAwMCBuIAowMDAwMDE1MzAyIDAwMDAwIG4gCjAwMDAwMTUzMzQgMDAwMDAgbiAKMDAwMDAxNTQ3NiAwMDAwMCBuIAowMDAwMDE1NDk3IDAwMDAwIG4gCjAwMDAwMTU1MTggMDAwMDAgbiAKMDAwMDAwMDA2NSAwMDAwMCBuIAowMDAwMDAwMzQ1IDAwMDAwIG4gCjAwMDAwMDQ2NDYgMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDA0NjI1IDAwMDAwIG4gCjAwMDAwMTM4NDcgMDAwMDAgbiAKMDAwMDAxMzY0MCAwMDAwMCBuIAowMDAwMDEzMTM3IDAwMDAwIG4gCjAwMDAwMTQ5MDAgMDAwMDAgbiAKMDAwMDAwNDY2NiAwMDAwMCBuIAowMDAwMDA0OTg2IDAwMDAwIG4gCjAwMDAwMDUxMzcgMDAwMDAgbiAKMDAwMDAwNTUxNyAwMDAwMCBuIAowMDAwMDA1NjYyIDAwMDAwIG4gCjAwMDAwMDU4MDMgMDAwMDAgbiAKMDAwMDAwNjEwOCAwMDAwMCBuIAowMDAwMDA2NDEyIDAwMDAwIG4gCjAwMDAwMDY3MzQgMDAwMDAgbiAKMDAwMDAwNzIwMiAwMDAwMCBuIAowMDAwMDA3NDExIDAwMDAwIG4gCjAwMDAwMDc3MzMgMDAwMDAgbiAKMDAwMDAwNzg5OSAwMDAwMCBuIAowMDAwMDA4MzEzIDAwMDAwIG4gCjAwMDAwMDg1NTAgMDAwMDAgbiAKMDAwMDAwODY5NCAwMDAwMCBuIAowMDAwMDA4ODQ5IDAwMDAwIG4gCjAwMDAwMDkwMjEgMDAwMDAgbiAKMDAwMDAwOTI1NyAwMDAwMCBuIAowMDAwMDA5NjUyIDAwMDAwIG4gCjAwMDAwMDk4MDcgMDAwMDAgbiAKMDAwMDAwOTkzMCAwMDAwMCBuIAowMDAwMDEwMjQ2IDAwMDAwIG4gCjAwMDAwMTA0NzkgMDAwMDAgbiAKMDAwMDAxMDg4NiAwMDAwMCBuIAowMDAwMDExMDI4IDAwMDAwIG4gCjAwMDAwMTE0MjEgMDAwMDAgbiAKMDAwMDAxMTUxMSAwMDAwMCBuIAowMDAwMDExNzE3IDAwMDAwIG4gCjAwMDAwMTIxMzAgMDAwMDAgbiAKMDAwMDAxMjQ1NCAwMDAwMCBuIAowMDAwMDEyNzAxIDAwMDAwIG4gCjAwMDAwMTI4NDkgMDAwMDAgbiAKMDAwMDAxNTYyNyAwMDAwMCBuIAp0cmFpbGVyCjw8IC9TaXplIDUxIC9Sb290IDEgMCBSIC9JbmZvIDUwIDAgUiA+PgpzdGFydHhyZWYKMTU3ODQKJSVFT0YK" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "for curve_idx, sty in zip(range(param_dict[\"num_states\"]), color_ls_alpha_cycler):\n", + " plt.plot(chi_ac_linspace, data_dict[\"all_quasienergies\"][omega_d_idx, :, curve_idx]/2/np.pi,\n", + " label=curve_idx, **sty)\n", + "ax.set_xlabel(\"induced ac stark shift [GHz]\")\n", + "ax.set_ylabel(\"quasienergies [GHz]\")\n", + "ax.legend(fontsize=12, ncol=2, loc='upper left', bbox_to_anchor=(1, 1))\n", + "ax.set_ylim(-27/2/np.pi, -20/2/np.pi)\n", + "plt.tight_layout()\n", + "plt.show()" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [], + "metadata": { + "collapsed": false + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/floquet.md b/docs/floquet.md new file mode 100644 index 0000000..090d566 --- /dev/null +++ b/docs/floquet.md @@ -0,0 +1,24 @@ +# Python API + +The **floquet** Python API consists of the **FloquetAnalysis** class function **optimize** which performs optimal control, and various models +that can be optimized including **SESolveModel** and **MESolveModel**. There are additionally various utility functions +and classes that help to define an optimization routine (cost functions, file input output, options, etc.) + +## Optimization + +::: floquet.floquet + options: + show_source: false + +## Amplitude conversion utilities + +::: floquet.amplitude_converters + options: + show_source: false + +## File utilities + +::: floquet.utils.file_io + options: + show_source: false + diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..fa27671 --- /dev/null +++ b/docs/index.md @@ -0,0 +1,31 @@ +# Getting started + +**floquet** is a python package for performing floquet simulations on quantum systems to identify nonlinear resonances. + +## Installation + +For now we support only installing directly from github +```bash +pip install git+https://github.com/dkweiss31/floquet +``` + +Requires Python 3.10+ + +## Jump in + +Please see the demo in the transmon notebook on the left for an example of how to set up a floquet simulation on a transmon. There we plot the "scar" plots which allow for "birds-eye view" of ionization. We then also plot the [Blais branch crossing results](https://arxiv.org/abs/2402.06615) to understand which states are responsible for the ionization. + +## Citation + +If you found this package useful in academic work, please cite + +```bibtex +@unpublished{floquet2024, + title = {Floquet: Identifying nonlinear resonances in quantum systems due to parametric drives}, + author = {Daniel K. Weiss}, + year = {2024}, + url = {https://github.com/dkweiss31/floquet} +} +``` + +Also please consider starring the project on [github](https://github.com/dkweiss31/floquet/)! \ No newline at end of file diff --git a/floquet/__init__.py b/floquet/__init__.py new file mode 100644 index 0000000..5008167 --- /dev/null +++ b/floquet/__init__.py @@ -0,0 +1,14 @@ +from importlib.metadata import version + +from .amplitude_converters import ChiacToAmp as ChiacToAmp +from .amplitude_converters import XiSqToAmp as XiSqToAmp +from .floquet import floquet_analysis as floquet_analysis +from .floquet import floquet_analysis_from_file as floquet_analysis_from_file +from .options import Options +from .utils.file_io import extract_info_from_h5 as extract_info_from_h5 +from .utils.file_io import generate_file_path as generate_file_path +from .utils.file_io import update_data_in_h5 as update_data_in_h5 +from .utils.file_io import write_to_h5 as write_to_h5 +from .utils.parallel import parallel_map as parallel_map + +__version__ = version(__package__) diff --git a/floquet/amplitude_converters.py b/floquet/amplitude_converters.py new file mode 100644 index 0000000..92fba9b --- /dev/null +++ b/floquet/amplitude_converters.py @@ -0,0 +1,113 @@ +import numpy as np +import qutip as qt + + +class ChiacToAmp: + r"""Convert given induced ac-stark shift values to drive amplitudes. + + Consider a qubit coupled to an oscillator with the interaction Hamiltonian + $H_I = g(ab^{\dagger} + a^{\dagger}b)$. If the oscillator is driven to + an average occupation number of $\bar{n}$, then the effective drive strength + seen by the qubit is $\Omega_d = 2 g \sqrt{\bar{n}}$. On the other hand based + on a Schrieffer-Wolff transformation, the interaction hamiltonian is + $H^{(2)} = \chi a^{\dagger}ab^{\dagger}b$. The average induced + ac-stark shift is then $\chi_{ac} = \chi \bar{n}$. Thus $\Omega_d = 2g\sqrt{\chi_{\rm ac}/\chi}$. + Observe that since $\chi \sim g^2$, $g$ effectively cancels out and can be set to 1 + (still need 2 pi out front since result needs to have units of radial frequency to + be a drive strength. So really we are taking $g = 2 \pi g'$ where $g$ has the right units + and $g'$ has units of GHz and we are setting $g'=1$) + """ # noqa E501 + + def __init__( + self, + H0: qt.Qobj, + H1: qt.Qobj, + state_indices: list, + omega_d_linspace: np.ndarray, + ): + self.H0 = H0 + self.H1 = H1 + self.state_indices = state_indices + self.omega_d_linspace = omega_d_linspace + + def amplitudes_for_omega_d(self, chi_ac_linspace: np.ndarray) -> np.ndarray: + if isinstance(chi_ac_linspace, float): + chi_ac_linspace = np.array([chi_ac_linspace]) + chis_for_omega_d = self.compute_chis_for_omega_d() + return np.einsum( + 'a,w->aw', + 2.0 * np.sqrt(chi_ac_linspace), + 1.0 / np.sqrt(chis_for_omega_d), + ) + + def compute_chis_for_omega_d(self) -> np.ndarray: + chi_0 = self.chi_ell( + np.diag(self.H0.full()), + self.H1.full(), + self.omega_d_linspace, + self.state_indices[0], + ) + chi_1 = self.chi_ell( + np.diag(self.H0.full()), + self.H1.full(), + self.omega_d_linspace, + self.state_indices[1], + ) + return np.abs(chi_1 - chi_0) + + @staticmethod + def chi_ell_ellp( + qubit_energies: np.ndarray, + op_mat: np.ndarray, + E_osc: float, + ell: int, + ellp: int, + ) -> np.ndarray: + E_ell_ellp = qubit_energies[ell] - qubit_energies[ellp] + return np.abs(op_mat[ell, ellp]) ** 2 / (E_ell_ellp - E_osc) + + def chi_ell( + self, qubit_energies: np.ndarray, op_mat: np.ndarray, E_osc: float, ell: int + ) -> np.ndarray: + n = len(qubit_energies) + return sum( + self.chi_ell_ellp(qubit_energies, op_mat, E_osc, ell, ellp) + - self.chi_ell_ellp(qubit_energies, op_mat, E_osc, ellp, ell) + for ellp in range(n) + ) + + +class XiSqToAmp: + r"""convert given $|\xi|^2$ value into a drive amplitude. + + This is based on the equivalence $\xi = 2 \Omega_d \omega_d / (\omega_d^2-\omega^2), + where in this definition $|\xi|^2= 2 \chi_{\rm ac} / \alpha$ where $\chi_{\rm ac}$ is + the induced ac stark shift, $\alpha$ is the anharmonicity and $\Omega_d$ is the + drive amplitude. + """ # noqa E501 + + def __init__( + self, + H0: qt.Qobj, + H1: qt.Qobj, + state_indices: list, + omega_d_linspace: np.ndarray, + ): + self.H0 = H0 + self.H1 = H1 + self.state_indices = state_indices + self.omega_d_linspace = omega_d_linspace + + def amplitudes_for_omega_d(self, xi_sq_linspace: np.ndarray) -> np.ndarray: + if isinstance(xi_sq_linspace, float): + xi_sq_linspace = np.array([xi_sq_linspace]) + idx_0 = self.state_indices[0] + idx_1 = self.state_indices[1] + drive_matelem = self.H1[idx_0, idx_1] + omega_01 = self.H0[idx_1, idx_1] - self.H0[idx_0, idx_0] + return np.einsum( + 'x,w->xw', + np.sqrt(xi_sq_linspace) / np.abs(drive_matelem), + np.abs(omega_01**2 - self.omega_d_linspace**2) + / (2 * self.omega_d_linspace), + ) diff --git a/floquet/floquet.py b/floquet/floquet.py new file mode 100644 index 0000000..b189b33 --- /dev/null +++ b/floquet/floquet.py @@ -0,0 +1,799 @@ +from __future__ import annotations + +import ast +import functools +import itertools +import time +import warnings +from itertools import chain, product + +import numpy as np +import qutip as qt +import scipy as sp + +from .options import Options +from .utils.file_io import extract_info_from_h5, update_data_in_h5, write_to_h5 +from .utils.parallel import parallel_map + + +def floquet_analysis( + H0: qt.Qobj | np.ndarray | list, + H1: qt.Qobj | np.ndarray | list, + state_indices: list | None = None, + omega_d_linspace: np.ndarray | list = 2.0 * np.pi * np.linspace(6.9, 13, 50), # noqa B008 + amp_linspace: np.ndarray | list = 2.0 * np.pi * np.linspace(0.0, 0.2, 51), # noqa B008 + options: Options = Options(), # noqa B008 + init_data_to_save: dict | None = None, +) -> FloquetAnalysis: + """Perform a floquet analysis to identify nonlinear resonances. + + Arguments: + H0: qt.Qobj | np.ndarray | list + Drift Hamiltonian (ideally diagonal) + H1: qt.Qobj | np.ndarray | list + drive operator + state_indices: list + state indices of interest + omega_d_linspace: ndarray | list + drive frequencies to scan over + amp_linspace: ndarray | list + amp values to scan over. Can be one dimensional in which case + these amplitudes are used for all omega_d, or it can be two dimensional + in which case the first dimension are the amplitudes to scan over + and the second are the amplitudes for respective drive frequencies + options: Options + Options for the Floquet analysis. + init_data_to_save: dict | None + Initial parameter metadata to save to file. Defaults to None. + """ + if state_indices is None: + state_indices = [0, 1] + if not isinstance(H0, qt.Qobj): + H0 = qt.Qobj(np.array(H0, dtype=complex)) + if not isinstance(H1, qt.Qobj): + H1 = qt.Qobj(np.array(H1, dtype=complex)) + if isinstance(omega_d_linspace, list): + omega_d_linspace = np.array(omega_d_linspace) + if isinstance(amp_linspace, list): + amp_linspace = np.array(amp_linspace) + if len(amp_linspace.shape) == 1: + amp_linspace = np.tile(amp_linspace, (len(omega_d_linspace), 1)).T + else: + assert len(amp_linspace.shape) == 2 + assert amp_linspace.shape[1] == len(omega_d_linspace) + return FloquetAnalysis( + H0, + H1, + state_indices, + omega_d_linspace, + amp_linspace, + options, + init_data_to_save=init_data_to_save, + ) + + +def floquet_analysis_from_file(filepath: str) -> FloquetAnalysis: + _, param_dict = extract_info_from_h5(filepath) + floquet_init = ast.literal_eval(param_dict['floquet_analysis_init']) + return floquet_analysis( + floquet_init['H0'], + floquet_init['H1'], + state_indices=floquet_init['state_indices'], + omega_d_linspace=floquet_init['omega_d_linspace'], + amp_linspace=floquet_init['amp_linspace'], + options=Options(**floquet_init['options']), + init_data_to_save=floquet_init['init_data_to_save'], + ) + + +class FloquetAnalysis: + def __init__( + self, + H0: qt.Qobj, + H1: qt.Qobj, + state_indices: list, + omega_d_linspace: np.ndarray, + amp_linspace: np.ndarray, + options: Options, + init_data_to_save: dict | None = None, + ): + self.H0 = H0 + self.H1 = H1 + self.state_indices = state_indices + self.omega_d_linspace = omega_d_linspace + self.amp_linspace = amp_linspace + self.options = options + self.init_data_to_save = init_data_to_save + # Save in _init_attrs for later re-initialization. Everything added to self + # after this is a derived quantity + self._init_attrs = set(self.__dict__.keys()) + self.num_states = H0.shape[0] + self.exponent_pair_idx_map = self._create_exponent_pair_idx_map() + self._num_fit_ranges = int(np.ceil(1 / options.fit_range_fraction)) + self._num_amp_pts_per_range = int( + np.floor(len(self.amp_linspace) / self._num_fit_ranges) + ) + array_shape = ( + len(self.omega_d_linspace), + len(self.amp_linspace), + len(self.state_indices), + ) + self.max_overlap_data = np.zeros(array_shape) + self.quasienergies = np.zeros(array_shape) + # fit over the full range, using states identified by the fit over + # intermediate ranges + self.coefficient_matrix = np.zeros( + (len(self.state_indices), self.num_states, len(self.exponent_pair_idx_map)), + dtype=complex, + ) + self.displaced_state_overlaps = np.zeros(array_shape) + self._displaced_state_overlaps = np.zeros(array_shape) + self.floquet_mode_idxs = np.zeros(array_shape, dtype=int) + self.floquet_mode_data = np.zeros( + (*array_shape, self.num_states), dtype=complex + ) + self.avg_excitation = np.zeros( + (len(self.omega_d_linspace), len(self.amp_linspace), self.num_states) + ) + self.all_quasienergies = np.zeros_like(self.avg_excitation) + + def _get_initdata(self) -> dict: + """Collect all attributes for writing to file.""" + init_dict = {k: v for k, v in self.__dict__.items() if k in self._init_attrs} + new_init_dict = {} + for k, v in init_dict.items(): + if isinstance(v, qt.Qobj): + new_init_dict[k] = v.data.toarray().tolist() + elif isinstance(v, np.ndarray): + new_init_dict[k] = v.tolist() + elif isinstance(v, Options): + new_init_dict['options'] = vars(v) + else: + new_init_dict[k] = v + return new_init_dict + + def param_dict(self) -> dict: + if self.init_data_to_save is None: + self.init_data_to_save = {} + param_dict = vars(self.options) | self.init_data_to_save + return param_dict | { + 'state_indices': self.state_indices, + 'omega_d_linspace': self.omega_d_linspace, + 'amp_linspace': self.amp_linspace, + 'num_states': self.num_states, + 'num_fit_ranges': self._num_fit_ranges, + 'num_amp_pts_per_range': self._num_amp_pts_per_range, + 'exp_pair_map': self.exponent_pair_idx_map, + 'floquet_analysis_init': self._get_initdata(), + } + + def __str__(self): + params = self.param_dict() + params.pop('floquet_analysis_init') + parts_str = '\n'.join(f'{k}: {v}' for k, v in params.items() if v is not None) + return f'Running floquet simulation with parameters: \n' + parts_str + + def run(self, filepath: str = 'tmp.h5py') -> None: + """Run the whole floquet simulation.""" + write_to_h5(filepath, {}, self.param_dict()) + print(self) + start_time = time.time() + self.floquet_main(filepath=filepath) + print(f'finished floquet main in {(time.time()-start_time)/60} minutes') + update_data_in_h5(filepath, self.assemble_data_dict()) + + def assemble_data_dict(self) -> dict: + """Collect all computed data in preparation for writing to file.""" + data_dict = { + 'max_overlap_data': self.max_overlap_data, + 'floquet_mode_idxs': self.floquet_mode_idxs, + 'fit_data': self.coefficient_matrix, + 'quasienergies': self.quasienergies, + 'displaced_state_overlaps': self.displaced_state_overlaps, + '_displaced_state_overlaps': self._displaced_state_overlaps, + 'all_quasienergies': self.all_quasienergies, + 'avg_excitation': self.avg_excitation, + } + if self.options.save_floquet_mode_data: + data_dict['floquet_mode_data'] = self.floquet_mode_data + return data_dict + + def hamiltonian(self, params: tuple[float, float]) -> list[qt.Qobj]: + """Return the Hamiltonian we actually simulate.""" + omega_d, amp = params + return [self.H0, [amp * self.H1, lambda t, _: np.cos(omega_d * t)]] + + def run_one_floquet( + self, params: tuple[float, float] + ) -> tuple[np.ndarray, qt.Qobj]: + """Run one instance of the problem for a pair of drive frequency and amp. + + Returns floquet modes as numpy column vectors, as well as the quasienergies. + Parameters. + ---------- + params: tuple + pair of drive frequency and amp + """ + omega_d, amp = params + T = 2.0 * np.pi / omega_d + f_modes_0, f_energies_0 = qt.floquet_modes( + self.hamiltonian(params), # type: ignore + T, + options=qt.Options(nsteps=self.options.nsteps), + ) + sampling_time = self.options.floquet_sampling_time_fraction * T % T + if sampling_time != 0: + f_modes_t = qt.floquet_modes_t( + f_modes_0, + f_energies_0, + sampling_time, + self.hamiltonian(params), # type: ignore + T, + options=qt.Options(nsteps=self.options.nsteps), + ) + else: + f_modes_t = f_modes_0 + return f_modes_t, f_energies_0 + + def _calculate_modes_quasies_ovlps( + self, + f_modes_energies: tuple[np.ndarray, qt.Qobj], + params_0: tuple[float, float], + disp_coeffs: np.ndarray, + ) -> np.ndarray: + """Return overlaps with "ideal" bare state at a given pair of (omega_d, amp). + + Parameters. + ---------- + f_modes_energies: tuple + output of self.run_one_floquet(params) + params_0: tuple + (omega_d_0, amp_0) to use for displaced fit + disp_coeffs: ndarray + matrix of coefficients for the displaced state + """ + f_modes_0, f_energies_0 = f_modes_energies + # construct column vectors to compute overlaps + f_modes_cols = np.array( + [f_modes_0[idx].data.toarray()[:, 0] for idx in range(self.num_states)], + dtype=complex, + ).T + # return overlap and floquet mode + modes_quasies_ovlps = np.zeros( + (len(self.state_indices), 3 + self.num_states), dtype=complex + ) + ideal_displaced_state_array = np.array( + [ + self.displaced_state(*params_0, disp_coeffs[array_idx], state_idx) + .dag() + .data.toarray()[0] + for array_idx, state_idx in enumerate(self.state_indices) + ] + ) + overlaps = np.einsum('ij,jk->ik', ideal_displaced_state_array, f_modes_cols) + # take the argmax along k + f_idxs = np.argmax(np.abs(overlaps), axis=1) + for array_idx, _state_idx in enumerate(self.state_indices): + f_idx = f_idxs[array_idx] + bare_state_overlap = overlaps[array_idx, f_idx] + modes_quasies_ovlps[array_idx, 0] = bare_state_overlap + modes_quasies_ovlps[array_idx, 1] = f_idx + modes_quasies_ovlps[array_idx, 2] = f_energies_0[f_idx] + modes_quasies_ovlps[array_idx, 3:] = ( + np.sign(bare_state_overlap) * f_modes_cols[:, f_idx] + ) + return modes_quasies_ovlps + + def bare_state_array(self) -> np.ndarray: + """Return array of bare states. + + Used to specify initial bare states for the Blais branch analysis. + """ + return np.squeeze( + np.array([qt.basis(self.num_states, idx) for idx in range(self.num_states)]) + ) + + def _step_in_amp( + self, f_modes_energies: tuple[np.ndarray, qt.Qobj], prev_f_modes: np.ndarray + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Perform Blais branch analysis. + + Gorgeous in its simplicity. Simply calculate overlaps of new floquet modes with + those from the previous amplitude step, and order the modes accordingly. So + ordered, compute the mean excitation number, yielding our branches. + """ + f_modes_0, f_energies_0 = f_modes_energies + f_modes_0 = np.squeeze(np.array(f_modes_0)) + all_overlaps = np.abs(np.einsum('ij,kj->ik', np.conj(prev_f_modes), f_modes_0)) + # assume that prev_f_modes_arr have been previously sorted. Question + # is which k index has max overlap? + max_idxs = np.argmax(all_overlaps, axis=1) + f_modes_ordered = f_modes_0[max_idxs] + avg_excitation = self._calculate_mean_excitation(f_modes_ordered) + sorted_quasi_es = f_energies_0[max_idxs] + return avg_excitation, sorted_quasi_es, f_modes_ordered + + def _calculate_mean_excitation(self, f_modes_ordered: np.ndarray) -> np.ndarray: + """Mean excitation number of ordered floquet modes. + + Based on Blais arXiv:2402.06615, specifically Eq. (12) but going without the + integral over floquet modes in one period. + """ + bare_states = self.bare_state_array() + overlaps_sq = np.abs(np.einsum('ij,kj->ik', bare_states, f_modes_ordered)) ** 2 + # sum over bare excitations weighted by excitation number + return np.einsum('ik,i->k', overlaps_sq, np.arange(0, self.num_states)) + + def _single_overlap_with_displaced_state( + self, omega_d_amp: tuple[float, float], disp_coeffs_for_new_amp: np.ndarray + ) -> np.ndarray: + """Calculate overlap of floquet mode with 'ideal' displaced state. + + This is done here for a given omega_d, amp pair. + """ + overlap = np.zeros(len(self.state_indices)) + omega_d, amp = omega_d_amp + for array_idx, state_idx in enumerate(self.state_indices): + floquet_data = self.floquet_mode_data[ + self.omega_d_to_idx(omega_d), self.amp_to_idx(amp, omega_d), array_idx + ] + disp_state = self.displaced_state( + omega_d, amp, disp_coeffs_for_new_amp[array_idx], state_idx=state_idx + ).dag() + overlap[array_idx] = np.abs(disp_state.data.toarray()[0] @ floquet_data) + return overlap + + def floquet_main(self, filepath: str = 'tmp.h5py') -> None: + """Perform floquet analysis over range of amplitudes and drive frequencies. + + This function largely performs two calculations. The first is the Xiao analysis + introduced in https://arxiv.org/abs/2304.13656, fitting the extracted Floquet + modes to the "ideal" displaced state which does not include resonances by design + (because we fit to a low order polynomial and ignore any floquet modes with + overlap with the bare state below a given threshold). This analysis produces the + "scar" plots. The second is the Blais branch analysis, which tracks the Floquet + modes by stepping in drive amplitude for a given drive frequency. For this + reason the code is structured to parallelize over drive frequency, but scans in + a loop over drive amplitude. This way the two calculations can be performed + simultaneously. + + A nice bonus is that both of the above mentioned calculations determine + essentially independently whether a resonance occurs. In the first, it is + deviation of the Floquet mode from the fitted displaced state. In the second, + it is branch swapping that indicates a resonance, independent of any fit. Thus + the two simulations can be used for cross validation of one another. + + We perform these simulations iteratively over the drive amplitudes as specified + by fit_range_fraction. This is to allow for simulations stretching to large + drive amplitudes, where the overlap with the bare eigenstate would fall below + the threshold (due to ac Stark shift) even in the absence of any resonances. + We thus use the fit from the previous range of drive amplitudes as our new bare + state. + """ + # for all omega_d, the bare states are identical at zero drive. We define + # two sets of bare modes (prev_f_modes_arr and disp_coeffs_for_prev_amp) + # because for the fit calculation, the bare modes are specified as fit + # coefficients, whereas for the Blais calculation, the bare modes are specified + # as actual kets. + prev_f_modes_arr = np.tile( + self.bare_state_array()[None, :, :], (len(self.omega_d_linspace), 1, 1) + ) + disp_coeffs_for_prev_amp = np.array( + [ + self.bare_state_coefficients(state_idx) + for state_idx in self.state_indices + ] + ) + for amp_range_idx in range(self._num_fit_ranges): + print(f'calculating for amp_range_idx={amp_range_idx}') + # edge case if range doesn't fit in neatly + if amp_range_idx == self._num_fit_ranges - 1: + amp_range_idx_final = len(self.amp_linspace) + else: + amp_range_idx_final = (amp_range_idx + 1) * self._num_amp_pts_per_range + amp_idxs = [ + amp_range_idx * self._num_amp_pts_per_range, + amp_range_idx_final, + ] + # now perform floquet mode calculation for amp_range_idx + # need to pass forward the floquet modes from the previous amp range + # which allow us to identify floquet modes that may have been displaced + # far from the origin + prev_f_modes_arr = self._floquet_main_for_amp_range( + amp_idxs, disp_coeffs_for_prev_amp, prev_f_modes_arr + ) + disp_coeffs_for_prev_amp = self._displaced_states_fit( + amp_idxs, disp_coeffs_for_prev_amp + ) + overlaps = self._overlap_with_displaced_states( + amp_idxs, disp_coeffs_for_prev_amp + ) + # save the extracted overlaps for later use when fitting over the whole + # shebang + self._displaced_state_overlaps[:, amp_idxs[0]: amp_idxs[1], :] = overlaps + update_data_in_h5(filepath, self.assemble_data_dict()) + # the previously extracted coefficients were valid for the amplitude ranges + # we asked for the fit over. Now armed with with correctly identified floquet + # modes, we recompute these coefficients over the whole sea of floquet mode data + # to get a plot that is free from numerical artifacts associated with + # the fits being slightly different at the boundary of ranges + full_displaced_fit = self._displaced_states_fit() + # try and be a little bit of a defensive programmer here, don't + # just e.g. overwrite self.coefficient_matrix + self.coefficient_matrix[:, :, :] = full_displaced_fit + true_overlaps = self._overlap_with_displaced_states( + [0, len(self.amp_linspace)], full_displaced_fit + ) + self.displaced_state_overlaps[:, :, :] = true_overlaps + + def _floquet_main_for_amp_range( + self, + amp_idxs: list, + disp_coeffs_for_prev_amp: np.ndarray, + prev_f_modes_arr: np.ndarray, + ) -> np.ndarray: + """Run the floquet simulation over a specific amplitude range.""" + amp_range_vals = self.amp_linspace[amp_idxs[0]: amp_idxs[1]] + + def _run_floquet_and_calculate( + omega_d: float, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + omega_d_idx = self.omega_d_to_idx(omega_d) + amps_for_omega_d = amp_range_vals[:, omega_d_idx] + avg_excitation_arr = np.zeros((len(amps_for_omega_d), self.num_states)) + quasienergies_arr = np.zeros_like(avg_excitation_arr) + modes_quasies_ovlps_arr = np.zeros( + (len(amps_for_omega_d), len(self.state_indices), 3 + self.num_states), + dtype=complex, + ) + prev_f_modes_for_omega_d = prev_f_modes_arr[omega_d_idx] + # for evaluating the displaced state (might be dangerous to evaluate + # outside of fitted window) + params_0 = (omega_d, amp_range_vals[0, omega_d_idx]) + for amp_idx, amp in enumerate(amps_for_omega_d): + params = (omega_d, amp) + f_modes_energies = self.run_one_floquet(params) + modes_quasies_ovlps = self._calculate_modes_quasies_ovlps( + f_modes_energies, params_0, disp_coeffs_for_prev_amp + ) + modes_quasies_ovlps_arr[amp_idx] = modes_quasies_ovlps + avg_excitation, quasi_es, new_f_modes_arr = self._step_in_amp( + f_modes_energies, prev_f_modes_for_omega_d + ) + avg_excitation_arr[amp_idx] = np.real(avg_excitation) + quasienergies_arr[amp_idx] = np.real(quasi_es) + prev_f_modes_for_omega_d = new_f_modes_arr + return ( + modes_quasies_ovlps_arr, + avg_excitation_arr, + quasienergies_arr, + prev_f_modes_for_omega_d, + ) + + floquet_data = list( + parallel_map( + self.options.num_cpus, _run_floquet_and_calculate, self.omega_d_linspace + ) + ) + ( + all_modes_quasies_ovlps, + all_avg_excitation, + all_quasienergies, + f_modes_last_amp, + ) = list(zip(*floquet_data, strict=False)) + floquet_mode_array = np.array(all_modes_quasies_ovlps, dtype=complex).reshape( + ( + len(self.omega_d_linspace), + len(amp_range_vals), + len(self.state_indices), + 3 + self.num_states, + ) + ) + f_modes_last_amp = np.array(f_modes_last_amp, dtype=complex).reshape( + (len(self.omega_d_linspace), self.num_states, self.num_states) + ) + self.max_overlap_data[:, amp_idxs[0]: amp_idxs[1]] = np.real( + floquet_mode_array[..., 0] + ) + self.floquet_mode_idxs[:, amp_idxs[0]: amp_idxs[1]] = np.real( + floquet_mode_array[..., 1] + ) + self.quasienergies[:, amp_idxs[0]: amp_idxs[1]] = np.real( + floquet_mode_array[..., 2] + ) + self.floquet_mode_data[:, amp_idxs[0]: amp_idxs[1]] = floquet_mode_array[ + ..., 3: + ] + self.avg_excitation[:, amp_idxs[0]: amp_idxs[1]] = np.array( + all_avg_excitation + ).reshape((len(self.omega_d_linspace), len(amp_range_vals), self.num_states)) + self.all_quasienergies[:, amp_idxs[0]: amp_idxs[1]] = np.array( + all_quasienergies + ).reshape((len(self.omega_d_linspace), len(amp_range_vals), self.num_states)) + return f_modes_last_amp + + def _omega_d_amp_params(self, amp_idxs: list | None = None) -> itertools.chain: + """Return ordered chain object of the specified omega_d and amplitude values.""" + if amp_idxs is None: + amp_range_vals = self.amp_linspace + else: + amp_range_vals = self.amp_linspace[amp_idxs[0]: amp_idxs[1]] + _omega_d_amp_params = [ + product([omega_d], amp_vals) + for omega_d, amp_vals in zip( + self.omega_d_linspace, amp_range_vals.T, strict=False + ) + ] + return chain(*_omega_d_amp_params) + + def _overlap_with_displaced_states( + self, amp_idxs: list, disp_coeffs_for_new_amp: np.ndarray + ) -> np.ndarray: + """Calculate overlap of floquet modes with 'ideal' displaced states. + + This is done here for a specific amplitude range. + """ + + def run_overlap_displaced(omega_d_amp: tuple[float, float]) -> np.ndarray: + return self._single_overlap_with_displaced_state( + omega_d_amp, disp_coeffs_for_new_amp + ) + + omega_d_amp_params = self._omega_d_amp_params(amp_idxs) + amp_range_vals = self.amp_linspace[amp_idxs[0]: amp_idxs[1]] + result = list( + parallel_map( + self.options.num_cpus, run_overlap_displaced, omega_d_amp_params + ) + ) + return np.array(result).reshape( + (len(self.omega_d_linspace), len(amp_range_vals), len(self.state_indices)) + ) + + def _displaced_states_fit( + self, amp_idxs: list | None = None, disp_coeffs_for_prev_amp: np.ndarray = None + ) -> np.ndarray: + """Loop over all states and perform the fit for a given amplitude range. + + If amp_idxs and disp_coeffs_for_prev_amp are None, then this indicates that we + are on the final pass and we should perform the fit over the whole range of + amplitudes. In this case we utilize the previously computed overlaps of the + floquet modes with the displaced states to obtain the mask with which we exclude + some data from the fit (because we suspect they've hit resonances). + """ + + def fit_for_state(array_state_idx: tuple[int, int]) -> np.ndarray: + return self._displaced_states_fit_for_amp_and_state( + amp_idxs, disp_coeffs_for_prev_amp, array_state_idx + ) + + array_idxs = np.arange(len(self.state_indices)) + array_state_idxs = zip(array_idxs, self.state_indices, strict=False) + fit_data = list( + parallel_map(self.options.num_cpus, fit_for_state, array_state_idxs) + ) + return np.array(fit_data, dtype=complex).reshape( + (len(self.state_indices), self.num_states, len(self.exponent_pair_idx_map)) + ) + + def _ravel_and_filter_params( + self, mask: np.ndarray, omega_d_amp_params: list + ) -> list: + mask = mask.ravel() + return [ + omega_d_amp_params[i] + for i in range(len(mask)) + if np.abs(mask[i]) > self.options.overlap_cutoff + ] + + def _ravel_and_filter_mode_data( + self, mask: np.ndarray, floquet_data: np.ndarray, state_idx_component: int + ) -> np.ndarray: + mask = mask.ravel() + floquet_idx_data_bare_component = floquet_data[ + :, :, state_idx_component + ].ravel() + # only fit states that we think haven't run into + # a nonlinear transition (same for omega_d_amp_filtered above) + return floquet_idx_data_bare_component[ + np.abs(mask) > self.options.overlap_cutoff + ] + + def _displaced_states_fit_for_amp_and_state( + self, + amp_idxs: list | None = None, + disp_coeffs_for_prev_amp: np.ndarray | None = None, + array_state_idx: tuple[int, int] = (0, 0), + ) -> np.ndarray: + """Find the displaced states. + + This is done here for all frequencies, a given amplitude window and a given + state. + """ + # For the bare state, we look at the amplitude at the leading edge of the + # window (that is, the smallest amplitude). This is the most natural choice, + # as it is most analogous to what is done in the first window when the overlap + # is computed against bare eigenstates (that obviously don't have amplitude + # dependence). Moreover, the fit coefficients for the previous window by + # definition were obtained in a window that does not include the one we are + # currently investigating. Asking for the state for amplitude values outside of + # the fit window should be done at your own peril. + array_idx, state_idx = array_state_idx + if amp_idxs is None: + floquet_idx_data = self.floquet_mode_data[:, :, array_idx, :] + else: + floquet_idx_data = self.floquet_mode_data[ + :, amp_idxs[0]: amp_idxs[1], array_idx, : + ] + if disp_coeffs_for_prev_amp is None: + # this means we are on the final lap, and want to compare with previously + # computed overlaps + ovlp_with_bare_state = self._displaced_state_overlaps[:, :, array_idx] + else: + + def _compute_bare_state(omega_d: float) -> np.ndarray: + omega_d_idx = self.omega_d_to_idx(omega_d) + return self.displaced_state( + omega_d, + self.amp_linspace[amp_idxs[0], omega_d_idx], + disp_coeffs_for_prev_amp[array_idx], + state_idx, + ).full()[:, 0] + + bare_states = np.array( + [_compute_bare_state(omega_d) for omega_d in self.omega_d_linspace], + dtype=complex, + ) + # bare states may differ as a function of omega_d, hence the bare states + # have an index of i that we don't sum over + # indices are i: omega_d, j: amp, k: components of state + # this serves as the mask that is passed to _disp_coeffs_fit + ovlp_with_bare_state = np.abs( + np.einsum('ijk,ik->ij', floquet_idx_data, np.conj(bare_states)) + ) + omega_d_amp_data_slice = list(self._omega_d_amp_params(amp_idxs)) + omega_d_amp_filtered = self._ravel_and_filter_params( + ovlp_with_bare_state, omega_d_amp_data_slice + ) + if len(omega_d_amp_filtered) < len(self.exponent_pair_idx_map): + warnings.warn( + 'Not enough data points to fit. Returning previous fit', stacklevel=3 + ) + return disp_coeffs_for_prev_amp + num_coeffs = len(self.exponent_pair_idx_map) + coefficient_matrix_for_amp_and_state = np.zeros( + (self.num_states, num_coeffs), dtype=complex + ) + for state_idx_component in range(self.num_states): + floquet_component_filtered = self._ravel_and_filter_mode_data( + ovlp_with_bare_state, floquet_idx_data, state_idx_component + ) + bare_same = state_idx_component == state_idx + bare_component_fit = self._disp_coeffs_fit( + omega_d_amp_filtered, floquet_component_filtered, bare_same + ) + coefficient_matrix_for_amp_and_state[state_idx_component, :] = ( + bare_component_fit + ) + return coefficient_matrix_for_amp_and_state + + def _disp_coeffs_fit( + self, + omega_d_amp_filtered: list, + floquet_component_filtered: np.ndarray, + bare_same: bool, + ) -> np.ndarray: + """Fit the floquet modes to an "ideal" displaced state based on a polynomial. + + This is done here over the grid specified by omega_d_amp_data_slice. We ignore + floquet data indicated by mask, where we suspect by looking at overlaps with the + bare state that we have hit a resonance. + """ + p0 = np.zeros(len(self.exponent_pair_idx_map)) + # fit the real and imaginary parts of the overlap separately + popt_r = self._fit_params( + omega_d_amp_filtered, np.real(floquet_component_filtered), p0, bare_same + ) + popt_i = self._fit_params( + omega_d_amp_filtered, + np.imag(floquet_component_filtered), + p0, + False, # for the imaginary part, constant term should always be zero + ) + return popt_r + 1j * popt_i + + def bare_state_coefficients(self, state_idx: int) -> np.ndarray: + """For bare state only component is itself.""" + coefficient_matrix_for_amp_and_state = np.zeros( + (self.num_states, len(self.exponent_pair_idx_map)), dtype=complex + ) + coefficient_matrix_for_amp_and_state[state_idx, 0] = 1.0 + return coefficient_matrix_for_amp_and_state + + def floquet_mode(self, omega_d: float, amp: float, array_idx: int) -> np.ndarray: + """Helper function for extracting the floquet mode.""" + omega_d_idx = self.omega_d_to_idx(omega_d) + amp_idx = self.amp_to_idx(amp, omega_d) + return self.floquet_mode_data[omega_d_idx, amp_idx, array_idx] + + def _fit_params( + self, XYdata: list, Zdata: np.ndarray, p0: tuple | np.ndarray, bare_same: bool + ) -> np.ndarray: + poly_fit = functools.partial(self._polynomial_fit, bare_same=bare_same) + try: + popt, pcov = sp.optimize.curve_fit(poly_fit, XYdata, Zdata, p0=p0) + except RuntimeError: + warnings.warn( + 'fit failed for a bare component, returning zeros for the fit', + stacklevel=3, + ) + popt = np.zeros(len(p0)) + return popt + + def _create_exponent_pair_idx_map(self) -> dict: + """Create dictionary of terms in polynomial that we fit. + + We truncate the fit if e.g. there is only a single frequency value to scan over + but the fit is nominally set to order four. We additionally eliminate the + constant term that should always be either zero or one. + """ + cutoff_omega_d = min(len(self.omega_d_linspace), self.options.fit_cutoff) + cutoff_amp = min(len(self.amp_linspace), self.options.fit_cutoff) + idx_exp_map = [ + (idx_1, idx_2) + for idx_1 in range(cutoff_omega_d) + for idx_2 in range(cutoff_amp) + ] + # kill constant term, which should always be 1 or 0 depending + # on if the component is the same as the state being fit + idx_exp_map.remove((0, 0)) + weighted_vals = [1.01 * idx_1 + idx_2 for (idx_1, idx_2) in idx_exp_map] + sorted_idxs = np.argsort(weighted_vals) + sorted_idx_exp_map = {} + counter = 0 + for sorted_idx in sorted_idxs: + exponents = idx_exp_map[sorted_idx] + if sum(exponents) <= self.options.fit_cutoff: + sorted_idx_exp_map[counter] = exponents + counter += 1 + return sorted_idx_exp_map + + def _polynomial_fit( + self, + xydata: np.ndarray, + *fit_params: np.ndarray | tuple, + bare_same: bool = False, + ) -> np.ndarray | float: + """Fit function to pass to curve fit, assume a 2D polynomial.""" + exp_pair_map = self.exponent_pair_idx_map + omega_d, amp = xydata.T + result = 1.0 if bare_same else 0.0 + for idx in exp_pair_map: + exp_pair = exp_pair_map[idx] + result += fit_params[idx] * omega_d ** exp_pair[0] * amp ** exp_pair[1] + return result + + def displaced_state( + self, omega_d: float, amp: float, coeffs: np.ndarray, state_idx: int + ) -> qt.Qobj: + """Construct the ideal displaced state, not including nonlinear transitions.""" + return sum( + self._polynomial_fit( + np.array([omega_d, amp]), + *coeffs[state_idx_component, :], + bare_same=state_idx == state_idx_component, + ) + * qt.basis(self.num_states, state_idx_component) + for state_idx_component in range(self.num_states) + ).unit() + + def amp_to_range_idx(self, amp: float, omega_d: float) -> int: + amp_idx = self.amp_to_idx(amp, omega_d) + return int(np.floor(amp_idx / self._num_amp_pts_per_range)) + + def omega_d_to_idx(self, omega_d: float) -> np.ndarray[int]: + return np.argmin(np.abs(self.omega_d_linspace - omega_d)) + + def amp_to_idx(self, amp: float, omega_d: float) -> np.ndarray[int]: + omega_d_idx = self.omega_d_to_idx(omega_d) + return np.argmin(np.abs(self.amp_linspace[:, omega_d_idx] - amp)) diff --git a/floquet/options.py b/floquet/options.py new file mode 100644 index 0000000..323222a --- /dev/null +++ b/floquet/options.py @@ -0,0 +1,54 @@ +from __future__ import annotations + + +class Options: + """Options for the floquet analysis. + + fit_range_fraction: float + Fraction of the amplitude range to sweep over before changing the definition of + the bare state to that of the fitted state from the previous range. For instance + if fit_range_fraction=0.4, then the amplitude range is split up into three + chunks: the first 40% of the amplitude linspace, then from 40% -> 80%, then + from 80% to the full range. For the first fraction of amplitudes, they are + compared to the bare eigenstates for identification. For the second range, they + are compared to the fitted state from the first range. And so on. Defaults to + 1.0, indicating that no iteration is performed. + floquet_sampling_time_fraction: float + What point of the drive period we want to sample the Floquet modes. Defaults to + 0.0, indicating the floquet modes at t=0*T where T is the drive period. + fit_cutoff: int + Cutoff for the fit polynomial of the displaced state. Defaults to 4. + overlap_cutoff: float + Cutoff for fitting overlaps. Floquet modes with overlap with the "bare" state + below this cutoff are not included in the fit (as they may be experiencing a + resonance). Defaults to 0.8. + nsteps: int + QuTiP integration parameter, number of steps the solver can take. Defaults + to 30_000. + num_cpus: int + number of cpus to use in parallel computation of Floquet modes over the + different values of omega_d, amp. Defaults to 1. + save_floquet_mode_data: bool + Indicating whether to save the extracted Floquet modes themselves. Such data is + often unnecesary and requires a fair amount of storage, so the default is False. + """ + + def __init__( + self, + fit_range_fraction: float = 1.0, + floquet_sampling_time_fraction: float = 0.0, + fit_cutoff: int = 4, + overlap_cutoff: float = 0.8, + nsteps: int = 30_000, + num_cpus: int = 1, + init_data_to_save: dict | None = None, + save_floquet_mode_data: bool = False, + ): + self.fit_range_fraction = fit_range_fraction + self.floquet_sampling_time_fraction = floquet_sampling_time_fraction + self.fit_cutoff = fit_cutoff + self.overlap_cutoff = overlap_cutoff + self.nsteps = nsteps + self.num_cpus = num_cpus + self.init_data_to_save = init_data_to_save + self.save_floquet_mode_data = save_floquet_mode_data diff --git a/floquet/utils/__init__.py b/floquet/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/floquet/utils/file_io.py b/floquet/utils/file_io.py new file mode 100644 index 0000000..a083d9a --- /dev/null +++ b/floquet/utils/file_io.py @@ -0,0 +1,53 @@ +import glob +import os +import re + +import h5py + + +def generate_file_path(extension: str, file_name: str, path: str) -> str: + # Ensure the path exists. + os.makedirs(path, exist_ok=True) # noqa: PTH103 + + # Create a save file name based on the one given; ensure it will + # not conflict with others in the directory. + max_numeric_prefix = -1 + for file_name_ in glob.glob(os.path.join(path, '*')): # noqa: PTH118, PTH207 + if f'_{file_name}.{extension}' in file_name_: + numeric_prefix = int( + re.match(r'(\d+)_', os.path.basename(file_name_)).group(1) # noqa: PTH119 + ) + max_numeric_prefix = max(numeric_prefix, max_numeric_prefix) + + # Generate the file path. + return os.path.join( # noqa: PTH118 + path, f'{str(max_numeric_prefix + 1).zfill(5)}_{file_name}.{extension}' + ) + + +def extract_info_from_h5(filepath: str) -> [dict, dict]: + data_dict = {} + with h5py.File(filepath, 'r') as f: + for key in f: + data_dict[key] = f[key][()] + param_dict = dict(f.attrs.items()) + return data_dict, param_dict + + +def update_data_in_h5(filepath: str, data_dict: dict) -> None: + with h5py.File(filepath, 'a') as f: + for key, val in data_dict.items(): + if key in f: + del f[key] + f.create_dataset(key, data=val, track_order=True) + + +def write_to_h5(filepath: str, data_dict: dict, param_dict: dict) -> None: + with h5py.File(filepath, 'a') as f: + for key, val in data_dict.items(): + f.create_dataset(key, data=[val], chunks=True, maxshape=(None, *val.shape)) + for kwarg in param_dict: + try: + f.attrs[kwarg] = param_dict[kwarg] + except TypeError: + f.attrs[kwarg] = str(param_dict[kwarg]) diff --git a/floquet/utils/parallel.py b/floquet/utils/parallel.py new file mode 100644 index 0000000..d1ab817 --- /dev/null +++ b/floquet/utils/parallel.py @@ -0,0 +1,11 @@ +from collections.abc import Iterable + +import pathos + + +def parallel_map(num_cpus: int, func: callable, parameters: Iterable) -> map: + if num_cpus == 1: + return map(func, parameters) + + with pathos.pools.ProcessPool(nodes=num_cpus) as pool: + return pool.map(func, parameters) diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..18ba3aa --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,72 @@ +site_name: floquet +site_description: Documentation for the floquet software library +site_author: Daniel Weiss +site_url: http://dkweiss.net/floquet/ + +repo_url: https://github.com/dkweiss31/floquet +repo_name: dkweiss31/floquet +edit_uri: "" + +theme: + name: material + features: + - navigation.sections # Sections are included in the navigation on the left. + - toc.integrate # Table of contents is integrated on the left; does not appear separately on the right. + palette: + - scheme: default + primary: indigo + accent: amber + toggle: + icon: material/weather-night + name: Switch to dark mode + - scheme: slate + primary: indigo + accent: amber + toggle: + icon: material/weather-sunny + name: Switch to light mode + icon: + repo: fontawesome/brands/github # GitHub logo in top right + +extra_javascript: + # To make MathJax work, see https://squidfunk.github.io/mkdocs-material/reference/mathjax/ + - https://polyfill.io/v3/polyfill.min.js?features=es6 + - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + +extra_css: + - _static/custom_css.css + +markdown_extensions: + - pymdownx.arithmatex: # Render LaTeX via MathJax + generic: true + - pymdownx.superfences # Seems to enable syntax highlighting when used with the Material theme. + - pymdownx.details # Allowing hidden expandable regions denoted by ??? + - pymdownx.snippets: # Include one Markdown file into another + base_path: docs + - admonition + - toc: + permalink: "ยค" # Adds a clickable permalink to each section heading + toc_depth: 4 + +plugins: + - search # default search plugin; needs manually re-enabling when using any other plugins + - autorefs # Cross-links to headings + - mknotebooks # Jupyter notebooks + - mkdocstrings: + handlers: + python: + rendering: + show_source: false + show_if_no_docstring: true + show_signature_annotations: true + members_order: source + heading_level: 4 + selection: + inherited_members: true # Allow looking up inherited methods + +nav: + - 'index.md' + - API: + - 'floquet.md' + - Examples: + - Transmon: 'examples/transmon.ipynb' diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..f415066 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,132 @@ +[project] +name = "floquet" +version = "0.1.0" +requires-python = ">=3.10" +description = "Floquet simulations for identifying resonances in quantum systems." +dependencies = [ + "numpy", + "scipy", + "scqubits", + "qutip>=4.7.6,<5.0", + "matplotlib", + "jupyter", + "h5py", + "ipython", + "pathos", +] + +[build-system] +requires = ["flit_core >=3.2,<4"] +build-backend = "flit_core.buildapi" + +[project.optional-dependencies] +dev = [ + "taskipy", + "ruff", + "codespell", + "pygments", + "pytest>=7.0,<8.0", + "pytest-xdist", + "pymdown-extensions", + "mkdocs", + "mkdocs-material", + "mkdocstrings[python]", + "mkdocs-gen-files", + "mkdocs-literate-nav", + "mkdocs-section-index", + "mkdocs-simple-hooks", + "mkdocs-glightbox", + "mkdocs-exclude", + "mknotebooks", + "nbconvert==6.5.0", + "sybil>=6", + "black", # needed by mkdocstrings to format function signatures +] + +[tool.ruff] +extend-include = ["*.ipynb"] +exclude = ["examples/*.py", "examples/*.ipynb"] + +[tool.ruff.format] +quote-style = "single" +docstring-code-format = true +skip-magic-trailing-comma = true + +[tool.ruff.lint] +select = [ + "F", "E", "W", "C90", "I", "D", "UP", "YTT", "ANN", "BLE", "B", "A", "C4", "FA", + "INP", "NPY201", "PIE", "T20", "PYI", "PT", "RSE", "RET", "SLF", "SIM", "INT", + "ARG", "PTH", "PL", "TRY", "FLY", "NPY", "RUF", +] +extend-select = ["D204", "D400", "D404", "D406", "D410"] +ignore = [ + "E741", + "C901", + "D100", "D101", "D102", "D103", "D104", "D105", "D106", "D107", "D417", + "ANN101", "ANN003", + "TRY003", + "PLR2004", "PLR0912", "PLR0913", + "T201" +] + +[tool.ruff.lint.isort] +split-on-trailing-comma = false + +[tool.ruff.lint.flake8-annotations] +suppress-none-returning = true + +[tool.ruff.lint.flake8-comprehensions] +allow-dict-calls-with-keyword-arguments = true + +[tool.ruff.lint.per-file-ignores] +"__init__.py" = ["F401", "PLC0414"] +"tests/**.py" = ["ANN"] +"examples/**.py" = ["INP001"] + +[tool.ruff.lint.pydocstyle] +convention = "google" + +[tool.codespell] +skip = ".git,*.ipynb" + +# === taskipy tasks definition === + +[tool.taskipy.tasks.lint] +cmd = 'echo "\n>>> ruff check --fix" && ruff check --fix' +help = "lint the code (ruff)" + +[tool.taskipy.tasks.format] +cmd = 'echo "\n>>> ruff format" && ruff format' +help = "auto-format the code (ruff)" + +[tool.taskipy.tasks.codespell] +cmd = 'echo "\n>>> codespell" && codespell tests floquet' +help = "check for misspellings (codespell)" + +[tool.taskipy.tasks.clean] +cmd = 'task lint && task format && task codespell' +help = "clean the code (ruff + codespell)" + +[tool.taskipy.tasks.test] +cmd = 'echo "\n>>> pytest -n=auto tests" && pytest -n=auto tests' +help = "run the unit tests suite (pytest)" + +[tool.taskipy.tasks.docbuild] +cmd = 'mkdocs build' +help = "build the documentation website" + +[tool.taskipy.tasks.docserve] +cmd = 'mkdocs serve' +help = "preview documentation website with hot-reloading" + +[tool.taskipy.tasks.all] +cmd = 'task clean && task test' +help = "run all tasks before a commit (ruff + codespell + pytest + doctest)" + +[tool.taskipy.tasks.ci] +cmd = '''echo "\n>>> ruff check" && ruff check && + echo "\n>>> ruff format --check" && ruff format --check && + task codespell && + task test && + task docbuild''' +help = "run all the CI checks" diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_floquet.py b/tests/test_floquet.py new file mode 100644 index 0000000..7ed34c9 --- /dev/null +++ b/tests/test_floquet.py @@ -0,0 +1,116 @@ +import itertools + +import numpy as np +import pytest +import qutip as qt +import scqubits as scq + +from floquet import ( + ChiacToAmp, + Options, + XiSqToAmp, + floquet_analysis, + floquet_analysis_from_file, +) + + +def _filepath(path): + d = path / 'sub' + d.mkdir() + return d / 'tmp.h5py' + + +@pytest.fixture(scope='session', autouse=True) +def setup_floquet(): + num_states = 20 + EJ = 29.0 + EC = 0.155 + ng = 0.0 + ncut = 21 + INIT_DATA_TO_SAVE = {'EJ': EJ, 'EC': EC, 'ng': ng, 'ncut': ncut} + tmon = scq.Transmon(EJ=29.0, EC=0.155, ng=0.0, ncut=21, truncated_dim=num_states) + omega_d_linspace = 2.0 * np.pi * np.linspace(6.9, 13, 39) + chi_ac_linspace = 2.0 * np.pi * np.linspace(0.0, 0.2, 40) + state_indices = [0, 1, 3] + + hilbert_space = scq.HilbertSpace([tmon]) + hilbert_space.generate_lookup() + evals = hilbert_space['evals'][0][0:num_states] + H0 = 2.0 * np.pi * qt.Qobj(np.diag(evals - evals[0])) + H1 = hilbert_space.op_in_dressed_eigenbasis(tmon.n_operator) + + options = Options(fit_range_fraction=0.5, num_cpus=6) + chi_to_amp = ChiacToAmp(H0, H1, state_indices, omega_d_linspace) + amp_linspace = chi_to_amp.amplitudes_for_omega_d(chi_ac_linspace) + floquet_transmon = floquet_analysis( + H0, + H1, + state_indices, + omega_d_linspace=omega_d_linspace, + amp_linspace=amp_linspace, + options=options, + init_data_to_save=INIT_DATA_TO_SAVE, + ) + return floquet_transmon, chi_to_amp, chi_ac_linspace + + +def test_chi_vs_xi(setup_floquet): + floquet_transmon, _, chi_ac_linspace = setup_floquet + amps_from_chi_ac = floquet_transmon.amp_linspace + EC = floquet_transmon.init_data_to_save['EC'] + xi_sq_linspace = 2.0 * chi_ac_linspace / EC / 2 / np.pi + xi_sq_to_amp = XiSqToAmp( + floquet_transmon.H0, + floquet_transmon.H1, + floquet_transmon.state_indices, + floquet_transmon.omega_d_linspace, + ) + amps_from_xi_sq = xi_sq_to_amp.amplitudes_for_omega_d(xi_sq_linspace) + rel_diff = np.abs( + (amps_from_xi_sq[1:] - amps_from_chi_ac[1:]) / amps_from_xi_sq[1:] + ) + assert np.max(rel_diff < 0.05) + + +def test_displaced_fit_and_reinit(setup_floquet, tmp_path): + floquet_transmon, chi_to_amp, _ = setup_floquet + filepath = _filepath(tmp_path) + floquet_transmon.run(filepath=filepath) + # for these random pairs, overlap should be near unity (most + # pairs don't correspond to a resonance!) + omega_d_vals = 2.0 * np.pi * np.array([7.5, 9.0, 11.8, 12.7]) + chi_ac_vals = np.array([0.01, 0.03, 0.05, 0.08, 0.1, 0.15, 0.19]) + omega_d_chi_ac = itertools.product(omega_d_vals, chi_ac_vals) + for omega_d, chi_ac in omega_d_chi_ac: + omega_d_idx = floquet_transmon.omega_d_to_idx(omega_d) + amp = chi_to_amp.amplitudes_for_omega_d(chi_ac)[0, omega_d_idx] + for arr_idx, state_idx in enumerate(floquet_transmon.state_indices): + disp_coeffs = floquet_transmon.coefficient_matrix + disp_gs = floquet_transmon.displaced_state( + omega_d, amp, disp_coeffs[arr_idx], state_idx=state_idx + ) + f_modes_energies = floquet_transmon.run_one_floquet((omega_d, amp)) + floquet_mode = floquet_transmon._calculate_modes_quasies_ovlps( + f_modes_energies, (omega_d, amp), disp_coeffs + ) + overlap = np.abs( + (qt.Qobj(floquet_mode[arr_idx, 3:]).dag() * disp_gs).data.toarray()[ + 0, 0 + ] + ) + assert 0.98 < overlap < 1.0 + # reinit + reinit_floquet_transmon = floquet_analysis_from_file(filepath) + assert reinit_floquet_transmon._get_initdata() == floquet_transmon._get_initdata() + + +def test_displaced_bare_state(setup_floquet): + floquet_transmon, chi_to_amp, _ = setup_floquet + for state_idx in floquet_transmon.state_indices: + ideal_state = qt.basis(floquet_transmon.num_states, state_idx) + disp_coeffs_bare = floquet_transmon.bare_state_coefficients(state_idx) + # omega_d and amp values shouldn't matter + calc_state = floquet_transmon.displaced_state( + 0.0, 0.0, disp_coeffs_bare, state_idx=state_idx + ) + assert calc_state == ideal_state