{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# CSNAnalysis Tutorial\n", "### A brief introduction to the use of the CSNAnalysis package\n", "---\n", "**Updated Aug 19, 2020**\n", "*Dickson Lab, Michigan State University*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "The CSNAnalysis package is a set of tools for network-based analysis of molecular dynamics trajectories.\n", " CSNAnalysis is an easy interface between enhanced sampling algorithms\n", " (e.g. WExplore implemented in `wepy`), molecular clustering programs (e.g. `MSMBuilder`), graph analysis packages (e.g. `networkX`) and graph visualization programs (e.g. `Gephi`).\n", "\n", "### What are conformation space networks?\n", "\n", "A conformation space network is a visualization of a free energy landscape, where each node is a cluster of molecular conformations, and the edges show which conformations can directly interconvert during a molecular dynamics simulation. A CSN can be thought of as a visual representation of a transition matrix, where the nodes represent the row / column indices and the edges show the off-diagonal elements. `CSNAnalysis` offers a concise set of tools for the creation, analysis and visualization of CSNs.\n", "\n", "**This tutorial will give quick examples for the following use cases:**\n", "\n", "1. Initializing CSN objects from count matrices\n", "2. Trimming CSNs\n", "2. Obtaining steady-state weights from a transition matrix\n", " * By eigenvalue\n", " * By iterative multiplication\n", "3. Computing committor probabilities to an arbitrary set of basins\n", "4. Exporting gexf files for visualization with the Gephi program" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting started\n", "\n", "Clone the CSNAnalysis repository:\n", "\n", "```\n", "git clone https://github.com/ADicksonLab/CSNAnalysis.git```\n", "\n", "Navigate to the examples directory and install using pip:\n", "\n", "```\n", "cd CSNAnalysis\n", "pip install --user -e\n", "```\n", "\n", "Go to the examples directory and open this notebook (`examples.ipynb`):\n", "\n", "```\n", "cd examples; jupyter notebook```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dependencies\n", "\n", "I highly recommend using Anaconda and working in a `python3` environment. CSNAnalysis uses the packages `numpy`, `scipy` and `networkx`. If these are installed then the following lines of code should run without error:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import networkx as nx\n", "import scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If `CSNAnalysis` was installed (i.e. added to your `sys.path`), then this should also work:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from csnanalysis.csn import CSN\n", "from csnanalysis.matrix import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook also uses `matplotlib`, to visualize output." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import matplotlib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great! Now let's load in the count matrix that we'll use for all the examples here:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "count_mat = scipy.sparse.load_npz('matrix.npz')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Background: Sparse matrices" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "It's worth knowing a little about sparse matrices before we start. If we have a huge $N$ by $N$ matrix, where $N > 1000$, but most of the elements are zero, it is more efficient to store the data as a sparse matrix." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "scipy.sparse.coo.coo_matrix" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(count_mat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`coo_matrix` refers to \"coordinate format\", where the matrix is essentially a set of lists of matrix \"coordinates\" (rows, columns) and data:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0 382.0\n", "0 651 2.0\n", "0 909 2.0\n", "0 920 2.0\n", "0 1363 1.0\n", "0 1445 2.0\n", "0 2021 5.0\n", "0 2022 7.0\n", "0 2085 4.0\n", "0 2131 1.0\n" ] } ], "source": [ "rows = count_mat.row\n", "cols = count_mat.col\n", "data = count_mat.data\n", "\n", "for r,c,d in zip(rows[0:10],cols[0:10],data[0:10]):\n", " print(r,c,d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although it can be treated like a normal matrix ($4000$ by $4000$ in this case):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(4000, 4000)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "count_mat.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It only needs to store non-zero elements, which are much fewer than $4000^2$:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "44163" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(rows)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**OK, let's get started building a Conformation Space Network!**\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1) Initializing CSN objects from count matrices\n", "\n", "To get started we need a count matrix, which can be a `numpy` array, or a `scipy.sparse` matrix, or a list of lists:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "our_csn = CSN(count_mat,symmetrize=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Any of the `CSNAnalysis` functions can be queried using \"?\"" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "CSN?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `our_csn` object now holds three different representations of our data. The original counts can now be found in `scipy.sparse` format:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "<4000x4000 sparse matrix of type ''\n", "\twith 62280 stored elements in COOrdinate format>" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.countmat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A transition matrix has been computed from this count matrix according to: \n", "\\begin{equation}\n", "t_{ij} = \\frac{c_{ij}}{\\sum_j c_{ij}}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "<4000x4000 sparse matrix of type ''\n", "\twith 62280 stored elements in COOrdinate format>" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.transmat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the elements in each column sum to one:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "matrix([[1., 1., 1., ..., 1., 1., 1.]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.transmat.sum(axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, the data has been stored in a `networkx` directed graph:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.graph" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "that holds the nodes and edges of our csn, and we can use in other `networkx` functions. For example, we can calculate the shortest path between nodes 0 and 10:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 1445, 2125, 2043, 247, 1780, 10]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nx.shortest_path(our_csn.graph,0,10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 2) Trimming CSNs\n", "\n", "A big benefit of coupling the count matrix, transition matrix and graph representations is that elements can be \"trimmed\" from all three simultaneously. The `trim` function will eliminate nodes that are not connected to the main component (by inflow, outflow, or both), and can also eliminate nodes that do not meet a minimum count requirement:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "our_csn.trim(by_inflow=True, by_outflow=True, min_count=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The trimmed graph, count matrix and transition matrix are stored as `our_csn.trim_graph`, `our_csn.trim_countmat` and `our_csn.trim_transmat`, respectively." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2282" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.trim_graph.number_of_nodes()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2282, 2282)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.trim_countmat.shape" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2282, 2282)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.trim_transmat.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3) Obtaining steady-state weights from the transition matrix\n", "\n", "Now that we've ensured that our transition matrix is fully-connected, we can compute its equilibrium weights. This is implemented in two ways.\n", "\n", "First, we can compute the eigenvector of the transition matrix with eigenvalue one:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "wt_eig = our_csn.calc_eig_weights()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This can exhibit some instability, especially for low-weight states, so we can also calculate weights by iterative multiplication of the transition matrix, which can take a little longer:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "wt_mult = our_csn.calc_mult_weights()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xUVfrH8c9DaFEpClgoCioWUJcSQeyrorAWUFlBWcv+QGzYRcGK2EBcXQEbinVtrAqyKkYUBBVBQhcQBUQpKiAdIiTh+f1xb3AYQzKRTGYm+b5fr3kxc+65N89hNA/nnnPPMXdHREQkVhUSHYCIiKQWJQ4RESkWJQ4RESkWJQ4RESkWJQ4RESmWiokOoDTUrl3bGzZsmOgwRERSytSpU1e5e53o8nKROBo2bEhWVlaiwxARSSlm9kNB5bpVJSIixaLEISIixaLEISIixaLEISIixaLEISIixVIuZlWJiJQnI6cvY2DmfJavzaZuzXR6nXEoHZvXK7HrK3GIiJQhI6cvo887s8nOyQNg2dps+rwzG6DEkkdcb1WZWTszm29mC8ysdwHHq5jZm+HxyWbWMCxvZWYzwtdMMzs31muKiJRnAzPnk52TxyErF9Nr/EvgTnZOHgMz55fYz4hb4jCzNOAJoD3QBLjQzJpEVesGrHH3g4HHgAFh+ddAhrs3A9oBz5hZxRivKSJSbq38dT3Xf/4a7714A11mZrLfhlUALF+bXWI/I563qloBC9x9EYCZvQF0AOZG1OkA9A3fvwUMMTNz980RdaoC+btNxXJNEZHyacoURr9yEwf98j0jm5xEv1N7sHq3GgDUrZleYj8mnreq6gFLIj4vDcsKrOPuucA6oBaAmbU2sznAbODK8Hgs1yQ8v4eZZZlZ1sqVK0ugOSIiSWrzZrjlFjjmGOp6Nld17ssNZ/fanjTSK6XR64xDS+zHJe10XHef7O5NgaOBPmZWtZjnD3X3DHfPqFPnD2t0iYiUDePGwZFHwr/+BZdfTvq333DGbd2pVzMdA+rVTOeh845MmVlVy4AGEZ/rh2UF1VlqZhWBGsCvkRXcfZ6ZbQSOiPGaIiJl37p1cOutMHQoHHRQkEBOPhmAjs1rlGiiiBbPHscUoLGZNTKzykAXYFRUnVHApeH7TsBYd/fwnIoAZnYAcBiwOMZrioiUbf/7HzRpAs89F9yimjVre9IoDXHrcbh7rpn1BDKBNOB5d59jZv2ALHcfBQwDXjGzBcBqgkQAcDzQ28xygG3A1e6+CqCga8arDSIiSWXlSrj+enj99eD21MiRcPTRpR6GuXvRtVJcRkaGaz8OEUlZ7kGyuO46WL8e7roLbrsNKleO6481s6nunhFdrifHRUSS2ZIlcNVV8P770Lo1DBsGTZsmNKSknVUlIlKubdsGzzwTJIlx4+Cxx+CLLxKeNEA9DhGR5PPdd3D55TB+PJx6ajBz6sADEx3VdupxiIgki9xceOQROOoomDEjmDU1ZkxSJQ1Qj0NEJDnMmgXdukFWFnToAE8+CXXrJjqqAqnHISKSSFu2wN13Q8uW8OOPMHw4jBiRtEkD1OMQEUmcSZOCXsbcuXDxxcEAeK1aiY6qSOpxiIiUtk2b4MYb4dhjYcMG+OADePnllEgaoB6HiEjp+uSTYMbU99/D1VfDQw9B9eqJjqpY1OMQESkNa9dC9+5w2mlQsWIw1faJJ1IuaYASh4hI/L37brAo4YsvBkuFzJwJJ56Y6Kj+NN2qEhGJl19+CdaXGj4c/vKXYFXbli0THdUuU+IQEdlFI6cvY2DmfJavzaZuzXR6nX4IHeeMgxtugI0b4YEHoFcvqFQp0aGWCCUOEZFdMHL6Mvq8M5vsnDwA/Mcf2OuC22BhFrRpEyxKePjhCY6yZClxiIjsgoGZ88nOycN8G12nj6b3+Bcxdx47uyc3jvg3pKUlOsQSp8QhIlJMkbemHGi0ehn9Rw+i9dI5TGjYnNvb9WRZjX24sQwmDVDiEBEplshbU2nb8rj8qxHc+Pmr/FaxMrf87QbeOuJUMKNezfREhxo3ShwiIsWQf2uqyS+LGDD6cY78ZSGjDzmWu9texco99gQgvVIavc44NMGRxo8Sh4hIMaxatY5bJr7BlZPeYs1u1bmyYx8+PPQ4AAyCWVVnHErH5vUSG2gcKXGIiMRq4kQyX76Bhit/5K0jTuW+U7qzLr0aAPVqpvNF71MSHGDpUOIQESnKxo1w++0wZAh771OX7hfez8f7N9t+uKzfmoqmJUdERArz0UdwxBEwZAj07Mlu387jrF6XUa9mOkbQ03jovCPL9K2paOpxiIgUZPVquPnmYH2pQw+Fzz6D44KxjI7Nq5WrRBFNPQ4RkWhvvx0sSvjKK8EtqhkzticNiXPiMLN2ZjbfzBaYWe8CjlcxszfD45PNrGFY3tbMpprZ7PDPUyLO+TS85ozwtXc82yAi5cjPP0OnTsGrbt1g/+8HHoCqVRMdWVKJW+IwszTgCaA90AS40MyaRFXrBqxx94OBx4ABYfkq4Gx3PxK4FHgl6ryu7t4sfK2IVxtEpJxwD25JNWkC770XbK40eTI0a1bkqeVRPHscrYAF7r7I3bcCbwAdoup0AF4K378FnGpm5u7T3X15WD4HSDezKnGMVUTKq8WLoV07+Oc/oWnTYK+M3r3LzEq28RDPxFEPWBLxeWlYVmAdd88F1gHRm+6eD0xz9y0RZS+Et6nuMjMr6IebWQ8zyzKzrJUrV+5KO0SkLNq2DQYPDmZMTZwY7MY3fnwwEC6FSurBcTNrSnD76oqI4q7hLawTwtfFBZ3r7kPdPcPdM+rUqRP/YEUkdcybByecEGyydMIJ8PXXwf7fFZL6V2LSiOff0jKgQcTn+mFZgXXMrCJQA/g1/FwfGAFc4u4L809w92XhnxuA1whuiYmIFC0nBx58MBi7+OYbePll+OADOOCAREeWUuKZOKYAjc2skZlVBroAo6LqjCIY/AboBIx1dzezmsD7QG93/yK/splVNLPa4ftKwFnA13Fsg4iUFdOmQatWcMcd0KEDzJ0LF18MBd/tlkLELXGEYxY9gUxgHjDc3eeYWT8zOyesNgyoZWYLgJuA/Cm7PYGDgbujpt1WATLNbBYwg6DH8my82iAiZUB2NvTpEySNn3+Gd94J9gDfZ59ER5ayzN0THUPcZWRkeFZWVqLDEJHS9vnn0K0bfPst/N//wSOPwJ57JjqqlGFmU909I7pcI0EiUvZs2AA9ewYD31u3wpgxwd7fSholQolDRMqW0aOD5zGefBJuuCGYMXXaaYmOqkzRIocikrIi9/4+rNJWnpn+Kvu//zYcfjh88QW0aZPoEMskJQ4RSUnb9/7emsvf5n/BvWOepuZvG/jm8us5bPAAqKLFJuJFiUNEUtLAzPnssWYl//7oSc74bhKz9j2YSzr3Y/2BTflCSSOulDhEJPW4c9yEd7lz7DAq5+Xw4Mn/ZNjRHcmrkIatzU50dGWeEoeIpJZFi6BHDx7+5BMmNziC29pdy+K9fl8Gr27N9AQGVz4ocYhIasjLCxYlvOMOSEtjxu0P8U87is25vz+LVt72/k4UTccVkeQ3dy4cfzzceCOcfDLMmUOzB3rz4Pl/Kdd7fyeKehwikry2boUBA+D++6FaNfjPf+Cii7avL9WxeT0ligRQ4hCR5JSVFSwXMmsWdOkCjz8Oe2un6GSgW1Uiklw2b4Zbb4XWrWHVKnj3XXj9dSWNJKIeh4gkj/HjoXt3WLAALr8cBg6EGjUSHZVEUY9DRBJv/Xq46qpg4HvbNvjkExg6VEkjSanHISKlKnJ9qbo103lkt6W0GXgHLF8ON90E990Hu+2W6DClEEocIlJqtq8vlZPHnpvXccv/HqHN3E9Zf9ChVJ84MRjXkKSnxCEipWZg5nyyt+Zy9rwJ9P34Gapt2cy/j7uQEe0uZbySRspQ4hCRUpO3ZAnPfvQUbRdMZsZ+jbmt/fXMr9MQ25ib6NCkGJQ4RCT+3OG55/h42I2k5eVy31+78ULGOWyrkAZofalUo8QhIvG1cGEwtXbcOLIzjqVrq258W22f7Ye1vlTq0XRcEYmPvDx49FE48kiYOhWGDqXOV59zdfcztL5UilOPQ0RK3tdfB8uFfPUVnH02PPUU1AuSg9aXSn3qcYhIydm6Ffr2hRYtgn0zXn89WDKknhJFWRLXxGFm7cxsvpktMLPeBRyvYmZvhscnm1nDsLytmU01s9nhn6dEnNMyLF9gZoPMwmUyRSSxvvoqSBj33gsXXADz5gWLE+p/0TInbonDzNKAJ4D2QBPgQjNrElWtG7DG3Q8GHgMGhOWrgLPd/UjgUuCViHOeAi4HGoevdvFqg4jEYPNmuPlmaNMG1q2D994Llj+vXTvRkUmcxLPH0QpY4O6L3H0r8AbQIapOB+Cl8P1bwKlmZu4+3d2Xh+VzgPSwd7IfUN3dJ7m7Ay8DHePYBhEpzLhxweD3o49Cjx4wZw6ceWaio5I4i2fiqAcsifi8NCwrsI675wLrgFpRdc4Hprn7lrD+0iKuKSLxtm5dkChOOQUqVIBPPw0GwKtXT3RkUgqSenDczJoS3L664k+c28PMsswsa+XKlSUfnEh59b//QZMmMGwY9OoFM2fCSSclOiopRfFMHMuABhGf64dlBdYxs4pADeDX8HN9YARwibsvjKhfv4hrAuDuQ909w90z6tSps4tNERFWrIALL4RzzoFatWDyZHj4Ya1kWw7FM3FMARqbWSMzqwx0AUZF1RlFMPgN0AkY6+5uZjWB94He7v5FfmV3/wlYb2bHhLOpLgHejWMbRMQdXn016GW8/Tb06xds65qRkejIJEGKTBxm9vdYyqKFYxY9gUxgHjDc3eeYWT8zOyesNgyoZWYLgJuA/Cm7PYGDgbvNbEb4yt838mrgOWABsBAYXVQsIvInLVkSPMD3j39A48YwYwbcdRdUrpzoyCSBLJicVEgFs2nu3qKosmSWkZHhWVlZiQ5DJHVs2xbswHfrrcHSIQ8+CD17QlpaoiOTUmRmU939D13LnS45Ymbtgb8B9cxsUMSh6oDWQBYpq777LliUcPx4OPXUIIEceGCio5IkUtitquVAFvAbMDXiNQo4I/6hiUipys2FgQPhqKOCW1LDhsGYMUoa8gc77XG4+0xgppm95u45pRiTiJS2mTODRQmnToWOHeGJJ6Bu3URHJUkqlllVrcxsjJl9a2aLzOx7M1sU98hEJP62bAkGuzMygoHw4cPhnXeUNKRQsSyrPgy4keA2VV58wxGRUvPll0EvY948uOSSYNmQWtELN4j8USyJY527a8qrSFmxaRPccQcMGgT168MHH0D79omOSlJIYbOq8qfbjjOzgcA7wJb84+4+Lc6xiUhJ+/jjYMbU4sVwzTXw0ENQrVqio5IUU1iP419RnyPn8jpwCiKSGtasgVtugeefDx7kmzABTjgh0VFJiipsVtVfSzMQEYmTESPg6qth5Uro3RvuvhvS0xMdlaSwIsc4zOymAorXAVPdfUbJhyQiJeKXX+Daa+G//4VmzeD994Md+kR2USzTcTOAKwn2vahHsMR5O+BZM7s1jrGJyJ/hDi+/DIcfHuz3/cADv2/rKlICYplVVR9o4e4bAczsHoKVa08kmKL7cPzCE5Fi+fFHuOIK+PBDOPbY4Onvww5LdFRSxsTS49ibiNlUQA6wj7tnR5WLSKJs2xY87d20KXz2WTDV9rPPlDQkLmLpcbwKTDaz/H0vzgZeM7Pdgblxi0xEYjN/PnTvDp9/Dm3bBosSNmyY6KikDCsycbj7fWY2GjguLLrS3fPXKO8at8hEZAcjpy9jYOZ8lq/Npm7NdG495UA6jH0D+vYNduF78cXgCXCzRIcqZVxhDwBWd/f1ZrYXsCh85R/by91Xl0aAIhIkjT7vzCY7J1j1p+b8r2n8eA/4eSGcfz4MGQL77pvgKKW8KKzH8RpwFsEAuAMW9afWWhYpJQMz55Odk0eV3K1cO/ENrpz0Fmt2q84dXe/hgf/0TXB0Ut4U9gDgWeGfjUovHBEpyPK12bRcOpeHRw/ioNVL+e8Rp3H/Kd1Yn16NBxIdnJQ7sTwAaARjGY3C8Y79gX3d/au4RycisHEjD08YxvlfjmR59TpcfEE/PmsUPJNRr6aeAJfSF8usqieBbQRrU90HbADeBo6OY1wiApCZCT160GnJEv5z9Nk8dPzFbK4cJIv0Smn0OuPQBAco5VEsz3G0dvdrCLaQxd3XAJXjGpVIebd6NVx2GbRrB7vthn32GdWeeZI9994LI+hpPHTekXRsXi/RkUo5FEuPI8fM0ggGxDGzOgQ9EBGJh7ffDpY8X7Uq2DfjzjuhalU6ghKFJIVYEscgYASwt5k9AHQC7oxrVCLl0U8/Qc+ewdatzZsHy4Y0a5boqET+IJYHAF81s6nAqQRTcTu6+7y4RyZSxm1/oG/NZrovmkCvzKFU3vob9O8PN98MFWP5d51I6YtlVtV9wATgRXffVJyLm1k74HEgDXjO3ftHHa8CvAy0BH4FOrv7YjOrBbxFMAD/orv3jDjnU2A/IDssOt3dVxQnLpFEy3+gr9aq5bz04RBOXDydrAZNWfv4k5x27omJDk+kULH8k2YRcCEwyMw2AJ8BE9z93cJOCsdFngDaAkuBKWY2yt0j17fqBqxx94PNrAswAOhMMBB/F3BE+IrWNWLZE5GU86/Rc+n85Qh6TXgZN+POtlfxavP21J2fy2mJDk6kCLHcqnoBeMHM9gUuAG4BegBFbVTcCljg7osAzOwNoAM7LozYAegbvn8LGGJmFvZsPjezg4vRFpHUMG8ejz15PRnL5vFpo5bc3u4allffGwge9BNJdkVOxzWz58xsIvAUQaLpBOwZw7XrAUsiPi8Nywqs4+65BDsL1orh2i+Y2Qwzuyt8QLGguHuYWZaZZa1cuTKGS4rEWU5OsKlSs2Y0Xr2UG8+8icv+3nd70gCoqwf6JAXE8hxHLYIxirXAamBV+Es+Ubq6+5HACeHr4oIquftQd89w94w6deqUaoAifzBtGhx9dDC1tmNHJo6awIfN2+6wkq0e6JNUUWTicPdz3b01wU5/NYFxZrY0hmsvAxpEfK4flhVYx8wqAjUIBskLi2dZ+OcGgoUYW8UQi0hiZGdD797QqlWwB/iIEfDmm7Q/rRkPnXck9Wqm64E+STmxzKo6i+Bf9icSJI6xBAPkRZkCNDazRgQJogtwUVSdUcClwJcEt8DGursXEktFoKa7rzKzSgSr934cQywipW/ChGCDpe++g27dYOBA2PP3u7wdm9dTopCUFMusqnYEieJxd18e64XdPdfMegKZBLe6nnf3OWbWD8hy91HAMOAVM1tAcBusS/75ZrYYqA5UNrOOwOnAD0BmmDTSCJLGs7HGJFIq1q+HPn3gySehUSMYMwZO01wpKTuskH/glxkZGRmelaXZu1IKRo+GK66ApUvh+uvh/vth990THZXIn2JmU909I7pcj6aKlIRff4Ubb4RXXoEmTWDiRDjmmERHJRIXscyqEpGdcYfhw+Hww+H11+Guu4IZVEoaUobF8hzH9bGUiZQ7y5fDuedC586w//4wdSr06wdVqiQ6MpG4iqXHcWkBZZeVcBwiqcMdhg0LbkllZsLDD8OkSXDUUYmOTKRU7HSMw8wuJJg+28jMRkUcqkYwA0qk/Fm0CC6/HMaOhZNOgmefhcaNEx2VSKkqbHB8IvATUBv4V0T5BmBWPIMSSTp5eTB4cLCxUloaPP10kEAqaJhQyp+dJg53/4HguYk2pReOSBKaMyd4gG/yZDjzzCBp1K+f6KhEEqawW1UbCLeLjT4EuLtXj1tUIslg69ZgU6X774fq1eHVV+HCC3dYX0qkPCqsx1HUsukiZdeUKUEvY/Zs6NIFBg0CLZYpAsS2VtX+BZW7+48lH45I6dq+fevabOrWTKf3iftz9ohn4NFHYd994d134ZxzEh2mSFKJ5cnx9yPeVwUaAfOBpnGJSKSU5G/fmp2TB0CDWZM56uGLYc1y6NEjmGZbo0aCoxRJPrHsAHhk5GczawFcHbeIRErByOnLuHn4TPLcqbZlE70/fYGuMz5kcc39uLb7Iwx+5uZEhyiStIq9VpW7TzOz1vEIRqQ05Pc08tz568IpPPjhEPbetIahR5/Loyd0ZUulqgxOdJAiSSyWMY6bIj5WAFoAMS+vLpJsBmbOJ33dah76ZCgd547nm9oHcOW5tzOzbrD7Xj1t3ypSqFh6HJGzq3IJxjzejk84InHmTsuJH3LPx89QbctmHjvuIp5s83dy0ioB2r5VJBaxjHHcWxqBiMTd0qVw1VUMeu89Zux3CLe2v45v6zTcfjjNTNu3isSgsAcAR+3sGIC7a46iJK3Iabb1qldh8OapNB/8IOTkMPume+ha9Wg25f1eP71SmpKGSIwK63G0AZYArwOTCZ4YF0l6kdNsD1iznP6vD6b5j7NZmXEsdd54mSMPOogHop7f6HXGoUoaIjEqLHHsC7QF8lfJfR943d3nlEZgIn/WwMz5bN2yle5Z73LzZ6+SUyGN29pdy+cnduCLgw4CoGPzekoUIn9SYUuO5AEfAh+aWRWCBPKpmd3r7kNKK0CR4qr23TyeGP04zX76jjEHt+LO06/ml2q1sXW/JTo0kTKh0MHxMGGcSZA0GgKDgBHxD0vkT9iyBR58kPdefIC1Vfeg5zm38t5hJ2xflLCuptmKlIjCBsdfBo4APgDudfevSy0qkeKaPDlYlHDOHH5qfy4XHN6Znyrtsf2wptmKlJzCdqH5B9AYuB6YaGbrw9cGM1tfOuGJFGHTJrjpJmjTBtatg/feo8EH73DbP46nXs10jOCBPs2YEik5hY1xaGszSW5jxwa78C1aBFddFeydUT3YJkaD3yLxE9fkYGbtzGy+mS0ws94FHK9iZm+GxyebWcOwvJaZjTOzjWY2JOqclmY2OzxnkJl21Sl31q4NEsappwZbt376KTz55PakISLxFbfEYWZpwBNAe6AJcKGZNYmq1g1Y4+4HA48BA8Ly34C7gFsKuPRTwOUEt9EaA+1KPnpJWu++C02awPPPw623wqxZcNJJiY5KpFyJZ4+jFbDA3Re5+1bgDaBDVJ0OwEvh+7eAU83M3H2Tu39OkEC2M7P9gOruPsndHXgZ6BjHNkiyWLEi2ImvY0eoXTsYDB8wANI1U0qktMUzcdQjePI839KwrMA67p4LrANqFXHNpUVcEwAz62FmWWaWtXLlymKGLknDHf7zHzj8cBgxAu67D7KyICMj0ZGJlFtldgDc3Ye6e4a7Z9TRXtGpackSOOssuPhiOOQQmD4d7rwTKldOdGQi5VqxN3IqhmVAg4jP9cOyguosNbOKQA3g1yKuWb+Ia0oKil6UcMjGKTQb8hDk5cG//w09e0JaWqLDFBHimzimAI3NrBHBL/cuBGteRRoFXAp8CXQCxoZjFwVy95/CZ0mOIVh48RLQZm2pLnJRwkarl9H/tcE0W/I1K1qfwN6vvwSNGiU6RBGJELfE4e65ZtYTyATSgOfdfY6Z9QOy3H0UMAx4xcwWAKsJkgsAZrYYqA5UNrOOwOnuPpdgv/MXgXRgdPiSFJa/KOEVU0Zw4+evsTWtEr3aX8fEE87hCyUNkaQTzx4H7v4BwZIlkWV3R7z/Dfj7Ts5tuJPyLIKlUKSMqPHtHJ7+4HGO/GUhmY2P4a62V7GiWi0tSiiSpOKaOESiRY5lHLBHGs/88CGjXhrCmqrVuKpDb0YfepwWJRRJckocUmoixzJaLJvHgNGDaPzrEqaffDZXt/gHP1XafXtdLUookryUOKTUDMycj23ayN0TXuGyqf9jefXaXPr3e1nQ4nhuO+NQ7cgnkiKUOKTUHDh9Ig9mDqHBul94qcWZPHzipWyqshu2NluLEoqkECUOib81a+Dmm3ll+Ass3Ksef7+oP1Ma/D6/QWMZIqlFiUNKVOTgd92a6TxaaRGtH7kLVq5k/j97csE+bVnnvz/Ip7EMkdSjxCElJnLwu87GNdw+8iFaz/+CtYc2peb773NoixbcG5VYNJYhknqUOGSX5fcylq3NBnfOmzOWuz95lvScLTx84iW8d3pXJrRoAWiDJZGyQIlDdklkL6PeuhU8mDmEk76fRla9w7mt/XUsrNUA25CT6DBFpAQpccguGZg5n9+25nDJtPe5bXywtcrdp13BKy3OxC1YfFmD3yJlixKH7JKqC79j+OhBHL1sLhMaNuf2dj1ZWmOf7cc1+C1S9ihxyJ+TkwOPPMIHL95DdsUq3Py3G3n7iFO2LxcCUE+D3yJlkhKHFN/06dCtG0yfzqrTzqTLERexpEqN7YfTK6Xx0HlHKmGIlFFKHFKoO0fO5vXJS8hzJz0vh6e/f4+TRr4Y7Pv99tvUO+88btYUW5FyRYlDdurOkbP5z6QfAchYOocBowdx0OplTD2lIy3feh723BPQFFuR8kaJQ3YQ+eS3A7tv2cytE17i0mnvs7T63lx8QT8mHtiShWHSEJHyR4lDtrtz5GxenfQj+Xv3nrhoKg9mDqHu+lW80PJsBp54CZsrp8POd/cVkXJAiUOAoKeRnzRqZG/g7rHPcv7XY1mwV306dX2YafUP3143LWLmlIiUP0ocAgQP8jnQ/pvP6TfmaWr+toHBbToz5NjObKlYeYe6F7ZukJggRSQpKHEIADlLl/H0mKdo9+2XzN7nIC69oB9z9zlwhzppZlzYugH3dzwyQVGKSDJQ4ijv3OHFF/l42HVUydlC/5Mu49lW55JXIVj63IDHOjfTrCkR2U6Jo5yJnDXVcttanpzwDHtP/oytzVvT+ZjuzKu+3/a6BnQ9Zn8lDRHZgRJHOTFy+jL6jprD2uwcKmzL49Jp73PrhJdwq8CMPg/S7P7buGLmT3qQT0SKFNfEYWbtgMeBNOA5d+8fdbwK8DLQEvgV6Ozui8NjfYBuQB5wnbtnhuWLgQ1hea67Z8SzDWVB5NLnB61awsOjH6fl8m8Yd2BL7jjjGqz6AXxRoYIe5BORmMQtcZhZGvAE0BZYCkwxs1HuPjeiWjdgjbsfbGZdgAFAZzNrAnQBmgJ1gY/N7BB3zwvP+6u7r4pX7HvuPSoAABBDSURBVGXNwMz55Py2hZ6T3+LaiW+wuVI6N5x1MyObnAxm2NrsRIcoIikknj2OVsACd18EYGZvAB2AyMTRAegbvn8LGGJmFpa/4e5bgO/NbEF4vS/jGG+Ztdc3s3nug39z+MrFvHfYCdxz2hX8unvN7ce1X4aIFEc8E0c9YEnE56VA653VcfdcM1sH1ArLJ0Wdm38PxYGPzMyBZ9x9aEE/3Mx6AD0A9t9//11rSarKzoa+fRn58iOs2r0mPc69g48OabNDFe2XISLFlYqD48e7+zIz2xsYY2bfuPuE6EphQhkKkJGRUf7WyJgwAbp3h+++Y0nHLlxw8PmsSNuxZ7HnbpW45+ymGtcQkWKpEMdrLwMiHzGuH5YVWMfMKgI1CAbJd3quu+f/uQIYQXALS/KtXw9XXw0nnQS5ufDxxzQc8Tq3X9SGejXTMYINlv7duRnT7z5dSUNEii2ePY4pQGMza0TwS78LcFFUnVHApQRjF52Ase7uZjYKeM3MHiUYHG8MfGVmuwMV3H1D+P50oF8c25BaPvgArrwSli6FG2+E++6D3XcHtPS5iJScuCWOcMyiJ5BJMB33eXefY2b9gCx3HwUMA14JB79XEyQXwnrDCQbSc4Fr3D3PzPYBRgTj51QEXnP3D+PVhmQW+SDf4ZW28sy0/9Dgg3egSROYOBGOOSbRIYpIGWVeDpbIzsjI8KysrESHUWK2P5exNZezvvmMvh8/Q43fNrKw27UcNrg/VKmS6BBFpAwws6kFPSuXioPj5VJkD6OCGbXWr2LQR0/SdsFkZu7bmH90vp8NBzbhCyUNEYkzJY4U0PXZL/li4erggzudZmZyx7jnqZyXwwMn/x/PH92BvAppepBPREqFEkcSi9zzG6DB2p/p/+EgjvthFpMaHMFt7a/jhz3rbj+uB/lEpDQocSSpto9+yncrNgFQYVse/5z6P26Z8Aq5FSrQ54yevPGX03H7fTa1HuQTkdKixJFkgoHvWWTnbAPgkJWLGTB6MM1/ms8nBx3NHadfw8/VawPBxkrb3LWSrYiUKiWOJHLnyNnb9/2ulJfD1V/+l2u+HM6GKrtx3dm9GHX4iRCx3/e/LviLkoWIlDoljiQxcvqy7UnjqJ++5eEPHuewVT/w7uEnce9pPVi9W40d6h930F5KGiKSEEocCTRy+jLuGDGbTVuD1eKr5vzGTZ+9Sresd1mx+550O/8uPjl4x3Uh83fl077fIpIoShwJssMUW+CYH2fRf/RgGq79iVebtaP/yf9kQ5Xddzin8d67M+amk0s5UhGRHSlxlLLoXka1LZvoM+4FLpr5IYtr7seFXR7kywOO2uGcCgYXtVYvQ0SSgxJHKYruZZy6YDIPZD5BnU1reabVeTx2/EX8Vqnq9uO6LSUiyUiJoxSMnL6MG96csf3zXpvXcc/HQ+kwbzzf1D6AK869g5l1d3wGo56m2IpIklLiiLPIB/lw55x54+n78VD22LKZR4/vylPHdCInrdIO5/xDvQwRSWJKHHESfVtq3/WruP+jJzht4RSm73cot7a/ju/qHPCH8447aC8lDRFJakocJSz6tpT5Ni6cmUmfcc9Tcds27julOy+0PJttFdJ2OC+9UgUeOu8o3ZoSkaSnxFGCDu7zPrkR25scsGY5/T8cTJsfZ/PFAUfRu911LKm57w7nHHfQXrx6eZtSjlRE5M9T4igB0avYpm3L4/+mvMvNn/+HrRUqclu7a3nzqNN3WC4EoHqVNCUNEUk5Shy7KLqXcdiK7xkwehB/+fk7xhzcmjtPv4pfqtX+w3kVDWbd264UIxURKRlKHH9SdMKonJvDNV8O5+pJw1lXdQ96nnMr7x12wh96GaAnwEUktSlx/AkNe7+/w+fmy75hwOhBHPLrj7zT9K/cd0p31kQtSgjBA33f9z+zlKIUEYkPJY5iiE4Y6Vt/4+bPXuH/skbxc7VaXNbpHj496OgCz92nWmUm39G2NMIUEYkrJY4YRSeNYxfPoP+Hg9l/3S+80vxvDDjpMjZW2a3AcxerlyEiZYgSRxGiE0b13zZy+7jn6TLrIxbtWZcLLurPVw2OKPBcjWWISFmkxFGI6KTR9rtJ3P/Rk9TetJanW5/PY8ddxJZKVQo8V70MESmrKsTz4mbWzszmm9kCM+tdwPEqZvZmeHyymTWMONYnLJ9vZmfEes2SEpk0am9aw5B3B/DsO/ezOr06HS/+F/1P/meBSWOfapWVNESkTItbj8PM0oAngLbAUmCKmY1y97kR1boBa9z9YDPrAgwAOptZE6AL0BSoC3xsZoeE5xR1zV22PWm403Hup9zz8VB2y8lm4AkX80zr88lNK/ivTQlDRMqDeN6qagUscPdFAGb2BtABiPwl3wHoG75/CxhiZhaWv+HuW4DvzWxBeD1iuGaJqJiXy9B37ueURVlMrXsYt7a/noW1G+y0vpKGiJQX8Uwc9YAlEZ+XAq13Vsfdc81sHVArLJ8UdW7+6n9FXRMAM+sB9ADYf//9ix18blpFFu1VjwmNWvByizP/sChhPiUMESlvyuzguLsPBYYCZGRkeBHVC3T/qZcXelxJQ0TKo3gmjmVA5L2d+mFZQXWWmllFoAbwaxHnFnXNuFPCEJHyLJ6zqqYAjc2skZlVJhjsHhVVZxRwafi+EzDW3T0s7xLOumoENAa+ivGau2xnieHfnZspaYhIuRe3Hkc4ZtETyATSgOfdfY6Z9QOy3H0UMAx4JRz8Xk2QCAjrDScY9M4FrnH3PICCrhmP+JUgREQKZsE/8Mu2jIwMz8rKSnQYIiIpxcymuntGdHlcHwAUEZGyR4lDRESKRYlDRESKRYlDRESKpVwMjpvZSuCHP3l6bWBVCYaTKGpH8ikrbVE7kktJtuMAd68TXVguEseuMLOsgmYVpBq1I/mUlbaoHcmlNNqhW1UiIlIsShwiIlIsShxFG5roAEqI2pF8ykpb1I7kEvd2aIxDRESKRT0OEREpFiUOEREplnKVOMysnZnNN7MFZta7gONVzOzN8PhkM2sYcaxPWD7fzM6I9ZrxEqe2LDaz2WY2w8xKZVXIP9sOM6tlZuPMbKOZDYk6p2XYjgVmNijcjjgV2/FpeM0Z4WvvJG5HWzObGv69TzWzUyLOSaXvo7B2lPr3sYttaRUR60wzOzfWaxbJ3cvFi2AZ9oXAgUBlYCbQJKrO1cDT4fsuwJvh+yZh/SpAo/A6abFcM1XaEh5bDNROke9kd+B44EpgSNQ5XwHHAAaMBtqnaDs+BTJS5PtoDtQN3x8BLEvR76OwdpTq91ECbdkNqBi+3w9YQbCVxi7/3ipPPY5WwAJ3X+TuW4E3gA5RdToAL4Xv3wJODf911AF4w923uPv3wILwerFcM1Xakgh/uh3uvsndPwd+i6xsZvsB1d19kgf/x7wMdIxrK+LQjgTZlXZMd/flYfkcID38l3CqfR8FtiPO8RZmV9qy2d1zw/KqQP5MqF3+vVWeEkc9YEnE56VhWYF1wr/wdUCtQs6N5ZrxEI+2QPAf1kdhF71HHOKOtivtKOyaS4u4ZkmLRzvyvRDearirFG7xlFQ7zgemufsWUvv7iGxHvtL8PnaIM1SstphZazObA8wGrgyP7/LvrfKUOKRox7t7C6A9cI2ZnZjogMq5ru5+JHBC+Lo4wfEUycyaAgOAKxIdy67YSTtS7vtw98nu3hQ4GuhjZlVL4rrlKXEsAxpEfK4flhVYx8wqAjWAXws5N5ZrxkM82oK75/+5AhhB/G9h7Uo7Crtm/SKuWdLi0Y7I72MD8BpJ/n2YWX2C/24ucfeFEfVT6vvYSTsS8X3sEGfoT/235e7zgI2E4zYxXLNwpTnQk8gXwaDQIoIB4fwBoaZRda5hx0Gm4eH7puw4oLyIYICpyGumUFt2B6qFdXYHJgLtkrUdEccvo+jB8b+lWjvCa9YO31ciuHd9ZbK2A6gZ1j+vgOumzPexs3Yk4vsogbY04vfB8QOA5QQr5+7y7624NjrZXsDfgG8JZhTcEZb1A84J31cF/kswYPwVcGDEuXeE580nYlZIQddMxbYQzLCYGb7mlFZbdrEdi4HVBP+SWko4MwTIAL4OrzmEcIWEVGoHQfKeCswKv4/HCWe/JWM7gDuBTcCMiNfeqfZ97Kwdifo+drEtF4exzgCmAR0Lu2ZxXlpyREREiqU8jXGIiEgJUOIQEZFiUeIQEZFiUeIQEZFiUeIQEZFiUeKQMsPM8iJWA52Rv+qnmT1nZk2SIL6NCf75GWY2qIg6Dc3s650cu8zM6sYnOkklFRMdgEgJynb3ZtGF7t49EcEkG3fPAnZlufzLCJ7HWF5EPSnj1OOQMi/cRyEjfN/NzL41s6/M7Nn8PTDMrI6ZvW1mU8LXcWF5XzN7PrzGIjO7Lizvb2bXRPyMvmZ2i5ntYWafmNm0cE+HP6w6amYnm9l7EZ+HmNll4fuWZjY+XGgyM1xdNvLcNDP73gI1w17WieGxCWbW2Mx2D2P+ysym58cQ+XPD9o4xszlhj+wHM6sd/pi08O9mjpl9ZGbpZtaJ4EG+V8PeXHrJfDuSipQ4pCxJj7pV1TnyYHib5S6C5S+OAw6LOPw48Ji7H02wKupzEccOA84gWJvoHjOrBLwJXBBR54Kw7DfgXA8Wi/wr8K9YV1ENrzsY6OTuLYHngQci67h7HsET/00I9vGYBpwQLv3dwN2/I1gZYKy7twpjGGhmu0f9uHvCOk0Jls/YP+JYY+CJ8Nha4Hx3f4ugt9LV3Zu5e3YsbZKySbeqpCwp8FZVhFbAeHdfDWBm/wUOCY+dBjSJ+B1f3cz2CN+/78HS2lvMbAWwj7tPN7O9w2RUB1jj7kvCX/4Phr2AbQTLVe8D/BxD/IcSLEI3JowjDfipgHqfAScSrDX0EHA5MB6YEh4/HTjHzG4JP1dlx8QAQdI5F8DdPzSzNRHHvnf3GeH7qUDDGGKXckSJQyRQATjG3aM3hgKI3I8hj9//v/kv0AnYl6C3AdCVIJG0dPccM1tM8Is7Ui479vbzjxswx93bFBHrBOAqoC5wN9ALOJkgoeRf53x3nx/Vln2KuG6+6PbqtpTsQLeqpDyZApxkZnuGy0+fH3HsI+Da/A9mVljPJd+bBKuRdiJIIhAsab0iTBp/JViVNNoPBL2bKmZWEzg1LJ8P1DGzNmEMlcJ9IaJ9BRwLbAsT3QyCfSMmhMczgWvzb5GZWfMCrvEF4a02Mzsd2DOG9m4AqsVQT8o4JQ4pS6LHOPpHHvRgP4UHCX7xfkGwKu268PB1QIaZzTKzuQR7gBfK3ecQ/CJd5u75t5ReDa8zG7gE+KaA85YAwwlmKA0HpoflWwmS0AAzm0mQEI4t4PwtBDu4TQqLPgvjmB1+vo9g6e9ZFuz+dl8B4d8LnB5Ovf07wa20DUU0+UXgaQ2Oi1bHlXLFzPZw941hj2ME8Ly7j0h0XKUtHEzPc/fcsIfzVBHjQyLbaYxDypu+ZnYawbjCR8DIBMeTKPsDw82sArCVYIBdJCbqcYiISLFojENERIpFiUNERIpFiUNERIpFiUNERIpFiUNERIrl/wF5Y9VTEeC7QgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "plt.scatter(wt_eig,wt_mult)\n", "plt.plot([0,wt_mult.max()],[0,wt_mult.max()],'r-')\n", "plt.xlabel(\"Eigenvalue weight\")\n", "plt.ylabel(\"Mult weight\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These weights are automatically added as attributes to the nodes in `our_csn.graph`:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'label': 0,\n", " 'count': 482,\n", " 'trim': 0.0,\n", " 'eig_weights': 0.002595528367725156,\n", " 'mult_weights': 0.0025955283677248217}" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.graph.node[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4) Committor probabilities to an arbitrary set of basins\n", "\n", "We are often doing simulations in the presence of one or more high probability \"basins\" of attraction. When there more than one basin, it can be useful to find the probability that a simulation started in a given state will visit (or \"commit to\") a given basin before the others.\n", "\n", "`CSNAnalysis` calculates committor probabilities by creating a sink matrix ($S$), where each column in the transition matrix that corresponds to a sink state is replaced by an identity vector. This turns each state into a \"black hole\" where probability can get in, but not out. \n", "\n", "By iteratively multiplying this matrix by itself, we can approximate $S^\\infty$. The elements of this matrix reveal the probability of transitioning to any of the sink states, upon starting in any non-sink state, $i$.\n", "\n", "Let's see this in action. We'll start by reading in a set of three basins: $A$, $B$ and $U$." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "Astates = [2031,596,1923,3223,2715]\n", "Bstates = [1550,3168,476,1616,2590]\n", "Ustates = list(np.loadtxt('state_U.dat',dtype=int))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then use the `calc_committors` function to calculate committors between this set of three basins. This will calculate $p_A$, $p_B$, and $p_U$ for each state, which sum to one." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "basins = [Astates,Bstates,Ustates]\n", "labels = ['pA','pB','pU']\n", "comms = our_csn.calc_committors(basins,labels=labels)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The committors can be interpreted as follows:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "comms[0] = [0.26406217 0.29477873 0.44115911]\n", "\n", "In other words, if you start in state 0:\n", "You will reach basin A first with probability 0.26, basin B with probability 0.29 and basin U with probability 0.44\n" ] } ], "source": [ "i = our_csn.trim_indices[0]\n", "print('comms['+str(i)+'] = ',comms[i])\n", "print('\\nIn other words, if you start in state {0:d}:'.format(i))\n", "print('You will reach basin A first with probability {0:.2f}, basin B with probability {1:.2f} and basin U with probability {2:.2f}'.format(comms[i,0],comms[i,1],comms[i,2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5) Exporting graph for visualization in Gephi\n", "\n", "`NetworkX` is great for doing graph-based analyses, but not stellar at greating graph layouts for large(r) networks. However, they do have excellent built-in support for exporting graph objects in a variety of formats. \n", "\n", "Here we'll use the `.gexf` format to save our network, as well as all of the attributes we've calculated, to a file that can be read into [Gephi](https://gephi.org/), a powerful graph visualization program. While support for Gephi has been spotty in the recent past, it is still one of the best available options for graph visualization.\n", "\n", "Before exporting to `.gexf`, let's use the committors we've calculated to add colors to the nodes:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "rgb = our_csn.colors_from_committors(comms)\n", "our_csn.set_colors(rgb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have added some properties to our nodes under 'viz', which will be interpreted by Gephi:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'label': 0,\n", " 'count': 482,\n", " 'trim': 0.0,\n", " 'eig_weights': 0.002595528367725156,\n", " 'mult_weights': 0.0025955283677248217,\n", " 'pA': 0.26406216543613925,\n", " 'pB': 0.2947787254045238,\n", " 'pU': 0.4411591091593356,\n", " 'viz': {'color': {'r': 152, 'g': 170, 'b': 255, 'a': 0}}}" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "our_csn.graph.node[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can use an internal `networkx` function to write all of this to a `.gexf` file:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "nx.readwrite.gexf.write_gexf(our_csn.graph.to_undirected(),'test.gexf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After opening this file in Gephi, I recommend creating a layout using the \"Force Atlas 2\" algorithm in the layout panel. I set the node sizes to the \"eig_weights\" variable, and after exporting to pdf and adding some labels, I get the following:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![Gephi graph export](committor_net_3state.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**That's the end of our tutorial!** I hope you enjoyed it and you find `CSNAnalysis` useful in your research. If you are having difficulties with the installation or running of the software, feel free to create an [issue on the Github page](https://github.com/ADicksonLab/CSNAnalysis)." ] }, { "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.7" } }, "nbformat": 4, "nbformat_minor": 1 }