Untitled.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import trig_poly\n",
    "import test_data\n",
    "\n",
    "import numpy as np\n",
    "from scipy import integrate as sp_integrate\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fc = 20\n",
    "min_separation = 0.55 / fc\n",
    "\n",
    "kernel = trig_poly.TrigPoly.multi_dirichlet(fc, [0.5, 0.5])\n",
    "N = int(1.0 / min_separation) - 8\n",
    "\n",
    "support = test_data.jittered_supports(N, min_separation=min_separation, jitter_factor=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shifted_supports = np.outer(np.ones(N), support) - np.outer(support, np.ones(N))\n",
    "K0 = kernel.eval(shifted_supports).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-62-4da310d22c2d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m         \u001b[0mM\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msp_integrate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mquad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msupport\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msupport\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/tkunisky/.virtualenvs/random_opt/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc\u001b[0m in \u001b[0;36mquad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)\u001b[0m\n\u001b[1;32m    339\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mweight\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    340\u001b[0m         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,\n\u001b[0;32m--> 341\u001b[0;31m                        points)\n\u001b[0m\u001b[1;32m    342\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    343\u001b[0m         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,\n",
      "\u001b[0;32m/Users/tkunisky/.virtualenvs/random_opt/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc\u001b[0m in \u001b[0;36m_quad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points)\u001b[0m\n\u001b[1;32m    446\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mpoints\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    447\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0minfbounds\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 448\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0m_quadpack\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_qagse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsabs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsrel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlimit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    449\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    450\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0m_quadpack\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_qagie\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mbound\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0minfbounds\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsabs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsrel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlimit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-62-4da310d22c2d>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(t)\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m         \u001b[0mM\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msp_integrate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mquad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msupport\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msupport\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/tkunisky/projects/sr_multisignals/trig_poly.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, t)\u001b[0m\n\u001b[1;32m    103\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    104\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 105\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    106\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    107\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__add__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/tkunisky/projects/sr_multisignals/trig_poly.pyc\u001b[0m in \u001b[0;36meval\u001b[0;34m(self, t)\u001b[0m\n\u001b[1;32m     46\u001b[0m             \u001b[0mreshaped_coeffs\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m             \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2.0\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m1j\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mrepeated_ts\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mreshaped_freqs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 48\u001b[0;31m             axis=0)\n\u001b[0m\u001b[1;32m     49\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mconjugate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/tkunisky/.virtualenvs/random_opt/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc\u001b[0m in \u001b[0;36msum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m   1776\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1777\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1778\u001b[0;31m \u001b[0;32mdef\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_NoValue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1779\u001b[0m     \"\"\"\n\u001b[1;32m   1780\u001b[0m     \u001b[0mSum\u001b[0m \u001b[0mof\u001b[0m \u001b[0marray\u001b[0m \u001b[0melements\u001b[0m \u001b[0mover\u001b[0m \u001b[0ma\u001b[0m \u001b[0mgiven\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "M = np.zeros((N, N))\n",
    "\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        M[i, j] = sp_integrate.quad(lambda t: (kernel(t - support[i]) * kernel(t - support[j])).real, 0.0, 1.0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10b1581d0>"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD/BJREFUeJzt3W2MXNV9x/Hv3+u1IYDBLg9xwY3TCEVFqDLpytSiihy5\nxM5DZRIJFF4gJyI4lYLUSEEq4k14UwlVeSgvqlQGuzFRgCQKBLd1cIhL6qaiDotDMWGTBiEnbGzZ\nJH5O8NPOvy/2ulrM3nNm5947d7z/30eKdnfOnJmTYX++M/s/D+buiEg8c9oegIi0Q+EXCUrhFwlK\n4RcJSuEXCUrhFwlK4RcJSuEXCUrhFwlqbj+f7PJFQ750yXBp+9j4Fcn+c4+8WdrmE52exyUyW5zg\nd5zyk9bNfSuF38zWAA8CQ8DD7v5A6v5Llwzz421LSttX3PPXyedb+K+vlLZNHD2a7It19XqUqzIN\nOvfcmmItM5H4fdrZ+UHXD9Pz234zGwL+EfgQcB1wu5ld1+vjiUh/VfnMvxx41d1fc/dTwOPA2nqG\nJSJNqxL+q4HXp/w8Xtz2Fma23sxGzWz0jd9OVHg6EalTlfBP98HjbR9e3X2Du4+4+8gVfzBU4elE\npE5Vwj8OTP3r3TXA3mrDEZF+qRL+54FrzezdZjYP+ASwpZ5hiUjTei71ufsZM7sb2MZkqW+Tu/80\n1Wds/IpkOe+5L/5T8jk/8l9/Vdo2NJT+SDFx+HCyPStVrsuV6nLtKgXKTNT0+1Cpzu/uW4GttYxE\nRPpK03tFglL4RYJS+EWCUvhFglL4RYJS+EWC6ut6/rlH3kwuy03V8QH+7bl/KW1bdcedyb7zd72W\nbJ84dCjZbnPLXyrvZOquHa1pkBmosvx8BlMAdOUXCUrhFwlK4RcJSuEXCUrhFwlK4RcJqq+lPp/o\nJHfZzS3LTZXztn99Y7Lvh/90VbJ96LJLk+0Th48k25PmZHYwypUCteQ3lj7999SVXyQohV8kKIVf\nJCiFXyQohV8kKIVfJCiFXySovtb5gWTNOre9dmpZbq6Ov/Wl7cn23AnBl357V2mbnz6V7Gtz0nV6\n94p1/CrbiktYuvKLBKXwiwSl8IsEpfCLBKXwiwSl8IsEpfCLBFWpzm9me4BjwARwxt1H6hhUmdT2\n2rn1+Lk6fu548BWU91+4dSzZN7cXQGpbcACfyKz3T9Xyq2wDnXtsOa/VMcnnA+7+mxoeR0T6SG/7\nRYKqGn4Hvm9mL5jZ+joGJCL9UfVt/03uvtfMrgSeMbOfufuOqXco/lFYD3AB76j4dCJSl0pXfnff\nW3w9ADwJLJ/mPhvcfcTdR4aZX+XpRKRGPYffzC4ys0vOfg98EHi5roGJSLOqvO2/CnjSJktJc4FH\n3f3pWkYlIo0z72Mdd4Et8hstse4+U5O2xL7+fuZMuu/wvGT7kVvfl2xPzQNYfc2fJfsOLbg42Z49\nEyC37793Em0V//vqzIDzyk7fzlE/2NXkDpX6RIJS+EWCUvhFglL4RYJS+EWCUvhFgur/1t0pmbKR\nd3ovK+W2184ty139eHk5b9v4C8m+VbYFB9KlPAAbLu+a+f8tcenKLxKUwi8SlMIvEpTCLxKUwi8S\nlMIvEpTCLxLUQB3RnV0e2klsYZ1Z9po7Jju3rDa1NXiT24IDXPbEi8n2zqnTpW25pczZbcFTrzlo\nye95TFd+kaAUfpGgFH6RoBR+kaAUfpGgFH6RoBR+kaD6X+dv6jjpTD3aPbMteOaY7NQ8gNx6/Fwd\nPzcPYOWhu5LtF+54pbSt8+aJZN/sXgG5bcOrzAPQHIBW6covEpTCLxKUwi8SlMIvEpTCLxKUwi8S\nlMIvElS2zm9mm4CPAgfc/fritkXAN4GlwB7gNnc/VHk0Veq+FdeVZ9e1p+rdmVp5bj1+ro7/w40P\npfvfWd7/wmd3J/t2Tp5Mtjc6D0B7AbSqmyv/14A159x2L7Dd3a8Fthc/i8h5JBt+d98BHDzn5rXA\n5uL7zcAtNY9LRBrW62f+q9x9H0Dx9cr6hiQi/dD43H4zWw+sB7iAdzT9dCLSpV6v/PvNbDFA8fVA\n2R3dfYO7j7j7yDDze3w6Ealbr+HfAqwrvl8HPFXPcESkX7LhN7PHgOeA95rZuJndCTwA3GxmvwBu\nLn4WkfOIeR9rqQtskd9oq/r2fG9Rtaac6G9zh9MPnZlDMOfCC5Ltb77/umR7ah5A7kyBhVvHku0T\nR44m2yup8Jp31T+gnb6do36wq40xNMNPJCiFXyQohV8kKIVfJCiFXyQohV8kqDilvpwGy0pVj8me\nMy9dSjz88WWlbbltwVf/YXlfgDmXXJJs7xw/nmyvpGopsMpjn6dU6hORLIVfJCiFXyQohV8kKIVf\nJCiFXyQohV8kqP4f0R1QdlvwzPbYue21U8tyVz+aruNv25veVvwDn/p0sv2C/3g52e5nziTbK/W1\nzLUrt+14cLryiwSl8IsEpfCLBKXwiwSl8IsEpfCLBKXwiwSlOv9ZTa4dTx1TDfljrjP16tT22rn1\n+Lk6/rP//HCyPbc1+KXf3pVsT7G56V/PKnMImtzKvav+A0BXfpGgFH6RoBR+kaAUfpGgFH6RoBR+\nkaAUfpGgsnV+M9sEfBQ44O7XF7fdD9wFvFHc7T5339rUIAdClbptribc4DyA3L76ufX4uTp+7lyA\nFZT3z84ByKzXz84DSO2jEKCOn9PNlf9rwJppbv+Kuy8r/je7gy8yC2XD7+47gIN9GIuI9FGVz/x3\nm9lLZrbJzBbWNiIR6Ytew/9V4D3AMmAf8KWyO5rZejMbNbPR06T3ohOR/ukp/O6+390n3L0DPAQs\nT9x3g7uPuPvIMPN7HaeI1Kyn8JvZ4ik/fgxI/8lYRAZON6W+x4CVwOVmNg58AVhpZssAB/YAn2lw\njCLSgGz43f32aW7e2MBYZq+qNeXcPIAKew3k1sTnavGpOj6k5wHk+laeB5CYHpE/S2H2zwPQDD+R\noBR+kaAUfpGgFH6RoBR+kaAUfpGgtHX3IGiyrFRly/EuVCkFVlkO3M1zp8qYlZYDw6woBerKLxKU\nwi8SlMIvEpTCLxKUwi8SlMIvEpTCLxKU6vzngyo15Uzf3JLeXD08t6w2VYuvshy4m/4Lt46Vtk0c\nPpLs2+o8gD7NAdCVXyQohV8kKIVfJCiFXyQohV8kKIVfJCiFXyQo1flngyrr+TN1+ibnATS5LTjA\n6sduKG0buuzSZN9W5wFU2YNhBlMEdOUXCUrhFwlK4RcJSuEXCUrhFwlK4RcJSuEXCSpb5zezJcAj\nwDuBDrDB3R80s0XAN4GlwB7gNnc/1NxQpRHeqdY9U89OHpOdmUOQWo8P6To+wLZf/6S0bcU9FY8H\nz7C5mfkTp09Vevw6dHPlPwN83t3/BPhz4LNmdh1wL7Dd3a8Fthc/i8h5Iht+d9/n7ruK748BY8DV\nwFpgc3G3zcAtTQ1SROo3o8/8ZrYUuAHYCVzl7vtg8h8I4Mq6Bycizek6/GZ2MfAd4HPufnQG/dab\n2aiZjZ7mZC9jFJEGdBV+MxtmMvjfcPcnipv3m9nion0xcGC6vu6+wd1H3H1kmPl1jFlEapANv5kZ\nsBEYc/cvT2naAqwrvl8HPFX/8ESkKd0s6b0JuAPYbWYvFrfdBzwAfMvM7gR+BdzazBClkqrbQFc8\najpVCswti80tq80ty02V8xo/HjxXAh2eV963T2XAbPjd/UdA2W/AqnqHIyL9ohl+IkEp/CJBKfwi\nQSn8IkEp/CJBKfwiQWnrbkmrctR0pn+2Fl5xHkCbx4Nf9tTuZHvn978vbUvNAQDwM6cTjcmub6Er\nv0hQCr9IUAq/SFAKv0hQCr9IUAq/SFAKv0hQqvNLWsX1/Mn+FfYCgC6OB09o/HjwV+9Itg/97Jel\nbRPHjiX71kVXfpGgFH6RoBR+kaAUfpGgFH6RoBR+kaAUfpGgVOeXtKr7/qf6N3gmAKSPyc71za3H\nz9Xxt33368n2lZ++q7Ttwn/P7AVw4kSyvVu68osEpfCLBKXwiwSl8IsEpfCLBKXwiwSl8IsEla3z\nm9kS4BHgnUAH2ODuD5rZ/cBdwBvFXe9z961NDVRmoQbPBID0Ofe5vfFT++pDej0+pOv4AD98+KHS\nthX3pPcSWPi9n5e22ZGhZN+pupnkcwb4vLvvMrNLgBfM7Jmi7Svu/sWun01EBkY2/O6+D9hXfH/M\nzMaAq5semIg0a0af+c1sKXADsLO46W4ze8nMNpnZwpI+681s1MxGT3Oy0mBFpD5dh9/MLga+A3zO\n3Y8CXwXeAyxj8p3Bl6br5+4b3H3E3UeGmV/DkEWkDl2F38yGmQz+N9z9CQB33+/uE+7eAR4Cljc3\nTBGpWzb8ZmbARmDM3b885fbFU+72MeDl+ocnIk0xz5RLzOwvgP8EdjNZ6gO4D7idybf8DuwBPlP8\ncbDUAlvkN9qqCqPtfRtomYVypcCU3O99lWOygTnz0x9xD398WWlbblvwNe8qf5P936ef5mjnt129\nMN38tf9HwHQPppq+yHlMM/xEglL4RYJS+EWCUvhFglL4RYJS+EWCytb567TAFvmNc/6y/A4NHvcs\ns1CDdf5Kj93F4w8tnHYpDACd479L9n36lz8ubVu++nVG/+dEV4PXlV8kKIVfJCiFXyQohV8kKIVf\nJCiFXyQohV8kqL7W+c3sDWDqnseXA7/p2wBmZlDHNqjjAo2tV3WO7V3ufkU3d+xr+N/25Gaj7j7S\n2gASBnVsgzou0Nh61dbY9LZfJCiFXySotsO/oeXnTxnUsQ3quEBj61UrY2v1M7+ItKftK7+ItKSV\n8JvZGjP7uZm9amb3tjGGMma2x8x2m9mLZjba8lg2mdkBM3t5ym2LzOwZM/tF8bV8bWj/x3a/mf26\neO1eNLMPtzS2JWb2rJmNmdlPzexvittbfe0S42rldev7234zGwL+F7gZGAeeB25391f6OpASZrYH\nGHH31mvCZvZ+4DjwiLtfX9z298BBd3+g+Idzobv/7YCM7X7geNsnNxcHyiyeerI0cAvwSVp87RLj\nuo0WXrc2rvzLgVfd/TV3PwU8DqxtYRwDz913AAfPuXktsLn4fjOTvzx9VzK2geDu+9x9V/H9MeDs\nydKtvnaJcbWijfBfDbw+5edxBuvIbwe+b2YvmNn6tgczjavOnoxUfL2y5fGcK3tycz+dc7L0wLx2\nvZx4Xbc2wj/dFkODVHK4yd3fB3wI+Gzx9la609XJzf0yzcnSA6HXE6/r1kb4x4ElU36+Btjbwjim\n5e57i68HgCcZvNOH9589JLX4eqDl8fy/QTq5ebqTpRmA126QTrxuI/zPA9ea2bvNbB7wCWBLC+N4\nGzO7qPhDDGZ2EfBBBu/04S3AuuL7dcBTLY7lLQbl5Oayk6Vp+bUbtBOvW5nkU5Qy/gEYAja5+9/1\nfRDTMLM/ZvJqD5OHmD7a5tjM7DFgJZOrvvYDXwC+C3wL+CPgV8Ct7t73P7yVjG0lMzy5uaGxlZ0s\nvZMWX7s6T7yuZTya4ScSk2b4iQSl8IsEpfCLBKXwiwSl8IsEpfCLBKXwiwSl8IsE9X+u1WMR5F57\n8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b0825d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(K0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10ae2ff10>"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADylJREFUeJzt3V2MXdV5xvHnsWObho8Wlw8hMIEaaxqrVZzKIpGIWqeU\nyOTG5CISaVVZrdOJElCDlBuUm+SmEjf5qpLQGOzaqgwIiRCI5BCIk5ZUailDhIKJQ6HIIcaWDXWD\nHVQ84zlvL+Y4HYy91/Kcfb7m/f8ka86cvWaf13vmmX3O2e+s5YgQgHyWDLsAAMNB+IGkCD+QFOEH\nkiL8QFKEH0iK8ANJEX4gKcIPJPWuQT7YJSuXxjWrljWO2fv6pY3bV7w2XX6gzmxxSHRqOhsLY2iO\nxIh5S29qOk64ZmxP4be9UdLXJC2VdG9E3NU0/ppVy/Qf31/VuM+JbZ9u3H7dt14p1hXHjhfHdE6c\nKI7RbPMvkShsr0aLdQ6uyKR7ezL+1Ozj1WMX/Ei2l0r6hqSbJa2V9Anbaxe6PwCD1cuvmeslvRQR\nL0fEtKQHJG1qpywA/dZL+K+U9Mt5nx/o3gdgDPQS/jO9gHnHi1fbk7anbE+99t8tvUYG0LNewn9A\n0vx3766SdPD0QRGxNSLWR8T6S393aQ8PB6BNvYT/aUlrbF9re7mkWyU92k5ZAPptwZf6IuKk7dsl\nfV9zl/q2R8TzrVUGoK96us4fEbsl7a4dv/f1S4vX8V/Ycnfj9gk1f71U1wuw5FhxSLkXoKZRKDrl\nMaXrv/QBLA5V38fCz0uPfQDz0d4LJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kNdDKPFa9NF6/Bl67j\nl/oAavYhtdML0Km4zh81f85Q6gWo+TtwegEWh+L3saJvpBJnfiApwg8kRfiBpAg/kBThB5Ii/EBS\nhB9IivADSQ20yUed2eKCGr02AUmDawTyyZPFfdQoNgK1MSGIRCPQYtDi95AzP5AU4QeSIvxAUoQf\nSIrwA0kRfiApwg8kRfiBpAba5BOdKK6CU5o9p2YGnkE1Aq3++kxxHzr+6/KYQhNPdCp+R9c0AgHz\ncOYHkiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5DUYGfyUUizzdPW9NoEJA2uEejPfvjXxX0s/9mB\n4pj41RuN263yml9VjUA1Sz0x208aPYXf9n5JxyXNSjoZEevbKApA/7Vx5v9wRLzewn4ADBCv+YGk\neg1/SHrc9jO2J880wPak7SnbUzPR/HoewOD0+rT/hog4aPsySU/Y/nlEPDl/QERslbRVki5aspJ3\nk4AR0dOZPyIOdj8ekfSwpOvbKApA/y04/LbPt33hqduSPiJpb1uFAeivXp72Xy7pYc+tFPMuSfdF\nxGONXxFSFK7zq9P8yqBT2C7VraRTMxFH6Tr+D3ZtL+5jYlu5n2D11/+rcXvUTAgyU/4/F1cGklTs\nBaAPYNFYcPgj4mVJ72uxFgADxKU+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0hqwJN5VCitXlPVqFKh\nonGmNBFHTQNPzcpAN/3grxq3L/v5q8V9tLEykFRxfOeaugo7oRFoHHDmB5Ii/EBShB9IivADSRF+\nICnCDyRF+IGkCD+QlGOADRkXeWV8wDc2D6ppIilx+Xeal1Q8ztKljZuX/M5vF3cx8/tXFsc8cf8/\nNm6vaSaqWaWoc/R/imNierp5e2kmplo0AvXFU7FHx+JoVYg48wNJEX4gKcIPJEX4gaQIP5AU4QeS\nIvxAUiM4mUfh+m/VZBIVk1Z0KnoB1HxNu2YlnZqJOErX8WsmBJlQzcpA5VWK4ldvNA+oWDGp5vhj\n+DjzA0kRfiApwg8kRfiBpAg/kBThB5Ii/EBShB9IavSafEpqJoEYVCPQzMny41Q0ApUm4qhp4Klp\nBPrTH3+yOOa8f//Pxu2diiafulWVCsefyT76rnjmt73d9hHbe+fdt9L2E7Zf7H68uL9lAmhbzdP+\nHZI2nnbfnZL2RMQaSXu6nwMYI8XwR8STko6edvcmSTu7t3dKuqXlugD02ULf8Ls8Ig5JUvfjZe2V\nBGAQ+v6Gn+1JSZOSdJ7e3e+HA1BpoWf+w7avkKTuxyNnGxgRWyNifUSsX6YVC3w4AG1baPgflbS5\ne3uzpEfaKQfAoNRc6rtf0r9JmrB9wPYWSXdJusn2i5Ju6n4OYIyM3oo9o6TULNTSykBevrx5+4UX\nFPfx1vuuLo754Y57i2Nu3nhrcy37Dxb30TlxojgmSg1SNbMB0Qj0DqzYA6CI8ANJEX4gKcIPJEX4\ngaQIP5AU4QeSIvxAUuM3k88gFZtIKmYDqpnVZnq6eR+lJbRUnoFHKjfwSNL3HnugcfuGT/5NcR/v\nfuYXxTGdwv+p2AQ0t5eKMRWSNgtx5geSIvxAUoQfSIrwA0kRfiApwg8kRfiBpLjO34uWVg+K2UIz\nQMUqOTUr6SypmIijdB3/n++9p7iPiW3lFYZWf62wMlDFSkcqHTdJUXFcsq4exJkfSIrwA0kRfiAp\nwg8kRfiBpAg/kBThB5Ii/EBSNPn0WxuNQBWr19RMGlKzkk5pIo6aBp4XttxdHHPzQ3/euH3J/vJk\nHp3/fas4xqppBCqdAxfn6kGc+YGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKcIPJOUYYHPCRV4ZH/CN\nA3u8VCpmDJLLv+u9rLnva8mFFxT3EVddXhzzvd33NW7/k8nJ4j7Of3p/uZY3jpXHzDY38RRnWpKq\nGrGq9JjHp2KPjsXRih+GijO/7e22j9jeO+++L9p+1faz3X8f7aVgAINX87R/h6SNZ7j/KxGxrvtv\nd7tlAei3Yvgj4klJRwdQC4AB6uUNv9tt/7T7suDisw2yPWl7yvbUjMp/WAJgMBYa/rslrZa0TtIh\nSV8628CI2BoR6yNi/TKtWODDAWjbgsIfEYcjYjYiOpLukXR9u2UB6LcFhd/2FfM+/ZikvWcbC2A0\nFSfzsH2/pA2SLrF9QNIXJG2wvU5SSNov6VN9rBFAH9Dkk0kLjUClJiBJWvJb5xXHvPmhicbt/7J1\na3Eff/jVzxTHXL1rf3FM59jxxu0xPV3cR8yUZx6qUmoWKuS11SYfAIsT4QeSIvxAUoQfSIrwA0kR\nfiApwg8kxYo9mVT1dBSuM1dMbFGzkk5pIo6aa/jP3fHN4pi1y8v7uXZH8ypFnYoJQapUHLviykvF\n1Z3qy+HMDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKZp8cE6iU+4isSqaWQqNMzWTcNQ08Pzs\nM+VGoIkVn27cft0/NDcBSVIc/3VxTOdExezVpePb1spA4swPpEX4gaQIP5AU4QeSIvxAUoQfSIrw\nA0kRfiApmnzwdsXZfspNJtGpOKfMNu+ntIqOVJ6BRyo38EjSC1vubtz+B2+Wm4ne80/7i2NqzrSd\nQpNPcaafc8CZH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeS4jo/zk0bq/6o4nr19HT5USpW0qmZiKN0\nHX/v31ZMCHJ+uZ/gum+9Uhzjiv93o05hRZ95imd+26ts/8j2PtvP2/5s9/6Vtp+w/WL348U9lAxg\nwGqe9p+U9LmIeK+kD0q6zfZaSXdK2hMRayTt6X4OYEwUwx8RhyLiJ93bxyXtk3SlpE2SdnaH7ZR0\nS7+KBNC+c3rDz/Y1kt4v6SlJl0fEIWnuF4Sky9ouDkD/VIff9gWSHpJ0R0RUr1lse9L2lO2pGVXM\nXgpgIKrCb3uZ5oK/KyK+3b37sO0rutuvkHTkTF8bEVsjYn1ErF+mFW3UDKAFNe/2W9I2Sfsi4svz\nNj0qaXP39mZJj7RfHoB+qbnOf4Okv5T0nO1nu/d9XtJdkh60vUXSK5I+3p8SAfRDMfwR8a+SztY5\ncGO75WBRaKERKGZOtlNKxUo6pYk4ahp4ShOCSNKEyvtZ/feFJp/jzZOceKb4EL9Bey+QFOEHkiL8\nQFKEH0iK8ANJEX4gKcIPJEX4gaQGP5OPCzONVDWIAJJmy8vXdE6U/5isdAasmYGnpoGnphFo43f+\nonH7klIpM0uLj/GbfVWPBLCoEH4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSGoITT6l3zeFpZ5oAloc\nit/HFpb8kqRO+eelUxhTs4RWcQYelRt4JOmx7+5q3D6xrbmZ6MQ36yfJ5cwPJEX4gaQIP5AU4QeS\nIvxAUoQfSIrwA0kN/jp/Sa99ABK9AItBzfewNDGMJEVL/QIlhZV0pIqJOFS+jl+aEOT6B18rP8ip\neqpHAlhUCD+QFOEHkiL8QFKEH0iK8ANJEX4gKcIPJOUYYEOM7dck/WLeXZdIen1gBfRunOodp1ql\n8ap3lGt9T0RcWjNwoOF/x4PbUxGxfmgFnKNxqnecapXGq95xqrUJT/uBpAg/kNSww791yI9/rsap\n3nGqVRqvesep1rMa6mt+AMMz7DM/gCEZWvhtb7T9gu2XbN85rDpq2N5v+znbz9qeGnY9p7O93fYR\n23vn3bfS9hO2X+x+vHiYNc53lnq/aPvV7jF+1vZHh1njKbZX2f6R7X22n7f92e79I3t8aw0l/LaX\nSvqGpJslrZX0Cdtrh1HLOfhwRKwb0Us8OyRtPO2+OyXtiYg1kvZ0Px8VO/TOeiXpK91jvC4idg+4\nprM5KelzEfFeSR+UdFv3Z3WUj2+VYZ35r5f0UkS8HBHTkh6QtGlItYy9iHhS0tHT7t4kaWf39k5J\ntwy0qAZnqXckRcShiPhJ9/ZxSfskXakRPr61hhX+KyX9ct7nB7r3jaqQ9LjtZ2xPDruYSpdHxCFp\n7gdY0mVDrqfG7bZ/2n1ZMHJPo21fI+n9kp7SeB7ftxlW+M80+dooX3a4ISL+SHMvU26z/cfDLmgR\nulvSaknrJB2S9KXhlvN2ti+Q9JCkOyLi2LDracOwwn9A0qp5n18l6eCQaimKiIPdj0ckPay5ly2j\n7rDtKySp+/HIkOtpFBGHI2I2IjqS7tEIHWPbyzQX/F0R8e3u3WN1fM9kWOF/WtIa29faXi7pVkmP\nDqmWRrbPt33hqduSPiJpb/NXjYRHJW3u3t4s6ZEh1lJ0KkhdH9OIHGPblrRN0r6I+PK8TWN1fM9k\naE0+3Us5X5W0VNL2iPi7oRRSYPv3NHe2l+amOr9v1Gq1fb+kDZr7a7PDkr4g6TuSHpR0taRXJH08\nIkbiTbaz1LtBc0/5Q9J+SZ869Zp6mGx/SNKPJT2n/583/vOae90/kse3Fh1+QFJ0+AFJEX4gKcIP\nJEX4gaQIP5AU4QeSIvxAUoQfSOr/AMn9gwjEj/HZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ac47510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(M)"
   ]
  },
  {
   "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": 0
}