diff --git a/.zenodo.json b/.zenodo.json index de384312..604945ea 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -31,6 +31,10 @@ "name": "Magnus Nord", "orcid": "0000-0001-7981-5293", "affiliation": "Norwegian University of Science and Technology" + }, + { + "name": "Zachary Varley", + "affiliation": "Carnegie Mellon University" } ] -} \ No newline at end of file +} diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 1af0a279..12a864d1 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -19,6 +19,8 @@ Unreleased Added ----- +Fast PCA Dictionary Indexing + Changed ------- diff --git a/doc/tutorials/pattern_matching_fast.ipynb b/doc/tutorials/pattern_matching_fast.ipynb new file mode 100644 index 00000000..2bffe829 --- /dev/null +++ b/doc/tutorials/pattern_matching_fast.ipynb @@ -0,0 +1,329 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "nbsphinx": "hidden" + }, + "source": [ + "This notebook is part of the `kikuchipy` documentation https://kikuchipy.org.\n", + "Links to the documentation won't work from the notebook." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Pattern matching\n", + "\n", + "Crystal orientations can be determined from experimental EBSD patterns by matching them to simulated patterns of known phases and orientations, see e.g. Chen et al. (2015), Nolze et al. (2016), Foden et al. (2019).\n", + "\n", + "In this tutorial, we will perform *dictionary indexing* (DI) using a small Ni EBSD data set and a dynamically simulated Ni master pattern from EMsoft, both of low resolution and found in the [kikuchipy.data](../reference/generated/kikuchipy.data.rst) module.\n", + "The pattern dictionary is generated from a uniform grid of orientations with a fixed projection center (PC) followng Singh and De Graef (2016).\n", + "The true orientation is likely to fall in between grid points, which means there is always a lower angular accuracy associated with DI.\n", + "We can improve upon each solution by letting the orientation deviate from the grid points.\n", + "We do this by maximizing the similarity between experimental and simulated patterns using numerical optimization algorithms.\n", + "This is here called *orientation refinement*.\n", + "We could instead keep the orientations fixed and let the PC parameters deviate from their fixed values used in the dictionary, here called *projection center refinement*.\n", + "Finally, we can also refine both at the same time, here called *orientation and projection center refinement*.\n", + "The need for orientation or orientation and PC refinement is discussed by e.g. Singh et al. (2017), Winkelmann et al. (2020), and Pang et al. (2020).\n", + "\n", + "The term *pattern matching* is here used for the combined approach of DI followed by refinement.\n", + "\n", + "Before we can generate a dictionary of simulated patterns, we need a master pattern containing all possible scattering vectors for a candidate phase.\n", + "This can be simulated using EMsoft (Callahan and De Graef (2013) and Jackson et al. (2014)) and then read into kikuchipy.\n", + "\n", + "First, we import libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2023-05-13T12:26:10.880617050Z", + "start_time": "2023-05-13T12:26:09.630506890Z" + } + }, + "outputs": [], + "source": [ + "# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting\n", + "%matplotlib inline\n", + "\n", + "import tempfile\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import hyperspy.api as hs\n", + "import kikuchipy as kp\n", + "from orix import sampling, plot, io\n", + "from orix.vector import Vector3d\n", + "\n", + "\n", + "plt.rcParams.update({\"figure.facecolor\": \"w\", \"font.size\": 15})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load the small experimental nickel test data" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [ + "nbval-ignore-output" + ], + "ExecuteTime": { + "end_time": "2023-05-13T12:26:18.263046637Z", + "start_time": "2023-05-13T12:26:10.882308954Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2023-05-13 08:26:12,186 - numexpr.utils - INFO - Note: NumExpr detected 16 cores but \"NUMEXPR_MAX_THREADS\" not set, so enforcing safe limit of 8.\n", + "2023-05-13 08:26:12,187 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.\n", + "2023-05-13 08:26:12,535 - numba.cuda.cudadrv.driver - INFO - init\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[########################################] | 100% Completed | 100.95 ms\n", + "[########################################] | 100% Completed | 403.06 ms\n", + "Dictionary size: 100347\n", + "[########################################] | 100% Completed | 1.61 ss\n" + ] + } + ], + "source": [ + "# Use kp.load(\"data.h5\") to load your own data\n", + "s = kp.data.nickel_ebsd_large(allow_download=True) # External download\n", + "s.remove_static_background()\n", + "s.remove_dynamic_background()\n", + "# load master pattern\n", + "energy = 20\n", + "mp = kp.data.nickel_ebsd_master_pattern_small(\n", + " projection=\"lambert\", energy=energy\n", + ")\n", + "ni = mp.phase\n", + "rot = sampling.get_sample_fundamental(\n", + " method=\"cubochoric\", resolution=2.0, point_group=ni.point_group\n", + ")\n", + "print(\"Dictionary size: \", rot.shape[0])\n", + "det = kp.detectors.EBSDDetector(\n", + " shape=s.axes_manager.signal_shape[::-1],\n", + " pc=[0.4198, 0.2136, 0.5015],\n", + " sample_tilt=70,\n", + ")\n", + "sim = mp.get_patterns(\n", + " rotations=rot,\n", + " detector=det,\n", + " energy=energy,\n", + " dtype_out=np.float32,\n", + " compute=True,\n", + ")\n", + "signal_mask = ~kp.filters.Window(\"circular\", det.shape).astype(bool)\n", + "p = s.inav[0, 0].data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's use the [pca_dictionary_indexing()](../reference/generated/kikuchipy.signals.EBSD.pca_dictionary_indexing.rst) method to match the simulated patterns to our experimental patterns. Let's keep the best matching orientation alone. The results are returned as a [orix.crystal_map.CrystalMap](https://orix.readthedocs.io/en/stable/reference/generated/orix.crystal_map.CrystalMap.html)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2023-05-13T12:27:16.660506988Z", + "start_time": "2023-05-13T12:27:09.423828612Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "Experiment loop: 0%| | 0/2 [00:00", + "image/png": "iVBORw0KGgoAAAANSUhEUgAABJ4AAAFCCAYAAACw1xcZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdIklEQVR4nO39ebyWZbn//x+LeZRJRgGZHCASJcXMVEyrrabuTDRTERQ1J8yxj7NuNdui4ZAjGiQOkZVaqduyNEwrNZQE0QRFQWYQZFqwgPX9Y/+2v+p4awf3eV73fYOv5+OxH4/P5+i8rvNai7Xue62r1euqqa+vrzcAAAAAAAAgswaVvgAAAAAAAABsnbjxBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAACsGNJwAAAAAAABSCG08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgEJw4wlV78orr7Samhqrqalx/9n/zf/x/xo0aGBt27a1wYMH24UXXmjvvfeePG+vXr3k8f/6fyNGjAhf69tvv22tW7e2mpoaO/fccz9x7fz5861Dhw5WU1Njo0aNCu8BAPCq/b3i6quvtpqaGmvSpIm9+uqr//bjue2226ympsYaNmxokydPjnwKAABCtb8//KsJEyaEzqv+D6hWjSp9AUAOLVu2tFatWpmZ2caNG23JkiX2yiuv2CuvvGJ33HGHTZo0yQ4++GB5bLNmzaxNmzYfe+5P+s/+VZ8+fWzMmDF22mmn2c0332xf//rXbZ999pFrR40aZcuWLbOePXvaD37wg/AeAIDSVPK94uKLL7Zf//rX9uKLL9rw4cPt5ZdftiZNmsi1s2bNsu9+97tmZnbOOefYvvvuG/nwAAAlqpbfJczMmjdvbp07dw6t/b9rNTNr2rTpZu0DlFU9UOWuuOKKejOrV1+u/ze/4oor/mn+4Ycf1t955531rVu3rjez+latWtUvWLDgn9Zsv/329WZWf8IJJ2S/5i9/+cv1Zlbfp0+f+lWrVrn//O677643s/qampr6p59+Ovv+APBpsyW8V7zxxhv1zZs3rzez+u9+97tyzcaNG+v33nvvejOrHzBgQP3atWuT9wWAT7Mt4f2hVGedddZHH8O9995bsesA/h3+p3bYKrVu3dpOPfVUGzt2rJmZrVq1yiZMmFC2/e+9917bZptt7O2337YLLrjgn/6z2bNn23nnnWdmZqeffrodcMABZbsuAMD/X7nfK3baaSf7/ve/b2ZmY8aMsRdeeMGtufHGG+3555+3Ro0a2cSJE61Zs2aFXQ8AQKv07xIREyZMsFtvvdXMzM4880w78cQTK3xFwMfjxhO2ascee6w1aPC/X+YvvfRS2fbt0aOH3XTTTWZmduedd9pvf/tbMzOrr6+3ESNG2MqVK61fv352/fXXl+2aAABaOd8rzjrrLNt///1t06ZNdsIJJ9jq1as/+s9ef/11u+yyy8zM7LLLLrPBgwcXei0AgE9Wqd8l/p2XXnrJTjvtNDMz22+//T66QQZUK248YavWrFkz69Chg5mZffjhh2Xde+TIkXbIIYdYfX29nXTSSbZixQq76aab7A9/+IM1aNDAJkyYYC1atCjrNQEAvHK+V9TU1NiECRNsm222sZkzZ9qFF15oZmYbNmyw4cOH27p162z33Xe3iy++uNDrAAD8e5X8XeLjLFy40I444girra21nj172sMPP2yNGpFuRnXjxhO2aqtXr/4ouNe+ffuy7z9u3Dhr3769zZkzx4455hi75JJLzMzs3HPPtb333rvs1wMA8Mr9XtGzZ0+7+eabzczsjjvusN/+9rd27bXX2l//+ldr1qyZTZw4kV8iAKAKVPp3iX9VV1dnw4YNs7lz51rz5s3tkUcesY4dO1b6soB/i59qsFW74447rL6+3szMPv/5z8s1kyZNsv/5n/+R/9nAgQPt6aefLnn/rl272q233mrHHnusPfnkk2ZmNmDAALvmmmtKPicAIK9KvFeMGDHCHnnkEfvlL39pxx9/vC1btszMzK677jrbeeedN+tcAIBiVPp3iX919tln23PPPWdm//tfcPM/ycaWghtP2Ops3LjR3nnnHXvggQfsuuuuM7P//W8oTjjhBLm+trbWamtr5X/WpUuX5OsZNmyYnXvuubZw4UIzMxs7diyPOwWACquG94q7777bXnjhhY/eH4YOHWpnn312SecCAORRDe8Pyr333mt33HGHmf3v/3ri2GOPzXZuoGj8T+2wVbjqqquspqbGampqrFGjRrbDDjvYlVdeaevWrbOOHTvaY489Zu3atZPHnnDCCVZfXy//79VXX02+tmuvvfajXyrMzO65557kcwIANl+1vVd07tzZLr300o/+/zfddJPV1NSUdC4AQOmq7f3hX/35z3+2008/3czMDjzwQB5QhC0Of/GErULLli2tVatWZmbWoEEDa9WqlfXp08cOOOAAO/HEEz+KAqY6++yzbdKkSfI/W7BggZtNmTLFrr32WjMzO/jgg+2JJ56whx9+2B5++GEbNmxYlmsCAMRU43tFmzZt5P8bAFA+1fj+8H/mz59v3/jGN2z9+vXWu3dvmzRpkjVs2DDL9QDlwo0nbBXOP/98u/LKKwvfZ8WKFf/010ufZN26dXbCCSfYhg0bbM899/yo4/HQQw/ZGWecYUOHDiUGCABlVI3vFQCAyqvW94f169fbkUceafPmzbOWLVvao48+WhWRc2Bz8T+1AzbDhAkTPvZPaf/VFVdcYdOmTbNmzZrZj3/8Y2vYsKH98Ic/tC5dutjixYs/+nNZAMDWZXPeKwAAnx6b+/5w5pln2gsvvGBmZuPHj7dddtmlnJcLZMONJ6AAf/7zn+2GG24wM7Pvfe97ttNOO5nZ/4YJ77zzTjMz+9nPfmY//elPK3aNAAAAAKrTnXfeaePGjTMzs4suuohMB7Zo3HgCMlu7dq2NGDHCNm7caPvss497QtHhhx/+0VMozjjjDFu0aFElLhMAAABAFXr++edt9OjRZmZ20EEH2TXXXFPhKwLScOMJyOziiy+2N99801q2bGnjx4+3Bg38t9ktt9xiXbt2tSVLlvA/uQMAAABgZmbvv/++HXnkkVZXV2c77LCDPfjgg/L3CWBLwlcwkNHkyZPt5ptvNjOz66+/3vr27SvXtW/f3u666y4zM/v5z3/+sU+3AAAAAPDpcfXVV3/0hLv58+fbzjvvbF26dAn93//1oIBqw1PtgExWr15tI0eOtPr6evvSl75kp5122ieuP/TQQ+3444+3iRMn2plnnmn777+/derUqUxXCwAAAKDa1NbWfvT/XrVqla1atSp87Pr164u4JCBZTT2PWAEAAAAAAEAB+J/aAQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIXgxhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAAChEo+jC42yom21nm9zsv2/8fOlXU3Oln81o4mf91/vZprf87KV1fnaE+JDnNPSzBjuEruW2ce/7dZltMnF9CWY+0d3N+h08N+seKb72LX99Ub2/9oWMV5JoxUY3mn7eX9ys2Zpxebe1Q92szha5WWPrVPIegy84uORjo6aMeaLwPaJqzX+fN7TGbqY+pz/d1MXNjmqwILTvtAH+2AZr/Lqe4Vdyb9+ZpR/r7Og/JybeJ/Z6w39vRP1p5AV+2GkbP1v0oRs9/M4YN6t5xx+6vrOfNW7hZ8N6x67lP4+7zK/LbI8f5D1fa/G2s7L0f7bsmvz4syUfe8Tw1zJeSZomA/ysx777uNm4Z5/Luu/Jbw33w6lL/WxQh5L3GDR+z5KPjZo60r+nVszy1X72bq2fic/pb66+z82+cpn4N1L2ONrP5oh1q/xrYtitx5R+7D840fy1dhP/Hfhd1rfkPRavEz9Trazzs9bi/aphLzda02h7N2tuf3Szmg09/Pk2zo5dy8zmfpbZV+fn/X3l9K/7N4TbH8n7+0qK9jeUfuxD956f70JStRFvvK+JH/paLsm67aS9n3ezZeZ/YGxv/Ure4wG7u+Rjo461UwrfI2q1rXWzdeL3C/U5PfoJ/8PCpINfD+37WfHv1tj869V21iZ0PqWF/S20jr94AgAAAAAAQCG48QQAAAAAAIBCcOMJAAAAAAAAheDGEwAAAAAAAApRU19fXx9ZqOLiUfffOCS28I3vxdbtND22bsEqP+vSys9USLxZ8J5crQ/nXnquDzh2NRE6D8odF692DSxWsD3jwTyxyyxESHzWaaPcrIO417tUxJfD24qQeDmoWPmeF4wofN9qCo5Hqbi4ooLjq/r5YxuLGHGd6Au+t8HPNolA9vBYmzBGxsVj9npDBFeFP3336tC6h/8SC3pvEB9/I/E5liHx7uLhF8pc/0CMXU72ofM+rWOnU3LHxavdS+fG1t14Y7HXsTlUSPx/xMv/qjf8rNXOpe8rQ+LlIGLlg6YU/2CKqgqOB6m4uCKD460O8rOFIiTeWTyEoZWYiUa2nX7Iv722CBUXj/qViZ/VhcVLbo+dsG3sfCsa+Z/b2mwQ37gqJN4gGHTfJP4dZosHwHywLHY+IXdcvNo91XW72MK2VfS7hAqJv73cz9qJuPgH4oe+IBUSLwcVK3/Sfl/4vtUUHI9ScXFFBccPtWZutkI8X66N+a+h922Fm9WJp1cMCP4uy188AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIUQdbL8jjvvRTcLB8cP/JufzRGx7fFiNjIYEldeCX5qdvMhrpSQOOLmfeurbrbWupV8vmj4OyUQHl2XO0KeW2Pz0UsV/h58gQjJNhLl0g0+VKeo86l9VfxcXbOijk0xW8TFeyXc8lchceW4r53kh0uPEiv991El/GlnHyaPBscffMaHxBuLkPr85/ys6z5+JkPiynofbJS6+1FKSBxxj4hnjLTZvfTzRcPfcl3wJTwaEs8dIc9uUAc3UuHvQeP3dLOGttLNNlrsm0adTwbHRfxcXbOkjk2wSLztdIq9ZWkqJK68utaN9mp/pJv9KeFScjnU3nKzaHC8flv/Wl0jWsx/bbSrm31uw6t+oQqJKyuax9a1ERHyD8ryKxrWiNeVtbNLP180/K3WvR0MhEdD4pkj5Lm1t35upsLfD9jdbtbAtnWzTbYktK86n9pXxc/VNSvq2CTb7uFnS1aXfDoVEld2GO9fm7qP9A+0mBvcl794AgAAAAAAQCG48QQAAAAAAIBCcOMJAAAAAAAAheDGEwAAAAAAAApRU19fXx9ZeJwNLfhSzO6/4WY/bPeBnz3Ux89OjgXF7P1YeXP+H32ZtmsHEeLqvz50vtvOXez3sKZ+DxEm39NOCO0R9Uu7x806WjCSWwYNbGNoXQfz0cY97DE3a27z3CwlQp4iJRq+wg7NfTlOSpQ7SgbHM1PB8Ur5qYiLzxb/5Bc2WuBmq/r5YxsP8MdGg+PKvjn7hzv6QHhuD3fxwfEuvkNsr4mX5rbBzuu3PhOMiy8S3xsDa8U6EY4V/vPcMW72tvjYVJj8Ov/Sl+T3Isq9fGHePVK8dG5sXYv5frb/FD9b8bKfpUTIU6REw09+a3jei1FSotxBKhCemwyOV8hvrr7PzVRw/Lhbxb9vKx96tYXiNScaHFduPab0Y//BiXZ0lvN8kh9tmO2Hi/3n4+2uX3Cz3vZkaI+adbGfU61hfz9bLh6a1Dr43qkefNSuvZ99sMyNrnvm/dgeQS+28bM1XbNukeSprtvFFnbwv4vZov38rHkvP0uJkKdIiIZP2vv5zBfjpUS5o1QgPDcVHK+Uo58QP/xv29KNJg15yc0OFb/jrxDPl4sGx5UWJl7XBP7iCQAAAAAAAIXgxhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAACiEqJOVx/33Xhlb+EE7P/sPERxXgiFxpeuBIhCron4zmoTOd8YPOrrZy+f6yHIj6yGO9rMNNid4rHeEXRVap/zRvuNm9da55POlWGprQ+teF0G73rYmdGw0Bh5dFw2JR6XEwGvNhyabWTDGmHC+v4yZ4GZ7XjCi5H2VB3bzUe5jX/Hx7qjc0fVeCbf8Z4l+3yb/LARph7xtx8I919KHxJUFIrbtX3HN6sS3fTgkrgxe7mfrxYMaOsWivo/+wF/LjL4+OL5urj/WP5bCrGl3sU4cq+ythqLBqjz2rp818t+SZRGN3W5UzzwQcWclGgOProuGxMNSYuDLV5d+bPR8bX0gdepg/4CIQVPyPphi1lgf9O57TkKcPXN0vVPKcz36im/+2I+IZs8GXySqwI/mizfEjuL1Vsz6bJjmZvXix/xwSFypVwF78ZSLlbH3OtvNj677wWtu1ki8DanZBvGrjlqnfEG96cyOHftr8WtDs+DDP7Jbqj4QYbV4iFX0Z7loDDy6LhgSj0qJga8Wv4ulhMTV+VqK75mD7Etu9qT9vuR9lV8v+7ubfa39jiWfL3t0fYl4Pw1aZf796gNxb0FpZ7u6WfDXEP7iCQAAAAAAAMXgxhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAAChETX19fX1k4XE2NOvG999wc2jdeQ8McrMbj53qFzbYofSLaSbuv9UmRKBFcPzVcSNLP5+QEhfPTV3Ln+3G0LEqTN7ASg85dhABurZ2npvtImJua61byfuWw0Ib7GYpMfBqouLduYPjKmquRKPh0eB4Spj89V6+yLxehJv79Q2dTtrXtw5Lt2PjjCcze7hLLLg6+Jd+NuUwPxvWOyEk3l08SGLu+tLPJ4Ljr664rPTzCSlx8dzUtfz0+dixKkz+0rmlX0uL+X721bf9bPHTftZm99L3LYeTX/qGH4p49xZJxLtzB8dV1FyKRsOjwfGUMHmXo/1sd/HEhWf9A0DCbj2m9GP/wYkmrjXBjzbMDq07/kMf+Z64zZ5+4cbY+aQGH/rZptiDJSQRHB8zcVnp5xNS4uK5qWv5n+1jx6ow+VNdE34+7iCepPE38abTW3zfrp1d+r5l8KO9/RubindviVS8O3dwXEXNlWg0PBocTwmTDxaPn9lG/O7e2gaGzqe0MPGgB4G/eAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIXgxhMAAAAAAAAK0Sj3CbczH+V+P3h/64gzfUj8F+ZD4nLd7WtCe6jwt+22IbaufzAkG10XFA2J5w6Op5zvi3aTm/3RvhM6dr75qF9XEUZTltpaN+tjPnisQuIdxNfpUvH1XClbS0hcica2U6hYuQqOqxi4kjskro7tqarKIiQ+c5afqeC4Wrdv4NryU99XsfeJh/9TDEVIXK0bJp5LIYnwt5kon6p1i0RgVlHrRM80KhoSzx0cTznf4SIc+5iI5ytvr/SzPqKnrKzp6mcrn/QzFRJf9Yaftdo5tm9ZbC0hcSUa207ZQsTKZXBcxcDlCTOHxNWxrcRrybNiNlT8DKGC42pdmXUT7wfzgj+Pdf6LD4kv3NOHxOW63f0DbyQR/rY2wXWtgw/hiK4LiobEcwfHU873tYV+9uvgP5G1a+9nHwTj7EvF7xzq30OFxNuJX60/EL9nVsjWEhJXorHtFCpWroLjKgau5A6Jq2O3k9Fw/4K10qa5mQqOq3UtQlfHXzwBAAAAAACgINx4AgAAAAAAQCG48QQAAAAAAIBCcOMJAAAAAAAAhaipr6+vjyw8zoaWvMma2mfd7Bc/FJXXl0TMbb9d/Kw2IfgcjYu/EuyuB499edyBbpYS/k5Rjli5Olb5k93vZg1sY+jYqFPsT262zGKhv1UiVl4OK+zQrOerNR8Q3RJj5YMv8OHX3KaMESHZzKIR8lX9fFx8tngJ65XwPIN9Y/3DmB1Lj6Ee9W0fYR32qF9X846fHTn8Ej+cm/BJUdHwJqKQuj5YSA0eO6P1ZW6WEv5OUY5YuTpW+YWI4r90buzYqEcX7+GHwSDsuLdeyXsxQSe/NTzvCZev9rMtMFY+aLwPSec2daQPU+ffJBghb3WQn/UXr2Ezgg8+UG49pvRj/8GJdnTJxz7+/CQ3U9HwNY38Uwta1P3Kn3CTeohEkIyLi58XVwRDzsFjrxvn49gp4e8U5YiVq2OVZ8UDHp7qmvdn3F/+xf8c3Tr4YIn9W43Mei1Rk/Z+Puv5VovfibbEWPkDdnfhexxrpxS+RzRCfqj5b7jZIhDeS0bIY1rY30Lr+IsnAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDBgnYaGRJfsMrP9tzNz1Tku7+Ixkaj4SkhcbXvK2JfsS4lyp07Ql6OkHh0j33sIjf7mY1xs64mwvNB0ZC40kpE88oRHM8dA98SQ+LKX8ZMcLM9L/JBZtsQi9srKmAeDY5Ho+Fqpo4183FxZaaIL/frGzq0aqiQ+IbX/eyYQy7ww/UqEiter6PRcEsIiS8SAV+1r1i3bkVsi5Sgd1Q5QuLRPQ5p6meTVvpZn9axa5GCIXHl5B38zy5lCY7njoFvgSFxZepg/3o9eMo+brbRSv+CUQHzcHA8Gg1XM3Xs3rFtbaj4OeBZ/7NGtVIh8RWN/H9/3mbdH/zBK0UAWf3zR6PhbcSx0ZB4a/EQDvXaL9alRLlzR8jLERKP7nHgbHHCAe397AMfZ4+KhsSVZ1aNd7NyBMdzx8C3xJC4cpB9yc2eEnHsTbak5D1UwDwaHI9Gw9VMHWvBaPhKERxvnRAcV/iLJwAAAAAAABSCG08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQ2ePi29mm3Kf0mon7ZSr8re6r1QavT0bIYwHzl0/3ceJoqDt3SDyqmuLnnRNi4B2qKHzXQXz9LRXfHyvsUDdrKEPTXjRmXe2iH8fYDbu42TnXXe1mlQqOp3zu1bEqqd9LvNSplGB1B8fL8D7RXbxezxXhb7lOhclF5VSFyYMB8xl9/UMULBjqzh0Sj8odEk9x9Dg/e+nc2LEt5ue9lhSr3vCzViJie/Jbw/3wXfH1pwLh0Zh1tQt+HPff87ibTRl8iJtVLDie8rmPHjtDvNZ1Fh9blQbHu5XjvxdvID5HrdWDKkRcfJNaFwyTq5C4WHfdDT6EHQ115w6JR1VT/Pyrr7/mZk91DT5op4N4okXw4R/ZtRO/qouHYUza+3k3WyceVNQyGKRW4epqF/04TnvxGH/wEP/7RaWC4ymf++ixvUQ0XMXFcwfH+YsnAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDhuPgaEZXexYa42ZU3/D9/8EsikbuNjzBaNxUIF1QgfEYs/C3D5DIkLj41wXUpse1ovDv3upRriYpeX4r97TdZzxeVEhJXopHqLTEkrqiPQwXHezXoEjo2JSQeFQ2OR8nA+gD/8da97o9t1N/PNswo+VJK11C8PjT0r/8Pt/fBy5p3/KHHjPTh+GB3/2MC4SoIKyqnKjiu1q0PBsfFunWLxRZB0Xh37nUp1xIVvT4T8fyoB1fuVvrBCZJC4ko0NL0lhsQV9XGI4Hgn9bYojk0JiUeFg+NRKrDeT4XExWvdwpViXfGfg3/U2Fq5WTtb4GbXbPAB3zWNtnezNkue9pu0bhy7GBUIX6kC4eJYFSZX61aIB97I4LhflxLbjsa7c69LuZao6PWleObtt/OeMCohJK5EQ9NbYkhcUR+HCo7bkgGhY1NC4lHR4HiU+nhX2K5u1kY8yEu8u1junx74iycAAAAAAAAUghtPAAAAAAAAKAQ3ngAAAAAAAFAIbjwBAAAAAACgEOG4uHLlxGAA843z3OjEQc+52Y9MxGBV0FvdL1PrVIRcUcHx/uJaXhHBWbEuJcqdsi53DLxSvmLfd7OJdqmbDRRfB6vMRxtbmYg7JlDnWyr2hVmt+Yh0M9sudGxKOF1FvlUMPLfBF53qr+W6u0LHqo93lQiJK73Ey5XIKdp76uW0YPue4r8GzvqZX7fnnKP98Pof+Nmim/xMBb2jgXAVIVfU+RaJwKwKmIt1Tfv6ZdEod8q63DHwShkt4ux7i3bywFV+Nu6tV9zs5B3yBsfV+caZ3xdmtny1n7VtGTs2IZyuIt8qBp7b4PE+bDtlZPDFXn28rdRDE4QZKkIu4uLR82WyzLqJqY/6NrAvutmQbQe52YvrbvGnU0Fv9XObCoSrCLmiguMqdL5CHCvWpUS5U9bljoFXyld/LYYXtvezBv71Z/+NI93smVXjM1zVJ59vf/P7wmy1+F5tGfzdLiWcriLfKgae20P2Czc7xo4IHas+XhUSV3rZQDdbadPc7H3xIrZDaAf+4gkAAAAAAAAF4cYTAAAAAAAACsGNJwAAAAAAABSCG08AAAAAAAAoRDgu/tRls/xwxyf97HFR9xzyjBv9qP+a2MavJPXPvWiEXAXH5bF+lDvoraLhSnTfaIQ8GjCPXkt037UiWtbV1rnZNw8d4WbnTPJfa2PP8P+Ws8b74HPuCPkKO9TNUmLbKepskZulxLuje6R8bOp8szd1cbOfmp8d1WCBm5UlOL7Bf42rjyNq5iz/sfX3X1ZWF+zS9sz8cuq8co6fiQ7jsEHf9cMLfNRQhsSV9Qnl007iWBUrVxHyaMBcyB30VtHwlH2jEfJowDx6LdF9Z4l1fUQn+ZAvTnSzY7v6B57sO+h8Nzv+exe6WfYI+VviIS0pse0UU5f6WUK8uyx7iPMtEi+5v7n6Pjf7ymX+c1+O4PhGE1+o6vMSdYl4n500w886ByPZq0QkO5Ol5t8j2tvNfuHjy9yo2d5XutmLKt6trEj5+a7Oj1SsXEXIVXBcBsz9x5E76K2i4Up032iEPBowj15LdN+N4kcK+8B/XY35/FFudk+3sW5W09JH69943j8wIneEfNLez7tZSmw7xTLx2JqUeHc59lDns233cKOjn/APfZh0sP/huhzB8U3i4Qry4wgSd2Zsk+3qZtEI+XbyRSyGv3gCAAAAAABAIbjxBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAACpGUml1/10FutkF0N1v0FiHxGSLK2l/EW1XQW4kGwnOvE7NoRDtFSvhbyR0NT6HON3PN9m6mQuLdxJfVObf5f8tG7/l1l/S8182Wmj82JUJejpC4slHEMYNJzrDcsfKo2eJbNXpL/S9jJrjZnheMSLmc0PlU6DxKhcTfEy+T/fqWvEVWkzc96oe9Putnf5vqZ51ErHWRiLWqGLgSDYTnXidmTZv6ZbmD4ynhbyV3NDyFOt+j6jVIfFnZQh+znjz1Br/uoEluNHasj8m22tkfmhQhL0dIXNk+c81YyR0rD1LB8aipg/3r9aApeR9Moc6nQudhKiTeSsye9Q88Kbe75/hosx0q3tTqRGB4pQh/q+C4ioEr0UB47nViFo1op0gJfyu5o+Ep1Pl2+qz4ehfPcujbcqObzVo92s36713vZjNeP9mf8AP/9ZwSIS9HSFxpWobfYXLHysO2Lf199yD7kps9ab9PuZrQ+VToXPNhchUSf1883Gs7U5X+0vEXTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDceAIAAAAAAEAhkuLiR/zxVTf7u8ViiH//wfGxTVTkW4kGwqOi8fMtUDQQnjuIHvXAmmdD6xotEcNusT02bOtn1753UmjdmS3GudkKOzS2cYVUKmqeQsXKe4lvcxkXT9ijHMHxwRf4kKzat1H/Lv5g8TK0qYXYpEperraf/J9u1kJdrx3nJjMWXRDbREW+lWggPErFz+eKcKxSJfH3jxMNhOcOokf1Pzz4MINN6n1MPPBEaeDfF88553OhdeMe9zXrk98aHtu3UioVNU8hYuWdxJdGSlxc7VGW4Pj4PUP72h4rYydU3wozNu+aivBBD/8EgHZzeouVf/OjTrvENlGRbyUaCI+S8fOE81WRaCA8dxA9aoJvgUsNWvpAeNSm1TVu1n/A3aF1b7460s0m7f18yddSDpWKmqeQsfIloiifEBdXe5QjOP6A+a81te9ScWxrMauTDwtrs9nX9Un4iycAAAAAAAAUghtPAAAAAAAAKAQ3ngAAAAAAAFAIbjwBAAAAAACgEDX19fWhqtoe9j03G2RD3WyyvVbyxfz99hP8MCUQniIYK1+4yldy14k2V1MRdZxz+VA3yx3+Vscq0X1ze9gmutn6NaVHzcee4f/dfrd33q+hjR39rOFiP2s1yocwVcy61t53sy0xBl4p12/wAe4LGy3Iuked+TJt7uC48qc3p/hreT3vHvvOzHeuFhf52fqFfrZDQs97xvaX+GFKIDxFMFa+68QxbtaklT90/So/e+hKP8sd/lbHKtF9c+u/g3jKw4DS/3uzfQed72bdb7iw5PMpq0XguqUIXD/44CF+KGLWtlzEULfEGHiF3H/WfW523K2ZY+9TfcI1d3BcbnuUeE5Q58wB61uPyXOeKdPcqH7wbDermbNH6Xt06eVnKYHwFMFY+SlN/XvJGvMvuC3Mv+D2vk78jpA5/K2OVaL75vabi/zThuatFu8bQTUtb3Gze2x8yedTOonney2yDW72mvn4vopZr7a1brYlxsAr5egX/WvOpCEvZd1jmfkfuHMHx5W77Fw3ayO+1lK0UA9/EPiLJwAAAAAAABSCG08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQ4bj45+xRN9vHnnGzJ2yAm7UxH1xcYT6eKePir4hoYv+EkKyKhgctXOKreSokrqi4ePTYzk397LUfDo0dHBSNlSvR+PnSe3xIXInGu8sh97Uc8OZgN5sy5onST1ghKratqJh6biou3kt8mx/VYMsLjk9+xMfFldmicd0r+DKZMy6+43F+1v13fjb/AD97V7wEbS9eI2VcfL0omi4SUdcoFQ0P2vW6a91MhcQVFRePHtvEfxvYOPHvkSIaK1ei8fPrd4ydLxrvLofc1zL9/NvcbOrIv5R+wkoRsW1JxdQzU3HxTuLf7SuXbXnB8alntY8t7C8C2zOCr5O54uIr3najYQ/2cbOHvyaeSrFY/PDVUfyQpuLiK0RkuXVjcYFBKhoedEpjH71WIXFFxcWjx643/yax23XzQ8dGRWPlSjR+vtc5u4XOF413l0Pua7nexrrZsXZKyeerFBXbVlRMPTcVF7cl/qEekw7O+4SfcgTHf2KjQ+tmm3/4Qy8bGDqWuDgAAAAAAAAqihtPAAAAAAAAKAQ3ngAAAAAAAFAIbjwBAAAAAACgEKLc/XELO5e8iQqJS7Wb/EyFxFUgXB0bDYlve50btX/pKjeb0Tp2OhUDX5gQIVfHdjr52dC+igqTRwPhap2y6DEfEm8YOjIt3p07Bp77WuxNPxp8Qd7Q6AtjxrlZM9su6x7liIZXExUSD38OGolv6g3BF4QtTEpYVIXEpbmqmi5mKhCujg2GxI852kfDD3sldKikYuAW7N+rCLk69oQbgvsKKkweDYSrdcoVX/SzaJQ7Jd6dOwae+1qUQeP3LH0TYerXRby0bcuse5QjGl5VVEw9+DloaCvdbKMFf+jckrR8Xwx9XFxSIXFlk4ioq0+lCoSrY4Mh8b6Nf+xmp1vpsXoVA7eECLk6dvpFNcF9PRUmjwbC1Trl8HN88Dn6cpsS784dA899LcoDdnfJeyj/ace7WUsTkf4E5YiGVxMVEo9+DhqYfyjBJluSfE3lxl88AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIUIx8VrbWrWjf8+6QM/7P7fsYM3iphWQx/dCh+70Y9eXS6OFXFCFQNXhwab3LZOrJOx8nWxmRINk6sIuQqOr3lUhMTVpzlz+FvJfb4UlbqWL1xwsh+KwPWU6+4qw9UUr5e4fT5bPGugHLfZp4x5ws1yx+PfE43KXmLdzFmx8+2bcjH/Ym2swRp27Ll+dkyjMaFjdxaN2DeC16eO3VmsG3tW7HwqBr7eNyatSavY+dQ6FQ1fL4LjaqZEw+QqQq6C4z9r7GctW/hZ7vC3kvt8KSp1LYMe+ZKbqcD1lJGvl+NyCtdJfF0tqtDnfurIv7hZ7ni8tRIvYjPEC+DQvA8e+bcWDsh6usN6/M3NptutoWN3Ej/Pvhn8YV0du5NY94r5H8zFy56Mgat1a/79pf3/1vkXYRUNbyKeSqFmSjRMriLkKji++JwVbqai3LnD30ru86Wo1LU8av53OxW4PsaOKMflFG/Jaj/bNvMDN4KOtVPcLHc8/n3z32+9bKCbrbRpofOp1yuFv3gCAAAAAABAIbjxBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAAChGOi681H5GL+vtDPv6XFAOPikbIp/pwW89Fp7rZwh4+xtxWxCmjkW8VJleWq2EwVh6lrllFyOce4tcte9/PmotYuYptlyM4XinqY/udTXGzA94cXPzFbPBfMCp6reLYUXXmvxkamyi6ZiZD4pmlfBwyOH6Rf31R/0ZKz/Crttevb+nHRjRrUvqx37zGz6YnxMCjohHyh/o3d7MrpnZ1s137vu1mrfr580Uj3ypMLtclxMqj1DWrCPlXjvez5uprT3zvqth2OYLjlaI+ts/ccIabTT//tsKvZaN4goqKXqs4dtjUpX42qEPp5wsqS0g84eNQn9PB432IW/0bSasSnvTwrPih7huln+6fbCo9lHxwD//a+qaJBxUJKgYeFY2Qz7r9MH9w3Qw3OuXsL7vZchGQjka+VZhcr/OisfIodc0qQt5N/GbTIfh7ZqWC45WiPrYL7Rw3u97GFn4tm8z/Xq2i1yqOHbXM/A8z7U38EJVbGULiKR+H+pw+ZL9wM/VvpGxnbUq+ltYiQh7FXzwBAAAAAACgENx4AgAAAAAAQCG48QQAAAAAAIBCcOMJAAAAAAAAhUjI1Gp//+m82EIV/v6JCMENDEbIpwVLrcFYn9J53SVutrDTtW4mutq2TnSD1UwFx2XAXOyRW2fxgTQVnfhlwZ5m7pD4Ac/7+6a/29vXaqP7qnVR6nxbYiQ9JThejpD41mLKdf5BBepzP1vEujf4Vmk4Gv5elbQ2d70wtk6Fv6+o7+OH282PnbDex8CtXqxT/d654ljh1XU+Qr5r07WhY1UMXM1UcDwlYJ6iSRc/+/UdfjZMRMiV3CHxuedf72bdb/BfgNF91boodb4tMZKeFBwvQ0h8azFl5Otupj731l+8UHYWL2IqGq60SnhaQyZ9e/zBzVTkW82euLmdP2HjXrGN68Sxklo3O3Tk3be/5WannL6DWOmj4WtsbmimguM6YF689ebfJJaLOHEHmx46X+6Q+Cgb6Wb32PiS91XrotT5tsRIekpwvCwh8a3EMXaEm6nP/Wyb5mbqnTgaDX/fVriZegVT+IsnAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBA19fX1Kq3q3HGIn60RffDzupxa+tWokPg0sUnu86Wsm7vejda1HuNmC9epC4xRwXEVJi/H+VRw/K8H+tnahH+2qC6d/WxBQnU9d/w86oA3Bxe/SWbR4LhSaz5w2sy2S7kc56ebfMzyqAZlKC0nqDNfGd630S5u9qxvkNtbM/2sl39pkvYVx5ZqifkXg1p7wc26vyOK1FHdRUg8GP5OOl/Kuk0+1vvFz/oXq5QYuAqOqzB5Oc6nguNDv+xnb/fevGsqRW9Rz30n+L2h5I6fR00//7biN8ksHBxXlq/2s7YtSz+f8Jur73Ozr1w2POse2U1d6kZnTrnIzX54lAg1D/XBaZvxYWzfW4+Jrft3LhY/pO0lXmzerit9DxUSr5td/PlS1tX6B+Mcda5/akiThAckqeC4CpOX43wqON7L3nCzw2zoZl/X5upqa9xsvrUo+Xy54+dR19vYwvfILRocV1abf2hLS/MPd0lx9BMD3GzSwf6hD9Vkmfkf6tuKr7+htrubfWCvulmvYHC8hYmnjwn8xRMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQiHBc/Jb9/Wz0NlUUEo/uERW9lrY+Gmvdfc30vUXXupkKdUcj5LmD47mvZVq30q8lasD2fvaabwRK5YiGp8TKP23BcSUlQn79Bh+uvLBRdcfFj9ohIbgd9MAv/GzPz+Y7/3zbyc26vpMQic0dEo/uERW9lm4ilNzAR32HdPCxTBXqjkbIcwfHc1/LgOtLv5aonv5lxGbsFju2HNHwlFj5py44riREyO8/y8fFj7u1uuPiI8aPL3yPfW2ym51o4peAUtwuflCtE/+GUblD4tE9oqLXsrGnnzXzfxtw3Ok+nqxC3dEIee7geO5rGW4Jv8cF1ZiPRXc0/2AXpRzR8JRY+actOK6kRMiPfnEPN5s05KXkayrSW/Z44Xtceo9/j7BRJ4WO5S+eAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgEL4YtnHqFPR5i+UId4dlTtWHj2fCInLdSIgulzMRONbRsNTQuIqBr4w4XzKwHl+lhIcf0B0zEae7mcN3y19j5QYuFKOgHk1GXzBwW72lzETQsc2Ft8gKiReZ7Hqbq8GxYe6t0THHuFnM9/Kd/6NJl5IuoeeX6GVIySeskf0fCIkHt131czYpahoeEpIXMXAg43YsNcv9LOU4PiDm3xJ/NCd/etI64RoeEoMXClHwLyaDBq/p5tNHRx8MMWgDn6mQuJTl4ZO10n8W8Jssu3rZifmOnmP9X42p1fp5ytHSDxlj+j5REg8um9b+atcLBqeEhJXMXBLOJ9yn/nfp1KC4+fakW72Y7vVzVKi4SkxcKUcAfNq8oDd7WYH2ZdCx7a3fm6mQuLLLPiD1ZIBsXWfMteM8u8RlwaP5S+eAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgEKE4+LnHXiJH0bj3SrUnZu4lvPeu8vNbux5aul7JHwcPTv5z997c64t/VqCVEhchcnVurYierpwXWxfdewOYt+3fGdUOtZ3zKzRcj9bGwyE5w6JQ9vzghFupoLjKhquguPKI5t22dzL+kQp1xI9376N/DU3qvHHdu1V8rYV0b1OhMSj8W4V6s5NXMusdm+7Wd8P+pS+R8LH8eJSH8Ec0mFt6dcSpELiKkyu1rXyLU9bH4yQq2Nn/pef9bs8dr5vNfAh8Zbix5TV4v1ORb5zh8ShDZriH0whg+MqGq6C48Jvfvn45l7WJ0u4luj5zpxykZu1M/969YElvF6V2xz/PRqOd6tQd27iWo4/fmc3mzjxjdL3SPg47r/dvx8cd7p/38hNhcRVmFytWy7i2E2CT6pQx95pr7rZt23X0Pl+YD9zMxXv7hoMhOcOiUN70n7vZio4rqLhKjiunPbEYZt/YZ8g5Vqi52srvtZ62EI3m2PqaTGVwV88AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIWoqa+vFzVY4fBglFsFuKMR8qho5Fvtm3Js9Fqi+4p1Cz/ng+hKSvhbhcSjomHy6DplWjc/yx0Dr/a4+AFvDq70JRSnkf9CmHJd7Os+6qebuoTWHdXABy5VDFxJCY7v1chfX4/eft2cd/xMrUsx862MJ3snGLlVAe5ohDwqGvlW+6YcG72W6L5i3a59fWBYSQl/q5B4VDRMHl2nDLjez3LHwKs9Lj79/NsqfQmFaWgr3WzKyNez7vGbq+8LrfvKZcP9UMXFlYTg+JnjT3GzVbaXm7WyP4XWpZhgDfOc6OZglFsFuKMR8qho5Fvtm3Js9Fqi+4p1p5z95dC2KeFvFRKPiobJo+uU4eZ/78odA6/2uPj1NrbSl1CYBuLf9xg7IuseRz8xILRu0sH+vUnFwJWU4Pgm86+ns2yWm/W1vqF1KS61HULr+IsnAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDxuPgll/hZ7mh4itxR83JEyIPeW3Stm0WD3ko08q0C5stFbDVl3+ixb+0ZW6dUe0g8aqsOjgtTxjxR8rHXb4jFxc9p9Dc3S4mGK0ftELuWqMNFkPmxhBB01rh4XXM/yx0NT5E7al6OCHnQkA5r3Swa9FaikW8VMF8lmpop+0aP7Xd5bJ1S7SHxqK05OK5MHfmXko+9/6xYXPy4UYf4YUI0XBkxfnzW81223MfAr267seTzZYuL317rZ7mj4SlyR83LESEPOu50//4cDXor0ci3Cpi3FVHulH2jx37bdg2tU6o9JB61NQfHlWPNP6Qh6ugX9witu2PIQ26WEg1X3rLHs57vzWU+Br5T+9J/ISAuDgAAAAAAgIrixhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAACiEL6V9HBXRjga9qyjUnaQcMXWxR0871c3es7uKvxZBxcA7N/Wzhetix0YNnOdn07rFjk0JiaswecoeW0vovBwGX3Cwm6ngeJ2pArAPevcSt9mjIXG1R0qEfM47ftajd+zYlJB44VREOxr0rqJQd5JyxNTFHi9O7eNmQwa9Xfy1CCoG3kQ09tcviB0b9fqFfjbg+tixKSFxFSZP2WNrCZ2Xw6Dx/skjMjg+dWnofJ3Uv2U0JK72SIiQt7I/udkq2yt0bEpIvFAqoh0NeldRqDtJOWLqYo/7b57hZsed3b/4axFUDHy9+Lmtifk3iWhIXLnP/O9Ywy32O2VKSFyFyVP22FpC5+XwgN3tZio4vszEU1FMxMWXrHajaEhc7ZESIe9rfd1sls0KHZsSEk/BXzwBAAAAAACgENx4AgAAAAAAQCG48QQAAAAAAIBCcOMJAAAAAAAAhaipr6+vD6285JLYGVMi5EruY5WUaHhK6Dzz5+C9Tj443rRHbIuU8Lei9lV7RNcpKmqu/HHX2Lpq98BkP/tR58Hlv5Aqo4LjP90kSsbCUQ1E3TjB4EHfDa3rt+bHWfdVogHzmTn7gnXNY+tSIuRK7mOVlGh4Sug88+dABcfXB4P1KeFvRe2r9oiuU1TUXOnpW6NbpAc3ve9mgy58tPwXUmVUcPw3V98XOvYrlw3PezGzHg0tGzH5g7z7CtGA+QRrmGfD22tj61Ii5EruY5WUaHhK6Dzz50AFx1tY99AWKeFvRe2r9oiuU1TUXBkdXFftzrUj3ayL7V2BK6kuKjh+9BMDQsdOOvj1rNdy9LDYL7jXPPyLrPsq0YD5pbZD6Hz8xRMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAACgEN54AAAAAAABQiHBcvHdjP3vn4FNju0QD4bnXRVVT/Dzh41j4OR8Xjwa9lZQYeEo0XImGxBeui617a8/SryXFxo5+1nBx7NirRl0fWjf5gp9sxhV9uqkweZ0tCh2756AbS963HHFxZf5sP1tdl+/833nez27q1id2cDQQnntdVDXFzxM+jl37+rh4NOitpMTAU6Lh8lqC7df1wWcK9Lu89GtJsbqTn7WMvSzZg9Nj/33ioPG3bsYVfbqpMLlNXRo7uNVzJe9bjri40s78a8TYYDj23znP7nSzG28eGjs4GgjPvS6qmuLnCR/HKWd/2c2iQW8lJQaeEg1XoiHxJhZ7k/i27VrytaToZI3cbJFtCB07ykaG1u1ou2zWNX2aqTD5MpsZOva0YZ8ped9yxMWVHuZ/UDvBYj+X8hdPAAAAAAAAKAQ3ngAAAAAAAFAIbjwBAAAAAACgENx4AgAAAAAAQCF8nawI0UC4Uu1hciX6sWX+ODrfcmlo24Wjr3GzlAi5khISV5aL2Gp0j5SPo1JUhPw984X/nubL0PuO+WZoj+j5cqum+PngQd8NrXvzHh8DV7OdRp0QOt/MFn5d0+n+fD16h05njZb72TuiSxs9X0VEA+FKtYfJlejHlvnjeNV+5deJyPeudqibpUTIlZSQuLJKtDyje6R8HJWiIuQ2dxs/W/GhG00deVZskzax8+VWVfHzWY/6mfi6Only7IeScfvGfiiZsG87Nztzsn8gxirbK3S+JWsbulmv5n90sw+C5yu7aCBcqfYwuRL92DJ/HHffd2Fo21OG+wfepETIlZSQuNJW/Nob3SPl46gUFSGfaW3crJ+tcLO/299Ce0TPl1s1xc+PHqaegOWj4aOXvuZn4shbOnw2tO+lw45ws4cenupms2xW6Hytba2bdbKBJZ9P4S+eAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgELU1NfX14dWXnKJn6VEwysVCM8dF0/5eFM8L+rTq0QdtVWsrPre131wPBrljobJU9Zts7efffj8v7+2zfHWnqUfq2LgDReXfr6rRvloo6IC4YqKhqfExVP2nfylb4WOlZqK77fuL8eOfXFc6fsKL1zlY+DKF66IBcfDdmjmRv2m3hU69NuP+9n5O6Ze0D+oa+5nKdHwSgXCc8fFUz7eFH1v9rMB4i3/9ZrQ6Ya08sHxaJQ7GiZPWVe3yc8aZ/6v1/pdXvqxKgbeclHp53twevCDU4FwRUXDU+LiCfue9Lz4xwxavN7PfvktFX8V/hhcF3RhMDh+fTA4Hia+Z0ZMEU+cEF60v7rZ6zYk9Yr+1+21fpYSDa9UIDx3XDzl403R7rd+1lw8d2rthtDpjhvuU8nRKHc0TJ6yrqN93s0W259D1xf1bdu15GNVDHyRxT73yigbGVqnAuGKioanxMVT9v2WxX7uVebZUje78/6jYgc/9kzJ+yqXdvLBceWaRbHgeFjT9/we98d+n/rh0B3cbMGzsW35iycAAAAAAAAUghtPAAAAAAAAKAQ3ngAAAAAAAFAIbjwBAAAAAACgEPnj4ko0wB2VO9SdOzgedOhjpR/7q23zxsWjxy4c7SPkKgbeWfQ5F66LXYpSu7ufNQs2pZVo1DwlOJ7ipFFjQ+uiMfCoaHBcrZv9pWFZryVJ39l+ljkurrx5Tyw4vtOozMFxod+a2LXMfCvjptG4uBINcEflDnXnDo4Hnbr07ZKPvWvIr/wwIS4ePXZX8xFyFQNv0sXP1i+IXYrStJufrZtX+vmiUfOU4HiKB+e2jS2MxsCjosFxse6kJ5bnvZYE9x4pflDJHBdXTg4Gx8flDo4LIybHguMTrGGeDaNxcSUa4I7KHerOHRyPevmJ0o/9snhQTEJcPHrsKcP9A3RUDHy9+TeJJlb6m8R2tr2bvW/vlny+aNQ8JTieYqh9J7QuGgOPigbH1brLLfZwpXK4/P4D/DBzXFwZvTQWHL+lQ+bguHDNw78IrbvUfHBc4S+eAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgELkj4tHQ+JliHdXKhqufPOxjW62OhhrbGn+2J80E4XTzCHxFO993UfIU6gYeJSKhkf1/NaZbva7WT8s+XwbRRO+4WI/u2qUj+ulhL+j677Z9xw3+/72PpS4RfqUBcebTvf79ujt11UkLh4NiZch3l2paLhyzjQfEl8TvJQW4sMYe9Ev/TBzSDzFkFY+Qp5CxcCjVDQ86sVO/d3sW9+eUfL5Vnfys5aL/OzB6eK/O0wIf4fXrXvdjU76nSi7b4E+bcHxMyf7MPUq28vNyh4Xj4bEyxHvrlQ0XHnqWT/rsCZ27NIWfjZSvM5nDomnOG746KznUzHwKBUNj7r/Zf+5/+PuO5Z8vk7mP/eLzH/uR9lIN0sJf0fX3ShmK0zcR9gCfdqC4w89PNXNZtksNyMuDgAAAAAAgIrixhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAAChEPC5++Kl+Fo135w6Op0TDyxAcP/SxrKeTftVMRHzLERdPOHbhaB8cbyuCqQvX+ZmKi6touFoX3UPpLJqiTY/8th8e/t3YCTPrteJRN1MhcWXyl77lZgd38eHFJxaIoGJT8X20rjLh/iRLi48dqrj40vdjx37hirzB8X5r/LVkjYu/08fPovHu3MHxlGh4GYLjpy71IfHc7hpRobh4wrG7mg+Ot+rn161fIGbi7UlFw9W66B5Kky5+9sca/xo5fnFlXiNHDm/rhyokLpz0/CY3m/53f77P7LjczRav9+fr2CS0bVW5t4sIrGem4uLtgsdenzk4PmLyB26WLS5+8xt+Fo135w6Op0TDyxEcf9mH37M7TPwuUY64eMKxpwz3D99ZLsLaTcy/gKu4uIqGq3XRPZT15t8kfvqiP9/fh4wJnS+3Aebfn1VIXPmW3eVml/f9vZv916wvudk8W+pm3axDaN9qcvmw4v+OR8XFt/HPNZGuWZQ3OH7Nw79wM+LiAAAAAAAAqChuPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgELE4+KXBCO8mUPdSQHzzMHxckTDW9pGN1stoo5rzVc7n952u9gm0UB4SoQ8YV8VIVch8dyiAfOeP32n+ItJ0dRH/cIGfjXfdVRSNH6eEhevEenXeh9lVV64yke+o4Zd6GfNxav4zBY+TN50ut93XbAFGlInQqVK5lB3UsA8c3C8HNHwFuJS1ogP464vi9jjih/ENokGwlMi5An7qgi5ConnFg2YnybaudVkZI146kbQSR06ZbySyonGz5Pi4q3FJitjhfULRXA8aq6tdrMm1tnNJuzr38fOnOyj1j+0L5Z8Lf/k9trYutyh7pSAee7geDmi4Utb+FmHNW70uYn/7WZ/ve/S2B7RQHhKhDxhXxUhVyHx3KIB88ttfOHXkmLHZ2IPKlL+a/+b8l1IBUXj50lx8dbi/WCleN8QLu3kg+NR29zmHzCy3hb6PYYd4WYPPTzVzY4J/i0Tf/EEAAAAAACAQnDjCQAAAAAAAIXgxhMAAAAAAAAKwY0nAAAAAAAAFCJ//lKFupVoIDwqc0i8HKIhcaWBdSx942hIXImuS9i38y3BsKHw3td9mFyJhsTVOmmjCLU23AIDrNOeynu+SsXKVUhc6XaXn807NXZsMCSufOEKH/5+855YcPz9AeeG1vUTIfEevUOHFk+FupVoIDwqc0i8HKIhcaltaz9bETw2GhJXousS9n319V/5dcG3pyGtfJhciYbE1Tpr60cb2vhZo+i/RxW5d2npYXKlUrFyFRJXDnvPR1h/2TMYHA+GxJXr9/U/gJwcDI4fZ31C60ZM9hHyqqBC3Uo0EB6VOyReDsGQuNLAEiLf0ZC4El2XsO/d94knsQQdN3x0aF00JK7WKYttgZt1tC6hY6vJ5c98J+v5KhUrVyFx5dsPL3azO4cFf08PhsSVaxZ91s1GL40Fx0+cod7D/Oyah38h1s0SM/EwG4G/eAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIXgxhMAAAAAAAAKEa+7pQS4KxUcD/rmY6VHvpVoNDxlj03mQ2Zm2/lRSgxcUYHwqJSAeTBM3vMRESYX6xZefJObqZD49pef52a1N57pF26JIfFyyB0rjzrg6Ni65YVexWbZaVTpwXFl5mdiEXKzH5S8h5MS4K5UcDzonGlvu1k48i1Eo+Epe9jylbF1KTFwRQXCo1IC5sEw+YsqTC7W7fr6YW6mQuKvThjnZuO/c7KbbYkh8XLIHSuPOnZ58H07+pCRMhiXEBxXjrKWwZW1Je/xT1IC3JUKjkc99ayfBSPfUjQanrDHpuivgSkxcEUFwqNSAubBMPn9990SWnfKsO+5mQqJPz3Nf4+eP3CCm22JIfFyyB0rj7p5/xtC6+YVfB2b45YOpQfHlW/NiEXDrX9sGX/xBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAACsGNJwAAAAAAABSipr6+PlbyPPxUP0uJhkepPcT5Dn2s9C3KQQXHlWhw/Mv73xVaN/pP3wmtyy4YAw8fq0Qj5GqLlde42YZ1PiSeoraNiJBvFBHV3GHyzzXOez6lUtHwKBVVbts6dmzXyXmvJcELV/m4eAfxDIGdro2GxIVHM8bF3+njZynR8Ci1hzjfqUt9ILyaqOC4Eg2OX7n7TqF1XQeOiZ0wt2AMPHysEo2QCw03ne9mG+/zsxQqQr6hjV+XO0z+x+Kf21KxaHhYU/G1sS74dfWZfnmvJcGFIi7+gVj3detc8h4H54qL3/yGn6VEw6PUHup8Lz9R+h7loILjSjQ4vtew2LqOv46tyy0YAw8fq0Qj5MLjt57lZmeNj20b9T8iQr7YFrhZ7jD5T6z4nwsqFQ0PayG+NtYEv65uF4H6Crm0k4+LfzjDr/vObaXv0Ye4OAAAAAAAACqJG08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQaXFxJRgDj/rmYz7KHQ1w56YC4SnXos73k8Nj57vlQx8XH72N/zeS61RwPBr0joqGxJWUaxH7tqq/KXRo7rh4lIyQ51aO4HiKaoqVqzB5/1fKfx0f4817fHC8quPiSjAGHnXONB8Njwa4c1OB8JRrUecbOzD2eZ7f23/fd32nLrZOBcejQe+oaEhcSbkWsW/DaReEDs0dF49SEfLcyhEcT1FVsXIVJu+3Q/mv42OcLILjVRsXV6Ix8KinnvWzaIA7NxUIT7kWdb6vDo0d27iln9Wtjq1TwfFo0DsqGhJXUq5F7Pv49aeFDs0dF49SEfLcyhEcT1FVsXIVJr8h48/biUYv9cFx4uIAAAAAAADYYnHjCQAAAAAAAIXgxhMAAAAAAAAKwY0nAAAAAAAAFCIeF7/kktgZRUj80Mf8smioO2WdUo4w+a8Oz3u+Gxf4QHhj0RKMxsWV0a9dutnX9ZFoSFxFw9Wx0bh48NhWK69xs0qFxKOyB8c3iijrkO3y7pGiabBquy74oIJorPxzj8bWzWsfW1chZ0z9rJvNr/UR5K7NbnCz2yxjLLmueWydCImfutQHwqOh7pR1SjnC5Hd1CIbYg+b2XuBmDa2Hm0Xj4krXBv7rJywaElfRcHVsNC4ePLbhJv99UKmQeFTu4PiGNn725yp6LsXi9bF1HZvE1oVj5Y8+HFt3yUWxdRXSo/5DN9vr637dnx7xsznPZbqI24ORchUSf/kJP4uGulPWKeUIk+9+cN7z9RHfzHPEN0s0Lq60enSzLumfREPiKhqujo3GxYPHPn7rWW5WqZB4VO7g+GLzP2f8ziZm3SPFPFsaWtfNOoTWhWPl28yLrfv+z2LrKuU2/yClzm/574+FO4jvrU6x71/+4gkAAAAAAACF4MYTAAAAAAAACsGNJwAAAAAAABSCG08AAAAAAAAoRDgu/k0b5WblCHVH5Y6VR/f4ySWd/UIRWI+KxsDfXe5n27cteVtp9J++E1uYEgOPni94bKv6m0Lrqj0uHpUSIW+24of+fAdd6xdGg96Kioa3FVG6hctLP5+w35PnuNkfDpob2yOqioLjKi4elTMufs6v/Kwcoe6o3LHy6B5j/0NE10VgPSoaA1+z/E03a9F2p5L3VboOHBNbmBIDj54veGzDaReE1lV7XDwqJUI+8qZxbnbY5f580aC3oqLhref4WbO+pZ9P+UZvPzvhv26LHRxVRcFxFRePyhYXP/hOPytHqDsqd6w8usfXRviZCqxHBWPgx3/TP4Bi4k/EN1+Kjr+OrUuJgUfPFzz28etPC62r9rh4VEqE/D+mjXCzrwzs4mbRoLeiouH9zT8NY4XF/n2jEfLn/sPH1Kdfe1Xo2LBqCo6LuHgYcXEAAAAAAABUEjeeAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIcJx8UPt1NAJU+LduX0o7qt1tTo3i17fry4RYeNoSHygP/aOV33IuS6hsaiObSwahlGXPLObm61stl/s4IRAeAoVF99aQuLVZNmIG92sxcLa2MEpsXKlaSwCWRZlCI5fOXUPN1tswc+9kDMufurLsXUp8e7cviA606+e62fR67trUEJIvLv/xCxpvLeb1Vnp8deN4tiG5gOzUV33e9IPl06LHZwQCE+h4uJbS0i8mnz1VB8hbyqex6KkxMqVe+/NHA1PUYbg+IAVPiS+cpvSz5ctLn78D2LrUuLdmX1msX9DmN4gIZK+1+l+Fg2JN+7lZ3NX+Fn3YGFf6SGOnVP6N2S/Kf79ZeZX74gdnBAIT6Hi4ltLSLyaHCl+N/4wGAhPiZUrl//1qKznS1KG4PjK+19zs9YrgoF/hbg4AAAAAAAAKokbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDZC22VCokr29gmN1PXJ6PhylwR3BNhNElEyKMh8dHb+LD7LR/eFTs4KHotrWv/EFq3cttD/LBCwXFsjkVuUiuDmd+InW7QzX7WVHzPqOB4NUXDhX2ePMXNnhuUNwh4xtTPimnpIfFqUamQuPLCmNg6GQ1XNomCr4iGSyJCXtc7FhLv+o5/cMb83o1j+wapMLnUYWDwjNP9qELBcWyGpeLf46DD3OipVxaETnfSYV3cbLH4cUsFx6sqGq50zRvAVXrU+5C4JYTEq0KFQuLK9I7BILqKhiu1/ncTa9YrdqyKkHdvGTx2tZ81Dh4bpcLkU/yo31Oxn+9m7nurH1YoOI64deJtfOZn9hcr1c+43hjxu8k8W+pmKjheVdFwof1vr3CzZZY5Ln7bK26UFBJPwF88AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIWoqa+vD1U7DzUfuG5pG90sGhdPOTbqV4eLYY+mfrbmQD+b/4yfDW6RfE3/6JYXrnWzcoTEFRUXbyw+3Eue2c3NGojA2yYReFvZbD9/wtb+fEk+9LG0Zuted7Pm4prXimuOrlNSji2PaEhcWCFipm0SaqZDRJCyiqiQeFRKcFzHxfO6zc7Pdq5TX/azFqKrHY2LpxwbdVeHPn7Y+D0/u+p9Pzv9CD8bKEqqCeY37ulm5QiJKyou3tB6uFnX/Z70B78rgrrbi//ua+k0P+sk3rdTbGjmZ6P+7mfvzvWz7buXvk5JObYcgiFxaVvxPrEk4X1i8c9LP7YcUkLiJ5b+HiPj4pnNeS7TiY4Xoe6l4ofN6M8iKcdG7X6wnzUVr11t7/ezp47zsz38a3qaxX5UjpC4ouLic/xTAfpN8e8ba8WvgM39r4o286t3+OFrmWv63Xwdu88Xv+VmLQb677010/y1RNcpKceWQzwk7tWaf29vZutKvpbmfx1a8rHloELiUcv+mvD7gIiLZ9cpFvjnL54AAAAAAABQCG48AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQsRKUJshGg2PhsRlIHzgtn42bUnofNamtZ/NedzPuh7iZ3NFXbG7j+apa1GBcBUSlx/bC34U9arv99muvg8rQ+LKykGt/LDuD362+BE/E5/63FRIXIlGvlNi4JULiftouFknNwmHxJWUkHiVSwqJH3SNH3YU39OL/ezKqXuIM9aWfC3VLBoNj4bEZSC8u9hkbvCEnUWo++K9/Ox7f/KzH+7iZw1E/FdciwqEq5C4/Nis9GBt1+t/76/lwi+5mQqJK3P+cJGb9djtJL+woQhX5w6JKyokrkQj3ykx8EqFxFU0vIOoxEZD4kpKSLzapYTEF7/mZ2df7Wc3X+ZGA9SDPbbGT3M0Gh79OUYFwhv38rO62bHzbRKvtwt9fNq+KoLj8/6fnzUTfwegrkUFwmVIvJefqQh50JDT73OzF28f7heKkLgy7Z5b3KzZ5C/6hU/P8LPcIXFBhcSVaOQ7JQZeqZC4ioY3FW8b0ZC4khISr3ZJIfEvj/XDzxztZ/dNcqOV9/v3l9Yr/IO3KoW/eAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIXgxhMAAAAAAAAKEY6LR6PhSjgQLoloeDQkrkSPXeqD4x2a+ED4UrvEzVRIPLeu/p/D5ot/DhUST1InQuxKx6/7mQqOl0FzW+pmlQt/p4hFw3VIfHbui8nrRRGkHLK+8G3TQuJ3u9l+Hf8SOvYPjQ91s8VbQUg8Gg1XwoFwRUXDoyHx6PkUERy/57aFbjbqrOZupkLiuXXd0Nbv22i5n4mQeAoZElc+EPXSdqJeWg7vzvWzSoW/U0Sj4Wo2xL8uVZWO3/CzxT8vft+UkPh8//PHPfs8EDp0VKfr3GzlxaVfSlWIRsOVaCBcUaHuaEg8ej5FBMd7nHunm825/XR/rAqJ57ZmmZ+1aO9GMiSeQIbElQP7+5kKjpdBi4E+7F+p8HeKaDRczaYPLD0kXg5rP/esmzX/69DC900LiV/lZi+uEedr60dD/p+/GVBNIXGFv3gCAAAAAABAIbjxBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAAChGOi//kks6xhdF4d0ogXFGx8pQ9xPlUSDxq9Danhtbd+PS1btZYtBhVSFzZXjSb3y2+2ayp4HhTEdncJGJ980qPpamQeHUFx0uPhitVFRJf4WOM1qYMMca6TW60z9PfLvl0KiSutF7hZyvb+Nl+HX/lZgPths29rKoz9j98RFuKxrtTAuGKipWn7CHOp0LiUV3fqQutm9vY10HVW4IKict9awf6Y5tNCx2bnQqOi+637brWzzY2K31fFRKvpuB4SjRcqaaQ+LbifWJJGd4nFomfKwYFf9ZVREhcWSb+O9/25t+z7tlnjDj6jM29qurytRGxddF4d0ogXFGx8pQ9xPlkSDyqbnVsXY/3/WyO+IVAhMSlv6/ysx1bxY7NTQXHX37Xz6aL94jtSn94hQqJV1NwPCUarlRTSLzWmrpZM1tX/Mbz/E9W7adfWvLpVEhcaSxuX9SJ2xwqQj7EiIsDAAAAAADgU4gbTwAAAAAAACgEN54AAAAAAABQCG48AQAAAAAAoBDhuPh5T/jotXJjz1hEO0xFw5VoSDx6vuAet3x4l5u9u1wc+4VY/FyFxFNULCQeNUAETn/zu8K3zR8SjwbCleg6r6pC4ko5QuKCCokvszVu1t78N9xznz/Ln7DdZ0P7/nr1IW62X5vHQ8duDWatFjFPoa/1ybuxioYr0ZB49HzBPeb3buxma5a/KfYVYXJxvobWo6RL+zgVC4lHdWvnZxtLj8SG5Q6JRwPhSnSdUk0hcaUcIXFFhcSX1/pZWxGtX/eYOOEXQ9te+POz3eyeb4wNHbulO/6bsQcATJyYeWMVDVeiIfHo+aJ7NG7pRsd/07/OT/yJ+PpU51Mh8RSVColHtRA/R283u/Btc4fEo4FwJbpOqaaQuFKWkLigQuItxK/za8QtiGWHfdcP18Z+yd+t0TVu9qKVHjWvJvzFEwAAAAAAAArBjScAAAAAAAAUghtPAAAAAAAAKAQ3ngAAAAAAAFCIcFz8xoMv8cNo0FtRke/o+dS66PkSrlmFxBUZWJ/rA3TR8yl1vpOcPUw+euGHeU+oNBIx4o5f97PFjxR/LVI0Gl56IDyq6kPiKVaIr7UX/ffMPktHlLyFDIkfdLdfGAyJS40/3ffy+7aMxbHDVOQ7ej61Lnq+hGtWIXGl7wcisN7cf81Hz6dstDluljtM3nW3W7OeT+rW0M8WiApruzIEx5VoNDwlEB5V7SHxFNuqn0m+4UeNni19DxUSn79ULIyFxKVO4R+9tzrhOHaUinxHzycj38HzJV2zD4krEye+4Ycbe/pZs9j5pB7iCUS5w+QDF+c9n/K+eJ/8cn8/e3pG8dciRKPhKYHwqGoPiaeotaZ++Lln3Wi77z9T8h4yJP7lq/wwGBKXum0UG5d+umry6f4tCQAAAAAAAIXhxhMAAAAAAAAKwY0nAAAAAAAAFIIbTwAAAAAAAChEvHAYjXKnRL5zB8fLYPQ2IiQu3PLeGDfrKtph80VDVYmGxMsRIU/ySuvYuuzB8dKj4c3Nh0bXWoeEa/G26pC40sZHlVMsExW+6Qfd7xfWbfKzzIHwPyz2sd/9Ov7KzZrtdr6b1b5yQ9ZrKVw0yp0S+c4dHC+Dru/UhdbNb9rWH7vBz+Y3Wh46XzQkXo4IeZKGH8TWfZA5OJ4SDX93rp9t3730a1G25pC4siTv+4QtF6Hrtav9bNEGP8scCB/13AVuds8+/ufGZ/1bhw3dkr4MolHulMh37uB4OdSJrzulmQjdr1nmZy3ax84XDYmXI0Ke4jPqWsTPdwfmDY6nRMNbDPQPS1gzLe9r3NYcElea2bqs52shbkvM3U2ExOeJX95VIDzBkBZ+3xfXXOFmxzfy700TN4j3sArhL54AAAAAAABQCG48AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQqTVEVNi4NFjU/ZIcMuHd7lZNCQur/kFP1qiutoiBp5ChcRVcHyL1H0fNzpjtxPc7LZXzhMHq5B48dqZDyPP75A3hhe2wocNc0e+K2X6gfdV+hI+kQyOf9ZXYweW42KKlhIDjx6bskeC+b0bu1k0JC6v2Xq6yZJG6o1ieWyPIBUSV8HxLVK/pm4067CT3azvL8f5Y1U0vBxeF7HgEf69rSy2Fe8TuSPflfLuikpfwSfSwfFeYuXsoi+lWCkx8OixKXukaNxS7BsMiatrtsV+tKyLn7UQMfAUKiSuguNboja7utHe3fyy5+f5mYqGl0N7cX3Pta9MSLzW/Hts7sh3pcztfI0fqpB4hajg+FnTfXDcdirDxQTxF08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQ3HgCAAAAAABAIdLi4ikx8GggvAwh8dxueeHa0LqUyLc6VoXEq95+zf3sD+KD2yA+uHX+a+3Gizu42W02ooQLS1fbYXZF9g3bSkLizx14px825p66cttZf/TDW88vdtOUGHg0EF6GkHhu8xv7kLhSlxD5VoFwFRKveq129rPt3vCz9zf62Us+gHtW2z+72ROlXFcOQ/yDBmxI+S/jY20tIfGpC/2sU9qPwFurOU1uFtOzi9swJQYeDYSXIySenQiJK90TIt8qEK5C4tVurnqdEg9G6LbJz/bwIeyn2i9ys1bzKvNQoukDKxMNj9paQuLLPiNC4t3EzxQw+++r/exGHzpX+O0MAAAAAAAAheDGEwAAAAAAAArBjScAAAAAAAAUghtPAAAAAAAAKERNfX19fWjl4afm3TkaIa/U+YJu+fCuwveolNELRZgvqlEwEJ7i4t5ZT9fsuBElH1vb6G9+uJXEu6vKkISI5hbojCf3c7PbRj+Td5NbHs93rnf65DuXWTxCXqnzBc3v3bjwPSql6263ln7wdg39TAXCExz8St7X4Sfue630gw/2389bTby7miz+eaWvoKx6XHyGm83Z/Sd5Nxn7zTznuVk8FCBFNEJeqfOF921Z/B6VMjAYSVfmib+PUIHwBJ327Z71fK2mlR4h/+vA/3CzrSXeXU2a/3VopS+hvHpc4me/uDfvHt+O3SfiL54AAAAAAABQCG48AQAAAAAAoBDceAIAAAAAAEAhuPEEAAAAAACAQuSPi6dEvisUCE8RjYvXida20jhzf1vtq/ZIComXQ+aQeJQKjtd2mF36CVeIz3O1R8jVNTdv7WdNakre4ryePgZ743vfcLPbp8aisacf2csP2312cy8rjzoRwmwcvOffT8SXlZq6+PX8q0rExVMi3xUKhKeIxsU32pzQuobWI+VyQvuqPZJC4mWQOyQeJYPjQw4t/YTbitfcao+Qq2tu18/P3lpU8hazunRxs74LFrjZ737h1ykHHH6MHy774mZfVxaLNvhZp0axY68Uxyq7dI5fz78qd1w8JfJdqUB4imhcvEfwAStzmpR+LdF91R4pIfEyyB0Sj1LB8ekD9y/5fLXW1M2qPUKurvldW+VmO1npD2M5ZZl/+M7d7f3nuflrwdezhm/72drMv6hHzRO/D3QLPozlwWCQf4du8ev5V8TFAQAAAAAAUEnceAIAAAAAAEAhuPEEAAAAAACAQnDjCQAAAAAAAIWIx8UBAAAAAACAzcBfPAEAAAAAAKAQ3HgCAAAAAABAIbjxBAAAAAAAgEJw4wkAAAAAAACF4MYTAAAAAAAACsGNJwAAAAAAABSCG08AAAAAAAAoBDeeAAAAAAAAUAhuPAEAAAAAAKAQ/x89R0J+ybYotwAAAABJRU5ErkJggg==" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pgm3m = ni.point_group\n", + "ckey = plot.IPFColorKeyTSL(pgm3m, direction=Vector3d.xvector())\n", + "ori = xmap.orientations\n", + "directions = Vector3d(((1, 0, 0), (0, 1, 0), (0, 0, 1)))\n", + "\n", + "fig, axes = plt.subplots(figsize=(15, 5), ncols=3)\n", + "for ax, v, title in zip(axes, directions, (\"X\", \"Y\", \"Z\")):\n", + " ckey.direction = v\n", + " rgb = ckey.orientation2color(ori).reshape(xmap.shape + (3,))\n", + " ax.imshow(rgb)\n", + " ax.axis(\"off\")\n", + " ax.set_title(f\"IPF-{title}\")\n", + "fig.subplots_adjust(wspace=0.05)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "nbsphinx": "hidden", + "tags": [], + "ExecuteTime": { + "end_time": "2023-05-13T12:27:23.209337535Z", + "start_time": "2023-05-13T12:27:23.205863969Z" + } + }, + "outputs": [], + "source": [ + "# Remove files written to disk in this tutorial\n", + "import os\n", + "\n", + "os.remove(temp_dir + \"ni.h5\")\n", + "os.remove(temp_dir + \"ni.ang\")\n", + "os.rmdir(temp_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "outputs": [], + "source": [], + "metadata": { + "collapsed": false, + "ExecuteTime": { + "end_time": "2023-05-13T12:26:31.757585100Z", + "start_time": "2023-05-13T12:26:31.690796926Z" + } + } + } + ], + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/kikuchipy/indexing/__init__.py b/kikuchipy/indexing/__init__.py index 1d8665f9..212fdf92 100644 --- a/kikuchipy/indexing/__init__.py +++ b/kikuchipy/indexing/__init__.py @@ -26,6 +26,8 @@ "NormalizedCrossCorrelationMetric", "NormalizedDotProductMetric", "SimilarityMetric", + "PCAIndexer", + "DIIndexer", "compute_refine_orientation_projection_center_results", "compute_refine_orientation_results", "compute_refine_projection_center_results", @@ -44,6 +46,8 @@ def __getattr__(name): "NormalizedCrossCorrelationMetric": "similarity_metrics", "NormalizedDotProductMetric": "similarity_metrics", "SimilarityMetric": "similarity_metrics", + "PCAIndexer": "di_indexers", + "DIIndexer": "di_indexers", "compute_refine_orientation_projection_center_results": "_refinement._refinement", "compute_refine_orientation_results": "_refinement._refinement", "compute_refine_projection_center_results": "_refinement._refinement", diff --git a/kikuchipy/indexing/_dictionary_indexing_contiguous.py b/kikuchipy/indexing/_dictionary_indexing_contiguous.py new file mode 100644 index 00000000..956be12e --- /dev/null +++ b/kikuchipy/indexing/_dictionary_indexing_contiguous.py @@ -0,0 +1,167 @@ +# Copyright 2019-2023 The kikuchipy developers +# +# This file is part of kikuchipy. +# +# kikuchipy is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# kikuchipy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with kikuchipy. If not, see . + +"""Private tools for approximate dictionary indexing of experimental + patterns to a dictionary of simulated patterns with known orientations. +""" + +from time import sleep, time +from typing import Union + +import dask +import dask.array as da +import numpy as np +from orix.crystal_map import create_coordinate_arrays, CrystalMap +from orix.quaternion import Rotation + +from kikuchipy.indexing.di_indexers import DIIndexer + +from tqdm import tqdm + + +def _pca_dictionary_indexing( + experimental: Union[np.ndarray, da.Array], + experimental_nav_shape: tuple, + dictionary: Union[np.ndarray, da.Array], + step_sizes: tuple, + dictionary_xmap: CrystalMap, + indexer: DIIndexer, + keep_n: int, + n_experimental_per_iteration: int, + navigation_mask: Union[np.ndarray, None], + signal_mask: Union[np.ndarray, None], +) -> CrystalMap: + """Dictionary indexing matching experimental patterns to a + dictionary of simulated patterns of known orientations. + + See :meth:`~kikuchipy.signals.EBSD.custom_di_indexing`. + + Parameters + ---------- + experimental + experimental_nav_shape + dictionary + step_sizes + dictionary_xmap + indexer + keep_n + n_dictionary_per_iteration + keep_dictionary_lazy + n_experimental_per_iteration + navigation_mask + signal_mask + + Returns + ------- + xmap + """ + # handle dictionary reshaping + dictionary_size = dictionary.shape[0] + dictionary = dictionary.reshape((dictionary_size, -1)) + + # handle experimental pattern reshaping + n_experimental_all = int(np.prod(experimental_nav_shape)) + experimental = experimental.reshape((n_experimental_all, -1)) + + # mask the dictionary and the experimental patterns + if signal_mask is not None: + dictionary = dictionary[:, ~signal_mask.ravel()] + with dask.config.set(**{"array.slicing.split_large_chunks": False}): + experimental = experimental[:, ~signal_mask.ravel()] + + # cap out the number of experimental patterns per iteration + n_experimental_per_iteration = min(n_experimental_per_iteration, n_experimental_all) + + # cap out the number of dictionary patterns to keep + if keep_n >= dictionary_size: + raise ValueError(f"keep_n of {keep_n} must be smaller than the dictionary size of {dictionary_size}") + + # prepare the output arrays + simulation_indices = np.zeros((n_experimental_all, keep_n), dtype=np.int32) + distances = np.full((n_experimental_all, keep_n), np.inf, dtype=np.float32) + + # dictionary must fit into memory + if isinstance(dictionary, da.Array): + dictionary = dictionary.compute() + + # calculate the number of loops for the experimental patterns and the start and end indices for each loop + n_experimental_iterations = int(np.ceil(n_experimental_all / n_experimental_per_iteration)) + experimental_chunk_starts = np.cumsum([0] + [n_experimental_per_iteration] * (n_experimental_iterations - 1)) + experimental_chunk_ends = np.cumsum([n_experimental_per_iteration] * n_experimental_iterations) + experimental_chunk_ends[-1] = max(experimental_chunk_ends[-1], n_experimental_all) + + # start timer + time_start = time() + + # set the indexer with the dictionary chunk + indexer(dictionary, keep_n) + + # inner loop over experimental pattern chunks + for exp_start, exp_end in tqdm(zip(experimental_chunk_starts, experimental_chunk_ends), + total=n_experimental_iterations, + desc="Experiment loop", + position=1, + leave=False): + experimental_chunk = experimental[exp_start:exp_end] + # if the experimental chunk is lazy, compute it + if isinstance(experimental_chunk, da.Array): + experimental_chunk = experimental_chunk.compute() + # query the experimental chunk + simulation_indices_mini, distances_mini = indexer.query(experimental_chunk) + # simulation_indices_mini, distances_mini = simulation_indices_mini, distances_mini + # fill the new indices and scores with the mini indices and scores + simulation_indices[exp_start:exp_end] = simulation_indices_mini + distances[exp_start:exp_end] = distances_mini + + # stop timer and calculate indexing speed + query_time = time() - time_start + patterns_per_second = n_experimental_all / query_time + # Without this pause, a part of the red tqdm progressbar background is displayed below this print + sleep(0.2) + print( + f" Indexing speed: {patterns_per_second:.5f} patterns/s, " + ) + + xmap_kw, _ = create_coordinate_arrays(experimental_nav_shape, step_sizes) + if navigation_mask is not None: + nav_mask = ~navigation_mask.ravel() + xmap_kw["is_in_data"] = nav_mask + + rot = Rotation.identity((n_experimental_all, keep_n)) + rot[nav_mask] = dictionary_xmap.rotations[simulation_indices].data + + scores_all = np.empty((n_experimental_all, keep_n), dtype=distances.dtype) + scores_all[nav_mask] = distances + simulation_indices_all = np.empty( + (n_experimental_all, keep_n), dtype=simulation_indices.dtype + ) + simulation_indices_all[nav_mask] = simulation_indices + if keep_n == 1: + rot = rot.flatten() + scores_all = scores_all.squeeze() + simulation_indices_all = simulation_indices_all.squeeze() + xmap_kw["rotations"] = rot + xmap_kw["prop"] = { + "scores": scores_all, + "simulation_indices": simulation_indices_all, + } + else: + xmap_kw["rotations"] = dictionary_xmap.rotations[simulation_indices.squeeze()] + xmap_kw["prop"] = {"scores": distances, "simulation_indices": simulation_indices} + xmap = CrystalMap(phase_list=dictionary_xmap.phases_in_data, **xmap_kw) + + return xmap diff --git a/kikuchipy/indexing/di_indexers/__init__.py b/kikuchipy/indexing/di_indexers/__init__.py new file mode 100644 index 00000000..b1fc9cb8 --- /dev/null +++ b/kikuchipy/indexing/di_indexers/__init__.py @@ -0,0 +1,26 @@ +# Copyright 2019-2023 The kikuchipy developers +# +# This file is part of kikuchipy. +# +# kikuchipy is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# kikuchipy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with kikuchipy. If not, see . + +from kikuchipy.indexing.di_indexers._di_indexer import DIIndexer +from kikuchipy.indexing.di_indexers._pca_indexer import ( + PCAIndexer, +) + +__all__ = [ + "PCAIndexer", + "DIIndexer" +] diff --git a/kikuchipy/indexing/di_indexers/_di_indexer.py b/kikuchipy/indexing/di_indexers/_di_indexer.py new file mode 100644 index 00000000..f1394671 --- /dev/null +++ b/kikuchipy/indexing/di_indexers/_di_indexer.py @@ -0,0 +1,125 @@ +# Copyright 2019-2023 The kikuchipy developers +# +# This file is part of kikuchipy. +# +# kikuchipy is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# kikuchipy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with kikuchipy. If not, see . + +import abc +from typing import List, Union, Tuple + +import numpy as np + + +class DIIndexer(abc.ABC): + """Abstract class implementing a function to match + experimental and simulated EBSD patterns in a dictionary when + the dictionary can be fit in memory. + + When writing a custom dictionary indexing class, the methods listed as + `abstract` below must be implemented. Any number of custom + parameters can be passed. Also listed are the attributes available + to the methods if set properly during initialization or after. + + Parameters + ---------- + datatype + Which data type to cast the patterns to before matching to. + """ + + _allowed_dtypes: List[type] = [] + + def __init__( + self, + datatype: Union[str, np.dtype, type] + ): + """Create a similarity metric matching experimental and + simulated EBSD patterns in a dictionary. + """ + self.dictionary_patterns = None + self._keep_n = None + self._datatype = np.dtype(datatype) + + def __repr__(self): + string = f"{self.__class__.__name__}: {np.dtype(self._datatype).name}" + return string + + def __call__(self, dictionary_patterns: np.ndarray, keep_n: int): + """Ingest a dictionary of simulated patterns to match to. + Must be 2D (1 navigation dimension, 1 signal dimension). + + Parameters + ---------- + dictionary_patterns + Dictionary of simulated patterns to match to. Must be 2D + (1 navigation dimension, 1 signal dimension). + keep_n + Number of best matches to keep. + """ + self._keep_n = keep_n + assert len(dictionary_patterns.shape) == 2 + assert dictionary_patterns.shape[0] > 0 and dictionary_patterns.shape[1] > 0 + self.prepare_dictionary(dictionary_patterns) + + @property + def allowed_dtypes(self) -> List[type]: + """Return the list of allowed array data types used during + matching. + """ + return self._allowed_dtypes + + @property + def datatype(self) -> np.dtype: + """Return or set which data type to cast the patterns to before + matching. + + Data type listed in :attr:`allowed_dtypes`. + """ + return self._datatype + + @datatype.setter + def datatype(self, value: Union[str, np.dtype, type]): + """Set which data type to cast the patterns to before + matching. + """ + self._datatype = np.dtype(value) + + @property + def keep_n(self) -> int: + """Return or set which data type to cast the patterns to before + matching. + + Data type listed in :attr:`allowed_dtypes`. + """ + return self._keep_n + + @keep_n.setter + def keep_n(self, value: int): + """Set which data type to cast the patterns to before + matching. + """ + self._keep_n = value + + @abc.abstractmethod + def prepare_dictionary(self, *args, **kwargs): + """Prepare all dictionary patterns before matching to experimental + patterns in :meth:`match`. + """ + return NotImplemented + + @abc.abstractmethod + def query(self, experimental_patterns: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Match all experimental patterns block to dictionary patterns + and return their similarities. + """ + return NotImplemented diff --git a/kikuchipy/indexing/di_indexers/_pca_indexer.py b/kikuchipy/indexing/di_indexers/_pca_indexer.py new file mode 100644 index 00000000..846cb9e5 --- /dev/null +++ b/kikuchipy/indexing/di_indexers/_pca_indexer.py @@ -0,0 +1,347 @@ +import numpy as np +from numba import njit, prange +import heapq +from scipy.linalg import blas, lapack +from typing import Tuple + +from kikuchipy.indexing.di_indexers import DIIndexer + + +@njit +def _fast_l2(query, dictionary_transposed, dictionary_half_norms): + """ + Compute the L2 squared distance between a query and a dictionary of vectors. + + Parameters + ---------- + query : np.ndarray + dictionary_transposed : np.ndarray + dictionary_half_norms : np.ndarray + + Returns + ------- + np.ndarray + + """ + pseudo_distances = dictionary_half_norms - np.dot(query, dictionary_transposed) + return pseudo_distances + + +@njit +def pseudo_to_real_l2_norm(pseudo_distances, query_pca): + """ + Convert pseudo L2 distances to real L2 distances. + + Parameters + ---------- + pseudo_distances : np.ndarray + query_pca : np.ndarray + + Returns + ------- + np.ndarray + + """ + return pseudo_distances + np.sum(query_pca ** 2, axis=1)[:, np.newaxis] + + +@njit(parallel=True) +def find_smallest_indices_values(matrix, k): + """ + Find the smallest k values in each row of a matrix. + + Parameters + ---------- + matrix + k + + Returns + ------- + smallest_indices : np.ndarray + smallest_values : np.ndarray + + """ + n_rows = matrix.shape[0] + smallest_indices = np.zeros((n_rows, k), dtype=np.int32) + smallest_values = np.zeros((n_rows, k), dtype=matrix.dtype) + + for i in prange(n_rows): + row = -1.0 * matrix[i] + heap = [(val, idx) for idx, val in enumerate(row[:k])] + heapq.heapify(heap) + + for idx in range(k, len(row)): + heapq.heappushpop(heap, (row[idx], idx)) + + sorted_heap = sorted(heap) + + for j in range(k): + smallest_values[i, j] = sorted_heap[j][0] + smallest_indices[i, j] = sorted_heap[j][1] + + return smallest_indices, smallest_values + + +def _eigendecomposition(matrix: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """ + Compute the eigenvector decomposition of a matrix. + Uses LAPACK routines `ssyrk` followed by `dsyevd`. + Parameters + ---------- + matrix + + Returns + ------- + eigenvalues : np.ndarray + eigenvectors : np.ndarray + + """ + """ + Compute the eigenvector decomposition of a matrix. + + Parameters + ---------- + matrix : np.ndarray + The matrix to compute the eigenvector decomposition of. + + Returns + ------- + np.ndarray + The eigenvector decomposition of the matrix. + + """ + if matrix.dtype == np.float32: + # Note: ssyrk computes alpha*A*A.T + beta*C, so alpha and beta should be set to 1 and 0 respectively + cov_matrix = blas.ssyrk(a=matrix.T, alpha=1.0) + + elif matrix.dtype == np.float64: + # Note: dsyrk computes alpha*A*A.T + beta*C, so alpha and beta should be set to 1 and 0 respectively + cov_matrix = blas.dsyrk(a=matrix.T, alpha=1.0) + # Compute eigenvector decomposition with dsyevd + else: + raise TypeError(f"Unsupported dtype {matrix.dtype}") + + w, v, info = lapack.dsyevd(cov_matrix.astype(np.float64), compute_v=True) + + # if the info is not 0, then something went wrong + if info != 0: + raise RuntimeError("dsyevd (symmetric eigenvalue decomposition call to LAPACK) failed with info=%d" % info) + + return w, v + + +class PCAKNNSearch: + """ + A PCA KNN search object. Look up nearest neighbors in a PCA component space. + + """ + + def __init__(self, + n_components: int = 100, + n_neighbors: int = 1, + whiten: bool = True, + datatype: str = 'float32'): + """ + + Parameters + ---------- + n_components : int + The number of components to keep in the PCA. The default is 100. + n_neighbors : int + The number of nearest neighbors to return. The default is 1. + whiten : bool + Whether to whiten the patterns. The default is True. + datatype : str + The dtype to use for matching. The default is 'float32'. + + """ + + self._n_components = n_components + self._n_neighbors = n_neighbors + self._whiten = whiten + self._datatype = datatype + + self._n_features = None + self._singular_values_ = None + self._explained_variance_ = None + self._explained_variance_ratio_ = None + + self._transformation_matrix = None + self._inverse_transformation_matrix = None + + self._dictionary_pca_T = None + self._dictionary_pca_half_norms = None + + def fit(self, dictionary: np.ndarray): + """ + Fit the PCA KNN search object to a dictionary of patterns. + + Parameters + ---------- + dictionary : np.ndarray + The dictionary of patterns to fit the PCA KNN search object to. + + """ + # check the data dimensons + if dictionary.ndim != 2: + raise ValueError("A must be a 2D array but got %d dimensions" % dictionary.ndim) + + # set the data to contiguous and C order as the correct dtype + data = np.ascontiguousarray(dictionary.astype(self._datatype)) + + # get the number of samples and features + n_entries, n_features = data.shape + self._n_features = n_features + + # check the n_components is less than the number of features + if self._n_components > n_features: + raise ValueError(f"n_components, {self._n_components} given, exceeds features in A, {n_features}") + + # compute the eigenvector decomposition + eigenvalues, eigenvectors = _eigendecomposition(data) + + # sort the eigenvalues and eigenvectors + idx = np.argsort(eigenvalues)[::-1] + eigenvalues = eigenvalues[idx] + eigenvectors = eigenvectors[:, idx] + + # keep only the first n_components and make sure the arrays are contiguous in C order + eigenvalues = np.ascontiguousarray(eigenvalues[:self._n_components]).astype(self._datatype) + eigenvectors = np.ascontiguousarray(eigenvectors[:, :self._n_components]).astype(self._datatype) + + # calculate the explained variance + self._explained_variance_ = eigenvalues / (n_entries - 1) + self._explained_variance_ratio_ = self._explained_variance_ / np.sum(self._explained_variance_) + + # calculate the singular values + self._singular_values_ = np.sqrt(eigenvalues) + + # components is an alias for the un-whitened eigenvectors + # calculate the transformation matrix + if self._whiten: + self._transformation_matrix = eigenvectors * np.sqrt(n_entries - 1) / self._singular_values_[ + :self._n_components] + self._inverse_transformation_matrix = np.linalg.pinv(self._transformation_matrix) + else: + self._transformation_matrix = eigenvectors + self._inverse_transformation_matrix = np.linalg.pinv(self._transformation_matrix) + + # make sure the transformation matrix is contiguous in C order + self._transformation_matrix = np.ascontiguousarray(self._transformation_matrix) + self._inverse_transformation_matrix = np.ascontiguousarray(self._inverse_transformation_matrix) + + # calculate the PCA of the dictionary using calls to BLAS + dictionary_pca = blas.dgemm(a=data, b=self._transformation_matrix, alpha=1.0) + + # calculate the norms of the PCA of the dictionary + self._dictionary_pca_half_norms = (np.sum(dictionary_pca ** 2, axis=1) / 2.0).astype(self._datatype) + # make sure the transpose of the dictionary PCA is contiguous in C order + self._dictionary_pca_T = np.ascontiguousarray(dictionary_pca.T).astype(self._datatype) + + def query(self, query: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """ + Query the dictionary. + + Parameters + ---------- + query : np.ndarray + The query to perform. + + Returns + ------- + np.ndarray + The indices of the nearest neighbors. + np.ndarray + The distances to the nearest neighbors. + + """ + query_pca = np.dot(query, self._transformation_matrix) + # calculate distance matrix + pseudo_distances = _fast_l2(query_pca, self._dictionary_pca_T, self._dictionary_pca_half_norms) + # get the indices of the nearest neighbors + indices, pseudo_distances = find_smallest_indices_values(pseudo_distances, self._n_neighbors) + # make the distances true distances by adding the query norms + distances = pseudo_distances + np.sum(query_pca ** 2, axis=1)[:, np.newaxis] + return indices, distances + + +class PCAIndexer(DIIndexer): + """ + Class for performing PCA KNN search using JAX. All nearest neighbor + lookups are performed in the truncated PCA space. This class is a + wrapper around the PCAKNNSearch class that allows for the use of + """ + + def __init__(self, + zero_mean: bool = True, + unit_var: bool = False, + n_components: int = 750, + n_neighbors: int = 1, + whiten: bool = False, + datatype: str = 'float32'): + """ + Parameters + ---------- + n_components : int, optional + The number of components to keep in the PCA. The default is 750. + n_neighbors : int, optional + The number of nearest neighbors to return. The default is 1. + whiten : bool, optional + Whether to whiten the patterns. The default is False. + datatype : str, optional + The datatype to use for matching. The default is 'float32'. + """ + super().__init__(datatype=datatype) + self._pca_knn_search = PCAKNNSearch(n_components=n_components, + n_neighbors=n_neighbors, + whiten=whiten, + datatype=datatype) + self._zero_mean = zero_mean + self._unit_var = unit_var + + def prepare_dictionary(self, dictionary: np.ndarray): + """ + Fit the PCA and KNN lookup objects. + + Parameters + ---------- + dictionary : numpy.ndarray + The dictionary to fit. + + Returns + ------- + None. + + """ + dictionary = dictionary.astype(self.datatype) + # preprocess the dictionary + if self._zero_mean: + dictionary -= np.mean(dictionary, axis=0)[None, :] + if self._unit_var: + dictionary /= np.std(dictionary, axis=0)[None, :] + self._pca_knn_search.fit(dictionary) + + def query(self, query: np.ndarray): + """ + Query the dictionary. + + Parameters + ---------- + query : np.ndarray + The query to perform. Must be of shape (n_samples, n_features). + + Returns + ------- + np.ndarray + The indices of the nearest neighbors. + np.ndarray + The distances to the nearest neighbors. + + """ + query = query.astype(self.datatype) + # preprocess the query + if self._zero_mean: + query -= np.mean(query, axis=0)[None, :] + if self._unit_var: + query /= np.std(query, axis=0)[None, :] + return self._pca_knn_search.query(query) diff --git a/kikuchipy/signals/ebsd.py b/kikuchipy/signals/ebsd.py index 76377f21..d0179276 100644 --- a/kikuchipy/signals/ebsd.py +++ b/kikuchipy/signals/ebsd.py @@ -44,6 +44,7 @@ from kikuchipy.filters.fft_barnes import _fft_filter, _fft_filter_setup from kikuchipy.filters.window import Window from kikuchipy.indexing._dictionary_indexing import _dictionary_indexing +from kikuchipy.indexing._dictionary_indexing_contiguous import _pca_dictionary_indexing from kikuchipy.indexing._hough_indexing import ( _get_pyebsdindex_phaselist, _indexer_is_compatible_with_kikuchipy, @@ -60,6 +61,8 @@ NormalizedCrossCorrelationMetric, NormalizedDotProductMetric, ) +from kikuchipy.indexing.di_indexers import DIIndexer + from kikuchipy.io._io import _save from kikuchipy.pattern import chunk from kikuchipy.pattern.chunk import _average_neighbour_patterns @@ -98,7 +101,6 @@ from kikuchipy.signals._kikuchipy_signal import KikuchipySignal2D, LazyKikuchipySignal2D from kikuchipy.signals.virtual_bse_image import VirtualBSEImage - _logger = logging.getLogger(__name__) @@ -198,10 +200,10 @@ def detector(self) -> EBSDDetector: @detector.setter def detector(self, value: EBSDDetector): if _detector_is_compatible_with_signal( - detector=value, - nav_shape=self._navigation_shape_rc, - sig_shape=self._signal_shape_rc, - raise_if_not=True, + detector=value, + nav_shape=self._navigation_shape_rc, + sig_shape=self._signal_shape_rc, + raise_if_not=True, ): self._detector = value @@ -221,7 +223,7 @@ def xmap(self) -> Union[CrystalMap, None]: @xmap.setter def xmap(self, value: CrystalMap): if _xmap_is_compatible_with_signal( - value, self.axes_manager.navigation_axes[::-1], raise_if_not=True + value, self.axes_manager.navigation_axes[::-1], raise_if_not=True ): self._xmap = value @@ -248,7 +250,7 @@ def static_background(self, value: Union[np.ndarray, da.Array]): # ------------------------ Custom methods ------------------------ # def extract_grid( - self, grid_shape: Union[Tuple[int, int], int], return_indices: bool = False + self, grid_shape: Union[Tuple[int, int], int], return_indices: bool = False ) -> Union[Union[EBSD, LazyEBSD], Tuple[Union[EBSD, LazyEBSD], np.ndarray]]: """Return a new signal with patterns from positions in a grid of shape ``grid_shape`` evenly spaced in navigation space. @@ -290,7 +292,7 @@ def extract_grid( nav_shape = self.axes_manager.navigation_shape if len(grid_shape) != len(nav_shape) or any( - [g > n for g, n in zip(grid_shape, nav_shape)] + [g > n for g, n in zip(grid_shape, nav_shape)] ): raise ValueError( f"grid_shape {grid_shape} must be compatible with navigation shape " @@ -361,7 +363,7 @@ def extract_grid( return out def set_scan_calibration( - self, step_x: Union[int, float] = 1.0, step_y: Union[int, float] = 1.0 + self, step_x: Union[int, float] = 1.0, step_y: Union[int, float] = 1.0 ) -> None: """Set the step size in microns. @@ -421,13 +423,13 @@ def set_detector_calibration(self, delta: Union[int, float]) -> None: dx.offset, dy.offset = -center def remove_static_background( - self, - operation: str = "subtract", - static_bg: Union[np.ndarray, da.Array, None] = None, - scale_bg: bool = False, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + operation: str = "subtract", + static_bg: Union[np.ndarray, da.Array, None] = None, + scale_bg: bool = False, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: """Remove the static background. @@ -505,6 +507,7 @@ def remove_static_background( # Get background pattern if static_bg is None: + static_bg = self.static_background try: if not isinstance(static_bg, (np.ndarray, da.Array)): @@ -555,15 +558,15 @@ def remove_static_background( return s_out def remove_dynamic_background( - self, - operation: str = "subtract", - filter_domain: str = "frequency", - std: Union[int, float, None] = None, - truncate: Union[int, float] = 4.0, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, - **kwargs, + self, + operation: str = "subtract", + filter_domain: str = "frequency", + std: Union[int, float, None] = None, + truncate: Union[int, float] = 4.0, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, + **kwargs, ) -> Union[None, EBSD, LazyEBSD]: """Remove the dynamic background. @@ -679,14 +682,14 @@ def remove_dynamic_background( return s_out def get_dynamic_background( - self, - filter_domain: str = "frequency", - std: Union[int, float, None] = None, - truncate: Union[int, float] = 4.0, - dtype_out: Union[str, np.dtype, type, None] = None, - show_progressbar: Optional[bool] = None, - lazy_output: Optional[bool] = None, - **kwargs, + self, + filter_domain: str = "frequency", + std: Union[int, float, None] = None, + truncate: Union[int, float] = 4.0, + dtype_out: Union[str, np.dtype, type, None] = None, + show_progressbar: Optional[bool] = None, + lazy_output: Optional[bool] = None, + **kwargs, ) -> Union[EBSD, LazyEBSD]: """Return the dynamic background per pattern in a new signal. @@ -771,7 +774,7 @@ def get_dynamic_background( pbar = ProgressBar() if show_progressbar or ( - show_progressbar is None and hs.preferences.General.show_progressbar + show_progressbar is None and hs.preferences.General.show_progressbar ): pbar.register() @@ -786,13 +789,13 @@ def get_dynamic_background( return s_out def fft_filter( - self, - transfer_function: Union[np.ndarray, Window], - function_domain: str, - shift: bool = False, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + transfer_function: Union[np.ndarray, Window], + function_domain: str, + shift: bool = False, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: """Filter patterns in the frequency domain. @@ -900,7 +903,7 @@ def fft_filter( return_lazy = lazy_output or (lazy_output is None and self._lazy) register_pbar = show_progressbar or ( - show_progressbar is None and hs.preferences.General.show_progressbar + show_progressbar is None and hs.preferences.General.show_progressbar ) if not return_lazy and register_pbar: pbar = ProgressBar() @@ -924,13 +927,13 @@ def fft_filter( return s_out def average_neighbour_patterns( - self, - window: Union[str, np.ndarray, da.Array, Window] = "circular", - window_shape: Tuple[int, ...] = (3, 3), - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, - **kwargs, + self, + window: Union[str, np.ndarray, da.Array, Window] = "circular", + window_shape: Tuple[int, ...] = (3, 3), + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, + **kwargs, ) -> Union[None, EBSD, LazyEBSD]: """Average patterns with its neighbours within a window. @@ -1066,7 +1069,7 @@ def average_neighbour_patterns( return_lazy = lazy_output or (lazy_output is None and self._lazy) register_pbar = show_progressbar or ( - show_progressbar is None and hs.preferences.General.show_progressbar + show_progressbar is None and hs.preferences.General.show_progressbar ) if not return_lazy and register_pbar: pbar = ProgressBar() @@ -1094,12 +1097,12 @@ def average_neighbour_patterns( return s_out def downsample( - self, - factor: int, - dtype_out: Optional[str] = None, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + factor: int, + dtype_out: Optional[str] = None, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: r"""Downsample the pattern shape by an integer factor and rescale intensities to fill the data type range. @@ -1202,12 +1205,12 @@ def downsample( return s_out def get_neighbour_dot_product_matrices( - self, - window: Optional[Window] = None, - zero_mean: bool = True, - normalize: bool = True, - dtype_out: Union[str, np.dtype, type] = "float32", - show_progressbar: Optional[bool] = None, + self, + window: Optional[Window] = None, + zero_mean: bool = True, + normalize: bool = True, + dtype_out: Union[str, np.dtype, type] = "float32", + show_progressbar: Optional[bool] = None, ) -> Union[np.ndarray, da.Array]: """Get an array with dot products of a pattern and its neighbours within a window. @@ -1279,7 +1282,7 @@ def get_neighbour_dot_product_matrices( if not self._lazy: pbar = ProgressBar() if show_progressbar or ( - show_progressbar is None and hs.preferences.General.show_progressbar + show_progressbar is None and hs.preferences.General.show_progressbar ): pbar.register() @@ -1293,9 +1296,9 @@ def get_neighbour_dot_product_matrices( return dp_matrices def get_image_quality( - self, - normalize: bool = True, - show_progressbar: Optional[bool] = None, + self, + normalize: bool = True, + show_progressbar: Optional[bool] = None, ) -> Union[np.ndarray, da.Array]: """Compute the image quality map of patterns in an EBSD scan. @@ -1359,13 +1362,13 @@ def get_image_quality( return image_quality_map.data def get_average_neighbour_dot_product_map( - self, - window: Optional[Window] = None, - zero_mean: bool = True, - normalize: bool = True, - dtype_out: Union[str, np.dtype, type] = "float32", - dp_matrices: Optional[np.ndarray] = None, - show_progressbar: Optional[bool] = None, + self, + window: Optional[Window] = None, + zero_mean: bool = True, + normalize: bool = True, + dtype_out: Union[str, np.dtype, type] = "float32", + dp_matrices: Optional[np.ndarray] = None, + show_progressbar: Optional[bool] = None, ) -> Union[np.ndarray, da.Array]: """Get a map of the average dot product between patterns and their neighbours within an averaging window. @@ -1461,7 +1464,7 @@ def get_average_neighbour_dot_product_map( if not self._lazy: pbar = ProgressBar() if show_progressbar or ( - show_progressbar is None and hs.preferences.General.show_progressbar + show_progressbar is None and hs.preferences.General.show_progressbar ): pbar.register() @@ -1475,10 +1478,10 @@ def get_average_neighbour_dot_product_map( return adp def plot_virtual_bse_intensity( - self, - roi: BaseInteractiveROI, - out_signal_axes: Union[Iterable[int], Iterable[str], None] = None, - **kwargs, + self, + roi: BaseInteractiveROI, + out_signal_axes: Union[Iterable[int], Iterable[str], None] = None, + **kwargs, ) -> None: """Plot an interactive virtual backscatter electron (VBSE) image formed from intensities within a specified and adjustable @@ -1537,9 +1540,9 @@ def plot_virtual_bse_intensity( out.plot(**kwargs) def get_virtual_bse_intensity( - self, - roi: BaseInteractiveROI, - out_signal_axes: Union[Iterable[int], Iterable[str], None] = None, + self, + roi: BaseInteractiveROI, + out_signal_axes: Union[Iterable[int], Iterable[str], None] = None, ) -> VirtualBSEImage: """Get a virtual backscatter electron (VBSE) image formed from intensities within a region of interest (ROI) on the detector. @@ -1582,13 +1585,13 @@ def get_virtual_bse_intensity( return vbse_sum def hough_indexing( - self, - phase_list: PhaseList, - indexer: "EBSDIndexer", - chunksize: int = 528, - verbose: int = 1, - return_index_data: bool = False, - return_band_data: bool = False, + self, + phase_list: PhaseList, + indexer: "EBSDIndexer", + chunksize: int = 528, + verbose: int = 1, + return_index_data: bool = False, + return_band_data: bool = False, ) -> Union[ CrystalMap, Tuple[CrystalMap, np.ndarray], @@ -1713,11 +1716,11 @@ def hough_indexing( return xmap def hough_indexing_optimize_pc( - self, - pc0: Union[list, tuple, np.ndarray], - indexer: "EBSDIndexer", - batch: bool = False, - method: str = "Nelder-Mead", + self, + pc0: Union[list, tuple, np.ndarray], + indexer: "EBSDIndexer", + batch: bool = False, + method: str = "Nelder-Mead", ) -> "EBSDDetector": """Return a detector with one projection center (PC) per pattern optimized using Hough indexing from :mod:`pyebsdindex`. @@ -1817,15 +1820,15 @@ def hough_indexing_optimize_pc( return new_detector def dictionary_indexing( - self, - dictionary: EBSD, - metric: Union[SimilarityMetric, str] = "ncc", - keep_n: int = 20, - n_per_iteration: Optional[int] = None, - navigation_mask: Optional[np.ndarray] = None, - signal_mask: Optional[np.ndarray] = None, - rechunk: bool = False, - dtype: Union[str, np.dtype, type, None] = None, + self, + dictionary: EBSD, + metric: Union[SimilarityMetric, str] = "ncc", + keep_n: int = 20, + n_per_iteration: Optional[int] = None, + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, + rechunk: bool = False, + dtype: Union[str, np.dtype, type, None] = None, ) -> CrystalMap: """Index patterns by matching each pattern to a dictionary of simulated patterns of known orientations @@ -1837,6 +1840,10 @@ def dictionary_indexing( One EBSD signal with dictionary patterns. The signal must have a 1D navigation axis, an :attr:`xmap` property with crystal orientations set, and equal detector shape. + approx + Whether to use the approximate indexing algorithm (hnswlib). + If ``True``, then the keep_n parameter is used to set the + number of neighbors to search for. metric Similarity metric, by default ``"ncc"`` (normalized cross-correlation). ``"ndp"`` (normalized dot product) is @@ -1854,7 +1861,7 @@ class methods. See Number of dictionary patterns to compare to all experimental patterns in each indexing iteration. If not given, and the dictionary is a ``LazyEBSD`` signal, it is equal to the - chunk size of the first pattern array axis, while if if is + chunk size of the first pattern array axis, while if it is an ``EBSD`` signal, it is set equal to the number of dictionary patterns, yielding only one iteration. This parameter can be increased to use less memory during @@ -1975,24 +1982,143 @@ class methods. See return xmap + def pca_dictionary_indexing( + self, + dictionary: EBSD, + indexer: DIIndexer, + keep_n: int = 5, + keep_dictionary_lazy: bool = False, + n_experimental_per_iteration: Optional[int] = None, + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, + ) -> CrystalMap: + """Index patterns by matching each pattern to a dictionary of + simulated patterns of known orientations + :cite:`chen2015dictionary,jackson2019dictionary`. + + Parameters + ---------- + dictionary + One EBSD signal with dictionary patterns. The signal must + have a 1D navigation axis, an :attr:`xmap` property with + crystal orientations set, and equal detector shape. + indexer + A custom indexer that implements the :class:`DIIndexer` + keep_n + Number of best matches to keep, by default 5 or the number + of dictionary patterns if fewer than 5 are available. + keep_dictionary_lazy + If ``True``, the dictionary will be kept in memory as a + :class:`dask.array.Array` with the same chunks as the + dictionary signal. + n_experimental_per_iteration + Number of experimental patterns to use per iteration. If not + given, all experimental patterns will be used. + navigation_mask + A boolean mask equal to the signal's navigation (map) shape, + where only patterns equal to ``False`` are indexed. If not + given, all patterns are indexed. + signal_mask + A boolean mask equal to the experimental patterns' detector + shape, where only pixels equal to ``False`` are matched. If + not given, all pixels are used. + + Returns + ------- + xmap + A crystal map with ``keep_n`` rotations per point with the + sorted best matching orientations in the dictionary. The + corresponding best scores and indices into the dictionary + are stored in the ``xmap.prop`` dictionary as ``"scores"`` + and ``"simulation_indices"``. + + See Also + -------- + refine_orientation + refine_projection_center + refine_orientation_projection_center + kikuchipy.indexing.merge_crystal_maps : + Merge multiple single phase crystal maps into one multiphase map. + kikuchipy.indexing.orientation_similarity_map : + Calculate an orientation similarity map. + """ + am_exp = self.axes_manager + am_dict = dictionary.axes_manager + dict_size = am_dict.navigation_size + + keep_n = min(keep_n, dict_size) + + nav_shape_exp = am_exp.navigation_shape[::-1] + if navigation_mask is not None: + if navigation_mask.shape != nav_shape_exp: + raise ValueError( + f"The navigation mask shape {navigation_mask.shape} and the " + f"signal's navigation shape {nav_shape_exp} must be identical" + ) + elif navigation_mask.all(): + raise ValueError( + "The navigation mask must allow for indexing of at least one " + "pattern (at least one value equal to `False`)" + ) + elif not isinstance(navigation_mask, np.ndarray): + raise ValueError("The navigation mask must be a NumPy array") + + if signal_mask is not None: + if not isinstance(signal_mask, np.ndarray): + raise ValueError("The signal mask must be a NumPy array") + + sig_shape_exp = am_exp.signal_shape[::-1] + sig_shape_dict = am_dict.signal_shape[::-1] + if sig_shape_exp != sig_shape_dict: + raise ValueError( + f"Experimental {sig_shape_exp} and dictionary {sig_shape_dict} signal " + "shapes must be identical" + ) + + dict_xmap = dictionary.xmap + if dict_xmap is None or dict_xmap.shape != (dict_size,): + raise ValueError( + "Dictionary signal must have a non-empty `EBSD.xmap` attribute of equal" + " size as the number of dictionary patterns, and both the signal and " + "crystal map must have only one navigation dimension" + ) + + with dask.config.set(**{"array.slicing.split_large_chunks": False}): + xmap = _pca_dictionary_indexing( + experimental=self.data, + experimental_nav_shape=am_exp.navigation_shape[::-1], + dictionary=dictionary.data, + step_sizes=tuple(a.scale for a in am_exp.navigation_axes[::-1]), + dictionary_xmap=dictionary.xmap, + indexer=indexer, + keep_n=keep_n, + n_experimental_per_iteration=n_experimental_per_iteration, + navigation_mask=navigation_mask, + signal_mask=signal_mask, + ) + + xmap.scan_unit = _get_navigation_axes_unit(am_exp) + + return xmap + def refine_orientation( - self, - xmap: CrystalMap, - detector: EBSDDetector, - master_pattern: "EBSDMasterPattern", - energy: Union[int, float], - navigation_mask: Optional[np.ndarray] = None, - signal_mask: Optional[np.ndarray] = None, - pseudo_symmetry_ops: Optional[Rotation] = None, - method: Optional[str] = "minimize", - method_kwargs: Optional[dict] = None, - trust_region: Union[tuple, list, np.ndarray, None] = None, - initial_step: Union[float] = None, - rtol: float = 1e-4, - maxeval: Optional[int] = None, - compute: bool = True, - rechunk: bool = True, - chunk_kwargs: Optional[dict] = None, + self, + xmap: CrystalMap, + detector: EBSDDetector, + master_pattern: "EBSDMasterPattern", + energy: Union[int, float], + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, + pseudo_symmetry_ops: Optional[Rotation] = None, + method: Optional[str] = "minimize", + method_kwargs: Optional[dict] = None, + trust_region: Union[tuple, list, np.ndarray, None] = None, + initial_step: Union[float] = None, + rtol: float = 1e-4, + maxeval: Optional[int] = None, + compute: bool = True, + rechunk: bool = True, + chunk_kwargs: Optional[dict] = None, ) -> Union[CrystalMap, da.Array]: r"""Refine orientations by searching orientation space around the best indexed solution using fixed projection centers. @@ -2177,22 +2303,22 @@ def refine_orientation( ) def refine_projection_center( - self, - xmap: CrystalMap, - detector: EBSDDetector, - master_pattern: "EBSDMasterPattern", - energy: Union[int, float], - navigation_mask: Optional[np.ndarray] = None, - signal_mask: Optional[np.ndarray] = None, - method: Optional[str] = "minimize", - method_kwargs: Optional[dict] = None, - trust_region: Union[tuple, list, np.ndarray, None] = None, - initial_step: Union[float] = None, - rtol: float = 1e-4, - maxeval: Optional[int] = None, - compute: bool = True, - rechunk: bool = True, - chunk_kwargs: Optional[dict] = None, + self, + xmap: CrystalMap, + detector: EBSDDetector, + master_pattern: "EBSDMasterPattern", + energy: Union[int, float], + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, + method: Optional[str] = "minimize", + method_kwargs: Optional[dict] = None, + trust_region: Union[tuple, list, np.ndarray, None] = None, + initial_step: Union[float] = None, + rtol: float = 1e-4, + maxeval: Optional[int] = None, + compute: bool = True, + rechunk: bool = True, + chunk_kwargs: Optional[dict] = None, ) -> Union[Tuple[np.ndarray, EBSDDetector, np.ndarray], da.Array]: """Refine projection centers by searching the parameter space using fixed orientations. @@ -2366,23 +2492,23 @@ def refine_projection_center( ) def refine_orientation_projection_center( - self, - xmap: CrystalMap, - detector: EBSDDetector, - master_pattern: "EBSDMasterPattern", - energy: Union[int, float], - navigation_mask: Optional[np.ndarray] = None, - signal_mask: Optional[np.ndarray] = None, - pseudo_symmetry_ops: Optional[Rotation] = None, - method: Optional[str] = "minimize", - method_kwargs: Optional[dict] = None, - trust_region: Union[tuple, list, np.ndarray, None] = None, - initial_step: Union[tuple, list, np.ndarray, None] = None, - rtol: Optional[float] = 1e-4, - maxeval: Optional[int] = None, - compute: bool = True, - rechunk: bool = True, - chunk_kwargs: Optional[dict] = None, + self, + xmap: CrystalMap, + detector: EBSDDetector, + master_pattern: "EBSDMasterPattern", + energy: Union[int, float], + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, + pseudo_symmetry_ops: Optional[Rotation] = None, + method: Optional[str] = "minimize", + method_kwargs: Optional[dict] = None, + trust_region: Union[tuple, list, np.ndarray, None] = None, + initial_step: Union[tuple, list, np.ndarray, None] = None, + rtol: Optional[float] = 1e-4, + maxeval: Optional[int] = None, + compute: bool = True, + rechunk: bool = True, + chunk_kwargs: Optional[dict] = None, ) -> Union[Tuple[CrystalMap, EBSDDetector], da.Array]: r"""Refine orientations and projection centers simultaneously by searching the orientation and PC parameter space. @@ -2586,11 +2712,11 @@ def refine_orientation_projection_center( # ------ Methods overwritten from hyperspy.signals.Signal2D ------ # def save( - self, - filename: Optional[str] = None, - overwrite: Optional[bool] = None, - extension: Optional[str] = None, - **kwargs, + self, + filename: Optional[str] = None, + overwrite: Optional[bool] = None, + extension: Optional[str] = None, + **kwargs, ) -> None: """Write the signal to file in the specified format. @@ -2646,9 +2772,9 @@ def save( _save(filename, self, overwrite=overwrite, **kwargs) def get_decomposition_model( - self, - components: Union[int, List[int], None] = None, - dtype_out: Union[str, np.dtype, type] = "float32", + self, + components: Union[int, List[int], None] = None, + dtype_out: Union[str, np.dtype, type] = "float32", ) -> Union[EBSD, LazyEBSD]: """Get the model signal generated with the selected number of principal components from a decomposition. @@ -2860,12 +2986,12 @@ def _slicer(self, *args, **kwargs): # ------------------------ Private methods ----------------------- # def _check_refinement_parameters( - self, - xmap: CrystalMap, - detector: EBSDDetector, - master_pattern: "EBSDMasterPattern", - navigation_mask: Optional[np.ndarray] = None, - signal_mask: Optional[np.ndarray] = None, + self, + xmap: CrystalMap, + detector: EBSDDetector, + master_pattern: "EBSDMasterPattern", + navigation_mask: Optional[np.ndarray] = None, + signal_mask: Optional[np.ndarray] = None, ) -> np.ndarray: """Check compatibility of refinement parameters with refinement. @@ -2952,11 +3078,11 @@ def _check_refinement_parameters( return points_to_refine def _prepare_patterns_for_refinement( - self, - points_to_refine: np.ndarray, - signal_mask: Union[np.ndarray, None], - rechunk: bool, - chunk_kwargs: Optional[dict] = None, + self, + points_to_refine: np.ndarray, + signal_mask: Union[np.ndarray, None], + rechunk: bool, + chunk_kwargs: Optional[dict] = None, ) -> Tuple[da.Array, np.ndarray]: """Prepare pattern array and mask for refinement. @@ -3029,13 +3155,13 @@ def _prepare_patterns_for_refinement( return patterns, signal_mask def _prepare_metric( - self, - metric: Union[SimilarityMetric, str], - navigation_mask: Union[np.ndarray, None], - signal_mask: Union[np.ndarray, None], - dtype: Union[str, np.dtype, type, None], - rechunk: bool, - n_dictionary_patterns: int, + self, + metric: Union[SimilarityMetric, str], + navigation_mask: Union[np.ndarray, None], + signal_mask: Union[np.ndarray, None], + dtype: Union[str, np.dtype, type, None], + rechunk: bool, + n_dictionary_patterns: int, ) -> SimilarityMetric: metrics = { "ncc": NormalizedCrossCorrelationMetric, @@ -3071,7 +3197,7 @@ def _prepare_metric( @staticmethod def _get_sum_signal( - signal, out_signal_axes: Optional[List] = None + signal, out_signal_axes: Optional[List] = None ) -> hs.signals.Signal2D: out = signal.nansum(signal.axes_manager.signal_axes) if out_signal_axes is None: @@ -3091,17 +3217,17 @@ def _get_sum_signal( # purposes def rescale_intensity( - self, - relative: bool = False, - in_range: Union[Tuple[int, int], Tuple[float, float], None] = None, - out_range: Union[Tuple[int, int], Tuple[float, float], None] = None, - dtype_out: Union[ - str, np.dtype, type, Tuple[int, int], Tuple[float, float], None - ] = None, - percentiles: Union[Tuple[int, int], Tuple[float, float], None] = None, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + relative: bool = False, + in_range: Union[Tuple[int, int], Tuple[float, float], None] = None, + out_range: Union[Tuple[int, int], Tuple[float, float], None] = None, + dtype_out: Union[ + str, np.dtype, type, Tuple[int, int], Tuple[float, float], None + ] = None, + percentiles: Union[Tuple[int, int], Tuple[float, float], None] = None, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: return super().rescale_intensity( relative, @@ -3115,13 +3241,13 @@ def rescale_intensity( ) def normalize_intensity( - self, - num_std: int = 1, - divide_by_square_root: bool = False, - dtype_out: Union[str, np.dtype, type, None] = None, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + num_std: int = 1, + divide_by_square_root: bool = False, + dtype_out: Union[str, np.dtype, type, None] = None, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: return super().normalize_intensity( num_std, @@ -3133,13 +3259,13 @@ def normalize_intensity( ) def adaptive_histogram_equalization( - self, - kernel_size: Optional[Union[Tuple[int, int], List[int]]] = None, - clip_limit: Union[int, float] = 0, - nbins: int = 128, - show_progressbar: Optional[bool] = None, - inplace: bool = True, - lazy_output: Optional[bool] = None, + self, + kernel_size: Optional[Union[Tuple[int, int], List[int]]] = None, + clip_limit: Union[int, float] = 0, + nbins: int = 128, + show_progressbar: Optional[bool] = None, + inplace: bool = True, + lazy_output: Optional[bool] = None, ) -> Union[None, EBSD, LazyEBSD]: return super().adaptive_histogram_equalization( kernel_size, @@ -3156,7 +3282,7 @@ def as_lazy(self, *args, **kwargs) -> LazyEBSD: def change_dtype(self, *args, **kwargs) -> None: super().change_dtype(*args, **kwargs) if isinstance(self, EBSD) and isinstance( - self._static_background, (np.ndarray, da.Array) + self._static_background, (np.ndarray, da.Array) ): self._static_background = self._static_background.astype(self.data.dtype) @@ -3180,12 +3306,12 @@ def compute(self, *args, **kwargs) -> None: super().compute(*args, **kwargs) def get_decomposition_model_write( - self, - components: Union[int, List[int], None] = None, - dtype_learn: Union[str, np.dtype, type] = "float32", - mbytes_chunk: int = 100, - dir_out: Optional[str] = None, - fname_out: Optional[str] = None, + self, + components: Union[int, List[int], None] = None, + dtype_learn: Union[str, np.dtype, type] = "float32", + mbytes_chunk: int = 100, + dir_out: Optional[str] = None, + fname_out: Optional[str] = None, ) -> None: """Write the model signal generated from the selected number of principal components directly to an ``.hspy`` file. @@ -3274,11 +3400,11 @@ def get_decomposition_model_write( def _update_custom_attributes( - attributes: dict, - nav_slices: Union[slice, tuple, None] = None, - sig_slices: Union[slice, tuple, None] = None, - new_nav_shape: Optional[tuple] = None, - new_sig_shape: Optional[tuple] = None, + attributes: dict, + nav_slices: Union[slice, tuple, None] = None, + sig_slices: Union[slice, tuple, None] = None, + new_nav_shape: Optional[tuple] = None, + new_sig_shape: Optional[tuple] = None, ) -> dict: """Update dictionary of custom attributes after slicing the signal data.