PGA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1b6fb7fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "from scipy.optimize import minimize, NonlinearConstraint, LinearConstraint\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ff163f34",
   "metadata": {},
   "outputs": [],
   "source": [
    "def arccosh(x):\n",
    "    if x < 1:\n",
    "        x = 1 + np.abs(x)\n",
    "    return np.arccosh(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ee2b63c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hyperboloid_to_ball(X):\n",
    "    x0 = X[:,0]\n",
    "    newX = []\n",
    "    for i in range(X.shape[0]):\n",
    "        newX.append(X[i]/(1+x0[i]))\n",
    "    \n",
    "    newX = np.array(newX)\n",
    "    return newX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e32befa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ball_to_hyperboloid(X):\n",
    "    dim = X.shape[1]\n",
    "    N = X.shape[0]\n",
    "    hyp = []\n",
    "    for ix in range(X.shape[0]):\n",
    "        hyp.append(list())\n",
    "        v = X[ix]\n",
    "        l2 = np.linalg.norm(v,2)\n",
    "        hyp[ix].append(np.sqrt(1 + l2**2)/(1-l2**2))\n",
    "        \n",
    "        for jx,x in enumerate(X[ix]):\n",
    "            hyp[ix].append(2*x/(1-l2**2))\n",
    "    hyp = np.asarray(hyp)\n",
    "    return hyp\n",
    "\n",
    "# def ball_to_hyperboloid(X):\n",
    "#     dim = X.shape[1]\n",
    "#     N = X.shape[0]\n",
    "#     hyp = np.zeros((N,dim))\n",
    "#     l2 = np.linalg.norm(X,2, axis=1)\n",
    "#     hyp[:,0] = np.sqrt(1 + l2**2)\n",
    "#     hyp[:, 1:] = X[:, 1:]\n",
    "#     return hyp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "54f24f78",
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj(gamma, w, dim='2d'):\n",
    "    p1 = np.zeros(gamma.shape[0])\n",
    "    p2 = gamma\n",
    "    \n",
    "    d = lambda p3 : np.abs(np.linalg.norm(np.array([np.cross(p2-p1,p1-p3)]),2))/np.linalg.norm(p2-p1,2)\n",
    "    f = (1/4)*np.sum([arccosh(1 + d(w_i)**2)**2 for w_i in w])\n",
    "    return f\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "910a2f6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isometry_map(X, fmean):\n",
    "    rho = lambda v: arccosh(v[0])\n",
    "    alpha = lambda jx,v: np.sinh(rho(v))/v[jx]\n",
    "    dim = X.shape[1]\n",
    "    N = X.shape[0]\n",
    "    \n",
    "    y = np.zeros((N,dim))\n",
    "    y[0,:] = np.cosh(rho(X[0,:]) - rho(fmean))\n",
    "    for ix in range(1,N):\n",
    "        v = X[ix]\n",
    "        for jx in range(len(v)):\n",
    "            y[ix,jx] = alpha(jx,v)*np.sinh(rho(v)-rho(fmean))\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b3a841c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_mean(points):\n",
    "    dim = points.shape[1]\n",
    "    Q = np.zeros((dim,dim))\n",
    "    Q[0,0] = 1\n",
    "    for i in range(1,dim):\n",
    "        Q[i,i] = -1\n",
    "    dL = lambda x,y: arccosh((x.T).dot(Q).dot(y)) # stability\n",
    "    phi = lambda z, points: np.sum([np.sinh(dL(z,p))**2 for p in points])\n",
    "    \n",
    "    z0 = np.array([1,1,1])\n",
    "    opt = {}\n",
    "    opt['disp'] = True\n",
    "    mean = minimize(phi, z0, args=(points),tol=1e-6, options=opt).x\n",
    "    print('Found Mean')\n",
    "    return mean\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5402ae26",
   "metadata": {},
   "outputs": [],
   "source": [
    "def PGA(points):\n",
    "    points_hyp = ball_to_hyperboloid(points)\n",
    "    fmean = find_mean(points_hyp)\n",
    "    points_hyp_centered = isometry_map(points, fmean)\n",
    "    points_ball_centered = hyperboloid_to_ball(points_hyp_centered)\n",
    "    \n",
    "    \n",
    "    X = points_ball_centered\n",
    "    w = np.zeros(X.shape)\n",
    "    for ix in range(X.shape[0]):\n",
    "        w[ix] = np.sqrt(8)*X[ix]/(1-np.linalg.norm(X[ix],2))\n",
    "    # phi in range [0,pi]\n",
    "    \n",
    "#   Code is given as proof of concept for 2D\n",
    "#   however this is how it would be parametrized for 3D\n",
    "#     if dim == '3d': \n",
    "#         rads = np.deg2rad(np.random.rand(2)*180)\n",
    "#         gamma0 = [np.cos(rads[0])*np.sin(rads[1]), \n",
    "#                  np.sin(rads[0])*np.cos(rads[1]),\n",
    "#                  np.cos(rads[1])]\n",
    "#     elif dim == '2d':\n",
    "    rad = np.deg2rad(np.random.rand(1)*180)\n",
    "    gamma0 = [np.cos(rad), \n",
    "             np.sin(rad)]\n",
    "    gamma0 = np.array([x[0] for x in gamma0])\n",
    "    \n",
    "    # and we do it for PG1 and PG2\n",
    "    pg1 = minimize(obj, gamma0, args=(w), tol=1e-11).x\n",
    "    con = lambda gamma: np.dot(gamma, pg1) \n",
    "    nlc1 = NonlinearConstraint(con, 0, 0)\n",
    "    opt = {}\n",
    "    opt['disp'] = True\n",
    "    pg2 = minimize(obj, gamma0, args=(w), constraints = (nlc1), tol=1e-11, options=opt).x\n",
    "    return pg1, pg2, fmean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0ca6145e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: 0.000084\n",
      "         Iterations: 26\n",
      "         Function evaluations: 413\n",
      "         Gradient evaluations: 102\n",
      "Found Mean\n",
      "[-4.29156971  2.89183681]\n",
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: 67.58389305357818\n",
      "            Iterations: 3\n",
      "            Function evaluations: 9\n",
      "            Gradient evaluations: 3\n",
      "[-4.29156971  2.89183681] [-0.50978785 -0.75653995] 1.509903313490213e-14\n"
     ]
    }
   ],
   "source": [
    "# Sample generated with sampling code on  https://github.com/lateral/frechet-mean-hyperboloid\n",
    "with open('sample.pkl', 'rb') as f:\n",
    "    ppoints = pickle.load(f)\n",
    "\n",
    "pg1, pg2, fmean = PGA(ppoints,'2d')\n",
    "print(pg1, pg2, np.dot(pg1, pg2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4d4d2093",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsG0lEQVR4nO3deXxU1f3/8dcnISFhDbIJhCWyCcg+oN8irXUpbi3YqrXlq/JVy9e6tdaq8FNbv9oFxda6VeqColXRimXRKqhUW7VqEoLsS2SRBNmiYQ2QhPP7497gBLNMYGbuJHk/H495ZO65d+793DOT+cw55y7mnENERCQSSUEHICIi9YeShoiIRExJQ0REIqakISIiEVPSEBGRiClpiIhIxJQ0pM7M7Gkz+03QcdSFmU0zsztqmH+nmf01njEFycycmfUKOo7GxMw2mNmZ/vN6+3lT0giI/wEqMbM9ZrbV/yJuETb/LDP7p5ntNrMiM1tsZreaWZo//yQzm29mO8wsYU62MbMJZlbu79cuP+7zw+a3NLM/+vu/18w+M7OXzexkf34HM3vBzDab2U4ze79i3rFwzl3tnLvb38ZpZlZwFPv2YzPL8fftczN73cxOPdbYEu0LxMzeMbP9/mdvl5nlmtkkM2tah3XEJSnFM/n5/6MH/fd/t18v34rHthOJkkawvuucawEMA0LA7QBmdhHwMvA80N051xb4IZAJdPVfWwq8BFwZ76Aj8B9/vzKAJ4GXzKyN/6WzEBgInA+0AvoBM4Fz/Ne2ALKB4cBxwAzgtfCEGgQz+wXwJ+B3QEegG/BnYGyAYQFgZk1isNrrnHMtgU7ATcAlwD/MzGKwrfrkXv+z3Qp4FHjFzJIDjim+nHN6BPAANgBnhk1PBV4FDNgE3BThenp5b2Otyz3gr3cXkAuMDpt3J14CegbYDSwHQmHzhwKL/Hkv4n3J/6aa7UwA3gubbg44vKR4FfA50LyOdbULGF5FeRpQArTzp28DyoBW/vTdwJ/8508Dv/HjKQEOAXv8R+ea6gBo7S93UQ0xJgGTgE+BIn9dx/nzevh1cDnwGbADuM2fdzZwEO9HwB7gk7BtPunXV6Efe3JYHb8P3O9v6zdAU+A+f/1bgWlAelh8N/vr2gxc4cfTq5p9eQe46oiybsA+4Hx/eiTwH6DYX+/DQKo/71/++vf6+/RDoA3e53s78KX/PPOIz806v+7XA+PD5l0BrPRfNx/vh1SV26nmfbkd2Ahs89/f1rW9L9XUy9OEfe6BZv7rO/vTPfF+FBX563oOyKjqfx7v8/bXoL+HjuahlkYCMLOuwLlAHtAXr0UxK8qbyQaG4P16fx74W0VXl+97eMkgA5iL9yWAmaUCs4Fn/df+DfhBJBv0fwFfhfcPvRY4E5jvnNsbadBmNgRIBfKPnOec2+/vV0UXwbfwvhxGhU2/e8Rr9uK1ajY751r4j83+7CrrAPgvvAT19xpCvR4Y52+zM94X3CNHLHMq3vt7BvArM+vnnHsDr/Xyoh/LYH/Zp/ESYC+8pP0dvLqscDLel2xH4LfAFKAP3nvcC+gC/ArAzM4GfgmcBfTGex/qxDn3GZADjPaLyoEbgXZ49XMGcI2/7Df9ZQb7+/Qi3pf3U0B3vARUwlefsebAg8A5zmvdfANY7M8bC/w/4PtAe+DfwAs1bOdIE/zHt4ET8FqyDx+xzNfel9rqw29dXIaX4LZWFAO/x3v/++H1CtxZ27rqnaCzVmN94P3q2IP3S20jXldHOt4H2AFpYcvO9JfbB1x6xHoiamlUsf0v8f7ZwPtgvxU2rz9Q4j//Jt6vUwub/wE1tzTK/Hh3AB/y1a+rt4ApYcsO8ZfbBayuYl2tgKXA5Br24268L5wmwBbgZ3hfoBWtkLb+ck9XxAycBhQcsZ6a6mA8sKWW+lwJnBE23Qmv9dCEr37Rhv+y/hi4JGzbfw2b1xE4QOWWwo+Af4bV8Wdh8wzv13bPsLL/Atb7z6cfUe99qGNLI+xz+Hg1r/k58Pew6WrXH/bef+k/b+5/Dn4Qvs/+vNeBK8Omk/D+D7pHuJ23gWvCpvtG+r5Usa6ngf1+rCX+8/E1bHsckBc2vQG1NOQYjXPOZTjnujvnrnHOleA1bcH70gHAOXeJcy4Dr4voqPpPzeyXZrbSH1wuxuv+aBe2yJaw5/uANL+l0BkodP4n3bexls196O9XO+fcKc65t/zyoiP2a7G/X9/H614JjzcdmOev6/c1bOtdvCQwDC/BvIn3a/8UIN85V1T9S7+mujooAtrVMnbQHfi7mRX79bsS79d4xxrWX904TXcgBfg8bH1/ATqELbMp7Hl7vK6S3LDl3/DLwXsPw5ev7f2rThfgCwAz62Nmr5rZFjPbhddaalfdC82smZn9xcw2+sv/C8gws2Tntf5+CFzt7/NrZnai/9LuwANh+/UFXpLsEmHMnam8vxvxEsbRvC8A9/mf2WZ4Xa5Tzewcfx87mtlMMyv09/Gv1FAn9ZWSRuJZjdeH/f1ordDMRgO3ABcDbfwP/U68f77afA50OWIAtNtRhvI28B2/O6Ja/oD5bKAA+N9a1vkB3q/HC4B3nXMr/PjO5YiuqTB1PdrsP3i//MfVsMwmvO6VjLBHmnOuMIL1HxnPJn977cLW1co5N6Ca1+zA++U7IGz51s4bsAXvPewatnyd3z+/C3U4XvcQeIPAq4DezrlWeF1INX2ebsJ7n072l6/oWjIA59x859xZeD8qVgGP+/M3Af97RL2mO+c+iDD0zXiJp0I3vJbw1qoXj4zzLMMbWzrPL/4d3vsy0N/H/yay/7F6RUkjwTjnDuH9g/3azH7iH3VkZtabsF9HflkaXn8/ZpZWwyGRLfH+UbYDTczsV3hdP5H4j//aG8wsxcy+jzcIejSewfsC+7t5hwwn+/sQqljAzFLwjhwrAS7366Nazrl9eAP71/JVkvgA71drdUljK9DWzFpHErRzbife+MAjZjbO/9WcYmbnmNm9/mLTgN+aWXd/P9r7/fGR2Ar0MLMkf3ufAwuAP5hZKzNLMrOe1R3e6dfR48D9ZtbB334XMxvjL/ISMMHM+ptZM+DXEcZV0UL4FjAHr+vmH/6slnjdinv8VsFPq9inE8KmW+K9p8Vmdlx4DP4v9LH+j4kDeN22Fe/7NGCymQ3wl23tH11Y3XaO9AJwo5ll+UfgVYwflUVWA9Xz9/tUvIMmKvZxD7DTzLrgHXzQ4ChpJCDnDehdjPdLZRPeL8mXgMfwBqLB+/VUwlcf2BK8VkpV5uN1V6zBa57vp3J3RU2xHMRr9UzA6xr4IfBKXfYnbF378QYkVwCv4Y9lACPw9he8QdDz8QZ+i/1j4vf4raXqvIvXnfNx2HRLvC6QquJYhfdlss7v9ugcQex/AH6BdyTOdrz6uw6vRQTe0WlzgQVmthtvLCfS80sq3tMiM1vkP78M7wfBCrzxp5cJ69qrwq14Bwt86HeNvIX3yx7n3Ot4hwsv9JdZGEFMD/v7sdV/7Szg7LAk/kvgx3hHOz2Od1RduDuBGX79XuyvI52vxrneCFs2Ca9uN+N9xr6Fn4Scc38H7gFm+vu1jK8Oz65qO0eajncQx7/wBq334x20cLRu8T+Pe/ES+1N4XYcA/4fXTboT7/N9VP8nic4qd1WLiIhUTy0NERGJmJKGiIhETElDREQipqQhIiIRi8WFzhJGu3btXI8ePYIOQ0SkXsnNzd3hnGtf1bwGnTR69OhBTk5O0GGIiNQrZlbtVQPUPSUiIhFT0hARkYgpaYiISMQCHdMwswzgCeAkvAt9XYF3WYkX8S5bvAG42Dn3pX/BvAfwLkS3D5jgnFv09bXWrLS0lIKCAvbv3x+NXah30tLSyMzMJCUlJehQRKQeCnog/AHgDefchebd7KcZ3tUy33bOTTGzSXh3Q7sV73ozvf3HyXhX2azzvaMLCgpo2bIlPXr0wBrZnSudcxQVFVFQUEBWVlbQ4YhIPRRY95R/hdFv4t3SEufcQedcMd49l2f4i83gq8tRjwWe8S9J/CHetfhruoBblfbv30/btm0bXcIAMDPatm3baFtZIkGYnVfIqCkLyZr0GqOmLGR2XiRXy09cQY5pZOFdLfQpM8szsyf8SyN39C8NDd7NUSouB96FyldmLaCKG7GY2UQzyzGznO3bt1e54caYMCo05n0XibfZeYVMfmUphcUlOKCwuITJryyt14kjyKTRBO8ywo8654bi3a5yUvgC/t3i6nQZXufcY865kHMu1L59leemiIjExdT5qykpLa9UVlJaztT51d3FIPEFOaZRgHef5o/86ZfxksZWM+vknPvc737a5s8vpPLdxzL9smPyjr1zrKuo5DR3Wq3LJCcnM3DgQMrKyujXrx8zZsygWbNmbN26lRtvvJEPP/yQNm3akJqayi233MIFF1xAUVERF154IdnZ2UyYMIGHH344qnGLSPRtLi6pU3l9EFhLwzm3BdhkZn39ojPwbjgzF7jcL7sc745h+OWX+XesOwXYGdaNVa+kp6ezePFili1bRmpqKtOmTcM5x7hx4/jmN7/JunXryM3NZebMmRQUFADeUU9333039913X8DRi0ikOmek16m8Pgj6PI3rgefMbAkwBO9WjFOAs8xsLXCmPw3ebSbX4d157HHgmrhHGwOjR48mPz+fhQsXkpqaytVXX314Xvfu3bn+eu8mY82bN+fUU08lLS0tqFBFpI5uHtOX9JTkSmXpKcncPKZvNa9IfIEecuucW0zY/aHDnFHFsg7vPtANRllZGa+//jpnn302y5cvZ9iwYUGHJCJRNG6od6zO1Pmr2VxcQueMdG4e0/dweX0U9HkajVJJSQlDhgwBvJbGlVdeybRp0yotc+211/Lee++RmppKdnZ2AFGKSDSMG9qlXieJIylpBKBiTCPcgAEDmDVr1uHpRx55hB07dhAKVdUQExEJRtBjGuI7/fTT2b9/P48++ujhsn379gUYkYjI1zX6lkYkh8jGg5kxe/ZsbrzxRu69917at29P8+bNueeeew4v06NHD3bt2sXBgweZPXs2CxYsoH///gFGLSKNTaNPGkHYs2dPleWdOnVi5syZ1b5uw4YNMYpIRCQy6p4SEZGIKWmIiEjElDRERCRiShoiIhIxJQ0REYmYkoaIiESs0R9ya/8X3ZsSuV/XfvuPo7k0+ptvvsmkSZM4ePAgqampTJ06ldNPPz2qsYuI1EYtjQAczaXR27Vrx7x581i6dCkzZszg0ksvDXgvRKQxavQtjaCNHj2aJUuW1Hpp9KFDhx4uHzBgACUlJRw4cICmTZvGPWYRabzU0ghQxaXRBw4cWKdLo8+aNYthw4YpYYhI3ClpBKDi0uihUIhu3bpx5ZVXfm2Za6+9lsGDBzNixIhK5cuXL+fWW2/lL3/5S7zCFRE5TN1TATjaS6MXFBRwwQUX8Mwzz9CzZ894hSsicphaGgmitkujFxcXc9555zFlyhRGjRoVRIgiImppRHKIbDzUdmn0hx9+mPz8fO666y7uuusuABYsWECHDh2CDFtEGhnzbr3dMIVCIZeTk1OpbOXKlfTr1y+giBKD6kBEamJmuc65Km8bqu4pERGJmJKGiIhETElDREQipqQhIiIRU9IQEZGIKWmIiEjElDRERCRiShoiIhIxJY3aLHkJ7j8J7szw/i556ZhXuWHDBk488UQmTJhAnz59GD9+PG+99RajRo2id+/efPzxx+zdu5crrriCkSNHMnToUObMmXP4taNHj2bYsGEMGzaMDz74AIB33nmH0047jQsvvJATTzyR8ePH05BP3BSRYDT6y4jUaMlLMO8GKC3xpndu8qYBBl18TKvOz8/nb3/7G9OnT2fEiBE8//zzvPfee8ydO5ff/e539O/fn9NPP53p06dTXFzMyJEjOfPMM+nQoQNvvvkmaWlprF27lh/96EdUnPWel5fH8uXL6dy5M6NGjeL999/n1FNPPaY4RUTCKWnU5O27vkoYFUpLvPJjTBpZWVkMHDgQ8K5we8YZZ2BmDBw4kA0bNlBQUMDcuXO57777ANi/fz+fffYZnTt35rrrrmPx4sUkJyezZs2aw+scOXIkmZmZAAwZMoQNGzYoaYhIVClp1GRnQd3K6yD8BkpJSUmHp5OSkigrKyM5OZlZs2bRt2/fSq+788476dixI5988gmHDh0iLS2tynUmJydTVlZ2zHGKiITTmEZNWmfWrTyKxowZw0MPPXR4XCIvLw+AnTt30qlTJ5KSknj22WcpLy+PeSwiIhWUNGpyxq8gJb1yWUq6Vx5jd9xxB6WlpQwaNIgBAwZwxx13AHDNNdcwY8YMBg8ezKpVq2jevHnMYxERqaBLo9dmyUveGMbOAq+Fccavjnk8I2i6NLqI1KSmS6NrTKM2gy6u90lCRCRa1D0lIiIRa5RJoyF3ydWmMe+7iBy7Rpc00tLSKCoqapRfns45ioqKKh2mKyJSF4GPaZhZMpADFDrnzjezLGAm0BbIBS51zh00s6bAM8BwoAj4oXNuQ123l5mZSUFBAdu3b4/aPtQnaWlph08AFBGpq8CTBvAzYCXQyp++B7jfOTfTzKYBVwKP+n+/dM71MrNL/OV+WNeNpaSkkJWVFZ3IRUQamUC7p8wsEzgPeMKfNuB04GV/kRnAOP/5WH8af/4Z/vIiIhInQY9p/Am4BTjkT7cFip1zFde/KAC6+M+7AJsA/Pk7/eUrMbOJZpZjZjmNtQtKRCRWAksaZnY+sM05lxvN9TrnHnPOhZxzofbt20dz1SIijV6QYxqjgO+Z2blAGt6YxgNAhpk18VsTmUChv3wh0BUoMLMmQGu8AXEREYmTwFoazrnJzrlM51wP4BJgoXNuPPBP4EJ/scuBOf7zuf40/vyFrjEeNysiEqCgxzSqcivwCzPLxxuzeNIvfxJo65f/ApgUUHwiIo1WIhxyi3PuHeAd//k6YGQVy+wHLoprYCIiUkkitjRERCRBKWmIiEjElDRERCRiShoiIhIxJQ0REYmYkoaIiERMSUNERCKmpCEiIhFT0hARkYgpaYiISMSUNEREJGJKGiIiEjElDRERiVhCXOVWRESiY3ZeIVPnr2ZzcQmdM9K5eUxfxg3tUvsLI6SkISLSQMzOK2TyK0spKS0HoLC4hMmvLAWIWuJQ95SISAMxdf7qwwmjQklpOVPnr47aNpQ0REQaiM3FJXUqPxpKGiIiDUTnjPQ6lR8NJQ0RkQbi5jF9SU9JrlSWnpLMzWP6Rm0bGggXEWkgKga7dfSUiIhEZNzQLlFNEkdS95SIiERMSUNERCKmpCEiIhFT0hARkYgpaYiISMSUNEREJGJKGhI1+/L3sXriara/sp3yveW1v0BE6h2dpyFR06xXM1qf2prlP1hOUloSbc5sQ7tx7Wj73bakdkgNOjwRiQIlDYmq4y87nj1L9lDwhwKKXi2i6NUiMGj1jVa0G9eOduPa0axXs6DDFJGjZM65oGOImVAo5HJycoIOo9Fx5Y6l5y/lize+qHJ+s/7NvARyQTtahVrFOToRqY2Z5TrnQlXN05iGRJ0lG/1e6Ed6n6qvrLlvxT52/ntnnKMSkWhQ0pCYSMlIYeDcgSS3Tv7avPRe6fR/sb9aGSL1kJKGxEyzvs3o/0L/r33KSvJLyD4pm60vbKUhd4+KNERKGhJTbc9pS897ex6ePu7c4wAo+6KMlT9eyfKLlnNw28GgwhOROlLSkJjL/EUmHS/tSFJaEgNfHUj/F/vTpK134N6OWTvIHpDNtpe3BRyliERCSUNizszo81gfOlzSATOjw8UdGLl8JO3GtQOgdEcpKy5awYofraC0qDTgaEWkJkoaEhfJacn0eazP4enUjqkMeGUA/f7ajyZtvFbHtpnb+HjAx+yYuyOoMEWkFkoaEjdJKZU/bmZGx/EdGbFsBMed5411lG4tZdnYZay8bCWlX6rVIZJoAksaZtbVzP5pZivMbLmZ/cwvP87M3jSztf7fNn65mdmDZpZvZkvMbFhQsUt0Ne3clIHzBtL3qb4kt/IO0d367FayT8qm6PWigKMTkXBBtjTKgJucc/2BU4Brzaw/MAl42znXG3jbnwY4B+jtPyYCj8Y/ZIkVM6PThE6MWDaCNt9pA8DBzQdZeu5SVl21irJdZQFHKCIQYNJwzn3unFvkP98NrAS6AGOBGf5iM4Bx/vOxwDPO8yGQYWad4hu1xFpa1zQGvTGIPn/pQ3ILr9Wx5cktZJ+UzRdvVX1ZEhGJn4QY0zCzHsBQ4COgo3Puc3/WFqCj/7wLsCnsZQV+2ZHrmmhmOWaWs3379tgFLTFjZnSe2JnQ0hAZ384A4MCmAyw5awlrrllD2R61OkSCEnjSMLMWwCzg5865XeHznHe6cJ1OGXbOPeacCznnQu3bt49ipBJv6T3SGfzWYHo91IukZt5HdfOjm8kZlEPxu8XBBifSSAWaNMwsBS9hPOece8Uv3lrR7eT/rTjrqxDoGvbyTL9MGjBLMjKvyyT0SYjWp7YGYP/6/Sw+bTFrf7aW8n262ZNIPAV59JQBTwIrnXN/DJs1F7jcf345MCes/DL/KKpTgJ1h3VjSwDXr1Ywh7wyh5x97kpTmfWwLHywkZ3AOO9/XFXNF4qXWpGFm11cc9hplo4BLgdPNbLH/OBeYApxlZmuBM/1pgH8A64B84HHgmhjEJAnMko2uN3YltDhEq1O8K+SW5JeQNzqP/F/mU16iVodIrNV6EyYz+w1wCbAImA7Md/Xk0qS6CVPD5codm/6wifV3rMcd9D6OzU5sxolPn0irk3XJdZFjcUw3YXLO3Y53bsSTwARgrZn9zsx61vhCkRiyZKPbLd0ILQrRYngLAPat2seibyxi3f9bx6EDhwKOUKRhimhMw29ZbPEfZUAb4GUzuzeGsYnUqvmA5gz7zzB63N0DSzE4BJ/9/jNyQ7nsXrQ76PBEGpxIxjR+Zma5wL3A+8BA59xPgeHAD2Icn0itklKS6HF7D4ZnD6f54OYA7F22l9yRuaz/9XoOHVSrQyRaImlpHAd83zk3xjn3N+dcKYBz7hBwfkyjE6mDFoNbMPzj4XS/ozskA+Ww8a6NLDp5EXuW7Ak6PJEGIZIxjV875zZWM29l9EMSOXpJqUlk3ZXF8I+G02xAMwD2LN5DbiiXjb/dyKEytTpEjkXgZ4SLxELL4S0J5YboNqkbJIErday/fT15/5XH3hV7gw5PpN5S0pAGK6lpEif8/gSGvj+U9L7pAOzO2U3OsBw+m/oZrrxeHDkuklCUNKTBa31Ka0J5ITJ/kQkG7oBj3S3ryBudx741+4IOT6ReUdKQRiE5PZlef+jFkH8NIa1nGgC7/rOLnME5bPrTJtwhtTpEIqGkIY1KxqkZjPhkBF2u966qf2j/IT698VMWn7aYkk9LAo5OJPEpaUijk9w8md4P9mbwPweT1sNrdez8906yB2VT+OdCtTpEaqCkIY1Wm9PaEFoSovPVnQE4tO8Qa69dyydnfULJBrU6RKqipCGNWpOWTejzaB8GLRhE065NASheWEzOwBw2P76ZenJtTpG4UdIQAY476zhGLB3B8VccD0D5nnLWTFzDknOWsL9gf8DRiSQOJQ0RX5PWTTjxyRMZ+OpAUjulAvDl/C/JPimbLTO2qNUhgpKGyNe0Pa8tI5aPoON/dwSgfGc5qyasYtnYZRz4/EDA0YkES0lDpAopbVLo92w/Bvx9ACkdUgAomldE9oBstj6/Va0OabSUNERq0H5ce0YsH0H7H7YHoOzLMlaOX8nyC5dzcNvBgKMTiT8lDZFapLZLZcDMAfR/qT9N2jYBYMcrO8gekM22l7cFHJ1IfClpiESow0UdGLl8JO0uaAdA6Y5SVly0guWXLKe0qDTg6ETiQ0lDpA5SO6YyYNYA+j3XjyZtvFbH9he38/GAj9kxZ0fA0YnEnpKGSB2ZGR1/3JERy0fQ9vy2AJRuLWXZuGWsvGwlpV+q1SENl5KGyFFq2qkpJ809iROfPpHk1skAbH12K9knZVP0j6KAoxOJDSUNkWNgZhx/+fGMWDaCNmPaAHBw80GWnreUVVeuomxnWcARikSXkoZIFKRlpjHo9UH0eawPyS28VseW6VvIHpjNF29+EXB0ItGjpCESJWZG5590JrQ0RMa3MwA4sOkAS76zhDU/XUPZHrU6pP5T0hCJsvQe6Qx+azC9H+5NUjPvX2zztM3kDMzhy3e+DDg6kWOjpCESA5ZkdLm2CyOWjKD16NYA7N+wn0++/Qlrb1hL+d7ygCMUOTpKGiIxlN4znSHvDKHn/T1JSvP+3QofKiRnSA47398ZcHQidaekIRJjlmR0/XlXQotDtDqlFQAl+SXkjc4j/5f5lJeo1SH1h5KGSJw069uMoe8N5YR7TsBSDRwU/KGA3GG57PpoV9DhNQiz8woZNWUhWZNeY9SUhczOKww6pAZHSUMkjizZ6HZLN0J5IVqGWgKwb9U+Fn1jEesmr+PQgUMBR1h/zc4rZPIrSyksLsEBhcUlTH5lqRJHlClpiASgef/mDP3PULJ+k4WlGByCz6Z8Rs7wHHbn7g46vHpp6vzVlJRW7uorKS1n6vzVAUXUMClpiAQkqUkS3W/rzvCc4bQY0gKAfcv3kXtyLut/vZ5DB9XqqIvNxSV1Kpejo6QhErAWg1ow7KNhdP91d6yJQTlsvGsji05exJ4le4IOr97onJFep3I5OkoaIgkgKTWJrDuzGPbRMJoNaAbAnsV7yA3lsvG3GzlUplZHbW4e05f0lORKZekpydw8pm9AETVMShoiCaTlsJaEckN0m9wNksCVOtbfvp68/8pj74q9QYeX0MYN7cLvvz+QLhnpGNAlI53ff38g44Z2CTq0BsWcc0HHEDOhUMjl5OQEHYbIUdn10S5WXr6SktVen7ylGll3Z9H1pq5YsgUcnTRkZpbrnAtVNa/etTTM7GwzW21m+WY2Keh4RGKl1cmtCOWFyLwpEwzcQce6W9eRd2oeB+b/Fe4/Ce7M8P4ueSnocAOnczTio14lDTNLBh4BzgH6Az8ys/7BRiUSO8npyfS6rxdD/z2U9F7egG7annk0+dcvYOcmwHl/593QqBOHztGIn3qVNICRQL5zbp1z7iAwExgbcEwiMdd6VGtCi0N0uaELJ5zxLMkpByovUFoCb98VTHAJQOdoxE99SxpdgE1h0wV+2WFmNtHMcswsZ/v27XENTiSWkpsn0/uB3jTN2FH1AjsL4htQAtE5GvFT35JGrZxzjznnQs65UPv27YMORyTq9jdvW/WM1pnxDSSB6ByN+KlvSaMQ6Bo2nemXiTQKuw7s4tZDe9nLEUc9pqTDGb8KJqgEoHM04qe+JY1soLeZZZlZKnAJMDfgmETi5qb5N/HQ/q38hBJ2pWUABq27wncfhEEXBx1eYHSORvw0CTqAunDOlZnZdcB8IBmY7pxbHnBYInHx+trXeSLvCQC2ZI2mxWVvgdW3332xM25oFyWJOKhXSQPAOfcP4B9BxyEST1+WfMlV864CoEVqC6aPnU6SEoYEQJ86kXrg5/N/zubdmwG476z76JHRI9iApNGqdy0NkcZm7uq5PPPJMwCcdcJZTBw+sdbXzM4rZOr81WwuLqFzRjo3j+mrrhuJCiUNkQRWtK+IifO8JNGqaSue/N6TmNV83amKs6MrTnarODsaUOKQY6buKZEEdv3r17N171YA/jTmT3Rt3bWWV+jsaIktJQ2RBDVrxSxeWPYCAOf2PpcJQyZE9DqdHS2xpKQhkoC27d3G1a9dDUBGWgaPf/fxWrulKujsaIklJQ2RBOOc45rXrmHHPu8aUw+d8xCdW3aO+PU6O1piSQPhIgnmxeUvMmvlLADGnTiO8QPH1+n1FYPdOnpKYkFJQySBbNmzhWv/cS0AbdPbMu28aRF3S4XT2dESK+qeEkkQzjkmzpvIFyVfAPDn8/5MxxYdA45KpDK1NEQSxLNLnmXemnkAXNT/Ii4eEL0LEOpkP4kWJQ2RBFC4q5AbXr8BgA7NO/Dn8/4ctXXrZD+JJnVPiQTMOcdV865i54GdAEw7bxrtmrWL2vp1sp9Ek5KGSMCm503njfw3ABg/cDwX9LsgquvXyX4STUoaIgHaWLyRG+ffCECnFp148JwHo74Nnewn0aSkIRIQ5xxXzr2S3Qd3A/DYdx/juPTjor4dnewn0aSBcJGATMuZxtvr3wZgwpAJnN/n/JhsRyf7STSZc672peqpUCjkcnJygg5D5GvWfbmOQY8OYm/pXjJbZbL0p0vJSMsIOiwRAMws1zkXqmqeuqdE4uyQO8QVc65gb+leAJ747hNKGFJvKGmIxNnDHz/MuxvfBWDisImM6TUm4IhEIqekIRJHa4rWMOmtSQB0b92d+75zX8ARidSNkoZIjG3d4915r/xQOf8z538oKfPOj3hq7FO0bNoyyNBE6kxJQySG8r/IZ+Kr3j2+7//wfj7Y9AEA1424jm9nfTvI0ESOipKGSAzNWTWHuavn8sjHj3D7wtsB6NmmJ1POnBJwZCJHR+dpiMTQnNVzALju9esAMIynxj5F89TmQYYlctTU0hCJkW17t/H+pvcrlXVq2Yl3N77Lc0ueo2hfUUCRiRw9JQ2RGHl1zasccocqlW3evZk7/nkHa4rW0Ca9TUCRiRw9dU+JxEhF11S4ZinNeHrs01w04KIAIhI5dkoaIjGw9+BeFny6oFJZt9bdmHPJHIYcPySYoESiQElDJAYWfLqA/WX7D0+f2u1UZl08iw7NOwQYlcix05iGSAyEd039ZNhPePuyt5UwpEFQS0MkysoOlTFvzTySLZkHzn6Aa0Zcg5kFHZZIVChpSIMyO68w8PtGvPfZewAsuHQBp2edHtdti8SakoY0GLPzCpn8ylJKSssBKCwuYfIrSwHimjg+/eJTsn+SzQltTojbNkXiRWMa0mBMnb/6cMKoUFJaztT5q+MaxxVDr1DCkAZLSUMajM3FJXUqjxWNX0hDpqQhDUbnjPQ6lYtI3SlpSEKbnVfIqCkLyZr0GqOmLGR2XmG1y948pi/pKcmVytJTkrl5TN9YhynSaGggXBJWXQe2K8qCPnpKpCELJGmY2VTgu8BB4FPgf5xzxf68ycCVQDlwg3Nuvl9+NvAAkAw84ZzTDQkauJoGtqtLBOOGdlGSEImhoLqn3gROcs4NAtYAkwHMrD9wCTAAOBv4s5klm1ky8AhwDtAf+JG/rDRgiTKwLSJfCSRpOOcWOOfK/MkPgUz/+VhgpnPugHNuPZAPjPQf+c65dc65g8BMf1lpwDSwLZJ4EmEg/Argdf95F2BT2LwCv6y68q8xs4lmlmNmOdu3b49BuBIvGtgWSTwxG9Mws7eA46uYdZtzbo6/zG1AGfBctLbrnHsMeAwgFAq5aK1X4k8D2yKJJ2ZJwzl3Zk3zzWwCcD5whnOu4su9EOgatlimX0YN5dKAaWBbJLEEdfTU2cAtwLecc/vCZs0FnjezPwKdgd7Ax4ABvc0sCy9ZXAL8OL5RNy6JcOE/EUk8QZ2n8TDQFHjTv+TCh865q51zy83sJWAFXrfVtc65cgAzuw6Yj3fI7XTn3PJgQm/4EuXCfyKSeOyrnqGGJxQKuZycnKDDqHdGTVlIYRWHtXbJSOf9SbrUt0hDZ2a5zrlQVfN0RngVGnvXjM6PEJHqJMIhtwmlomumsLgEx1ddMzVd86ih0fkRIlIdJY0jJMo9GYKk8yNEpDrqnjqCumZ0foSIVE9J4widM9KrHARubF0zOj9CRKqi7qkjqGtGRKR6amkcQV0zIiLVU9KogrpmRESqpu4pERGJmJKGiIhETElDREQipqQhIiIRU9IQEZGINeir3JrZdmBj0HEco3bAjqCDSCCqj8pUH5WpPio72vro7pxrX9WMBp00GgIzy6nuEsWNkeqjMtVHZaqPymJRH+qeEhGRiClpiIhIxJQ0Et9jQQeQYFQflak+KlN9VBb1+tCYhoiIREwtDRERiZiShoiIRExJI4GY2VQzW2VmS8zs72aWETZvspnlm9lqMxsTVn62X5ZvZpMCCTwGzOwiM1tuZofMLHTEvEZVF1VpTPsazsymm9k2M1sWVnacmb1pZmv9v238cjOzB/06WmJmw4KLPPrMrKuZ/dPMVvj/Kz/zy2NbH845PRLkAXwHaOI/vwe4x3/eH/gEaApkAZ8Cyf7jU+AEINVfpn/Q+xGluugH9AXeAUJh5Y2uLqqom0azr1Xs+zeBYcCysLJ7gUn+80lh/zfnAq8DBpwCfBR0/FGui07AMP95S2CN//8R0/pQSyOBOOcWOOfK/MkPgUz/+VhgpnPugHNuPZAPjPQf+c65dc65g8BMf9l6zzm30jm3uopZja4uqtCY9rUS59y/gC+OKB4LzPCfzwDGhZU/4zwfAhlm1ikugcaBc+5z59wi//luYCXQhRjXh5JG4roC71cBeB+ETWHzCvyy6sobMtVF49rXSHR0zn3uP98CdPSfN5p6MrMewFDgI2JcH7pzX5yZ2VvA8VXMus05N8df5jagDHgunrHFWyR1IVIXzjlnZo3qPAIzawHMAn7unNtlZofnxaI+lDTizDl3Zk3zzWwCcD5whvM7IoFCoGvYYpl+GTWUJ7za6qIaDbIu6qimOmiMtppZJ+fc5353yza/vMHXk5ml4CWM55xzr/jFMa0PdU8lEDM7G7gF+J5zbl/YrLnAJWbW1MyygN7Ax0A20NvMsswsFbjEX7YhU100rn2NxFzgcv/55cCcsPLL/KOGTgF2hnXb1HvmNSmeBFY65/4YNiu29RH0EQB6VDoaIh+vz3Gx/5gWNu82vCNmVgPnhJWfi3fUxKd43TqB70eU6uICvD7XA8BWYH5jrYtq6qfR7OsR+/0C8DlQ6n8+rgTaAm8Da4G3gOP8ZQ14xK+jpYQdhdcQHsCpgAOWhH1nnBvr+tBlREREJGLqnhIRkYgpaYiISMSUNEREJGJKGiIiEjElDRERiZiShoiIRExJQ0REIqakIRJHZjbCv5dBmpk19++DcFLQcYlESif3icSZmf0GSAPSgQLn3O8DDkkkYkoaInHmXy8qG9gPfMM5Vx5wSCIRU/eUSPy1BVrg3W0tLeBYROpELQ2RODOzuXh328sCOjnnrgs4JJGI6X4aInFkZpcBpc65580sGfjAzE53zi0MOjaRSKilISIiEdOYhoiIRExJQ0REIqakISIiEVPSEBGRiClpiIhIxJQ0REQkYkoaIiISsf8PK52jimKZeTMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mn = np.mean(ppoints, axis=0)\n",
    "pg1 = pg1/np.abs(np.linalg.norm(pg1))\n",
    "pg2 = pg2/np.abs(np.linalg.norm(pg2))\n",
    "plt.quiver(*pg1, scale=3, color=['m'], label = 'PG1')\n",
    "plt.quiver(*pg2,scale=7, color=['g'], label='PG2')\n",
    "plt.scatter(ppoints[:,0] - mn[0], ppoints[:,1] - mn[1])\n",
    "plt.scatter(fmean[0], fmean[1], label='mean')\n",
    "plt.title('PG1 and PG2 withCentered Dataset on Ball')\n",
    "plt.legend()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0cc6e14",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}