{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Merging with a Robust Error Model \n", "\n", "In the previous example, we computed the common merging statistics $CC_{1/2}$ and $CC_{anom}$ to explore a dataset with significant anomalous signal from native sulfur atoms. One assumption that we made while merging is that the scaled reflection observations are normally distributed about the mean. This assumption is consistent with the merging strategy used by [AIMLESS](https://doi.org/10.1107/S0907444913000061), which was used to scale the data in the first place. In this example, we will explore whether the scaled reflection observations are normally distributed, and whether we can improve the anomalous signal ($CC_{anom}$) by using a different error model.\n", "\n", "_Note:_ See [pytorch.org](https://pytorch.org/) for customizable PyTorch installation instructions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import seaborn as sns\n", "sns.set_context(\"notebook\", font_scale=1.3)\n", "from tqdm.notebook import tqdm" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: torch in /Users/jgreisman/miniconda3/envs/test/lib/python3.9/site-packages (1.10.1)\r\n", "Requirement already satisfied: typing-extensions in /Users/jgreisman/miniconda3/envs/test/lib/python3.9/site-packages (from torch) (4.0.1)\r\n" ] } ], "source": [ "# Install PyTorch to running kernel\n", "import sys\n", "!{sys.executable} -m pip install torch" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import torch" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import reciprocalspaceship as rs" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9.18\n" ] } ], "source": [ "print(rs.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Normal Error Model\n", "\n", "It is common in merging to assume that scaled intensities are normally distributed about the true mean. We can assess the validity of this assumption by looking at the residuals between scaled intensities and the maximum likelihood esimate of their true intensity." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "ds = rs.read_mtz(\"data/HEWL_unmerged.mtz\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "merged_normal = rs.algorithms.merge(ds)\n", "merged_normal = merged_normal.stack_anomalous()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "ds.hkl_to_asu(anomalous=True, inplace=True)\n", "ds[\"IML\"] = merged_normal.loc[ds.index, \"I\"]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAFUCAYAAAAH7mfYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABLyElEQVR4nO3deXhU5fn/8fedjUDCvoRAhLAJalgl4i6uUMW60FZxrVZt1WpdqvVb/SFYrFqtVqu27lqtW8Xiilut4gpUlEUEQYiIEBJIAiRkz/P740wwhAkkYTJnls/rus6VyTkn59xDyMw9z3I/5pxDREREJJIl+B2AiIiIyO4oYREREZGIp4RFREREIp4SFhEREYl4SX4H4AczawfkAuuBWp/DERGRHSUCmcB851yl38FIZIibhMXMpgE3+h2HiIg022HAh34HIZHB4nFas5kNAlZ+8MEHZGVl+R2OiIg0sHbtWg477DCAwc65b/yORyJD3LSwNFILkJWVRXZ2ts+hiIhIE9RlL9tp0K2IiIhEPCUsIiIiEvGUsIiIiEjEU8IiIiIiEU8Ji4iIiEQ8JSwiIiIS8eImYTGzaWbmzMwBq/2OR0RERJovbhIW59w055w55wwY4Hc8IiIi0nxxk7CIiIhI9FLCIiIiIhFPCYuIiIhEPCUsIiIiEvGUsIiIiEjEU8IiIiIiES/J7wBay8zygK14y4/XOOfG+huRiIiItJWoTVgCjnTObfQ7CBEREWlb6hISERGRiBf2hMXMDjezV8xsXaBU/ulNnPdjM1tsZpVmtsLMzm90igPeMrPPzOyito9cRERE/OJHl1AasBB4FHgx2AlmNi5wbAbwLHAM8KCZFTnnZgVOO9Q5972Z9QLeNrNlzrk5Qa7VBejSaHdWCJ6HiIiIhEnYExbn3GxgNoCZNXXalcAHzrlpge+XBZKYa4BZget8H/haYGb/Bg4AdkpYgCuAG0MTvUh0yM3NZf369Tvtz8zMZP78+T5EJCKyZyJ10O3BwN8a7XsDeMzMkoEUIME5t9XM0oDjgJuauNZfgMcb7csCPghZtCI+aCopASgoKOCQQw7Zaf+KFSvaOiwRkTYRqQlLb2BDo335QDLQA2gP/DvQQpMEPO2ceyPYhZxzJUCJmU1DLS0SQ9avX8+QIUOaPBZMQUEBWVnBe0TV+iIikSxSE5Zdcs6tAka28GemAdMAzCwbWB3quETCpraW9nV1pNbW7nTIATjXxI/VNpnkqPVFRCJZpCYs+UBGo30ZQA2guisSN3Jzc6lcu5Z9qqsZVlPjfa2uZmhNDSucgyZaUgqB7774glVpaaxOT+ebtDTy0tLCG7yISAhFasLyMXAscEuDfROBec656tZcUF1CEjVqauA//4F//pNXFywgo64u6GnbAEvYuTJBonP0dI6eJSWMKSnZvr8O+Br4bNUq3unVi7z09DYJX0SkLYQ9YTGzdGBwg139zWwUUOqcWxnYdxfwkZlNBZ4DjgbOAH7S2vuqS0gimnMwfz7885/w3HOwwRvClQFsS0xkVVraDq0lq9PSeP2TTzj88MODXmvVnDmckZPDgNJSBpWVMbC0lL3KyxnmHMPWrOHMNWv4Ji2NdzIy+E+vXhSmpmp8i4hEND9aWMYC/23w/a2B7X1gPIBzbq6ZTQZuBm4AvgN+2aAGS4uphUUi0qZNcN998NRT0HAMyZAhcOaZjL//fhL22QfXdAmAnZmx1oxPu3fn0+7dt+9Orquj5oMPuKp3b44oLGRQWRmDVq3iwlWrWNS5M3fX1LBm0CDqgrTaaHyLiPjNjzos7wG7ffV1zr0EvBTC+05DLSwSKUpL4S9/gdtvhy1bvH29e8Ppp8OZZ8L++4MZKx96iCEtSVZ2oTohgTlmuKFDuWfIEA4oKuKYDRs4aNMmRm3ezGPAmvnzeWTAAOb07Akhuq+ISChE6hgWkdhUVQUPPQR/+MP2bh8mTuTSb77h1bIyav/1L/jXv7afXlBQ0OSsnj1RnZDARz168FGPHnSoqeGIwkJOWb6cIeXlTF+6lGUdO/LQwIEs6No15PcWEWmNuFn80MymBdYucqh1RcKtrg6efhr22Qd+/WsvWRk3Dv77X5g9m5e2bWPg3nszZMiQHbbaINOWQ21bUhKzMzPJMePOIUPYlJLCsK1b+fPChfxp4UL23rq1zWMQEdmduElYnHPTnHPmnDNggN/xSBz58kvIzfW6elatgmHD4MUX4ZNPYPx4v6PbrsaMV/r25axx43howABKExPJLS7mgc8+487i4h+6rkREfBA3CYtI2DkH994LY8fCggWQlQWPPAKLF8Mpp0TsGJGKxESe7t+fMw88kGf32osqM362bRuMGgUff+x3eCISp5SwiLSFDRtg0iS47DKoqIDzz4evvvK+JkXH0LEtyck8MGgQF4wdy+cAq1dTe8gh3NmpE/379iUrK4usrCxyc3P9DlVE4kDcJCwawyJh89prMHw4vP46dO3qDaJ95BGI0kJt36WlcZAZz+y1FwZctXUrr5WWcnhWFkOGDGly3SIRkVCKm4RFY1ikzZWXewNqJ02CwkI48khYtAh+4tU7zM3N3d4q0XgrKCjwOfhdqzbjwUGD+O3IkRSmpJCzZQsP/e9/HJef3+S6RSIioRQdbdMikS4/H044wRurkpwMN98MV18NDYqwtWZ15Ujzedeu/CI3l6uXL+eIjRv5v2XL6JOWBrW1kJjod3giEsPipoVFpM0sXw4HHeQlK4MGwaefwjXX7JCsxJKtyclM228//jR0KFVmnFtW5rUilZf7HZqIxLDYfEUVCZePP4aDD4a8PDjgAG+q8pgxfkfV9syYnZnJ1SNHUmIGs2bB0UfDRi2mLiJtI24SFg26lZCrf5MuKvLGrbz7LvTs6XdUYbWkSxdO7tkT+vXzkrVDDoHV+vMSkdCLm4RFg24lpO6/HyZP9qYsX3gh/PvfkJbmd1S+WJmc7CUrI0bA11//0D0mIhJCcZOwiISEc/D738Oll3rl9m+6CR54IGpqq7SZPn1gzhyvxWnDBjjiCHjzTb+jEpEYEuevsiIt4BxcfrlXvTYxER580CsE10Bubm6TM37aaiFDvxUUFJCVlQVAsnP8uX17Ti0tpeZHPyLppZfgxBN9jlBEYoESFpHmmjrVS1ZSUrwuoOOP3+mUWJi63FK1tbU7POd7naP2m2/46dq18NOfwhtvRNSaSSISneKmS0iDbmWP3HEHzJjhtaw891zQZEU8zoz7Bw3iyQ4doLLSa2GZP9/vsEQkysVNwqJBt9JqDz3k1VUBeOwxOPlkX8OJCmZc36ULTJkCpaUwcaK3arWISCvFTcIi0irPPQe//KX3+K9/hbPP9jeeKFJnBk884U35LiqCY4+FVav8DktEopQSFpGmvP46nHWWN9h2xgxvnSBpmeRkeP55bwzL+vVwzDGwbp3fUYlIFFLCIhLMnDlenZWaGvjtb72pzNI67dvDyy9Dbq5XVO7YY2HTJr+jEpEoo1lCIo19+aU3ULSiAi64AP70JzDbfjgepy7vsY4dYfZsOPxwWLrUWyjyvfcgNdXvyEQkSihhEWmouNgbVLtli9fC8ve/75CsQHxOXQ6J7t3h7be9tZfmzoWLL4ZHH93p31dEJBh1CYnUq62FM8+ElSth1Cj4xz+8acwSOn36eGswtW8Pjz8O993nd0QiEiXiJmFRHRbZrf/3/7xui+7dvcJwHTr4HVFsGjXKa1kBuOIKr2tIRGQ34qZLyDk3DZgGYGbZKGmRhv71L7jlFq9F5fnnITvb74iiXsOS/Y1lZmYy/5pr4PbbvWq4n33mrfgsItKEuElYRJq0aBH8/Ofe4z//GY46ytdwYkXjkv0NrVixwksQFy6Et96CU06BDz/0uopERIJQwiLxrajIG2S7bRucc463uGFAU7OBNBMoRBIT4ZlnvOnOCxbARRd544Y0CFdEglDCIvGrpgZOP92rDbL//jvNCGpqNpBmAoVQt27eINyDDoKnnvJ+D1dc4XdUIhKB4mbQrchOrr/em2bbs6c3yFbdEf4YPtybMQRekb533/U1HBGJTEpYJD699ZZXEC4pCV54Afbay++I4ttPfuJVE66t9ZZDUCVcEWlECYvEn02bfhhke9NNXvVV8d9NN8Ehh3hrDl10kbeGk4hIgBIWiS/OeW+G69fDoYfCtdf6HZHUS0yEJ5/0yvi/+KK30rOISEBUJyxmlmhmn5vZq37HIlHiiSe8N8OOHb03R1WyjSwDBsC993qPL7sMVq3yNx4RiRjRPkvoN8BXQCe/A5EosGqV9yYI3ptidrYWMoxEZ58Nr77qFfM7+2x4/31vrJGIxLWofRUwsyzgBOBm4Cqfw5FIV1PjvfmVlnqVVc8+G9BChn7ZbRXcN9+Ejz6Cjz+GW2+FG24Ic4QiEmnC3iVkZoeb2Stmti6wts/pTZz3YzNbbGaVZrbCzM5vdMpfgGuBuraOWWLArbd6b359+gRdgVnCq74KbrBt/fr1Xn2W+jEs06bBvHm+xisi/vNjDEsasBC4tKkTzGwc8CIwExgJ3A08aGYnB45PAgqcc5/t7mZm1sXMshtuQPCPdhKb5s3z3vTAexPs1s3XcKSZjjkGrrzyh6nOZWV+RyQiPgp7wuKcm+2cu8E59+9dnHYl8IFzbppzbplz7l7gGeCawPFDgB+bWR7wLHCUmT3VxLWuwFvosOH2wZ4/E4kKZWXem11trffmd8wxfkckLfHHP0JODqxYAVdf7Xc0IuKjSJ0ldDDwVqN9bwC5ZpbsnPs/51yWcy4bOB141zl3VhPX+gswoNF2WJtELZHnuuu8N7ucHO/NT6JLair885+QkgIPPABvvOF3RCLik0hNWHoDGxrtyweSgR4tuZBzrsQ5lwf8HLWwxJe5c+G++6gx49iCArIGDyYrK2uHraCgwO8oZXdGjIA//MF7fPHF3kKVIhJ3IjVhaTbn3HvOuUnNOG+ac86cc4bXyiKxrLp6e7XUB9LSqNl336ADPGtra/2OVJrjyith5EjIy/Mq4opI3InUhCUfyGi0LwOoATaGPxyJOnffDYsWQXY2d3Xs6Hc0sqeSk70uITP4859h8WK/IxKRMIvUhOVj4NhG+yYC85xz1a25oJlNC0yjdnjdQhKr8vLgxhu9x/ffT0VCpP43lxYZNw4uucSrqXPRRVCnigYi8cSPOizpZjbKzEYFdvUPfD+4wWl3AYeb2VQzG2pmlwBnALe39r7qEooTzsGll3rjHE47DX70I78jklC6+WbIzIRPP4UHH/Q7GhEJIz8+eo4FPg9sALcGHj9cf4Jzbi4wGfgZsBi4Gvilc25Wa2+qFpY4MXMmvP46dO4Md93ldzQSap07wz33eI+vu85bxFJE4oIfdVjeq2/paLSNb3TeS865HOdcinNukHPukT28r1pYYtz40aPJP+00AK4zIys3VzOBYtHkyXDCCbB5szcYV0TiQtSuJSTS2M9XrKB3XR1fdurEvJEjGRIov681gaJLU+sMZWZmMn/+fG/g7b33wn//C889B+eeq64/kTgQNwmLmU0DbvQ7Dmkjn37KOWVl1Jjx5733xmmtoKhVv85QYytWrPjhm+xsmD4drrnGG4j75ZfQoUP4ghSRsIub6RPqEoph1dXwy1+SADyflcXq9HS/I5Jw+M1vVJtFJI7ETcIiMexvf4NFi1iTmMg/srP9jkbCpXFtluXL/Y5IRNpQ3HQJSWzIzc3dYUxKl7o6PszPpwtwFVCZmOhXaOKHcePgF7+Ahx+Ga6+Fl17yOyIRaSNx08Kiac2xYf369TuU1r8xMZEuzvF5ly78W4XE4tMf/gBpafDyy/Duu35HIyJtJG4SFo1hiT19t23j5O+/pw64f9Agr2tA4k/v3vB//+c9vvpq0PpQIjEpbhIWiT0XrVpFsnO82bs3K7VeUHy76irYay/44gv4xz/8jkZE2oASFolKI0pKOHzjRsoTEnhkgBrM4l779nDLLd7j66+H0lJ/4xGRkIubhEVjWGKHOcclK1cC8Gy/fmxq187niCQiTJkCubleuf7bW73smIhEqLhJWDSGJXrkBkrqB9sKCgo4ZsMGhpaWUpiSwvN77eV3uBIpEhLgzju9x7ffDmvX+huPiIRUs6Y1m1ktcJBzbl6QY/sD85xzmk8qIVE/EyiYknXruHDVKgAeHjiQCk1jloYOPRR+8hN44QWva+iJJ/yOSERCpLktLLuafpEEaD6phMVVQM+qKr5OT+ftjAy/w5Ewql9jKNiWm5v7w4m33gopKd7g2//9z7+ARSSkmmxhMbN0oFODXT3MrE+j01KB0wC1vUqb615Zye+cA+D+wYO1XlCcaWqNIWi0ztCgQXD55XDHHd405/fe05R3kRiwqxaWq/ESke8AB7wSeNxwWwH8BniwbcMUgfNXryYNmNOjBwu7dPE7HIlk118P3bvDnDkwa5bf0YhICOxqDMssIA+vO+hRYAbwTaNzKoGvnHML2yK4UNJqzdGtf1kZE/PzqQYeHDjQ73Ak0nXp4q3m/Otfw3XXwYknQpJWIhGJZk3+BQeSkIUAganArznnNoYrsFBzzk0DpgGYWTaa2uy7xusC1SsoKNip6f/neXkkAH8Hvu/QITwBSnS76CJv1tDXX8PTT8M55/gdkYjsgWYNunXOPRHNyYpEpsbrAtVvtY1Kqw8qLWV8YSFVZtyisQjSXMnJMHWq93j6dKiu9jceEdkjzUpYzCzFzKaa2RIz22JmVY22yrYOVOLXz1d7jWEv9+nDOiUs0hJnngl77w2rVmmKs0iUa26n7p+Bi4HXgJfxxq6ItLm9t2zh0E2bqEhI4Ol+/bwqpiLNlZQEN97oJS5/+AOcfTaoMrJIVGpuwvIT4Hrn3G1tGYxIY+fl5QEwq29fivVGI61x2mlw882wdCk88ghcconfEYlIKzS3cFwqoApMElb7bt7MgUVFbEtM5BmV4JfWSkz0xrCAl7iUl/sbj4i0SnMTln8Bx7VlIBKbdrcu0K6cHxi7MrNvX7akpIQjXIlVp54KI0fCunXwwAN+RyMirdDcLqEXgXvNrAPwBlDc+ATn3MehDCzUVIfFH7taFyjYlOZ6I0tK2L+khNLERC1wKLtVX7Y/mMzMTObPnw833QQnnQS33AIXXghpaWGOUkT2RHMTltcDXy8NbK7BMQt8H9Gr0KkOSxRxjvMCrSv/2msvSpOTfQ5IIl2zyvafeCLk5sL8+XD//XDNNWGMUET2VHMTliPbNAqRBo4GRm7ezJakJGY28alZpMXMvFaWH/0IbrsNfvUr6NjR76iklRYsWDAhKSnpRudcb5o/vEEiV52Z5dfU1EwfM2bMm8FOaFbC4px7P7RxiTTBOaYHFjh8bq+9KFM5dQmlCRPg4IPh44/hnnu8NYck6ixYsGBCu3bt7s3Ozq5q3759cUJCgtv9T0kkq6urs/Ly8s55eXn3Lliw4NfBkhZlpRJRxhUVcRBQnJzMi337+h2OxBozrx4LeKs5l5T4Go60TlJS0o3Z2dlVaWlp5UpWYkNCQoJLS0srz87OrkpKSgo63rRZH1/NrJodx63sxDmnaRxxqqk1gSD4ukBNco5zA3VXnunXjwq1rkhbOOooGD8e3nvPa2WpL98vUcM517t9+/Y7Tf6Q6Ne+ffuKQDffTpr7jnAzOycs3fGmOrcDHm91dBL1WjsTqLFRJSXss3UrhXhl+EXazNSpPyQsV1+tGUPRJ0EtK7Ep8HsN2vvT3DEs04LtN7NE4BWgtLXBidQ7Y80aAP5qRmViRE86k2g3fjwccADMmwePPgqXXeZ3RCKyG3s0hsU5Vwv8DfhNaMJpHjNLNbN5ZrbQzL40s+nhvL+E3pCtW8ktLqY8IYH7/Q5GYp8ZXHed9/iOO7SSs0gUCMWg2xSgUwiu0xKVwFHOuZHAKGCimR0Y5hgkhKYEWlde6dOHYq3ILOFw0kkwdCisWQPPPut3NCIx68EHH+xqZvvv6XWalbCYWb8g22AzOxm4FfhoTwNpCeep74ZKDmzqz4xSfbdt4/DCQqrN+Jfqrki4JCTA737nPb7tNqir8zceiQuTJ0/ONrP9G28ZGRkj/I6toQEDBux31VVXRdRgwuYOus0jeEJgwIfAL5t7QzM7HLgG2B/IBKY453b6eGNmP8Yb7Ls3sAa4xTn3aIPjicBnwGDgPufc3ObGIC0XsplAQZz23XckAm9mZLAxNbXV1xEJpqmy/ZmZmcz/6CP4f/8PvvwSXnvNq4Yr0sbGjRu39fnnn1/VcF+SZkXuVnO7hI4Ejmq0HQz0dc4d7pxb24J7pgEL8Ur8B2Vm4/DWL5oJjATuBh4MtOgA3vgZ59woIAs4wMxyWhCDtFD9TKBgW21tbauv262ykgn5+dQBz2rNIGkD9WX7G2/r16+HlBRvlhB4aww5NdRK20tOTnb9+vWrabj16dOnZsmSJe3S0tJG33bbbT3rz/3www87pKSkjHn44Ye7Alx11VV9BgwYsN99993XrW/fvsPbtWs35tBDDx2yYsWKHUqLvPDCC51GjBgxLDU1dUxmZubwc845p19xcfH29/y6ujpmzJjRa+DAgfulpKSM6dGjx8jJkydnAxxwwAFD8/LyUu+6667M+hag5cuXpwAsWrSo3dFHHz2oY8eOo7p06TLqqKOOGvzll1+2a3jvadOmZfTq1WtE+/btRx977LGDCgsLQ7K+Stgr3TrnZgOzAazpsQpXAh80mJ20LJDEXAPManS9EjP7LzARWNL4QmbWBejSaLf6HSLET9auJcU55vTowXeaWip+uPBCmDEDPvkEPvwQDjvM74ikNUIwRqJVnPssVJfKycmpvO2229b89re/7X/UUUdt3XvvvavOOuusgaeccsqmCy64YHvdmfz8/JTHHnus57PPPvtNTU0Nl112Wf/JkycPXLRo0TKAWbNmdTz77LMHTZ8+fe3xxx+/ev369UlXXnllv7POOiv7tddeWwVwxRVX9Hn00Uczbrzxxu8mTpy4tbi4OPGtt97qCPDKK6+sPOCAA/aZMGFCyQ033JAP0KdPn5o1a9YkjR8/fthJJ51UdOutty5LSEhgxowZvY877ri9ly9fviQ9Pd09/vjjXWbMmNF3+vTpa0888cTNr7zySuc//elPIelaalEblJkNBQ4DugJFeEnF16EIpJGD8WYfNfQG8JiZJeMlINWBZKU9cCxwWxPXugKt0hyROjvHj9etA+Dpfv18jkbiVnq6N615+nS49VYlLNLmPvnkk04dOnQY3XDfcccdVzJr1qzVv/71rze9/fbbnaZMmTJo33333ZaQkOAefvjh7xqeW1lZmfDUU0+tHjZsWBXAE088sXrs2LH7vfnmm+kTJkwovfnmm/tcdNFFG6699tpC8BKh++6779vx48fvs27duqS0tLS6Bx98sPcNN9yw9pprrtlYf93DDz98G0BGRkZtQkKCS09Pr+vXr19N/fG77rqr14ABAyoee+yx7fE8//zzed26dRv1/PPPdzn//POL77nnnt4//vGPi66//voCgBEjRhT873//6zBr1qzue/rv1txKt8nAo8AZeONW6jkzewr4hXOuJugPt05vYEOjffl4g2t7AD2BJwLjWBKA551zrzZxrb+wc2G7LOCDUAUrrfMrIK22ls+6dGF5p3BPNBNp4Ne/httvh9dfh0WLYEREjX+U5ghhS0dbGzVqVOnjjz+e13Bfp06dtvetP/7449/us88+ObNnz+46Z86crzp27LjDiPAePXpU1ycrAPvvv39Fx44daxcvXpw6YcKE0iVLlnRYtGhR2kMPPZRRf44LdHcuXbq0XWJioqusrLTjjz9+S0vi/vzzzzssXLgwrXGyVVFRkbBixYp2ACtXrkydMmXKxobHDz744NKwJSzANOBnwPXAM3jJQ2/gdOAmvEGx/29Pg2ku59wiYPRuT/TOLQFKzGwaamnZraYG1+7pwNrGUmpr+U3gD0itK+K7Hj28rqG77/ZmDP3zn35HJDGsffv2dTk5OZVNHV+0aFHqpk2bkpxztnLlynYHHXRQeUuu75yzyy+/fP25555b1PhYdnZ21fz589u3Ju66ujo77LDDtvz1r3/9rvGxnj17hrLRIqjmJixnAdOdc7c22PctcFugleMiQpuw5AMZjfZlADXAxp1P373AeJhpAGaWDaxudXQxrKky+y0psd8cE/PzyQCWp6ezoGvXkF5bpFWuugruu8+ryfKHP8DAgX5HJHGouLg44eyzzx7405/+dFO/fv0qL7vssuyDDz74y0GDBm2vbrhx48bk5cuXpwwdOrQKYMGCBalbt25NzMnJqQDYb7/9ypYtW5baVFI0ZsyYipSUFPf66693GjNmTEWwc5KTk13jCRWjRo0qmzlzZvfBgwdXpaamBh2hPnjw4IpPPvkknQbv1YHv91hzZwn1BuY1cWxe4HgofYw3LqWhicA851yrSlKa2TQzc2bmULLiq4S6Ok77zkvQn+7Xz6s6KuK3fv3gzDO9eix//rPf0UgMq66utjVr1iQ13gDOP//8fsnJye6BBx747uabb84fOnTotjPOOGNAw+ShXbt2dWefffaADz74oMP777/f4dxzzx2w3377bZs4cWIpwNSpU9fNnj2768UXX9x37ty57RcvXtzu6aef7jxlypT+AJ07d6678MILN9x66619b7/99h5Llixp99FHH7WfOnXq9oaCvfbaq3Lu3LnpK1asSFm/fn1SbW0t11xzTUFVVZX96Ec/GvT222+nLVu2LOW1115Lv+CCC7IWL17cDuCyyy7b8PLLL3f74x//2HPx4sXt/vjHP/Z89913u4Ti3625Ccsa4Jgmjh0D7NQ81BQzSzezUWY2KrCrf+D7wQ1Ouws43MymmtlQM7sEb/zM7c29T2POuWnOOXPOGTCgtdeRPXdEYSF9Kir4GviwZ8/dni8SNtde63199FEoLPQ3FolZc+fO7di/f/+Rjbe//vWv3V9//fVuTz311KqOHTvWJSYm8uyzz67++uuv2//+97/f3jDQu3fvqrPPPnvjaaedNui4444b1rFjx5oXX3zxm/rjJ5100tZZs2at+Oyzz9LHjx8/bNy4cftOmzatb+/evbd/4L/nnnu+v+qqq9bde++9vUePHr3fySefPOSrr77a3lV00003rSstLU0cPnx4Tp8+fUauXLkyZa+99qr54IMPlqWlpdX99Kc/HTJy5MicX/3qV9lbt25N7NGjRy3AeeedV/z73//++7vuuivzgAMO2Pedd97pfM0116wLxb9bc7uEHgJuMbMO7DyG5VLguhbccyzw3wbf3xrY3gfGAzjn5prZZLzCcTfgJUS/dM7NasF9JEL9dK1XtucuM+rUuiKRZN99YdIkePVVeOABuOEGvyOSGDNz5sw8vGKsQV122WWbGn4/YMCA6uLi4oXBzmt8bkOTJk3aOmnSpOVNHU9ISOCmm27acNNNNzWe4ALAIYccUr5kyZKvGu8fNmxY1auvvroq2M/UC3bd+llDe6K5dVj+ZGbd8RY5bFjwrRq4wzl3R3Nv6Jx7jx1nGjV13kvAS8297u5o0G1k2GfLFvbZupUtSUk8WVtLrt8BSdxqqgLuoRUVPAtw//1e6f7kkNS8EpE91Ow6LM6535nZbcCBeHVYioFPnXM7jUKORBp0GxkmB1pXXsvMpPz7732ORuJZfQXcxjY4x7KPPmLY+vXwwgswZYoP0YlIYy1ardk5V+Sce90598/A16hIViQy9Kis5IjCQmqBWX37+h2OSHBmPJYemNRw993+xiLSyJ133rlu9erVX/odhx+aTFjMbIyZbTKzJlcDM7NJZrbRzIa3TXiho1lCP8jNzSUrKyvoVlCwx92MTfrx99+T5Bwf9OxJgRY5lAg2s3176NoV5s71NhHx3a5aWC4H5jvnXmnqhEB12U/x1v6JaJol9IO2WshwV5JrazkxUMvlRbWuSISrSEjwCskB3HOPv8GICLDrhOVI4MlmXOMZvNWbRZp0dEEBXaqr+To9ncWdO/sdjsjuXXopJCbC88/DupDMyhSRPbCrhKU3XjXb3fkWyAxNOBKTnNs+2PbFrCwVipPo0K8fnHIK1NTA3xqvxSoi4barhKUUaM5iRd0D50Y0jWHxz8jNmxlcVkZRcjLv9urldzgizfeb33hfH3gAKoJWMBeRMNlVwvIZcGozrjE5cG5E0xgW/5waaF15tU8fqhNaNDFNxF+HHAKjR3tVb595xu9oROLaruqwPAg8Z2YfO+ceCHaCmV2MVzL/tLYITlqvqVWXIfQrL+9Kf+c4ZONGqs14qU+fsNxTZE81LCr3k7Iy/gIs+eUvmXjDDWT26cP8+fN9jU+Cy8nJ2aewsDBslf569uxZHawarJ+WL1+eMmzYsOHvvPPOsqOPPros2Dlr1qxJ6t+//8hXXnnl60mTJm0N1b1fffXVjieeeOLe33777cJ+/fqFfPXmJhMW59wLZvYI8Dcz+yXwGt6aQgD9gOOBUcAjzrmZoQ5M9kxTqy7XHwuXS5wjEXi3Z0+K2rUL231F9kTDonJL6+oo+uQTcqqrmdyrFzPD+PcjLVNYWJicnZ0dtr67vLy8FtdnmDx5cvaLL77YHbzy+D169Kg++OCDt9x5553fDxgwoFWL+zY0aNCgqm+//XZhRkZG20z59NEu2+edcxcBFwIdgOuBBwLb9UAacFHgnIinMSzhl1pbyy8Cj2cGKYEuEg2qExJ4JdA6WD94XGRPjBs3buu33367cOXKlYseffTRVV9++WWHU089dVAorp2UlES/fv1q2rVr50JxvUiy2wEFzrlHnHPDgCy8svwHAns554Y55x5p6wBDRWNYwu+4/Hy6Aks6dWJ5p05+hyPSai/36UO1GYds3EhWTchbuiXOJCcnu379+tUMGDCg+oQTTig999xzN37xxRdpRUVFCQAvvPBCpxEjRgxLTU0dk5mZOfycc87pV1xcvP39evbs2emjRo0a1qFDh9EdO3YclZOTs8+cOXM6gNclZGb7/+c//0mrP//f//53p8GDB+/Xrl27MTk5Oft8/PHHaQ3jefXVVzua2f5r1qzZodfFzPZ/8MEHu9Z/f9lll/UdOHDgfu3btx/du3fvEWeccUa/TZs2JTb1PDdt2pR4yimnZHfr1m1ku3btxvTt23d4w1WnW6olawmtA1SMQJrHOU4NrBWkQnES7YrateO9nj05tqCAc8uCDgsQaZW8vLzkWbNmdU1MTCQpKYlZs2Z1PPvsswdNnz597fHHH796/fr1SVdeeWW/s846K/u1115bVV1dzWmnnTb49NNP3/jUU0+trqurY+7cuR2Sk5ODtqisXr06ecqUKYNPOeWUTc8///w3y5cvb/e73/1ur9bE2r59+7r777//2wEDBlQtX7683ZVXXtnvwgsv3OvFF1/MC3b+VVdd1eerr77qMHPmzBV9+/atWbFiRcp3332X0pp7QwsSFpGWGF1SQv9t2/gemNOzp9/hiOyxmVlZHFtQwJSyMm+Ks5aXkFb65JNPOnXo0GG0c46KiooEgPPOO6+gU6dOdTfffHOfiy66aMO1115bCJCTk1N53333fTt+/Ph91q1bl5SYmOi2bt2aeNJJJ5Xk5ORUAowYMaKyqXvdeeedvXr06FH95JNPfpuUlMSYMWMqNmzYkHTllVdmtzTuP/3pT9sHcA0dOrSqtLR07c9//vNBtbW1eYmJOze0fPfdd+1ycnK2HXnkkdsA9t5776qW3rMhzTGVNnFioDLow0CtpjJLDFjeqRPL09Pp4py3irNIK40aNap03rx5S+fMmfPV5Zdfvn748OFl99xzz1qAJUuWdPj73//eu0OHDqPrt4kTJw4FWLp0abuMjIzayZMnbzr11FP3Hj9+/OCpU6dmLF++vMlWi2XLlqWOHj26NCnph/aJI488slW105544okuY8eOHdqrV68RHTp0GP2LX/xiYHV1tX333XdBZ2b96le/Knjttde67r333vtecMEFWS+//HLH1ty3nt5JJOS6VFVx6MaN1AKPqKqtxJD6wbdzf/GLoIuH5ubm+hyhRIP27dvX5eTkVObm5lbcfffd6/r161d13nnn9Qdwztnll1++ft68eUvrt/nz5y9dvHjxkgMPPHAbwAsvvJA3Z86crw499NCtb775ZucRI0bkPPXUU11aG09iYqKrv3f9vsrKyh1evN999920888/f9Chhx669bnnnlv56aefLr3zzju/DXZuvZ/97GdbVq1atfiKK67ILykpSZoyZcrgSZMmDWxtnHHTJWRm04Ab/Y4jHkzMzyfZOT7q3p3vi4sJydB3kQjwn169uOjrrxlXVcWxffqQl7bD2EVWrFjhU2QSzf7whz+sGzt27H5z5swp2G+//cqWLVuWWt/d05Rx48aVjxs3rnzGjBkbjjnmmEGPPfZY97POOquk8XnDhg2r+Pe//92tpqaG+laW9957L73hOb17964BWLNmTXL//v2rAT766KMODc95//3307t37179l7/8ZftY1ueee64ru5GZmVlzySWXFF1yySVFTz75ZMk555wzqLCwMLFnz54tnnYdNy0smiUUHuYckwLdQa+qUJzEmIqkJJ4OPJ6kBRElREaPHl0xfvz4kt///vd9p06dum727NldL7744r5z585tv3jx4nZPP/105ylTpvQHWLJkSbtLL72079tvv522YsWKlNdeey196dKlHYYOHRq0/swVV1xRUFhYmHzOOef0W7BgQepzzz3X+a677tph/b+cnJyKPn36VE2dOrXP559/nvrGG2+kX3PNNTvUohg2bFjFpk2bku+5557uS5cuTbn33nu7P/zww7tca+Xyyy/v8+STT3ZZvHhxu88//zx15syZXXv27FndvXv3VtWIaVYLi5lN3cXhOmAz8IVz7oPWBCGt11RF23BWs21odHExfSsqyG/XjnnduoX9/iJt7SEzfuUcx23YwIMDB1IVZLCh+KNnz57VrSnmtif3C9W1fve73+VPmDBh2HXXXZc/a9asFTfddFOfxx9/vJeZkZWVVXnCCSeUAHTs2LHu66+/Tj3zzDMHlZSUJHXr1q3m+OOPL77jjjuCZtCDBg2qfvrpp1dee+21ex144IE9Bg8eXH7LLbd8d8YZZwyuPyc5OZknn3xy1WWXXdbvoIMO2jc7O7vinnvu+XbChAnD6s+ZMmXK5k8//XT99OnTs373u98ljB07tnTGjBnfXXTRRU128aSkpLjp06f3/f7771OSk5Pd8OHDy1566aUVCa0c19jcLqH/F/ga7C+zNrDfmdn/gOOdc5taFY20WFMVbcNZzbahEwP3fT0zkzqNX5EY9IUZX6Wns8/WrRxZWMibvVtdVkJCLNLK5Aczc+bMvGD7jzvuuDLn3PZ1+SZNmrQ82Hn9+/evfvvtt79p6vpDhw6tangdgMmTJ2+ZPHnylw33TZkyZYdzjjrqqLIvv/xyh3+/xte5++6719199907JEYXXnhhw5i3NvyZO+64Y/0dd9wRsjej5qY5Y4BvgIuBPkBy4OulwCrgUOBkIBv4U6iCk+jStbJy+2Db2ZmZuz1fJFrVd3eqW0gkfJqbsNwP3Ouce8A5l++cqw18/RtwH/An59wrwB/x1hiSOPSj/HySnOOT7t3ZqHWDJIa926sXpYmJ5GzZwoDSVs0QFZEWam7Csj/wdRPHluG1wAAsBDRwIcRyc3ODTqHMysqioKDA7/AAb7DtCYHuoFc02FZiXEViIu9kZAAwSYshioRFc8ew5ANnAG8FOXZW4DhAZ6AoBHFJA5Gy8vKujC0upk9gsO3/NNhW4sCrffpw8rp1HJefz4MDB1Kpwbcibaq5CcutwN/NbG9gFlAI9MQbt3IgUL9i89HAvNCGGBqqw9K2Gk5l1mBbiQffpKeztGNH9t26lSMLCnhD47bCqa6urs4SEhJibkXieFdXV2d4s4930qyExTn3oJnlA9cDMwI/VwssAE4KjF8BuAEI2TSvUHLOTQOmAZhZNrDax3BiSm/nOGTjRmrMmK0ZExJHXunTh32XL+fEdet4IzOTgoICsrKydjovMzOT+fPn+xBhbDKz/PLy8s5paWnlfscioVVeXp4ayDd20pLVml8GXjazBLzWlY3OudpG52zZo0glKp2HN699TvfuFGmwrcSR//bqxaUrV7Lv1q0MKi1lTm1t0O5bVcANrZqamul5eXn3Zmdn0759+wq1tES/uro6Ky8vT83Ly0upqamZHuyc1pTmT8JrrknEa2WROJbgHBc477VCg20l3lQmJvJW796c+v33nLhuHY/5HVCcGDNmzJsLFiz49TfffHOjc643cVS1PYbVmVl+TU3N9DFjxrwZ7IRmJyxmdiwwHcjF+89RZ2bzgBudc++EJFyJOmOLisgG1qWm8lnX3S4rIRJzXs3M5NTvv+eYDRtIc/qgHy6BN7Wgb2wSm5qVlZrZRGA23iygGcAlwM1AF2C2mU1oqwAlstVP6XwtMxOnwbYSh1anp7OkUyfSams53e9gRGJYc5vRpgOvAznOuemBAnLTgBzgjcBxiTNdq6o4eONGakCDbSWuvRqYIXSuWlhE2kxzu4SGA9Oc2/Gv0TnnzOxvwAshjyzONLWIIfi3kOHuHL1hA4nAK0CxBttKHHu/Z09+s2IFh9TVkbVtG2s7dNjheFOzh0AziESaq7kJSznQo4ljPQLHw8rM9gL+AWQADnjQOXd3uOMIlWgoDrcD55iY7808e0JdQRLnKpKSeK9XL36Un8+E/HweGbjjAra1TcweAs0gEmmu5nYJvQXcYmZjGu40s9F4Y1neCHVgzVADXO2c2xeveN2lZravD3HEpcGlpQwqK2NzUhKv+h2MSASoX7X5uA0bSFDXkEjINTdh+S1QAcw3s+Vm9r6ZLQf+F9h/TVsF2BTn3Hrn3ILA463AV0DfcMcRr+pbV/6TkUG1WlhEWNS5M6uAXpWVjC4u9jsckZjTrITFOfc9MAq4Am+Bw6rA198Ao5xzLVpj3cwON7NXzGydmTkzCzq43sx+bGaLzazSzFaY2flNnJcNjAbmtiQOaZ2kujqO2bABgDc02FYEAGfGPwLJe31CLyKh05JKt6XAXwPbdma2n5md6Zz7fQvum4aX8DwKvBjsBDMbFzg2A3gWOAZ40MyKnHOzGpyXDswErlCl3fA4cNMmOtfU8E1aGivS0/0ORyRiPIm3/sdhGzeSVlNDWdLuX2I1IFekeVpT6baxYcDvgGYnLM652Xh1XbCmuxOuBD4ITJ8GWBZIYq7BW4ARM0vGS1b+6ZxrKvHpglcvpqHgrw7SLPWfHt/s3RvUHSSyXZ4Zn3fuzOiSEsYXFPBaM6o/a0CuSPOEImFpKwcDf2u07w3gsUCiUgM8AnzlnLtzF9e5gghapbmp6cuROnW5sa5VVRy4aRO1wDsZGX6HIxJx3uzdm9ElJUzMz29WwiIizRPJCUtvYEOjfflAMt5U6kHA2cBiM/sicPz3zrnXG/3MX4DHG+3LAj4IYazN1tT05YicuhxEfe2Vj7p3pzglxe9wRCLO+z168JuEBHK2bAlak0VEWidqF4xyzn3onDPn3Ajn3KjA1jhZwTlX4pzLA34OrA5sviQrUa9B7ZU3NdhWJKj6miwAEzT4ViRkIjlhyccrCtdQBl5X0MaWXsw5Ny2Q4BgwIATxxZ2GtVc+6d7d73BEIlaoarLUD8htvOXm5oYqVJGo0WSXkJlV41WQ3Z22GnX5MXAscEuDfROBec656pZezMymEUFjWaJRw9orNQmRnOuK+GtR586sS02lT0UFo4uLea+V12lqQK4G40o82tUYlptpXsLSYoGpyIMb7OpvZqOAUufcysC+u4CPzGwq8BxwNHAG8JPW3DMw22ha4P7ZeF1D0kzJzqn2ikgzOTPe7N2b8/LymJifz11+ByQSA5pMWBpMJ24LY4H/Nvj+1sD2PjA+cP+5ZjYZL3G6AfgO+GXDGiwSPieAaq+ItMBbGRmcl5fHYRs30kml+kX2mC+zhJxz79GMriTn3EvAS6G4p7qE9sy5gRdc1V4RaZ789u35vEsXRpeU8DMzvvY7IJEoFzcDETTotvW6VlVxPKj2ikgL1Q++PUctLCJ7LG4SFmm9ozZsIAn4VLVXRFrk/R49KE9I4BAga9s2v8MRiWpxk7CY2bTAQosODbhtkWMKCgCvT15Emq8iKYk5PXsCcFTg7ygUmprurCnPEssiudJtSGmWUOtkbdvGsK1b2YzXwiIiLfNORgYTNmzg6A0b+Ef//iEZA6b1hyQexU0Li7RO/afCfwNViYn+BiMShRZ06cIGoF95OXuXlvodjkjUipsWFs0SagXnODpQe+UZzQwSaZW6hASeBy7DW4vr644d2/R+9d1FwWRmZjJ//vw2vb9IW4mbhEVdQi23d2kp/crLKUpO5r81NRzid0AiUeoZMy5zjiMLCnhg0CDq2vADgLqLJFapS0iaVN+68t9evahVC4tIq80F1qWm0rOqihElJX6HIxKVlLBIUAmBT4PgrR0kInvAjP8EVnA+OoSzhVpKiylKNIubLiFpmRElJfSsquL71FS+auM+d5F48J+MDM5es4YjCgu5Z8gQqn1YQFSLKUo0i5uEJZyDbnNzc1m/fn3QYwUFBU32L0eSoxu2rqg7SGSPfZuWxsq0NAaXlXFAUREf9ejhd0giUSVuEpZwDrpdv359k0lJU4lMJEmuq+OIwkKA7c3YIrLn/pORweBVqzh6wwYlLCItpDEsspMDioroWFPDivR01qSl+R2OSMx4N/AB4OBNm2hfU+NzNCLRJW5aWKT56mcHqXVFJLQKUlNZ2LkzIzdv5tCNG3k7sDii31S7RaKBEhbZQbpzHLxpE/DDp0ERCZ3/9OrFyM2bObqgIGISFtVukWigLiHZwclAu7o6FnbuTGFqqt/hiMSc93v2pMaMsUVFdKmq8jsckagRNwmLVmtuninOAaq9ItJWtqSkML9rVxKB8YHB7SKye3GTsDjnpjnnzDlnwAC/44lEXaqqOAaoMeP9nj39DkckZtV/IKgfLyYiuxc3CYvs3vjCQpKAed26sSU52e9wRGLWR927U56QQM6WLWQHWjVFZNeUsMh2mh0kEh4VSUl8HKjDcprPsYhECyUsAkBGeTk5W7ZQBttfSEWk7bwT+GAwRS0sIs2ihEWAHwb/vQJUJCb6G4xIHPhft25sSUpiONCvrMzvcEQinhIWAdheiv8FrRskEhY1CQnby/MfodlCIrulhEXIqKhgn61bKU9IYLbfwYjEkfcCs/GUsIjsXtwkLKrD0rTDAy+Wn3bvToVaWETCZkHXrhQDg8rK2GvbNr/DEYlocZOwqA5L08YXFAA/fNoTkfCoSUjg5cDjIwJ/hyISnNYSinO9KirYN9AdNLd7d7/DEYk7L5hxrnMcUVjIU9nZfoezEy2MKJFCCUucq+8Omtu9O5WaHSQSdu8ApYmJDC4ro++2bXzfoYPfIe1ACyNKpIibLiEJrn6wn0rxi/ijymx77SOtLSTSNCUscaxnRQU5W7ZQkZDAp926+R2OSNzSbCGR3VPCEscO37gR8NYOqkhS76CIX/7XtStliYkMKS2lj2YLiQSld6k4doRmB4lEhOrERD7u3p1jCwo4orCQZ/r39zukZtGAXAmnqE1YzOxRYBJQ4JzL8TueaNOjspLhW7ZQmZDAp5odJOK793v14tiCAsZHUcKiAbkSTtHcJfQ4MNHvIKJV/eyged26Ua7uIBHfze/alW2JiexdWkpmebnf4YhEnKhNWJxzc4Aiv+OIVodrdpBIRKlKTOSTQGunBt+K7CzsCYuZHW5mr5jZukCp/NObOO/HZrbYzCrNbIWZnR/uWHclNzeXrKysoFtBhFes7O0cwzdv9qZTqjtIJGJotpBI0/zoC0gDFgKPAi8GO8HMxgWOzQCeBY4BHjSzIufcrJbczMy6AF0a7Q4+SqwF1q9f32Tf7fr16/f08m3qVLxMVd1BIpFlXrdulCckMGzrVnqrW0hkB2F/t3LOzQZvUWBreqG9K4EPnHPTAt8vCyQx1wCzWnjLK4AbWxpnLPuJc4C6g0QiTX230FGFhRxRWMjzfgckEkEidQzLwcBbjfa9AeSaWXILr/UXvMUOG26H7WmA0aprZSWH4VXX/CRQXVNEIsf7vXoB6hYSaSxSE5bewIZG+/KBZKAHgJk9A3wCDDWztWb2i2AXcs6VOOfygJ8DqwPbB20TduQ7fONGEoD53bpRpu4gkYgzN9AttM/WrfQLtIaKSOQmLLvlnJvinMt0ziU757Kcc4/s5vxpzjlzzhleK0tcqv/UNkfdQSIRqTIxcfvK6ZN9jkUkkkRqwpIPZDTalwHUABvDH05s6FxVxYiSEqqAjzQ7SCRi1c8WOkUtLCLbRWrC8jFwbKN9E4F5zrnq1lzQzKYFplE7vG6huHPQpk0kAv8FypJbOhRIRMJlXrduVJlxEN64MxHxpw5LupmNMrNRgV39A98PbnDaXcDhZjbVzIaa2SXAGcDtrb2vuoTg0MBihy81PTtLRCJAeVISn3XtSgJw8KZNfocjEhH8aGEZC3we2ABuDTx+uP4E59xcvO7bnwGLgauBX7a0BktD8d7Cklpby9jiYgBe9jkWEdm9DwOz+Oo/aIjEu7AnLM659+pbOhpt4xud95JzLsc5l+KcG7S7QbXNuG9ct7CMLSqiXV0dSzt2ZL1aWEQi3ic9elAHjCkupn1Njd/hiPguUsewSIjVf0r7ULVXRKJCcUoKHwMpznFAkZZNE4mbhCWeu4QS6uo4KNAProRFJHrUjzdTt5BIHCUs8dwlNGLzZjrV1LCmfXu+S0vzOxwRaab68WYHbtpEUl2dr7GI+C1uEpZ4pu4gkei00ozVHTqQXlvLyJISv8MR8ZUSlljnHIcoYRGJWpotJOKJm4QlXsewDC4tpXdlJZtSUljWqZPf4YhIC9UnLIds3AiqfCtxLG4Slngdw1LfuvJR9+44TWcWiTpfd+xIYUoKPauq2HvrVr/DEfFN3CQs8UrjV0SinJm6hURQwhLTepeXM7isjLLERL7o2tXvcESklZSwiChhiWn13UFzu3WjOkG/apFotbBLF0oTExmwbRt9t23zOxwRX8TNu1g8DrpVd5BIbKhNSOCT7t2BHz6IiMSbuElY4m3QbXfnGL55M9VmzAu80IlI9FK3kMS7uElY4s0kIBH4vEsXypKS/A5HRPbQ/G7dqDJjvy1b6KXpzRKHlLDEqJMCL2jqDhKJDeVJSXzWtSsJwIl+ByPiAyUsMahdbS3HBh5/rIRFJGbUfwA5SS0sEofiJmGJp0G3Y4uK6AAs7diRTe3a+R2OiITIJz16UAccDbSvqfE7HJGwipuEJZ4G3Wp2kEhsKk5J4ctOnUgFcouK/A5HJKziJmGJFwnOMS7wQvaJEhaRmFPfzXvQpk0+RyISXkpYYszQLVvoWl3NaiCvQwe/wxGREPs0UKZgXFERprEsEkeUsMSYAwOtK68BaLFDkZiT16EDq4Gu1dUM02KIEkeUsMSY+mbi15SsiMQmM+8DCXCguoUkjihhiSE9KioYUlpKeUIC7/sdjIi0mfoPJBrHIvFECUsMOSjQHfS/bt2oVAuLSMx6HyhPSGBIaSk9Kir8DkckLOImYYmHOiz1zcOfduvmcyQi0pYqzfisa1fghw8qIrEubhKWWK/DklJby5jiYgDmarFDkZhXv3rzOHULSZyIm4Ql1o0uKSG1ro7l6emqbisSB+o/mOxfXExKba3P0Yi0PSUsMaJ+8N0nal0RiQub2rVjeXo6qXV1jC4p8TsckTanhCUWOPfD+BUlLCJxo/7vXbOFJB4oYYkBA8rKyKispCg5ma87dvQ7HBEJkx3GsajqrcQ4JSwx4KAGrStO05lF4sbXHTtSlJxM78pKBpSV+R2OSJtSwhIDNH5FJD45M3ULSdxQwhLlujvHvlu2UNWgLoOIxI9PlLBInIjahMXMJprZcjNbaWbX+R2PX36E90tc2KUL5UlJfocjImH2WdeuVJux75YtdKqq8jsckTYTlQmLmSUC9+G9X+8LTDGzff2Nyh/HBwbaaXaQSHwqT0riiy5dSAAOUNVbiWFRmbAABwArnXOrnHNVwLPAST7HFHaJdXVMCDxWwiISvzSOReJB2BMWMzvczF4xs3WBtX1Ob+K8H5vZYjOrNLMVZnZ+g8N9ge8afL82sC+u5GzeTBcgr0MH1rVv73c4IuKT+oTlgKIikjS9WWKUH4Me0oCFwKPAi8FOMLNxgWMz8FpPjgEeNLMi59ysltzMzLoAXRrt7g+wdu3allxqBzU1NZSXlwc95pxr8bHW/MygDRvIA17p3Hmn4625Xqjji5R7Rfr1FHtsXM/P2L8BPkxNJauigjHge3zgvUbm5eUFPbY7DV6bE1t1AYlJ5nzMxgMrJ09xzj3baP+zQIZz7sgG+54EBjrnDjGzg4BpzrkJgWP/B+CcuyXIPaYBN7bdsxARkTZymHPuQ7+DkMgQqdNKDgb+1mjfG8BjZpYMzAeGmNkA4HvgdOCMJq71F+DxRvtSgIHACqA1q4ZlAR8Ah+F1R8USPbfoE6vPC/TcotWePrdEIBPvtV4EiNyEpTewodG+fCAZ6OGcW29mvwbexPuP/ahz7stgF3LOlQAlQQ593drg7Idqsmudc3mtvU4k0nOLPrH6vEDPLVqF6Ll9E5poJFZEasKyW86514HX/Y5DRERE2l6kTmvOBzIa7csAaoCN4Q9HRERE/BSpCcvHwLGN9k0E5jnnqn2IR0RERHwU9i4hM0sHBjfY1d/MRgGlzrmVgX13AR+Z2VTgOeBovEG1PwlnrLtQAkwn+NiYaFeCnlu0KSE2nxfouUWrEmL3uYlPwj6t2czGA/8Ncuh959z4BuedBNwM7I1XJO6PzrlHwhCiiIiIRBhf67CIiIiINEekjmERERER2U4Ji4iIiEQ8JSwhYGYZZvawmW1osFjjZL/jCiUzOyewWOUbfscSCmZ2rZnNNbPNZlZkZu8ElnyIKrtZJDRqxcrvZ3di8O8q5l8LxT9KWPaQmXUCPsRbLXoyMBQ4H/jWz7hCycyGArcCc/yOJYTGAw/ilQ4/BFgDvG1mg3f1Q5GkwSKhM4GRwN14i4Se7GdcITKeKP/97E6s/V3Fw2uh+EuDbveQmd0MTAGGOeeq/I4n1MysHTAXuANv1ezezrmJ/kYVemaWABQA051zf/U7nubY3SKh/kUWetH4+9mVWPy7ivXXQvGfWlj23MnAR8BfzSzfzJaa2bTAIo2x4E5gsXPuKb8DaWOpQDug2O9AWuBg4K1G+94AcmPo/1+9aPz97Eos/l2dTGy/ForPonYtoQgyKLA9DZwADADuB9KB3/oY1x4zs1OBCcBov2MJg9vxily95HMcLbHLRUKB9WGPqO1E4+8nqBj+u4rZ10KJDGphCSLwqcDtZvt1/el4TdUXOuc+c869AEwDLrYGS5ZGiuY+NzPbC/g7cKZzbqvfcTdHC39vDX/uBrxKyidHy3ONJ7H0+4nGv6sWiKrXQok+amEJ7l7g2d2ckx/4uh5Y6ZyrbXDsK6AD3qfcwtCHt0ea+9zGAz3xlkio358AYGY1wDjn3GdtFGNrteT3BoCZTQcuB46NwOezOzG/SGiU/36C2Z/o+7tqrmh7LZQoo4QlCOfcRpr/gv8RcIiZJTjn6gL79gbKWnCNsGnuczOz/wDDG+2egffC8ytg5U4/5LMW/t4ws9uAC4FjovRNon6R0Fsa7IuZRUJj4PcTTNT9XbVAVL0WSvRRwrLn/gz8FLjbzO7F67e9EbjXRfEUrEBz9ZKG+8ysBEh1zi0J+kNRxMz+CpyHt6Dm92bWO3Co3Dm32b/IWiTSFwlttRj5/ewkxv+uYvK1UCKHpjWHgJlNwPuUuy/wPfAEcEssfMptyMweJwamXwKYWVP/8Z9wzv08nLHsiVhdJDRWfj/NEWN/V3HxWij+UMIiIiIiEU+zhERERCTiKWERERGRiKeERURERCKeEhYRERGJeEpYREREJOIpYREREZGIp4RFREREIp4SFhEREYl4SlhEBDPLCqxm/Vmj/YcH9r/TYF+2mdWZWcfwRyoi8UoJi4gAjAK+AYZYg2WEgdvxFuT7osG+kcCqwLo4IiJhoYRFJETM7OhAa8QFfseyK03EORL4FFgHZAfO+xlQAhSyc8KyMByxiojU02rNIqEzMvA10t/Mg8U5CpgHtAP2NbO1wAzgNGBOo3OVsIhI2ClhEQmdkUAtsMTvQHYjWJwjgYeAVLyVdgfhtbhsBZKBrxqd+4+wRCoiEqCERSR0RgIrnHPlfgeyGzvEaWbpeAnKQqADcA6QCxwS+LrUOVfT4NyBqIVFRMJMY1hEQsDMkoF9gEV+x7IrTcQ5HCh0zm0AFgOnAM8459bgJTdfNDh3BLDZOZcXloBFRALUwiISGsOAFCK/5SFYnKMafL8KOBaYG/h+JPBug3NHAl+aWWqDfXXOuao2iVZEJEAtLCKhsdNAVjMrMrP+LbmImSWb2Xoz6xHkmJlZajO3pv62gw24HUmgxcV53mkwZblxC8tIvK6i8gbbmy15jiIiraEWFpHQ2CERMLO+eDNu1rTkIs65aiCzicOHAB8081JHAu8F2b9TwuKc+9Uu4slu9P2vgCbPFxFpK0pYREJjJFDknFsb+D4HWOaccyG8x9fAec08d1kT+xvHCXitQcBo59y3exCfiEibUcIiEhoj2Hkg6xozexiYDGwCJjvndjnGxcwuBQ5wzp3b+JhzrgB4PMRxtro1SEQknDSGRWQPmVkGkMGO40JygAPx6pV0B94HLmnG5XZKKEKliTihbVqDRERCSgmLyJ4LNpA1B5jqnJvjnKvD66JpTkLQZgkLTVfi3d4aZGbFZrbSzEYiIhJBlLCI7Ln6N/dFAIEZOvsCrzY4Zz9g6a4uElh0MIe2T1gaX7+1rUEiImGjhEVkDznnbnfOmXPus8CuQUCZc259g9O2Tx3ehQFAeaCAWzjirNfa1iARkbBRwiISesNpkJw0qC7bsEbL42b2eKOfa8vuoKBa2xokIhJuSlhEQq9xt84wvNL3xQ32ZQEfNfq5HRKdMGlta5CISFiZJgaIhJeZJeElBCMDheL8jOVU4FLn3NGB75PxVmjObJRgiYj4SnVYRMIssPLxvn7HEdCc1iAREd+phUVEREQinsawiIiISMRTwiIiIiIRTwmLiIiIRDwlLCIiIhLxNEtIZA+ZmUauN+CcM79jEJHYo1lCIiIiEvHUJSQiIiIRTwmLSAQzs1+Y2XwzKzMz12jbGljhWUQk5ilhEYlQZnY78DCwBfg/4GagInD4TeAhpz5dEYkTGsMiEoHM7GC8xREfBS6oT0zM7HLgbmCic+5NH0MUEQkrtbCIRKYrgM3AlY1aUd4PfB0W9ohERHykhEUkwphZAnAc8IZzbkujwymBr1sC52YFxrN0CWOIIiJhp4RFJPL0AzoDS4IcGxv4+lXg6yjgW+dcSduHJSLiHyUsIpGnY+BrVZBjpwPfA/MD348EvghDTCIivlLCIhJ5vgMccEzDnWZ2AXA4cKtzrjawexSwMKzRiYj4QAmLSIQJdO/MBI41s+fM7BIzewJ4CHgGuK/B6Tu0sJjZPWa2366ub2a/MbP+IQ9cRKQNaVqzSAQys87An4FTgDRgOfB34O8Npjin480kGuKcWxXY9zZwgnMuWHeSiEjUUsIiEqXM7CDgDaBLgyRmAV79lmOAW5xz/wjyc/91zh0Z1mBFRPaQuoREotcoYFGDZKU30AG4DpgETGj8A2aWAeSHMUYRkZBQwiISvRrPEBoO/NM5VwZ0B1YF+ZnhwKK2D01EJLSS/A5ARFrHOferRrty+CGBGUkgMTGzLOfc2sD+4cDisAQoIhJCamERiR0NW09GAovMLAl4solzRESihgbdisQwMzsAGOGce9jvWERE9oQSFhEREYl46hISERGRiKeERURERCKeEhYRERGJeEpYREREJOIpYREREZGIp4RFREREIp4SFhEREYl4SlhEREQk4v1/1UxG6UWZiQUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Expected residuals for normally distributed data\n", "x = np.linspace(-6., 6., 50)\n", "bin_width = x[1] - x[0]\n", "normal = bin_width*ds.shape[0]*(1/np.sqrt(2*np.pi))*np.exp(-0.5*x**2)\n", "\n", "# Histogram Residuals\n", "fig = plt.figure(figsize=(8, 5))\n", "sns.histplot((ds.I - ds.IML)/ds.SIGI, bins=x, color=\"k\", label=\"Residuals\")\n", "plt.plot(x, normal, 'r-', lw=2, label=\"Expected\")\n", "plt.ylabel(\"Log Count\")\n", "plt.yscale(\"log\")\n", "plt.ylim(1e0, 3e5)\n", "plt.xlabel(r\"$\\frac{I_{h,i} - I^{ML}_{h}}{\\sigma_{I_{h,i}}}$\", fontsize=24)\n", "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These residuals are not symmetric about $0$, which would have been expected for truly normally distributed data. In addition, the tails are much \"heavier\" than for a normal distribution. This suggests that it may be possible to do a better job merging these scaled observations by using a different error model that is more robust to outliers." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Generalized Merging with Flexible Error Model\n", "\n", "The inverse-variance weighting scheme implemented in `rs.algorithms.merge()` is the maximum likelihood estimator for the true mean if we assume the observations are normally distributed. \n", "However, we can write a more general form for the maximum likelihood estimator for the mean of each intensity distribution, $\\mu \\in\\mathbb{R}^{|\\mathbf{H}|}$, without assuming a specific distribution for the error model. \n", "Therefore, we will maximize the probability of the data given the model\n", "\\begin{align*}\n", "P(data | model) &= \\prod_{h,i}P(I_{h,i}|\\mu_h, \\sigma_{I_{h,i}}) \n", "\\end{align*}\n", "by minimizing the negative log likelihood\n", "\\begin{align*}\n", "\\mathcal{L} &\\triangleq -\\log P(data | model) \\\\\n", "&= -\\sum_{h,i}\\log P(I_{h,i}|\\mu_h, \\sigma_{I_{h,i}}) \\\\\n", "I^{ML} &= \\underset{\\mu}{\\mathrm{argmin}}\\ -\\sum_{h,i} \\log P(I_{h,i} | \\mu_h, \\sigma_{I_{h,i}}) \n", "\\end{align*}\n", "with respect to $\\mu$. With this formulation, it is possible to supply any parametric form for the error model belonging to the [location-scale family](https://en.wikipedia.org/wiki/Location%E2%80%93scale_family) of distributions. This maximum likelihood estimator is implemented in the function below." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def merge_mle(ds, distribution, *args, lr=.5, progress_bar=True, return_loss=False, **kwargs):\n", " \"\"\"\n", " Merge observations using the provided distribution as an error model.\n", " Additional arguments or keyword arguments will be passed to the PyTorch \n", " distribution constructor.\n", " \n", " Parameters\n", " ----------\n", " ds : rs.DataSet\n", " Scaled, unmerged observations\n", " distribution : torch.distributions.Distribution\n", " PyTorch distribution to use as error model\n", " lr : float\n", " Learning rate for Adam optimizer\n", " progress_bar : bool\n", " Whether to display a progress bar for optimization\n", " return_loss : bool\n", " Whether to return the loss function values during optimization\n", " \n", " Returns\n", " -------\n", " rs.DataSet or (rs.DataSet, losses)\n", " Merged DataSet with or without list of loss function values\n", " \"\"\"\n", " # Compute MLE with normal error model\n", " ds = ds.copy()\n", " mle_norm = rs.algorithms.merge(ds).stack_anomalous()\n", " ds[\"IML\"] = mle_norm.loc[ds.index, \"I\"]\n", " \n", " # Observed intensities and error estimates\n", " groupby = ds.groupby(['H', 'K', 'L'])\n", " idx = groupby.ngroup().to_numpy()\n", " I = torch.as_tensor(ds.I.to_numpy())\n", " SigI = torch.as_tensor(ds.SIGI.to_numpy())\n", " \n", " # Initialize optimization at MLE with normal error model\n", " mle = groupby.first()[\"IML\"].to_numpy()\n", " mean = torch.tensor(mle, requires_grad=True)\n", " \n", " # Define loss function\n", " def _evaluate_loss():\n", " return -torch.sum(distribution(*args, loc=mean[idx], scale=SigI, **kwargs).log_prob(I))\n", "\n", " # Setup and fit model\n", " losses = []\n", " opt = torch.optim.Adam([mean], lr=lr) \n", " for _ in tqdm(range(300), disable=not progress_bar):\n", " opt.zero_grad()\n", " loss = _evaluate_loss()\n", " losses.append(loss.detach())\n", " loss.backward()\n", " opt.step()\n", " grad = torch.autograd.grad(_evaluate_loss(), mean, create_graph=True)[0]\n", " hess = torch.autograd.grad(grad.sum(), mean, create_graph=True)[0]\n", " \n", " # Package results\n", " results = rs.DataSet({'I': mean.detach().numpy(), \n", " 'SIGI': np.sqrt(1./hess.detach().numpy())},\n", " index=groupby.first().index, \n", " spacegroup=ds.spacegroup,\n", " cell=ds.cell, \n", " merged=True)\n", " results.infer_mtz_dtypes(inplace=True)\n", "\n", " if return_loss:\n", " return results.unstack_anomalous(), losses\n", " return results.unstack_anomalous()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function must be passed a `torch.distribution.Distribution` to use as the error model. It then uses the [Adam optimizer](https://pytorch.org/docs/stable/optim.html#torch.optim.Adam) to minimize the negative log-likelihood and fit the merged intensities, $I^{ML}$. For stability, the model is initialized using the mean intensity values from `rs.algorithms.merge()`. The following cell fits the model using a [Student's _t_-distributed error model](https://pytorch.org/docs/stable/distributions.html#studentt) with `df=4.0`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "6cd5820834fa44d5b66d870ad7bd24e1", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/300 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "result1, loss1 = merge_mle(ds, torch.distributions.StudentT, 4.0, \n", " progress_bar=True, return_loss=True)\n", "\n", "# Plot loss function\n", "plt.plot(loss1)\n", "plt.xlabel(\"Iteration\")\n", "plt.ylabel(\"Loss\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
I(+)SIGI(+)I(-)SIGI(-)
HKL
1814913.3830950.5903978312.2444470.5472484
25241481.6251516.697405487.1359616.214811
3913762.359973.602170761.640683.1317017
29295235.866648.579554235.866648.579554
30145174.98365.340581189.37324.8809543
\n", "
" ], "text/plain": [ " I(+) SIGI(+) I(-) SIGI(-)\n", "H K L \n", "18 14 9 13.383095 0.59039783 12.244447 0.5472484\n", "25 24 1 481.62515 16.697405 487.13596 16.214811\n", "39 13 7 62.35997 3.6021707 61.64068 3.1317017\n", "29 29 5 235.86664 8.579554 235.86664 8.579554\n", "30 14 5 174.9836 5.340581 189.3732 4.8809543" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result1.sample(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a sanity check we can test this implementation by passing the normal distribution as the error model. In this case, the initial values should match the maximum likelihood estimate and optimization should not change the estimates. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "df1d4a7dda514a0d84aef6c962136263", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/300 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "result2, loss2 = merge_mle(ds, torch.distributions.Normal, lr=1e-3,\n", " progress_bar=True, return_loss=True)\n", "result2 = result2[[\"I(+)\", \"SIGI(+)\", \"I(-)\", \"SIGI(-)\"]].stack_anomalous()\n", "\n", "# Plot loss function\n", "fig, ax = plt.subplots(ncols=2, figsize=(12, 5.5))\n", "ax[0].plot(loss2)\n", "ax[0].set_ylim(4.15e6, max(loss1))\n", "ax[0].set_xlabel(\"Iteration\")\n", "ax[0].set_ylabel(\"Loss\")\n", "ax[0].set_title(\"Loss Function\")\n", "ax[1].loglog(merged_normal[\"I\"].to_numpy(), result2.loc[merged_normal.index, \"I\"].to_numpy(), \"k.\", alpha=0.5)\n", "ax[1].set_xlabel(r\"$I^{ML}$ (AIMLESS)\")\n", "ax[1].set_ylabel(r\"$I^{ML}$ (PyTorch)\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see here that the loss function does not change during the optimization, and that the final maximum likelihood estimates for the intensities are equivalent to the input. This validates that the merging function is working as expected." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Assess Student's _t_-Distributed Error Model\n", "\n", "A Student's _t_-distribution is useful in modeling data that contain outliers. This distribution places more density in the tails than a normal which makes maximum likelihood estimates robust to outlying measurements. It is parameterized by a degree of freedom, $\\nu$, which can take a value between $[0, \\infty)$ and controls the heaviness of the tails. As $\\nu\\to\\infty$, the probability density function approaches a normal distribution.\n", "\n", "Below, we will set up a few helper functions for merging our data within randomly partitioned half-datasets. We will do this for the Student's t-distribution scanning several degrees of freedom and for the normal distribution so that we can compare the results. Based on the previous example, we will compute $CC_{anom}$ using a Spearman correlation coefficient and repeated 2-fold cross-validation to compare the different error models." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def sample_halfdatasets(data):\n", " \"\"\"Randomly split DataSet into two equal halves by BATCH\"\"\"\n", " batch = data.BATCH.unique().to_numpy(dtype=int)\n", " np.random.shuffle(batch)\n", " halfbatch1, halfbatch2 = np.array_split(batch, 2)\n", " half1 = data.loc[data.BATCH.isin(halfbatch1)]\n", " half2 = data.loc[data.BATCH.isin(halfbatch2)]\n", " return half1, half2" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def merge_dataset(dataset, nsamples, distribution, *args, **kwargs):\n", " \"\"\"\n", " Merge dataset with repeated 2-fold cross-validation using `distribution`\n", " as error model.\n", " \"\"\"\n", " dataset = dataset.copy()\n", " samples = []\n", " for n in tqdm(range(nsamples)):\n", " half1, half2 = sample_halfdatasets(dataset)\n", " mergedhalf1 = merge_mle(half1, distribution, *args, progress_bar=False, **kwargs)\n", " mergedhalf2 = merge_mle(half2, distribution, *args, progress_bar=False, **kwargs)\n", " result = mergedhalf1.merge(mergedhalf2, on=[\"H\", \"K\", \"L\"], suffixes=(1, 2))\n", " result[\"sample\"] = n\n", " samples.append(result)\n", " return rs.concat(samples).sort_index()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def merge_dataset_normal(dataset, nsamples):\n", " \"\"\"\n", " Merge dataset with repeated 2-fold cross-validation using normal distribution\n", " as error model.\n", " \"\"\"\n", " dataset = dataset.copy()\n", " samples = []\n", " for n in tqdm(range(nsamples)):\n", " half1, half2 = sample_halfdatasets(dataset)\n", " mergedhalf1 = rs.algorithms.merge(half1)\n", " mergedhalf2 = rs.algorithms.merge(half2)\n", " result = mergedhalf1.merge(mergedhalf2, on=[\"H\", \"K\", \"L\"], suffixes=(1, 2))\n", " result[\"sample\"] = n\n", " samples.append(result)\n", " return rs.concat(samples).sort_index()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Note:* Using these settings, the following will take ~15 min to complete. For the pre-print, this was run with `nsamples=15`. It is possible to run into occasional numerical instabilities when computing the Hessian for `df=4.0`, but these will only impact $\\sigma_I^{ML}$ estimates which are not used here. " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "e90b8157e9544b92993a1bf9fdbd16de", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def plot(results, label, ax, color=None):\n", " ax.errorbar(results.index, results[\"mean\"], yerr=results[\"std\"], color=color, label=label)\n", " return\n", "\n", "fig = plt.figure(figsize=(10, 6))\n", "ax = fig.gca()\n", "with sns.color_palette(\"viridis\", 5) as palette:\n", " plot(results4, r\"Student-T ($d.f.=4$)\", ax, color=palette[4])\n", " plot(results8, r\"Student-T ($d.f.=8$)\", ax, color=palette[3])\n", " plot(results16, r\"Student-T ($d.f.=16$)\", ax, color=palette[2])\n", " plot(results32, r\"Student-T ($d.f.=32$)\", ax, color=palette[1])\n", " plot(results64, r\"Student-T ($d.f.=64$)\", ax, color=palette[0])\n", "plot(resultsinf,\"Normal\", ax, color=\"k\")\n", "plt.ylabel(r\"$CC_{anom}$ (Spearman)\")\n", "plt.xlabel(r\"Resolution Bin ($\\AA$)\")\n", "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), title=\"Error Model\")\n", "plt.ylim(0, 0.72)\n", "\n", "plt.xticks(resultsinf.index, labels, rotation=45, ha='right', rotation_mode='anchor')\n", "plt.grid(axis=\"y\", linestyle='--')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Summary \n", "\n", "This dataset was scaled and merged in AIMLESS, which involves rounds of outlier rejection and implicitly assumes that the intensities are normally distributed about the true mean. In the [first section](3_mergingerrormodel.ipynb#Normal-Error-Model), we observed that the residuals from merging are not normally distributed, suggesting that the normal error model may have been suboptimal. We implemented a more general maximum likelihood-based approach for merging data using different probability distributions and used it to evaluate the performance of Student's _t_-distributed error models. This seemed like a reasonable starting point because it is often more robust to outliers than a normal distribution. \n", "\n", "Using repeated 2-fold cross-validation, we saw that a Student's _t_-distribution with a low degree of freedom ($\\nu=4$) outperforms the normally distributed error model when assessed using $CC_{anom}$. Furthermore, the performance seems to approach that of the normal distribution when the degree of freedom is increased, which is expected since the _t_-distribution approaches a normal disitrbution as $\\nu\\to\\infty$.\n", "\n", "Although this dataset is quite high quality and was used to phase and refine a model from the [native sulfur SAD signal](http://doi.org/10.2210/pdb7L84/pdb), it still shows that there can be an incremental improvement in $CC_{anom}$ from revisiting assumptions regarding error models during merging. This model is implemented in ~40 lines of code using `PyTorch`, and can be quickly applied to any dataset of interest. By lowering the barrier to implementing such models, `reciprocalspaceship` makes it easy to try new analyses and to revisit some of the assumptions made in crystallographic data reduction." ] } ], "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.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }