diff --git a/notebooks/be/BE_Processing_sidpy.ipynb b/notebooks/be/BE_Processing_sidpy.ipynb index b48518a6..0876b0d1 100644 --- a/notebooks/be/BE_Processing_sidpy.ipynb +++ b/notebooks/be/BE_Processing_sidpy.ipynb @@ -6,7 +6,9 @@ "metadata": {}, "outputs": [], "source": [ - "%matplotlib widget\n" + "# !pip install dask-ml\n", + "# !pip install git+https://github.com/pycroscopy/sidpy.git@main\n", + "# !pip install \"dask==2025.1.0\" \"distributed==2025.1.0\"" ] }, { @@ -15,38 +17,15 @@ "metadata": {}, "outputs": [], "source": [ - "import sys\n", - "sys.path.insert(0, r'/Users/rvv/Github/BGlib/')\n", - "sys.path.insert(0, r'/Users/rvv/Github/sidpy/')\n", - "\n", - "import BGlib.be as belib\n", - "import numpy as np\n", "import os\n", + "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import h5py\n", "import sidpy\n", "import SciFiReaders as sr\n", - "from BGlib.be.analysis.utils.sidpy_sho_fitter import SHOestimateGuess, SHOestimateGuess, SHO_fit_flattened\n", - "\n", - "\n", - "def load_data(file_path):\n", - " \"\"\"\n", - " \n", - " Given the path to a .h5 dataset, reads, patches and returns the BEPS dataset, \n", - " the frequency vector, and the index of the independent dimensions\n", - " \n", - " \"\"\"\n", - " \n", - " patcher = belib.translators.LabViewH5Patcher()\n", - " patcher.translate(file_path)\n", - " reader = sr.Usid_reader(file_path)\n", - " beps_raw = reader.read() \n", - " freq_axis = beps_raw.labels.index('Frequency (Hz)')\n", - " freq_vec = beps_raw._axes[freq_axis].values\n", - " all_dims = np.arange(len(data.shape))\n", - " ind_dims = np.delete(all_dims, freq_axis)\n", - " \n", - " return beps_raw, freq_vec, ind_dims\n" + "import BGlib.be as belib\n", + "from BGlib.be.analysis.utils.sidpy_sho_fitter import SHOestimateGuess, SHO_fit_flattened\n", + "from sidpy.proc.fitter_refactor import SidpyFitterRefactor" ] }, { @@ -55,65 +34,37 @@ "metadata": {}, "outputs": [], "source": [ - "folder_path = r'/Users/rvv/Dropbox (ORNL)/PTO LSMO JCY/031822'\n", - "file_name = r'60x60_0deg_0004.h5'\n", - "path_to_file = os.path.join(folder_path, file_name)" + "%matplotlib widget" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, - "outputs": [], - "source": [ - "#path_gwy = r'/Users/rvv/Dropbox (ORNL)/AE Related Stuff/RL for Walls/Paper/Resubmission/12202023/PTO_110_Virgin0001.gwy'\n", - "#gwy_reader = sr.GwyddionReader(path_gwy)\n", - "#data = gwy_reader.read()" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "#fig = data[0].plot()\n", - "\n", - "#threadpoolctl version. We need to insist on version 3.2.0 \n" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "/Users/rvv/Github/sidpy/sidpy/sid/translator.py:42: FutureWarning: Consider using sidpy.Reader instead of sidpy.Translator if possible and contribute your reader to ScopeReaders\n", - " warn('Consider using sidpy.Reader instead of sidpy.Translator if '\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "File is already Pycroscopy ready.\n" + "C:\\jawad_pc\\2_code\\20250912_code_pyc_bglib\\env_pb2\\lib\\site-packages\\sidpy\\sid\\translator.py:42: FutureWarning: Consider using sidpy.Reader instead of sidpy.Translator if possible and contribute your reader to ScopeReaders\n", + " warn('Consider using sidpy.Reader instead of sidpy.Translator if '\n", + "2026-06-12 10:50:07,949 - BGlib.be.translators.labview_h5_patcher - INFO - File is already Pycroscopy ready.\n" ] } ], "source": [ + "folder_path = r'../../../inputs/'\n", + "file_name = r'60x60_0deg_0004.h5'\n", + "path_to_file = os.path.join(folder_path, file_name)\n", "patcher = belib.translators.LabViewH5Patcher()\n", "patcher.translate(path_to_file)\n", - "\n", "reader = sr.Usid_reader(path_to_file)\n", "data = reader.read()" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -261,7 +212,7 @@ "Cycle: Cycle (generic) of size (2,)" ] }, - "execution_count": 7, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } @@ -272,7 +223,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -285,166 +236,165 @@ ], "source": [ "beps_raw = data\n", - "\n", "freq_axis = beps_raw.labels.index('Frequency (Hz)')\n", "freq_vec = beps_raw._axes[freq_axis].values\n", "all_dims = np.arange(len(data.shape))\n", "ind_dims = np.delete(all_dims, freq_axis)\n", - "print(ind_dims, freq_vec.shape)\n" + "print(ind_dims, freq_vec.shape)" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Setting Ind_dims from given values from user\n" + "Setup Complete. Params: 4 | Spatial Dims: [0, 1, 3, 4, 5]\n" ] } ], "source": [ "lb = [1E-6, freq_vec.min(), 50, -2*np.pi]\n", "ub = [1E-3, freq_vec.max(), 500, 2*np.pi]\n", - "beps_small = beps_raw[:5, :5, :]\n", - "\n", - "fitter = sidpy.proc.fitter.SidFitter(beps_small, SHO_fit_flattened,num_workers=16,\n", - " guess_fn = SHOestimateGuess,ind_dims=ind_dims,\n", - " threads=1, return_cov=False, return_fit=True, return_std=False,\n", - " km_guess=True,num_fit_parms = 4, n_clus = 5)\n" + "beps_small = beps_raw[:2, :2, :, :, :, :]\n", + "fitter = SidpyFitterRefactor(beps_small, SHO_fit_flattened, SHOestimateGuess, ind_dims=(freq_axis,), num_params=4, lower_bounds=lb, upper_bounds=ub)\n", + "fitter.setup_calc()" ] }, { "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "fitter.do_guess()\n", - "#fit_dset = fitter.do_fit(bounds = (lb,ub))" - ] - }, - { - "cell_type": "code", - "execution_count": 11, + "execution_count": 8, "metadata": {}, "outputs": [ { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "Warning: complex dataset detected. For Kmeans priors, we will treat real part only\n" + "C:\\jawad_pc\\2_code\\20250912_code_pyc_bglib\\env_pb2\\lib\\site-packages\\sidpy\\sid\\dataset.py:1517: UserWarning: Dimensional information will be lost. Please use fold, unfold to combine dimensions\n", + " warnings.warn('Dimensional information will be lost.\\\n", + "2026-06-12 10:50:24,015 - root - INFO - Starting _check_array\n", + "2026-06-12 10:50:24,064 - root - INFO - Finished _check_array in 0:00:00.050186\n", + "2026-06-12 10:50:24,069 - root - INFO - Starting init_scalable\n", + "2026-06-12 10:50:24,071 - dask_ml.cluster.k_means - INFO - Initializing with k-means||\n", + "2026-06-12 10:50:24,141 - dask_ml.cluster.k_means - INFO - Starting init iteration 1/ 9 , 1 centers\n" ] }, { - "name": "stderr", + "name": "stdout", "output_type": "stream", "text": [ - "/Users/rvv/pycroscopy_env/lib/python3.12/site-packages/distributed/client.py:3371: UserWarning: Sending large graph of size 864.85 MiB.\n", - "This may cause some slowdown.\n", - "Consider loading the data with Dask directly\n", - " or using futures or delayed objects to embed the data into the graph without repetition.\n", - "See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.\n", - " warnings.warn(\n", - "/Users/rvv/pycroscopy_env/lib/python3.12/site-packages/distributed/client.py:3371: UserWarning: Sending large graph of size 864.85 MiB.\n", - "This may cause some slowdown.\n", - "Consider loading the data with Dask directly\n", - " or using futures or delayed objects to embed the data into the graph without repetition.\n", - "See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.\n", - " warnings.warn(\n" + "Starting Dask K-Means Guess with 4 clusters...\n" ] }, { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "---Finished KMeans, onto fiting each KM Center---\n", - "Fitting center 0\n", - "Fitting center 1\n", - "Fitting center 2\n", - "Fitting center 3\n", - "Fitting center 4\n", - "Shapes of output of fitting function is 246 and original data is 123 Reshaping output dataset. You are responsible for reshaping\n", - "using generic parameters for dimension 0\n", - "using generic parameters for dimension 1\n" + "2026-06-12 10:50:24,213 - dask_ml.cluster.k_means - INFO - Finished init iteration 1/ 9 , 1 centers in 0:00:00.071166\n", + "2026-06-12 10:50:24,238 - dask_ml.cluster.k_means - INFO - Starting init iteration 2/ 9 , 1 centers\n", + "2026-06-12 10:50:24,310 - dask_ml.cluster.k_means - INFO - Finished init iteration 2/ 9 , 1 centers in 0:00:00.072084\n", + "2026-06-12 10:50:24,337 - dask_ml.cluster.k_means - INFO - Starting init iteration 3/ 9 , 5 centers\n", + "2026-06-12 10:50:24,406 - dask_ml.cluster.k_means - INFO - Finished init iteration 3/ 9 , 5 centers in 0:00:00.069253\n", + "2026-06-12 10:50:24,433 - dask_ml.cluster.k_means - INFO - Starting init iteration 4/ 9 , 7 centers\n", + "2026-06-12 10:50:24,512 - dask_ml.cluster.k_means - INFO - Finished init iteration 4/ 9 , 7 centers in 0:00:00.077562\n", + "2026-06-12 10:50:24,536 - dask_ml.cluster.k_means - INFO - Starting init iteration 5/ 9 , 13 centers\n", + "2026-06-12 10:50:24,602 - dask_ml.cluster.k_means - INFO - Finished init iteration 5/ 9 , 13 centers in 0:00:00.066245\n", + "2026-06-12 10:50:24,628 - dask_ml.cluster.k_means - INFO - Starting init iteration 6/ 9 , 15 centers\n", + "2026-06-12 10:50:24,694 - dask_ml.cluster.k_means - INFO - Finished init iteration 6/ 9 , 15 centers in 0:00:00.065630\n", + "2026-06-12 10:50:24,717 - dask_ml.cluster.k_means - INFO - Starting init iteration 7/ 9 , 17 centers\n", + "2026-06-12 10:50:24,781 - dask_ml.cluster.k_means - INFO - Finished init iteration 7/ 9 , 17 centers in 0:00:00.063115\n", + "2026-06-12 10:50:24,810 - dask_ml.cluster.k_means - INFO - Starting init iteration 8/ 9 , 18 centers\n", + "2026-06-12 10:50:24,870 - dask_ml.cluster.k_means - INFO - Finished init iteration 8/ 9 , 18 centers in 0:00:00.060724\n", + "2026-06-12 10:50:24,902 - dask_ml.cluster.k_means - INFO - Starting init iteration 9/ 9 , 20 centers\n", + "2026-06-12 10:50:24,966 - dask_ml.cluster.k_means - INFO - Finished init iteration 9/ 9 , 20 centers in 0:00:00.064427\n", + "2026-06-12 10:50:28,057 - root - INFO - Finished init_scalable in 0:00:03.988505\n", + "2026-06-12 10:50:28,060 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 0.\n", + "2026-06-12 10:50:30,476 - dask_ml.cluster.k_means - INFO - Shift: 0.1536\n", + "2026-06-12 10:50:30,478 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 0. in 0:00:02.418565\n", + "2026-06-12 10:50:30,479 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 1.\n", + "2026-06-12 10:50:30,524 - dask_ml.cluster.k_means - INFO - Shift: 0.0095\n", + "2026-06-12 10:50:30,525 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 1. in 0:00:00.045118\n", + "2026-06-12 10:50:30,526 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 2.\n", + "2026-06-12 10:50:30,567 - dask_ml.cluster.k_means - INFO - Shift: 0.0052\n", + "2026-06-12 10:50:30,569 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 2. in 0:00:00.042578\n", + "2026-06-12 10:50:30,570 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 3.\n", + "2026-06-12 10:50:30,606 - dask_ml.cluster.k_means - INFO - Shift: 0.0019\n", + "2026-06-12 10:50:30,607 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 3. in 0:00:00.037589\n", + "2026-06-12 10:50:30,609 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 4.\n", + "2026-06-12 10:50:30,654 - dask_ml.cluster.k_means - INFO - Shift: 0.0013\n", + "2026-06-12 10:50:30,657 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 4. in 0:00:00.047890\n", + "2026-06-12 10:50:30,659 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 5.\n", + "2026-06-12 10:50:30,706 - dask_ml.cluster.k_means - INFO - Shift: 0.0015\n", + "2026-06-12 10:50:30,708 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 5. in 0:00:00.049308\n", + "2026-06-12 10:50:30,710 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 6.\n", + "2026-06-12 10:50:30,760 - dask_ml.cluster.k_means - INFO - Shift: 0.0008\n", + "2026-06-12 10:50:30,762 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 6. in 0:00:00.051734\n", + "2026-06-12 10:50:30,764 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 7.\n", + "2026-06-12 10:50:30,810 - dask_ml.cluster.k_means - INFO - Shift: 0.0002\n", + "2026-06-12 10:50:30,812 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 7. in 0:00:00.048638\n", + "2026-06-12 10:50:30,815 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 8.\n", + "2026-06-12 10:50:30,861 - dask_ml.cluster.k_means - INFO - Shift: 0.0002\n", + "2026-06-12 10:50:30,864 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 8. in 0:00:00.049983\n", + "2026-06-12 10:50:30,866 - dask_ml.cluster.k_means - INFO - Starting Lloyd loop 9.\n", + "2026-06-12 10:50:30,920 - dask_ml.cluster.k_means - INFO - Shift: 0.0000\n", + "2026-06-12 10:50:30,922 - dask_ml.cluster.k_means - INFO - Finished Lloyd loop 9. in 0:00:00.055872\n" ] }, { - "name": "stderr", + "name": "stdout", "output_type": "stream", "text": [ - "/Users/rvv/Github/sidpy/sidpy/base/num_utils.py:54: RuntimeWarning: invalid value encountered in divide\n", - " if var / step_avg < tol:\n" + "Calculating cluster means and fitting priors...\n" ] } ], "source": [ - "output = fitter.do_fit(bounds = (lb,ub))" + "output = fitter.do_fit(use_kmeans=True, n_clusters=4, fit_parameter_labels=['Amplitude', 'Resonant Frequency', 'Quality Factor', 'Phase'])" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "{0: X: X (nm) of size (5,),\n", - " 1: Y: Y (nm) of size (5,),\n", + "{0: X: X (nm) of size (2,),\n", + " 1: Y: Y (nm) of size (2,),\n", " 2: DC_Offset: DC_Offset (V) of size (64,),\n", " 3: Field: Field (generic) of size (2,),\n", " 4: Cycle: Cycle (generic) of size (2,),\n", - " 5: fit_parms: fit_parameters (a.u.) of size (4,)}" + " 5: fit_parameters: Label (generic) of size (4,)}" ] }, - "execution_count": 12, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "#Let's check the fitted dataset.\n", - "output[0].data_type = 'spectral_image'\n", - "\n", - "output[0]._axes" + "output.data_type = 'spectral_image'\n", + "output._axes" ] }, { "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "http://127.0.0.1:8787/status\n" - ] - } - ], - "source": [ - "address = fitter.client.dashboard_link\n", - "print(address)" - ] - }, - { - "cell_type": "code", - "execution_count": 14, + "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "2d2efed20806492db2bb444defb0b723", + "model_id": "5f14a9641eb74b8d93a850332241e0f4", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "HBox(children=(VBox(children=(Output(),), layout=Layout(width='45%')), VBox(children=(Output(), VBox(children=…" + "HBox(children=(VBox(children=(Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': \"Canvas(t…" ] }, "metadata": {}, @@ -461,14 +411,14 @@ "# Example synthetic dataset\n", "# ==============================================\n", "\n", - "data = output[0]\n", + "data = output\n", "\n", "# Axes\n", - "X_vals = output[0]._axes[0].values\n", - "Y_vals = output[0]._axes[1].values\n", - "DC_vals = output[0]._axes[2].values\n", - "field_vals = output[0]._axes[3].values\n", - "cycle_vals = output[0]._axes[4].values\n", + "X_vals = output._axes[0].values\n", + "Y_vals = output._axes[1].values\n", + "DC_vals = output._axes[2].values\n", + "field_vals = output._axes[3].values\n", + "cycle_vals = output._axes[4].values\n", "fit_labels = ['Amplitude', 'Resonant Frequency', 'Quality Factor', 'Phase'] #need to add into the fit labels\n", "\n", "# ==============================================\n", @@ -563,7 +513,7 @@ "# Wrapper to convert voltage to index\n", "def update(dc_value, field_value, cycle_value, fit_value):\n", " dc_index = int(dc_value)\n", - " cycle_index = int(cycle_value)-1\n", + " cycle_index = int(cycle_value) ## -1\n", " plot_visualizer(dc_index, field_value, cycle_index, fit_value)\n", "\n", "# ==============================================\n", @@ -592,204 +542,131 @@ }, { "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(5, 5, 64, 2, 2, 4)" - ] - }, - "execution_count": 15, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "output[0].shape" - ] - }, - { - "cell_type": "code", - "execution_count": 15, + "execution_count": 11, "metadata": {}, "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 15, - "metadata": {}, - "output_type": "execute_result" - }, { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "803243fd92da42b7bc59ad46653be3a4", + "model_id": "d6019616f87b4179b800ea36bddea5d7", "version_major": 2, "version_minor": 0 }, - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATRJJREFUeJzt3Qd4VVW6//E3vReSQEIgAQIJJQ3LgCioIAKK9DJ3RmeccRxnbDNKFQsIUgQFy1zr3GuZ8h+lCAIqCoigAiJoCoRASCihhBAgnfTzf9aK5IICUpKsc87+fp5ny0rOIedlG3J+7LXXel1sNptNAAAAYBmupgsAAABA8yIAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMe6mC3BkdXV1cvjwYQkICBAXFxfT5QAAgItgs9mkpKREIiMjxdXVotfCbE5q/fr1tjvuuMPWunVrm/pjLl269KzH7777bv35M4+BAwde0mvk5ub+5GtwcHBwcHBwOMaRm5trsyqnvQJYVlYmycnJcs8998jIkSPP+ZxBgwbJ22+/3fCxl5fXJb2GuvKn5ObmSmBg4BVWDAAAmkNxcbFERUU1vI9bkdMGwNtuu00fF6ICX0RExGW/xulpXxX+CIAAADgWFwvfvmXRie96X3zxhbRq1Uo6d+4s999/vxw/fvyCz6+srNT/ajjzAAAAcDSWDYBq+vcf//iHrF27VubOnSvr16/XVwxra2vP+3vmzJkjQUFBDYe6fAwAAOBoXNSNgGKBS7xLly6V4cOHn/c5OTk50rFjR1mzZo3ccsst570CqI4f30NQVFTEFDAAAA5CvX+rCzlWfv+27BXAH4uJiZGwsDDZs2fPBe8ZPH2/H/f9AQAAR0UA/MHBgwf1PYCtW7c2XQoAAECTctpVwKWlpWddzdu7d6+kpKRISEiIPqZPny6jRo3Sq4Czs7Nl0qRJ0qlTJxk4cKDRugEAAJqa0wbArVu3St++fRs+HjdunP717rvvltdee03S0tLk3XfflcLCQr0T+IABA+SZZ5655L0AAQAAHI0lFoE0FW4iBQDA8RTz/s09gAAAAFZDAAQAALAYAiAAAIDFEAABAAAshgAINKFFW3Plf77MEdZaAQDsidNuAwOYtmH3MZm4OE2P/b3c5b96RJsuCQAAjSuAQBMoOlUtk34If8ozKzMk90S50ZoAADiNAAg0genLd0hecYW0D/WVa9u1kLKqWpmwKFXq6pgKBgCYRwAEGtmq7XnywfeHxNVFZP7Y7jJ/bLL4errJN3tPyDsb95kuDwAAAiDQmApKK+WJpel6/KebOso17VpIu1A/efz2rvpzc1dlyp78UsNVAgCsjgAINBK10leFv+NlVdIlIkAe6R/b8NidPaOlT2yYVNbUyfhFqVJTW2e0VgCAtREAgUayLOWQfLrjqLi7uuhpXy93t4bHXFxcZN7oJAnwdpfU3EJ5fX220VoBANZGAAQawZGiUzL1wx16/NdbYiU+Mugnz2kd5CPTh8br8Utrs2TH4aJmrxMAAIUACDTC1K/a8qWkokaS2wbJ/Td3PO9zR1zVRgbGh0t1rU3GL0yVypraZq0VAACFAAhcoX9/c0C+zCoQL3dXverX3e38f63UVPCsEYkS6ucpmXkl8tKarGatFQAAhQAIXIH9x8tk9sc79XjSoC7SqZX/z/6eMH8vmTUiQY/VvYDb9p9s8joBADgTARC4TLV1Nr25c3lVrfTsECK/v779Rf/eQQmt9XSw2hdafY1TVUwFAwCaDwEQuExvfbVXvt13Uvw83eT5McniqnZ+vgRPD4mXiEBv2VtQpvcHBACguRAAgcuw+2iJPPfZLj1+8o5uEhXie8lfI8jXQ+aOTtJj1SFk456CRq8TAIBzIQACl6i6tk6v4K2qqZObO7eU//pF1GV/rZviWupNopWJi9OkuKK6ESsFAODcCIDAJXpl3R5JP1QkQT4eMndUkl7ZeyVUm7joEF85VHhKZq7MaLQ6AQA4HwIgcAnSDxbJf3++R49nDIuX8EDvK/6afl7u+h5ClSMXbj0oazKONkKlAACcHwEQuEgV1bUybmGK1NTZ5PbECBmaHNloX7tHhxC5t3cHPX7sg3Q5UVbVaF8bAIAfIwACF+mF1bslK79Uwvw9ZebwxCue+v2x8QM6630EC0or5akPtzfq1wYA4EwEQOAifLvvhLz5ZY4ezxmZJCF+no3+Gt4ebrJgbLK4ubrIR2lHZHnq4UZ/DQAAFAIg8DPKq2r0Zs02m8joa9rKrd3Cm+y1ktoGy0N9O+nxU8u2y9HiiiZ7LQCAdREAgZ+hWr3tP14ukUHeMnVItyZ/vYf6dZKENoFSdKpaHluSJjaVPAEAaEQEQOACNuw+Jv/afECP541OlkBvjyZ/TQ83V1kwtrt4urnKul3HZOHW3CZ/TQCAtRAAgfNQV+AmLU7T49/2aie9Y8Oa7bXjwgNk/IA4PZ6xIkNyT5Q322sDAJwfARA4j+krdkhecYW0D/WVx27r0uyvf2+fGLm2XQspq6qViYtTpa6OqWAAQOMgAALn8OmOPPngu0Pi6iIyf2yy+Hq6N3sNajWw2iDax8NNNuec0P2CAQBoDARA4EeOl1bK4x+k6/F9N3aUa9qFGKulfZifPD64qx7PXZUp2cdKjdUCAHAeBEDgDGrF7RNLt8vxsirpHB4gj94aa7okuatntPSJDZPKmjoZtzBVamrrTJcEAHBwBEDgDB+mHJZVO/LE3dVFT/16ubuZLkl3HJk7KkkCvN0lNbdQ3thQvyE1AACXiwAI/CCvqEKm/tCC7S+3xEpCmyCxF5HBPvL0kHg9fnHNbsk4XGy6JACAAyMAAj9M/U5ekibFFTWS3DZIHri5o9ibkVe3kQHdwqW61ibjFqZIZU2t6ZIAAA6KAAiIyH+25Mr63cfE091VT/26u9nfXw01FTx7ZKLuQ5yZVyIvrckyXRIAwEHZ37sc0MwOHC+XmR9l6PGkgZ2lU6sAsVdh/l4ye0SCHr++Plu+O3DSdEkAAAdEAISl1dbZZMKiVCmvqpUeHULknhs6iL0blNBahnePFLUv9ISFqXKqiqlgAMClIQDC0t7+eq9s2XdCfD3dZP6YZHFVOz87gOlDEyQ80EtyCsr0/oAAAFwKAiAsK+toicz7dJcePzm4m0SF+IqjCPL10FvDKKpDyMbsAtMlAQAcCAEQllRdWyfjF6VKVU2d3BTXUn7VI0oczc2dW8mve0br8cRFaVJSUW26JACAgyAAwpJeXZctaQeLJNDbXV9JUytsHdHjt3eVqBAfOVR4Smau3Gm6HACAgyAAwnK2HyqSv31ev4XKM8MTJCLIWxyVv5e7PD86WVR+fX9rrnyeedR0SQAAB0AAhKVUVNfqTZRr6mxyW0KEDE2OFEfXMyZU/vDD6uXJS9LlZFmV6ZIAAHaOAAhLeWHNbtl9tFTC/D1l5vAEh536/bEJev9CfzlWUilTl+8wXQ4AwM4RAGEZW/edkDc35OjxnJFJEurvJc7C28NNFoxNFjdXF1mRelhWph02XRIAwI4RAGEJ5VU1etWvzSYy6uq2cmu3cHE2SW2D5cG+nfT4yWXbJb+kwnRJAAA7RQCEJcz5OFP2Hy+X1kHeMnVIN3FWD/frJPGRgVJYXi1TlqSLTSVeAACsEgA3bNggQ4YMkcjISH2f17Jly856XL0xTp06VVq3bi0+Pj7Sv39/ycqqXxkK5/Jl1jH55+b9evzc6GQJ8vEQZ+Xh5ioLxnYXTzdXWZuZL4u2HTRdEgDADjltACwrK5Pk5GR55ZVXzvn4vHnz5OWXX5bXX39dvvnmG/Hz85OBAwdKRQXTZs6k6FS1TFqcpse/7dVOeseGibPrHBEg4wbE6fGMFRly8GS56ZIAAHbGaQPgbbfdJjNnzpQRI0b85DF19e/FF1+UJ598UoYNGyZJSUnyj3/8Qw4fPvyTK4VwbCoAHSmqkPahvvLYbV3EKv7YJ0auaddCSitrdACuq2MqGABggQB4IXv37pW8vDw97XtaUFCQ9OzZUzZt2nTe31dZWSnFxcVnHbBfn+3IkyXfHRRXF5H5Y5PF19NdrEKtBp4/Jll8PNxkY/Zx+cemfaZLAgDYEUsGQBX+lPDws1eCqo9PP3Yuc+bM0UHx9BEV5Xj9Y63ieGmlPL40XY//eKO6GhYiVtM+zE8ev73+quezqzIl51ip6ZIAAHbCkgHwck2ZMkWKiooajtzcXNMl4RzUFL/aBqWgtEriwv3l0f7198NZ0Z0920nvTmFSUV2nt8Gpqa0zXRIAwA5YMgBGREToX48ePbtvqvr49GPn4uXlJYGBgWcdsD/LUw/LJ9vzxN3VRa+IVZskW5Wrq4vMG50kAV7u8v2BQnnjh42wAQDWZskA2KFDBx301q5d2/A5dT+fWg3cq1cvo7XhyuQVVchTy7br8cP9YiWhTZBYXWSwj0wbGq/HL67ZLTuPcO8qAFid0wbA0tJSSUlJ0cfphR9qfODAAb0v4COPPKJXCS9fvlzS09Plt7/9rd4zcPjw4aZLxxVM/U5ekibFFTWS1DZIHujb0XRJdmPU1W1095PqWps8+n6KVNUwFQwAVua0yyK3bt0qffv2bfh43Lhx+te7775b3nnnHZk0aZLeK/C+++6TwsJC6d27t6xatUq8vb0NVo0r8Z8tubJ+9zHxdHfVK2DVpsiop/7RM3tEomzbf1Iy80rk5bVZMmFgZ9NlAQAMcbHRK+qyqWljtRpYLQjhfkCzDhwvl0EvbZDyqlp5cnBXubdPjOmS7NLH6UfkgX9/p7fGWXL/9XJVdAvTJQFAsyvm/dt5p4BhHWqT4wmLU3X469E+RH5/QwfTJdmt2xNby7DukaL2hR6/MFVOVdWaLgkAYAABEA7vra/3ypa9J8TX002eH5OsN0HG+c0YmiDhgV6SU1Am8z7NNF0OAMAAAiAc2p78Epn36S49fmJwV4kO9TVdkt0L8vWQZ0cl6fHbX++TjdkFpksCADQzAiAcltrUeNzCVL2i9ca4lvLrHtGmS3IYfTu3kl/9cL4mLkqTkopq0yUBAJoRARAO69UvsiXtYJEEervLvFFJeqUrLp66YhoV4iOHCk/JrI92mi4HANCMCIBwSNsPFemtTJQZwxIkIojtey6Vv5e7PDc6WVRufu/bXPk88+zOOAAA50UAhMOprKnVK1hr6mwyKD5Cr2rF5bkuJlTu+WHV9OQl6XKyrMp0SQCAZkAAhMN5YXWW7DpaIqF+njJrRAJTv1do4sDO0rGlnxwrqZSpy3eYLgcA0AwIgHAo2/afkDc3ZOvx7JGJEurvZbokh+ft4SYLxnbX2+esSD0sK9MOmy4JANDECIBwGOVVNXrqV21iPPLqNjIwPsJ0SU4jOSpYHry5vnfyU8u2S35JhemSAABNiAAIhzH3k0zZd7xcWgd5y7Qh8abLcToP9YuVbq0D5WR5tUxZki50iQQA50UAhEP4ek+BvLtpvx7PHZUkQT4epktyOp7urrLgl8ni6eYqazPzZdG2g6ZLAgA0EQIg7F5xRbVMXJSqx3ddF603fUbT6BIRKI/eGqfHM1ZkyMGT5aZLAgA0AQIg7J4KIoeLKqRdqK9Mua2r6XKc3n03xsjV0cFSWlkjkxanSZ266RIA4FQIgLBrqzOOyuJtB/Vmxc+PSRY/L3fTJTk9tRp4/tju4u3hKhuzj8s/N9dPvQMAnAcBEHbrRFmVTPkgXY/v6xMjv2gfYroky+gQ5tdwtXXOJzsl51ip6ZIAAI2IAAi7pFagPrksXQpKKyUu3L/hvjQ0n99c105u6BQqFdV1Mn5RqtQyFQwAToMACLu0PPWwfJyeJ+6uLnqTYrVZMZqXq6uLzBudLAFe7vL9gUJ544cNuAEAjo8ACLtztLhCpn5Y35Ls4X6xktAmyHRJltUm2EemDummxy+s3i07jxSbLgkA0AgIgLC7qd/JS9Kk6FS1JLYJkgf61nengDmjr2kr/buGS3WtTcYtTJWqmjrTJQEArhABEHblvW9z5Ytdx+o3JR6bLB5ufIua5uLiIrNHJkgLXw99BfBvn2eZLgkAcIV4d4XdyD1RLjNXZujxxAGdJTY8wHRJ+EGrAG+ZOTxRj1/9IltScgtNlwQAuAIEQNgFtdnwhEWpUlZVKz3ah8g9vTuYLgk/MjiptQxNjtSrgcctTJGK6lrTJQEALhMBEHbh7Y375Ju9J8TX001v+Kw2I4b9mTEsXloFeEnOsTKZt2qX6XIAAJeJAAjj9uSXyrxVmXr8xOCuEh3qa7oknEewr6fMHZWkx299vVc2ZR83XRIA4DIQAGFUTW2djF+YIpU1dXJjXEv5dY9o0yXhZ/Tt0kp+1SNKjycuTtU9gwEAjoUACKNe+yJbUg8WSaC3u8wblaRXnML+PTG4m7Rt4SMHT56SWR/VL9wBADgOAiCM2X6oSF5aW7+lyPRh8RIR5G26JFwkfy93fa+m8p8tubIuM990SQCAS0AAhBGVNbUyfmGq1NTZZFB8hAzv3sZ0SbhE18WEyj031K/WVpt3F5ZXmS4JAHCRCIAw4oXVWbLraImE+nnKzBEJTP06qEmDOkvHln6SX1LZ0L4PAGD/CIBodtv2n5A3N2Tr8eyRiRLm72W6JFwmbw83mT+2u962Z3nqYfko7YjpkgAAF4EAiGZVXlWjp37rbCIjr2ojA+MjTJeEK9Q9KlgeuLm+Z/OTy9Ilv6TCdEkAgJ9BAESzmvtJpuw7Xi4Rgd4ybWi86XLQSB7uFyvdWgfKyfJqefyDdLHZbKZLAgBcAAEQzebrPQXy7qb9ejxvdJIE+XiYLgmNxNPdVRb8Mlk83Fxkzc58WbztoOmSAAAXQABEsyiuqJaJi1L1+M6e0XrTZziXLhGB8uitcXo8Y0WGHCo8ZbokAMB5EADRLFQgOFxUIdEhvvL47V1Nl4Mm8qcbO8pV0cFSUlmjA3+dutkTAGB3CIBocqszjuopQbXTi9o82M/L3XRJaCJqNfCCsd3F28NVNmYfl39urp/yBwDYFwIgmtSJsiqZ8kG6Hv+xT4z06BBiuiQ0sQ5hfjLltvqrvHM+2Sl7C8pMlwQA+BECIJqMWgn61LLtUlBaKbGt/GXcD/eHwfn95rp2cn3HUKmorpPxC1OklqlgALArBEA0mRVpR+Sj9CNnTAu6mS4JzcTV1UWeG5OsewZ/d6BQ3vhh428AgH0gAKJJHC2u0Ff/lIf6dpLEtkGmS0IzaxPsI1OHdNPjF1dnSWZesemSAAA/IACiSaZ+H1uSJkWnqiWhTaA81K+T6ZJgyJhr2kr/rq2kqrZOxr2fKlU1daZLAgAQANEU3v82V9btOla/OfDY7uLhxreZVbm4uOh+zy18PSTjSLH87fMs0yUBAAiAaGy5J8rlmZUZejxhQJzEhQeYLgmGtQrwlpnDE/X41S+yJSW30HRJAGB5BEA0GrXp74RFqVJWVSu/aN9C/tA7xnRJsBODk1rLkORIvRpYrQquqK41XRIAWBoBEI3mnY375Ju9J8THw01v+KxW/wKnPTMsXloGeEn2sTJ57tNdpssBAEsjAKJR7MkvlbmrMvX48cFdpV2on+mSYGeCfT1l3qgkPX7r672yOee46ZIAwLIIgLhiNbV1Mn5RqlTW1Emf2DC5q2e06ZJgp/p2aSX/9YsosdlE3y5QWlljuiQAsCQCIK7Y6+uzJTW3UAK83WXe6CS98hM4nycGd9V7BB48eUpmfVS/YAgA0LwsHQCffvppHVbOPLp06WK6LIey43CRvLS2fmuP6UPjpXWQj+mSYOcCvD30PaLKf7bkyrrMfNMlAYDlWDoAKvHx8XLkyJGG46uvvjJdksOorKmV8QtTpbrWJgO6hcuIq9qYLgkOolfHUPn9De31ePKSNCksrzJdEgBYiuUDoLu7u0RERDQcYWFhpktyGC+tUe29SiTEz1Nv9svULy7F5EFdJKaln+SXVMq05TtMlwMAlmL5AJiVlSWRkZESExMjd955pxw4cOC8z62srJTi4uKzDqvatv+kvvdPmT0iQcL8vUyXBAfj7eEm88cki9ot6MOUw/Jx+hHTJQGAZVg6APbs2VPeeecdWbVqlbz22muyd+9e6dOnj5SUlJzz+XPmzJGgoKCGIyoqSqzoVFWtXsFZZxM97TsoobXpkuCgropuIfff3FGPn1iaLsdKKk2XBACW4GKzqQ0ZoBQWFkq7du1kwYIF8oc//OGcVwDVcZq6AqhCYFFRkQQGBopVPL18h970OSLQWz595EYJ8vUwXRIcWFVNnQz976/07QT9u4bL3397DbcTAGhSxcXF+kKO1d6/z2TpK4A/FhwcLHFxcbJnz55zPu7l5aW/Uc48rGbjngId/pS5o5MIf7hinu6u8sIvu4uHm4us2XlUlnx3yHRJAOD0CIBnKC0tlezsbGndminNcymuqJaJi9P0+Nc9o+WmuJamS4KT6No6UB7pH6fH05fvkEOFp0yXBABOzdIBcMKECbJ+/XrZt2+fbNy4UUaMGCFubm7yq1/9ynRpdmnmygz9xhwd4itP3N7VdDlwMn+6MUauig6Wksoambw4TerUTaYAgCZh6QB48OBBHfY6d+4sY8eOldDQUNm8ebO0bMmVrR9bu/OoLNx6UNStWWoTXz8vd9Mlwcm4u7nqVcHeHq7y1Z4C+dc3+02XBABOy9Lv4u+9957pEhzCybIqeeyDdD2+t3cH6dEhxHRJcFIxLf31/oDTV2TInI8zpU9sS+kQ5me6LABwOpa+AoiL8+SH2/X2HJ1a+cv4AZ1NlwMnd3ev9tIrJlROVddvN1TLVDAANDoCIC5oReph+SjtiLi5usiCsWp6zs10SXByrq4u8tyYJPH3ctcbjv/9yxzTJQGA0yEA4rzyiyvkqQ+36/FDfTtJUttg0yXBItq28JWpd3TT4wWf7ZZdeefenB0AcHkIgDgntT+4uu+vsLxaEtoEykP9OpkuCRYz5tq20q9LK6mqrZNxC1P0htEAgMZBAMQ5LdyaK59n5ounm6ssGKs26eVbBc1LdQN5dmSiBPt6yI7DxfLf6869QTsA4NLxro6fyD1RLjNWZOjx+AFxEhceYLokWFSrQG+ZOTxBj19Zt0dScwtNlwQAToEAiLOozXcnLk6VsqpaubZdC7m3T4zpkmBxdyRFyh1JrfVq4PGLUqWiutZ0SQDg8AiAOIvq87s554T4eLjpDZ/V6l/AtGeGJUjLAC/Zk18qz3+6y3Q5AODwCIBokH2sVOauytTjxwd3lfZswAs70cLPU+aOStTj//16r3yTc9x0SQDg0AiA0Gr0SstUqaypkz6xYXJXz2jTJQFn6dclXH55bZTYbCITFqdKaWWN6ZIAwGERAKG9sSFH32Af4O0uc0cl6RWYgL158o6u0ibYR3JPnJJZH+00XQ4AOCwCICTjcLG8uGa3Hj89JF4ig31MlwScU4C3h+4SovxnywH5Yle+6ZIAwCERAC2usqZWb7JbXWuTAd3CZeTVbUyXBFzQ9R3D5Pc3tNfjyUvSpKi82nRJAOBwCIAW99KaLMnMK5EQP0+ZNSKRqV84hEkDu0hMmJ8cLa6Uacvr2xUCAC4eAdDCvjtwUl5fn63Hs4bXb7MBOAIfTzd5fmyyqF2KlqUclk/Sj5guCQAcCgHQok5V1cqEhalSZxMZ3j1Sbktsbbok4JJcHd1C7r+5ox4/sWy7HCupNF0SADgMAqBFqf3+cgrKJDzQS6YPrW+1BTiav9wSK10iAuREWZU8vjRdbGqPGADAzyIAWtDGPQW644eitnwJ8vUwXRJwWbzc3WTB2O7i4eYiqzOOygffHTJdEgA4BAKgxZRUVMvExWl6/Kse0XJz51amSwKuSLfIQHmkf5weP718hxwuPGW6JACwewRAi5m5cqccKjwlUSE+8sTgrqbLARrFn26Mke5RwVJSWSOTFqcxFQwAP4MAaCGfZx6V97fmitrp5fnRyeLv5W66JKBRuLu5yvyxyeLt4Spf7SmQf23eb7okALBrBECLOFlWJZOXpOvxH27oID1jQk2XBDSqji39ZfKgLno8++NM2VdQZrokALBbBECLeOrD+m0yOrb0kwkDO5suB2gSd/dqL71iQuVUda1MWJQqtWqfIwDATxAALWBF6mFZmXZE3Fxd9IpJbw830yUBTcLV1UXmjU7Stzds3X9S/ufLHNMlAYBdIgA6ufziCn31T3nw5o6SHBVsuiSgSUWF+MpTd9QvcJr/2W7ZlVdiuiQAsDsEQCemVkJO+SBdCsurJT4yUB7qF2u6JKBZjL02Svp1aSVVtXUybmGKVNfWmS4JAOwKAdCJLdp6UNZm5ounm6ue+vV05383rMHFxUWeHZkowb4esuNwsfzt8z2mSwIAu0IicFIHT5bLjJUZejxuQJx0jggwXRLQrFoFesszw+rbHL6ybo+kHSw0XRIA2A0CoBOqq7PJxEVpUlpZI9e0ayF/7BNjuiTAiCHJkTI4qbVeDTxuYapUVNeaLgkA7AIB0An9Y9M+2ZRzXHw83GT+mGS9+hewqpnDEiTM30v25JfK/M92mS4HAOwCAdDJ5BwrlWdXZerxlNu7SPswP9MlAUa18POUuaMS9fh/vtor3+QcN10SABhHAHQiNbV1Mn6Rmuaqk96dwuSunu1MlwTYhVu6hsvYa9uKahE8YXGqlFXWmC4JAIwiADqRNzbkyPcHCiXAy13mjk7Sm+ICqPfUHd2kTbCP5J44JbM+3mm6HAAwigDoJHYeKZYX1+zW42lD4/UbHYD/E+DtIc+NTtLj//fNAVm/+5jpkgDAGAKgE6iqUZvdpkp1rU36dw2XUVe3MV0SYJeu7xQmv7u+vR5PWpwqReXVpksCACMIgE7g5bVZ+gpgC18PmTMyUW+CC+DcJg/qIh3C/ORocaU8vWKH6XIAwAgCoIP7/sBJefWL+i4Hs0YkSssAL9MlAXbNx9NN5o9NFnWL7NLvD8mq7UdMlwQAzY4A6MBOVdXK+IWpUmcTGdY9Um5PbG26JMAhXB3dQv58U0c9fnzpdikorTRdEgA0KwKgA5v3aabkFJRJqwAvmT403nQ5gEP5a/9Y6RIRICfKquTxD9LFpvaIAQCLIAA6qI3ZBfL21/v0WG35EuzrabokwKF4ubvJgrHdxcPNRT7LOKqngwHAKgiADkj1+FW9fpVf9YiWvp1bmS4JcEjdIgPlkf5xejxt+Q45XHjKdEkA0CwIgA5o5soMOVR4SqJCfOSJwV1NlwM4tD/dGCPdo4KlpKJGJi9JYyoYgCUQAB3Musx8ee/bXFE7vTw3Oln8vdxNlwQ4NHc3V70q2MvdVb7MKpB/fXPAdEkA0OQIgA6ksLxKX6FQ7rmhg1wXE2q6JMApdGzpr/cHVOZ8vFP2Hy8zXRIANCkCoAOZ+uEOyS+plI4t/WTiwM6mywGciuoQcl1MiJRX1cqERalSq/ZXAgAnRQB0EB+lHZHlqYfFzdVFr1z09nAzXRLgVFxdXRpuq/h230n5369yTJcEAE2GAOgA8ksq5Mll6Xr84M0dJTkq2HRJgFOKCvGVp+6oX1j1/Ke7ZffREtMlAUCTIADaObUiccqSdDlZXi3xkYHyUL9Y0yUBTm3stVHSr0srqaqtk3ELU6S6ts50SQDQ6AiAdm7RtoOyNjNfPH9Yqejpzv8yoCm5uLjIsyMTJcjHQ7YfKpZX1tX32gYAZ2L5NPHKK69I+/btxdvbW3r27ClbtmwRe3HwZLnMWJGhx4/eGiddIgJNlwRYQqtAb3lmeIIe//fneyT9YJHpkgCgUVk6AL7//vsybtw4mTZtmnz33XeSnJwsAwcOlPz8fNOlSV2dTSYtTtNdP66ODpb7bowxXRJgKUOTI2VwUmupqbPpqeCK6lrTJQFAo7F0AFywYIH88Y9/lN///vfSrVs3ef3118XX11feeust06XJPzfvl43Zx8XHw03mj+2uV/8CaF7PDEuQMH8vycovlQWrd5suBwAajWUDYFVVlWzbtk369+/f8DlXV1f98aZNm875eyorK6W4uPisoynkHCuVOZ/s1OMpt3eRDmF+TfI6AC4sxM9T3w+o/P3LHNmy94TpkgCgUVg2ABYUFEhtba2Eh4ef9Xn1cV5e3jl/z5w5cyQoKKjhiIqKapLanvt0l1RU18kNnULlrp7tmuQ1AFyc/t3CZcw1bUW1CFYbRJdV1pguCQCumGUD4OWYMmWKFBUVNRy5ublN8jrPjkqSX/WIlnmjk/XmtADMmjqkm7QJ9pEDJ8obrs4DgCOzbAAMCwsTNzc3OXr06FmfVx9HRESc8/d4eXlJYGDgWUdTUNtPzBmZqN9wAJgX4O0hz41O0uN/bT4gG3YfM10SAFwRywZAT09Pueaaa2Tt2rUNn6urq9Mf9+rVy2htAOzP9Z3CdL9gRa3QLzpVbbokALhslg2AitoC5u9//7u8++67snPnTrn//vulrKxMrwoGgB+bPKh+UVZecYVMX77DdDkAcNksHQB/+ctfyvPPPy9Tp06V7t27S0pKiqxateonC0MAQPHxdJPnxySLujX3g+8Pyart514wBgD2zsWmms3isqhtYNRqYLUgpKnuBwRgf+auypTXvsiWUD9P+fTRG/VegQAcRzHv39a+AggAl+OR/rHSJSJAjpdVyZNLtwv/jgbgaAiAAHCJvNxVh55k8XBzkVU78mRZyiHTJQHAJSEAAsBliI8Mkr/eEqvHUz/cIUeKTpkuCXAqqv/22p1nb9WGxkMABIDL9OebOkpyVLCUVNTorWGYCgYazwtrdssf3t0qT7PivkkQAAHgMrm7ucr8Mcni5e4qX2YVyP/bcsB0SYBT2LrvhLy5IUePb+gUZrocp0QABIAr0KmVv0wa1EWPZ320U/YfLzNdEuDQyqtqZPyiVN1/e9TVbeXWbmzN1hQIgABwhX5/fXvp2SFEyqtqZeKiNKmtYyoYuFxzPs6U/cfLJTLIW6YN7Wa6HKdFAASAK+Tq6qI3iPbzdJMt+07IW1/tNV0S4JC+zDom/9y8X4/njU6WQG8P0yU5LQIgADSCqBBfeeqO+qsVz322S7KOlpguCXAoqr+2Wkyl/LZXO+kdy71/TYkACACN5Je/iJKbO7eUqpo6GbcwVapr60yXBDiMGSsy5EhRhbQP9ZXHbqu/rxZNhwAIAI3ExcVF5o5KkiAfD0k/VCSvrss2XRLgED7bkSdLvjuo+2yrTdZ9Pd1Nl+T0CIAA0IjCA71lxrB4Pf7b51mSfrDIdEmAXTteWimPL03X4z/eGCPXtAsxXZIlEAABoJENTY6UwYmtpabOJuMXpeiOBgB+Sm2e/uSy7VJQWiVx4f7yaP840yVZBgEQAJpgKviZ4QkS5u8lu4+Wygurd5suCbBLy1MPyyfb88Td1UUWjO0u3h5upkuyDAIgADSBED9PmTMyUY/f/DJHdzYA8H/yiirkqWXb9fjhfrGS0CbIdEmWQgAEgCaiOhiMvqat7migOhuUVdaYLgmwm6nfyUvSpLiiRpLaBskDfTuaLslyCIAA0ISmDummOxqozgbPfpJpuhzALvxnS66s331MPN3r+2l7uBFHmhtnHACakOpk8NyYZD1WHQ5UpwPAyg4cL5eZH2Xo8aSBnSU2PMB0SZZEAASAJnZDpzC5u1c7PVadDlTHA8CK6upsMmFxqu6b3aN9iPz+hg6mS7IsAiAANIPJt3XRHQ5Up4PpK3aYLgcw4q2v98qWvSfE19NN9892Uzs/wwgCIAA0A9XZQHU4UO93H3x3SD7dkWe6JKBZ7ckvkXmf7tLjJwZ3lehQX9MlWRoBEACaiepwcN+N9asdn1iarjsgAFag+mKr/tiqT/ZNcS3l1z2iTZdkeQRAAGhGj94aK53DA3TngyeWbtfbYQDO7rUvsiXtYJEEervrftlqs3SYRQAEgGbk5e6mp4JV54NVO/Lkw5TDpksCmtT2Q0Xy8tosPZ4xLEEigrxNlwQCIAA0P9Xx4K+3xOrx1A+3y5GiU6ZLAppEZU2tjF+Yqvti35YQIcO6R5ouCT8gAAKAAfff3FGS2wbpTgiTl6QzFQyntGD1btl1tETC/D1l5vAEpn7tCAEQAAxwd3OV+WO7i5e7q2zYfUz+35YDpksCGtW2/SfkzQ05ejxrRKKE+nuZLglnIAACgCGdWvnLxIGd9XjWRzt1hwTAGZRX1ehVv+rC9sir28jA+AjTJeFHCIAAYNA9N3SQnh1CdGeECYtSpbaOqWA4PtX3WvW/bh3kLdOGxJsuB+dAAAQAg1xdXXRHBD9PN9my74S8/fVe0yUBV+SrrAL5x6b9ejxvdJIE+XiYLgnnQAAEAMOiQnzlyTu66bHqlJB1tMR0ScBlKa6olomLU/X4N9e1kz6xLU2XhPMgAAKAHfivX0TJzZ1b6k4J6t4p1TkBcDTTl2foftftQn1lyu1dTJeDCyAAAoAdUNtjqA4Jaros/VCRvLou23RJwCX5bEeeLPnuoKidXuaPSdb9r2G/CIAAYCfCA71lxrD6G+b/9nmWpB8sMl0ScFFUX+vHl6br8X19YuTa9iGmS8LPIAACgB0ZmhwptydG6M4J4xamSEV1remSgAtSm5g/9eF23d86LtxfHr01znRJuAgEQACws6ngZ4Yl6M4JWfml8sLq3aZLAi5oeeph+Tg9T/e3nj+mu3h7uJkuCReBAAgAdkZ1TJgzMkmP3/wyR7buO2G6JOCcjhZXyFPLtuvxw/1iJbFtkOmScJEIgABgh27tFi6jrm6rOymMX5QqZZU1pksCfjL1O2lxmu5nndgmSB7o29F0SbgEBEAAsFNTh3TTnRRURwXVWQGwJ+99myvrdx8TT3dXWTA2WTzciBSOhP9bAGCn1JYwz41O1uN/bt4vX2YdM10SoOWeKJeZKzP0eOKAzhIbHmC6JFwiAiAA2LHesWG6o4KiptuKTlWbLgkWV1dnq78toapWerQPkXt6dzBdEi4DARAA7JzqqKA6K6gOC9NX7DBdDizura/3ypa9J8TX002eG5Mkbq4upkvCZSAAAoCdUx0VVGcF1WHhg+8Oyac78kyXBIvak1+i+1Urj9/eVdqF+pkuCZeJAAgADkB1Vrjvxhg9fmJpuu68ADSnmto6Gb8wVfer7hMbJnf2jDZdEq4AARAAHMSj/eN0pwXVceGJpdv1NhxAc3nti2xJPVgkAd7uMm90kt60HI6LAAgADkJ1WFgwtrvuuLBqR558mHLYdEmwiB2Hi+SltVl6rPpVtw7yMV0SrhABEAAcSEKbIN1xQZn64XbJK6owXRKcXGVNrYx7P1X3px4YHy7Du7cxXRIaAQEQAByM6riQ1DZId2CYvCSNqWA0qRfXZMmuoyUS6ucps0YkMvXrJAiAAOBgVMcFtSpYdWBQnRj+syXXdElwUtv2n5A31mfrsQp/Yf5epktCI7FsAGzfvr3+V8yZx7PPPmu6LAC4KKrzwqSBnfV45kcZcuB4uemS4GTKq2r0qt86m8jIq9rIoIQI0yWhEVk2ACozZsyQI0eONBwPP/yw6ZIA4KL9/oYOuhNDeVWtTFicqjs0AI1l7ieZsu94uUQEesu0ofGmy0Ejs3QADAgIkIiIiIbDz48NLQE4DtWB4fkxybojg+rMoDo0AI3h6z0F8u6m/XqstnxRfanhXCwdANWUb2hoqFx11VXy3HPPSU1NzQWfX1lZKcXFxWcdAGBSdKivPDG4qx6rDg2qUwNwJYorqmXiolQ9vuu6aLkxrqXpktAELBsA//KXv8h7770n69atkz/96U8ye/ZsmTRp0gV/z5w5cyQoKKjhiIqKarZ6AeB8ft2j/k1adWgYtzBVd2wALtczKzLkcFGFRIf4ypTb6v9xAefjYnOi/QMee+wxmTt37gWfs3PnTunSpctPPv/WW2/pIFhaWipeXl7nvQKojtPUFUAVAouKiiQwMLAR/gQAcHnUfoADXlivt4YZd2uc/OWW+r0CgUuxJuOo3PuPrbrv9MI/9ZJftA8RZ1RcXKwv5Fj5/dupAuCxY8fk+PHjF3xOTEyMeHp6/uTzO3bskISEBMnMzJTOnetX1v0cvoEA2JNl3x+SR95P0Z1Clj14g940GrhYJ8qqZMALG6SgtFL3nX78due9+lfM+7e4ixNp2bKlPi5HSkqKuLq6SqtWrRq9LgBoDsO6R8qq7Xm6TZzavmP5wzeIl7ub6bLgANS1oKeWbdfhL7aVv76KDOdmyXsAN23aJC+++KKkpqZKTk6O/Pvf/5ZHH31U7rrrLmnRooXp8gDgsqj9TGeNSNAdG1TnhhdW1/duBX7OirQj8lH6EX31WPWbVn2n4dwsGQDVPX5qAchNN90k8fHxMmvWLB0A33zzTdOlAcAVCfX3ktkjE/X4zQ3ZupMDcCFHiyv01T/loX6dJLEttw5YgVPdA9jcuIcAgL0atzBFPvjukLQP9ZWP/9pHfD2d6o4fNBIVAX7/zrfyxa5jktgmSD544HrdatDZFfP+bc0rgADg7KYNiZfWQd66k8Ozn2SaLgd26v1vc3X4U32l549NtkT4Qz3+TwOAE1KdG+aOStLjf2zaL19lFZguCXYm90S5PLMyQ48nDIiTuPAA0yWhGREAAcBJqc2hVScHZdLiVN3hAVBU3+gJi1KlrKpWftG+hfyhd4zpktDMCIAA4MTUXm7tQn11Z4cZK+qv9gBvb9wn3+w9oftIq37Sqq80rIUACABOTC3+mD8mWXd2WLztoKzOOGq6JBi2J79U5q3KPOMfCH6mS4IBBEAAcHLXtg+R+/rUT/FN+SBNd3yANak+0eMXpUplTZ30iQ2TO3vW3yIA6yEAAoAFPHqrusnfXwpKq+TJZel6+w9Yz+vrsyU1t1ACvN1l3ugkvXk4rIkACAAWoDo7qA4PqtPDx+l5sjz1sOmS0Mx2HC6Sl9bWd4eZPlRtE+RjuiQYRAAEAItIaBMkD/eL1eOpH+7QHSBgDZU1tbo/dHWtTQbGh8uIq9qYLgmGEQABwEIe6NtRd3woOlUtk5ekMRVsES+uyZLMvBLdJ3rWiESmfkEABAArUZ0eFoxN1p0fVAeI977NNV0Smti2/SfljfXZejxrRIKE+XuZLgl2gAAIABYTGx4gEwd01uOZKzN0Rwg4p/KqGr3hc51N9LTvoITWpkuCnSAAAoAF3dO7g/RoH6I7QeiAoBICnM68Vbtkb0GZRAR6y9ND4k2XAztCAAQAC1KdH1QHCNUJQnWEeOvrvaZLQiP7ek+BvLNxnx7PHZ0kQb4epkuCHSEAAoBFRYf6yhODu+rxvE93yZ78EtMloZGovs+TFqfpsdrs+aa4lqZLgp0hAAKAhf26R7TcGNdSqmrq9DYhqlMEHN8zKzLkUOEpiQ7x1e3egB8jAAKAhantQOaNSpJAb3dJPVgkr35Rv1oUjmtNxlFZtO2g7v+spvn9vNxNlwQ7RAAEAIuLCPKW6cPqFwi8vDZLth8qMl0SLpPq8/zYB+l6fK9a6NMhxHRJsFMEQACADO/eRgbFR0hNnU1PBavOEXA8T324XQpKK6VTK38Z/8NWP8C5EAABAHoqeOaIBN0pYtfREnlhdX3PWDgO1d/5o7QjeoW32uxb9X8GzocACADQVIcI1SZMeXNDtmzbf8J0SbhI+cUV8tSy7Xr8YN9OktQ22HRJsHMEQABAg0EJETLyqja6c4SaCladJGDfVD9n1ddZ9XeOjwyUh/t1Ml0SHAABEABwlmlD43XniH3Hy2XuJ5mmy8HPWLg1V9btOiaeus9zd93vGfg5fJcAAM4S5OMh80Yn6fG7m/brjhKwT6qP84wVGXo8fkCcdI4IMF0SHAQBEADwE2pzaNVBQpm4KFV3loB9Uf2bJy5O1f2cr23XQu7tE2O6JDgQAiAA4JxUBwnVSeJwUUXDVSbYj3c37ZPNOSfEx8NNb/isVv8CF4sACAA4J9VBQgUL1VFi8baDsjrjqOmS8IPsY6Xy7A/3Zz5+exdpH+ZnuiQ4GAIgAOC8VCeJP/4wtTjlg3TdaQJmqX7N9Zt110mf2DC567p2pkuCAyIAAgAuaNytcRLbyl93mFB7zaltR2DOGxtyJCW3UAK83WXuqCS9iTdwqQiAAIALUh0l1PYi6h6zj9KPyIq0I6ZLsqyMw8Xy4prdevz0kHiJDPYxXRIcFAEQAPCzEtsGyUN96zcYVlcBjxZXmC7JclR/5nELU6S61ia3dguXkVe3MV0SHBgBEABwUR7q10kS2gTqjhOPLUljKriZvbw2SzLzSiTEz1PmjExk6hdXhAAIALgoqsOEmgr2dHfVnSfe/zbXdEmW8d2Bk/LaF9l6PHtEgu7bDFwJAiAA4KLFhQfIhAFxevzMygzdiQJN61RVrUxYmKr7Mw/vHimDElqbLglOgAAIALgkf+gdI79o30J3oJiwKFV3pEDTmbsqU3IKyiQ80EumD00wXQ6cBAEQAHBJ1GpgtUG06kDxzd4T8s7GfaZLclobswsazq/a8iXI18N0SXASBEAAwCVrF+onjw/u2nCFak9+qemSnE5JRbVMXJSmx7/uGS03d25luiQ4EQIgAOCy3NUzWneiUB0pxi9K1R0q0HhmrtwphwpPSVSIj+7LDDQmAiAA4LKobUjmjU7SHSlScwvl9fX1q1Rx5T7PPCrvb83VfZifH50s/l7upkuCkyEAAgAuW+sgH5k+NF6PX1qbJTsOF5kuyeGdLKuSyUvS9fje3h2kZ0yo6ZLghAiAAIArMuKqNjKgW7juUDF+YaruWIHL99SH2+VYSaV0auUv4wd0Nl0OnBQBEABwxVPBs0cm6g4VqlPFS2uyTJfksFakHpaVaUf0SusFY5N1H2agKRAAAQBXTHWmUB0qFHUv4Lb9J02X5HDyiyv01T/lwb6dJKltsOmS4MQIgACARqE6VKjpYLUvtNogWnWwwMVRfZUf+yBdCsurJT4yUB7q28l0SXByBEAAQKN5eki8RAR6y96CMr0/IC7Ooq0H5fPMfPE8o98y0JT4DgMANBrVqWLu6CQ9Vh0sNu4pMF2S3VP9lGeszNDjcQPipHNEgOmSYAEEQABAo7oprqXc2TNajycuTpPiimrTJdkt1Ud50uI0Ka2skWvatZA/9okxXRIsggAIAGh0qnNFdIiv7mQx84erW/ipdzftk005x3Vf5fljkvXqX6A5EAABAI3Oz8tdnh+TrDtZLNx6UNbuPGq6JLuTc6y04T7Jx2/vIu3D/EyXBAtxygA4a9Ysuf7668XX11eCg8+9jP7AgQMyePBg/ZxWrVrJxIkTpaamptlrBQBn1aNDiO5koagVrqrDBeqpvsnjFqZKRXWd9O4UJnf2bGe6JFiMUwbAqqoqGTNmjNx///3nfLy2tlaHP/W8jRs3yrvvvivvvPOOTJ06tdlrBQBnpjpZqI4WqrPF6T3uIPLGhhxJyS2UAC933U/ZlalfNDOnDIDTp0+XRx99VBITE8/5+GeffSYZGRnyr3/9S7p37y633XabPPPMM/LKK6/oUAgAaByqk4XqaKHubVMdLlSnC6vLOFwsL67ZrcfThsZLZLCP6ZJgQU4ZAH/Opk2bdDgMDw9v+NzAgQOluLhYduzYcd7fV1lZqZ9z5gEAuDDV0eL0xsbqKqDqeGFVVTVq6jdF902+tVu4jLq6jemSYFGWDIB5eXlnhT/l9MfqsfOZM2eOBAUFNRxRUVFNXisAOIOH+nWShDaButOFuh9Qdb6wopfXZul+yapv8uwRibqPMmCCwwTAxx57TP9FudCRmdm0u85PmTJFioqKGo7c3NwmfT0AcBYepztcuLnqjheq84XVfH/gpLz6xR49njU8QVoGeJkuCRbmLg5i/Pjx8rvf/e6Cz4mJubgNNCMiImTLli1nfe7o0aMNj52Pl5eXPgAAly4uPEDGD4iTOZ9k6s4XvTqGSlSIr1iB6os8fmGq7pM8vHuk3JbY2nRJsDiHCYAtW7bUR2Po1auX3iomPz9fbwGjrF69WgIDA6Vbt26N8hoAgJ+6t0+MrM44Klv3n9QdMP59b09LrICd92mm5BSUSXigl0wfmmC6HMBxpoAvhdrjLyUlRf+qtnxRY3WUlpbqxwcMGKCD3m9+8xtJTU2VTz/9VJ588kl58MEHucIHAE1IrQZWG0SrzheqA4bqhOHsNmYXyNtf1/85545K0v2SAdOcMgCq/fyuuuoqmTZtmg59aqyOrVu36sfd3Nxk5cqV+ld1NfCuu+6S3/72tzJjxgzTpQOA01MdL1TnC+XZTzIl+1j9P86dkerxO3FRmh7/qke03Ny5ftYJMM3FZtWlWI1AbQOjVgOrBSFq+hgAcHHUW89v39oiX2YVSPeoYFn8517i7uZ81yQeW5Im732bK1EhPvLJX28Ufy+HufPKqRXz/u2cVwABAPZN7dygpkMDvN11RwzVGcPZrMvM1+FP7fTy/Ohkwh/sCgEQAGCE6oDx9JB4PVadMVSHDGdRWF4lk5fUT/3+4YYO0jMm1HRJwFkIgAAAY0Ze3UYGdAvXnTFUhwzVKcMZPPXhDskvqdR9kCcM7Gy6HOAnCIAAAKNTwbNHJurOGKpDxktr63vkOrKVaYd1z2O14nn+mGTdDxmwNwRAAIBRYf5eMntE/d54r32RrTtmOKr8kgp5atl2PX7w5o6SHBVsuiTgnAiAAADjBiW01h0yVKcM1TFDdc5wxJXNU5aky8nyaomPDJSH+sWaLgk4LwIgAMAuqA4ZqlOG6pihOmc4mkXbDsrazHzd71j3PXbnLRb2i+9OAIBdUB0y1NYwiuqcoTpoOIqDJ8tlxooMPR43IE46RwSYLgm4IAIgAMBuqE4Zv+4Zrceqg0ZJRbXYu7o6m+5rrLp+XNOuhfyxT4zpkoCfRQAEANiVJ27vqjtnHCo8JTNX7hR798/N+2Vj9nHd31it+lWrfwF7RwAEANgVPy93mT+mu+6g8f7WXPk886jYq5xjpTLnk/qQOuX2LrrPMeAICIAAALvTo0OI3Nu7gx5PVitry6rE3tTW2WT8olSpqK6T3p3C5K6e7UyXBFw0AiAAwC6NH9BZd9I4VlIpU5fvEHvzxga1Z2GhBHi5y7zRSeLK1C8cCAEQAGCXVAeNBWPr76lTnTVUhw17kZlXLC+sru9aMm1ovO5rDDgSAiAAwG4ltQ2WB/t20uMnl23XnTZMU/2KH30/Vfcv7t81XEZd3cZ0ScAlIwACAOzaw/066c4aheXVutOG6rhh0t8+z5KdR4qlha+HzBmZqPsZA46GAAgAsGsepztruLnqThuLth40VktKbqG8+kW2Hs8akSgtA7yM1QJcCQIgAMDuqc4a4wfE6fGMlRm680Zzq6iulXELU/Tq32HdI+X2xNbNXgPQWAiAAACHcG+fGLm2XQvdcUN1CVEdOJrTvFW7JOdYmbQK8JLpQ+Ob9bWBxkYABAA4BLUa+Pkxybrjxqac4/Lupn3N9tqbc47LW1/v1eO5o5Mk2Nez2V4baAoEQACAw1CdNh6/vYsez12VqTtxNDV1xXHColQ9/lWPKOnbuVWTvybQ1AiAAACHctd17aRPbJjuwKE6cdTU1jXp6836SN1zeEratvCRJwZ3a9LXApoLARAA4FDUtitzRyXpDhyqE8cbG3Ka7LXWZebLf7bk6rGafvb3cm+y1wKaEwEQAOBwVOcN1YFDeXHNbr0vX2MrLK+SyUvS9PieGzrIdTGhjf4agCkEQACAQ1IdOG7tFq47coxbmKo7dDSmqR/ukPySSolp6SeTBnVu1K8NmEYABAA47FTw7BGJEuLnqa8Avrw2q9G+9kdpR2R56mG98lhtQq36EgPOhAAIAHBYqhPHzOEJevzqF3vk+wMnr/hrqn7DTy5L1+MHbu4o3aOCr/hrAvaGAAgAcGiqI4fqzKH2hR6/MFVOVdVe9tdSfYYf/yBdTpZXS7fWgfJwv9hGrRWwFwRAAIDDmzE0QcIDvSSnoEzmfZp52V9n8baDsmZnvni4uciCXyaLpztvk3BOfGcDABxekK+HPDsqSY/f/nqfbMwuuOSvcajwlMxYkaHHj94aJ10iAhu9TsBeEAABAE5Bdej4VY9oPVa9gksqqi/696q+wpMWp0pJZY1cFR0sf7qxYxNWCphHAAQAOI0nBneVqBAffTVv1kc7L/r3/XPzfvl6z3Hx9nCV+WOS9epfwJkRAAEATkN16nhudLK4uIi8922u7uTxc/YWlMmcT+rD4pTbukpMS/9mqBQwiwAIAHAqqmOH6tyhqE4eqqPH+dTW2WT8whTdV/j6jqHym+vaNWOlgDkEQACA05k4sLN0bOmnO3mojh7n8+aGHPnuQGH9lcMxyeLK1C8sggAIAHA6qnOH6uCh7uVTHT1UZ48fy8wrlhdW79bjqUO6SZtgHwOVAmYQAAEATik5KlgevLl+Na/q7KE6fJym+gaPez9VqmrrpH/XVjLmmrYGKwWaHwEQAOC0HuoXK/GRgbqzh+rwoTp9KH/7PEsyjhRLC18PmT0yUfcVBqyEAAgAcFqqk8f8scni6eaqO3yoTh8puYXy6hfZ+vGZwxOlVYC36TKBZkcABAA4NdXRQ3X2UFSnj0ffT9Grf4ckR8rgpNamywOMIAACAJzefTfGyNXRwbrTh9r3r2WAlzwzLN50WYAxBEAAgNNTq4Hnj+0uPh5u+uO5oxIl2NfTdFmAMe7mXhoAgObTIcxPFv25lxSWV0vv2DDT5QBGEQABAJaR0CbIdAmAXWAKGAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAsxikD4KxZs+T6668XX19fCQ4OPudzVN/HHx/vvfdes9cKAADQ3JxyG5iqqioZM2aM9OrVS/73f//3vM97++23ZdCgQQ0fny8sAgAAOBOnDIDTp0/Xv77zzjsXfJ4KfBEREc1UFQAAgH1wyingi/Xggw9KWFiY9OjRQ9566y2x2WwXfH5lZaUUFxefdQAAADgap7wCeDFmzJgh/fr10/cJfvbZZ/LAAw9IaWmp/OUvfznv75kzZ07D1UUAAABH5WL7ucteduKxxx6TuXPnXvA5O3fulC5dujR8rKaAH3nkESksLPzZrz916lR9T2Bubu4FrwCq4zR1BTAqKkqKiookMDDwov8sAADAnOLiYgkKCrL0+7fDXAEcP368/O53v7vgc2JiYi776/fs2VOeeeYZHfC8vLzO+Rz1+fM9BgAA4CgcJgC2bNlSH00lJSVFWrRoQcADAABOz2EC4KU4cOCAnDhxQv9aW1urw53SqVMn8ff3lxUrVsjRo0fluuuuE29vb1m9erXMnj1bJkyYcEmvc3r2nMUgAAA4juIf3rcd5C64pmFzQnfffbf6P/qTY926dfrxTz75xNa9e3ebv7+/zc/Pz5acnGx7/fXXbbW1tZf0Orm5ued8HQ4ODg4ODg77P3Jzc21W5TCLQOxRXV2dHD58WAICAnQnkcZ0eoGJWpRi1RtUz4dzc2Gcn/Pj3FwY5+f8ODfOdX5sNpuUlJRIZGSkuLpac0c8p5wCbi7qm6Zt27ZN+hrqL5Ij/GUygXNzYZyf8+PcXBjn5/w4N85zfoKCgsTKrBl7AQAALIwACAAAYDEEQDultqOZNm0a29KcA+fmwjg/58e5uTDOz/lxbi6M8+N4WAQCAABgMVwBBAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxBEA79Morr0j79u11n+KePXvKli1bxIo2bNggQ4YM0Tu1q04ry5YtO+txtX5p6tSp0rp1a/Hx8ZH+/ftLVlaWWMGcOXPkF7/4he5C06pVKxk+fLjs2rXrrOdUVFTIgw8+KKGhoboH9qhRo3QPbCt47bXXJCkpqWFT2l69esknn3zS8LiVz82PPfvss/rv1yOPPNLwOSufn6efflqfjzOPLl26NDxu5XOjHDp0SO666y7951c/dxMTE2Xr1q0Nj1v557KjIQDamffff1/GjRunl9N/9913kpycLAMHDpT8/HyxmrKyMv3nV4H4XObNmycvv/yyvP766/LNN9+In5+fPlfqB7SzW79+vX4T2rx5s6xevVqqq6tlwIAB+pyd9uijj8qKFStk0aJF+vmqbeHIkSPFClSHHhVstm3bpt+c+vXrJ8OGDZMdO3aI1c/Nmb799lt54403dFg+k9XPT3x8vBw5cqTh+Oqrrxoes/K5OXnypNxwww3i4eGh/0GVkZEh8+fPlxYtWjQ8x8o/lx2O6WbEOFuPHj1sDz74YMPHtbW1tsjISNucOXNsVqa+VZcuXdrwcV1dnS0iIsL23HPPNXyusLDQ5uXlZfvPf/5js5r8/Hx9jtavX99wLjw8PGyLFi1qeM7OnTv1czZt2mSzohYtWtj+53/+h3Pzg5KSEltsbKxt9erVtptuusn217/+VX/e6udn2rRptuTk5HM+ZvVzM3nyZFvv3r3P+zg/lx0LVwDtSFVVlb5ioS6Zn9lvWH28adMmo7XZm71790peXt5Z50r1dVRT5lY8V0VFRfrXkJAQ/av6PlJXBc88P2oaKzo62nLnp7a2Vt577z19dVRNBXNu6qkryIMHDz7rPCicH9FTlurWk5iYGLnzzjvlwIED+vNWPzfLly+Xa6+9VsaMGaNvPbnqqqvk73//e8Pj/Fx2LARAO1JQUKDfrMLDw8/6vPpY/aXC/zl9PjhXInV1dfr+LTU1k5CQoD+nzoGnp6cEBwdb9vykp6fre7RUZ4I///nPsnTpUunWrRvnRkQHYnWLibqX9Mesfn5UWHnnnXdk1apV+l5SFWr69OkjJSUllj83OTk5+pzExsbKp59+Kvfff7/85S9/kXfffVc/zs9lx+JuugAAV34lZ/v27WfdpwSRzp07S0pKir46unjxYrn77rv1PVtWl5ubK3/961/1vaNqoRnOdttttzWM1b2RKhC2a9dOFi5cqBc1WJn6x6a6Ajh79mz9sboCqH72qPv91N8vOBauANqRsLAwcXNz+8mKMvVxRESEsbrs0enzYfVz9dBDD8nKlStl3bp1euHDaeocqFsKCgsLLXt+1JWaTp06yTXXXKOvdKkFRS+99JLlz42axlSLyq6++mpxd3fXhwrG6sZ9NVZXa6x8fn5MXe2Li4uTPXv2WP57R63sVVfRz9S1a9eGKXJ+LjsWAqCdvWGpN6u1a9ee9S8u9bG6dwn/p0OHDvoHypnnqri4WK86s8K5UutiVPhT05qff/65Ph9nUt9HaqXemedHbROjflBb4fyci/q7VFlZaflzc8stt+jpcXV19PShruqoe91Oj618fn6stLRUsrOzdfix+veOus3kx9tN7d69W18hVaz+c9nhmF6FgrO99957esXUO++8Y8vIyLDdd999tuDgYFteXp7NatQqxe+//14f6lt1wYIFerx//379+LPPPqvPzYcffmhLS0uzDRs2zNahQwfbqVOnbM7u/vvvtwUFBdm++OIL25EjRxqO8vLyhuf8+c9/tkVHR9s+//xz29atW229evXShxU89thjekX03r179feG+tjFxcX22Wef2ax+bs7lzFXAVj8/48eP13+v1PfO119/bevfv78tLCxMr7S3+rnZsmWLzd3d3TZr1ixbVlaW7d///rfN19fX9q9//avhOVb+uexoCIB26G9/+5v+AePp6am3hdm8ebPNitatW6eD34+Pu+++u2HLgaeeesoWHh6uQ/Mtt9xi27Vrl80KznVe1PH22283PEf9wH3ggQf09ifqh/SIESN0SLSCe+65x9auXTv9d6hly5b6e+N0+LP6ubmYAGjl8/PLX/7S1rp1a/2906ZNG/3xnj17Gh638rlRVqxYYUtISNA/c7t06WJ78803z3rcyj+XHY2L+o/pq5AAAABoPtwDCAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAWMv/B8ze10AwkqEHAAAAAElFTkSuQmCC", - "text/html": [ - "\n", - "
\n", - "
\n", - " Figure\n", - "
\n", - " \n", - "
\n", - " " - ], "text/plain": [ - "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + "VBox(children=(IntSlider(value=0, description='Freq Slice:', max=122), Canvas(header_visible=False, toolbar=To…" ] }, "metadata": {}, "output_type": "display_data" } ], - "source": [ - "plt.figure()\n", - "plt.plot(output[0]._axes[2])" - ] - }, - { - "cell_type": "code", - "execution_count": 60, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[('Amplitude', 0),\n", - " ('Resonant Frequency', 1),\n", - " ('Quality Factor', 2),\n", - " ('Phase', 3)]" - ] - }, - "execution_count": 60, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "[(name, idx) for idx, name in enumerate(fit_labels)]" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [ - { - "ename": "NameError", - "evalue": "name 'beline_sidpy' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[31m---------------------------------------------------------------------------\u001b[39m", - "\u001b[31mNameError\u001b[39m Traceback (most recent call last)", - "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[14]\u001b[39m\u001b[32m, line 88\u001b[39m\n\u001b[32m 85\u001b[39m \u001b[38;5;28mself\u001b[39m.fig.tight_layout()\n\u001b[32m 86\u001b[39m \u001b[38;5;66;03m# Example usage\u001b[39;00m\n\u001b[32m 87\u001b[39m \u001b[38;5;66;03m# Assume raw_data and fit_data are your input matrices\u001b[39;00m\n\u001b[32m---> \u001b[39m\u001b[32m88\u001b[39m raw_data = np.array(\u001b[43mbeline_sidpy\u001b[49m) \u001b[38;5;66;03m# Replace this with your actual raw data\u001b[39;00m\n\u001b[32m 89\u001b[39m fit_data = np.array(fit_results[\u001b[32m0\u001b[39m]) \u001b[38;5;66;03m# Replace this with your actual fit data\u001b[39;00m\n\u001b[32m 91\u001b[39m \u001b[38;5;66;03m# Instantiate the visualizer\u001b[39;00m\n", - "\u001b[31mNameError\u001b[39m: name 'beline_sidpy' is not defined" - ] - } - ], "source": [ "\"\"\"\n", - "THis visualizer works for BELIne data. Just input a BEline sidpy dataset and an SHO fit dataset from sidpy fitter\n", - "Some things need work\n", - "Get X,Y units from sidpy dataset\n", - "Plot a scale bar and add a checkbox for whether you want it\n", - "Add export plot button\n", - "Make a simple viz for the plot parameters only, and enable figures to be output.\n", - "Get the frequency vector from the sidpy dataset\n", - "Then, start work on a BEPS visualizer.\n", + "Interactive per-pixel SHO fit-quality visualizer for BE-line / BEPS-slice data.\n", + "Shows: raw 2D heatmap at one frequency | raw vs fitted amplitude | raw vs fitted phase.\n", + "Click on the heatmap to pick a pixel; move the slider to change the frequency slice.\n", "\"\"\"\n", "\n", - "%matplotlib notebook\n", "import matplotlib.pyplot as plt\n", - "from mpl_toolkits.mplot3d import Axes3D\n", "import numpy as np\n", - "from ipywidgets import interact, widgets\n", + "from ipywidgets import widgets, VBox\n", + "from matplotlib.patches import Rectangle\n", + "from IPython.display import display\n", "\n", "\n", "class InteractiveVisualizer:\n", " def __init__(self, raw_data, fit_data, freq_vec):\n", - " self.raw_data = raw_data\n", - " self.fit_data = fit_data\n", - " self.freq_vec = freq_vec\n", - " self.fig, (self.ax1, self.ax2, self.ax3) = plt.subplots(1, 3, figsize=(12, 4))\n", - " \n", + " self.raw_data = raw_data # (X, Y, freq) complex\n", + " self.fit_data = fit_data # (X, Y, 4) real (Amp, w_0, Q, phi)\n", + " self.freq_vec = freq_vec # (freq,) real\n", + "\n", " self.x = 0\n", " self.y = 0\n", " self.freq_slice = 0\n", "\n", - " # Connect the button press event handler to the 2D image plot\n", + " # Build figure with ipympl-friendly pattern: ioff() so plt.subplots\n", + " # doesn't auto-display in a separate output area.\n", + " plt.ioff()\n", + " self.fig, (self.ax1, self.ax2, self.ax3) = plt.subplots(1, 3, figsize=(12, 4))\n", + " self.fig.canvas.header_visible = False\n", + " plt.ion()\n", + "\n", + " # Click handler on the heatmap\n", " self.fig.canvas.mpl_connect('button_press_event', self.on_click)\n", "\n", - " self.freq_slice_slider = widgets.IntSlider(min=0, max=len(freq_vec), value=0, description='Freq Slice:')\n", - " interact(self.update_plots, freq_slice=self.freq_slice_slider, x=widgets.fixed(0), y=widgets.fixed(0))\n", + " # Slider (max is INCLUSIVE in IntSlider, so last valid index = len - 1)\n", + " self.freq_slice_slider = widgets.IntSlider(\n", + " min=0, max=len(freq_vec) - 1, value=0, description='Freq Slice:'\n", + " )\n", + " self.freq_slice_slider.observe(self._on_slider_change, names='value')\n", + "\n", + " # First draw, then display slider + figure together\n", + " self.update_plots(self.freq_slice, self.x, self.y)\n", + " display(VBox([self.freq_slice_slider, self.fig.canvas]))\n", + "\n", + " def _on_slider_change(self, change):\n", + " self.update_plots(change['new'], self.x, self.y)\n", "\n", - " self.update_plots(self.freq_slice, self.x,self.y)\n", - " \n", " def on_click(self, event):\n", - " if event.inaxes == self.ax1:\n", - " self.x, self.y = int(event.xdata + 0.5), int(event.ydata + 0.5)\n", + " if event.inaxes == self.ax1 and event.xdata is not None and event.ydata is not None:\n", + " nx, ny = self.raw_data.shape[0], self.raw_data.shape[1]\n", + " self.x = int(np.clip(round(event.xdata), 0, nx - 1))\n", + " self.y = int(np.clip(round(event.ydata), 0, ny - 1))\n", " self.update_plots(self.freq_slice, self.x, self.y)\n", "\n", " def update_plots(self, freq_slice, x, y):\n", " self.freq_slice = freq_slice\n", - " # Plot the slice of the raw data\n", + "\n", + " # --- Raw 2D heatmap at this frequency slice ---\n", " self.ax1.clear()\n", - " self.ax1.imshow(np.abs(self.raw_data[:, :, freq_slice]))\n", - " self.ax1.set_title('Raw Data')\n", + " self.ax1.imshow(np.abs(self.raw_data[:, :, freq_slice]), origin='lower')\n", + " self.ax1.set_title(f'Raw |amp| @ freq idx {freq_slice}')\n", " self.ax1.set_xlabel('X')\n", " self.ax1.set_ylabel('Y')\n", - " self.ax1.add_patch(Rectangle((self.x-0.5, self.y-0.5), 1, 1, linewidth=2, edgecolor='red', facecolor='none'))\n", - "\n", - " \n", - " amp_data = np.abs(self.raw_data[x,y,:])\n", - " phase_data = np.angle(self.raw_data[x,y,:])\n", - " \n", - " fit_parms = self.fit_data[x,y,:]\n", - " \n", - " sho_fit = SHO_fit_flattened(freq_vec, *fit_parms)\n", - " \n", - " sho_fit_complex = sho_fit[:len(sho_fit)//2] + 1j*sho_fit[len(sho_fit)//2 :]\n", + " self.ax1.add_patch(Rectangle(\n", + " (x - 0.5, y - 0.5), 1, 1, linewidth=2, edgecolor='red', facecolor='none'\n", + " ))\n", + "\n", + " # --- Raw vs fitted amplitude/phase at selected pixel ---\n", + " amp_data = np.abs(self.raw_data[x, y, :])\n", + " phase_data = np.angle(self.raw_data[x, y, :])\n", + "\n", + " fit_parms = self.fit_data[x, y, :]\n", + " sho_fit = SHO_fit_flattened(self.freq_vec, *fit_parms)\n", + " n = len(sho_fit) // 2\n", + " sho_fit_complex = sho_fit[:n] + 1j * sho_fit[n:]\n", " amp_fit = np.abs(sho_fit_complex)\n", " phase_fit = np.angle(sho_fit_complex)\n", - " \n", - " # Plot the spectrum at the selected X, Y position from raw data and fit data\n", - " #Plot amplitude\n", - " \n", + "\n", " self.ax2.clear()\n", - " self.ax2.plot(self.freq_vec, amp_data,'ro', label = 'Raw Data')\n", - " self.ax2.plot(self.freq_vec, amp_fit,'r-', label = 'Fit')\n", - " self.ax2.set_title('Amplitude')\n", - " self.ax2.set_xlabel('Frequency')\n", + " self.ax2.plot(self.freq_vec, amp_data, 'ro', label='Raw')\n", + " self.ax2.plot(self.freq_vec, amp_fit, 'r-', label='Fit')\n", + " self.ax2.set_title(f'Amplitude @ ({x},{y})')\n", + " self.ax2.set_xlabel('Frequency (Hz)')\n", " self.ax2.set_ylabel('Amplitude (a.u.)')\n", - " \n", - " #Plot phase\n", + " self.ax2.legend()\n", + "\n", " self.ax3.clear()\n", - " self.ax3.plot(self.freq_vec, phase_data,'ro', label = 'Raw Data')\n", - " self.ax3.plot(self.freq_vec, phase_fit,'r-', label = 'Fit')\n", - " self.ax3.set_title('Phase')\n", - " self.ax3.set_xlabel('Frequency')\n", - " self.ax3.set_ylabel('Amplitude (a.u.)')\n", - "\n", - " # Update the figure\n", - " self.fig.canvas.draw()\n", - " self.fig.tight_layout()\n", + " self.ax3.plot(self.freq_vec, phase_data, 'ro', label='Raw')\n", + " self.ax3.plot(self.freq_vec, phase_fit, 'r-', label='Fit')\n", + " self.ax3.set_title(f'Phase @ ({x},{y})')\n", + " self.ax3.set_xlabel('Frequency (Hz)')\n", + " self.ax3.set_ylabel('Phase (rad)')\n", + " self.ax3.legend()\n", + "\n", + " try:\n", + " self.fig.canvas.draw_idle()\n", + " except AttributeError:\n", + " pass\n", + "\n", "# Example usage\n", - "# Assume raw_data and fit_data are your input matrices\n", - "raw_data = np.array(beline_sidpy) # Replace this with your actual raw data\n", - "fit_data = np.array(fit_results[0]) # Replace this with your actual fit data\n", + "# InteractiveVisualizer is BE-line-only (expects 3D inputs).\n", + "# For BEPS data, pick one (DC, field, cycle) slice.\n", + "raw_data = np.array(beps_small[:, :, :, 0, 0, 0]) # (X, Y, freq) at DC=0, field=0, cycle=0\n", + "fit_data = np.array(output[:, :, 0, 0, 0, :]) # (X, Y, 4 params) same slice\n", "\n", "# Instantiate the visualizer\n", - "visualizer = InteractiveVisualizer(raw_data, fit_data, freq_vec)\n", - "\n" + "visualizer = InteractiveVisualizer(raw_data, fit_data, freq_vec)" ] }, { @@ -816,7 +693,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.4" + "version": "3.10.0" } }, "nbformat": 4,