diff --git a/examples/example_measurement.ipynb b/examples/example_measurement.ipynb new file mode 100644 index 00000000..178c60ff --- /dev/null +++ b/examples/example_measurement.ipynb @@ -0,0 +1,291 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0a432940", + "metadata": {}, + "source": [ + "-This notebook is an example of how to properly set up the ScreenBeamProfileMeasurement Class:\n", + " -Since ScreenBeamProfileMeasurement takes a screen object as a pydantic field to use for live images and this example is mean't to be \n", + "done offline we subclass Screen and overwrite the the image property to a static 2d array. \n", + " -The array passed in this example is automatically" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "80146949", + "metadata": {}, + "outputs": [], + "source": [ + "from lcls_tools.common.measurements.screen_beam_profile_measurement import ScreenBeamProfileMeasurement\n", + "from lcls_tools.common.data_analysis.fit.projection import ProjectionFit\n", + "from lcls_tools.common.data_analysis.fit.methods import GaussianModel\n", + "from lcls_tools.common.image_processing.image_processing import ImageProcessor\n", + "from lcls_tools.common.image_processing.roi import ROI, RectangularROI, CircularROI\n", + "from lcls_tools.common.devices.device import Device, PVSet, ControlInformation, Metadata\n", + "from lcls_tools.common.devices.screen import ScreenControlInformation, ScreenPVSet, Screen\n", + "from matplotlib import pyplot as plt\n", + "from epics import PV, caget\n", + "import numpy as np\n", + "import h5py\n", + "import random" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7051ec06-01f3-4dc3-b86a-c732fd1c3df5", + "metadata": {}, + "outputs": [], + "source": [ + "class ScreenTest(Screen):\n", + " @property\n", + " def image(self) -> np.ndarray:\n", + " return self._image\n", + " \n", + " @image.setter\n", + " def image(self,image):\n", + " self._image = image" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f736c366-6c88-4187-b9a4-a4dcab9c7730", + "metadata": {}, + "outputs": [], + "source": [ + "def create_test_image(size:tuple,center:list,radius:int):\n", + " # make img that is a circle in the center of the image with known standard dev and mean. no imports, no calls to external or\n", + " # internal files. \n", + " image = np.zeros(size)\n", + " for y in range(image.shape[0]):\n", + " for x in range(image.shape[1]):\n", + " distance = np.sqrt((x - center[0])**2 + (y - center[1])**2)\n", + " if distance < radius:\n", + " image[y,x]= 1\n", + " return image" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "9fa172e9", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# creating projection fit class\n", + "x = np.linspace(0,800,800)\n", + "\n", + "\n", + "gauss_model = GaussianModel()\n", + "projection_fit = ProjectionFit(model = gauss_model,use_priors =True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "4289cca5", + "metadata": {}, + "outputs": [], + "source": [ + "# creating image processing class\n", + "# would be nice to have default way of finding the center of image\n", + "# also maybe have rectangularROI have default roi_type 'Rectangular', and change roi_type to roi\n", + "roi = RectangularROI(center =[400,400],xwidth=300,ywidth=300)\n", + "image_processor = ImageProcessor(roi = roi)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "8260d6b6", + "metadata": {}, + "outputs": [], + "source": [ + "PV_strings = {\n", + " 'image': 'ArrayData',\n", + " 'n_bits': 'N_OF_BITS',\n", + " 'n_col': 'ArraySize1_RBV',\n", + " 'n_row': 'ArraySize0_RBV',\n", + " 'resolution': 'RESOLUTION'\n", + "}\n", + "metadata = {\n", + " 'area': 'TEST',\n", + " 'beam_path': ['SC_TEST'],\n", + " 'sum_l_meters': 99.99\n", + " }\n", + "control_name = \"OTRS:TEST:650:\"" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "eb543534", + "metadata": {}, + "outputs": [], + "source": [ + "# creating Screen Device\n", + "# inconsistent use of Controls_information and control_information between classes recommend change.\n", + "screen_pvs = ScreenPVSet(**PV_strings)\n", + "meta_data = Metadata(**metadata)\n", + "controls_information = ScreenControlInformation(control_name = control_name, PVs = screen_pvs )\n", + "screen_test = ScreenTest(controls_information = controls_information, metadata = meta_data)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "d5d1add9", + "metadata": {}, + "outputs": [], + "source": [ + "test = ScreenBeamProfileMeasurement(name = control_name,device = screen_test, image_processor = image_processor, fitting_tool = projection_fit)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "71c61cf3-bf4a-4d56-beb7-07e55b61cce0", + "metadata": {}, + "outputs": [], + "source": [ + "test.device.image = create_test_image((800,800),[400,400],50)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "725f6fb5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['mean', 'sigma', 'amplitude', 'offset']\n" + ] + }, + { + "data": { + "text/plain": [ + "[{'raw_image': array([[0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " ...,\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.]]),\n", + " 'processed_image': array([[0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " ...,\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.],\n", + " [0., 0., 0., ..., 0., 0., 0.]]),\n", + " 'mean_x': 150.50167233387222,\n", + " 'sigma_x': 31.481275393806968,\n", + " 'amplitude_x': 98.871501132453,\n", + " 'offset_x': 0.9887150113245301,\n", + " 'mean_y': 150.50167233081507,\n", + " 'sigma_y': 31.481275399702522,\n", + " 'amplitude_y': 98.87150113245299,\n", + " 'offset_y': 0.9887150113245299}]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print(test.fitting_tool.model.param_names)\n", + "test.measure()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "29d792d7-dac2-45ad-888e-3cb29921ce92", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(
, )" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwiklEQVR4nO3deXiU1fn/8fdMVsgGSUhIQghh3xeDyKq4QVGxLlVaW0SLtmhFgWor1Wpdvlq1UlQKFTfwVxfqgrWKS4qyI/smoOyEJSEkkH3PPL8/nsxAyDZJZsnyeV3XXDN55jwzdyYDc8859znHYhiGgYiIiIiXWL0dgIiIiLRuSkZERETEq5SMiIiIiFcpGRERERGvUjIiIiIiXqVkRERERLxKyYiIiIh4lZIRERER8SpfbwfgDJvNxsmTJwkJCcFisXg7HBEREXGCYRjk5uYSGxuL1Vpz/0ezSEZOnjxJfHy8t8MQERGRBjh27BidOnWq8f5mkYyEhIQA5i8TGhrq5WhERETEGTk5OcTHxzs+x2vSLJIR+9BMaGiokhEREZFmpq4SCxWwioiIiFcpGRERERGvUjIiIiIiXqVkRERERLxKyYiIiIh4lZIRERER8SolIyIiIuJVSkZERETEq5SMiIiIiFcpGRERERGvUjIiIiIiXqVkRERERLyqWWyUJyItkK0cjq6D1O2QdwratIcOvaHr5eDf1tvRiYgH1btnZNWqVUycOJHY2FgsFguffPJJneesXLmSpKQkAgMD6dq1K//85z8bEquItARlJbB+PrzYGxZfB18/CutegeVPwvu3wQvd4atHoPCstyMVEQ+pdzKSn5/PoEGDmDdvnlPtDx8+zDXXXMOYMWPYtm0bf/rTn7j//vv56KOP6h2siDRzp/fBq5fCV7MhP93sDel3Iwy/Fwb9AsI6Q2k+rJ8HryTBwW+8HbGIeIDFMAyjwSdbLCxdupQbbrihxjZ//OMf+fTTT9m7d6/j2LRp09ixYwfr16936nlycnIICwsjOzub0NDQhoYrIt508BtYMhlK8iCoA0WX/omMbjeTW2qhoKSMAF8fQgN8CU9bSdDKJ7Cc/gGwwE/+CsOneTt6EWkAZz+/3V7Aun79esaNG1fp2Pjx49m8eTOlpaXVnlNcXExOTk6li7jW/PnzSUxMJDAwkKSkJFavXl1r+7/85S9YLJZKl44dOzb6caV1KN/3P2zvTIKSPH4IGMi1pX+l99JoRv9tDRNeWs3NC9Zz3StruPRvK+j/L4PhGX/mm7Y/AQz48o+UrH7Z27+CiLiR25ORtLQ0oqOjKx2Ljo6mrKyMjIyMas959tlnCQsLc1zi4+PdHWarsmTJEmbMmMEjjzzCtm3bGDNmDBMmTCAlJaXW8/r160dqaqrjsmvXLpc8rrRcx88W8MYHn1D67i+w2kr4ujyJidkPsjunDQABvlYig/3pHN6WqJAA2vr7AHCq0MKvz0zmpbIbAfBf/mfeee1v7Dqe7bXfRUTcyGgEwFi6dGmtbXr06GE888wzlY6tWbPGAIzU1NRqzykqKjKys7Mdl2PHjhmAkZ2d3ZhwnXbZZZcZ9913n/HAAw8Y7dq1M6KiooxXX33VyMvLM+644w4jODjY6Nq1q7Fs2TLHOTabzXjuueeMxMREIzAw0Bg4cKDxwQcfVHrcL774whg1apQRFhZmhIeHG9dee61x4MCBSs87ffp046GHHjLat29vREdHG48//rjLf79hw4YZ06ZNq3Ssd+/exsMPP1zjOY8//rgxaNAglz+utEyncgqNP3ywwxg2+x3j+GOJhvF4qLH28THG79/baCxed9jYeDjTyC4sqfbcwpIyY9fxLGPJphTjkY93GO88OdkwHg81Ch+LMCY+/JJx+xsbjB9Sczz8G4lIQ2RnZzv1+e32npGOHTuSlpZW6Vh6ejq+vr5ERERUe05AQAChoaGVLp62ePFiIiMj2bhxI9OnT+eee+7hlltuYeTIkWzdupXx48czefJkCgoKAHj00Ud56623WLBgAbt372bmzJn86le/YuXKlY7HzM/PZ9asWWzatInly5djtVq58cYbsdlslZ43KCiIDRs28Pzzz/Pkk0+SnJxcbYzPPPMMwcHBtV4uHCYpKSlhy5YtVYbOxo0bx7p162p9Tfbv309sbCyJiYn8/Oc/59ChQy55XGk5DMPg/313lCv+tpIlm1N43mcBcZZM8kMSufih//C3n1/M7SO6cHGXcEID/ap9jEA/H/rHhXHr0HievnEgv/jTm2TFX0GgpZR5/q+wZd9Rrnl5Nc8s20tRabmHf0MRcQePFLD+97//Zc+ePY5j99xzD9u3b2+yBaxjx46lvLzc8UFeXl5OWFgYN910E2+//TZgDj/FxMSwfv16BgwYQGRkJN988w0jRoxwPM5dd91FQUEB7777brXPc/r0aaKioti1axf9+/ev8rwAw4YN44orruCvf/1rlfPPnDnDmTNnav1d4uLiaNOmjePnkydPEhcXx9q1axk5cqTj+DPPPMPixYv58ccfq32cL774goKCAnr27MmpU6d4+umn+eGHH9i9ezcRERENflxpOc7mlzDz39tZ8eNpAGZHrua3eQvANxB+sxKiejf8wYuyYcFoyE5hbchP+OXp2wHoGR3M/F8m0T0q2BW/goi4mLOf3/Ve9CwvL48DBw44fj58+DDbt28nPDyczp07M3v2bE6cOOH40J42bRrz5s1j1qxZ3H333axfv5433niD9957rwG/lucMHDjQcdvHx4eIiAgGDBjgOGavg0lPT2fPnj0UFRVx9dVXV3qMkpIShgwZ4vj54MGD/PnPf+a7774jIyPD0SOSkpJC//79qzwvQExMDOnp6dXGGB4eTnh4eIN+P4vFUulnwzCqHDvfhAkTHLcHDBjAiBEj6NatG4sXL2bWrFkNflxpGX5Iy+GuxZs5fraQAF8rT10Wwi0bFpl3Xv1k4xIRgMAwuOlVeOsaRuV+ySfjbuaudeHsO5XHDf9Yy9xJg7mqb3TdjyMiTVK9h2k2b97MkCFDHB+ys2bNYsiQITz22GMApKamVipYTExMZNmyZaxYsYLBgwfz1FNP8fLLL3PzzTe76FdwDz+/yl3IFoul0jH7B6zNZnMkFZ9//jnbt293XPbs2cOHH37oOGfixIlkZmby2muvsWHDBjZs2ACYSUttz3v+MM75GjJMExkZiY+PT7VDZxcWGtcmKCiIAQMGsH//fpc+rjQ/Gw5lcss/13P8bCEJEW355N6R3Hr6FSxlhZAwGi6+2zVPlDASRvwOgMG7/o8v7h3KsMRw8orL+M3/28y7G1QoLdJc1btnZOzYsdQ2srNo0aIqxy677DK2bt1a36dqNvr27UtAQAApKSlcdtll1bbJzMxk7969vPrqq4wZMwaANWvWNOp5p02bxq233lprm7i4uEo/+/v7k5SURHJyMjfeeKPjeHJyMj/96U+dfu7i4mL27t3r+F1c9bjSvGw4lMntb26kuMzGxV3a8/rtFxOWkgz7vgSrH1w3B6wuLE0b+zB8/zGcPUKHHfN5566HeXTp9yzZfIw/Ld1FcVk5d45KdN3ziYhHaG8aFwgJCeHBBx9k5syZ2Gw2Ro8eTU5ODuvWrSM4OJgpU6bQvn17IiIiWLhwITExMaSkpPDwww836nkbOkwza9YsJk+ezNChQxkxYgQLFy4kJSWFadPOLSw1b948li5dyvLlywF48MEHmThxIp07dyY9PZ2nn36anJwcpkyZUq/HlZZj5/Espi7eTHGZjct7dWDBr5IItNog+c9mg5H3QYdern3SgBD4yTPwwR2w9mX8ku7krzcPIDzYnwUrDvLEf/fg62Nl8vAE1z6viLiVkhEXeeqpp4iKiuLZZ5/l0KFDtGvXjosuuog//elPAFitVt5//33uv/9++vfvT69evXj55ZcZO3asx2OdNGkSmZmZPPnkk6SmptK/f3+WLVtGQsK5/8AzMjI4ePCg4+fjx4/zi1/8goyMDDp06MDw4cP57rvvKp3jzONKy7D/VC5T3txIXnEZw7uGm4mInw9sWgSZB6BtBIyeVefjNEjfGyD+Eji2AVY8i+X6l/nDeDPpWbDiII//53vi2gVyRW8ND4o0F42aTeMpWg5epOlIzy3i+lfWkpZTxKBOYbxz93CCA3yhpABeHmzuwDvhBbjkN+4LImUDvDkOLFa4dwN06IlhGDz80S6WbD5GkL8PH0wbSd9Y/X8h4k1NZjl4EWk5ysptTH93G2k5RXTrEMSiO4eZiQjA1sVmItKuMyTd4d5AOl8Cva4FwwarXwTMYu+nb+zPyG4R5JeUM3XxJk7lFLk3DhFxCSUjIuK057/6kQ2HzxDk78Ork4fSPsjfvKO0CNbMNW+P+T34+rs/mMseMq93/RsyzSFFPx8rC36ZRLcOQaRmF3HPv7ZQVl79bDQRaTqUjIiIU5btSmXhKnPV3b/dMqjyQmPb/h/kpUFoJxh0m2cCih0CPcZX9I7McRwOa+vHW3cMIyTQl60pWbz8zYFaHkREmgIlIyJSp+NnC/jDhzsB+M2lXZkwIObcnTYbbPineXvU/Z7pFbG7tKJ3ZOcSyD23xk3niLb8343mIoXzvtnPxsO1r1QsIt6lZEREamUYBrM/3kVecRlJCe0dM1ccDi43Z9AEhMJgD/WK2MVfbM6ssZXCxoWV7rp+UCw3X9QJmwEzl2wnu7DUs7GJiNOUjIhIrZZsOsbq/RkE+Fp54WcD8fW54L8Ne6/IkMnmOiCeNuI+83rTG1CSX+muJ37aj4SItpzIKuQvn+72fGwi4hQlIyJSo5NZhfzf53sBeHBcL7p2uGBDutP74MD/AAsMc9Gy7/XV+1ponwhFWfD9R5XuCg7wZe6kwVgtsHTbCdYeyPBOjCJSKyUj9TR27FhmzJjh7TBE3M4wDP60dBe5xWUM6dyOX4+uZpl1+9BIrwkQ7qVl2K0+MPRO8/bmt6rcPaRze24f0QWAP3/yPcVl5R4MTkScoWSknj7++GOeeuopb4fhMvPnzycxMZHAwECSkpKqbKxXnVWrVjFx4kRiY2OxWCx88sknLntsaTq+3nOKFT+ext/Xygs/G4SP9YLdlwuzYPu75u1Lfuvx+CoZdJu5F87JrZC6o8rds8b1pENIAIcy8lm48pAXAhSR2igZqafw8HBCQrwwLu4GS5YsYcaMGTzyyCNs27aNMWPGMGHChEq7LlcnPz+fQYMGMW/ePJc/tjQNxWXlPLPMHJ75zZiulafx2m1/F0rzoUMfSKx+g0iPCe4AfSaat7csqnJ3aKAff76uLwDzvj3A0cz8Km1ExHuUjFTjww8/ZMCAAbRp04aIiAiuuuoq8vPN/7wuHKbJzc3ll7/8JUFBQcTExPD3v/+9UpuxY8cyffp0ZsyYQfv27YmOjmbhwoXk5+dz5513EhISQrdu3fjiiy8qxfDll18yevRo2rVrR0REBNddd12lvWJcYc6cOUydOpW77rqLPn36MHfuXOLj41mwYEGt502YMIGnn36am266yeWPLU3DorVHOJpZQFRIAPeM7Va1gWHAtn+Zty+eChZL1TaeZl/1decHUJxX5e6JA2MY3T2S4jIbT/53j2djE5FaKRm5QGpqKr/4xS/49a9/zd69e1mxYgU33XQTNW3hM2vWLNauXcunn35KcnIyq1evZuvWrZXaLF68mMjISDZu3Mj06dO55557uOWWWxg5ciRbt25l/PjxTJ48mYKCAsc5+fn5zJo1i02bNrF8+XKsVis33ngjNlvV1SSfeeYZgoODa71cOERSUlLCli1bGDduXKXj48aNY926dQ19+dz+2OJ+p3OLeaViobA//KQ3QQHV7KeZthPSd4NPAAz4mYcjrEHipRDeFUpyqxSygrlc/JM/7Yev1cLyH9K19ohIE6Jdey+QmppKWVkZN910k2O32QEDBlTbNjc3l8WLF/Puu+9y5ZVXAvDWW28RGxtbqd2gQYN49NFHAZg9ezZ//etfiYyM5O67zdkHjz32GAsWLGDnzp0MHz4cgJtvvrnSY7zxxhtERUWxZ88e+vfvX+m+adOmceutt9b6e8XFxVX6OSMjg/LycqKjK+9sGh0dTVpaGo3hzscW95uT/CN5xWUM7BTGTUPiqm9krxXpfS20ae+54GpjsZi9I8mPmUM1SVOqNOnaIZhJF8fzzoYU/vrFXj66ZySWptCrI9LKqWfkAoMGDeLKK69kwIAB3HLLLbz22mucPXu22raHDh2itLSUYcOGOY6FhYXRq1flRaEGDhzouO3j40NERESlBMf+oZ2enu44dvDgQW677Ta6du1KaGgoiYnmTIXqai7Cw8Pp3r17rZc2bdpU+ztc+B+xYRgu+8/ZnY8t7nEgPZf3Nx0D4LHr+mK9sGgVoKwEdv7bvD34lx6MzgmDf3mukDVtV7VNHriyB238fNiakkXynlMeDlBEqqNk5AI+Pj4kJyfzxRdf0LdvX1555RV69erF4cOHq7S1D91U96F7Pj8/v0o/WyyWSsfs558/BDNx4kQyMzN57bXX2LBhAxs2bADMIZALNWSYJjIyEh8fnyo9Fenp6VV6NOrLnY8t7vXy8gMYBozrG83QLuHVN9r/FRSegeCO0O1yzwZYl6BIc5oxnEuYLhAVGsivR3cB4IWvfqTcVv0QrIh4jpKRalgsFkaNGsUTTzzBtm3b8Pf3Z+nSpVXadevWDT8/PzZu3Og4lpOTw/79+xv1/JmZmezdu5dHH32UK6+8kj59+tTYOwPmMM327dtrvQwdOrTSOf7+/iQlJZGcnFzpeHJyMiNHjmxU/O58bHGf/ady+e/OkwA8cFWPmhtuf8+8HjTJXOOjqRlYMWT5/UfmvjnV+O1l3WjX1o/96Xl8tPW4B4MTkeqoZuQCGzZsYPny5YwbN46oqCg2bNjA6dOn6dOnT5W2ISEhTJkyhYceeojw8HCioqJ4/PHHsVqtjRqOaN++PRERESxcuJCYmBhSUlJ4+OGHa2wfHh5OeHgN32JrMWvWLCZPnszQoUMZMWIECxcuJCUlhWnTpjnazJs3j6VLl7J8+XLHsby8PA4cOLcT6uHDh9m+fTvh4eF07tzZ6ceWpuXlb8xekfH9oukXG1Z9o7zTZs8IeG533vrqfjUEhEHOCUhZB11GV2kSGujH78Z25/+W7eXl5fu5aUhc1WXuRcRjlIxcIDQ0lFWrVjF37lxycnJISEjgxRdfZMKECdW2nzNnDtOmTeO6664jNDSUP/zhDxw7dozAwMAGx2C1Wnn//fe5//776d+/P7169eLll19m7NixDX7M6kyaNInMzEyefPJJUlNT6d+/P8uWLXMU7oJZjHrhlOLNmzdz+eXnuudnzZoFwJQpU1i0aJHTjy1Nx75TuXxm7xW5smfNDb//EGxlEJcEUb09FF09+QVC34nm1ONdH1SbjABMHpHAq6sOcvxsIf/deZIbh3TycKAiYmcxapqz2oTk5OQQFhZGdnY2oaGh3g6nVvn5+cTFxfHiiy8ydepUb4cj4pT73t3KZztT+Um/jvxzclLNDV+/Co5vggkvwCW/8VyA9XVoJbx9PQS2gwf3g69/tc3+8e0BXvjqR3pGB/PlA5dWX7ArIg3m7Oe3+iUbadu2bbz33nscPHiQrVu38stfmrMLfvrTn3o5MhHnHDqdx+e7UgG4/8paakWyjpmJCBbo28Tf311GmwW2RVkVG/lV71fDEwgO8GXfqTyW/5BeYzsRcS8lIy7wt7/9jUGDBjlWal29ejWRkZHeDkvEKW+uPYxhwJW9o+gbW0vP457/mNcJoyCkic+KsvpA/4q1enZ9UGOzsDZ+/Gq4OXQ4f8WBGhc3FBH3UjLSSEOGDGHLli3k5eVx5swZkpOTa1wkTaSpOZtfwodbzNkkd43pWnvj3RUzyvrd4N6gXMW+MuyPX0Bxbo3Nfj26C/6+VralZLFBq7KKeIWSEZFW7N2NKRSV2ugXG8rwrrXMyMo6Bic2Axboc73H4muU2CEQ3g3KCmHfVzU2iwoJ5NahZvHqghWu3f9JRJyjZESklSops7F43REA7hqTWPt09OY0RGNnOa+2Ze9/a236mzHdsFhg5b7THDpddZM9EXEvJSMirdRnO0+SnltMdGgA1w6Irb1xcxuisetznXm9PxlKi2ps1jmiLVf0igLgX99V3XJBRNxLyYhIK2QYBq+vNrc4mDLSrJmoUVZK8xuisYu9CELjoDQfDn1ba9PJI8xC1g+2HKOgpMwT0YlIBSUjIq3QlqNn2ZOaQxs/H24b1rn2xs1xiMbOYoHeFb0jez+rtemlPTrQJaItuUVlfLLtpAeCExE7JSMirdC7G82hiImDYmjXtvoFwRx+WGZeN/W1RWpiH6r5cRmU19zjYbVaHNN8315/RNN8RTxIyUgTMnbsWGbMmFHjz+58Lmk9sgtK+XynucjZL+rqFcnPhGPfmbd7Vb8lQpPXeSS0CTd3Gk5ZV2vTW5LiaePnww9puWw6UvPmlCLiWkpGmrCPP/6Yp556yvFzS0ggVq1axcSJE4mNjcVisfDJJ584dd78+fNJTEwkMDCQpKQkVq9e3aA2Ap9sP0FxmY3eHUMYHN+u9sb7vwbDBtEDoF28R+JzOR9f6HWNebuOoZqwtn7cMMQs5l28/oibAxMROyUjTVh4eDghISHeDsOl8vPzGTRoEPPmzXP6nCVLljBjxgweeeQRtm3bxpgxY5gwYQIpKSn1aiNm4ep7FUM0vxjWue7dpX+sGKJprr0idvahmh8+gzqGX+xDNcm7T3E2v8TdkYkISkaq9eWXXzJ69GjatWtHREQE1113XaWda8eOHcv06dOZMWMG7du3Jzo6moULF5Kfn8+dd95JSEgI3bp144svvqj0uGPHjuW+++7jvvvuczz2o48+WuPY9Pk9IXfccQcrV67kpZdewmKxYLFYOHLkCABdunRh7ty5lc4dPHgwf/nLXwAzAbj99tsJDg4mJiaGF198scpzGYbB888/T9euXWnTpg2DBg3iww8/bNgLWIsJEybw9NNPc9NNNzl9zpw5c5g6dSp33XUXffr0Ye7cucTHx7NgwYJ6tRHYdiyLH9JyCfC1csOQuNoblxbBgeXm7eaejHS9HPyCIOcEnNxWa9N+sWH0iw2lpNzGpztUyCriCUpGqpGfn8+sWbPYtGkTy5cvx2q1cuONN2Kz2RxtFi9eTGRkJBs3bmT69Oncc8893HLLLYwcOZKtW7cyfvx4Jk+eTEFBQaXHXrx4Mb6+vmzYsIGXX36Zv//977z++ut1xvTSSy8xYsQI7r77blJTU0lNTSU+3rlu84ceeohvv/2WpUuX8vXXX7NixQq2bNlSqc2jjz7KW2+9xYIFC9i9ezczZ87kV7/6FStXrqzyeM888wzBwcG1Xlw1RFJSUsKWLVsYN25cpePjxo1j3bp1TrcR03sbzF6R6wbGEtbGr/bGR9aYU2JDYiBmsPuDcye/QOh+pXm7ltVY7W5JMldktS+VLyLu5evtAJqim2++udLPb7zxBlFRUezZs4f+/fsDMGjQIB599FEAZs+ezV//+lciIyO5++67AXjsscdYsGABO3fuZPjw4Y7Hio+P5+9//zsWi4VevXqxa9cu/v73vzvOq0lYWBj+/v60bduWjh07Ov275OXl8cYbb/D2229z9dVXA2ZC1KlTJ0eb/Px85syZwzfffMOIESMA6Nq1K2vWrOHVV1/lsssuq/SY06ZN49Zbb631eePi6vjW7aSMjAzKy8uJjq48pTQ6Opq0tDSn2wjkFpXymaNw1YlEdl9Fz17Pn4C1BXxv6Tke9n4K+7+Cy2fX2vT6wXH837K97DqRzQ9pOfTuWMsGgiLSaEpGqnHw4EH+/Oc/891335GRkeHoEUlJSXEkIwMHDnS09/HxISIiotIGefYPxvT0ytuSDx8+vNI4/YgRI3jxxRcpLy932+9SUlLiSDLArEXp1auX4+c9e/ZQVFTkSFbsSkpKGDJkSJXHDA8PJzy8ln1M3ODC2gbDMKocc6ZNa/bl92kUlpbTrUMQSQnta29sGOYGc3Cu+LO5617x/j65DXJP1bpmSniQP1f2jubL3Wl8uPk4j17X10NBirROLeDrjutNnDiRzMxMXnvtNTZs2MCGDRsA88PZzs+vche3xWKpdMz+IXj+0I67WK3WKnUnpaWlAE6tlWCP8fPPP2f79u2Oy549e6qtG/HkME1kZCQ+Pj5VejjS09MdCZ8zbcScRQNw45C4upO0tJ1mfYVfW0i81APReUBItLl5HsCB5Dqb/6xiqOaT7ScoLXf/v2OR1kw9IxfIzMxk7969vPrqq4wZMwaANWvWuOzxv/vuuyo/9+jRAx8fnzrP9ff3r7YHpUOHDqSmpjp+zsnJ4fBhc6nv7t274+fnx3fffUfnzuaaEmfPnmXfvn2O4Ze+ffsSEBBASkpKlSGZ6nhymMbf35+kpCSSk5O58cYbHceTk5P56U9/6nSb1u5UThHrDmYC8NPBTvxt7L0i3a4w6y1aih7jzJ6RfV/BkF/V2vSyXh2IDPYnI6+ElT+e5qq+SmxF3EXJyAXat29PREQECxcuJCYmhpSUFB5++GGXPf6xY8eYNWsWv/3tb9m6dSuvvPJKtbNbqtOlSxc2bNjAkSNHCA4OJjw8HKvVyhVXXMGiRYuYOHEi7du3589//rMjuQkODmbq1Kk89NBDREREEB0dzSOPPIL1vBqAkJAQHnzwQWbOnInNZmP06NHk5OSwbt06goODmTJlSqU4GjNMk5eXx4EDBxw/Hz58mO3btxMeHu5IlubNm8fSpUtZvtycyTFr1iwmT57M0KFDGTFiBAsXLiQlJYVp06Y5HseZNq3Zf3ecxDAgKaE98eFt6z5hf0XPQc/x7g3M03qMh5XPwcFvoawEfGtefdbPx8oNg+N4fc1hPthyTMmIiBspGbmA1Wrl/fff5/7776d///706tWLl19+mbFjx7rk8W+//XYKCwsZNmwYPj4+TJ8+nd/85jdOnfvggw8yZcoU+vbtS2FhIYcPH6ZLly7Mnj2bQ4cOcd111xEWFsZTTz3l6BkBeOGFF8jLy+P6668nJCSE3//+92RnZ1d67KeeeoqoqCieffZZDh06RLt27bjooov405/+5JLf227z5s1cfvnljp9nzZoFwJQpU1i0aBFgFqSeP5V60qRJZGZm8uSTT5Kamkr//v1ZtmwZCQkJ9WrTmtmHaG4YXMfuvGCuunqiYrZV96vcGJUXxA6BoA6QfxpS1kPX2nsCb07qxOtrDvPtD6fJLSolJLCOGUgi0iAWoxlswJCTk0NYWBjZ2dmEhjbfqvaxY8cyePDgKmuCiLjTgfQ8rpqzEl+rhY2PXEV4UB170ez6ED6aClH94N4WODV66T2w410YcR+M/79amxqGwVVzVnLwdD5zbh3ETRd1qrW9iFTm7Oe3ClhFWrj/VPSKXNqzQ92JCMCB/5nXPVpYr4hdz4r1aPZ/XWdTi8XCtQPN3iT7fj4i4npKRkRaMMMw+M92cxXRnzozRGOznUtGWtoQjV23K8DqCxn74MyhOptfNzAGgFX7T5NdWOru6ERaJSUjHrRixQoN0YhH7TieTcqZAtr6+3C1MwWYaTvNegq/IIgfXnf75igw7NzvZl/uvhY9o0PoGR1MablB8p5Tbg5OpHVSMiLSgn2121x75fLeUbT1d6Je3d4r0vWyWmeaNHvdrzCvD37rVPNrB9iHarRXjYg7KBkRaaEMw+Cr781kZHw/J7cQaOlDNHbdKpKRw6ugvO6hl2sHmq/f6v0ZZBdoqEbE1ZSMiLRQB9LzOJSRj7+Plct7daj7hMIsOLbRvG3fVK6l6jgI2oRDSS4c31xn8+5RIfTuGEKZzeCrPdrvSMTVlIyItFD2IZpR3SOcWx/j8EowyiGiB7Tv4t7gvM1qhW4V690c/MapU+yFrJ9pVo2IyykZaQDDMPjNb35DeHg4FouF7du3V3tMxJu+rEhGftJfQzTVsg/VOJmMTBhgJiPrD2aQW6ShGhFXUjLSAF9++SWLFi3is88+c6z2Wd2xxho7diwzZsxofMDVmD9/PomJiQQGBpKUlOTUxnZ1nbNq1SomTpxIbGwsFouFTz75xC2xS92Ony3g+xM5WC1wVR8nZtEYBhxcYd62f0i3dF0rekZOboWCM3U279YhmK6RQZSWG6zen+Hm4ERaFyUjDXDw4EFiYmIYOXIkHTt2xNfXt9pjTdWSJUuYMWMGjzzyCNu2bWPMmDFMmDCBlJSURp2Tn5/PoEGDmDdvnid+DanF17vNKagXdwknIjig7hPOHobsFLD6QcJIN0fXRITFQYfeYNjMQlYnXNknCoD/7dUUXxGXMpqB7OxsAzCys7M98nxFRUXG9OnTjQ4dOhgBAQHGqFGjjI0bNxqGYRhTpkwxAMclISGh2mOGYRgffPCB0b9/fyMwMNAIDw83rrzySiMvL88wDMOw2WzGc889ZyQmJhqBgYHGwIEDjQ8++MARw4WPCRiHDx92ye83bNgwY9q0aZWO9e7d23j44Ydddg5gLF26tNGxSsPc8s91RsIfPzPeWH3IuRM2vWEYj4caxpsT3BtYU/PFw+bv/Z/pTjVffzDDSPjjZ8aQJ782ysptbg5OpPlz9vNbPSPV+MMf/sBHH33E4sWL2bp1K927d2f8+PGcOXOGl156iSeffJJOnTqRmprKpk2bqj2WmprKL37xC37961+zd+9eVqxYwU033YRRsRXQo48+yltvvcWCBQvYvXs3M2fO5Fe/+hUrV64E4KWXXmLEiBHcfffdpKamkpqaSnx8fKU4n3nmGYKDg2u9XDiUUlJSwpYtWxg3blyl4+PGjWPduur3IWnIOeI9mXnFbD5iDjuMd7Ze5NAK87rrWLfE1GR1O2+9ESe26Rqa0J6wNn6cyS9hW8pZNwcn0no0aCxh/vz5vPDCC6SmptKvXz/mzp3LmDFjamz/zjvv8Pzzz7N//37CwsL4yU9+wt/+9jciIiIaHLi75Ofns2DBAhYtWsSECRMAeO2110hOTuaNN97goYceIiQkBB8fHzp2PPcf/YXHtm7dSllZGTfddJNj59gBAwY4nmPOnDl88803jBgxAoCuXbuyZs0aXn31VS677DLCwsLw9/enbdu2lZ7nfNOmTePWW2+t9feJi4ur9HNGRgbl5eVER1euI4iOjiYtrfopiw05R7xn5b7T2AzoFxtKXLs2dZ9gKz83TNHakpGEkeDjbw5RZR6AyB61Nvf1sTK2Vwf+s/0k/9ubztAu4R4KVKRlq3fPSH3rDdasWcPtt9/O1KlT2b17Nx988AGbNm3irrvuanTw7nDw4EFKS0sZNWqU45ifnx/Dhg1j7969Tj/OoEGDuPLKKxkwYAC33HILr732GmfPmt+k9uzZQ1FREVdffXWlXoy3336bgwcPOv0c4eHhdO/evdZLmzbVfxhZLJZKPxuGUeWYK84Rz1vx42kAxjqztgiYS8AXngX/EIi9yI2RNUH+QRB/iXnb3jtUhysrCoKXq25ExGXqnYzMmTOHqVOnctddd9GnTx/mzp1LfHw8CxYsqLb9d999R5cuXbj//vtJTExk9OjR/Pa3v2Xz5roXGvIG+zBKYz94fXx8SE5O5osvvqBv37688sor9OrVi8OHD2Oz2QD4/PPP2b59u+OyZ88ePvzwQ6efoyHDNJGRkfj4+FTp0UhPT6/S89GYc8Q7ym0Gq/bbk5Eo506yfwgnjgGfplt47TaJl5nXR+qeUQZwWc8O+Fot7E/PIyWzwI2BibQe9UpGGlI7MHLkSI4fP86yZcswDINTp07x4Ycfcu2119b4PMXFxeTk5FS6eEr37t3x9/dnzZo1jmOlpaVs3ryZPn361OuxLBYLo0aN4oknnmDbtm34+/uzdOlS+vbtS0BAACkpKVV6Ms6vC/H396e8vLzGx582bVqlZKa6y9ChQyud4+/vT1JSEsnJyZWOJycnM3Jk9bMoGnKOeMfO41lkFZQSEujLkPh2zp3UWutF7BIvNa8PrzZ3La5DWBs/Lq4YntGsGhHXqNfXoIbUDowcOZJ33nmHSZMmUVRURFlZGddffz2vvPJKjc/z7LPP8sQTT9QnNJcJCgrinnvu4aGHHiI8PJzOnTvz/PPPU1BQwNSpU51+nA0bNrB8+XLGjRtHVFQUGzZs4PTp0/Tp04eQkBAefPBBZs6cic1mY/To0eTk5LBu3TqCg4OZMmUKAF26dGHDhg0cOXKE4OBgwsPDsVrP5Y/h4eGEh9d/zHrWrFlMnjyZoUOHMmLECBYuXEhKSgrTpk1ztJk3bx5Lly5l+fLlTp+Tl5fHgQMHHD8fPnyY7du3O15HcT/7EM2YHpH4+jjxXaO0EI6uN2+31mQk7iJzl+LCM5C+GzoOqPOUK/tEsf5QJst/OMWvRyd6IEiRFq4+U3ROnDhhAMa6desqHX/66aeNXr16VXvO7t27jZiYGOP55583duzYYXz55ZfGgAEDjF//+tc1Pk9RUZGRnZ3tuBw7dsyjU3sLCwuN6dOnG5GRkVWm9hqGYfz97393TN+t6diePXuM8ePHO6YH9+zZ03jllVcc99tsNuOll14yevXqZfj5+RkdOnQwxo8fb6xcudLR5scffzSGDx9utGnTxqVTew3DMP7xj38YCQkJhr+/v3HRRRdVel7DMIzHH3+8yu9Y1znffvttlenIgDFlyhSXxS21u37eGiPhj58ZSzamOHfCgW/Mqa1/62UYtlY8VfX/3WS+Duv+4VTzg+m5RsIfPzO6/+lzI7+41M3BiTRfzk7ttRiGE/PZKpSUlNC2bVs++OADbrzxRsfxBx54gO3btzumpZ5v8uTJFBUV8cEHHziOrVmzhjFjxnDy5EliYmLqfN6cnBzCwsLIzs4mNDTU2XBFWpXMvGKG/t//MAzY8KcriQ4NrPuk5Mdh7VwY9Au48Z9uj7HJWvsSJD8GPSfAbe/X2dwwDEY/9y0nsgp5686LudzZ+hyRVsbZz+961Yw0pHagoKCg0tACmMWdcK5YVEQab/X+DAwDencMcS4RgXNFm/YiztbKXjdydC2Ul9XZ3GKxMLp7JABrtTS8SKPVezbNrFmzeP3113nzzTfZu3cvM2fOrFQ7MHv2bG6//XZH+4kTJ/Lxxx+zYMECDh06xNq1a7n//vsZNmwYsbGxrvtNRFq5FT+mA/WYRVOUAye3m7e7jHZPUM1Fx4EQGAbFOZC6w6lTRvcwk5E1B5SMiDRWvefxTZo0iczMTJ588knHhnDLli1zLOyVmppaac2RO+64g9zcXObNm8fvf/972rVrxxVXXMFzzz3nut9CpJWz2QxWVXxDd3p9kZTvwCiH9l2gXXydzVs0qw8kjIYfP4cjq6BTUp2njOoeicUCP6Tlkp5bRFSIk71RIlJFgxYVuPfee7n33nurvW/RokVVjk2fPp3p06c35KlExAnfn8zmTH4JwQG+XNS5vXMn2YdoWnuviF3ipWYycngVjJ5ZZ/PwIH/6xYby/Ykc1h7I4MYhnTwQpEjLpL1pRFqA9QczAbgkMRx/Xyf/WTuSkUvdFFUz46gbWQ9lJU6dMqqibmTN/kx3RSXSKigZEWkB1h8yPwxHdHNyv6ei7HO1EeoZMUX1gbaRUFYIJ5xbIXpMd3NIbM2B0yrIF2kEJSMizVxpuY1Nh81deod3dTIZOboeDBuEd4WwuLrbtwYWS+XVWJ0wtEt7AnytnMop5kB6nhuDE2nZlIyINHO7TmSTX1JOWBs/+sY4uQ6P6kWql1ix+7h9F+M6BPr5MCzRXAV5tab4ijSYkhGRZu78ehGr1cnNHI9U7L2kepHK7OutHN8IJc5tgudYb0RTfEUaTMmISDP3XX3rRQqzIG2neVs9I5WFd4XQOCgvgWMbnDrFXsS64fAZysrr3mhPRKpSMiLSjJWU2dh85CxQj2QkpaJeJKI7hNa9HUOrcn7dyBHn6kb6xIQSEuBLXnEZP6TlujE4kZZLyYhIM7bzeBaFpeWEB/nTMyrEuZMOq16kVl3qVzfiY7UwtIu5tsuGikJiEakfJSMizZi9XmR41/rUi9iTkTFuiqqZsxexntwGJflOnXJxRRHrJiUjIg2iZESkGXOsL+LslN7Cs5C2y7ytnpHqtesMYfFgK4NjG506ZViXimTkyBmtNyLSAEpGRJqp4rJythytZ73I0XWAARE9IKSj+4Jr7hIqdiE/utap5gM6heHvayUzv4RDGc71pojIOUpGRJqpHceyKS6zERkcQLcOwc6dZJ/Sm6ghmloljDKvj65zqnmArw+D49sBGqoRaQglIyLN1NYUs1dkaEJ7LBYn60VUvOocezJyfDOUFjl1in2oZuMRJSMi9aVkRKSZsg/RXJTQzrkTCs7Aqe/N2yperV1ENwiOhvJiOLHFqVPsK7FuUjIiUm9KRkSaIcMw2FbRM5KU0N65k1K+AwyI7AnBUe4LriWwWM6rG3FuqOaihPZYLXDsTCGp2YVuDE6k5VEyItIMHTtTSEZeCX4+FvrFhjl3Usp687rzCPcF1pI46kbWONU8OMDX8bfYqLoRkXpRMiLSDG1JMT/s+sWGEejn49xJSkbqx56MHNsI5aVOnXJxFw3ViDSEkhGRZmjr0SygHkM0JQXmIl4ACUpGnNKhN7RpD6UFkLrDqVOGJZp/j02Hz7ozMpEWR8mISDPkKF7t7GQycmKLuYhXSCy0S3BjZC2I1QqdK+pGjjg3VJOUYPaM7EvPJa+4zF2RibQ4SkZEmpn84jJ+SMsB6jGTxjFEM9wszhTndKnfeiMdQgKIa9cGw4Bdx7PdGJhIy6JkRKSZ2XE8C5sBsWGBxIS1ce4k+4epfYaIOMf+eqWsB1u5U6cMijeLWHccz3JTUCItj5IRkWZma8UQzRBn60XKy+D4JvO2ilfrJ3oA+IdAcc65NVrqMKhTOwB2HMtyX1wiLYySEZFmZmtKFgBJztaLnNoFJXkQEAZRfdwXWEvk42sObYHTQzWDKpaFVzIi4jwlIyLNiGEYjmXgL3K2Z+SovV7kErA6OQ1YzkmoXxHrgLgwrBY4mV1Eeo5zS8mLtHZKRkSakUMZ+WQVlBLga6VvTKhzJ6VUfKPXEE3D2PfxOboODKPO5kEBvvSICgFgh4pYRZyiZESkGbHP0OgXG4q/rxP/fA2jYhl4VLzaUDGDwbcNFJ6B0z84dcrATmYR604VsYo4RcmISDPy/QkzGRkQ5+QS8JkHIf80+ARA7BA3RtaC+fpD/MXm7aNrnTrFXjeyXXUjIk5RMiLSjHx/sqJnxNlkxD5EE5cEvgFuiqoVsC9+Zu9lqsPg84pYDSeGdkRaOyUjIs2EzWaw+4S52Fl/pzfHsw/RqF6kUewzapxMRnp1DMHf10pOURlHMgvcGJhIy6BkRKSZSDlTQG5xGf6+VnpEBzt30lEVr7pEp4vB4gPZxyDrWJ3N/Xys9I81C4w1xVekbkpGRJoJ+xBNn44h+Pk48U83Nw3OHgYsED/MvcG1dAHBEDPQvO1k74jqRkScp2REpJn4vmKIxvl6kYr1RTr2h0Anz5Ga2XuX7K9rHewrsWpGjUjdlIyINBO7K3pGnK4XcSx2piEal3AkI871jNin9+4+mUNZuc1dUYm0CEpGRJoBwzDYVd9pvVrszLXsRazpe6DwbJ3Nu0QEEeTvQ3GZjUMZ+W4OTqR5UzIi0gycyCokq6AUX6uFnh2dKF4tyoa0io3dlIy4RnAURHQHDDi2sc7mVquFPhWr5Np7tUSkekpGRJoBe71Iz+gQAnyd2F/m2CbAgPZdIDTGrbG1Ko4pvs7VjfSrmFGz52SOuyISaRGUjIg0A456kbj67kejJeBdqp51I31j7T0jSkZEaqNkRKQZqHe9iL14VYuduZY9GTmxBUrr3pG3X0Wx8Z7UHK3EKlILJSMiTZxhGI49aZya1ltWbH5YgupFXC28KwRFQXkJnNxWZ/Me0cH4Wi1kFZRyMrvu5EWktVIyItLEpecWk5FXgtUCfTo6MUxzchuUF0NQh4qCS3EZi6VedSMBvj50jzILjlU3IlIzJSMiTZz9Q6xbh2Da+DtRvOpYAn64+eEprpVg3zTPuSLWvipiFamTkhGRJm7fqVwAenYMce4Ee3Glilfdw9EzsgFsdS9m1rvi7/bjKSUjIjVRMiLSxP1YkYz0inYiGbHZ4Jg9GRnuxqhasegB4B8Mxdlwem+dzXtVDK39kJbr7shEmi0lIyJN3P5TeQD0dGan3vQ95oJn/sHQcaCbI2ulfHzNXXzh3JBYLew9I0cy8ikqLXdnZCLNlpIRkSbMZjPYn25+o+7hTM+IvY6h08Xmh6a4Rz3WG4kKCaBdWz9sBhxIz3NzYCLNk5IRkSbs2NkCikpt+PtaSQhvW/cJKdoczyMcdSN1JyMWi8UxxKahGpHqKRkRacL2VQzRdOsQjK9PHf9cDUOLnXlKp6Fg9YWc45CVUmdzRxFrmopYRaqjZESkCXPMpHGmXiTrKOSeND8k44a6ObJWzj8IYgaZt53oHVERq0jtlIyINGH7HcmIM/UiFR+KMYPB34khHWkcR91I3euN9HL0jCgZEamOkhGRJuxHx0waJ5IR+8wODdF4hj0ZOVp3MmLv2UrPLSa7sNSdUYk0S0pGRJqocpvBwdP1mNbrKF7VYmceYS9iPb0XCs7U2jQk0I+YsEBAM2pEqqNkRKSJOpqZT0mZjUA/K/Ht6xh2yc+AjH3mbS125hlBkRDZ07x9bGOdze171BxI11CNyIWUjIg0Ufbi1R5RIVitdewxY68X6dAb2oa7OTJxcEzxrXvxM3syYl/ETkTOUTIi0kTZp/X2qNcQjepFPKoei5/1iDLrfvZrmEakCiUjIk3UvnrNpLGvL6J6EY+yJyMntkJpYa1N7UmlakZEqlIyItJE7XN2g7ySfEjdYd5WvYhnte8CwR3BVmomJLXo3sFMRk5kFZJfXOaB4ESajwYlI/PnzycxMZHAwECSkpJYvXp1re2Li4t55JFHSEhIICAggG7duvHmm282KGCR1qC03MbhjHzAiWGa45vAVgahnaBdZw9EJw4Wy3l1I7VP8W0f5E9ksD+AY5aUiJjqnYwsWbKEGTNm8Mgjj7Bt2zbGjBnDhAkTSEmpeUnkW2+9leXLl/PGG2/w448/8t5779G7d+9GBS7Skh3JyKe03CDI34e4dm1qb2yvV9D6It5hHxpzom5ERawi1av3tp5z5sxh6tSp3HXXXQDMnTuXr776igULFvDss89Waf/ll1+ycuVKDh06RHi4WeXfpUuXxkUt0sLZ6wq6RwVjsdQxk8a+2JmGaLzD/rof2wC2crD61Ni0e1Qw3x06o54RkQvUq2ekpKSELVu2MG7cuErHx40bx7p11U9t+/TTTxk6dCjPP/88cXFx9OzZkwcffJDCwpqLvYqLi8nJyal0EWlNDlUM0XTtUMcQTXmpOUwDWuzMW6L7g38IFOdA+p5amyZGmn9P+xCciJjqlYxkZGRQXl5OdHR0pePR0dGkpaVVe86hQ4dYs2YN33//PUuXLmXu3Ll8+OGH/O53v6vxeZ599lnCwsIcl/j4+PqEKdLsHTpdkYxEBtXeMG0nlBZAYDtzjRHxPKsPxA8zb9cxVNO1g/n3tP99RcTUoALWC7uNDcOosSvZZrNhsVh45513GDZsGNdccw1z5sxh0aJFNfaOzJ49m+zsbMfl2LFjDQlTpNk6nGF24yd2qCMZse+L0nk4WDU5zmsc+9TUvviZPbk8nJmPzWa4OyqRZqNe/3tFRkbi4+NTpRckPT29Sm+JXUxMDHFxcYSFhTmO9enTB8MwOH78eLXnBAQEEBoaWuki0po4hmki6xim0WJnTUPCeTv4GjUnGZ3at8Xfx0pJmY0TWbWvSyLSmtQrGfH39ycpKYnk5ORKx5OTkxk5svrx6lGjRnHy5Eny8s4VbO3btw+r1UqnTp0aELJIy3Ymv4SsAnNn1y6RtexJYxha7KypiL0IrH6QmwpZNc8s9LFaSIgw/6aqGxE5p979urNmzeL111/nzTffZO/evcycOZOUlBSmTZsGmEMst99+u6P9bbfdRkREBHfeeSd79uxh1apVPPTQQ/z617+mTZs6piyKtEL2IZrYsEDa+tcy4S1jPxRkgm8gxAz2THBSPf+2EDvYvF3HeiOJ9qEaJSMiDvWe2jtp0iQyMzN58sknSU1NpX///ixbtoyEhAQAUlNTK605EhwcTHJyMtOnT2fo0KFERERw66238vTTT7vutxBpQQ6ednImjX1ztrih4Ovv5qikTp1HmDObjq6DQT+vsZn5dz3FIU3vFXGodzICcO+993LvvfdWe9+iRYuqHOvdu3eVoR0RqZ79G3NiXTNp7MWrWuysaUgYCeterrNnxF7Eekg9IyIOKr8XaWLs35i71jWTJuW8mTTiffGXmNcZ+yA/o8Zmmt4rUpWSEZEm5mhmAQBdausZyTkJWUfBYoVOwzwUmdSqbThE9TVv19I7Yu/xOpldSFFpuSciE2nylIyINCGGYZByxkxGEsJrmUljX8+i4wAI1NT3JsOJ9UbCg/wJDvDFMOD4WU3vFQElIyJNSmZ+CQUl5VgsENe+ltlm9pU+tQR802KfYl1LMmKxnJveezRTQzUioGREpEmx94rEhAYS4FvzhmuqF2mi7D0jaTuhOLfGZvZk5EjFkJxIa6dkRKQJSan4cIqvbYim8Cyc2m3e1mJnTUtYHLRLAMMGxzbW2CwhwqwbSVHPiAigZESkSXHUi0TUkoykbAAMiOgOwVGeCUycZ08QaylitdcDqWdExKRkRKQJsScjnWvrGbEvdqb9aJomRxFrLcmIvWfkjJIREVAyItKk2D+cah2msRdHJozyQERSb/aekROboay4+iYVPV/HzhRQVm7zVGQiTZaSEZEm5FhdPSMlBXBym3lbK682TRHdIagDlBWd+1tdoGNoIP6+VspsBiezijwcoEjTo2REpIkoKi0nLcf8YKoxGTmxGWxlEBJrFkpK02OxnJvlVMMUX6vV4vgbHz2jIlYRJSMiTcTxs4UYBgQH+BIeVMPGd44hmpHmh540TZ3rLmLtoum9Ig5KRkSaiGPn1YtYako0HMmIhmiaNMeMmg1gq37Jd3td0HEVsYooGRFpKs7NpKlh5dXyUnOLetDKq01dxwHgHwLF2ZC+p9omndpXJCNaEl5EyYhIU1HntN7UHVBaAG3aQ4feHoxM6s3qA/EVGxjWMMW3U8Vy/8fPqmdERMmISBNRZzJy9Lz1Raz6p9vk2YfSjq6t9m57MnIiSz0jIvofTaSJsC8F37liQayqDez70ahepFk4v4jVMKrcbR+mycgrobCk+roSkdZCyYhIE2AYRu09IzabFjtrbuKSwMcf8k7BmUNV7g5r40dIgC8AJ7I0VCOtm5IRkSYgI6+EwtJyLBaIa1dNAevpH6AoC/zaQsxAj8cnDeAXCLEXmbdrmOIbVzFUc0xFrNLKKRkRaQLsvSKxYW3w963mn6V9P5pOF4OPnwcjk0ZJqH2fGs2oETEpGRFpAs6tMVLDtF4N0TRP9r9XSvUrsTqKWJWMSCunZESkCai1XsQwzn2z1mJnzUv8MMBi1ozkplW5W9N7RUxKRkSagKMVM2kSqptJk3UUck+C1Q/ihno4MmmUwDDo2N+8Xc0+NeeSEfWMSOumZESkCbAP09g/nCqxf4jFDgH/GtYgkaarln1qVDMiYlIyItIE2Be+sn84VaL9aJq3WopY7clnRl4xRaVaa0RaLyUjIl5WbjNIyykCapjW61jsTPvRNEv2v9up76Ewq9JdYW38CHasNaLeEWm9lIyIeFl6bhHlNgNfq4UOIQGV78w9BZkHAAt0vsQr8UkjhURDeFfAgGMbK91lsVhUNyKCkhERrztZ8Y24Y1ggPlZL5TvtvSLR/cwN8qR5ctSNVC1itfeGaUaNtGZKRkS87GSWOUQTG1bbEI3qRZq1hIpkpJa6Ea01Iq2ZkhERL7P3jMS2C6x6p33HVxWvNm/2v9+JLVBSuQdEM2pElIyIeN25ZOSCnpGCM5D2vXm7yxgPRyUu1T4RQuPAVgrHN1W6SwufiSgZEfG6E/ZhmguTkZT1gAGRvSA4yvOBietYLNBltHn7yJpKd6lnRETJiIjX2XtGqkzrtX9oddF+NC2CfZ+aKsmI+XdPz9VaI9J6KRkR8bLU7BqGaY6sNq/t36ilebP/HU9srlQ30q6tH239fQBIzS7yRmQiXqdkRMSLCkrKOFtQCkDM+QWs59eLJCgZaRHCu0JILJSXVKobqbzWiOpGpHVSMiLiRfZpvSEBvoQG+p27w1Ev0tNcNEuaP9WNiNRIyYiIF9U4k8ZRL6JekRalhmTEXi+ktUaktVIyIuJFNa4xonqRlqmGuhH7EJ1qRqS1UjIi4kUns6uZ1qt6kZarhrqRmDAzGUnLUc+ItE5KRkS8qNphGtWLtFw11I10DDX//uoZkdZKyYiIF1U7TKN6kZatmmTE3jOSmlWEYRjeiErEq5SMiHiRIxk5f5M81Yu0bNXUjXSsSEYKS8vJKSzzVmQiXqNkRMRLDMOoWjOiepGWr5q6kUA/H8KD/AFIVd2ItEJKRkS8JDO/hJIyGxbLuW/GqhdpBWqsG9GMGmm9lIyIeIl9iCYqJAA/n4p/ivYPpwTtR9Oidam6T835dSMirY2SEREvqXYmjepFWocuY8zr8+pG7GuNpGVrmEZaHyUjIl5yIquWehElIy1beFcIialUNxITpum90nopGRHxEnvPiH0pcI6uAwyI6A4hHb0XmLhfpboRszfMXjOSlqNkRFofJSMiXnJuWm9F8erhleZ14mVeikg8KvFS8/qQ+Xe314zY3xcirYmSEREvsU/rjbH3jFR8KNF1rHcCEs+y/51PbIGiHMf7IDVbC59J66NkRMRL7IWKMWGBkHMSMn4ELKoXaS3adYb2iWCUw9G1jmGagpJycou18Jm0LkpGRLygrNzG6dxioKJW4PAq846YQdA23IuRiUd1rRiSO7SSNv4+tGvrB0CailillVEyIuIFmfkl2AzwsVqICA7QEE1rZf97V9QL2XtHVDcirY2SEREvsH/z7RAcgI+Fc8WrXVW82qp0qShiTd8Duacc07zVMyKtjZIRES+wT9+MDguEzAOQcwJ8/CF+uJcjE48KioCOA8zbh1c5tgXQWiPS2jQoGZk/fz6JiYkEBgaSlJTE6tWrnTpv7dq1+Pr6Mnjw4IY8rUiLcaoiGekYGgCHVpgH4y8B/7beC0q8wz6V+/AKYuxrjSgZkVam3snIkiVLmDFjBo888gjbtm1jzJgxTJgwgZSUlFrPy87O5vbbb+fKK69scLAiLYX9w8YsXtUQTavW9XLz+tBKMzkFTmpJeGll6p2MzJkzh6lTp3LXXXfRp08f5s6dS3x8PAsWLKj1vN/+9rfcdtttjBgxosHBirQUaY6eEb9zM2kSx3otHvGihBFg9YPsYyT6pAPqGZHWp17JSElJCVu2bGHcuHGVjo8bN45169bVeN5bb73FwYMHefzxx516nuLiYnJycipdRFoS+zBNL+MwFGVDQCjEDvFyVOIV/kHQ6WIAErLNfWqUjEhrU69kJCMjg/LycqKjoysdj46OJi0trdpz9u/fz8MPP8w777yDr6+vU8/z7LPPEhYW5rjEx8fXJ0yRJs/+YdMtd4t5oMto8HHu34e0QBVTfNufWg9AbnEZuUWlXgxIxLMaVMBqsVgq/WwYRpVjAOXl5dx222088cQT9OzZ0+nHnz17NtnZ2Y7LsWPHGhKmSJN1Ksdc8KxDxnfmAe1H07pV1Av5Hl1NWKD537J6R6Q1qddXscjISHx8fKr0gqSnp1fpLQHIzc1l8+bNbNu2jfvuuw8Am82GYRj4+vry9ddfc8UVV1Q5LyAggICAgPqEJtJs5BWXkVdcRgAltEndaB5U8WrrFpcE/sFQeIZRwWksK4oiNbuIHtEh3o5MxCPq1TPi7+9PUlISycnJlY4nJyczcuTIKu1DQ0PZtWsX27dvd1ymTZtGr1692L59O5dccknjohdphuz1IqMCDmEpK4LgaOjQ28tRiVf5+EGC+X/opb57APWMSOtS70HqWbNmMXnyZIYOHcqIESNYuHAhKSkpTJs2DTCHWE6cOMHbb7+N1Wqlf//+lc6PiooiMDCwynGR1uJUxYfMFQF7oQRziKaaYU5pZbqOhf1fM6R8BzBW03ulVal3MjJp0iQyMzN58sknSU1NpX///ixbtoyEhAQAUlNT61xzRKQ1s0/rHc4u84CGaAQcdUNd83fgT6mjB02kNbAYhmF4O4i65OTkEBYWRnZ2NqGhod4OR6RR5q84wMIvN7M1cBpWDJi5B8LivB2WeJthwIu9IO8Ut5X8iYAel/PWncO8HZVIozj7+a29aUQ87FR2EaOt35uJSFRfJSJisligm7lC9aXWnaRVzLgSaQ2UjIh4WFpOEZdZd5g/dNf2CHKeivfDZdYdGqaRVkXJiIiHpWUXcZnPTvOH7ld5NxhpWrpdgYGFPtZj+OWnUVxW7u2IRDxCyYiIh4Vk/UCUJYty3zbQWXs1yXnahptrjgCX+uwkXUM10kooGRHxoHKbwYCizQCUdR4NvlrcTyqzVPSWXWbd4Zh5JdLSKRkR8aCMvGIutZj1In69xtXRWlqlimRkjHUXp87meTkYEc9QMiLiQaczMhhq/REAaw/Vi0g14i4i3xpCmKUA2/HN3o5GxCOUjIh4UOmBFfhZykm1xkB4V2+HI02R1YejYeb6Iu1PrvJyMCKeoWRExIPaHlsBwA/BWsxKana64xgA4s+s93IkIp6hZETEUwyD6PQ1AJyIHOXlYKQpK064HIDOxT9CfqaXoxFxPyUjIp6SeZB2xakUG74UxGpKr9SsfcfO7LV1NlfpPfStt8MRcTslIyKecuB/AGyy9SKifbiXg5GmrGNoICttgwAwDiR7ORoR91MyIuIpFcnIStsgOoYFejkYacqiQgNYaRsIgLF/OdhsXo5IxL2UjIh4QmkhHDHrRVbaBhEdqmREahbg68PBwP7kGwFYC07DqV3eDknErZSMiHjC0bVQVkiqEc4+o5N6RqROEWEhrLP1N3/Yr6EaadmUjIh4wr6vAFhRPojgAD+CA3y9HJA0dR1DA1hRUTdif/+ItFRKRkTczTBg35cALLddRHSo9qORunUMC2R5+RDzh+ObID/DuwGJuJGSERF3S98LWSmUWwNYY+uvehFxSnRoIGlEcCKwB2DA/q+9HZKI2ygZEXG3fV8AcLzdxRQRQEclI+IE+/tkc0DFar0/fuHFaETcS8mIiLv9aA7R7AoyFzqLVvGqOMH+PvnGlmQeOPgNlBV7MSIR91EyIuJO+RnmeD+wxjoUQD0j4hT7+2RNficIjoaSPHNWlkgLpGRExJ32fw0Y0HEAPxSEAKhmRJxiT0YyC8oo7361ebCil02kpVEyIuJO9nH+nhM4lVMEoDVGxCnt2vrh72v+F30m7krz4L4vzNlZIi2MkhERdykrMcf5gfIePyE91xzv1zCNOMNisTjeKynthoFPAGSlwOkfvByZiOspGRFxl6NrzHH+4GgyQ/tQbjOwWiAy2N/bkUkzYU9GThZYIfFS8+CPy7wYkYh7KBkRcZe9n5nXPceTllsCQIeQAHx99M9OnGOfUXMqpwh6X2Me/OFzL0Yk4h76X1HEHWy2cx8afa7nVI6GaKT+Olas1puWXQS9rgUscGILZJ/wbmAiLqZkRMQdTmyGvDTwD4HES0mrKF6NUjIi9WCfeZWWUwQh0RBfsQCaekekhVEyIuIOe/9rXvccD74BnMqumEmjZETqoeP5wzQAfSaa1z/810sRibiHkhERVzOMc8lIn+sAHD0jmtYr9dHx/J4RgN7m+4kja6HgjJeiEnE9JSMirpa+B84eNqdiVixWZf9mqwXPpD7s75dTOcUYhgHhiRA9AIxy7VUjLYqSERFXs/eKdLsCAoKBigJENEwj9WNPRkrKbJwtKDUPVvS28cNnXopKxPWUjIi4mn1Kr/1Dg/OHaQK8EZE0U/6+ViKCzHVp7AmtY6jm4DdQku+lyERcS8mIiCudOQyndoHFB3pOAKCgpIzcojJAwzRSf+eGaiqSkeh+0D4RyorgwP+8GJmI6ygZEXGlvZ+a1wkjISgCOPeNNsjfh5BAP29FJs2UvejZUcRqsZybVbPnP16KSsS1lIyIuNL3H5vX/W50HLJ/iERrJo00gGOtEfswDUC/G8zrH7+EkgLPByXiYkpGRFwl8yCkbgeLFfr+1HE4XauvSiN0vHCYBiD2ImiXAKX5sP8rL0Um4jpKRkRcZc8n5nXipRAU6Ticpmm90gj2oue085MRi+Vc75u9N06kGVMyIuIq3y81r/vdVOmwvXtdyYg0RLXDNAD9K95n+7+G4lwPRyXiWkpGRFwhY785i8bqe664sIK9e92+6ZlIfVRZEt5xx0AI72bOqtmnoRpp3pSMiLjC7opeka6XQ9vwSndpKXhpDHvNyNmCUopKy8/dYbGc6x3RUI00c0pGRFyhmlk0do7VV8PaeDIiaSHC2vgR4Gv+V20vhnawv98OJENRtocjE3EdJSMijZW+F07vBR9/6H1tpbvKym2O7vVY9YxIA1gslnN1IxcO1UT1hcheUF4CPyzzQnQirqFkRKSxdn1gXne7Etq0q3RXem4xNgN8rRYiglUzIg1TZfdeO4sF+t9s3ra/D0WaISUjIo1hs8HOf5u3B95a5e7U82bS+FgtnoxMWhD7gnnpFyYjAAN+Zl4f+hZy0zwYlYjrKBkRaYyUdZB9DAJCodeEKnenZhcCENtOQzTScPaZWFWm9wJEdINOw8Cwwa4PPRyZiGsoGRFpjB3vm9d9fwp+VQtUU7PMD48YFa9KI9RYM2I3aJJ5vfN9D0Uk4lpKRkQaqrTw3EZlg35ebZOTFT0jMSpelUaoca0Ru343gdUP0nbBqd0ejEzENZSMiDTUj8ugOAfC4qHzyGqbnOsZUTIiDVdjAatd23DoOd68vUO9I9L8KBkRaagdS8zrgbeCtfp/SqkVHx4x7TRMIw0X7dgsrxjDMKpvNLBiqGbXB2Arr76NSBOlZESkIfJOw4H/mbcHVj9EA5CaVVHAqpoRaQR7MlJSZuNsQWn1jXqOh8B2kJsKh1d5LjgRF1AyItIQuz4Aoxxih0CHntU2KSmzcTrPXDEzRrNppBH8fa1EBPkDNcyoAfANOLci6473PBSZiGsoGRGpL8OArW+btwf/ssZmp3KKMAzw97ES3tbfQ8FJS3VuqKaGZARgyK/M6z3/gcKzHohKxDWUjIjU1/HN5vLvvm1gwC01Nkt17EkTiFULnkkj2WfU1FjEChCXBFH9zJ18teaINCNKRkTqa+si87rfDVWWfz9fqqb1igs51hqpaZgGzOXhL7rdvL1lsdmLJ9IMKBkRqY+inHM79Nr/06+BvWckVjNpxAXsGy3ak9waDbwVfALg1C5I3e7+wERcoEHJyPz580lMTCQwMJCkpCRWr15dY9uPP/6Yq6++mg4dOhAaGsqIESP46quvGhywiFft/hhKCyCiB3QeUWtT+0yajuoZERewJ7Uns2rpGQFzzZG+15u3tyx2c1QirlHvZGTJkiXMmDGDRx55hG3btjFmzBgmTJhASkpKte1XrVrF1VdfzbJly9iyZQuXX345EydOZNu2bY0OXsTj7P+5X3S72SVei5P2nhElI+IC9hlZJ+vqGYFzvXa7PoSSfDdGJeIa9U5G5syZw9SpU7nrrrvo06cPc+fOJT4+ngULFlTbfu7cufzhD3/g4osvpkePHjzzzDP06NGD//73v40OXsSj0r6Hk1vNZbcH/aLO5udqRjRMI40X5+gZKax54TO7LmMgvCuU5MLupR6ITqRx6pWMlJSUsGXLFsaNG1fp+Lhx41i3bp1Tj2Gz2cjNzSU8PLzGNsXFxeTk5FS6iHjdptfN697XQHCHOpvbCw21xoi4gn24r6i0loXP7CwWGDLZvL3pDTdHJtJ49UpGMjIyKC8vJzo6utLx6Oho0tLSnHqMF198kfz8fG699dYa2zz77LOEhYU5LvHx8fUJU8T1Cs/Czorl3y++u87mxWXlZOSVAFp9VVwjwNeHDiEBgNk7UqeLbjcLWU9uNaejizRhDSpgtVwwVm4YRpVj1Xnvvff4y1/+wpIlS4iKiqqx3ezZs8nOznZcjh071pAwRVxn+7tm4WpUX+gyus7m9l6RAF8r7dr6uTs6aSXs9UcnnElGgiKh/83m7Y0L3RiVSOPVKxmJjIzEx8enSi9Ienp6ld6SCy1ZsoSpU6fy73//m6uuuqrWtgEBAYSGhla6iHiNzQYbXzNvD7u7zsJVODfjIbZdG6cSdRFn2GfUpDqTjID5fgVzOnpeupuiEmm8eiUj/v7+JCUlkZycXOl4cnIyI0dWv4U6mD0id9xxB++++y7XXnttwyIV8ZYDyXD2MASGndsZtQ5a8EzcwTG9t7aFz84XdxF0uhhspbBlkfsCE2mkeg/TzJo1i9dff50333yTvXv3MnPmTFJSUpg2bRpgDrHcfvu5xaDee+89br/9dl588UWGDx9OWloaaWlpZGdnu+63EHGnDa+a10Mmg3+QU6fYFzzTTBpxpZj6DNPYDfuteb35TSivo/BVxEvqnYxMmjSJuXPn8uSTTzJ48GBWrVrFsmXLSEhIACA1NbXSmiOvvvoqZWVl/O53vyMmJsZxeeCBB1z3W4i4S8YBOLgcsMDFU50+TT0j4g7nT+91Wt+fQlAU5KbC3k/dFJlI4/g25KR7772Xe++9t9r7Fi1aVOnnFStWNOQpRJqG7+ab1z3Gmes2OCk1S9N6xfViG5KM+PrD0Dth5XOw/h/Q7yan6p5EPEl704jUJO80bH/HvD1yer1OPbf6qoZpxHXsyUh6bjElZTbnT7z4bvANhBNb4OhaN0Un0nBKRkRqsvFVcyv2uCSnpvOe7/jZAgDi2isZEdeJDPYnwNeKYTixYd75gjvA4F+at9fMdUtsIo2hZESkOsV556bzjnqgXt3a2YWl5BaVAdBJyYi4kMVicbynjp2pRzICMPI+sFjN2WFp37shOpGGUzIiUp2tb0NRFoR3g97X1evUY2fMXpGIIH/a+jeoLEukRp3atwXO9b45LbyrWcwKsO5lF0cl0jhKRkQuVF5qFvqBWSti9anX6fYPiU7hbV0dmQjx4RU9I/VNRsDs5QNzN9+s6ndaF/EGJSMiF9r1AeQcN6dDOrE774Xs3efxGqIRN4iv6Bmp9zANQOwQSLwMjHJYN8/FkYk0nJIRkfOVl8HK583bI+4Fv/pPzbX3jMSrZ0TcwP6+qvcwjd2YWeb1lkWQc9I1QYk0kpIRkfPtfN9c+r1tpFO781bn2Fl7z4iSEXE9R8/I2Qb0jIDZM9J5JJQXw+o5LoxMpOGUjIjYlZeaC0MBjJ4BAcENehh7Aat9bF/EleyzaU7nFlNUWl7/B7BY4PLZ5u2tiyFLu6KL9ykZEbHb/o5Z1BcUBUOdX/r9fIZhcLziG2sn9YyIG7Rr60dwgDlLq8FDNYmXQpcxUF4Cq190YXQiDaNkRASgrARW/c28PWYW+DcskcjML6GwtByLBWK1FLy4QaW1Rho6VAMwtqJ3ZNu/4OxRF0Qm0nBKRkTA7K7OPgbBHSHpjgY/jH2IpmNoIAG+9ZsSLOIsRxHrmQb2jAB0GWXWj9hKzxVti3iJkhGRomxY8ax5+9IHwa/htR4qXhVPsL+/UhqTjABc8Wfzevs7kLarkVGJNJySEZHVc6AgEyJ7QtKdjXqoIxn5ACREKBkR9+kSab6/jmQ2MhmJvxj63QgY8PWjYBiND06kAZSMSOt29ih8t8C8ffVT4NO45dvtyUiXyKDGRiZSoy4R5vvL/n5rlKv+Aj7+cGgF7E9u/OOJNICSEWndlj9prreQeCn0HN/ohztU8eGQqGRE3Mj+/jqaWUC5rZG9Ge27wCW/NW9//ai58J+IhykZkdbr+Gb4/kPAAuOertfOvDU5klnRMxKhZETcJ7ZdG/x9rJSU2ziZ1YgZNXZjHoQ24ZDxI2xd1PjHE6knJSPSOpWXwecVy2IP+gXEDGr0Q57NLyGroBQ4N6Yv4g4+VgudI+x1Iy4YqmnT7txU3+VPQd7pxj+mSD0oGZHWadPrkLoDAsPg6idc8pCHKz4UOoYG0ta/cbUnInWx974ddkXdCMDQX0PHAVCUBcl/ds1jijhJyYi0Pjkn4ZunzdtX/QWCo1zysOeKV9UrIu6XWPE+c1ky4uML180FLLDjPTi8yjWPK+IEJSPS+nw5G0pyodPFcNEdLnvYIypeFQ+yz9hyyYwau05DzR4SgM9mQVmx6x5bpBZKRqR12fc17PkELD5w3d/B6rp/Aocr1nxQMiKeYH+fNXqtkQtd+Zi5P1Pmflgz17WPLVIDJSPSehScgU+nm7eH32OOj7vQ4Yw8QDNpxDPsyUjKmQJKy22ue+A27eAnFSsSr3rerK0ScTMlI9J6LHsQ8tLMlVaveNSlD22zGRxMN7vLu3YIdulji1QnOiSQtv4+lNsMjrpiRs35+t8MfSaCrQyWToPSItc+vsgFlIxI67DrQ/j+I3N45sZ/Nmr/meocP1tIYWk5/j5WumgpePEAq9VCj+gQAH5My3Ptg1ssZjFrUAdI3wPf/p9rH1/kAkpGpOXLSYXPf2/evvQhiEty+VPsO5ULQLeoYHx99M9KPKNXtNkLZ3//uVRQJEx8yby97hU4us71zyFSQf9rSstWXgYf3WWunRAz2NyV1w1+rPgw6BmtIRrxnJ4VPSNuSUYAel8Lg38FGPDR3ZCf6Z7nkVZPyYi0bN8+DUfXgH8w3Pw6+Pi55Wn2O5KRELc8vkh13J6MgFnMGtEdco7Dx3eBrdx9zyWtlpIRabl+WAZr/m7e/uk8iOzhtqf68ZQ5Zq9kRDzJ/n47kllAcZmbkoTAULj1bfBtAwe/gVUvuOd5pFVTMiIt05nD8Mk08/Yl90C/G932VGXlNg6eNpORXkpGxIOiQwMIDfSl3GZw6LSLZ9RUeqJ+5ro8ACv+Cgf+577nklZJyYi0PIVn4d1boSgbOg2Dq59069MdPVNASZmNNn4+dGrv2lk6IrWxWCz06uiBoRqAwb+ApDsAAz74NaT/4N7nk1ZFyYi0LGUlsGQyZOyD0LiK7mV/tz7lvjTzQ6BHdDBWq8WtzyVyIfv03h/S3JyMAPzkOeg8Aoqz4d1bIC/d/c8prYKSEWk5DAP++wAcWW0WrN62BEJj3P60u05kA9A3JtTtzyVyIfv77vuK96Fb+QXCpHcgvCtkpcC7k6DExcvRS6ukZERajm+egh3vmgub3bLY5cu912TncfNDYGCndh55PpHzDap43+06kY1hGO5/wqAI+OWH0KY9nNwKH95p9kiKNIKSEWkZVj4Pq180b1/7N+hxlUee1jAMdh7PAmBgpzCPPKfI+Xp1DMHfx0pWQSkpZzzUSxHRDX7+HvgGwr4v4aOp5po+Ig2kZESavzVzzy1XPe7pc1uge8CRzAJyisoI8LU6CglFPMnf10qfWHOoZsdxDwzV2CWMgJ+/Az7+sPdTWPobrUEiDaZkRJovw4DVc+B/j5s/X/FnGDndoyHYe0X6xobip2XgxUsGVfTK7TyW5dkn7n6VWSRu9TX3fvr4NxqykQbR/57SPNls8OVsWP6E+fNlD7ttqffa7DhmfhMdpHoR8SJ7vdJOT/aM2PWaAD97qyIh+RDe+zkUu3jjPmnxlIxI81NWAh/fDRsWmD+PfwYun+2VUFQvIk2BvWfk+5PZlNs8UMR6ob7Xwy+WgF9bOLgcFk+E/AzPxyHNlpIRaV5yT8Hb15vfwKy+cNNrMOJ3XgmlqLScnRXTKYd0bu+VGEQAunYIJiTQl4KScnaf9ELvCJhF41P+C23CzVk2r10Babu8E4s0O0pGpPk4tgkWXgYp6yEg1FxHZOCtXgtn85GzlJTZiAkLpEtEW6/FIeJjtTC8awQAaw94cWfdTkNh6tfQvgtkHYXXr4adH3gvHmk2lIxI02ezwXf/hEXXQG4qRPaCu781i+e8aN1Bsxt6RLcILBatvCreNbKbmYzY35deE9nj3L/PskJzp99lD0FpoXfjkiZNyYg0bdkn4F83wpd/hPIS6DMR7l4Okd29HRlrD5rfQEd2i/RyJCLn3oebjpxx3w6+zmobDrf9G8b83vx540J49VI4uc27cUmTpWREmiabDbb9CxaMgEMrzO3Lr/kb3Pr/IMD763nkFJWyq6J41f6NVMSbekYHExnsT1Gpje0pWd4OB6w+cOVj8MuPILijuV/U61fBt8+ol0SqUDIiTU/qTnjrJ/Cf35k778ZeBNNWw7C7oYkMh3x3MBObAYmRQcS200694n0Wi4URFb0jaw80oZksPa6Ce9dDn+vBVgYrn4P5w2Hf196OTJoQJSPSdGSfMDe6W3gZHNsAfkFw9VNmQVxkD29HV8mXu9MAuKxnBy9HInKO/f1of382GW3DzcXRblkEIbFw9oi56+87t5hfPqTV8/V2ACLkpcPal2Dja1BebB7re4O5fkhYnFdDq05RaTnJu08BcO1A9+8KLOKsq/tG4+9jZd+pPPadyqVntPeHNB0sFuh3o1nYuvI5+G4B7P/avPS7CcY+DB16eTtK8RL1jIj3pP8A/7kP/t4P1s8zE5HOI+HOL+HWxU0yEQFYvT+D3OIyOoYGkqT1RaQJCWvjx6U9zaGaz3amejmaGgSEmHtI/W4j9P+ZeWz3x/CPYfDuJDi82tzqQVoVJSPiWaWFsPPf5gqN8y+Bbf/PnCXT6WL41Udw5zJzA64m7LOdJwG4ZkAMVmvTqGERsbP31n228yRGU/5Qj+gGP3sDpq2B3tcBFnMH4MXXwYJR5nT+gjPejlI8RMM04n6lRXDoW9j7GfzwX7MoFQAL9LkORkyHzpd4NURnnckv4auK8fjrBmmIRpqeq/pEE+Br5dDpfDYdOcuwxHBvh1S7jgPM3X8zDsB382H7O5C+25zOn/xn6DHOnNLfczy0UU9kS6VkRNwj+wQcWQM/fg77/wel+efuC+sMQ34Fg2+DdvHei7EB/t/6oxSV2ugfF8qQ+HbeDkekipBAP266KI73Nh5j4apDTT8ZsYvsDtfNgSv/DLs+hK1vQ9pO+OEz82L1hS6jzV6UxMvMovYmMrtOGk/JiDSezQZnDsKJLXBkNRxZC2cPV24T2gl6X2t+w0kYBdbmN0JYVFrO2+uPAPCbS7tp1VVpsu4a05X3Nx3jf3tPcfB0Ht06BHs7JOe1aW9O4x92t7m3zd7/mpf0PeaaQ4dWmO2CoiBhpJmgxA+DDn3A19+bkUsjKBkR5xkG5KaZiUfmAUj73vzmkvZ95Z4PAIsVYgZBtyvMbzKxQ5r9t5hF646QmV9CXLs2XNO/o7fDEalRtw7BXNUnmuQ9p/h78j7m3XaRt0NqmI4DzMvlf4LMg2YPyf5kOL4J8tNhzyfmBcDqB1G9oeMg85zI7hDRHcLizQXYpEmzGE26wsmUk5NDWFgY2dnZhIaGejuclslWbtZy5GdA7knISYWcE+ZeMDmpkJ0CmYeqJh12vm2gY/+KbypjIP4SCGw5f6uUzALGzV1JUamN5382kFuHNq/hJWl9vj+RzfXz1mAz4K07Luby3lHeDsl1yooremLXwtE15jLzRTXsVuzjD+0ToV1nCI01LyEx566DoyCwnXpV3MTZz28lI02JYZyb0maUm6sVlpea1xferva+UvMfaWmBOWulynUhlORDcQ4Unq24ZJmXYie3HbdYoV2CWQnfobfZ+9FxoDl+20K/fZSU2Zj8xgY2HD7DyG4RvHPXJRqikWbhmWV7WbjqELFhgfznvtF0CAnwdkjuYRiQlWL21KbuhFO7zR7cM4fPrV1UF/9gMylp0x7atDMvge3APwj82oBf24rr82+3Bd9AM+Gx+oKPr9lDY/UFnwuuHbf9zF5iixWouLZYmn3PcU3cmozMnz+fF154gdTUVPr168fcuXMZM2ZMje1XrlzJrFmz2L17N7GxsfzhD39g2rRpTj+f25KRz2bBD58D9gTAMG87XhKjhmM04ByjhnOaWC4YEHruW4PjG0SMWfMR0c1MRFrRN4hym8ED72/js52ptPX34fP7x5AYGeTtsEScUlhSzoSXVnEks4C+MaG895vhhLXx83ZYnmMrh+zjZmKSfRxyTpoXe49v7smK6cNN5f9hSw2JirWa++wJzHntqnu8Sj9e2OaCn6+bY9b2uZCzn9/1rhlZsmQJM2bMYP78+YwaNYpXX32VCRMmsGfPHjp37lyl/eHDh7nmmmu4++67+de//sXatWu599576dChAzfffHN9n961irIgr4ktm1wb6/lZt2/1P1+YtVfK5Ct+Dgwzs/9K3wLam8d9WtF/VHU4mVXIrH9v57tDZ/DzsfDPXyUpEZFmpY2/D4vuHMbP/rmOPak53PCPtcydNJhBrWUmmNUH2ieYl5rYh6iLss7rLT5b8XPWeb3L1fQ022+Xl1T0UJeZPdTn3y4vNXu6nVLxJdWwNfpXbxAvbmBY756RSy65hIsuuogFCxY4jvXp04cbbriBZ599tkr7P/7xj3z66afs3bvXcWzatGns2LGD9evXO/WcbusZOXvUfBPas0s4L9O0Z4yWWo7RgHOqe55qjlks5xINq6/5j6qFduM1JQUlZWw+cpbPd6aydNsJSspttPX3Yc6tg/mJilalmdpzMoe7Fm/iZHYRAFf1ieLGIZ0Y3jWc8CB/DTu6m2FUHk4/v3fckXzUdttW9Zzzj1dJXozKz+3sfe3iXb6Wi1t6RkpKStiyZQsPP/xwpePjxo1j3bp11Z6zfv16xo0bV+nY+PHjeeONNygtLcXPr+o38eLiYoqLz43z5eTk1CdMp310yIddJ6qrc6h++KSmvK2mbK62NM+o4ayazqn/c9Tv8Wt/bvfGWtPj19ZzWvNz1O9vVG4zyCkqI6ughMy8Ek5mF1aKc1iXcJ772UD1iEiz1jc2lC8euJTHPv2eT3ec5H970/nf3nQAggN8iQ9vS1RIAIF+Vvx9fQjwtRLga8V6QZJS6TvYBc9RW0KjXMeb6hiaOc/NF1np76VdOOqVjGRkZFBeXk50dHSl49HR0aSlVT/ckZaWVm37srIyMjIyiImpuorls88+yxNPPFGf0Bpk5b7TfLrjpNufR5qXmLBAxvSI5Jah8VzcpZksGCVSh7C2frz08yHcf2UP3tuQwop9pzmQnkdecRl7U3PY20S3shHPGdK5Pf3jwrzy3A1aZ+TCDNgwjDqy4qrtqztuN3v2bGbNmuX4OScnh/h410+lHNcvms7hbau9r7ZMvsa7ansN6n8KlhrOakhsNZ3TkO5Zb8fc0Oepvr2F0EBf2rf1p32QP4mRQYQHtZ4CXWl9unUI5tHr+vIo5kJ+x88WcuxMARl5xZSU2ygutVFcZqOotLxqj+IFPY/V9The2DlZY8+nNDk9ory3OF69kpHIyEh8fHyq9IKkp6dX6f2w69ixY7XtfX19iYiIqPacgIAAAgLcPwXtuoGxXDfQ7U8jItIkBfr50D0qmO5e/BASgXru2uvv709SUhLJycmVjicnJzNy5MhqzxkxYkSV9l9//TVDhw6ttl5EREREWpd6bxAya9YsXn/9dd5880327t3LzJkzSUlJcawbMnv2bG6//XZH+2nTpnH06FFmzZrF3r17efPNN3njjTd48MEHXfdbiIiISLNV75qRSZMmkZmZyZNPPklqair9+/dn2bJlJCSY87hTU1NJSUlxtE9MTGTZsmXMnDmTf/zjH8TGxvLyyy97f40RERERaRK0HLyIiIi4hbOf381vH3cRERFpUZSMiIiIiFcpGRERERGvUjIiIiIiXqVkRERERLxKyYiIiIh4lZIRERER8SolIyIiIuJVSkZERETEq+q9HLw32BeJzcnJ8XIkIiIi4iz753Zdi703i2QkNzcXgPj4eC9HIiIiIvWVm5tLWFhYjfc3i71pbDYbJ0+eJCQkBIvF4rLHzcnJIT4+nmPHjmnPGzfTa+0Zep09Q6+zZ+h19gx3vs6GYZCbm0tsbCxWa82VIc2iZ8RqtdKpUye3PX5oaKje6B6i19oz9Dp7hl5nz9Dr7Bnuep1r6xGxUwGriIiIeJWSEREREfGqVp2MBAQE8PjjjxMQEODtUFo8vdaeodfZM/Q6e4ZeZ89oCq9zsyhgFRERkZarVfeMiIiIiPcpGRERERGvUjIiIiIiXqVkRERERLyqxScj8+fPJzExkcDAQJKSkli9enWt7VeuXElSUhKBgYF07dqVf/7znx6KtHmrz+v88ccfc/XVV9OhQwdCQ0MZMWIEX331lQejbd7q+562W7t2Lb6+vgwePNi9AbYQ9X2di4uLeeSRR0hISCAgIIBu3brx5ptveija5qu+r/M777zDoEGDaNu2LTExMdx5551kZmZ6KNrmadWqVUycOJHY2FgsFguffPJJned4/LPQaMHef/99w8/Pz3jttdeMPXv2GA888IARFBRkHD16tNr2hw4dMtq2bWs88MADxp49e4zXXnvN8PPzMz788EMPR9681Pd1fuCBB4znnnvO2Lhxo7Fv3z5j9uzZhp+fn7F161YPR9781Pe1tsvKyjK6du1qjBs3zhg0aJBngm3GGvI6X3/99cYll1xiJCcnG4cPHzY2bNhgrF271oNRNz/1fZ1Xr15tWK1W46WXXjIOHTpkrF692ujXr59xww03eDjy5mXZsmXGI488Ynz00UcGYCxdurTW9t74LGzRyciwYcOMadOmVTrWu3dv4+GHH662/R/+8Aejd+/elY799re/NYYPH+62GFuC+r7O1enbt6/xxBNPuDq0Fqehr/WkSZOMRx991Hj88ceVjDihvq/zF198YYSFhRmZmZmeCK/FqO/r/MILLxhdu3atdOzll182OnXq5LYYWxpnkhFvfBa22GGakpIStmzZwrhx4yodHzduHOvWrav2nPXr11dpP378eDZv3kxpaanbYm3OGvI6X8hms5Gbm0t4eLg7QmwxGvpav/XWWxw8eJDHH3/c3SG2CA15nT/99FOGDh3K888/T1xcHD179uTBBx+ksLDQEyE3Sw15nUeOHMnx48dZtmwZhmFw6tQpPvzwQ6699lpPhNxqeOOzsFlslNcQGRkZlJeXEx0dXel4dHQ0aWlp1Z6TlpZWbfuysjIyMjKIiYlxW7zNVUNe5wu9+OKL5Ofnc+utt7ojxBajIa/1/v37efjhh1m9ejW+vi32n7tLNeR1PnToEGvWrCEwMJClS5eSkZHBvffey5kzZ1Q3UoOGvM4jR47knXfeYdKkSRQVFVFWVsb111/PK6+84omQWw1vfBa22J4RO4vFUulnwzCqHKurfXXHpbL6vs527733Hn/5y19YsmQJUVFR7gqvRXH2tS4vL+e2227jiSeeoGfPnp4Kr8Woz3vaZrNhsVh45513GDZsGNdccw1z5sxh0aJF6h2pQ31e5z179nD//ffz2GOPsWXLFr788ksOHz7MtGnTPBFqq+Lpz8IW+1UpMjISHx+fKhl2enp6lYzPrmPHjtW29/X1JSIiwm2xNmcNeZ3tlixZwtSpU/nggw+46qqr3Blmi1Df1zo3N5fNmzezbds27rvvPsD80DQMA19fX77++muuuOIKj8TenDTkPR0TE0NcXFylrdL79OmDYRgcP36cHj16uDXm5qghr/Ozzz7LqFGjeOihhwAYOHAgQUFBjBkzhqefflq91y7ijc/CFtsz4u/vT1JSEsnJyZWOJycnM3LkyGrPGTFiRJX2X3/9NUOHDsXPz89tsTZnDXmdwewRueOOO3j33Xc13uuk+r7WoaGh7Nq1i+3btzsu06ZNo1evXmzfvp1LLrnEU6E3Kw15T48aNYqTJ0+Sl5fnOLZv3z6sViudOnVya7zNVUNe54KCAqzWyh9bPj4+wLlv7tJ4XvksdFtpbBNgnzb2xhtvGHv27DFmzJhhBAUFGUeOHDEMwzAefvhhY/LkyY729ulMM2fONPbs2WO88cYbmtrrhPq+zu+++67h6+tr/OMf/zBSU1Mdl6ysLG/9Cs1GfV/rC2k2jXPq+zrn5uYanTp1Mn72s58Zu3fvNlauXGn06NHDuOuuu7z1KzQL9X2d33rrLcPX19eYP3++cfDgQWPNmjXG0KFDjWHDhnnrV2gWcnNzjW3bthnbtm0zAGPOnDnGtm3bHFOom8JnYYtORgzDMP7xj38YCQkJhr+/v3HRRRcZK1eudNw3ZcoU47LLLqvUfsWKFcaQIUMMf39/o0uXLsaCBQs8HHHzVJ/X+bLLLjOAKpcpU6Z4PvBmqL7v6fMpGXFefV/nvXv3GldddZXRpk0bo1OnTsasWbOMgoICD0fd/NT3dX755ZeNvn37Gm3atDFiYmKMX/7yl8bx48c9HHXz8u2339b6f25T+Cy0GIb6tkRERMR7WmzNiIiIiDQPSkZERETEq5SMiIiIiFcpGRERERGvUjIiIiIiXqVkRERERLxKyYiIiIh4lZIRERER8SolIyIiIuJVSkZERETEq5SMiIiIiFcpGRERERGv+v+mAMyzzd/PcAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "test.fitting_tool.plot_fit()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "79b0871e-66ba-4508-8ba4-5e21ea1d237c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.0" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "test.fitting_tool.model._forward(0,[0,10000,1,0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d32b918-caba-4e16-9840-e8393bf80761", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/lcls_tools/common/data_analysis/fit/method_base.py b/lcls_tools/common/data_analysis/fit/method_base.py new file mode 100644 index 00000000..58567fee --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/method_base.py @@ -0,0 +1,134 @@ +from abc import ABC, abstractmethod +import numpy as np +from matplotlib import pyplot as plt + + +class MethodBase(ABC): + """ + Base abstract class for all fit methods, which serves as the bare minimum skeleton code needed. + Should be used only as a parent class to all method models. + --------------------------- + Arguments: + param_names: list (list all of param names that the model will contain) + param_guesses: np.ndarray (array that contains a guess as to what + each param value is organized in the same order as param_names) + param_bounds: np.ndarray (array that contains the lower + and upper bound on for acceptable values of each parameter) + """ + + param_names: list = None + param_bounds: np.ndarray = None + + def __init__(self): + + self.init_values: dict = None + self.fitted_params_dict: dict = None + + @abstractmethod + def find_init_values(self) -> list: + ... + + @abstractmethod + def find_priors(self, data: np.ndarray) -> dict: + ... + + # TODO: move to plotting file + def plot_init_values(self): + init_values = np.array(list(self.init_values.values())) + """Plots init values as a function of forward and visually compares it to the initial distribution""" + fig, axs = plt.subplots(1, 1) + x = np.linspace(0, 1, len(self.profile_data)) + y_fit = self._forward(x, init_values) + axs.plot(x, self.profile_data, label="Projection Data") + axs.plot(x, y_fit, label="Initial Guess Fit Data") + axs.set_xlabel("x") + axs.set_ylabel("Forward(x)") + axs.set_title("Initial Fit Guess") + return fig, axs + + # TODO: move to plotting file + def plot_priors(self): + """Plots prior distributions for each param in param_names""" + num_plots = len(self.priors) + fig, axs = plt.subplots(num_plots, 1) + for i, (param, prior) in enumerate(self.priors.items()): + x = np.linspace(0, self.param_bounds[i][-1], len(self.profile_data)) + axs[i].plot(x, prior.pdf(x)) + axs[i].axvline( + self.param_bounds[i, 0], + ls="--", + c="k", + ) + axs[i].axvline(self.param_bounds[i, 1], ls="--", c="k", label="bounds") + axs[i].set_title(param + " prior") + axs[i].set_ylabel("Density") + axs[i].set_xlabel(param) + fig.tight_layout() + return fig, axs + + def forward(self, x: np.ndarray, params: dict) -> np.ndarray: + # TODO:test new usage + + params_list = np.array([params[name] for name in self.param_names]) + return self._forward(x, params_list) + + @staticmethod + @abstractmethod + def _forward(x: np.ndarray, params: np.ndarray) -> np.ndarray: + ... + + def log_prior(self, params: dict): + # TODO:test new usage + print(params) + params_list = np.array([params[name] for name in self.param_names]) + return self._log_prior(params_list) + + @abstractmethod + def _log_prior(self, params: np.ndarray): + ... + + def log_likelihood(self, x: np.ndarray, y: np.ndarray, params: dict): + # TODO:test new usage + params_list = np.array([params[name] for name in self.param_names]) + return self._log_likelihood(x, y, params_list) + + def _log_likelihood(self, x: np.ndarray, y: np.ndarray, params: np.ndarray): + # reducing error between data and fit + return -np.sum((y - self._forward(x, params)) ** 2) + + def loss(self, params, x, y, use_priors=False): + # TODO:implement using private functions _log_likelihood and _log_prior + # ML group way of iterating over priors/ reducing difference in fit and data + loss_temp = -self._log_likelihood(x, y, params) + if use_priors: + loss_temp = loss_temp - self._log_prior(params) + return loss_temp + + @property + def priors(self): + """ + Initial Priors store in a dictionary where the keys are the + complete set of parameters of the Model + """ + return self._priors + + @priors.setter + def priors(self, priors): + if not isinstance(priors, dict): + raise TypeError("Input must be a dictionary") + self._priors = priors + + @property + def profile_data(self): + """1D array typically projection data""" + return self._profile_data + + @profile_data.setter + def profile_data(self, profile_data): + if not isinstance(profile_data, np.ndarray): + raise TypeError("Input must be ndarray") + self._profile_data = profile_data + # should change these to private methods only called when profile data is set + self.find_init_values() + self.find_priors() + self.fitted_params_dict = {} diff --git a/lcls_tools/common/data_analysis/fit/methods.py b/lcls_tools/common/data_analysis/fit/methods.py new file mode 100644 index 00000000..fe4c832b --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/methods.py @@ -0,0 +1,101 @@ +import numpy as np +from scipy.stats import norm, gamma +from lcls_tools.common.data_analysis.fit.method_base import MethodBase + + +class GaussianModel(MethodBase): + """ + GaussianModel Class that finds initial parameter values for gaussian distribution + and builds probability density functions for the likelyhood a parameter + to be that value based on those initial parameter values. + Passing this class the variable profile_data automatically updates + the initial values and and probability density functions to match that data. + """ + + param_names: list = ["mean", "sigma", "amplitude", "offset"] + param_bounds: np.ndarray = np.array( + [ + [0.01, 1.0], + [0.01, 5.0], + [0.01, 1.0], + [0.01, 1.0], + ] + ) + + def find_init_values(self) -> dict: + """Fit data without optimization, return values.""" + + data = self._profile_data + x = np.linspace(0, 1, len(data)) + # init_fit = norm.pdf(data) + offset = data.min() + amplitude = data.max() - offset + + weighted_mean = np.average(x, weights=data) + weighted_sigma = np.sqrt(np.cov(x, aweights=data)) + + self.init_values = { + self.param_names[0]: weighted_mean, # data.mean() + self.param_names[1]: weighted_sigma, # data.std() + self.param_names[2]: amplitude, + self.param_names[3]: offset, + } + return self.init_values + + def find_priors(self) -> dict: + """Do initial guesses based on data and make distribution from that guess.""" + # Creating a gamma distribution around the inital amplitude. + # TODO: add to comments on why gamma vs. normal dist used for priors. + amplitude_mean = self.init_values["amplitude"] + amplitude_var = 0.05 + amplitude_alpha = (amplitude_mean**2) / amplitude_var + amplitude_beta = amplitude_mean / amplitude_var + amplitude_prior = gamma(amplitude_alpha, loc=0, scale=1 / amplitude_beta) + + # Creating a normal distribution of points around the inital mean. + mean_prior = norm(self.init_values["mean"], 0.1) + # TODO: remove hard coded numbers? + sigma_alpha = 2.5 + sigma_beta = 5.0 + sigma_prior = gamma(sigma_alpha, loc=0, scale=1 / sigma_beta) + + # Creating a normal distribution of points around initial offset. + offset_prior = norm(self.init_values["offset"], 0.5) + self.priors = { + self.param_names[0]: mean_prior, + self.param_names[1]: sigma_prior, + self.param_names[2]: amplitude_prior, + self.param_names[3]: offset_prior, + } + return self.priors + + def _forward(self, x: np.array, params: np.array): + # Load distribution parameters + # needs to be array for scipy.minimize + mean = params[0] + sigma = params[1] + amplitude = params[2] + offset = params[3] + return ( + np.sqrt(2 * np.pi) * amplitude * norm.pdf((x - mean) / sigma) + offset + ) # norm.pdf(x, loc=mean, scale=sigma) #( + + # TODO: remove when above is confirmed the same/below not needed. + # @staticmethod + # def _forward(x: np.ndarray, params_list: np.ndarray): + # amplitude = params_list[0] + # mean = params_list[1] + # sigma = params_list[2] + # offset = params_list[3] + # normal = norm() + # return ((np.sqrt(2 * np.pi)) * amplitude) * normal.pdf( + # (x - mean) / sigma + # ) + offset + + def _log_prior(self, params: np.ndarray) -> float: + return np.sum( + [ + prior.logpdf(params[i]) + for i, (key, prior) in enumerate(self.priors.items()) + ] + ) diff --git a/lcls_tools/common/data_analysis/fit/projection.py b/lcls_tools/common/data_analysis/fit/projection.py new file mode 100644 index 00000000..898125bf --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/projection.py @@ -0,0 +1,146 @@ +import numpy as np +import scipy.optimize +import scipy.signal +from matplotlib import pyplot as plt +from pydantic import BaseModel, ConfigDict +from lcls_tools.common.data_analysis.fit.method_base import MethodBase + + +class ProjectionFit(BaseModel): + """ + 1d fitting class that allows users to choose the model with which the fit + is performed, and if prior assumptions (bayesian regression) about + the data should be used when performing the fit. + Additionally there is an option to visualize the fitted data and priors. + -To perform a 1d fit, call fit_projection(projection_data={*data_to_fit*}) + ------------------------ + Arguments: + model: MethodBase (this argument is a child class object of method base + e.g GaussianModel & DoubleGaussianModel) + visualize_priors: bool (shows plots of the priors and init guess + distribution before fit) + use_priors: bool (incorporates prior distribution information into fit) + visualize_fit: bool (visualize the parameters as a function of the + forward function + from our model compared to distribution data) + """ + + # TODO: come up with better name + model_config = ConfigDict(arbitrary_types_allowed=True) + model: MethodBase + use_priors: bool = False + + def normalize(self, data: np.ndarray) -> np.ndarray: + """ + Normalize a 1d array by scaling and shifting data + s.t. data is between 0 and 1 + """ + data_copy = data.copy() + normalized_data = (data_copy - np.min(data)) / (np.max(data) - np.min(data)) + return normalized_data + + def unnormalize_model_params( + self, params_dict: dict, projection_data: np.ndarray + ) -> np.ndarray: + """ + Takes fitted and normalized params and returns them + to unnormalized values i.e the true fitted values of the distribution + """ + + projection_data_range = np.max(projection_data) - np.min(projection_data) + length = len(projection_data) + for key, val in params_dict.items(): + if "sigma" in key or "mean" in key: + true_fitted_val = val * length + elif "offset" in key: + true_fitted_val = val * projection_data_range + np.min(projection_data) + else: + true_fitted_val = val * projection_data_range + temp = {key: true_fitted_val} + params_dict.update(temp) + return params_dict + + def model_setup(self, projection_data=np.ndarray) -> None: + """sets up the model and init_values/priors""" + self.model.profile_data = projection_data + + def fit_model(self) -> scipy.optimize._optimize.OptimizeResult: + """ + Fits model params to distribution data and plots the fitted params + as a function of the model. + Returns optimizeResult object + """ + x = np.linspace(0, 1, len(self.model.profile_data)) + y = self.model.profile_data + + init_values = np.array( + [self.model.init_values[name] for name in self.model.param_names] + ) + + res = scipy.optimize.minimize( + self.model.loss, + init_values, + args=(x, y, self.use_priors), + bounds=self.model.param_bounds, + method="Powell" + ) + return res # TODO:optional argument to return fig,ax + + def fit_projection(self, projection_data: np.ndarray) -> dict: + """ + type is dict[str, float] + Wrapper function that does all necessary steps to fit 1d array. + Returns a dictionary where the keys are the model params and their + values are the params fitted to the data + """ + assert len(projection_data.shape) == 1 + fitted_params_dict = {} + normalized_data = self.normalize(projection_data) + self.model_setup(projection_data=normalized_data) + res = self.fit_model() + for i, param in enumerate(self.model.param_names): + fitted_params_dict[param] = (res.x)[i] + self.model.fitted_params_dict = fitted_params_dict.copy() + params_dict = self.unnormalize_model_params(fitted_params_dict, projection_data) + return params_dict + + def plot_fit(self, **kwargs): + """ + Plots the normalized fitted forward function against the normed data. + To plot unnormalized projection data, pass the unnormalized forward + function params and data as kwargs. + """ + fitted_params_dict = {} + projection_data = [] + if kwargs: + for key, val in kwargs.items(): + if key in self.model.fitted_params_dict: + fitted_params_dict[key] = val + else: + projection_data = val + x = np.linspace(0, len(projection_data), len(projection_data)) + else: + fitted_params_dict.update(self.model.fitted_params_dict) + projection_data = self.model.profile_data + x = np.linspace(0, 1, len(self.model.profile_data)) + + y = projection_data + fig, ax = plt.subplots() + y_fit = self.model.forward(x, fitted_params_dict) + ax.plot(x, y, label="data") + ax.plot(x, y_fit, label="fit") + + textstr = "\n".join( + [ + r"$\mathrm{%s}=%.2f$" % (key, val) + for key, val in fitted_params_dict.items() + ] + ) + ax.text( + 0.05, + 0.95, + textstr, + transform=ax.transAxes, + verticalalignment="top", + ) + return fig, ax diff --git a/lcls_tools/common/devices/screen.py b/lcls_tools/common/devices/screen.py index 868f00d6..121cc979 100644 --- a/lcls_tools/common/devices/screen.py +++ b/lcls_tools/common/devices/screen.py @@ -242,9 +242,11 @@ def _write_image_to_hdf5( # we update with original key if it isn't in our normal screen metadata # otherwise, prepend user_ to the key to retain all information. [ - dset.attrs.update({key: value}) - if key not in self.metadata - else dset.attrs.update({"user_" + key: value}) + ( + dset.attrs.update({key: value}) + if key not in self.metadata + else dset.attrs.update({"user_" + key: value}) + ) for key, value in extra_metadata.items() ] diff --git a/lcls_tools/common/devices/yaml/generate.py b/lcls_tools/common/devices/yaml/generate.py index 136b6327..1376a651 100644 --- a/lcls_tools/common/devices/yaml/generate.py +++ b/lcls_tools/common/devices/yaml/generate.py @@ -88,11 +88,11 @@ def _construct_information_from_element( item.strip() for item in element["Beampath"].split(",") if item ], "area": element["Area"], - "sum_l_meters": float( - np.format_float_positional(sum_l_meters, precision=3) - ) - if sum_l_meters - else None, + "sum_l_meters": ( + float(np.format_float_positional(sum_l_meters, precision=3)) + if sum_l_meters + else None + ), }, } diff --git a/lcls_tools/common/image_processing/image_processing.py b/lcls_tools/common/image_processing/image_processing.py index dc5f810b..22ed8f75 100644 --- a/lcls_tools/common/image_processing/image_processing.py +++ b/lcls_tools/common/image_processing/image_processing.py @@ -1,45 +1,65 @@ import numpy as np -import scipy.ndimage as snd - - -def fliplr(image): - """Flip over vertical axis""" - return np.fliplr(image) - - -def flipud(image): - """Flip over horizontal axis""" - return np.flipud(image) - - -def center_of_mass(image, sigma=5): - """Find center of mass, sigma sets threshold""" - return snd.center_of_mass(image > image.mean() + sigma * image.std()) - - -def average_image(images): - """If we can do things with an average image, do it!""" - return sum(images) / len(images) - - -def shape_image(image, x_size, y_size): - """Shape typical returned array, rows x columns so y x x""" - return image.reshape(y_size, x_size) - - -def x_projection(image, axis=0, subtract_baseline=True): - """Expects ndarray, return x projection""" - proj = np.sum(image, axis=0) - if subtract_baseline: - return proj - min(proj) - - return proj - - -def y_projection(image, subtract_baseline=True): - """Expects ndarray, return y projection""" - proj = np.sum(image, axis=1) - if subtract_baseline: - return proj - min(proj) - - return proj +from matplotlib import pyplot as plt +from scipy.ndimage import gaussian_filter +from pydantic import BaseModel, PositiveFloat, ConfigDict +from lcls_tools.common.image_processing.roi import ROI + + +class ImageProcessor(BaseModel): + """ + Image Processing class that allows for background subtraction and roi cropping + ------------------------ + Arguments: + roi: ROI (roi object either Circular or Rectangular), + background_image: np.ndarray (optional image that will be used in + background subtraction if passed), + threshold: Positive Float (value of pixel intensity to be subtracted + if background_image is None, default value = 0.0) + visualize: bool (plots processed image) + ------------------------ + Methods: + subtract_background: takes a raw image and does pixel intensity subtraction + process: takes raw image and calls subtract_background, passes to result + to the roi object for cropping. + """ + + model_config = ConfigDict(arbitrary_types_allowed=True) + roi: ROI = None + background_image: np.ndarray = None + threshold: PositiveFloat = 0.0 + + def subtract_background(self, raw_image: np.ndarray) -> np.ndarray: + """Subtract background pixel intensity from a raw image""" + if self.background_image is not None: + image = raw_image - self.background_image + else: + image = raw_image - self.threshold + return image + + def clip_image(self, image): + return np.clip(image, 0, None) + + def process(self, raw_image: np.ndarray) -> np.ndarray: + """Process image by subtracting background pixel intensity + from a raw image, crop, and filter""" + image = self.subtract_background(raw_image) + clipped_image = self.clip_image(image) + if self.roi is not None: + cropped_image = self.roi.crop_image(clipped_image) + else: + cropped_image = clipped_image + processed_image = self.filter(cropped_image) + return processed_image + + def filter(self, unfiltered_image: np.ndarray, sigma=5) -> np.ndarray: + # TODO: extend to other types of filters? Change the way we pass sigma? + filtered_data = gaussian_filter(unfiltered_image, sigma) + return filtered_data + + def plot_raw_and_processed_image(self, raw_image: np.ndarray, processed_image): + fig, ax = plt.subplots(2, 1) + c = ax[0].imshow(raw_image > 0, origin="lower") + rect = self.roi.get_patch() + ax[0].add_patch(rect) + ax[1].imshow(processed_image > 0, origin="lower") + fig.colorbar(c) diff --git a/lcls_tools/common/image_processing/roi.py b/lcls_tools/common/image_processing/roi.py new file mode 100644 index 00000000..c046f326 --- /dev/null +++ b/lcls_tools/common/image_processing/roi.py @@ -0,0 +1,109 @@ +from abc import ABC, abstractmethod +from typing import ( + List, +) +import numpy as np +from matplotlib import patches +from pydantic import BaseModel, PositiveFloat, Field + + +class ROI(BaseModel, ABC): + roi_type: str + center: List[PositiveFloat] + + @abstractmethod + def crop_image(self, img, **kwargs) -> np.ndarray: + """crop image using ROI""" + pass + + @abstractmethod + def get_patch(self): + pass + + +class CircularROI(ROI): + """ + Define a circular region of interest (ROI) for an image, cropping pixels outside a + bounding box around the ROI and setting pixels outside the boundary to a fill + value (usually zero). + """ + + roi_type: str = Field(default="circular", frozen=True) + radius: PositiveFloat + + @property + def bounding_box(self): + return [ + int(self.center[0] - int(self.radius)), + int(self.center[1] - int(self.radius)), + int(self.radius * 2), + int(self.radius * 2), + ] + + def crop_image(self, img, **kwargs) -> np.ndarray: + fill_value = kwargs.get("fill_value", 0.0) + img = self.fill_value_outside_circle(img, self.center, self.radius, fill_value) + bbox = self.bounding_box + img = img[bbox[1]: bbox[1] + bbox[3], bbox[0]: bbox[0] + bbox[2]] + return img + + def fill_value_outside_circle( + self, + img: np.ndarray, + center: List[PositiveFloat], + radius: PositiveFloat, + fill_value: float, + ): + height, width = img.shape + for y in range(height): + for x in range(width): + distance = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2) + if distance > radius: + img[y, x] = fill_value + return img + + def get_patch(self): + return patches.Circle( + tuple(self.center), self.radius, facecolor="none", edgecolor="r" + ) + + +class RectangularROI(ROI): + """ + Define a rectangular region of interest (ROI) for an image, cropping pixels outside + the ROI. + """ + + # set roi_type : str = 'rectangular' and freeze it. + roi_type: str = Field(default="rectangular", frozen=True) + xwidth: int + ywidth: int + + @property + def bounding_box(self): + return [ + int(self.center[0] - int(self.xwidth / 2)), + int(self.center[1] - int(self.ywidth / 2)), + int(self.xwidth), + int(self.ywidth), + ] + + def crop_image(self, img, **kwargs) -> np.ndarray: + x_size, y_size = img.shape + if self.xwidth > x_size or self.ywidth > y_size: + raise ValueError( + f"must specify ROI that is smaller than the image, " + f"image size is {img.shape}" + ) + bbox = self.bounding_box + img = img[bbox[1]: bbox[1] + bbox[3], bbox[0]: bbox[0] + bbox[2]] + return img + + def get_patch(self): + return patches.Rectangle( + (self.bounding_box[0], self.bounding_box[1]), + self.bounding_box[2], + self.bounding_box[3], + facecolor="none", + edgecolor="r", + ) diff --git a/lcls_tools/common/io.py b/lcls_tools/common/io.py index 73fb61d6..a88b10af 100644 --- a/lcls_tools/common/io.py +++ b/lcls_tools/common/io.py @@ -8,20 +8,17 @@ def __init__(self): pass def write( - self, - measurement_data: dict, - measurement_obj: Measurement, - filename: Optional[str] = None + self, + measurement_data: dict, + measurement_obj: Measurement, + filename: Optional[str] = None, ): """ Write data to h5file """ raise NotImplementedError - def read( - self, - filename: Optional[str] = None - ): + def read(self, filename: Optional[str] = None): """ Read data from h5file """ diff --git a/lcls_tools/common/measurements/measurement.py b/lcls_tools/common/measurements/measurement.py index 3917c3d1..cb99a50f 100644 --- a/lcls_tools/common/measurements/measurement.py +++ b/lcls_tools/common/measurements/measurement.py @@ -1,15 +1,28 @@ from abc import ABC, abstractmethod from typing import Optional - from pydantic import BaseModel, DirectoryPath class Measurement(BaseModel, ABC): + """ + Base abstract class for all meaurements, which serves as the bare minimum skeleton code needed. + Should be used only as a parent class to all performable measurements. + --------------------------- + Arguments: + name: str (name of measurement performed) + save_data: bool = True (specifies whether or not to dump data and meta data to h5py file) + save_location: Optional[DirectoryPath] = None + (optional argument that is a path to the directory you wish to save the data) + --------------------------- + Methods: + measure: abstractmethod to be implemented in all children classes wear measurement is performed + """ + name: str save_data: bool = True save_location: Optional[DirectoryPath] = None @abstractmethod def measure(self, **kwargs) -> dict: - """ Implements a measurement and returns a dictionary with the results""" + """Implements a measurement and returns a dictionary with the results""" raise NotImplementedError diff --git a/lcls_tools/common/measurements/screen_beam_profile_measurement.py b/lcls_tools/common/measurements/screen_beam_profile_measurement.py new file mode 100644 index 00000000..22e4c44a --- /dev/null +++ b/lcls_tools/common/measurements/screen_beam_profile_measurement.py @@ -0,0 +1,99 @@ +from lcls_tools.common.devices.screen import Screen +from lcls_tools.common.image_processing.image_processing import ImageProcessor +from lcls_tools.common.data_analysis.fit.projection import ProjectionFit +from lcls_tools.common.measurements.measurement import Measurement +import numpy as np + + +class ScreenBeamProfileMeasurement(Measurement): + """ + ScreenBeamProfileMeasurement class that allows for background subtraction and roi cropping + ------------------------ + Arguments: + name: str (name of measurement default is beam_profile), + device: Screen (device that will be performing the measurement), + image_processor: ImageProcessor () + fitting_tool: ProjectionFit () + fit_profile: bool = True () + ------------------------ + Methods: + single_measure: measures device and returns raw and processed image + measure: does multiple measurements and has an option to fit the image profiles + """ + + name: str = "beam_profile" + device: Screen + image_processor: ImageProcessor + fitting_tool: ProjectionFit + fit_profile: bool = True + # charge_window: Optional[ChargeWindow] = None + + def measure(self, n_shots: int = 1) -> dict: + """ + Measurement function that takes in n_shots as argumentment, where n_shots is the number of + image profiles (is this what I want to call it?) + we would like to measure. Invokes single_measure per shot + Then if fit_profile = True, fits the x and y projections of each image using the provided + fitting_tool. Results are stored in list of dictionaries where each element of the list + is results for the image fitting process stored as a dictionary. + ---- + Currently under work, what do we want to do with the images? + Needs dump_controller + + """ + images = [] + while len(images) < n_shots: + measurement = self.single_measure() + if len(measurement): + images += [measurement] + + results = {"images": images} + if self.fit_profile: + final_results = [] + for i, ele in enumerate(images): + temp = {} + temp["raw_image"] = ele["raw_image"] + temp["processed_image"] = ele["processed_image"] + + projection_x = np.array(np.sum(ele["processed_image"], axis=0)) + projection_y = np.array(np.sum(ele["processed_image"], axis=1)) + + for key, param in self.fitting_tool.fit_projection( + projection_x + ).items(): + key_x = key + "_x" + temp[key_x] = param + + for key, param in self.fitting_tool.fit_projection( + projection_y + ).items(): + key_y = key + "_y" + temp[key_y] = param + final_results += [temp] + + # what should I do with results now? + results["fits"] = final_results + + # no attribute dump controller + if self.save_data: + pass + # self.dump_controller.dump_data_to_file(final_results, self) + + return final_results + + def single_measure(self) -> dict: + """ + Function that grabs a single image from the device class + (typically live beam images) and passes it to the + image processing class for processing (subtraction and cropping) + returns a dictionary with both the raw and processed dictionary + """ + # measure profiles + # get raw data + raw_image = self.device.image + # get ICT measurements and return None if not in window + # if self.charge_window is not None: + # if not self.charge_window.in_window(): + # return {} + processed_image = self.image_processor.process(raw_image) + return {"raw_image": raw_image, "processed_image": processed_image} diff --git a/tests/datasets/fit/make_test_guassian.py b/tests/datasets/fit/make_test_guassian.py new file mode 100644 index 00000000..7ad4a63e --- /dev/null +++ b/tests/datasets/fit/make_test_guassian.py @@ -0,0 +1,15 @@ +import numpy as np +from scipy.stats import norm + +# import matplotlib.pyplot as plt + +# Generating 1000 point gaussian distribution +data = np.random.normal(size=1000) +fit = norm.pdf(data) + +# Plotting distribution +# plt.plot(data,fit, '.') +# plt.show() + +# Saving distribution +np.save("test_gaussian.npy", data) diff --git a/tests/datasets/fit/test_gaussian.npy b/tests/datasets/fit/test_gaussian.npy new file mode 100644 index 00000000..a2026819 Binary files /dev/null and b/tests/datasets/fit/test_gaussian.npy differ diff --git a/tests/datasets/h5py/test_image.h5 b/tests/datasets/h5py/test_image.h5 new file mode 100644 index 00000000..601a01fc Binary files /dev/null and b/tests/datasets/h5py/test_image.h5 differ diff --git a/tests/datasets/images/numpy/test_roi_fill_value_image.npy b/tests/datasets/images/numpy/test_roi_fill_value_image.npy new file mode 100644 index 00000000..9a28bf7f Binary files /dev/null and b/tests/datasets/images/numpy/test_roi_fill_value_image.npy differ diff --git a/tests/datasets/images/numpy/test_roi_image.npy b/tests/datasets/images/numpy/test_roi_image.npy new file mode 100644 index 00000000..d9fe6e67 Binary files /dev/null and b/tests/datasets/images/numpy/test_roi_image.npy differ diff --git a/tests/unit_tests/lcls_tools/common/data_analysis/fit/test_methods.py b/tests/unit_tests/lcls_tools/common/data_analysis/fit/test_methods.py new file mode 100644 index 00000000..4d5b3bf2 --- /dev/null +++ b/tests/unit_tests/lcls_tools/common/data_analysis/fit/test_methods.py @@ -0,0 +1,97 @@ +from lcls_tools.common.data_analysis.fit.methods import GaussianModel +import numpy as np +import unittest +import os + + +class TestGaussianModel(unittest.TestCase): + def setUp(self) -> None: + self.data_location = "./tests/datasets/fit/" + self.data_filename = os.path.join(self.data_location, "test_gaussian.npy") + self.data = np.load(self.data_filename) + self.gaussian_model = GaussianModel() + self.gaussian_model.profile_data = self.data + # Decimal place to check calculation errors to + self.decimals = 2 + return super().setUp() + + @unittest.skip + def test_find_init_values(self): + init_dict = self.gaussian_model.find_init_values() + self.assertIsNotNone(init_dict) + # TODO: pick up test here, the init dict is not as expected + self.assertAlmostEqual(init_dict["mean"], -0.07874, places=self.decimals) + self.assertAlmostEqual(init_dict["sigma"], 1.0, places=self.decimals) + return + + @unittest.skip + def test_find_priors(self): + priors_dict = self.gaussian_model.find_priors() + self.assertIsNotNone(priors_dict) + return # TODO: flesh out this test + + @unittest.skip + def test_forward(self): + init_dict = self.gaussian_model.find_init_values() + fit = self.gaussian_model.forward(self.data, init_dict) + # TODO: better way to check this? w/o hardcode? + self.assertAlmostEqual(fit.max(), 0.4, places=self.decimals) + self.assertAlmostEqual(fit.min(), 0.0, places=self.decimals) + return + + # def super_gaussian(self, x, amp, mu, sig, P, offset): + # """Super Gaussian Function""" + # """Degree of P related to flatness of curve at peak""" + # return amp * np.exp((-abs(x - mu) ** (P)) / (2 * sig ** (P))) + offset + # def double_gaussian(self, x, amp, mu, sig, amp2, nu, rho, offset): + # return ( + # amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) + # + amp2 * np.exp(-np.power(x - nu, 2.0) / (2 * np.power(rho, 2.0))) + # ) + offset + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_gaussian(self): + # # Test that the fitting tool can fit each type of gaussian distribution + # x_data = np.arange(500) + # # generated data for pure gaussian + # test_params = [3, 125, 45, 1.5] + # y_data = self.gaussian(x_data, *test_params) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # y_test = y_data + y_noise + # fitting_tool = FittingTool(data=y_test) + # fits = fitting_tool.get_fit() + # self.assertIsInstance(fits, dict) + # for key, val in fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # self.assertLessEqual(val["rmse"], 0.4) + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_super_gaussian(self): + # x_data = np.arange(500) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # test_params_super_gauss = [4, 215, 75, 4, 1] + # y_data_super_gauss = self.super_gaussian(x_data, *test_params_super_gauss) + # y_test_super_gauss = y_data_super_gauss + y_noise + # super_gauss_fitting_tool = FittingTool(data=y_test_super_gauss) + # s_fits = super_gauss_fitting_tool.get_fit() + # self.assertIsInstance(s_fits, dict) + # for key, val in s_fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # if key == "super_gaussian": + # self.assertLessEqual(val["rmse"], 0.4) + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_double_gaussian(self): + # # generated data for super gaussian + # x_data = np.arange(500) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # test_params_double_gauss = [2, 100, 25, 10, 240, 25, 2] + # y_data_double_gauss = self.double_gaussian(x_data, *test_params_double_gauss) + # y_test_double_gauss = y_data_double_gauss + 3 * y_noise + # double_gauss_fitting_tool = FittingTool(data=y_test_double_gauss) + # d_fits = double_gauss_fitting_tool.get_fit() + # self.assertIsInstance(d_fits, dict) + # for key, val in d_fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # if key == "double_gaussian": + # self.assertLessEqual(val["rmse"], 0.8) diff --git a/tests/unit_tests/lcls_tools/common/data_analysis/fit/test_projection.py b/tests/unit_tests/lcls_tools/common/data_analysis/fit/test_projection.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py b/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py deleted file mode 100644 index 028f250b..00000000 --- a/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py +++ /dev/null @@ -1,71 +0,0 @@ -from lcls_tools.common.data_analysis.fitting_tool import FittingTool -import numpy as np -import unittest - - -class TestFitTool(unittest.TestCase): - def gaussian(self, x, amp, mu, sig, offset): - return amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) + offset - - def super_gaussian(self, x, amp, mu, sig, P, offset): - """Super Gaussian Function""" - """Degree of P related to flatness of curve at peak""" - return amp * np.exp((-abs(x - mu) ** (P)) / (2 * sig ** (P))) + offset - - def double_gaussian(self, x, amp, mu, sig, amp2, nu, rho, offset): - return ( - amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) - + amp2 * np.exp(-np.power(x - nu, 2.0) / (2 * np.power(rho, 2.0))) - ) + offset - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_gaussian(self): - # Test that the fitting tool can fit each type of gaussian distribution - x_data = np.arange(500) - - # generated data for pure gaussian - test_params = [3, 125, 45, 1.5] - y_data = self.gaussian(x_data, *test_params) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - y_test = y_data + y_noise - - fitting_tool = FittingTool(data=y_test) - fits = fitting_tool.get_fit() - self.assertIsInstance(fits, dict) - for key, val in fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - self.assertLessEqual(val["rmse"], 0.4) - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_super_gaussian(self): - x_data = np.arange(500) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - test_params_super_gauss = [4, 215, 75, 4, 1] - y_data_super_gauss = self.super_gaussian(x_data, *test_params_super_gauss) - y_test_super_gauss = y_data_super_gauss + y_noise - super_gauss_fitting_tool = FittingTool(data=y_test_super_gauss) - s_fits = super_gauss_fitting_tool.get_fit() - self.assertIsInstance(s_fits, dict) - for key, val in s_fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - if key == "super_gaussian": - self.assertLessEqual(val["rmse"], 0.4) - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_double_gaussian(self): - # generated data for super gaussian - x_data = np.arange(500) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - test_params_double_gauss = [2, 100, 25, 10, 240, 25, 2] - y_data_double_gauss = self.double_gaussian(x_data, *test_params_double_gauss) - y_test_double_gauss = y_data_double_gauss + 3 * y_noise - double_gauss_fitting_tool = FittingTool(data=y_test_double_gauss) - d_fits = double_gauss_fitting_tool.get_fit() - self.assertIsInstance(d_fits, dict) - for key, val in d_fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - if key == "double_gaussian": - self.assertLessEqual(val["rmse"], 0.8) diff --git a/tests/unit_tests/lcls_tools/common/image_processing/test_image_processing.py b/tests/unit_tests/lcls_tools/common/image_processing/test_image_processing.py index 5f6b8e57..3276574e 100644 --- a/tests/unit_tests/lcls_tools/common/image_processing/test_image_processing.py +++ b/tests/unit_tests/lcls_tools/common/image_processing/test_image_processing.py @@ -1,88 +1,63 @@ import unittest -import os import numpy as np -import lcls_tools.common.image_processing as ip -from lcls_tools.common.matlab2py.mat_image import MatImage - -CAMERA = "CAMR:LGUN:210" - - -class ImageProcessingTest(unittest.TestCase): - data_location: str = "/tests/datasets/images/matlab/" - - def setUp(self): - self.file = os.path.join(self.data_location, "test_image.mat") - try: - if not os.path.isfile(self.file): - raise FileNotFoundError(f"Could not find {self.file}, aborting test.") - except FileNotFoundError: - self.skipTest("Invalid dataset location") - """Use test image""" - self.MI = MatImage() - self.MI.load_mat_image("test_image.mat") - - def test_fliplr(self): - """Test that fliplr does the right thing""" - col_init = self.MI.image[:, 0] - col_final = ip.fliplr(self.MI.image)[:, -1] - self.assertEqual(np.array_equal(col_init, col_final), True) - - def test_flipud(self): - """Test that flipud does the right thing""" - row_init = self.MI.image[0] - row_final = ip.flipud(self.MI.image)[-1] - self.assertEqual(np.array_equal(row_init, row_final), True) - - def test_center_of_mass(self): - """Test that we get correct x and y centroids""" - (x1, y1) = ip.center_of_mass(self.MI.image) - self.assertEqual((int(x1), int(y1)), (522, 669)) - (x2, y2) = ip.center_of_mass(self.MI.image, sigma=2) - self.assertEqual((int(x2), int(y2)), (540, 660)) - - def test_average_image(self): - """Test that we can average a number of images""" - images = [] - while len(images) < 10: - images.append(self.MI.image) - - ave_image = ip.average_image(images) - self.assertEqual(np.array_equal(ave_image, self.MI.image), True) - - def test_shape_image(self): - """Test that we can reshape our ndarray""" - self.assertEqual(self.MI.image.shape, (1024, 1392)) - image = ip.shape_image(self.MI.image, 16, 89088) - self.assertEqual(image.shape, (89088, 16)) - - def test_x_projection(self): - """Test we get expected value for x projection""" - x_proj = ip.x_projection(self.MI.image) - self.assertEqual(x_proj.sum(), 9792279) - self.assertEqual(int(x_proj.mean()), 7034) - self.assertEqual(int(x_proj.std()), 10005) - - def test_y_projection(self): - """Test that we get expected value for y projection""" - y_proj = ip.y_projection(self.MI.image) - self.assertEqual(y_proj.sum(), 9623943) - self.assertEqual(int(y_proj.mean()), 9398) - self.assertEqual(int(y_proj.std()), 10398) - - def test_gauss_func(self): - """Test we get correct value for a gaussian evaluation""" - ans = ip.gauss_func(1.0, 2.0, 3.0, 4.0) - self.assertEqual(round(ans, 2), 1.76) - - def test_gauss_fit(self): - """Test that we get the correct gaussian fit parameters""" - x_proj = ip.x_projection(self.MI.image) - y_proj = ip.y_projection(self.MI.image) - _, a_x, x0_x, sigma_x = ip.gauss_fit(x_proj) - _, a_y, y0_y, sigma_y = ip.gauss_fit(y_proj) - self.assertEqual(int(a_x), 30373) - self.assertEqual(int(x0_x), 660) - self.assertEqual(int(sigma_x), 124) - self.assertEqual(int(a_y), 29528) - self.assertEqual(int(y0_y), 541) - self.assertEqual(int(sigma_y), 127) +from lcls_tools.common.image_processing.image_processing import ImageProcessor +from lcls_tools.common.image_processing.roi import RectangularROI + + +class TestImageProcessing(unittest.TestCase): + data_location: str = "tests/datasets/images/numpy/" + + def __init__(self, methodName: str = "runTest") -> None: + super().__init__(methodName) + self.center = [400, 400] + self.size = (800, 800) + self.widths = (300, 300) + self.radius = 50 + self.image = np.load(self.data_location + "test_roi_image.npy") + + def test_process(self): + """ + Given an np.ndarray and roi process + and assert the return in an np.ndarray + """ + image_processor = ImageProcessor() + image = image_processor.process(self.image) + assert isinstance(image, np.ndarray) + roi = RectangularROI( + center=self.center, xwidth=self.widths[0], ywidth=self.widths[1] + ) + image_processor = ImageProcessor(roi=roi) + image = image_processor.process(self.image) + assert isinstance(image, np.ndarray) + # TODO:run coverage + + def test_subtract_background(self): + """ + Given an np.ndarray, check that when the image_processor + is passed a background_image. the subtract_background function + call subtracts the returns an np.ndarray + that is the difference between the two np.ndarrays + """ + background_image = np.ones(self.size) + image_processor = ImageProcessor(background_image=background_image) + image = image_processor.subtract_background(self.image) + assert image.all() == (self.image - 1).all() + + """ + Given an np.ndarray check that when the image_processor + is passed a threshold check that subtraction occurs correctly + """ + image_processor = ImageProcessor(threshold=1) + image = image_processor.subtract_background(self.image) + assert image.all() == (self.image - 1).all() + + def test_clip(self): + """ + Given an np.ndarray check that when the image_processor + is passed a threshold check that the np.ndarray elements + are clipped at to not drop below zero + """ + image_processor = ImageProcessor(threshold=100) + image = image_processor.subtract_background(self.image) + clipped_image = image_processor.clip_image(image) + assert clipped_image.all() >= np.zeros(self.size).all() diff --git a/tests/unit_tests/lcls_tools/common/image_processing/test_roi.py b/tests/unit_tests/lcls_tools/common/image_processing/test_roi.py new file mode 100644 index 00000000..9acc01e2 --- /dev/null +++ b/tests/unit_tests/lcls_tools/common/image_processing/test_roi.py @@ -0,0 +1,54 @@ +import unittest +import numpy as np +from lcls_tools.common.image_processing.roi import RectangularROI, CircularROI + + +class TestROI(unittest.TestCase): + data_location: str = "tests/datasets/images/numpy/" + + def __init__(self, methodName: str = "runTest") -> None: + super().__init__(methodName) + self.center = [400, 400] + self.size = (800, 800) + self.widths = (300, 300) + self.radius = 50 + + self.image = np.load(self.data_location + "test_roi_image.npy") + self.filled_value_test_image = np.load( + self.data_location + "test_roi_fill_value_image.npy" + ) + + def test_circular_roi_crop_image(self): + """ + Given a circular ROI and image, + test that image has correct size after cropping + (size of roi) + """ + circular = CircularROI(center=self.center, radius=self.radius) + cropped_image = circular.crop_image(self.image) + assert cropped_image.shape[0] == 2 * self.radius + assert cropped_image.shape[1] == 2 * self.radius + + def test_fill_value_outside_circle(self): + """ + Given a test image change the pixel values + of all elements of the image outside a radius of 10 to zero + Test that this image is now exactly equal to the test filled value image + """ + circular = CircularROI(center=self.center, radius=self.radius) + image = circular.fill_value_outside_circle( + img=self.image, center=self.center, radius=10, fill_value=0 + ) + assert image.all() == self.filled_value_test_image.all() + + def test_rectangular_roi_crop_image(self): + """ + Given a rectangular ROI and image, + test that image has correct size after cropping + (size of roi) + """ + rectangular = RectangularROI( + center=self.center, xwidth=self.widths[0], ywidth=self.widths[1] + ) + cropped_image = rectangular.crop_image(self.image) + assert cropped_image.shape == self.widths diff --git a/tests/unit_tests/lcls_tools/common/measurements/__init__.py b/tests/unit_tests/lcls_tools/common/measurements/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/unit_tests/lcls_tools/common/measurements/test_screen_beam_profile.py b/tests/unit_tests/lcls_tools/common/measurements/test_screen_beam_profile.py new file mode 100644 index 00000000..6820d7ab --- /dev/null +++ b/tests/unit_tests/lcls_tools/common/measurements/test_screen_beam_profile.py @@ -0,0 +1,120 @@ +from lcls_tools.common.measurements.screen_beam_profile_measurement import ( + ScreenBeamProfileMeasurement, +) +from lcls_tools.common.data_analysis.fit.projection import ProjectionFit +from lcls_tools.common.data_analysis.fit.methods import GaussianModel +from lcls_tools.common.image_processing.image_processing import ImageProcessor +from lcls_tools.common.image_processing.roi import RectangularROI +from lcls_tools.common.devices.device import Metadata +from lcls_tools.common.devices.screen import ( + ScreenControlInformation, + ScreenPVSet, + Screen, +) +import numpy as np +import unittest + + +class ScreenTest(Screen): + @property + def image(self) -> np.ndarray: + return self._image + + @image.setter + def image(self, image): + self._image = image + + +class TestScreenBeamProfileMeasurement(unittest.TestCase): + def __init__(self, methodName: str = "runTest") -> None: + super().__init__(methodName) + self.test_instantiate_pydantic_objects() + + def create_test_image(self, size: tuple, center: list, radius: int): + # make img that is a circle in the center of the image with known + # standard dev and mean. no imports, no calls to external or + # internal files. + image = np.zeros(size) + for y in range(image.shape[0]): + for x in range(image.shape[1]): + distance = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2) + if distance < radius: + image[y, x] = 1 + return image + + def test_instantiate_pydantic_objects(self): + # creating image processing class + self.radius = 50 + self.size = (800, 800) + self.center = [400, 400] + self.xwidth = 300 + self.ywidth = 300 + self.means = [150, 150] + self.sigmas = [30, 30] + self.amplitude = [99, 99] + self.offsets = [1, 1] + self.roi = RectangularROI(center=[400, 400], xwidth=300, ywidth=300) + self.image_processor = ImageProcessor(roi=self.roi) + + self.pvs = { + "image": "ArrayData", + "n_bits": "N_OF_BITS", + "n_col": "Image:ArraySize1_RBV", + "n_row": "Image:ArraySize0_RBV", + "resolution": ":RESOLUTION", + } + self.metadata = { + "area": "TEST", + "beam_path": ["SC_TEST"], + "sum_l_meters": 99.99, + } + self.control_name = "OTRS:TEST:650:" + + self.screen_pvs = ScreenPVSet(**self.pvs) + self.meta_data = Metadata(**self.metadata) + self.controls_information = ScreenControlInformation( + control_name=self.control_name, PVs=self.screen_pvs + ) + self.screen_test = ScreenTest( + controls_information=self.controls_information, metadata=self.metadata + ) + self.screen_test.image = self.create_test_image( + size=self.size, center=self.center, radius=self.radius + ) + self.gauss_model = GaussianModel() + self.projection = ProjectionFit(model=self.gauss_model) + + self.screen_beam_profile_measurement = ScreenBeamProfileMeasurement( + name=self.control_name, + device=self.screen_test, + image_processor=self.image_processor, + fitting_tool=self.projection, + ) + + def test_single_measure(self): + perform_single_measure = self.screen_beam_profile_measurement.single_measure() + self.assertIsInstance(perform_single_measure, dict) + assert len(perform_single_measure) == 2 + + def test_measure(self): + perform_measure = self.screen_beam_profile_measurement.measure() + self.assertIsInstance(perform_measure, list) + self.assertIsInstance(perform_measure[0], dict) + if isinstance(self.gauss_model, GaussianModel): + for key, val in perform_measure[0].items(): + if key == "amplitude_x": + self.assertTrue(90 <= val <= 100) + elif key == "amplitude_y": + self.assertTrue(90 <= val <= 100) + elif key == "mean_x": + self.assertTrue(140 <= val <= 160) + elif key == "mean_y": + self.assertTrue(140 <= val <= 160) + elif key == "sigma_x": + self.assertTrue(20 <= val <= 40) + elif key == "sigma_y": + self.assertTrue(20 <= val <= 40) + elif key == "offset_x": + self.assertTrue(0 <= val <= 1) + elif key == "offset_y": + self.assertTrue(0 <= val <= 1)