2_beads.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9ff51ed7",
   "metadata": {},
   "source": [
    "# Laser power increase - beads\n",
    "\n",
    "In order to download the data, run the [download_data notebook](0_download_data.ipynb) first. \n",
    "\n",
    "We immobilized TetraSpeck beads (ThermoFischer) onto a coverslip and placed it on a microscope. We cropped a small region around a single bead and acquired time-series for different values of a laser pulse duration. The time-series were acquired using a [Micro-Manager script](https://github.com/mufpga/analysis/blob/main/acquisition_scripts/beads_pulse_acquisition.bsh).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ffbea704",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, re\n",
    "import os.path as op\n",
    "import numpy as np\n",
    "import read_roi as rr\n",
    "import tifffile as tff\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# define some matplotlib style for all figures\n",
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['axes.spines.right'] = False\n",
    "mpl.rcParams['axes.spines.top'] = False\n",
    "mpl.rcParams['xtick.direction'] = 'in'\n",
    "mpl.rcParams['ytick.direction'] = 'in'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eafebd71",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_image(path, norm = False):\n",
    "    tif = tff.TiffFile(path)\n",
    "    \n",
    "    pages = []\n",
    "    for page in tif.pages:\n",
    "        array = page.asarray()\n",
    "        if norm:\n",
    "            array = array / np.max(array)\n",
    "            \n",
    "        pages.append(array)\n",
    "        \n",
    "    return np.array(pages)\n",
    "\n",
    "def read_roi(path):\n",
    "    return rr.read_roi_file(path)\n",
    "\n",
    "def extract_roi(roi):\n",
    "    y0 = roi['bead']['left']\n",
    "    x0 = roi['bead']['top']\n",
    "    w = roi['bead']['width']\n",
    "    h = roi['bead']['height']\n",
    "    \n",
    "    return x0, y0, w, h\n",
    "\n",
    "def average_img(img):\n",
    "    return np.mean(img, 0)\n",
    "\n",
    "def average_over_roi(img, roi):\n",
    "    if roi:\n",
    "        x0, y0, w, h = extract_roi(roi)\n",
    "    else:\n",
    "        x0 = 0\n",
    "        y0 = 0\n",
    "        w = img.shape[1]\n",
    "        h = img.shape[2]\n",
    "    \n",
    "    av_img = average_img(img)\n",
    "    \n",
    "    return np.mean(av_img[x0:x0+w, y0:y0+h])\n",
    "\n",
    "def get_roi_mask(img, roi):\n",
    "    mask = np.ones(img.shape)\n",
    "    \n",
    "    if roi:\n",
    "        x0, y0, w, h = extract_roi(roi)\n",
    "        mask[x0:x0+w, y0:y0+h] = 0\n",
    "    \n",
    "    return mask \n",
    "\n",
    "def overlay_roi_im(img, roi, value = 1):\n",
    "    arr = np.copy(img)\n",
    "    x0, y0, w, h = extract_roi(roi)\n",
    "    \n",
    "    for x in range(w+1):\n",
    "        arr[x0+x,y0] = value\n",
    "        arr[x0+x,y0+h] = value\n",
    "        \n",
    "    for y in range(h):\n",
    "        arr[x0,y0+y] = value\n",
    "        arr[x0+w,y0+y] = value\n",
    "    \n",
    "    return arr\n",
    "\n",
    "def extract_pulse(name):\n",
    "    rev = name[::-1]\n",
    "    ind = rev.find('_')\n",
    "    \n",
    "    return int(name[len(name)-ind:])\n",
    "\n",
    "def generate_filename(file):\n",
    "    rev = file[::-1]\n",
    "    ind = rev.find(os.sep)\n",
    "    \n",
    "    return file[len(file)-ind:]+'_MMStack.ome.tif'\n",
    "\n",
    "def run_all(root, roi):\n",
    "    dirs = [x[0] for x in os.walk(root)]\n",
    "    if root in dirs:\n",
    "        dirs.remove(root)\n",
    "        \n",
    "    pulses = []\n",
    "    diff_means = []\n",
    "    inners = []\n",
    "    outers = []\n",
    "    stds = []\n",
    "    for f in dirs:\n",
    "        pulse = extract_pulse(f)\n",
    "        file_path = op.join(f, generate_filename(f))\n",
    "    \n",
    "        img = read_image(file_path)\n",
    "        inner = np.mean(average_over_roi(img, roi))\n",
    "        \n",
    "        inners.append(inner)\n",
    "        pulses.append(pulse)\n",
    "    \n",
    "    return pulses, inners\n",
    "\n",
    "def plot_fig(t, v, reg, name, s_x, s_y, file_name = None):\n",
    "    ax = plt.subplot(111)\n",
    "    \n",
    "    ind_min = np.argmin(t)\n",
    "    ind_max = np.argmax(t)\n",
    "    yt0 = reg.intercept_+t[ind_min]*reg.coef_[0]\n",
    "    yt1 = reg.intercept_+t[ind_max]*reg.coef_[0]\n",
    "    plt.plot([t[ind_min], t[ind_max]], [yt0, yt1], '-', linewidth = 0.7)\n",
    "    \n",
    "    plt.plot(t, v, '.')\n",
    "\n",
    "    #plt.ylim([-0.1, 1])\n",
    "    #plt.xlim([-0.1, 1])\n",
    "\n",
    "    plt.title(name)\n",
    "    plt.xlabel(s_x)\n",
    "    plt.ylabel(s_y)\n",
    "\n",
    "    if file_name:\n",
    "        dir_path = 'figs'\n",
    "        if not os.path.isdir(dir_path): os.mkdir(dir_path)\n",
    "        plt.savefig(op.join(dir_path,file_name+'.eps'), format='eps')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a6328a00",
   "metadata": {},
   "outputs": [],
   "source": [
    "root = op.join('data','beads','640_5-500')\n",
    "roi = read_roi(op.join(root,'bead.roi'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5c1128a9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 432x288 with 1 Axes>,\n",
       " <matplotlib.image.AxesImage at 0x2096d29ac40>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD1CAYAAAB9TzjVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAANGUlEQVR4nO3dT4xd5XnH8d9vJoNp3NCAEMi1qYItKypd1GlTXImqBZFQhw3JolJYNEiNlCyC1ErZoGwSKaqUTZpuoiyiWrBoQEgtDQtKcaxKbtXWglSoOKIRlkuIsWUH0T+kJMaeebqYa3cYPD6P73nvuec+/n4ka+beOfOeZ87c35xz7338vo4IAahpad4FAJgdAg4URsCBwgg4UBgBBwoj4EBh75vl4Dd7R3xk+feuvFGstduhk3+vMvtsOVZL8/gZx2qsx76lxM94aPVJb/W1XhXbPmD7h7aP235k89ff0bk+wwPoaeqA216W9E1Jn5B0h6QHbd/RqjAA/fU5g98p6XhEnIiIdyQ9IemBNmUBaKFPwHdK+vGG2ycn9wEYiT4vsl3uif27GtvP65yOrj536fZO79GupT09dgngavQJ+ElJt224vUvSqY0brGib9i/f12MXAProc4n+vKS9tm+3fZ2kT0t6uk1ZAFqY+gweERdsPyzp7yQtSzoYET9474YDv1fZyjzqHuv70tfC+80tjagHoVejS0Q8I+mZRrUAaGzB/1QCuBICDhRGwIHCCDhQGAEHCiPgQGEEHCiMgAOFzXRGl6YynT9j7Zob6+wwYzXWLsK11dxYS8vtxvKWk7XkSun13QBGjYADhRFwoDACDhRGwIHCCDhQGAEHCiPgQGGzb3TpaiBo2dQQ0b2NlGsemEdzSstmnpZTAo11CqVMs0im6URKTrPUr+lkHkb6mwPQAgEHCiPgQGEEHCiMgAOFEXCgMAIOFEbAgcIIOFDY4kzZlJHtWhpaw+4zL+d+xlhLdPW17IobuqNPGm9nWeZYtOywu9Juen03gFEj4EBhBBwojIADhRFwoDACDhRGwIHCCDhQ2OwbXbreqG85NdJcxko0lGSbGhJTEEUkmzsy9c9jWqfMsWeNtmY4gwOF9TqD235V0luSViVdiIiPtigKQBstLtHviYg3GowDoDEu0YHC+gY8JD1n+/u2P9eiIADt9L1EvysiTtm+RdIh2/8eEUcufvG8zuno2vcubbzTu7XLu3vuEkBWr4BHxKnJx7O2n5J0p6RLAV/RNu1f+li/CgFMbepLdNvbbX/g4ueS7pN0rFVhAPrrcwa/VdJTXp9V432SvhMRzzapCkATUwc8Ik5I+vXeFcxj2qCWY7WcJmqsU04leCnXYRdrA79xc413xfE2GVAYAQcKI+BAYQQcKIyAA4URcKAwAg4URsCBwhZnbbKRNixkGjxS64TNYax5aPkzjvUxkZrGS8PUzhkcKIyAA4URcKAwAg4URsCBwgg4UBgBBwoj4EBhBBwobPadbF3TI2W7kVou8tewAyrddZXw7GsvNBvrWvD7u36z3WAtu+LcsIsw1RW3Nc7gQGEEHCiMgAOFEXCgMAIOFEbAgcIIOFAYAQcKm32jS6sGgpbNA4m1yZpOeZRZCy3p/l+7J7fLG3+pc5u1G96fG+v0G53bxFs/TY0Vq6ud2zz7H0dTYzXV8HeUMtB0U5zBgcIIOFAYAQcKI+BAYQQcKIyAA4URcKAwAg4URsCBwhZnyqaMptPurOR2mejMWrq+3WFe+9Avp7b7yW/d0LnNT2/L7XPHv2zv3Ob9R0/kBvvZz3PbJbTsNsz8HhdR5xnc9kHbZ20f23DfTbYP2X5l8vHG2ZYJYBqZS/RHJR3YdN8jkg5HxF5Jhye3AYxMZ8Aj4oikNzfd/YCkxyafPybpk23LAtDCtC+y3RoRpyVp8vGWdiUBaGWmL7Kd1zkdXX3u0u2d3qNdS3tmuUsAG0wb8DO2d0TEads7JJ293EYr2qb9y/dNXx2AXqa9RH9a0kOTzx+S9N025QBoKfM22eOS/lnSh22ftP1ZSV+T9HHbr0j6+OQ2gJHpvESPiAe3+NK9qT00m7IpcbGR3VfLsTLW2o21uj3XgPPfe7u38c6fpcZ6+/gvdG6zfSXZGNSw0SWzLtxYp95KP76y6+1t9e29vhvAqBFwoDACDhRGwIHCCDhQGAEHCiPgQGEEHCiMgAOFLc6UTQMt1na1vNzdaRSr7Wq/7tR/pba79fnrO7d5+0e5xQdveK27+yz+9+3UWC2ljv2F87nBorsrThr+Mdi3E48zOFAYAQcKI+BAYQQcKIyAA4URcKAwAg4URsCBwmbf6DLklE3zkPj5MlMLpXf35n+mtvvgke7pmD54/bbcTn9+rnOT9E/YcA2w1Hpi2cdNpp+k6dp3ubr6rpk20tQAaIGAA4URcKAwAg4URsCBwgg4UBgBBwoj4EBhBBwobPadbB0y0+5I/Tt63r3T7r9rLbvPWi6AF+8kpyBK1J+eZqnh4omRmhopO1iirmwn29BjDTQFGWdwoDACDhRGwIHCCDhQGAEHCiPgQGEEHCiMgAOFzb3RpWkDS0vZRoSBp5L62+P/NOj+Ri1z7Ofxe8w087hd89OVcAYHCusMuO2Dts/aPrbhvq/Yft32i5N/98+2TADTyJzBH5V04DL3fyMi9k3+PdO2LAAtdAY8Io5IenOAWgA01udFtodtf0bSC5K+GBHvmbD7vM7p6Nr3Lt3e6d3a5d09dgngakwb8G9J+qrW57v/qqSvS/qjzRutaJv2L31s+uoA9DLVq+gRcSYiViNiTdK3Jd3ZtiwALUwVcNs7Ntz8lKRjW20LYH46L9FtPy7pbkk32z4p6cuS7ra9T+uX6K9K+vzsSgQwLTedQmeTG3xT7Pe9HRU07Ohp2o00fAdU06mdElM2eSX3Ekycv9A9VrL21FRYA01n9C4tu+JaStR16MITWx58OtmAwgg4UBgBBwoj4EBhBBwojIADhRFwoDACDhQ2+ymbuhpZso02S4k1zEY6zVK2rlgbdgqiTANLVnott3k0iwwt8/hay05V1u94cQYHCiPgQGEEHCiMgAOFEXCgMAIOFEbAgcIIOFAYAQcKm/vig+kpm1p2QC1yN1W2Cy/zMybHykzHlO5ky2jZ3TgPI3p8cQYHCiPgQGEEHCiMgAOFEXCgMAIOFEbAgcIIOFDY7BtdupoW0o0uieaHlmPNo4li6AaJeUwlldFyvbpsXUM/vlqOdQWcwYHCCDhQGAEHCiPgQGEEHCiMgAOFEXCgMAIOFEbAgcLmP2VTVqbzJ9v1kxprDgsZtuyAymg5/VNLLetq+ZjIajlWT51H0vZttv/e9su2f2D7jyf332T7kO1XJh9vnH25AK5G5k/lBUlfjIhflfTbkr5g+w5Jj0g6HBF7JR2e3AYwIp0Bj4jTEfGvk8/fkvSypJ2SHpD02GSzxyR9ckY1ApjSVT2BtP0hSR+RdFTSrRFxWlr/IyDplubVAegl/SKb7V+U9FeS/iQi/seJFxLO65yOxuFLt3fqdu3y7mnqBDCFVMBtr2g93H8ZEX89ufuM7R0Rcdr2DklnN3/firZpv+9tVy2Aq5J5Fd2S/kLSyxHxZxu+9LSkhyafPyTpu+3LA9BH5gx+l6Q/lPSS7Rcn931J0tckPWn7s5Jek/QHM6kQwNQ6Ax4R/yhpqyfc3dffXc/Vs00Na6vd22SnWWrZIKGRTlWUkTmm0njryvy+IzlWRss107LNQz2nDqNVFSiMgAOFEXCgMAIOFEbAgcIIOFAYAQcKI+BAYQQcKGz+Uza1nA7oWlhobkTTAc1M+nglft8DLfL37rEaPqazXX1b4AwOFEbAgcIIOFAYAQcKI+BAYYMH/GScGHqXzSxy7dJi10/t0xk84K8v8C9qkWuXFrt+ap8Ol+hAYY6Wb/BvHtz+iaQfbbr7ZklvzGyns7XItUuLXT+1b+2NiDhwuS/MNOAA5otLdKAwAg4UNmjAbR+w/UPbx20v1Gqktl+1/ZLtF22/MO96rsT2QdtnbR/bcN/CLPe8Rf1fsf365Pi/aPv+eda4lbEttz1YwG0vS/qmpE9IukPSg5NliBfJPRGxLyI+Ou9COjwqafOLLou03POjem/9kvSNyfHfFxHPDFxT1qiW2x7yDH6npOMRcSIi3pH0hNaXIEZjEXFE0pub7l6Y5Z63qH8hjG257SEDvlPSjzfcPjm5b1GEpOdsf9/25+ZdzBQqLPf8sO1/m1zCj/YpxkVjWG57yIBf7n/eL9J7dHdFxG9o/SnGF2z/7rwLusZ8S9IeSfsknZb09blW02HzctvzqmPIgJ+UdNuG27sknRpw/71ExKnJx7OSntL6U45FcmayzLO2Wu55zCLiTESsRsSapG9rxMf/SsttT74+2PEfMuDPS9pr+3bb10n6tNaXIB4929ttf+Di55Luk3Tsyt81Ogu93PPFcEx8SiM9/mNbbnvQTrbJWxt/LmlZ0sGI+NPBdt6D7d1aP2tL6/PYfWfMtdt+XNLdWm+RPCPpy5L+RtKTkn5Fk+WeI2KUL2RtUf/dWr88D0mvSvr8xee0Y2L7dyT9g6SX9P9Lz35J68/DBz/+tKoChdHJBhRGwIHCCDhQGAEHCiPgQGEEHCiMgAOFEXCgsP8D+sqgI9jD0xoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# example showing a single frame with the roi and the outer mask\n",
    "img = read_image(op.join(root,'Beads_pulse_length_100','Beads_pulse_length_100_MMStack.ome.tif'))\n",
    "plt.figure(), plt.imshow(overlay_roi_im(average_img(img), roi, 1_000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "630e4c16",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear fit: y=0.7928435863994558*x+202.52089647495367, R2 = 0.9994466237475355\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAETCAYAAADQ97psAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAupElEQVR4nO3de5xVdb3/8dd775lBU1FJJBQMNajQFGtCUylvpZWX7HasVDRyPCfLOtlF9JxTWZ6fld3OqV+/7EqWGWkGkqKCImjKiNe4qKGAQxCMKAapM7P3/vz+WGtgMw7DBmbPzN77/Xw85jF7r73WXt81yvp81/f2UURgZma1KdPfBTAzs/7jIGBmVsMcBMzMapiDgJlZDXMQMDOrYQ4CZmY1zEHArMJIulDS93bw2AMkbZSUTd//QdIpvVpAqyjyPAEbKCQtBz4REbP6uywDlaQG4CngqIj4Wy9833jgRxHxlp0unFUkPwlYTVOikv4dnAE8vrUAIKlue74sIpqBwZIae6NwVnkq6X9+qyGSzpN0j6SrJT0vaZmkdxd9PkfS1yTdK2mDpNsl7VP0+VGS/ixpvaRHJR3X5dgrJd0LvAgc1M35l0v6gqTHJP1T0s8kDZN0a3q+WZL2LvF850takh73tKQLiz47TtJKSZdIWitptaTze/jTvBu4u+j4UZJC0iRJzwB3SspI+g9JK9Lv/JWkPbvsXxws5gDv7em/h1UvBwEbyI4EngD2Ab4J/EySij7/KHA+sC/QAHweQNL+wJ+ArwND0u03ShpadOw5QBOwB7BiK+f/APBOYAxwGnArcFlangxwcYnnWwucCgxOy/tdSW8uOs9rgD2B/YFJwA+LA0wXb0r/Jl29A3gjcDJwXvpzPEmA2x34wVa+D2AJcHgPn1sVcxCwgWxFRPwkIvLAFGA4MKzo819ExJMR8RIwFRiXbj8buCUibomIQkTcASwA3lN07C8jYlFE5CKiYyvn/9+IWJM2vcwD5kfEwxHRBtwEHFHK+SLiTxHxVCTuBm4HJhSdpwO4IiI6IuIWYCPw+q2UaS9gQzfbvxIR/0z/Fh8DvhMRT0fERmAycFYPTUUb0u+1GuQgYAPZ3ztfRMSL6cvdu/ucpFmn87PXAh9Km2bWS1oPHEsSRDq1lHD+NUWvX+rmfUnnk/RuSfdLei797D0kTxOd1kVEbivX0tXzJE8vXRVfz35s+XSzAqhjywBabA9g/VY+syq3XZ1IZhWiBbg2Ii7oYZ/eHBa31fNJGgTcCJwLTIuIDkl/BNR13xI9RtI81VXx9awiCUydDgByJEFsRDfHvhF4dAfLYxXOTwJWjX4NnCbpZElZSbukHbDd3QDLfb4GYBDQCuTSzu137cS5biFp/+/Jb4F/l3SgpN2B/wZ+1+Vpo9g7SPo7rAY5CFjViYgWkqGUl5HcfFuAL1Cm/997Ol9EbCDpQJ5K0pTzUWD6TpzuZuANkvbrYZ+fA9cCc4FlwMvAp7vbUdJbgX+mQ0WtBnmymFmFkdQEjI2Iz/bCd90I/CztkLYa5CBgZlbD3BxkZlbDHATMzGqYg4CZWQ1zEDAzq2EVN1nslFNOiZkzZ/Z3MczMKk23ExQr7kng2Wef7e8imJlVjYoLAmZm1nscBMzMapiDgJlZDXMQMDOrYWUNApL2knSDpMfT9HpvkzRE0h2S/pr+Lk7RN1nSUklPSDq5nGUzM7PyPwl8H5gZEW8gSV+3BLgUmB0Ro4HZ6XskjQXOAg4BTgH+r6RsmctnZlbTyhYEJA0G3g78DCAi2iNiPcmSu1PS3aYA70tfnwFcHxFtEbEMWAqML1f5zMysvE8CB5Gsrf4LSQ9L+qmk3YBhEbEaIP29b7r//myZIm9lum0Lra2tNDY2bvq55ppryngJZmb975l1L/LTeU+X5bvLOWO4Dngz8OmImC/p+6RNP1vR3Wy2V6xzPXToUBYsWNBLRTQzG0BammH5PBg1AUaOJ5cv8It7l3P/0+v4r9PGluWU5QwCK4GVETE/fX8DSRBYI2l4RKyWNBxYW7T/yKLjR5DkSjUzq34tzTDldMi3Q7aB1W/7MvMee4IDDzmBT5wwDC3+f5uCQ28qWxCIiL9LapH0+oh4AjgRWJz+TASuSn9PSw+ZDlwn6TvAfsBowCnvzKw2LJ+XBIDIk8+1se+8y/mQQPN/A/MFhRxkG2Di9F4NBOVeQO7TwG8kNQBPA+eT9ENMlTQJeAb4EEBELJI0lSRI5ICLIiJf5vKZmQ0MoyaQz9RDPpAyZChAFCBfSHeIJEgsn1c5QSAiHgEau/noxK3sfyVwZTnLZGY20Gx4uYNvPfQqho/6Hufvv5JdBg+FmZcmN/1MFih6Ehg1oVfPXXFLSZuZVZPZS9bwk3lPc/GJozn64EM3fzBs7OZOYtiiw7g3OQiYmfWDZze2ceWflnCEnuTaMc9Q35AF9tm8w8jxW97we/nm38lBwMysD0UEd8/+E2sem8Ulh41hxP1XwJJ2uOfqXu/0LYWDgJlZX2hp5oUldzLtyZf5yHP/l7rIoT8r6fyNQlk6fUvhIGBmVmb5FfOJKaexe6GDszMZMp03/shAJgOoLJ2+pXAQMDPrTcWzfoHWhbN5bNFCji90JMM+C2x54z/lKnhpXVk6fUvhIGBm1luKZv1GJku+EAyJPCdk61C2fvMwz36+8RdzEDAz6y1Fs34jXyAZ4R9QyMNbzoU9Rw6IG38xBwEzs17y0n5HU0cdGYJMtg4VT/I6/KMD6ubfyUHAzGxntTTz9AMz+cWq/fnAO69lXH5h2Sd59RYHATOznfCPJ+9ll9+eyajo4Iq6QWjkdBh5yeYdBujNv5MTzZuZ7YCIYNojf+PWGb+nnhwZCqhzrH8F8ZOAmdl2WrX+Ja79/e+ZUP844489DN3esCkPQH+M9d8ZDgJmZiUqFIJfz1/B8kfm8B/rvkQm3wErB9aQz+3lIGBm1pN08tfKPd/C5QtexYlv3Jf/OGQdmbs6IPLJE8BL62DCJdv+rgHIQcDMbGtamokppxO5NvZVPd87/uvsndsAr3p10vRToU1AxcoaBCQtBzYAeSAXEY2SvgJcALSmu10WEbek+08GJqX7XxwRt5WzfGZm3Upr/61/e4ohuTayFKiPDvaeMxkiBtys353RF08Cx0fEs122fTciri7eIGkscBZwCEmO4VmSxjjFpJmVXZf1fjpr/3uRJZOtg0IeqcuKnxXcBFRsIDUHnQFcHxFtwDJJS4HxwH39Wywzq2pF6/2QbWD1qDPZN639ZyQ4Il3uYddXb075WOFNQMXKHQQCuF1SAD+OiGvS7Z+SdC6wALgkIp4H9gfuLzp2ZbptC62trTQ2bk5b3NTURFNTU7nKb2bVrmi9n3yujZXrX+Q1dYM23+yLl3soTvlYwU1AxcodBI6JiFWS9gXukPQ48CPgayQB4mvAt4GPA+rm+Oi6YejQoSxYsKCMRTazWhKvPZZCph7ygbINvPWMi4CLur/Zd035WAXKGgQiYlX6e62km4DxETG383NJPwFmpG9XAiOLDh8BrCpn+cyshrU0s+Hxu/jxiv0YMfaHfHCf5WQPevvmm3yV3ey3pmxBQNJuQCYiNqSv3wVcIWl4RKxOdzsTWJi+ng5cJ+k7JB3Do4HmcpXPzGpXYcV8ClNO51WFdj5X10Dm5Jth5Af7u1j9opxPAsOAmyR1nue6iJgp6VpJ40iaepYDFwJExCJJU4HFQA64yCODzKxXtTTz3KLZPLJwIccV2pNMX/mOfsntO1Ao4hXN7gNaY2NjuE/AzEqWDv/M7bI3ceulZAodKFtHpnit/4nTayEIdNfvOqCGiJqZ9a50+Gfk2yGgjkAUBnSmr77mIGBmVavjqblkc21kKJAlgzIZCA3oTF99zUHAzKrSn596lpkLh/DlbAMUOlAVLfXQmxwEzKw6pG3/G4e/jSsf2536bIYvfOIcsq2HVt0Er97kIGBmlS9t+y/k2qinnnPecx1jx5+UfFaFE7x6k9NLmlnF2/jEXeTTtv8G5Ri7dgbM+3YSHKxHfhIws4oVEUxd0MLjT+7Lf3a2/Wey8PB1tTb8c4f5ScDMKtKKdf/kgl8tYMPLOS6/cCKZ826GEy6HI85OAkBn1q8KS/ze1/wkYGYVJZcv8LN7lvHA8uf4+lte4jXP3wSrJmxu+29phkd+W3VLPpeLg4CZVYxFq17ghj/+gTOHLKPpkAPRHydvvtl3NvuMHJ+89oigkjgImNmA93JHnv+Z/VfqVi3gv567FLV2wJNdMn0Vr//jEUElcxAwswFt/tPr+NMt0zj7NS2M2Xc9tHQk7f2RgUwGkJt9doKDgJkNSP94uYNvznycERsX8tX1k9G6DshkIVMHBaoq2Xt/chAwswHnjsVrmHvnDJpGruKAvZ5LlnuOfHLz98JvvcpBwMwGjOefuId7Zv2R7G77cMX676F17a+s/Xvht15V1iAgaTmwAcgDuYholDQE+B0wiiSpzIfTRPNImgxMSve/OCJuK2f5zGxgiAjumvUnjvnzxzmVHFpX1Onr2n9Z9cWTwPER8WzR+0uB2RFxlaRL0/dfkjQWOAs4hCS95CxJY5xdzKwKpYu9MWoCLbsdyhUzFvMJ7qOBHOqu09e1/7Lpj+agM4Dj0tdTgDnAl9Lt10dEG7BM0lJgPHBfP5TRzMqlKNFLXnXM3uuTfGPMLgzZZzQ807B53L87fftEuYNAALdLCuDHEXENMKwz0XxErJa0b7rv/sD9RceuTLeZWTVZPo/It6PIk4lg4vr/RfPDN/5+Uu4gcExErEpv9HdIeryHfbvLf/mKBMitra00NjZuet/U1ERTU9POl9TMyidt/mkfeTQ3rj2AD1BHvSCjLhO+XloHEy7p79LWlLIGgYhYlf5eK+kmkuadNZKGp08Bw4G16e4rgZFFh48AVnX9zqFDh+JE82YVZFPzTxtBPaPe8Uvqz78ZrbgHdn01zLzU6/z0o7IFAUm7AZmI2JC+fhdwBTAdmAhclf6elh4yHbhO0ndIOoZHA14M3KzCtS29m7pcG9l0rf+3ZZfAAZfAAUcmOwwb63V++lE5nwSGATdJ6jzPdRExU9IDwFRJk4BngA8BRMQiSVOBxUAOuMgjg8wqQNFIny1u4i3NPPXATGav6GBScZ7frrV9r/PTrxTximb3Aa2xsTHcHGQ2QKRNPV1X8nzhyXvY9bfvpy46UN0g5A7fgaC7flfPGDaznbB8XhIA0gQusWwe057dn7a7buDD5BDu8B3onFnMzHbcqAnJE4CyFLL1/PqxDQx77Id8YMJhSdOPsu7wHeD8JGBm26+oHyB/zjQemXczjz6X5bz1PyKzrgNaPOa/UpQcBCTtFhH/LGdhzGwA67zxFw3rLGTrufLV3+DAIy7gvPYbyNzVsTm3r5uAKsI2g4Cko4GfArsDB0g6HLgwIj5Z7sKZ2QBR3AEsEYUCokDkgs++bg2Dj3ottEyAuQ0e819hSnkS+C5wMsk4fiLiUUlvL2upzGxgKeoAjsiQR2TIkqlrYPAbjk/2cW7filRSc1BEtKTj/Tt5/L5ZLShqAopsA4VcG3nV03HSf7Nb4YVX3uw95r/ilBIEWtImoZDUAFwMLClvscys3xU1AeUz9fxs9yZOGlXPQY2n0OAbfdUoJQj8K/B9khU9VwK3AxeVs1BmNgAUrfZJPjhv3GAajv9Cf5fKetk2g0CaEOZjfVAWMxsgIoJ7O97AW6mjQZDNNpB93Tv6u1hWBqWMDvoF3SzpHBEfL0uJzKzvFY37Xz34TXxtxmJeP2wMR06cjlrudUdvFSulOWhG0etdgDPpZolnM6tQW2T6qufWPf+VK8bswj6vPxFGHgWjjurvEloZldIcdGPxe0m/BWaVrURm1re6ZPo6/4UfoOaAB7+/aUE4q147snbQaOCA3i6ImfWxlmbyd1/Nnc/k6aCOUJZMJoOisHnW7/J5/V1KK7NS+gQ2kPQJKP39d5LE8GZWCYrX+4dN4/4Lt34J8u28PdNA9r3fQC8950xfNaiU5qA9+qIgZlYGxcs9ZLKAiEKOQoAokCUgOuCl5zav8+NMXzVlq0FA0pt7OjAiHirlBJKywALgbxFxqqSvABcArekul0XELem+k4FJJDOSL46I20o5h5ltxRbr/RfSR/pAZMhkshDxyhq/Z/3WlJ6eBL7dw2cBnFDiOT5DMsN4cNG270bE1cU7SRoLnAUcQpJjeJakMU4xabYT0vX+I99OLjJIkKVAJuulni2x1SAQEcfv7JdLGgG8F7gS+Nw2dj8DuD4i2oBlkpYC44H7drYcZjWnqB/gvmN/zooHb+fwCafyxuGD3dRjWyhpATlJhwJjSeYJABARvyrh0O8BXwS69it8StK5JM1El0TE8yTLUtxftM/KdJuZbY+icf8d1LH80B/y/s98m0F12eRz3/ytyDaHiEr6MvC/6c/xwDeB00s47lRgbUQ82OWjHwEHA+OA1WxuduouCfIrZiq3trbS2Ni46eeaa67ZVlHMakph2TwKuTYUeerJ8ZF9n9kcAMy6KOVJ4IPA4cDDEXG+pGEkSWa25RjgdEnvIXmCGCzp1xFxducOkn7C5hnJK4GRRcePoJuZyUOHDmXBggUlnN6sxrQ08/yiO5n25MucnalHkUvy/HqYp/WglCDwUkQUJOUkDQbWAgdt66CImAxMBpB0HPD5iDhb0vCIWJ3udiawMH09HbhO0ndIOoZHA83bczFmNaWo3T9XKBBTTmdwoYNz6xrIvOeb7vS1kpQSBBZI2gv4CfAgsJGduzl/U9I4kqae5cCFABGxSNJUYDGQAy7yyCCzrSga/1/I1DNnlxM5odBBhgLkO5zf10pWymSxzlzC/0/STGBwRDy2PSeJiDnAnPT1OT3sdyXJSCIz605n7f+FlZvW+4l8cPTB+5BZPMgzfW27lbJsxDTgd8C0iFhe9hKZ2ZaKUjx2LumQV5Z8ZKhL1/p/1VvPhree7eGftt1KaQ76DvAvwP+R1EwSEGZExMtlLZmZbbnsg0QUCogCRBDjziXz6gO2vOn75m/bqZTmoLuBu9PlH04gWfLh52w5A9jMyqFo2YdCZCggssommb7e8jHf9G2nlTpZbFfgNJIngjcDU8pZKDNLjZpAZOsp5IJCpp44+f9Q177eTT7Wa0rpE/gdcCQwE/ghMCciCuUumFlNa2kmls3jzrYxzNv7Kv5t1CqGHfZO3/it15XyJPAL4KMermnWR1qaKUw5jci18/ZMPe+YOJ26157f36WyKlVKn8DMviiIWc1raSb/9FyeeGIJr8+1k6VANnLwzL3wWuf5tfIoqU/AzMqspZnCL0+DfDtjVEcmWw+FnMf8W9k5CJj1p5ZmOp6ay6Ili3hTPqn9Qx6OOBf2HOkOYCu7smcWM7OtaGkm/8vTyOTbeVOmjmxx7f/wj/rmb32iLzKLmVkXG17u4P5bbuCEztp/5OHNrv1b3ytrZjEzSxWt+PnA8udY2nwrbx37OrKtRev9uPZv/aCUeQKvIkkNeUBENEkaDbw+ImZs41Azgy0yfeUiwxGCRgroQef5tf5X6jyBB4Gj0/crgd+zORmMmXUnrf3H+hYi10aGAnUUUABE8gTgJZ+tn5USBA6OiH+R9BGAiHhJUnepIM2sU5fav1SHyKNMFpCHf9qAUUoQaE/XDgoASQcDbWUtlVmlSmv/hfUt0Fn7F+gtRZ2+4CWfbcAoJQh8hWTdoJGSfkOSO7jkOezp6qMLgL9FxKmShpAsRz2KJLPYhyPi+XTfycAkIA9cHBG3lXwlZv2tS+0/m6kjIp/k+e3a6eubvw0QpSwbcbukB4GjAAGfiYhnt+McnwGWsHnp6UuB2RFxlaRL0/dfkjQWOAs4hCTH8CxJY7xmkVWK3NNzyaS1/3qBPOTTKkBmWztIuhbIRcSf0hFBu0maXcqXSxoBvBf4adHmM9i8FPUU4H1F26+PiLaIWAYsBfwvxwa+lmb+dvPX+cVD/yCyDYSym2v/Ey5xALABrZTmoHuA+ZI+B+wPfAEodTjD94AvAnsUbRsWEasBImK1pH3T7fsD9xfttzLdZjbwpG3/bfV7kbltMq+JDj5R14De/Q0P+bSKUkpz0I8lLQLuAp4FjoiIv2/rOEmnAmsj4kFJx5VQlu5GHEXXDa2trTQ2Nm5639TURFNTUwlfb9ZLOtv+c21kEVkiSfmY7/CQT6s4pUwWOwf4T+Bc4DDgFknnR8Sj2zj0GOB0Se8BdgEGS/o1sEbS8PQpYDiwNt1/JTCy6PgRwKquXzp06FAWLFiwrWKb9b609v/ysyuoz7UlSz2TQZkMhDzk0yqSIl5R2d5yB+mPQFNErE3fjweuiYhxJZ8keRL4fDo66FvAuqKO4SER8UVJhwDXkfQD7AfMBkZ37RhubGwMBwHrcy3NRFr7z5GlLiMykU9u/J71a5Wh2/ldpTQHva/L++Y0EOyoq4CpkiYBzwAfSr93kaSpwGIgB1zkkUHWr4rW+/nHkrvYLa3910se+WNVo6elpL8YEd+U9D9b2eXiUk8SEXOAOenrdcCJW9nvSuDKUr/XrGyKxvznVccNe/4b59U1QL6j+3H/ZhWqpyeBJenvB/uiIGYDQmft/4WVRL4dRZ5MwHlHDCZz4M2e6WtVp6elpG9Of08BkDQ4eRsb+qhsZn2j88a/66th5qVJ7Z8MQYY6QSbbAAemN37f/K3KlDI6qJFkJdE9krdaD3w8IvyEYJUvbfYh3w4SUSggCmRI1vvRXm73t+pWymSxnwOfjIh5AJKOJQkKh5WzYGZ9Yvm8JABEnnwIyJBRNqn9j3O7v1W/UoLAhs4AABAR90hyk5BVtqImoHymHvIB2QaynvFrNaaUINAs6cfAb0lm8P4LMKczEb0TzlvF6Kbtv4M6Zo74DCcf1MCg173DN36rOaUEgXHp7y932X40TjhvlaKo7T8kSNv+6wWnj9nFSz1YzSplspgTzlvlK2r7L4RAWTIoGfPvpR6shpXyJGBWWYpm+gJJpq9dhlBQPYpA2QYybvs3AxwErNoUD/lM8/lGIUcHdTzw+i9y9H4ic6Bv/GadHASsOhTN9O1s9ol8AQARNAiO3V9u+zfroqQgIOlokpzAm/aPiF+VqUxm26dr7T9TRxSgIzJk09U+3fZv1r1SZgxfCxwMPEKSAB6SUUEOAta/uqv9F+DhfU5nFftw7Dvfx167Nni9H7MelPIk0AiMjW0lHjDrS93U/gsF6Ig68oedxakTTtm8r2/+ZltVShBYCLwGWF3mspiVrmjIZxRg/t6nsmHQazj6pPfx1oOP7u/SmVWMUoLAPsBiSc1AW+fGiDi9bKUy25riWb/ZBgq5djqijsFHnctR40/q79KZVZxSgsBXduSLJe0CzAUGpee5ISK+LOkrwAVAa7rrZRFxS3rMZGASSd/DxRFx246c26pU8azfbD2/HfJJRuzyMkcefzpjD3xbf5fOrCKVMmP47h387jbghIjYKKkeuEfSreln342Iq4t3ljQWOAs4hCTH8CxJY5xi0rpL9FLIBSe+tp5h7+26momZbY9SRgcdBfwv8EagAcgC/4yIwT0dl3Ykb0zf1qc/PXUunwFcHxFtwDJJS0mSzt+3rTJaFeqy2Bv5dgqZLPnIkAUydQ0MO8zNP2Y7q5TmoB+Q1NB/TzJS6FxgdClfLilLkp7ydcAPI2K+pHcDn5J0LrAAuCQingf2B+4vOnxlus1qTddEL1FAUSDywUuHfozBww70kE+zXpIpZaeIWApkIyIfEb8AjivxuHxEjANGAOMlHQr8iGTewTiSEUffTndXd1/RdUNrayuNjY2bfq655ppSimKVZIuRPwXyIQpkydQNYvCR5yazfh0AzHpFKU8CL0pqAB6R9E2SG/du23OSiFgvaQ5wSnFfgKSfADPStyuBkUWHjQBWdf2uoUOHsmDBgu05vVWKoiagSEf+5FXHyyd9ncGFDa79m5VBKUHgHJInhk8B/05yo/7Atg6SNBToSAPArsBJwDckDY+IzjkHZ5LMQwCYDlwn6TskHcOjgebtuRirYEVNQPlMPT/f/ULeMTLLmCPfTYNv/GZlU8rooBXpTXx4RHx1O757ODAl7RfIAFMjYoakayWNI2nqWQ5cmJ5nkaSpwGIgB1zkkUFVrnjJ5+XzNo38IR+cc/ge7HLCF/q7hGZVr5TRQacBV5OMDDowvYFfsa3JYhHxGHBEN9vP6eGYK4Ert1UmqwLFnb/ZBhYffhkHRx0Ngmy2gezod/R3Cc1qQqmTxcYDcwAi4hFJo8pXJKsJRZ2/+Vwba9es4uCJ09DKP7vt36wPlRIEchHxgtTd4B2z7ZQ2AcWuQ8gXZfo67uT3Jzd+z/w161MlLSAn6aNAVtJo4GLgz+UtllWVLhO/It9OB3X8efTnmbB/hsxBb3fN36yflBIEPg1cTrIMxG+B24CvlbNQVkWK1/uRiEKeDEG94LiRWWf6MutnpYwOepEkCFxe/uJY1Slq+y+EkLIE4UxfZgNEKaODGoHLeGV6ycPKVyyreGkTUMegvUB1ZCLI1DWgU74BL61z56/ZAFFKc9BvgC8AfwEK5S2OVYW0CSjybRSinieOuIw37Z1HB/rGbzbQlBIEWiNietlLYlWjbend1OXayFKgQTkOG5J327/ZAFVKEPiypJ8Cs9kys9gfylYqq0wtzTz1wExmr+hgUrYBCh1u+zcb4EoJAucDbyDJB9DZHBSAg4BtavvfmBlMw6zLOTA6uKBuEHq32/7NKkEpQeDwiHhT2UtilaN4tc+ZlxK5NnZBZAlEIRkN9NI6NwGZVYBSgsD9ksZGxOKyl8YGvqJx/wUJ0nH/IoMyGQiBm4DMKkYpQeBYYKKkZSR9AiLJHukhorWoaLXPCJHJZCFd+oFTrnITkFmFKSUInFL2UtjAlzYBrenYjb2oox7Scf++8ZtVspLyCfRFQWwAa2kmppxO5NrYW/X88/ivs7ec6cusGpTyJGC1Kq39t/7tKYak4/7rySUBwJ2+ZlWhbEFA0i7AXGBQep4bIuLLkoYAvyNZhmI58OGIeD49ZjIwCcgDF0fEbeUqn21DUe1/L7JksnVQyHvcv1mVKeeTQBtwQkRslFQP3CPpVuD9wOyIuErSpcClwJckjQXOAg4hyTE8S9IYp5jsY2ntf9WKvzIsrf1nJDjiXNhzpJuAzKpM2YJARASwMX1bn/4EcAZwXLp9CknGsi+l26+PiDZgmaSlJBnN7itXGa2LlmZiymkUcu3sS5ZMth4KuWTI5+Ef9c3frAqVtU8gTTL/IPA64IcRMV/SsIhYDRARqyXtm+6+P3B/0eEr021Wbi3NxLJ5rHj6CUbm2slSAAmOONu1f7MqV9YgkDbljJO0F3CTpEN72L27/JXRdUNrayuNjY2b3jc1NdHU1LSzRa1dLc0UppxG5NoZIdf+zWpNn4wOioj1kuaQzDlYI2l4+hQwHFib7rYSGFl02AhgVdfvGjp0KAsWLCh3kWtCoRA8OvdmDuus/ePav1mtyZTriyUNTZ8AkLQrcBLwODAdmJjuNhGYlr6eDpwlaZCkA4HRQHO5ylfrVv5lDjf+z+d4sW5PMnWDQNnNtf8JlzgAmNWIcj4JDAempP0CGWBqRMyQdB8wVdIk4BngQwARsUjSVGAxkAMu8sigXtbSTP7pudy9ssCxS6/mg+TQRi/3YFbLlAziqRyNjY3h5qAd0NJM4ZenEfl2UIYMBRSF5AnghMs9+cus+nXX7+oZw1WvpZmOp+by2KKFjMunbf8BZDKAV/w0q3UOAtWspZn8L08jk29nXKaObPHIHzcBmRkOAlXrhRc7uG/G73nnptp/Ht7sWb9mtiUHgWrQmekrbdZ54v5b+M2aA/jI+JPJ3n5tkunL4/7NrBsOApWuKNNXZLLk8sFo8ny1bhAaPh0mTt8cIBwAzKwLB4FKV5zpK1+gDhCR1P6Xz/OYfzPrUdkmi1nfWL1XI+3UUSCLsg0oW7954pdH/ZjZNvhJoELlVtzPgjk3c+fLY5j0wRsY9twDm2/6bv4xsxI5CFSgpx66kxE3n8X46ODIukFo8HQ4pGiyl2/+ZlYiB4EK0rbsPu6bPY26jX/jIHKIwua2f9/4zWwHOAhUiIXz72D0rR/j7XQkqR4zdVDAbf9mtlMcBAa4jUv/zL2z/she7WtoUC5Z76eQh7d44peZ7TwHgQHs/rm3csSdE3kXHahr7d8Tv8ysFzgIDEDPPTGPe++YxojMOtf+zaysHAQGkIjgzlkzOPbPkziVHMpkXfs3s7JyEBgg/r5wLvfOvok37PoCDeRQ5JObv2v/ZlZGZQsCkkYCvwJeQ3I7uyYivi/pK8AFQGu662URcUt6zGRgEpAHLo6I28pVvgEhzfR176rgyCe/xfvJoX+49m9mfaecTwI54JKIeEjSHsCDku5IP/tuRFxdvLOkscBZwCHAfsAsSWOqNsVkmumLfDvHFGf6cu3fzPpQ2YJARKwGVqevN0haAuzfwyFnANdHRBuwTNJSYDxwX7nK2F9e7sjTfMdNHLO1TF+u/ZtZH+mTPgFJo4AjgPnAMcCnJJ0LLCB5WnieJEDcX3TYSnoOGhVpcfMsHps3g0NHjyK7atDmtf6d6cvM+kHZg4Ck3YEbgc9GxD8k/Qj4Gkn992vAt4GP030S5Oi6obW1lcbGxk3vm5qaaGpqKkfRe09LM21L72bm0+28e+X3eSM59Jhv/GbW/8oaBCTVkwSA30TEHwAiYk3R5z8BZqRvVwIjiw4fAazq+p1Dhw5lwYIFZStzr0vz/Nbl2zktbfsn0jV/XlqXrPdvZtZPypZPQJKAnwFLIuI7RduHF+12JrAwfT0dOEvSIEkHAqOB5nKVry88u7GNGdOnorTtPxMFUMbr/ZvZgFHOJ4FjgHOAv0h6JN12GfARSeNImnqWAxcCRMQiSVOBxSQjiy6q1JFBEcFND/+NGx9ayRVHvYfMrb9227+ZDUiKeEWz+4DW2NgYA7k5qOW5F7lixmJOHvwMZ+79NNmD3p584EQvZta/uut3dRDoLflCcMst03h56d28Y9zr2feer2yu/U+c7pu/mfW3boOAl43oBU+u2cD1f7iRy1q/RDY60Fwlnb/hpC9mNrA5COyEtlyem6bdxK6r7uOS/dqpW9sBkYfIbDn5yx3AZjZAOQjsoAdXPM+0m2/iy89PJlPoQOu7rPnjDmAzqwAOAtvppafu497Zf+TxQYcx+Y3Pkr0nrf17zR8zq0AOAtvhoXtv45A7zuFEOjipbhAcclVS6+/sAPaaP2ZWYRwESvDcP9u58k9LeO8Ld27O9NU543fidA//NLOK5SDQg4hg+qOreOie27hk1Gr2O+iNsKao5t954/fN38wqlIPAVrQunsc9s/7I4CH78pX130UPesavmVUfB4EuCoVg5m3TObH5At5HDq3vMubfi76ZWRVxECiy8rE53Dd7GoftsWFznl+P+TezKuYgALTnCky7+SbOeOzf+CA5tNFj/s2sNtR8EHhywWwemnszRw55kfrO2r/H/JtZjajZIPBie47r/3AjZz/xaUaTQ//sUvv3mH8zqwE1GQTu+euz/OCuv3LlPk/QQM4zfs2sZtVUEFj/Yjv/fcsSxrQvYcroFgbtcQAs9oxfM6tdZQsCkkYCvwJeQ1LPviYivi9pCPA7YBRJZrEPR8Tz6TGTgUlAHrg4Im7rrfIse/gu5t85jU++aTSjmr8Gf/W4fzOzcj4J5IBLIuIhSXsAD0q6AzgPmB0RV0m6FLgU+JKkscBZwCHAfsAsSWN6JcVkSzOj/vQRRuXb0f2CQgHwuH8zs7Ilmo+I1RHxUPp6A7AE2B84A5iS7jYFeF/6+gzg+ohoi4hlwFKgd6rmy+ehfHs67r+QjPt3snczs77pE5A0CjgCmA8Mi4jVkAQKSfumu+0P3F902Mp0284bNWHL1T7dBGRmBvRBEJC0O3Aj8NmI+IfUbZpL6D7/5SsSILe2ttLY2LjpfVNTE01NTT0XYuR4r/ZpZtaNsgYBSfUkAeA3EfGHdPMaScPTp4DhwNp0+0pgZNHhI4BVXb9z6NCh7FCiea/2aWb2CmXrE1BS5f8ZsCQivlP00XRgYvp6IjCtaPtZkgZJOhAYDTSXq3xmZlbeJ4FjgHOAv0h6JN12GXAVMFXSJOAZ4EMAEbFI0lRgMcnIoot6ZWSQmZltlSJe0ew+oDU2NsYONQeZmdW2bjtky9YcNFBcc801/V2EPleL1wy1ed21eM1Qm9ddrmt2EKhCtXjNUJvXXYvXDLV53Q4CZmbW6yquT0BSK7BiOw7ZB3i2TMUZqGrxmqE2r7sWrxlq87p39pqfjYhTum6suCBgZma9x81BZmY1zEHAzKyGVW0QkHSKpCckLU2XrK4akn4uaa2khUXbhki6Q9Jf0997F302Of07PCHp5P4p9c6RNFLSXZKWSFok6TPp9mq/7l0kNUt6NL3ur6bbq/q6ASRlJT0saUb6vhauebmkv0h6RNKCdFt5rzsiqu4HyAJPAQcBDcCjwNj+LlcvXt/bgTcDC4u2fRO4NH19KfCN9PXY9PoHAQemf5dsf1/DDlzzcODN6es9gCfTa6v26xawe/q6nmQl3qOq/brTa/kccB0wI31fC9e8HNiny7ayXne1PgmMB5ZGxNMR0Q5cT5KvoCpExFzguS6b+z5PQx+KgZSfog9FYmP6tj79Car8uiWNAN4L/LRoc1Vfcw/Ket3VGgT2B1qK3vdeboKBa4s8DUBxnoaq+lv0lJ+CKrzutFnkEZIVd++IiFq47u8BXyRJTdup2q8ZkgB/u6QHJXWukV/W667WRPMl5SaoEVX1t+jt/BSVIJKFFMdJ2gu4SdKhPexe8dct6VRgbUQ8KOm4Ug7pZltFXXORYyJiVZps6w5Jj/ewb69cd7U+CZSUm6DKrEnzM7AjeRoqQU/5KdLPq/K6O0XEemAOcArVfd3HAKdLWk7SlHuCpF9T3dcMQESsSn+vBW4iad4p63VXaxB4ABgt6UBJDSQJ7Kf3c5nKrarzNNRqfgpJQ9MnACTtCpwEPE4VX3dETI6IERExiuTf7p0RcTZVfM0AknaTtEfna+BdwELKfd393Rtexl7295CMIHkKuLy/y9PL1/ZbYDXQQVIbmAS8GpgN/DX9PaRo/8vTv8MTwLv7u/w7eM3HkjzqPgY8kv68pwau+zDg4fS6FwL/lW6v6usuupbj2Dw6qKqvmWQ046Ppz6LO+1a5r9vLRpiZ1bBqbQ4yM7MSOAiYmdUwBwEzsxrmIGBmVsMcBMzMapiDgFkPJJ0n6Qe98D27SrpbUnYHjv2UpPN3tgxm3XEQMOsbHwf+EMkSENvr58DFvVweM8BBwGqIpFGSHpc0RdJjkm6Q9Kr0s+WS9klfN0qa083xH5K0MF3bf266LSvpW5IeSL/zwq2c/mOkMz0lHde5Rn76/geSzktfXyVpcfpdVwNExIvAcknVtDKmDRAOAlZrXg9cExGHAf8APrkdx/4XcHJEHA6cnm6bBLwQEW8F3gpckE7h3yRduuSgiFje05dLGgKcCRySlu/rRR8vACZsR1nNSuIgYLWmJSLuTV//mmQ5ilLdC/xS0gUkiYsgWd/l3HSp5/kkU/xHdzluH2B9Cd//D+Bl4KeS3g+8WPTZWmC/7SirWUmqdSlps63puk5K5/scmytFu3R7YMS/SjqSJNnJI5LGkSzn++mIuK2Hc77UzXcWLwNcn35/Lm3yOZFk4bRPAScUlemlHs5htkP8JGC15gBJb0tffwS4J329HHhL+voD3R0o6eCImB8R/wU8S7KM723Av6XLXCNpTLoC5CYR8TyQlVQcCA5NRww1AG9LP98d2DMibgE+C4wr2n8MyQJyZr3KQcBqzRJgoqTHgCHAj9LtXwW+L2kesLURPN9Kk4AvBOaSrPb4U2Ax8FC6/cd0/4R9O1s2PT0PzATuI1kZcjKwJzAjLdvdwL8X7X8MMGs7r9Vsm7yKqNWMNC3ljIjoKTNXuc59BPC5iDgnzZb1+Yg4dXuPLWMRrUb5ScCsD0TEw8BdOzJZjKRj+T97uUhmgJ8EzMxqmp8EzMxqmIOAmVkNcxAwM6thDgJmZjXMQcDMrIY5CJiZ1bD/D1GZaoOcMIjEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute values\n",
    "pulse_roi, inner_roi = run_all(root, roi)\n",
    "\n",
    "# linear fit\n",
    "x = np.array(pulse_roi).reshape(-1, 1)\n",
    "reg = LinearRegression().fit(x, inner_roi)\n",
    "print(f'Linear fit: y={reg.coef_[0]}*x+{reg.intercept_}, R2 = {reg.score(x, inner_roi)}')\n",
    "\n",
    "plot_fig(pulse_roi, inner_roi, reg, 'Inner mean (roi)', 'pulse (\\u03bcs)', 'mean pixel value', 'fig_beads_5-500')\n",
    "\n",
    "dir_path = 'figs'\n",
    "if not os.path.isdir(dir_path): os.mkdir(dir_path)\n",
    "plt.savefig(op.join(dir_path,'fig_beads_pulse.eps'), format='eps')"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}