notebook_demo/Sim_with_cascade.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib ipympl\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "from pylab import *\n",
    "import pandas as pd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a function to load the data\n",
    "def load_data(start,end,ch,name=\"OKSeq\",root=\"../data_hela/\"):\n",
    "    #Start and end are in kb\n",
    "    #return x in kb and signalvalue\n",
    "    \n",
    "    data = {\"OKSeq\":\"OKSeq_5kb.csv\",\n",
    "            \"MRT\":\"MRT_10kb.csv\",\n",
    "            \"Ori\":\"ori_position.tsv\"}\n",
    "    \n",
    "    if \"tsv\" in data[name]:\n",
    "        sep=\"\\t\"\n",
    "    else:\n",
    "        sep=\",\"\n",
    "    \n",
    "    read = pd.read_csv(root+data[name],sep=sep)\n",
    "    #print(read.head())\n",
    "    \n",
    "    read = read[read.chrom==\"chr%i\"%ch]\n",
    "    read = read[ (start<read.chromStart/1000) &  (read.chromStart/1000< end)]\n",
    "    return np.array(read.chromStart)/1000,np.array(read.signalValue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation de reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.4e-04 ,NDiff 49\n",
      "EXperimental\n",
      "Fs cte , dt 3.3 (min)\n",
      "Average introduced 0.0\n",
      "Maxiiiiiiiiiiiiiii 126\n"
     ]
    }
   ],
   "source": [
    "from repli1d.fast_sim_break import get_fast_MRT_RFDs\n",
    "\n",
    "fork_speed= 1.5  # in kb / minute\n",
    "ndiff=1.1  # Number of firing factors per Mb\n",
    "Nsim = 200\n",
    "binsize = 5 # in kb\n",
    "dori = 10 # Average distance between origins in kb\n",
    "noise_ratio = 0.1 # 10 percent of random activations \n",
    "\n",
    "startp=5000 # in kb\n",
    "endp = 50000 # in kb\n",
    "ch=1\n",
    "genomic,DNase=load_data(startp,endp,ch=1,name=\"Ori\")  #Generate random profile\n",
    "RandomActivation = np.ones_like(DNase) * noise_ratio * np.sum(DNase) / len(DNase)  # Random activation\n",
    "profile = DNase+RandomActivation\n",
    "\n",
    "\n",
    "#Rescaling parameters to size of simulated chromosome\n",
    "kon = 1e-5 * 2875000 / (len(DNase) * binsize ) #whole genome ( 2875000 == whole genome I am working with)\n",
    "ndiff = int(ndiff * len(DNase) * binsize / 1000)\n",
    "print(\"%.1e ,NDiff %i\"%(kon,ndiff))\n",
    "\n",
    "data = get_fast_MRT_RFDs(Nsim,\n",
    "                         profile,\n",
    "                         ndiff,\n",
    "                         kon=kon,\n",
    "                         fork_speed=fork_speed / binsize, #Fork speed is rescaled in bin/minute\n",
    "                         single_mol_exp=True,\n",
    "                         pulse_size=5,\n",
    "                         dori=dori,\n",
    "                         return_dict=True,n_jobs=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation avec cascade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EXperimental\n",
      "Fs cte , dt 3.3 (min)\n",
      "Average introduced 0.0\n",
      "Maxiiiiiiiiiiiiiii 158\n"
     ]
    }
   ],
   "source": [
    "cascade = {\"infl\":50,\n",
    "           \"amount\":100,\n",
    "           \"down\":1,\n",
    "          \"damount\":0.001}\n",
    "\n",
    "\n",
    "data_cascade = get_fast_MRT_RFDs(200,\n",
    "                         profile,\n",
    "                         ndiff,\n",
    "                         kon=kon,\n",
    "                         fork_speed=fork_speed / binsize, #Fork speed is rescaled in bin/minute\n",
    "                         single_mol_exp=True,\n",
    "                         pulse_size=5,\n",
    "                         dori=dori,\n",
    "                         cascade=cascade,\n",
    "                         return_dict=True,n_jobs=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the cascade function is the followwing:\n",
    "\n",
    "def comupute_enhanced(proba, actual_pos,cascade={}):\n",
    "\n",
    "    #assert()\n",
    "\n",
    "    infl = cascade[\"infl\"] #,50)\n",
    "    amount = cascade[\"amount\"] #,4)\n",
    "    down = cascade[\"down\"]#,10)\n",
    "    damount = cascade[\"damount\"]#,1 / 10000.)\n",
    "\n",
    "    cproba = np.zeros_like(proba) + 1\n",
    "    for position, direction, init_position in actual_pos:\n",
    "        # print(position,direction,cproba,max(position-infl),position)\n",
    "        if direction == \"L\":\n",
    "            cproba[position:min(position + infl, len(cproba) - 1)] *= amount\n",
    "        if direction == \"R\":\n",
    "            cproba[max(position - infl, 0):position] *= amount\n",
    "\n",
    "        if direction == \"L\":\n",
    "            cproba[position:min(position + down, len(cproba) - 1)] *= damount\n",
    "        if direction == \"R\":\n",
    "            cproba[max(position - down, 0):position] *= damount\n",
    "\n",
    "    return cproba\n",
    "\n",
    "# then the proba i use to draw the origin is the current proba * cproba\n",
    "\n",
    "# you could replace by one of your own .\n",
    "# by doing:\n",
    "\n",
    "#from repli1d import fast_sim_break\n",
    "\n",
    "#def custom_comupute_enhanced(proba, actual_pos,cascade={}):\n",
    "#     blabla\n",
    "#     return cproba\n",
    "\n",
    "#fast_sim_break.compute_enhanced = your_own\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single mol analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Without correlation\n",
      "\n",
      "Mean inter ori distance 87.9833836858006\n",
      "Mean real inter ori distance 233.8429903125571\n",
      "Nfork 0, percent 0.61\n",
      "Nfork 1, percent 0.18\n",
      "Nfork 2, percent 0.15\n",
      "Nfork 3, percent 0.02\n",
      "Nfork 4, percent 0.03\n",
      "Mean number of fork / fiber having forks 1.8764493950693035\n",
      "\n",
      "With correlation\n",
      "\n",
      "Mean inter ori distance 45.784747527329515\n",
      "Mean real inter ori distance 64.05875776664588\n",
      "Nfork 0, percent 0.70\n",
      "Nfork 1, percent 0.17\n",
      "Nfork 2, percent 0.06\n",
      "Nfork 3, percent 0.01\n",
      "Nfork 4, percent 0.01\n",
      "Mean number of fork / fiber having forks 2.8777227775245646\n"
     ]
    }
   ],
   "source": [
    "from repli1d.single_mol_analysis import compute_info, compute_real_inter_ori ,compute_mean_activated_ori\n",
    "\n",
    "def compute_corre(datav,fiber_size=200): #fiber size in kb\n",
    "\n",
    "    resolution = binsize\n",
    "    Deltas0 = compute_info(Sme=datav[\"single_mol_exp\"], size_single_mol=int(fiber_size/resolution), n_sample=500)\n",
    "\n",
    "    print(\"Mean inter ori distance\", np.mean(np.concatenate(Deltas0[\"dist_ori\"])*resolution))\n",
    "\n",
    "    Deltas = compute_real_inter_ori(datav[\"position_time_activated_oris\"])\n",
    "    print(\"Mean real inter ori distance\", np.mean(Deltas)*resolution)\n",
    "\n",
    "    nfork_per_fiber = np.array([len(di) for di in Deltas0[\"size_forks\"]])\n",
    "    # hist(nfork_per_fiber,range=[0,5],bins=6,normed=True)\n",
    "    for i in range(5):\n",
    "        print(\"Nfork %i, percent %.2f\" % (i, np.sum(nfork_per_fiber == i)/len(nfork_per_fiber)))\n",
    "\n",
    "    print(\"Mean number of fork / fiber having forks\",\n",
    "          np.mean(nfork_per_fiber[nfork_per_fiber != 0]))\n",
    "    return Deltas0,Deltas\n",
    "    \n",
    "print(\"Without correlation\\n\")\n",
    "Delta = compute_corre(data)\n",
    "\n",
    "print(\"\\nWith correlation\\n\")\n",
    "\n",
    "Delta_corre = compute_corre(data_cascade)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd9218b9e90>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAu4AAAHwCAYAAAD5Keq8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5xWVb3//dcH0SEQUETDhBotFcM0wdsMyx90H4PyRyaW+bU00xL17lhi347pcUxLTSkTT3aO3qnl3Vc9ePIcUiyPyo9TmCGlDxMli9HwwDFFfsgIqKz7j2tfNgzz62KumWvWzOv5eMxje+211t7r2q6HvmfN3mtHSglJkiRJvduAWndAkiRJUscM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgbWugO9RUQsA4YBjTXuiiRJkvq2emBtSmnPShoZ3P9m2Nve9rYR++2334had0SSJEl915IlS3jttdcqbmdw/5vG/fbbb8Rjjz1W635IkiSpD5swYQKLFy9urLSd97hLkiRJGTC4S5IkSRkwuEuSJEkZMLhLkiRJGTC4S5IkSRkwuEuSJEkZMLhLkiRJGXAdd0mS1Kds3ryZVatWsW7dOjZu3EhKqdZdUh8WEdTV1TF06FBGjBjBgAHdNy9ucJckSX3G5s2b+ctf/kJTU1Otu6J+IqXEhg0b2LBhA+vXr2fMmDHdFt4N7pIkqc9YtWoVTU1NDBw4kFGjRjFkyJBunQGVNm/ezPr161m5ciVNTU2sWrWKkSNHdsu5HMmSJKnPWLduHQCjRo1i6NChhnZ1uwEDBjB06FBGjRoF/G0Mdsu5uu3IkiRJPWzjxo0ADBkypMY9UX9THnPlMdgdDO6SJKnPKD+I6ky7elpEAHTrw9COakmSJKmLysG9OxncJUmSpAwY3CVJkqQMuBykJEnqN+q/fm+tu9ApjVd9vNZdUC9kcK+1huG17kHPalhT6x5IkiRlyeAuSZL6nd46o53LXwRUGwb33qKvz0T3t78sSJIkVZkPp0qSJPVBjY2NRASnn346jY2NnHzyyYwcOZJBgwZx8MEH8/Of/3yrNhs3buSqq67igAMOYPDgwQwbNowPf/jD3HXXXV3uz9NPP80ZZ5xBfX09dXV17Lbbbnz4wx/mxhtv3Krugw8+yOTJkxkxYgSDBg1in3324etf/zpr1mw90XnkkUcSEWzatIlvfvOb7LvvvtTV1XH66acDcOuttxIR3Hrrrdx///0ceeSRDB8+fKvlG59++mlOP/10xowZQ11dHW9/+9s55ZRTeOaZZ7r83avFGXdJkqQ+7LnnnuOQQw5hr7324rOf/SyrVq3izjvv5Pjjj+c///M/OeqoowDYtGkTH/3oR5k3bx5jx47l3HPPpampiVmzZvHpT3+a3//+93z729/epj7ce++9nHTSSWzcuJHJkyfzmc98htWrV/P444/zne98h2nTpr1V95//+Z+ZNm0aQ4YM4aSTTmK33XZj7ty5XH311cyePZtf/epX7LTTTlud48QTT+S3v/0tU6ZM4ROf+AS77bbbFuWzZs3i/vvvZ8qUKZx99tk0Nja+VXb//ffzyU9+ktdff51jjz2W97znPSxfvpx/+7d/49577+Xhhx9m/Pjx2/Tdq8ngLkmS1IfNnTuXhoYGLr300rf2nXLKKUyePJlrrrnmreA+Y8YM5s2bx5QpU/iP//gPBg4sxcRLL72UQw45hCuvvJJjjjmGiRMnVnT+l156iVNOOYU33niDhx56iCOOOGKL8uXLl7/1z8899xxf/vKX2XHHHXn00UcZO3bsW2XnnHMON954I1/72tf4l3/5l63O89xzz/Hkk08ycuTIVvtx3333cd999zF58uQt9r/yyit85jOfYfDgwcyfP5/3vve9b5X94Q9/4AMf+ABnnnkmixcvruh7dwdvlZEkSerD3vWud3HxxRdvse+jH/0o73znO3n00Uff2vejH/2IiOC73/3uW6EdYLfdduOSSy4B4Oabb674/Lfddhtr165l2rRpW4V2gNGjR7/1z7fffjubNm3ivPPO2yK0A3zrW99i6NCh/OQnP2Hjxo1bHefyyy9vM7QDHH/88VuFdoAf//jHrF69mssuu2yL0A4wbtw4zjrrLH73u9/x1FNPdfhdu5sz7pIkSX3Y+9//frbbbrut9o8ZM4aFCxcCsG7dOp599ln22GOPrQIzwKRJkwD43e9+V/H5H3nkEQCmTJnSYd3yrHb5fM3tvPPOHHTQQcyfP5+nn36aAw88cIvyQw45pN1jt1VevgaPP/44DQ0NW5UvXboUgCVLlmwV7HuawV2SJKkPa+1+cICBAweyefNmgLce+tx9991brVvev3r16orPX26zxx57dFi3K/0YNWpUu8duq/zll18G4Kabbmq3/auvvtpueU/wVhlJkqR+bvjw0rLNK1eubLV8xYoVW9SrRPkXhxdeeKFb+9FylZjOlpeP9fjjj5NSavPntNNO67D/3c3gLkmS1M8NHTqUd7/73bzwwgv88Y9/3Kr84YcfBtimlVUOPfRQAObMmdNh3YMOOggoPVDb0urVq/n973/PoEGD2G+//SruR0f9W7BgQdWO2V0M7pIkSeKMM84gpcSFF17Im2+++db+l156icsvv/ytOpU67bTTGDZsGDfeeCPz58/fqrz5qjKnnnoq22+/PTNnzuTZZ5/dot4ll1zC2rVrOfXUU6mrq6u4H235/Oc/z0477cRll122xcO6ZZs3b271F4la8B53SZIkMX36dObMmcO///u/c+CBB/Kxj32MpqYm/vVf/5UXX3yRr33ta3zoQx+q+LgjR47kpz/9KVOnTuWoo45iypQpHHDAAaxdu5YnnniCv/zlLyxbtgyA+vp6rrvuOs4991zGjx/Ppz71KXbddVfmzZvHwoULGTt2LFdffXVVv/cuu+zCrFmzOOGEEzj00EP5yEc+wrhx4xgwYADPP/88Cxcu5OWXX2bDhg1VPe+2MLhLkqR+p/7r99a6C73ODjvswAMPPMB3v/tdfvrTnzJz5kwGDhzIgQceyHXXXcdnPvOZbT72xz/+cRYtWsTVV1/Ngw8+yC9/+Ut23nlnxo4dyz/8wz9sUfecc87hPe95D9deey133303TU1NjBkzhgsvvJCLLrqozYdtu+IjH/kITzzxBNdeey2/+MUvWLBgATvssAPveMc7mDRpEieeeGLVz7ktIqVU6z70ChHx2Pjx48c/9thjPXvihuLhioatX+Hbp/SX7ylJqqklS5YAtHkPdC6BvfGqj9e6C9oGHY2/sgkTJrB48eLFKaUJlRzfGXdJktRvGIiVMx9OlSRJkjLgjLskSZI6bfXq1Vx33XWdqnv66adTX1/fvR3qRwzukiRJ6rTVq1dz2WWXdarukUceaXCvIoO7JEmSOq2+vh4XN6kN73GXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIy4JtTJUlS/9EwvNY96JyGNbXugXohZ9wlSZKkDDjjLkmS+p/eOqOdy18EVBPOuEuSJEkZqGpwj4gPR8TdEbEiIjYW219GxMdaqTsxIu6LiFUR0RQRT0TE+RGxXTvHPyYi5kbEmoh4NSJ+ExGnVfM7SJIk9RWPPvoon/70p9ljjz2oq6tj99135+ijj+auu+7aot5dd93F4YcfzvDhw3nb297G+973Pq688ko2bty41THr6+upr69n7dq1fPWrX6W+vp7tt9+ehoYGABoaGogI5s6dy09/+lM+8IEPsOOOO1JfX7/FcX7zm98wdepURo0axQ477MCYMWP40pe+xH//93931+XIXtVulYmIi4HLgZeAnwMrgJHAQcCRwH3N6h4P3A1sAO4EVgHHAt8DDgNOauX45wEzgZeB24FNwFTg1oh4X0pperW+iyRJUu5uuukmpk2bxnbbbcdxxx3H3nvvzYsvvsiiRYv4wQ9+wKc+9SkALrroIq688kpGjhzJKaecwo477sicOXO46KKL+MUvfsEDDzzA9ttvv8WxN23axKRJk1i1ahVHH300w4YNY88999yizowZM3jggQc49thjOeqoo1iz5m+3J91yyy2cddZZ1NXVcdxxxzFmzBj++Mc/cvPNNzN79mweeeQR3vnOd3b/RcpMVYJ7RJxEKbT/J/DJlNK6FuXbN/vnYcBNwJvAkSmlRcX+S4CHgKkRcXJK6Y5mbeqBaykF/INTSo3F/m8CvwUuiIi7U0oLq/F9JEmScvbUU09xzjnnMGzYMBYsWMC4ceO2KF++fDkACxcu5Morr2TMmDE8+uijjBo1CoArr7ySE044gZ///Odcc801XHTRRVu0X7FiBe9973uZN28eQ4YMabUPDz30EAsXLuSggw7aYv/SpUv50pe+RH19PfPmzWOPPfbYos3f/d3f8fd///f87Gc/6/J16Gu6fKtMRAwArgaagFNahnaAlNLrzT5OBXYF7iiH9qLOBuDi4uO0Foc4A6gDbiiH9qLNK8C3i49nd+2bSJIk9Q033ngjb7zxBpdccslWoR1g9OjRAPzoRz8C4OKLL34rtAMMHDiQGTNmMGDAAG6++eZWzzFjxow2QzvAF7/4xa1Ce7lvr7/+Ot///ve3CO0AkyZN4rjjjmP27NmsW7dVpOz3qjHjPhHYE5gFvBIRHwf2p3QbzKOtzIJPKrb3t3Ks+ZR+AZgYEXUppY2daDOnRR1JkqR+7ZFHHgFgypQp7dZbvHgxUArMLe2zzz6MHj2aZcuWsXr1anbaaae3ygYNGsQBBxzQ7rEPOeSQVvcvXFiKhvPmzeO3v/3tVuUvvvgib775JkuXLmXChAntnqO/qUZw/7+K7f8Ai4H3NS+MiPnA1JTSX4td+xbbpS0PlFJ6IyKWAeOAvYAlnWizIiLWA6MjYnBKqam9zkbEY20UjW2vnSRJUi5Wr14NsNWMdkvl+8533333Vst33313nn/+edasWbNFcN9tt92IiHaP3XwGv7mXX34ZgGuuuabd9q+++mq75f1RNVaV2a3Yng28Dfi/gaGUZt1/ARwO/Guz+uUFSttaQLW8f6dm+zrbxsVPJUlSv1cO2S+88EK79YYPL0WnlStXtlq+YsWKLeqVdRTa26tTPtaaNWtIKbX5c8QRR3R4jv6mGsG9vHxjUJpZfzCl9GpK6Q/ACcBy4IiI+GAnj1f+t5wq6EOn26SUJrT2AzxdwfkkSZJ6rUMPPRSAOXPmtFuvfA/63Llztyp79tlnWb58OXvuuecWs+3V6tuCBQuqdsz+ohrB/ZVi++eU0uPNC1JKr1GadQco3+jU0ez4sBb1KmmztsPeSpIk9XHTpk1j4MCBXH755Tz11FNblZdXlTnjjDMAuOKKK/jrX//6Vvmbb77J9OnT2bx5M1/4wheq2rfzzjuP7bffnq985SssXbrVXdBs2rTJUN+Gatzj/kyxXd1GeTnYv61Z/YOBfYAt7jePiIGUHnR9A/hzi3OMLNosbNFmd2AIsLyj+9slSZL6g/e+97384Ac/4Oyzz+aggw7i+OOPZ++99+bll19m0aJFDB06lIcffpiJEyfyta99je985zvsv//+TJ06lSFDhjBnzhyefPJJPvShD3HhhRdWtW9jx47lRz/6EWeccQbjxo1j8uTJ7LPPPrz++us8//zzLFiwgF133ZWnn/ZmiJaqEdznUwrae0fEDimlTS3K9y+2jcX2IeB/AZOB/9Oi7uHAYGB+sxVlym0OK9q0XKVmSrM6kiRJHWvo+4/FnXXWWey///5ce+21zJ07l3vuuYeRI0dywAEHcOaZZ75V7+qrr+aggw7ihhtu4Mc//jGvv/467373u7niiiu44IIL2GGHHaret1NPPZUDDzyQGTNm8PDDD/PLX/6SIUOG8I53vIOpU6fy6U9/uurn7AsipUpuJW/jIBG3Uwrj30opXdxs/99RulVmLVCfUlpdvIDpT5Rubzms2QuYBlEK3x8EPtPiBUx7UlphZj0wodkLmHam9AKmdwMTu/ICpoh4bPz48eMfe6ytRWe6Sfk/HA1tPXfbR/SX7ylJqqklS0oL0u23336tV8glsPv/yyx1OP4KEyZMYPHixYuL5yw7rSpvTgW+CnwA+EZEHA48CryL0sOpbwJnpZRWA6SU1kbEWZTWfZ8bEXdQeiPqcZSWfZwF3Nn84CmlZRFxIXA9sCgi7gQ2UXqZ02hghm9NlSRJHTIQK2NVCe4ppRcj4gOU3nx6AnAosA64F7gypfRIi/r3RMQRwDeAE4FBwLOUfgG4PrXyZ4CU0syIaASmA5+j9GDtU8DFKaXbqvE9JEmSpN6qWjPupJRWUQreX+1k/V8BH6vwHLOB2ZX3TpIkScpbNZaDlCRJktTNDO6SJElSBgzukiRJUgYM7pIkSVIXVWOJ9Y4Y3CVJUp8REQBs3ry5xj1Rf1MO7uUx2B0M7pIkqc+oq6sDYP369TXuifqb8pgrj8HuYHCXJEl9xtChQwFYuXIl69atY/PmzT1yC4P6p5QSmzdvZt26daxcuRL42xjsDlVbx12SJKnWRowYwfr162lqamL58uW17o76mcGDBzNixIhuO77BXZIk9RkDBgxgzJgxrFq1inXr1rFx40Zn3NWtIoK6ujqGDh3KiBEjGDCg+25oMbhLkqQ+ZcCAAYwcOZKRI0fWuitSVXmPuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUAYO7JEmSlAGDuyRJkpQBg7skSZKUgaoE94hojIjUxs/KNtpMjIj7ImJVRDRFxBMRcX5EbNfOeY6JiLkRsSYiXo2I30TEadX4DpIkSVJvNrCKx1oDXNfK/ldb7oiI44G7gQ3AncAq4Fjge8BhwEmttDkPmAm8DNwObAKmArdGxPtSStOr8zUkSZKk3qeawX11Sqmho0oRMQy4CXgTODKltKjYfwnwEDA1Ik5OKd3RrE09cC2lgH9wSqmx2P9N4LfABRFxd0ppYRW/jyRJktRr1OIe96nArsAd5dAOkFLaAFxcfJzWos0ZQB1wQzm0F21eAb5dfDy7uzosSZIk1Vo1Z9zrIuJU4J3AeuAJYH5K6c0W9SYV2/tbOcZ8oAmYGBF1KaWNnWgzp0WddkXEY20Uje1Me0mSJKkWqhncRwE/abFvWUR8PqU0r9m+fYvt0pYHSCm9ERHLgHHAXsCSTrRZERHrgdERMTil1NSVLyFJkiT1RtUK7rcAC4A/AOsohe7zgC8CcyLigymlx4u6w4vtmjaOVd6/U7N9nWkzpKjXbnBPKU1obX8xEz++vbaSJElSrVQluKeULmux60ng7Ih4FbgAaABO6OThonzYCrqwLW0kSZKkbHT3w6k/LLaHN9tXnjUfTuuGtahXSZu1FfVOkiRJykR3B/cXi+2QZvueKbb7tKwcEQOBPYE3gD93ss3uxfGXe3+7JEmS+qruDu4fLLbNQ/hDxXZyK/UPBwYDv262okxHbaa0qCNJkiT1OV0O7hExLiJGtLL/XcANxcfbmxXNAl4CTo6Ig5vVHwRcUXy8scXhbgE2AucVL2Mqt9kZuKj4+EMkSZKkPqoaD6eeBHw9Ih4GllFaVebdwMeBQcB9lN56CkBKaW1EnEUpwM+NiDsovRH1OErLPs4C7mx+gpTSsoi4ELgeWBQRdwKbKL3MaTQww7emSpIkqS+rRnB/mFLgPojSrTFDgNXAf1Fa1/0nKaUtVntJKd0TEUcA3wBOpBTwnwW+Clzfsn7RZmZENALTgc9R+mvBU8DFKaXbqvA9JEmSpF6ry8G9eLnSvA4rbt3uV8DHKmwzG5hd6bkkSZKk3HX3w6mSJEmSqsDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGXA4C5JkiRlwOAuSZIkZcDgLkmSJGWgW4J7RHw2IlLxc2YbdY6JiLkRsSYiXo2I30TEaR0c97SIeLSov6Zof0x3fAdJkiSpN6l6cI+IMcBM4NV26pwHzAb2B24HbgLeAdwaEde20eZa4FZg96L+7cD7gNnF8SRJkqQ+q6rBPSICuAV4GfhhG3XqgWuBVcDBKaVzU0pfAQ4A/gRcEBEfbNFmInBBUX5ASukrKaVzgQnFca4tjitJkiT1SdWecf8yMAn4PLC+jTpnAHXADSmlxvLOlNIrwLeLj2e3aFP+/K2iXrlNI/BPxfE+38W+S5IkSb1W1YJ7ROwHXAV8P6U0v52qk4rt/a2UzWlRpyttJEmSpD5jYDUOEhEDgZ8AzwMXdVB932K7tGVBSmlFRKwHRkfE4JRSU0QMAfYAXk0prWjleH8stvt0sq+PtVE0tjPtJUmSpFqoSnAH/hE4CPhQSum1DuoOL7Zr2ihfAwwp6jV1sj7ATp3rqiRJkpSfLgf3iDiE0iz7jJTSwq53iSi2qcJ2naqfUprQ6klLM/HjKzynJEmS1CO6dI97s1tklgKXdLJZeYZ8eBvlw4rt2k7W72hGXpIkScpeV2fcd+Rv95ZvKK0GuZWbIuImSg+tng88A4ws2m0xQx8Ru1O6TWZ5SqkJIKW0PiJeAPaIiN1buc9972K71T3z6oUa2vr9qw9p8HdISZJUfV0N7huB/7eNsvGU7nv/L0phvRzSHwIOAybTIrgDU5rVae4h4LNFm1s62UaSJEnqM7oU3IsHUc9srSwiGigF99tSSjc3K7oF+BpwXkTcUl7LPSJ25m8r0rR8edMPKQX3b0TEPeW13IuXLp1L6ReIloFevUl/mIXuD39NkCRJNVOtVWU6LaW0LCIuBK4HFkXEncAmYCowmlYeck0p/Toivgt8FXgiImYBOwCfBkYA/0/zlzlJkiRJfU2PB3eAlNLMiGgEpgOfo/SQ7FPAxSml29poc0FEPAGcB3wR2AwsBq5JKf28RzouSZIk1Ui3BfeUUgPQ0E75bGB2hce8DWg12EuSJEl9WZeWg5QkSZLUMwzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgbWugNSn9MwvNY96H4Na2rdA0mS+h1n3CVJkqQMOOMuVUt/mIXuD39NkCSpl3LGXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKQFWCe0RcHREPRsRfIuK1iFgVEb+LiEsjYpc22kyMiPuKuk0R8UREnB8R27VznmMiYm5ErImIVyPiNxFxWjW+gyRJktSbVWvG/SvAEOAB4PvA/we8ATQAT0TEmOaVI+J4YD5wOPAz4J+AHYDvAXe0doKIOA+YDewP3A7cBLwDuDUirq3S95AkSZJ6pYFVOs6wlNKGljsj4lvARcA/AOcU+4ZRCt1vAkemlBYV+y8BHgKmRsTJKaU7mh2nHrgWWAUcnFJqLPZ/E/gtcEFE3J1SWlil7yNJkiT1KlWZcW8ttBfuKrZ7N9s3FdgVuKMc2psd4+Li47QWxzkDqANuKIf2os0rwLeLj2dvU+clSZKkDHT3w6nHFtsnmu2bVGzvb6X+fKAJmBgRdZ1sM6dFHUmSJKnPqdatMgBExHRgR2A4cDDwIUqh/apm1fYttktbtk8pvRERy4BxwF7Akk60WRER64HRETE4pdTUQR8fa6NobHvtJEmSpFqqanAHpgNvb/b5fuD0lNJfm+0bXmzXtHGM8v6dKmwzpKjXbnCXJEmSclTV4J5SGgUQEW8HJlKaaf9dRByTUlrcycNE+XAVnLrTbVJKE1o9QGkmfnwF55QkSZJ6TLfc455S+p+U0s+Ao4FdgB83Ky7Pmg/fqmHJsBb1KmmztsKuSpIkSVno1odTU0rPAU8B4yJiZLH7mWK7T8v6ETEQ2JPSGvB/blbUXpvdKd0ms7yj+9slSZKkXHX3qjJQekkSlNZth9Ja7QCTW6l7ODAY+HVKaWOz/e21mdKijiRJktTndDm4R8TYiBjVyv4BxQuYdqMUxF8pimYBLwEnR8TBzeoPAq4oPt7Y4nC3ABuB84qXMZXb7EzpBU8AP+zqd5EkSZJ6q2o8nDoZuCYi5gN/Al6mtLLMEZSWdFwJnFWunFJaGxFnUQrwcyPiDkpvRD2O0rKPs4A7m58gpbQsIp069xMAABR8SURBVC4ErgcWRcSdwCZKL3MaDczwramSJEnqy6oR3P8T+BfgMOBASss4rqe05vpPgOtTSquaN0gp3RMRRwDfAE4EBgHPAl8t6m+1OkxKaWZENFJacvJzlP5a8BRwcUrptip8D0mSJKnX6nJwTyk9CZy7De1+BXyswjazgdmVnkuSJEnKXU88nCpJkiSpiwzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGDO6SJElSBgzukiRJUgYM7pIkSVIGuhzcI2KXiDgzIn4WEc9GxGsRsSYi/isivhARrZ4jIiZGxH0RsSoimiLiiYg4PyK2a+dcx0TE3OL4r0bEbyLitK5+B0mSJKm3G1iFY5wE3AisAB4GngfeDnwSuBmYEhEnpZRSuUFEHA/cDWwA7gRWAccC3wMOK465hYg4D5gJvAzcDmwCpgK3RsT7UkrTq/BdJEmSpF6pGsF9KXAccG9KaXN5Z0RcBDwKnEgpxN9d7B8G3AS8CRyZUlpU7L8EeAiYGhEnp5TuaHaseuBaSgH/4JRSY7H/m8BvgQsi4u6U0sIqfB9JkiSp1+nyrTIppYdSSrObh/Zi/0rgh8XHI5sVTQV2Be4oh/ai/gbg4uLjtBanOQOoA24oh/aizSvAt4uPZ3ftm0iSJEm9V3c/nPp6sX2j2b5Jxfb+VurPB5qAiRFR18k2c1rUkSRJkvqcatwq06qIGAh8rvjYPHDvW2yXtmyTUnojIpYB44C9gCWdaLMiItYDoyNicEqpqYN+PdZG0dj22kmSJEm11J0z7lcB+wP3pZR+0Wz/8GK7po125f07bUOb4W2US5IkSVnrlhn3iPgycAHwNPDZSpsX29RurW1sk1Ka0OoBSjPx4ys4pyRJktRjqj7jHhHnAt8HngKOSimtalGlo9nxYS3qVdJmbQVdlSRJkrJR1eAeEecDNwBPUgrtK1up9kyx3aeV9gOBPSk9zPrnTrbZHRgCLO/o/nZJkiQpV1UL7hHxvym9QOn3lEL7i21UfajYTm6l7HBgMPDrlNLGTraZ0qKOJEmS1OdUJbgXL0+6CngM+EhK6aV2qs8CXgJOjoiDmx1jEHBF8fHGFm1uATYC5xUvYyq32Rm4qPj4QyRJkqQ+qssPp0bEacA3Kb0JdQHw5YhoWa0xpXQrQEppbUScRSnAz42IOyi9EfU4Sss+zgLubN44pbQsIi4ErgcWRcSdwCZKL3MaDczwramSJEnqy6qxqsyexXY74Pw26swDbi1/SCndExFHAN8ATgQGAc8CXwWuTylttTpMSmlmRDQC0ymtDz+A0gOwF6eUbqvC95AkSZJ6rS4H95RSA9CwDe1+BXyswjazgdmVnkuSJEnKXXe+gEmSJElSlRjcJUmSpAwY3CVJkqQMGNwlSZKkDBjcJUmSpAwY3CVJkqQMGNwlSZKkDBjcJUmSpAwY3CVJkqQMGNwlSZKkDBjcJUmSpAwY3CVJkqQMDKx1BySp12kYXuse9JyGNbXugSSpk5xxlyRJkjLgjLukyvWXGem+PBvdX/4dSlIf4oy7JEmSlAFn3CV1Xl+egZYkqZdzxl2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScrAwFp3QJJUQw3Da92D7tewptY9kKSqcMZdkiRJyoAz7pLUH/WHWej+8NcESf2KM+6SJElSBqoS3CNiakTMjIgFEbE2IlJE3N5Bm4kRcV9ErIqIpoh4IiLOj4jt2mlzTETMjYg1EfFqRPwmIk6rxneQJEmSerNq3SpzMXAg8CqwHBjbXuWIOB64G9gA3AmsAo4FvgccBpzUSpvzgJnAy8DtwCZgKnBrRLwvpTS9St9FkiRJ6nWqdavMV4B9gGHAtPYqRsQw4CbgTeDIlNIXUkoXAu8HFgJTI+LkFm3qgWspBfyDU0rnppS+AhwA/Am4ICI+WKXvIkmSJPU6VQnuKaWHU0p/TCmlTlSfCuwK3JFSWtTsGBsozdzD1uH/DKAOuCGl1NiszSvAt4uPZ29j9yVJkqRerxYPp04qtve3UjYfaAImRkRdJ9vMaVFHkiRJ6nNqsRzkvsV2acuClNIbEbEMGAfsBSzpRJsVEbEeGB0Rg1NKTe2dPCIea6Oo3fvyJUmSpFqqxYx7eWHdthYRLu/faRvauGivJEmS+qTe+AKmKLaduV++4jYppQmtHqA0Ez++gnNKkiRJPaYWM+4dzY4Pa1GvkjZru9AvSZIkqdeqRXB/ptju07IgIgYCewJvAH/uZJvdgSHA8o7ub5ckSZJyVYvg/lCxndxK2eHAYODXKaWNnWwzpUUdSZIkqc+pRXCfBbwEnBwRB5d3RsQg4Iri440t2twCbATOK17GVG6zM3BR8fGH3dRfSZIkqeaq8nBqRHwC+ETxcVSx/WBE3Fr880sppekAKaW1EXEWpQA/NyLuoPRG1OMoLfs4C7iz+fFTSssi4kLgemBRRNwJbKL0MqfRwIyU0sJqfBdJkiSpN6rWqjLvB05rsW+v4gfgOWB6uSCldE9EHAF8AzgRGAQ8C3wVuL61N7CmlGZGRGNxnM9R+mvBU8DFKaXbqvQ9JEmSpF6pKsE9pdQANFTY5lfAxypsMxuYXUkbSZIkqS+oxT3ukiRJkipkcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyYHCXJEmSMmBwlyRJkjJgcJckSZIyMLDWHZAk5av+6/fWugttahxU2nZXHxuv+ni3HFeS2uKMuyRJkpQBZ9wlqUp68+xzd+uVs88NpU21+9af/z1Lqi1n3CVJkqQMOOMuqVv051nJXjn7LEnKnsFdkiS1yV/Cpd7D4C6pW/k/PkmSqsPgLvWA/jxjJalv6E+/hPvfbPVWPpwqSZIkZcAZd/W4/jyT0Z9mrCT1gIbh3X6K8ousystr1lTDmlr3QKopZ9wlSZKkDDjjrppx9llSznrDXw/Ls+H1G37a/eeq5X+ze+AvC1IOnHGXJEmSMuCMuyRJFehVfy1sKG16VZ8kdRuDuySpb/M2C0l9hLfKSJIkSRlwxl2S1De5dKCkPsYZd0mSJCkDBndJkiQpAwZ3SZIkKQMGd0mSJCkDBndJkiQpAwZ3SZIkKQMGd0mSJCkDruPeS9R//d5ad0GSpN6th96C2ziofL4eOd2WfP+A2uGMuyRJkpQBZ9x7icarPl7rLkiS1Dv18Cx0Lf4K3jjolJqde4t+mEd6NWfcJUmSpAxkNeMeEaOBbwKTgV2AFcA9wGUppVdq2TdJktQ31GTWuaE4dzHzXjMNPXEO7+PfVtkE94h4N/BrYDfg34GngUOAvwcmR8RhKaWXa9hFSZIkdaDWtwO1lNPtQdkEd+AHlEL7l1NKM8s7I+K7wFeAbwFn16hvkiRJ267Gs9A9EaZr/teEPiCL4B4RewFHA43AP7UovhT4IvDZiLggpbS+h7snSZKUtR6ZdW7owXN1Qm+b+e+MXB5OnVRsf5lS2ty8IKW0DvgVMBg4tKc7JkmSJPWELGbcgX2L7dI2yv9IaUZ+H+DB9g4UEY+1UXTgkiVLmDBhwrb1cFuteLW0nd3D55UkSepJvSzzrHihdHvShAf+scfPvWTJEoD6StvlEtzLr0pr6waw8v6dunCON1977bU1ixcvbuzCMbbdisU1OW03Gltsn65pL/o+r3P38xr3DK9z9/Ma9wyvc0e6nnmqeo0X/081jlKxemBtpY1yCe4diWKbOqqYUuodv+b1ceW/bHi9u5fXuft5jXuG17n7eY17hte5+/Xna5zLPe7lGfXhbZQPa1FPkiRJ6lNyCe7PFNt92ijfu9i2dQ+8JEmSlLVcgvvDxfboiNiizxExFDgMeA14pKc7JkmSJPWELIJ7SulPwC8p3ch/boviy4AhwI9dw12SJEl9VU4Pp54D/Bq4PiI+AiwBPgAcRekWmW/UsG+SJElSt4qUOlyIpdeIiDHAN4HJwC7ACuAe4LKU0qpa9k2SJEnqTlkFd0mSJKm/yuIed0mSJKm/M7hLkiRJGTC4S5IkSRkwuEuSJEkZMLhLkiRJGTC4S5IkSRkwuKtqIqIxIlIbPytr3b+cRMTUiJgZEQsiYm1xDW/voM3EiLgvIlZFRFNEPBER50fEdj3V79xUcp0jor6d8Z0i4o6e7n9vFxG7RMSZEfGziHg2Il6LiDUR8V8R8YWIaPX/QY7lylR6nR3L2yYiro6IByPiL8U1XhURv4uISyNilzbaOJYrVMl17o9jOac3pyoPa4DrWtn/ak93JHMXAwdSum7LgbHtVY6I44G7gQ3AncAq4Fjge8BhwEnd2dmMVXSdC49TevFbS09WsV99xUnAjZRelvcw8DzwduCTwM3AlIg4KTV7oYhjeZtUfJ0LjuXKfAVYDDwAvAgMAQ4FGoAvRsShKaW/lCs7lrdZRde50H/GckrJH3+q8gM0Ao217kdf+AGOAvYGAjgSSMDtbdQdRuk/bhuBg5vtHwT8umh7cq2/U2/8qfA61xflt9a637n8AJMoBZUBLfaPohQuE3Bis/2O5Z65zo7lbbvOg9rY/63iev6g2T7Hcs9c5343lr1VRuqFUkoPp5T+mIr/MnVgKrArcEdKaVGzY2ygNKMMMK0bupm9Cq+zKpRSeiilNDultLnF/pXAD4uPRzYrcixvg224ztoGxThszV3Fdu9m+xzL26jC69zveKuMqq0uIk4F3gmsB54A5qeU3qxtt/q0ScX2/lbK5gNNwMSIqEspbey5bvVZ74iILwG7AC8DC1NKT9S4Tzl6vdi+0WyfY7n6WrvOZY7l6ji22Da/do7l6mvtOpf1m7FscFe1jQJ+0mLfsoj4fEppXi061A/sW2yXtixIKb0REcuAccBewJKe7Fgf9XfFz1siYi5wWkrp+Zr0KDMRMRD4XPGxebBxLFdRO9e5zLG8DSJiOrAjMBw4GPgQpTB5VbNqjuUu6uR1Lus3Y9lbZVRNtwAfoRTehwDvA/6Z0j1ocyLiwNp1rU8bXmzXtFFe3r9TD/SlL2sCLgcmADsXP0dQehjwSODBiBhSs97l5Spgf+C+lNIvmu13LFdXW9fZsdw104FLgfMphcn7gaNTSn9tVsex3HWduc79biwb3FU1KaXLinst/yel1JRSejKldDbwXeBtlJ4IV8+LYut93F2QUnoxpfSPKaXFKaXVxc984GjgN8B7gDNr28veLyK+DFwAPA18ttLmxdax3IH2rrNjuWtSSqNSSkFpkuqTlGbNfxcR4ys4jGO5A525zv1xLBvc1RPKD0cdXtNe9F3lmZvhbZQPa1FPVZRSeoPSknvgGG9XRJwLfB94CjgqpbSqRRXHchV04jq3yrFcmWKS6meUQuIuwI+bFTuWq6SD69xWmz47lg3u6gkvFts+9eeqXuSZYrtPy4LiHtc9KT2Y9uee7FQ/U/7TrWO8DRFxPnADpXWVjypWPGnJsdxFnbzO7XEsVyil9BylX5LGRcTIYrdjucrauM7t6ZNj2eCunvDBYut/oLrHQ8V2citlhwODgV+7ckG3OrTYOsZbERH/m9JLZ35PKUy+2EZVx3IXVHCd2+NY3jbvKLblFdQcy92j5XVuT58cywZ3VUVEjIuIEa3sfxel2R+AVl8lry6bBbwEnBwRB5d3RsQg4Iri44216FhfEhEfiIgdWtk/idKb/sAxvpWIuITSQ5KPAR9JKb3UTnXH8jaq5Do7lisXEWMjYlQr+wdExLeA3SgF8VeKIsfyNqj0OvfHsRy+d0TVEBENwNcpPcm9DFgHvBv4OKU3xd0HnJBS2lSrPuYkIj4BfKL4OAr4KKVZgwXFvpdSStNb1J9F6dXad1B6tfZxlJYkmwV8ypcMba2S61wsLTYOmAssL8oP4G/rNV+SUir/D1lARJwG3Eppdmwmrd/P25hSurVZG8dyhSq9zo7lyhW3IF1DaQ32P1FaK/ztlFYw2QtYSekXpqeatXEsV6jS69wfx7LBXVUREUcAZwMH8bflIFdT+pPtT4Cf+B+ozit+Ebq0nSrPpZTqW7Q5DPgGpVuTBgHPAj8CrvcFWK2r5DpHxBeAEygtrzcS2B74H2AhcENKaUFbB+mvOnF9AeallI5s0c6xXIFKr7NjuXIRsT+lN50eBoymtIzjekrrtN9LaWxu9RCwY7kylV7n/jiWDe6SJElSBrzHXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKgMFdkiRJyoDBXZIkScqAwV2SJEnKwP8Px+jggVdrSYIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist(np.concatenate(Delta[0][\"dist_ori\"]),histtype=\"step\",label=\"no_corre\");\n",
    "hist(np.concatenate(Delta_corre[0][\"dist_ori\"]),histtype=\"step\",label=\"corre\");\n",
    "legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}