{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from tqdm.auto import tqdm\n", "import noctiluca as nl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MSD analysis\n", "The mean squared displacement—MSD—is one of the most frequently used statistics on SPT data. ``noctiluca`` provides some functionality to calculate these. Note, however, that MSD analysis is statistically non-trivial (e.g. [Vestergaard et al.](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.89.022726) note that estimating diffusion coefficients from linear MSD fits is suboptimal). Thus, for production MSD analysis, we recommend [bayesmsd](https://github.com/OpenTrajectoryAnalysis/bayesmsd), a downstream analysis library building on ``noctiluca``.\n", "\n", "In the language of ``bayesmsd``, the MSDs calculated by ``noctiluca`` are *empirical MSDs*, i.e. directly calculated from the raw data, using the formula\n", "\\begin{equation}\n", "\\text{MSD}(\\Delta t) := \\left\\langle (x(t+\\Delta t) - x(t))^2 \\right\\rangle_t\\,.\n", "\\end{equation}\n", "Let's see how this is done." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Step 1: generate an example data set that we will use for the remainder of this example\n", "# We use 100 standard random walk (i.e. simple diffusion) trajectories in 2D, with 100 data points per trajectory\n", "data = nl.make_TaggedSet(np.cumsum(np.random.normal(size=(100, 100, 2)), axis=1))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "50431a3321fd417c95f46486871824b0", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/100 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Step 2: look at MSDs\n", "nl.plot.msd_overview(tqdm(data), color='tab:blue', alpha=0.2)\n", "\n", "plt.xlabel('time [frames]')\n", "plt.ylabel('MSD [a.u.]')\n", "plt.ylim([1, None])\n", "plt.title('Single trajectory & ensemble MSDs')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's dissect what just happened:\n", "\n", "+ in the first step we created an example data set. For more details on this, check the previous tutorials on [TaggedSet](02_TaggedSet.ipynb) and [I/O](03_IO.ipynb).\n", "+ [msd_overview()](../noctiluca.rst#noctiluca.plot.msd_overview) is a function that generates plots like the one above. It is generally useful to get a quick overview over your dataset.\n", "+ the actual MSD calculation is done by the ``nl.analysis`` module, which we will describe in more detail below; ``msd_overview()`` internally simply calls ``analysis.MSD()``.\n", "+ ``msd_overview()`` allows some customization of the plot by forwarding all keyword arguments to the ``pyplot.plot()`` calls for the single trajectory MSDs. For more detailed customization, use the returned list of ``matplotlib.Line2D`` objects\n", "+ note that we can use the nice progress bar provided by [tqdm](https://github.com/tqdm/tqdm) to keep tabs on the calculation while it is running; this comes in handy when working with large datasets\n", "\n", "## More on [analysis.MSD()](../noctiluca.analysis.rst#noctiluca.analysis.p2.MSD)\n", "As mentioned, the actual MSD calculation is performed by the ``nl.analysis`` module, which provides the ``MSD()`` function that can be used to calculate MSDs for single trajectories, or whole datasets:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvy0lEQVR4nO3dd3hUZfbA8e+ZNAgJJQQIzQwlIQkoYFyKZUV01zKKK0XBioggKCtFZUTEujpihwUbIrZFUFfUzK7+RFdFRV0QCyhiG0SaSAsdkry/P94bNyIJ6Xdmcj7Pk0fyzs2dE7mcvHnvuecVYwxKKaWii8ftAJRSSlU/Te5KKRWFNLkrpVQU0uSulFJRSJO7UkpFIU3uSikVhWLdDgBARDYBq92OQymlIky6MabZoV4Ii+QOrDbGHON2EEopFUlEZElpr+myjFJKRSFN7kopFYU0uSulVBQKlzX331m6dGnz2NjYWUAX9IdQNCgClhcUFAzPzc392e1glIp2YZvcY2NjZ6WlpWU3a9Zsq8fj0e5mEa6oqEg2bdqUs2HDhllAP7fjUSrahfOMuEuzZs3yNbFHB4/HY5o1a7Yd+5uYUnWeWCNEJKEmzh/Oyd2jiT26OH+f4XzNKVWbrgUeAV4VEanuk+s/tAo677zz0pcuXVqvMl/79ddfx2dkZHQu7/G//PJLTCAQOOQDCodz4okndvzll19iKvp1eXl5yW+88UaDyrynUqp8RGQoEHA+nWFqYGMNTe4VNG/evNW5ubl7a+O9Nm/eHPP44483P9RrBQUFZX7tO++8821qamphRd/zrbfeSl60aFFSRb7mwIEDFX0bpeokZylmCvAEIIDfGPNyTbyXq8ldRM4SkUeBRm7GcSj5+fmePn36dOzUqVNORkZG58cee6wJQI8ePTq9++67iQCJiYndx4wZ07pTp045Xbt2zVqzZk0swIoVKxK6du2a1aVLl+yxY8e2SkxM7H7w+QsKChg5cmSbLl26ZGdmZubcfffdqQcfM2HChDZr1qxJyMrKyhk5cmSbvLy85J49e2aeddZZ7Tp16tQZ4JRTTunQuXPn7I4dO3a+5557fj1H69atj1y/fn0swMyZM1OOPPLI7KysrJzzzz8/vfgHwwsvvNAwJycnu1OnTjm9e/fO/Prrr+OfeuqpZg8//HCLrKysnNdeey1p1apV8b17987MzMzM6d27d+Y333wTDzBgwADv8OHD2/Ts2TPziiuuaJuent5l3bp1sQCFhYUcccQRXYrfXykFItIOeAa4BVs9NsYYc1dNvZ+r//iMMa9i15tKfYQW4NoXPmu7asOOxOp878y05N13D+y6prTX//nPfzZMS0s78Pbbb38LdhZ98DF79uzx9O7de+f06dPXXnHFFW2mT5/ebOrUqeuvuuqqtqNHj/555MiRW6ZOnXrIZZUHHnggtVGjRoXLly//as+ePfKHP/wh66yzzsrPysraX3zMvffe+9OZZ55Zf+XKlV+CXTL5/PPPGyxbtmxF8XHPPvtsqEWLFoU7d+6U7t2751x44YVb09LSfp2xf/LJJ/VeeOGFlCVLlqxMSEgwF1544REPP/xw0/79+2+/6qqrvG+//fbKrKys/Rs3boxp0aJF4cUXX7wpKSmp8NZbb90I0Ldv347nn3/+5jFjxmx+4IEHmo4aNartwoULvwP47rvv6r3//vurYmNjmTBhQsGsWbNSpkyZ8vPLL7/cMDs7e0/Lli3L/vVCqSgmIh7gOOBs56Oj89I+YIgx5qWafH9dlinF0UcfvWfRokUNR40a1fq1115Latq06e+WOOLi4szgwYO3A+Tm5u5avXp1PMCyZcuShg0btgVg+PDhmw91/oULFzacP39+06ysrJzu3btnb926NfbLL7887Fr+UUcdtavkD4C77rqrRadOnXJyc3OzN2zYELdixYrfnOO1115LXr58eWLXrl2zs7Kyct57772G33//fcLbb7/doEePHjuKz9WiRYtDLuEsW7aswYgRI7YAjBo1asvSpUt/XbLp37//1thYOz8YNWrUL88991xTgNmzZ6cOHTr0l8N9L0pFIxGJcdbUVwLvAhOwiX0nduZ+ZE0ndgjjOveSypph15Sjjjpq3yeffPLliy++2OiGG25ovXDhwvx77rlnfcljYmNjjcfjKf4zBQUF5b7jbYyRe++998cBAwbkVySuxMTEouI/5+XlJb/zzjvJS5YsWZmcnFzUo0ePTnv27PnND2xjjAwaNGjzjBkz1pYcf/bZZxtV9QZ9UlLSr7F07NjxQGpqasErr7ySvGzZsgYLFiz4vkonVyrCOCWN5wJ+IMcZ/gl4HngR+MgYU2u/zerMvRShUCguOTm5aPTo0VvGjh278dNPPy33slC3bt12zpkzpwnA7NmzUw51zJ/+9KftDz30ULN9+/YJwOeff56Qn5//m7+PRo0aFe7atavUv6Nt27bFNGrUqDA5Oblo2bJl9T777LPfVbmcdtpp+Xl5eU3Wrl0bC7Bx48aYVatWxZ900km7Pvroo+SVK1fGF48DJCcnF+7YsePXJaju3bvvmjVrVhOARx55JOWYY47ZWVo8w4YN2zR8+PB2/fr121I8o1cq2olIBxG5G1gHPIVN7KuBi4F2xpjxxpj3azOxgyb3Ui1durR+t27dsrOysnLuuuuullOmTFl/+K+ypk+fvmb69OktjjzyyOz169fHJSUl/W7JY9y4cb9kZWXtPfLII7MzMjI6X3755ekHDhz4zVQ6LS2tMDc3d2dGRkbnkSNHtjn4HAMGDNheUFAgmZmZOZMmTWrVtWvXXSVfFxGTm5u7d/LkyWtPPvnkzMzMzJy+fftmrlmzJq5Vq1YF06ZNC51zzjkdO3XqlHPOOee0d865LRgMNi6+ofrQQw/9+PTTT6dmZmbmzJ07t+nMmTNL/S1qyJAh23fv3h0zYsSIQy5FKRUtRCRBRM4TkbeAb4FrgBTgM+AyoJMx5unaTui/ibEGyisrHoTIkoP7uX/22Wehrl27RuS67Y4dOzwNGjQo8ng8PProo03mzZuX8uabb35XW+9fUFBAampqt40bN36WkJBQa3/B7777buK4cePaLl269OvSjvnss89Su3bt6q2tmJSqTiLSGbgEGAoUF0vsAeYDM4H/1kTNehnx/C53FtPfnWvA+++/n3j11VcfYYyhYcOGhXPmzAnV5vtnZGR0HjJkyKbaTOyTJk1KmzNnTrMnnnjih9p6T6Vqg4jEAucBVwN/KPHS58CjwDPGmO1uxFYWnbmrWqUzdxUpnFLGc4FbgQxneDt2lj4be4PU1QSqM3ellKoAEekD3APkOkPfAXcC/zDG7HEprArR5K6UUg4R6QDci33oCGwFzM3AHGNMRPXZ0OSulKrznBr167E16gnALuAu4D5jzK6yvraY1x+sD3QA8kMB3481FWt5aSmkUqpOE5FjgKXATdjE/jSQYYy5rQKJ/QIgH/gCWO31B5/2+oM10qe9vDS5V8GhGoKVZfz48a2mTJnSorLvV7IZWGn8fn9aZc9fmqeffrpxZdscKxWuRCRWRG4CPgQ6A98AfYwxFxtjyv1ci9cfbAX8HfgEOB/byvdC7PKOazS5R5lp06a1rO5zLliwoPHnn39ev7rPq5RbnLX1Rdj1dA9wP9DVGPNORc7j9QfjgReAeODiUMA3NxTwXY9N7Fd6/cHTqzXwCtDkXorJkye3uP3225sDXHbZZW179eqVCfDyyy8nn3322e2KjztUy9/S2uSWtGLFioQTTjgho3Pnztm5ubmdli1b9ruZ8YYNG2KOO+64jOzs7Jzzzz8/vWTV1aFa/Y4ePbr1vn37PFlZWTn9+vVrV9pxBQUFDBgwwJuRkdE5MzMz55ZbbmleWkxvvPFGg4ULFzaePHlym6ysrJwVK1a4+qumUlVRvLUd9knSXsBa4BSnRUBlqmDuB3oDQ0MBX8mH924AvgRmurY8Y4xx/QNYcvDYp59+GjLGLDHGLAFMTXwUn/9QHwsXLvzqtNNO22KMWXL00Ufv6NKly669e/cuHT9+/NqpU6eGiuN69tlnvzHGLBk5cuSGa6+9dq0xZslJJ520bdq0aT8YY5bcf//9P5x88slbjTFLxo0bt+7GG29cY4xZ0qtXr/zPP//8C2PMkjfffPOrnj175h8cwyWXXLJxwoQJa40xS+bOnfsNYNatW/epMWbJhg0blhljluzYsWNpx44d96xfv36ZMWZJ/fr1C0ue41DHvfvuu1/27t17e/ExmzZtWlZWTP379/9l9uzZ35X1/6u8H87fq+vXnH7UvQ+gOfBKiRzwDyClsudLn5h3TvrEPJM+Me+eUl7/k/P6X2vwe/pd7iz+0Jl7KY4//vjdX3zxRYOtW7d6EhISzDHHHLNz0aJFiYsXL07u27fvTiiz5W+pbXIBtm/f7lm2bFnSoEGDOmRlZeWMHj06/eeff447OIYPP/wwediwYZsBBg8evL1hw4a/9qg5XKvfso7Lysrat2bNmoRLLrmk7QsvvNCwSZMmheWNSalIJCL9sDc7zwK2Yfupn2+M2VKZ83n9wQ7YB5k+ASYd6phQwPcGsBgY7fUHq32P1MOJiFJIY8zS2n7PhIQE06ZNm30zZsxI7dGjx86uXbvuWbhwYfLq1asTunfvvhcq3/K3sLCQ5OTkguJNOMpSfP6SytPqt6zjmjVrVrh8+fIvX3rppYYzZ85sPm/evJRHHnnkx/LGpFSkEJEU4EHsDU6A/wCXGGMq3UbcWWaZj539DwoFfPvLOHwW8Dh26eaDg86Thv2BMzoU8D1f2XhKozP3Mhx77LE7Z8yY0aJPnz47TjnllB1PPvlks5ycnN2HSrglHa5NbkpKSlGbNm32z549uwlAUVERixcv/t0Ny169eu2YPXt2U4D58+c3zM/Pj4GyW/3Gxsaa4jbCpR23fv362MLCQoYOHbrt9ttvX/vFF18klhVTUlJS4cHtiJUKd85sfQU2se8BxmPX16u6P8RE4Gjg0lDAd7h9C57H7rx07iFeywRSsS0Nqp3+gy3DiSeeuGPTpk1xffv23dW2bduChIQEc9xxx5Xaz7xYedrkzp079/snnngitXiP1hdffLHxwccEAoF177//flJOTk7266+/3qhly5b7oexWvxdccMGm7OzsnH79+rUr7bhQKBR3/PHHd8rKysoZNmxYu1tvvfWnsmK64IILtkybNi0tOztbb6iqsCciKSLyD+BlIA14H1sJc78xpqjsry6b1x9MBsYCC0IB32E3tg4FfDuA14GBTmVNSZnOf1dVJabSaOMwVau0cZiqSSJyJrZTY0tgN3Y9/O/GmENuI1lRXn/wRmwjsV6hgO+jcn7NGUAQeDAU8I0tMT4V+CvQIBTwVSq+shqH6cxdKRXxRKSJiDwJvIpN7MWz9QerMbHnYEsc55c3sQOEAr5/AQ8BY7z+YOcSL2UC31Q2sR9OjSR3EfmLiDwmIi+LyJ9r4j2UUgpARHzAcuy2dnuxuyKdaIz5trrew+sPNgXysOvjYytxihux7Qke9vqDxXm3A7bbZI0od3IXkdki8rOILD9o/DQR+VpEvhURP4AxZoEx5nLsbiXnVWvESikFiEhDEXkMm3RbYcsOuxlj7q2u2XoJfwfaAP1CAV+5WxMUCwV8m4HrgOOdD7B19xurLcKDVGTmPgc4reSAiMQAM4DTsZvCDhGRnBKHTHZer4yioqKiWq8NVTXH+fus0g0tpQCcFYHlwHBsNco1wAnGmFK3eKwsrz94MTAYuLUiyzGH8Krz3+7O7D0V+Lmq8ZWm3MndGPMucHDBfw/gW2PM98aY/cBzwNnOI753Af82xnxSydiWb9q0qZEm+OhQVFQkmzZtaoT9B6lUpYhIYxGZha1AaQssAXJraLaO1x/si10vfwfbEKzSQgHfBmwy74rdTNsDbKpqjKWp6kNMrYGSZX4/AT2BMcApQCMR6WiMefjgL3T6O4xwPk09+PWCgoLhGzZsmLVhw4Yu6I3faFAELC8oKBjudiAqMjl16w9hl2D2Y1v03mOMKaiJ9/P6gz2Af2G7RZ4XCviq430+A7phl2SgBmfuVU3uh5pVG2PMNGBaWV9ojHkUW7KEiCw5+PXc3NyfgX5VjE8pFeFEpDl2zXuQM/QhMMwY81VNvafXH2wBLMDuxNQ3FPBV1wz7U+xG262dz8N25v4T9lejYm2w/zOUUqpKRESwT3bOAJpid0e6gWqsWz8Urz/YDNteoDHQsxoTO9hlpHigr/N52M7c/wtkiEg7bOvMwdhm9UopVWki0gJ4GPiLM7QQGG6MWQ3gNOJqA/wB6A4cBbTDLuG+j32AqR3wA3a2HAwFfId9KNLrD3YF/g9ogm0v8EW1fVNW8SpFcZ9392fuIjIX6AOkishPwE3GmMdF5CrszY0YYLYxZkUFznkWtktbowpFrZSKWiJyLjATO1vfga2EeSx9Yl6c1x+8ELgUyOV/eaMQ+wj/KmweuhCb5H8A/gQkAvu9/uA6bM77EVv8sQr4HLt+PwS77NMbW56YWwOJHSemLdibqvuAGnsKP2zbDyil6hYRaYZdWy9usvUGMDx9Yt4a7IrA7YAX++DP68BX2NWDz0MB368bbXj9wXqACQV8+7z+YAx2Vn8htnCjEFv0UbJku9gXwL+BGTW5wbXXHxyP3anpllDAd3NVzlVW7oyIlr9KqejlrK1fiN3VqHht/RrgkfSJeR2xZYgnYHunXwm8Fgr4Sn1eIhTw7S3x50JgmfMB/Lqkk479QXEkkOSccxm1437sD65yr3JUhs7clVKuEREvdm39VGdoITDCGPOD1x88D9sP/QBwLfBEWUm9Liord7qa3EusuZ9kjMlwLRClVK1ynm6/ErgDaABsBSYAc9In5sUAd2H7r38ADA4FfFXtwR6VwnZZxhjzKvDqoerclVLRSUQ6YbeoO9YZeh4YY4zZ6PUHU4B52IcgpwMTQgHfAXcijWy65q6UqhXObH0c9sZoArABGGWMWQDg9QczsE+EHoEtQ5zjTqTRQZO7UqrGiUg2drbeyxl6EhhnjNkK4PUHe2K7OwKcFAr4Pvj9WVRFaHJXStUYEYnDtrqdgn0ycx1wuTHmX/Br5coV2AqStcBpoYDvG5fCjSrakEspVSNEpBvwMXYZJh5b+dK5RGJvBbyEfWDpP9it6zSxVxNXZ+76hKpS0UdE6mNn6tdinxgNYVsHvFl8jNcfHILt8JjgHHefljlWL62WUUpVGxH5I3aGngEYbMXLJGPMTgCvP1gfO1Mfit05aWgo4FvlTrTRTdfclVJVJiLJ2Nr0Uc7Ql9jZ+uLiY5xuiy9jb6reht3ZqEZ6sStN7kqpKhKRU4DHsSWMB7APJt1pjNlXfIzXH8wEgthOjoNCAd+LbsRal2hyV0pViog0Aqbyvx3VlgKXGmN+003R6w/+EXvjtAg4Wcsca4dWyyilKkxETsfuhzsC2zJ3MtD7EIn9Umy/mE3YahhN7LVEq2WUUuUmIk2xNekXOUMfY7e8+02HQ68/GIvdUHoCNrmfGwr4ttZmrHWdqzN3Y8yrxpgRwHY341BKHZ6IDMDeKL0I2IstYTz2EIm9Gbbf+gTsFnlnaGKvfbrmrpQqk7Pl3d+Bgc7Qu9hKmN88cOT1B+OB4cCN/G+bujm1GKoqQdfclVKHJNal2B2PBmI30bgS26L74MTeDTurn4HdKelYTezu0pm7Uup3RKQ98ChwsjP0OnCFMSZ08LFef3AQMAfbk/0M7K5G7u8CVMdpcldK/UpEYoGrsQ8Z1Qc2A2OBZ81BO/t4/cG22FLIwdhNNQaEAr4NtRqwKpUmd6UUACJyFPZhpOKdff4BjDXGbCp5nNcfTMDeLL0Bu7R7K3BHKODbhwobWgqpVB0nIvWwidqPzQk/YZdggiWP8/qDcdh/rwFs75h/YndKCtVqwKpctHGYUnWYiPQBHgEynaGZwPXGmPziY7z+YCfgMuBioAXwDbbv+uu1GqyqEF2WUaoOEpFU7Hr5pc7QSmCEMWYRgNcf9AD9gb8CJwAF2J2SngD+rfuahj9N7krVISLiwSb0qUAKtnXA7cBUY8w+Z2ekM7E3VLtiyxr9wJN6szSyaHJXqo4QkS7Aw8BxztCbwGhjzCoArz+Yjt1A43RsUr8QeC4U8BW6EK6qIk3uSkU5Z2ekG7HtAmKBDcA4YJ4xxjjVL6OwM3ic12bo0ktkk4NKV90JQmSJMeaYwx+plKoIETkDuxtSe+zOSA9hd0ba7vUH62FvlPqxfdZfB67Q6pfIUVbu1PYDSkUhEUkXkZewG2S0B74AjjPGXJk+MW+X1x+8HPgW2zMmBPwJOF0Te/RwdeZeos79JGNMhmuBKBUlRCQeGI/doLo+sBO4CZiePjEP4Fzn8wzgQ2wf9re0XUBkKmvmrnXuSkUJETkRW6ee4wzNB8anT8zbjm34NQ67Fd5y4GzgVU3q0UtvqCoV4USkOXA39iEjsA8ZXZk+MW8htu/L/diHjxYBVwHBUMBX5EasqvZoclcqQjk165cDd2L7p+/Dbk49NX1iXmvgNeDPwH+Bc0IB32K3YlW1T5O7UhFIRI7GLsH0dIZeB65Kn5i3EVv2OAH7gNIY4CGtVa97NLkrFUFEpDG2Hn0UttptPXB16tkTX2yQdcIl2Jl7GvAsMDEU8K11K1blLk3uSkUAZwnmQuzaenOgELuWfnP6xLxs4CNsq94Pgb+EAr6P3IpVhQdN7kqFOWcJZjpwrDO0CHvDdAM2wQ8D1mGT/z+0AkaBJnelwpaINMUuwYwEBNgIXJd69sRnGmSdMAJ7IzUJ2wTs9lDAt8O1YFXY0eSuVJgp0bnxLqApdgnmQeDW9Il57YHFQA/gLeDKUMC30q1YVfjS5K5UGBGRrtj+L72dof8AY5wlmNuwDyP9gi7BqMPQ5K5UGBCRZOAW7OYYMdjOjeNbXf7oi3Epra7EthNoiG3Ze0Mo4NvmVqwqMugeqkq5SEQEOA+4F2gFFAHTJSHxxiPGzj8J+BLogK1jvyYU8C13LVgVUbS3jFIuEZEsbFfGk52hj4FRTi+YudhNM77Edmt8zZ0oVaTSZRmlapmIJGK7MV4DxAFbAH/aJfc/k5CWcQ0wCbtn6XhgeijgK3AtWBWxNLkrVYtE5HRgBtDOGZqV0LbzTWnn33U6tud6B2AeMEGfLlVVocldqVogIq2AB4BBztDncalHTGh12cxcYAnQEvgU+HMo4HvDlSBVVNHkrlQNEpEYbB+Yv2GrXXbHNGx2X6vhDyd64hL+CSQDb2Db9b6ppY2qumhyV6qGiEgu8AiQC+Cpn/x2iyF3bolv5vVjnzh9HpgaCviWuRimilKa3JWqZiLSiP89cOSR2Phfmpxyxdrkrn/ug932bjrwYCjgW+1imCrKaXJXqpo4NesDsa0CWoIUJXU/fWuTPpemeuLrHwD8wCP6AJKqDZrclaoGItIeWwVzGkB8y8yipqdd5Ylv3n4ldqb+Yijg2+9mjKpu0eSuVBWISDyemGsRzxRMUbwnoQGNT7ykIKnrqc+KJ2Z6KOBb6naMqm7S5K5UJdX3djvH06DJI0W7tjYDSMw8dlej489/ML6Z98FQwPez2/Gpuk2Tu1IV1Owv/l67lr/1zN7Vn3UAiElqujsx+4Q7U/oOv1P3KlXhQpO7UuWUfu2CTts/enHOjiWv9Crakw/iKYpJajqzcMem6/I/fmmP2/EpVZImd6UOw+sPZuz98Yt7ty16+qx9P31pB2Ni36ew4PKC/J+/cjc6pQ5NW/4qVQqvP+gt3LPjtvwPn78gf8nLQlEhiGzCmPEUFjxrjNGnSVXY0pa/Sh3E6w+mGVM0efeqxSO3vjUrtjB/E4ABHsaYG4wxW10OUanD0mUZpRxef7AJcN2BzWvGbln4aMLe0DJxXloGXGGM+djF8JSqEE3uqs7z+oNJwF+L9u+5bvsH8xrlf/zPIkyRANuAG4BHjDFaBaMiiiZ3VWd5/cF6wEhjzKTdK99rvuWNh/YW7ckH8ACPA9cbYza5GqRSlaTJXdU5Xn8wFhgKTNn/y49tN//r/q37138DUA/bW/1KXYJRkU6Tu6ozvP6gBzgXuLVo786MLW/N2rDri4VFQBPsVnfXA4/rEoyKBprcVdTz+oMCnAncbkzRUTuW5v209e0ntlF4IA0oAh4CbjTGbHY1UKWqkSZ3FdW8/mAW8Chwwt41y9dsWnDnd0W7t3dwXn4PGGOM+dS1AJWqIZrcVVTy+oNxwHXAlMKdW3b//MIti/dv/K638/I64Fpgrj6IpKKVJncVdbz+4DHA46aw4Kht7z71af7HL3UA0xs4ANwL/M0Ys9PdKJWqWZrcVdTw+oOJwK3AuD3fL92y6eW71pj9u7s5L/8LGGuM+ca1AJWqRZrcVVTw+oN9gccK8je1/+WVqd/vW/tVe+el77BJPc/F8JSqdZrcVUTz+oONgXtMwYHLtn/w3ObtH87fizHtgT3A34B7jTF7XQ1SKRdoclcRy+sPHgvM3f3tx202/+uBbUV78ps6Lz0PXGOM+dHF8JRylSZ3FXGch5EmHtiy9rYt//fQvr2rP/UAjYEvsaWNb7kaoFJhQJO7iihefzCtcM+Of+R/9MJJ+R+/VIQpSgR2ADcBfzfGHHA5RKXCgiZ3FTHSr3np1J0r/jN/26KnGxbt2ga2wddsYJIxZqOrwSkVZjS5q7Dn9Qfjdq1877H8/750yf51XxcPfwT8VRt8KXVomtxVWEvpOzx337qV/9r99QfN7WZIbAAmAs8YY4rcjU6p8KXJXYUlEYmr1z73kX0/fXmp2b8HkELgPuB2Y0y+y+EpFfY0uauwE5ucenZMUsqcvd8vbQwg8fX/Y/bvucIYs8rl0JSKGJrcVdgQkYyYBk2eLNy1tTeAJ7HRVlOw/9Kifbtfdjs2pSKNx+0AlBKRhp64evchnpWFu7b2lrh6RfEtO80s2r09TRO7UpVT7TN3EWmP3VS4kTFmYHWfX0UPEfEAFxMTe68p2JcCkHDEUSvEE3PWnh8++cHl8JSKaOWauYvIbBH5WUSWHzR+moh8LSLfiogfwBjzvTHmspoIVkUPEemFeD4GnqCwICW+RYe9DXsNunzv6s+6aGJXqurKuywzBzit5ICIxAAzgNOBHGCIiORUa3Qq6ohIKxF5CliMKcqNSUqh0XFD/t30zPGp2xfPn+V2fEpFi3Ityxhj3hUR70HDPYBvjTHfA4jIc8DZ2P4ehyUiI4ARzqep5YpWRSwRqQeMBSYDDYiJJbnbGVuSuv558LrHr3zD3eiUij5VWXNvDawp8flPQE8RaYpttdpdRK43xtx5qC82xjyK3dsSEVlShThUGBMRAfpha9TbA9TP6EWjngNnJ7TOGhMK+Ha7GqBSUaoqyV0OMWacHeSvqMJ5VZQQkc7A/cCfAOKaHkHjPkPXJnbsMTgU8L3nbnRKRbeqJPefgLYlPm+D3XhY1XEikgLcDIwGYiQ+sajxHy+S5K6nPiix8TfobF2pmleV5P5fIENE2gFrgcHA+RU5gYicBZwFNKpCHCpMiEgscDlwG9AUxCR1P53Gxw35IaZBk4tDAd8HLoeoVJ1RruQuInOBPkCqiPwE3GSMeVxErgJeB2KA2caYFRV5c2PMq8CruuYe+UTkJOAB4CiAhNbZ+1L+PDouvnm7e4EpoYBvj5vxKVXXiDHG7RgQkSXGmGPcjkNVnFNFdQ8wAMBTv+HOlFOvTErMPHaliFwaCvg+dDVApaJYWblTe8uoShGRBoAfuBZIQDz7GvUatK9h73OTPHEJdwE3hwI+3ZhaKZdoclcV4pQ2DgGmYsthSTjiyO9SfRM6xDZM/R4YGgr4dAMNpVzmanLXG6qRRURygWnAsQCeeknfNjtnUqN6RxzVDrgTuFVn60qFB1eTu95QjQwi0gK4A7gUEJBNjY4bsqrRcYOPE/GsAHyhgO+/7kaplCpJl2VUqUQkHhgDTAEaAgfiW3R4tcXgO3p76jXohX0S+bZQwLfPzTiVUr+nyV0dkoicjn26tBMAMXH/l3Z+YE9Cq079gS+As0IB31IXQ1RKlUGTu/oNEcnE9oHxOUOrGvYcuKBJn6HDsPdGbgX+Fgr49rsVo1Lq8PSGqgLsbkjAjcDVQBywI7ZRi2kthz+U64mNvw74CBgeCviWl3UepVR40BuqdZyzG9Il2GqXFoABeaLl0Ae/i2/R/nrnsKuBGaGAr9CtOJVSFaPLMnWYiPTCljb+wRlanNT11AeanjZmHLYy5t/AqFDAt9qtGJVSlaPJvQ4SkVZAALjIGVoncfUmtR03v72I5xlgO3ABMDcU8Lnfn0IpVWGa3OsQEUkAxlG8GxLsB+5pOezv/4lv5p0GZAPPAONCAd8v7kWqlKoqTe51QIndkO4FOjjDC+qld72pxeC/jQTeAFYDp4cCvtdcClMpVY20WibKiUg28CDObkjACmBs+sS8ekAQ2x/mAeDGUMC305UglVLVTqtlopSzG9JNwJXYfvvbgCmtLn/4xbiUNvcB5wHLgYGhgO8j1wJVStUIj9sBqOolIrEiMhr4Bvgrdq/bhyUuISN9Yt6OuJQ2y4FzsOvuuZrYlYpOuuYeRUTkZOwSSxdn6D/YJZidwFzgFOA94PJQwLfSlSCVUrVCk3sUEJGO2N2QznaGfgCuaX3lk/+OTWp6NbbxVwEwCng0FPAVuROpUqq2aHKPYE7LgBuAsUA8sAu4Pb5FhwdbDn1wIPA10BZ4GbgqFPD95FasSqnapck9AjktA4Zie6y3cIbnAJPSJ+Z1AT4AugFLgYtDAd/btR6kUspVWgoZYUTkOGxpY64ztBi4On1i3n7gCeBUIITdCm++LsEoVTdpKWSEEJEjgLuAwc7QWuC6thNeXOSJTbgNuBhb7jgemKkbaChVt+myTJgTkUTgOuejPrAXmNrk5MsfbnjM2X8FHseWO94D3BkK+La6FqxSKmxocg9TTsuAwcBUoI0zPC+2SavJrUc8egZ2N6SmwNPYp0u1c6NS6lea3MOQiByDXVc/1hlahniuTr/ulTTgNWx/mDeBa0MB3zKXwlRKhTFN7mFERFpiK2CGOkM/A5Pajnv+W098/buBntgZ++nA69qOVylVGk3uYcBpxTsW2xIgCTgAPND83Nvm1W/X/Ubsw0lrgWHAU7ojklLqcDS5u8hZVz8b24q3vTP8SoOcPoHUs665GLtv6W5gEvBgKODb7U6kSqlIo8ndJSJyJLYPTF9n6MuY5FR/m9Fzjsb2V08AZgK3hQK+Te5EqZSKVPoQUy0TkabArcAV2K6cW/HE3Nx27Pz9nriER4E04AVgUijg+8bFUJVSEUwfYqolIhILjARuA5oAhcDf0y66b3FCq8zJ2C3u3gf6hwK+xe5FqpSKBrosUwtEpA8wDTjSGXqz0XFDHmt8/AWjgKuAVdge6y9rBYxSqjpocq9BIpIO3A0McoZC8Wkd70i7+L6TRDzPYUsdRwOzQgHfAbfiVEpFH03uNcBpGXAt4AfqAXskrt49rUfPiY+plzTNOex2YGoo4NvhVpxKqeilyb0aOaWNA7CljUc4w8+1GHzHknrpR12DvVn6LHB9KOBb41KYSqk6QJN7NXFKG6cBfZyhz5Jz+z2ecsqIS7E9Yj4E/qJ7liqlaoMm9yoSkRRsaeMobGnj5rimbR9oeen0oyUmdhqwBttbfZ7eLFVK1RZN7pUkIjHACOzaeQpQKLEJj7YaPnN/bKMWU4D92HYC94UCvj0uhqqUqoM0uVeCiJyIXYI5yhl6O7Xfde81yP7jFdg2vE8Ak0MB33q3YlRK1W2a3CtARNpiSxvPc4ZWN+h80lNNfeP7i8hk4B1gnLbhVUq5TdsPlIOI1Od/pY31gT2xjdNmtbx0eoYnvv6NwHdAf2CBrqsrpcKBth8og1Pa2B9b2pgOILHxC9IuvGdrfIv2o4Fd2KQ/XfcsVUqFE12WKYWIdMHuhlTctfGLJqeMfKdh7lkXYH/TeAS4STs2KqXCkSb3g4hIE+AWbFuAGGBL/Q495jU7Z9LJEhN7FfB/wIRQwLfczTiVUqosmtwdTmnjcOBv2IqXopjk1HlpF9/XMjYpZRTwNeAD/q3r6kqpcKfJHRCRE7Cljd0AiIl7v/nAKRvqe7sPArYBfwUe1uZeSqlIUaeTu1PaOBXbHgDgx4a9z13U+ISL+olIT2zCvzUU8G11LUillKqEOpncRaQecA1wPZAI7E1o2+XV5gNv+oMnvv4FwKvANaGAb5WbcSqlVGXVqeTulDb+BVva2A7Ak9h4Ydr5gYZxTdsMAr4ATgkFfG+6F6VSSlVdnUnuIpKDLW08BQBPzMpU37jVDXL6nIrdNGMEMDsU8BW6F6VSSlWPqE/uItIYuBm7nV0MsDWp2+kfpPzpipPEE9MeuAu4IxTw5bsXpVJKVa+oTe5OaeNl2NLGVKAoPq3jm80H3pwd06CxD3gemBgK+H5wM06llKoJUZncReQ4YDrQHcBTL+nT5gNvik1onX0ysBQYHAr4FrkZo1JK1aSoSu4i0ga7zHK+HfCsa3LSsFDyMWcfKyLrgEuAZ0IBX5GLYSqlVI2LiuTulDaOB27AljbuS8w6/uOmZ4zt4Ymr1wS7U9LUUMC3y804lVKqtkR0cndKG8/Glja2B4hNabO0+cCbjohr0vIE4Blgkm5GrZSqayK2n7uIZGFLG/8MIHH1Qqn9ri1I7NgzF1gMnBkK+D6uzniVUipSRFw/dxFpCEwBrgZiEclv2GtQqPHxFxwlnpgfsa0E5mtzL6VUXRYxyzIi4gEuwt4wbQGYeuldP0vtd21OTGLj9tj19vt1M2qllIqQ5C4ix2BLG3sBeBo0+bZ5/8mNE1p16ordjHpSKODb4GaMSikVTsI6uYtIc+AOYBggeGJ+adJ3+Lbko8/sKCIfAz5dV1dKqd8Ly+QuIrHYnZBuxd5sPZCYdcKXTU8bc5QnIbEIuBR4SuvVlVLq0MIuuYvISdglmM4AsY1bftVswI2t41OP6Aw8ANwSCvi2uxiiUkqFvbBJ7iJyBHAPMAhAYuPXNT396gOJ2X/MFpGFwNWhgO9LV4NUSqkIES7JvSWwEqgP7Enqfsa3KX2HHymx8SGgP7BASxuVUqr8wiW5twKIT+v4RbNzJmXENmyeAdwE3K2ljUopVXHhkdzFs7/5wJs3129/9JHAi8CEUMC32u2wlFIqUoVFco9v3i6+fvujtwIX6RZ3SilVdWGR3At3blkDdAsFfAfcjkUppaKBGOP+fUoRWWKMOcbtOJRSKpKUlTs9tR2MUkqpmqfJXSmlopAmd6WUikKa3JVSKgppcldKqSikyV0ppaKQJnellIpCmtyVUioKaXJXSqkopMldKaWiULX3lhGRBsBMYD/wtjHm2ep+D6WUUmUr18xdRGaLyM8isvyg8dNE5GsR+VZE/M5wf+AFY8zlQL9qjlcppVQ5lHdZZg5wWskBEYkBZgCnAznAEBHJAdoAa5zDCqsnTKWUUhVRrmUZY8y7IuI9aLgH8K0x5nsAEXkOOBv4CZvgP6WMHx4iMgIY4XzaWUSWlBFCI6C0TbFTgV8O8y2Eo7K+p3B+r6qcq6JfW97jy3Pc4Y7Rayx83quy56qp66s8x7p1faWX+ooxplwfgBdYXuLzgcCsEp9fBPwdaAA8ATwEXFDOcz9a2deBJeX9HsLp43Dfc7i+V1XOVdGvLe/x5TlOr7HIea/Knqumrq/yHBuO11dVbqjKIcaMMWYXcGkFz/VqFV+PRLX5PVXne1XlXBX92vIeX57j9BqLnPeq7Llq6voqz7Fhd32Ve7MOZ1kmzxjTxfm8N3CzMeZU5/PrAYwxd9ZMqKXGpRt9qBql15iqSTV1fVWlzv2/QIaItBOReGAw8Er1hFUhj7rwnqpu0WtM1aQaub7KNXMXkblAH+zC/0bgJmPM4yJyBvAAEAPMNsb8rSaCVEopVTFhsYeqUkqp6qXtB5RSKgppcldKqSgUdcldRBqIyJMi8piIXOB2PCq6iEh7EXlcRF5wOxYVnUTkL07+ellE/lzZ80REctfeNqomVeT6MsZ8b4y5zJ1IVaSq4DW2wMlfQ4HzKvueEZHc0d42qmbNofzXl1KVMYeKX2OTndcrJSKSuzHmXWDLQcO/9rYxxuwHDu5tAxHy/Sl3VfD6UqrCKnKNiXUX8G9jzCeVfc9ITn6t+d8MHWxSbw38ExggIg8RnY+Uq9pxyOtLRJqKyMNA9+KnspWqpNJy2BjgFGCgiFxR2ZNX+2Ydtag6e9sodbDSrq/NQKX/wSlVQmnX2DRgWlVPHskz95+AtiU+bwOscykWFX30+lI1rUavsUhO7uHS20ZFJ72+VE2r0WssIpK709tmMdBJRH4SkcuMMQXAVcDrwFfAfGPMCjfjVJFJry9V09y4xrS3jFJKRaGImLkrpZSqGE3uSikVhTS5K6VUFNLkrpRSUUiTu1JKRSFN7kopFYU0uSulVBTS5K6UUlFIk7tSSkWh/wcz/cKHPynnMAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "traj = data[0] # pick some trajectory\n", "msd_traj = nl.analysis.MSD(traj)\n", "msd_data = nl.analysis.MSD(data)\n", "\n", "plt.loglog(msd_traj, color='tab:blue', label='single trajectory')\n", "plt.loglog(msd_data, color='k', linewidth=2, label='whole dataset')\n", "plt.ylim([1, None])\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are running this notebook live, you might have noticed that this time the calculation of the ensemble MSD was almost instant, while for the first plot it took about a second. What's different here?\n", "\n", "``analysis.MSD()`` caches the MSD calculations it does for any ``Trajectory``—in that ``Trajectory``'s ``meta``-dict:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "traj.meta: dict_keys(['MSD'])\n", "traj.meta['MSD']: dict_keys(['N', 'data'])\n", "\n", "nl.analysis.MSD(traj) is traj.meta['MSD']['data']: True\n" ] } ], "source": [ "print(\"traj.meta: \", traj.meta.keys())\n", "print(\"traj.meta['MSD']:\", traj.meta['MSD'].keys())\n", "print()\n", "print(\"nl.analysis.MSD(traj) is traj.meta['MSD']['data']:\", nl.analysis.MSD(traj) is traj.meta['MSD']['data'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When ``nl.analysis.MSD()`` is called on a ``Trajectory``, it checks whether the ``meta['MSD']`` entry exists; if so, it skips the whole calculation and just returns the stored value. This means that once you calculate MSDs on all your trajectories, running downstream analyses (e.g. MSD grouped by different conditions) is fast. So don't feel bad about calling ``nl.analysis.MSD()`` on the same data 5 times in a row! It will do the actual calculation only once." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Time averaging\n", "By default, ``analysis.MSD()`` calculates a time-averaged MSD (TA-MSD). To prevent this:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqQklEQVR4nO3dfbyUc/7H8dd35pw63V6VTuKIC0WlhPx2sdZq5a7BUiGnrdZNsWTdMyhdhYxFLMImN5FYcpOM211Le2ORVkqKcKHoThk3dTo3c/3++M6x3TvVOeeamfN+Ph4eHjNzzcxndi+frr7X5/v5mCAIEBGR/BIJOwAREal9Su4iInlIyV1EJA8puYuI5CEldxGRPKTkLiKShwrCDgDAGLMc+CzsOEREcsxuQRAUb+qFrEjuwGdBEBwYdhAiIrnEGDNzc69pWUZEJA8puYuI5CEldxGRPJQta+4beeedd9oVFBRMBLqhP4Q2JQ3MraysPKtnz57Lwg5GRLJL1ib3goKCie3bt+9SXFy8KhKJqLvZBtLptFm+fHnXJUuWTAROCDseEcku2XxF3K24uPhbJfZNi0QiQXFxcQr7NxsRyUHjRwzqUVefnc3JPaLEvmWZ/32y+f9DEdmM0VcPv/D06IvvPjHi+Hvq4vOzdlkmm5WUlHRv3759+TvvvLOg+rnOnTt3raqqMh999NH73333XWTgwIG7zZ8/v0kQBKZly5aVr7766keO46Sj0WjPTp06ramsrDTRaDQ47bTTvh45cuTSaDQa5k8SkXo055oe+5xTkBr3A00qPg/a3VkX36Hkvo1++OGH6MKFCws7duxYMWvWrKJ1Xxs7dmy7du3aVTz77LOfAsyePbtxo0aNAoDGjRun58+fPw9g8eLFBSeffPIeqVQqeuutt35Z/79CROqd5xzUyRS+sobGZkJl7Owrrp8wty6+JtS/0htjjjfGTACcMOPYkrvuuqtN9+7du3Tu3LlraWnpbpWVlQCceOKJKx966KE2AA899FCbfv36rax+z1dffVVYUlJSUf24R48ea5s0abLRElNJSUnlxIkT/QceeKBdOp2uh18jIvXJjSdjbjw5xo0nn+gYnzYIzxmUDszflwatm5WWX/3m3VW/ebCuvjvUK/cgCKYD07e0hRbgsqmzO3y45Lumtfnde7Vvsfqm/j2+2NIxs2bNKpo6dWqbmTNnzm/cuHHw29/+dtd77rlnB4DS0tJVQ4YM2X3MmDFLX3rppVaTJ0/+5PHHH98BYNiwYSuOO+64vaZNm9b6sMMO+3bo0KFfd+/efe2mvqNr167l6XSaxYsXF3To0KGyNn+jiNQfN55sBpT5iVhV5nG/ItZOPTQyl2Ojb5b/PDK/P8DcwK04q/zSJcto3d9PxOrsvqKWZbbgxRdfbDF37tymPXr06AJQVlYWadeuXSVAcXFxleM4lRMmTGjdsWPHNc2bN//x0vuQQw5Z8+mnn8555plnWr7yyistDznkkC6vv/76/AMOOKBsU9+jObYiuc2NJ0cAY4Dv3Hjy3b3N58ENBS8ddmL0X+kmpjwSBPwwO9hj8Y0VA3ZPpn/+TZrokX4itqguY8qJ5P5TV9h1JQgCc/LJJ389fvz4xes+/5e//GUHgP79+6+6/PLLd7vrrrs+3fC9juOkhwwZ8s2QIUO+GTx4MNOmTXM2ldznzZvXKBqNUlJSoqt2kSznxpMdgNZ+IvZe5rEDPAocC7zWxXy29OyC6b+MRd7cCahKYx4BHjWGV/cb/d+K6fHkbsDXfiL2fV3HmhPJPSzHHHPMt3379u141VVXLS0pKalcunRpNJVK/VjWMnDgwFVfffVVYd++fb/97LPPCquff/nll5vtv//+ZcXFxVVlZWXmww8/LOrVq9d3G37+l19+WTB06NDdTj/99GWRiCoaRbKVG0/uAjwL7J95PB3oAnSMkK48K5qcfEXBY22iJjgZu3v8XmA0Xmq9Qgk/Eau31uZK7lvQs2fPshEjRiw+4ogj9kqn0xQWFga3337759Wvt27dOn399dcv2fB9H374YdHw4cN3A7uTtHfv3qkhQ4asAli7dm2kc+fOXatLIU899dSvR40atbT+fpWIbIOrsBsGLwH2Avq7Zsm7Z0afr+gb/UdxM7P2t8AS4EbgbrxUKKsN6zLZsN5rjJm5YT/32bNn+z169FgRVky5Yvbs2W179Ojhhh2HSD5y48koMAC4H3jALyodDhwJnA0cj604/Cv2Sv1pvFTF5j6rLmwqd1bTlbuIyCa48WRb4B9AZ9csmf1Uo2uqgEXAjsAy4I/ARLzUxyGGuVlK7iIi63DjSQMcA1y2Myv2mNjo5le7mM8PM4ZuwHPAQ8BzeKnyUAP9CUruIiIZbjzZHni4JT/0PrdgWnpo9Pl01KQPBe4B/pgNa+k1peQuIgK48WTHKFV//230rzteWTDlh8ZUNDWGR4GReKl6q3KpLUruItLgufFkm5+bef8cU/hg270ji6LAP4FL8FL/DTu2baXkLiINm+c0u6fw//51TPTtHcuDgqXAecBTeKnwSwm3g3bObMaSJUuinTt37tq5c+eubdu27dGuXbt9qx9/8cUXBQUFBQfcdNNNbTf3/osvvnhnY0zPuXPnNq5+bvTo0e2MMT1nzJjRFOC2227bYa+99uq61157de3UqdM+kydPbgXQr18/t6SkpPvee+/d1XXdbieddJL76aefFm7mq0RkG1WOan3wyqDFF8dE3+48o6r7W41MZUe81JO5nthByX2z2rdvXzV//vx58+fPnzd48ODl55xzztLqx4888kjrHj16/PDEE0/ssKXP6NSp05rqzpEA06ZNa7PnnnuWAXz88ceFt9xyy05vvPHGgg8//HDezJkzPzjwwANXVx973XXXLVqwYMG8Tz75ZO5+++23ulevXnuXlZWZuvvFIg2I5xTiOaMjpP+1msatLyg/79nBFVcegpeq87YA9UXJfRs88cQTbW6++eYvlixZUrilK+o+ffp88/zzz7cC20OmRYsWlW3atKkE2xa4WbNmacdxqsD2ouncufNGpVWRSIRRo0Yta9u2bcXUqVOztjWySK7498iDD1watPoIuOaZ9KGcsPa65LT0L06q7uaYL3Jjzf2Z8zqwbF6ttvylXdfVnDh+q8uaFi5cWLhixYrCXr16rT7hhBNWTZo0qY3neZtsH9CyZcuqnXfeufztt98umjp1aqv+/fuvevjhh9sCHHTQQavbtm1b0aFDh+6/+MUvvuvbt++q0tLS1Oa+d9999139wQcfFG3udRH5CZ5jyoLC83pGgj+tpigyvHx42XPpQ+4H4n4ilncDFXTlvpUmTZrU5oQTTlgFMGjQoJVTp05ts6XjTznllJUPP/xwm2Qy2XrgwIGrqp8vKChgxowZH02ZMuXjTp06lcXj8Q4XX3zxzpv7nGxoEyGSszynfTowySJTcceb6S6Rs8svOu+59CFN/UTsPD8R26ipXz7IjSv3bbjCritPPvlkmxUrVhQ+9dRTbQCWLVtWOGfOnMYTJkxo+8orrzgA1WP0AAYMGPDNNddcs0v37t1Xt2nTZr2rg0gkQq9evVb36tVr9bHHHvvtWWed5Y4bN26T4/bmzJnTtHfv3hs1KRORLfAcA5QGAbdXEm11fcVAJlf1vuzjxAl3hR1aXcuN5J4lZs+e3Xj16tXRZcuWvVf93EUXXbTzQw891OaOO+5YDCze8D3NmzcPPM9b1LVr1/UmMfm+X7ho0aLCQw89dDXAzJkzm5aUlGy05p5Opxk7dmy75cuXF/br1+/bOvhZIvnJc9r/EDR+oJlZe8zcYPe1F1acG/k4KLnAT8RuDzu0+qDkvhUmTZq0Q58+fVat+9yAAQNWlZaW7nHTTTd9tbn3DRs2bNWGz5WXl5tLL710l6VLlxY2btw4aNOmTcW99977YzvhESNG7JJIJHYqKyuL7L///j+8+uqrC4qKirQ2I/JTPMdUBNFBEB1fSFWzRMWAtfdWxf5WRXSqn4g9EHZ49UUtf3OcWv6KrMNzdgTuBk56O70X8YqhZR8HJb/xE7GXww6tLqjlr4jkN7u2fipwZxDQ/KbKU1b/uer4v1URLa2PkXbZqE6qZYwxJxpj7jXGTDPGHFUX3yEiAlRfrU/FzjL9+OKK3195V9WJTauIjm+oiR22IrkbY+43xiwzxszd4PljjDELjDELjTFxgCAIngmCYCjwO+yfpiIitc9z+gLvA8cB8YvLf9/r6fQvTwM+x05IarC2ZlnmQeBObKN6AIwxUWA8duzUIuBtY8yzQRBUlwKOyLy+LdLpdNpEIpHwbwpkqXQ6bbDDeEUaFs9pC9yBHYE3Cxjklk1ZDDwA/B9wWr7tON1aNb5yD4JgBrByg6d/BiwMguCTIAjKgceA3xjrRuCFIAhmbWNsc5cvX+5kEphsIJ1Om+XLlzvA3J88WCSfeE5/YB7QD7gGOMgtm1IMfAacgt1x+liIEWaF7b2hWgKsu8FoEfBz4HygN+AYYzoGQXDPhm80xgwDhmUebtRdsbKy8qwlS5ZMXLJkSTe0k3ZT0sDcysrKs8IORKReeE577EpAX+zV+hF4qTluPNkcu7KwHPi1n4ht6wVlXtne5L6pq+ogCILbgS1uFAiCYAIwAWw5z4av9+zZcxlwwnbGJyK5zlbCnB4E3BJgmgJXRkxws1s2pR3x5G3AoYALHKbE/j/bm9wXAR3WebwLsMnt8yIiW81zOgD3Ake/G3Rce0nFOY0+CXbuD3QDjgeKgPeAC/1E7B8hRpp1tje5vw10Msbsjt16PwAo3e6oRKRhs1frQ4GbqwJTOLpycHpy1ZEfpYlMAs4BDgNeAEb4idjCMEPNVjVO7saYR4HDgbbGmEXAqCAI7jPGDAdeAqLA/UEQvL8Vn3k89k9f9SkXEcterU8EjloZtJh1Uvlo97OgvQ8c6idiKeDmUOPLEVnbfkBEGhh7tT4E+FMQEL2tst+nf6rq2w3MV8DhfiL2YcgRZh21HxCR7OY5LvBn4ChgxqnlI996K+hyKXAVcJufiK0JM7xcpOQuIuHxnCjwB+A6bHnv+V3KHvjnGhq/BfzFT8RuCDW+HBZqcteau0gD5jndgPuwmyGTwLlu2ZTWwMvACmB4iNHlvFCTexAE04Hpm6pzF5E85TmNgDi2PUkKW2H3mFs25RfAdOB7oLefiKnl93bQzk8RqT+esz+2hHo0tpNjV7zUo27ZlBOBV4Bl2KqYBeEFmR+U3EWk7nlOIzxnNPAW0A74DV6q1C2bssKNJ+PAk8Bs4Bd+IvZZmKHmC91QFZG65Tk/w66tdwMmAxfgpVa68aTBdpo9F9t08AxVxdQeJXcRqRue0wQYA1wMfAUch5dKAmQS+x+xif1m4HI/EQt/000eUbWMiNQ+z/kVtidMJ2yDwMvxUikAN57cEXvF3h/b5VGJvQ6oWkZEao/nONgr8mHAJ9i2vK9Wv+zGk/sCrwLNsRuUblRirxtalhGR2uE5JwB3A+2xSy2j8FKrq19248kSbD37WuCXfiL2QShxNhBK7iKyfTxnB+z8hlJgDrYS5se/jbvxZDPsEsyVQCuU2OuFSiFFZNt5zknYAdWnAB5w4AaJfU9gPnZSUgQ40U/E3q33OBsgXbmLyNbznGLsTdFTgHeBo/FSs+HH5ZeWQCV2Y1IT4Ajg71pfrz+qlhGRmrNteU/DLsO0wLYQ+CNeqgLAjSf7AE8DjTLv+BY71/SdEKJt0FQtIyI14zkl2Bumx2N3mp6Bl/pxOI8bT8aAp7DLNOOwYzef9xOx90KItsHTsoyIbFlmQDU2YTcGLgVuw0tVVR/ixpOl2HX194Aj/URsVQiRyjqU3EVk8zxnT+wQjSOAGcBZeKmPql/O7DS9BLgJeB17w/SbECKVDSi5i8jGPKcAuAjbvbES2ybgz3ipdPUhbjzZHpv4TwAeBwb7idjaEKKVTVByF5H1eU4PbKOvnsA04Dy81GI3njyIePJCYA9gDdAdWwlzKXCrn4ilN/OJEgIldxGxPKcIGAlcDqwETgaedMum7EQ8+TpwGPAN8CZQhF2GudJPxOaHE7BsiUohRQQ851BgIrA39sboJZm2vK2AF7BX6xcC9/mJ2PchRSlbQaWQIg2Z57QCbsQ2+vKxm5FeBnDjySLsskwXIOYnYq+EFKVsAxME4W8YM8bMDILgwLDjEMlHbjy5M1DmJ2Irf3zSljeeDPwpCGi3MCh5YlD5lW8soU0H7FX6LoALFAOn+YnYY/UfufyULeVOJXeRPOTGkwcCJwHHAfsCaeDfQDIW+c938YIp53SIrOjmp3dcdUHF8ILZwZ4tMm9dC3wKfAEsApJ+IvZkCD9BamBLuVM3VEXyiBtPOtjNRmcAVcA/sDdIWzSm/Lizo8/dcF7BNMopwKsYzMNVR1ZUEX0GeA47uHqxql7yg5K7SJ5w48mjsTdFdwYS2EEY3wDVk5FOAVgRtHxpYmWfUQ9WHfMxsFLJPD8puYvkODee3B0YCwzAttc92E/E3gKquzf+Efgddrnl2Lajv3gxDsRDiVbqi5K7SI5y48kdgKuB4dhdpNcCY/1ErAzPiQJDsUm/BbYiZsy6k5Ekv6nOXSQHufHk/the6a2B+wHPT8QWA+A5PbHdG/8PeA27w3ReOJFKWFTnLpJj3HiyOzax/wD08hOxOUB1zfq12D4wy4CBwKN4qfBL4qTeaVlGJIe48WRX4G9AGXYIxsfrDNAYh61LHw+MxEulwotUwqbkLpIj3HiyM/Aqdn29OrHvBdyFbcn7NhDDS2nqkSi5i2SzzO7SE7EbknoBXwNH+EWlX+AxBrgC26HxPGxL3qrNfZY0LEruIlkos64+DuideWoBdiDGBL+odB8gCewOTME2+VoSSqCStZTcRbJIZofpaGx54zfANcBUPxH7AM/ZHfgTtsLsA+DXeKm/hxWrZDcld5EskBlX91vs1Xk77ISjEX4i9jWe0wSPUdh9R1XYdgJ/wkuVhxawZD0ld5GQufFkF2xd+q+wgzCO8xOxmXiOweNE4FZsh8a/AJfipRaFFavkDiV3kZC48WQT7A7Ty4HvsT3V7/MTsTSe0xG4AzgGeB8twchWUnIXCYEbT/bCNvnaA3gIuMxPxJbhOUV4XAFcCZRjh1SPx0tVhBet5CIld5F65saTQ7HLMB9j69XtFbnnHIXdgNQReAxbBfNlWHFKblNvGZF64saTEeB67I3RF4FT/UTsWzxnF+y6en/gQ+BIvNRfw4tU8oEmMYnUg8z6+oPYnup/Bob7RaUGO3R6FBAFrgNuxkutDSlMyTGaxCQSIjeeLMYOmj4YuAy4xS8qPQzbNqAr8CxwAV7KDy1IyTtK7iJ1yI0n+2GTeEugv19U+jq2Re/vAB84AS81PbQAJW8puYvUATee3BF7c7Qf8N8CKo9aWDT4AOykJAe4AbhOwzOkrii5i9SizE7Tgdg2Ac2Bq95u/Pvnik3qTuAw4J/AOXip90MMUxqASNgBiOQLN548BFsF8zCwoLv55CC/qLR5sUm9A3QHzgJ+pcQu9UFX7iLbIXOlfjR209Fh2Ja8Fy5oPHhhY1P5JLZz4yTgMrzU8vAilYZGyV1kG7jxZBToi03q+wOLgYsmF4597tDo3BuA27Dr673wUq+FFac0XEruIlvBjSf3BIZk/tkVu+nozFjkP1PGN7p9GDALKMT2jLlZnRslLEruIj/BjSdbACdjyxd/CQTAy8AlwNN+Uel+2BulPYGXgPPwUh+HEqxIhpK7yCZk1tIPwXZq7A80xV6lXwlM9hOxRXhOC+AW4HxgGTAAeBwvFf62b2nwlNxF1uHGk62AQdik3g34DngEeAD4j5+IBZk+6/2w5Y47Y5uAXY2X+iaUoEU2QcldGrzMVfrPgLOxV99NgJnAUOAxPxH7/seDPWcP4E7gWGA20A8v9WZ9xyzyU5TcpUFz48l2wFPAL4AfsDXqf/YTsVnrHeg5jbF9Ya4GKoGLgTvwUpX1GrBIDanlrzRYbjxZAvwV2A07kPphPxH7dqMDPacXdullb2AqcJFG3Um2CzW5B0EwHZhujJkZZhzS8LjxpAv8DSgGjvYTsX9sdJDn7IQdWD0Q+ATog5d6oR7DFNlmWpaRBseNJzthE3sLoLefiL213gGeU4itgPGAxtg+62PxUmvqN1KRbafkLg2KG092wy7FRIDD/URs9noHeM7h2Bum+wAvYPusf1TPYYpsNyV3aTDcePIA7OajtUAvPxH74McXPacEuBlbLeMDJwLPqmZdcpWSuzQIbjx5MPZK/BvgCD8RsztIPacRcAFwDbZtwGjgRi3BSK5Tcpe858aTRwJPA19hE/vnAHjO0diNSHsD04EL8VKfhBWnSG1Scpe8lalhvxHbE2Ye9ubpV5mNSLcCJwAfATG81POhBSpSB5TcJe9k2vGeg61yaQYkgOv9otIAj2uxm5EqgCuA29S5UfKRkrvklcza+nhsj/W/Auf7RaULsF0dbwY6AJOBK/BSX4YWqEgdU3KXvODGk8XYK/QzsIMzTgGm+kWl+wCvAocD/wVOw0v9K6w4ReqLkrvktMwSzFBgLHZT0k3AGL+otAC7rj4cSAG/B+7FS1WFFatIfVJyl5yVqYK5CegB/B0Y7heVzsdOSUpgWwvcA4zES30dWqAiIVByl5zjxpP7YatgjgI+JTMkwy8q7Qn8G/h55t/H4qVmbe5zRPKZkrvkDDee3BW4FjtMYxVwEXC3X1TaAvgzcBZ2ItJgYLJ2l0pDpuQuWS8zHelK7E5SsEsxN/hFpd9jSx7HYNfbbwVG46U2btsr0sAouUvWcuPJxsC5wAigNXaQxkg/Efs802P9duwovFeB8/FS80ILViTLKLlL1nHjySLgVGAUsDvwCnC5n4i9i+fshscT2KHVPtAPeFpLMCLrU3KXrOHGk7tj55ieCbTFzig92k/EXsZzmuLhYXeVBthGXzerwZfIpim5S6gyderHYOvQ+2AT9zTgLuDVTMuA6t2luwJ/AS7HS30eUsgiOUHJXUKR2VF6BvaGqAsswfaCuddPxL4AwHP2xXZtPBx7FT8ILzUjhHBFco6Su9QbN540wEHYm6SnAI2A17BLLc/4iZht4OU5O2D7qv8e239du0tFtpKSu9Q5N56MYGvTLwT2A74DJgD3+InY+z8e6DkF2DX3MUAr4G7gGrzUynoNWCQPKLlLnXLjyZ7Y9fOfAXOwyzCP+InY9+sd6Dm/xi7BVJc2XoiXmlO/0Yrkj1pP7saYPYCrAScIgv61/fmSG9x4sg12Df0c7K7RQdikvn7Joufsjr1Z2heVNorUmhold2PM/cBxwLIgCLqt8/wx2KutKDAxCIJEEASfAGcaY6bWRcCS3TJLML/D9n5pg91oNMpPxFLrHeg5zYA4dnBGFXaj0jiVNorUjppeuT8I3Ak8VP2EMSaKHYpwJLAIeNsY82wQBNol2EC58eT+2CWYg4B/Aef5idjs9Q7yHAOcBvwRKAGmYAdnLKrfaEXyW42SexAEM4wx7gZP/wxYmLlSxxjzGPAb7KzKn2SMGQYMyzxsW6NoJSu58WRrbEOv3wMrsC13H97EEswB2Cv5XwCzgFM1OEOkbmzPmnsJ8MU6jxcBPzfG7ABcD+xvjLkyCIIbNvXmIAgmYCsmMMbM3I44JCSZJZjB2KvwHbB/k7vGT8S+We9Az2mHPSfOxCb/s4AHVdooUne2J7mbTTwXBEHwNfYmmuSxTE/18cAh2N7pR/mJ2LvrHeQ5jbCTkEYBTbFdG8fgpdZffxeRWrc9yX0RdthwtV0ADRzOc5n2u2OA84CvgdOBh/xELL3egZ5zLDaZ7w28AFyEl1pQr8GKNGDbk9zfBjoZY3bHDiQeAJRuzQcYY44Hjgec7YhD6oEbTxZgl2DGYsfX3Y1tv7tqvQM9Zy9sUu8DfAjE8FLP12+0IlLTUshHsf092hpjFgGjgiC4zxgzHHgJWwp5fxAE72/hYzYSBMF0YLrW3LNXpmXASdia9S7Af4A+fiK2/vg6z3GAkdiBGquBS4E78FLl9RqwiAA1r5Y5bTPPPw/oqixPufHkr4EbsJVRH2A3Gj2zXhWM50Sxde3VV/T3ASPwUkvrO14R+R+1H5CNuPHkgdhkfSS2IuoMbGlj5XoHes6h2E1sB2Dr2vvgpd6p32hFZFOU3OVHbjy5N3b5pT+2ZPEibHOvsvUO9JwO2B2op2FvrJcCj6llgEj2CDW564ZqdnDjyV2w5YqnA2uw7XbH+YnY+oOmPacJtl1AHFsKey1wI17qh3oNWER+UqjJXTdUw+XGkzsAV2Jr0Q1wBzDWT8SWr3egbRnQD9vgazfgCew0JL8+4xWRmtOyTAPkxpPNsb3VLwOaY3sGeX4i9tlGB68/Dek9oBde6rV6ClVEtpGSewPixpONsf18RgDtgGeAEesNzKhmpyGNwe42/gY7PelevFTlRseKSNZRcm8AMkOoS7HJ2sWOtvuNn4j9Z6OD7TSkczLHtsR2eRylaUgiuUU3VPOcG08eCYzDTjiahR1j98pGHRthU9OQLsBLza2/aEWktuiGap5y48kibLniH4CPgFOBqRv1gIENpyF9mvn3MyptFMldWpbJQ2482Q07BKM79ko8vlGtOmxuGtIteKmNjxWRnKLknkcyfWDOw16FfwMc6ydiL250oC1tHIDtw74LmoYkkneU3POEG0+2Ax7AdmN8HjjdT8SWbXTgxtOQBmgakkj+UXLPA248eSx2zq0DnA+M38SIu3WnIS1H05BE8pqqZXLYBjdN5wBH+InY+tUtnlPI/6YhNUPTkEQaBFXL5Kga3TS105DGAZ2BF7HTkObXc6giEgIty+SYdW6a3gR8ix2c8cJ6B3nO3tik3gdbBnk8kFRpo0jDoeSeQzI3Te8HYtibpmf4idj/hmJ4TivgGuy6+2rgEuBOTUMSaXiU3HOEG08eg71p2ooNb5raaUhnYXux7wDcC4zES21cLSMiDYKSe5bL3DRNYGeTzgV6r3fT1HN6AbcB+wKvAxfipd6t90BFJKsouWcxN57cB3vTdF9sbfoVP940Xb9lwGfAycCTWlcXEVApZFba4k1Tz2mBHbBxMf9rGTAOL7UmnGhFJBupFDLLZG6a3gccx7o3TT0nAgwCbgB2Ah4GrsRLLQ4tWBHJWlqWySIb3DT9A3Cnn4gFeM7B2Fr2/wPeBPripTbuxS4ikqHkngU2cdP0SD8Rm4Pn7IJHAhgIfAkMBh7BS23ctldEZB1K7iHb4KbpHcAVflEpeIzEtuONYnvCJPBS34cXqYjkEiX3kGRump6LrXj5Foj5RaUvYKtebgJ2BaYCl+OlPg0tUBHJSUruIdjgpukLwOl+UenO2Dr1XwKzgcF4qdfDi1JEcpmSez1z48l+wN3Y4dN/+KDx7x5vYsqvB84AVmBnnN6nVrwisj2U3OuJG0+2wa6plwKzdjNLznq98cW9sY29mmAbfV2rVrwiUhu0iakeuPFkDNvvpdiQHvV+4zPfa2rWPg50BJLAJXipBaEGKSJ5RZuY6pAbTzrY4RinA3P+EH3qvIsLp54NjAbmA8fipTaecSoisp20LFNH3HjyKOxN053as/LmfzS+IFpoqh7HtuK9CBiPl6oINUgRyVtK7rXMjSdbYEsZz46Qnv9I4dgbDo7OGwYUY5P91WrFKyJ1Tcm9Frnx5K+ABwA3Fnljyu2Fd+4VNcEI4A2gD17qnXAjFJGGQsm9FrjxZFNgLHDBTnztP9F49Eu7mBWlwFfAb4EpasUrIvVJyX07ufHkwcCkRlR08gom/fu06Kv7GsOvsb1ixuKlvgs5RBFpgJTct1Gm2ddo4NKjI2+tuLXwrsVNTfkhwHTgYrzUwnAjFJGGTMl9G7jx5IHApD3Ml13/VHjnF90jfgdgASptFJEsoeS+Fdx4shEwojmrr7q04Ik1g6MvV0ZM0Aq4BLgTL1UeboQiIpaSew258eS+EdKTTor8c79RhQ+VtTSrmwP3A1fhpZaGHZ+IyLrUfuAnuPFkAXDFfmahd23h/UH3iA/wLvAHvNTbYcYmIrI5aj+wBW482WUnVky5uGDqficXzCAdmKXA5cBkTUMSkWymZZlNcOPJaDPWXHp29K/XnV/wdLQpayuBcRETXKfSRhHJBUruG3DjyU7HRd545qKCqV33jHzFmqDRKxETnIeX+ijs2EREakrJPcONJyO/irzrPVD40tW9orMj3wZNl1QF5swmo5c/H3ZsIiJbS8kdOPOqa/dJFLzzfN/ojF2riFZ+GbS5fmez8jqVNopIrmrQyT058qjIqqD5xEThzN8Vm5R5P73bv/YwX/XfefSnS8KOTURkezTY5P63kYfHOprlk/eOLmr1Ybrk23fTew4+8tq/Tws7LhGR2tDgkvuaUcU7fhLs9HivyOeHfU2LYHLlEQ+vwDn9wuse0EBqEckbDSe5e06jJUHrK1vCiE5mUcHjVb9a9HL6wOPuHztydtihiYjUtgaR3KtGtTr6B5re396s2vm1qn3TD1YdfcNr6f1H+ImYNiKJSF7K7+TuOR3LgsI7i0xw9NfpFlxdccaC6elDTvQTsflhhyYiUpfyM7l7TnPg6nRgLq0iEk1UDKiaXNV71Pc0vdFPxCrDDk9EpK7lV3L3HAMMTAfmjxET7PR0+lDGVfSfu5ji0/xEbG7Y4YmI1Jf8Se6e0xO4HTjkg2DXihHlZ1T9N+h0LTDWT8QqQo5ORKRe5X5y95x2wPVBwJnf02TtmMpBPFl12II0kSF+IjYr7PBERMKQu/3cPacQOBcYnQ5oPrnqyO9vqjy12Xc0vQEY7Sdia2s3WhGR3JGb/dw95wjsEkzXeeldvzi/4nzn46DkS2CIn4i9WQehiojklNxalvGc3YBbgH7fB0XLLqs4e+UL6Z/tAmYcMMJPxNaEHKGISFbIjeTuOU2wE5Di6cAwserYz26pPGW3tTR6B4j5idh/Qo5QRCSrZHdyt6WNJwLjAHdmeq+F55cP3+Mr2jYDhgL3a5epiMjGsje5e04X7Lp676+DFl/+oeL87/6V7rY7cCfg+YnYqnADFBHJXtmX3D3HAUYB51cGkbJbKk9ePKHquJIqon8HLvATsTkhRygikvWyJ7l7TgQYAiSCgOK/pg/4+IqKYR1X0nIVcAow1U/EgnCDFBHJDVmR3Js3oinwb+Dni4K2n59bfkHZe8GeHYBrgYSfiK0ON0IRkdySFcl97x0iXdYGhau8ysHLHqvqtWtA5BngEj8R+yTs2EREclFWJPflVc3W9lx7d+vvaboAGOQnYi+HHZOISC7LiuS+xBQXGppeCtzhJ2LlYccjIpLrTBCEf4/SRAvfDaoq9gs7DhGRXGKMmRkEwYGbei1S38FsUrpSAzRERGpRdiR3ERGpVUruIiJ5SMldRCQPKbmLiOQhJXcRkTyk5C4ikoeU3EVE8pCSu4hIHlJyFxHJQ0ruIiJ5qNYbhxljmgF3AeXAa0EQPFLb3yEiIltWoyt3Y8z9xphlxpi5Gzx/jDFmgTFmoTEmnnm6LzA1CIKhwAm1HK+IiNRATZdlHgSOWfcJY0wUGA8cC3QFTjPGdAV2Ab7IHFZVO2GKiMjWqNGyTBAEM4wx7gZP/wxYGATBJwDGmMeA3wCLsAn+Xbbwh4cxZhgwLPNwH2PMzC2E4ACpzbzWFljxEz8hG23pN2Xzd23PZ23te2t6fE2O+6ljdI5lz3dt62fV1flVk2PDOr922+wrQRDU6B/ABeau87g/MHGdx4OAO4FmwAPA3cDAGn72hG19HZhZ09+QTf/81G/O1u/ans/a2vfW9PiaHKdzLHe+a1s/q67Or5ocm43n1/bcUDWbeC4IguAH4PSt/Kzp2/l6LqrP31Sb37U9n7W1763p8TU5TudY7nzXtn5WXZ1fNTk2686vGk9iyizLPBcEQbfM44MBLwiCozOPrwQIguCGugl1s3FtdhKJSG3QOSZ1qa7Or+2pc38b6GSM2d0Y0wgYADxbO2FtlQkhfKc0LDrHpC7VyflVoyt3Y8yjwOHYhf+lwKggCO4zxvQBbgOiwP1BEFxfF0GKiMjWyYoB2SIiUrvUfkBEJA8puYuI5KG8S+7GmGbGmEnGmHuNMQPDjkfyizFmD2PMfcaYqWHHIvnJGHNiJn9NM8Ycta2fkxPJXb1tpC5tzfkVBMEnQRCcGU6kkqu28hx7JpO/fgecuq3fmRPJHfW2kbr1IDU/v0S2xYNs/Tk2IvP6NsmJ5B4EwQxg5QZP/9jbJgiCcmDD3jaQI79PwrWV55fIVtuac8xYNwIvBEEwa1u/M5eTXwn/u0IHm9RLgKeAfsaYu8nPLeVSPzZ5fhljdjDG3APsX70rW2QbbS6HnQ/0BvobY87Z1g+v9WEd9ag2e9uIbGhz59fXwDb/Byeyjs2dY7cDt2/vh+fylfsioMM6j3cBvgwpFsk/Or+krtXpOZbLyT1bettIftL5JXWtTs+xnEjumd42bwB7G2MWGWPODIKgEhgOvAR8ADweBMH7YcYpuUnnl9S1MM4x9ZYREclDOXHlLiIiW0fJXUQkDym5i4jkISV3EZE8pOQuIpKHlNxFRPKQkruISB5SchcRyUNK7iIieej/AZvRTuVwBivaAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "e_msd = nl.analysis.MSD(data, TA=False, recalculate=True)\n", "\n", "plt.loglog(e_msd, label='eMSD')\n", "plt.loglog(msd_data, label='TA-MSD') # this was calculated above\n", "\n", "plt.ylim([1, None])\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition to ``TA=False``, note the use of ``recalculate=True``; with this switch ``MSD()`` will ignore cached values and just run the calculation in any case. Remember to use this when calculating MSDs with different settings for ``TA``!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parallelization\n", "If the dataset in question is very large, calculating MSDs can take some time; we might want to run those calculations on multiple cores in parallel. ``noctiluca`` has the following mechanism for this:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# For this example, we only use n=1 processes, so not really parallelized\n", "# Also note that chunksize can be an important parameter to actually\n", "# achieve a speedup!\n", "with nl.Parallelize(n=1, chunksize=10):\n", " msd = nl.analysis.MSD(data, recalculate=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [nl.Parallelize](../noctiluca.rst#noctiluca.parallel.Parallelize) context manager provides a simple way to parallelize computations in certain ``noctiluca`` functions, without having to change anything else about your code; parallelization is taken care of under the hood.\n", "\n", "Functions wrapped by this context manager need to be \"aware\" of this mechanism; note the corresponding statement in the [analysis.MSD()](../noctiluca.analysis.rst#noctiluca.analysis.p2.MSD) docstring.\n", "\n", "The function ``TaggedSet.apply()``—demonstrated in the ``TaggedSet`` tutorial—is another parallel-aware function in ``noctiluca``." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Other two-point functions\n", "Using ``nl.analysis.MSD()`` in the examples above is actually a short-cut: under the hood, ``MSD()`` lives in ``nl.analysis.p2`` and is in fact just a simple call to ``nl.analysis.p2.P2()`` (which then in turn dispatches to ``p2.P2dataset()`` or ``p2.P2traj()``). What's all this?\n", "\n", "The MSD is just one of various \"two-point functions\" that may be of interest for SPT analysis. Others include the autocovariance\n", "\\begin{equation}\n", "\\text{ACov}(\\Delta t) := \\left\\langle x(t+\\Delta t)x(t) \\right\\rangle_t\\,,\n", "\\end{equation}\n", "the velocity autocorrelation\n", "\\begin{equation}\n", "\\text{VACov}(\\Delta t) := \\left\\langle v(t+\\Delta t)v(t) \\right\\rangle_t \\qquad v(t) \\equiv x(t+\\tau_\\text{frame}) - x(t) \\,,\n", "\\end{equation}\n", "\n", "or other user-defined versions of these. ``noctiluca.analysis`` provides [MSD()](../noctiluca.analysis.rst#noctiluca.analysis.p2.MSD), [ACov()](../noctiluca.analysis.rst#noctiluca.analysis.p2.ACov), [VACov()](../noctiluca.analysis.rst#noctiluca.analysis.p2.VACov), [ACorr()](../noctiluca.analysis.rst#noctiluca.analysis.p2.ACorr) (auto*correlation*, i.e. normalized), and [VACorr()](../noctiluca.analysis.rst#noctiluca.analysis.p2.VACorr) (velocity auto*correlation*), all of which work pretty much exactly like demonstrated above for ``MSD()``. As an example, here are two ways to calculate the velocity autocorrelation for our diffusive dataset (expected to be delta-correlated):" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAACMCAYAAAAEJtCrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0RElEQVR4nO2deVxUVf/HP+feGfZFEWRVUGGAAQVFQVLDfcMlRUtzL9fKUtOkR7Is/EWZpWSmVoLbk/ZYueGSlqHG4zKgoCigJG6AAiIMO3Pv+f3BjE3KMiCGM895v17zYu4959z7Pefe+XCWe79fQikFg8FgMACupQ1gMBiMZwUmiAwGg6GGCSKDwWCoYYLIYDAYapggMhgMhhomiAwGg6FG0tIG1AYhJA/AjZa2g8FgGByulFK7uhIbFERCyGYAIwDco5T61pJOAKwFMBxAGYDplNIkddpQdRoP4FtKaZSORt+glHbXMS+DwWDoBCFEUV+6LkPmWABD60kfBsBD/ZkN4Gv1iXkAX6nT5QAmEkLkOpyPwWAwWoQGBZFSegLA/XqyjAawldZwGkArQogjgEAA1yilf1JKqwDsVOdlMBiMZ5LmmEN0BnBLa/u2el9t+4Oa4Xx/wy08bgqAaQAGZUWFsvcQGQxGk2kOQSS17KP17K/9IITMRs2QGwBsdT35F9KvJo/iEgasE8aYAKHl9eVNTExsK5FIvgXgC7bCzmAYKiKASyqVamZAQMC9xhRsDkG8DaCd1rYLgGwARnXsrxVK6SYAm4CGJz614UCreUJBQK0A1CuIEonkWwcHB287O7tCjuNYb5LBMEBEUSR5eXny3NzcbwGMakzZ5ugl7QMwldTQE0ARpTQHwDkAHoSQDoQQIwAT1HmbFQFcBQBwEC11yO5rZ2dXzMSQwTBcOI6jdnZ2RagZCTYKXR67+R5AXwC2hJDbAN4HIAUASukGAAdR88jNNdQ8djNDnaYihLwB4AhqHrvZTClNbayBDSGqBZGA6iKIHBNDBsPwUf/OG93h02WVeSKl1JFSKqWUulBKv6OUblCLIdSry69TSjtRSjtTShVaZQ9SSmXqtJWNNU4XBErUPURq8TSO39KYmZl1bUq5BQsWOO3Zs8cSAD788MO2SqXyqc6ZHjhwwPLo0aPmT/McjaGhdsvPz+ejoqIePqCblZUlHTp0aMenadPhw4ct3N3dfby8vOQlJSW1zbE3SEhIiHt+fj5fXx7ta99YDhw4YNmvXz/3hvIFBgZ6njhxwqy+PE/7vktISDDdtWuXdXMeU+8XFq5Q1+sxqiHIpE5CS9vyLLFmzZrsF154QQkAGzdutC8pKXmq1/q3336zPHny5FP7pySKIgRBqHO7sRQUFPDfffddW822m5tb9eHDh/98MivrZ+vWrTbz58/PTUtLu2xhYdGokYqmvvHx8ddsbW3rrbj2tW9JnvZ9p1AozOLi4pggahMrDE1boZqG3UJIaUvb0hDz5s1z1u6VLFq0yOn999+3B4D33nvP3tfX11smk8kXLlzo9GhZURQxZ84cFw8PDx+ZTCb/5ptvWmvSIiIi7GUymdzT01P+2muvOQNAWFiYW0xMTOvIyMi29+7dk4aEhMiCgoJkX3zxhe2rr776cLFr9erVtjNnznR59HyTJk1q7+vr6+3u7u6jbY+zs3PnnJwcCQCcOHHCLDAw0DM9Pd1o69atdhs2bLD38vKSHz582CIjI8MoODhYJpPJ5MHBwbKrV68aAcCtW7ckgwYN6uTp6Sn39PSUa3qVH3zwgb2Hh4ePh4eHz4cfftgWANLT0406duzoM3ny5PY+Pj7yw4cPW2hvZ2ZmGjXUbkVFRVxwcLBMLpd7y2Qy+fbt21sBwNtvv+1y69YtYy8vL/mcOXNc0tPTjTw8PHwAoKysjIwbN85NJpPJvb295fv377cEgOjo6DaDBw/u1KdPHw9XV1ffuXPnPtZuALB3715Lb29vuUwmk48fP96tvLycfP7557ZxcXE2n376qdOoUaM6PFpGl/pnZmYaabf/kiVLHDt06ODz3HPPeYwcObLD8uXL7bWvveZ6LVy40ElT//Pnz5sAwPHjx826du3q5e3tLe/atatXcnKycW110VBSUkJGjBjRUSaTyUNDQztWVFQ87OHWdq88et/Vd09pExkZ2bZTp04+MplMPmLEiI4AUFxczI0fP97N19fX29vbW759+/ZWFRUV5OOPP3bav39/ay8vr7/9Hp4ISukz9wGg0DWv+9KfR8qW/kR9lv4nqKG8Fy5cyKKUKlrqc+rUqdTu3bsrNdsdO3Ysz8jISPnxxx8zJkyYkCcIgkKlUin69u374ODBg2mUUoWpqalAKVXExMRcCw4OLqqurlbcvHnzgoODQ2VWVlbyrl27Mvz9/UuKi4uTKKWK3Nzc85RSxdixY/M3b96cSSlVODk5VWZnZ1+glCqKioqSXFxcKioqKhIppQp/f/+SM2fOpD5qq+Y41dXVih49eihPnz6d+uix4uPjL/fo0UNJKVUsXLgw+7333rulKd+vX78H0dHR1ymlii+++OL6gAEDCimliuHDh99fsWLFTc2x8/Pzz584ceKyh4dHWVFRUdKDBw+SOnXqVH7q1KnUtLS0FEIIPXbs2BVKqeLRbV3araqqSlFQUJBEKVVkZ2dfaNeuXYUgCIq0tLQUd3f3co292tvLly+/FRYWlk8pVSQlJV1ycHCoLC0tTVy7du11Z2fnyvz8/POlpaWJjo6OlVevXk3WbrfS0tJEe3v7quTk5IuUUsULL7yQr6mv9jXR/uhaf+32j4+Pv+zp6VmmVCoT79+/n9S+ffsKTfs/eu0jIyNvUkoVH3/88Y0XX3wxj1KqKCgoSKqqqlJQShU///xz+uDBgwsppYr9+/en9+3b98GjNr7//vu3xo0bl08pVZw+fTqV53kaHx9/Wdd7pb582h87O7uqsrKyREqpIi8v7zylVPH666/nfPXVV39q9rm6ulYUFRUlrV279vqUKVPu1fV7U//eG6Utz6Rzh8YQIdkhnyb5BZtVQ7oB487oWm7J7uR2GbnKeudAGovMwbJs1Ti/W3Wl9+rVq7ygoECSlZUlzcnJkVhbWwseHh5Vq1atanvixAkruVwuB4CysjIuLS3NZNiwYSWasidPnrR88cUX70skErRr104VFBRUcurUKbPff//dcvLkyfmWlpYiANjb29c7nLKyshJ79eql3LVrl3Xnzp0rqqurSWBg4GOPK23ZssUmNjbWVqVSkby8PGlycrJJUFBQvY81aXP+/HnzQ4cOZQLAvHnz7q9YscIFABISEix37959HQAkEgnatGkj/P777xbDhw9/YGVlJQJAaGho4fHjxy3Hjx//wNHRsWrAgAEPe//a24cPH7ZqqN1EUSQLFixwOX36tAXHcbh3757R7du3673vExISLObPn38PALp27Vrh5ORUdfHiRRMA6N27d3GbNm0EAHB3d6/IzMw0dnd3r9aUTU5ONnFxcans0qVLJQBMnz694KuvvmoLoM7n4RpTf+0yw4YNe6AeetNBgwY9qOv4L7/8ciEABAYGlu3bt681ANy/f59/6aWXOmRlZZkQQmh1dXW9c5qnTp2yePPNN+8BQFBQULlMJivTpOl6r+iSz9PTs3zMmDEdRo0a9WDSpEkP1HW1OnLkSKvo6GgHAKisrCTXrl0zqs/epqL3gkhrVrbBgT4zE/r1MXLkyMLt27e3zs3NlYaFhd0HanrpCxYsyFmyZEl+XeVoHcHAKKWo8a+hO7Nnz85fuXKlg0wmq5g8efJj50xLSzNat26dfWJi4hU7OzshLCzMraKiggMAnuepKIoAgPLy8maZcqmrbgBgZmYm1rWtS7tt3LjRpqCgQHLx4sUrxsbG1NnZuXNDdtdnj5GR0cNEnucfE5L6yjblfI/WvynnMTExoQAgkUioSqUiALB06VLnkJAQ5dGjRzPT09ON+vfv79nQcWq7z+q7V5qS7/jx41cPHTpkuWfPnlaffvqp09WrVy9RSrF79+5rfn5+ldp5T5061ey/eb0XRBFcKQAQUNPGlKuvJ/c0mTJlyv1Zs2a5FRYWSuLj49MBYNiwYcUffPCB0+zZs+9bW1uL169flxoZGVFnZ2eVplxISIjym2++sXvjjTcK7t27Jzl79qxFdHT0LWNjY7py5UqnWbNm3be0tBTv3r3LP9pLNDc3F4qKijhHR0cAQP/+/UvfeOMNo9TUVPOLFy8+9ihUYWEhb2pqKtrY2Ai3bt2S/P7779YhISFKAHBxcan6448/zF588cXiH3744eG8jaWlpVBcXPxw9bNr166l3377bevXX3/9/saNG226d+9eAgC9evVSrlq1ym758uX3VCoViouLuf79+5e88sorbh999FEupRQHDx5sHRsb2+AChy7tVlRUxNva2lYbGxvT/fv3W2ZnZxsBgLW1tVBaWlqrMPbu3btk+/btNqNGjVKmpKQY5+TkGHXp0qXizJkzDY4o/P39K+7cuWN06dIlY19f38qtW7e26dOnT70LHE2pf9++fUvmzZvnWlZWllNdXU2OHTvWaurUqXkN2aehuLiYd3FxqQKAjRs3NvhmmKZNRo4cqTx37pxJRkaGGVD/vaJ939WXT4MgCMjMzDQaOXKkcvDgwSVOTk42RUVFfL9+/YpXr15tHxsbe5PjOPzxxx+mvXr1KreyshKae9HGAASRlAIAD9qsw9+nRffu3StKS0s5e3v7KldX12oAGDt2bHFqaqpJjx49vICaXsGOHTuua/+wp0yZ8iAhIcHC29vbhxBCV6xYcbt9+/aq9u3bFyclJZn5+/t7S6VSOnDgwKJ169bd0T7ntGnT8ocNG+bRtm3b6jNnzmQAwAsvvFCYkpJiZmdn99gQOzg4uNzX17fMw8PDp3379pUBAQEPh6DLly/Pnjt3rtsnn3xSHRAQ8HAoFxYW9mDcuHGdDh061GrNmjU3v/7665vTpk1zW7t2rUObNm1UW7duzQKAr7/++ub06dNdZTKZLcdxWLdu3Y2BAweWvvzyywXdunXzVtc1r1evXuXp6en1Dot0abeZM2feHzZsmLuvr6+3j49PWYcOHSoAwMHBQQgICCjx8PDw6d+/f9GiRYseDmnfeeede1OmTHGVyWRynuexcePGLFNTU526ZGZmZnTDhg1Z48eP7yQIAvz8/MoWL15cr1D17t27rLH1DwkJKRs6dGiRXC73cXZ2ruzSpUuptbW1zsvuS5cuzZ05c2aH6Ohohz59+hQ3lH/x4sX3JkyY0EEmk8l9fHzKOnfuXArUf688et/VlU+DSqUiL7/8cgelUslTSsmcOXPu2traClFRUdmzZ89u7+XlJaeUEhcXl8rjx49fGzZsmPKzzz5z9PLykr/99ts5s2bNKtS1/nVBmtLFf9oQQhRUR3+I30RMHDFLcnD/v1X9o1+O/Pmt+vImJydn+fn51Tm8+l+iX79+7gsWLLg7evToFn88g9E0ioqKOGtra1GpVHLBwcGeGzZsuNG7d++yhkv+b5CcnGzr5+fnpr2vIW3R+8duMqhL7jrVaCRRD52HC//L5Ofn825ubr4mJiYiE0P9ZvLkya5eXl7yLl26eI8cObKQieGTo/dD5v8IfXPVX+991qKW6Ae2trZCVlbWpZa2g/Hk7N+//3pL22Bo6L0gtkFRFQWBFCqDfHWPwWD8c+j9kHm65Ig0yWQu3pX+u29L28JgMPQbvRdEEaQYAHjQel89YjAYjIbQe0GkICUAwEFkgshgMJ4IvRfEBZExtIry4CEYpCBu27atVWJioklL2wHUOBvQOECoL8+GDRtsNNsnTpwwmz59erv6yjwpmzdvbt2xY0cfjRMBDf+ES6+nRX2urZrq9upJ2iM6OrpNVlaWtLHlnsQVWUugkyASQoYSQtIJIdcIIeG1pC8hhFxQfy4RQgRCiI06LYsQclGdpnNogMZQDSk40KfybmNLs2fPnlYpKSmNegunMahUqnq3G8vVq1eNd+3a9VAQn3/++bLY2Nin+lZQTEyM7dq1a29qHjrX0FwuvaqrqxvO1MzU59qqvrT6bH2S9ti+fbvtzZs3GyWIKpXqmXFFpisNCqIu8ZUppasopf6UUn8A7wKIp5Rqhy7tp05/KsHn16rGlv4idr/9NI7d3AwcOLCTj4+Pt7u7u89nn3328JUpbYemMTExrcPCwtyOHj1qfuzYsVYREREuXl5e8tTUVOOEhARTPz8/L5lMJh80aFCnvLw8HgAuXbpk/Nxzz8k8PT3lcrncOzU11bgul2EHDhywDAoKko0cObKDp6enz6PbKpUKc+bMcdG41Vq1atVjr3alp6cbBQQEeMrlcm+5XO6tceO1bNkyZ4VCYeHl5SVfsWJFW22Ho3fv3uUHDhzYSSaTyf38/LzOnDljCtS4QRs/frxbYGCgp4uLS+fIyMi2j54PqHkvWSaTyT08PHzmzZvnDACLFy92TExMtJg/f77rnDlz/uaOS7tHW5/rrt27d1vJ5XJvT09PeXBwsExj08SJE1179erlMXbs2A7Z2dmSIUOGdPL19fX29fX1/uWXX8w1+caOHevWq1cvD2dn585btmxpNXfuXBeZTCbv06ePR2VlJQGAkydPmvXo0cPTx8fHu3fv3h43btyQAjWOVufNm+fcuXNnbzc3N9/Dhw9b1Ofaqra0R22t69pot0d91/hRd3IxMTGtL126ZDZ16tSOGue2tbk4A2rcjS1evNgxICDAc/Pmza21XZHV1Qa1ufxqMepzhaN+iyUYwBGt7XcBvFtP/n8DmKW1nQXAtqHzNMZFz6Mf16UH7rguPfBtQ/la2v0X1XKBpFQqE93d3ctzcnLOUy13VZRSxebNmzPHjh2bT2txGeXh4VF24MCBNEqp4q233sqeMWPGXUqponPnziVbtmy5RtUuqIqLi5Pqchm2f//+dBMTE+HKlSspVO3ySXt71apVWUuWLLlDKVWUlZUl+vj4lF65ciVF20VWcXFxUmlpaSKlVJGSknLRx8enVHMsbfdR2ttTp069u2jRojuUUsXevXvTPT09y6jadZi/v39JWVlZYnZ29gVra2uVxj2Z5nP9+vVkBweHyjt37lyoqqpSBAUFFW/duvUapVTRo0cPpcYVlfZH2966XHfduXPngr29fZWm7prrs3Dhwmy5XF6qVCoTKaWKESNGFBw+fDiNUqrIyMhI6dChQ7kmX7du3ZQVFRWJCQkJqSYmJsKuXbsyKKWKgQMHFm7duvVaRUVFor+/f8mdO3cuUEoVmzZtytS40urRo4dy5syZuZRSxc6dO68GBwcXa+yty7XVo2mP2lrXtdFuj7qucV3u5LTbuD4XZ05OTpXLli176AZOc//W1wa1ufxqjs/Tcv+lc3xlQogZgKEA3tDWXAC/EEIogI20Jrpes+JC7qkkEFo1uuCmfo97+PAeeR99FuWhsoTDlpEej6X7vZSPoLkFUOZK8P3ETn9Lm308vaFTfvLJJ/ZxcXGtACA3N1eamppq4uDgoJNz24KCAl6pVPKhoaElADBr1qyC8ePHdywsLOTu3r1rNHXq1AdAzfu0AGhdLsOsra3FLl26lHp5eVVpjq29fezYMau0tDQzjasopVLJX7582cTHx6dCk7+qqoq8+uqrrpcvXzblOA43btxocA737Nmzlj/++OM1ABg1apRy9uzZkoKCAh4ABg8e/MDU1JSampqqbGxsqm/fvi3p1KnTw/HfqVOnzHv27Kl0cnJSAcBLL710Pz4+3mLKlCkPdGk7oHbXXQUFBXxgYKBSU3dtxxhDhw7VuNfCH3/8YXX16tWHUxclJSV8YWEhBwADBw4sMjY2poGBgeWCIJBx48YVA4CPj0/59evXjVJSUoyvXr1q2r9/fxlQ4+zXzs7uYd3Gjx9fCADPPfdc6ZIlS5o09aNtqy7Xpq5rfPToUauG3Mk15OJs6tSpj71TXF8b1Obyq6XQRRAbE195JIA/6N+Hy70opdmEkLYAjhJC0iilJx47SRPjMgPADun/OdyhtoHAzMYU+8c5cOCAZXx8vKVCoUiztLQUAwMDPTWuqLRdK2mGH7pC63ENVhcNuNUiq1evvhkWFva3l/61nQ2sXLnSvm3bttU//vjjdVEUYWpqGtAUO9X/KGFsbKztVgsaN1W61EVXanPdRetxn2Zubv43V2MKheJKba7/NbbzPA+JREI5rmYmiuM4qFQqQikl7u7u5RcuXEir7Txa7rkgCEKTYq1o26rLtanrGh88eNCqIXdyDV0LjZg+er662qA2l19SaaPXb5oFXRZV6oq7XBsTAHyvvYNSmq3+ew/AzwACaytIKd1EKe1Oa+YZG+WAQQVelBCh8S04+3j6Y58+i2reiTa2EGtND5pbAACwdFA9ltYADx484K2trQVLS0vx/PnzJsnJyQ/9ubVp06Y6KSnJRBAE7N279+GckYWFhVBcXMyp8whWVlbC4cOHLQDgu+++axMcHFxiY2MjOjg4VG3btq0VUCOoSqWSCwkJUe7evdtGpVIhOztbcvbsWYs+ffo02BsdNGhQ0ddff22nmf9KSUkx1tigoaioiHd0dKzmeR7r169vo4lvYm1tLZSUlNQaBKlnz57KmJiYNkDNP4fWrVurbGxsavX39yjPP/986ZkzZyxzcnIkKpUK//nPf2z69u37mMeUxtKvX7/SM2fOWKalpRkBNfOcteXr3bt38SeffPJwbjMhIUHnha4uXbpU3L9/X3Ls2DFzoMbBqUKhqPfJgfpcWzXk9qqua6NNXdd46NChxdu2bbPVBIfStIeFhYVQVFTEA393cQYAurg4q6sNtF1+rV+//rZSqeQ152kJdBFEneIrE0KsAYQA2Ku1z5wQYqn5DmAwgGZ/j1YFXpBAeOZfQwwLCytSqVREJpPJ//Wvfzn5+fk9FKcVK1bcGT16tHtwcLCnvb39w+HUpEmT7kdHRzt4e3vLU1NTjWNiYq4vXbrURSaTyVNSUkyjoqKyAWD79u3Xv/rqq7YymUzevXt3r1u3bkmmTJnywMfHp9zb29unb9++Mo3LsIbsXLhwYb6Xl1dF586dvT08PHxmzZrl+qgj1AULFtz7/vvv2/j5+XllZGSYmJqaigAQGBhYLpFIqKenp3zFihV/Wxz55JNPspOSksxkMpl82bJlzrGxsTq/i+vq6lq9fPnyOyEhITJvb2+fLl26lE2ePPmBruXrwsnJSRUdHZ01ZswYd09PT/mYMWNqndTftGnTraSkJHOZTCbv1KmTz7p16+xqy1cbJiYmdOfOnZnh4eEunp6ech8fH3l8fHy9r5oOGzZMmZGRYVpbvJD60oC6rw3wV4+8rms8bty44mHDhj3w9/f39vLykn/00UcOADB16tT8+fPnu3p5eclFUYTGxZlMJpNzHIeGXJzV1QYal18ymUzu6+sr17j80rVtmxud3H8RQoYDWIO/4iuvJITMBR7GZgYhZDqAoZTSCVrlOqKmVwjUDM//TXUIR9oY918AcGl5lyIBvOj34fl6A80w91+M/2VOnjxptmjRonbnzp1rcDRjCDTF/ZdOvSpK6UHUBKTX3rfhke1YALGP7PsTgJ8u53gSVOBVUjRhyMxg/I9w4sQJs6lTp3ZcsWKFXjye1lI888NMXdguDPpTpFzrz1vaEAbjGeX5558vY27fGsYgBHG3EHILgAkTRAaD8SQYhCB6kFu8EVS6vC8piqJIOI579uImMBiMZkMURQJApycYtNF75w4A8J5ku8+3Rp+5NJwTl/Ly8qzVjcVgMAwQURRJXl6eNZrwRItB9BAFcFVSNPxAq0qlmpmbm/ttbm6uLwzknwGDwXgMEcAllUrV6Dc1DEIQRXDVUqgaFMSAgIB7AEb9AyYxGAw9xCB6STU9xIYFkcFgMOrDUASx0ghP5sePwWAwDEIQ9wi9UxZVz8OEf61qsXcgGQyG/mMQc4hHxB431F+NAbBg3QwGo0kYhCB2JpnG5qQS1ZQ3B0KZIDIYjCZhEEPm6ZIj3juNIhHCJ9u3tC0MBkN/MQhBFChfAQA8qN5E92IwGM8eBiGIIkg5AHAQ6/Uxx2AwGPVhIILIlQMAD9G8obwMBoNRF80Vl7kvIaRIKzbzcl3LNgcaQeSYIDIYjCegwVVmrbjMg1ATX+UcIWQfpfTyI1lPUkpHNLHsE3FY7H7p9yo/3KeWec92mCkGg/Eso8tjN4EArqm9X4MQshPAaAC6iNqTlNWZU2KXXPXX6nozMhgMRj3oMmSuLS6zcy35ggkhyYSQQ4QQn0aWfSICSLq0P5eE3txFh+Y+NoPB+N+hueIyJwFwpZSWqANS7QHgoWPZmpM8QVzmYfzZtjMlh7BdNcC3MeUYDAZDm2aJy0wpLaaUlqi/HwQgJYTY6lJW6xhNjsssgisFAAKqc6xcBoPBeJRmictMCHEghBD190D1cQt0KdsciCAlAMBDrDf4N4PBYNRHg0NmSqmKEPIGgCP4Ky5z6iNxmccBmEcIUQEoBzCB1gR8rrVsc1dCI4gc6yEyGIwnoFniMlNK1wFYp2vZ5oaCUwIAT0QmiAwGo8kYhLebeKHL3VOiL8qoSVpYSxvDYDD0FoMQxLPUu0y9dl3RwqYwGAw9xiAEEUDlSC4BFTBqD4S2tC0MBkNPMQhBzIoKFarfn4KjYoBfS9vCYDD0F4PwdgMAVZBAAsGope1gMBj6i8EIogo85SEyQWQwGE3GYASxGhLKQZS2tB0MBkN/MShBZD1EBoPxJBjEogoAzK5adLsS0su/tLQhDAZDbzEYQbxIOyoBCC1tB4PB0F8MRhDHcCdNOVAX9hwig8FoKgYjiLMkB9tWgW8NfNLSpjAYDD3FYBZVVOBVEggGI/AMBuOfx4AEkRMkEPiWtoPBYOgvBiSIvErKBJHBYDwBzRWXeRIhJEX9SSCE+GmlZRFCLqrjNSua03htVJSvlkAwGIFnMBj/PM0Vl/k6gBBKaSEhZBiATQCCtNL7UUobFSelsUSoXjmrAh944mmehMFgGDTNEpeZUpqglf80aoJJ/aP8SZ0ewIBWzRkMxj9Pc8Zl1vAqgENa2xTAL4SQRHWo0afCWO6E7Uw+zuppHZ/BYBg+zRWXuSYjIf1QI4i9tXb3opRmE0LaAjhKCEmjlD42sn2SuMwAMJQ/164Xl2oOrG9sUQaDwQDQTHGZAYAQ0gXAtwBGU0oLNPsppdnqv/cA/IyaIfhjPElcZgAQwVVKUd3YYgwGg/GQ5orL3B7ATwCmUEoztPabE0IsNd8BDAZwqbmM10YAV2VEBKyJmMEevWEwGE2iueIyLwfQBsB6dbx6lbqnZw/gZ/U+CYB/U0oPP42KCOAqAYCAWgAoehrnYDAYhk1zxWWeCWBmLeX+BPCPxDkRwFUAAAdqBSaIDAajCRjMg8xR1RPju1d8jRjV0KqWtoXBYOgnBvPc3l3YKNVfmddsBoPRJAxGEMO4+LZu3F3cpG3tgdBbDZdgMBiMv2MwQ+YA7qrdfMkeeJA7ji1tC4PB0E8MRhBFkHIAIKDmLW0Lg8HQTwxGEAVwZQDAQWSCyGAwmoTBCCJV9xA51kNkMBhNxGAEUQQpBQACatbStjAYDP3EYFaZ16rCzq9SvQQV+MvTW9oYBoOhlxiMIBbBokz9VdqihjAYDL3FYATxBe6UaWfuOm7Stt4sNjODwWgKBjOH6Mndkr4qOQQfktWhpW1hMBj6icEIogp8CQBwRDRpaVsYDIZ+YjCCKIJTAgAPygSRwWA0CYOZQ1TRGkHkwHqIDMazzpqIGcQIqoAqSBIXRMZQANgQMWmgA7n/1i3a9rX5kVtaxB9Bc8VlJoSQaHV6CiGkm65lm4tyGBfXfKPM280zxNqIGaYtbQPj2aMdubfrNcm+c2aoyHMLj/u/1cteVczm446+wCeM6EquxbeUXc0Vl3kYAA/1JwjA1wCCdCzbLFiQ8kq3ih0UIKdHN+Nx3cLjuKyoULEZD/k/w/qIyT0n8fF/nFkelHta9Ja/FRn7xI5710dM7ulJbm9Ko+0nvB65tdnvI31hTcQMnzWqcb4AugPYnxUV+jBwm1t4HAHQJ4BkLNls9GnoAaHn6kmRPy95WrZsiJjUrxPJ/vwqdZn5WuS2RM3+4PCtvtMlh20rYJT4VmSsxj0f1kZMt5zIXwzLFB0rfhW6XQPwzl6xV5lcvJHAQzQezCcG7IgYGzUp8qfHOlCRy14zaUseDC+HUeAXqvG/ZUWF/tKcdSGU1hpA768MhAQD+IBSOkS9/S4AUEo/1sqzEcDvlNLv1dvpAPoCcGuobB3nVKhDEDQKt/C4cgBfZkWFvqNjfgLAGoBXlOSbOaakwjGHtvl8buSOX95f9qb9FmHIEgDzX+EPFg/kksry0OpSMTVLLIXJiQoYHV8QGSM0cHyZBCrza1Gjz2vtkwIIeJ3f09OF5LmXwfi8EmZbF0TGPIyQ5R++y3KBZPdMc1R0u4fWH2v/8NdEzDC2QPm4NkQ51hSVjpnUcc7rkdsuap93TcQMaQeSu8cCZe0v0Q69HxWi7yJeerUnd+Xz29Tu9pzqRd2yokIrG2qrDRGTBvpzmTG5tPX5LOowWjPMqYu1ETPMB3Pnct3IXQtTUoXfBP+7r1S/45EVFaqsrxwAfBkxzVkAl/to+66NmG4Zyp3JdeeyzcZWfrDvp48XjgYetulCAL9nRYWe1eSPjphuBwBvRsbmNXTOughf9rbjTqH/UACFAPb9E/8c10TMIDzE9vMjt9x4NG1jxMuh3bmMmADuql3Pii+RCxuskmykIrhvXorcO+fbiAnTbtG2H2wRhrgBtOiC8WxLE1RxO4SBM16N3BmrfayvIqb4tULpxDxYv9vQ9awL7/Cf7HcafZTlx/1pkk+txH3Cc/M+VE3dDyDag9wed9S45qd4RWxfckzs5js/csuN3e+N3D6OPzEpVjXknemRP6xyC49rBaAqKyq0bG3EDNPh3JkCc1JuPKvqbZe4j9/I0ZzLLTxuoivJjY03XmR0n1ogoHLjT9ejRoQ1xt6GtEUXQRwHYKg6TAAIIVMABFFK39DKcwBAFKX0lHr7VwBLUSOI9ZZtitF1sTViXGURzP8rgBvAgVpJoVrxiWpiTwDW0/gjf8rJjVvmpNzdiRR0diIFrdPEdmR6dTgHAL8ZLUJHLhcAcEl0K3Ig9y0XV88lv4v+378n2dp5OH9W7kjuPwxgdZvaCv0qV8dWQ3phk3T16zJyu2M+rAuTxU7/+lnotfUS7bjEFBUfbJZ+Jv2TOq6eFPnTkon/+tT+v6LPjwB6rZNGYwR/GgBQSC3oNep8N1N0VIWrZqsAuJ4yfpO4kHzco9ZinNDztRWqad8AmOKI/FX/NXnTDgBESqCEKT0kBH2di9ZvLIiMoesjJvfsw1083JnLsgaAC2LH+3uFXp7vr1yXvzFiUrAXuRkbwqfIAGB21UL8IvY4DmDsa/xe0YHcf8MEVe5VkOT8SR1TNwvDM01Qhc+l66MGcOf7chAhJQJ+Fbpe+kYI9e/FXXJyIgWf2eFBNzNSaWWCSpPr1PHMd6phM+ZKDuwbyp/rtkPV/zMAZluEIXMyaLvLw7nTv/bg0jMrIc3gIdoQUL4Epjs14hcb8eKCF/njX+RSm8r/ivJ1eWi1VJP263t9kwfw57ssrZ51cZfQrzOAkKyo0BNx7w1O3C8EdzssBmIglxj3PJeS0IHkvNqDS+9oQqqxXjWq/FPVhMsmqLy4UrrZjIdg6kzyA93IXbsyaqz6VDXhUpzY88+BXGK7RZLd8kpIqkqpabGECNIuJNOlf+Vq5KINZvCHbnuQO7urIL0hgWBLQKXZtM2fXwkvnOcgtP1MunE6D9E0m7b5eF7k9lMA4BYe5wKgbwBJNx7J/9daAqFNK1IaZINiDwLKrVBNO5tG298OIOlXB/OK6zZQvtSVuzbencs2OSt63p9Q9d5cN5KbNpM/+IojKRjSh7voXQZjnBC7/PaT0OcTS5RfniE5lOzP/WmTJLrn+5NM2zvUVhhW9fHrJTDbtkSyy20MfzLZCCruByGk12uR208DwJaI8UtH8v/92IYoyUnBNzOJyvzeiowp1dzfX0ZMc/Mgt7/nQCVXqGtfTdpP74X+R0ZuD7osuv70q9ht+RGxx8/tyF3fefz+2D5cyqvtuHzpUaFb9azqxaIryV3zGr/X1pKUywZyiX1SaYeCY0I390mSY3lF1KLiiNjdqjYhXhcxddTPQu/dmdQ5DkBYVlSo6B++0+kBLC9zEDMiJd+dKYHZ0f8IIYeOfjynUaE2G9KW5orLXFeexsR0fqK4zAAwmFfwDqQwpIryKhV4mJEq7BeCr12mbulDuHNDn+Mv8wIlyKIOFVdF5xt3qF0WgEMArm4WhmU7igVV7UhelB/J7J9LW5cEcVdei/2/ZTs0D3p/GTHN1QyVAyxR1qccxgHVkL4I4FUeYlUBrIo6key23SUZ3w3lz256tWoJn0et91iTkgETuN8W73lvmOf/Sf4cEq0aw/0s9nntD9En5w5tU22FssEuJG9Ee3LP5QGxKAQQDyAmVjUk34nk08F84tpp/C8b7Enhx69VL2idA9vTMaohW8thvJOHYPc8d/GHMfzJ1wZVrXpxTXicarP0VtuOJJfbpeq7Sd0ss7YIQw62iphuNp+P8xHA4YjQ/Ww6dXnhF7HHAACb35b8cOMV/pCVOfmro7hdNQAAUA0ew/mzSBDkNy9Q97EdSfZXSaJH0GlRfm4Ul9A5jD8puUetxSJqUSmAiO7kzqAL1P36atV4aQG1PKkZqi0LjzsI4It+3IUF4yV/D8t9Qez4ZXD4lu4zJQeHTOB//yKPtqriIZJx/Im3+1V+PmZNeNzvn0g2ur8kOd/lFyEgaZfQr48tHqTNk+zfvfu9r38dx5/pxkM8f1QMOPIvyY6lHbncUCU1RYLok1oC09sZokslADN3cmd0GH+yNQDcEm2r08R2WVIiGN2jraQAfMtgLFaDV5miysSBK2zFQyD/FX0u9uEvvvuj8Lz1eD7+Ozl3c4G27UeE7oAAiODxPJcCW1IMkZLh55d3K5BCJQ0iU63OUG84kvuYLqkZ3YmUIAc2KhXlxWrwPQAMCeKumM+RxAEArolOFYeF7mcyqZOnCO6HW7QtxvE102rxYpdLV6hr2BuRWzNGqG1YG5HTXknNUnpxqR2Pi/4pl6nriEtR49WLEqGXN0W8PHYi/9u+yfyvCYHh238vh/Hd34zPTCii5pXnRM+UIbyih5lYebNT+L5QL3LTca7kwJypfPJgS5QTjlC0EYtvBIVvk62Qxh4cyyuCCqkF9ZVkzfAUb804WhUg3qL2Y16O/HnflxHTPgnClQRbUmzkT6713/PxWyk1YdqBHyJGxbwoiZ9+VXTOOCL0kABYVVev9I3Irfs+C497B6BfTOSP/7QmYveY9dK0M5epq0WkasrklyP3ZAB/iUVzYlBD5m8iJoa1Jsowc1S4GkFleZO2/bYI5l8uiIyhayJmGJuhckgFpOfejNyS0/DRGsYtPM4IgCOAm1lRoTQ6YpqjK7n7XSeS89xJsfP6eZHb/xUdMd3ueS45zZ/706aYmtEfhT4LZkT+EK3rOb6MmOYayKUldCTZDi9XRbx5lbp8rT1sWxMxQ3qXtl75vTCgFQAsluzqIIHw8dzIHb+pbRwDYKcPuV70umRv2g1qHz4vcnuCpvzod9eO/kgasyWfWt++Q213l8D0VylU7dNoe8luISRPimpuqWSn9czInTu06j0fwPudSeZvo/mE/TMjd27TpI1594uu56lsBoBWAGZmRYX+LcbNuoipMmuUviyBYCuAKzIjlZ451GbkKtWEymjpl6b+5BrdLwYHVEOSVkCtlm4TBvcGEPCD0QpbC1Qoj4oB9m9FxpRviwiLnCI5tgwAkkT3/BNiF6cFkTHVGyIm9bNE+dACWH326FBZs7JJgerXI7cl63oNtMobt0LJ6wCoCvw9CqK6Te0kW4QhRQDyZvJx+Tak2K49ufd/XblrfSqpVPhe6P/DN8KI1UO4s9Y9uSv+IkhpKUwPvBkZe1erPfm5/L7BHUnOuBKYKIphvmFBZAx1C4/jAQwBYLdQ8h+lCaqPzYncUVyHbUQCocMbkVv/rC39u4gJ09uTu0tnV79dTkE8e3MXfwnk0ma+GRlbsDNi9AYF9Zy5WwjhO5E7+NV4CZJE9/xEUTbXCmVDTover+wXg0s2S1dZq8CnptCOAa1Q8poHufPWb6L/lvdWrn9fl/bbEjHu+EbVyL7ZsE0G0DUrKrRO8XELjyOv8gd/fk+6ffQ50TO3B5fusF/o+cvIj44M0eVcddEcQ2YJgAwAAwDcQU2c5pcppalaeUIBvAFgOGoWVaIppYG6lG2K0frG2ojp1h1I7q4cavPlnMh/xzXlGGsiZpCmzvO4hcfZAXiQFRXaqOHFP4VbeJwbgG+NUO00mz8wcfHKb5MfSSeLJbsCRHCZb0bGFgI17RFErty0J4V2cWJPr/mRW7JawHS9xC08jjwqRm7hcXIAHR1QUDJdckQyN3LHMa20fgBirVC6+RXJoQ+f4D7kAUQAiMuKClU0lH9NxAw+gGSk9+EvdcoUHSsOiME2b0XGlDfl3BqeWBDVBxkOYA3+isu8UjsuM6kJvLwOwFAAZQBmUEoVdZV9UqMZhkltP9T6WBMxgwfAaS9IMQyLtREzTN1Izo+5tE30nMgdTxzTvVkE8Z+GCSKDwXgaNKQtBvPqHoPBYDwpTBAZDAZDDRNEBoPBUPOsziHmAXjsKf16sAWQ/5TMaQlYfZ5tDKk+hlQXoOH6uFJK7epKfCYFsbEY2iIMq8+zjSHVx5DqAjx5fdiQmcFgMNQwQWQwGAw1hiKImxrOolew+jzbGFJ9DKkuwBPWxyDmEBkMBqM5MJQeIoPBYDwxei+I/1SIgqcBIaQdIeQ4IeQKISSVEPKWer8NIeQoIeSq+m/rlra1MRBCeELIebWfTL2uDyGkFSFkNyEkTX2dgvW8PgvV99olQsj3hBATfaoPIWQzIeQeIeSS1r467SeEvKvWhnRCSIOecvRaELVCFAwDIAcwkRAib1mrGoUKwNuUUm8APQG8rrY/HMCvlFIPAL+qt/WJtwBc0drW5/qsBXCYUuoFwA819dLL+hBCnAG8CaA7pdQXNQ5XJkC/6hOLGicy2tRqv/q3NAGAj7rMerVm1A2lVG8/AIIBHNHafhfAuy1t1xPUZy9q4s+kA3BU73MEkN7StjWiDi7qm7I/gAPqfXpZHwBWAK5DPdeutV9f6+MM4BYAG9Q4hz4AYLC+1Qc1flYvNXQ9HtUDAEcABNd3bL3uIeKvC6zhtnqf3kEIcQPQFcAZAPaU0hwAUP9t24KmNZY1AN4BoB17RF/r0xFAHoAY9RTAt4QQc+hpfSildwB8BuAmgBwARZTSX6Cn9dGiLvsbrQ/6Log6hyh4liGEWAD4EcACSmmtHpH1AULICAD3KKWJDWbWDyQAugH4mlLaFUApnu3hZL2o59ZGA+gAwAmAOSFkcsta9VRptD7ouyDeBtBOa9sFQHYL2dIkCCFS1IjhDkrpT+rddwkhjup0RwD3Wsq+RtILwChCSBaAnQD6E0K2Q3/rcxvAbUrpGfX2btQIpL7WZyCA65TSPEppNYCfADwH/a2Phrrsb7Q+6LsgngPgQQjpQAgxQs0E6r4Wtkln1J7GvwNwhVL6uVbSPgDT1N+noWZu8ZmHUvoupdSFUuqGmmvxG6V0MvS3PrkAbhFCPNW7BgC4DD2tD2qGyj0JIWbqe28AahaJ9LU+Guqyfx+ACYQQY0JIB9TEjT9bS/m/aOkJ0maYYB2OmrgtmQCWtbQ9jbS9N2q68CkALqg/wwG0Qc3CxFX1X5uWtrUJdeuLvxZV9LY+APwBKNTXaA+A1npenxUA0gBcArANgLE+1QfA96iZ/6xGTQ/w1frsB7BMrQ3pAIY1dHz2pgqDwWCo0fchM4PBYDQbTBAZDAZDDRNEBoPBUMMEkcFgMNQwQWQwGAw1TBAZDAZDDRNEBoPBUMMEkcFgMNT8P+SZpcmV2+eyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "vacf0 = nl.analysis.VACorr(data)\n", "vacf1 = nl.analysis.ACorr(data.apply(lambda traj: traj.diff()))\n", "\n", "plt.figure(figsize=[5, 2])\n", "plt.plot(vacf0, label='velocity autocorrelation of original data set')\n", "plt.plot(vacf1, linestyle='--',\n", " label='autocorrelation of increment trajectories',\n", " )\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All of the functions mentioned above are implemented as special cases of a generic \"two-point function\" (the \"two-point\" referring to the fact that these are averages including two different time points along the trajectory):\n", "\\begin{equation}\n", "\\text{P2}(\\Delta t) = \\left\\langle f\\left(x(t+\\Delta t), x(t)\\right) \\right\\rangle_t\\,,\n", "\\end{equation}\n", "with different functions $f$ (and in the case of correlation functions also some normalization on each trajectory before taking the ensemble average).\n", "\n", "You can define your own two-point function and take advantage of the full functionality implemented in ``P2()`` (such as parallelization of ensemble averages, time averaging or no, caching, etc.) by just specifying the function $f$. For example, here is the full definition of ``analysis.MSD()``:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def SD(x, y):\n", " return np.sum((x-y)**2, axis=-1)\n", "\n", "def MSD(*args, **kwargs):\n", " return nl.analysis.p2.P2(*args, **kwargs, function=SD, writeto='MSD')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The last keyword ``writeto`` specifies which ``meta`` entry of the trajectories to use for caching the calculation.\n", "\n", "For more details, read the documentation of [analysis.p2.P2traj()](../noctiluca.analysis.rst#noctiluca.analysis.p2.P2traj), [analysis.p2.P2dataset()](../noctiluca.analysis.rst#noctiluca.analysis.p2.P2dataset), and [analysis.p2.P2()](../noctiluca.analysis.rst#noctiluca.analysis.p2.P2)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }