Metric - Manifold Alighnment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.spatial.distance as dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Loading + Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PBMC immune cells from Seurat/CCA alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dat = pd.read_csv('immune_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "meta = pd.read_csv('immune_meta.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dat1 = meta['x'] == 'IMMUNE_CTRL'\n",
    "dat2 = meta['x'] == 'IMMUNE_STIM'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dat_array = dat.values[:,1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14053, 14039)\n"
     ]
    }
   ],
   "source": [
    "print(dat_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14053, 6573)\n"
     ]
    }
   ],
   "source": [
    "ctrl = dat_array[:,dat1.values]\n",
    "print(ctrl.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14053, 7466)\n"
     ]
    }
   ],
   "source": [
    "stim = dat_array[:,dat2.values]\n",
    "print(stim.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>AAACATACATTTCC.1</th>\n",
       "      <th>AAACATACCAGAAA.1</th>\n",
       "      <th>AAACATACCTCGCT.1</th>\n",
       "      <th>AAACATACCTGGTA.1</th>\n",
       "      <th>AAACATACGATGAA.1</th>\n",
       "      <th>AAACATACGGCATT.1</th>\n",
       "      <th>AAACATACTGCGTA.1</th>\n",
       "      <th>AAACATACTGCTGA.1</th>\n",
       "      <th>AAACATTGAGTGTC.1</th>\n",
       "      <th>...</th>\n",
       "      <th>TTTGACTGCCCTAC.1</th>\n",
       "      <th>TTTGACTGCTACCC.1</th>\n",
       "      <th>TTTGACTGCTTGGA.1</th>\n",
       "      <th>TTTGACTGGCGAAG.1</th>\n",
       "      <th>TTTGACTGGGTACT.1</th>\n",
       "      <th>TTTGCATGAACGAA.1</th>\n",
       "      <th>TTTGCATGACGTAC.1</th>\n",
       "      <th>TTTGCATGCCTGTC.1</th>\n",
       "      <th>TTTGCATGCTAAGC.1</th>\n",
       "      <th>TTTGCATGGGACGA.1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AL627309.1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>RP11-206L10.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LINC00115</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NOC2L</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>KLHL17</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 14040 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Unnamed: 0  AAACATACATTTCC.1  AAACATACCAGAAA.1  AAACATACCTCGCT.1  \\\n",
       "0     AL627309.1               0.0               0.0               0.0   \n",
       "1  RP11-206L10.2               0.0               0.0               0.0   \n",
       "2      LINC00115               0.0               0.0               0.0   \n",
       "3          NOC2L               0.0               0.0               0.0   \n",
       "4         KLHL17               0.0               0.0               0.0   \n",
       "\n",
       "   AAACATACCTGGTA.1  AAACATACGATGAA.1  AAACATACGGCATT.1  AAACATACTGCGTA.1  \\\n",
       "0               0.0               0.0               0.0               0.0   \n",
       "1               0.0               0.0               0.0               0.0   \n",
       "2               0.0               0.0               0.0               0.0   \n",
       "3               0.0               0.0               0.0               0.0   \n",
       "4               0.0               0.0               0.0               0.0   \n",
       "\n",
       "   AAACATACTGCTGA.1  AAACATTGAGTGTC.1        ...         TTTGACTGCCCTAC.1  \\\n",
       "0               0.0               0.0        ...                      0.0   \n",
       "1               0.0               0.0        ...                      0.0   \n",
       "2               0.0               0.0        ...                      0.0   \n",
       "3               0.0               0.0        ...                      0.0   \n",
       "4               0.0               0.0        ...                      0.0   \n",
       "\n",
       "   TTTGACTGCTACCC.1  TTTGACTGCTTGGA.1  TTTGACTGGCGAAG.1  TTTGACTGGGTACT.1  \\\n",
       "0               0.0               0.0               0.0               0.0   \n",
       "1               0.0               0.0               0.0               0.0   \n",
       "2               0.0               0.0               0.0               0.0   \n",
       "3               0.0               0.0               0.0               0.0   \n",
       "4               0.0               0.0               0.0               0.0   \n",
       "\n",
       "   TTTGCATGAACGAA.1  TTTGCATGACGTAC.1  TTTGCATGCCTGTC.1  TTTGCATGCTAAGC.1  \\\n",
       "0               0.0               0.0               0.0               0.0   \n",
       "1               0.0               0.0               0.0               0.0   \n",
       "2               0.0               0.0               0.0               0.0   \n",
       "3               0.0               0.0               0.0               0.0   \n",
       "4               0.0               0.0               0.0               0.0   \n",
       "\n",
       "   TTTGCATGGGACGA.1  \n",
       "0               0.0  \n",
       "1               0.0  \n",
       "2               0.0  \n",
       "3               0.0  \n",
       "4               0.0  \n",
       "\n",
       "[5 rows x 14040 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat_pre = dat = pd.read_csv('immune_pre.csv')\n",
    "dat_pre.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14053, 6573)\n"
     ]
    }
   ],
   "source": [
    "dat_pre_array = dat_pre.values[:,1:]\n",
    "ctrl_pre = dat_pre_array[:,dat1.values]\n",
    "stim_pre = dat_pre_array[:,dat2.values]\n",
    "print(ctrl_pre.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Topography Preservation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def markov_matrix(dat,length):\n",
    "    x = dist.cdist(dat,dat)/length\n",
    "    y = np.exp(x) - np.eye(len(x))\n",
    "    row_sum = np.sum(x,axis = 0)\n",
    "    z = np.divide(x,row_sum)\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def KNN(dat,K):\n",
    "    x = dist.cdist(dat,dat)\n",
    "    z = np.zeros((dat.shape[0],K))\n",
    "    #print(z)\n",
    "    for i in range(0,len(dat)):\n",
    "        y = np.argsort(x[i,:])\n",
    "        #print(i)\n",
    "        #print(y[1:K+1])\n",
    "        z[i] = y[1:K+1]\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def markov_RMSE(dat1,dat2,len1,len2):\n",
    "    a = markov_matrix(dat1,len1)\n",
    "    b = markov_matrix(dat2,len2)\n",
    "    deviation = a - b\n",
    "    RMSE = np.mean(np.linalg.norm(deviation,axis = 0))\n",
    "    MSTD = np.std(np.linalg.norm(deviation,axis = 0))\n",
    "    return RMSE, MSTD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def jaccard_score(dat1,dat2,K):\n",
    "    a = KNN(dat1,K)\n",
    "    b = KNN(dat2,K)\n",
    "    jaccard = np.zeros(len(dat1))\n",
    "    for i in range(0,len(a)):\n",
    "        #print(a[i])\n",
    "        #print(b[i])\n",
    "        intersect = len(np.intersect1d(a[i],b[i]))\n",
    "        union = len(np.union1d(a[i],b[i]))\n",
    "        jaccard[i] = float(intersect)/float(union)\n",
    "        #print(jaccard)\n",
    "    return jaccard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "kBET in R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cross Entropy - Fang Wang and Ken Chen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "658\n"
     ]
    }
   ],
   "source": [
    "# Subsample data\n",
    "sample = range(0,ctrl.shape[1],10)\n",
    "print(len(sample))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### KNN Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "5\n",
      "10\n",
      "25\n",
      "50\n"
     ]
    }
   ],
   "source": [
    "ks = [3,5,10,25,50]\n",
    "jacs = np.zeros(len(ks))\n",
    "jacs_var = np.zeros(len(ks))\n",
    "for i in range(0,len(ks)):\n",
    "    print(ks[i])\n",
    "    #print(np.mean(jaccard_score(ctrl_pre.T[sample],ctrl.T[sample],ks[i])))\n",
    "    jacs[i] = np.mean(jaccard_score(ctrl_pre.T[sample],ctrl.T[sample],ks[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEYCAYAAAANjbKIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHrFJREFUeJzt3XtwVOX9x/HPbgLEJSGwGyQTCFoSKhev4ab5IZIG0hmn\naEZRLqEdQJ1aoQiOloDcvKC2QhFKuAiIg61D1HoBWzuhKokQIkFMARnUgAS5NSQLAdIkkJz9/UHZ\nskKSA+7ZnJD3a4Zh9+zJOd98BT+cc548j8Pn8/kEAIBNOZu6AAAAGkJQAQBsjaACANgaQQUAsDWC\nCgBgawQVAMDWwkN1oqKiIq1evVqGYSg1NVXp6ekBn3/11Vf6wx/+oGuvvVaSNGDAAA0fPjxU5QEA\nbCokQWUYhlatWqUZM2bI4/Fo2rRp6tu3r7p06RKwX8+ePZWZmRmKkgAAzURIbv0VFxcrNjZWnTp1\nUnh4uJKTk1VYWBiKUwMAmrmQXFF5vV55PB7/e4/Ho2+//fai/b7++ms9+eSTcrvd+uUvf6n4+PhL\nHu/w4cOW1drcxMXF0Y8G0J/G0aPG0aOGBas/cXFxl9wesmdUjfnJT36ipUuXKiIiQtu3b9fLL7+s\nRYsWXXLf+r6Zlop+NIz+NI4eNY4eNczK/oQkqNxut8rLy/3vy8vL5Xa7A/ZxuVz+10lJSVq1apVO\nnjypdu3aXXQ8/mXzP/xLr2H0p3H0qHH0qGFWX1GF5BlVQkKCjhw5otLSUtXW1io/P199+/YN2OfE\niRM6Pz9ucXGxDMNQVFRUKMoDANhYSK6owsLCNH78eM2dO1eGYSglJUXx8fHKycmRJKWlpamgoEA5\nOTkKCwtT69atNXnyZDkcjlCUBwCwMUdzXOaDS/D/4ZZEw+hP4+hR4+hRw66KW38AAFwpggoAYGsE\nFQDA1ggqAICttcigMta92dQlAABMapFB5Vu/tqlLAACY1CKDCgDQfBBUAABbI6gAALZGUAEAbI2g\nAgDYGkEFALA1ggoAYGsEFQDA1ggqAICtEVQAAFsjqAAAtkZQAQBsjaACANgaQQUAsDWCCgBgawQV\nAMDWCCoAgK0RVAAAWyOoAAC2RlABAGyNoAIA2BpBBQCwNYIKAGBrBBUAwNYIKgCArRFUAABbI6gA\nALZGUAEAbI2gAgDYGkEFALA1ggoAYGsEFQDA1ggqAICtEVQAAFsLWVAVFRXp8ccf129/+1u9//77\n9e5XXFyskSNHqqCgIFSlAQBsLCRBZRiGVq1apenTp2vBggXavHmzDh48eMn9/vKXv+iWW24JRVkA\ngGYgJEFVXFys2NhYderUSeHh4UpOTlZhYeFF+3300UcaMGCA2rVrF4qyAADNwGUFVVlZmb755pvL\nPonX65XH4/G/93g88nq9F+2zdetWpaWlXfbxAQBXr3AzO5WVlWnhwoXav3+/JOmNN95QQUGBioqK\n9OijjwalkNdff10ZGRlyOhvPzri4uB91ru+DcAw7uZq+FyvQn8bRo8bRo4ZZ2R9TQfXqq6/qtttu\n0zPPPKOHHnpIknTzzTdrzZo1pk7idrtVXl7uf19eXi632x2wz969e7Vw4UJJ0smTJ/Xll1/K6XSq\nf//+Fx3v8OHDps7bkGAcww7i4uKumu/FCvSncfSocfSoYcHqT31hZyqoiouLlZmZGXC143K59J//\n/MfUyRMSEnTkyBGVlpbK7XYrPz9fkyZNCtgnKysr4HWfPn0uGVIAgJbFVFBFR0fr6NGjAWl38OBB\nxcTEmDpJWFiYxo8fr7lz58owDKWkpCg+Pl45OTmSxHMpAEC9TAXVsGHD9Pvf/17p6ekyDEObNm3S\ne++9p/T0dNMnSkpKUlJSUsC2+gJqwoQJpo8LALi6mQqqn/3sZ4qKitI///lPeTwe5eXlacSIEdya\nAwBYrtGgMgxDb7/9tu677z7169cvFDUBAODX6Fhwp9OpnJwchYWFhaIeAAACmPqB30GDBmnDhg1W\n1wIAwEVMD0//xz/+oXXr1snj8cjhcPg/e+aZZywrDgAAU0GVmpqq1NRUq2sBAOAipoJq8ODBFpcB\nAMClmQoqSfr000+Vl5cnr9crt9utQYMGKSUlxcraAAAwF1TvvvuucnNzNWzYMMXExKisrEzr1q3T\n8ePHdd9991ldIwCgBTMVVB9//LHmzJmjjh07+rfdcsstmj17NkEFALCUqeHpNTU1Fy1mGBUVpTNn\nzlhSFAAA55kKqltvvVWLFi3S4cOHdebMGR06dEiLFy9myXgAgOVM3fobP368XnvtNT355JOqq6tT\neHi47rjjDo0bN87q+gAALZypoHK5XJo4caIee+wxnTp1SlFRUaZW4gUA4McylTa5ubkqKSmR0+lU\ndHS0nE6n9u/fr7y8PKvrAwC0cKaCKjs7Wx6PJ2BbTEyM1q5da0lRAACcZyqoqqqq5HK5Ara5XC5V\nVlZaUhQAAOeZCqouXbqooKAgYNvWrVvVpUsXS4oCAOA8U4MpMjIy9OKLLyo/P1+xsbE6evSodu7c\nqWnTplldHwCghTMVVD169ND8+fO1adMmlZWVKTExUWPHjlVMTIzV9QEAWjjTk9LGxMQoPT1dknT6\n9GlFRkZaVhQAAOc1GFS5ubmKjo7WrbfeKknat2+fXn75ZXm9XsXGxmrq1KmKi4sLSaEAgJapwcEU\n69evV/v27f3vly1bpptuuknz5s3TTTfdpDfeeMPyAgEA9lbxl+WWHr/BoCovL1fXrl0lSWVlZfr+\n++/1q1/9SvHx8crIyFBxcbGlxQEA7O/kmyssPX6DQeV0OlVbWytJ+uabbxQXF+d/NtWmTRtmTwcA\nWK7BoOrVq5fWrl2rkpISffTRR+rTp4//s0OHDgXcFgQAwAoNBtW4ceP03XffaebMmWrTpo1/1J8k\n5eXlscwHAMByDY76c7vdmj179iU/y8jIsKQgAAAuxFodAABbI6gAALZGUAEAbI2gugzGujebugQA\naHHqHUyxa9cuUwe48cYbg1aM3fnWr5XuGd3UZQBAi1JvUC1dujTgvdfrlcPhUFRUlE6dOiWfzyeP\nx6PFixdbXiQAoOWqN6iysrL8r999912dPn1aI0aMUJs2bVRTU6Ps7GxFRUWFpEgAQMtl6hnV3/72\nN40ePVpt2rSRdG76pNGjR+vDDz+0tDgAAEwFVURExEUT0O7du9cfXAAAWMXUwokjRozQCy+8oD59\n+sjj8ai8vFzbt2/XQw89ZHV9AIAWzlRQDRo0SN26dVNBQYGOHz+uzp076/7771eXLl2srg8A0MI1\nGlSGYejZZ5/V008/reHDh1/xiYqKirR69WoZhqHU1NSACW4lqbCwUNnZ2XI4HAoLC9PYsWPVo0eP\nKz4fAODq0GhQOZ1OlZaWyufzXfFJDMPQqlWrNGPGDHk8Hk2bNk19+/YNuCK76aab1LdvXzkcDpWU\nlGjBggV65ZVXrvic9fGd8J77veK4HNEdgn58AEBwmRpMMXz4cK1YsULHjh2TYRgBv8woLi5WbGys\nOnXqpPDwcCUnJ6uwsDBgn4iICDkcDklSTU2N/3WwGR9m//f3tZYcHwAQXKaeUS1fvlzSuTWofig7\nO7vRr/d6vfJ4PP73Ho9H33777UX7bd26VW+++aYqKio0bdo0M6VdFt8Jr5T/8bk3mz+W7xcjuaoC\nAJszFVShmn2if//+6t+/v3bv3q3s7GzNnDnzkvvFxcVd0fG9761R5flbmD6frvl0vdyPZZr++u9/\nxLmtZMea7IT+NI4eNY4e1c/q/zeaCqqOHTv+qJO43W6Vl5f735eXl8vtdte7f69evbRkyRKdPHlS\n7dq1u+jzw4cPX3YNvhNeGRvWSbVnz22oPavKnHWqShl2WVdVV3JuK8XFxdmuJjuhP42jR42jR40L\nRn/qCztTQSVJ27Zt0+7du3Xy5MmA7RMnTmz0axMSEnTkyBGVlpbK7XYrPz9fkyZNCtjn6NGj6tSp\nkxwOh/bt26ezZ88GdYom48Ns6YfP1AxDxodrFZbxm6CdBwAQXKaC6u2339aGDRuUnJysgoICDRky\nRJs3b9Ydd9xh6iRhYWEaP3685s6dK8MwlJKSovj4eOXk5EiS0tLSVFBQoLy8PIWFhal169aaMmVK\ncAdU7Nsj1dUGbqurlfbuCd45AABBZyqoPv30U82YMUNdu3bVxo0bNXbsWA0cOFB//etfTZ8oKSlJ\nSUlJAdvS0tL8r9PT0y/62apgCpu10P+67pF7FLZinWXnAgAEj6nh6ZWVlerataskKTw8XLW1tUpM\nTNTu3bstLQ4AAFNXVLGxsfr+++8VHx/vv2UXGRmpyMhIq+sDALRwpielPXXqlCQpIyNDCxcuVHV1\nNZPSAgAsZyqoLny2lJiYqD/96U+WFQRc7Yx1b8p5z+imLgNoNkw9o8rNzVVJSUnAtv37919ypgoA\nDfOtZ/ou4HKYCqrs7OyAKZAkKSYmRmvX8hcOAGAtU0FVVVUll8sVsM3lcqmystKSogAAOM9UUHXp\n0kUFBQUB27Zu3crCiQAAy5kaTJGRkaEXX3xR+fn5io2N1dGjR7Vz505LZjgHAOBCpoKqR48emj9/\nvjZt2qSysjIlJiZq7NixiomJsbo+AEALZyqozp49q/bt2wdMcVRbW6uzZ8+qVatWlhUHAICpZ1TP\nP/+89u3bF7Bt3759mjt3riVFAQBwnqmgOnDggLp37x6wLTEx8aKfrQIAINhMBZXL5VJFRUXAtoqK\nCrVp08aSogAAOM9UUA0YMEALFy7UgQMHVFNTowMHDmjx4sWm16MCAOBKmRpMMXLkSK1Zs0bTp0/X\n2bNn1bp1aw0ePFijRo2yuj4AQAtnKqhat26thx9+WA899JBOnTqlqKio4K6+CwBAPUwF1XnV1dWq\nrq5WVVWVf1unTp2CXhQAAOeZCqqDBw9q0aJFlxzll52dHfSiYF7FX5ZLKcOaugwAsIypwRQrV65U\n79699dprr8nlcmn16tUaOnSoJkyYYHV9aMTJN1c0dQkAYClTQVVSUqKMjAy1bdtWPp9PLpdLY8aM\n4WoKAGA5U0HVqlUr1dXVSZKioqJUVlYmn8+n06dPW1ocAACmJ6XdsmWLBg8erNtvv10vvPCCWrVq\npd69e1tdHwCghTMVVE888YT/9ahRoxQfH6/q6moNGjTIssIAAJAuc3i6JDmdTgIKABAyDQZVY4Ml\nHA6HIiMjlZCQoBtuuCGohQEAIDUSVOXl5Y0e4MCBA3rvvfeUkpKi0aNHB60wAACkRoLqscceM3UQ\nr9erqVOnElQAgKAzNTy9MW63Ww888EAwDgUAaEZ8J7znfq84btk5ghJUkpSWlhasQwEAmgnjw+z/\n/r7WsnMELagAAC2L74RXyv/43JvNH1t2VUVQAQCuiPFhtmQY/31jWHZVZTqoDh06pHfeeUcrV670\nv7/UbOoAgKuf/2qqrvbchrpay66qTAXVli1bNHv2bHm9Xn322WeSzq1NtWbNmqAXBACwv4CrKf9G\na66qTM1M8dZbb2nGjBm6/vrrtWXLFknSddddp/379we9IABAM7Bvz/+ups6rq5X27gn6qUwFVUVF\nha677rqAbQ6Hg+XoAaCFCpu10P+67pF7FLZinWXnMnXrr1u3bsrLywvYtnnzZiUmJlpSFAAA55m6\noho3bpyef/55ffLJJ6qpqdHcuXN1+PBhzZgxw+r6AAAtnKmg6ty5s1555RV98cUX6tOnjzwej/r0\n6aOIiAir6wMAtHCml/lo06aNkpOTrazF1i6cJsQR3aGJqwGAlsNUUM2aNeuSAyfCw8Pl8XjUv39/\n9e3bt8FjFBUVafXq1TIMQ6mpqUpPTw/4/LPPPtMHH3wgn8+na665Rg8//LCuv/5689+JxS6cJiQs\n4zdNXA0AtBymBlP06tVLpaWl6tmzp+6880717NlTx44dU0JCgqKjo7V06VJ98MEH9X69YRhatWqV\npk+frgULFmjz5s06ePBgwD7XXnut5syZo/nz5+v+++/Xq6+++uO+syAK1TQhAICLmbqi2rFjh55+\n+ml16dLFv+3OO+9UVlaWXnjhBQ0YMEALFy7Uvffee8mvLy4uVmxsrDp16iRJSk5OVmFhYcDxLlx4\nsXv37qbWwgqVS00TwlUVAISGqSuqQ4cO+UPmvI4dO+rw4cOSpMTERJ04caLer/d6vfJ4PP73Ho9H\nXq+33v0/+eQT3XbbbWZKs1wopwkBAFzM1BVVz549tWTJEo0YMUJut1ter1dvvfWWevToIencKr8d\nOgRngMGuXbv06aef6tlnn613n7i4uB91ju8v4xje99ao0ucL3Ojz6ZpP18v9WOaPqiMYLud7aans\n1h87/jezWz12RI/qZ/WfaVNBNXHiRK1cuVJTpkyRYRgKCwtT//79/SsAh4eH6/HHH6/3691ud8Ct\nvPLycrnd7ov2Kykp0fLlyzVt2jRFRUXVe7zzV3I/htlj1O38Qqo9G7ix9qwqd3yh6iDUEQzB6MfV\nKi4uzpb9sVNNdu2RndCjxgWjP/WFnamgioyM1OTJk2UYhk6ePKl27drJ6fzfXcPGkjQhIUFHjhxR\naWmp3G638vPzNWnSpIB9ysrKNG/ePE2cONFW/3IJ5TQhAICLmf45KkmqqanRmTNndOzYMf+2Hz67\nupSwsDCNHz9ec+fOlWEYSklJUXx8vHJyciSdWx34nXfe0enTp/3LiISFhemll166nPIAAFchU0F1\n8OBBLVq06JLrT2VnZ5s6UVJSkpKSkgK2Xbh8/aOPPqpHH33U1LEAAC2HqVF/K1euVO/evfXaa6/J\n5XJp9erVGjp0qCZMmGB1fQCAFs5UUJWUlCgjI0Nt27aVz+eTy+XSmDFjTF9NAQBwpUwFVatWrVRX\nVydJioqKUllZmXw+n06fPm1pcQAAmHpG1aNHD23ZskWDBw/W7bffrhdeeEGtWrVS7969ra4PANDC\nmQqqJ554wv961KhRio+PV3V1tQYNGmRZYQAASCZu/RmGoTlz5ujs2XM/9Op0OjVo0CClpaWxHhUA\nwHKNBpXT6VRpaal8P5xGCMBlu3BdMwDmmBpMMXz4cK1YsULHjh2TYRgBvwCYd+G6ZgDMMfWMavny\n5ZKkvLy8iz5jiDpgzkXrmv1iJKtFAyaYCqrFixdbXQdw1WNdM+DKmAqqjh07Sjo3sKKioiJoS3oA\nLUW965pxVQU0ytQzqsrKSi1cuFAZGRn+Wc+3bdumtWu5zw6YEXA15d9o8KwKMMFUUK1YsUIul0tL\nlixRePi5i7Cf/vSnys/Pt7Q44Kqxb8//rqbOq6uV9u5pmnqAZsTUrb+dO3dq+fLl/pCSpHbt2qmi\nosKywoCrCeuaAVfO1BWVy+XSqVOnAraVlZXxrAoAYDlTQZWamqr58+dr165d8vl8+uabb5SVlaWh\nQ4daXR8AoIUzdevv3nvvVevWrbVq1SrV1dVp6dKlGjJkiO6++26r60MDLpzlgJFjAK5WpoLK4XDo\n7rvvJphs5sJZDvh5HABXK1O3/p566imtW7dO5eXlVtcDky6a5YC54wBcpUwF1QMPPKDi4mJNnjxZ\ns2fP1oYNG1g0sYldapYDALgambr1179/f/Xv319VVVX6/PPPtXnzZq1Zs0Y33nijpk6danWN+AFm\nOQDQkpgKqvOuueYaDRw4UG3btlVtba2+/PJLq+pCAxqa5YBnVQCuNqaCyufzadeuXdq0aZO2bt2q\njh07auDAgZowYYLV9eFSmOUAQAtiKqh+/etfKyIiQsnJyXruuefUpUsXq+tCA5jlAEBLYiqofve7\n3ykxMfGi7YZhyOk0NR4DAIArYiplfhhSBw4c0Jo1a/Sb3/A8BABgLdODKU6ePKlNmzYpNzdX+/fv\nV8+ePTV27FgLSwMAoJGgqq2t1bZt27Rx40b961//UmxsrP7v//5Px44d05QpUxQdHR2qOgEALVSD\nQfXII4/I6XTqrrvu0oMPPqhu3bpJknJyckJSHAAADT6juu6661RZWani4mLt3buX2SgAACHX4BXV\nnDlzdOzYMeXm5mr9+vVavXq1br75ZtXU1Kiuri5UNQIAWrBGB1N07NhRw4cP1/Dhw7Vnzx7l5ubK\n4XDoqaeeUkpKisaMGROKOgEALdRlTaHUo0cP9ejRQ+PGjdPWrVuVl5dnVV0AAEi6zKA6r3Xr1ho4\ncKAGDhwY7HoAAAjAtBIAAFtrkUHlGDayqUsAAJjUIoPKec/opi4BAGBSiwwqAEDzQVABAGztikb9\nXYmioiKtXr1ahmEoNTVV6enpAZ8fOnRIS5Ys0XfffaeRI0fqnnvuCVVpAAAbC8kVlWEYWrVqlaZP\nn64FCxZo8+bNOnjwYMA+kZGRGjdunIYNGxaKkgAAQdJu9COWHj8kQVVcXKzY2Fh16tRJ4eHhSk5O\nVmFhYcA+0dHRSkxMVFhYWChKAgAESXTGry09fkiCyuv1yuPx+N97PB55vd5QnBoA0MyF7BlVMMXF\nxTXJeb9vwnPXx4412Y3d+mPH/2Z2q8eO6FHDrOxPSILK7XarvLzc/768vFxut/uKj3f48OFglNXs\nzl0fO9ZkF3Fxcbbsj51qsmuP7IQeNSxY/akv7EJy6y8hIUFHjhxRaWmpamtrlZ+fr759+4bi1ACA\nZi4kV1RhYWEaP3685s6dK8MwlJKSovj4eP9KwWlpaTpx4oQyMzNVVVUlh8Ohv//97/rjH/8ol8sV\nihIBADYVsmdUSUlJSkpKCtiWlpbmf92+fXstW7YsVOUAAJoJZqYAANgaQQUAsDWCCgBgawQVAMDW\nCKrLwIKLABB6BNVlYMFFAAg9ggoAYGsEFQDA1ggqAICtEVQAAFsjqAAAtkZQAQBsjaACANgaQQUA\nsDWCCgBgawQVAMDWCCoAgK0RVECIMbkxcHkIKiDEmNwYuDwEFQDA1ggqAICtEVQAAFsjqJq5dqMf\naeoSAMBSBFUzF53x66YuAQAsRVABAGyNoAIA2BpBBQCwNYIKAGBrBBUAwNYIKgCArRFUAABbI6gA\nALbm8Pl8vqYuAgCA+nBFBQCwNYIKAGBrBBUAwNbCm7oAmLNkyRJt375d0dHRmj9/viTp9OnTWrBg\ngY4dO6aOHTtqypQpioyMbOJKm05ZWZmysrJ04sQJORwODRkyRHfffTd9+q8zZ85o9uzZqq2tVV1d\nnW6//XY9+OCD9OcSDMNQZmam3G63MjMz6dEPTJgwQREREXI6nQoLC9NLL71kaY8YTNFM7N69WxER\nEcrKyvIH1Z///GdFRkYqPT1d77//vk6fPq0xY8Y0caVN5/jx4zp+/Li6deumqqoqZWZm6qmnntLG\njRvpkySfz6eamhpFRESotrZWs2bN0tixY7V161b68wMffvih9u7d6/9zxN+1QBMmTNCLL76odu3a\n+bdZ2SNu/TUTvXr1uuhfJ4WFhbrrrrskSXfddZcKCwubojTb6NChg7p16yZJuuaaa9S5c2d5vV76\n9F8Oh0MRERGSpLq6OtXV1cnhcNCfHygvL9f27duVmprq30aPGmdlj7j114xVVFSoQ4cOkqT27dur\noqKiiSuyj9LSUn333XdKTEykTxcwDENTp07V0aNH9fOf/1zdu3enPz/w+uuva8yYMaqqqvJvo0cX\ne+655+R0OjV06FANGTLE0h4RVFcJh8Mhh8PR1GXYQnV1tebPn6+xY8fK5XIFfNbS++R0OvXyyy+r\nsrJS8+bN04EDBwI+b+n9+eKLLxQdHa1u3brpq6++uuQ+Lb1H0rmQcrvdqqio0PPPP6+4uLiAz4Pd\nI4KqGYuOjtbx48fVoUMHHT9+POB+cUtVW1ur+fPn684779SAAQMk0adLadu2rXr37q2ioiL6c4Gv\nv/5a27Zt05dffqkzZ86oqqpKixYtokc/4Ha7JZ37u9WvXz8VFxdb2iOeUTVjffv2VW5uriQpNzdX\n/fr1a+KKmpbP59OyZcvUuXNn/eIXv/Bvp0/nnDx5UpWVlZLOjQDcsWOHOnfuTH8uMHr0aC1btkxZ\nWVmaPHmybrzxRk2aNIkeXaC6utp/W7S6ulo7duxQ165dLe0Ro/6aiVdeeUW7d+/WqVOnFB0drQcf\nfFD9+vXTggULVFZWxpBZSXv27NGsWbPUtWtX/22HUaNGqXv37vRJUklJibKysmQYhnw+n+644w4N\nHz5cp06doj+X8NVXX2n9+vXKzMykRxf497//rXnz5kk6Nyhn4MCBuu+++yztEUEFALA1bv0BAGyN\noAIA2BpBBQCwNYIKAGBrBBUAwNYIKgCArRFUgM1MmDBBO3bs8L/fvHmzxo0bp927dzdhVUDTYQol\nwMY2btyoNWvWKDMzUzfccENTlwM0CYIKsKkNGzZo7dq1evrpp5WQkNDU5QBNhqACbCgnJ0dff/21\nZs6cqeuvv76pywGaFM+oABvauXOnunfvrq5duzZ1KUCTI6gAG3r44Yd15MgRLVu2TEzHiZaOoAJs\nKDo6WjNnztSePXu0cuXKpi4HaFIEFWBTbrdbs2bNUlFRkV5//fWmLgdoMgQVYGMxMTGaPXu2Pv/8\nc7355ptNXQ7QJFiPCgBga1xRAQBsjaACANgaQQUAsDWCCgBgawQVAMDWCCoAgK0RVAAAWyOoAAC2\nRlABAGzt/wFUUlT5y+B82QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x120dca890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use('ggplot')\n",
    "plt.errorbar(ks,jacs,jacs_var, linestyle='None', marker='^', elinewidth = 1.0)\n",
    "plt.xlabel('K')\n",
    "plt.ylabel('Average Jaccard Score')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Diffusion Distance/Markov Transition Matrix Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/architverma/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:3: RuntimeWarning: overflow encountered in exp\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n",
      "1.0\n",
      "10.0\n"
     ]
    }
   ],
   "source": [
    "ls = [0.1,0.5,1.0,10.0]\n",
    "mars = np.zeros(len(ls))\n",
    "mars_var = np.zeros(len(ls))\n",
    "for i in range(0,len(ls)):\n",
    "    print(ls[i])\n",
    "    [mars[i], mars_var[i]] = markov_RMSE(ctrl_pre.T[sample],ctrl.T[sample],ls[i],ls[i])\n",
    "    #print(np.mean(jaccard_score(ctrl_pre.T[sample],ctrl.T[sample],ks[i])))\n",
    "    #mars[i] = np.mean(t)\n",
    "    #mars_var[i] = np.std(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAEWCAYAAADy9UlpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X90VPWd//HXZIaEjTQhM9HEmGCaYKwUxWTHqNGExnhy\natWW2kN12bqLLOA6ujnqwYMCVU/dRCylQCyISDJSz3f7A3ehP7a75VQPsBK6JCahusKSCMqPhB2S\nCQRDY5y58/2DddZhQhKSyczAfT7O8ZA783l/7vt+znFeuXcmcy2BQCAgAABMICHWDQAAEC2EHgDA\nNAg9AIBpEHoAANMg9AAApkHoAQBMg9ADAJiGLdYNxFpHR8cFjc/KyrrgmksdazI41iUcaxKONQkX\niTXJysoa9HHO9AAApkHoAQBMg9ADAJgGoQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmEbVv\nZGltbZXb7ZZhGKqoqNCsWbNCng8EAnK73WppaVFSUpJcLpfy8vKGrH3jjTf07rvvymazKSMjQy6X\nS5dddpkkacuWLXr77beVkJCghx56SDfeeOO4HJfx639SwjfnRK0u0nMAQLw59f9elcrvHZe5o3Km\nZxiG6urqtGTJEq1atUq7du3S0aNHQ8a0tLTo+PHjqq2t1cKFC7Vx48Zha2+44QatXLlSP/rRj3Tl\nlVdqy5YtkqSjR4+qoaFBP/7xj7V06VLV1dXJMIxxObbAb34e1bpIzwEA8ab3n14bt7mjEnrt7e3K\nzMxURkaGbDabSkpK1NjYGDKmqalJZWVlslgsKigoUF9fn3p6eoasnTFjhqxWqySpoKBAXq9XktTY\n2KiSkhJNmDBBV1xxhTIzM9Xe3h6NQwUAxLGohJ7X65XD4QhuOxyOYEB9cUx6enrYmJHUStLbb78d\nvIR5bo3dbh+0BgBgLpfEXRb+5V/+RVarVaWlpRdce75v4h5pzZFRzjHaukjPESnx0ke8YV3CsSbh\nWJNQ4/naFpXQs9vt6u7uDm53d3fLbreHjenq6gob4/f7h6zdvn273n33XT377LOyWCyD7s/r9Ybt\n73ORuLXQaG+BEYnbicTDLUm4NcrgWJdwrEk41mRwF/WthfLz89XZ2SmPxyOfz6eGhgY5nc6QMU6n\nUzt37lQgENCBAweUnJystLS0IWtbW1v1q1/9SosXL1ZSUlLIXA0NDfrss8/k8XjU2dmpqVOnRuNQ\nAQBxLCpnelarVfPmzVN1dbUMw1B5eblycnK0bds2SVJlZaUKCwvV3NysqqoqJSYmyuVyDVkrSXV1\ndfL5fHrhhRckSddcc40WLlyonJwc3XrrrXryySeVkJCgv/u7v1NCAn+SCABmF7X39IqKilRUVBTy\nWGVlZfBni8Wi+fPnj7hWkl5++eXz7u++++7TfffdN8puAQCXIk5/AACmQegBAEyD0AMAmAahBwAw\nDUIPAGAahB4AwDQIPQCAaRB6AADTIPQAAKZB6AEATIPQAwCYBqEHADANQg8AYBqEHgDANAg9AIBp\nEHoAANMg9AAApkHoAQBMg9ADAJgGoQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmQegBAEyD\n0AMAmAahBwAwDUIPAGAahB4AwDRs0dpRa2ur3G63DMNQRUWFZs2aFfJ8IBCQ2+1WS0uLkpKS5HK5\nlJeXN2Tt7t27tXnzZh07dkw1NTXKz8+XJPl8Pm3YsEEffvihEhISNHfuXH31q1+N1qECAOJUVM70\nDMNQXV2dlixZolWrVmnXrl06evRoyJiWlhYdP35ctbW1WrhwoTZu3DhsbU5OjhYtWqTrrrsuZK4/\n/OEPkqSVK1dq2bJl+ulPfyrDMKJwpACAeBaV0Gtvb1dmZqYyMjJks9lUUlKixsbGkDFNTU0qKyuT\nxWJRQUGB+vr61NPTM2Rtdna2srKywvZ39OhRTZ8+XZKUmpqqyy67TAcPHhz/AwUAxLWohJ7X65XD\n4QhuOxwOeb3esDHp6elhY0ZSe67c3Fw1NTXJ7/fL4/Ho4MGD6urqitDRAAAuVlF7Ty+aysvLdfTo\nUT399NO6/PLLde211yohYfB8H+xMcThfrDkyyjlGWxfpOSIlXvqIN6xLONYkHGsSajxf26ISena7\nXd3d3cHt7u5u2e32sDFfPBv7fIzf7x+29lxWq1Vz584Nbi9btuy8C9jR0XEhh6KsrKywmgudY6x1\nkZ5jrAZbE7Aug2FNwrEmgxvrmpzvNT8qlzfz8/PV2dkpj8cjn8+nhoYGOZ3OkDFOp1M7d+5UIBDQ\ngQMHlJycrLS0tBHVnuvTTz9Vf3+/JOlPf/qTrFarsrOzx+34AAAXh6ic6VmtVs2bN0/V1dUyDEPl\n5eXKycnRtm3bJEmVlZUqLCxUc3OzqqqqlJiYKJfLNWStJO3Zs0f19fXq7e3V8uXLlZubq6VLl+rU\nqVOqrq5WQkKC7Ha7HnvssWgcJgAgzkXtPb2ioiIVFRWFPFZZWRn82WKxaP78+SOulaTi4mIVFxeH\nPX7FFVdozZo1Y+wYAHCp4RtZAACmQegBAEyD0AMAmAahBwAwDUIPAGAahB4AwDQIPQCAaRB6AADT\nIPQAAKZB6AEATIPQAwCYBqEHADANQg8AYBqEHgDANAg9AIBpEHoAANMg9AAApkHoAQBMg9ADAJgG\noQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmQegBAEyD0AMAmAahBwAwDUIPAGAahB4AwDRs\n0dpRa2ur3G63DMNQRUWFZs2aFfJ8IBCQ2+1WS0uLkpKS5HK5lJeXN2Tt7t27tXnzZh07dkw1NTXK\nz8+XJPl8Pq1fv16HDh2SYRgqKyvTt7/97WgdKgAgTkXlTM8wDNXV1WnJkiVatWqVdu3apaNHj4aM\naWlp0fHjx1VbW6uFCxdq48aNw9bm5ORo0aJFuu6660Lm+uMf/yifz6eVK1dq+fLl+sMf/iCPxxON\nQwUAxLGohF57e7syMzOVkZEhm82mkpISNTY2hoxpampSWVmZLBaLCgoK1NfXp56eniFrs7OzlZWV\nNeg++/v75ff7NTAwIJvNpuTk5HE/TgBAfItK6Hm9XjkcjuC2w+GQ1+sNG5Oenh42ZiS157rllls0\nceJELVy4UC6XS/fee68mTZoUoaMBAFysovaeXjS1t7crISFBr776qvr6+vTss8/q+uuvV0ZGRtjY\n850pDuWLNUdGOcdo6yI9R6TESx/xhnUJx5qEY01CjedrW1RCz263q7u7O7jd3d0tu90eNqarqyts\njN/vH7b2XO+8845uvPFG2Ww2paam6tprr9WHH344aOh1dHRc0LFkZWWF1VzoHGOti/QcYzXYmoB1\nGQxrEo41GdxY1+R8oRmVy5v5+fnq7OyUx+ORz+dTQ0ODnE5nyBin06mdO3cqEAjowIEDSk5OVlpa\n2ohqz5Wenq73339f0tn39tra2nTVVVeN2/EBAC4OUTnTs1qtmjdvnqqrq2UYhsrLy5WTk6Nt27ZJ\nkiorK1VYWKjm5mZVVVUpMTFRLpdryFpJ2rNnj+rr69Xb26vly5crNzdXS5cu1de//nWtW7dOTz75\npAKBgMrLy3X11VdH41ABAHEsau/pFRUVqaioKOSxysrK4M8Wi0Xz588fca0kFRcXq7i4OOzxiRMn\n6sknnxxjxwCASw3fyAIAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmMWzo1dfXh2y//fbbIds/\n+tGPItsRAADjZNjQ27FjR8j2G2+8EbL93nvvRbYjAADGybChFwgEotEHAADjbtjQs1gs0egDAIBx\nN+zXkPn9/uCXN0tn72R+7jYAABeDYUMvNTVVr7zySnB70qRJIdspKSnj0xkAABE2bOitXbs2Gn0A\nADDuRvV3eh0dHdqzZ49OnDgR6X4AABg3w57pbdq0SV/+8pdVVlYm6eyfMLzyyiu67LLL1N/fr0WL\nFqmwsHDcGwUAYKyGPdNrbGzUtGnTgts/+9nP9NBDD6murk4LFizQm2++Oa4NAgAQKcOG3unTp5We\nni5JOnz4sE6fPq077rhDklRWVqaOjo7x7TCOBU56z/57qicqdZGeAwDizXi/tg0besnJyTp58qQk\naf/+/crPz9eECRMkST6fb1yaulgYv/3F//7786jURXoOAIg34/3aNmzo3XrrrVqzZo1+97vfaevW\nrbr99tuDz7W3tysjI2NcGot3gZNeqeGtsxu73hrxbyWjrYv0HAAQb6Lx2jZs6M2ZM0fTpk3Tn/70\nJ91555268847g8999NFHIdtmYvz2F9Lnf5hvGCP+rWS0dZGeAwDiTTRe24b99KbNZtPs2bMHfe4b\n3/hGxBu6GAR/G/H/7+Vdv+/sbyX3PCBLalrE6yI9BwDEm2i9tg0beufeZWEwM2fOjEgzF4uQ30aC\nD579rcT6149EvC7ScwBAvInWa9uwobdu3TplZmZq8uTJg95xwWKxmC70dHD///028jm/T/pw//jU\nRXoOAIg3UXptGzb07rrrLv3xj3/UxIkTNXPmTN10003BT2+alfXZNcGf/Qu+Ketrvx7XukjPAQDx\nJlqvbcOG3ty5c/U3f/M3am1t1Y4dO/T666+rqKhIX/va1/SVr3xlXJoCAGA8jOi7NxMSElRUVKQn\nnnhCq1ev1qRJk/T888+H3GIIAIB4N+yZ3ufOnDmjXbt2aceOHert7dV3vvMd5ebmjmNrAABE1rCh\n19TUpJ07d2r//v1yOp363ve+x2VNAMBFadjQW7FihbKyslRaWqrExETt3btXe/fuDRlz//33j1uD\nAABEyrChV1ZWJovFotOnT0ejHwAAxs2woffoo4+e97mPP/5Y//zP/zyiHbW2tsrtdsswDFVUVGjW\nrFkhzwcCAbndbrW0tCgpKUkul0t5eXlD1u7evVubN2/WsWPHVFNTo/z8fEnSf/zHf+jXv/6/j7se\nPnxYL730Eu9BAoDJDRt6n376qbZs2aKPPvpIV155pWbPnq3Tp0/rpz/9qd57773gzWWHYhiG6urq\ntGzZMjkcDj3zzDNyOp3Kzs4OjmlpadHx48dVW1urtrY2bdy4UTU1NUPW5uTkaNGiRdqwYUPI/kpL\nS1VaWirpbOCtWLGCwAMADB96dXV1OnTokGbMmKHW1lYdPnxYHR0dmjlzph5++GGlpKQMu5P29nZl\nZmYG78hQUlKixsbGkNBramoKXkotKChQX1+fenp6dOLEifPWfrH+fN555x2VlJQMOw4AcOkbNvT2\n7t2rH/7wh0pNTdVdd90ll8ul559/Xtddd92Id+L1euVwOILbDodDbW1tYWM+v1nt52O8Xu+Iaoey\ne/duPfXUUyMeDwC4dA0bev39/UpNTZV0NnAmTpx4QYEXS21tbUpMTNSUKVPOOyYrK+uC5/1izZFR\nzjHaukjPESnx0ke8YV3CsSbhWJNQ4/naNmzo+f3+sG9eOXd7+vTpQ85ht9vV3d0d3O7u7pbdbg8b\n09XVFTbG7/cPW3s+u3bt0m233TbkmI6OjhHN9bmsrKywmgudY6x1kZ5jrAZbE7Aug2FNwrEmgxvr\nmpwvNIcNvdTUVL3yyivB7UmTJoVsWywW/eQnPxlyjvz8fHV2dsrj8chut6uhoUFVVVUhY5xOp/79\n3/9dt912m9ra2pScnKy0tDSlpKQMWzsYwzC0e/du/eAHPxh2LADAHIYNvbVr1455J1arVfPmzVN1\ndbUMw1B5eblycnK0bds2SVJlZaUKCwvV3NysqqoqJSYmyuVyDVkrSXv27FF9fb16e3u1fPly5ebm\naunSpZKkffv2KT09PfgBGAAARvzdm2NVVFSkoqKikMcqKyuDP1ssFs2fP3/EtZJUXFys4uLiQWu+\n+tWvqrq6egwdAwAuNSO6ywIAAJcCQg8AYBqEHgDANAg9AIBpEHoAANMg9AAApkHoAQBMg9ADAJgG\noQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmQegBAEyD0AMAmAahBwAwDUIPAGAahB4AwDQI\nPQCAaRB6AADTIPQAAKZB6AEATIPQAwCYBqEHADANQg8AYBqEHgDANAg9AIBpEHoAANMg9AAApmGL\n1o5aW1vldrtlGIYqKio0a9askOcDgYDcbrdaWlqUlJQkl8ulvLy8IWt3796tzZs369ixY6qpqVF+\nfn5wvo8//lgbNmzQn//8Z1ksFr344otKTEyM1uECAOJQVELPMAzV1dVp2bJlcjgceuaZZ+R0OpWd\nnR0c09LSouPHj6u2tlZtbW3auHGjampqhqzNycnRokWLtGHDhpD9+f1+vfzyy3rssceUm5ur06dP\ny2aLWr4DAOJUVC5vtre3KzMzUxkZGbLZbCopKVFjY2PImKamJpWVlclisaigoEB9fX3q6ekZsjY7\nO1tZWVlh+9u7d6+mTJmi3NxcSdKXvvQlJSRwJRcAzC4qpz9er1cOhyO47XA41NbWFjYmPT09ZIzX\n6x1R7bk6OztlsVhUXV2t3t5elZSU6Fvf+laEjgYAcLG6JK/5+f1+7d+/Xy+++KKSkpL0gx/8QHl5\nebr++uvDxg52pjicL9YcGeUco62L9ByREi99xBvWJRxrEo41CTWer21RCT273a7u7u7gdnd3t+x2\ne9iYrq6usDF+v3/Y2nM5HA5dd911SklJkSQVFhbq0KFDg4ZeR0fHBR1LVlZWWM2FzjHWukjPMVaD\nrQlYl8GwJuFYk8GNdU3OF5pReaMrPz9fnZ2d8ng88vl8amhokNPpDBnjdDq1c+dOBQIBHThwQMnJ\nyUpLSxtR7blmzJihI0eO6NNPP5Xf79e+fftCPjQDADCnqJzpWa1WzZs3T9XV1TIMQ+Xl5crJydG2\nbdskSZWVlSosLFRzc7OqqqqUmJgol8s1ZK0k7dmzR/X19ert7dXy5cuVm5urpUuXatKkSbr77rv1\nzDPPyGKxqLCwUEVFRdE4VABAHIvae3pFRUVhwVNZWRn82WKxaP78+SOulaTi4mIVFxcPWlNWVqay\nsrIxdAwAuNTwOX4AgGkQegAA0yD0AACmQegBAEyD0AMAmAahBwAwDUIPAGAahB4AwDQIPQCAaRB6\nAADTIPQAAKZB6AEATIPQAwCYBqEHADANQg8AYBqEHgDANAg9AIBpEHoAANMg9AAApkHoAQBMg9AD\nAJgGoQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA0yD0AACmQegBAEyD0AMAmAahBwAwDVu0dtTa2iq3\n2y3DMFRRUaFZs2aFPB8IBOR2u9XS0qKkpCS5XC7l5eUNWbt7925t3rxZx44dU01NjfLz8yVJHo9H\nTzzxhLKysiRJ11xzjRYuXBitQwUAxKmohJ5hGKqrq9OyZcvkcDj0zDPPyOl0Kjs7OzimpaVFx48f\nV21trdra2rRx40bV1NQMWZuTk6NFixZpw4YNYfvMzMzUihUronF4AICLRFQub7a3tyszM1MZGRmy\n2WwqKSlRY2NjyJimpiaVlZXJYrGooKBAfX196unpGbI2Ozs7eDYHAMBwohJ6Xq9XDocjuO1wOOT1\nesPGpKenh40ZSe1gPB6PnnrqKT333HPat29fBI4CAHCxi9p7etGUlpamdevW6Utf+pIOHjyoFStW\naOXKlUpOTg4bO5ozxS/WHBnlHKOti/QckRIvfcQb1iUcaxKONQk1nq9tUQk9u92u7u7u4HZ3d7fs\ndnvYmK6urrAxfr9/2NpzTZgwQRMmTJAk5eXlKSMjQ52dncEPunxRR0fHBR1LVlZWWM2FzjHWukjP\nMVaDrQlYl8GwJuFYk8GNdU3OF5pRubyZn5+vzs5OeTwe+Xw+NTQ0yOl0hoxxOp3auXOnAoGADhw4\noOTkZKWlpY2o9ly9vb0yDEOS9D//8z/q7OxURkbGuB0fAODiEJUzPavVqnnz5qm6ulqGYai8vFw5\nOTnatm2bJKmyslKFhYVqbm5WVVWVEhMT5XK5hqyVpD179qi+vl69vb1avny5cnNztXTpUn3wwQf6\n5S9/KavVqoSEBC1YsECTJk2KxqECAOJY1N7TKyoqUlFRUchjlZWVwZ8tFovmz58/4lpJKi4uVnFx\ncdjjt9xyi2655ZYxdgwAuNTwjSwAANMg9AAApkHoAQBMg9ADAJgGoQcAMA1CDwBgGoQeAMA0CD0A\ngGkQemNkufeBqNZFeg4AiDcpcxaM29yE3hglfHNOVOsiPQcAxJvUv3543OYm9AAApkHoAQBMg9AD\nAJgGoQcAMA1CDwBgGoQeAMA0CD0AgGkQegAA07AEAoFArJsAACAaONMDAJgGoQcAMA1CDwBgGrZY\nN3AxaW1tldvtlmEYqqio0KxZs2LdUkx1dXVp7dq1OnnypCwWi+6880594xvfiHVbccEwDD399NOy\n2+16+umnY91OzPX19Wn9+vU6cuSILBaLHnnkERUUFMS6rZj67W9/q7ffflsWi0U5OTlyuVxKTEyM\ndVtRtW7dOjU3Nys1NVUrV66UJH3yySdatWqVTpw4ocsvv1xPPPGEJk2aFLF9cqY3QoZhqK6uTkuW\nLNGqVau0a9cuHT16NNZtxZTVatWDDz6oVatWqbq6Wr///e9Nvyaf+93vfqerrroq1m3EDbfbrRtv\nvFGrV6/WihUrTL82Xq9X//Zv/6bly5dr5cqVMgxDDQ0NsW4r6r72ta9pyZIlIY9t3bpV119/vWpr\na3X99ddr69atEd0noTdC7e3tyszMVEZGhmw2m0pKStTY2BjrtmIqLS1NeXl5kqS/+Iu/0FVXXSWv\n1xvjrmKvu7tbzc3NqqioiHUrceHMmTPat2+f7rjjDkmSzWbTZZddFuOuYs8wDA0MDMjv92tgYEBp\naWmxbinqpk2bFnYW19jYqJkzZ0qSZs6cGfHXWS5vjpDX65XD4QhuOxwOtbW1xbCj+OLxeHTo0CFN\nnTo11q3E3Ouvv67vfe97+vOf/xzrVuKCx+NRSkqK1q1bp48//lh5eXmaO3euJk6cGOvWYsZut+ve\ne+/VI488osTERM2YMUMzZsyIdVtx4dSpU8FfACZPnqxTp05FdH7O9DBm/f39WrlypebOnavk5ORY\ntxNT7777rlJTU4NnwJD8fr8OHTqkyspK/fCHP1RSUlLEL1ldbD755BM1NjZq7dq1evXVV9Xf36+d\nO3fGuq24Y7FYZLFYIjonoTdCdrtd3d3dwe3u7m7Z7fYYdhQffD6fVq5cqdLSUt18882xbifm/vu/\n/1tNTU169NFHtXr1ar3//vuqra2NdVsx5XA45HA4dM0110iSbrnlFh06dCjGXcXWe++9pyuuuEIp\nKSmy2Wy6+eabdeDAgVi3FRdSU1PV09MjSerp6VFKSkpE5yf0Rig/P1+dnZ3yeDzy+XxqaGiQ0+mM\ndVsxFQgEtH79el111VW65557Yt1OXJgzZ47Wr1+vtWvX6vHHH9f06dNVVVUV67ZiavLkyXI4HOro\n6JB09gU/Ozs7xl3FVnp6utra2vTpp58qEAjovffeM/2Hez7ndDq1Y8cOSdKOHTt00003RXR+vobs\nAjQ3N2vTpk0yDEPl5eW67777Yt1STO3fv1/PPvuspkyZErwE8Vd/9VcqKiqKcWfx4b/+67/0m9/8\nhj9ZkPTRRx9p/fr18vl8uuKKK+RyuSL6MfSL0S9/+Us1NDTIarUqNzdXf//3f68JEybEuq2oWr16\ntT744AOdPn1aqamp+u53v6ubbrpJq1atUldX17j8yQKhBwAwDS5vAgBMg9ADAJgGoQcAMA1CDwBg\nGoQeAMA0CD0gQjZs2KA333wzuL1t2zYtWLBADz74oE6fPq39+/erqqpKDz74oPbs2TPq/dTU1Gj7\n9u0R6Dg+rV27Vj//+c9j3QYuUXz3JjACjz76qE6ePCmr1aqEhARlZ2errKxMd955pxISzv7uuHDh\nwuB4n8+nTZs2qbq6Wrm5uZLO/l3W17/+9THffuncb6WPlLVr18rhcOiBBx4Yl/kHs337dr311lt6\n4YUXorZPmBuhB4zQ4sWLdcMNN+jMmTP64IMP5Ha71d7eLpfLFTb21KlT+uyzz5STkxN87MSJE6b/\nJhIg1gg94AIlJyfL6XRq8uTJWrp0qe655x5NmTIleKZUVlamxYsXS5Lmzp2rqVOnqqurSx6PRy+9\n9JISEhJUX1+vxx9/XA8//LBuuOEGSWfPBI8fP66qqioNDAxo/fr1am1tlWEYuvLKK7V48WJNnjxZ\nzz//vEpLS1VRUSHDMLRlyxa99dZbGhgY0I033qh58+YpOTlZHo9Hjz32mFwul37xi19oYGBAd999\n96i+SejYsWOqr6/XwYMHlZKSovvvv18lJSWSzp4hJiUl6cSJE9q3b5+ys7NVVVWlzMxMSdLevXtV\nX1+vkydPqrS0VEeOHFFZWZmuvfZavfbaa/L5fHrwwQdltVr1+uuvSzr7hcwvvvjioPMBY8F7esAo\nTZ06VXa7Xfv37w95PCsrK3gX6Ndff13PPfecXn75ZaWnp2vx4sV64403hv26qR07dujMmTN65ZVX\nVF9frwULFgx6V+3t27dr+/bteu655/STn/xE/f39qqurCxmzf/9+rVmzRt///vf15ptvXvCNfvv7\n+/WP//iPuv3227Vx40Y9/vjjqqurC5mnoaFBs2fPltvtVmZmZvA9ud7eXv34xz/WnDlzVF9fr6ys\nrOAXK2dnZ2vBggUqKCjQG2+8EQy8oeYDxorQA8bAbrfrk08+ifi8VqtVn3zyiY4fP66EhATl5eUN\netumd955R/fcc48yMjI0ceJEzZkzRw0NDfL7/cExs2fPVmJionJzc3X11Vfr448/vqBempubdfnl\nl6u8vFxWq1Vf/vKXdfPNN2v37t3BMcXFxZo6daqsVqtuv/12ffTRR5KklpYWZWdn6+abb5bVatVd\nd92lyZMnD7vP880HjBWXN4Ex8Hq94/LFyWVlZeru7tbq1at15swZlZaW6oEHHpDNFvq/bE9Pjy6/\n/PLgdnp6uvx+f8iNN78YMklJServ77+gXk6cOKG2tjbNnTs3+Jjf71dZWdmw++jp6Qm5+bLFYhnR\nLbnG2jNBriRJAAACQ0lEQVRwPoQeMErt7e3yer36yle+Mqr6pKQkDQwMBLdPnjwZ/Nlms2n27Nma\nPXu2PB6PXnzxRWVlZemOO+4ImSMtLU0nTpwIbnd1dclqtSo1NTXk/o9j4XA4NG3aNH3/+9+/4NrJ\nkyfL6/UGtwOBQMg2EG1c3gQu0JkzZ/Tuu+9qzZo1Ki0t1ZQpU0Y1T25urnbt2iWfz6cPP/xQ//mf\n/xl87v3339fhw4dlGIaSk5Nls9kGvYP0bbfdpn/913+Vx+NRf3+/fvazn+nWW2+V1WodVU+GYWhg\nYCD4n8/n01/+5V+qs7NTO3fulM/nk8/nU3t7+4jeGywqKtLhw4e1Z88e+f1+/f73vw8J989D0efz\njapf4EJxpgeM0EsvvSSr1SqLxaLs7GzdfffdqqysHPV8999/v9asWaOHHnpI06ZN02233RZ8f/Dk\nyZN67bXX5PV6NXHiRN16660hlxM/V15erp6eHj333HMaGBjQjBkzNG/evFH3tHXrVm3dujW4fe21\n1+qFF17QsmXLtGnTJm3atEmBQEBXX321/vZv/3bY+VJSUvTkk0/K7XZr7dq1Ki0tVV5eXvCDPNOn\nTw9+oCUhISHsQzhApHE/PQBRYxiGHnnkEf3DP/yDpk+fHut2YEJc3gQwrlpbW9XX16fPPvtMW7Zs\nUSAQUEFBQazbgklxeRPAuDpw4IBqa2vl8/mUnZ2tp556atC/OQSigcubAADT4PImAMA0CD0AgGkQ\negAA0yD0AACmQegBAEyD0AMAmMb/BwLNziK2CQX5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f940890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(ls,mars,mars_var, linestyle='None', marker='^', elinewidth = 1.0)\n",
    "plt.xlabel('Diffusion Length')\n",
    "plt.ylabel('RMSE')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}