Untitled2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate information theory measures from free-energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 810,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic IT class and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 869,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def islinspace(a,tol=1e-12): \n",
    "    \"\"\"check that a vector is linearly spaced\"\"\"\n",
    "    diffs = abs(a-np.linspace(min(a),max(a),len(a)))<tol\n",
    "    if not np.all(diffs):\n",
    "            print(a-np.linspace(min(a),max(a),len(a)))\n",
    "    return np.all(diffs)\n",
    "\n",
    "class free_energy(object):\n",
    "    \"\"\"\n",
    "    free-energy for arbirary dimensions\n",
    "    \"\"\"\n",
    "    \n",
    "    def read(self, filename=None, table=None, kT=1, verbose=True):\n",
    "        \"\"\"\n",
    "        reads a file with cv1 cv2 cv3 ... fes\n",
    "        or uses a table in the same format\n",
    "        \"\"\"\n",
    "        ## read file, id dimensions, identify values of colvars \n",
    "        ## check that they conform to standard format\n",
    "        ## reshape fes to array (n1,n2,n3, ...)\n",
    "        assert filename or table\n",
    "        if not filename==None:\n",
    "            tt = np.loadtxt(filename)\n",
    "        elif not table==None:\n",
    "            tt = np.array(table)\n",
    "            \n",
    "        nvars = (tt.shape)[1] - 1\n",
    "        edges_ = [np.unique(tt[:,i]) for i in range(nvars)]\n",
    "        nbins_ = [len(edge) for edge in edges_]\n",
    "\n",
    "        for edge in edges_:\n",
    "            assert islinspace(edge)\n",
    "        fes_ = (tt[:,nvars]).reshape(*nbins_)\n",
    "        bin_sizes = np.array([edge[1]-edge[0] for edge in edges_])\n",
    "        bin_volume = bin_sizes.prod()\n",
    "        \n",
    "        self.dim = nvars\n",
    "        self.fes = fes_\n",
    "        self.volume = bin_volume\n",
    "        self.edges = edges_\n",
    "        #calculates the normalised probability from fes\n",
    "        # note that normalisation is prob.sum()*self.volume = 1\n",
    "        f0 = -(self.fes)/kT\n",
    "        f0max = f0.max()\n",
    "        f1 = f0 - f0max\n",
    "        p0 = (np.exp(f1))\n",
    "        p1 = p0/(p0.sum()*self.volume)\n",
    "        self.prob = p1\n",
    "        self.kT = kT\n",
    "        \n",
    "        if verbose: self.describe()\n",
    "\n",
    "    def describe(self):\n",
    "        print(\"free-energy has %d dimensions\"%self.dim)\n",
    "        print(\"binning of free-energy: %s\"%str(self.fes.shape))\n",
    "        print(\"bin volume: %f\"%self.volume)\n",
    "        nor = (self.prob).sum()*self.volume\n",
    "        print(\"nomalization: %f\"%nor )\n",
    "        assert np.abs(nor-1.0) < 1e-8\n",
    "        \n",
    "    def H(self):\n",
    "        ## mask and calculate entropy\n",
    "        sum1 =  -(np.ma.log(self.prob)*self.prob).sum()\n",
    "        return(sum1*self.volume)\n",
    "    \n",
    "    def marginalise(self,setvars):\n",
    "        \"\"\"\n",
    "        marginal on set of variables\n",
    "        the order of setvars does not count anywhere,\n",
    "        order in original fes is preserved\n",
    "        \"\"\"\n",
    "        setvars.sort()\n",
    "        nvars = len(setvars)\n",
    "        bin_volume = self.volume\n",
    "        new_edges = [self.edges[i] for i in setvars]\n",
    "        old_bin_sizes = np.array([edge[1]-edge[0] for edge in self.edges])\n",
    "        \n",
    "        integration_vars = np.setdiff1d(range(self.dim),setvars)\n",
    "        bin_integration_vol = [old_bin_sizes[i] for i in integration_vars]\n",
    "        bin_integration_vol = np.array(bin_integration_vol).prod()\n",
    "        new_bin_vol = self.volume/bin_integration_vol\n",
    "        prob_ = (self.prob).sum(axis=tuple(integration_vars))*bin_integration_vol\n",
    "        \n",
    "        fes_ = (self.kT)*np.log(prob_)\n",
    "        fmarg = free_energy()\n",
    "        fmarg.dim = nvars\n",
    "        fmarg.fes = fes_\n",
    "        fmarg.prob = prob_\n",
    "        fmarg.volume = new_bin_vol\n",
    "        fmarg.edges = new_edges\n",
    "        return fmarg\n",
    "    \n",
    "    def plot(self,dims, show=True, **kwargs):\n",
    "        \"\"\"\n",
    "        plot 1d or 2d free-energyies\n",
    "        \"\"\"\n",
    "        if isinstance(dims,int): dims = [dims]\n",
    "        assert len(dims)<3\n",
    "        if len(dims)==1:\n",
    "            f1 = self.marginalise(dims)\n",
    "            plt.plot(f1.edges[0],f1.fes,**kwargs)\n",
    "            if(show): plt.show()\n",
    "        if len(dims)==2:\n",
    "            f2 = self.marginalise(dims)\n",
    "            #plt.imshow(f2.fes,interpolation='None', **kwargs)\n",
    "            plt.contourf(f2.edges[0], \n",
    "                       f2.edges[1], \n",
    "                       f2.fes, **kwargs)\n",
    "            if(show): plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### define some information measures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 870,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def conditional_entropy(fes, set1, set2):\n",
    "    #H(set1|set2) = H(set1, set2) - H(set2)\n",
    "    #print(set2)\n",
    "    f1 = fes.marginalise(set2)\n",
    "    f2 = fes.marginalise(np.union1d(set1,set2))\n",
    "    centropy = f2.H() - f1.H()\n",
    "    return centropy\n",
    "\n",
    "def MI(fes, set1, set2):\n",
    "    f1=fes.marginalise(set1) \n",
    "    f2=fes.marginalise(set2)\n",
    "    f3 = fes.marginalise(np.union1d(set1,set2))\n",
    "    mi = f1.H() + f2.H() - f3.H()\n",
    "    return mi\n",
    "\n",
    "def CI(fes,set1,set2,set3):\n",
    "    \"\"\"\n",
    "    calculates the 3-body co-infirmation \n",
    "    \"\"\"\n",
    "    # mutual info - conditional mutual info\n",
    "    i2 = MI(fes,set1,set2)\n",
    "    i2c = conditional_entropy(fes, set1, set3) + \\\n",
    "        conditional_entropy(fes, set2, set3) - \\\n",
    "        conditional_entropy(fes, np.union1d(set1,set2), set3) \n",
    "    return i2 - i2c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### generate multinormal free-energy for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 871,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def multinormal(cov_,nsample_=10000,nbins_=10):\n",
    "    ### generate a sample gaussian multivariate free-energy\n",
    "    #cov_ = np.array([[10,1,1],[4,5,2],[1,2,8]])\n",
    "    mean_ = 0.0*np.diag(cov_)\n",
    "    pts = np.random.multivariate_normal(mean_,cov_,nsample_);\n",
    "    probs, edges = np.histogramdd(pts, bins=nbins_)\n",
    "\n",
    "    # calculate bin volume and normalise prob\n",
    "    bin_sizes = np.array([edge[1]-edge[0] for edge in edges])\n",
    "    bin_volume = bin_sizes.prod()\n",
    "    probs = probs/(probs.sum()*bin_volume)\n",
    "\n",
    "    # calculate free-energy (kT units)\n",
    "    fes_ = ((-1)*np.log(probs))\n",
    "    it = np.nditer(fes_, flags=['multi_index'])\n",
    "    tab = []\n",
    "    while not it.finished:\n",
    "        #print(it.multi_index, fes_[it.multi_index])\n",
    "        ind = it.multi_index\n",
    "        pos = [edge[ind[i]] for i,edge in enumerate(edges)] \n",
    "        tab.append(pos + [fes_[ind]])\n",
    "        it.iternext()\n",
    "    #np.savetxt('fes1',tab)\n",
    "    return tab\n",
    "\n",
    "def gaussian_entropy(cov_):\n",
    "    \"\"\"\n",
    "    calculates gaussian entropy for a given correlation matrix\n",
    "    \"\"\"\n",
    "    if isinstance(cov_,np.ndarray):\n",
    "        det = np.linalg.det(cov_)\n",
    "        k = cov_.shape[0]\n",
    "    else: #for some reason numpy does not calc dets of 1d matrix\n",
    "        det = cov_\n",
    "        k= 1\n",
    "    return(0.5*np.log(det*(2*np.pi*np.exp(1))**k))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (1) calculate entropy for a 2D normal and compare with exact value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 872,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "free-energy has 2 dimensions\n",
      "binning of free-energy: (60, 60)\n",
      "bin volume: 0.103012\n",
      "nomalization: 1.000000\n"
     ]
    }
   ],
   "source": [
    "## genrate sample free-energy for testing\n",
    "cov_ = np.array([[10,4],[4,2]])\n",
    "np.linalg.matrix_rank(cov_)\n",
    "tab1 = multinormal(cov_,nsample_=300000,nbins_=60)\n",
    "\n",
    "f0 =free_energy()\n",
    "f0.read(table=tab1)\n",
    "\n",
    "# or read from file: np.savetxt('fes1',tab1); f0.read('fes1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 873,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numerical entropy: 3.540991\n",
      "theoretical entropy for gaussian 3.531024\n"
     ]
    }
   ],
   "source": [
    "print(\"numerical entropy: %f\"% f0.H())\n",
    "print(\"theoretical entropy for gaussian %f\" % gaussian_entropy(cov_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (1.1) same as above, convergence with bins and sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 886,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNXd+PHPNxsEQgibQFiSGGsRfiK41a1t0IIUl7o/\nYNCi1oUq2sfaRx9pDDFuVbS0tepji0VKEJda94VUSVVE5XEBRFyekCBkYZUlQNb5/v6YSZgkk8lM\nMnu+79drXnPn3nPPPZcJ5zv3nnPPEVXFGGOM6UhcuAtgjDEmslmgMMYY45UFCmOMMV5ZoDDGGOOV\nBQpjjDFeWaAwxhjjVaeBQkR6iciHIvKpiKwTkfwO0uW40nwuIivc1k8VkS9F5GsRuTWQhTfGGBN8\n4stzFCLSR1UPiEg8sBK4UVU/ctveH3gfmKKqFSIyWFV3iEgc8DVwBlAJrAamq+qXwTgZY4wxgefT\nrSdVPeBa7AUkAG2jy6XAP1S1wpV+h2v9icA3qrpJVRuAZcDPul1qY4wxIeNToBCROBH5FKgGilV1\ndZskRwIDRWSFiKwWkctc60cAm93SbXGtM8YYEyUSfEmkqg5gooikAi+IyFhV/aJNPscCpwN9gVUi\nsirgpTXGGBNyPgWKZqq619VQPRVwDxRbgB2qWgvUisg7wDFABTDaLd1I17p2RMQGnTLGGD+pqgT7\nGL70ehrsaqxGRJKByUDbxugXgdNEJF5E+gA/ADbgbLw+QkQyRCQJmA681NGxVDUmX/n5+WEvg52f\nnZ+dX+y9QsWXK4rhwJOuHkxxwNOq+pqIXAuoqj6uql+KyJvAWqAJeFxdt6ZE5AZguWvfhaq6IShn\nYowxJig6DRSqug5n+0Pb9f/T5vN8YL6HdG8A3+9GGY0xxoSRPZkdAjk5OeEuQlDZ+UU3Oz/TGZ8e\nuAsFEdFIKYsxxkQDEUEjoTHbGGNMz2aBwhhjjFcWKIwxxnhlgcIYY4xXFiiMMcZ4ZYHCGGOMVxYo\njDHGeGWBwhhjjFcWKIwxxnhlgcIYY4xXFiiMMcZ4ZYHCGGOMVxYojDHGeGWBwhhjjFcWKIwxxnhl\ngcIYY4xXFiiMMcZ4ZYHCGGOMVxYojDHGeGWBwhhjjFcJ4S6AMSY8yso2kZe3iIoKByNGxFFYOIus\nrIxwF8v4oPm7CxVR1ZAdzBsR0UgpizGxrqxsE5Mn/4nS0gKgL7Cf7Ox8iovnWLCIcK2/uxRUVYJ9\nTAsUxsSIxkbYvx/27YOamvbv7svLlhXw5Ze34AwSzfaTkTGfU0/NB0Bc1Y9I6+WO3n1JE6y0PenY\nixcX8L//2/zdSUgChd16MiYMVOHAgY4rcm/rOtpWXw8pKc5Xv37tl93XNTQ4aB0kAPrSt6+DadOc\n5Wsup/tyZ++BShNpx+5s31CeS1WVp+8uuCxQGOOD+vruV+Tu6/bvh169Oq7I3df16wfp6d4r/n79\nIDn50K/OzmzcGEdp6X7aXlFMnBhHbm4w/gVNoOzeHUdRUdvvLrjs1pPpskhtDHU4DlXKgarcm5oO\nVdq+VO6dbevbFxLC+DPN2iiil7VRREhZTOcCVdGoQm1tYH6lN78fPOisiLtSgXe0LSnJ91/r0aI5\n0FdWOkhPj5xAbzrX/N0VFc2zQBHtIvUXd1eoQkOD8776gQNw3XUFvPxy+8bQE06YzyWX5PtV8Sck\ndO/XedttffpAnD0hZHoAEWvMjmqefnF/8EHgL+3bVuDNr4MH26/r7noRZyXcpw/s3u25MXTLFgdV\nVc5Ke/hw3yr3xMSA/XMYY4LAAkWQ5OUtcgsSAH0pLS3guuvmc9tt+QGtxN0r8Lav5GTP6/v2hSFD\nfE+fnNy6Qp8501OD2n5OPz2OBx8M3b+zMSb4LFAESUWF51/cq1Y5uPNO/yvwjirxthV4qBQWzuKD\nD/LbtVEUFs4JfWGMMUFlgSIIGhth3744oP0v7nPPjWPJkjAVLICysjIoLp5DXt58t8ZQ6zFjTCyy\nxuwAW7UKZs+GPn02sWXLn9i82bofGmOCI1SN2RYoAmTHDrjtNnj9dZg/H6ZPh/Jy635ojAmeiAkU\nItILeAdIwnmr6jlVLWiT5sfAi8BG16rnVfUu17ZyYA/gABpU9cQOjhOVgcLhgCeegLlzYcYMKCiA\n/v3DXSpjTE8QMd1jVbVORCap6gERiQdWisjrqvpRm6TvqOq5HrJwADmq+l0gChxJPv0UfvlL5/Kb\nb8KECeEtjzHGBINPjyWp6gHXYi+cwcXTT/+Oopr4epxosWcP3HQTTJ0Kv/gFrFxpQcIYE7t8qsBF\nJE5EPgWqgWJVXe0h2cki8pmIvCoiY93WK1AsIqtF5OoAlDlsVGHpUjjqKOfzC+vXw1VX2VPAxpjY\n5lP3WFV1ABNFJBV4QUTGquoXbkk+Bka7bk/9FHgBONK17VRVrRKRITgDxgZVfc/TcebNm9eynJOT\nQ05Ojt8nFCwbNsD118N338E//gEnnxzuEhljepqSkhJKSkpCfly/ez2JSB6wX1Uf8pKmDDhOVXe1\nWZ8P7PO0b6Q2Zu/fD3fdBX/9K+TlOdskwjnqpzHGNAtVY3anN01EZLCI9HctJwOTgS/bpBnqtnwi\nzgC0S0T6iEiKa31fYArweQDLHzSq8MILMG4cfPstrF0LN95oQcIY0/P4Uu0NB54UkTicgeVpVX1N\nRK4FVFUfBy4SkdlAA3AQ+A/XvkOBf4qIuo5VpKrLA34WAbZxozMolJY6u76efnq4S2SMMeFjD9y5\nqauDBx6ABQvgllvg5pud8xAYY0wkKSsvI++hPIr+VBQZz1H0FMXFzsbqcePg448hwx6g7lTzH2vF\n3gpGpI6g8OZCsjKzwl0sY2JaWXkZk2+YTOkxpSE7Zo8PFBUVziuH1avhj3+Es88Od4miQ6s/1kFA\nPXxwwwcUP1xswSJKWKDvGlWl0dFIg6PB+d7UQIOjgYamhpb1bZfbpu3Ofq//9XXn/7sQ3u3osYGi\noQH+9Ce45x5nT6ZFi5xDdhvf5D2U1/qPNQlKjynl5t/dzMO/e5g4iUNEEAQRcX52LXe0zpd9TGCE\nItA71BGyyrPT/bqRZ9v8m7SJeIknMT6RxLhEEuISWpYT412fPSy3TdtqPy/beif0JiUppSUfQaAG\njlgB/xeQb6pzPbKN4r33nMFh2DB4+GE48sjO9+mJVJWdB3eyafcmNu3ZdOh9zyaKFxaz/9T97fZJ\neieJQdMGoSgOdaCqKIqq67NruaN1bfdpXtess2DT1QDU2T6+5NOVfYJVXo/7cCjNikUrKD2qza/S\nehi9fjTHX3p8QCpyhzq6V3l2tF8w8vQjj4S4hLD+aDn/5+ehL71I0W5IAWujCLTt2+G//svZHvH7\n38NFFzlnh+upHOqguqaaTbs3Ub67vF0w2LR7EwlxCWSkZZDR3/nKTMvktNGncTDjIG/Wv9muorl4\n3MUs+XVwJtzoKJh0JQB1J2gFep9wlPcdeaf9rYskSElMYcb/mxGQSjde4u0q0BeqUF8PtbXOHjXN\nr7afXevGfLGF3+5uPy1aMPWIQNHUdOiBucsucz5l3a9fuEsVfA1NDWzZu6VVAHAPCFv2bqF/7/7O\nIOAKBmOHjGXa96a1fO7f2/NQuMfNPe7QrYskoB6y12RT+HBh0M7H/dez6Z53RrzD19u+5oh3Ydg+\nqO4H//dDmDh8IheNvSjcxQu+xkafKuVO1wUiTX29s3tlr16HXr17t/7sti5p07chDRLQA249ffyx\ncyKhpCR49FE4+uiAHyJsDjYcbH0V0OZqoLqmmmEpw1pdETQvZ6ZlMrr/aJITu94w09wYWrm3kvTU\ndGsMjSLvvfsOf5h6BosONLqm1YJZfRK46Y23OO2HPwr8AZt/NQezwvVnP4ejdWXspWIOepqkJL8G\njCuYOZNbioroi3PE1YiYjyJUAh0odu+G3/4WnnsO7rsPLr889IP3dbdXyZ7aPS2VfsuVgFtA2FO7\nh5GpI8lMy2wVBJrfR6aOJDE+DBNqm4hXMGMGtyxb1maiXph/0knk//znga/M6+udk7tHQsXcq5dz\niIUovS22qayMP02eTEFpqbVRdJUqLFnibIs47zz44gsYODD05eisV4mqsv3A9nZXA+63hhodja2u\nADL6Z3B8+vEtwWBYyjC7DWPaa2iAykrYvBm2bHG+t3k5tm1rd/uiL+AoL3dOtOJeySYnQ1pa9yrm\nXr1smOUAycjKYk5xMfPz8qCoKCTHjKkrivXrnb2Z9u933mY64YQAFa4LZt44k6J+Re0ae4evHU7/\nqf3ZtHsTyYnJ7W4HuX8emDzQGgNNa01NUF3dvvJ3Dwjbt8Nhh8GoUc7XyJGHll2vgltu4Zannmp/\nRZGbS/6S4HRGMIEXMTPcRYOaGigshL/9DebNg2uvhfj48JapYm+F80rCXRIM6TOEpRcvZXT/0fTr\n1QNa1I3vHA7Yts1z5d/8qq52XiK3qfg56aRDAWH4cOdtHi9m3X03+R99REFpaUsbRX52NnMKg9cZ\nwUSvqA4UqvDPf8KvfgWTJsG6dTB0aOf7hcKI1BFQT7sriqOHHs24w8aFq1gmXFRh506Pt4FaAkJF\nBaSmtr8SmDDh0LoRIwIyAJn77QtHZSVx6enMKSwkI8s6IxgPVDUiXoA6/ze1fuXna4uNG8s1N3ee\n5uTcoUceuaLT9O7y89unDUV6TstX5qHcjmafla0byzaGtTyWPgjpHQ7VXbtU16xRfeUV1cce0/wf\nvuU5/WGPqE6bpnrttap33aX65JOqb72l+TfsiJ7ztfQRk95ZhQe/fo6aNoqysk1MnvwnSksLwHWx\nPHBgPqtWzeHIIzNCVk5ffFTxEWctPYuFP1rIM0XPWPfRaLdvn/crgc2bnQ21XtoEGDkSUlLCfSYm\nRpSVl5P36KMU3X8/at1jD5k5s4CiolugTfNbbu58lizJD3r5fFX2XRmnPnEq/3P2/3DO988Jd3FM\nZw4c8N4msHmzsxdR24q/bUDo7/nBxEjWXNlU1NYyondvCmfPJiszM9zFMp0oKy9ncn4+pdOnw7Rp\nIQkUUdNGUVHhoP1D632prHR4Sh4W3x38jmlLp/Hfp/23BYlIUFvrvO/v7Wpg//5DFX7z+8SJcO65\nh4LAgAFR2+e+I60qm+RkOHiQD/LzKS4osGDhgariABzdfA9EPvf84Q+HvrcQiZpAMWJEHM6+Ga2v\nKNLTI6Nvdn1TPRc8cwFnZp/JnB/MAZwPxizKy8NRUUHciBHMssbCwGloaB0EPF0N7N4N6emtrwLG\njYMzzzz0eciQmAsCbTU4HOxpbGRPU5PzvbGReQsWtK5skpMpnT6d8+fP5+Lf/MZjBaV0v6Jsl18w\n8gxCGcE5vWecSJffJQB5xAGf7d0b8qGuoyZQFBbO4oMP8lu1UWRn51NYOCfMJXP+SvjFS7+gf6/+\nPDjlQaD105Mt3Q8/+IA5xcUWLDrT2Oj5WQH3gLBjh7OLm/ttoOxsyMk5tG7o0Kh/yKvR4WhVwe92\nvbuv6+xzncNBakIC/RMS6B8fT/+EBNbX1LSvbJKT+a6+ngNNTa0qp3jXe3cruOZ3CVA+waiE2+Xp\nVuZIMXPwYIoOHrQrCk+ysjIoLp5DXt58KisdpKfHUVg4h6ys8Ddk3/nvO9mwYwMlPy8hPs75AMei\nvLyWIAHO0FZQWsr8q68m/847Dz296uk9SoYX6NIVk8MBW7d6vxLYuhUGDWrfJnDKKa2fFUiI7D/f\ntpV82wrcl0rfUyXfv83noUlJHOklTd/49qO4zhw0qH1lc/AgPxwwgLsPPzzE/1LGH4WzZ/NB823D\nEImaxuxItXjNYvJL8ll11SqGpQxrWZ8/aRIFJSXt0uenpVEwZkzrcXHavjscHQcRX98Dsa+XgOXx\niik7mzlPPUVGfHzHVwKVlc6GX289hNLTwz5ZeWeVvC+Vvi+VfNvPaT5U8oHgqY0ie9kya6OIEtbr\nKYqUlJdwybOXUDKrhLFDxrba5j7CYzOfh0hwHwK5o2Diy3t39u0kYBVs2cItO3e2P7+kJPLHjOm4\nh9DIkc58gqjR4WBv20rcz0q/trmSj49vXXn7UemnBKmSD5TmyqaytpZ06/UUlUI1hIcFii7asH0D\nOU/msPSCpZxx+Bnttm8qK+NP48dTUFPT+hd3NLVReAlY+VddxaxPPiFvxAgqBg5kxK5dFFZUsGjS\nJArefrvrh2xTybet0H2p9GsdDvp1UIH7WulHeiVvDNhYTxFta81Wzlp6Fr/7ye88BgmAjG+/ZU7f\nvsyfNg3H9u3ROURCQoLz1bf9NCm7R43ijORkyubObbl18f7ddzN5yBDKDh70qcHVU6XvrZJvXh6c\nmMgRyclWyRsTInZF4acDDQeY9OQkzsw+kzsn3ek5UUODsy/+vHnO+VajWL3Dwdb6eirr66mqq6PK\ntbz4nnv49oIL2jWGJjzzDCOuucbnWzVpVskb02V2RRGBHOrgsn9expGDjqQgp6DjhH/4g/Ne/IUX\nhq5wfqpzOFpV/B0t72ls5LDERIb36sXwpCTSk5IY3qsXvePjPXav/GFqKm+ffHJ4TsoYExQWKPzw\nX8X/xY4DO1g+c3nHv3q3bHFOqbdqVVi6uB5oaqLKS8XfvFzT1MRQt4q/OQic1r9/q4AwODGReA/n\n8U3//nztoXtlepAbqo0xoWe3nnz0yOpH+OOHf+T9q95nYLKXKfMuuQTGjIE77wzoWDo1jY2d/vqv\nqquj1uFo9+vf0/LAxETiuhHIrHulMeFnvZ4iyKtfv8ovXv4FK69cyeEDvDyMtHw5XHcdrF9P2dat\nnVakqsrepiafbgE1qXqt+JuXByQkhOwev3WvNCa8LFBEiE+qPuHMJWfy8oyXOWnkSR0nrKuDo4+G\nhx6Cs89m5q23UpST0+7WzOgXX2T0tde2BAEB0psrfNe7p+VUa+Q1xrRhjdkRYPOezZz71Lk8dtZj\n3oMEwAMPwNixcPbZAFTU1nps7E0R4a6srJYg0C/Ch6EwxhirpTqwt24vZy09i1+d9CsuHNtJ76Wy\nMliwAD7+uGXViN69wUNj78TUVH6clhakUhtjTOBF99CaQdLQ1MDFz17MqaNO5dcn/7rzHW66CW6+\nGTIODVBYOHs2mcuWOYMFtLRRFM6eHaRSG2NMcFgbRRuqyrWvXMuWvVt4acZLJMR1ctH10kvwm9/A\n2rXOcZDc3LpyJc8tWUJGUpI19hpjAs7aKMLk/pX3s7pyNe/MeqfzIHHggPNq4i9/aRckVJUXEhL4\n+333cUoUTpNpjDHNLFC4efrzp/nz6j+z6qpV9OvVr/Md7rkHfvAD+MlP2m1auWcPcSKcnJoahJIa\nY0zodNpGISK9RORDEflURNaJSL6HND8Wkd0i8onr9Vu3bVNF5EsR+VpEbg30CQTKym9XMuf1Obw8\n42VGpI7ofIevvoLHHoMHH/S4eWF1NVcNG2ZdWo0xUa/TKwpVrRORSap6QETigZUi8rqqftQm6Tuq\neq77ChGJAx4GzgAqgdUi8qKqfhmoEwiE/9v1f1z4zIUsPn8xxww7pvMdVOGGG+D222FE+6Cyt7GR\nF3bs4HcnnhiE0hpjTGj51OtJVQ+4FnvhDC6eWp09/XQ+EfhGVTepagOwDPhZVwoaLDsP7GRa0TTu\nnHQnU4+Y6ttOzz7rnNN5juf5up/eto1JaWkcFuZZ2owxJhB8ChQiEicinwLVQLGqrvaQ7GQR+UxE\nXhWR5uneRgCb3dJsca2LCLWNtZz39HmcP+Z8rjnuGt922rfP2RX20UchMdFjkoVVVVw1fHgAS2qM\nMeHj6xWFQ1UnAiOBH7gFgmYfA6NVdQLOW00vBLaYgedQB1e8eAXDU4Zz70/u9X3HggJn4/Vpp3nc\n/HlNDVvq6jhzwIAAldQYY8LLr15PqrpXRFYAU4Ev3NbXuC2/LiKPiMhAoAIY7ZbFSNc6j+bNm9ey\nnJOTQ05Ojj/F80ve23ls2r2Jty5/izjx8bnDzz+HxYud7x1YWF3NrGHDSIizZxmNMYFVUlJCSUlJ\nyI/b6QN3IjIYaFDVPSKSDLwJ3Keqr7mlGaqqW13LJwLPqGqmq/H7K5yN2VXAR8AMVd3g4Tghe+Bu\n4ScLufe9e1l11SqG9B3i206q8OMfw/Tp8MtfekxS53AwatUqPjj2WA5vO86TMcYEWCQ9cDcceNLV\ngykOeFpVXxORawFV1ceBi0RkNtAAHAT+A+fGJhG5AVju2nehpyARSsWlxcx9ey7vXPGO70EC4O9/\ndz5gd+21HSZ5accO/l/fvhYkjDExpUcN4fH5ts85/cnT+ccl/+CHGT/0fcfdu+Goo5zDdZxwQofJ\npq5Zw2XDhpE7dGgASmuMMd6F6oqix9xIr9xXyVlLz2LB1AX+BQmA3/4WfvYzr0Hi29paVu/bxwWD\nB3ezpMYYE1l6xBAeNfU1nPPUOVxz7DVcevSl/u388cfw3HPwxRdeky2qrmb6YYeRHB/fjZIaY0zk\nifkriiZHE5f+41ImDJ3A7T+83b+dHQ5nw/U998DAjufJdqjyRFUVv7BnJ4wxMSimA4Wq8qs3fsXB\nxoM8dvZj/o+79Ne/QkICzJrlNdlb333HwMREJvbzYSBBY4yJMjF96+kPH/6BFeUrWHnlShLjPT9F\n3aEdOyAvD5Yvh06eibAnsY0xsSxmryhe+PIF5r8/n1cvfZX+vbswH8Rtt8GMGXCM90ECdzY08Mau\nXVx62GFdLKkxxkS2mLyi+KjiI65++WreyH2DjLSMzndo6/334fXXO23ABijaupWzBg1iQAfjPhlj\nTLSLuSuK8t3lnP/0+Txx7hMcl36c/xk0NjobsB94ADqZmU5V7baTMSbmxVSg2F27m2lF07jt1Ns4\n5/vndC2TRx5x9nCaMaPTpB/v28e+piZy0tK6dixjjIkCMfNkdn1TPVOXTGX80PEsmLqga5lUVcH4\n8fDOO84nsTtx3VdfMap3b+ZmdOH2ljHGdJM9me0HVeWal68htVcqD07xPDWpT37zG7jqKp+CxIGm\nJp7Zvp1Zw4Z1/XjGGBMFYqIx+6537mL99vWU/LyE+LguPhm9YgW8+65PDdgAz23fzsmpqYzo1atr\nxzPGmCgR9YFiydolPPHZE6y6ahV9k/p2LZP6erj+evj976Gvb3ksrKrippEju3Y8Y4yJIlEdKP5d\n/m9+vfzXrPj5CoaldOMW0IIFkJkJ55/vU/JvDhzgywMHOHvQoK4f0xhjokTUBoovd3zJJc9dwtIL\nljJ2SNuZWf2weTPcfz98+CH4OMTHE9XVXDZ0KEk2i50xpgeIqkBRVl5G3kN5lH1Xxtpta8m7MY8z\nDj+je5n+53/CDTdAdrZPyRsdDp6sruatTp7YNsaYWBE1gaKsvIzJN0ym9JhSGASMhscffZyLx11M\nVmZW1zJ94w349FPn7HU+em3XLrJ69+YoH9syjDEm2kXNvZO8h/KcQSLJtSIJSo8pJe+hvK5lWFvr\nvJJ4+GHwY+pSexLbGNPTRE2gqNhbcShINEuCyr2VXcvw/vudD9f99Kc+71JVV8c7e/ZwyRA/5to2\nxpgoFzW3nkakjoB6WgeLekhPTfc/s9JS+OMf4ZNP/Npt8datXDh4MCkJUfPPZowx3RY1VxSFNxeS\nvSbbGSwA6iF7TTaFNxf6l5Eq3Hgj3HILjB7tx27OWezstpMxpqeJmp/GWZlZFD9cTN5DeVTurSQ9\nNZ3Chwv9b8h+8UXYuBH++U+/dntvzx7iRTgpNdW/4xljTJSLmUEBfbJ/P4wbB088Aaef7teuszZs\n4OiUFH49alSQCmeMMf6xQQGD4e674ZRT/A4SexobeWHHDi4bOjRIBTPGmMgVNbeeuu3LL+Hxx2Hd\nOr93XbZtGz8ZMIDDktp2uzLGmNjXM64oVJ2D/v32t9CFxmh7dsIY05P1jEDx9NOwY4fzATs/raup\noaq+nikDBwahYMYYE/li/9bT3r3OrrBPPw1deP5hYVUVs4YNI97HAQONMSbWxP4Vxbx5MGUKnHqq\n37vWORwUbdvGFTaLnTGmB4vtK4q1a2HJEli/vku7v7hjB8f07cvhfowFZYwxsSZ2rygcDvjlL+HO\nO6GLYzP91RqxjTEmhgPF4sXOKU6vvrpLu2+qreWTffs4f/DgABfMGGOiS2zeevruO7jtNnjlFYiP\n71IWf6uqYsbQofTu4v7GGBMrYvOKYu5cuOACOP74Lu3epMrfqqu5yhqxjTEmBq8oVq+G55+HDRu6\nnMVb333H4MREJvTrF8CCGWNMdIqtK4qmJmcD9n33wYABXc7GnsQ2xphDOg0UItJLRD4UkU9FZJ2I\n5HtJe4KINIjIBW7rykVkjWv/jwJVcI/+8hfo3Rsuv7zLWeyor+fNXbu49LDDAlgwY4yJXp3eelLV\nOhGZpKoHRCQeWCkir6tqq0pfROKA+4A322ThAHJU9buAldqT7dvhjjvgrbcgrusXSku2buWcwYNJ\nS0wMYOGMMSZ6+VSjquoB12IvnMHF08QRc4DngG1t1ouvx+mWW2+FmTPh6KO7nIWqstAasY0xphWf\nGrNdVwsfA9nAn1V1dZvt6cB5qjpJRE5ss7sCxSLSBDyuqn8JQLlbW7kSli+HL77oVjar9+3jYFMT\nP05LC1DBjDEm+vkUKFTVAUwUkVTgBREZq6rutfIC4Fa3z+4j6J2qqlUiMgRnwNigqu91u+TNGhud\nDdjz50M3pyldWFXFlcOHIzYAoDHGtPCre6yq7hWRFcBUwD1QHA8sE2cNOxj4qYg0qOpLqlrl2ne7\niPwTOBHwGCjmzZvXspyTk0NOTk7nhXr4YecQHf/xH/6cSjv7m5p4dvt21p1wQrfyMcaYYCkpKaGk\npCTkx+10zmwRGQw0qOoeEUnG2Vh9n6q+1kH6vwEvq+rzItIHiFPVGhHpCywHClR1uYf9/J8zu7IS\nxo933nr6/vf927eNJ6ureXbbNl4ZP75b+RhjTKiEas5sX64ohgNPutop4oCnVfU1EbkWUFV9vE16\n99p+KPBPEVHXsYo8BYkuu+UWuOaabgcJcN52unnkyAAUyhhjYkunVxSh4vcVxVtvwVVXOYcQ79u3\nW8f+6sBwVzHSAAAZiElEQVQBfvzpp2w++WQSu9G11phYkpmZyaZNm8JdDANkZGRQXl7ebn0kXVFE\nnvp657Smf/hDt4MEwBNVVVw+bJgFCWPcbNq0iUj5IdnThbuDTVQFik1lZSzKy8OxahVxtbXMOvpo\nMrqZZ4PDweKtW1lxzDEBKaMxxsSaqAkUm8rK+NPkyRSUltIX2A/kT5nCnOJiMrKyupzva7t2kd27\nN2MCcGVijDGxKGrutSzKy2sJEgB9gYLSUhbl5XUrXxsA0BhjvIuaQOGoqKDtb/6+gKOysst5VtXV\n8e6ePVzcxalSjTGmJ4iaQBE3YgT726zbD8Slp3c5zyerq7l4yBBSEqLmDpwxxiUrK4u3337b47b3\n3nuPo446KsQlil1REyhmFRaSn53dEiz2A/nZ2cwqLOxSfi0DANptJ2NizmmnncaGbkxe5qvTTz+d\nww47jLS0NCZOnMhLL73UavvSpUvJzMykX79+XHDBBezevTvoZQqGqHqOoqXXU2UlcenpzCos7HJD\n9r937+b6r79m3QknhL3rmTGRyNVH3+O2srJN5OUtoqLCwYgRcRQWziIry78+iN3NIysri4ULF3L6\n6af7ddxAWrduHWPGjCExMZGPPvqIn/zkJ3zzzTcMHTqU9evXc/LJJ/P6668zceJErr76ahwOB089\n9ZTfx+nouwjVcxSoakS8nEUJncu++EIf+vbbkB7TmGjS0f/JjRvLNTv71wo1CqpQo9nZv9aNG8t9\nzjsQeWRmZuq9996rY8eO1YEDB+qVV16pdXV1qqpaUlKiI0eObJV2/vz5On78eE1LS9Pp06e3pN2x\nY4eeffbZmpaWpgMHDtQf/ehHPpfB3YcffqjJycm6evVqVVW9/fbbNTc3t2V7aWmpJiUlaU1Njd95\nd/RduNYHvX6OmltPgbSnsZGXduzgsqFDw10UY6JOXt4iSksLwK0PYmlpAXl5i0KaBzhv7RQXF1Na\nWspXX33FXXfd1bKt7Z2CZ599luXLl1NWVsaaNWtYtMh5rAcffJBRo0axc+dOtm3bxj333NOyz/XX\nX88NN9zgtQznnHMOycnJnHTSSeTk5HD88ccDsH79eo5xez7r8MMPp1evXnz99dd+nWMk6JGtuE9t\n3crkgQMZnJQU7qIYE3UqKhzgoQ9iUZGDoiJfc/GcR2Wlw6+yzJkzh3RXh5a5c+dy4403cuedd3pM\ne9NNNzHU9ePwnHPO4bPPPgMgMTGRqqoqysrKyM7O5tRTT23Z589//nOnZXj55ZdpamriX//6V6t2\nkZqaGvr3798qbWpqKvv27fPrHCNBj7yisFnsjOm6ESPiwEMfxNzcONR1I6mzV26u5zzS0/2rkka6\nDeSZkZFBpZfu8kPd7iD06dOHmpoaAH7zm9+QnZ3NlClTOOKII/jd737nVxkA4uPjOfPMM3nzzTd5\n5ZVXAEhJSWHv3r2t0u3Zs4d+/fr5nX+49bhAsbamhq319UweODDcRTEmKhUWziI7Ox/c+iBmZ+dT\nWDgrpHkAbN68uWV506ZNLVcX/khJSWH+/PmUlpby0ksv8dBDD7FixQq/8wFobGyktLQUgHHjxrFm\nzZqWbaWlpTQ0NHDkkUd2Ke9w6nG3nhZWVXHFsGHEW08nY7okKyuD4uI55OXNp7LSQXp6HIWFc/zs\nsdT9PMB5a+iss84iOTmZe+65h+nTp/t7Orz66quMGTOG7Oxs+vXrR0JCAnE+DBD61VdfUVZWRk5O\nDgkJCSxbtox3332XBx54AIDc3FxOOeUUVq5cyYQJE7jjjju48MIL6RuFwwX1qEBR29RE0dat/O9x\nx4W7KMZEtaysDJYsyQ9rHiLCpZdeypQpU6iqquK8885j7ty5HabtyDfffMMNN9zAjh07GDBgANdf\nfz0//vGPAZg9ezYiwiOPPNJuP1Vl3rx5bNiwgfj4eL73ve/xzDPPMGHCBADGjh3LY489xqWXXsqu\nXbuYPHkyTzzxRJfPN5yi6jmK7lq2dSsLq6sptpFijemUt+coTGiF+zmKHtVGYY3Yxhjjvx4TKMoP\nHuTTffs4b/DgcBfFGGOiSo8JFH+rrubSoUPpHR8f7qIYY0xU6RGBokmVv1VX8wsbANAYY/zWIwLF\nv777jqFJSYxPSQl3UYwxJur0iEDx16oqa8Q2xpguivlAsb2+nuJdu5hhAwAaY0yXxHygWLJ1K+cO\nHkx/m8XOGGO6JKYDhaqysKrKZrEzJgbZVKihE9OB4qN9+6hT5Udthvo1xsS2UE2FmpmZSZ8+fUhN\nTSU1NZWpU6e22h4rU6HGdKBY6GrEtqlOjQmssvIyZt44k0mzJjHzxpmUlZeFJY9wExFeffVV9u7d\ny969e3njjTdatq1fv57rrruOoqIitm7dSnJyMrNnzw5jabsuZgNFTWMjz23fzs+tt5MxAVVWXsbk\nGyZT1K+IkqwSivoVMfmGyX5V9IHIA+Cjjz5i3LhxDBo0iKuuuor6+noA/v3vfzNq1KiWdFlZWTz4\n4IMcc8wxDBgwgBkzZrSk3blzJ+eccw4DBgxg0KBBLQMC+qqj8bCWLl3Kueeey6mnnkqfPn0oLCzk\n+eefZ//+tvNwRL6YDRTPbt/Oaf37M7xXr3AXxZiYkvdQHqXHlELzBJFJUHpMKXkP5YU0D4iMqVBz\nc3MZOnQoU6dOZe3atS3rbSrUKLCwqorfjB4d7mIYE3Mq9lbAoDYrk6BobRFFBT7OhboWmNQ+j8q9\nHc9Q50m4p0JdunQpxx57LKrKggULOPPMM/nqq69ITU21qVAj3Zf791NaW8s0m8XOmIAbkToC6tus\nrIfc8blovvr0yh2f6zGP9FT/ZqgL91SoJ598Mr169aJ3797cdtttpKWl8e677wI2FWrEe6K6msuH\nDiXRh1mqjDH+Kby5kOw12Ycq+nrIXpNN4c2FIc0DIm8qVPd5I2JpKtSYq0kbHA4WV1fbsxPGBElW\nZhbFDxeTuy+XSWWTyN2XS/HDxWRlZoU0D3DeGqqoqGDXrl3dmgq1eZ5rf6ZC3bx5M++//z4NDQ3U\n1dXxwAMPsHPnzpZbV7m5ubz88susXLmS/fv321SokeTVnTs5sk8fjuzTJ9xFMSZmZWVmseSPS8Ka\nR7inQt23bx+zZ89m48aN9O7dmwkTJvDGG28wYMAAwKZCDYpATYV6zrp1XDRkiHWLNaabbCrUyGFT\noQZQRV0dK/fs4aIhQ8JdFGOMiRmdBgoR6SUiH4rIpyKyTkTyvaQ9QUQaROQCt3VTReRLEflaRG4N\nVME9ebK6mouHDKGvzWJnjDEB02mgUNU6YJKqTgQmAD8VkRPbphOROOA+4M026x4GzgTGATNEZEyA\nyt6KQ5UnbABAY4wJOJ9uPanqAddiL5wN4J5uXM4BngO2ua07EfhGVTepagOwDPhZ14vbsXd27yY5\nPp4TorCPsjHGRDKfAoWIxInIp0A1UKyqq9tsTwfOU9VHAfeGlRHAZrfPW1zrAm6ha05sGwDQGGMC\ny6fusarqACaKSCrwgoiMVdUv3JIsALrd/jBv3ryW5ZycHHJycnzab3dDA6/s3MmCI47obhGMMSZi\nlZSUUFJSEvLj+t09VkTygP2q+pDbuo3Ni8BgYD9wDc7bUPNUdaor3W2Aqmq7Z+S70z320YoKSnbv\n5ulx47q0vzGmPeseGzkivnusiAwWkf6u5WRgMvClexpVPdz1ysLZTvFLVX0JWA0cISIZIpIETAde\nCvRJ/NUasY0xJmh8aaMYDqwQkc+AD4E3VfU1EblWRK7xkL4l7KlqE3ADsBxYDyxT1YBOO/XZvn3s\naGjgJ66nIY0xPYNNhRo6vnSPXaeqx6rqBFUdr6p3u9b/j6o+7iH9lar6vNvnN1T1+6r6PVW9L7DF\ndzZiXzFsGHHWiG2McQnFVKibN2+mX79+LdOg9uvXj7i4OH7/+9+3pLGpUCNAbVMTT23dyhV228mY\nkNpUVkbBzJnkT5pEwcyZbCrzfxrTQOQRTqNGjWLfvn0t06CuW7eO+Ph4LrroIiC2pkJFVSPi5SyK\nf5ZWV+uUzz7zez9jTOc6+j9ZvnGj/jo7W2tAFbQG9NfZ2Vq+caPPeQcij8zMTL333nt17NixOnDg\nQL3yyiu1rq5OVVVLSkp05MiRrdLOnz9fx48fr2lpaTp9+vSWtDt27NCzzz5b09LSdODAgfqjH/3I\n5zK4mzdvnp5++uktn2+//XbNzc1t+VxaWqpJSUlaU1Pjd94dfReu9UGvn6P6imKhNWIbE3KL8vIo\nKC2lebDsvkBBaSmL8nyfxjQQeUBkTIXa7O9//zuzZs1q+RxLU6FGbaAoO3iQNfv387PBg8NdFGN6\nFEdFBW1nVOgLOIqKQMSnl6OoyHMeXmao86R5KtS0tDTmzp3LU0891WHa5qlQ09LSOpwKNT4+vt1U\nqA8//HCn5Xj33XfZtm0bF154Ycs6mwo1AjxRXU3uYYfRy2axMyak4kaMYH+bdfuBuNxc142kzl9x\nubme8/BzhrpwT4XabPHixVx44YX0cZsHx6ZCDbMmVRbZLHbGhMWswkLys7NbKvr9QH52NrMKfZ/G\nNBB5QGRMhVpbW8uzzz7b6rYTxNZUqFE5w93yXbsYnpTE0Skp4S6KMT1ORlYWc4qLmZ+Xh6Oykrj0\ndOYUFpKR5fs0poHIA5y3hs466yySk5O7NRXqmDFjyM7O9msq1GbPP/88AwcObJkVr1lubi6nnHIK\nK1euZMKECTYVaqiUlZeT9+ijvLlzJ1l9+lA2aBBZmZnhLpYxPU5GVhb5S7o3FWp38wj3VKjNFi9e\nzOWXX95uvU2FGgSdjfVUVl7O5Px8SqdPh+RkOHiQ7GXLKC4osGBhTBDYWE+RI+LHeooUeY8+eihI\nACQnUzp9OnmPPhreghljTIyLmkBRUVt7KEg0S06msrY2PAUyxpgeImoCxYjeveHgwdYrDx4kvXfv\n8BTIGGN6iKgJFIWzZ5O9bNmhYOFqoyiM1rFTjDEmSkRNYzYc6vVUWVtLeu/eFM6ebQ3ZxgSJNWZH\njnA3ZkdVoDDGhI4FisgR7kARNbeejDHGhIcFCmOMMV5ZoDDGRCWbCjV0LFAYY2JOKKZCBbjjjjsY\nP348iYmJ3Hnnne22e5sKtb6+niuvvJL+/fuTnp7eagrVSGOBwhjjt7LycmbeeiuTbrqJmbfeSll5\neVjyCLfvfe97PPDAA5x99tnttnU2FWp+fj6lpaVs3ryZt99+m/vvv5/ly5eHsvi+C8U0er686MJU\nqMaY4Ono/+TGsjLNvvxy5bXXlBUrlNde0+zLL9eNZWU+5x2IPCJpKtSZM2dqQUFBq3WdTYWanp6u\n//rXv1q233HHHTpjxgyP+Xf0XWBToRpjIlEgxl0L1NhtkTQValvepkLdvXs3VVVVjB8/vmX7Mccc\nw/r167t0rGCLqmHGjTHh19G4a0XV1RSVlPiWSXV1QMZua54KFWDu3LnceOONHtsK4NBUqECHU6Fm\nZ2e3mwq1q7xNhVpTU4OItNoeydOkWqAwxvilZdw194r+4EFyhw1jSU6OT3nMfP11ijzk4e/Ybd2Z\nCrWqqgpwToU6b948pkyZgohw9dVXc+utt/pVDk+8TYWa4pp0be/evQwePLjVtkhkt56MMX4JxLhr\ngRq7LRKmQu2It6lQ09LSGD58eKvta9asYdy4cd0+bjDYFYUxxi9ZmZkUFxS0HnfNzwnEApEHhH8q\n1MbGRhobG3E4HDQ0NFBXV0diYiJxcXGdToV62WWXcdddd3HcccdRVVXFX/7yFxYvXux3+UMiFC3m\nvrywXk/GRJRI/z+ZlZWl9913n44dO1YHDBigV1xxhR48eFBVnb2eRo0a1SrtW2+91fJ53rx5etll\nl6mq6u9//3vNzMzUlJQUHTVqlN59990t6a677jqdPXt2h2WYNWuWiojGxcW1vJ588smW7U899ZSO\nHj1aU1JS9Pzzz9fvvvuuZVtdXZ1eeeWVmpqaqsOGDdMFCxZ0eJyOvgtC1OvJBgU0xnhkgwJGDhsU\n0BhjTESzQGGMMcYrCxTGGGO8skBhjDHGKwsUxhhjvLJAYYwxxit74M4Y41FGRka7gfVMeGRkZIT1\n+J0+RyEivYB3gCScgeU5VS1ok+ZcoBBwAA3Af6rqSte2cmBP8zZVPbGD49hzFMYY44eIeY5CVeuA\nSao6EZgA/FRE2lb2/1LVY1xprgL+6rbNAeSo6sSOgkSsK/F1RM0oZecX3ez8TGd8aqNQ1QOuxV44\nryq0g+0AKTiDQzPx9TixKtb/UO38opudn+mMTxW4iMSJyKdANVCsqqs9pDlPRDYALwNXum1SoFhE\nVovI1YEotDHGmNDx9YrC4bqtNBL4gYiM9ZDmBVU9CjgPuMtt06mqeiwwDbheRE4LQLmNMcaEiN+D\nAopIHrBfVR/ykqYUOEFVd7VZnw/s87SviFhLtjHG+CkUjdmddo8VkcE4eyvtEZFkYDJwX5s02apa\n6lo+FkhS1V0i0geIU9UaEekLTAEK8CAUJ2uMMcZ/vjxHMRx4UkTicN6qelpVXxORa3GOhf44cKGI\nXA7UAweBS1z7DgX+6bpaSACKVHV5wM/CGGNM0ETMfBTGGGMiU9i7rYrIVBH5UkS+FpHuz2geRCJS\nLiJrRORTEfnItW6AiCwXka9E5E0R6e+W/r9F5BsR2SAiU9zWHysia13nvMBtfZKILHPts0pERgf5\nfBaKyFYRWeu2LiTnIyI/d6X/ynU1GqrzyxeRLSLyies1NRrPT0RGisjbIrJeRNaJyI2u9THx/Xk4\nvzmu9bHy/fUSkQ9ddck6cbbfRu73F4pp9Dp64QxU/wdkAInAZ8CYcJapk/JuBAa0Wfc74L9cy7cC\n97mWxwKf4rzlluk6z+YruA9xNvYDvAac6VqeDTziWv4PYFmQz+c0nA9Rrg3l+QADgFKgP5DWvByi\n88sHbvaQ9qhoOj9gGDDBtZwCfAWMiZXvz8v5xcT35zpOH9d7PPABcGKkfn/hvqI4EfhGVTepagOw\nDPhZmMvkjaeHB38GPOlafhJn92CAc3F+MY2qWg58A5woIsOAfnroWZTFbvu45/UccEbAz8CNqr4H\nfNdmdTDP53TX8pnAclXdo6q7geVAyy/DQOng/MD5Pbb1M6Lo/FS1WlU/cy3XABtwdl+Pie+vg/Mb\n4doc9d8fdPggc0R+f+EOFCOAzW6ft3DojyESuT88+AvXuqGquhWcf9zAYa71bc+twrVuBM7zbOZ+\nzi37qGoTsFtEBgbjRLw4LIjns8d1Ph3lFSo3iMhnIvJXt0v7qD0/EcnEeeX0AcH9ewz3+X3oWhUT\n3594fpA5Ir+/cAeKaNP24cEf0mY4Ew+fuyMSugzH2vk8AhyuqhNw/gd9MIB5h/z8RCQF56/Fm1y/\nvGPq79HD+cXM96etH2Q+UUTGEaHfX7gDRQXg3mA70rUuIqlqlet9O/ACzltnW0VkKIDrMnCbK3kF\nMMpt9+Zz62h9q31EJB5I1TYPLYZAKM4nbN+7qm5X141a4C84v8NWZW1Tpog9PxFJwFmJ/l1VX3St\njpnvz9P5xdL310xV9wIlOG//ROb3F+gGGj8bc+I51JidhLMx+6hwlslLWfsAKa7lvsBKnA8Q/g64\nVTtufEoCsmjd+NTccCU4G5+mutb/kkONT9MJcmO26ziZwDq3z0E/H1o3pjUvp4Xo/Ia5Lf8nsDRa\nzw/n/eiH2qyLme+vg/OLie8PGIyrARlIxjmVw7RI/f6CWgn5+A82FWePhm+A28JdHi/lzMIZyD4F\n1jWXFRgI/Mt1Dsvd/8GB/3Z9oRuAKW7rj3Pl8Q3wB7f1vYBnXOs/ADKDfE5LgUqgDvgWuML1hxP0\n8wFmudZ/DVwewvNbDKx1fZcv4LwnHHXnB5wKNLn9TX7i+r8Ukr/HMJ5frHx/R7vO6TPX+cx1rY/I\n788euDPGGONVuNsojDHGRDgLFMYYY7yyQGGMMcYrCxTGGGO8skBhjDHGKwsUxhhjvLJAYUw3iMgK\ncc7qaEzMskBhjDHGKwsUJuaISB8RecU1KcxaEblYRPJE5CPX58fc0q4QkYdcIwKvF5HjReQfrgld\nCl1pMlyTxSwRkS9E5BkR6e3huJNF5H0R+V8ReVqcc8YjIveJyOeuEU/vD92/hDGBYYHCxKKpQIWq\nTlTV8cAbwJ9U9UTX5z4icpZb+jpVPQH4H+BFnBO+HA3MEpEBrjTfBx5W1bHAPpzj6LQQkUHAb4Ez\nVPV44GPgZtewzuep6v9T54indwXrpI0JFgsUJhatAyaLyL0icpqq7gPOEJEPxDkt6iRgnFv6l9z2\n+1xVt6lqPc7B0ppH5vxWVT9wLS/BOXueu5NwDty20jXHwOU4R+jcAxx0zZ1wPnAwsKdqTPAlhLsA\nxgSaqn7jamCeBhSKyNvA9cCxqlrpmp/Y/dZRnevd4bYMzrkAOvo/0naQNME5a1hu24QiciLO2Qov\nBm4gyDMXGhNodkVhYo6IDAcOqupSYD5wLM6KfZdrIpyLupDtaBH5gWv5UuDdNts/AE4VkWxXGfqI\nyPdEpC/OEUDfAG4Gxnfh2MaElV1RmFh0NPCAiDiAepxtDucBnwNVwEduab0Nn+y+7Sucsxr+DVgP\nPOaeRlV3iMgs4CkR6eVa/1uc7RkvujV+/2c3zsuYsLBhxo3phIhkAK+o6tHhLosx4WC3nozxjf2i\nMj2WXVEYY4zxyq4ojDHGeGWBwhhjjFcWKIwxxnhlgcIYY4xXFiiMMcZ4ZYHCGGOMV/8fIdWcEkmi\n1jkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3f7e662b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## try for different bin and sample size\n",
    "nbins_list = [30,50,70,100]\n",
    "nsamples_list = np.array([10,50,100,200,300])*1000\n",
    "for nbins in nbins_list:\n",
    "    nh=[]\n",
    "    for nsamples in nsamples_list:\n",
    "        fs =free_energy()\n",
    "        fs.read(\n",
    "            table=multinormal(cov_,nsample_=nsamples,nbins_=nbins),\n",
    "            verbose=False\n",
    "        )\n",
    "        nh.append(fs.H())\n",
    "    plt.plot(nsamples_list,nh,label='bins: %d'%nbins,marker='o')\n",
    "plt.xlabel(\"samples\")\n",
    "plt.axhline(y = gaussian_entropy(cov_), linestyle='--')\n",
    "plt.legend(loc=4)\n",
    "plt.show()\n",
    "        \n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (2), calculate marginal distribution and calculate entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 798,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "free-energy has 1 dimensions\n",
      "binning of free-energy: (50,)\n",
      "bin volume: 0.565828\n",
      "nomalization: 1.000000\n"
     ]
    }
   ],
   "source": [
    "## keep only variable 0\n",
    "f1 = f0.marginalise([0])\n",
    "f1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 800,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVNUbwPHvQUDEDXEHFRB3s7DMNRVU3MV9BZPMtMUl\nl35aRkhUtliW2WKluZGamqm5oojlbq65K5sKruGKGzjn98cMiKSyzcydGc7neeZ5uHfu3Pu+3pnX\nmXfunCOklCiKoigFg53WASiKoijmo4q+oihKAaKKvqIoSgGiir6iKEoBooq+oihKAaKKvqIoSgGS\no6IvhGgvhDgmhDghhBj/iPtrCiG2CSHuCCHGPOJ+OyHEXiHECmMErSiKouRNtkVfCGEHTAfaAXWB\n/kKIWlk2+xcYAXz2mN2MAo7kI05FURTFCHLyTr8hcFJKmSClTAUWAl0zbyClvCyl3AOkZX2wEKIS\n0BH4yQjxKoqiKPmQk6LvDpzJtHzWsC6npgJvAeqnv4qiKBoz6Re5QohOwAUp5X5AGG6KoiiKRuxz\nsE0iUCXTciXDupxoBgQIIToCRYDiQoi5UsoXs24ohFCfBBRFUXJJSpmrN9M5eae/G6gmhPAQQjgC\n/YAnXYWTEYCU8h0pZRUpZVXD46IeVfAzbW+Tt9DQUM1jUPmp/FR+tnfLi2zf6Usp7wshhgPr0f8n\nMVNKeVQIMUx/t/xBCFEe+BsoDuiEEKOAOlLKm3mKysbEx8drHYJJqfysm8qvYMlJewcp5VqgZpZ1\nMzL9fQGonM0+NgOb8xCjoiiKYiTqF7lmEBwcrHUIJqXys24qv4JF5LUvZGxCCGkpsSiKolgDIQTS\nBF/kKvkUHR2tdQgmpfKzbpnz8/T0RAihbhZ28/T0NNr5zlFPX1GUgiEhISHPV4UopiOE8X7ipNo7\niqJkMLQLtA5DyeJx50W1dxRFUZQnUkXfDApST9gWqfwUW6KKvqIoVuHEiRPUr1+fkiVLMn36dLMf\n/+LFi7Ro0YKSJUvy1ltvmf34xqJ6+oqiZLDknv6QIUMoWbIkn3/+uSbH/+CDD9i/fz9Lliwx+7GN\n2dNXV+8oipKtuLgEQkJmk5iow93djvDwYLy8PMy6j4SEBPr37//Y+3U6HXZ2pmteJCQkUKdOHZPt\n32y0HjAo08BB0lZt2rRJ6xBMSuVn3TLn96jXYWxsvPT2HivhpgQp4ab09h4rY2Pjc3yM/O6jVatW\nslChQtLJyUkWL15cnjx5UgYHB8vXXntNduzYURYrVkxu3LhRXrt2TQ4cOFCWLVtWenp6yg8++CBj\nH7Nnz5bNmjWTo0ePli4uLtLb21tu27ZNzp49W1auXFmWL19ezpkz55HHDw4Olg4ODtLR0VEWL15c\nbty4Uep0Ojl58mTp7e0ty5QpI/v27SuvXLkipZTyzp07MigoSJYuXVq6uLjIhg0byosXL+b43yur\nx9VHw/rc1drcPsBUN1X0rZfKz7plV/QDAydlKtYyo2gHBk7K8TGMsQ9fX185c+bMjOXg4GDp4uIi\nt2/fLqXUF9qBAwfKbt26yZSUFBkfHy9r1KghZ82aJaXUF30HBwc5Z84cqdPp5LvvviurVKkihw8f\nLu/duyfXr18vixcvLlNSUh55/ODgYBkSEpKx/OWXX8omTZrIpKQkee/ePfnqq6/K/v37SymlnDFj\nhgwICJB37tyROp1O7t27V964cSPHuWalir6iKCbxqNehr+97WYp1+u1x63O+rZ/fezmO7VFFf9Cg\nQRnL9+/fl46OjvLYsWMZ62bMmCH9/PyklPqiX6NGjYz7/vnnH2lnZycvXbqUsa506dLywIEDjzx+\n1qJfu3ZtGRUVlbGclJQkHRwc5P379+WsWbNks2bN5MGDB3Oc35MYs+irq3cURXkid3c7ICXL2hQC\nA+1yXPIDAx+9Dze3/JWgypUfDO57+fJl0tLSqFLlwZxPHh4eJCY+mPOpfPnyGX8XKVIEgDJlyjy0\n7ubNnI0In5CQQPfu3XF1dcXV1ZU6derg4ODAhQsXGDhwIO3ataNfv35UqlSJCRMmcP/+/TznaUyq\n6JuBrV8Hbc35xcUlEBQUhp9fKEFBYcTFJfxnG2vOLyeyyy88PBhv71AeFO0UvL1DCQ8PzvExjLGP\nR8k8PEGZMmVwcHAgIeHBOUxISMDdPTdTeudclSpVWLNmDcnJySQnJ3PlyhVSUlKoWLEi9vb2hISE\ncPjwYbZt28bKlSuZO3euSeLILXX1jlJgxcUl4O//NTExYUBRIIUdO0KJjBwBkHGlib19Ah4eXrm+\nWsVWeHl5EBk5gpCQKSQl6XBzsyM8fESu/j2MsY/s2NnZ0adPHyZOnMicOXP4999/mTp1Kv/73/8e\n+xh9hyRvhg0bxjvvvMOcOXOoUqUKly5dYvv27QQEBBAdHU2ZMmWoU6cOxYoVw8HBwaRXFuWGKvpm\n4Ovrq3UIJmWt+YWEzDYUfGcQ90EWJSYmjNGj3+XQoUIP/Wfg76//z8AWC39Ozp+Xlwfz54fm6zj5\n3UfWQcceNQjZtGnTGDFiBFWrVqVIkSIMHTqUl156KV/7fNx9o0aNAqBt27acO3eOcuXK0bdvXwIC\nAjh//jyvvvoqiYmJFCtWjH79+jFw4MBsczQH9eMspcDy8wslOnoSBLwCPj9DSnm47o64mYi82omB\nMVfodzaBxbdeZxGd6RH4bb4Ln6Wz5B9nFWRqwDUrU9B7wpbK3d0Onv8S3HfBJ8nw4y5Y/RnP73fj\nz7+iGHF+K4ubxVHDOZzT1Kb75jVw7JjWYRudtZ4/JW9U0VcKrDrt6oHvRFg0H+6WpMR1O2bdHMu6\nmCNE3BxJ4+tnmb1nO+90vkYD7xDsSzqCry/4+cHixaDTaZ2CouSaau8oBdLB2PM8+93zvFYljCs7\nTlN730Fej4/Evksnkt/6H637/fKgp19hGwxswzs1ZvLhiz3h99/h/fchOBjGjdM6FaNS7R3LpMbe\nUZR8uHUnlRbT+tLMZTBfjxgMiRPgUAxsXA9NmlAciIws/dCVJl41Z/HhiZHU2VKTwD59oGFD/c3X\nFxo00DolRckx9U7fDKKjo632CpecsIb8Mg/2daDCekRZJy5M3Yj9vr3QqRMcOgRlyz7ysen5vT3n\ndz49/DqzWkQQufBPqu/7h2FnNnFn65941qtr5oyMJ/P5U+/0LZPZv8gVQrQXQhwTQpwQQox/xP01\nhRDbhBB3hBBjMq2vJISIEkIcFkL8I4QYmZvgFMUY0q/Hj4gYR/Tl2lwpd4Gi6+pwJiYWhgyBKVMe\nW/AzmzyoG91Kvknwxm5ErBjApCNL+ONGF/Y27fPIH3UpiiXK9p2+EMIOOAG0BpKA3UA/KeWxTNuU\nATyAbsAVKeUXhvUVgApSyv1CiGLAHqBr5sdm2ofNvtNXtBUUFEZExDgoFweD/GBeJJyvzsJnutK3\ngj2sWQM5nHg6KCiMiHg7eGYRzNyG891C7KE+UU3r8frWpSbOxPTUO33LZO53+g2Bk1LKBCllKrAQ\n6Jp5AynlZSnlHiAty/rzUsr9hr9vAkcB0/wmWlEeIzFRB6II9AiE9VPgvA/VSaT9kW3w/fc5LvgZ\n+9r6Lpx/Bpp8zi2K0p9FDNi9FmJjTZiFohhHToq+O3Am0/JZ8lC4hRCegA+wM7ePtXa2fh20pefn\n7m4HNX8FnQMceBGBjh94mRVPtwBPz2wfnzk//eBjtyDqA2g4HZwvsZ8arHy6OfTvD6mpJsvDVCz9\n/KXTerpEY/Py8iIqKsrsxzXL1TuG1s4SYJThHf8jBQcH42l4Ebq4uODj45PxBVP6E1Mtq+XcLoeF\nDSJi9LOwZjggGMJ0TjkcpujoB4Ujp/sLDw9mx45QYmJaw+YW0Px9vE8Vxn5MMNHTvsT3vfdg8mSL\nyj83y5bs008/pVWrVuzbt0/rUDSRfo6io6OJj4/P+46yG3sZaAyszbQ8ARj/mG1DgTFZ1tkDa9EX\nfDWevmJ2Yb+slo5v1pT9B7wnezUdLa8VdpZnVq/N8/5iY+NlYOAk2cB3jOR/ReTSDdv0d1y4IKWb\nm5SZxli3Npb8OmzTps1D4+lndf/+fTNG82RpaWnZbuPp6Sk3btyYo/097rxgiklUgELAKfRf1DoC\n+4Haj9k2FBibZd1c4IscHCdHyStKbty/r5PFRjWRI2YskFKnk7JbNynffddo+2/y7tuyxriXH6xY\nuFDKZs2Mtn9zs9TXodbTJUqpn8Tl7bfflg0bNpQlSpSQ3bp1y5geMT4+Xgoh5MyZM2WVKlVky5Yt\npZRSLl++XNatW1eWKlVK+vn5yaNHj2bsz9PTU06ePFnWqVNHurq6ysGDB8u7d+8+8thmLfr6/dIe\nOA6cBCYY1g0Dhhr+Lo++738VSAZOA8WAZsB9w38U+4C9QPvHHOOx/9jWriBNt2dpPv8tSjqMqSHv\n3kuTcskSKWvVkvLOnVzt40n5xSYlSzG+jFy9yzBbU2qqlJ6eUu7YkY+ozSu76RIthdbTJfr6+spK\nlSrJI0eOyFu3bsmePXvKoKAgKeWDoj9o0CB569YteefOHXnixAlZtGhRuXHjRpmWliY//fRTWa1a\nNZmamiql1Bf9evXqycTERHnlyhXZrFmzh2bmyszsRd8cN0t+suWXJRdFY7Dk/EqNaiWHTJ8tZVqa\nlF5eUkZH53of2eXXNvwjWXl0nwcrvvxSyt69c30creSm6DMJo9zyQuvpEtPf6ac7cuSIdHR0lDqd\nTsbHx0s7OzsZHx+fcX94eLjs27dvxrJOp5Pu7u5y8+bNUkp90f/hhx8y7l+9erWsVq3aI49tzKKv\nhmEwA0v/tWp+WWp+M1Zv44Z9LF+9PABWr9b/AKtly1zvJ7v85r0xkoqTq7Mgeh/9fevD4MEQHg5x\nceDllcfozSc350+GWtY1/OaeLjHz8Tw8PEhNTeXy5csZ6ypVqpTxd1JSEh4eD+ZfEEJQuXLlh+LJ\nvL2HhwdJSUlPyNY41Cibis0KifyAPu4TcHZygOnTYfhwkxynXKmi9Cj7DqOWT9SvKF4cXn4ZvvrK\nJMdTHjD3dIlnzjy4ej0hIQFHR8eH/tPIHI+bm9tDsaQ/PnOhz7o/Nzc3o8X6OKrom4E1XA6XH5aY\n3/yNe0i2/4fvhgbD8eOwfz/07p2nfeUkv5+HDyW50D806DkYP79QRp5M5f7s2XDlSp6OaU6WeP7y\nIvN0iTdv3iQhIYGpU6c+ccYqfYck5+bPn8+xY8e4desWoaGh9O7dO6PQZ91Xnz59WLVqFZs2bSIt\nLY0pU6bg5OREkyZNMrb55ptvSExMJDk5mY8++oh+/frlKp68UEVfsUn/++MDupZ9ixJFC8O33+rH\n2HFyMtnxLp0/h/Oup9hT8hjR0ZP4elk4K++7kfzxJyY7ZkGT0+kSnZ2dqVq1Ki1atCAoKMho0yUC\nDBw4kEGDBuHm5sa9e/f4KtOnuayPrVGjBvPnz2f48OGULVuWVatWsXLlSuzt7TO2HzBgAG3btqVa\ntWpUr16diRMnPvH4xqBG2VRsRvpImvvPnePws/P5e+DfPOflDh4ecOAAZOrHGltQUBgRv4yG15pA\n5GdwsiPPsJ2oIu1wvXoZHB1NdmxjUmPvPJ6fnx8DBw5k8ODBZj+2mi5RUbLIPJLm4VLXYPsE+vb4\nictffqWf6cqEBR8MY/LIEhAdBi3DAMkBmhBfuDQsWmTSYytKbqiibwa20jN9HEvILyRktn6mK5dL\n4B0Ju0cREzOJO59Py/cXuDnJTz8mTwoc7QGFb+hjIIWNPk3h88/Bgt89W8L5swbZtX6shSr6ik1I\nTNQBRcFnNvwzAO4Vx5fdpKahn93KxMLDg/H2DgV5G/6cCC3C8PYOpddPH8C9e6DBwFqKcUVFRWnS\n2jE2VfTNwFKvYzcWS8jP3d0OxA3w+Rn26V+YI5jKpjoNczV08qPkJD8vLw8iI0cQGDiFFmWOIUoc\nY8Dbz+Pl7QVjxujf7VsoSzh/ivmoL3IVmxAXl0CDPuNIfu4EzDhAZY5xwK4+1w7swfOpOmaP56Vp\ns1gRu4B/v4yEO3f0QzhHRUEd88eSG+qLXMukvsi1MrbeM7WE/Ly8PCjSNJXKl93x8wtlet1XsHux\nv1EKfl7y+2boQK7Zn+SntTv0l4q+/jp88UW+YzEFSzh/ivmooq/YhLhzV0h0jmbjl/OIWv02AReP\nU/LtCZrF4+zkQO+K45m47gP9ildfhSVL4NYtzWJSFFDtHcVG9P/8O7YmRnP6i0Uwdy788gusXatp\nTFdv3qH0+97Mbb+SwFbPQrt2+nF5+vbVNK4nUe0dy6TaO4qSxYqzsxjWyHBlhQnH2ckNl2JOBJR+\ni//98aF+Rf/++v+MFEVDquibga33TLXOb8lfB7lrf563erSBgwfhwgXo0MFo+89Pfj++OpTzDltZ\ntvUQdO8O0dEWNx6P1ucvp7SaU9bWqKKvWL2PVv9M4yKDcHQoBL/+qm+fFCqkdVgAlCnpzAv2wQTN\n6I9fty/YVdKdSzN+1DosJZP79+9rHYJ55XYAflPdsOFJVBTTuXHrrhT/Kys37D2pnw6xenUpd+/W\nOqwMsbHx0qP6CMlbpSWlj8uezJNbilSRsbHx2T9YA5b6Ohw4cKC0s7OTzs7Osnjx4vKzzz6TO3bs\nkE2bNpUuLi7Sx8dHRmeaICcpKUkGBARIV1dXWb16dfnjjz9m3Ddp0iTZq1cvGRQUJEuWLPnEeXct\nxePOC2rmLKWgeWvWUllyVEv9wr59+tmxdDpNY8osMHCShJuSlmGSrsHSiVsyGRf5eo+xWof2SJb8\nOvT09JRRhknnExMTZenSpeXatfoJ7jds2CBLly4tL1++LKWUsnnz5hnTIO7fv1+WLVs2Y4awSZMm\nSUdHR7lixQoppX6aRUtnzKKv2jtmYC0907zSMr85B2fRp7rhC9xff4U+ffL9C9ys8pNfxvAQu1+H\n2su44yBZSRdq/XPIaPHlV67yE8I4tzyShitY5s+fT6dOnWjXrh0ArVu3pkGDBqxevZqzZ8+yfft2\nPvnkExwcHHjmmWcYMmQIc+fOzdhPkyZN6NKlCwCFCxfOczzWSBV9xWrtPZnEJadtfDCgp35As/Si\nb0EyBmK7VQaSGoD3ehbQnfZXDmsdWt7o2wP5v+VTQkICv/76K66urri6ulKqVCm2bt3KuXPnSEpK\nwtXVFWdn54zts06bWNnEo65aMjVHrhnY+tgm5s4vfdz8NTf+pFhxL1KuXobYy/p3kPXrG/14+ckv\nPDyYHTtC9SOAHu0BtRcRd68CVa/espg5dK3l+Zl5lMvKlSvz4osvMmPGjP9sd/bsWZKTk0lJSaFo\n0aIAnD59+qFpE21lxMy8UO/0FavyYNz8sSR7nOHG5s/w9/+aqz/8aJLWTn5lHoitQbFTUON3Fq98\nnUJ9+sDChVqHZ1UqVKhAbGwsAEFBQaxcuZL169ej0+m4c+cOmzdvJikpiUqVKtG0aVPefvtt7t69\ny8GDB5k5c+YTp00sUHL7JYCpbljwF0j5lf4Fkq0yZ34ZX4xW+UvyRi0JOgk35IWiLlLu32+SYxoz\nv+KjmsoPFq6VcvNmKevVM9p+8yNzfpb8Oly+fLmsUqWKLFWqlPz888/lrl27ZMuWLaWrq6ssV66c\n7Ny5szxz5oyUUv9Fb+fOnaWrq6usVq2a/OGHHzL2M2nSJDlw4ECt0siTx50X8vBFbo7aO0KI9sCX\n6D8ZzJRSfpLl/prAz8CzwDtSyi9y+lhFyY2ML0afWggHgwBBA45yV1cInn5a6/Cy1bJ8D+bvWcrE\nj7/X/0jr0CF46imtw7IKAQEBBAQEPLTucV9Cu7m5sXLlykfeFxoaauzQrEq27R0hhB0wHWgH1AX6\nCyFqZdnsX2AE8FkeHmvzrKVnmlfmzC9j3Pxay+BILwD6EMEe77oma+0YM7+3OvXguFjOvftS/yOy\nBQuMtu+8svXnp/KwnPT0GwInpZQJUspUYCHQNfMGUsrLUso9QFpuH6souREeHkwZn2FwpwT8WxO4\nyQD7WTT49B2tQ8uRFk974XTPne9Wb9GPxbNwoUVPpajYnpwUfXfgTKbls4Z1OZGfx9oMdZ2+8Xh5\neeDWqjhlL5fDzy+U99q9SVmPclRq39ZkxzR2fk1L9eDnHb/Bs8/qh4vYtcuo+88tW39+Kg+zqEs2\ng4OD8fT0BMDFxQUfH5+Mj57pT0y1XLCXW7RoyREiebvJO7TyqYbvihXQOIjozZstIr6cLL/Zridd\nv2jBxqiutO7fHxYsIPr2bYuIT7Fc6ecoOjqa+Pj4PO8n2/H0hRCNgUlSyvaG5QnovzH+zxeyQohQ\n4Eb6F7m5fKzMLhZFWbR5P0Ere3D30xjskODhoR83v25drUPLlcJjavN9+zm8VKUEtGoFZ85YxCBx\najx9y2Tu8fR3A9WEEB5CCEegH7DiSfHl47GK8kTTo5biU7gndnYCduyAkiWtruADNCjWgxl//ga1\nakHZsrBtm9YhKQVEtkVfSnkfGA6sBw4DC6WUR4UQw4QQQwGEEOWFEGeA0cBEIcRpIUSxxz3WVMlY\nKlv/6GzO/Hbf+I1XXuihXzDTsAumyO91357svb0UnU5C167wmMsLzSFzfh4eHggh1M3Cbh4eHkY7\n3znq6Usp1wI1s6ybkenvC8AjB7N41GMVJS9W7zpGqv1VBrdtBDodLF4MGzZoHVae9Petz6DVqSzb\ndoieAQEQFASffqp1WPnqFVuq6OhodVlqJmqOXMVqtA3/iPMpSRz8eDr89Re88YZ+piwr9dw7Yyju\nWILo996DypUhKgpqqvdHSs6ZqqevKBZhS/JSghv11C8sWwY9e2obUD4NadaTndd/Azs76NJF0xaP\nUnCoom8Gqqeff1sOxXOn8Gle79Rcv+KPP/SF0gxMld8r7Ztwz/4SkXtOQkAArNDmGgf1/CxYVNFX\nrMKUP5ZR/X5XnBzt4fhxuHXLJMMom5N9ITuqprZm0Cev0/bjrdzasYuEPfu0Dkuxcaqnr1iFEm++\nwNiG7xA6oCNMmQKnTsH332sdVr7ExSXQoM84kn3i4afdLCWAbeVSeWPH93h5Ge9qDcV2qZ6+YpP2\nx5zjZpHDjO7aWr9i5UqztXZMKSRkNsl7fwTXWCieyAp60PhiYUJCZmsdmmLDVNE3A1vvKZo6v09+\n/50q9zpSomhhSE6Gffv0v2I1E1Pll5ioA50LxLWCqhtYRSf82cSlM/dMcrzHUc/PgkUVfcXirT/z\nG32eMlyps2YN+PpCkSKaxmQMGfPnxraBqhu5TFn+oS6tCyVoHZpiw1TRNwNb/2GIqfKLi0ugS9/x\nJBf5i9j1e4iLS9CktWOq/MLDg/H2DoW4JlB1A3CTra72vOpm3ikf1fOzYFFf5CoWKX0u3Jji1aHG\nevh1LjWrvsvhy7ModOwYVKyodYhGEReXwMR3f2ZB+S9pcTqIiFe7Uyl4kH4ANgub71exPOqLXAtl\n6z1FU+QXEjKbmJgwqLUajnUHilIx1p8E+2JmL/imPH9eXh78EjGJ6oV6Uq5RLSq1bgXOzvrvLcxE\nPT8LFlX0FYuUmKiDQvbgGQ0nOwDQhQ1sLlFd28BMxN+7NVuSNujf3Wv4Qy3F9qmibwa23lM0RX7u\n7nZQZQNcrg23SwOSLqzgVG3zT7FsjvP3WrvWXHDazJ17aWYv+ur5WbCooq9YpPDwYArXDYNTfgDU\nZB/FC11gyPQJ2gZmIk95lcfxbiUWRO+Fpk3h9Gl9X19RjEwVfTOw9Z6iKfLz8vIA7xs0dL2En18o\nofUn4NwnAK+qnkY/VnbMdf5qF27Ngp0bwN4eOnUy2wBs6vlZsKiir1ikvSeTuFfkApsXfE9UVBj9\ni92hROAArcMyqS512/D3vxv1C6qvr5iIumRTsUhDps9mXexqznzxq/5XuJ6ecOGCTfwo63HOXrpO\n5S/d+Hf8JVxFGri7w9mzUKKE1qEpFkpdsqnYjPWxa2nj2V6/sGYN+PnZdMEHqFS2BMVSnuan9Vuh\neHFo1gzWrdM6LMXGqKJvBrbeUzR2fvdS73PWMZLX27XVr9B4gDVznr/6Jduw7IChxdO5M6xebfJj\nqudnwaKKvmJx5kX9jeM9N56vWQlSU/Xvdjt10joss+j5bGsO3jTM+9uxo/5Tjk6nbVCKTVE9fcXi\n+IWFcfPeTXZ/+Jl+3tjx42H3bq3DMovrKXcp+WFZYkck4FWxFNSpA3PnQoMGWoemWCDV01dswt9X\n19K7fjv9go2MnZ9TJYoWpvTtpny3dpN+RceOZmnxKAWHKvpmYOs9RWPmF5OUzM0ihxna/gX9ilWr\n9L1tDZn7/DUq25o1xwx9/Y4d9f8GJqSenwVLjoq+EKK9EOKYEOKEEGL8Y7aZJoQ4KYTYL4TwybR+\ntBDikBDioBAiQgjhaKzgFdszffUGyt1ugUsxJzh5Em7etPq5cHMrqGkbTqQZiv4LL+jnBL54Udug\nFJuRbdEXQtgB04F2QF2gvxCiVpZtOgDeUsrqwDDge8N6N2AE8KyU8mnAHuhn1AysgK2P/WHM/FYd\nX0dzN0NrZ80a6NBB8yGGzX3+ejd/hlSHy+w+fhYcHaF1a1i71mTHU8/PgiUn7/QbAiellAlSylRg\nIdA1yzZdgbkAUsqdQEkhRHnDfYWAokIIe8AZSDJK5IrN0ekkMWItw1oZrs9ftUrf3ihg7AvZ4X6v\nFT9syNTiUX19xUhyUvTdgcwjP501rHvSNomAu5QyCfgcOG1Yd1VKuSHv4VonW+8pGiu/ZdsOYacr\nTOv61SAlBbZtA39/o+w7P7Q4fy0qtyYqzlD0O3SA9eshLc0kx1LPz4LF3pQ7F0K4oP8U4AFcA5YI\nIQZIKX951PbBwcF4enoC4OLigo+PT8ZHs/QTp5Ztd3lKxCJqurbHzk4Q/dVXUK0avoYhCCwhPnMu\ne969zS9nFuPr50kl90K8WLIkjt9+i+/IkRYRn1rWZjn97/j4ePIq2+v0hRCNgUlSyvaG5QmAlFJ+\nkmmb74FNUspFhuVjQEugOdBOSvmKYf1AoJGUcvgjjqOu0y/gXN9sw2vPjeTDgQHw2mvg7Q3jxmkd\nltnFxSWA6CYoAAAgAElEQVTQxn8asV2WwPy1cLkK011aE9jvWVy++1br8BQLYqrr9HcD1YQQHoYr\nb/oBWYf/WwG8aAiiMfo2zgX0bZ3GQggnIYQAWgNHcxOgUjBcvJLClaI7eaOjH0hZYPv5oJ8qMjbm\nfYj1B+/1QFEirn7ItYW/aR2aYgOyLfpSyvvAcGA9cBhYKKU8KoQYJoQYathmNRAnhDgFzABeN6zf\nBSwB9gEHAAH8YIpELFnmj2a2KD/5xcUlEBQURsO+Q3C4VIa715Ph8GEoVAhq1zZekPlg7vOXmKgD\nisLxAKi9DICd+OKSct0kE6uo52fBkqOevpRyLVAzy7oZWZb/07IxrA8DwvIaoGK74uIS8Pf/Wj8B\neocJcGQw/v5fs6uXA64dO2p+qaZW3N3tgBSIaQvdX4SiF9ClFOOQezWarVkDQ4dqHaJixdTYO4pm\ngoLCiIgYBxSFEdVh8WI4X50j5epSe9Y3BWaQtawe+s+w5xCIb4L31dNse6My5TZvgt9/1zpExUKo\nsXcUq5LRxnCJg8I34MLTlCCNKv+e04+fX0B5eXkQGTmCwMApVL55H8f6XxEZOYJyLwbBpk1w967W\nISpWTBV9M7D1nmJe88toY3hHQmwbkHa0ZQUny1cBZ2ejxpgfWpw/Ly8P5s8P5e8FP3Ov7GXSHIpB\n6dJQty78+adRj6WenwWLKvqKZsLDg/H2DoWqayDGH0ihb/FPqDR0oNahWYxypYrifsefj35brl9h\nhgHYFNumevqKpo6fiKXWzKd4bvcwalcswc/rv8V+506oWlXr0CzGiBkLWHRsPhenroJ9+6BvXzhx\nQuuwFAugevqK1dl5OhnHVE/+jprKvDEB2Jcpowp+FuN7dOJSkS2cvngNfHz0I4+ePKl1WIqVUkXf\nDGy9p5if/CJ2RFK7sGF8HQv9QZbW569S2RKUv92SyUtX6i9j7dDBqAOwaZ2fqdl6frmlir6iqV2X\nIwmoayj6q1dbZNG3BF2q9WT5iaX6hU6dVF9fyTPV01c0c/naLcp+Wo7EMedw092B6tX1k4U4qnl2\nsoo7d4WqX3twblwSFex14O4O585BsWJah6ZoSPX0FasyY+1flEipj1vp4rBuHbRqpQr+Y3hVLEXp\n20355LfVUKIENG4MGwrcKOWKEaiibwa23lPMa37LDkbyXCnL7ueD5Zy/Tp69WHrE0OLp3Bn++MMo\n+7WU/EzF1vPLLVX0Fc0cvhVJv4b++slB1q/Xf0GpPNaE7l0547SW5Ou3H/T1dTqtw1KsjOrpK5o4\nFHeBej/U5HbYZZx2boeRI/XXoCtP5Ppma4bVH8HkQd2gVi2YPx8aNNA6LEUjqqevWI3v1m2kwh1f\nnBztYeVK6NJF65CsQtvKPVl4cIl+oXNndRWPkmuq6JuBrfcU85Lf+lORNHc39PP/+MOii74lnb+B\nDZ4n3mEpzX0n8uGBy9xdmv+JVSwpP1Ow9fxySxV9xex0Okmc2MDLvv4QEwNXrsBzz2kdlsWLi0tg\n1MuL4GJ9tiQ1Y9KGL7lz+Dind+7WOjTFiqievmJ2q3cdI2BxW+59koDd19Pgn3/gp5+0DsviZcw/\n0PhHKH8Qls9iIT253iiNV3Ys1zo8RQOqp69YhVmbI6mKP3Z2Qt/P79xZ65CswkPTKFZbA0j+oBu1\nYtTga0rOqaJvBrbeU8xtfluSImlf3R+uX4edO6FNG9MEZiSWcv4y5h+4UhVSnaHcIdbQgueuxeVr\nYhVLyc9UbD2/3FJFXzGrW3dSuVBkM6+1b63/Fe4LL6ihBHIoY/4BUiDWH7xX4eL9NXb16sLmzVqH\np1gJ1dNXzOq7VVsZu2E4t6bug0GDoFEjeP11rcOyGnFxCYSEzObPy4e56L6Po+9uwGtBBJw/D9Om\naR2eYmZ56emroq+Yle+kSdxKvcWu9ydDhQqwZw9UqaJ1WFYnfQC2axMvUSLmOHTrpr8SSuTq9a9Y\nOfVFroWy9Z5iTvKLi0sgKCiMLUkzST12nqRly8HNzSoKviWeP6+KpSh6uzY/rd8G9epBaiocO5an\nfVlifsZk6/nlVo6KvhCivRDimBDihBBi/GO2mSaEOCmE2C+E8Mm0vqQQYrEQ4qgQ4rAQopGxgles\nQ1xcAv7+XxOxZCj3y1xl/4ovWPHKJ1xt3kLr0KyaT3F/luyN1L+7N+IAbIpty7boCyHsgOlAO6Au\n0F8IUSvLNh0AbylldWAY8H2mu78CVkspawPPAEeNFLvV8PX11ToEk8ouv5CQ2cTEhIHXTkhsBGll\naHb1JlNP3TZPgPlkqeev93P+HLgZqV/IR9G31PyMxdbzy62cvNNvCJyUUiZIKVOBhUDXLNt0BeYC\nSCl3AiWFEOWFECWA5lLKnw33pUkprxsvfMUaZFxfXmcxHO2BB/GU4zJb7rppHZpVe7ltE24VOc7J\ns//q5yLYt0//62ZFeYKcFH134Eym5bOGdU/aJtGwzgu4LIT4WQixVwjxgxCiSH4Ctka23lPMLj93\ndztwuAQ1VsGRXnRhJatpS0X3QuYJMJ8s9fwVK+JI2dvN+W5dFBQpAi1a6C+DzSVLzc9YbD2/3LI3\nw/6fBd6QUv4thPgSmACEPmrj4OBgPD09AXBxccHHxyfjo1n6iVPL1rccHh7M7/sCSNlTFVLK0Znl\nfFvqNv06dSKdJcVrTcvNKviz5ngkAdFloUYNfP/4A/r1s5j41LJxl9P/jo+PJ6+yvWRTCNEYmCSl\nbG9YngBIKeUnmbb5HtgkpVxkWD4GtDTcvV1KWdWw/gVgvJTyP0Mqqks2bVv5kR0omuRI7Yu1WbJ9\nKhf27sWzXl2tw7J6K3YcoeeSTtz9NBa7xLPg4wMXLoC9qd/PKZbAVJds7gaqCSE8hBCOQD9gRZZt\nVgAvGoJoDFyVUl6QUl4Azgghahi2aw0cyU2AivU7e+k6F4ttY+PXs1n1ZkOKtPJVBd9IOjesjc7u\nHpsOxEDlyuDhAVu3ah2WYsGyLfpSyvvAcGA9cBhYKKU8KoQYJoQYathmNRAnhDgFzAAy/8RyJBAh\nhNiP/uqdj4ycg8XL/NHMFmWX3weLl1P+dku8KpayygHWLPn82dkJvHT+/BRluIqne3f4LXdj7Fty\nfsZg6/nlVo4+A0op1wI1s6ybkWV5+GMeewB4Pq8BKtbvt5ML6FtroH4+19Wr4b33tA7JprSt5s/v\nx5cCr0GPHvq5hr/8Uv06V3kkNQyDYlLHz1ym1nfenBuXSIUTh2DIEDh0SOuwbMqhuAs8/UMtboVd\nwsmhkH7u3IgINXduAaCGYVAsTvjS36h8pz0VXIvBsmXQNetPPJT8esqrPIXvVGF+1N/6d/c9euS6\nxaMUHKrom4Gt9xSflN+q+IX0f7ofSAlLl0KvXuYLzEis4fzVdfZnwa689fWtIb/8sPX8cksVfcVk\n9sec41qRfbzdqwMcOKDv6fv4ZP9AJde6P+3P38mGot+gAaSkwNECN+KJkgOqp6+YTM9Pp7H/4h5i\npsyBkBD97E6ffqp1WDbp8rVblP20PIljknArXRxGjoTy5WHiRK1DU0xI9fQVi7I+aSGDnuunX1i6\nFHr21DYgG1ampDMlbjxNu6HD8PML5cMj17m7cKHWYSkWSBV9M7D1nuKj8ttyKJ4UpxOM694GjhyB\nGzfgeeu8ctcazl9cXAKpx+w4dMuF6OgwQjd+RcqRGM78tSXbx1pDfvlh6/nllir6ikl8tPxXaup6\n4uzk8OBdvp16uplKSMhsbh+aDN7RANynJMt0vYke/YG2gSkWR70KzSB90CRblTm/9Bmy1iVOxTnm\nLnFxCVbf2rGG85eYqIPzTaHIv1AqBoBl9OaZU9l/mWsN+eWHreeXW6roK0aTMUPWui7onAV7l01j\niO/7pCUmQtOmWodn09zd7UDehsN94ekIADbQBO9bF/QDsCmKgSr6ZmDrPcX0/DJmyKq7Eg73AVmC\nBqc9iS7lCYWsY+z8R7GG8xceHoy3dygc6AXPzAVuUsl7MrKdP6zIOj7iw6whv/yw9fxySxV9xWj0\nM2Q5wzPz4J8BAPRkJWucaz35gUq+eXl5EBk5gv6+UQh5lYY93iQycgTFBg5Uv85VHqKu01eMJigo\njIi/mkKXUfDNYapwmr95jrH9X2fuL+9rHV6B4R/+IeduJnLok2/1V025u8OZM1CypNahKUamrtNX\nNBUeHozD82NgXyAg6MFCNhWvSNiHL2sdWoEyqUcQR8SvXE+5C8WLQ8uWsGqV1mEpFkIVfTOw9Z5i\nen5FXcqQWv00Hdyv4OcXyqtlv6HFV2Px8vLQNsB8srbz16yuByXuPMWHi1frV2QzAJu15Zdbtp5f\nbqmirxjNewt+o+ydF1j96xSi5g+jZtpNKgQO0DqsAqlb1YHMOzBPv9ClC0RGwu3b2galWATV01eM\nxvXN1gTXe5UvXu4N33wDO3fC3Llah1Ugnb54DY+pHpwaEYu3myu0agUjRuhH4FRshurpK5rZciie\nq04HmNQvQL9iyRKr/kGWtatSriSV77bn3QWL9Cv69oVfftE2KMUiqKJvBrbeU4yOjua9pXN5in6U\nKFoYLl2CffugbVutQzMKaz1/g58byB9nDC2ePn1g/Xq4evU/21lrfjll6/nllir6Sr7dv6/jrxuz\nGdcmWL/i99+hfXsoUkTTuAq6//VsS0rhGDbuOwWlSkGbNvpPYEqBpoq+Gdj62B+HbhaikM6ZoFbP\n6VcsWmSVM2Q9jrWeP2cnB56x68/7yw3v9oOCYP78/2xnrfnllK3nl1uq6Cv5Nv2v2fiXfQk7OwGn\nT8P+/dC5s9ZhKcDYNgPZdnM+Op2Ejh31k9InJGgdlqIhVfTNwJZ7iueTb3IqaTEf9g3Ur5g3T98/\ndnLSNjAjsubzN8DvWex0hZmxZhsULgy9e0NExEPbWHN+OWHr+eVWjoq+EKK9EOKYEOKEEGL8Y7aZ\nJoQ4KYTYL4TwyXKfnRBirxDiySM/KVYnZMFSXO7V4+mqFfSTn8+ZA4MGaR2WYmBnJ2jk1IUJv4zC\nzy+UsJh73Js1S3+ulAIp2+v0hRB2wAmgNZAE7Ab6SSmPZdqmAzBcStlJCNEI+EpK2TjT/aOB54AS\nUsqAxxxHXadvhVze9GXoMyP59KUesG0bDB6sn5Bb5OrSYcVE4uISaN75QxI7L4UvEiEtjQT7ShRa\nOh/3ANWCs3amuk6/IXBSSpkgpUwFFgJds2zTFZgLIKXcCZQUQpQ3BFUJ6Aj8lJvAFMsWF5dAm96j\nuOa4m/j1+/STpcyZA8HBquBbkJCQ2SQemQoXnoHqq4BizE57lYPjP9Q6NEUjOSn67sCZTMtnDeue\ntE1ipm2mAm8BBfZtvK31FNMnS9l4uSgcepnFC5vSufUX3P/1V/0VIjbGms+ffrjrovBPf6j7KwDz\neJlGsf9AWhpg3fnlhK3nl1v2pty5EKITcEFKuV8I4Qs88S1gcHAwnp6eALi4uODj45NxuVX6iVPL\n2i/rJ0tpBT6DYdtq4Drl44ow29WZlytV0jw+tfxg2d3dDkiBw+Wg5h/gkMKpVDfWOBXB/fPP8R0/\n3qLiVctPXk7/Oz4+nrzKSU+/MTBJStnesDwBkFLKTzJt8z2wSUq5yLB8DGgJjAKCgDSgCFAc+E1K\n+eIjjqN6+lbCzy+U6Hg/6DASvjsACFbTgT11ivDuYTVhhyVJ/1QWExMGQT1gbxDedw+wY2AJypw4\n/p8reRTrYqqe/m6gmhDCQwjhCPQDsl6FswJ40RBEY+CqlPKClPIdKWUVKWVVw+OiHlXwFevi7m4H\nPj/B/mBAUJEkGrGT2Hq1tQ5NySJ9Rq3AwCmU/7cwRZ7/mMjIEZR543X9GPs3bmgdomJm2RZ9KeV9\nYDiwHjgMLJRSHhVCDBNCDDVssxqIE0KcAmYAr5swZquT+aOZLRg5rjfUWgIHuwHQghAii1ciZPJQ\njSMzDWs/f15eHsyfH8pfP/zM7YqJFHUpA2XKQIsWsGyZ1eeXHVvPL7dy1NOXUq4FamZZNyPL8vBs\n9rEZ2JzbABXL8+PW3ZROeYH23eaRlHifnruX0uynn3Cz8slSbF31SqUpc6cJk5euYuqQPvov3X/8\nESZO1Do0xYzUePpKrpV6sxVDnn6Dzwb3hN27oX9/OHlSXappBV6aNovI+NWc/WKJflIVd3f90Axu\nblqHpuSBGk9fMbkth+K55nSQkL6GH/ak/wJXFXyr8G6vbiQ6RXI++aZ+FNTu3dU4+wWMKvpmYEs9\nxdDf5j0YN//uXVi4kOjq1bUOy6Rs6fx5u7lS9nYzPlryh35FcDDR06aBTqdtYCZkS+fPGFTRV3JM\np5P8eW0OY1sbxtb54w+oVw8qVNA2MCVXArz7sOSo/odavPACODjAhg3aBqWYjerpKzn27R9bGBM1\njFtTDumHUQ4I0E+JqAZYsypx565Q9WtPEseexa10cZg5E5Yt0/8nrlgV1dNXTOqrzbNpXSb4wbj5\nW7aoeXCtkFfFUpS73ZzJ6S2eAQP0k9ifPKltYIpZqKJvBrbQU7x87RYn7ZfyYV/D2Dpff60fXK1Y\nMZvI70lsMb+u1fqw9Ji+xRO9cycMGQLffKNxVKZhi+cvP1TRV3LkvQXLKHOnCT7eFfW/4pw1C0aN\n0josJY8m9grgXJEozl66rl/x2mswd676hW4BoHr6So6UftOfoLpD+OqVvvDVV/qx8xct0josJR/K\nj+5Cz5r9+PZVw6xnvXtDy5Yw/Im/s1QsiOrpKyax8+gZrjjtJax/V/1wvF9+CWPGaB2Wkk/dq/dh\n2fHFD1aMHKlv29nw5ZuKKvpmYa09xbi4BIKCwugwYRDFz1TlyqUL8Pvv+l9xNmqUsZ215pdTtprf\nwEY+nHdaS12fAQQFhRHnXhmcnWH9eq1DMypbPX95pYq+8kjpQ/JGRIzlSpVErv/5Kf7+X3Nn8mT1\nLt8GxMUlMKjvHIhvw5EbXkREjMO/7XQu9R8A06ZpHZ5iQqqnrzxSUFAYERHjoMo+CBgC04/SmE0s\nL9aTclcvQ6FCWoeo5EPG+a23HJ79CeZsBG4R3G8yP2/8QX85bo0aWoepZEP19BWjyZhmr9knsHMk\nIBjDdyyq2EgVfBuQcX4P94biSVB9DVCUhAuF4JVXYPp0rUNUTEQVfTOwxp6iu7sdVNgGFffCvsF4\nEocfUeyv/+x/trXG/HLDFvPLmEZR5wC/D4K2Y8HuKm5udvrLN+fPh+vXtQ7TKGzx/OWHKvrKI4WH\nB1PI7yXYNhLSnBjFFJa6ePPux8O0Dk0xgvDwYLy9Q4EUONsYblSgaItehIcHQ6VK4O8Ps2drGqNi\nGqqnrzzSih1H6LasJT1OD+Xu2TQW7JpG8qZNVGnaWOvQFCOJi0sgJGQ2SUk67rleZJv3EuLGnMSj\nvAts3ar/xfWxY6qdZ8Hy0tNXRV95JK+xQdQoVZd1774NU6bA/v36j/yKzar11isUd3Rh94efgZTQ\nvDkMHQovqmmtLZX6ItdCWVtPceO+UyQ4rGXmsDcgNVV/Cd/o0Y/d3tryy62Ckt+vr4Wz5/7PRO2P\n0U+K8/HH8N57+nkTrJitn7/cUkVf+Y/Xf/mY5k5vUKlsCf27e29veO45rcNSTOzpqhVoXXQ0L86d\noF/xwgvw1FPw/ffaBqYYlWrvKA/ZfuQ0zeb5cPyNk1Qv5QQ1a8KSJdBY9fILguTrtykXXpNpLX/h\n9c4vwMGD0Latftjl4sW1Dk/JQrV3lHwbNvdTni/0CtUrlYbPP9e/21MFv8BwLVGEod6TGRc5nAGB\nofiNWsqWIuW58m6I1qEpRqKKvhlYS09xf8w5DolfmDVkLJw7px9Nc/LkbB9nLfnlVUHLb3Sbpty5\nc4EF/3gQHR3GwPhfYPoMEnbv0SbAfLL185dbOSr6Qoj2QohjQogTQojxj9lmmhDipBBivxDCx7Cu\nkhAiSghxWAjxjxBipDGDV4wjfWC1FuN74JJQA2d5G0JC4OWXwctL6/AUMwubNBe5Zh60mQSON4mn\nLvN0wRwd+IbWoSnGIKV84g39fwynAA/AAdgP1MqyTQdgleHvRsAOw98VAB/D38WA41kfm2kfUjG/\n2Nh46e09VuIcJxlfSlLiuOxY6UWZVrqMlFeuaB2eogFf3/ckSEnXYEmH4RKkLMd5edW+iJRxcVqH\np2RiqJvZ1vHMt5y8028InJRSJkgpU4GFQNcs23QF5hoq906gpBCivJTyvJRyv2H9TeAo4J77/5oU\nUwkJmU1MTBg0/Q4O94HrNRh5NpH5Xs+Di4vW4SkayBiiYd0XUPs38PiTixQjqtbzEBqqdXhKPuWk\n6LsDZzItn+W/hTvrNolZtxFCeAI+wM7cBmntLLmnmJiog/KnoP7PsPk92rEWT84wv+h/x9h5HEvO\nzxgKWn4ZQzTccYRV30LAYDyqvUP9iG9g7Vr45x9N4swrWz9/uWVvjoMIIYoBS4BRhnf8jxQcHIyn\npycALi4u+Pj44OvrCzw4cWrZuMvlK0qoEQy/vYTdjaNMYRz/431wWEN0dLTm8all8y97eXkQHt6A\nWbPe4P59D3ZeK4LdszHEJ1/Gc8IEmDiRaMOcCpYQb0FaTv87Pj6evMr2On0hRGNgkpSyvWF5Avo+\n0ieZtvke2CSlXGRYPga0lFJeEELYA38Aa6SUXz3hODK7WBTja/HOeLYmLUI35xCv8AsDmMeQqg2J\n3DASLy8PrcNTLMDxM5ep83U9ZrRaxhBfH/1vN+bO1c+nq2jKVNfp7waqCSE8hBCOQD9gRZZtVgAv\nGoJoDFyVUl4w3DcLOPKkgq9oY83u42zRzWTBgPm83OcjJjuOZWP72qrgKw+pWbkMI2t8xRvrB3M1\nDf2wHIMHw40bWoem5EG2RV9KeR8YDqwHDgMLpZRHhRDDhBBDDdusBuKEEKeAGcBrAEKIZkAg0EoI\nsU8IsVcI0d5EuViszB/NLEXafR19I16mR+lQ+rR9gZ8871O6TzfC1/yQ64JvifkZk8oPPh/cm9Ky\nFq3fH0fQ4v2suVWCqGf9iItLMH2A+WTr5y+3ctTTl1KuBWpmWTcjy/LwRzxuK6DGZbVA/b74BoCF\nY96A6GiYNw/27tU2KMVi2dkJZvZ4m44rfNkbtZ4V50dzkHqENx3Nu9umqk+GVkSNvVOApI+ffjgp\nmf3P/8isZn/wUtOnoX59+PFHaF/gPoQpuRAUFEbEP+7Q9Gv4cTct728jggGE9xrI94s/yX4HitGp\nsXeUx4qLS8Df/2siIsayv9IR2DaBD0evIqV/f+jbVxV8JVuJiTo4+DJcqwKt32EzviyiH/02/6Ef\nf1+xCqrom4El9BQzfoTlsxiKJMO2d2gX686FvYfho4/ytW9LyM+UVH56+h9t3YLfZ+t/tPXUQiby\nDtXvXoSFC00ZYr7Y+vnLLVX0C4jERB2UiwP//8GKmdTTHWESkwmt3gMcHbUOT7ECGT/auu0Ei36D\nDsMp8fRYmPczjBoFZ89qHaKSA6qnX0B06DOWte6LYeNknP/pxm6e52NGowtMYv589dN6JWcyz6t7\nttQh4jz3cGLcPrx+/Bq2bIF16/SzbilmoebIVR4p6d8beIc3QRwtxe31a5nBaJy5zqSq7uqafCVf\n6r/9Jol3TnD2o2U4+raAwEAYqQbTNRf1Ra6F0rKneOdeGs980AdPx6b88+08vnrhJTo7LWFzHy+j\nFXxb75mq/B5va+hn3JO3aDbpLcZWaMS/b03ky5b9LOr6fVs/f7mlir4N0+kkPhNfA2DfB9/inXyJ\nkcc24bZ5LT8umqze4Sv55uzkwMJeU/g7dTZfnGiK/73NDPhzIyHN3rKowq88oNo7Nia955qYqONY\n2S1cdbtATMh23OJPQMeO8NNP0KWL1mEqNiQoKIyIaF8Y0Bt+/pMWly+wmF5Mb9+d99f8oHV4Nk21\ndwq4B9fijyP63xqcdz9J6XUvIHf8pS/433+vCr5idImJOkhsCZGfQFA7/ixTniHM5M2oX+DoUa3D\nU7JQRd8MzNVTzLgWv+YGaD8aItbgeuwlivbqA998A927m+S4tt4zVfk9WcakK/tfgugwCPZlpZsL\nvz7nD+3awenTRokzr2z9/OWWKvo2JDFRBw3mQufXIGIVdS9J1tGNaVXbQ69eWoen2KiM6/dJgf3B\nsPIrCGxPQnBvGDMG/P3h4kWNo1TSqZ6+jdDpJOUGtODfcudh/jpqXbnLRlozjg+xCzytrsVXTCrz\n9ftubnZU9avFBydGMLh8GJ1nLaBuwkm+9u/P6M9HqQsIjEhdp1+AZP7CtmwFHdvLHSaZeEqva0zZ\n4wNYQV/eZhLbvI8TGTlCvdAUs/vq11W8uasvbPqI8XtvMYov+V/F5ry/dYp6PhqJ+iLXQhm7p/jQ\nF7bbx7DYcQvn/j3Iny/NZ1+PUkQVbsePdZuiC0wyS8G39Z6pyi9vdq/4G37eCs2/5JMWaQQxn0/P\n/cWu7i+DTmeSYz6KrZ+/3DLLHLmKcWV8YVviCgzoBGeaUXbVDBw3dqB01Qpw7DCTDHMNK4pWEhN1\nkPwMzNoCfbsT5bGZ539fzrq4XtCtG8yZA6VKaR1mgaPaO1bI1+89Nt+oBK3ehW3j6Li1Dj/xCus9\nazDo5EawV/+XK9oLCgojImIcUBTs0qDZp9D4C2qeasrRSp6kLV9OaN0ObL9dHnd3O8LDg1XbJ5fy\n0t5R1cHCZe7du7vb0bxnU7bWmA33K1J49mo+uTSPbnxDH+bg0Ww7g1TBVyxEeHgwO3aE6j+V6orC\nX6Mod+0Y8X5/U1Z3iz5p1QlbswR7RvIFw9ix4zP1/ZMZqJ6+GeS1p/hQ737LBCLOpfDq9q74Fe/J\nhOUlOHqpF24k4cNWznmvJzw82Jhh55it90xVfnnj5eVBZOQIAgOn4OcXSmDgFHb8Hs7lyfuxu5jK\ndwMO0qj6e9TkOCfxoUNMBSa985PR47D185db6m2hBdP37ieB9xboOBzHxHoMnjae8EJTcX7uaT5z\naiPJB9kAAAnWSURBVM3mtEp0cptJeLh6h6RYHi8vj0deLlz3gi/Rf00hrlswQX5FeGbzUD46EYnP\n77u4OLU0Y/++wtkkVNvHBFRP3wJkbeGEhwdzQxbG780gkj3OUDhVMOQ3X8ZfWsVBnmb1sxX4Zs9M\nrcNWlDzL6PcLJ6i2Dhp8B5W30WGXK+/vvIbD7fJ8xLusoDXu3h+rts9jqOv0LdyjijuAv//X+r6n\nnSNU/41CDSYi3C7SYU85OuysQ4+UneykEeGE8Dd1CAycon5spVi19NZlTEwYUBRIwa32aG7XPcgV\nj1gC9pZn+G4dDW4mspzOnG9lR98ZkwiZNO+h109B/48gL0UfKWW2N6A9cAw4AYx/zDbTgJPAfsAn\nN481bCdtRWxsvAwMnCR9fd+TgYGT5C+/LJCxsfHS23ushJtSP4v0TentPVa26TxcUmeepMsQ6Tiq\nvOzgX1vOdG0irzkWkbefeUZ+7NpcVuXgQ4+JjY3XOsWHbNq0SesQTErlZxrprxM/P/3rJDY2Xvr6\nvicpdFdSd6Gkdy9Z4Y2S8s2mFeXfLs7ynIOTnCqGykZsl/Zceei1kPU1l/k1Ysvnz1A3c1TH0285\nKfh2wCnAA3AwFPVaWbbpAKwy/N0I2JHTx2bahxn+iXLvcU+mJ63PWtxLl/aVAQHj9OvsUiXuOyQt\nw2SF/vVlQE87GV63qlxbsqZMFiXkXzSTo5gq+zYZ9dBxMr8wLM3UqVO1DsGkVH7mExg4KdNrR0pE\nmsR9k3RqW1vW6OEjJzVzkAddisrrhRxllFM9ubhOU3nup1nyOc/X//OGKjY2XsbGxstnn22X49ev\ntclL0c/JF7kNgZNSygTDx4mFQFfDu/d0XYG5hsq9UwhRUghRHvDKwWOf6FEtkfSPdI+7L7frH7cv\n4D8fQXfsCGXWrO4MHrzsP+sjI0fov3yNC4HiV8F1D4Vcj5B2+Qo3Umbzkt98vFIv8VSiM8/vkBRJ\nFRxwKM3Wu734hv+3d/axVZV3HP98zy20lpKCpVQ3XgY4UBdRyKYG3dQ5FZ2ubo4No0zUOCIjG9HE\niZspYyZTs/iHbmRza6x/7MUtCxOY2UBxdWQxc0G3ERHFibyKSmSOMuD23u/+uAe4lLb0tvf29LbP\nJznpeZ5znpPvL09/v/Oc5+1exNe4gHdpANq4afIPga4HwgYS+/btS1pCSQn29R/HTfNkBPggU6pW\nU//hdby45iGWVv6HpRNfoPbc57ig8jlmHXiRUUv/zrp3zK7KFWyMzmRLdhpb3ppK87y7+fOOOja8\nPQP4Hj3x397Ej2LGqH7hZG8F4Abg8bz0zcCjHe5ZBczKS68FZvakbN61E95iXXWJHHmDd3attXV9\nQfldPWvS5Lt89XWLTeVOc8oeV9Vuct3YdZ4w7gl/6pxz/Mmp9/iSs+/0tdPn+JbzPu27Zk71gzNG\n+KfTK/y7aXLruGF+q6bShxR5carSf6kZ6Zaqa9yk+/xVfuVJvGn4rxsbF3epq1xoampKWkJJCfb1\nL5193Z7wBRD7ypy53/W0y29z6sKHPf2Kz3nOZ6d4yaxaN58buXWcvLM65fuEX6+p8frRp3lF3cfd\nfPqpfqD+8/5W3Z2+qfb7/sLIH/uy6hbfffk8b376GV888Q6PYrurOGDxYVHjysnKFAq9aOmfdCBX\n0g3AVba/HqdvBs63/c28e1YBP7D91zj9LHAPuZZ+t2XznuGOWo6t6DuFP3ykIXcfJooOIEM2W82R\nEQxhRBbRhqg+dq9NRJaIg4jhRIYIk7KpyGZJOZM7t0gZhmXN8IwZnuHYkYWDKWirSNFWUUFbKkOb\nammLqtivavZSx/vZBqJTdxCNGs0LL89jr8ezjQlsZzxpbqWxcSwbN6aOa11MmZJrXQDH7VBYbgNU\n8+fPp6WlJWkZJSPYlzydDfwe8Z/77285tvL3KPupG3c7ew8uJfW/RUw5/VrGVmyjPtpFQ6qVMZl6\n6rP7GNO+n5pMmpHpNDXtaWrSZuRhqDkMVe1QlYF0BIdScCgVcSiKSEeiPYJ2iUzUTruGkZFojyKy\nQFbCSmNBxlVkEZYw5A61YY/EiMeihazetQxo69UEjZLM3pF0IbDU9uw4fS+5t8tDeff8BHje9lNx\n+jXgEnJBv9uyec8Y3FN3AoFAoAQUGvR70qf/EnCGpInAbmAucGOHe1YC3wCeil8S+2zvkfR+D8r2\nSnggEAgECuekQd92RtIiYA252TjNtjdJWpC77MdtPyPpGklbgDbg1u7KlsyaQCAQCHTLgFmcFQgE\nAoHSk+iGa5K+LGmjpIykmXn5EyUdkLQhPpYnqbO3dGVffG2JpDckbZJ0ZVIai4WkJkk78upsdtKa\n+oqk2ZJek/S6pG8nrafYSNoq6R+SXpb0t6T19BVJzZL2SPpnXt5oSWskbZb0J0m1SWrsC13YV7Df\nJb3L5r+ALwKtnVzbYntmfCzsZ13FolP7JJ0FfAU4i9zCtuWSBsOYxiN5dfbHpMX0BUkR8CPgKuAT\nwI2SzkxWVdHJApfanmH7/KTFFIEnyNVXPvcCz9qeBqwDlvS7quLRmX1QoN8lGvRtb7b9BtBZwCv7\nINiNfY3Ar223295KbvuKweB0ZV9neRxdlGg7DRxZWDiYEMk3/IqG7fXABx2yG4En4/Mngev7VVQR\n6cI+KNDvBnKFfyz+XHle0sVJiykyHwW256V3xnnlziJJr0j6eTl/Rsd0rKMdDI46ysfAWkkvSboj\naTElYqztPQC23wHGJqynFBTkdyXfT1/SWqAhP4vcP9t3bK/qotguYILtD+K+8N9LOtv2/hLLLZhe\n2leWdGcrsBxYZtuSHgAeAW7vf5WBArjI9m5J9eSC/6a4NTmYGWwzVwr2u5IHfdtX9KJMmvgzxvYG\nSW8CU4ENRZbXZ3pjH7mW/fi89Lg4b0BTgK0/I7c1RzmzE5iQly6LOioE27vjv+9JWkGuS2uwBf09\nkhridUOnAe8mLaiY2H4vL9kjvxtI3TtH+6UkjYkH0pA0GTgD+HdSwopEfr/bSmCupOGSJpGzr6xn\nT8QOdYQvARuT0lIkji5KlDSc3MLClQlrKhqSqiXVxOcjgCsp/zqDnJ919LX58fktwNP9LajIHGdf\nb/wu0Z9LlHQ98BgwBlgt6RXbVwOfAZZJOkxuhsEC2wNnK8Ae0pV9tl+V9BvgVSANLDxh46Hy42FJ\n55Grr63AgmTl9I0hsLCwAVgRb39SAfzC9pqENfUJSb8ELgXqJG0DmoAHgd9Kug14m9ysubKkC/su\nK9TvwuKsQCAQGEIMpO6dQCAQCJSYEPQDgUBgCBGCfiAQCAwhQtAPBAKBIUQI+oFAIDCECEE/EAgE\nhhAh6AcCgcAQIgT9QCAQGEL8H+lwZ4RCTyWZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3f7c2a3cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(f1.edges[0],np.exp(f1.fes),marker='o', label='from fes')\n",
    "plt.plot(f1.edges[0],f1.prob,label='from prob')\n",
    "\n",
    "theo_prob = np.exp(-(f1.edges[0])**2/(2*10))/sum(np.exp(-(f1.edges[0])**2/(2*10)))/f1.volume\n",
    "plt.plot(f1.edges[0], theo_prob,label='teor')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 801,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theo sum:  2.5701360853\n",
      "theor:     2.5702310797\n",
      "sum:       2.57056375488\n",
      "method:    2.57056375488\n"
     ]
    }
   ],
   "source": [
    "#plt.plot(f1.edges[0], theo_prob*np.log(theo_prob))\n",
    "#plt.plot(f1.edges[0], f1.prob*np.log(f1.prob))\n",
    "#plt.show()\n",
    "\n",
    "print(\"theo sum: \",sum(-theo_prob*np.log(theo_prob))*f1.volume)\n",
    "print(\"theor:    \", gaussian_entropy(cov_[0,0]))\n",
    "print(\"sum:      \", sum(-f1.prob*np.log(f1.prob))*f1.volume)\n",
    "print(\"method:   \",f1.H())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 802,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0029385319146779665"
      ]
     },
     "execution_count": 802,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## check Mutual information \n",
    "MI(f0,[0],[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (3) start from 3D gaussian and marginalise to 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 803,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cov2_ = np.array([[10,2,1],[2,5,5],[1,2,4]])\n",
    "tab1 = multinormal(cov2_,nsample_=500000,nbins_=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 804,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "free-energy has 3 dimensions\n",
      "binning of free-energy: (50, 50, 50)\n",
      "bin volume: 0.110069\n",
      "nomalization: 1.000000\n",
      "numer:  6.5076798134\n",
      "theor:  6.52311534619\n"
     ]
    }
   ],
   "source": [
    "f3 = free_energy()\n",
    "f3.read(table=tab1)\n",
    "\n",
    "print(\"numer: \" ,f3.H())\n",
    "print(\"theor: \", gaussian_entropy(cov2_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 805,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numer:  4.01692136252\n",
      "theor:  3.98916961291\n"
     ]
    }
   ],
   "source": [
    "f3_1d= f3.marginalise([1,2])\n",
    "print(\"numer: \", f3_1d.H())\n",
    "print(\"theor: \",gaussian_entropy(cov2_[1:3,1:3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (4) calculate conditional entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 859,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7725900870187348"
      ]
     },
     "execution_count": 859,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_entropy(f0,[0],[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEST (5) three-body co-information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 860,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "free-energy has 3 dimensions\n",
      "binning of free-energy: (50, 50, 50)\n",
      "bin volume: 0.110069\n",
      "nomalization: 1.000000\n"
     ]
    }
   ],
   "source": [
    "f3.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 861,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.017592154669062054"
      ]
     },
     "execution_count": 861,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CI(f3,[0],[1],[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:mylab3]",
   "language": "python",
   "name": "conda-env-mylab3-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}