GP_multichannel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'o2')"
      ]
     },
     "execution_count": 247,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAFACAYAAADasjVjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFPNJREFUeJzt3V9sZHd5xvHnTewdY0fCbNZxliTt5sICRVQ0rFtSoFFEuAglIukNhYoWUeiqEuAEaGnKBbNz0SqoFCUXFcpogaQiCqAQCkKIQtNSqIhcnCwSCVtqFELYsOOYGgfsLbO2eHvh2XC8sHlnxufPzG++Hyla++z4PL+fSR5mfMbnNXcXAOD8Lqh6AQAw6ChKAAhQlAAQoCgBIEBRAkCAogSAAEUJAAGKEgACFCUABMaqXkA3Dhw44IcOHap6GQAS8/DDD//Y3Weixw1FUR46dEhLS0tVLwNAYszsB908jpfeABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKEoACFCUABCgKAEgQFECQKCw3/U2s49JulHS0+7+ks6x/ZI+JemQpCckvcHdf1JE/trampaXl7W+vq7p6WnNzc1p//79ZFSQkcIeyOhNs7mhen1cKyv7NDt7Ro3Glo4cuSjXjDL2cVaRzyjvlnTDOcduk/Sgu89JerDzee7W1ta0uLiodrutiy++WO12W4uLi1pbWyOj5IwU9kBGb5rNDS0sTKjVqsnd1GrVtLAwoWZzI7eMMvaRVVhRuvvXJJ276psk3dP5+B5JNxeRvby8rKmpKU1NTcnMnv14eXmZjJIzUtgDGb2p18fVbu9+sdpuj6leH88to4x9ZJX9M8pZdz/V+bglafZ8DzSzI2a2ZGZLq6urPYWsr69rcnJy17HJyUmtr6/3ul4yBvz8ZAxexsrKvp6O96OMfWRVdjHH3V2SP8ffN9193t3nZ2bC+2ruMj09rdOnT+86dvr0aU1PT/e1VjIG9/xkDF7G7OyZno73o4x9ZJVdlCtmdlCSOn8+XUTI3NycNjc3tbm5KXd/9uO5uTkySs5IYQ9k9KbR2FKttr3rWK22rUZjK7eMMvaRZTtP7IphZockfSFz1fvvJf2vu99uZrdJ2u/u74vOMz8/773e4TyVq4cpZKSwBzJ6MyxXvc3sYXefDx9XVFGa2X2SrpN0QNKKpLqkf5b0aUm/IekH2nl7UHiZqp+iBIBIt0VZ2Pso3f1N5/mr64vKBIAi8Js5ABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKEoACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAECrsfZdVSuVN0Chkp7CGljGG5+/ggZJyV5DPKVOYjp5CRwh5Sykhl5nYyc72rlMp85BQyUthDShmpzNxOfa53KVKZj5xCRgp7SCkjlZnbIzPXu0ipzEdOISOFPaSUkcrM7dTnepcilfnIKWSksIeUMlKZuZ3UXO+8MNd7uDNS2ENKGVz1/qXK53rnibneAIrQbVEm+dIbAPJEUQJAgKIEgABFCQABihIAAhQlAAQoSgAIUJQAEKAoASBAUQJAgKIEgABFCQABihIAAhQlAAQoSgAIUJQAEKikKM3s3Wb2mJk9amb3mdlEFesAgG6MxQ/Jl5ldJmlB0lXu/n9m9mlJb5R0d545w3Ir+lHISGEPZWUwpmGwMs6q6qX3mKTnmdmYpElJP8rz5KkMYE8hI4U9lJXRbG5oYWFCrVZN7qZWq6aFhQk1mxu5ZaTyvSojI6v0onT3pyR9SNKTkk5Jesbdv5xnRioD2FPISGEPZWXU6+Nqt3e/yGu3x1Svj+eWkcr3qoyMrNKL0sxeIOkmSVdKeqGkKTN786953BEzWzKzpdXV1Z4yUhnAnkJGCnsoK2NlZV9Px/uRyveqjIysKl56v0bS99191d23JD0g6RXnPsjdm+4+7+7zMzMzPQWkMoA9hYwU9lBWxuzsmZ6O9yOV71UZGVlVFOWTkq4xs0kzM0nXSzqRZ0AqA9hTyEhhD2VlNBpbqtW2dx2r1bbVaGzllpHK96qMjKxK5nqbWUPSH0nalnRc0tvdvX2+x/cz1zuVq24pZKSwh7IyuOpdbka3c70rKcpe9VOUABDptij5zRwACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKEoACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAESp/rXZZhucPyKGSksAeJu4+PYsZZST6jTGWucAoZKexBYub2KGZkJVmUqcwVTiEjhT1IzNwexYysJIsylbnCKWSksAeJmdujmJGVZFGmMlc4hYwU9iAxc3sUM7KSLMpU5gqnkJHCHiRmbo9iRlay42pTueqWQkYKe5C46p1iBnO9ASDAXG8AyAlFCQABihIAAhQlAAQoSgAIUJQAEKAoASBAUQJAgKIEgABFCQABihIAAhQlAAQoSgAIUJQAEKAoASBAUQJAoJK53mY2LemYpJdIckl/5u4P5ZlRxh2WN5pNjdfr2reyojOzs9pqNHTRkSO5ZgzLnaKrPL/E3cfJSHOu952SvuTuL5b0Ukkn8jx5GTN/N5pNTSwsqNZqydxVa7U0sbCgjWYzt4wU5iMzc5uMYc3IKr0ozez5kq6V9FFJcvcz7p7rjMkyZv6O1+saa7d3HRtrtzVer+eWkcJ8ZGZukzGsGVlVPKO8UtKqpI+b2XEzO2ZmU+c+yMyOmNmSmS2trq72FFDGzN99Kys9He9HCvORmblNxrBmZFVRlGOSXibpI+5+taRNSbed+yB3b7r7vLvPz8zM9BRQxszfM7OzPR3vRwrzkZm5TcawZmRVUZQnJZ1098XO5/drpzhzU8bM361GQ9u12q5j27WathqN3DJSmI/MzG0yhjUjq5JxtWb2dUlvd/fvmtlRSVPu/lfne/ygzvXmqvdgnF/iqjcZCc71NrPf1s7bg/ZJelzSW939J+d7PHO9ARSh26Ks5H2U7v4tSeHiAGAQ8Js5ABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKEoACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAJhUZrZi83sejO76JzjNxS3LAAYHM9ZlGa2IOlzkt4l6VEzuynz139X5MIAYFBEN+79c0mH3X3DzA5Jut/MDrn7nZKs6MXtBaMgBieDMQ1kDGvGWdFL7wvcfUOS3P0JSddJeq2ZfVgDXJRlDEffaDY1sbCgWqslc1et1dLEwoI2ms3cMlIYJN9sbmhhYUKtVk3uplarpoWFCTWbG7mcX0rj+0TG4GVkRUW50plvI0nqlOaNkg5I+q1CVpSDMoajj9frGmu3dx0ba7c1Xq/nlpHCIPl6fVzt9u4XLu32mOr18VzOL6XxfSJj8DKyoqL8U0ktSTKzl5rZOyX9haR/kHRtISvKQRnD0fetrPR0vB8pDJJfWdnX0/F+pPB9ImPwMrKesyjd/aS7t8zsFkn3Srqk888nlPMs7jyVMRz9zOxsT8f7kcIg+dnZMz0d70cK3ycyBi8jq9v3Ub5N0svd/QPu/gFJ12jnQs9AKmM4+lajoe1abdex7VpNW41GbhkpDJJvNLZUq23vOlarbavR2Mrl/FIa3ycyBi8jq6u53mb2bUm/4+4/73w+Iemb7l7Kzyn7mevNVe/ByeCqNxmDmtHtXO9ui/I9kt4i6bOdQzdLutvd7+hpVX3qpygBINJtUUbvo5QkufuHzeyrkl7VOfRWdz++h/UBwNDoqiglyd0fkfRIgWsBgIHETTEAIEBRAkCAogSAAEUJAAGKEgACFCUABChKAAhQlAAQoCgBIEBRAkCAogSAAEUJAAGKEgACFCUABLq+zVrezOxCSUuSnnL3G/M+P3c4717RdyAflrtdk0HG+VT5jPIWSSeKODFzvbtX9NztVGY8kzF6GVmVFKWZXS7pdZKOFXF+5np3r+i526nMeCZj9DKyqnpGeYek90n6xfkeYGZHzGzJzJZWV1d7OjlzvbtX9NztVGY8kzF6GVmlF6WZ3SjpaXd/+Lke5+5Nd5939/mZmZmeMpjr3b2i526nMuOZjNHLyKriGeUrJb3ezJ6Q9ElJrzazT+QZwFzv7hU9dzuVGc9kjF5GVlfjaotiZtdJ+svoqjdzvbnqTQYZRWTkOte7KEUWJQBEcp3rXRR3/6qkr1a5BgCI8Js5ABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKEoACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAIUJQAEKr0fZZG4w/ngZKSwBzJGM+OsJJ9RMtd7cDJS2AMZo5mRlWRRMtd7cDJS2AMZo5mRlWRRMtd7cDJS2AMZo5mRlWRRMtd7cDJS2AMZo5mRlWRRMtd7cDJS2AMZo5mRVem42m4x13u4rx6msAcy0swYirne3WKuN4AidFuUSb70BoA8UZQAEKAoASBAUQJAgKIEgABFCQABihIAAhQlAAQoSgAIUJQAEKAoASBAUQJAgKIEgABFCQABihIAAhQlAARKn+ttZldI+idJs5JcUtPd78w7hzucD05GCnsgYzQzzqriGeW2pPe6+1WSrpH0DjO7Ks8A5noPTkYKeyBjNDOySi9Kdz/l7o90Pv6ZpBOSLsszg7neg5ORwh7IGM2MrEp/RmlmhyRdLWnx1/zdETNbMrOl1dXVns7LXO/ByUhhD2SMZkZWZUVpZhdJ+oykW939p+f+vbs33X3e3ednZmZ6OjdzvQcnI4U9kDGaGVmVFKWZjWunJO919wfyPj9zvQcnI4U9kDGaGVmlj6s1M5N0j6Q1d7+1m69hrvdwXz1MYQ9kpJkxsHO9zexVkr4u6duSftE5/H53/+L5voa53gCK0G1Rlv4+Snf/T0lWdi4A9IvfzAGAAEUJAAGKEgACFCUABChKAAhQlAAQoCgBIEBRAkCAogSAAEUJAAGKEgACFCUABChKAAhQlAAQoCgBIEBRAkAg2aJsPtTUwQ8e1AWNC3TwgwfVfCi/edtkDNb5ySCjqIyzLjx69GhhJ89Ls9k8eqSHWTTNh5pa+NcFPbP1jCRpY3tDX/n+V3TJvkt0+IrD+ayJjIE4Pxlk7CWj0WicOnr0aNiwpc/M6UevM3MOfvCgWj9v/crxSycu1am/PpXLmsgYjPOTQcZeMrqdmZPkS++Vn6/0dJyM4jJS2AMZo5mRlWRRzk7M9nScjOIyUtgDGaOZkZVkUTaubah2QW3XsdoFNTWubZBRckYKeyBjNDN2cfeB/+fw4cPeq7u+cZdfevulbkfNL739Ur/rG3f1fA4yhuP8ZJDRb4akJe+ig5K8mAMA3RjpizkAkCeKEgACFCUABChKAAhQlAAQoCgBIEBRAkCAogSAAEUJAAGKEgACFCUABChKAAhQlAAQoCgBIEBRAkCAogSAwFgVoWZ2g6Q7JV0o6Zi73553xtrampaXl7W+vq7p6WnNzc1p//79ZFSQkcIeyBjNjLNKf0ZpZhdK+kdJr5V0laQ3mdlVeWasra1pcXFR7XZbF198sdrtthYXF7W2tkZGyRkp7IGM0czIquKl9+9K+p67P+7uZyR9UtJNeQYsLy9rampKU1NTMrNnP15eXiaj5IwU9kDGaGZkVVGUl0n6Yebzk51ju5jZETNbMrOl1dXVngLW19c1OTm569jk5KTW19f7WC4Zg3x+MsgoKiNrYC/muHvT3efdfX5mZqanr52entbp06d3HTt9+rSmp6dzWx8Zg3F+MsgoKiOriqJ8StIVmc8v7xzLzdzcnDY3N7W5uSl3f/bjubk5MkrOSGEPZIxmRlbp42rNbEzS/0i6XjsF+U1Jf+zuj53va/oZV5vKVbcUMlLYAxlpZnQ7rraSud5m9geS7tDO24M+5u5/+1yPZ643gCJ0W5SVvI/S3b8o6YtVZANArwb2Yg4ADAqKEgACFCUABChKAAhQlAAQoCgBIEBRAkCgkjec98rMViX9oM8vPyDpxzkupyop7COFPUjsY9DsZR+/6e7hzSSGoij3wsyWunnn/aBLYR8p7EFiH4OmjH3w0hsAAhQlAARGoSibVS8gJynsI4U9SOxj0BS+j+R/RgkAezUKzygBYE8oSgAIJFuUZnaDmX3XzL5nZrdVvZ5+mNkVZvbvZvYdM3vMzG6pek17YWYXmtlxM/tC1Wvpl5lNm9n9ZvbfZnbCzH6v6jX1ysze3fn36VEzu8/MJqpeUzfM7GNm9rSZPZo5tt/MvmJmy50/X1BEdpJFWcbs8JJsS3qvu18l6RpJ7xjSfZx1i6QTVS9ij+6U9CV3f7Gkl2rI9mNml0lakDTv7i/RzpSBN1a7qq7dLemGc47dJulBd5+T9GDn89wlWZQqYXZ4Gdz9lLs/0vn4Z9r5j/JXRvsOAzO7XNLrJB2rei39MrPnS7pW0kclyd3PuHsx81GLNSbpeZ35VZOSflTxerri7l+TtHbO4Zsk3dP5+B5JNxeRnWpRdjU7fJiY2SFJV0tarHYlfbtD0vsk/aLqhezBlZJWJX288yOEY2Y2VfWieuHuT0n6kKQnJZ2S9Iy7f7naVe3JrLuf6nzckjRbREiqRZkUM7tI0mck3eruP616Pb0ysxslPe3uD1e9lj0ak/QySR9x96slbaqgl3pF6fwM7ybtlP4LJU2Z2ZurXVU+fOe9joW83zHVoix8dnhZzGxcOyV5r7s/UPV6+vRKSa83sye082OQV5vZJ6pdUl9OSjrp7mef1d+vneIcJq+R9H13X3X3LUkPSHpFxWvaixUzOyhJnT+fLiIk1aL8pqQ5M7vSzPZp54fVn694TT0zM9POz8NOuPuHq15Pv9z9b9z9cnc/pJ3/Lf7N3YfuWYy7tyT90Mxe1Dl0vaTvVLikfjwp6Rozm+z8+3W9huyC1Dk+L+ktnY/fIulzRYRUMq62aO6+bWbvlPQv+uXs8McqXlY/XinpTyR928y+1Tn2/s64X1TjXZLu7fwf8OOS3lrxenri7otmdr+kR7TzrorjGpJfZTSz+yRdJ+mAmZ2UVJd0u6RPm9nbtHMrxjcUks2vMALAc0v1pTcA5IaiBIAARQkAAYoSAAIUJQAEKEok59fdZQbYC4oSKbpbv3qXGaBvFCWGmpm9p3NfxUfN7FbpvHeZAfqW5G/mYDSY2WHt/GbMyyWZpEUz+w93P17typAanlFimL1K0mfdfdPdN7Rzg4ffr3hNSBBFCQABihLD7OuSbu7cCWdK0h92jgG5oigxtDpjMu6W9F/aufP7MXc/3rnLzEOSXmRmJzt3lgH6xt2DACDAM0oACFCUABCgKAEgQFECQICiBIAARQkAAYoSAAL/D7Yr+hpw7yDqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.mgrid[0:10:11j,0:10:11j]\n",
    "X = np.einsum('abc->bca',X).reshape(11**2,2)\n",
    "plt.scatter(X[:,0],X[:,1], alpha = 0.2, c = 'k')\n",
    "\n",
    "\n",
    "\n",
    "o1 = np.linspace(0,10,11)\n",
    "o2 = np.ones_like(o1)*1\n",
    "\n",
    "chan1 = np.concatenate([o1.reshape(-1,1),o2.reshape(-1,1)],axis=1)\n",
    "plt.scatter(chan1[:,0],chan1[:,1], c = 'g')\n",
    "\n",
    "\n",
    "o1 = np.linspace(3,9,7)\n",
    "o2 = o1+1\n",
    "\n",
    "chan2 = np.concatenate([o1.reshape(-1,1),o2.reshape(-1,1)],axis=1)\n",
    "plt.scatter(chan2[:,0],chan2[:,1], c = 'b')\n",
    "\n",
    "o1 = np.ones(5,)*1\n",
    "o2 = np.linspace(2,6,5)\n",
    "chan3 = np.concatenate([o1.reshape(-1,1),o2.reshape(-1,1)],axis=1)\n",
    "plt.scatter(chan3[:,0],chan3[:,1], c = 'r')\n",
    "\n",
    "plt.gcf().set_size_inches(5,5)\n",
    "plt.xlabel('o1')\n",
    "plt.ylabel('o2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "observX = np.concatenate([chan1,chan2,chan3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEzCAYAAABHZATQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAEFpJREFUeJzt3XGIpPV9x/HPt3MjnbtU1+Bi3T1BoTKHWMolQzE5CKLCWJLGK5RgwWAl5f5pownpBDfQaqF/CJuGSCmBxZikRAzhMmxExI1o2pL8cXTWSRnjZZpgot7smdsQ5mLSse5tvv1jZ4/bq3e7M88z88w83/frn9v93c78vrN7+3aefXZ8zN0FANH8TtYDAEAWiB+AkIgfgJCIH4CQiB+AkIgfgJB2jZ+ZPWFmZ8zs5QvW3mtmz5vZj/t/Xj3aMQEgXXt55vdVSXddtPaQpBfc/SZJL/TfB4CpYXv5JWczu0HSM+5+S//9tqTb3P20mV0n6d/cvTzKQQEgTcP+zO9adz/df/tNSdemNA8AjMW+pHfg7m5ml3z6aGbHJB2TpAMHDrz/0KFDSbcEMIFanbOX/Ls/nL9qpHuvrq7+wt1nB7nNsPH7uZldd8Fh75lLfaC7L0lakqRKpeKNRmPILQFMsiOPvqhOt/f/1udnSvr+Q7ePdG8ze23Q2wx72Pu0pPv6b98n6dtD3g+AnKhVyyoVCzvWSsWCatXJPB2w6zM/M3tK0m2SrjGzU5IelvSopG+a2SckvSbpY6McEsDkO3p4XpK0uNLWWrenuZmSatXy+fVJs6ezvWnhsBfAKJjZqrtXBrkNr/AAEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEFLi6/YCmFzLzc7UXFBo3IgfkFPLzY4W6i31NjYlSZ1uTwv1liQRQHHYC+TW4kr7fPi29TY2tbjSzmiiyUL8gJxa6/YGWo+G+AE5NTdTGmg9GuIH5FStWlapWNixVioWVKuWM5posnDCA8ip7ZManO19d8QPyLGjh+eJ3SVw2AsgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgpETX8DCzT0v6K0kuqSXpfnd/O43BgDxabna4oNCEGPqZn5nNS3pAUsXdb5FUkHRPWoMBebPc7Gih3lKn25NL6nR7Wqi3tNzsZD1aSEkPe/dJKpnZPkn7Ja0lHwnIp8WVtnobmzvWehubWlxpZzRRbEPHz907kj4v6XVJpyWddffvXPxxZnbMzBpm1lhfXx9+UmDKrXV7A61jtJIc9l4t6W5JN0qak3TAzO69+OPcfcndK+5emZ2dHX5SYMrNzZQGWsdoJTnsvVPST9193d03JNUlfTCdsYD8qVXLKhULO9ZKxYJq1XJGE8WW5Gzv65JuNbP9knqS7pDUSGUqIIe2z+pytncyDB0/dz9hZsclvSTpnKSmpKW0BgPy6OjheWI3IRL9np+7Pyzp4ZRmAYCx4RUeAEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEJKdA0PYNotNztcTS0o4oewlpsdLdRb6m1sSpI63Z4W6i1JIoABcNiLsBZX2ufDt623sanFlXZGE2GciB/CWuv2BlpHvhA/hDU3UxpoHflC/BBWrVpWqVjYsVYqFlSrljOaCOPECQ+EtX1Sg7O9MRE/hHb08DyxC4rDXgAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACElip+ZzZjZcTP7kZmdNLMPpDUYAIxS0mt4PCbpOXf/czO7QtL+FGZKbLnZ4aI0U4qvHcZl6PiZ2VWSPiTpLyXJ3d+R9E46Yw1vudnRQr2l3samJKnT7Wmh3pIkvokmHF87jFOSw94bJa1L+oqZNc3scTM7kNJcQ1tcaZ//5tnW29jU4ko7o4mwV3ztME5J4rdP0vskfcndD0v6jaSHLv4gMztmZg0za6yvryfYbm/Wur2B1jE5+NphnJLE75SkU+5+ov/+cW3FcAd3X3L3irtXZmdnE2y3N3MzpYHWMTn42mGcho6fu78p6Q0zK/eX7pD0SipTJVCrllUqFnaslYoF1arlS9wCk4KvHcYp6dneT0p6sn+m91VJ9ycfKZntH4xzxnD68LXDOJm7j22zSqXijUZjbPsBiMHMVt29MshteIUHgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCSXsBoIi03O1wEJyV8LpFXuYvfcrOjhXpLvY1NSVKn29NCvSVJfNMOiM8l8ix3h72LK+3z36zbehubWlxpZzTR9OJziTzLXfzWur2B1nFpfC6RZ7mL39xMaaB1XBqfS+RZ7uJXq5ZVKhZ2rJWKBdWq5Ywmml58LpFnuTvhsf2DeM5QJsfnEnlm7j62zSqVijcajbHtByAGM1t198ogt8ndYS8A7AXxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQEvEDEBLxAxAS8QMQUuL4mVnBzJpm9kwaAwHAOKRxAaMHJZ2UdGUK95WK5WYn1xfdyfvjA8Yh0TM/Mzso6cOSHk9nnOSWmx0t1FvqdHtySZ1uTwv1lpabnaxHS0XeHx8wLkkPe78o6bOSfpvCLKlYXGmrt7G5Y623sanFlXZGE6Ur748PGJeh42dmH5F0xt1Xd/m4Y2bWMLPG+vr6sNvt2Vq3N9D6tMn74wPGJckzvyOSPmpmP5P0DUm3m9nXL/4gd19y94q7V2ZnZxNstzdzM6WB1qdN3h8fMC5Dx8/dF9z9oLvfIOkeSS+6+72pTTakWrWsUrGwY61ULKhWLWc0Ubry/viAcUnjbO9E2T7rmdezoXl/fMC4mLuPbbNKpeKNRmNs+wGIwcxW3b0yyG14hQeAkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkIgfgJCIH4CQiB+AkHJ3ASNJWm52uMAPgMvKXfyWmx0t1FvqbWxKkjrdnhbqLUkigADOy91h7+JK+3z4tvU2NrW40s5oIgCTKHfxW+v2BloHEFPu4jc3UxpoHUBMuYtfrVpWqVjYsVYqFlSrljOaCMAkyt0Jj+2TGpztBXA5uYuftBVAYgfgcnJ32AsAe0H8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhET8AIRE/ACERPwAhDT0NTzM7HpJ/yrpWkkuacndH0trsCSWmx0uYATgspJcwOicpM+4+0tm9nuSVs3seXd/JaXZhrLc7Gih3lJvY1OS1On2tFBvSRIBBHDe0Ie97n7a3V/qv/2WpJOSMq/L4kr7fPi29TY2tbjSzmgiAJMolZ/5mdkNkg5LOvEuf3fMzBpm1lhfX09ju8ta6/YGWgcQU+L4mdl7JH1L0qfc/VcX/727L7l7xd0rs7OzSbfb1dxMaaB1ADElip+ZFbUVvifdvZ7OSMnUqmWVioUda6ViQbVqOaOJAEyiJGd7TdKXJZ109y+kN1Iy2yc1ONsL4HKSnO09Iunjklpm9oP+2ufc/dnkYyVz9PA8sQNwWUPHz92/J8lSnAUAxoZXeAAIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIKckFjAbW6pzVkUdfHPnV1JabnbFevS3P++X5sbFffva74vf/4P2D3rbwyCOPjGCkd/eP//TPj/ihO/Xv/72ug1eXdOi6K1PfY7nZ0UK9pV/+zzuSpLfePsd+U7AX+7Ffkv1+/V/P6e/+9oF/GOT2mRz29jY2tbjSHsl9L6601dvYZL8p24v92C+N/QaR2c/81rq9sd4v+032XuzHfmntt1eZxW9upjTW+2W/yd6L/dgvrf32KpP4lYoF1arlkdx3rVpWqVhgvynbi/3YL439BjH2Ex6Hbvsz/f2f3jyyM0CHrrtSB68uqdU5q1+/fU7zMyX2m4K92I/9kuy3duKZgU94mLuPZLB3U6lUvNFojG0/ADGY2aq7Vwa5Db/kDCAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHICTiByAk4gcgJOIHIKRE8TOzu8ysbWY/MbOH0hoKAEZt6PiZWUHSv0j6E0k3S/oLM7s5rcEAYJSSPPP7Y0k/cfdX3f0dSd+QdHc6YwHAaCWJ37ykNy54/1R/DQAm3r5Rb2BmxyQd67/7v2b28qj3zNA1kn6R9RAjkufHJvH4pt3AV0ZPEr+OpOsveP9gf20Hd1+StCRJZtYY9Nqa0yTPjy/Pj03i8U07Mxv4guBJDnv/U9JNZnajmV0h6R5JTye4PwAYm6Gf+bn7OTP7G0krkgqSnnD3H6Y2GQCMUKKf+bn7s5KeHeAmS0n2mwJ5fnx5fmwSj2/aDfz4zN1HMQgATDRe3gYgpLHEL88vgzOz683su2b2ipn90MwezHqmUTCzgpk1zeyZrGdJm5nNmNlxM/uRmZ00sw9kPVOazOzT/X+bL5vZU2b2u1nPlISZPWFmZy78tTkze6+ZPW9mP+7/efVu9zPy+AV4Gdw5SZ9x95sl3Srpr3P2+LY9KOlk1kOMyGOSnnP3Q5L+SDl6nGY2L+kBSRV3v0VbJyfvyXaqxL4q6a6L1h6S9IK73yTphf77lzWOZ365fhmcu59295f6b7+lrW+cXL3SxcwOSvqwpMezniVtZnaVpA9J+rIkufs77t7NdqrU7ZNUMrN9kvZLWst4nkTc/T8k/fKi5bslfa3/9tckHd3tfsYRvzAvgzOzGyQdlnQi20lS90VJn5X026wHGYEbJa1L+kr/sP5xMzuQ9VBpcfeOpM9Lel3SaUln3f072U41Ete6++n+229Kuna3G3DCIyVm9h5J35L0KXf/VdbzpMXMPiLpjLuvZj3LiOyT9D5JX3L3w5J+oz0cMk2L/s++7tZW5OckHTCze7OdarR861dYdv01lnHEb08vg5tmZlbUVviedPd61vOk7Iikj5rZz7T1I4vbzezr2Y6UqlOSTrn79rP149qKYV7cKemn7r7u7huS6pI+mPFMo/BzM7tOkvp/ntntBuOIX65fBmdmpq2fF5109y9kPU/a3H3B3Q+6+w3a+tq96O65eebg7m9KesPMtl8Yf4ekVzIcKW2vS7rVzPb3/63eoRyd0LnA05Lu6799n6Rv73aDkf9fXQK8DO6IpI9LapnZD/prn+u/+gXT4ZOSnuz/x/lVSfdnPE9q3P2EmR2X9JK2fjOhqSl/tYeZPSXpNknXmNkpSQ9LelTSN83sE5Jek/SxXe+HV3gAiIgTHgBCIn4AQiJ+AEIifgBCIn4AQiJ+AEIifgBCIn4AQvo/Bged3y05j8YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(observX[:,0],observX[:,1])\n",
    "plt.xlim(0,10)\n",
    "plt.ylim(0,10)\n",
    "plt.gcf().set_size_inches(5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x113fec438>]"
      ]
     },
     "execution_count": 250,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl4VPW9x/H3l4SELECAQICEEDZB9iUsSrWtSy9qFZfWiooICNrFpbZ16+1je629tre1SrcrBRSQRVwqXosLoi2uQCJI2AlhCRCyEEhiQsgyv/tHRh9KUSCTyZmZfF7Pk2dmTmbmfAbCh5PfOed3zDmHiIhErlZeBxARkeBS0YuIRDgVvYhIhFPRi4hEOBW9iEiEU9GLiEQ4Fb2ISIRT0YuIRDgVvYhIhIv2OgBAcnKyy8jI8DqGiEhYyc7OLnHOdT7d80Ki6DMyMsjKyvI6hohIWDGzvWfyPA3diIhEOBW9iEiEU9GLiEQ4Fb2ISIRT0YuIRLjTFr2ZzTOzIjPbdMKyjma20sx2+m87+Jebmc0ys1wz22hmI4MZXkRETu9MtuifASactOwBYJVzrh+wyv8Y4DKgn/9rJvCXpokpIiKNddqid86tBkpPWjwRmO+/Px+4+oTlC1yDj4AkM+vWVGFFRCKFz+d49O9byC+tCvq6GjtGn+KcK/DfPwSk+O+nAvknPG+/f9m/MbOZZpZlZlnFxcWNjCEiEp5mvb2Tv767m3d3lgR9XQHvjHUNVxc/6yuMO+dmO+cynXOZnTuf9gxeEZGI8c72Ip5ctZNrR6QyaUyPoK+vsUVf+NmQjP+2yL/8AHBi6jT/MhERAfJLq7hn6Qb6p7Tl0WuGYGZBX2dji/4VYIr//hRg+QnLb/EffTMOKDthiEdEpEWrrq3nu4uy8TnHU5NHERcT1SzrPe2kZma2BPgakGxm+4GHgceAZWY2HdgLXO9/+grgciAXqAKmBiGziEhYenj5ZjYdKGfOLZn07JTQbOs9bdE75yZ9wbcuPsVzHfD9QEOJiESa59bt47msfL7/9T5cMjDl9C9oQjozVkQkyHL2l/Gz5Zv5St9k7r20f7OvX0UvIhJER6tq+O6ibJITYpg1aQRRrYK/8/VkIXHhERGRSFTvc9y5ZD1F5cdZdsd5dEyI8SSHil5EJEh+9+Z23t1Zwn9fO4ThPZI8y6GhGxGRIHgtp4A//2MXk8akM2lMuqdZVPQiIk1sZ2EFP37+E0akJ/HzqwZ6HUdFLyLSlMqra5m5MJu4mGj+ctMoYqOb56SoL6OiFxFpIj6f44dLN5BfWsWfbxpJ1/ZtvI4EqOhFRJrMrLd3smpbET/75kDG9OrodZzPqehFRJrAqq2FPPHWTq4dmcot5/X0Os6/UNGLiAQot6iCu5duYHBqO37VTDNSng0VvYhIAMqO1TJjQTZtWrfiqcmZtGnt/c7Xk+mEKRGRRqr3Oe5asp79R6pYPGMcqUlxXkc6JRW9iEgj/eaNbfxzRzGPXjOY0Rmhs/P1ZBq6ERFphOUbDvDUP/O4aWw6N40NrZ2vJ1PRi4icpZz9Zdz3wkbGZHTk4SsHeR3ntFT0IiJnobjiODMXZtEpIYY/3zySmOjQr1GN0YuInKHjdfV8b1E2R6pqeOGO80lOjPU60hlR0YuInAHnHP/5t02s23OEWZNGMDi1vdeRzljo/84hIhIC5r63m+ez93PXRX25alh3r+OcFRW9iMhpvLO9iF+t2MqEQV2555JzvI5z1lT0IiJfIreogrsWr2dA13Y8/p1htPLgmq+BUtGLiHyBI5U1TJ+fRWzrKP46JZP4mPDcramiFxE5hdp6H99b9DEFR6t5avKokJ3e4EyE539PIiJB5Jzj569s5sO8w/zu28MY1bOD15ECoi16EZGTPPPBHhat2cftX+3NdaPSvI4TMBW9iMgJ3t5WyCOvbuEbA1O4/z8GeB2nSajoRUT8thaUc+fi9Qzs3o4nbhgelkfYnIqKXkQEKKqoZvoz60hsE82cW0aH7RE2pxI5n0REpJGqa+uZsSCbI1W1PH/HeXRt38brSE0qoC16M/uhmW02s01mtsTM2phZLzNbY2a5ZvacmcU0VVgRkabm8zl+tOwTNu4/ypM3DA+rOWzOVKOL3sxSgbuATOfcYCAKuAH4NfB751xf4AgwvSmCiogEw+Mrd/D3nAIevGwA3xjU1es4QRHoGH00EGdm0UA8UABcBLzg//584OoA1yEiEhTPZ+Xzx3dyuWF0D2Zc0NvrOEHT6KJ3zh0Afgvso6Hgy4Bs4Khzrs7/tP1AaqAhRUSa2vu5JTz4Ug5f6ZvMI1cPxiwyjrA5lUCGbjoAE4FeQHcgAZhwFq+faWZZZpZVXFzc2BgiImdt+6EK7liYTZ/Oifz55pG0jorsAxAD+XSXALudc8XOuVrgJWA8kOQfygFIAw6c6sXOudnOuUznXGbnzp0DiCEicuaKyquZ9sw64mKimDd1NO3atPY6UtAFUvT7gHFmFm8Nv/NcDGwB3gG+5X/OFGB5YBFFRJpG5fE6ps1fx5GqGubdOjqsJyo7G4GM0a+hYafrx0CO/71mA/cD95pZLtAJmNsEOUVEAlLvc9y1ZD1bDpbzxxvD61KAgQrohCnn3MPAwyctzgPGBPK+IiJNyTnHL/5vM6u2FfHIxEFcNCDF60jNKrL3QIiIAHPe3c2CD/cy88LeTD4vw+s4zU5FLyIR7ZVPDvLoiq1cMaQbD0yIjNkoz5aKXkQi1ge7SvjRsg2M6dWR310fntd7bQoqehGJSNsOlXP7gmwyOiXw18mZtGkd5XUkz6joRSTiHDx6jFvnrSM+Nor508bQPj7yj5X/Mip6EYkoZcdqufXptVQer+OZqWPo3kKOlf8ymo9eRCLG8bp6Zi7IYndJJfOnjuHcbu28jhQSVPQiEhHqfY57n/uENbtLefKG4ZzfN9nrSCFDQzciEvY+OyHq7zkF/PTyc5k4XJPmnkhFLyJh7w9v57Lgw73cfmFvZlwYufPKN5aKXkTC2qI1e3l85Q6uHZnK/S30hKjTUdGLSNh6fVMBP3t5ExcN6MKvrxvaYk+IOh0VvYiEpY/yDnPX0g0M75HEn26M/IuHBEJ/MiISdjYfLGPG/CzSO8Yz79bRxMW03LNez4SKXkTCyp6SSqbMW0dim2gWTBtDUnyM15FCnopeRMLGobJqbpqzBp9zLJyus17PlIpeRMJCaWUNN89dQ9mxWuZPHUPfLm29jhQ2dGasiIS8T4/XMfXptewrrWLBtDEMSWs5lwFsCtqiF5GQVl1bz4z5WWw6WM6fbxzJuN6dvI4UdlT0IhKy6up93LlkPR/mHea33x7KJQNb1rVem4qKXkRCks/nuO/FjazcUsgvrhrENSPSvI4UtlT0IhJynHM8/MpmXvr4APdeeg5Tzs/wOlJYU9GLSEhxzvHYa9tY+FHDJGV3XtTX60hhT0UvIiFl1qpcnlqdx+RxPXngsgGYaf6aQKnoRSRk/HV1Hr9/awffGpXGL64apJJvIip6EQkJCz/ay6MrtnLF0G6aibKJqehFxHMvZu/nZy9v4uIBXfj99cOJUsk3KRW9iHjq1Y0H+ckLnzC+byf+dNNIYqJVS01Nf6Ii4pnXNxVw99INjOrZgb/ekkmb1ppuOBhU9CLiiZVbCvnB4vUMS2vP01PHEB+jqbeCRUUvIs3unW1FfG9RNoNS2/PMtDEkxqrkgymgojezJDN7wcy2mdlWMzvPzDqa2Uoz2+m/7dBUYUUk/K3eUcztz2bTv2tbFkwbQ7s2rb2OFPEC3aJ/EnjdOTcAGAZsBR4AVjnn+gGr/I9FRHg/t4QZC7Lo0zmRZ6ePpX2cSr45NLrozaw9cCEwF8A5V+OcOwpMBOb7nzYfuDrQkCIS/j7cdZjp89eR0SmBRbeN1SUAm1EgW/S9gGLgaTNbb2ZzzCwBSHHOFfifcwg45byiZjbTzLLMLKu4uDiAGCIS6j7YVcK0Z9aR1iGeZ28bS8cElXxzCqToo4GRwF+ccyOASk4apnHOOcCd6sXOudnOuUznXGbnzp0DiCEioeyD3M9KPo4lM8bRuW2s15FanECKfj+w3zm3xv/4BRqKv9DMugH4b4sCiygi4er93BKmzV9Hesd4lsxUyXul0UXvnDsE5JtZf/+ii4EtwCvAFP+yKcDygBKKSFh6b2fDlnzPjgksnjGO5ESVvFcCPXj1TmCRmcUAecBUGv7zWGZm04G9wPUBrkNEwsy7O4u5bX4WvZIbdrx2Usl7KqCid85tADJP8a2LA3lfEQlfq3cUM2NBQ8kvnjFOO15DgM6MFZEms2prIbfNz6J350SVfAhR0YtIk3gtp4DbF2YzoFtblszQIZShRBNMiEjAlm84wL3LPmF4jySenjpa0xqEGG3Ri0hAlmXlc89zGxid0UFz14QobdGLSKM9+9Fe/vPlTVzQL5nZkzOJi9F88qFIRS8ijTL3vd088uoWLh7QhT/dNFIXDQlhKnoROSvOOWatyuX3b+1gwqCuzJo0Qpf/C3EqehE5Y845Hv37Vua8t5trR6bym+uGEh2lkg91KnoROSP1PsdP/5bD0nX5TDmvJw9fOYhWrczrWHIGVPQiclo1dT7uXbaBVzcWcOdFfbn30nMwU8mHCxW9iHyp6tp6vvtsNu9sL+ahywcw88I+XkeSs6SiF5EvVF5dy4z5WazdU8qvrhnCjWPTvY4kjaCiF5FTKq44zpR5a9lRWMET3xnOxOGpXkeSRlLRi8i/2Xe4isnz1lBUfpw5UzL5Wv8uXkeSAKjoReRfbDlYzpSn11Jb72PxjLGMSO/gdSQJkIpeRD63Ju8wt83PIrFNNItvO49+KW29jiRNQEUvIgC8ufkQP1iynrQOcSycPpbUpDivI0kTUdGLCEvX7uOhv+UwJLU9T08do7nkI4yKXqQFc87x+7d2MmvVTi48pzN/uWkkCbGqhUijv1GRFqq23sdP/5bDsqz9fHtUGr+6dgitNW9NRFLRi7RAlcfr+P7ij/nH9mLuuqgvP9SUBhFNRS/SwhRXHGfaM+vYfLBMZ7u2ECp6kRZkd0klU+atpaiimtmTM7lkYIrXkaQZqOhFWoi1u0uZuTCLVmYsmTFOJ0K1ICp6kRZg+YYD/OT5jaR1jOPpW0fTs1OC15GkGanoRSKYc44/vp3L71buYGyvjjw1eRRJ8TpGvqVR0YtEqJo6Hw++lMOLH+/n2hGp/Pd1Q4iN1gW8WyIVvUgEKquq5Y5ns/kw7zD3XNKPuy/up8MnWzAVvUiE2VNSybT568gvreLx64dx7cg0ryOJx1T0IhHkg10lfPfZj2llsHD6WMb17uR1JAkBAZ/vbGZRZrbezF71P+5lZmvMLNfMnjMz7fkRaQaL1+zjlrlr6dw2lpe/P14lL59riokt7ga2nvD418DvnXN9gSPA9CZYh4h8gbp6H7/4v8089LccxvdN5qXvna/DJ+VfBFT0ZpYGXAHM8T824CLgBf9T5gNXB7IOEfli5dW13LYgi6ff38O08b2YOyWTdm1aex1LQkygY/RPAPcBn12GphNw1DlX53+8H9AVhUWCYHdJJTMWZLGnpFJz1siXanTRm9k3gSLnXLaZfa0Rr58JzARIT9cPqMjZ+Mf2Iu5asp6oVsaC6WM4v0+y15EkhAWyRT8euMrMLgfaAO2AJ4EkM4v2b9WnAQdO9WLn3GxgNkBmZqYLIIdIi+Gc46nVefzm9W3079qO2ZNH0aNjvNexJMQ1eozeOfegcy7NOZcB3AC87Zy7CXgH+Jb/aVOA5QGnFBGO1dRz99INPPbaNi4f0o0Xv3ueSl7OSDCOo78fWGpmvwTWA3ODsA6RFmX/kSpuX5jNloJy7pvQn+9+tY/OdJUz1iRF75z7B/AP//08YExTvK+IwAe5JfxgyXpq633MmzKarw/o4nUkCTM6M1YkRJ04Ht+7cyJPTR5Fn86JXseSMKSiFwlBFdW1/Pj5T3hjcyFXDO3Gb64bSkKs/rlK4+gnRyTE7Cis4I6F2ewtreI/rziX6V/ppfF4CYiKXiSE/N8nB7n/xY3Ex0Sz+LaxjNV8NdIEVPQiIaCmzsevVmzlmQ/2MKpnB/5800hS2rXxOpZECBW9iMfyS6v4weKP+WR/GVPHZ/DgZecSE90U8w2KNFDRi3ho5ZZCfrRsA87B/948kgmDu3kdSSKQil7EA7X1Pv7nje3MXp3H4NR2/OnGkZpaWIJGRS/SzA4ePcadS9aTvfcIk8f15KdXnEub1rpotwSPil6kGb2x+RD3v7iR2joff5g0giuHdfc6krQAKnqRZlBdW8+jf9/Kwo/2MiS1PbMmjaBXsoZqpHmo6EWCbGdhBXcuWc+2QxXMuKAXP/mPATqqRpqVil4kSJxzLFmbz3+9upmEmGienjqar/fXhGTS/FT0IkFwpLKGn76cw4qcQ3ylbzKPf2cYXdrqBCjxhopepIm9u7OYHz//CaWVNdw/YQC3X9ibVq00V414R0Uv0kSqa+v5zevbmff+bvp2SWTulNEMTm3vdSwRFb1IU9hysJx7nlvPjsJPufX8DB64bICOjZeQoaIXCUC9zzH3vTx++8YO2se35pmpo/madrhKiFHRizTSnpJKfvz8J2TtPcJ/DErhv68dSseEGK9jifwbFb3IWfL5HAs/2stjr20jOsp4/PphXDMiVRcHkZClohc5C/mlVdz/4kY+2HWYr57TmV9fN5Su7XXYpIQ2Fb3IGXDOsXRdPr98dQsAj107hO+M7qGteAkLKnqR08gvreLBl3J4L7eE83p34jffGkqPjvFexxI5Yyp6kS9Q73PM/2AP//PGdqJaGY9cPZibxqTr5CcJOyp6kVPYWVjBfS9uZP2+o3y9f2cevWYI3ZPivI4l0igqepET1NT5+N9/7uKPb+eSEBvFE98ZzsTh3TUWL2FNRS/it25PKQ+9lMPOok+5clh3Hr5yIMmJsV7HEgmYil5avKNVNTz22jaWrssnNSmOuVMyufjcFK9jiTQZFb20WM45Xt5wgF++upWjx2qZeWFv7rmkH/Ex+mchkUU/0dIi5RV/ys+Wb+L93MMM75HEwmuGMLB7O69jiQSFil5alKqaOv7wdi5z3s2jTesofnn1YG7UIZMS4Rpd9GbWA1gApAAOmO2ce9LMOgLPARnAHuB659yRwKOKNJ5zjhU5h/jl37dQUFbNdSPTuP+y/rrqk7QIgWzR1wE/cs59bGZtgWwzWwncCqxyzj1mZg8ADwD3Bx5VpHFyiyp4+JXNvJ97mIHd2vGHSSPIzOjodSyRZtPoonfOFQAF/vsVZrYVSAUmAl/zP20+8A9U9OKB8upa/vh2LvPe2018TBT/NXEQN45JJzqqldfRRJpVk4zRm1kGMAJYA6T4/xMAOETD0M6pXjMTmAmQnp7eFDFEAKir9/FcVj6Pv7mD0qoavj0qjfsmDNAx8dJiBVz0ZpYIvAjc45wrP/EMQuecMzN3qtc552YDswEyMzNP+RyRs/V+bgmPvLqFbYcqGJPRkflXDtR1W6XFC6jozaw1DSW/yDn3kn9xoZl1c84VmFk3oCjQkCKnk1f8Kb9asY23thbSo2Mcf7lpJBMGd9XUBSIEdtSNAXOBrc65x0/41ivAFOAx/+3ygBKKfImST48za9VOFq/ZR5vWUdw/YQBTx2fowtwiJwhki348MBnIMbMN/mUP0VDwy8xsOrAXuD6wiCL/rvJ4HXPe3c3s1buorvNxw+ge3H1JPx0uKXIKgRx18x7wRb8XX9zY9xX5MrX1Pp5bl88Tb+2k5NPjTBjUlZ9M6E+fzoleRxMJWTozVsKCz+dYsamAx1fuIK+4ksyeHXhq8khG9dTx8CKno6KXkOacY9XWIn63cgdbC8rp1yWR2ZNHcenAFO1oFTlDKnoJSc453s89zG/f3M6G/KP07BTPE98ZzpXDuhOleWlEzoqKXkLOmrzDPL5yB2t2l9K9fRseu3YI141Ko7XOaBVpFBW9hATnHB/sOsyTq3aydncpyYmx/PzKgUwam05stA6VFAmEil485ZzjnzuKmbVqJx/vO0pKu1gevnIgk8ak61h4kSaiohdP+HyOt7YW8sd3ctm4v4zUpDgeuXow3x6VpoIXaWIqemlWx+vqWb7+IE+t3sWu4krSO8bz6+uGcM2INGKiNQYvEgwqemkWFdW1LFm7j7nv7aaw/DgDu7Vj1qQRXD64q6YNFgkyFb0EVUHZMeZ/sJdFa/ZSUV3H+X068T/fGsYF/ZJ1HLxIM1HRS1B8vO8I897bzWubDuGcY8Lgrtzx1T4MTUvyOppIi6OilyZTW+/jtU2HmPfebjbkH6VtbDTTxmdwy3kZ9OgY73U8kRZLRS8BKyyvZunafJas3ceh8moyOsXzi6sGcd2oNBJj9SMm4jX9K5RG+ewEp2c/2subWwqp9zku6JfMo9cM5uv9u9BK0xSIhAwVvZyVI5U1vLT+AIvW7CWvuJKk+NZM/0ovbhyTTkZygtfxROQUVPRyWvU+x/u5JTyXlc/KzYXU1PsYkZ7E7749jCuGdtMJTiIhTkUvXyi/tIrns/fzQlY+B8uqSYpvzY1j07k+swcDu7fzOp6InCEVvfyLsmO1vJZTwN/WH2DN7lLM4IJ+nXnoinO5dGCKJhgTCUMqeuF4XT3vbCvi5fUHeXtbETX1PnonJ3Dvpedw3ag0UpPivI4oIgFQ0bdQtfU+3s8tYUVOAa9vOkR5dR3JibHcPK4nV4/ozpDU9jpzVSRCqOhbkJo6H+/vKmHFxgLe3FJI2bFaEmOj+cbAFCaOSGV8n06ad0YkAqnoI1zl8Tre3VnCyi2FrNzSsOXeNjaaSwemcNmQblzQL1lHzYhEOBV9BDpUVs1bWwt5a2shH+w6TE2dj3ZtorlkYApXDOnGV/ola6eqSAuioo8AtfU+1u87yuodxfxzRzE5B8oASO8Yz+RxPbnk3BQyMzromqsiLZSKPkzll1axemcx/9xezIe7DlNxvI6oVsbwHkncN6E/l56bQt8uidqhKiIq+nBRUHaMj/IO8+Guw3yYd5j80mMApCbF8c1h3biwX2fO75tM+7jWHicVkVCjog9Bzjn2Hq4ie+8R1u0p5aO8w+w5XAVAuzbRjO3diann9+LCc5Lp01lb7SLy5VT0IeBYTT05B8rI3nuE7L1HWL/vCIcrawBo2yaasb06cfO4npzXpxMDurYjSjNDishZUNE3s6qaOrYcLCfnQBk5B8rYfKCcnUUV+FzD93snJ/D1AV0Ymd6BUT070K9Loqb8FZGAqOiDpN7n2FdaxfZDFeworGB7YQXbD1WQV/zp56WenBjLkNR2fGNQCkPTkhiZnkSnxFhvg4tIxAlK0ZvZBOBJIAqY45x7LBjrCQVHq2rIK6lkj/8rr6SS3SWV5BZ9yvE6HwBm0KNDPOektOXyId0YktqeIantSWkXq/F1EQm6Ji96M4sC/gRcCuwH1pnZK865LU29rmCrq/dRWlnDofJqDhw5xoGjDV8H/bf7jxzjaFXt589vZZDWIZ6M5ATO79OJc1La0r9rW/p2SSQ+Rr88iYg3gtE+Y4Bc51wegJktBSYCnhS9c47jdT6O1dRTVVvPsZp6yqtrKTtWS/mxhtujVQ23pZU1lHx6nOKKhq/Sqhqc+9f3i4+JIjUpju5JcQxNS6J3cgIZnRLISE4gvWM8MdE6KUlEQkswij4VyD/h8X5gbBDWw7J1+Ty1ehc+B3U+Hz5fw9h4nc9R7/NRXevjWG39Gb1XfEwUHRNiSE6MpUfHeEb27EDnxFiS28bSpW0sqUlxpHWIo31caw23iEhY8Ww8wcxmAjMB0tPTG/UeHRJiPj/c8PMvM6KiGm7jYqJo0zqK+Jgo4lo3fMW2bkW7uNa0P+GrXZvW2hIXkYgVjKI/APQ44XGaf9m/cM7NBmYDZGZmupO/fyYuHZjCpQNTGvNSEZEWIxibseuAfmbWy8xigBuAV4KwHhEROQNNvkXvnKszsx8Ab9BweOU859zmpl6PiIicmaCM0TvnVgArgvHeIiJydrQHUkQkwqnoRUQinIpeRCTCqehFRCKcil5EJMKZO3kyFy9CmBUDexv58mSgpAnjhAN95pZBn7llCOQz93TOdT7dk0Ki6ANhZlnOuUyvczQnfeaWQZ+5ZWiOz6yhGxGRCKeiFxGJcJFQ9LO9DuABfeaWQZ+5ZQj6Zw77MXoREflykbBFLyIiXyKsi97MJpjZdjPLNbMHvM4TbGbWw8zeMbMtZrbZzO72OlNzMLMoM1tvZq96naU5mFmSmb1gZtvMbKuZned1pmAzsx/6f6Y3mdkSM2vjdaamZmbzzKzIzDadsKyjma00s53+2w7BWHfYFv0JFyG/DBgITDKzgd6mCro64EfOuYHAOOD7LeAzA9wNbPU6RDN6EnjdOTcAGEaEf3YzSwXuAjKdc4NpmN78Bm9TBcUzwISTlj0ArHLO9QNW+R83ubAtek64CLlzrgb47CLkEcs5V+Cc+9h/v4KGAkj1NlVwmVkacAUwx+sszcHM2gMXAnMBnHM1zrmj3qZqFtFAnJlFA/HAQY/zNDnn3Gqg9KTFE4H5/vvzgauDse5wLvpTXYQ8okvvRGaWAYwA1nibJOieAO4DfF4HaSa9gGLgaf9w1RwzS/A6VDA55w4AvwX2AQVAmXPuTW9TNZsU51yB//4hICjXRg3nom+xzCwReBG4xzlX7nWeYDGzbwJFzrlsr7M0o2hgJPAX59wIoJIg/TofKvzj0hNp+E+uO5BgZjd7m6r5uYZDIINyGGQ4F/0ZXYQ80phZaxpKfpFz7iWv8wTZeOAqM9tDw9DcRWb2rLeRgm4/sN8599lvai/QUPyR7BJgt3Ou2DlXC7wEnO9xpuZSaGbdAPy3RcFYSTgXfYu7CLmZGQ1jt1udc497nSfYnHMPOufSnHMZNPz9vu2ci+gtPefcISDfzPr7F10MbPEwUnPYB4wzs3j/z/jFRPgO6BO8Akzx358CLA/GSoJyzdjm0EIvQj4emAzkmNkG/7KH/NfolchxJ7DIvwGTB0z1OE9QOefWmNkLwMc0HFm2ngg8Q9bMlgBfA5LNbD/wMPAYsMzMptMwg+/1QVm3zoz2WX1qAAAAOElEQVQVEYls4Tx0IyIiZ0BFLyIS4VT0IiIRTkUvIhLhVPQiIhFORS8iEuFU9CIiEU5FLyIS4f4f4YsPilZdPtAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,10)\n",
    "y = x**2\n",
    "plt.plot(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "def truth(X,alphas):\n",
    "    al0, al1 = alphas[:,0], alphas[:,1]\n",
    "    x0,x1 = X[:,0], X[:,1]\n",
    "\n",
    "    norm1 = 1.5-al1\n",
    "    norm2 = 3.0-norm1\n",
    "    f0 = 3-al0\n",
    "    f1 = 3-al1\n",
    "    len0 = 5*f0\n",
    "    len1 = 7*f1\n",
    "    \n",
    "    return 20 * (norm1 * np.exp(-x0/len0) + norm2 * np.exp(-x1/len1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 252,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sign(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.random.randint(len())\n",
    "# observX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22]\n"
     ]
    }
   ],
   "source": [
    "global_chan_numbers = np.arange(len(observX))\n",
    "print(global_chan_numbers)\n",
    "\n",
    "ind = [0] + np.cumsum([len(chan1),len(chan2),len(chan3)]).tolist()\n",
    "sli = [slice(ind[i],ind[i+1]) for i in range(len(ind)-1)]\n",
    "\n",
    "\n",
    "meas_chan = [0,1,2]\n",
    "\n",
    "chan_datas = []\n",
    "for chan_num in meas_chan:\n",
    "    global_chan_numbers[sli[1]]\n",
    "    n_evals = 20\n",
    "    bin_nums  = np.random.choice(global_chan_numbers[sli[chan_num]], size = n_evals)\n",
    "    alphas    = np.random.uniform(-1.3,1.3, size = (n_evals,2))\n",
    "    chan_data = np.concatenate([bin_nums.reshape(-1,1),alphas], axis=1)\n",
    "    chan_datas.append(chan_data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x116954eb8>"
      ]
     },
     "execution_count": 255,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAD+dJREFUeJzt3W+IXXedx/HPJzNddCokIxmKSZu5eVCUGcWte1mqBZHOPOhu1fpIKlMty8JAp6utCFI3sKWFLPtApD7oBC5ttcVLRWJZu3ZxNaMi+6TspC22SRTFNtM0aTPFporzoG387oNz02RmMn/u3HPu797fvF9wOfd8c+b+voc5+eTknHPPcUQIAND/dqRuAABQDgIdADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkInBbg62e/fuqNVq3RwSAPre0aNHX4+IkY2W62qg12o1zc/Pd3NIAOh7tk9uZjkOuQBAJgh0AMgEgQ4AmSDQASATBDoAZGLDQLf9iO2ztl+4pPZ+2z+z/bvWdLjaNgGgA82mVKtJO3YU02YzdUeV2Mwe+ncl3bSido+kuYi4VtJcax4Aek+zKU1PSydPShHFdHo6y1DfMNAj4leS/riifIukR1vvH5X0uZL7AoByHDggLS0try0tFfXMbPUY+lURcab1/lVJV621oO1p2/O25xcXF7c4HABs0cJCe/U+1vFJ0SieMr3mk6YjohER9Yioj4xs+M1VACjXvn3t1fvYVgP9NdsfkKTW9Gx5LQFAiQ4elIaGlteGhop6ZrYa6E9Kur31/nZJPyqnHQAo2dSU1GhIo6OSXUwbjaKeGRdHTNZZwH5c0qck7Zb0mqR7Jf2npB9I2ifppKTPR8TKE6er1Ov14OZcANAe20cjor7RchvebTEivrDGH0203RUAoDJ8UxQAMkGgA0AmCHQAyASBDgCZINABIBMEOgBkgkAHgEwQ6ACQCQIdADJBoANAJgh0AMgEgQ4AmSDQASATBDqA7mg2pVpN2rGjmGb4kObUNrx9LgB0rNmUpqcvPqz55MliXsryQROpsIcOoHoHDlwM8wuWloo6SkOgA6jewkJ7dWwJgQ6gevv2tVfHlhDoAKp38KA0NLS8NjRU1FEaAh1A9aampEZDGh2V7GLaaHBCtGRc5QKgO6amCPCKsYcOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBkgkAHgEwQ6ACQiY4C3fZXbR+z/YLtx22/p6zGAADt2XKg294r6SuS6hHxYUkDkm4tqzEAFRkeLu5JfuE1PJy6I5Sk00Mug5Lea3tQ0pCk0523BKAyw8PSuXPLa+fOEeqZ2HKgR8Qrkr4paUHSGUlvRsRPy2oMQAVWhvlGdfSVTg65DEu6RdJ+SXskXWn7tsssN2173vb84uLi1jsFAKyrk0Muk5JejIjFiHhb0hOSPrFyoYhoREQ9IuojIyMdDAcAWE8ngb4g6XrbQ7YtaULSiXLaAlCJXbvaq6OvdHIM/WlJhyU9I+n51mc1SuoLQBXeeGN1eO/aVdTR9wY7+eGIuFfSvSX1AqAbCO9s8U1RAMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh1IZWZGGhyU7GI6M5O6I/S5ju6HDmCLZmakQ4cuzp8/f3F+djZNT+h77KEDKTTWeLjXWnVgEwh0IIXz59urA5tAoAMpDAy0Vwc2gUAHUpiebq8ObAInRYEULpz4bDSKwywDA0WYc0IUHSDQgVRmZwlwlIpDLgCQCQIdADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJnoKNBt77J92PZvbJ+w/fGyGgMAtKfTPfRvS/pJRHxI0kclnei8pd4xOVk8kP3Ca3IydUco3fDw8l/y8HDqjoAt23Kg294p6ZOSHpakiHgrIs6V1Vhqk5PS3Nzy2twcoZ6V4WHp3IpN9tw5Qh19q5M99P2SFiV9x/azth+yfWVJfSW3Msw3qqMPrQzzjepAj+sk0AclfUzSoYi4TtJfJN2zciHb07bnbc8vLi52MBwAYD2dBPopSaci4unW/GEVAb9MRDQioh4R9ZGRkQ6GAwCsZ8uBHhGvSnrZ9gdbpQlJx0vpqgdMTLRXRx/atau9OtDjOr3K5cuSmrZ/LelvJf175y31hiNHVof3xERRRybeeGN1eO/aVdSBPjTYyQ9HxHOS6iX10nMI722A8EZG+KYoAGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBkgkAHgEwQ6ACQCQIdADJBoANAJgh0AMgEgb4OHgjfRc2mVKtJO3YU02YzdUdA3+nofug5W++B8NxCu2TNpjQ9LS0tFfMnTxbzkjQ1la4voM+wh74GHgjfRQcOXAzzC5aWijqATSPQkd7CQnt1AJdFoCO9ffvaqwO4LAJ9DTwQvosOHpSGhpbXhoaKOoBNI9DXwAPhu2hqSmo0pNHR4nKi0dFinhOiQFu4ymUdhHcXTU0R4ECH2EMHgEwQ6ACQCQIdADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIRMeBbnvA9rO2f1xGQwCArSljD/0uSSdK+JyeMz5e3J77wmt8PHVHFZucXL7Ck5OpOwLQho4C3fbVkm6W9FA57fSO8XHp+PHltePHMw71yUlpbm55bW6OUAf6SKd76A9I+rqkv5bQS09ZGeYb1fveyjDfqA6g52w50G1/WtLZiDi6wXLTtudtzy8uLm51OADABjrZQ79B0mdtvyTp+5JutP29lQtFRCMi6hFRHxkZ6WA4AMB6thzoEfGNiLg6ImqSbpX084i4rbTOEhsba6/e9yYm2qsD6Dlch76GY8dWh/fYWFHP0pEjq8N7YqKoA+gLg2V8SET8UtIvy/isXpJteK+F8Ab6GnvoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBkgkAHgEwQ6ACQCQIdADJBoANAJgh0AMgEgQ4AmSDQASATBPo6ZmakwUHJLqYzM6k7AoC1lXI/9BzNzEiHDl2cP3/+4vzsbJqeAGA97KGvodForw4AqRHoazh/vr06AKRGoK9hYKC9OgCkRqCvYXq6vToApMZJ0TVcOPHZaBSHWQYGijDnhCiAXkWgr2N2lgAH0D845AIAmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJCJLQe67Wts/8L2cdvHbN9VZmMAgPZ0sof+jqSvRcSYpOsl3Wl7rJy2esPevZJ98bV3b+qOAGBtWw70iDgTEc+03v9Z0glJ2UTe3r3S6dPLa6dPE+oAelcpx9Bt1yRdJ+npMj6vF6wM843qAJBax4Fu+32Sfijp7oj402X+fNr2vO35xcXFTocDAKyho0C3fYWKMG9GxBOXWyYiGhFRj4j6yMhIJ8MBANbRyVUulvSwpBMR8a3yWuoNe/a0VweA1DrZQ79B0hcl3Wj7udbrH0vqK7lXXlkd3nv2FHUA6EWDW/3BiPhfSS6xl55DeAPoJ3xTFAAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBkoucDvfl8U7UHatpx3w7VHqip+Xwz+7HHHxyX7/O7r/EHx7sy7sxTMxq8f1C+zxq8f1AzT810ZdzUY6eSavtKNe7kY5PLtuvJxya7Mq6Ubvvq9jo7Iiod4FL1ej3m5+c3vXzz+aam/2taS28vvVsbumJIjc80NPWRqSpaTD72+IPjOv768VX1sd1jOnbnscrGnXlqRofmD62q31G/Q7M3z1Y2buqxU0m1faUad/KxSc29OLeqPrF/Qke+dKSycaV021eZ62z7aETUN1yulwO99kBNJ988uao+unNUL939Uomd9c7Yvm/tZ4bEvdX9rgbvH9T5OL+qPuABvfNv71Q2buqxU0m1fW237VpKt32Vuc6bDfSePuSy8OZCW/Vcxk7hchv8evVcxk4l1fa13bZraXttXz0d6Pt27murnsvYKQx4oK16LmOnkmr72m7btbS9tq+eDvSDEwc1dMXQstrQFUM6OHEw27HHdo+1VS/L9N9Nt1XPZexUUm1fqcad2D/RVr1MqbavFOvc04E+9ZEpNT7T0OjOUVnW6M7RrpwQTTn2sTuPrQrvqk+IStLszbO6o37Hu3stAx7o2knJlGOnkmr7SjXukS8dWRVk3TghKqXbvlKsc0+fFAUAZHJSFACweQQ6AGSCQAeATBDoAJAJAh0AMtHVq1xsL0pa/b3jzdkt6fUS2+kHrPP2sN3Webutr9T5Oo9GxMhGC3U10Dthe34zl+3khHXeHrbbOm+39ZW6t84ccgGATBDoAJCJfgr0RuoGEmCdt4ftts7bbX2lLq1z3xxDBwCsr5/20AEA6+iLQLd9k+3f2v697XtS91M129fY/oXt47aP2b4rdU/dYHvA9rO2f5y6l26wvcv2Ydu/sX3C9sdT91Q1219tbdMv2H7c9ntS91Q224/YPmv7hUtq77f9M9u/a02Hqxi75wPd9oCkByX9g6QxSV+wXe3NwdN7R9LXImJM0vWS7twG6yxJd0k6kbqJLvq2pJ9ExIckfVSZr7vtvZK+IqkeER+WNCDp1rRdVeK7km5aUbtH0lxEXCtprjVfup4PdEl/L+n3EfGHiHhL0vcl3ZK4p0pFxJmIeKb1/s8q/qLvTdtVtWxfLelmSQ+l7qUbbO+U9ElJD0tSRLwVEefSdtUVg5Lea3tQ0pCk04n7KV1E/ErSH1eUb5H0aOv9o5I+V8XY/RDoeyW9fMn8KWUebpeyXZN0naSn03ZSuQckfV3SX1M30iX7JS1K+k7rMNNDtq9M3VSVIuIVSd+UtCDpjKQ3I+Knabvqmqsi4kzr/auSrqpikH4I9G3L9vsk/VDS3RHxp9T9VMX2pyWdjYijqXvpokFJH5N0KCKuk/QXVfTf8F7ROm58i4p/zPZIutL2bWm76r4oLi2s5PLCfgj0VyRdc8n81a1a1mxfoSLMmxHxROp+KnaDpM/afknFIbUbbX8vbUuVOyXpVERc+J/XYRUBn7NJSS9GxGJEvC3pCUmfSNxTt7xm+wOS1JqerWKQfgj0/5N0re39tv9GxUmUJxP3VCnbVnFs9UREfCt1P1WLiG9ExNURUVPx+/15RGS95xYRr0p62fYHW6UJSccTttQNC5Kutz3U2sYnlPmJ4Es8Ken21vvbJf2oikEGq/jQMkXEO7b/RdL/qDgr/khEVPvE5PRukPRFSc/bfq5V+9eI+O+EPaF8X5bUbO2o/EHSPyXup1IR8bTtw5KeUXEl17PK8Fujth+X9ClJu22fknSvpP+Q9APb/6zijrOfr2RsvikKAHnoh0MuAIBNINABIBMEOgBkgkAHgEwQ6ACQCQIdADJBoANAJgh0AMjE/wNvKnI5YgSGHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "this = observX[np.asarray(chan_datas[0][:,0], dtype='int')]\n",
    "plt.scatter(this[:,0],this[:,1], c = 'g')\n",
    "\n",
    "\n",
    "this = observX[np.asarray(chan_datas[1][:,0], dtype='int')]\n",
    "plt.scatter(this[:,0],this[:,1], c = 'r')\n",
    "\n",
    "this = observX[np.asarray(chan_datas[2][:,0], dtype='int')]\n",
    "plt.scatter(this[:,0],this[:,1], c = 'b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.gaussian_process.kernels import RBF\n",
    "from sklearn.gaussian_process.kernels import ConstantKernel\n",
    "from sklearn.gaussian_process.kernels import WhiteKernel\n",
    "from sklearn.gaussian_process.kernels import Matern\n",
    "from sklearn.gaussian_process.kernels import RationalQuadratic\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "\n",
    "def get_gp(X, y, alpha=10**-3, kernel_name='const_rbf', seed = 1234):\n",
    "    kernel = ConstantKernel(\n",
    "#             constant_value_bounds=[1.,100.0]\n",
    "            ) * Matern(\n",
    "                length_scale_bounds=[2.,10.], length_scale = [10.,10.,1.,1.],\n",
    "            ) + WhiteKernel(noise_level_bounds=[1.,3.]\n",
    "#                 constant_value_bounds=[0,0.001]\n",
    "#             ) * Matern(\n",
    "#                 length_scale_bounds=[.001,100]\n",
    "            )\n",
    "    gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=alpha, random_state=seed)\n",
    "    gp.fit(X, y.ravel())\n",
    "    return gp\n",
    "# ConstantKernel?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(60, 5)"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for_gp_list = []\n",
    "eval_list   = []\n",
    "for nchan in meas_chan:\n",
    "    thisX = observX[np.asarray(chan_datas[nchan][:,0], dtype='int')]\n",
    "    thisA = chan_datas[0][:,1:]\n",
    "    v = truth(thisX,thisA)\n",
    "    full_chan_data = np.concatenate([chan_datas[nchan],v.reshape(-1,1)],axis=1)\n",
    "    all_in_one = np.concatenate([observX[np.asarray(full_chan_data[:,0],dtype='int')],full_chan_data[:,1:]],axis=1)\n",
    "    for_gp_list.append(all_in_one)\n",
    "\n",
    "#     eval_this = np.concatenate([np.asarray(full_chan_data[:,0],dtype='int'),full_chan_data[:,1:]],axis=1)\n",
    "    eval_list.append(full_chan_data)\n",
    "    \n",
    "for_gp = np.concatenate(for_gp_list)\n",
    "# len(for_gp),len(observX)\n",
    "# for_gp[0].shape\n",
    "for_gp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# full_chan_data = np.concatenate([chan_datas[1],v.reshape(-1,1)],axis=1)\n",
    "# thisX = observX[np.asarray(chan_datas[0][:,0], dtype='int')]\n",
    "# thisA = chan_datas[0][:,1:]\n",
    "# v = truth(thisX,thisA)\n",
    "\n",
    "# np.random.choice(obvervX)\n",
    "\n",
    "\n",
    "# for_gp = np.concatenate([observX[np.asarray(full_chan_data[:,0],dtype='int')],full_chan_data[:,1:]],axis=1)\n",
    "for_gp_X = for_gp[:,:4]\n",
    "for_gp_y = for_gp[:,4]\n",
    "gp = get_gp(for_gp_X,for_gp_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(a0,a1):\n",
    "    alphas = [a0,a1]\n",
    "\n",
    "    f,axarr = plt.subplots(1,3)\n",
    "\n",
    "    f.set_size_inches(15,5)\n",
    "\n",
    "    \n",
    "    vmin,vmax = 40,60 \n",
    "    \n",
    "    ax = axarr[0]\n",
    "    X = np.mgrid[0:10:31j,0:10:31j]\n",
    "    evalX = np.einsum('abc->bca',X).reshape(31**2,2)\n",
    "    evalA = np.einsum('x,a->xa',np.ones((len(evalX),)),alphas)\n",
    "    ax.scatter(evalX[:,0],evalX[:,1], c = truth(evalX,evalA), vmin = vmin,vmax = vmax)\n",
    "    ax.set_ylim(0,10)\n",
    "    ax.set_xlim(0,10)\n",
    "    \n",
    "#     evalX  = observX[global_chan_numbers]\n",
    "#     evalA = np.einsum('x,a->xa',np.ones((len(evalX),)),alphas)\n",
    "    w,werr = gp.predict(np.concatenate([evalX,evalA],axis=1), return_std = True)\n",
    "    \n",
    "    \n",
    "    ax = axarr[1]\n",
    "    ax.set_xlim(0,10)\n",
    "    ax.set_ylim(0,10)\n",
    "    ax.scatter(evalX[:,0],evalX[:,1], c = w, vmin = vmin,vmax = vmax)\n",
    "\n",
    "\n",
    "    evalX  = observX[global_chan_numbers]\n",
    "    evalA = np.einsum('x,a->xa',np.ones((len(evalX),)),alphas)\n",
    "\n",
    "    ax.scatter(evalX[:,0],evalX[:,1], c = 'r')\n",
    "\n",
    "\n",
    "    w,werr = gp.predict(np.concatenate([evalX,evalA],axis=1), return_std = True)\n",
    "    t = truth(evalX,evalA)\n",
    "    \n",
    "    ax = axarr[2]\n",
    "    for sl in sli:\n",
    "        bin_num = global_chan_numbers[sl]\n",
    "        ch_y = w[sl]\n",
    "        ch_t = t[sl]\n",
    "        y_m_err = (w-werr)[sl]\n",
    "        y_p_err = (w+werr)[sl]\n",
    "        ax.fill_between(bin_num,y_m_err,y_p_err, facecolor = 'grey', alpha = 0.2)\n",
    "        ax.scatter(bin_num,ch_y)\n",
    "        ax.plot(bin_num,ch_y)\n",
    "        ax.plot(bin_num,ch_t, c = 'k', linestyle = 'dotted', alpha = 0.5)\n",
    "    ax.set_ylim(0,1.2*vmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e8bc9251cb3e426cbc9e8551ea93c3a8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='a0', max=1.0, min=-1.0), FloatSlider(value=0.0, desc…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot(a0, a1)>"
      ]
     },
     "execution_count": 260,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ipywidgets import interact\n",
    "interact(plot, a0 = (-1,1,0.1), a1 = (-1,1,0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = {'channels': []}\n",
    "for i,(c,cX) in enumerate(zip(eval_list,[chan1,chan2,chan3])):\n",
    "    spec['channels'].append(\n",
    "        {\n",
    "            'name': 'channel{}'.format(i),\n",
    "            'samples': [\n",
    "                {\n",
    "                    'name': 'someprocess'.format(i),\n",
    "                    'data': [0.]*len(cX),\n",
    "                    'modifiers': [\n",
    "                        {\n",
    "                            \"name\": \"mygp\",\n",
    "                            \"type\": \"gpsys\",\n",
    "                            \"data\": {\n",
    "                                \"kernel_spec\": {\n",
    "                                    'type': \"c_matern_white\"\n",
    "                                },\n",
    "                                \"evals\": c.tolist(),\n",
    "                                \"plcmt\": cX.tolist()\n",
    "                            }\n",
    "                        }\n",
    "                    ]\n",
    "                },\n",
    "            ]\n",
    "        }\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "json.dump(spec,open('gpsys.json','w'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 4)"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eval_list[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}