two_way_shooting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-Way Shooting\n",
    "\n",
    "This example shows how use two-way shooting in OpenPathSampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import openpathsampling as paths\n",
    "import numpy as np\n",
    "\n",
    "from openpathsampling.engines import toy as toys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start out by creating a simple toy potential energy surface, defining a couple of states, and creating a trajectory that connects the two states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a simple toy PES\n",
    "pes = (toys.OuterWalls([1.0,1.0], [0.0,0.0]) + \n",
    "       toys.Gaussian(-0.7, [12.0, 0.5], [-0.5, 0.0]) +\n",
    "       toys.Gaussian(-0.7, [12.0, 0.5], [0.5, 0.0]))\n",
    "topology = toys.Topology(n_spatial=2, masses=np.array([1.0]), pes=pes)\n",
    "engine = toys.Engine({'integ': toys.LangevinBAOABIntegrator(dt=0.02, temperature=0.005, gamma=0.000001),\n",
    "                      'n_frames_max': 5000,\n",
    "                      'n_steps_per_frame': 10}, topology)\n",
    "template = toys.Snapshot(coordinates=np.array([[0.0, 0.0]]),\n",
    "                         velocities=np.array([[0.1, 0.0]]),\n",
    "                         engine=engine)\n",
    "\n",
    "cv_x = paths.CoordinateFunctionCV(\"x\", lambda snap: snap.xyz[0][0])\n",
    "cv_y = paths.CoordinateFunctionCV(\"y\", lambda snap: snap.xyz[0][1])\n",
    "left_state = paths.CVDefinedVolume(cv_x, lambda_min=float(\"-inf\"), lambda_max=-0.3)\n",
    "right_state = paths.CVDefinedVolume(cv_x, lambda_min=0.3, lambda_max=float(\"inf\"))\n",
    "network = paths.TPSNetwork(left_state, right_state)\n",
    "\n",
    "# a fake trajectory\n",
    "fw_shot = engine.generate(template, paths.AllOutXEnsemble(left_state | right_state))\n",
    "bw_shot = engine.generate(template.reversed, paths.AllOutXEnsemble(left_state | right_state))\n",
    "trajectory = bw_shot.reversed + fw_shot[1:]\n",
    "\n",
    "initial_conditions = paths.SampleSet([\n",
    "    paths.Sample(replica=0,\n",
    "                 ensemble=network.sampling_ensembles[0],\n",
    "                 trajectory=trajectory)\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADuNJREFUeJzt3V+MnXWdx/H3h9a6JIKolD8CY0ms\nJGgMxBN2zQZwd2msNy0XRDG6loRNLwhxV3dJajBmgzeIMXohFza6CXqDQvzTKApa/2SzEcIQjC4a\naGHXpZTYEv8kLqCyfPeip3aYnplzOs/Tc2b4vV9JM+c5z+88v19/DO8+PTOTpqqQJLXllFkvQJI0\nfcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQetnvYClnHnmmbVp06ZZL0OS1pSH\nHnromaraOG7cqo3/pk2bmJ+fn/UyJGlNSfLLScb5to8kNcj4S1KDjL8kNcj4S1KDjL8kNcj4S1KD\njL8kNcj4S1KDjL8kNcj4S1KDjL8kNcj4S1KDjL8kNcj4S1KDjL8kNcj4S1KDeol/kq1JHk2yP8mu\nZcZdk6SSDPqYV5K0Mp3jn2QdcDvwLuBi4L1JLh4x7jTgg8ADXeeUJHXTx53/ZcD+qnqiqv4I3Als\nHzHu48BtwPM9zClJ6qCP+J8HPLng+MDwuT9LcilwQVV9s4f5JEkd9RH/jHiu/nwyOQX4NPDPYy+U\n7Ewyn2T+8OHDPSxNkjRKH/E/AFyw4Ph84OCC49OAtwA/TPLfwF8Be0Z90beqdlfVoKoGGzdu7GFp\nkqRR+oj/g8DmJBcm2QBcC+w5erKqfldVZ1bVpqraBNwPbKuq+R7mliStQOf4V9ULwI3AvcAvgK9U\n1SNJbkmyrev1JUn9W9/HRarqHuCeRc99bImx7+hjTknSyvkTvpLUIOMvSQ0y/pLUIOMvSQ0y/pLU\nIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMv\nSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUIOMvSQ0y\n/pLUIOMvSQ0y/pLUIOMvSQ0y/pLUoF7in2RrkkeT7E+ya8T5Dyf5eZKfJtmb5A19zCtJWpnO8U+y\nDrgdeBdwMfDeJBcvGvYwMKiqtwJ3A7d1nVeStHJ93PlfBuyvqieq6o/AncD2hQOq6gdV9ezw8H7g\n/B7mlSStUB/xPw94csHxgeFzS7ke+HYP80qSVmh9D9fIiOdq5MDk/cAAuHKJ8zuBnQBzc3M9LE2S\nNEofd/4HgAsWHJ8PHFw8KMlVwM3Atqr6w6gLVdXuqhpU1WDjxo09LE2SNEof8X8Q2JzkwiQbgGuB\nPQsHJLkU+BxHwn+ohzklSR10jn9VvQDcCNwL/AL4SlU9kuSWJNuGwz4JvAq4K8lPkuxZ4nKSpCno\n4z1/quoe4J5Fz31sweOr+phHktQPf8JXkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWp\nQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZf\nkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk/CWpQcZfkhpk\n/CWpQb3EP8nWJI8m2Z9k14jzr0zy5eH5B5Js6mNeSdLKrO96gSTrgNuBLcAB4MEke6rq5wuGXQ/8\npqremORa4BPAe7rO/XLz9Yef4pP3PsrB3z7H6884lZveeRFXX3rerJc1kWmtvcs8k752uXErPdfH\n+vt4/Wqbp29rdd2z0Dn+wGXA/qp6AiDJncB2YGH8twP/Onx8N/DZJKmq6mH+l4WvP/wUH/nqz3ju\nT/8HwFO/fY6PfPVnAKv+k3daa+8yz6SvXW4csKJzC/9w6LJPa2GfZ2mtrntW0rW/Sa4BtlbVPwyP\n/x74y6q6ccGY/xyOOTA8fnw45pmlrjsYDGp+fr7T2taSv771+zz12+eOe/4v1p/C5W/aOIMVTe7f\nHzvM8y+8eNzzfa99uXmuWDRP8tIxP3rsMM//acRrX3EK73jTWX8+/uFjh5YcB6zo3N9cdOT6P3h0\n6WsvXMNSllvbJK+f1LTm6dtS6z7vjFP5j11/O4MVzUaSh6pqMG5cH3f+GfHc4j9RJhlDkp3AToC5\nubnuK1tDDo4IP8DzL7zIgd+MPrdajAry0ef7XPty8/zPr59d/rUjonD0+f965n8nGnei1z567vHD\nvz+hNZzoPJO+flLTmqdvS617qf+3WtdH/A8AFyw4Ph84uMSYA0nWA68Gfr34QlW1G9gNR+78e1jb\nmvH6M04deed/3hmn8u1/vHwGK5rcUn9r6Xvty83znX+6YsWvvfdDV0w0DljRufs+dOUJraHr76Gr\nac3Tt6XW/frhfx+9VB/f7fMgsDnJhUk2ANcCexaN2QPsGD6+Bvi+7/e/1E3vvIhTX7HuJc+d+op1\n3PTOi2a0oslNa+1d5pn0tcuNW+m5Ptbfx+sntVY/F9fqumel851/Vb2Q5EbgXmAd8G9V9UiSW4D5\nqtoDfAH4UpL9HLnjv7brvC83R78gtRa/U2Faa+8yz6SvnWTcSs913ae1sM+ztFbXPSudv+B7srT2\nBV9J6sOkX/D1J3wlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwl\nqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHG\nX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUHGX5IaZPwlqUGd4p/k\ntUm+m2Tf8ONrRoy5JMmPkzyS5KdJ3tNlTklSd13v/HcBe6tqM7B3eLzYs8AHqurNwFbgM0nO6Div\nJKmDrvHfDtwxfHwHcPXiAVX1WFXtGz4+CBwCNnacV5LUQdf4n11VTwMMP5613OAklwEbgMc7zitJ\n6mD9uAFJvgecM+LUzScyUZJzgS8BO6rqxSXG7AR2AszNzZ3I5SVJJ2Bs/KvqqqXOJflVknOr6ulh\n3A8tMe504FvAR6vq/mXm2g3sBhgMBjVubZKklen6ts8eYMfw8Q7gG4sHJNkAfA34YlXd1XE+SVIP\nusb/VmBLkn3AluExSQZJPj8c827gCuC6JD8Z/rqk47ySpA5StTrfXRkMBjU/Pz/rZUjSmpLkoaoa\njBvnT/hKUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1\nyPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhL\nUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1yPhLUoOMvyQ1qFP8k7w2yXeT7Bt+fM0y\nY09P8lSSz3aZU5LUXdc7/13A3qraDOwdHi/l48CPOs4nSepB1/hvB+4YPr4DuHrUoCRvA84G7us4\nnySpB13jf3ZVPQ0w/HjW4gFJTgE+BdzUcS5JUk/WjxuQ5HvAOSNO3TzhHDcA91TVk0nGzbUT2Akw\nNzc34eUlSSdqbPyr6qqlziX5VZJzq+rpJOcCh0YMeztweZIbgFcBG5L8vqqO+/pAVe0GdgMMBoOa\n9DchSToxY+M/xh5gB3Dr8OM3Fg+oqvcdfZzkOmAwKvySpOnp+p7/rcCWJPuALcNjkgySfL7r4iRJ\nJ0eqVue7K4PBoObn52e9DElaU5I8VFWDceP8CV9JapDxl6QGGX9JapDxl6QGGX9JapDxl6QGGX9J\napDxl6QGGX9JapDxl6QGGX9JapDxl6QGGX9JapDxl6QGGX9JapDxl6QGGX9JatCq/Ze8khwGfjnr\ndQBnAs/MehGrhHtxjHtxjHtxzGrYizdU1cZxg1Zt/FeLJPOT/JNoLXAvjnEvjnEvjllLe+HbPpLU\nIOMvSQ0y/uPtnvUCVhH34hj34hj34pg1sxe+5y9JDfLOX5IaZPwXSfLaJN9Nsm/48TXLjD09yVNJ\nPjvNNU7LJHuR5JIkP07ySJKfJnnPLNZ6siTZmuTRJPuT7Bpx/pVJvjw8/0CSTdNf5ck3wT58OMnP\nh58De5O8YRbrnJZx+7Fg3DVJKsmq+w4g43+8XcDeqtoM7B0eL+XjwI+msqrZmGQvngU+UFVvBrYC\nn0lyxhTXeNIkWQfcDrwLuBh4b5KLFw27HvhNVb0R+DTwiemu8uSbcB8eBgZV9VbgbuC26a5yeibc\nD5KcBnwQeGC6K5yM8T/eduCO4eM7gKtHDUryNuBs4L4prWsWxu5FVT1WVfuGjw8Ch4CxP2CyRlwG\n7K+qJ6rqj8CdHNmThRbu0d3A3yXJFNc4DWP3oap+UFXPDg/vB86f8hqnaZLPCzhyc3gb8Pw0Fzcp\n43+8s6vqaYDhx7MWD0hyCvAp4KYpr23axu7FQkkuAzYAj09hbdNwHvDkguMDw+dGjqmqF4DfAa+b\nyuqmZ5J9WOh64NsndUWzNXY/klwKXFBV35zmwk7E+lkvYBaSfA84Z8Spmye8xA3APVX15Fq/yeth\nL45e51zgS8COqnqxj7WtAqP+4y7+9rhJxqx1E/8ek7wfGABXntQVzday+zG8Ofw0cN20FrQSTca/\nqq5a6lySXyU5t6qeHgbt0IhhbwcuT3ID8CpgQ5LfV9VyXx9YlXrYC5KcDnwL+GhV3X+SljoLB4AL\nFhyfDxxcYsyBJOuBVwO/ns7ypmaSfSDJVRy5abiyqv4wpbXNwrj9OA14C/DD4c3hOcCeJNuqan5q\nqxzDt32OtwfYMXy8A/jG4gFV9b6qmquqTcC/AF9ci+GfwNi9SLIB+BpH9uCuKa5tGh4ENie5cPj7\nvJYje7LQwj26Bvh+vfx+eGbsPgzf5vgcsK2qRt4kvIwsux9V9buqOrOqNg0bcT9H9mXVhB+M/yi3\nAluS7AO2DI9JMkjy+ZmubPom2Yt3A1cA1yX5yfDXJbNZbr+G7+HfCNwL/AL4SlU9kuSWJNuGw74A\nvC7JfuDDLP/dYWvShPvwSY78Lfiu4efA4j8kXzYm3I9Vz5/wlaQGeecvSQ0y/pLUIOMvSQ0y/pLU\nIOMvSQ0y/pLUIOMvSQ0y/pLUoP8HrTVcaSdBRfUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1112bbf10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(cv_x(trajectory), cv_y(trajectory), '-o')\n",
    "plt.xlim(-0.5, 0.5);\n",
    "plt.ylim(-0.5, 0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shooting point modifiers\n",
    "\n",
    "One of the things that distinguishes two-way shooting from one-way shooting is that two-way shooting involves modifying the shooting point snapshot. There are a lot of options for how you might change that snapshot; here we have implemented a few, and we'll show how they work.\n",
    "\n",
    "All `SnapshotModifier`s have an initialization option `subset_mask`, which allows you to specify a subset of the total atoms that this applies to. This allows you to, for example, only change the velocities of the solute atoms. This is not illustrated for this simple model, since the 2D toy has only one (pseudo-)atom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"Direction\" changing modifiers\n",
    "\n",
    "Usually, when people are first introduced to path sampling (and to the shooting move in particular), they are told that you perturb the velocity at the shooting point. That's what the `VelocityDirectionModifier`s do. However, it is very hard to do this correctly in the presence of constraints, so if your system involves constraints, please use the thermal randomization approach.\n",
    "\n",
    "There are a couple versions of velocity direction modification, although in this toy model they do the same thing:\n",
    "\n",
    "* `VelocityDirectionModifer`: changes the direction of all atoms from the subset\n",
    "* `SingleAtomVelocityDirectionModifier`: changes the direction of a randomly selected atom from the subset\n",
    "\n",
    "In this example, we'll only use the first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in a real system, you remove linear momentum. in this toy, doing so gives velocity=0!\n",
    "# the default is remove_linear_momentum=True\n",
    "direction_modifier = paths.VelocityDirectionModifier(delta_v=0.04, remove_linear_momentum=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thermal randomization modifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = template.instantaneous_temperature  # use the temperature of the initial snapshot\n",
    "thermal_randomizer = paths.RandomVelocities(beta=1.0/temperature, engine=engine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manually building a two-way shooting mover\n",
    "\n",
    "Here we will build the two-way shooting mover, and show its output, as well as how different modifiers sample differently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "tps_ensemble = network.sampling_ensembles[0]\n",
    "dir_change_mover = paths.TwoWayShootingMover(ensemble=tps_ensemble,\n",
    "                                             selector=paths.UniformSelector(),\n",
    "                                             modifier=direction_modifier,\n",
    "                                             engine=engine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to shoot several times from the same initial trajectory, to show the variety of shots we can get with this mover."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "change = dir_change_mover.move(initial_conditions)\n",
    "changes = [change] + [dir_change_mover.move(initial_conditions) for i in range(4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8W9Wd9/HPT6tlyZJtebfjLBAC\nIdBAwpISoGVJadlKS9tAFzoPHTpPy1Napp2203aaobN0VkqnDCXTdqYt7VBKF6BpgbJDCGkIhCUJ\ngWwk3pJ4t2zLsqTz/CFZ1rWk2I7lLf69Xy+9pHt1r3SuEp/vvefce64YY1BKKaWG2Ka7AEoppWYW\nDQallFIWGgxKKaUsNBiUUkpZaDAopZSy0GBQSillkZdgEJHLRGSXiOwWka9kef8CEXlJRKIicu2I\n92Iisi35eDAf5VFKKXXsZKLXMYiIHXgTuBRoALYA1xljdqQtswDwA18EHjTG3J/2XsgY45tQIZRS\nSuWNIw+fcTaw2xizF0BE7gWuBlLBYIzZn3wvnofvU0opNYnyEQy1wMG06QbgnHGsXyAiLwJR4NvG\nmN9mW0hEbgJuAvB6vStOPvnkYyyuUkrNTVu3bm01xpSPtlw+gkGyzBtP+1S9MaZJRBYBT4jIa8aY\nPRkfaMx6YD3AypUrzYsvvnhspVVKqTlKRN4ey3L56HxuAOalTdcBTWNd2RjTlHzeCzwFnJGHMiml\nlDpG+QiGLcBiEVkoIi5gLTCms4tEpERE3MnXZcB5pPVNKKWUmnoTDgZjTBS4GXgE2AncZ4zZLiK3\nichVACJylog0AB8C7haR7cnVTwFeFJFXgCdJ9DFoMCil1DSa8Omq00H7GJRSavxEZKsxZuVoy+mV\nz0oppSw0GJRSSlloMCillLLQYFBKKWWhwaCUUspCg0EppZSFBoNSSikLDQallFIWGgxKKaUsNBiU\nUkpZaDAopZSy0GBQSillocGglFLKQoNBKaWUhQaDUkopCw0GpZRSFhoMSimlLDQYlFJKWWgwKKWU\nstBgUEopZaHBoJRSykKDQSmllIUGg1JKKQsNBqWUUhYaDEoppSw0GJRSSlloMCillLLISzCIyGUi\nsktEdovIV7K8f4GIvCQiURG5dsR7N4jIW8nHDfkoj1JKqWM34WAQETtwJ/BeYClwnYgsHbHYAeCT\nwM9HrFsKfBM4Bzgb+KaIlEy0TEoppY5dPo4YzgZ2G2P2GmMiwL3A1ekLGGP2G2NeBeIj1n0P8Edj\nTLsxpgP4I3BZHsqklFLqGOUjGGqBg2nTDcl5k72uUkqpSZCPYJAs80y+1xWRm0TkRRF58ciRI2Mu\nnFJKqfHJRzA0APPSpuuApnyva4xZb4xZaYxZWV5efkwFVUopNbp8BMMWYLGILBQRF7AWeHCM6z4C\nrBGRkmSn85rkPKWUUtNkwsFgjIkCN5Oo0HcC9xljtovIbSJyFYCInCUiDcCHgLtFZHty3XbgWyTC\nZQtwW3KeUkqpaSLGjLU7YOZYuXKlefHFF6e7GEopNauIyFZjzMrRltMrn5VSSlloMCillLLQYFBK\nKWWhwaCUUspCg0EppZSFBoNSSikLDQallFIWGgxKKaUsNBiUUkpZaDAopZSy0GBQSill4ZjuAiil\nlIKdzz7Js/f+hJ62VoqCZZy/9hOccv67p6UsGgxKKTXNdj77JI+u/x7RyAAAPa1HeHT99wCmJRy0\nKUkppaaBMYZIfx/drYd5+p4fpUJhSDQywLP3/mRayqZHDEopNQHxWIxwb4hwqIdwKES4N/kcSs5L\nTg/0hugPDb8Oh3qIx2JH/eyettYp2gorDQal1Iw1Ve3uxhgGB8LDlbulgu8h3BtiwFLR9yafe4j0\n9x/1s92FXtxeHwU+HwW+IopKy1KvC3xFFHh9PPvz/6G/pztj3aJgWd63dSw0GJRSM9KxtLsP770n\nKvGhPfP+tL30oYo+vdIPh0LEY9GcZbHZHZbK3FdaSln9/FTFnnrPO7SMLxEGXh82u33UbXU4nZZt\nBXC43Jy/9hPj+cnyRoNBKTUjPXvvT7K2uz/2w7to2Ll9RPNNKLn33nfUz3R5Ci2VeVn9gkSl7rXu\nwQ9V9G6vD4+vCIfbjYhM2rYOBd1MOStJb+2plJo2g+EwXYdb6Dx8iK5DLXQdTjw6D7XQ3ngw53qF\ngWLL3nl6xZ6ozIcr9vRlxrL3fjwb66099YhBKTVp4vEYofb2RIWfqvgP0Zmc7uvqtCzv8ngIVFRR\nWlNHqL01a/t9UVk5N93531O1CXOSBoNSCjj2jt6Bvj5Lxd95+FBquvvIIWLR4bZ7ERtFZeUUV1Zy\nwoqzCVRUEaisojj5XOArSjXZjOxjgOltd59LNBiUUkft6F3yzgvoaTuS2NMf2utPC4HwiLNpCrw+\nApVVlM9fyIlnr0pU+hVVBCoqKSorx+4YW7Uz09rd5xLtY1BKsf6zf0ZP65GM+WJLXANr4vHUPJvd\njr+8gkBFFcWVw5V+IBkABT7flJVbjY/2MSiljmpwIEzzW7s4uOO1rKEAiUA455oPpyr94soqfMEg\nNtvc7sQ93mkwKDVHDIbDNL35Bgd3vEbDztdofutN4rEoIjZsDgfxaOZ5/EVl5azWNv05R4NBqVlq\ntM7iSLifpl07adj5Oge3v0bLnrcSQWCzUbnoRFZcfjXzlp5GzZKl7N26WTt6VYoGg1KzUPbO4v/g\nyIF9iNg4uOM1Du3dTTwWw2a3U7noRFZe8X7qlp5G7ZJTcHkKLZ+nHb0qnXY+KzUL5eoshsTwDVUn\nLGbeqadRd8oyapacgqvAM8UlVDPRlHY+i8hlwB2AHfiBMebbI953Az8BVgBtwEeMMftFZAGwE9iV\nXPQFY8xf5KNMSh2vIv19OUMB4OYf3YuzoGAKS6Qm24anvsEde39Diw2q4nDLomu4/F3fmrTvm3Aw\niIgduBO4FGgAtojIg8aYHWmL3Qh0GGNOFJG1wD8BH0m+t8cYs3yi5VDqeNd1+BAvP/wQrz3xaM5l\nisrKNRSOMxue+gbr9v2GsD1x4V+zHdbt+w3ApIVDPo4YzgZ2G2P2AojIvcDVQHowXA2sS76+H/ie\nTOaIVEodJ4wxNO3aydbf/5bdf3oBBE46dzWltXVseeBX2lk8GxkDkRCR3iO0dx+kraeRtt5DtPUd\npi3cRttAJ22RbtqjvbTFwuyN9xO3WavLsE24Y+9vZnQw1ALpo101AOfkWsYYExWRLiCYfG+hiLwM\ndANfN8Y8m+1LROQm4CaA+vr6PBRbqZll5FlGi848i5bdb3Fo71sUeH2svOoDnPGeK1Jj9JdUVmtn\n8XSIx2CgG8JdyUfitenvpL+/LVHB97fSNtBBW6SbtsEQbbF+2uMR2kyUdpuhzWanx579BpqF8TjB\nmCGIMB8Hu3PU0i2TeP/NfARDtj3/kT3auZZpBuqNMW0isgL4rYicaozJuGOFMWY9sB4Snc8TLLNS\nM0q2s4xeefT3FBaXcPGNn+HUCy7KaCI65fx3z80gePU+ePw26GqAQB1c/Ddw+ofHvn40MqJiT3uM\nmG/6u+ge6EzsxQ+GaIv10WYitNnttNvttNlsiWe7nTa7jbAtS20tEHDYCYqHoKOAJQ4fQZefoLuY\noKeM0sJygr5qgkV1lBbNw+OrAIc7tfqaHy2jOcv1hFXxzHn5ko9gaADmpU3XAU05lmkQEQcQANpN\n4pSoAQBjzFYR2QOcBOgpR2pOyXbvAQC7w8HyNe+bhhLNUK/ex4bHvsQd/kJaSuqoisa45Y9f5PKD\nm6Hy1BEVfZbKf6Cb2GAfHXYbbbbhCj29cm+z22l3uhLPAlE7UAhgAxLDfdgRShyFlDqLCLqLqfcE\nCXrKKfVWEvTVEPRWECwIEvQEKXGX4LQ7j3mTb1l0TaKPIa05qSBuuGXRNRP6KY8mH8GwBVgsIguB\nRmAtcP2IZR4EbgA2AdcCTxhjjIiUkwiImIgsAhYDe/NQJqVmjb6uzpxnGU3XPX+njTHQ3wHdjdDV\nmHhOe72hdRvrgoHUnnmz08G60iKiO37OOS8NJCp3p5M2t482t4c2h4s2p512t522QAnt+OiIRzKa\nNACcNkeyMi+j3BPk5IIgpQWlBD3BVCU/NF3sLsYmk9iWk2aoH2FWnZWU7DO4GXiExOmqPzLGbBeR\n24AXjTEPAj8Efioiu4F2EuEBcAFwm4hEgRjwF8aY9omWSanZYt/LL/LwXd/J+f503fN3UqQq/aYR\nFX4TdDcMv46OuAeD2OkvqubtQAXfLvVnNNeEbTa+XpHtd4pT6LAR9JQSLAhSX1DKGZ5gqqIfWen7\nnL5JvUvbRFz+rm9NahCMpBe4KTVFLJ3LpUFKauo48No2yubN5+TV7+KFX92bcZbRmptunp5+hPG2\n4xsD4c5ExZ6+p9/dlPiMoTAYHHHrTbFBUTX4azFF1RwpKmef281+G+yL97Mv0sH+3iaaepuPXl5j\n+OY716Uq+aG9e49DL+xLp6OrKjWDZHQut7XS09bK/Hecyfu/+HUcLhf+YNnMOMvo1fvgoc/BYHLP\nvesgPPi5RIVftSxzT38oDAZ7rZ8jNvBVQaAWKpfC4jXgr2GgqIK3HU72m4FExd9zgH1d+9jf/Sa9\nR15OrV7oKGRBYAFnVJ7JNf4FLAws5Nsb19EaDWUUudpVzLUnXTuZv8qcMmeCoeuhhzh8+3eINjfj\nqK6m4gufJ3DlldNdLDVH5Opcbm88iMPlAmbIWUbhLnjkr4dDYUi0Hx5flzZDoKgK/LVQfjKceAn4\naxLT/lqMv4Y2h4t9oaFKf3/iufllGt9qxKS18ld7q1kYWMjVJ1zNwsBCFgQWsNC/kIrCioymnWg8\nyrrnvkHYDKbmFYiTW8796iT8GHPXnAiGroceovkbf4MJhwGINjXR/I2/AdBwUFMiVyfytHYuR/qg\n5TVoegkaX0o8t+0GYIO3kDtKimlx2BNn/nR0cnlvH/yfRxKVf1EV2J0MxgY50HOA/V372de9j32d\nW9n/9q/Y17WPnsGe1Fd5HB7m++dzWtlpXHXCVYnKP7CQ+qJ6Cp2FuUqY4fJFlwNwx0t30NLbQpW3\nilvOvCU1X+XHnOhjeOuii4k2jTyDFhDBXlKCuN3YXC7E7U4+XNhc7tS0ze1CXGnvud05ppOvc0zb\n3MnvcLkQ+9Td6ESPlqZXPB7je5/8CIMD4Yz3Jnxj+7H2BUQjcHh7MgBeTjwO7wQTSxakGmrOgJoz\n2fDKf7GuyGnp5HXH43wyLNRcchv7uvaljgIaehqIDX0GUFFYkdjrTzb9LAwsZKF/IZXeyik7i0fl\npn0MaaLNOTqujMF/2XuIDwxgBiKYgQHikcTreH8/prMz+d5A8r1I6jUTDVSn0xJGkxVMvc9vovV7\n30uUGT1ammrRwUH+8L1/Y3AgjM1uJx4brkQnPIRFtr6Ahz4HJg7V7xg+Cmh6OXFkEIsklvOUQM2Z\nsOS9qTDAX5362O80/pbwoPUa0wGbjbsLgee/icvmYn5gPktKlnDZgstSe/8L/AvwOr3Hvj1qxpjT\nRwyOmhoWP/H4uL/fGAODg5agiA8MYCzTEUwk7b1k8JjIQNYgGvW9fAeTzYZr4ULsPh+2oiJsRT7s\nviJsRUXYi3zYfMl5RUXYfMl5RUXYfIl5MsYbuo/V8XRUkzr7qLUVu9NBbHCQCz9+I95AcX47l29f\nlgiDDEJq8AGXD6qXQ20yAGrOgJIFMKLt/mDPQZ5rfI6NjRt5uuHpnF/5hw/8gWpvNXa9teespEcM\naSq+8HlLHwOAFBRQ8YXPH9PniQi4XNhdLpiGG5+PJ5gaPvPZ7B8Sj+M+4QTioR5iXV0MNjQQC4WI\n9/Skji6ORjwea6h4RwbMcKhkCxi7z4ckO12Ppz6gkWcfxQYHsTkceAPF+e1cHgznCAUAA+//PtSe\nCcETIUsl3h/tZ0vLFjY2bmRj00be7n4bgDpfHYWOQvqifRnrVHurqSuqy0/51Yw2J4JhqHI5XvZI\nxxNMjpqanEdLdd+9I+s6JhJJhUSsJ5QIj54e4umvQ73J18nlQj0MtrQkX4cw/f1ZP9uyHW43tqIi\nYh0dkNbEAmDCYQ79wz/irKvDWVWFo7w870cpkyHb2UfxaJRn7/3JxEMhHocDm+DVe2H7A7mXC8yD\n5ddZZhlj2Ne1L3FU0LSRF1teJBKPUGAv4Kyqs7j+5OtZXbuaen89G/ZuYN3z6wjHhnekCuwF3HLm\nLRMrv5o1Zv5fWp4Errxy1gbBRBzL0ZK4XDhKS6G09Ji/1wwOJsLlqAGTeK/zl7/M+hmxjg7evi45\nuorNhqO8HEdVJc6qapxVlTiGniurcFYfPTymqqlqQkNb5OpIbn0LXrkXXrsPOg+A0wunXAlFVWx4\n7X8S4wYNnT3U3cflFyeOtkKREJtbNieOCho30tSb2EFYFFjE2pPXcl7teayoXIHb7rYUQ8/8UXOi\nj2Gum+nt97n6gOzl5dT8/d8x2NxC9FALgy2HiLY0M9hyiMGWFkzfiOYOmw1HWRmO6iqcQ2FRWUWk\nsYGuX96PiURSi0pBAdXfui2vv0PjGzu4d92Xs/b/jHr20ciOZACbM3FtQOfbiYvFFr0LTl8LJ18O\nbl9iz37EOf0ucXBR/SW0hlvZdngbUROl0FHIudXnsrpuNefVnEeNryZv26xml7H2MWgwqGk3so8B\nRq+4jTHEexLNV9GWluTzoeHpQ4eINjcTHxkeaezFxcz/2T24FizIevrwaIGaPsSFp6iIcKgXT1ER\nA319xAaHQ2hMQ1vk6ki2OeGSdXDatYlrB9KsuX8NzTmGilhSsoTVtas5r/Y8lpcvn9DonrPJm5tb\n2PTAHkLtA/hK3ay6+gROOqdq9BXnCO18VrPGsfQBiQh2vx+73w8nnZR1GWMM8VCIN88+J+tefKyz\nk72XX4EUFlJw8skULF2aeJx6KuGdO2lZty5nh/jITub+7m4QYdW11+P2eMZ/9lFXQ/b58Si88+as\nb+UKBUG4/6r7j/59x6E3N7fw5M/eIBpJ3Kgg1D7Akz97A0DDYZw0GNSMMBl9QCKCvagIR3V1zqaq\nir+8lfD2HYR37KDz17/G3HNPzs8z4TCHb/8OgSuvzD7EhTH86YFfctOd/z32jubYIPzpv3K/H7Ce\nBWSMYUvLFr7/6vdzrlLlPb4rwXgsTn/PIL1dA/R2DtDbFaG3a4Btjx1MhcKQaCTOpgf2aDCMkwaD\nOu7l6oCv/KsvJcLo/e8HwMRiRN5+m/D2HTR96UtZPyva1MTgoUO5O5lzzAcyO5eXfRB2/QFad0HF\nUmjfC9G0q6OdnkQHNIlA2Nyymbu23cVLh1+i3FPOFQuv4LEDjx03Zw/lqvD7ks+9nYnXfT2R7PeI\nzNEqHmof/fRrZaXBoI57Y22qErsd96JFuBct4vDtt2c9yoiJsOEjH4DKkqzf5RmM0vXQQ5lHP9mu\nUt74HSgsg+vuhZMug9d+yYZnb+MOdyxxlpGrmM95CyltfJ67XrmLbUe2UVFYwVfP/iofPOmDuO1u\nVu9dPSVnD02k7d5S4XdFkpX+cIXfl5zX3xPJbPET8BS58AZceANuKuqLKCx24w248QZcFAYSrwv9\nTn76jU1ZQ8BX6s6Yp45OO5+VyiJbh3hHiZ/tJy+kuy9EMC50ECOeNp6QLR7ntINHmF/oz7yiPlfn\nsr8Obt0OkPX6AUEwGKq8VXxq2ae4ZvE1uOyu/G7sKEa23QM4XDYuvG4J804pzUuFn6rki60Vvsfv\nxG4f2xhLucr57o+erE1JSdr5rNQENBX7eObMJYRCPRQMRvEZG61uO36vlw9+4cssOP0MHlt1Fruq\nSgk7HRQMRlnS3E5tZ4hoV2/mGU31rQQWZPmi7sbUyzteusMSCgAGQ8AVYMM1G6Y8EABisTgbf707\na9v94z/embnCiAq/bJ4vUfEXuyn0u5IV//gq/LEaqvz1rKSJ02BQagTLGUcihF1OwsD808/gqr/8\na1wFibuCzff4qd15IPMDjKHpy19JXKlM8oymQ8UABBaMuCI8rXO5pbcla3m6I91TEgp93RHaGkK0\nNoZoawzR2hCio7mXeCx3q8KF1y9Jq/BdePyuvFf443HSOVUaBHmgwaDUCDlvqtPUkAoFyNGp7XaD\niWEiUcu6JiYcftVvDYa0zuXDfYdx2BwMxgcZaaxnGY21HyA2GKe9pZe2xhBtDckQaOylv3v42gtv\nwEWwzkf90lJ2bmwi3BvN+BxfqZtlF9SOqWxqdtFgUCqp6/AhXnvikTEPa5G1U/uD59L0H78icZqM\nVbTPTtvbNbS/EiXa58BRFqDidA97q17m1qduBQNOm9MSDmM9yyjXOfzhvkECFYWJI4FkCHS29BGP\nJ44C7A4bpTVe5i8LUlbrI1jrJVjrw1M0fIRSVuvL2na/6uoTRi2Xmp00GNScFo/F2PPSn3j1sYfZ\n/8pLCILD5SKaNnzGkKJgWca8jOsvbl/G4cIY0b7sf1qHN8HQn120tYuGr3+NH1xmKDynnvVXrufN\njjeP6SyjTQ/sydoP8Owv3kpN+0rcBOt8LDi9jLI6H8FaH8UVHmyjNP1o2/3co8Ggjmvpw1akX4Xc\n3XqY1554lNefeJRQRzu+0iCrPriWZe9eQ+PO1y1XNcM4bqrT1UDF6QU0bwlgYsMVrtjjiLeYeLf1\nBji2gUH+7Akni9fdQ8BTwuKSxazeHufwf36HaHMDjup/p+sL8ZwX/3W39bNvW+tRz9W/5i/PoLTG\nR4H32IfF0Lb7uUWDQR23Rg5b0dN6hEe+/122PPRrWg+8jcGw8B1ncvGnPsuiM1ZiS46X5E9etTzq\nsBaWC9ZqoWYFYFL9CIdfLSLaZ8dRGKPiXCeNT3ZnaWACb88ghy65kr6LL8bm99Nxzz05h+IwxtDe\n1MvebUfYu+0IrQdDANjskrWT2FfqpmZx9msulMpFg0Edt7J1Iseig7QeeJuz3/8hTrtoDYGKyqzr\njnpTnYwL1hoSj5IToKeJwIL+4Y5mpweu/C6vb72N0q5YxkeFCm1Un30WXb/7Haavj5aKlexZfhUD\n7lLcA+0s2vsQPXf9ksHIKezddoSuI/0gULUwwKoPnMCi5eUc3tet/QAqbzQY1HElEu6n+a1dNL6x\nI2cnsjGG1Ws/PrEvevw26xDZQ+IRuOq7We+rcM+Ff8tNv4eCtBN8wg740SVw1z/fTjwc5uk1f84b\nS64nnrxHwkBBkJ2n3AAiyB/3U10pnHZFDSesXoivuCD1OcUVhfS9/BJbt4QJOwIURLtYsbxAm3/U\nMdFgULNaqL2Nxl07aNy1g6ZdOzm8fy8mHgcRbHY78VjmHnpRWWYn8rjlGg21qyFxc53TP2yZ3dbf\nxqZlTgwRrn/KEOyGNj/8/F3C3rMTp3zaCgrYu/gDqVBIEcExGGLVlm/hjITgF9BUVkbBqUvxnLqM\ngmXLiDQ1UfD9f+Wd6afObi2gqzY2o+69oWYHDYbjTO9AlLuf2cM9m96mo2+QkkInH1s1n09fcAJe\n98z65x5vWU08TlvDgWQQ7KTxjR10HzkEgMPtpvrEJZxzzYepPekUqk86mZ2bX+Dx//peYgTTIXYn\nZ1/7sWMv28oA3uf/FYPJ2l8Q99fSP2LdYm8U/6IfALDp1AI2njpcHru4+MZpNxOLxtn5fDNhpz9r\neaJOL8s2P0P4jTcSo8G+/jrh7dtpffa51IV0Gb9XOEzLP/8LRRdfjK2wcNRtPubfZIb935ot5ZzJ\n8jJWkohcBtwB2IEfGGO+PeJ9N/ATYAXQBnzEGLM/+d5XgRuBGPA5Y8wjo32fjpWUXe9AlGv+cyNv\nt/UxEB2uLNwOG/ODhfzmM+fNmD+MobI69m3jrNZNFMVC9Nh9bClbRXThcn7zmfNwSYyW3W/StGtn\n4ojgzZ0M9PYC4C0uoWbJKdQuOZXaJadQvmAR9rTbeo7l83P9FkPrntb+KF+Qe6mRVppNkJc5iXfb\nXqFQBnjJtoyl0Z14ZPi01n7j4vaCm3nCfSGNkY3Ygg8jzk4wdpA4g41/xpJ+D6u6KyiKBOhxdfFC\n0WHKnUs4P+Ii1D5AjMQf0Uh9DvjMv12QUeZ4Xx+dr26n5ZOfyBpUQ+zFxThqqnFW1+Csrk48aqpx\nVFXhrKnBUVaW9WZFI3+T+pef5WOvbaC8v5MjnmLuOe1yDpxx/oz5vzWb/gamw5SNlSQiduBO4FKg\nAdgiIg8aY3akLXYj0GGMOVFE1gL/BHxERJYCa4FTgRrgMRE5yRiTefyvRnX3M3sy/iAABqJx3m7r\n4+5n9nDrpUumqXRWdz+zB8e+bVxw+EmcJtHo7o+FuPDQ4zR27+L7f/kLbO1NxGOJ94J19Zx07mpq\nlyyldslSApVViOSuClO/ReGJ7Kg/0fKee5Tf4u5n9nBa+6N8y7aewmTFXytt1LKJ7fEF/KDiq/y+\nJcB74s/wV477qJE2mkyQf45+mA2Rldj9z+Gq+jViSx4ZSAwTt3Nir4eL2hfjTFbh/kgxl7YFEAzd\nJXFi5wV5fFsLF/c6UssADGJ4piCKPUuZbYWF/E+PnxWFJVT0dWRsS7erkOb3fJBzfINEm5oZPHiQ\nvj/9iXhPj3VBhwNnZSXO6mprgNQkQuSHO0PUv/wsn916HwXJI7DK/k4+u/U+7gTufqZqRvzfmk1/\nAzPZhI8YRGQVsM4Y857k9FcBjDH/mLbMI8llNomIA2gByoGvpC+bvtzRvlOPGLI787ZHae/LHFIh\nnd12tP3KhNGXOMq6o6xsMzG80RDegW7ec/hRPPHM8+8N0Oyu4pCnhhZPNYcKqhiwF4xasPS3R1YM\n2XicdkuZh9bvjcR43nUzNbb2jHUa4mWsjnz3qJ/rPeHb2FydGfOv3/q3+CPFGfP7xPCf/jBIYttP\nHrBzQdiB3wjdkgiFN9wxBCgqyNyX6wlHufDgVm7Zdn+q0gYI253csfxanp63ImM9T6Sfsr5Ogn0d\nlPV1EOzrpKw38bqsr4PSvk4cxvobxhFsWW560OP0cPdp7yfm9dHr8hByFdLr8tDrKmRwim8p2hOO\n5rotAwClXhcvfePSKSvPTDOVo6vWAunjCTcA5+RaxhgTFZEuIJic/8KIdbMOviIiNwE3AdTX1+eh\n2MefjlFCASAWn+Rh1o2hIB77oYLiAAAZJklEQVQmEO3GP9iDP9pNYLAbfzTxKIqGslYuI/2q5hrr\njNjoFf149Q9aD0y99HOJbStXOjdRLZmhAFAjbaN+rjgzQwGgKBLIOt9jwKSl2hvuGG+4Mw+aDdAd\nzhyzCOCpeSsA+OSOP6Saef5n6XtT80eu142TQwXlUFAOpZmfZzNxSsI9lPd3UNHXSXl/Bzdu35B9\nuwb7+eJL/5v1vQGbg5CrkJDTk3r0pE+PeN2T9joygVB518GtWX+Ljr7MK9pVpnwEQ7b9uGz3V8q2\nzFjWTcw0Zj2wHhJHDOMp4FxRUug86hFDqdfFn/764qN+xlh+2GgkQveRw3QfOUTX4Ra6D7fQdfgQ\n3YcT04Nh62mcHn+AQF0V/ooTCVRU4q+o4iuPNrDq4CP4Yr0Zn99j91HqdbHpqxcNl2uc/+Kr/vHx\nowZlSaGTjV+5CBPpx77njzh2/gb77keRaJhmgoSMhyIyT0dtMsGj3SwsIe7mxPZlnHPgCnyREkKu\nDrZXPkecOPYsPQjdYigpdGKAzlHK/NQXM6+tuPBfn6Szb5Cn5q1IBcFY1huPC//1Sa7cu5HK/szQ\nO1IQ4B8u/X/89MNLiXd1YXq6iXcnHt7ubkrSpuM9PZjuZuId3ZhQ6Ohf6nJh8/uxFfmR5LPNn3hI\nUVHqdWI68fyhn7/OCftf57Ov/dbS5HXLtsQ9sF89edWEfoe5Ih/B0ADMS5uuA0be+mpomYZkU1IA\naB/jumqMPrZqPnc/vTdrM4rbYeNj59bjGMOQyMYYejs76DrUQtfhFjoPtdB95BCdyelQu3Wv2eF0\nEaisIlBRybxTlxEoryJQWUVxRSX+ikrLiKRD1sR38dhDPZY+BoBBcbClbBUfO7cetyN3Z+hoPp72\nW1xley7ZF9BKkynju+ZDrFqymMLffQbe2ACREHgr4MwbYNkH+d83/DQ+e4+ljwGgz7i43azlrIUl\nvHKwK/vvHHyaxe3LuHDvWpzxxEB0RZFSzjl4JYPEMRgcI/oPNnljfHzVAoCj/vt9fNV8AoWZe9Gf\nGOXfPdd64/GJVfO5Z/fllj4GSDRX/eT0K1izZiXlZ42v7d7EYsSSgRHr7ibW1U28uyv1OtbdlXiv\nK/l++xGi+/Yk1hnZR5L0nzm+qyA2yJ/tfJitn/zIuMo4V+UjGLYAi0VkIdBIojP5+hHLPAjcAGwC\nrgWeMMYYEXkQ+LmI/DuJzufFwJ/yUKY56dMXnMDDr7fkPBPn0xcMXwUbCffTdfhQqvLvSu7tdx5q\nofvwIaKDaYfcIvhKgxRXVDH/tOUEKqqSQVBFcWUVhYHio3YE5y7rcp6BUcs6kd8i0Yn8g1QFXyet\n/BN3ITuBgmJY9oHEvZfnrwZ74s/h05VRrtm+hm+0kzwrKdG5fLtZy2ula7jnujP52A8383ZbH7HC\nrbjLH0k0H8U9iL2fc1/8+1QopH5ChAGbsLEwyqpee6r/YJM3xkBtQWp7H369JecZNbl+k6FtHe96\n4/09r3n9fO6ErGcl3XEM3yF2O46SEigZ/5AdJhYj3tOTESJ9bR20/923sjZFlPd15OW3mAvydbrq\n+4DvkDjT7kfGmL8XkduAF40xD4pIAfBT4AwSRwprjTF7k+t+Dfg/QBT4vDHmD6N9n3Y+57btyccz\nz9232Zm/fAUej4euQy10Hm6hv7vLsp7L40lU+BVDe/tVqaMAf1kFDlf+bxSTOt/8hQN09EUoKXTx\nsXPrJ36+uTHQ1UD4wFbkgc/gztJcRWEZ3LoTHNm3a7Sy9Q5E+auH/5tn2u8C2/BvbcPGTZv+nVw9\n5QMfrD3q9h7rbzJpv+UUf0c+vPnui4g1N2fMt1dXc9KTT0xDiWaOsXY+6z2fjzPrP/tnOYeCCFRU\nWvb2AxWVqQAo8BWNe69/RjAGOvZD8zZofgWaks/92TuPhwmsy95JPFZr7l9Dc6+1AnJFPXxi6204\n4pmB4yt1c8M/nDeh71Sjy3a/bikooPpbt835q8D1ns9z1MibyaSI8Kn/+OHUFibf4nFo35sMgW2J\nEGh5FcLJox+bEypOgZMvh5rlUL0c7rsBurMMX5F2S81j1dLbwolHVqQ6mfsdPdiMDVvckTHaqQ5o\nN3Wy3kDpC5+f86EwHhoMx5miYFnWI4ZsN5mZMSzDVycHnVv2QWh9K7H3P3Q00PwqRJKdjnY3VJ4K\np34gGQLvgIql4BgxztAl37SOggqWW2oeq3A0zCnt57Jq7wdS/QmFUT+GODsXPcstF/5fvbHNNMq4\ngZIaFw2G48z5az9x7DeZmQ6v3gcPfg6iQ8NXH4TffBp++xkYusWlwwNVp8E71g6HQPnJMJbz3IcG\ns8sy2umxago18fknP8+Z+z6cpZPZxmmt5+uNbdSspsFwnDllrDeZmQrGQH8HdDdCd1Pac9rr1rfI\nuCrAxMFZCO/7j0QQBBenzhg6JllGOz1Wm5s388Wnv0iwrZ6iSJYrw4BY9+inBCs1k2kwHIdGvclM\nPhgDfW3DFXxXQ2al3900fCQwRGxQVA3+mkR/QOub2T8/0gvLr5vcbRiDX/zuYQ481o8n7KfP3cnr\n5c9zaf+NBNvqEVsiw0bylbozZyo1i2gwHK+ytduPda85HofeI1n28NP3/JshNmKcI5tjuNKvfgcs\neS/4axPTQ8++Suve/+3LEs1HI+Whc3iifvG7h2n5PRTGE0NZeAdKOKfhcnDEOe/aE3EXOnjm3jf1\nrmnquKPBcDzKuO3kwcQ0JDp1Q4dG7NmPDIDm4fb9ITbncAVfuxJOSb4OpFX83nKwjfNq5Yv/ZlI6\nh/PhwB/7U6GQrs8eYvklifG67HabdjKr444Gw/Eo220nB/vhN3+ReIwc1dzuTlTugTqoX2Xdw/fX\ngL8OCoNgm4S280noHJ6IWDTOvh2HefLpF/EMZL9pjmegKPVaO5nV8UiD4XiU67aTJgbnf3FExV8L\nhaWjj5c9mfLYOXw0b25uybp3HxuMc3BnO7u2NrN72yEYsBGx2xBbBGc8s7+gv6B70suq1HTSYDge\nBepytNvPg4u/MfXlmQHe3NzCkz97I9UfEGof4PGf7uS1pxtoa+plMBwj4uhnX8mr2E/oZ+0lV7B/\nWzstv49YrmKO2iLUX5I5KKBSxxMNhuPRDG63ny6bHthj6SQGiEcNLXu72Fv1MjuLNzPv5CA3nfHn\nrKhMDF19Th38goc58FgXnrCf/oJu6i/x8JErLpuOTVBqymgwHI9mWLv9VItF47Q399LWEKL1YIjW\nhh5C7WGyDWxngL7Ve/iH07/G8orlGe9/5IrL4IrJL7NSM4kGw/Fqitrtp1u4dzARAA0hWg/2cKQh\nREdzb2qcIofTRrDOx6A9gjOW2V/Q5+7irkvumupiKzWjaTCoGSVXB7Exhp62MK0HQxxp6KH1YIi2\nhhA97cMjaBb6XZTN8zH/1CAltR7afc28HH6BXzdvpN/ustw8B2DQFmHTvAeBD0zDlio1c2kwqBkj\nWwfxYz/eyZbf76OvK0IknDjNVgSKKwupWuRn2YW1lNX5CNb5CDk6eb7pee5vfI4X9rxAd6Qbm9g4\nrew0Wmp28TT3Wm63ubn+d/Qu0BsGKjWSBoOaMrFonL7uCL2dA/R2DdDbGaG3a4C+5HTjm52WoaoB\nTNzQ3RZm6TtrKJvno6yuiNJaL06XnUgswkuHX+L+xud57unneKvjLQAqPBVcVH8R59Wex6rqVQTc\nATbs3cC659exu3xr6rML7AWsO3PdVP4ESs0KGgxqwuJxQ3/PUIU/XPH3JadDnQP0dQ3Q35N5o3ub\nTSgMuPAWuzNCIfX50TgXXp+4n/DB7oPcv3cDGxs38qeWP9Ef7cdhc7CiYgW3rriV82rPY3Hx4oyb\nDl2+6HIA7njpDlp6W6jyVnHLmbek5iulhukd3OaYXG342RhjCPcOpvbsezsHLJV/X3JeX3eEjP9G\nAoVFiQrfG3BRWOzGG0i89g69Lnbj8TkRW6IS//FfbyTUPpBRjh53O3ve9whNvU0c6DkAQJ2vjtW1\nq1ldu5qzqs6i0FmY199JqeOR3sFNZcjWhv/ET3fSuLsTf7CA3q5IqlmntzNCb/cA8WjmjkOB14m3\n2IU34Ka01oc34MJX7KYwWdl7A24K/U5s9rEPoWGMYdl7K3juf/dYLigbtEXYPO937Gl5ifPrzuej\np3yU1bWrqffXT/wHUUplpcEwyxhjGByIEemPMtAfJdKfeB0JR9PmDc8f6I8yGE48tzWGMoaJjkUN\nO55NdMC6CuyJir3YTc3iYrzFrkRln6rwXRQGXDic4xwoL00kFuFA9wH2de9jf9d+9nXtY3/3fvZ3\n7adnsIcTF63I6CDeXb4VQbjz4jsn8tMppcZozgbDeJpU8iUeN6lKOlWhj6i80yv0oco+VdEnp0dr\n/RMBl8eBq8CRePbY8RW7aT0YyrnOn3/nAlwF+fnvYIyhLdyWqPhHBEBjqJF4WjpVFFawMLCQ9y16\nHwsDC1n/6np+Vv63GZ9Z5dWB6pSaKnMyGLI1qTz5szcAcoZDLBZnsD+WtkeevfIesExHGUjbox8M\nx7J+djqbXZKVuh2Xx4Hb48Bf5klW8IlpZ4Edt8dhmZceAk63PaPzFXK34ftK3ccUCpFYhIM9B1OV\n/r6u4RDoGexJLee2u5nvn8/S4FLetzARAAsCC1jgX4DX6bV8ZrG7mHXPryMcG74+ocBewC1n3jLu\n8imljs2c7HzOVUE63Dbmnxq0VujJR3Qwy626RrA7bWkVtT1Vcbs8DtwFiUo7V2U+NM/utGWt1PNh\nZCACRG2DPLXof+ld0JT1LB1jDO3hdmvFnzwKaAg1WPf+PRWpSn9hYCEL/AtYEFhAtbcam4y9v2HD\n3g169pBSk2Csnc9zMhju/Isncr5XUlU4XHEPPQpGVOg59t7tjpl/r9/hJrQwIXcnL8x7KHVuv9vu\n5kMnfYigJ2g5CuiJWPf+6/31LPQPB8BC/0Lm++fjc/mma7OUUmOgZyUdha/UnbNJ5fp1505DiSZP\nLB6jY6CDtv422sJttJe3w0fb+Okrd9E72GtZdiA2wD077wGg3FPOwsBC3rvgvZYjgGpvNfbx3qVN\nKTWrzMlgWHX1CRlNKrPpXr2RWIT2cHuqsh/5PPRee7idjnAHhrEfFQrC89c9r3v/Ss1hczIYhjqY\nZ8q9eo0x9Ef7Myv65Ov0ir6tv83SsZvO4/AQLAhS6illXtE8llcsT0wXlBL0BAkWBAl6EtMfeuhD\nNPc2Z3xGlbdKQ0GpOW5CwSAipcAvgAXAfuDDxpiOLMvdAHw9Ofl3xpgfJ+c/BVQDQ3eUWWOMOTyR\nMo1V+r16N+zdwM0vfYKWN/LX2WmMoTvSnXWvvj3cntiz729PzUs/Cyed3+VPVeonlZxEsCaYqvyH\nKvqhyn88V//ecuYtevaPUiqriR4xfAV43BjzbRH5SnL6y+kLJMPjm8BKEvdF2SoiD6YFyEeNMdM2\nvsXQ4GpDFWRzbzPrnl8HkBEO0XiUzoHORCU/otlmaA+/vX947z5qohnfZxMbJe6S1J77PP88y578\nyMreaXdOynbr2EFKqVwmdFaSiOwC3mWMaRaRauApY8ySEctcl1zm08npu5PL/W/yiOGL4w2GfI6V\ntOb+NVmbVLwOLxfUXZDau2/rb6NzoDNre73T5rRU5iObbdKnA66Adt4qpabFVJ2VVGmMaQZIhkNF\nlmVqgfQ70zck5w35bxGJAb8i0cyUNalE5CbgJoD6+vyNk9PS25J1fm+0l9fbXidYEKS+qJ4zKs7I\nWfn7nL5Ju/ZAKaWm2qjBICKPAdl6Zb82xu/IVmMOVf4fNcY0ikgRiWD4OPCTbB9ijFkPrIfEEcMY\nv3tUVd6qrEcM1d5qfv+B3+fra5RSatYY9YosY8wlxphlWR4PAIeSTUgkn7N1HDcA89Km64Cm5Gc3\nJp97gJ8DZ09sc8bvljNvocBeYJmnnbBKqblsopfqPgjckHx9A/BAlmUeAdaISImIlABrgEdExCEi\nZQAi4gSuAF6fYHnG7fJFl7Puneuo9lYjCNXeata9c512wiql5qyJ9jF8G7hPRG4EDgAfAhCRlcBf\nGGM+ZYxpF5FvAVuS69yWnOclERBOwA48BvzXBMtzTC5fdLkGgVJKJc3JsZKUUmouGutZSTN/1Del\nlFJTSoNBKaWUhQaDUkopCw0GpZRSFhoMSimlLDQYlFJKWWgwKKWUstBgUEopZaHBoJRSykKDQSml\nlIUGg1JKKQsNBqWUUhYaDEoppSw0GJRSSlloMCillLLQYFBKKWWhwaCUUspCg0EppZSFBoNSSikL\nDQallFIWGgxKKaUsNBiUUkpZaDAopZSy0GBQSillocGglFLKQoNBKaWUhQaDUkopiwkFg4iUisgf\nReSt5HNJjuUeFpFOEfndiPkLRWRzcv1fiIhrIuVRSik1cRM9YvgK8LgxZjHweHI6m38BPp5l/j8B\ntyfX7wBunGB5lFJKTdBEg+Fq4MfJ1z8G3p9tIWPM40BP+jwREeAi4P7R1ldKKTV1JhoMlcaYZoDk\nc8U41g0CncaYaHK6AaidYHmUUkpNkGO0BUTkMaAqy1tfm+B3S5Z55ijluAm4CaC+vn6CX62UUiqX\nUYPBGHNJrvdE5JCIVBtjmkWkGjg8ju9uBYpFxJE8aqgDmo5SjvXAeoCVK1fmDBCllFITM9GmpAeB\nG5KvbwAeGOuKxhgDPAlceyzrK6WUmhwTDYZvA5eKyFvApclpRGSliPxgaCEReRb4JXCxiDSIyHuS\nb30ZuFVEdpPoc/jhBMujlFJqgkZtSjoaY0wbcHGW+S8Cn0qbPj/H+nuBsydSBqWUUvmlVz4rpZSy\n0GBQSillocGglFLKQoNBKaWUhQaDUkopCw0GpZRSFhoMSimlLDQYlFJKWWgwKKWUstBgUEopZaHB\noJRSykKDQSmllIUGg1JKKQsNBqWUUhYaDEoppSw0GJRSSlloMCillLLQYFBKKWWhwaCUUspCg0Ep\npZSFBoNSSikLDQallFIWGgxKKaUsNBiUUkpZaDAopZSy0GBQSillocGglFLKQoNBKaWUxYSCQURK\nReSPIvJW8rkkx3IPi0iniPxuxPz/EZF9IrIt+Vg+kfIopZSauIkeMXwFeNwYsxh4PDmdzb8AH8/x\n3peMMcuTj20TLI9SSqkJmmgwXA38OPn6x8D7sy1kjHkc6JngdymllJoCjgmuX2mMaQYwxjSLSMUx\nfMbfi8jfkDziMMYMZFtIRG4CbkpOhkRk1zGVePqVAa3TXYhpoNs9t+h2z0zzx7LQqMEgIo8BVVne\n+tp4S5TFV4EWwAWsB74M3JZtQWPM+uQys5qIvGiMWTnd5Zhqut1zi2737DZqMBhjLsn1nogcEpHq\n5NFCNXB4PF8+dLQBDIjIfwNfHM/6Siml8m+ifQwPAjckX98APDCelZNhgogIif6J1ydYHqWUUhM0\n0WD4NnCpiLwFXJqcRkRWisgPhhYSkWeBXwIXi0iDiLwn+dbPROQ14DUSbXN/N8HyzAazvjnsGOl2\nzy263bOYGGOmuwxKKaVmEL3yWSmllIUGg1JKKQsNhkk21mFDksv6RaRRRL43lWWcDGPZbhFZLiKb\nRGS7iLwqIh+ZjrLmg4hcJiK7RGS3iGSMACAibhH5RfL9zSKyYOpLmV9j2OZbRWRH8t/2cREZ0zn0\nM91o25223LUiYkRk1p2+qsEw+cY6bAjAt4Cnp6RUk28s290HfMIYcypwGfAdESmewjLmhYjYgTuB\n9wJLgetEZOmIxW4EOowxJwK3A/80taXMrzFu88vASmPM6cD9wD9PbSnzb4zbjYgUAZ8DNk9tCfND\ng2HyjWnYEBFZAVQCj05RuSbbqNttjHnTGPNW8nUTietgyqeshPlzNrDbGLPXGBMB7iWx/enSf4/7\nSZyhJ1NYxnwbdZuNMU8aY/qSky8AdVNcxskwln9rSOzk/TMQnsrC5YsGw+SzDBsCZAwbIiI24N+A\nL01x2SbTqNudTkTOJnEF/J4pKFu+1QIH06YbkvOyLmOMiQJdQHBKSjc5xrLN6W4E/jCpJZoao263\niJwBzDPGWEaTnk0mOlaSIi/DhnwG+L0x5uBs2onM13ApyQsdfwrcYIyJ56NsUyzbP9rI88DHssxs\nMubtEZGPASuBCye1RFPjqNud3Mm7HfjkVBVoMmgw5EEehg1ZBZwvIp8BfIBLRELGmKP1R0y7fAyX\nIiJ+YAPwdWPMC5NU1MnWAMxLm64DmnIs0yAiDiAAtE9N8SbFWLYZEbmExI7ChbkGyJxlRtvuImAZ\n8FRyJ68KeFBErjLGvDhlpZwgbUqafKMOG2KM+agxpt4Ys4DEeFE/memhMAajbreIuIDfkNjeX05h\n2fJtC7BYRBYmt2ktie1Pl/57XAs8YWb31aWjbnOySeVu4CpjzLjGUZvBjrrdxpguY0yZMWZB8u/5\nBRLbP2tCATQYpsKYhg05Do1luz8MXAB8cjbfxS/ZZ3Az8AiwE7jPGLNdRG4TkauSi/0QCIrIbuBW\njn522ow3xm3+FxJHwL9M/tuODMtZZ4zbPevpkBhKKaUs9IhBKaWUhQaDUkopCw0GpZRSFhoMSiml\nLDQYlFJKWWgwKKWUstBgUEopZfH/AUNn8FKPNPzGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11132ff90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(cv_x(trajectory), cv_y(trajectory), '-o', lw=2, markersize=8)\n",
    "for change in changes:\n",
    "    plt.plot(cv_x(change.trials[0].trajectory), cv_y(change.trials[0].trajectory), '-o')\n",
    "plt.ylim(-0.15, 0.15);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the magnitude of the change is rather small, and because the total kinetic energy in the direction of the transition is large in the original trajectory, note that nearly all the paths make a full crossing. Now let's compare that with the thermal randomizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "thermal_mover = paths.TwoWayShootingMover(ensemble=tps_ensemble,\n",
    "                                          selector=paths.UniformSelector(),\n",
    "                                          modifier=thermal_randomizer,\n",
    "                                          engine=engine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcXGWd7/HPr9be09VL9j0Ewioh\nLTGGIIhs44I4CGEEMooTnJE7ztXxjuK4DDgMXvWq85I7Q66KoCKD2xBEQQg7hiXshDV7muy9pfeu\n5bl/nOpOne6q9FK9Jt/361WvqjrnOXWeOumcb53zPOc55pxDRESkW2CsKyAiIuOLgkFERHwUDCIi\n4qNgEBERHwWDiIj4KBhERMRnWILBzC4wszfNbJOZfSnL/DPN7HkzS5jZJb3mJc3sxfRj7XDUR0RE\nhs7yvY7BzILAW8C5QC3wLHC5c+61jDJzgTLgH4G1zrlfZ8xrcc6V5FUJEREZNqFh+IzTgU3OuS0A\nZnYncBHQEwzOuW3pealhWJ+IiIyg4QiGGcDOjPe1wNJBLF9gZhuABHCTc+6/sxUys9XAaoDi4uIl\nixYtGmJ1RUSOTs8999wB51x1f+WGIxgsy7TBnJ+a7ZzbZWbzgYfM7BXn3OY+H+jcGmANQE1Njduw\nYcPQaisicpQys+0DKTccjc+1wKyM9zOBXQNd2Dm3K/28BXgEWDwMdRIRkSEajmB4FlhoZvPMLAKs\nBAbUu8jMYmYWTb+uApaT0TYhIiKjL+9gcM4lgGuB+4HXgbuccxvN7Hoz+wiAmb3bzGqBjwO3mNnG\n9OLHAxvM7CXgYbw2BgWDiMgYyru76lhQG4OIyOCZ2XPOuZr+yunKZxER8VEwiIiIj4JBRER8FAwi\nIuKjYBARER8Fg4iI+CgYRETER8EgIiI+CgYREfFRMIiIiI+CQUREfBQMIiLio2AQEREfBYOIiPgo\nGERExEfBICIiPgoGERHxUTCIiIiPgkFERHwUDCIi4qNgEBERHwWDiIj4KBhERMRHwSAiIj4KBhER\n8VEwiIiIj4JBRER8hiUYzOwCM3vTzDaZ2ZeyzD/TzJ43s4SZXdJr3iozezv9WDUc9RERkaHLOxjM\nLAjcDFwInABcbmYn9Cq2A/hr4I5ey1YAXweWAqcDXzezWL51EhGRoRuOI4bTgU3OuS3OuS7gTuCi\nzALOuW3OuZeBVK9lzwcecM7VO+cagAeAC4ahTiIiMkTDEQwzgJ0Z72vT00Z6WRERGQHDEQyWZZob\n7mXNbLWZbTCzDfv37x9w5UREZHCGIxhqgVkZ72cCu4Z7WefcGudcjXOuprq6ekgVFRGR/g1HMDwL\nLDSzeWYWAVYCawe47P3AeWYWSzc6n5eeJiIiYyTvYHDOJYBr8XborwN3Oec2mtn1ZvYRADN7t5nV\nAh8HbjGzjell64Eb8MLlWeD69DQRERkj5txAmwPGj5qaGrdhw4axroaIyIRiZs8552r6K6crn0VE\nxEfBICIiPgoGERHxUTCIiIiPgkFERHwUDCIi4qNgEBERHwWDiIj4KBhERMRHwSAiIj4KBhER8VEw\niIiIj4JBRER8FAwiIuKjYBARER8Fg4iI+CgYRETER8EgIiI+CgYREfFRMIiIiI+CQUREfBQMIiLi\no2AQEREfBYOIiPiExroCIiJyeG89vYf1d2+mpb6Tkoooyy5awLFLp47Y+hQMIiLj2FtP7+HhX7xB\noisFQEt9Jw//4g2AEQsHBYOIyBhLpRwdLXE6WuK0t3Sln+N0tHTx/P07ekKhW6Irxfq7NysYREQm\ninhXkvbmjB18cxftPTv7uDevNU57sxcEnW0JcINbR0t958hUnmEKBjO7APgBEAR+5Jy7qdf8KHA7\nsASoAy5zzm0zs7nA68Cb6aJPOec+Mxx1EpGjy0idh3cpR0db9w494xd9c99f+O0tXXQ0x0nEU1k/\nywJGQUmYwvSjckYJhaXh9LQIhSVhCkq750coKA7z86+vzxoCJRXRvL9bLnkHg5kFgZuBc4Fa4Fkz\nW+ucey2j2NVAg3PuGDNbCXwLuCw9b7Nz7tR86yEiR6/BnIdPdCUP/XJv6fLt4H2/6NM7+87WOC7H\nr/lwNOjt2IvDFJVFqJhe7O3cS8IUlqZ39N07/JIw0aIQZjao77bsogW+7wYQigRYdtGCQX3OYAzH\nEcPpwCbn3BYAM7sTuAjIDIaLgG+kX/8a+KENduuIiOSw/u7NWc/DP/yLN3jzmT2Hfu23xkl0JrN+\nhhkUZOzIK6YVU1B6aKfe8yu+JNzzKz8UDo74d+sOtonWK2kGsDPjfS2wNFcZ51zCzJqAyvS8eWb2\nAnAQ+Gfn3OPZVmJmq4HVALNnzx6GaovIRNPREqdxXxsNe9po3Os9Gva25TzfnuhK0d4cp7A0TGxq\n8aGderH3iz5zhx8tCmGB8fl79dilU0c0CHobjmDItiV7H3jlKrMbmO2cqzOzJcB/m9mJzrmDfQo7\ntwZYA1BTUzPIZhoRmSiSiRRN+9t9O/7GdBB0tMZ7ygUCxqTJhUyaXERLQwfxjr5HAiUVUS697t2j\nWf0jwnAEQy0wK+P9TGBXjjK1ZhYCJgH1zjkHdAI4554zs83AscCGYaiXiIyiwTT+OudoO9hF4570\njn9f+ghgTxsH6zpwqUO//QrLIsSmFDF/cTXlU4qITSmifEoRZVUFBIKBnnWP9nn4I9lwBMOzwEIz\nmwe8A6wE/qpXmbXAKmA9cAnwkHPOmVk1XkAkzWw+sBDYMgx1ysvrjz/M43feTnPdAUorq1ix8iqO\nX3H2WFdLZNzK1fibSKSYPKfUd+qn+9GV8Qs/GA5QPrmIqlmlLHz3FMqnFFE+uYjyqUVEC/vfTY3F\nefgjWd7BkG4zuBa4H6+76k+ccxvN7Hpgg3NuLfBj4GdmtgmoxwsPgDOB680sASSBzzjn6vOtUz5e\nf/xh/rTmhyS6vHOWzQf286c1PwRQOIjkkLPx92dv+KaVxKKUTyniuKVTKZ/q/fIvn1JEaawg7/P7\no30e/kg2LNcxOOf+APyh17SvZbzuAD6eZbnfAL8ZjjoMl8fvvL0nFLolujp57I6fsmj5+7CAxh0U\nAUglU+zdepAdr9Uf9mKr8z59IrGpRUyaXEQ4MvK9eCR/uvK5l+a6A1mnt9TX8f0rLqa4vIKSigpK\nKiq9R6yS0u7X6Uc4WjDKtRYZHc31Hex8rZ4dG+vY+UYDXe0JzCAQMlKJvn1CSiqiLKyZMgY1lXwo\nGHoprayi+cD+PtOjxSW869wLaamvo6W+jrqdO9j+8gt0tbdnKVtMSexQUPiCIz29qGySjj5kXDhc\no3EinmTX243seK2eHRvradjdCninhI45rZpZJ1Qyc1GMHa/WqfH3CKJg6GXFyqt8bQwAoUiUcz55\nTdY2hq72NprTYdHzaKjLCJDttDY24pz//GsgGKI4FvOCIx0WxbEKHX3IqMrWaPzQz99g+8Y6OloT\n7HqrgUQ8RSBkzFhYzgnLpzHrhAoqphX7ruBV4++RxVyua73HsZqaGrdhw8j1aB3uXkmpZJLWpgZ/\neGSESHN9PS31dcQ7dPQho+u2657M2T5QPqWI2SdUMOuECmYcGyMcVfvARGdmzznnavorpyOGLI5f\ncfaw9kAKBIOUVlRRWlF12HKdbW19jzqGePSR7RGODH7QLXXdPbIdrtH4E//ynlGsiYwnCoZxJFpU\nRLSoiMqZs3KW6e/oY//O7Wx96fmsRx8FxSV9A6PX0UhhaVnP0Ye67h659u9o5pnfb805fyRH7pTx\nT8EwweRz9JHZFrJ/xzZaGxvoPWxkIBjyel3FKtm3bUvWrruP33m7gmECyNaoXD61iGd/v5Vtr9QR\nLQoxf3EV21+tJxlXo7EcojaGo1gqmaS18dDRR3OvU1c7N76cc9lTPnABk+fOp3rOfKpnzyVcoEby\n8STbEBFm3u+AaFGIUz8wi5PPnkW0MDTq9xOWsTPQNgYFg+S05rOfzNp1NxgOE45E6Wht8SaYEZs6\nneq585k8Z54XGHPnU1weG/TY8zI8cjUqRwpDrLrxvUQGMMyEHHnU+Cx5y9V197zV17LojLNortvP\nvm1b2b9tC/u2bWHv5rd4a/2hUdOLJpVTPWce1emwmDx3PrFpMwgE1btlpOVqVO5qTygUpF/6C5Gc\nutsRcvVKKquaTFnVZI6pOXT7jc62VvZv3+oFxnYvMF7441qSiQQAoXCEqtlz0kcX3pFF9ew5RAqL\nRv8LHqGcc0QKQ3S1J/rMU6OyDIROJcmISyYS1O+q7Tmy8AJjKx0tzV4BM8qnTO0JCu9U1DxKYpU6\nFTUAvjaCWJSSWJQ9Ww5iAcjs2RyKBDj7E4vUfnAUUxuDjGvOOVrq67yg2LaFfdu3sH/bVhr37u4p\nU1ha5h1RdJ+KmjOP2PSZBEOHDnSP9usssjUyAyxYUs28U6p46u4talSWHmpjkHHNzCitrKK0sooF\nS07vmd7Z1sb+Hd3tFlvZv30rL97/e5Jx785dwXCYqllzqJ4zn1QywZvrH++ZdzReZ5FtuGuAvVsP\ncsHfnMxxS6eNQa1kolMwyLgSLSpi5qITmbnoxJ5pqWTy0Kmo7VvZt20Lmzc8RXtznzvAHnXXWeRq\nZD7cFc0i/VEwyLgXCAapmjWHqllzenb4zjn+z8oPZy3ffGA/2156njknn3rEjyFVEovS0tA3BNTI\nLPlQMMiEZGaUVlVnvc4CM35z49coq57MSWefy0lnnUtp5eGvFJ+IUsmU1/W0VzDoymXJl4JBJqyc\nQ6Rf/beEIhFeWXc/f77rF6z/1S+Zt3gJJ7//fOaf9u4Jfx1FZi8kgOkLJ3GwrkONzDJs1CtJJrT+\neiU17t3Dqw//iVcfeZDWhnqKYxWc+L5zOPns8yifOvEaZrP1QlI3VBkodVcVyZBKJtnywgZeWXcf\nW194DudSzD7pFE5+//kcc/p7CYXDE6Lra66hLkoqoqy6cfkY1EgmEnVXFckQCAY5pmYpx9Qspbn+\nABsffpBXHn6Ae//92xSUlDJlwUJqX3uVZLwLGL9dX9ULSUaDgkGOOqUVVbznL1ey9OJL2f7qS7yy\n7n7eeuqJPuXGW9dXl3IEQwGSib7XLagXkgwnBYMctSwQYO4pi5l7ymK+e9mHspZprjswyrXqq3dj\nc7ahLtQLSYbTkd3JW2SASquqs08f426u3Y3NvlNFBgXF3m+6koqoGp5l2OmIQYTsXV8xY/mlV4xd\npcg+5IVLQiga5LPfPXOMaiVHOgWDCH2HGC8sKaW9+SBN+/aMab3U2CxjQcFwhGvtTHDLY5v5+frt\nNLTFiRWFuWLZHK45cwHF0Yn1zz/S3+X4FWf7Gpr/8MPv8tRv/4v18cn8YpMb0joHU+feZacUhLnC\nQliWHuVNluK06//U81nAkNczkt9puEz0v+N7t9zLD57/AXta9zClaCrHRj7O+pfmjNvvMizXMZjZ\nBcAPgCDwI+fcTb3mR4HbgSVAHXCZc25bet6XgauBJPD3zrn7+1ufrmMYmNbOBBf/3yfZXtdGZ0ZP\nlmgowJzKIn73d8vHzR9if8biu9TXN/Ifn/sMiUSKlAUoTbbQHCzh2aplJOad2u86B1Pn7rLRdzpY\n1hqkzBlJvEZAZ0Yw479pHMd9hXHeiCaJhgLMihWCwc769gGvZ6jbsXv52S88zhWv3Et1eyP7C8v5\n+ckfZMfiFb7lm+65h33f+z6J3bsJTZvG5P/5D0z6cPbxrQ6ne52VXT9hX+UG9oeM6oRjcl0NdZFP\nZa3zvVvu5aZnbqKxsxGASZFJfHnpl/ng/A8Oev35unfLvXzjz9+gI9nRM82lwnTs/hiJg4uB0fs/\nOWrXMZhZELgZOBeoBZ41s7XOudcyil0NNDjnjjGzlcC3gMvM7ARgJXAiMB140MyOdc4l862XeL8g\ne+8AADoTKbbXtXHLY5v5/LnHjVHtBmcsvstPn9vL60XHsrjx+Z5pZckWztz3MI8Btzw29bDrHEyd\nb3lsM9F3OjinJUQY7+ZEISCBY2MkwXGpENG446A5HitI8EY02fNZWw+0ghnJlBvQevLZjrc8tpnZ\nLzzOZ5+7i4KkN9z5lPZGPvvcXdycsU2a7rmH3V/9Gq7D2xkmdu1i91e/BjDocLjlsc1Udv2ETZM3\n0JEeFHFf2Dg4eQMXrT/Am2c1U9h4oCd8njgxwFef/CrxVLznM5q6mvjnJ/6ZjkQH5809j5RL9Twc\njmQqicP5pqdcihQpnHMkXRLnXM55Kee9TpEilfLmdU/71jPf8oUCgAXifKDhl1z+0O+oam3rCdf+\n/qZGS95HDGa2DPiGc+789PsvAzjn/i2jzP3pMuvNLATsAaqBL2WWzSx3uHXqiGFgTrv+T9S3xQ9b\nJhgY3B3ShuN+akO5KVs82f/faSSYpZPdANeVrVhnIsWqHT+jLNnSZ97BYAm3zb6SwrA37lLmd+p+\n2dp1+N83BpQUhDCguSPB3zRFmeT6focmS7Fm0tDbFAwoKwx79W6Pc7gtaUBpQe7fi80dCW69/5tM\naW/sM68tGOW+uUshGuX8tx+nON7Rp8zBSDG3LvlLUhYgZZZ+zvE64L1u6XJMm/3vHIgESBmkApAy\nOG1TissfdUQz7mDaGYLbzjGeXuSVBa8sBg5wlvHI8n5If5xDsHxjkmv+4CjIqHtHMMyPl67khz+5\nbsTWO5pXPs8Adma8rwWW5irjnEuYWRNQmZ7+VK9lZ2RbiZmtBlYDzJ49exiqfeRr6CcUgD6/Miey\nrmTfC7/yVZolFDKnt8eHfnDr8Ha03cpc9p1SrumDWU9Te/9/C91lD3b0vVd0puosoQBQmOzkwm1P\nEUnGCeaIn7KuVj63/vYB1aWv/rd1NAGr73esvn/o/y4pAKNPsOQMGvoPnN6vq5sg1OvPtSAZ55IX\n7gFGLhgGajiCIdtfbe+/ilxlBrKsN9G5NcAa8I4YBlPBo1WsKHzYI4aK4gjPXHfOgD9vODb6UA9Q\nl9744GGDrqI4wvovv39Y1tVt2Y0PkrAQYdd3R9kcLCFWFObJL73ft57MVa741kOHrXOsKMwj/3g2\nia4kX7nhCSzH4c1BcxhD3/6xojAPfeEszOCs7zxC4wDqlMv7vvMw+wvLsx4x7Css539e/C88/IWz\n2PvB80nt6dujK1BdTeV/rIFkElIpXDIFqSQkU7iUN41kEpd+JpXiut+8SFnsxzQHjEAKAs57XHtP\nKucO5NZzA5jDe3DouThYzIfmfgTnnLc+53CpFM6lvHWmn51LZczznsG/DM7hnPO+R3d5l4LUoWdc\niqaOBna3voNzjoDzKji1Ifu/Zq7QHW3DEQy1wKyM9zOBXTnK1KZPJU0C6ge4rAzRFcvmcMujW/qc\nTwavseuK98wmlO30yzh05QC+SzQ0fMNpO+e4ghfAJUgSIMih9cYtxLNVy7hy2RyKIrn/C2Wr86LO\nIGd2hLyjgM4Qbz9Sy1vP7OXk1gDbQymmJ6ynjQG8hub1xUnePS/GSzubsn7/oJG1jQG8bXPlsjnE\niiMAXNXPdrxy2RwmFYVzfqerls3h55s+6GtjAO80yM9P/iBXLptDeXEE+8LnfW0MAFZQwNT/9UUm\nnbQo5+dnc2JgJs8+/7KvjQHg8kdTVPW9iR+JyeWsO73T18YAELIQ3zzj67xrjBqgr3viOlLpS9YX\n1SaozlL3zlj2Cy1H23DsFZ4FFprZPDOL4DUmr+1VZi2wKv36EuAh5zVurAVWmlnUzOYBC4FnhqFO\nAlxz5gLmVBYRDfn/mbt7QHR3c5wIRvu7PP27u+C1J9g2dQmPTH4/B4Ml3mmWYAmPTT6bxLxT+11n\n7zov6gxyQXuYSS6AYVh7kqfXbqW1qZPz/vZkXlwQZl1JgibzGkSbLMW6kgSdMwr44eWn5fz+86qK\nmVc1sG2T73a85swF7Fi8gpuXXMrewnJSwN7Ccm5ecik7Fq/oWX7Shz/MtBuuJzR9OpgRmj6daTdc\nP6ReSdecuYC6yKc4Zl8Nk+MpzDkmx1NsX1yIBf0BZwUFzPniddyw/AbKo+U90ydFJvHNM745Jr2S\nAD44/4PceMaNFAQLALjjLKOj12+KrlCEGV/8whjUrq/h6q76F8D38bqr/sQ5969mdj2wwTm31swK\ngJ8Bi/GOFFY657akl/0K8CkgAfyDc+6P/a1Pjc8D19P/+6kdNLR1ESuKcMV7Zo+rPtMDNVrf5aUH\n/siDP7qZE1aczYqr/wdrntg65HVm1vmS3cakVN/fYsWxKH/9b8v7/X6Hmw8MeNvkux3H4m8q6zqX\nzmLVk1+l/p4XSbSFCE2fPuQusaPl3i338v3nfsCetj0sf7mAyx9JUtXaRmesmhlf/AJTPnbRiK5f\n92MQGYSeezGkbxVaNWceV9z4PYKh4dnRxTuTrPncoznnf/Y/359znhxGKgm/uwZe+RVc+G1Yunqs\nazSu6X4MIgP0+uMP9xknqXH3O7y1/vG8h9yu29XCxsd28eZTu3OW0ZDZeQgE4aP/AV1t8McvQrQE\nTv2rsa7VhKdgkKNaV0c7D922xj94HpDo6hryvRiS8RSbX9jHq4+9w+5NTQRCxjGnTaasqoAXH9zZ\n57acGjI7T8EwXPIT+OVlcPdnIVwEJ350rGs1oSkY5KiTiMfZ9tLzvPHEI2x+7pk+odCtv3sxZN4n\noaQiyrveP5O2g3Fe//NuOlrilFUX8t6PHcOi906lsMTrFRSbUuxbZtlFCzRk9nAIF8DKO+BnF8Nv\nPg2RYlh47ljXasJSG4McFVKpJLWvvcobTz7KW08/SWdrK4WlZRy7bAVvP/0kbU19+4+XVlWz+uZb\ns35e930Seg+JDTB/cTUnrZjBzEUxbJBXlkue2hvhtg/Dgbfgit/A3DPGukbjitoY5KjnnGPv5rd5\n/clHeXP947Q21BMuKGThu9/DojPOYvZJ7yIYCjHj2EV92hhCkSgrVl7V5zNbGjrZvbmRR+54M2so\nFJdHufCak0f0e8lhFJbDlb+DW/8C7rgMrloLM5eMda0mHAWDHHHqanfyxp8f5Y0nH6Vxz26CoRDz\nFtewaPlZzD+thnC0wFe+970YSiurWLHyKo5bfhb7dzazZ3MTuzc3sWdzE831fcf/ydTaqPskjLni\nKrjqv+EnF8DPPwaf/ANMOXGsazWh6FSSTCg93UozduDHrzibgwf28+afH+P1Jx9l/7YtmAWYdeLJ\nLDrjfSw8/b0UFJf0+9ld7Qn2bD0UAnu3HiTe6Y25UzQpwrQFk5i2oJypCyZx3y2v0NLQNwRKKqKs\nunH5sH9vGYKGbV44pJLwqfugUo38uo5BjjjZupUGgiHKJk+hcfc7AEw75jgWLT+TY5etoCRWAfRt\nJF520QIWnj6F5rqOnhDYvbmJul0t4LwBNitmlKSDYBJT50+itLIAyxh5M1sbQygS0P2Xx5v9b8Kt\nF3o9lT75Ryif1f8yRzAFgxxxbvnbVbTU1/WZHggGWfaXl7No+fsonzrNNy/bDtwMQgVB4u3e0UC4\nIMjUeWVMne8dEUyZV0aksP+zrNkCR6EwDu1+CX76YQhFIRCC5t0waSac8zU45dKxrt2oUuOzTFgd\nrS3U7dzBgZ3bObBzO3U7t3OgdgftB5uylk+lUiy9+DJam7rYvamRg3UdNNe1c7Cug7ee3kuy14Bx\nzoFLOs5ceSzTjplExfQSAkPoPXTs0qkKgolg2rvgPZ+BR791aFrTTrjn773XR1k4DISCQcZMvKOD\nutp0ANTu8AJg53bfUUG4oJCqWbNZsOR03n76z3S2tfb5nECwlP/8+0dIJfxHv0VlkT6h0C3RleLk\ns2YO7xeS8evFO/pOi7fDuusVDFkoGGTEJeJxGnbV9hwBHNi5nbraHTTtPTRefzAcJjZtJtVzjmfW\nSVMJRapwVNDZVkhzfSdbX+0gSQh4AG+8RY9ZmOnHXcisk2ZRVlVIaWUBZZUFlFYUEIoEue26J2mp\nz95ILEeRptoc03fCrhdg+uLRrc84p2CQvGSeZy+OhTn5zBJKY60c2HnoCKBhz670jU7AAkFKKqZQ\nWDqd6YtOxVFBIl5O68ECWpqhpfnQZxeUJCmrTFA5o4R5p1RRWnkMDbtn8dqjv6W1oY7SqkO9knJZ\ndtGCrI3EGobiKDNpphcC2aw5C6aeAktWwckfh4JJo1q18UiNzzJoLpXi4IF9vPTgi7y07iUSXftx\nyTpcsp5Dt180osWVhAuqwSqJx8txrgILxDDzbqgTLQ5RVlno/cKvLKC0+3WV94s/cph7Dw+GGomF\nl+/y2hTi7YemhQvh/JsgFYfnboO9r0CoEE682AuJWUtH7R7Qo0W9kmTAel8bsPzSK5m/5Aw6WuI0\n7N7Hvu3bqH9nB417ammue4e2xj2kkl2HPiBQSiBQiQWrCASrsGAlFqygoLj71E76uSpj519RMKCe\nPyLD5uW7vDaFptq+vZKc804pPX8bvPJr6GqB6kVw2lVwykoorhzbug8TBYOQSqbobEvQ0Rrvee5o\njdPZmp7WGmf3pmd457Xf4ny3QTSwSUA7uIzz81ZEIFhFKFpNtGQKnW1lXghY3/P1n/7emUS145eJ\nqLMFNv7WO4p4ZwMEI7DoQ95RxNwzITAxboebjbqrHkGSyRSdrQk62+J0tCboTO/gM3f4h6ZllGuL\ng2vHpVpxrgWXaoWM12atJLv20Pc28w6zFmaftJzY9NlUz57L5HlzKJ9SSaQw1NO183ANuwoFmbCi\nJd6RwmlXwd6NXkC8fKcXFrF5cNqVcOonoPTIPR2pI4YhGOo56+4dfM+OvC1BR0s8vSNP/5Jvi6en\nHSrX1ZH0fY5zrmeHj2shGO4gEGzDaAPXSirRTKKrmXhnMy6V7FOPaHEJJeUxiisq2fHKi9kra8YX\n7ryn3+2gq3/lqBBvh9fv8UJi+xNgQTjuQjhtFRxzjnfDoAlARwwjpPfOsKW+k3U/e5133m6gfEpx\nzw7f96s+vcOPd/TdSXczg0hRiEhBnHC4HQu2ES1oJxJpIZVsIdl1kK72g3S2NdHZ0pR1h19QXEJx\nrIKSiiqKyxdSEqvw3scqKC6voKSigqLyGOHIoVM/az77yZ7bWWYqrazqd1t07/zVsCtHvHCh1x5x\nyqVwYJPXFvHiHfDG76FsBiwGNnp5AAAKzklEQVS+AhZfecQMuaEjhkHKdfqkmwWMguIQ0aIwBcVh\nCopDRIqChIJdYK24VCupZAuJroPEO7wdfUdzI21NDbQ2NpBKZtnhl5RSXB6jpKIyvZOPURyrpCSW\n8VxeQSgSGfT3yTb+UCgS5bzV1+Z9W0uRI1qiC976o3cUsfkhb9ox53hHEcdd6N1ZbpzREcMIyQyF\nQ6d0WnCuhRUfn0FHSyOtjQ20NNTT2ljPgW31tB1mh9/9q75q1myKM37ZF5dX9ITAUHb4A5VryGmF\ngkg/QhE44SLv0bAdXvi597jrSiie7N17+rSrJuSorkfNEUOu4ZqzcakUbQebenbuLfXec2tDPa8/\nuYl458H0+f02oO+QCwWlZd45/FgFJbFKimOxPqd1RnqHLyJjIJmATQ/C87fDW/eBS8LcFd5RxPEf\n9m5BOobUXTVDttMlwVCY41ecTVlVtbfzb/B2/C0N9bQ1NWb/hV9aRjhaStvBMFAMgWLMiglFSlly\n4YmccMZCispjhMLj7xBSREbZwd3w4i+8kGjcDoUx75qIJatg8vFjUiUFQ4ZcDazdCkvLDv2az3wu\nP/Q+c4evK2lFZMBSKdj6qBcQr9/jXWk9893eUcRJH4NI8ahVRcGQ4bsrP+xd2ZjFP/zidwRD+oUv\nIqOgtQ5e+qXXq+nAWxAphZMv8Y4iRmEgPzU+ZyitrMreJbOqWqEgIqOnuBLeey0s+yzseMoLiJfu\nhOduPTSQXyAEj30n+9Ado2TiXts9CCtWXkUo4h+2IRSJsmLlVWNUIxE5qpnBnGVw8X/CF96Av/gO\n4ODeL8A9n0uPBOsO3VDo5btGtXp5BYOZVZjZA2b2dvo5lqPcqnSZt81sVcb0R8zsTTN7Mf2YnE99\ncjl+xdmct/paSquqwYzSqmr10xeR8aGwHE7/G7jmcSiZ0nd+9w2FRlG+p5K+BKxzzt1kZl9Kv/+n\nzAJmVgF8HajBG5TnOTNb65xrSBf5hHNuxK9WO37F2QoCERm/zKBlX/Z5uW40NELyPZV0EXBb+vVt\nwEezlDkfeMA5V58OgweAC/Jcr4jIkWdSjtvN5po+QvINhinOud0A6edsp4JmAJm3TqpNT+t2a/o0\n0lfNct8Vw8xWm9kGM9uwf3/urqciIhPWOV/zxmXKFC70po+ifk8lmdmDQLZO+l8Z4Dqy7ey7+45+\nwjn3jpmVAr8BrgRuz/Yhzrk1wBrwuqsOcN0iIhNHd++jXDcUGiX9BoNz7gO55pnZXjOb5pzbbWbT\ngGwnyGqBszLezwQeSX/2O+nnZjO7AzidHMEgInJU6B7FdQzleyppLdDdy2gVcHeWMvcD55lZLN1r\n6TzgfjMLmVkVgJmFgQ8Br+ZZHxERyVO+wXATcK6ZvQ2cm36PmdWY2Y8AnHP1wA3As+nH9elpUbyA\neBl4EXgH+H951kdERPJ0VAyJISIiAx8S46i48llERAZOwSAiIj4KBhER8VEwiIiIj4JBRER8FAwi\nIuKjYBARER8Fg4iI+CgYRETER8EgIiI+CgYREfFRMIiIiI+CQUREfBQMIiLio2AQEREfBYOIiPgo\nGERExEfBICIiPgoGERHxUTCIiIiPgkFERHwUDCIi4qNgEBERHwWDiIj4KBhERMRHwSAiIj4KBhER\n8ckrGMyswsweMLO308+xHOXuM7NGM/t9r+nzzOzp9PL/ZWaRfOojIiL5y/eI4UvAOufcQmBd+n02\n3wauzDL9W8D30ss3AFfnWR8REclTvsFwEXBb+vVtwEezFXLOrQOaM6eZmQHvB37d3/IiIjJ68g2G\nKc653QDp58mDWLYSaHTOJdLva4EZedZHRETyFOqvgJk9CEzNMusrea7bskxzh6nHamA1wOzZs/Nc\ntYiI5NJvMDjnPpBrnpntNbNpzrndZjYN2DeIdR8Ays0slD5qmAnsOkw91gBrAGpqanIGiIiI5Cff\nU0lrgVXp16uAuwe6oHPOAQ8DlwxleRERGRn5BsNNwLlm9jZwbvo9ZlZjZj/qLmRmjwO/As4xs1oz\nOz8965+Az5vZJrw2hx/nWR8REclTv6eSDsc5Vweck2X6BuDTGe9X5Fh+C3B6PnUQEZHhpSufRUTE\nR8EgIiI+CgYREfFRMIiIiI+CQUREfBQMIiLio2AQEREfBYOIiPgoGERExEfBICIiPgoGERHxUTCI\niIiPgkFERHwUDCIi4qNgEBERHwWDiIj4KBhERMRHwSAiIj4KBhER8VEwiIiIj4JBRER8FAwiIuKj\nYBARER8Fg4iI+CgYRETER8EgIiI+CgYREfFRMIiIiE9ewWBmFWb2gJm9nX6O5Sh3n5k1mtnve03/\nqZltNbMX049T86mPiIjkL98jhi8B65xzC4F16ffZfBu4Mse8LzrnTk0/XsyzPiIikqd8g+Ei4Lb0\n69uAj2Yr5JxbBzTnuS4RERkFoTyXn+Kc2w3gnNttZpOH8Bn/amZfI33E4ZzrzFbIzFYDq9NvW8zs\nzSHVeGxVAQfGuhLjnLbRwGg79U/bqK85AynUbzCY2YPA1CyzvjLYGmXxZWAPEAHWAP8EXJ+toHNu\nTbrMhGVmG5xzNWNdj/FM22hgtJ36p200dP0Gg3PuA7nmmdleM5uWPlqYBuwbzMq7jzaATjO7FfjH\nwSwvIiLDL982hrXAqvTrVcDdg1k4HSaYmeG1T7yaZ31ERCRP+QbDTcC5ZvY2cG76PWZWY2Y/6i5k\nZo8DvwLOMbNaMzs/PesXZvYK8Are+cBv5lmf8W5CnwobJdpGA6Pt1D9toyEy59xY10FERMYRXfks\nIiI+CgYREfFRMIyggQ4Zki5bZmbvmNkPR7OOY20g28jMTjWz9Wa20cxeNrPLxqKuo83MLjCzN81s\nk5n1GVXAzKJm9l/p+U+b2dzRr+XYGsA2+ryZvZb+u1lnZgPqx3+0UzCMrIEOGQJwA/DoqNRqfBnI\nNmoDrnLOnQhcAHzfzMpHsY6jzsyCwM3AhcAJwOVmdkKvYlcDDc65Y4DvAd8a3VqOrQFuoxeAGufc\nKcCvgf89urWcmBQMI2tAQ4aY2RJgCvCnUarXeNLvNnLOveWcezv9ehfe9TLVo1bDsXE6sMk5t8U5\n1wXcibetMmVuu1/j9fqzUazjWOt3GznnHnbOtaXfPgXMHOU6TkgKhpHlGzIE6DNkiJkFgO8CXxzl\nuo0X/W6jTGZ2Ot6V8ptHoW5jaQawM+N9bXpa1jLOuQTQBFSOSu3Gh4Fso0xXA38c0RodIfIdK+mo\nNwxDhvwd8Afn3M4j9cfecA2rkr4g8mfAKudcajjqNo5l+2Po3bd8IGWOZAP+/mZ2BVADvG9Ea3SE\nUDDkaRiGDFkGrDCzvwNKgIiZtTjnDtceMaEMx7AqZlYG3Av8s3PuqRGq6nhSC8zKeD8T2JWjTK2Z\nhYBJQP3oVG9cGMg2wsw+gPcj5H25BukUP51KGln9DhninPuEc262c24u3lhRtx9JoTAA/W4jM4sA\nv8PbNr8axbqNpWeBhWY2L/39V+Jtq0yZ2+4S4CF3dF2x2u82MrPFwC3AR5xzgxrL7WimYBhZAxoy\n5Cg3kG10KXAm8NdHy93+0m0G1wL3A68DdznnNprZ9Wb2kXSxHwOVZrYJ+DyH7/V2xBngNvo23pH4\nr9J/N73DVbLQkBgiIuKjIwYREfFRMIiIiI+CQUREfBQMIiLio2AQEREfBYOIiPgoGERExOf/AylW\nKvTvrMSgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e539050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "change = thermal_mover.move(initial_conditions)\n",
    "changes = [change] + [thermal_mover.move(initial_conditions) for i in range(4)]\n",
    "\n",
    "plt.plot(cv_x(trajectory), cv_y(trajectory), '-o', lw=2, markersize=8)\n",
    "for change in changes:\n",
    "    plt.plot(cv_x(change.trials[0].trajectory), cv_y(change.trials[0].trajectory), '-o')\n",
    "plt.ylim(-0.15, 0.15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#[len(change.trials[0].trajectory) for change in [changes[0]]]\n",
    "modified = changes[0].canonical.details.modified_shooting_snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "trial_trajectory = changes[3].trials[0].trajectory\n",
    "modified = changes[3].canonical.details.modified_shooting_snapshot\n",
    "index = trial_trajectory.index(modified)\n",
    "forward_part = trial_trajectory[index:]\n",
    "backward_part = trial_trajectory[:index]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, this doesn't even get a crossing. Here you'll often see paths that only go one direction: this is because the are forward-first shots that end in the left state, or backward-first shots that end in the right state. These trajectories could never satisfy the left-to-right ensemble. You also see far more left-to-left or right-to-right paths here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the two-way shooting move strategy\n",
    "\n",
    "As with most path movers in OpenPathSampling, the preferred way to select two-way shooting is to use its associated move strategy. This is also the easiest way. All you need to do is select an approach for the shooting point selection (`selector`) and for modifying the shooting point (`modifier`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy = paths.strategies.TwoWayShootingStrategy(modifier=thermal_randomizer,\n",
    "                                                   selector=paths.UniformSelector(),\n",
    "                                                   engine=engine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "scheme_thermal = paths.OneWayShootingMoveScheme(\n",
    "    network=network,\n",
    "    selector=paths.UniformSelector(),\n",
    "    engine=engine\n",
    ")\n",
    "scheme_thermal.append(strategy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "storage = paths.Storage(\"thermal.nc\", \"w\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on Monte Carlo cycle number 100\n",
      "Running for 10 seconds -  0.11 seconds per step\n",
      "Estimated time remaining: 0 seconds\n",
      "DONE! Completed 100 Monte Carlo cycles.\n"
     ]
    }
   ],
   "source": [
    "simulation = paths.PathSampling(storage=storage,\n",
    "                                move_scheme=scheme_thermal,\n",
    "                                sample_set=initial_conditions)\n",
    "simulation.run(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shooting ran 100.000% (expected 100.00%) of the cycles with acceptance 10/100 (10.00%)\n"
     ]
    }
   ],
   "source": [
    "scheme_thermal.move_summary(storage.steps)\n",
    "storage.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a pretty low acceptance probability. Let's try the same, but with the direction based modifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "scheme_direction = paths.OneWayShootingMoveScheme(\n",
    "    network=network,\n",
    "    selector=paths.UniformSelector(),\n",
    "    engine=engine\n",
    ")\n",
    "scheme_direction.append(\n",
    "    paths.strategies.TwoWayShootingStrategy(modifier=direction_modifier,\n",
    "                                            selector=paths.UniformSelector(),\n",
    "                                            engine=engine)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "storage = paths.Storage(\"direction.nc\", \"w\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on Monte Carlo cycle number 100\n",
      "Running for 14 seconds -  0.14 seconds per step\n",
      "Estimated time remaining: 0 seconds\n",
      "DONE! Completed 100 Monte Carlo cycles.\n"
     ]
    }
   ],
   "source": [
    "simulation = paths.PathSampling(storage=storage,\n",
    "                                move_scheme=scheme_direction,\n",
    "                                sample_set=initial_conditions)\n",
    "simulation.run(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shooting ran 100.000% (expected 100.00%) of the cycles with acceptance 80/100 (80.00%)\n"
     ]
    }
   ],
   "source": [
    "scheme_direction.move_summary(storage.steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This acceptance probability is actually too *high*. It indicates that each trajectory is probability too much like the previous trajectory. The acceptance can be decreased by increasing the delta-v.\n",
    "\n",
    "Why is the direction-based modifier better than thermal randomization for this system? Fundamentally, the issue is that this model is very ballistic. The thermostat friction is so low that this is effectively deterministic. If the motion is ballistic, then completely randomizing your velocity is likely to do a lot of harm to your chances of getting another crossing. On the other hand, if you have diffusive dynamics, then the system will largely lose memory of the direction of your velocity in some relatively short time frame. In that case, thermal randomization should be just as good as direction modification (and it works in more situations, e.g., with constraints.)"
   ]
  },
  {
   "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.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}