tram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "504c2d7f",
   "metadata": {},
   "source": [
    "# TRAM\n",
    "This example shows how the transition-based reweighting analysis method (TRAM) <cite data-footcite=\"wu2016multiensemble\">(Wu et al, 2016)</cite> can be used to estimate a Markov state model of a simple double-well potential.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In the TRAM problem, we have a system which is governed by some unknown equilibrium distribution, $\\mu(x) = e^{f - u(x)}$, where $u(x)$ is the dimensionless potential function. $f$ is the free energy, which ensures normalization of $\\mu(x)$, and is unknown.\n",
    "\n",
    "We assume that we have sampled our system at $K$ different thermodynamic states, which are related to our reference state through the addition of some bias potential, $b^k(x)$ (at state $k$, $k \\in \\{1, ... K\\}$). This means each thermodynamic state has an equilibrium distribution $\\mu^k(x) = e^{f^k(x) - b^k(x)} \\mu(x)$. Again, the free energy $f^k$ for state $k$ ensures normalization.\n",
    "\n",
    "The shape of the bias potentials depends on the chosen sampling method. In replica exchange for example, it depends on the temperature difference between states. For umbrella sampling, the bias potential is usually a quadratic well around some center coordinate. In the following example we will use umbrella sampling.\n",
    "\n",
    "For each thermodynamic state, samples are _not_ assumed to be taken from the global equilibrium distribution of that state, i.e. from $\\mu^k(x)$. In stead, each thermodynamic state $k$ is governed by a Markov State Model (MSM) consisting of $M$ states. Samples are assumed to be taken from the _local equilibrium_ that governs each markov state $S_i$, $i \\in \\{1, ... M\\}$, given by $\\mu^k_i = e^{f^k_i - b^k(x)} \\mu(x)$.\n",
    "\n",
    "TRAM is thus a combination of a thermodynamic reweighting estimator, that combines data from multiple thermodynamic states, and a MSM, allowing use of time-lagged data. It produces a maximum-likelihood estimate for the underlying Multi Ensemble Markov Model of the biased thermodynamic states.\n",
    "\n",
    "## Application to a double well\n",
    "To apply TRAM to a double well we use the deeptime MarkovStateModel for sampling, and KMeans for clustering the data into the $M$ Markov states. We will then use TRAM to estimate the TRAM model, which holds the transition matrices and free energies from the data. Finally, we use the model to re-obtain the PMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b22f2497",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from deeptime.data import tmatrix_metropolis1d\n",
    "from deeptime.markov.msm import MarkovStateModel, TRAM, TRAMDataset\n",
    "from deeptime.clustering import KMeans"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74ac4c8f",
   "metadata": {},
   "source": [
    "Define the potential we will be sampling from. This is a simple double well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2fabade2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApnUlEQVR4nO3deXycdbX48c+ZyUy2Zm+Sbkk3uu8ldEHBqghYloKAgiKKFxEu+nPBBVHRq8JFXK6y1npFFrnsCAhFBGUpSLFpaUpouqT7kq3ZJntm+f7+yExaStpM25l5nmdy3q9XXp1kns6cJ5OcnDnPdxFjDEoppZzPZXUASimlYkMTulJKJQlN6EoplSQ0oSulVJLQhK6UUkkixaonHj58uBk3bpxVT6+UUo60du3aA8aYwoHusyyhjxs3jvLycqueXimlHElEdh3pvkFbLiJSIiKviEiViLwnIl8f4JglItIqIuvDHzedaNBKKaWOTTQVegC43hizTkSygLUi8pIxZuNhx60yxpwb+xCVUkpFY9AK3RhTY4xZF77dBlQBo+MdmFJKqWNzTKNcRGQcMA94e4C7F4tIhYi8ICIzjvD/rxaRchEpb2hoOPZolVJKHVHUCV1EhgFPAt8wxvgOu3sdMNYYMwe4A3h6oMcwxqwwxpQZY8oKCwe8SKuUUuo4RZXQRcRDXzJ/yBjz1OH3G2N8xpj28O2VgEdEhsc0UqWUUkcVzSgXAf4IVBljfnOEY0aEj0NEFoQftzGWgSqllDq6aCr0DwGfBz52yLDEpSJyjYhcEz7mYqBSRCqA24FLTZzW5d1c28avXtxMU0dvPB5eKaXi6rcvb2HV1vhcQxx02KIx5g1ABjnmTuDOWAV1NDsOtHPnK9UsnTWS/ExvIp5SKaViwhjDHf+s5tqPTOS0SbG/jui4tVyy0jwAtHX7LY5EKaWOTWdvkGDIkJUWn0n6jkvo2eGE7usOWByJUkodG1+4EM1O98Tl8R2X0CN/2bRCV0o5TVu4ENUKPSzyl83XpQldKeUskbwV6TTEmuMS+sEKXVsuSiln0Qr9MB63i3SPu78XpZRSThHJW1laoR+UlZaiFbpSynEigzmy07VC76cJXSnlRJHBHNpDP0R2ukdbLkopx/F1BfC6XaSmxCf1OjKhZ6V5dBy6Uspx2rr9ZKWlEF76KuYcmdCz01Jo02GLSimHaesOxG2ECzg0oWuFrpRyIl+3P26zRMGhCT07LUVniiqlHEcr9AFkp3voCYToCQStDkUppaLm6/LHbYQLODSh62xRpZQTaYU+gP4VF/XCqFLKQfpGuWiF/j5aoSulnCYQDNHRG9SWy+H6V1zUC6NKKYdo74nvwlzg0ISuFbpSyml8XZF1XLRCfx/dhk4p5TQHV1rUCv19ssPfkMhfPKWUsrt4r4UODk3omd4URLRCV0o5hy/OKy2CQxO6yyVkpabo9H+llGNEKnRN6APoW89FK3SllDP07ycap80twNEJPUV76Eopx4hU6MNSNaF/QHa6R3voSinHaOv2k+F1k+KOX9p1bkLXbeiUUg7i647vwlzg6ISuPXSllHPEe2EucHBC142ilVJOEu/NLcDBCT3SQzfGWB2KUkoNSiv0o8hKSyFkoKNXN7lQStlfX0LXCn1AWbomulLKQfp2K7K4QheREhF5RUSqROQ9Efn6AMeIiNwuItUiskFE5scn3IOy+xfo0j66UsrejDEJqdCj+XMRAK43xqwTkSxgrYi8ZIzZeMgxnwQmhT8WAveE/42bg0voaoWulLK3nkCI3mAorrNEIYoK3RhTY4xZF77dBlQBow87bBnwgOmzGsgVkZExj/YQusmFUsopDi6da6MeuoiMA+YBbx9212hgzyGf7+WDSR8RuVpEykWkvKGh4RhDfT/d5EIp5RQHF+ayySgXERkGPAl8wxjjO/zuAf7LB8YTGmNWGGPKjDFlhYWFxxbpYbL610TXCl0pZW/9C3PZoUIXEQ99yfwhY8xTAxyyFyg55PMxwP4TD+/IIt8YXUJXKWV3idjcAqIb5SLAH4EqY8xvjnDYs8AV4dEui4BWY0xNDOP8gDSPG6/bpS0XpZTt9W9uEeeZotH8ufgQ8HngXRFZH/7ajUApgDFmObASWApUA53AlTGPdADZ6Sl6UVQpZXuJqtAHfXRjzBsM3CM/9BgDXBeroKKVlebRCl0pZXuRHrqtRrnYTXZail4UVUrZXlt3AJdAptcd1+dxdELvq9A1oSul7K2t209Wmoe+S5Lx4/CErhtFK6Xsz9cdiPssUXB4Qs/WCl0p5QBt3X6yUuPbPweHJ3Td5EIp5QS+rvivhQ4OT+jZ6R46e4P4gyGrQ1FKqSNKxG5F4PCEHvmL165VulLKxhKxWxE4PqHriotKKfvzdfvjvo4LODyhZ+uKi0opmwuFDO09gbivtAhOT+jhnlSrTi5SStlUW08AY+K/jgs4PKHnZ3oBaO7stTgSpZQaWEs4P+VleOP+XI5O6JFvUHOHJnSllD01hfNTpACNJ0cn9NyMvrcwTR3aclFK2VOkg5CnCf3oPG4XWWkp2nJRStlWpODMy9Ae+qDyM72a0JVSttWiFXr08jK8/T0qpZSym6aOXlJcQlaqDlsclFboSik7a+7sJS/TG/elcyEJEnpuhodmvSiqlLKppo5e8hMwZBGSIKHnZ2iFrpSyr+YOP3mZ8b8gCkmQ0PMyvXT2Bun2B60ORSmlPqCpszchk4ogCRK6zhZVStlZS7iHngiOT+iRv3w60kUpZTehkKG506899GhFBuvrhVGllN20dQcIhoxW6NHSlotSyq6aOiPruOhF0ajkaUJXStlUpBWcqy2X6OSmRxbo0oSulLKXyLR/7aFHKcXtIifdo0voKqVsJ5FL50ISJHTouzDa1KkXRZVS9pLIpXMhWRJ6prf/rY1SStlFU4cfr9tFptedkOdLioSerysuKqVsqLmjl7xMT0IW5oIkSeh5mV7toSulbKc5gdP+IYqELiL3iki9iFQe4f4lItIqIuvDHzfFPsyjy8/09o/3VEopu7BdQgfuA84e5JhVxpi54Y+fnnhYxyY3w0O3P0RXry7QpZSyj6aO3oSNcIEoErox5nWgKQGxHLfIGE+dXKSUspPmzsQtnQux66EvFpEKEXlBRGYc6SARuVpEykWkvKGhIUZPfXBIkF4YVUrZRTBkaOlM3OYWEJuEvg4Ya4yZA9wBPH2kA40xK4wxZcaYssLCwhg8dR9dz0UpZTe+Lj8hk7gx6BCDhG6M8Rlj2sO3VwIeERl+wpEdA11CVyllN/2TipxUoYvICAkPshSRBeHHbDzRxz0WkSV0W3S2qFLKJhI9SxQgZbADRORhYAkwXET2Aj8GPADGmOXAxcC1IhIAuoBLjTEmbhEPICfdg4hW6Eop+2gK79GQyB76oAndGHPZIPffCdwZs4iOQ/8CXdpDV0rZRGSyoxNHuVhOp/8rpezk4OYWDuqh20VeplcrdKWUbTR39uJNcZHuSczCXJBMCT3Do/uKKqVso7mjbwx6ohbmgqRK6FqhK6Xso6nDn9ARLpBECT0/s6+HnuABNkopNaDmzt6EbQ4dkTQJPS/TS08gRJdfF+hSSlmvuSOxKy1CEiX0fJ0tqpSykb4KXRP6ccnV2aJKKZsIhgwtXX5ytUI/PpG/hI1aoSulLNba5ccYyM/QHvpxKRiWCkBje4/FkSilhroD4TwUyUuJkjQJvSir7xtX59OErpSyVp2vG4Di7LSEPm/SJPTM1BSyUlP6v5FKKWWVSGFZnK0V+nEryk6lvk0TulLKWpHCsihLK/TjVpydRr22XJRSFmto6yE7LYV0b+LWcYEkS+hFWanUaYWulLJYna+bogT3zyHJEnpxdhp1vh6d/q+UslSdrzvh/XNIsoRelJ1GbyBEa5dOLlJKWafO10NxgvvnkGQJPfIXUYcuKqWsYoyhvk1bLicsMuZThy4qpazS3OnHHzTacjlRkclF9W1aoSulrBEZOp3oIYuQdAldK3SllLWsmlQESZbQ071ustNSqNeErpSyiFXT/iHJEjocHLqolFJWiBSUhVlaoZ+w4uw0nVyklLJMna+H3AwPaZ7EzhKFJEzoRVmpOv1fKWWZOl93/wCNREu+hJ6dRn1bt84WVUpZor6tx5L+OSRhQi/OTsUfNDTrVnRKKQvU+7otGbIISZnQdeiiUsoaoZAJV+jacomJg9P/NaErpRKrqbOXQMhoyyVWIm919MKoUirRDm5soRV6TBT2T//XCl0plViRZUesWJgLkjChp3nc5GZ4dHKRUirh6vtnidq0QheRe0WkXkQqj3C/iMjtIlItIhtEZH7swzw2xVlp2kNXSiVcpJC0YpYoRFeh3wecfZT7PwlMCn9cDdxz4mGdmKLsVOp0xUWlVILV+brJy/CQmpL4WaIAKYMdYIx5XUTGHeWQZcADpm8mz2oRyRWRkcaYmlgFeayKstKorj9g1dMrB+nsDVCxp5X39reycb+PhvaDhUBWWgrTR2YzY1QOc0pyyc/0WhipcoI6n3WTiiCKhB6F0cCeQz7fG/7aBxK6iFxNXxVPaWlpDJ56YMXZqTS09RAKGVwuidvzKGcyxrBmZzOPl+/h+Xdr6OwNAn0/N6Ny04n8xOxu6mTlu7UAuF3CR6cUcUnZGD42tQiPO+kuP6kYsGqnoohYJPSBMuaA8+6NMSuAFQBlZWVxm5tfnJ1GIGRo6uxl+DBrelnKnl7f0sAtK6vYVNtGptfN+XNGcdbMEcwclTNg39PX7adqv49XNjfw1Lq9vFxVR1FWKtefOZmLTy7BrQWDOkS9r4cpxVmWPX8sEvpeoOSQz8cA+2PwuMft0MlFmtAVQHV9Oz9/fiOvbm6gND+DX148m3NmjyTDe/Rfgew0DwsnFLBwQgHfPnMyr21p4K5Xqvnek+9y/7928cNzp3HqxOEJOgtlZ8GQoaHd2pZLLN43PgtcER7tsghotbJ/DlCok4tUmDGGB9/aydLbV7F2ZzM3Lp3KS986nUvKSgZN5odLcbv4+LRinrz2VO64bB6tXX4++4e3+elfN9IbCMXpDJRTNHb0EAwZiiwasghRVOgi8jCwBBguInuBHwMeAGPMcmAlsBSoBjqBK+MVbLQiFXpNqw5dHMpau/zc8OQGXqisZcmUQn51yZyYvGMTEc6bM4pPTC/m1hc2ce+bOyjf1cQdl81jbEFmDCJXTlTXGp5UZNHCXBDdKJfLBrnfANfFLKIYGJGdRopL2NPcaXUoyiLbG9r54p/WsL+lixuXTuWqD0+I+QXyNI+bn5w/g8UTC/juExs49/Y3uPvy+Zw2qTCmz6OcYXdTX74pyU+3LIakvFSf4nYxJi+9/xushpb1e1q4ePlbdPQEePQri7n69IlxHe101owRrPz6aYzOS+dL963hmfX74vZcyr4i+aY0P8OyGJIyoQOU5GewRxP6kPPK5nouW7GazFQ3T1x7KiePzUvI847OTefRryxmfmkeX39kPf+7antCnlfZx+6mTvIzvWSleSyLIWkTeml+hlboQ8w/qur48v3ljB+eyZPXnsr44YntZ+eke7j/Swv45MwR/Pz5Ku56pTqhz6+staepkxILq3NI8oTe0umntUt3LhoK3qw+wLUPrWP6qGwe+coiyy5MpXnc3PnZ+SybO4pfvriZ+97cYUkcKvF2N3Va2m6BJE/ogLZdhoDynU1cdX854wsyuf/KBWRb+JYX+maV/uqSOZw5vZif/HUjj63ZM/h/Uo7mD4bY19JFqYUXRCGZE3pBX0LXtkty21zbxpV/WsOInDQevGoBeTZZb8XjdnHHZ+dx+uRCvvfUBv5WWWt1SCqOalq6CYYMY/OtHbaatAk90svShJ68Gtp6+NJ9a0j3uvnzVQstHf87kNQUN7+//GTmluTyjUffYcPeFqtDUnFycMiitlziIjvNQ16GRxN6kur2B7n6wXIaO3r44xdOYXSutW91jyTd62bF58soyEzlqvvLqWntsjokFQf9QxYLNKHHTakOXUxKxhi+/XgF7+xu4befmcusMTlWh3RUhVmp3PvFU+jsDfKl+8rp6AlYHZKKsd1NnXjcwggL13GBJE/oJfkZ7GrUhJ5slr+2nec21PC9s6dy9syRVocTlSkjsrjzs/PYXOvju09soG+CtUoWu5s6GJOXYfnqm0md0McWZLCvpYtAUBdOShb/qj7AL1/cxLmzR3LNRyZYHc4xWTKliO+ePZXn363h3jd3Wh2OiiE7DFmEJE/opfkZBENGF+lKEjWtXXzt4XeYUDiMX1w0GxHnrUX+ldMncOb0Yv57ZRVrdjZZHY6Kkd2NmtDjTke6JI/eQIjrHlpHtz/I8stPJjM1Fkv5J56I8KtPz6EkP4PrHlpHfZsWG07X2unH1x3QhB5vpZrQk8av/76ZdbtbuO3iOZxUNMzqcE5IdpqHey6fj6/bz7cerSAU0n66k9llyCIkeUIfmZOOxy16YdThXt/SwO9f387li0o5Z7YzLoIOZuqIbH5y3gzeqD7ACl3Iy9F2NXUAfdfsrJbUCd3tEsbk6dBFJzvQ3sO3HqtgcvEwfnjOdKvDianPnFLC0lkj+NWLm6nY02J1OOo4aYWeQCW66qJjhUKG6x+roK3bzx2XzSfN47Y6pJgSEf77wtkUZ6fx/x55h7ZuXUjOifY0dVKQ6WWYDa7rJH1CL83XjS6c6v63dvLalgZ+eM40poywbif1eMrJ8PC7S+eyp6mTnzy70epw1HHYbYNlcyOGQELPoLXLT2unVj9OUl3fxq0vbOJjU4u4fNFYq8OJq7Jx+Xz1oyfx5Lq9uoiXA9llDDoMiYTet/qZVunO4Q+G+NZjFWR43dx60SxHjjc/Vl/7+CRmjs7mxr+8S0Nbj9XhqCj5gyH2t3Tb4oIoDImErkMXnebOf1azYW8rt1w4y3YrKMaLx+3ifz49l/aeAN9/SpcGcIr9LV0EQ0ZbLokSWf1sZ2OHxZGoaFTsaeHOV6r51LzRfHJWcgxRjNak4iy+d/ZUXq6q57Fy3RTDCXY2Wr8x9KGSPqEPS01hdG46m2vbrA5FDaLbH+T6xysoykrlx+fPsDocS1x56jgWTyjgZ89Vsa9Fl9q1u821PgAmF9vjon3SJ3SAqSOy2BT+xiv7+u3LW6mub+fWi2aTk27tNnJWcbmE2y6eTcgYbnhSWy92t6mmjeLsVPJtslPWkEjo00Zms62hg55A0OpQ1BGs293Mite3cekpJXxkcqHV4ViqJD+DG5dOY9XWAzz8b2292FlVbRvTRmZbHUa/IZHQp47MIhgyVNe3Wx2KGkC3P8h3Hq9gRHYaPzhnmtXh2MLnFpbyoZMKuPn5jext1gv6duQPhqiub2PqCE3oCRX5hlfVaB/djv7npS1sa+jg1otmk5U2NFsthxMRfnHRbAC+p60XW9rW0I4/aJg20h79cxgiCX1cQQapKS421Wgf3W7e2d3MH1Zt57IFJZw+xFsthxuTl8GN50zjzepGHlmjrRe72RQuELVCT7AUt4vJxVls0pEuttITCPLdJzZQnJ3G95dqq2Ugl51SyuIJBdz8fBX7ddSLrVTV+vC4hQmFmVaH0m9IJHTQkS52dPs/trK1vp1bPjWLbG21DMjl6mu9BEOG7z/1rrZebGRTTRsnFWXhcdsnjdonkjibNjKbA+29Oq3aJir3tbL8te1cNH8MH51SZHU4tlZakMH3zp7Ca1saeGLtXqvDUWGban226p9DlAldRM4Wkc0iUi0iNwxw/xIRaRWR9eGPm2If6omZGv7Ga5Vuvd5AiG8/XkFBppebzk2uNc7j5YrF4zhlXB4/e24jdT7dts5qTR291Pl6mGaj/jlEkdBFxA3cBXwSmA5cJiID/RauMsbMDX/8NMZxnrCDI100oVvtrleq2VTbxs0XziInQ1st0eibcDSHnkCIH/xFWy9WiwywmOrACn0BUG2M2W6M6QUeAZbFN6zYy8/0Upyd2n9lWlmjqsbHXa9Us2zuKD4xvdjqcBxl/PBMvn3mFF6uqueZ9futDmdIq6q13wgXiC6hjwYOHTO1N/y1wy0WkQoReUFEBlyIQ0SuFpFyESlvaGg4jnBPzNQR2f0vhEo8fzDEd56oIDfDw0/OG5prtZyoL314PPNKc/nJX9+jvk1bL1bZVONj+DAvhVmpVofyPtEk9IEWoz78/d46YKwxZg5wB/D0QA9kjFlhjCkzxpQVFiZ+zPHUkVlU17fhD4YS/twKVry+ncp9Pn62bCZ5Nln7wmncLuGXF8+mszfIj56u1NaLRTbV2muGaEQ0CX0vUHLI52OA973fM8b4jDHt4dsrAY+IDI9ZlDEyfWQ2/qBhe4MupZtom2p9/PblLZwze+SQWxY31k4qyuKbZ0zmxffq+OuGGqvDGXICwRBb6tpsN8IFokvoa4BJIjJeRLzApcCzhx4gIiMkvK2MiCwIP25jrIM9UZG/qDrSJbH8wb5RLTnpHn62bKbV4SSFL582nrkludz0TKW2XhJsZ2MnPYGQMyt0Y0wA+CrwIlAFPGaMeU9ErhGRa8KHXQxUikgFcDtwqbHhe8EJhZl43S4q97VaHcqQsvzVbVTu8/HzC2baZplRp0txu/jVJXPo7A3yw79o6yWR3tvflz/sNsIFICWag8JtlJWHfW35IbfvBO6MbWix53G7mFOSw793NFkdypBRVePj9n9u5fw5ozh7prZaYumkomF8+8zJ3LJyE89W7GfZ3IHGKqhYe3tHE8NSU5hik00tDjVkZopGLJpQQOV+H23dfqtDSXq9gb7NnnPSvfzXEN2BKN7+48MTmF+ay03PvKcTjhLk7e2NnDIujxQbTfmPsF9EcbZoQgHBkKF8V7PVoSS93768haoaH7d+apaOaokTt0v41SVz6AkE+c4TusxuvNW3dbOtoYNFEwqsDmVAQy6hzy/Nw+MWVm+33TXbpLJ2VxPLX9vGZ8pKOEMnEMXVhMJh/GDpNF7f0sCf395tdThJ7e3tfe3ahZrQ7SHd62bOmNz+F0bFXkdPgG89VsGo3HR+dJ6u1ZIIly8ay+mTC7nl+Sp2HNBhufHy9o5GMr1uZo6y3wgXGIIJHfraLu/ua6W9J2B1KEnp5pVV7G7q5Defnsuw1Kiuu6sTJCLcdtFsvCkuvvnoegI6eS4uVm9v4pTx+bbsn8MQTegLJ+QTDBnWah895v5WWcv/vb2bq0+fwILx+VaHM6SMyEnj5xfMZP2eFn73j61Wh5N0DrT3UF3fzsLx9my3wBBN6CePzSPFpX30WKtp7eKGpzYwe0wO139iitXhDEnnzRnFxSeP4c5XqvXnO8YibdpFE+xbqAzJhJ7hTWFOSa7+wMdQMGT4xiPr6Q2E+N2l8/CmDMkfLVv4r/NnMK4gk28+up6Wzl6rw0kaq7eH++ejc6wO5YiG7G/dogn5vLu3lQ7to8fEPa9W8/aOJn66bCbjh9tnj8WhKDM1hdsvnceB9h5ueFLXTo+Vt3c0UjYu31Zbzh3OvpHF2cLxBQS0jx4Tb21r5DcvbeH8OaO4aL7OVrSDWWNy+M5ZU/jbe7Xc/6+dVofjeAfae9hS185CG7dbYAgn9Egf/S1tu5yQel83X3v4HcYNz+SWT80ivEabsoGrPjyBM6YVcfPKKtbt1sLlRBzsn9v3gigM4YSemZrCgvH5/K2yVt+SHqdAMMRXH36Hjp4A93zuZB2iaDMul/DrS+YyIieN6x5aR1OH9tOP1wuVNeRnepll4/45DOGEDrBs7ih2HOjgXV198bj88u+b+feOJm6+cCZTRthvoSIFORke7vncyTR29PL1R94hGNLi5Vi19wR4uaqOc2aNtHX/HIZ4Qj97xki8bpfuz3gcnlm/j9+/tp3PLizlU/PHWB2OOoqZo3P4r/NnsGrrAW57cZPV4TjOSxtr6faHWDZ3lNWhDGpIJ/ScDA9LphTy14r9Wrkcg4o9LXz3iQ0sGJeve4M6xGULSvn8orH8/rXtPLl2r9XhOMoz6/czOjed+aV5VocyqCGd0AGWzR1NfVsPb+vF0ajUtnbz5QfKGT4slXsun6/jzR3kpvOmc+rEAr7/1Ls6uitKje09rNp6gPPnjsLlsv8F/yH/2/jxaUVket08vX6f1aHYXkdPgKsfLKejJ8Afv1hGwTB77Xiujs7jdnH35+YzMjeNrzxYzp6mTqtDsr3n360hGDKOaLeAJnTSPG7OmjmCFypr6fYHrQ7HtnoDIa59aB2V+1r53aXzbLmfohpcboaXP36hjN5AiCvu/TcH2nusDsnWnlm/nynFWY75eR/yCR3ggrmjaesO8OrmBqtDsaVQyPCdJyp4fUsDt1w4S9c3d7iTirK494unUNPaxZV/WqOrjh7BnqZO1u5qZtk8Z1TnoAkdgFMnFjB8mJfHyvdYHYrtGGP4+fNVPLN+P985awqXLii1OiQVA2Xj8rn7c/PZWOPjKw+W0xPQd6eHezx88fi82ZrQHSXF7eILi8fxz031VOxpsToc2zDGcOsLm7j3zR188dRx/OeSiVaHpGLoY1OLue2i2bxZ3cg1D67VluMhWjp7+dMbO/jE9GJK8jOsDidqmtDDrvzwePIyPPzmpS1Wh2ILxhh++txGfv/6di5fVMpN507Xaf1J6KKTx3DLhbN4ZXMDX36gXJN62B9Wbae9N8D1Z062OpRjogk9bFhqCtd8ZCKvbWlgzc6hvT1dKGT48bPv8ac3d/LFU8fxs2UzHTFkSx2fzy4s5baLZ/NG9QG+dN8aOnuHdk/9QHsPf3pzJ+fMGumYi6ERmtAPccXicQwflsqv/77Z6lAs0+0P8rWH3+GBt3bx5dPG8+PztDIfCj5dVsKvL5nD6u2NfOb3q6nzdVsdkmWWv7qNbn+Qb5zhrOocNKG/T7rXzVc/OpHV25v4V/UBq8NJuAPtPVy6YjUrK2u4celUblw6TZP5EPKp+WNY8fkytjW0c+Fdb1JV47M6pISr83Xz4OpdXDhvDCcVDbM6nGOmCf0wly0sZVROGj99buOQ6idW7mvlgrveZFOtj3s+dzJXnz5Rk/kQdMb0Yh77ymKCxnDJ8rf4W2Wt1SElTGREVzBk+PrHJ1kdznHRhH6Y1BQ3N184i021bfzw6cqkX1rXGMMf39jBhXe/SSBoePTqxZw9c4TVYSkLzRydwzPXfZgJhZlc8+e1/OjpyiFR3Dzw1i7+WrGfb5wxidIC54xsOZQm9AF8dGoR/+/jk3hi7V4e/nfyjk2v93XzH/eX87PnNvKRyUW88PXTmFOSa3VYygZG5KTxxDWn8uXTxvPg6l1ccNebbNyfvC2Ytbua+NlzGzljWhH/ueQkq8M5bmJVBVpWVmbKy8stee5oBEOGK+9bw+ptjTx+zeKkSnSBYIj739rF/7y0hd5giB8sncYVi8dqi0UN6NXN9Xz78QqaO/18YfE4vvmJSWSleawOK2Ya2no4945VpHncPPvVD5OTbu9zE5G1xpiyAe/ThH5kzR29nHvHG3T5g9z9ufm2335qMMYYXt3SwC9e2MSm2jY+Mrmwb4d43dRZDaKls5fbXtzMw//ezfBhqXz7zMlcOG+M41fb3NXYwZcfKGd3UydPXfshpo+y/zBFTegnYHtDO1c9UM7uxk5+dO50R1ayoZDh7xvruPOVrVTu8zE6N50fnTuNs2aMcNy5KGut39PCj5+ppGJvK6Nz07lmyUQuOXkMaR631aEds9e2NPC1/1uHyyXc9dn5fOik4VaHFBVN6CfI1+3nW4+u5+WqepbNHcX1n5jiiIsmuxo7eHLdPp5cu5d9LV2MLcjguiUnccG80Y6vrJR1jDG8urmB2/+5lXd2t5CT7mHZ3FFcfPIYZo3OsX2R0NzRyx9WbWf5a9uYXJzFH64oc9T0/hNO6CJyNvA7wA38rzHm1sPul/D9S4FO4IvGmHVHe0wnJXToq3Jv/+dW7n5lG4FQiPPnjOKq0yYwY1S2bX6Au/1B1u1uZtXWA6za2kDlPh8icNqkQj5dNoazZ4wgxeZ7IirnMMbw1vZGHl2zh79V1tITCDG2IIPTJxVy2qThLJxQYKt+9J6mTh54aycPvb2bzt4gF84bzc0XziTD66zNzU8ooYuIG9gCfALYC6wBLjPGbDzkmKXA1+hL6AuB3xljFh7tcZ2W0CPqfN388Y0d/Hn1Ljp7gwwf5mXxxOGcMi6PsQWZjMlLZ1ROOmkeV8wTfSAYwtcdoLXLT72vm30tXexv6aK6vp2NNT62NXQQDBlSXML8sXl8dEoRF8wbxcic9JjGodThfN1+nt9Qw8sb63hreyOdvX3DHEvy05kxMofJxcMYnZfO6NwMRuSkkpPuJSfdE5d3ir2BEHW+bvY2d7GnqZN39jTzr22N7GrsxO0Szp8zimuXTGRysTM3Nj/RhL4Y+Ikx5qzw598HMMb89yHH/B541RjzcPjzzcASY0zNkR7XqQk9ormjl5fCP7xvVh+gvu39GwWkuIQMr5sMbwopbiHFJbhcQiTFixy8HXkFjDEYwBgIGUMwZPAHQ/QEQvT4Q3QdYSzwiOw0ZozKZvqobOaMyWXRxAKGpTqr6lDJozcQYu2uZtbtbmZjjY+N+33sbOxgoFSTmuIizePGm+LC63bhdglulyACAv1F0ZF+V0IhQ9AYgkFDpz9IZ0+Q3mDofc+RlZrCwgkFnDqxwHGrJw7kaAk9mt/60cChg7H30leFD3bMaOB9CV1ErgauBigtdfa62nmZXj59SgmfPqUEYwy14Ypgb3Mnta09tPf46egJ0tETIBj+oQtENqI2YHj/T3d/ehdwieAO/+txu0j19P2wD0tLISfdQ26Gh4LM1HDFk+7IC1IqeXlTXCyeWMDiiQdHhfmDIWpb+95V1vm68XX5aen009YToDcQoicQpCcQIhQyhAwEjenP3kf7XXHLwWIpUkANS3VTlJXGmLz0/t+RodJqjCahD9Q3OPxvbTTHYIxZAayAvgo9iud2BBFhZE46I3PSOWVcvtXhKGU7HreLkvwMx1fHdhfNn629QMkhn48B9h/HMUoppeIomoS+BpgkIuNFxAtcCjx72DHPAldIn0VA69H650oppWJv0JaLMSYgIl8FXqRv2OK9xpj3ROSa8P3LgZX0jXCppm/Y4pXxC1kppdRAohoKYYxZSV/SPvRryw+5bYDrYhuaUkqpYzE0Lv0qpdQQoAldKaWShCZ0pZRKEprQlVIqSVi22qKINAC7jvO/DweSZRdnPRd7SpZzSZbzAD2XiLHGmMKB7rAsoZ8IESk/0loGTqPnYk/Jci7Jch6g5xINbbkopVSS0ISulFJJwqkJfYXVAcSQnos9Jcu5JMt5gJ7LoBzZQ1dKKfVBTq3QlVJKHUYTulJKJQlHJHQRuURE3hORkIgccaiPiOwUkXdFZL2I2HJ/u2M4l7NFZLOIVIvIDYmMMVoiki8iL4nI1vC/eUc4zpavy2Df4/By0LeH798gIvOtiDMaUZzLEhFpDb8G60XkJiviHIyI3Csi9SJSeYT7nfSaDHYusX9NjDG2/wCmAVOAV4Gyoxy3Exhudbwnei70LVO8DZgAeIEKYLrVsQ8Q523ADeHbNwC/cMrrEs33mL4loV+gb0euRcDbVsd9AueyBHjO6lijOJfTgflA5RHud8RrEuW5xPw1cUSFboypMsZstjqOWIjyXBYA1caY7caYXuARYFn8oztmy4D7w7fvBy6wLpRjFs33eBnwgOmzGsgVkZGJDjQKTvl5GZQx5nWg6SiHOOU1ieZcYs4RCf0YGODvIrI2vCG1Ux1p0227KTbhnanC/xYd4Tg7vi7RfI+d8jpEG+diEakQkRdEZEZiQos5p7wm0YrpaxLVBheJICIvAyMGuOsHxphnonyYDxlj9otIEfCSiGwK/5VMqBicS1SbbifC0c7lGB7GFq/LYWK2+bkNRBPnOvrWAGkXkaXA08CkeAcWB055TaIR89fENgndGHNGDB5jf/jfehH5C31vRROeOGJwLrbZdPto5yIidSIy0hhTE37bW3+Ex7DF63KYZNr8fNA4jTG+Q26vFJG7RWS4McZpi1055TUZVDxek6RpuYhIpohkRW4DZwIDXl12gGg25raDZ4EvhG9/AfjAuw8bvy7JtPn5oOciIiNERMK3F9D3u9+Y8EhPnFNek0HF5TWx+kpwlFeLL6TvL3MPUAe8GP76KGBl+PYE+q7uVwDv0dfesDz24zmX8OdLgS30jV6w67kUAP8Atob/zXfS6zLQ9xi4BrgmfFuAu8L3v8tRRlhZ/RHFuXw1/P2vAFYDp1od8xHO42GgBvCHf0/+w8GvyWDnEvPXRKf+K6VUkkialotSSg11mtCVUipJaEJXSqkkoQldKaWShCZ0pZRKEprQlVIqSWhCV0qpJPH/AUFjsqsnSa/eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def potential(x):\n",
    "    return (x - 1) ** 4 * (x + 1) ** 4\n",
    "\n",
    "xs = np.linspace(-1.5, 1.5, num=100) # We will sample 100 discrete bins\n",
    "plt.plot(xs, potential(xs));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c8bb85b",
   "metadata": {},
   "source": [
    "We will use umbrella sampling to sample the potential. Define a sampling axis and the umbrella centers. These umbrella centers define our thermodynamic states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2c070f24",
   "metadata": {},
   "outputs": [],
   "source": [
    "bias_centers = np.linspace(-1.5, 1.5, 20) # 20 umbrella's\n",
    "\n",
    "# we use a harmonic bias potential\n",
    "def harmonic(x0, x):\n",
    "    return 2 * (x - x0) ** 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97c77830",
   "metadata": {},
   "source": [
    "For each of the bias centers, contruct the bias potential as a sum of the original potential and the harmonic bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d20aaf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "biased_potentials = []\n",
    "for bias_center in bias_centers:\n",
    "    bias = lambda x, x0=bias_center: harmonic(x0, x) + potential(x)\n",
    "    biased_potentials.append(bias)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebd15f5f",
   "metadata": {},
   "source": [
    "Now we will sample a discrete trajectory for each of the thermodynamic states. We use a deeptime MSM to perform the sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cc53313b",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 10000  # trajectory length\n",
    "\n",
    "trajectories = np.zeros((len(bias_centers), n_samples), dtype=np.int32)\n",
    "\n",
    "for i, U in enumerate(biased_potentials):\n",
    "    biased_energies = U(xs)  # get the potential at each of the bins\n",
    "\n",
    "    transition_matrix = tmatrix_metropolis1d(biased_energies)\n",
    "\n",
    "    msm = MarkovStateModel(transition_matrix)\n",
    "    trajectories[i] = msm.simulate(n_steps=n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "584a8429",
   "metadata": {},
   "source": [
    "Plot the sampled trajectories and biased potentials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8831eb32",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAO2klEQVR4nO3db4hd+V3H8ffHSJ4UrbZJ3ZpsmlSDdQsVyjRttWgVV7KLkC5WTBVbtSVEiOIDwYBQVny0PhLp1hAkWJ8YBG0NdZZYC1KhVpItu9tmbdrpWsmYtZtuy5bF0t20Xx/MLb2dvXfumcm9c8/93fcLhrnnnF/OfH9zZj7zu7/zJ6kqJEmL7/vmXYAkaToMdElqhIEuSY0w0CWpEQa6JDXi++f1hfft21eHDx+e15eXpIX06KOPfqWq9o/aNrdAP3z4MFevXp3Xl5ekhZTkv8dtc8pFkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaMbc7RSU16sGXj1n/3O7WsYQMdEl3blyIa1c55SJJjTDQJakRBrokNcI5dEk747x57xjokrozxHvNKRdJaoSBLkmNMNAlqRHOoUva2rTmzYf3412jM9FphJ7keJLrSdaSnB2x/e1Jnkvy2ODj/dMvVZK0lYkj9CR7gIeBe4F14EqSS1X15Kam/1ZVvzyDGiVJHXQZoR8D1qrqqap6AbgInJhtWZKk7eoS6AeAG0PL64N1m701yeNJHkny+lE7SnIqydUkV2/durWDciVJ43Q5KZoR62rT8qeB11TV80nuBz4CHH3JP6o6D5wHWFlZ2bwP7YQnmiQNdAn0deDuoeWDwM3hBlX19aHXq0k+mGRfVX1lOmXqe3i3nmbNn7GF1GXK5QpwNMmRJHuBk8Cl4QZJ7kqSwetjg/0+O+1iJUnjTRyhV9XtJGeAy8Ae4EJVXUtyerD9HPBO4HeT3Aa+AZysKqdUJE3mtOHUdLqxqKpWgdVN684Nvf4A8IHplqZt8xdDWmreKSpp9zlHPxMG+qLwF0DLwHeZd8RA7xt/oCXtkIHeqs0jev84SM0z0PvMaRZJ2+Dz0CWpEY7Q+2A3RuLOzUvNc4QuSY1whL6Mxr0jcOSuPtnuu0rfhRrokhaYFw58D6dcJKkRjtAlOdJthIEuqT1LOp/ulIskNcIRuqT+c0qoEwN9Xvr4A7qkb1OlVjjlIkmNMNAlqREGuiQ1wjl0jeZ8urRwDHRNZrhrkS3Rz6+Bvpv6eGWLpGY4hy5JjTDQJakRTrloe5ZoPlJaNI7QJakRjtBnzROhknaJI3RJaoQjdGlZLeO7x8bPAXUaoSc5nuR6krUkZ7do96Yk30ryzumVKEnqYmKgJ9kDPAzcB9wDvCvJPWPaPQRcnnaRkqTJuozQjwFrVfVUVb0AXAROjGj3e8DfA89MsT5JUkdd5tAPADeGlteBNw83SHIAeAD4BeBN43aU5BRwCuDQoUPbrVWSpqfB+fQugZ4R62rT8p8Df1RV30pGNR/8o6rzwHmAlZWVzftoxzKebJI0d10CfR24e2j5IHBzU5sV4OIgzPcB9ye5XVUfmUaRc9fgX3JJ7ekS6FeAo0mOAP8DnAR+fbhBVR35zuskfw18tJkw32xcuDsqlzRnEwO9qm4nOcPG1St7gAtVdS3J6cH2czOusb8McUk90unGoqpaBVY3rRsZ5FX1W3deVg8Y1pM5FSX1irf+S1IjvPVfWia+82yaI3RJaoSBLkmNcMpF0+fJUmkuDPRhzi9KWmAGuqbDP4bS3DmHLkmNcIQuSY2c93GELkmNcISu2Wpk5CMtAkfoktQIA12SGmGgS1IjnEOXWuc9AkvDEbokNcJAl6RGGOiS1AgDXZIaYaBLUiO8ykW7x7tGpZlyhC5JjTDQJakRBrokNcI5dEkatsDnegx0b4uW1AinXCSpEY7QpRb5znMpOUKXpEY4Qtd8LPCJJ6mvOo3QkxxPcj3JWpKzI7afSPJEkseSXE3ytumXKknaysQRepI9wMPAvcA6cCXJpap6cqjZx4FLVVVJ3gD8HfC6WRQsSRqtywj9GLBWVU9V1QvAReDEcIOqer6qarD4MqCQJO2qLoF+ALgxtLw+WPc9kjyQ5HPAPwG/M2pHSU4NpmSu3rp1ayf1SpLG6BLoGbHuJSPwqvpwVb0OeAfwp6N2VFXnq2qlqlb279+/rUIlSVvrEujrwN1DyweBm+MaV9UngB9Lsu8Oa5MkbUOXQL8CHE1yJMle4CRwabhBkh9PksHrNwJ7gWenXawkabyJV7lU1e0kZ4DLwB7gQlVdS3J6sP0c8CvAu5O8CHwD+LWhk6SSpF3Q6caiqloFVjetOzf0+iHgoemWJklztvkRCj2/Cc47RTV/3jUqTYWBLrXCB3ItPR/OJUmNMNAlqREGuiQ1Yjnn0J1rlNQgR+iS1AgDXZIaYaBLUiMMdElqxHKeFFV/edfo9niCX0McoUtSIwx0SWqEgS5JjTDQJakRnhSVpK56ftLeEbokNcIRurRovFRRYxjo6q+ev72V+sYpF0lqhIEuSY1wykVaBM6bqwMDXYvB+XRpouUJdEc4khq3PIEuLQLfiegOGOhaPMsSer6r1DYZ6NK8GdyaEi9blKRGGOiS1IhOgZ7keJLrSdaSnB2x/TeSPDH4+GSSn5p+qdIID778ux/Skps4h55kD/AwcC+wDlxJcqmqnhxq9l/Az1XV15LcB5wH3jyLgqUm+Ado8fXw5HyXEfoxYK2qnqqqF4CLwInhBlX1yar62mDxU8DB6ZYpSZqkS6AfAG4MLa8P1o3zXuCROylKkrR9XS5bzIh1NbJh8vNsBPrbxmw/BZwCOHToUMcSJUlddBmhrwN3Dy0fBG5ubpTkDcBfASeq6tlRO6qq81W1UlUr+/fv30m9kqQxugT6FeBokiNJ9gIngUvDDZIcAv4B+M2q+vz0y5QkTTJxyqWqbic5A1wG9gAXqupaktOD7eeA9wOvBD6YBOB2Va3MrmxpAXlli2as063/VbUKrG5ad27o9fuA9023tCnwF0jSEvFZLmpHD68LlnaTgS5Nm39YNCcGujRLTvtpFxnoWi6bA9YRtBri0xYlqRGO0LXcnO9WQwx0tc95bC0Jp1wkqRGO0NUmR+VaQga6NA3+AVEPGOjSd3iCVAvOOXRJaoSBLkmNcMpFku5UT6brDHRpkp78skqTGOjSKOOuWjHc1WMGurRTXqqonvGkqCQ1wkCXpEYY6JLUiPbm0J3XlLSkHKFLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjOgV6kuNJridZS3J2xPbXJfn3JN9M8ofTL1OSFsSDL//uxy6beOt/kj3Aw8C9wDpwJcmlqnpyqNlXgd8H3jGLIiVJk3UZoR8D1qrqqap6AbgInBhuUFXPVNUV4MUZ1ChJ6qBLoB8Abgwtrw/WbVuSU0muJrl669atnexCkjRGl0DPiHW1ky9WVeeraqWqVvbv37+TXUiSxugS6OvA3UPLB4GbsylHkrRTXQL9CnA0yZEke4GTwKXZliVJ2q6JV7lU1e0kZ4DLwB7gQlVdS3J6sP1ckruAq8APAt9O8gfAPVX19dmVLkka1ul/LKqqVWB107pzQ6//l42pGEnSnHinqCQ1wkCXpEYY6JLUCANdkhphoEtSIzpd5dJrc3iimST1kSN0SWqEgS5JjTDQJakRBrokNcJAl6RGLP5VLpLUV8NX4T343My/nCN0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxGLeKeoz0CXpJRyhS1IjFnOELkmLZhee6+IIXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDWiU6AnOZ7kepK1JGdHbE+SvxhsfyLJG6dfqiRpKxMDPcke4GHgPuAe4F1J7tnU7D7g6ODjFPCXU65TkjRBlxH6MWCtqp6qqheAi8CJTW1OAH9TGz4F/FCSV0+5VknSFrrcKXoAuDG0vA68uUObA8DTw42SnGJjBA/wfJLr26r2u/YBX9nhv+0b+9JPrfSllX5AS335k9xJX14zbkOXQM+IdbWDNlTVeeB8h6+5dUHJ1apaudP99IF96adW+tJKP8C+dNFlymUduHto+SBwcwdtJEkz1CXQrwBHkxxJshc4CVza1OYS8O7B1S5vAZ6rqqc370iSNDsTp1yq6naSM8BlYA9woaquJTk92H4OWAXuB9aA/wN+e3YlA1OYtukR+9JPrfSllX6AfZkoVS+Z6pYkLSDvFJWkRhjoktSIhQj0JL+a5FqSbycZe6lPki8l+UySx5Jc3c0au9pGX7Z83EIfJHlFko8l+cLg8w+PadfL49LSIy069OXtSZ4bHIPHkrx/HnVOkuRCkmeSfHbM9kU6JpP6Mv1jUlW9/wB+EvgJ4F+BlS3afQnYN+9677QvbJx8/iLwWmAv8Dhwz7xrH1HnnwFnB6/PAg8tynHp8j1m40T/I2zcZ/EW4D/mXfcd9OXtwEfnXWuHvvws8Ebgs2O2L8Qx6diXqR+ThRihV9V/VtVO7yrtlY596fK4hT44AXxo8PpDwDvmV8q2tfRIi0X5eZmoqj4BfHWLJotyTLr0ZeoWItC3oYB/TvLo4DEDi2rcoxT65kdqcL/B4POrxrTr43Hp8j1elOPQtc63Jnk8ySNJXr87pU3dohyTrqZ6TLrc+r8rkvwLcNeITX9cVf/YcTc/U1U3k7wK+FiSzw3+Su6qKfSl06MUdsNWfdnGbnpxXDaZ2iMteqBLnZ8GXlNVzye5H/gIG09HXTSLcky6mPox6U2gV9UvTmEfNwefn0nyYTbeiu56cEyhL715lMJWfUny5SSvrqqnB297nxmzj14cl01aeqTFxDqr6utDr1eTfDDJvqpatIddLcoxmWgWx6SZKZckL0vyA995DfwSMPLs8gLo8riFPrgEvGfw+j3AS9599Pi4tPRIi4l9SXJXkgxeH2Pjd//ZXa/0zi3KMZloJsdk3meCO54tfoCNv8zfBL4MXB6s/1FgdfD6tWyc3X8cuMbG9Mbca99JXwbL9wOfZ+Pqhb725ZXAx4EvDD6/YpGOy6jvMXAaOD14HTb+c5cvAp9hiyus5v3RoS9nBt//x4FPAT8975rH9ONv2Xjs9ouD35P3LvAxmdSXqR8Tb/2XpEY0M+UiScvOQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmN+H9N41Qh3C6geQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(xs[trajectories.flatten()], bins=100, density=True, color='C1',\n",
    "         label='Histogram over visited states');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35fe092c",
   "metadata": {},
   "source": [
    "Now use the deeptime KMeans module to cluster the data into 2 states (using the knowledge that we are working with a double-well potential). For each sample, the Markov state that sample was in is stored in `dtrajs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e402203f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[23.82456604]\n",
      " [76.19136811]]\n"
     ]
    }
   ],
   "source": [
    "estimator = KMeans(\n",
    "        n_clusters=2,  # place 2 cluster centers\n",
    "        init_strategy='uniform',  # uniform initialization strategy\n",
    "        max_iter=10,\n",
    "        fixed_seed=13,\n",
    "        n_jobs=8\n",
    "    )\n",
    "\n",
    "clustering = estimator.fit_fetch(trajectories.flatten())\n",
    "dtrajs = clustering.transform(trajectories.flatten()).reshape(\n",
    "        (len(trajectories), n_samples))\n",
    "\n",
    "print(clustering.cluster_centers)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb31c79",
   "metadata": {},
   "source": [
    "The final piece of input we need are the bias matrices, $K$ matrices, one for each thermodynamic state. \n",
    "\n",
    "The bias matrix of thermodynamic state $k$ is of shape $(N_k, K)$, $N_k$ being the number of samples taken at state $k$. `bias_matrices[k, i, l]` is the potential of the $i$-th sample taken at state $k$, **evaluated at** state $l$. \n",
    "\n",
    "In symbols, when $x_i^k$ is the $i$'th sample taken at thermodynamic state $k$, `bias_matrices[k, i, l]` $= b^l(x_i^k)$.\n",
    "\n",
    "This means that for each sample taken __at any given state__, the bias potentials of that sample __in all states__ need to be stored!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "57200274",
   "metadata": {},
   "outputs": [],
   "source": [
    "# move the range of the trajectories back to the range of xs: (-1.5, 1.5)\n",
    "trajectories_shifted = trajectories / 100 * 3 - 1.5\n",
    "\n",
    "bias_matrices = np.zeros((len(bias_centers), n_samples, len(bias_centers)))\n",
    "for i, traj in enumerate(trajectories_shifted):  # enumerate the states (=trajectories)\n",
    "    for j, bias_function in enumerate(biased_potentials):  # for each term. state...\n",
    "        # ... store the bias potential of all samples as if they were sampled at that thermodynamic state\n",
    "        bias_matrices[i, :, j] = bias_function(traj) - potential(traj) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d31289e0",
   "metadata": {},
   "source": [
    "Construct a `TRAMDataset` object from the discrete trajectories and the bias matrices. This object will construct the state and transition count matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "049722d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 10000)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset = TRAMDataset(dtrajs=dtrajs, bias_matrices=bias_matrices)\n",
    "bias_matrices.shape\n",
    "dtrajs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d3a8357",
   "metadata": {},
   "source": [
    "We can use the dataset to reduce our data to a chosen submodel, or to the largest connected set according to chosen criteria. One could also choose to reduce the data to a chosen submodel, for example if one wants to compute the connected set by a different method than the ones available in the TRAM module. In this case, we want to reduce our dataset to the largest connected set by the BAR variance connectivity criterium. This may take a while, so we pass a progress bar to show the progress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "64170f26",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1fa67e5139b24baf9b8165200e52b4af",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/40 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm.notebook import tqdm\n",
    "dataset.restrict_to_largest_connected_set(connectivity='BAR_variance', progress=tqdm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef201496",
   "metadata": {},
   "source": [
    "Now we have gathered all our input, we can run a TRAM estimate. Because we set maxiter to 100, this might give us a warning that TRAM did not converge (100 iterations is enough for demonstration purposes). The model we obtain contains the data from the connected set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "597f8ea7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "985ca6cd703c46ee8ebd46a2de8c7665",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b397a4787934d90b0fe1f1421b7e34e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\users\\maaike\\documents\\phd\\deeptime\\deeptime\\deeptime\\markov\\msm\\tram\\_tram.py:256: ConvergenceWarning: TRAM did not converge after 100 iteration(s). Last increment: 0.0005462764459391867\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "tram = TRAM(lagtime=50, maxiter=100, progress=tqdm, callback_interval=3, init_strategy=\"MBAR\")\n",
    "model = tram.fit_fetch((dtrajs, bias_matrices))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "903cf917",
   "metadata": {},
   "source": [
    "We can plot the `biased_conf_energies` in a contour plot to see what they look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9ca4f380",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAADCCAYAAACVO1RNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfn0lEQVR4nO2de7AdRZ3HPz/uDfL2AiGBBSLRjQhLCSriA12DrBgiwrqFCj4Q1FUsUSylRLRESootLHyQXQRBwEjJwy0FYRV5lFsIlgZ5FE+DGBEhhE14v0tzwm//6DnJuXPnzPT0dM/0nDPfqluZmX78+pycz/x6erp/LapKp06dRk8bNd2ATp06hVEHd6dOI6oO7k6dRlQd3J06jag6uDt1GlF1cHfqNKKabLoBWdpGRHdOXZuc2HAsE+nEIcfpc9ty6TomctJy6l83KeuPewMZ16UyDp6n0/LK2dbZSzU6z960cuuml9PeQN4e2cdNp5F+tbs2p2A/bcWjqrodBdpVRJ/LSX8IrlbVRUX11KUo4Z4/Acu2MMeztk0lbjNwbJu2TSqfS1rKlqbKPbnNphuOmZqelnP+xIy0rYeWeyzViMH0dN4npqVtPS0trz1PPr3heBJ44dGBvI9suFklDZquR4ccZ50/4ljPI6nzx4B+s3RtKnH1wPGanLTFf8VCzwGfzUk/Hmbb1FOXooRbJhKoQ8Kal1YCZMiHLA86k+4Gc/o8r96isoNAQwpoyIe6CGJXqIvKpsvngj2eihJuNgEWJMcBYPUFss15KJiL6rZqax7UZby0zXkemGXLp8sWgp3nte01C5jjVLIZxQn3JAZAD7DGAnJWfogMaIgb6qzyncceqjjhfgkbPHdNIIN/mLPK+IbZpo7SQJuGTldVqF3qLALbSkVee3RvBnHCPYmB0+WZm5kgQ3mvnHUtBMy+7BQBDQ5eGspDDeW9ddY1G7BLd8eraRKYW6G8iJwPHASsUdU9kmsnAf/Ohk/8ZVW9MqPsImAJ5p3Muap6qk1749Og5w4Acta1rDw+YPZly+ozuAANbvC5QO2zHqfueONeeilwBnBB6vp3VPWbwwqJyATwXeAdwErgJhG5QlX/kGcsTriznrmxAxn8eGWTxw/MrvZs6gGLbjfUC7Vr3Vl1OYNt47XrhV1VrxeRXRyK7gOsUNX7AETkEuAQoBrcQ7oSPwZ2TbJMAU+q6l4ZZe8HngHWAT1V3dvqo0yCLph52RfI4OaVwa2bnWXP1qYNzBAY6GHXQkKdVRdUALs6yLMo7JbPFpGbB87PUdVzLKo+RkSOAG4GvqCqT6TSdwQeHDhfCbyhqFIbz72UVFdCVd/fPxaRbwFP5ZTfT1Wz/juHa8Kta511rW6PXMXm0Gs+YYbwQIeobwbUUA3sIF77UWsHtkFnASdjptadDHwL+GgqT8Z/7IypeDNUCHdeV0JEBHgf8Paiesrob5OzouxaD7vuG2aoCehh131DXbVOa7DbJ1Vd/0FE5PvAzzOyrQQGZ2TvBKwqqrvqM/dbgdWq+qch6QpcIyIKnG3ZRWEdE6VfEZk8zYBsazvXli3M0CzQPmyVqbcU2GG64yElIjuo6sPJ6XuAuzKy3QQsEJH5wEPAYcAHiuquCvfhwMU56fuq6ioRmQNcKyL3qOr1WRlF5BPAJwDmzts4d2oluE0QGZYv77rtTWRYHcFghjBAQxxQQ0Cw3WeoVXwVdjGwEPNsvhL4GrBQRPbCOML7gU8mef8B88prsar2ROQY4GrMq7DzVfXuInvOcIvIJPBvwOuG5VHVVcm/a0TkMsyoXybciVc/B+Af936plln4AOFAzrKfV08ZmMGDd4ZygOVdt+0m+6p7WP2ZUEM5sOOTqh6ecfm8IXlXAYsHzq8EZrz/zlMVz/0vwD2qujIrUUQ2BzZS1WeS4wOAr9tU3GNyPVAxg5x73QfM0AzQPu168dZQHmy/XruNsnkVNqMroarnYfr9F6fyru9KYHowl5kxNyaBi1T1KptGDT5zlwE5L60pmKFBoPPSQkPtYiNyjz05AXO2yMmQ986oAdmMlmd1JVDVIzOure9KJC/c93RpVI/J9VDXAXJeXbn2fXlnaBZo3/bLQg2ewe68NkQ6Q81mtLyvaEGG8p4Z3KB1AS2vHS5tyLPlBDX4BXv81Bq4oR6QnbrX4OaVwR0k30DXbc8JanAHuzr064OIDFPbuuVNqWmQoUUwQ1xAF9msHezxVJRw9wP6tQJkcIenKD0UXFXaFAxqcAe7SOMJfpRwm2758HBEhdergAzuXhmqgVOUXgWsorqr9AxsbDt7aygGu0p3vMRNo79acZjus6+qDkUM99S0a42ADM3CDGGBrsN+JW8N1cAeb7UGbvAAMjQPs02e0EDX0YZCqCE82OMNfpRwQzHI4MErQz2g2NRRFaaY2lHZW0M9Hrvkc3w//JejhsRGOBkTeOHFpEFH9qdtD5TbGbPkevsk3zmqusSmudGpt25yRsNqARnq8cxQDBHEA7RtW2rx1jZ12OapXUuZGWbpNFX9KoCIfBY4ETg6Va6HCeJwq4hsCdwiIte2MsyS9ibCPyuDP1Bjgtk2j017oEaoYcTBzoyNoKpPD5xuTkYQhmRJ6MPJ8TMishwTnaV9cM9Q20AGfzDb2myiXd6ghvpjnjm8WisaLXcMsyQipwBHYKbB7FeQdxfgNcCNRfXGCXePsK+jbPOUyVc3zGXy+fTSYAk1+PPWtnWVyRdELmGWUNWvAF8RkROAYzDrvGdIRLYAfgp8LuXxMxUn3INqK8h9jQrQ4BlqGDGwfegi4BdkwC0iszBgX6iql9pU5hr99CQCBVIHjOd2WdRgk+6St0mYy+ZtBdTgH+wyiie4g4gsGAhTdjBwT0YewQR1WK6q37at2yn6aaIggdSnqQmQIQzMZdsQQ3sbhbpsvTV47QmqvgrLCrO0WER2xbzi+ivJSHkqNsK+wIeBO0XktqS6TIc6qErRTwvkFEgdMJ67TMyvLMUAR9l2xNJma6hhbMD2INcwS6r6G7LDG+eqyjO310DqgwES2XZeOCigHBgQzjuXzR+y3aWAhvJAhQK7rOLpkoeWK9zeA6kPBkiUnffOD7geCxAQFmYI3/6ooHapv0avXfwqLCo5wR0ykPoMhYQhNMguZcrCDDUADeGhdrHRju54U3KCO2QgdWD4M3dasXll1zJ1AA01QQ31gN2pSE7RTwkYSH2oQoMM9cEMkQMN9UHtaquBm0HFhSN1yzX6abBA6usVsnsN7lDWCTO4fTaoGWpoB9jjM5gGsc5Q6xWktwFmaBHQ0A6oO5VRnHAPyvWH3iaYwf1zQsugrmKzatmKGofR8uDKm36apSpAVi3fSqChGiRjCHYLFSfcRaoCY9UbQRWYoRrQ0JCX7qvqM2sHZ52KH+62euW+GoW5r6ah9tGGCG4M1eeWZy3Cei9wErAbsI+q3pxTfgIzI/QhVT2oyF6ccNu+506raa8M1WGGSIAGf6PLEYAZh5YycxHWXZitsM+2KH8ssBzYysZYnHDbqmmvDH5gBk9Agz+QYvDWvutpVkPCLC0HSHbDHSoR2Ql4F3AK8Hkbe+2BOwavDP5ghgiBhs5bR6vTgS8CW9oWiBfuWGCGSIGGOKEGv+2K6CYRKIZakUSk/5x+i4gstC0XJ9xFk1jS8gky+IUZIgYa/M/aigjGGZpD4FlqTjHULLQvcLCILAY2AbYSkR+p6ofyCsUJd5HGCua+Yoca/Lcx5htFfVLVE4ATABLPfVwR2NAGuH2DDP5hhpYADe2AOlJV33EkaxHW48B/AdsBvxCR21T1nakwS87NLWpQ1ru504B3A38H/gwcpapPZpS9H3gGWAf0rLssPeJ9Zu4rCMwQDpRQ3dFQ7Q1Vb/Cu+VANWYQFcFlG3mmLsAauXwdcZ2NvI4s8S4FFqWvXAnuo6quBe0m6DEO0n6ruFehZZKYey/irKl0788+rVg/8+dYa2gd2aM1pugG1yClAoqpeM3C6DDjUc7vsFMIjQ0Cv3FcdUIT2Tm0Fuy8HDz4BOmYLRz4K/HhImgLXiIgCZ1d6LdBakPsaBaD7ajvYfTXXRa9DleAWka9gnpAvHJJlX1VdJSJzgGtF5B5VvX5IXRuin240r4O5lOr8gY4K2H2NbhfdGW4R+QhmoG1/Vc2MVtrfZ1hV14jIZZhY5plwT4t+Kntr+SjNWZXWBXJfdf/wO6jr1LpJ4cltNsnJ8UJtbbGRa4DERcDxwNtU9fkheTYHNkq2HN0cOAD4unNLbVQ7zDDaQPfVgd1GuQZIPAF4CaarDbBMVY9OvZubC1yWpE8CF6nqVd5a3gjI0NwPvalnww7stipYgMRkG6E9K7VuWuXjBnNfTQ74NPnZ50bQhunqMcmTTOXkGIFueXjpmHSxsxTL6G0M3wUYyGNpS7sUKdx1KJYfTCwwDyqW76avDnAXjQHcMf4oYgS6rxi/L4ixm15WInIsZl97Ab6vqqdn5FmIWbs9C7PK7G2u9kYM7lj/42OGua9Yv7u0mvPi65goeOZ+eGiKiOyBAXsfzJqMq0TkF6r6p4E8U8CZwCJVfSCZH+Ism7nlkWp1xl9MWkPYed0+Fdt3V6S5xVni026Yt0rPq2oP+DVmn71BfQC4VFUfADM/pIrBFnjutvzw2gBxltry/abVum76XcApIrItZlh9MSaS6aBeCcwSkesw4ZSWqOoFOCpSuNcS/39aW2EeVOzfsY3q66YXd8uHh1lS1eUi8g3MispngduZGXNoEngdsD+wKfA7EVmmqve6tDdSuGPVKADd1yiA3Ve6m97YZ8sNs6Sq55HMERGR/8DsYT+olUkdzwHPicj1mLkiHdz+NEoQZ2mUwM5SNLBPk4jMSdZZzMPEKn9TKsvlwBkiMglsDLwB+I6rvTGGe9QBHqY4fuj1KmsArvz3sI4Jnsjvlhfpp8kz91rg06r6hIgcDaCq30u67lcBdwAvYqZy3+VqbIThHld4h2kcoc5T/SPuqvrWjGvfS52fBpzmw17L4e4AtlMH9jjKNUDiNpjoK7sA9wPvU9UnMsouApZgtlA7V1VPLd/EDuBq6sD2JbNwZOumm2Et1wCJXwJ+paoLgF8l59OU7Ej4XeBAYHfgcBHZ3a5ZPdozASRmdWCPswrhTsIiPZ66fAjww+T4h8C/ZhTdB1ihqvep6t+BS5JynWpRB/a4y/WZe66qPgygqg8PmQO7I/DgwPlKzNB+p6DqoA4li0ksUSnkgFpWFLTMWGuQCpDI9mFaNBbyNQrc3STaLle4V4vIDonX3oHsh+OVwM4D5zsBq4ZVOD1A4u5DbwKd6lJ3k2i7XOG+AvgIcGry7+UZeW4CFojIfOAh4DDMqpdOY6XuJtGUXAMkngr8t4h8DHgAeG+Sd32ARFXticgxwNWYV2Hnq+rdYT5Gp9GXz0knbjeKHpM8VmUnwJrlGiARzMqVdN5pm5ep6pXAlc6t69QpiFq5Hry0Ip2hVrTkczz+czp1qqJI4S6Sr+ev7ibRyV5VX4WJyK5M31fv5cCJg7HUROSlwI+AeRg+v6mqP3Cx11K4fcnnIE13o+iUL1X9I7AXrJ/B+RAz9+b+NPAHVX23iGwH/FFELkwmgpXSmMPtU11volMp7Q/8WVX/mrquwJZiturZAjM7NB2xxUod3NGp6CbRwd+UqoRZytBhwMUZ18/AvGpehYmj9n5VfdGhuR3c7VIHduTKDbPUl4hsDByM2XMvrXcCtwFvB16B2Y/vBlV9umxjWhzaeJw0lw7sEJpDQ/tzHwjcqqpZ3bSjMOGNVVVXAH8BXuVipPPc0auD2l1+we0xWTXMUl+Hk90lBzMpbH/gBhGZC+wK3OdipIM7anVgF6sRz+ssEdkMeAfwyYFr6+OoAScDS0XkTsziq+NV9VEXWx3c0aoD26hd8BZJVZ+H6XNYB+OoJbM8D/BhK1K4nUb+R0jjBHZ74DWj5e0JsxQp3JD9SmgcfvSj9hnbA++oKWK4s5QGftRAGIXP08Eci5xfhYnIriJy28Df0yLyuVSehSLy1ECeEyu3eJpi3N3TVW0Gew4NvlYqqbkZf6MpZ89tOU8W4AZVPcjVjp36gLf1P6pt7W4DxOD7ex3XGGrD5snWrNW0D5S2tDdmoNvyHdYrX3APmycL8CYRuR0zV/a4YdFYpgdIrBLtok1ePOY2xghzzN9XfKoMd8E82VuBl6nqsyKyGPgZsCCrnukBEud7CJAYuxePrW2xwRzb9zOe3fKh82QHJ7ur6pUicqaIzHadcVNesXrxWNoTA9ANfRcya8PxiMba9QH30HmyIrI9sFpVVUT2wYzOP+bBZknF5MWbbEfTMEcA8hipEtwW82QPBT4lIj3gBeAwVW3oPtm0F2/KblNAN/B5A0PcWzfJk09POZcXkSngXGAPTH/ho6r6u4F0wWycuRh4HjhSVW91tVcJbot5smdgFp9HpCa8eBNgNwF1zZ+zfR55CXCVqh6ajFVtlko/EDMmtQCz9dZZVNiCq2Uz1NqousGuE+oOZluJyFbAPwNHAiQx0dJx0Q4BLkh6t8tEZKq/s4+LzQ7uoKrzx18X1DV9ppAgp9+0Wg7vam+CFx6dysuSF2bp5cAjwA9EZE/gFuBYVX1uIH/W5pk7Ah3ccakusOuAuobPEgrmejcIyQuzNAm8FviMqt4oIksw+9p/dSBPqc0zi9TBHUShYRgBoEPAHPdOPyuBlap6Y3L+Ewzc6TzWm2cWKVK41zbdgAoKCUVIqAO2ezS8sgkz8EiWcy2Wqv6fiDwoIrsm6zL2B/6QynYFcIyIXIIZSHvK9XkbooUbzK7ATb+XLatQgIT6HgK1tw1eebuB45qmVAGfAS5MRsrvA45KvTq+EvMabAXmVdhRVYxFDDe0C/AQoIT47AHa6RvmkCA3KFW9DUg/kw++OlbMjiNeFDnc0A7AfQPj+/N6bl8HcyvUArhjl09wIobaJ9BthblHI5OnXTWGcPuEMUawRxxoHyDP9lBHC9QSuGPsmvuCKDKofQEdC8xjAnKWWgK3L/kC0kc9EUHtA+hxgLlHnSPrlVV1Vdj9wDPAOqCXnp3je5VLHKoKkw+oRwToJkEeA4/uw3PvlxN8wesql+bVNNgV7VcFukmYO5BLK3S33Osql2qqCmaV8lWgbhDoDubpGqduOWZS+zUiosDZGRuNe1zl0tSgWguhbgLoNoIc91z0yqoK976qukpE5mA2Cb9HVa8fSLde5TI9+ulUxWal5QpKE2A72hwXoF3LjTjIWaoaiWVV8u8aEbkM2AcYhNt6lcv06Kc7ZdwAaobFuVxLoK4b6FGA2UO3PNnA42bgofRmHSLyQeD45PRZ4FOqerurLWe4RWRzYCNVfSY5PgD4eiqb11Uu9alOsGuEuk6gXaCsA+Tmp6YeCywHtspI+wvwNlV9QkQOxDi7RsIszQUuM2+7mAQuUtWrQq5ycWtiHWVqgrouoGOGuV0gT5OI7AS8CzgF+Hw6XVV/O3C6DNPTdVaVvcLuA/bMuB5slUs5xQp2DVDXAXTbYXZpSw8TKCmn1pwwSwCnA18EtrSw9jHglyVbOE1jNkMtT2WhCwz1KABdNn9MILtpaJglETkIWKOqt4jIwrxKRGQ/DNxvqdKYlsAd2juGBDsiqGOCOQaQ653csi9wcLKt1ibAViLyI1X90GAmEXk1Jrb5gapaaQ1aS+AOqTLwBYI6FqBDgtHGNnuUqp5Asp9e4rmPywB7HnAp8GFVvbeqzRGEuwysIcAOBHWb4GhTW8vUGWA9d2oA+kTMt3dmMlA9Y71GGY0Y3CHADuCtQ0AdApIQ3ewm2xnJHHNVvQ64LjkeHID+OPBxX3ZaAHeTccQ8e2tbqEN4Pt8AxA5zUZ4xmLHWArhtZQusT7AbgjpmoJtom00eHzCP2cKRSOQTbI/e2ifUPqGJFWZfwBe1e7sR3ZA7pcjh9tkl9wV2hFDX6d3qhjkwyJvOfpIXLJrRRkUOt41soK0JbBuo2wZ0nW2pAeRK6rrldcoH2JFA7QOiWGBuAchTW21I7zx3K1UV7AigrgpSDDBH4JEHYR4XRQx3RfCaBrtJqJu+oUTklWekMTPNeg1y8cKRqFRlPffOwAXA9sCLmBUwS1J5FgKXY9apAlyqquk13wEUEOyQUIeEJiTQAWH2DXLe9VFTFc/dA76gqreKyJbALSJyraqmtyW9IR1xorpChT+q4K1DQd0U0MHa5O6VfYI8xRO57RgFVVnP/TBJjyaJxrIcE/wwDbeDqoTgzQPUEezYoHYFOpR3doS5aZC3LuvBK8wtt+npDuR9PSZYw/tV9SduFj09c4vILsBrgBszkt8kIrdjYqcdp6p3D6ljIEBi3q8lD95IwG4L0C2AOTTINXbRrXq6SYy1bwBXVzVYGW4R2QL4KfA5VX06lXwr8DJVfTZZx/ozzAYFMzQ9QOL8Ib8Q32BHDrVvoH3D3IFsrRI93c9geHp9VZuV4BaRWUlDLlTVS9Ppg7Cr6pUicqaIzM7ZocRRNYHtE2onmBzq8mo/G+Y2gRzDYNqwnq6I7Ai8B3g7TcKd7AN2HrBcVb89JM/2wGpVVRHZB9gI56cWn/HNhtTVJNS+gG4Q5raAvK3zomwFXZuXoSiGWlFP93TgeFVdl6znrqQqnntf4MPAnSJyW3Lty8A8WL9O9VDgUyLSw0wEOiwJmuhRw6APCLYv4JoCOiKYmwI5kAcfGkMNinu6wN7AJQnYs4HFItJT1Z+5NKbKaPlvyN5RZDDPGcAZrjY2aBjAHsD25a2zQAoJdUCgY4fZB8h1vwqz6emq6vyB/EuBn7uCDVHPUHNVILDrhroMvGXslfDOdcNcpl5bkLPsl34Ftl5rgdWOZa16ul7VYrizvHZFsEPBF+JmUSPQVWGu6pV9gtzUgJpNTzeV/8iqNlsAdxbEFWOl2YJdF9RVu92WQDcJc9XutQ+Qy7RhFNQCuG2V5bVrALsOqBsCug6YXUC2sWdjf+rxsos9e5itpNuhyOG29dqOYIeGOjTQDcA8KiDL4zOKjJwih9tGFmC7euuQUNcMtBUAjl7ZBmRfNw4rsG1A9hx/PEa1DO60N/YEtguMvm4ODkDXCXPrQc6C2NlrVxotr10Rw+1hl04XsG2A9FFHYZ3FHtoH0DajyelrPmAO1rVOw5wH+4h3zSOGO60ir21xMygLZVmobeooCXUR0GDj8aqB5uNZuaxXzhrsKvTKZUDuuuUtVtHgmW+wa4a6bqDrhtkbyHnlSgPeo+uWV1a6WSW9dhmwm4S6QaCbhLlxkCsB3h5FCneeKoBdxVtXgbqCl64CdBnv7PNGUcorNwXyiD9vQ/X13IuAJcAEcK6qnppKlyR9MfA8cKSq3lrOSoUtcV3hrOLpLaEu46XzwCvjQfO8sy+YS3llF5ArllmbuimseYraVMRLkmchZunnLMwqs7e52quynnsC+C7wDmAlcJOIXJEKG3MgJvLKAuANwFnJv47K8dquYLt6a+s6hnvpEEDbemdXj5/nldsEstuT81pcZ6jZ8CIiU8CZwCJVfUBEKu2nVcVz7wOsUNX7koZdAhzC9LAxhwAXJGu4l4nIlIjskIScsZBlWCVbsG29rGeobbvdtkDbemcf9Vl7ZZ8gO3S5bUEedlyDbHj5ACb89wMAqlpprmsVuHcEHhw4X8lMr5yVZ0dKxIHfIMubmE+w80AeWr4a1C4euirQTp45IpjTXWsbmBuYIW7DyyuBWSJyHbAlsERVL3A1WAXurOVr6ZkXNnlMxmnRT/kbLL7LuiWDNT465Nhes51LppR+Uzt4PuTu5s22g5q03bT9Xe2yPXQ1HD9svR7AJjlhlmxYmAReB+wPbAr8TkSWqeq9du2bWZmrVgI7D5zvhAlfXDYPkI5+KjfnhasJqc52M2r6s9vkU9VFFczY8vKoqj4HPCci1wN7Ak5wb+RSKNFNwAIRmS8iGwOHAVek8lwBHCFGbwSesn/e7tRppGTDy+XAW0VkUkQ2w3Tbl7sadPbcqtoTkWMwwdMngPNV9W4ROTpJ/x5wJeY12ArMq7CjXO116tRm2fCiqstF5CrgDsyuJOeqqv3jaUriPRipB4nIJ9IhYTvbo2u7aftNf/ZQihLuTp06VVeVZ+5OnTpFrMbgFpFFIvJHEVkhIl/KSBcR+c8k/Q4ReW2Ntj+Y2LxDRH4rInv6sm1jfyDf60VknYgcWqdtEVkoIreJyN0i8uu6bIvIS0Xkf0Tk9sS2tzEaETlfRNaISOYzbMjfW2NS1dr/MAMKfwZeDmwM3A7snsqzGPgl5v3gG4Eba7T9ZmDr5PhAX7Zt7Q/k+1/MoOShNX72KcysqXnJ+ZwabX8Z+EZyvB1m2srGnuz/M/Ba4K4h6UF+b03+NeW510/FU9W/A/2peINaP3VVVZcBUyKyQx22VfW3qtqf2rUM807Sl2w+O2zY7dHnZCob216nQJa0rcCWyYKjLTBw93wYV9XryV8LFur31piagnvYtNSyeULZHtTHMHd0Xyq0P7Dbo+9dKGw++yuBrUXkOhG5RUSOqNH2GcBumMkddwLHquqLnuwXKdTvrTE1tZ7b69TVALZNRpH9MHC/xYPdMvZPx+NujyVte50CWdL2O4HbMFvYvgK4VkRu0Jm7YYZQqN9bY2oKbq9TVwPYRkReDZwLHKiqPuN12NjfG4+7PZa07XUKZEnbRwGnqnkIXiEifwFeBfy+om1f7WuXmnjQx9xU7gPms2Fw5Z9Sed7F9AGO39doex5mVt2bm/jsqfxL8TegZvPZdwN+leTdDLgL2KMm22cBJyXHc4GHgNkev/tdGD6gFuT31uRfI55bG5y6amn7RMxCzjMT79lTT4saLO0HkY1t9TwFsoxt4GRgqYjciYHseFX1slJMRC4GFgKzRWQl8DVMtJOgv7cm1c1Q69RpRNXNUOvUaUTVwd2p04iqg7tTpxFVB3enTiOqDu5OnUZUHdydOo2oOrg7dRpRdXB36jSi+n/3PyUepJJXigAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 252x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1) = plt.subplots(1, figsize=(3.5, 3))\n",
    "im = ax1.contourf(model._biased_conf_energies, levels=50, cmap='jet')\n",
    "plt.colorbar(im);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fbb607c",
   "metadata": {},
   "source": [
    "## Recovering the PMF\n",
    "\n",
    "Now we want to recover the potential function (PMF). To do this, we define a number of bins (in this case 25) that we want to calculate the unbiased PMF over. We discretize all samples into these 25 bins, and use the model to calculate the PMF over the bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "eccbf4e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+wUlEQVR4nO3dd3gU5fbA8e+76Q2SAAkQSgiErlRpooKotHBV9IINvBaEH1jQq4B6BctVEbsCIoJ6KYogUgxNVJAiIKFKkQChhZJOSG/7/v5IEcgmJJDd2d2cz/PkyWZnMnNmNzmZvHPmvEprjRBCCMdnMjoAIYQQVUMSuhBCOAlJ6EII4SQkoQshhJOQhC6EEE7C1agd165dW4eGhhq1eyGEcEg7duxI1FrXsbTMsIQeGhpKVFSUUbsXQgiHpJQ6UdayKw65KKUaKqXWKaUOKqX2K6WesbBOL6VUqlJqd9HHxGsNWgghROVU5Aw9H/i31nqnUsoP2KGUWqu1PnDZehu11hFVH6IQQoiKuOIZutb6rNZ6Z9HjNOAgEGLtwIQQQlROpapclFKhQAdgm4XF3ZVSe5RSq5RSbcr4/ieUUlFKqaiEhITKRyuEEKJMFU7oSilfYDEwVmt94bLFO4HGWut2wKfAUkvb0FrP1Fp31lp3rlPH4kVaIYQQV6lCCV0p5UZhMp+vtf7h8uVa6wta6/SixysBN6VU7SqNVAghRLkqUuWigNnAQa31B2WsU7doPZRSXYq2m1SVgQohhChfRapcbgSGAX8qpXYXPfcS0AhAaz0DuBf4P6VUPpAF3KelL68QQthURapcNmmtldb6eq11+6KPlVrrGUXJHK31VK11G611O611N63179YKeNu2bfTs2ZMDBy6vmhRCCPtWUFBA9+7dmT9/vlW273C9XAoKCti8eTOnTp0yOhQhhKiU5ORktm7dSkpKilW273AJPSgoCID4+HiDIxFCiMopzlvFeayqOWxCj4uLMzgSIYSoHEnol/Hz88PDw0PO0IUQDkcS+mWUUvTu3dtqL4gQQlhLu3bteP/992nYsKFVtq+Mqi7s3Lmzlva5QghROUqpHVrrzpaWOdwZuhBCOKojR44QExNjte07ZEJ//fXX6dGjh9FhCCFEpTz//PPcddddVtu+Qyb09PR0du7cidyMKoRwJPHx8Va9/ueQCT04OJicnBwuXLi86aMQQtgvSegWyM1FQghHFB8fT3BwsNW2LwldCCFsICsri7S0NKueoVek26LdCQ0NJSIiAi8vL6NDEUKICjGZTCxYsIC2bdtabR9Shy6EEA5E6tCFEMJgZ86cYd26dWRmZlptHw6b0Nu0acNzzz1ndBhCCFEha9as4dZbb7XqtT+HTej5+fmcPn3a6DCEEKJCihN5nTp1rLYPh03owcHBUuUihHAY8fHx+Pj44OPjY7V9OGxCDwoKkoQuhHAY1r6pCCShCyGETdgioTtkHTrALbfcAoDWGqWUwdEIIUT53nnnHbKysqy6D6lDF0IIB+K0degFBQUUFBQYHYYQQpRLa83cuXOJjo626n4cNqFv2bIFNzc31q1bZ3QoQghRrpSUFIYPH86qVausuh+HTeiBgYForeXCqBDC7ll7cuhiDpvQpeOiEMJRSEK/An9/f1xdXSWhCyHsniT0K1BKSS26EMIhxMXFAdZP6A5bhw7w9NNPExoaanQYQghRrvvvv58bbriB2rVrW3U/UocuhBAOxGnr0HNzczl79qzRYQghRLl++OEHlixZYvX9OHRCf+mll2jWrJnRYQghRLk+/PBDPvnkE6vv54oJXSnVUCm1Til1UCm1Xyn1jIV1lFLqE6XUEaXUXqVUR+uEe6mgoCAyMzPJyMiwxe6EEOKq2KIxF1TsDD0f+LfWuhXQDRijlGp92Tr9gfCijyeAz6o0yjJILboQwhHYTULXWp/VWu8sepwGHARCLlvtTmCOLrQV8FdK1avyaC9T/AIVlwQJIYS9yc3N5fz58/aR0C+mlAoFOgDbLlsUApy66OtYSid9lFJPKKWilFJRCQkJlQy1NDlDF0LYu+JcZ4uEXuE6dKWUL7AYGKu1vnD5YgvfUqoeUms9E5gJhWWLlYjToiZNmjB58mRatGhxrZsSQgirqFu3LrGxsfj6+lp9XxVK6EopNwqT+Xyt9Q8WVokFGl70dQPgzLWHV75atWoxfvx4a+9GCCGumouLCyEhpQYsrKIiVS4KmA0c1Fp/UMZqy4HhRdUu3YBUrbVNCsRPnjzJyZMnbbErIYSotD/++IM33niDCxcuH9ioehUZQ78RGAbcqpTaXfQxQCk1Sik1qmidlUAMcAT4AhhtnXBLu/XWW5kwYYKtdieEEJWyYcMGJk6caJN9XXHIRWu9Cctj5Bevo4ExVRVUZQQHB8tFUSGE3YqPj8fDwwM/Pz+r78uh7xQFpOOiEMKuxcfHExwcbJPJ7J0ioUsduhDCXtnqpiJwkoSemJgok0ULIexSQkKCzRK6Q/dDBxg8eDAtWrTAbDbj4uJidDhCCHGJrVu3kpWVZZN9OXxC79ChAx06dDA6DCGEsMjFxcUmNxWBEwy5ZGZmsnnzZqqilYAQQlSl9PR0Ro4cyebNm22yP4dP6MeOHaNnz578+uuvRocihBCXOHv2LDNnzuTYsWM22Z/DJ3Rp0CWEsFfFeUmqXCqoVq1amEwmSehCCLsjCb2STCYTderUkYQuhLA7ktCvgtwtKoSwR5mZmbi7u1O7dm2b7E8VtmGxvc6dO+uoqKgq2db69eupUaMGHTvaZCpTIYSoMK11ld72r5TaobXubGmZw9ehA/Tq1cvoEIQQwiJb9HAp5hRDLocPH2bRokVGhyGEEJd46aWXmDJlis325xQJ/fvvv2fIkCE2u71WCCEqYunSpWzfvt1m+3OKhF58BVnuFhVC2BNbdloEJ0voUukihLAX+fn5JCUlSUKvLEnoQgh7UzxiIAm9kiShCyHsTXp6Ok2aNCEkJMRm+3SKssWQkBA2bdpEy5YtjQ5FCCEACA8PJyYmxqb7dIqE7u7uzo033mh0GEIIYSinGHIBWLZsGatWrTI6DCGEAOCbb76hX79+ZGZm2myfTnGGDvD2229Ts2ZN+vfvb3QoQgjBn3/+yS+//IKXl5fN9uk0Z+jSoEsIYU+Ka9Dl1v+rIAldCGFPbH1TEThhQjebzUaHIoQQktCvRVBQEPn5+Zw/f97oUIQQgsaNG9O+fXub7tMp+qEDJCUlkZaWRsOGDXFxcamy7QohhD1x+n7oUDi3aK1atYwOQwghDOM0Qy4pKSm888477Nmzx+hQhBDV3IkTJ2jbti2rV6+26X6dJqHn5OQwYcIENm/ebHQoQohq7uzZs+zfv9/mRRpXTOhKqS+VUvFKqX1lLO+llEpVSu0u+phY9WFeWe3atVFKSemiEMJwxXnI1lUuFRlD/xqYCswpZ52NWuuIKonoKrm6ulKrVi3i4uKMDEMIIUoSenBwsE33e8UzdK31BiDZBrFcM7m5SAhhD4rzUJ06dWy636oaQ++ulNqjlFqllGpT1kpKqSeUUlFKqShrTBcnCV0IYQ9CQkIYMGAAnp6eNt1vherQlVKhQKTWuq2FZTUAs9Y6XSk1APhYax1+pW1WdR06FFa6+Pj44O7uXqXbFUIIe1FeHfo1n6FrrS9ordOLHq8E3JRSta91u1cjICBAkrkQwq59t/0kf8amWmXb15zQlVJ1VVE7MaVUl6JtJl3rdq/Gxo0beeqpp8jNzTVi90IIAUDPnj0ZPXp0qefzCsy8tGQfq/eftcp+K1K2+C2wBWihlIpVSj2mlBqllBpVtMq9wD6l1B7gE+A+bVA/gQMHDjB16lSsMT4vhBAVdfjwYQoKCko9fzI5kwKzJqy2r1X2e8WyRa31/VdYPpXCskbDFZcIxcfH23RiViGEKFZQUEBiYqLFGvSYhAwAwur4WGXfTnOnKPxdxC+VLkIIoyQnJ2M2m8tI6OkAVjtDl4QuhBBVqLybio4lZlDLx52a3m5W2bfTJXSTyURqqnWuIAshxJV4enoybNgwWrRoUWpZTEKG1YZbwIna5wL4+fmRm5sr/dCFEIZp2rQpc+ZY7pQSk5hOn5bWawfgVGfoSilJ5kIIQ1mqbgFIzcojMT2XJlY8Q3eqhA4wefJk/vvf/xodhhCimpo0aRI1a9Ys1Tr37wuiktArbOPGjSxdutToMIQQ1VR8fDze3t6YTJem12OJxSWL1qlwASdM6NKgSwhhpPj4+DJr0F1MikaB3lbbt9MmdKMmvxZCVG9lJvTEdBoFeuPuar2065QJPScnh7S0NKNDEUJUQ3FxcWWeoTex4vg5OFnZIkD9+vVp0KABqamp1KhRw+hwhBDVzPDhw2nevPklz5nNmmOJGfRsZt1GtE6X0O+//37uv7/c9jNCCGE1kyZNKvXcmdQscvLNVr0gCk445CKEEEbJy8vj/Pnzpa7hWbspVzGnS+ipqakMGjSIJUuWGB2KEKKa2b17NwEBAURGRl7yvC1q0MEJE7qXlxeRkZHs27fP6FCEENVMccn05RdFjyVm4OvhSh0/D6vu3+kSuru7OwEBAVKLLoSwubISekxiYVOuosndrMbpEjrIzUVCCGOUmdATMqw+3AJOnNDj4uKMDkMIUc0U3/bv4/N38s7KLeD0+SyrV7iAA5Yt/hadwNsrDzLn0S4E1fC0uM71119PYmKijSMTQlR3/fr1o379+pc8V9zDxdo3FYEDJnQ3F8Vf59I4FJdWZkKfOtUupjgVQlQzt99+O7fffvslz/3dlEuGXEppEewHQHRcusGRCCHEpY4ePUpKSsolzxWXLNriDN3hEnotXw9q+bhzOK7sXi0//vgjHTp0kHF0IYRN3XTTTbzwwguXPBeTmEH9mp54u1t/QMThEjpAeLAvh8pJ6Lm5uezevVsSuhDCZsxmMwkJCRYqXNKtOkvRxRwyobcI9uNwXHqZLXKLX1ApXRSVlZKRy9ebj3HgzAWjQxEO5vz58+Tn51+S0LXWhTXota1f4QIOeFEUIDzYj/ScfM6kZhPi71VqefELKmfooqLMZs33O2N5e+VBUjLzALijdTBP9wmnbUhNg6MTjqD4BDI4+O9JoBPTc0nLzrfJBVFw0ITevOTCaFq5CV3O0EVFHDqXxn+W/sn24yl0bhzAhP4t2XwkidmbYvjpQBx9WgbxVJ9w2jf0NzpUYccs3VRU0sPFBjXo4LAJvfDFORyXRu8WpRvJ+/v7c/PNN1OnTh1bhyYcSEZOPp/8cphZm45Rw9OVKfdcz72dGmAyKTqHBvJoz1DmbDnBFxtjuGvaZm5pXoen+4TTqXGA0aELO9S0aVNmzJhBmzZtSp6LKS5ZtEGFCzhoQvf3difIz4ND5yyXLiql+O2332wclXAUWmvW7I/j9R/3cyY1m6GdGzKhf0sCfNwvWc/P040xvZvxcI9Q5m09wRcbYrjns9/p2aw2T/cJp0uTQIOOQNijkJAQRo4ceclzxxIzcHc1Ud/CSII1OGRCh8Jhl8PxMs2cqJxTyZlMWr6fX/+Kp2VdPz65vwOdQ8tPzL4eroy6pSnDuzdm/taTfL4hhiGfb6FbWCBP9wmne1gtqzddEvbv6NGjpKWl0b59+5LnYhLSaVLLBxeTbX4+HLLKBYoSelw6ZrPlSpcnnniC/v372zgqYa9y881MW3eE2z/8ja0xSbw8oBU/PtXzisn8Yt7uroy4OYyN43ozMaI1MQkZPPDFNoZ8voUtR5OsGL1wBO+9916pu0RjEjJsdkEUHDqh+5KVV0BsSpbF5VlZWfz11182jkrYo9+PJtL/4w28u+YQvZoH8fNztzDi5jDcXK7ux9/L3YVHezZhw7jevH5nG2JTsnhg1lb+OJZcxZELRxIfH3/JBdG8AjMnkzMloVdE+EWVLpYEBwdLlYvg1eX7eeCLbeQWmPnqXzcwY1inKhvP9HRzYXj3UNY+dwsNA7x59rvdpGXnVcm2heO5PKGfTM4k36xpYqMadKhAQldKfamUildKWZwCSBX6RCl1RCm1VynVserDLK240iW6jHH0oKAgMjMzycjIsEU4wg79uOcMX/9+nAe7NuKnsbfQu2Xpiqiq4OvhyodD23M2NYvXfjxglX0I+xcfH39JDfoxG80jerGKnKF/DfQrZ3l/ILzo4wngs2sP68r8PN2oX9OT6HNlJ3SQWvTqKj4tm1eW7aNdg5q89o82eLm7WHV/nRoHMKZ3M77fEcvqfWetui9hny4/Q49JLKzCa2pPZ+ha6w1AeYODdwJzdKGtgL9Sql5VBVie8GC/MrsutmrViqFDh2IyOeyokrhKWmteXPwnWbkFvD+kPa4VHCtPTk7myJEjJV+fO3euUvt9uk841zeoyYs//En8hexKfa9wfHPnzuXRRx8t+TomIYNaPu7U9HazWQxVke1CgFMXfR1b9FwpSqknlFJRSqmohISEa95xi7p+HElIp8BCpUvXrl1ZsGABjRs3vub9CMfy/Y5Yfvkrnhf6tqBZUPlnR2azmV9//ZUHHniA+vXr8/TTTwMQHR1NSEgId9xxBwsXLiQnJ+eK+3VzMfHh0PZk5RXwwvd7y+w1JJxTRETEZSWLtq1wgapJ6JYKLC3+JGutZ2qtO2utO1fFXZzhQb7k5ps5kVT2OLn8UlUvp89n8fqPB+jSJJBHb2xS7rqzZs2iWbNm9OnTh1WrVjFixAjeeustoPBu40mTJnHo0CGGDh1KSEgIzz333BWH8JrW8eWlAa34LTqBuVtPVNlxCfuWmJhIZGQkycl/D2bEJGbYpAf6xaoioccCDS/6ugFwpgq2e0XNy5nsIjc3F39/f95++21bhCLsgNmsGf/9Xgq05r1722G67GaOvLw8li1bRnZ24XBIUlISoaGhzJ8/nzNnzvDpp5+WnGEFBQUxceJEYmJiWLNmDbfeeiuzZ8/Gza3w3+dDhw6Rnm55uG9Yt8bc0rwOb644yJF4mYilOoiKimLQoEElpdIXsvNITM+xWQ+XYlWR0JcDw4uqXboBqVprm1wVCr+op8vl3N3dMZvNclG0Gpm/7QSbjiTy8sBWNKrlXfL84cOHmTBhAg0bNuSuu+5i+fLlAIwbN65kuMXLy3Ipo4uLS8mwy9mzZwkIKOzj8vDDD1OvXj1GjBjB1q1bL/lPUCnFu/dej7e7C89+t5vcfLMVj1rYg8s7LcYk2LaHS7GKlC1+C2wBWiilYpVSjymlRimlRhWtshKIAY4AXwCjrRbtZbzdXWkY6FXmZBdBQUGS0KuJ44kZvLXyL24Kr80DXRoBkJOTw5AhQ2jevDnvvfceXbt2Zfny5QwePBig0rfre3v//Ufigw8+YMiQIXzzzTd0796dLl26XHIRNaiGJ28Pvo4/T6fy6a+Hq+AIhT27vNOirbssFqtIlcv9Wut6Wms3rXUDrfVsrfUMrfWMouVaaz1Ga91Ua32d1jrK+mH/rXlQYQsAS+TmouqhwKx5ftEeXF0UU+69viRR5+fnk5KSwn/+8x9OnTrFsmXLGDRoEK6u197CqEePHsyePZuzZ88yY8YMDh48yIoVKy5Zp1/betzbqQHT1h1hxwm5i9SZxcfH4+npia9vYQKPScjAxaRoFOh9he+sWg7bnKtY87p+bDicQF6BudSt3EFBQRw9etSgyIStzN4UQ9SJFD4Y0o56Nb0wm81kZ2fj4+PDmjVrrFq6WqNGDUaOHMnAgQNp0KBBqeWTBrVma0wSz363h5XP3ISvh8P/ygkL4uLiCAoKKjmZOJaYQcMAL9xdbVs27fBF2s2Dfckr0BxPLF3pMmjQIO655x4DohK2cjgujfd+iuaO1sHc3SEErTVjx47l1ltvJTMz02b3IRQn861bt3LvvfeWXHj183Tjw6HtOZWSyRtyF6nTeuWVV1iwYEHJ10cT0m0+3AJOkNDDg8qudHn00UeZNGmSrUMSNpJXYObfi/bg4+7Cm3dfh1KKiRMn8umnn9KjR48yL3Ra05EjR1i8eDFDhw4lL6+wr8sNoYGMuqUp30WdYs3+yt2sJBxDs2bN6N69O1BYbXU8KcPmF0TBCRJ6syBfTKrsJl0ZGRkUFBTYOCphC5+tP8re2FTevPs66vh58N577/Hf//6Xxx57jPfff9+QHuUPPfQQU6dOZfny5TzyyCOYzYUVLs/e1pw29WsU3kWaJneROpu5c+eyY8cOAM6kZpGdZ5Yz9Kvh6eZC41o+FhP6okWL8PX1lTa6Tmjf6VQ++eUw/2hXnwHX1WPOnDm88MILDBkyhM8//9zQCSfGjBnDm2++yfz58xkzZgxaa9xdTXw0tD0ZOfmMl7tInYrWmscff5yFCxcChePngM1vKgInSOhQeMeopYRefDeqVLo4l5z8Av69cA8BPu68fmfh/I09e/Zk5MiRzJ07FxcX6zbiqogXX3yRcePGkZCQQH5+PlDYe2hC/5asO5TAN3+cNDhCUVUuXLhAbm7uRSWLhQm9qY1v+wcnSegt6vpxPCmTnPxLh1ak46Jz+ujnwxyKS+Ode67j9LHDaK0JCwtjxowZuLu7X3kDNqCUYvLkyXz33Xe4ubmV9IJ5uHsoN4XX5r+RB0tqlYVjs1SD7uvhSh0/D5vH4hQJPTzYjwKzLvnLWKz4BY6LizMiLGEFO06k8PlvRxnauSEucX/RuXNnJk+ebHRYFimlcHFxISEhgU6dOjF9+nRMJsW797bD3dXEsNl/EHVc6tMdXamEnljYlMuIYT+nSOglk11cNuwSGBiIyWSSM3QnkZVbwPOL9lCvphf962YQERFBWFgYI0aMMDq0cvn7+xMWFsaYMWOYN28edWt6MufRLriYFEM+38IHa6PJL5D2AI7K0m3/Royfg5Mk9LDavriaVKk7Rk0mE6+++iq33HKLQZGJqvTB2kMcS8zg/6534547I6hTpw5r166ldu3aRodWLjc3NxYuXEjv3r3517/+xdKlS2nX0J8VT/fk7g4N+OSXw/zz8y3ldg0V9uuOO+5g3759tGzZkuy8As6kZhFmw0ktLuYUCd3d1URobR+LPV1eeeWVUjNxC8dz8OwFvtx8nCEd6vLKmOF4enry888/U79+faNDqxBPT0+WLVtG586dGTp0KOvXr8fP0433h7Tj0/s7cCQ+nQEfb2TxjlipgHEwPj4+tGnTBk9PT44lZqC1baedu5hTJHQoHHax1HXxwoULnDp1ysJ3CEdhNmteXvInNb3ceGnQdcycOZO1a9cSFhZmdGiV4ufnx8qVK7n77rtp0aJFyfOD2tVn9dibaRNSk38v2sNT3+4iNVMmm3YUkZGRzJo1C7ioy6Ik9GvTPNiPE8mZZOVeWukyatQoevfubVBUoiosjDrFjuNJDK6fjr+3O7fffjutW7c2OqyrEhgYyIIFC6hXrx5aa7KysgAI8ffi2xHdeKFvC1bvO0f/jzewLSbJ4GhFRcyZM4f33nsP+LvLooyhX6PmwX5oXdhD4WLScdGxJWfkMnn1X/gf+YlXHh/Mtm3bjA6pSmitufPOOy+Zg9LFpAonmv6/Hri7mrjvi628u+Yv8uSCqV27eHLoY4kZ1Kvpibe7MU3YnCihW650CQoKIi0treRMSDiWt1ceJOnkYQ6u+ILBgwfTpUsXo0OqEkqpknlvL27qBNC+oT8rnr6Jf3ZqwLR1R7n3s99L7j4U9ic+Pr6kwuVoou3nEb2Y0yT0xrV8cHcxlbowKjcXOa4/jiWzcNsx8n/9lMCAAGbMmGHoLf1Vbfz48XTr1o3Ro0dz+vTpS5b5eLgy5d52TH+wI8eTMhn4yUYWbj8lF0ztUPEZutaamIR0wypcwIkSupuLibA6PqVKFyWhO6a8AjP/WfonBVELOXfsELNmzaIqJha3J66ursyZM4ecnBwee+wxi8l6wHX1WPXMTVzfoCbjFu9l9qZjBkQqypKfn09SUhJBQUEkpueSlp0vZ+hVJTzYr9SQS/v27fn0008JCQkxKCpxNWZvOkZ0XDpDbu3MM888Q0REhNEhWUV4eDjvvfce+/btK7Maq76/F/Mf78ZtrYJ4d80haRlgR1xcXEhJSWHs2LGGNuUq5lQJvUWwL7EpWWTk5Jc817BhQ5588kmHqVcWEJuSycc/H+a2VsF8MPE5PvroI6NDsqpRo0axf/9+GjVqVOY6LibFW3dfh4erifGL92I2y9CLPVBK4e/vT82aNUv+0DY1oG1uMadK6OHBhZNdHI7/+wxGa82+ffs4eVK62zmK1348QPzPs2idvsPoUGxCKUXNmjXJzc1l+vTpJd0ZLxdUw5OJg9qw/XgK/9ty3LZBCov279/PhAkTOH36NDGJGbi7mqjvb/uJVYo5VUJvHlw8e9Hfwy7F1QQff/yxUWGJSvj5QBzLf4wkeetizp2oXvPB/vTTT4wZM6bcZmP3dAyhV4s6TFl9SFoF2IHdu3fzzjvvkJ6eTkxCOk1q+eBiMu7CvVMl9EaB3ni4mog+V7rSRS6K2r/M3Hxe+vZ3zq/5lDZt2vL6668bHZJNRUREcN999/Haa6+xc+dOi+sopXh78HW4mpQMvdiBizstxhhcsghOltBdTIpmQb5Ex5eudJGEbv8++eUw+xd9gDk7jXnz5uLp6Wl0SDY3bdo0goKCGDZsWMlE05erV9OLlwe2YmtMMvNlogxDxcfH4+rqio9fDU4mZRp6QRScLKFD4bDL5T1dgoODOXdOJue1Z4fj0pi6YBWZhzbx+muv0b59e6NDMkRgYCBffvklBw4c4D//+U+Z6w29oSE3hddm8sqDxKZk2jBCcbG4uDiCgoKITcki36wNmUf0Yk6Z0M+mZpOa9Xdzo9atW3PgwAFSU1MNjEyURWvNf5buo1ZYW5asWM24ceOMDslQffv25Z133uGhhx4qc53ioReACYv/lBuODHL69GnatGljeFOuYk6Y0Av/Qh6J//ss/dFHH2X58uV4eRl39VmU7fuoU2yK2sP4fi25a0Bfu5gT1Gjjxo0r+S+loKDA4joNAryZMKAVm44ksmC7dBQ1wurVq1m0aBExiYXDvGEy5FK1/q50+XscvXnz5vTv399u5psUfzufmctzr77D2S+fpIVrotHh2BWtNSNHjix3RqYHuzSie1gt3lxxkDPnpV+RrRWXnB5LzCDQxx1/b2NzjNMl9BB/L7zdXTh0WaVLdHQ0b7zxRplnO8IY42ev5vRPs+jZqw8dO3YwOhy7opSiTp06fPXVVyxbtsziOiaT4p17rqfArHnxBxl6saUHHniAiRMnAnA0IcPws3NwwoRuMinCg3w5HH9pQt+5cycTJ07kjz/+MCgycbk/jsbz9dvP4+nlxXfzvnaqxltVZeLEiXTo0IERI0aUWanVqJY34/u14LfoBL7fEWvjCKuntLQ0vv/++5JKpJgE40sWwQkTOhT3dLm0dLFv38Kx2RUrVhgUlbhYfoGZB598kdyzh/l8xgzq1atndEh2yd3dnblz53LhwgWeeOKJMs/Ah3cPpUtoIG9EHiDuguVyR1F11q5dS15eHhEREVzIziMxPcfwChdw0oTeItiPhLQcUjJyS54LCAigZ8+eREZGGhiZKDZ51V8kpOdw2z+GMOyB+4wOx661adOGt956i3Xr1hETE2NxHZNJ8c6915OTb+blJTL0Ym2RkZH4+/vTo0cPjhVXuDjKkItSqp9S6pBS6ohSaoKF5b2UUqlKqd1FHxOrPtSKCy9jsouBAweyZ88emWPUYIt3xDJr0zHGvvAiPy1dcOVvEIwdO5YDBw7QtGnTMtdpUtuHF/q24OeD8SzbfcaG0VUvZrOZFStW0K9fP1xdXf+ucHGEIRellAswDegPtAbuV0pZmtBxo9a6fdGHofdsl1S6XHbHaEREBN7e3uzbt8+IsASw59R5Rr7wCqFZh3l5YCsZN68gk8lESEgIWmuWLFmC2Wx5WrpHbmxCx0b+vPrjfuLTZOjFGjIzM7nvvvt44IEHgMLxcxeTolGgAyR0oAtwRGsdo7XOBRYAd1o3rGtTr6Ynfh6upXq6tGzZkqSkJPr3729QZNVbfFo297/+JUnr/0dY5kHcXJxyxM+qVq5cyeDBg5k6darF5S4mxZR725GZW8ArS/fJ0IsV+Pr68vHHHzNo0CCgMKE3DPDC3dX4n+eKRBACXDxGEVv03OW6K6X2KKVWKaXaWNqQUuoJpVSUUioqISHhKsKtGKUU4cG+pYZclFLVsj+IPcjJL+DxWRs5smgKDRuH8ulHHxgdkkMaMGAAAwcOZPz48Rw6dMjiOs2CfHn2tuas2R9H5N6zNo7Q+e3ateuS8ufCplzGXxCFiiV0S/8TX/5nfyfQWGvdDvgUWGppQ1rrmVrrzlrrztaeTqx50exFl5+hHD16lPbt20u1iw1prZm0bD/r/vc++Rfi+XbeXHx8jP/31BEppfjiiy/w9vZm+PDhZfZOH3FTE9o1qMmk5fu5kJ1ncR1ReadPn6Zjx458+OGHAJjNmmOJ6YY35SpWkYQeCzS86OsGwCVXXLTWF7TW6UWPVwJuSqnaVRblVWge7EdKZh6J6bmXPN+gQQOOHDki1S42NG/rCf639CfS96zhhRde4MYbbzQ6JIdWr149pk+fzh9//ME777xjcR1XFxP/ves6kjNymbvlhI0jdF4rV64EoF+/fgCcvZBNdp7ZLi6IQsUS+nYgXCnVRCnlDtwHLL94BaVUXVV0dUsp1aVou0lVHWxlFF8YvbzzooeHB3fccQeRkZEyvmgDW2OSeO3HA/S79WbmzJlb7XqcW8vQoUN55pln6Ny5c5nrXNegJrc0r8OXm46RlSt3SFeFyMhIGjduTJs2haPKxdPOhdV2kCEXrXU+8CSwBjgILNRa71dKjVJKjSpa7V5gn1JqD/AJcJ82OFs2L6N0EQqrXWJjY9m7d6+tw6pWYlMyGT1/J3U98/j4/g4MG/YQHh4eRoflND766CP69u1b7jpjejcjKSOX77ZL3/RrlZ2dzc8//0xERERJdVZxl8WmDnSGjtZ6pda6uda6qdb6zaLnZmitZxQ9nqq1bqO1bqe17qa1/t2aQVdEHT8P/L3dOBRXeob0AQMGAMg4uhVl5RYwcu4OEvesY897w4n5a7/RITklrTWTJk3i5Zdftri8S5NAbggNYOaGGHLzLZc6iopZv349mZmZDBw4sOS5mIR0fD1cqeNnHycqxtfZWIlSiuZBpSe7AKhbty7PP/88119/vQGROT+tNeMW72Xv4eOk/jyD1q1a0rZtW6PDckpKKc6ePcvbb7/N5s2bLa4zunczzqRms3TXaRtH51x69+7NTz/9RK9evQDIzitg1b5zdGjkbzf3UzhtQofCO0YPWah0AXj33XeJiIgwICrnN+O3GJbvPk3N7bPJzcnif//7H66urkaH5bTef/99QkNDefjhh0lPL/0faa/mdWhTvwaf/XaUApmD9Kp5eHhw++23l8yrsHhnLPFpOYy6pey7d23NqRN6i7p+pGXnE3chx+LyU6dOER0dbeOonNu6v+KZsuYvmiZvY/+235g8eTItW7Y0Oiyn5ufnx9dff01MTIzF2Z6UUozp3YxjiRms2id16VcjOjqaF198kbNnC1+//AIzM347SruG/vRoWsvg6P7m1Ak9PKh4sovSwy5aa7p27corr7xi67CcVkxCOk8v2EXLujVo6XGe3r1789RTTxkdVrVw88038+yzzzJr1ixOnChdpti3TV3C6vgwbd1Rqe66CkuWLGHy5MklLRd+3HuGU8lZPNm7md0Mt4CTJ/TyKl2UUgwYMIDVq1eTlyc3XlyrtOw8RsyJws3FxMxhnZg+9VNWrlyJyeTUP2J25c033yQqKorGjRuXWuZiUozu1YyDZy+w7pDlvuqibJGRkXTo0IGQkBDMZs30dUdpEexHn5ZBRod2Caf+bavl60FtX3eLCR0KyxcvXLjApk2bbByZc0nPyWf0/J0cT8rkTt8Ykk4WDmNJmwXb8vT0LLnQv2vXrlLL72xfnxB/L6b+ekTO0ishKSmJ33//veSa29qDcRyOT2d076aYTPZzdg5OntChcNjl8skuit122224u7vLXaPX4FhiBndP28zmI4mMbufO2y8+wxtvvGF0WNXa8uXL6dixI4sXL77keTcXEyNvCWPnyfNsO5ZsUHSOZ/Xq1ZjNZgYOHIjWmunrjtAo0JuB19nfpCxOn9CbB/tyuIxKF19fX3r16iX16Fdp3V/x/GPqJhLTc/hqeEe+fXcCvr6+TJs2zejQqrX+/fvTqVMnRo4cSVxc3CXLhnRuSG1fD6atO2JQdI4nISGBZs2accMNN7D5SBJ7YlMZdUtTXO2wW6j9RVTFmtf1IyO3gNNlzIj+4Ycfsn79etsG5eC01kxbd4RH/7edhgHeLH+yJ+u++5wdO3YwY8YM6tata3SI1Zqbmxtz5swhPT291LR1nm4uPH5TEzYeTmTPqfPGBelAxo4dS3R0NCaTianrDhNcw4N7OllqOGs850/oJT1dLA+7tG7dWhJQJWQUjZe/u+YQg66vz+L/68GJAzt5/fXXefDBB7n33nuNDlFQ+HP99ttvs3z5cl599dVLlj3YtRE1PF2Zvl7O0q+kuJulUoodJ5LZGpPMiJvC8HB1MTgyy5w/oReVLh4q48IowKJFi8q8dVr87XhiBoOn/86a/ed4eUArPr6vPV7uLnTr1o0pU6Ywc+ZMo0MUFxk7diyjRo2iVatWlzzv5+nGv3qEsmZ/nMU7qcXfXnrpJTp27EhBQQHT1x0lwNuN+7s0MjqsMjl9Qq/p7UZwDY8yK10AoqKimDJlChcuXLBhZI5l/aHC8fK4tGzmPNqVETeHERMTw5kzZ3B1deX555/H29vb6DDFRZRSfPbZZ9x3X+Ek3Glpf/8OPHJjE7zdXZi+/qhR4TmEFStWUKtWLQ7FZfDLX/E8cmMTfDzs965np0/oUDjssuvk+TKbE0VERJCfn8/atWttHJn901ozff0RHvl6O/X9vfjxyZ70DK9NQkICffv2ZdCgQVIC5wB+/PFHmjRpwu7duwEI8HHngS6NWL7nDCeTMo0Nzk7FxMRw4MABIiIi+Oy3o/h6uPJw91CjwypXtUjo93dpxLHEDF784U+Lyad79+4EBARI+eJlMnLyefKbXUxZfYiB19Xjh9E9aBjoTWZmJoMGDeL06dNMmzbNru6UE5Z17NgRLy8vBgwYwMmTha10R9wchotSfL5BztItKa5+u757b1bsPcND3RpT09vN4KjKVy0S+oDr6vFMn3AW74y1+C+mq6sr/fr1Y+XKlWXOpl7dnEzK5J7PfmfVvrO82L8ln97fAW93VwoKCnjwwQf5448/+Oabb+jWrZvRoYoKCAkJYdWqVWRkZDBgwADOnz9PcA1P7u3cgEVRscRfyDY6RLuzYsUKWrRoweqThTX8j/VsYnRIV1QtEjrA2NvC+Ue7+ry75hAr/yzdoGjQoEGEhoYSH1+9b4suMGuW7T7NoKmbOJuazdePdGHkLU1LzsKnTJnC0qVL+eijj7j77rsNjlZURtu2bVmyZAnR0dEMHjyY3NxcRt3clHyzmS82xhgdnt0ZPnw4o8e+wA+7Yhl6Q0O76XleHmXU+Gfnzp11VFSUTfeZnVfAg7O2se90KgtHdqddQ3+b7t+eFZg1kXvP8MkvhzmakEHbkBpMf6ATjWpdeqEzJSWFhQsXMnLkSIMiFddq3rx57N+/nzfffBOTycTYBbv46UAcm8ffSoCPu9Hh2ZXXftzP3C0nWP9CLxoE2MdFf6XUDq21xbkHq1VCB0hMz+GuaZvJyTezdMyNhPh7XbI8Ozu7WvUgKTBrftxzhk9+PUxMQgYt6/rxdJ9w+rWpe0mfis2bN9OpU6dq9dpUB+fPnyc+x4U7PtzAM33Cefb25kaHZBd++eUXAus24P4FMURcX5/3/tnO6JBKlJfQq82QS7Havh589a8byM4t4LGvt5Oek1+ybOHChQQEBHD6tPPP7JJfYGbJrlhu/+A3xn63G3cXE5892JGVT9/EgOvqXZLMN23aRJ8+fZgwYYKBEYuqdvToUVq0aMH6Zd9yR+tgvv79+CW/D9WV2WzmoYce4pExz5GTb7arCSyupNoldIDwYD+mP9SRw/HpPP3trpJZXNq2bUt2drZT93bJLzDzw85Ybv9wA89+twd3VxMzHipM5P0vS+QAhw4d4s4776Rx48bSO97JNG7cmM6dOzN69Giu0zGkZuUxf2vpXurVzc6dOzl37hyJAa3p37YuzYJ8jQ6pwqplQge4KbwOr/2jDb/+Fc9/VxwAoFWrVoSGhjpl+WJ+gZnvd8Ry2we/8dzCPXi6uTDjoU6sfPom+rUtncgB4uLi6N+/P66urqxatYpatexnZhZx7VxdXfnuu+9o3749L455lFZuiXyx8RjZeQVGh2aoyMhIlFKohh0Y3auZ0eFUSrVN6AAPdWvMozc24avNx5m79QRKKSIiIvj555/JyrLczMvR5BWYWRR1ij4f/Mbzi/bg7e7K58M6seKpnvRrW7fcfs4PP/ww586dIzIykrCwMBtGLWzF19eXyMhIateuzR8zxnHu9CkWRZ0yOixD/Ri5Au8GrejToRltQ2oaHU6lVOuEDvDywFb0aRnEq8v381t0AhEREWRlZTl8B8b0nHxmbYzhlinreOH7vfh6uDJzWCdWPN2Tvm3KT+TFPvroI3744QduuOEGG0QsjFK3bl1WrVpF/763c0OLhkxe9Rfbj1fPfunJycns3r0LtyadGdPbsc7OoRpWuViSkZPPvTO2EJucyTePdmT5vJk8+OCDhIaGXvF7s3IL2HA4gZ/2x7HvdCo9w2szuGMIrevVMOQOyvgL2Xz1+3HmbT1BWnY+XZsEMvKWMHq3CKpQPLGxscybN4/x48fLHaDVUPyFbO58Zynn4hNZ8OL9dA2rXsNsuflmbpy0hPoBPix7vp/R4VgkZYsVcOZ8FndN24y7q4mlY26ktm/ZNxEkpefwy1/x/LQ/jo2HE8jJN1PD05XW9Wuw40QKeQWalnX9GNwxhLvahxBUw/qlfkfi0/liQwxLdp0m32ymf9t6PHFzWIVr7c1mMzNnzmTcuHHk5+ezZ88ewsPDrRu0sEv3DLmPJYsXEdDtHr6f+T692zQ0OiSbWbj9FOMW7+XrR26gVwv7mi+0mCT0Ctobe54hn2+hRS0PHmxwni6dO5acpZ9IyuCn/XGsPRBH1IlkzBrq1/TkjjZ1uaN1MDc0CcTNxURKRi6Rf57lh52x7Dp5HpMqvAA7uGMId7Sui5d71fVR1loTdSKFz387ys8H4/F0M/HPTg15/KYmNK7lU+HtREdHM2LECDZs2ECfPn2YOXOmjJlXYykpKTw19jnmz/kat8D6fPjpZ4x54B9Gh2V1mVnZNOp0K836DGXLJ0/b7X+oktArYfW+szz+2VpOT3+Ym2/rS6/HXmHr2YKSfuqt6tXg9tbB3NE6mDb1yx9WiUlI54edp1my6zSnz2fh6+HKgOvqMrhjA7qEBl71BLMFZs3aA+f4fEMMu06eJ8DbjeHdQxnevTG1yvnPwpL8/HyaNWtGamoqH3zwAf/617/s9gdZ2NaSyFU8+K/HyUo6w8QPPue1Z58wOqQqZzZrdp06z4qdx5n10duc2rCIV6fNZdLoh4wOrUyS0Cvps/VHeenVNzi/cT4mNw/a3z2K/xs1ir7X1adhYOVv/zWbNduOJfPDzlhW/nmWjNwCGgR4MbhDCD2a1UZryDebyS/Q5BWYySvQ5JuLPheYyTMXfs4v0GTk5rNs9xmOJWbQKNCbx29qwj87Naz0mf++ffto2bIlrq6ubNq0iaZNm1Kvnv1NeiuMFZtwnl7DnkVfF8Hsx2+iYz1PatSoYfM4cvILyM41V0m3w+y8An4/msjaA3H8fDCeE7s2kPLzTPJT47hpwD38svRb3Nzst6uiJPRK0lozb9tJkk8fY8nU1/lt/Truuecevv/++2vedlZuAT8dOMfinafZdDgB81W8/Nc3qMnIm5vSr21dXCp5lp+VlcXrr7/Ou+++y/vvv88zzzxT+QBEtZKSkctDs7cRfSaF/MXj6NqpPZ988glBQdYfYz6RlMG8rSdYGBVLalYegT7uNK3jQ9M6vjSt40tY0eMGAV7lTtp8PjOXX4uue204nEBmbgE+7i60dInjh1f/RctWrfhs+nR69epl9WO6VpLQr4HWmoULFxIYGMjtt99OVlYWWVlZBAYGXvO24y5kEx2XhqvJhJuLwtXFhKtJ4eZiwtVF4WYq/HzxYzcXE55uVzcOv3HjRh5//HGio6N59NFHee+99wgICLjm4xDO73xmLg/O/J2tS2Zz4ffvqFHDj48++oiHHnqoyofozGbNb9EJzNlynPXRCZiUol+bulzfoCbHkzI4Gp/B0YR0kjJyS77H3cVE41rehYk+qDDJNwz0Zm9sKmsPnGP78RQKzJrgGh70ahZASN5pRg4diIerC4sWLeLOO+/E3d0xGpNJQq9CEydO5LPPPmPKlCk8/PDDmEyOUco/efJkXnzxRZo0acLMmTO57bbbjA5JOJjUrDyGz97G7r378PljNgf3RNGvXz/mz59fJSc45zNzWRh1inlbT3IyOZM6fh480KURD3RtRLCFSrGUjFxiEtM5mlCY4I/GZxCTkM6J5MySdh4AzYN9i6571SXu4HaeeupJTpw4wfHjxx1ygnhJ6FVo7969jB49ms2bN9OjRw+mT59Ou3b204kNoKCggIMHD7Jt2zY6depE+/bt2bJlC4sWLeKNN97Ax6fiFTBCXOxCdh7DZ//Bn6eSuUPtZv/mtfzyyy+4ubmxaNEiGjZsSIcOHfDwqPjF+X2nU5mz5TjLdp8hJ99Ml9BAhnVvTN82dXF3rfwJU26+mZPJmZxIyqBZkC+Na/lw5swZ/v3vf7NgwQKaNm3K1KlT6dfPPuvMr6S8hI7W+oofQD/gEHAEmGBhuQI+KVq+F+h4pW126tRJO6qCggL91Vdf6dq1a2uTyaSnTZtmdEg6IyNDjx8/Xvfq1Uv7+vpqQAN60qRJRocmnMyFrFx997RNOuzFFfrH3ae11lrn5eVpb29vDWg3Nzd9ww036CeffFL/+uuvFreRnZevf9h5St81bZNuPD5St/zPKv3iD3v1gTOpVR5vSkqK9vf31x4eHvrVV1/VWVlZVb4PWwKidBl59Ypn6EopFyAauB2IBbYD92utD1y0zgDgKWAA0BX4WGvdtbztOuoZ+sWSk5N5+eWXGTVqFO3atSMpKYmYmBg8PT1LPry9vausqVV6ejpZWVklZ9/btm2jRYsWvPnmm5jNZurWrUtoaChdu3alS5cudO3alfDwcClDFFUuPSefR776g50nz/PR0PYMalefU6dO8dvmLWz6fStR2/9g/95d3Dns/7jtwdGcjktkzlvP4d+4Nd4hLTjvFUJqvomwurUYfmMY93RqQE2vq68sycrKIjk5mezsbLKzs8nKyiIvL4/u3bsDMGvWLHr37k3Tpo7TCrcs1zTkopTqDryqte5b9PWLAFrrty9a53Ngvdb626KvDwG9tNal53or4gwJ/XLDhg1j3rx5lzwXGBhIUlISAEOHDmXFihV4enri5eWFp6cnYWFhrFmzBoB//vOfrFu3jry8vJKPtm3bsmvXLgC6dOnC9u3bS7YdFhbG0KFDeeutt4DCmnJXV1dbHKoQZOTk88jX24k6nkwtXw9SMnLJv2jsWpsL0Pl5mNw9MSefInH5O2TFn4CLcs533y1kyJB/snbtWgYMGICbmxuurq64ubnh5ubGggUL6NWrF2vXruXJJ58sSdbFiXvjxo107dqV2bNn8/jjj5eKcePGjfTs2dMmr4etlJfQK/LbHwJc3H4tlsKz8CutEwJcktCVUk8ATwA0atSoArt2HFprRo4cyX333XfJD9zFF00jIiJo0KBBybLs7Gz8/f1Llvfs2ZOgoKCSH2Y3Nzfq169fsvzZZ58lISGBsLAwunbtSp06dS6JQZK5sCUfD1e+fuQGPvgpmvScfAJ83An0di/87ONGgLc7gT6FX/t5uKJmjiItLY2oqCgOHDhATk4O7dsXXn8KDQ1l3Lhxl5zM5OXllZRG+vr60qFDB7y8vEpOhjw9PUsuat544418/vnnlyzz8vKiQ4cOhr0+RqjIGfo/gb5a68eLvh4GdNFaP3XROiuAt7XWm4q+/gUYp7XeUdZ2nfEMXQghrO1ap6CLBS7uztMAOHMV6wghhLCiiiT07UC4UqqJUsoduA9Yftk6y4HhqlA3ILW88XMhhBBV74qDrlrrfKXUk8AawAX4Umu9Xyk1qmj5DGAlhRUuR4BM4BHrhSyEEMKSCl1F01qvpDBpX/zcjIsea2BM1YYmhBCiMhzjvnUhhBBXJAldCCGchCR0IYRwEpLQhRDCSRjWbVEplQCcuMpvrw0kVmE4RpJjsU/OcizOchwgx1Kssda6jqUFhiX0a6GUiirrTilHI8din5zlWJzlOECOpSJkyEUIIZyEJHQhhHASjprQZxodQBWSY7FPznIsznIcIMdyRQ45hi6EEKI0Rz1DF0IIcRlJ6EII4SQcIqErpf6plNqvlDIrpcos9VFKHVdK/amU2q2UssvZMypxLP2UUoeUUkeUUhNsGWNFKaUClVJrlVKHiz4HlLGeXb4vV3qNi9pBf1K0fK9SqqMRcVZEBY6ll1Iqteg92K2UmmhEnFeilPpSKRWvlNpXxnJHek+udCxV/56UNXu0PX0ArYAWwHqgcznrHQdqGx3vtR4LhW2KjwJhgDuwB2htdOwW4pwCTCh6PAF4x1Hel4q8xhS2hF4FKKAbsM3ouK/hWHoBkUbHWoFjuRnoCOwrY7lDvCcVPJYqf08c4gxda31Qa33I6DiqQgWPpQtwRGsdo7XOBRYAd1o/ukq7E/hf0eP/AXcZF0qlVeQ1vhOYowttBfyVUvVsHWgFOMrPyxVprTcAyeWs4ijvSUWOpco5REKvBA38pJTaUTQhtaMqa9JtexOsi2amKvocVMZ69vi+VOQ1dpT3oaJxdldK7VFKrVJKtbFNaFXOUd6TiqrS98RupolXSv0M1LWw6GWt9bIKbuZGrfUZpVQQsFYp9VfRX0mbqoJjURaeM6S+tLxjqcRm7OJ9uUxFXmO7eR+uoCJx7qSwB0i6UmoAsBQIt3ZgVuAo70lFVPl7YjcJXWt9WxVs40zR53il1BIK/xW1eeKogmOxm0m3yzsWpVScUqqe1vps0b+98WVswy7el8s40+TnV4xTa33hoscrlVLTlVK1tdaO1uzKUd6TK7LGe+I0Qy5KKR+llF/xY+AOwOLVZQdQkYm57cFy4OGixw8Dpf77sOP3xZkmP7/isSil6iqlVNHjLhT+7ifZPNJr5yjvyRVZ5T0x+kpwBa8W303hX+YcIA5YU/R8fWBl0eMwCq/u7wH2Uzi8YXjsV3MsRV8PAKIprF6w12OpBfwCHC76HOhI74ul1xgYBYwqeqyAaUXL/6ScCiujPypwLE8Wvf57gK1AD6NjLuM4vgXOAnlFvyePOfB7cqVjqfL3RG79F0IIJ+E0Qy5CCFHdSUIXQggnIQldCCGchCR0IYRwEpLQhRDCSUhCF0IIJyEJXQghnMT/A94tnuHsuFZ8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_bins = 25\n",
    "bins = np.linspace(-1.5, 1.5, n_bins)\n",
    "\n",
    "# the bin number for each sample. Samples are binned into 25 bins\n",
    "bin_numbers = np.floor(trajectories/(100/n_bins)).astype(np.int32)\n",
    "\n",
    "# pass samples to the model with their corresponding bin numbers\n",
    "pmf = model.compute_PMF(dtrajs, bias_matrices, bin_numbers)\n",
    "pmf -= pmf.min()\n",
    "\n",
    "plt.plot(bins, pmf)\n",
    "\n",
    "# compare with the original distribution, shifted so that the minimum lies at zero.\n",
    "plt.plot(bins, potential(bins) - potential(bins).min(), 'k--');"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3a4872e3",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. rubric:: References\n",
    "\n",
    ".. footbibliography::"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}