Gaussian_kernels_spam.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy \n",
    "from scipy.io import loadmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of X:  (4000, 1899)\n",
      "Shape of y_label:  (4000, 1)\n",
      "Shape of y_label after reshape:  (4000,)\n"
     ]
    }
   ],
   "source": [
    "spamTrain_sample = loadmat('spamTrain.mat')\n",
    "X = spamTrain_sample['X'] # (4000, 1899)\n",
    "y_label = spamTrain_sample['y'] # (4000, 1)\n",
    "y_label = np.int_((y_label - 0.5)*2)\n",
    "\n",
    "num_train = 4000\n",
    "\n",
    "X = X[0:num_train, :]\n",
    "y_label = y_label[0:num_train, :]\n",
    "\n",
    "print(\"Shape of X: \", X.shape)\n",
    "print(\"Shape of y_label: \", y_label.shape)\n",
    "y_label = y_label.reshape(y_label.shape[0])\n",
    "print(\"Shape of y_label after reshape: \", y_label.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4000,)\n"
     ]
    }
   ],
   "source": [
    "labels = y_label\n",
    "print(labels.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1899, 4000)\n",
      "[[0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 1 0]\n",
      " ...\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]]\n"
     ]
    }
   ],
   "source": [
    "points = X.T\n",
    "print(points.shape)\n",
    "print(points)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4000, 4000)\n",
      "[0. 0. 0. ... 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "def gaussian_kernel(M):\n",
    "    pass\n",
    "\n",
    "# points_in_6d = np.apply_along_axis(polynomial_kernel, 0, points)\n",
    "# print(points_in_6d.shape)\n",
    "# print(points_in_6d[:,0])\n",
    "# #for checking\n",
    "# print(2**0.5*points[0,0])\n",
    "# # d = yiyje^(xi - xj)\n",
    "from scipy.spatial import distance_matrix\n",
    "print(distance_matrix(points.T,points.T).shape)\n",
    "norm_matrix = (distance_matrix(points.T,points.T))**2\n",
    "print(np.diag(norm_matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# dataset = np.vstack([points_in_6d, labels])\n",
    "# print(dataset.shape)\n",
    "# print(dataset[:,51])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YiXi = np.multiply(points_in_6d,labels.T)\n",
    "# print(YiXi.shape)\n",
    "# print(YiXi[:,51])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "True\n",
      "[20.87140953+0.j 14.80490113+0.j 12.72439438+0.j ...  1.00001   +0.j\n",
      "  1.00001   +0.j  1.00001   +0.j]\n"
     ]
    }
   ],
   "source": [
    "# D = YiXi.T @YiXi + np.eye(YiXi.shape[1])*1e-5\n",
    "#print(labels.shape)\n",
    "YiYj = labels[:, np.newaxis] @ labels[:, np.newaxis].T\n",
    "#print(YiYj)\n",
    "D = np.multiply(YiYj, np.exp(-1*norm_matrix/2)) + np.eye(norm_matrix.shape[0])*1e-5\n",
    "print(np.argmax(D))\n",
    "print(np.all(np.linalg.eigvals(D) >= 0))\n",
    "print(np.linalg.eigvals(D))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-237.85680994828616"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = 2\n",
    "C = 0.1\n",
    "print(cp.settings.EIGVAL_TOL)\n",
    "cp.settings.EIGVAL_TOL = 1e-08\n",
    "#cp.solvers.options['abstol'] = 1e-10\n",
    "from cvxpy.atoms.affine.wraps import psd_wrap\n",
    "\n",
    "n_feat = X.shape[1]\n",
    "n_train = X.shape[0]\n",
    "\n",
    "\n",
    "lambd = cp.Variable(shape=(n_train,))\n",
    "gamma = cp.Variable(shape=(n_train,))\n",
    "objective = cp.Minimize(-cp.sum(lambd) + 0.5*cp.quad_form(lambd, psd_wrap(D)))\n",
    "\n",
    "constraints= [lambd >= 0, lambd.T @ labels == 0 ,lambd <= C*np.ones(n_train,), lambd+gamma==C*np.ones(n_train)]\n",
    "problem = cp.Problem(objective, constraints)\n",
    "problem.solve(solver = cp.ECOS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3.405054016525355e-13\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "print(lambd.value.T @labels)\n",
    "print(lambd.value[np.where(lambd.value <=0)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([   2,    3,    4, ..., 3996, 3998, 3999], dtype=int64),)\n",
      "[0.1   0.1   0.048 ... 0.1   0.048 0.048]\n"
     ]
    }
   ],
   "source": [
    "lamb2 = np.round(lambd.value,3)\n",
    "k = np.where((lamb2 < C) & (lamb2 > 0))\n",
    "print(k)\n",
    "print(np.round(lambd.value,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2692"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(k[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4000,)\n",
      "(1899, 4000)\n",
      "(5000, 4000)\n",
      "(5000, 4000)\n",
      "(5000,)\n",
      "[-0.0042648  -0.06131075 -0.06978035 ...  0.01313748 -0.00433603\n",
      "  0.01354559]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.kernel_approximation import RBFSampler\n",
    "\n",
    "b = 0\n",
    "rbf_feature = RBFSampler(gamma = 0.5 , random_state = 1, n_components = 5000)\n",
    "print(labels.shape)\n",
    "print(points.shape)\n",
    "points_in_infd = rbf_feature.fit_transform(points.T).T\n",
    "print(points_in_infd.shape)\n",
    "YiXi = np.multiply(points_in_infd,labels.T)\n",
    "j = np.multiply(YiXi, lambd.value.T)\n",
    "print(j.shape)\n",
    "w = np.sum(j,axis = 1)\n",
    "print(w.shape)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.85521481 -0.84396599 -0.92154752 ... -0.99572641 -0.90395484\n",
      " -1.02893131]\n"
     ]
    }
   ],
   "source": [
    "some_num = 100\n",
    "for i in k:\n",
    "     print(labels[i] - w.T @ points_in_infd[:,i])\n",
    "b = labels[k[0]] - w.T @ points_in_infd[:,k[0]]\n",
    "b= max(b)\n",
    "# print(max(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([  46,  152,  153,  188,  242,  244,  248,  257,  357,  401,  430,\n",
      "        463,  528,  573,  606,  711,  732,  792,  834,  859,  884,  907,\n",
      "        922,  951,  955,  958,  979, 1003, 1016, 1082, 1194, 1287, 1333,\n",
      "       1374, 1453, 1465, 1468, 1473, 1549, 1602, 1629, 1630, 1702, 1803,\n",
      "       1820, 1884, 2025, 2079, 2089, 2130, 2156, 2261, 2262, 2283, 2433,\n",
      "       2498, 2528, 2602, 2605, 2616, 2625, 2673, 2744, 2762, 2903, 2936,\n",
      "       2959, 3013, 3022, 3047, 3290, 3291, 3359, 3490, 3576, 3591, 3620,\n",
      "       3653, 3673, 3758, 3782, 3797, 3799, 3800, 3814, 3897], dtype=int64),)\n",
      "86\n"
     ]
    }
   ],
   "source": [
    "print(np.where(w.T @ points_in_infd + b >=0))\n",
    "print(len(np.where(w.T @ points_in_infd + b >=0)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([   0,    1,    2, ..., 3997, 3998, 3999], dtype=int64),)\n",
      "3940\n"
     ]
    }
   ],
   "source": [
    "print(np.where(w.T @ points_in_infd + b <=0))\n",
    "print(len(np.where(w.T @ points_in_infd + b <=0)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3.2212575036992632 3.5730260083918064\n",
      "-3.405434719073713 3.654481663898312\n",
      "(353, 340)\n"
     ]
    }
   ],
   "source": [
    "# h = 0.02\n",
    "# x_min, x_max = points[0, :].min() - 1, points[0, :].max() + 1\n",
    "# y_min, y_max = points[1, :].min() - 1, points[1, :].max() + 1\n",
    "# print(x_min, x_max)\n",
    "# print(y_min,y_max)\n",
    "# xx, yy = np.meshgrid(np.arange(x_min, x_max, h),\n",
    "#                      np.arange(y_min, y_max, h))\n",
    "# print(xx.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def gaussian_hyperplane(x,y,w,b):\n",
    "#     x_feat = rbf_feature.transform([[x,y], [x, y]])[0, :]\n",
    "#     return (w.T @ x_feat + b)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Z = np.ones((xx.shape))\n",
    "# for i in range(xx.shape[0]):\n",
    "#     for j in range(xx.shape[1]):\n",
    "#         Z[i][j] = gaussian_hyperplane(xx[i][j],yy[i][j],w,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlMklEQVR4nO3deXgUZbYG8PekE0ISCGvYCYGArLJIiAZUdkVcQBH3fUauDsyI46iDqKhcBXXcBnUGBByZQXHFbZCrIBGFsIRN9kUQCIZdSCAhS+fcPzrRELN0d1V3VXW/v+fJE5J0qk8l5K2vz/dVlagqiIjIuSKsLoCIiIxhkBMRORyDnIjI4RjkREQOxyAnInK4SCuetHHjxpqUlGTFUxMROdaaNWuOqmpCxc9bEuRJSUnIzMy04qmJiBxLRPZW9nm2VoiIHI5BTkTkcAxyIiKHY5ATETkcg5yIyOEY5EREDscgJyJyOAY5EZHDMciJiByOQU5E5HAMciIih2OQExE5HIOciMjhGORERA7HICcicjjDQS4itUVklYhsEJHNIvKkGYUREZF3zLixRAGAQap6SkSiAHwnIl+o6goTtk1ERDUwHOSqqgBOlX4YVfqmRrdLRETeMaVHLiIuEVkP4DCAr1R1ZSWPGSMimSKSeeTIETOeloiIYFKQq6pbVXsCaAUgVUS6VfKYGaqaoqopCQm/uXcoERH5ydRVK6p6AsASAMPM3C4REVXNjFUrCSJSv/TfMQCGAthmdLtEROQdM1atNAfwloi44DkwvKeqn5uwXSIi8oIZq1a+B9DLhFqIiMgPPLOTiMjhGORERA7HICcicjgGOVEIy8gApkzxvKfQZcaqFSKyoYwMYPBgoLAQqFULWLwYSEuzuioKBI7IiUJUeronxN1uz/v0dKsrokBhkBOFqAEDPCNxl8vzfsAAqyuiQGFrhShEpaV52inp6Z4QZ1sldDHIiUJYWhoDPBywtUJE5HAMciIih2OQE5mI67bJCuyRE5mE67bJKhyRE5mE67bJKgxyohp42y7hum2yClsrRNXwpV3CddtkFQY52V5GhnXhWFm7pLoauG6brMAgJ1uzegKxrF1S9vxsl5AdMcjJ1nwdEZuN7RJyAgY5+SVY7Q47jIjZLiG7Y5CTz4LZ7uCImKhmDHLyWbDbHRwRE1XP8DpyEWktIktEZIuIbBaR+8wojOyL66WJ7MWMEXkxgAdUda2I1AWwRkS+UtUtJmybbIjtDiJ7MRzkqpoNILv037kishVASwAM8hDGdgeRfZh6ir6IJAHoBWBlJV8bIyKZIpJ55MgRM5+WiCismRbkIlIHwIcAxqtqTsWvq+oMVU1R1ZSEhASznpaIKOyZEuQiEgVPiM9V1Y/M2CaRnfA642RnhnvkIiIAZgHYqqovGi+JyF6svkwAUU3MGJH3A3ArgEEisr70bbgJ2yWyBV5nnOzOjFUr3wEQE2ohsiU7XCaAqDo8s5MMsfISs8HCdfNkdwxy8ls49Y6rWjcfDgcysj8GOfnN6kvMWi2cDmRkb7xnJ/kt3K+5wklQsguOyMlv4d475iQo2QWDPEQFq3cbztdcCfcDGdkHgzwEsXcbPOF8ICP7YI88BLF3SxReGOQhKNwnIYnCDVsrIYi9W6LwwiAPUezdkhV4gpQ1GOREZApOsluHPXKiCnjtcf9wkt06HJGHOae+FA5U3RxV+o8nSFmHQR7GnBpagaw73K8fYwQn2a3DIA9jTg2tQNRdUlKCvJx89Oqcj3qRBXCXnEG0nEHS8T1YvyQJqgoRgUQIompFIio6CtGx0YipUxuxdWsjpm4MIiLYqeQkuzUY5GHMVi+FfeiV+FJ3SUkJjh44joN7DuPQ3iM4sv8YjmYdw/GDP+P4wRM4eSQHJ4/mIi8nH6oKADiv7JvdwOwXALxQc/kigjoN4lA/IR71m9ZDw2b10bhlIzRJbIxmSU3QvF0TNE9uhtqx0TVvjMhHUvafN5hSUlI0MzMz6M9Lv2WLHrkfvZKKdRcVFmH/tp+wZ+M+7N2yH/u2HUDW9p/w0w+HUFRQdNb31m1YB41aNECDpvVRv0k84hvWRd2GdRBXLxax8bGovWghar8/D1ElRXBFCFx33Qm59RZoicLtLoG7qBiFZ4pQkFeAM6cLkJebj1M/n0bOsVycPJaLnw+dwPHsEziadQwF+YVnPXfTNglI7NIKbbslIrlnEs7p3Q4t2jfjaJ68IiJrVDWl4uc5Ig9zQX0pXNVRw8deibvYjYTa+9CtwU4sn/UD3hq7G3s370dxkRsA4Ip0oUVyU7Tu1BKpl/VCi/bN0axtEzRLSkBC60aIjqlhVNw1Fvj0tV8PLHddCaR18Xl3VRU5x3JxcM9h/PTDIWTt+An7tx/A3s1ZWLfo+1/qrVM/Dh1T26Nbv07oMaArOp3fHlG1onx+PgpfDHIKjupG3TX0SvJPn8HmZduxcekWbFq2DTtW/4AzeQUAgPhGddGhdzv0HtoDyT3aoO25iWjVsYWxIDRp1k5EUK9xPOo1jkfHPu3P+lpRYRH2bT2AHZk/YPuqXdicsR1znngPqorasdHoPqALUi87D2lXpaBJ68b+7wuFBbZWKDjuvReYPh1Q9VwEZvJkYMKEX79ebrTu7pOKrSt3Ys2XG7Du643YtnIX3MVuRLgikNwzCV37dkTnC85B56hcNNuxDjJwYEjMsOUcz8XGpVuxdtH3WL1wPbJ3HwIAdOyTjP7X9cPAG/qicctGFldJVqqqtcIg94Mt+spOkpHh+WEVlvaLo6OBJUvO+uH9fPgkVi1Yi5UL1mLtV9/j9Mk8REQIzklJRs+B3dB9QFd07dsRsXVjft2mE9dOeklVsX/7T1j+yWos/SADO9fshogg5dIeGH73EKRdmQJXpMvqMinIAtojF5HZAK4AcFhVu5mxTbsK8fwIjPR0T/8bAESAO+8E0tKQvecQvvtoFZZ9vBJblu+AqqJRiwa4aNQF6DOsJ3oNPhd1G9SpeptOXDvp5ShARJDYqSUSO7XEDQ+PRNbObCya8w2+fCsdT476GxJaN8LVf7ocl48Z8uvBjcKWWT3yfwF4FcAck7ZnW07ND0uV64EfjopHuqsz0vs8jJ1rdgMAknsm4dbHR+OCK3ujfa+2EBGftmn92kkvGRgFtOrQHHdMvgG3ThqNFZ+vwfy/L8CMB+fgnWc+xLUPXIWr7xuOmLjaAd4BsitTglxVl4pIkhnbsjsn5ofVTnXujqX3vYRF767Axh9PAa9/i459knH3s7fgolEXoHm7pr5v1ImnEZowCnBFutBvZCr6jUzFtlU78Z/JH+DNR9/BJ68txO+euQlDb+vv3YGQQoppPfLSIP+8qtaKiIwBMAYAEhMTe+/du9eU57WC03vkwai/pKQE67/ehIVvfo1l81eh8EwRWndqicE3XYSBN/ZDi+RmgXliOwtQX27Tsm2Y/sBb2LZqF869uDMemHkvWrZvbkLBZDcBn+ysKcjLc/pkp5MFusd/4shJLJy9BAveWITs3YdQp34cBt10IS65fQDOSUl25mjRzCOfwW1V9e0lJSVYOHsJZjw4B+4iN/7wyp0YdtcgZ/68qUo8ISiMlf/jD1SPf/vqXZg/bQGWvpeBosJidO/fBXc8dT0uvOZ81Kpdy/gTVCYYLy3MPvIZOAOrulIiIiIw/PeDkXpZTzx3x2t48e5/Ysvy7fjj63ejVjRPLgp1DPIQUzHbKv7xv/yygR5/hY273W5kfJqJD178DJuXbUds3RgMv3sIrrz3ErTp0joAe1ehlmAsH7LR7LY3pTRu2QhTFk7Ev594H3Of/hDZew7jyY8fQlx87FmPc3p7kM5m1vLDdwAMANBYRLIATFLVWWZsm7xXWbZV/OM/dszPOcJyGy+Kisai+1/Cux9uwoGd2WjWtgnufekOXHrnwN8ERsAEK2BtNLvtbSkulwt3TL4BrTu1xPN3voa/XjIZU7987JffDZfQhh6zVq3caMZ2yJjKsq2yP36/Xt2np6OwoBgLS9pgXkknHJn6FTqc1xaPzrsfF446Hy5XkE9OCVbA2mh1jK+lDL75IsTUqY2nRr+ASSOfwzNfTESt6Cg7vcggk/DMzhBS1UjL6Mtod7Eb//fYTPzn2c9wBDHoKsdx899+h5TxN5ozmeZvgcHuDzi0H7F47reYeuvfMfz3g3H/jHs4IncwnqIfgirLFTOzRlWx7ONVmDVhLrJ2ZKNTl2a48/w66PW7KyF9+xrbeBmnpIpT6qzCrEfexryp8/HI2+Mx8IZ+Tj0mhT2uWgkxVeWKWZel3bV+D14f/yY2Lt2KxM4t8eT8h5B2VYr5y9mC9TrfaHI5vB9xx1PXY0P6JkwbNxO9BndDWlo9J5VPNeDV7B0qUHcsP3XiNKaNm4mxKQ9j35Ys/On1uzFjwwvoO6JPYNYkl/W6Xa7A9brLjnqPPeZ5n5Hh+zaCUWcAuSJd+PMb9+L0yTzMmfSe1eWQyTgid6hAzPV9+9FKvDpuJk4cPokr770Utz91fdUXrTJLMCYTzRhN22jS019JXVvj8jFDseCNRbjhryPRJDHB6pLIJOyRO5hZfc6c47mYNm4W0uctQ3LPJPz5jXtwTu9ks8q0nsP722Y6vP8obkseh2vuG44xz99mdTnkI052UqU2fLMZU2/5O34+dBK3Pj4a1z88ApFRIfhCjbN7v3jy2r9h03fbMC9rOq9p7jCc7KSzqCreffZjvPnoO2jRvhmmrXgGHc5rZ3VZgRPUm5Pa26AbL8R3H63Epu+2oceArlaXQybgZGcYKsgvwDM3vYxZj7yNi0en4bXVz4Z2iDtdRgYwZYp/k7SV6H1JD7giXVjz1QZTtkfW44g8zOT+fAqPXTUVW5bvwO+n3oLrHryq+tUobElYKwD9/di6MWjXow22rtxpUpFkNQZ5GMk5louHhj6FfVuyMHHe/eg/uoZA4CSh9QK0fj25exusXLDW8HbIHthaCRP5p/LxyPCnsW/rATz58UM1hzgQuMXq5L0ArV9v1q4pfj50EoVnCk3ZHlmLI/IwUFJSgim3/B071+zGpI8eRJ9hvbz7RjMWq4dTayYQ+xqg9ev1GscD8LxKa9yykSnbJOswyMPAe899goxPMzH2lbvQ96o+3n+j0RAJp9ZM+X11uYC77gJuu82c/Q3AipvoGM/NPgrPFJm6XbIGgzzE7dm0D29NehcXj07DiHHDfN+AkRBx+PVJvFI2Ct+379d9dbuB6dOBt96y7cHL7S4BAES42F0NBQzyEPeP+/+FmLox+NNrvw/+/RvLt2YiIz1hl5FhXrBZ3bapOAqPjARKSgBVz5uND155OXkAPCtYyPkY5CFs8/LtWLd4I+554fZfeqJBVdaamTMHmD0beOON345SfQnj8o8FrG/blH/FAQB33+15/+abQHGxrS+udTz7Z0RGuVCnQZzVpZAJGOQh7PPpXyKuXiyGjxliXRFpaZ7AK2s5lB+l+tJDr/jY228HCgo8I+CCAmtGvhUng8t64rfdFtgLxZvgpx8OomlSE0REsLUSChjkIcpd7Lkx8oVXn4+YuNrWFlPV6hdfeugVH3vwoCfEAc/7RhasvKhqMrjivIINJ333bNyHpK6tLK2BzMMgD1F7t2Th9Mk89Bp8rtWlVB14vixvrPjYZs2AiAhPiEdEeO4qbYWaJoMzMoAnnvj11YORvrlJo/qc47nI2pGNobcN8HsbZC8M8hB1YNdBAECbLjYZdVUWeL4sb6z4WMDTb7fB3e2rVDYSLwvxiAjva60Y2iaO6jekbwEAnHtRZ7++n+zHlCAXkWEAXgHgAjBTVaeasV3yX9mqhLj6sRZXUgNfljdWfKydb/RQcSQeEQEMGeL5nDeTuhVD28SlnKv+uwax8THofEEHv76f7MdwkIuIC8BrAIYCyAKwWkQ+VdUtRrdN/ouKjgIAFOaH8CnYdr00bWUj8eho70IcqDy0TbolVHFRMZZ/monzLz/P7+vO22zelmDOiDwVwC5V3Q0AIjIPwAgADHILNWntmfw7+OMRtOnS2uJqTFJZuyGQieLv9suC2NeReJnKQtukU/VXfbEOOcdyMfCGC/36fhvO2xLMCfKWAPaX+zgLwPkmbJcMaNu9DSIiBFszduD84edZXY5xFRPk5ZeB8eMDlyhGEqtiEPsS4oD3q2H88Pn0r9CweQOkXubl9XYqCIeTdZ0oaItIRWSMiGSKSOaRI0eC9bRhKy4+Fp3TzsHyT1dbXYo5KibIhx8G9sqMRq78WBbEkyf7f4BJSwMmTDA1JfdvP4DMhetx+d1D/L7FW4AuxkgGmRHkBwCUf+3eqvRzZ1HVGaqaoqopCQm8e3cwDLrxIuzZuM9+NxDw5443FRNk1KjAJorRxApAEBv17rOfICo6Elf+4VK/t2HGMYrMZ/jmyyISCWAHgMHwBPhqADep6uaqvoc3Xw6OvNx83NzmXnS7qBMmf/JXq8vxMNKycEqP3Iaydmbjd13GY8TYYfjDy3daXQ75KWA3X1bVYhEZB+D/4Fl+OLu6EKfgia0bg+seHIHZE9/G2sUbcZ4dTg4y0mSt2COu7AxKM4PXrqti/DD7kbmoVTsKN0642upSKABM6ZGr6gJVPUdVk1X1aTO2Sea4ZvxwtEhuilfumY78U/lWlxO4JmvZSP+xxzzvTbpRselMvpGyN9Z9vRHffrgS1z80Eg2a1g/a81Lw8Io5IS46Jhp/nnkvsncfxrQ/zoLRVpphgWqyOuG2dDUdbAIQ8oVnCjFt7Ew0a9sEo/9ypWnbJXvhKfohprLuQo/+XXHLY9fi30+9j3N6J2PkuMusLDEwLQuTTpgJqOraSgFaoP3vpz7A/u0/YcrCRxEdE214e04VQtMdlWKQh5DqsuDmx0bhhw0/4h/j30RCq0boNzLV2mLNFqB7W5qquoNNABZob1q2De899zGG3TkQKZf0MLQtfy8bb4dfQ1icxKSqQX/r3bu3kvmeeUbV5fLcnsbl8nxcXt6pfB13wQS9LPoGXbVwnSU1hr3lyz2/mOXLf/v5mBjPLy4m5rdf99HJYzl6U5t79NbksXrq5GlD25o+XTUqSjUioubSTN4NU9T0d+EkADK1kkxljzyE1DSPGBNXG0//dwISu7TCpJHPIeMzLgE1hS+97arWl5s4d+B2u/HsbdNwPPtnTHxnPOLi/b9wWkYGMHYsUFR09j08qmLHqYqwOImpsnQP9BtH5IFT1YCvvJPHcnRs6sN6SeR1+sXsr4NXXCjyZrjqzS/FRDP/+h8dItfqp//4P8PbeuYZz66V3Yg0Ksp5I3LVoP8KAgZVjMgZ5GHqdE6ePnTJUzpErtWZE+aq2+22uiTnWb5cNTLy15SLiPjt6/YgJ9sXs7/WIXKtvvQ/07WkpMTw9srKj4jw7Or06d59TyiEph1VFeSc7AxxVU08xdaNwdOfT8C0sTMxb+p87Nm4Fw/P+SPqNqhjVanOk57+6+3mAM9r94qv24N4lamV/12Dl8b8E70v6YFx0+6CiBjepj9zyCF0HpVzVJbugX7jiDw4vBkMlpSU6MevfqGXRl2vt7T7g25btTP4hTqVN8PVII3I1y/ZpMNjbtR7Ux7S0zl5AXkOsh442Rl+vJl4EhGMGDsMLy59Cu5iN+7r9yjmPfsx3G53sMu1p+omMsuGq//7v8DSpcCYMVU/JoBXmdrwzWY8euUUNG/XFFO+mIjYujGmPwfZm+GLZvmDF80KjurWz1bWcsk5nouX75mBbz9Yga79OuIvs8eiVYfmwV8YbJeFyA5YgJz55QY8cfVzaJqUgOcXT0LDZg2sLokCqKqLZrG1EuIqm3iq7tV+SUmJfvXvb3Rkg9t1eMyNOm/sK1pUO9a/1oA/s152WvZg8wXIi9/+VofVul7H9HxAjx86YXU5FARgayU8VbZsubqWi4hgyC0X441NL6LPsJ6Y+fq3GFvQD5vc9X1bGOzvRazstBDZpguQVRXzps7HlJtfQee0c/DCkifRoEk9q8siCzHIw5A3+dS4RUM88dFDeGLKVchFLdwvAzEVfXC0a2/vnsTfQLZTeNrwLgoF+QV47o5XMeuRtzHghn6YuvBR1KkfZ3VZZDH2yMOUL23o/Gn/wDsvf4EP9sfAVSsKo/58Ba57cMSvk2qVbczMG0gEkl368V7I3nMIk0e/gJ1r9+D2J6/HzY+OMmWJITlHVT1yBjlVr1wgZ0fWw+yL7kX64u2onxCPGydcgyt6xqPW8Eu9n1E1ox6ztumAycwy37yfgZfG/BMA8PCcPyLtyt/Od9mVg46VthewOwRRiCvXImmOk5g4KA6jnn4Gsx55G//487/wfnwUbihoictKdqNWxRNezD4zxOzgDeTJOial1+mTp/H6/f/Cl/9KR6fU9njknfFo3rapOTUGgYOOlY7GHjlVr5KedafUDnh+0SQ8t+hxNE1KwKvaE7fKZXg/oiPy+vQNXC1mT4Ta/G5FKxesxd3dH8CiOd/gpkeuwUvfTnZUiAP2mrsOZRyRh4CqBn+mDAqrOUe716Bz0XPdy1j/2nt4++WFmLE7Bm9f9yYuH7MfI8ZdhoRWjfzdpcqZffOIQF3D3OBI/+iBY/jnA2/hm/cy0KZLKzy+/Gl0Su1gTm1B5oT7fYQC9sgdrqqXrla8pN2+ehfeff4TLPtoJSQiAhdek4qR4y5D136djE/KlR2VGjUCjh2zd8PVzx9+QX4BPnjxc8ybOh/u4hLc9Mg1uO6hEagVHRWEogOHPXLzsEceoqoa/FX8/Jw5gf9j6tinPR5/7wFk7zmET15diIWzv8Y372UgqVtrXH73UAy6+ULEN6zr+4ad1mj1caRfXFSML9/6Bv956n0cyTqGfiP74H/+djuat3NWG6UqvIhW4HFE7nDejMgjIz3XWXW7g5uD+afPYMnb3+HzGV9h55rdiIqOQr+rUzH01v7oPbQ7XJEu7zY0ZYqn3+x2e/rZkyd7znJyuKLCIiz691K8M2U+sncfQqfzO+D3U29Gj/5drS6NbCogyw9FZDSAJwB0BpCqql6lM4PcXDX1yPftA954w9oc3LVuD76YtRhL5i1D7vFTqN+kHvqPTsOA6/uiS9+OiIioZt7daSPyGuQcy8WCmYvx8bQFOPbTz+jQux1ufXw0LriiN9eFU7UCFeSdAZQAmA7gLwxye7JTDhYWFGHVgrX4+p3vsPLzNSg8U4RGLRqg38hU9B3RB937d0FUrUp6wg5vtKoqtmTswIKZi5A+bxkKzxSh56BuuO7BEUi5pAcDnLwS0BOCRCQdDHJbs2MO5uXmY8VnmVj64QpkLlyPgvxCxNaNwXlDu6PPsF5IuaQ7miQmWF2mIfu3H0D6u8uxeO63OLAzG7XjojH4poswYtwwtD23jdXlkcMwyMnWzuQVYN3ijVjxWSZWfbEORw8cBwC07NAcPQd0Rff+XdDtwk5IaN3Y1qPX4qJibFu5EysXrEPGp6uxd0sWRATd+3fBkFv7o//oCxBTh9cLJ//4HeQisghAs0q+NFFVPyl9TDpqCHIRGQNgDAAkJib23rt3r/fVU1hRVfy4eT/WLdqIdV9vxPdLtyAvJx8A0LhlQ3Q6vwPO6Z2MDr3bIblnkqVX/juTV4Bdcz7D5k++xcbcaGzceBB5ufmIcEWg+8Wd0XdEKi4adT4atzR5TT2FJY7IyWd2ace43W7s+X4fNn23DVtWbMf2Vbvw0w+Hfvl6/Sb10KZLKyR2aomWHZqjeXJTNEtqgoTWjVCnfpzhEbyq4vTJPBz88TCyfziEfdsOYN/WLOzesBf7tmWhxO35G2otp9D9qr447+YhOG9Id16VkEzHdeTkEztNkLpcLrTv1Rbte7XFyD9eBsBzN6Mf1v+I3Rv24sdN+7B3axaWzFuGUydOn/W90TG10KBZfdRPiEfdhnUQVz8OMXG1ER1bC7WioxAR6QJUUVKicBcVoyC/EGfyCnD6ZB5yjuXixOEcHM/+Gfmnzpy13SaJjdGuexv0SyhCx2/mo1PJUTSIKAZSLwautckkBOxzMKbAMhTkInI1gGkAEgD8V0TWq+qlplRGlgrizd/9Et+wLnoNOhe9Bp37y+dUFTnHcpG9+xAO7T2KI/uP4uiB4zhx+CROHMnByaO5OLDrIM6cOoOC/EIUFxbDXey5N6lERCAyyoXomFqoHReN2HqxiG9YBx3Oa4sGTXshoVUjNGmTgBbJTdHqnOa/9rkzMoDBs4DCYtudg26ngzEFlqEgV9X5AOabVAvZiK/XyLDDyE9EUK9xPOo1jg/etUkCdb0WE3h7MLbD746MYWuFKuVLPoX8yK+mpLPpOejeHIx9vUE32RODnKrkbT7ZvQ1jiIOPUt4cjKv63Tl4t38RTgciBjkZFtKXKjXxKGV2sHizvZoOxlX97px+cA6FA5EvGORkmI3bxMaZdJQyO1jM2l5VvzunH5ydfiDyFYM8QMLpZR3gf5vY9j8nk45SZgeLmdur7Hfn9IOz0w9EvmKQB0C4vazzl2N+TiZMZpodLMEIKpvO4XrF6QciXzHIAyDcXtb5K5x+TmYHS7gFlT+cfCDyFYM8AMLtZZ2/wu3nVF2w+NNiCqegouoxyAOAoyXv8Ofk4ZgWE9kWgzxAnDZasmrS0Wk/p0AIpxYTBQaDPEQYCWKOCM8W7INauLWYyHwM8hBgNIhDeUToayhbcVBji4mMYpCHAKNBHKojQn9C2aqDGltMZASDPAQYDeJQHRH6E8qhelCj0MYgDwFmBHEojgj9CeVQPahRaDPlVm++4q3eKFhsfwkAIh/wVm8UlkLxlQZRRRFWF0BERMYwyImIHI5BTmEvIwOYMsXznsiJ2COnsMazWikUcEROYa2yteZETmMoyEXkeRHZJiLfi8h8EalvUl1EQVG21tzl4glA5FxGR+RfAeimqt0B7AAwwXhJRMFTdgLQ5Mlsq5BzGeqRq+qX5T5cAeBaY+VQuLDTiTpca05OZ+Zk510A3q3qiyIyBsAYAEhMTDTxaclpOMFIZK4aWysiskhENlXyNqLcYyYCKAYwt6rtqOoMVU1R1ZSEhARzqidH4gQjkblqHJGr6pDqvi4idwC4AsBgteLCLeQ4vMIgkbkMtVZEZBiAhwD0V9U8c0qiUMcrDBKZy2iP/FUA0QC+EhEAWKGq9xiuikIeJxiJzGN01Up7swohIiL/8MxOIiKHY5ATETkcg5yIyOEY5EREDscgJyJyOAY5hTXeVIJCAW8sQWGL13yhUMEROYUtXvOFQgWDnMIWbypBoYKtFQpbvOYLhQoGOYU1XvOFQgFbK0REDscgJyJyOAY5EZHDMciJiByOQU5E5HAMciIih2OQExE5HIOciMjhGORERA7HICcicjgGORGRwxkKchGZLCLfi8h6EflSRFqYVRgREXnH6Ij8eVXtrqo9AXwO4HHjJRERkS8MBbmq5pT7MA6AGiuHiIh8ZfgytiLyNIDbAJwEMLCax40BMAYAEhMTjT4tERGVEtXqB9EisghAs0q+NFFVPyn3uAkAaqvqpJqeNCUlRTMzM32tlYgorInIGlVNqfj5GkfkqjrEy+eYC2ABgBqDnIiIzGN01UqHch+OALDNWDlEROQroz3yqSLSEUAJgL0A7jFeEhER+cJQkKvqKLMKISIi//DMTiIih2OQExE5HIOciMjhGORERA7HICcicjgGORGRwzHIiYgcrsZrrQTkSUWOwHMCkTcaAzgawHKCIRT2AQiN/eA+2Eco7Eew96GNqiZU/KQlQe4LEcms7CIxThIK+wCExn5wH+wjFPbDLvvA1goRkcMxyImIHM4JQT7D6gJMEAr7AITGfnAf7CMU9sMW+2D7HjkREVXPCSNyIiKqBoOciMjhHBHkIjJZRL4XkfUi8qWItLC6Jl+JyPMisq10P+aLSH2ra/KViIwWkc0iUiIili+58pWIDBOR7SKyS0T+anU9vhKR2SJyWEQ2WV2Lv0SktYgsEZEtpf+X7rO6Jn+ISG0RWSUiG0r340lL63FCj1xE4lU1p/TffwLQRVUddTciEbkEwNeqWiwizwKAqj5scVk+EZHO8NwNajqAv6iqY+6gLSIuADsADAWQBWA1gBtVdYulhflARC4GcArAHFXtZnU9/hCR5gCaq+paEakLYA2AkU76PQCAiAiAOFU9JSJRAL4DcJ+qrrCiHkeMyMtCvFQcAPsffSpQ1S9Vtbj0wxUAWllZjz9Udauqbre6Dj+lAtilqrtVtRDAPHjuM+sYqroUwHGr6zBCVbNVdW3pv3MBbAXQ0tqqfKcep0o/jCp9syyXHBHkACAiT4vIfgA3A3jc6noMugvAF1YXEWZaAthf7uMsODBAQomIJAHoBWClxaX4RURcIrIewGEAX6mqZfthmyAXkUUisqmStxEAoKoTVbU1gLkAxllbbeVq2ofSx0wEUAzPftiON/tAZJSI1AHwIYDxFV5xO4aqulW1JzyvrlNFxLJ2l6GbL5tJVYd4+dC5ABYAmBTAcvxS0z6IyB0ArgAwWG06OeHD78FpDgBoXe7jVqWfoyAr7Sl/CGCuqn5kdT1GqeoJEVkCYBgASyaibTMir46IdCj34QgA26yqxV8iMgzAQwCuUtU8q+sJQ6sBdBCRtiJSC8ANAD61uKawUzpJOAvAVlV90ep6/CUiCWUrz0QkBp5JdMtyySmrVj4E0BGeFRN7Adyjqo4aTYnILgDRAI6VfmqFA1feXA1gGoAEACcArFfVSy0tygciMhzAywBcAGar6tPWVuQbEXkHwAB4Lp16CMAkVZ1laVE+EpELAXwLYCM8f88A8IiqLrCuKt+JSHcAb8HzfykCwHuq+pRl9TghyImIqGqOaK0QEVHVGORERA7HICcicjgGORGRwzHIiYgcjkFORORwDHIiIof7f0gEhUT6hzfHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.plot(Pos[0,:], Pos[1,:], 'r.')\n",
    "# plt.plot(Neg[0,:], Neg[1,:], 'b.')\n",
    "# plt.contour(xx, yy, Z,[0])\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}