phsamuel/kpca.ipynb

## kpca.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def gauss_kernel(x,y,sigma2):\n",
    "    return np.exp(-(x-y).transpose().dot(x-y)/sigma2)\n",
    "\n",
    "def gram(x,y,kernel):\n",
    "    g=np.zeros((x.shape[1],y.shape[1]))\n",
    "    for i in range(x.shape[1]):\n",
    "        for j in range(y.shape[1]):\n",
    "            g[i,j]=kernel(x[:,i],y[:,j])\n",
    "    return g\n",
    "\n",
    "def rand_circle(N,r,rsigma):\n",
    "    angles=np.random.uniform(0,2*np.pi,N)\n",
    "    rs=r+rsigma*np.random.randn(N)\n",
    "    return np.concatenate(([rs*np.cos(angles)],[rs*np.sin(angles)]),axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f25df6e6470>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dYYwcZ3kH8P+ze/Eho0Jg44QUfHVUpR+CTonpKdVKRFxkl7hRpTRYoFDUsxJkR8EBHfRDuSKEVQsdQkU5lFBkp3LwfQgokhWIwDKJraxyH7YCp9glCSCixNixEhKuDSAFr3M7Tz/Mzt3c3szuzM68M/PO/H+Sdb5de3dubvc/zzzzvu+KqoKIiMqplvcGEBGROQx5IqISY8gTEZUYQ56IqMQY8kREJTaW9wb4XXXVVbpt27a8N4OIyCrPPvvs71R1S9B9hQr5bdu24fTp03lvBhGRVUTkN2H3sV1DRFRiDHkiohJjyBMRlRhDnoioxBjyREQlxpAnIioxhjxRntptYH7e/UpkQKHGyRNVSrsN7NgBXL4MbNoEnDoFNJt5bxWVTCqVvIgcEZHXReQ5323vFZGnROTXva/vSeO5iEqj1XIDvtt1v7ZaeW8RlVBa7ZrvANjVd9sXAZxS1esBnOp9T1Rsptsn/sefnnYr+Hrd/To9beY5qdJSadeo6jMisq3v5jsATPf+fhRAC8C/pPF8RJG02251PD0drQ1iun0S9PinTsXbxrDHjfIYcfcHlYLJnvw1qvpq7++vAbgm6B+JyD4A+wBgYmLC4OZQpYwS2EHtkzTDMOjx5+aSPUfUn5P9/8rKZHSNuh8kG/hhsqp6WFWnVHVqy5bARdTIFkUaKTJKv9t0+8TE40f9Odn/ryyTlfxvReRaVX1VRK4F8LrB56K8ZVEpxmk3eIHqbU+UQG0202mfZPn4UX/OUfYHlYLJkH8CwB4AX+t9/YHB56K8mW51xD2IDAvUsANGs2m2jZH240c9cJg+gFFhpRLyIvJduBdZrxKRVwB8BW64PyYinwbwGwCfSOO5qKCSVIpRKvRRDiJhgVq2/nTUA4fpAxgVUlqjaz4ZcteONB6fLDBqpRg1cNNsN5g+6ygajqqpNM54pfT4K8WowRI1cNNsN1SpP122sxaKjSFP6YsTLHECN612Q5X601U7a6ENGPKUvjjBklfgVqU/XaWzFgrEkKf0xQ2WqgRuHqp01kKBGPKUPgZLsfAgWmkMeTKDwZI9jqKhAAx5ItsEhTlH0VAIhjyRTcLCnKNoKAQ//q/MBi0YVqTFxCi6sIXGuDY9hWAlX1aDTt95am+vsJFLZbjYzWsKRjDky2rQ6TtP7e01KMxtvtjNwsMYhnxZDRqrzgkydrM5zMMkKTx4BjAQQ76shlV8tp/aU/EkCdtRCw+eAQzFkC+zQRVfGatByk/SsB218GDrcSiGPBEll0bYjlJ4sPU4FEOeiKIZ1I7JK2zZehyKIU9Eww1rx+QZtmw9DsSQJ6LhorRjGLaFxBmvRDRcmjNqOds6U6zki4zjf6ko0mrHcMhj5hjyeYgS3nwzUNGk0Y7JcsgjiyQADPnsRQ1vjv+lMspqFA6LpFXsyWctbBXBflxVkMrIa/scPJg8eAf19qO+zyqAlXzWolYyHP9LZZVG22dYpc5JUqsY8lmLE94ckkYUbFg7k0XSKoZ8HhjeRMlEqdT5PgPAkCciG7FSj4whT0R2YqUeCUfXEBGVGEOeiKjEGPJERMNYvN4Oe/JERINYPnuWlTwR0SCWz55lyBNRdUVpw1i+xAjbNURUTVHbMM0msLAAHDsG7N5tVasGYMgTUVVFXem13QZmZ91/s7QETE5aFfRs1xClqH2hjfmlebQv2DcKo3KitmEs78mzkifyaV9oo3Wuhelt02hubYbeFvZ/dyzuwOXuZWyqb8KpmVMD/z3lLOrSCJavaMmQH4afLlMZQSENIHJwt861cLl7GV3t4nL3MlrnWgz5oouyNILl6+QYD3kROQfgjwC6AFZUdcr0c6bG8vGxVRe1AvcEhTSAyME9vW0am+qbVg8I09umU/xpKFcWr5OTVSV/q6r+LqPnSg8/gq/wwoJ8lNZJWEhHDe7m1iZOzZyKdWAhMo3tmkEs78XZblglPijIR2mdhIV0nOBubm1uONgM+79xzziI4sgi5BXAkyKiAA6p6mH/nSKyD8A+AJiYmMhgc2KwvBdnsyiV+KAgH7V10h/SYbel9TPwYi2ZlkXIf1hVL4rI1QCeEpFfquoz3p290D8MAFNTU5rB9qw37MKqxb0423gVbWNzA8deOIbOSgcOHFxauYTFs4sbwm9QkHtV+eLZRQDAz1//eebVcpSzCV6sJdOMh7yqXux9fV1EHgdwM4BnBv+vjPDCamF4FW2n24GjDgQChXvMVygeOfMIZm6cWReAUXrgR88eXT1Y1KSG8fp4ZtVylLMJXqwl04yGvIi8E0BNVf/Y+/tHAfybyeeMhRdWC8OraB11ALjB7g/6FWclsMod1EpZfUy4j+mog85KBwdaB3Bg+oDxoB90EPL34XmxlkwyXclfA+BxEfGe61FVPWH4OaPjhdXC8Cpar5KvoYax+hgEghVnZaQqd/Uxe5W8QODAwcmXT2Lp/FImFX3QQSioDz93y5zR7aDqMhryqvoSgBtNPkcivLBaKHtu3AMA2H7tdiy/tbwa6qNWuf5K2uvzn3z5JBx1cu1/x+nDDxoiyuqfouAQSl5YzV1/ZRvUex+Vv5KevHoSS+eXMu1/B4Vx1D582MgbjsihOBjylLusRphkPVkpLIyjbkfYfuGIHAvluDxKtUKe69BkKmpLIcsRJqOOeR/FoDCOsh1h+4UjciyT8yi+6oT8KDuaB4WRxWkplHU5gKRh3L9fAGB+aT50RA779AWV8yi+6oR83B3NMfSJxG0pZFlhZyWNg5e3X4aNyInbp+cBIUM5j+KrTsjH3dEcQ58IWwqutA5eww6aYfeHrY/PC7cZynkUX3VCftiO7m/NcAx9bF6gvNl5E2dePYPP/s1nceX4lawWUzDsoBl0f1iY88KtAQVeHqU6IQ+E7+iw1gzH0EfmBcqllUurs1SffOlJHPr7QwyQFAxr/QTdP780HxjmPMtKWcFbu9UK+TBhrRmOoY/Mqw69gPcce+EY9v31vpy2qlyGtX767w+r7lvnWljYtbA64YwH4YQK3tqtZsizNZM6L1D8lTwA7L5hd45bVW1Bo3OCPt7QG7HDsB9RwfOjeiHP1owR/kDxevK7b9jNKj5H/Rdd+9s3i2cXcfTsUV6ATarg+VG9kGdrJpFBQ+/KOAzSVkEXXfvbN0D0z6+lIQqcH9UL+YKfWhUZh97ZI2gEzdwtcxvaN/5KvrG5MbB1w7H1dqpeyBf81KrIOPTOHmEjaPrPtvwttv3H98NxHIyPbfxgFR7g7VW9kAcKfWpVVO0LbZz//XmM1cYABxx6V3BRZ9t6t3/kOx/BirMCAOisdCJPtqLiq2bIUyz+Kq5eq2Pvh/ZuWA6YiifqNZLWuRa62l39vlarRZpsRXZgyNNQ/ioODjDx7gkGfIlMb5vGeH0cnW4HdanjodsfijTZiiIowCKHDHkailVcuTW3NrGwawHHXjg2cNgrR0/FVJCZsAx5WidoBIU/BG669ia0zrVWbyf7tS+0MXtiFpe7l7F0fgmTV0/yd5uGgsyEZcjTqkEfNzd7YhadbgdPvvQkaqgFjsAgO/GiqiH+4dr1OnD+vFvdZxz0tUyfjQot6M3uv91RBwDgwFl3P9nNa8fVpb6hHde+0Mb80jzaF9r5baCtvOHae/cCIsDDDwO33grcd58b9hlhJU+rhn3cXKfbgaMOaqixN18iYRdVOTY+Bc2m26ZZWXHbNt0ucOgQcPTo+h69wQu0DHlaFfZm99/e2Nzg6oUlFHRRlW2clHhtm0uXAFX3j79Hb/gCLUOeVnFdGvLjqKqUeG2bxUXgyBG3mvcvqWL4Ai1DngDw1Jw24tj4FHmz7Gdm1rdl2m33gmy97v47A+tpMeQJAE/NKRjP4FLmX1LF36YZG3Mv0M7MpN6T5+gaAjB4hAURGeBv06ysABMTRoZXspInADw1J8pcRsueM+RpFU/NiTKU0bLnDHkiorxksOw5e/JEI2i3gfn5TCcuEo2ElTxRRN6kxEYDmJ0NnrtSgJVlidZhyFcEP58zGf9oNxHAcdw/GU5cJBoJQ74CONFpvVGqbf9ot1rNnbsisn5QxOLi2sz1HFeWJVqHIV8BnOi0JqjaBtbaMMvL6ycjegcDb7Rbp+OG/Oc/D1x55fp/+8gjbsAD7kHA0Ig4olgY8hXANUjW9C8TsrjoLgjY6bjtl1rNnXx4++3A8eNry4ycOgUsLACf+Yx72ze/CTz99FqlvrjoPh7gVvj33FO+Kp4tPzsx5CuAE53W9M8/ee21tRYLsNZn//731/6P13o5f94NeMA9KCwurlXxR46sPcamTe7s9DJhy89eDPmK4EQnl3/+SaMB7N+/Fs5hul3g8GF31nmQVmst/EWAu+8uXxUf1vJjdV98xkNeRHYB+CaAOoD/VNWvmX5Ooih+9rO1cB7m3Dn3j4h7UKjX3bOA++4Dtm9f36/fvt3kVucjqOXH6t4ORkNeROoAvgXgbwG8AuCnIvKEqr5g8nlpuP4KrCoVWbvtfgKbt/BfXF7V3+2utXTGxoAvfAF44AH39s99zj2AGFhQMDdBLb/5pXle0A9ToAkTpiv5mwG8qKovAYCIfA/AHQAY8jnqr8AWdi1g9sRsJSqyxUW34gaAt99O5zFXVoBHH3UD3nHcxz90yO3T33NPecK+v+XHC/ohCjZhwvSyBu8HcMH3/Su921aJyD4ROS0ip9944w3Dm0PAxv7qsReOBX6AN0X3yituwHu8sfKHDrnvd1uWP4jzwd1edX/w1oOlLgxi8yZM+D/pKUe5r12jqodVdUpVp7Zs2ZL35lRC/9rxu2/YXZm15Gdm3OJKxMzjb9u2vg3knxhVdN4Z3pef/jJ2LO6IHPRzt8wx4D0FnDBhul1zEcBW3/cf6N1GOQrqr05ePVmJnnyzCTz4IHD//W6bZdjImrj8n+QGuBdiDS4Vnqo0J81V5RrPBq2W+8ICCjNhwnTI/xTA9SJyHdxwvwvAPxp+zkqL+uby+qve6fn0tmnM3TKX4ZbmZ3nZba2kHfDA+set1YCdO4EDB3J/n0eSVo+90qNu+idiFGDChNGQV9UVEbkfwI/hDqE8oqrPm3zODQp0ldu0uG+uNN+MNlVu/vchEDyM8l3vAv7wh7Xvr7wSePPNwY8r4lbxqu6f8XF7Ah6IN2lu0O+70stoZPRBIHEYHyevqscBHDf9PIEKdpXbtLhvrrTejDZVbt4xf2HBrei9CVHeGTbgVuB/+tPa9zffDJw5s/5xxsbcIB8bcyc/bd/uDpt85BF31E6t5j6HbS+3KJPmhv2+Kz/qJoMPAomj3DNe+xcqKdmygP3VVNw3V1pvRlsqt0HHfH+P3ltGGHAr8yuvXD+j9d573bPw/mJtft59DMdx/93yctY/YTqGnZUN+31zGY1iKXfIZ/RBuXkIq6bivLnSejPaUrmFHfP37QMmJ93WysmTawHvLSW8ezewtLS+zRpUrJXh5RblrCzK75vLaBRHuUO+2XTPmY8dc9+pJariw6qpuG+uNN6MtlRug0K42XRD3gvzen39RKbJyeFt1gK2Y2OLclY2yu/bpms2ZVOukO+/yNpur31O29KS+0618Z0XoGjVsw2V27AQHnR/1DZrwdqxsXmvq063g5rU0NjcCPx3cX7fUc4OeBAwpzwhH9RwLXFP3pbquWiGhbDtIZ1Uc2sTC7sWsP/4fnSdLmZPzGLy6slEr69hZwc2Xbi3UXlCPijQy9AkHcCG6pnss/zWMlQVDpxULqQPO+u05cK9rcoT8kGBXoYmKVHGorQC47RXhp11Fq31WDaiJqb9jWhqakpPnz49+gNUaOITkUmDQnyUSXfDDgjsyScjIs+q6lTQfeWp5AE2VFPGN151DWoFxmmvRD0gsPVoTrlCnlLDi2EUJk57hf32/DHkK6h9oY3Fs4sAgJkbZ0aa1UjVFWdkF/vtERhuMzPkK6Z9oY3po9O43HVX53rkzCN4es/Tq29Ur0XT2Nzgm5NCRW2vcKjvEBmsr8WQr5jWuRbe7q597p2/Sg/6WMDlt5b55qRE2G8fIIO5PAz5ipneNo0r6lesVvL+Kr2/RbP81vK6NeZ5IbYa+HvOUAZzeRjyFdPc2kRrTyuwJz+of8oLsdXA33PGMpjLw5CvoLDT50H9U16IrQb+nnNgeOg3Q57WCToAtC+0cf7351Gv1QEHvBBbYhwNUz4MeRrIf/o+VhvD3g/tDR12yV6u/TgapnwY8jSQ//QdDjDx7onQgGcv135JDtQ8yBdTNUOea9xEFvX0nb1c+wUdqAFE/mBvHuSLqXohX7EP904q6uk7e7nZMFkt9x+oF88u4ujZo5GCmwf54qpeyJf4g0RMiTKZhb1c88Kq5WHBH/XA0H+gBhA5uHmQL67qhXzJP0gkT5zZaFZQtQxgYJskThul/0ANYF0lPyi4eZD38beDgdxbw9ULeX6QCFmqv1pubG7gQOsAOt0OHA3+FKe4bZT+A3Wc4OZBHuvbwfU6IAKsrOTaGq5eyAMbJx/wQixZwF8tNzY3MHtiFp2VDhw4qEktsNpO2kZhcMfkbwc7jnubaq6t4WqGvB8vxJJFvNCdX5rH5e5lN+BRw87rduLA9IENgcw2Ssb87eD+Sj6n1jBDfnERuHQp96MtURz9FXpQwHtYjWfA3w3wt4OB3LsE1Q75dhs4csQNeAAYG+OFWLKCqQqdE5pGENQNmFtbvTXvorHaId9qub0zwD2tuvvu3H8hRFGlXaFzQtOICj4su5b3BuTK65/V68A73gHMzOS9RdZpX2hjfmke7QvtvDelVPLYr2FDNGkIf44UcFh2tSt5DqdMJK3Kjy2C9cI+oauxuTHSJ3WNOhmq8hOaoo66K3iOVDvkAeNrOZdZGlPZ2SLYyL9fO90O9h/fD8dxVodKjtfHI++nQfu3P/w5EsenRKPuyhHyHOeeizQqP655spF/v9akhq7ThQN3zHXYpKcwYfs3LPw5EqcnTp+94AcE+0O+4Du4zOJWfkFtA7YINooz6WlYKyZs//LgOkSc5U8KfuHV/pAv+A4uu6iV36DKMemBooz8+3Xy6snVwPf35KO0usL2Lw+uQ8Tpsxd8PSz7Q77gO5hcgyrHpAeKsgvbP1Gr8aD/z/57gP62b9TrdbzwaljBdzC52L9Pn3+fjtXGcP7359G+0I68T9h/90na9i3wAI5yjJNvNt0ZZiZ3crsNzM+7Xyk2r3I8eOvBkStwL9TqUmeLAWv7dO+H9kKhePi/H8aOxR2cszCKoLZvSdhfyWeBF3dTkbRyrHqLIeh6RHNrE61zLXSdLs9wkhjW9rV4BJ+xkBeRAwD2Anijd9O/qupxU89nFC/uFsYoBwovHEedTFQEg65H8CJqCga1fS0v8kxX8g+o6r8bfg7zeHHXGv3VrheO3gdr1FDD+Fj0yURFMezCdZXPcCIbVo2H9dUtL/LYromCF3etEFTteuHoaG8yERx0Vjo40DowcHneohlWrfvPcKoyzHRVlFZKkmrc8iLPdMjfLyIzAE4D+GdV/b/+fyAi+wDsA4CJiQnDm5NAga+ekyuo2vXC0avkATfon3rpKSydX7Kmoo9arVdumGnU8E5SjVte5CUaXSMiJ0XkuYA/dwD4NoC/BHATgFcBfCPoMVT1sKpOqerUli1bkmwOVVzQ6Jvm1iYWdi1g53U78anJT6Em7kteoeisdKxaabG5tYm5W+YGhnblVpKMOiom6UqRWYzgMyRRJa+qO6P8OxF5GMAPkzwXVVfU9kNQtdu+0MbsidnVdWDU+4AYALVarXQXKSt3ETZqK8XyajwJk6NrrlXVV3vf3gngOVPPlSmLh1LZKG77oX/0jb+yVSjqtTocdVCXOh66/aGRWxlF7XtX7iJsnPCuaMvVZE/+6yJyEwAFcA7AvQafKxuWD6WyUdJZrv2Vrbc2e5IALHrfu3IzWSsa3lEZC3lV/SdTj50by4dS2Shp+8FEZTvswJOkyi/qGQLZi0Mo47B8KJWN0gjp/so2aZAOOvAkqfKLfoZAdmLIx7Vnj/t1Zsat4tmjNy7N9kMaQTrowJOkvcQF2MgEhnxU/f34mRn26C2UVpCGHXiStJcqNzKGMsGQjypsPC579FYxHaRJ2kuVGxljEs+wV4l/3HDepqam9PTp03lvRrCgqh1wX0Rvvw1ccUV4yPMFVyi8uFlyFTzDFpFnVXUq6D5W8lEFjcdttwER937va78KvuCKrnJDDKum1QI6HcBx3K8VP8NmyMfRPx631QJWVgBV92vQi4nDLomy1Wi4AQ+4XxuNfLcnZ+X4ZKi8RFkPI+maGf34CVVEgy0vA7VetNVq7vcVxko+iShTqtNcM4OtH6LhpqeB8XHOZ+lhyEcx6MJplCnVaU279rd+Oh3gwAH3D4OeaE2FFyMLwpAfpkjVs9f68S4qnTwJLC2xoifqx/VsVrEnP0yRPsXdq1B27nR7jY6T/zYRUaEx5IdJ48JpmhdLm023RTM+nt7FXKIscNBALtiuGSZpf89Eu4c9R7KFdz2r0QBmZ4vR9qwYhnwUSfp7psbJs+dIRecvcETc9qK/xcjXbybYrjEt7XHyRLbwFziO474H4rwP2N5JBSt509haoarq//yFhQV3YlKjsTZYIOz9UKRRbZZjyGeBrRWqorD1nqKEN5cDSQ1DnoJx5UxKQ9B6T1HCm5/ClhqGPG3EU2UyJWp4s82ZGoY8bcRTZTIlrIUTFOZsc6aCIU8b8VSZTPKHN88ajeMQStrIq7YOHuSbjsxKsmwIh1hGwkreRllcFOWpMmVh1LNGngFExpC3gT/UAb64qTxGvcDK60aRMeSLrr9i2bNn/Yt7cTHdqp5DJylro5w18rpRZAz5ouuvWIC1F/fYGHDkiHtfGlU9T4HJFhxiGRlDvuj6K5aZGfdPqwWcPw88/HB6p6w8BSab8LpRJAz5ogurWLzxxUePpnfKylNgotIRVc17G1ZNTU3p6dOn894Mu8TpoUf5t+zJlxt/v6UkIs+q6lTgfQz5imC/nfgaKK1BIc/JUFVRpM+qpXzwNVBJDPmqsP3DSzi7MTnbXwM0El54rYqoQ86K2LNlmyEdHHZYSQz5Khk25CyvMB12YOHQzvRw2GHlMORpTR5hGuXAYuvQziKeFVHlMORpTR5hGuXAEmcN8qIY9ayo6D8XWYchT2vy6NnG+aQgm9YgH+WsyIafi6yTaHSNiHxcRJ4XEUdEpvrumxORF0XkVyJyW7LNpMw0m8DcXHqLnQ0bETPK2vU2DAUcZSTLoJ+Lo4toREkr+ecAfAzAIf+NInIDgLsAfBDAnwM4KSJ/pardhM9HtohTlca9GGhDj36Us6L+n6vRcIO90QBmZ1nh00gShbyq/gIARKT/rjsAfE9VOwBeFpEXAdwMgGVIVZi8iGvLUMC4By//z+UP9lrN3Y+Ow9FFFJupnvz7AfyX7/tXerdtICL7AOwDgImJCUObQ5kzXW2bGgqY94VP7+ean187SKq6QS9S3DMXKqyhIS8iJwG8L+CuL6nqD5JugKoeBnAYcNeuSfp4ZFCcALSl2vZL68JnGgeK/oPkwgKwvGzPvqTCGBryqrpzhMe9CGCr7/sP9G4jW40SgLZNvEmjxZTWgcLGgyQVkqm1a54AcJeIjIvIdQCuB/ATQ89FWbBhREtSjYbbEqnVRm+LpLmf0hjpxFE5lZeoJy8idwJ4EMAWAD8SkTOqepuqPi8ijwF4AcAKgP0cWWM5G0a0JNFuuxc6Hccd9riwMFq4Fmk/cdw9IfnomscBPB5y31cBfDXJ41OBlL194FXgjuNW88vLoz1OkfYT1/whcMYrxWFbjz2ONCvwouynIp1VUG4Y8kRAsSrwtJTxZ6LY+PF/RESW48f/ERFVFEOeiKjEGPJERCXGkCciKjGGPBFRiTHkiYhKrFBDKEXkDQC/SenhrgLwu5Qeq6y4j6LhfhqO+ygaU/vpL1R1S9AdhQr5NInI6bBxo+TiPoqG+2k47qNo8thPbNcQEZUYQ56IqMTKHPKH894AC3AfRcP9NBz3UTSZ76fS9uSJiKjclTwRUeUx5ImISqxUIS8iHxeR50XEEZGpvvvmRORFEfmViNyW1zYWjYgcEJGLInKm9+f2vLepKERkV+/18qKIfDHv7SkqETknIj/vvX64VniPiBwRkddF5Dnfbe8VkadE5Ne9r+8xvR2lCnkAzwH4GIBn/DeKyA0A7gLwQQC7APyHiNSz37zCekBVb+r9OZ73xhRB7/XxLQB/B+AGAJ/svY4o2K291w/Hyq/5Dty88fsigFOqej2AU73vjSpVyKvqL1T1VwF33QHge6raUdWXAbwI4OZst44sczOAF1X1JVW9DOB7cF9HRJGo6jMA/rfv5jsAHO39/SiAfzC9HaUK+QHeD+CC7/tXereR634R+Z/e6aXx00dL8DUTnQJ4UkSeFZF9eW9MwV2jqq/2/v4agGtMP6F1n/EqIicBvC/gri+p6g+y3h4bDNpnAL4N4CDcN+pBAN8AcE92W0cl8GFVvSgiVwN4SkR+2atiaQBVVRExPobdupBX1Z0j/LeLALb6vv9A77ZKiLrPRORhAD80vDm2qPRrJg5Vvdj7+rqIPA631cWQD/ZbEblWVV8VkWsBvG76CavSrnkCwF0iMi4i1wG4HsBPct6mQui90Dx3wr14TcBPAVwvIteJyCa4F+6fyHmbCkdE3ikif+b9HcBHwdfQIE8A2NP7+x4AxrsP1lXyg4jInQAeBLAFwI9E5Iyq3qaqz4vIYwBeALACYL+qdvPc1gL5uojcBLddcw7AvfluTjGo6oqI3A/gxwDqAI6o6vM5b1YRXQPgcREB3Dx5VFVP5LtJxSAi3wUwDeAqEXkFwFcAfA3AYyLyabjLqn/C+HZwWQMiovKqSruGiKiSGPJERCXGkCciKjGGPBFRiTHkiYhKjCFPRFRiDHkiohL7fyo4yRhgoR0AAAABSURBVI27MYJ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline \n",
    "from matplotlib import pyplot as plt\n",
    "circ=rand_circle(100,10,0.5)\n",
    "circ=np.concatenate((circ,rand_circle(100,5,0.5)),axis=1)\n",
    "circ=np.concatenate((circ,rand_circle(100,0,0.5)),axis=1)\n",
    "\n",
    "plt.plot(circ[0,:99],circ[1,:99],'r.')\n",
    "plt.plot(circ[0,100:199],circ[1,100:199],'g.')\n",
    "plt.plot(circ[0,200:],circ[1,200:],'b.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute normalized Gram matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "G=gram(circ,circ, lambda x,y: gauss_kernel(x,y,4.))\n",
    "C=np.eye(G.shape[0])-np.ones(G.shape)/G.shape[0] # centering\n",
    "Gt=C.dot(G).dot(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find the first two principal components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "w,v=np.linalg.eig(Gt)\n",
    "\n",
    "ind=w.argsort()[::-1]\n",
    "w=w[ind]\n",
    "v=v[:,ind]\n",
    "\n",
    "a0=v[:,0]\n",
    "a1=v[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project data to principal component and visualize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f25df6bea90>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAbh0lEQVR4nO3df3Bc5X3v8fd315KAQsKN7AQwFkpvPEnc+toEjRNBJ2xjSIEwsS+kCcncLBAakRSGuAMDpJmkTDODCUwzzgATULHBapjQ3vIjJnWA4LLY1JsU+VeJzeXGl6tYBqc4ouXHBevH7vf+cXallbSyZO2RdrXP5zWzs3vOebTn8Y7nu89+z/c8j7k7IiJS/xLV7oCIiMwOBXwRkUAo4IuIBEIBX0QkEAr4IiKBmFftDkxk/vz53traWu1uiIjMKTt27Piduy8od6xmA35rayvd3d3V7oaIyJxiZr+Z6JhSOiIigVDAFxEJhAK+iEggFPBFRAJRccA3s0Vm9oyZ7TOzvWb2jTJtUmb2hpntLjy+U+l5RUTk2MRRpTMEXO/uO83sJGCHmf3c3feNabfN3S+O4XwiIjINFY/w3f2Qu+8svH4LeBFYWOn7ViLbm2XttrVke7PV7IaISE2JtQ7fzFqBM4FfljncbmZ7gFeBG9x9b5m/7wA6AFpaWqbVh2xvlpVdKxnIDdCYbGRLegvti9qn9V4iIvUktou2ZnYi8DCwxt3fHHN4J3CGuy8D7gQeK/ce7t7p7m3u3rZgQdkbxSaV6ckwkBsg5zkGcgNkejLTeh8RkXoTS8A3swaiYP+guz8y9ri7v+nubxdebwYazGx+HOceK9WaojHZSMISJCxB8wnNM3EaEZE5J44qHQPWAy+6+/cnaHNKoR1mtqJw3r5Kz11O+6J21l2wjoQlyOVzrHlijXL5IiLEk8M/B/gy8IKZ7S7s+0ugBcDd7wE+B3zdzIaAd4HLfAbXVux7pw93J09+OK2jPL6IhK7igO/uzwE2SZu7gLsqPddUNZ/QjJlhbsNpnWxvlkxPhlRrSsFfRIJUs7NlTle2N8uaJ9aQy+dwnKH8ENf97DocJ5fPqXJHRIJVd1MrFKt0nChj5DgDuQEGc4PkPMeRoSN07emqci9FRGZf3QX80iodgAQJGpONJBNJIPoCuH/3/bqQKyLBqbuUTvuidrakt5DpydB8QjN97/SRak3RtaeLe3fcO5zm0YVcEQlN3QV8iIJ+uWC+cc/G4TtwU62p2e+YiEgV1WXALxpbmVMc+atSR0RCVLcBf6I5dRToRSRUdXfRtqhrTxdHho5oTh0RkYK6DPjZ3iz3775/uDQzmUgqZy8iwavLgJ/pyTCUHwLAML6y/CtK5YhI8Ooy4Bdr8ZOW5Lh5x5Felq52l0REqq4uL9oWK3KKd9S+8NoLE1bnaI4dEQlFXQb8oo17NtI/1E+ePAlL0JRsGjWPjlbHEpGQ1GVKB6I8fjHYA+Q9P65aR6tjiUhI6jbgN5/QPBzsYWROndJqndJcv+6+FZF6V7cpnV2HdmEYjmMY5/3+edySumVUykZ334pISOoy4Gd7s2zYvWG4Dr8x2Tgu2Bfp7lsRCUVdpnQyPRly+RwQ1eFfufxKBXURCV5dBnzV4YuIjFdxSsfMFgFdwAcABzrd/Qdj2hjwA+Ai4B3gCnffWem5J6LcvIjIeHHk8IeA6919p5mdBOwws5+7+76SNhcCiwuPjwM/LDzPGOXmRURGqzjgu/sh4FDh9Vtm9iKwECgN+KuALnd34BdmdrKZnVr42xmT7c0O322bXpbWF4CIBC3WKh0zawXOBH455tBCoLdk+2Bh36iAb2YdQAdAS0tLRX3J9mb5441/TH+uH4ANuzeQuVzLGopIuGK7aGtmJwIPA2vc/c3pvIe7d7p7m7u3LViwoKL+FO+iLRrMDepOWhEJWiwB38waiIL9g+7+SJkmrwCLSrZPL+ybMcVKnaKGZIPupBWRoMVRpWPAeuBFd//+BM02Adea2UNEF2vfmOn8ffuidp65/Bnl8EVECuLI4Z8DfBl4wcx2F/b9JdAC4O73AJuJSjL3E5VlXhnDeSelSh0RkRFxVOk8B9gkbRy4ptJziYjI9NXlnbYiIjKeAr6ISCDqPuBne7Os3baWbG+22l0REamqupweuUhLGIqIjKjrEb6WMBQRGVHXAV9LGIqIjKjrlI6mSRYRGVHXAR9085WISFFdp3RERGREMAFf5ZkiErq6T+mAyjNFRCCQEb7KM0VEAgn4Ks8UEQkkpaPyTBGRQAI+qDxTRCSIlI6IiCjgi4gEQwFfRCQQsQR8M9tgZq+Z2a8mOJ4yszfMbHfh8Z04zisiIlMX10XbB4C7gK6jtNnm7hfHdD4RETlGsYzw3X0r8Hoc7yUiIjNjNnP47Wa2x8x+ZmZ/UK6BmXWYWbeZdR8+fHgWuyYiUv9mK+DvBM5w92XAncBj5Rq5e6e7t7l724IFC2apayIiYZiVgO/ub7r724XXm4EGM5s/G+cWEZHIrAR8MzvFzKzwekXhvH2zcW4REYnEUqVjZj8GUsB8MzsI/BXQAODu9wCfA75uZkPAu8Bl7u5xnFtERKYmloDv7l+c5PhdRGWbVZPtzWryNBEJWhCTp2kBFBGRQKZW0AIoIiIBjPCzvVkOvHGAZCIJebQAiogEq64DfmkqZ15iHl/92FdJL0srnSMiQarrgF+ayiEPLe9tUbAXkWDVdQ5fa9mKiIyo6xG+1rIVERlR1wEftJatiEhRXad0ysn2Zlm7bS3Z3my1uyIiMqvqfoRfSjdgiUjIghrh6wYsEQlZUAFfVTsiErKgUjqq2hGRkAUV8EFVOyISrqBSOiIiIVPAFxEJhAK+iEggggz4uvlKREIU15q2G4CLgdfc/Q/LHDfgB8BFwDvAFe6+M45zT1VxicPmE5pZ88Qa+of6SSQS3H3R3XSc1TGbXRERqYq4qnQeIFqztmuC4xcCiwuPjwM/LDzPitI7bBOWYCg/hOPk83mu3XwtS9+/VJU7IlL3YknpuPtW4PWjNFkFdHnkF8DJZnZqHOeeitI7bHOeI/rBEcl5TnfcikgQZiuHvxDoLdk+WNg3K0rvsG1KNnHD2TfQkGggYQmakk2641ZEglBTN16ZWQfQAdDS0hLb+5a7w3b1h1frjlsRCYq5ezxvZNYK/HSCi7b3Ahl3/3Fh+yUg5e6HJnq/trY27+7ujqVvIiKhMLMd7t5W7thspXQ2AWmLfAJ442jBfiaoFFNEQhdXWeaPgRQw38wOAn8FNAC4+z3AZqKSzP1EZZlXxnHeqdI8+CIiMQV8d//iJMcduCaOc01HuXnwFfBFJDRB3GmrefBFRGqsSmemFKt0uvZMdF+YiEj9C2KEX7Rxz0b+duffsrJrpS7eikhwggn4Ws9WREIXTMBXHl9EQhdEDh+0nq2ISDABH7SerYiELZiUDuhuWxEJWzAjfN1tKyKhC2aEf7QqHY38RSQEwYzwi1U6xRF+sUpHI38RCUUwAX+iKh3NsyMioQgm4MPoKp3SRc3LjfxFROpNUAG/KNubJbUxxWBukIZkA3deeCd97/SpPl9E6lqQAb9rTxcDuQEABnID7Dq0ix9e/MMq90pEZGYFU6UjIhK6IAN+elmapmQThtGUbCK9LF3tLomIzLggUzoAVy6PVllML0srby8iQQgu4I+tu9foXkRCEUtKx8wuMLOXzGy/md1c5vgVZnbYzHYXHn8Wx3mnQ/Pii0ioKh7hm1kSuBs4HzgIPG9mm9x935imf+/u11Z6vkpNdMetiEi9iyOlswLY7+4vA5jZQ8AqYGzArwmaF19EQhVHSmch0FuyfbCwb6xLzezfzOwfzWxRuTcysw4z6zaz7sOHD8fQNRERKZqti7aPAz92934zuxrYCHxqbCN37wQ6Adra2nwmOqLJ0kQkVHGM8F8BSkfspxf2DXP3PnfvL2zeB5wVw3mnRRdtRSRUcQT854HFZvZBM2sELgM2lTYws1NLNj8LvBjDeadFi5mLSKgqTum4+5CZXQs8CSSBDe6+18z+Guh2903AdWb2WWAIeB24otLzTpcu2opIqMx9RlLlFWtra/Pu7u5qd0NEZE4xsx3u3lbuWJBz6YiIhCjYgK91bEUkNMHNpQMqzRSRMAU5wldppoiEKMiAr9JMEQlRkCkdlWaKSIiCCvjZ3uyoIK9ALyIhCSbgZ3uznPvAuQzmB2lINPDsFc8q4ItIUILJ4d/+L7czmB8EYDA/yM1Pj1unRUSkrgUT8F9969VR288deE41+CISlGAC/uLmxaO28+RVjilS47JZWLs2epbKBZPDP/z/Ri+okrCEyjFFalg2CytXwsAANDbCli3QrstuFQlmhH/pkktHbd9w9g26aCtSwzKZKNjnctFzJlPtHs19wQT8pe9fyrxE9INmXmIeqz+8uso9EpGjaW4GM0gkohF+KlW+XSVpn9BSRsGkdDI9GYpTQbs7mZ6MRvgiNSqbhTVrIJ+HZBLWrSufzqkk7RNiyiiYgF+cTqE4YZry9yK1q5jOyeejUX5f39HblaZ92tujYJ7JRL8KJgriE/3tWNksdHVFr888M+pLKgUvvAAPPwyXXgodHZX9e2dLMAG/dDqF5hOa6drTRdeeLtLL0hrpi9SYVCoadRdH3xOlc8q1m8rIPZuFAweiXw8w8Tmy2Wj/wMDIvkQi+hLK5aLtp56KnpcuHf0lU/zSaW6GXbuiNul0lX9FuHtNPs466yyfCdsPbPfG7zY6t+Dcgjd9t8m3H9g+I+cSkenbvt391luj52Npd+ut7smkO0TPt946vv3xx0fHmprcv/a18ufYvt390592N4ve62iPFStG3vP4493vvTd6Hvu3jY3jzzXVf+dUES0tWzauxjLCN7MLgB8QrWl7n7vfNuZ4E9AFnAX0AV9w9544zn2sMj0ZBnODw9vF6ZE1yhepLe3tUxsNj2032a+D0lRO6b7ie8HIr4T+/ihUlzIbv++002DHjug9+/vhjjvgyJHx7QYG4Pbb4fXXYd8+WLgweh4shKSTT47SQ9/73uT/7umoOOCbWRK4GzgfOAg8b2ab3H1fSbOrgP9w9w+Z2WXA94AvVHru6Ui1pmhINjCQi36jKZ8vUl/a26M0zkQ5/FQK5s2Lrg8kErBhQxSoS9M/Y68hfOQjcO65UQ5/1y64//7ouBnccAOsXg1PPhkF+Xwe9u+fuH+PPTby+ne/G33sP/8z+kK4/XY46SR48814PpOiOEb4K4D97v4ygJk9BKwCSgP+KuCWwut/BO4yMyv8/JhV7YvayVyeoWtPdBUmvSxN+12PwSNpuOSSmftqFZFZM9mvg2Lkyeej1/n86Au3qVSU38/louP798P69SPvmU6P/0JZtw6+/vXRo/pyvwam6q234D3viTfoxxHwFwK9JdsHgY9P1Mbdh8zsDaAZGPX9ZmYdQAdAS0tLDF0rb9TUyDfdFH2dwsizgr5I1ZSrsJlK1c1UZTIjgdx95CJs6UXfTAY+8QnYujX6m8HBKDysWDG6D5lMVK3T1wdPPBF9cRSZwZe+BL/+Nbz6Khw8eOx9feutiv6p49RUlY67dwKdAG1tbbMy+s9ufZDMH0GqB9oPAo88Ev0+K15eL9Zg1XuBrkiMphugy1XYQLz18qU5/mQSLroITjklStd0dUXpmqGh8X/3+OPRo7ExGs2vWVM+T1/kDg8+OP1+QnSuOMUR8F8BFpVsn17YV67NQTObB7yX6OJtVWV7s6TOe4XBBDTkIfMAtO/fD2efPbrh8ceHcVeGSAwquaGpXG38gQMjgbVcvfyxfrkUc/zF4P7441FO/777Rgf64sjffeS5mPp5+OHyF3TjFvd0EnFMrfA8sNjMPmhmjcBlwKYxbTYBlxdefw7452rk78fq+n6agSR4AgaS0LVsgobvvquJPESmqJI5cIqj72Qyem5uji6qlkaL5uaR18Uvl29/O3qe6hQJ7e3Q0hIF+GJlzdhRfTI5ct5i6qfYr0svHanhn47ilBENDbBkCSxfDu97Hxx3HCxeDF/7GmzfHv8Ys+IRfiEnfy3wJFFZ5gZ332tmf01UD7oJWA/8nZntB14n+lKovpf/L5w5xbYT3fkhIqNM9aapcsZW2BTz7UW5XJRKWbp0dDXNZHfLTtbP4nsXrVgRPf/rv47s+9jHomxv6S+Ja66J/q74a6D4pVH8cviLv4hKLffujW7Q+r3fg29+c/xNWrMllhy+u28GNo/Z952S10eAP43jXLG56SbSO3OsX8ZwSie9Z4K2n/600jkiUzRZWeRU/r70bxobR+fKx1bTxPHl0twM11038j7r1kUXY0sD/lVXjZ5CoaNjdOCGY7v0V42QUlMXbWdNsTLndDAfeYxTvMz+ox/NehdF5rKp3jQ1lfcp5ttL6+WLATbOL5exo+7i/qPNlzP231nr40KrgVR6WW1tbd7d3T0zb754Mezfz9o/gm9/CnIJSObgu8/AN58rtLnxRpVnitSQOEsz65mZ7XD3tnLHwhzhX3IJ3H47qR5ozMGAQ2MeUpwBt16t/1EiNSiuXw0hCzPgF0bu7XfeyZaud8l8aB6pj3+B9m1K3YhI/QozpSMiUqeOltIJZonDiWR7s6zdtpZsbyBrnIlIsMJM6RRke7Os7Fo5vArWlvQWTZMsInUr6BF+pifDQG6AnOeG58UXEalXQY7ws71ZuvZ08du3f0sykYS85sUXkfoXXMDP9mY594FzGcxHS8w0JBr46se+qrVtRaTuBZXS6dzRyef/5+eHgz3AUH6Ilve2TBjsdVFXROpFMCP8zh2dXP3Tq8ftT1hiwlSOLuqKSD0JZoS/fuf6svuvP/v6CYO4LuqKSD0JJuAf13DcuH0JS3By08kT/k2qNUVjspGkJXVRV0TmvGBSOkvmL2Hrb7aO2teQaDhqEG9f1M6W9BYyPRlSrSmlc0RkTgsm4J956viVTobyZRauHGPUguciInNYMCmdvnfGL6Gb85zy8iISjGACfqo1RUOiYdS+eYl5ysuLSDCCCfjti9p59opnWf3h1SyZv4TVH17N1iu2Kl0jIsGoKIdvZu8D/h5oBXqAz7v7f5RplwNeKGwecPfPVnLe6Wpf1M6jlz1ajVOLiFRdpSP8m4Et7r4Y2FLYLuddd19eeFQl2IuIhK7SgL8K2Fh4vRFYXeH7zYrOHZ38yd/9CZ07OqvdFRGRWVNpWeYH3P1Q4fVvgQ9M0O44M+sGhoDb3P2xco3MrAPoAGhpaamwa+WVTrHw1MtPAdBxVpnl6EVE6sykI3wze9rMflXmsaq0nUdrJU60XuIZhSW3vgSsM7P/Wq6Ru3e6e5u7ty1YsOBY/y1T8vC+h4+6LSJSryYd4bv7eRMdM7N/N7NT3f2QmZ0KvDbBe7xSeH7ZzDLAmcD/mV6XK3PpkkuHR/bFbRGREFSaw98EXF54fTnwk7ENzOy/mFlT4fV84BxgX4Xnnbal7186XI/fkGhg6fuXVqsrIiKzqtKAfxtwvpn9GjivsI2ZtZnZfYU2HwW6zWwP8AxRDr9qAT/TkyHveQDyntedtiISjIou2rp7H7CyzP5u4M8Kr7cDNTOMTrWmSCaS5HN5komk7rQVkWAEc6dtKXfHcaLrzCIiYQgu4Hft6Rpe4nAwP0jXnq4q90hEZHYEFfCzvVl2HtpZ7W6IiFRFMPPhF9en7R/qB8AwGpINpJelq9wzEZHZEcwIP9OToX+onzx5DOP83z+fzOUZzZYpIsEIJuA3n9BMnqgc03GWn7pcwV5EghJMwN91aNeo7b/Z/jdke7NV6o2IyOwLJuCPpeUNRSQ0wQT89LL0qCUOG5ONuulKRIISTJVOcYnDYt19ellaOXwRCUowAR+ioK8gLyKhCialIyISumBG+NneLJmeDM0nNNP3Th+p1pRG+yISlCACfrY3yycf+CRD+SEAEpagKdnElvQWBX0RCUYQKZ0//6c/Hw72EM2DP5AbUFmmiASl7gN+545Odv/77nH7VZYpIqGp+4C/fuf6cftWf2S10jkiEpy6D/invee0UdsnNpzIhR+6UMFeRIJTUcA3sz81s71mljeztqO0u8DMXjKz/WZ2cyXnPFY3nn3jqDts3x58m6t/ejU3PX0Ta7et1Xw6IhKMSqt0fgVcAtw7UQMzSwJ3A+cDB4HnzWzTbC5k/pnFn+Gxlx4bte+Of7mDhCVoTDYqvSMiQah0EfMXAczsaM1WAPvd/eVC24eAVcCMB/zioidHho6MO+Y4Oc8NV+so4ItIvZuNHP5CoLdk+2Bh3zhm1mFm3WbWffjw4YpPXFz0xBm/WHlDooGkJVWtIyLBmHSEb2ZPA6eUOfQtd/9JnJ1x906gE6CtrW18lD5GpYuelEqQ4Kozr6LlvS2641ZEgjFpwHf38yo8xyvAopLt0wv7ZlzfO31l9ycSCc2WKSLBmY2UzvPAYjP7oJk1ApcBm2bhvOw9vLfsfveKfzyIiMw5lZZl/nczOwi0A/9kZk8W9p9mZpsB3H0IuBZ4EngR+Ad3Lx+JY/Zsz7Nl97u7plUQkeBUWqXzKPBomf2vAheVbG8GNldyrjg1JBt0oVZEglPXd9r25/rL7r9y+ZXK34tIcOo64H90wUfH7ZuXmEd6WboKvRERqa66Dvi3rbyNpCWHtz/Z8km2XrFVo3sRCVJdL4DSvqidbVduI9OTUb29iASvrgM+aOFyEZGiuk7piIjICAV8EZFAKOCLiARCAV9EJBAK+CIigVDAFxEJhNXqzJFmdhj4TbX7cQzmA7+rdifmAH1OU6PPaXL6jMo7w90XlDtQswF/rjGzbnefcCF3iehzmhp9TpPTZ3TslNIREQmEAr6ISCAU8OPTWe0OzBH6nKZGn9Pk9BkdI+XwRUQCoRG+iEggFPBFRAKhgB8jM7vDzP6Xmf2bmT1qZidXu0+1xMwuMLOXzGy/md1c7f7UGjNbZGbPmNk+M9trZt+odp9qmZklzWyXmf202n2ZKxTw4/Vz4A/d/b8B/xv4ZpX7UzPMLAncDVwILAG+aGZLqturmjMEXO/uS4BPANfoMzqqbwAvVrsTc4kCfozc/Sl3Hyps/gI4vZr9qTErgP3u/rK7DwAPAauq3Kea4u6H3H1n4fVbRMFsYXV7VZvM7HTgM8B91e7LXKKAP3O+Avys2p2oIQuB3pLtgyiYTcjMWoEzgV9Wtyc1ax1wI5Cvdkfmkrpf4jBuZvY0cEqZQ99y958U2nyL6Of5g7PZN6kPZnYi8DCwxt3frHZ/ao2ZXQy85u47zCxV7f7MJQr4x8jdzzvacTO7ArgYWOm6yaHUK8Ciku3TC/ukhJk1EAX7B939kWr3p0adA3zWzC4CjgPeY2Y/cvf/UeV+1TzdeBUjM7sA+D5wrrsfrnZ/aomZzSO6kL2SKNA/D3zJ3fdWtWM1xMwM2Ai87u5rqt2fuaAwwr/B3S+udl/mAuXw43UXcBLwczPbbWb3VLtDtaJwMfta4Emii5H/oGA/zjnAl4FPFf7/7C6MYkVioRG+iEggNMIXEQmEAr6ISCAU8EVEAqGALyISCAV8EZFAKOCLiARCAV9EJBD/H6MiiH7UUoHXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs=np.real((a0*Gt).sum(axis=1))\n",
    "ys=np.real((a1*Gt).sum(axis=1))\n",
    "\n",
    "plt.plot(xs[:99],ys[:99],'r.')\n",
    "plt.plot(xs[100:199],ys[100:199],'g.')\n",
    "plt.plot(xs[200:],ys[200:],'b.')"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}