E10_understand the T square_kpca_experiment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "78fe4589",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.decomposition import KernelPCA\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ccebe7ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_normal = pd.read_csv('./data/normal_training.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d496c37",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_faulty = pd.read_csv('./data/faulty_testing.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5b2fa46f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = train_normal[(train_normal.simulationRun>=1)&\n",
    "                                 (train_normal.simulationRun<5)].iloc[:,3:]\n",
    "\n",
    "df_test = train_normal[(train_normal.simulationRun>5)&\n",
    "                                (train_normal.simulationRun<10)].iloc[:,3:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "114d11c2",
   "metadata": {},
   "source": [
    "### Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e1625cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the mean and std\n",
    "mu_train = df_train.mean(axis=0)\n",
    "std_train = df_train.std(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "71290934",
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalizing training data\n",
    "X = np.array(((df_train - mu_train)/std_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1fdf05d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 52)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8edca9c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "kpca = KernelPCA() \n",
    "kpcat = kpca.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "42134e15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 52)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kpcat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "113886e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = kpca.eigenvalues_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "115b8d4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52,)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "63897508",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = kpca.eigenvectors_.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d764c24e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 2000)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "53a867b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "conf_T2 = 0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5abb4ca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# limit of the T^2 statistic\n",
    "from scipy.stats import f\n",
    "F = f.ppf(conf_T2, 30, X.shape[0]-30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "45e969a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "T2_lim = ((30*(X.shape[0]**2-1))/(X.shape[0]*(X.shape[0]-30)))*F"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c824d1",
   "metadata": {},
   "source": [
    "### Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "01a22bcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "IDV = 1\n",
    "df_test = test_faulty[(test_faulty.faultNumber==IDV) & \n",
    "                               (test_faulty.simulationRun==1)].iloc[:,3:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "edb6e3da",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xt = np.array((df_test - mu_train)/std_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9f4cca71",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960, 52)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xt.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "eb5a7d28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-6.20769948e-05, -7.61391924e-04,  7.57745152e-04, ...,\n",
       "         1.78223809e-02,  2.39561614e-02,  2.28408337e-02],\n",
       "       [-4.65880913e-03,  7.21150055e-04, -6.01273347e-04, ...,\n",
       "        -3.33695021e-04,  2.22174848e-03,  2.97980204e-03],\n",
       "       [ 8.58241009e-04, -9.23134036e-03, -1.17745970e-02, ...,\n",
       "         2.30615950e-02,  6.82544410e-03,  1.50406119e-02],\n",
       "       ...,\n",
       "       [ 5.53672185e-03,  4.81433502e-03, -1.24390158e-02, ...,\n",
       "         2.34059981e-02,  4.10398486e-02,  3.35663080e-02],\n",
       "       [ 3.05457324e-02,  1.56026507e-03,  2.92493357e-02, ...,\n",
       "         1.64957748e-02,  3.19910137e-02,  1.44186629e-02],\n",
       "       [ 7.36538976e-03, -2.03821522e-02, -2.22145646e-02, ...,\n",
       "         1.58277009e-02,  2.98979380e-02, -4.86720506e-03]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[:,:52]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0ee3e34c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 30)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[:,:30].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "86cea1f4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 52)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[:,:52].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "14780782",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = Xt@P[:,:30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "86d991eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960, 30)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d2faa2d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "T2 = np.array([T[i,:]@np.linalg.inv(np.diag(L[:30]))@T[i,:].T for i in range (Xt.shape[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "af357b51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960,)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7a5a0ee1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD+CAYAAADRRMnDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq0klEQVR4nO3deZhU9Z3v8fe3qmkURGwB2ZpuQAjD4o2hW0BjEjMuI45ER4wImUkyMSKZODPem7k3JjHGa3IT8zyZO4nPmCBRrk8yirgvDI4L4xoVoYkZQETa1m4akc1iEZBe6nf/qIVTa1dXVdep7v68nqcfun7n1Dm/U6c53/rt5pxDRETEK+B3BkREpPwoOIiISAoFBxERSaHgICIiKRQcREQkhYKDiIikUHAQEZEUCg4iIpKiolQnMrPLgb8ETgbuds49U6pzi4hI9xRUcjCz5Wa228w2JaVfbGZbzazRzG4EcM495py7FlgCLCjkvCIi0rMKrVa6B7jYm2BmQeAOYC4wDVhoZtM8u9wU3S4iImWqoGol59xLZjY+KXkW0OicawIws/uBy8xsC3Ab8JRzbkMuxx8+fLgbPz758CIikk1DQ8Ne59yIQo7RE20OY4HtntetwGzg74ELgKFmNsk5tzTdm81sMbAYoKamhvXr1/dAFkVE+i4zay70GCVrkHbO3Q7cnsN+y8xsJzCvsrKyrudzJiIiyXqiK+sOYJzndXU0LWfOuSedc4uHDh1a1IyJiEhueiI4rAMmm9kEM6sErgae6M4BzGyemS07cOBAD2RPRES6UmhX1hXAa8AUM2s1s2uccx3A9cDTwBbgAefc5u4cVyUHERF/FdpbaWGG9NXA6nyPa2bzgHmTJk3K9xAiIlKAspw+QyUHERF/lWVwUJuDiIi/StaVtTucc08CT9YPGXIt552XuPGqq+Dv/g6OHIFLLkl989e/HvnZuxeuvDJ1+7e+BQsWwPbt8Dd/k7r9O9+BefNg61a47rrU7TfdBBdcAG++CTfckLr9pz+Fc86BV1+F738/dfsvfwlnngnPPQc/+Unq9jvvhClT4Mkn4Z//OXX7738P48bBypXwm9+kbn/oIRg+HO65J/KTbPVqGDQIfv1reOCB1O0vvBD59xe/gFWrEredeCI89VTk9x//GNasSdw+bBg8/HDk9+99D157LXF7dTX8279Ffr/hhshn6PWpT8GyZZHfFy+Gd95J3H7mmZHPD+Cv/xpaWxO3n302/Oxnkd/nz4d9+xK3n38+/PCHkd/nzoWjRxO3X3op/NM/RX5P/rsD/e3pby/ye2/828tDWZYcRETEX+ac8zsPKTwN0tdu27bN7+yIiPQqZtbgnKsv5BhlWXJQg7SIiL/KMjiIiIi/yjI4qLeSiIi/yjI4qFpJRMRfZRkcRETEXwoOIiKSoiyDg9ocRET8VZbBQW0OIiL+KsvgICIi/lJwEBGRFAoOIiKSQsFBRERSlGVwUG8lERF/lWVwUG8lERF/lWVwEBERfyk4iIhICgUHERFJoeAgIiIpFBxERCRFyYKDmU00s7vN7KFSnVNERPJTUHAws+VmttvMNiWlX2xmW82s0cxuBHDONTnnrinkfCIiUhqFlhzuAS72JphZELgDmAtMAxaa2bQCzyMiIiVUUHBwzr0EfJSUPAtojJYU2oD7gctyPaaZLTaz9Wa2fs+ePYVkT0RE8tQTbQ5jge2e163AWDMbZmZLgc+Y2fcyvdk5t8w5V++cqx8xYkQPZE9ERLpSUaoTOef2AUty2dfM5gHzJk2a1LOZEhGRtHqi5LADGOd5XR1NExGRXqIngsM6YLKZTTCzSuBq4InuHEAT74mI+KvQrqwrgNeAKWbWambXOOc6gOuBp4EtwAPOuc3dPK6m7BYR8ZE55/zOQ0b19fVu/fr1fmdDRKRXMbMG51x9IcfQ9BkiIpKiLIODqpVERPxVlsFBDdIiIv4qy+CgkoOIiL/KMjio5CAi4q+yDA4iIuIvBQcREUlRlsFBbQ4iIv4qy+CgNgcREX+VZXAQERF/KTiIiEiKsgwOanMQEfFXWQYHtTmIiPirLIODiIj4S8FBRERSKDiIiEgKBQcREUlRlsFBvZVERPxVlsFBvZVERPxVlsFBRET8peAgIiIpFBxERCSFgoOIiKSoKNWJzGww8GugDXjBOXdvqc4tIiLdU1DJwcyWm9luM9uUlH6xmW01s0YzuzGafAXwkHPuWuBLhZxXRER6VqHVSvcAF3sTzCwI3AHMBaYBC81sGlANbI/u1lngeUVEpAcVFByccy8BHyUlzwIanXNNzrk24H7gMqCVSIAo+LwiItKzeuIhPZbjJQSIBIWxwCPAfDP7DfBkpjeb2WIzW29m6/fs2dMD2RMRka6UrEHaOXcY+Nsc9ltmZjuBeZWVlXU9nzMREUnWEyWHHcA4z+vqaFrONH2GiIi/eiI4rAMmm9kEM6sErgae6M4BNPGeiIi/Cu3KugJ4DZhiZq1mdo1zrgO4Hnga2AI84Jzb3J3jquQgIuKvgtocnHMLM6SvBlbne1wzmwfMmzRpUr6HEBGRApRll1KVHERE/FWWwUFtDiIi/irL4KCSg4iIv8oyOIiIiL/KMjioWklExF9lGRxUrSQi4q+yDA4iIuKvsgwOqlYSEfFXWQYHVSuJiPirLIODiIj4S8FBRERSlGVwUJuDiIi/yjI4qM1BRMRfZRkcRETEXwoOIiKSQsFBRERSKDiIiEiKsgwO6q0kIuKvsgwO6q0kIuKvsgwOIiLiLwUHERFJoeAgIiIpFBxERCSFgoOIiKQoWXAws4lmdreZPVSqc4qISH5yCg5mttzMdpvZpqT0i81sq5k1mtmN2Y7hnGtyzl1TSGZFRKQ0KnLc7x7gX4HfxRLMLAjcAVwItALrzOwJIAj8LOn933DO7S44tyIiUhI5BQfn3EtmNj4peRbQ6JxrAjCz+4HLnHM/Ay7NN0NmthhYDFBTU5PvYUREpACFtDmMBbZ7XrdG09Iys2FmthT4jJl9L9N+zrllzrl651z9iBEjCsieiIjkK9dqpYI55/YBS3LZ18zmAfMmTZrUs5kSEZG0Cik57ADGeV5XR9NERKSXKyQ4rAMmm9kEM6sErgaeKEamNPGeiIi/cu3KugJ4DZhiZq1mdo1zrgO4Hnga2AI84JzbXIxMacpuERF/mXPO7zxkVF9f79avX+93NkREehUza3DO1RdyDE2fISIiKcoyOKhaSUTEX2UZHNQgLSLir7IMDio5iIj4qyyDg0oOIiL+KsvgICIi/lJwEBGRFGUZHNTmICLir7IMDmpzEBHxV1kGBxER8ZeCg4iIpCjL4KA2BxERf5VlcFCbg4iIv8oyOIiIiL8UHEREJIWCg4iIpFBwEBGRFGUZHNRbSUTEX2UZHNRbSUTEX2UZHERExF8VfmdAiquhOcTrTfuYM3EYdbVVadMamkPc9tQW3t55kBMGBPlMTRXXfeH0+P4iIgoOfUhDc4iv3PU6bR1hKisC3PvNOQAs/G0kLWAwrmoQzR8dib/n0LFOnnlrF89v3c39i8/OKUCkC0Ai0rcoOPQhrzfto60jTNhBW0eYmx7dyLbdH9MRdgCEHQmBwau90/F6074uH/Y33P9HHn/zAwAGDogEIAUIkb6npMHBzC4H/hI4GbjbOfdMKc/f1x062k40DhB2sOXDQ91+fzaX/+srvNl6vAfZJ+3hLgNKrJRx6Gg7m3ceZO6M0SyaXdOtfIlI6eUcHMxsOXApsNs5N8OTfjHwKyAI3OWcuy3TMZxzjwGPmVkV8AtAwaEI7lvbwvJXmmjcc7ig4yx7uYkLp49K+7C/bfWWhMAQUzWoMuPxGppD8SqtmJe37aVl32FuvGRqQXkVkZ7Vnd5K9wAXexPMLAjcAcwFpgELzWyamZ1hZquSfk7zvPWm6PukQPetbeH7j27MKzCMOnlgwuuwg6Uvvpt23wfWb0+b/i/PbqWhOZR2250vvpsQGGKWvtTEfWtbuplbESmlnIODc+4l4KOk5FlAo3OuyTnXBtwPXOac2+icuzTpZ7dF/Bx4yjm3oXiX0X+tXNf9h2z1KSfw8LfO4Y6v1GFJ29Zs2ZXysL9vbQsfHUlf5bTn4zauuvO1lPc0NId45q1dGfNw8+ObMgYVEfFfoeMcxgLer5St0bRM/h64ALjSzJak28HMFpvZejNbv2fPngKz1/cNrMh+C8+sHsqkEYMJGBhQWRHgVwtnUldbRV1tFRdOG5mwf9hFGrZjblu9he8/ujHrOTrDjoc3tMZfNzSH+IcV2WN/Z9glnEdEyktJG6Sdc7cDt3exzzIz2wnMq6ysrCtNznqnhuYQG1pSv31/fvJwZk8clnWsQ8x1Xzid57bsijdkw/GG6dtWb2HpS00pxx9cGeRwW2dC2v1vtDB/ZjUAVy19lU6X8rYEgYAxZ+KwnK5TREqv0OCwAxjneV0dTSuIc+5J4Mn6+vprCz1WX/bwhla8VfqnDq7kny6akrY3UKykkC79jLFD+ZOnsfneN5rZ0BLijffTV/v84C+n8evnt9G6/5N4WtgRLz2kCwxVgyoIHemIv55Zc4q6wIqUsUKDwzpgsplNIBIUrgYWFZopM5sHzJs0aVKhh+rT9h46lvC6vrYqr26iC86q4U+tx6uODn3SmTEwXH7mmPg5kqub9h46xv4jbWnfN/eMMTywriUezBqaQyz+3fqijMyOlYqqBlUSOtLW491mG5pDPLyhFQOumFmtICd9Une6sq4AzgOGm1kr8CPn3N1mdj3wNJGurMudc5sLzZRKDvkZPmRg1zulsWh2TUpJIJ0ln58Y74K6aHYNj/2xNSGIbP3wUNpBdpUVAebPrMaAe6O9lMIOnnlrFy+8s4cV1+Y/kC42KvxYe5jkAsvL2/ay7KV3+eyk4UV7iDc0h1i47DXaosWj+9dt59pzJzDkxAEaMS59Ss7BwTm3MEP6amB10XKESg65aGgOsWZLYm+gGWPyn8V22pihWYPDRdNGpoxN+O7cqSy489V4aSDT6Osr6yIP5q1pBuW1dXQ9kC6dWGnhg/1H0waGmPf3HeH9fS082NBaUBCKeb1pH+2eerPOsIu3y1QEjQX145g+ZiihI229NlhoehSBMp0+QyWHrj28oTWlbj+UoUonF9d94fSMXU8DFtmerK62igVn1cRLA+nESg3Z8pduZHa2B5T327tBxsDg1Z5nEEpWNaiSgKVvV+nodPHPwih8epFsn8F9a1t4atPOoledxQYutneEGVARKEpAld6pLIODSg5dS25vCBbY+6eutopZ46vStjVcMHVkxgfE9CyllYumjUxoU5gzcRjBNA/W377yXsLI7HQTCHrPf+eL78ardXIJDAADKgIF945qaA5xy5Obu+yJFcvXJ+1hlr74Lr/9an1e58r0kI4NfIRI1RlQtADxyIbW+MDFto4wj2xoVXDop8pyPQct9pNdQ3OINW8nfsu/9twJBf8n/u7cqVQEE4fFVVYE0pYaYjZ9kH61viWfn8iyr9Yn5KmutoofX35Gyr7J4ySSJxD0joe4b20Lz2YZXJeOGdwyb3pB3+DveL4x44jvbJ59axe3rd4SP0auA/9iD2nH8Yd0TPLAx+WvpHY3zldy3MsWB7t7Tdne39Ac4stLX2Xmrc9w1dJX48eM7RMLiD94dGNJB08Weo29WVmWHCS7Rza00pn0jBpy4oCCj1tXW8XKxWcn9Pzpqt45uQQDiQ3XyRbNruGFrbtTqrBWrjs+TmLVnz5ImEDw0NF2vv/oRvYeOsaat3fnXFqIcQ5uX/NO/Pzd4S3FdPvEUXe+1MTdf3iP9mg12OmnncQ3PjuBKaOGZKw2Sj7Vph0HuG9tC6EjbbQn3fzGPYe5bfWWosxXNWPMUIIW+dwHBI35M6sTqrBieT50tJ27XnmPsHNpS3cx6Xp2xdIeamilvSOMAd4reuP9EPN/8yqzxlfxx+37E9p4gG61H6VbyySX9pTYfGVNew8TdhAMRD6bBWfVZL1vfYk5l+dffA/yVCtdu23bNr+zU3a+/+jGhLmJggHjgetyW4uh2H7w6MaUNodFs2v46V+llhBiYt8Sw0l/emeNr+KPLSGSv5wHog+rYpg1vorvzp2a8UGW/CD55XPv8IfGvUU7v9eAoNEZTv9wbWgOcdWdr6Z8Cchm6qghVFYEWHBWTUIQzPWBmNwTa+qoIRxp60zoaBCwSLD1fhwBg+9cNIU5E4fFzwORebrWeAZYVlYEuGXedG5dtTlrJ4JcfG7ycG644FMpn1ny/fOuZVJfW8WbrQfirxd/biIXTh8Vf8/WDw+xcl0L+4+0Z+xcAZH71tHpCBicP3Vkt7pjl6obtJk1OOe6X5/pPUY5BoeY+vp6t379er+zUXYamkMsWPZa/A/0J5ef4ds02Oke9F+ZXcP/yRIcABb/bn3WuZe6YkRGWXdmeGpnCyjez+y+tS2sXNfCwIoAb7YeoKMz0s5x86WRh9gn7d2rRsqHAedOHs7cGaPjpbVnN3+YdnR6rn76V5Hr66r9JqahOcStT25OGAyZq9jnecsTm2iP/k06Uj//2HUWK9jGHvA3XjKV21ZvYdnLTTh3vCPAIxtas3aWgOMBGpdYesknH0NOHJB1rE1y8K3swQb/YgQHVSv1UgEi/9kqAsaUUUN8y0ddbRU/ufwMfvj4JsJhx4CKAFdEq4eySTdtR3c4SAgMwYBx7bkTOHisI/6tLNPDLuwipa/b17zDhwdTq8WOtYdZ/of3cgoMY6tOZProkxkxZCCHj3Xw+JsfdL/ai0jDcqxxOdceWNncvuYdNkfbg2LtN+l6bMW+yT6wfjsdubS0p2Fm/P619+MPvUyHGVARiFQhFen7aNhFZvh9a+dBXop+dhDpCPDdh/7E9tDRLo+RXGVVSD7S3TfvFPXJ3aCL1YOup5RlyUHVStl5v3UHDf7HRVP49hf97dmVT994b6+bYviff5H4ORT7+Mkqg8aKpKVVY1VRL3seVn4LRuuCBlQE+PrZ49m88yDTR5/MwWMdPNTQ2u1G9nyMHzaIM8edwmPRVQT7m0mnncQFf3Yay//wXkLJ4ZZ50xPa9hqaQ9z06Ea2h45wwdSR/PLqz+R1PlUr9UPJk+ENCFrOaz+Xo0Krl7xiVSlemSYPLETAYOGsmox1xsnVB+VgxEmVDD9pYLdXByy1YpSaytnYqhMZO/QEAHbsP8rOA58QdpHrHjakkr2HEscCXX7mmLwChKqV+pmG5hB3Jj3oak8d1GsDA0Sql9a8vTtj20GuApZ+kN2Nl0ylZthgbnp0Y951ysnn6aqNp662ihWLz867Dr8n7Pm4jT0f5z9IslQWzq5h/szqeI+oZS839UhnAL/sCB1lR5rqLgcpgQHghXf8W7agLMc5SHqvN+1L+VY1YcRJvuSlWOpqq/jxZTOoCETGVwQs8m3pc5OHZ31f0CKD7CorAgQtUkTPNMht0ewaHvzWOUwaMbjg/IZd5rEdXnW1Vdw8bzonDAgQgPh6Gt2V/J5gwAganDAgwOVnjsnjiLlZNLuGh791Dotm1+SV73xUBCJdZ+tqq/j2Fycx5MQBdLdiI1DkzFoPHLM7zvvUCN/OXZYlB42QTi95mgkDlmQZoNZbLJqd2ne8oTnEa+/uo8PztTFW5RD0dCG87gun59TWUVdbxeyJwwpeZzuWj1zU1VZx7zfnxMeN3Lpqc6Rfv8H4YYM5cLS9y2/zl505hv/Y/GF8pPTNlybWUc+aMIwfPrYxUjVRpC6/sSlPYtO8nzywouhVc8kqAsatl81IuIdzJg5j4IAAbe1hAgFjZs0pGWcLjjlhQJAjSWuN5NIVeuqoIRxt7+T9fYldWM3g6lk1bN5xIG0pMN2I/2L5/OThebc5FIPaHHqRv7l7bUJD56erh/L49ef6mKOedd/aFm5+fFN8oNXNl05n0wcHeDDas6a7c/94+73nYuigCg4e7Yh/ezUoaL6hdP3wvRMXJgsA3/mLxPED2cZnxAJQId1vjUjVTvI4lauWvprwYI41sG7eeZB9Hx/jrZ3da8v4/OThOGD66JOzzmjr/cyA+Ay8MV09vT5dHRm49sPHNmZ9iF80bSTnTTktoQODd34sIKEdydvu9OzmD7l3bTOHjnWmO3ReAgYPLjkn7ypjtTn0M3NnjE4IDgvO8mdsQ6mkK1H84NGN8e6A3Z37p662ihXXzkkY4dsZdikPmEC0mmr512YBdGvEeFfnT55OZOV157D0xXd5bsuuhCoUAyoHBOLn66pU5N1+U7QkkasLp43k5W174qWT+Wm6In937tSEuZ5+Pv+/JcyFFduWbeyJ93y5zjeVfG2xktgH+492OYYhYHCzZ9qUWAkr3cSJI4YMjLcjPbVpZ9qgtWLx2WkHsNXVVnHwWEfCwNRscm1097ubq4JDL+L94+2phWzKTfLDoTtz/3R1vNjo2Be37k74VnzG2KEJD5We/A9aV1vFb79an1BKqggYX64fl9cI2tCRtpR6+mzVKpUVAZZ84XSWdFE95w2syft4t1UNquTmxzcmlIa8568IWkFVobH719AcYsXalqydDM73TBiZ/EXj2c0fcmd00FylZ2zOotk1Gf9fZQvSyVWNAYtUSaUb4T4m2mNpw/b9KWNLjMj7srWhlYqCQy/S0BwidKQtZdqA/mT+zGoeWr+d9k4Xn/snX7H/7B/sP5oQHGaMHVryzzddKSkfyfX0t142I37cWOmnalAlmz44kPYbcDbZHo7ebVNGDWHpi++y++An8bmIij1lRF1tFRdMG5nQDdr7jTxdEPLmsa62KmHqjELzdMXMah6MjhkJBowfRz/32HUPGVjBc1t20bjncLzH0pLPT+TgsY6EatLkcQ9+UptDL5HrNAj9QbEXo+lraxj0l8V6kquznHPxMQPp2k1KkZ9sn3tym+HnJg/n99fM7pH71WfbHNRbKZV3GutyH3bf07qqg8/neJmqTHqjYn8+5Sq5OivWGyxTu0kp8pPtc09uM5w7Y3RO7/NLWQYHrQSXas7EYVRWBOJ//H7XR/Y15fofVLJLrs4q5wDf29oMVa3US5Rqql8R6f36bLWSJEpub8hl1lMRkUJo+oxewNvecKw9nLCkpohIT1Bw6AXmTBwWn3vIAQ+u394v17QVkdJRcOgF6mqrOG/KafHX7Z0uYcF5EZFiK1lwMLOpZrbUzB4ys2+V6rx9xfAhAxNel283AhHpC3IKDma23Mx2m9mmpPSLzWyrmTWa2Y3ZjuGc2+KcWwJcBXw2/yz3TzPGDCUYnfa5ssCRwSIiXcm15HAPcLE3wcyCwB3AXGAasNDMppnZGWa2KunntOh7vgT8O7C6aFfQDzQ0h7h11WbCLjKf/y1fmqGurCLSo3Lqyuqce8nMxiclzwIanXNNAGZ2P3CZc+5nwKUZjvME8ISZ/TtwX7p9zGwxsBigpqa8B4mUyiMbWuPTMIfDLu2KZyIixVTIOIexwHbP61Zgdqadzew84ApgIFlKDs65ZcAyiAyCKyB/fUJDc4iV645PBRwGqgZV+pchEekXSjYIzjn3AvBCLvv257mVkifher1pX8K0v0b6tZJFRIqpkOCwAxjneV0dTROP7sy4mG7m1TkThzEgaPEVqDSvkoiUQiHBYR0w2cwmEAkKVwOLipGpvjLxXnen2U438+q3vzgp4wpUIiI9JafgYGYrgPOA4WbWCvzIOXe3mV0PPA0EgeXOuc3FyFRvqVbqqlSQbtqLTA/2huYQO/YfpSIYoLMzceZVzRgqIqWmWVnzlEupoKE5lLAoeWyFqOSper3HSl4isr8s3CIixVOMWVn73fQZDc0h7ni+MWVuokzpmTyyoZVj7YlVQMnqaqv4cv3xZpnOsOPmxzclnKOhOcQvn3snXsLoDDvGnHJiPDB85a7X+edntvKVu17XfEoiUjJlOWV3T1UrZfq2n8sykd5v8BCZ/C5W5goG0zcSNzSH2H3oGGbEF30POxdfxS2Wn2PtYRyRRcm91UmxAOSAT9rDLH3xXX771YK+DIiI5KQsg0NPNUhnWmrzkQ2RhcEB2jrCPJLUNpAcVObPrKYjfLw67gufGpGy/8MbWnkgunC4V0XA4g//WH4ckSLcGWOHMmPsULZ+eIilL77Lmi27EuZQevatXdy3tqXsV5ASkd6vLINDsUoOyfX1samv2zsdQc9DOrnVJfl1clBxRB7ysbaEF9/ZQ0NzKG1pIJk3zbv0ZzBgbPrgAH9qPZD1mp7atFPBQUR6XFm2OTjnnnTOLR46dGjex8hUX++iP50uMu11Q3OI+TOrqQxaxkntYg/xYLTaZ8aYoUwdfXJ8e2fn8TYHb1VQOu2djoej5329aR83XzqdBbNqGHrigITBbplM95xXRKSnlGXJoRjSVSF9sP8o7dFv+51huHdtCw82tLLi2jnc8qUZ8YW/k9sb6mqruPebc3i9aR9Vgyq5ddXm+FxHxvF2gobmUGJbRIC0D/wVa1tYua6FznA0OhuEc+w0NuTEAXl9HiIi3dFng0Nylc0H+4+y59CxlP3aOsLc+eK7vLRtD20dYda9/xEQmaLC2zYwZ+Iwvv3FSfzg0Y3xwABgBjdfOp262irueL4x3hZhwIThJ9G4++OUczqOB41wLCGNUScPZNfBY/HN3vYKEZGeVJbBoVhtDlfMrGbvoWO88M4eVrzRgmXYb9fBT+KljLb2MDc/vonOsCNggBnhsGNARYBvnDOe+99oSXhv2MHKdS1MGTWEqkGVBMwAR0XAeH/f4bzybcDAAQH+4fxPceuqzbS1hwkEjFsv01TdIlIafXIQnLd3UcCMsHMZq20MuO7zE7nrlfcIRz+LbPtm+rSCAQgGjpdU/vzPTuO5Lbtyri6KmTpqCJUVARacVcOi2TXct7YlXt2lhmgRyUUxBsGVZcmhUN72BpwjELD4gz+ZGSx/9X06ww4jWs2TQbbnfGcYOsORdzvncOTejuC15cND0X8jM5HcumpzvLpryqghKjmISEmUZW+lQnl7F1UEI9/ip44aknbfsCM+1iCHzkI5CQYDnDZkYLway4gMcEsnU1VXW0eY5X94r8tR2CIiPaEsg4OZzTOzZQcOZO/zn0msd9GCWTVgxnNbdrEtTcNwynnTpE0bnT6oZHNm9VCmjxlKMGAELNJ+sPhzEwmmOUG2wkXj7o+7HIUtItITyjI4FGOcQ11tFWNPOZGOzsg3745wZOBb0CJjGS6aNpJTByeuqDbkhNRatgNH27t97o8Ot3Hrqs3RRm3j5kunc+MlU3lgyTlcOG0kwYAlBKJYycJIH6AMuLJOU3WLSOn0yTaHmKpBlYkPYYMFZ9VwRXSQ23/e+WrC/gc/6Ug5xgf7P0lJC0TnSgpYJAitez+UUAI4dXAlTXsPRwbcueNrPtfVVrHkC6czYshA9h46xvNbd9PRGWkTibVTpAsODpgxJv9AKSLSXX02ONy3tiXSJdXz1HaeGU/veL4xpwbj5F0MOH/qSJ5/ezdh59jQsj9hYj2I9JYyM8y5hOqg5Km5XfyYkf06O8OYGZ1hl3DeAFoaVERKq08Gh4bmEDc/vilhcjzvSGYgPs9SW2fmXkzeB36sQbkiYOw++Em8e6xzDrPEbq6d0V5SQMJBEkZtR88b2/XKumrGnnJifAR2bAoOAyoHqL1BREqrTwaH15v20ekJDAFg4ojBTBxxUjwtttbCvWtbUt5vkFJkOHfScKpPHcRDDa38V+uB+IPbRfcNGGm7r7Z3Hp+iO3nUNmbxVd/mJy3/GRuIFwxYfAS2iEiplGVwKHSEdNWgSoIBi49dcEDjnsM07jkcGS0dXa/hipnVPLyhlbb2MGZExkOEXbxqx+ulbXuZNb6Kjs7jk+qNPHkguw8diweE86eO5Lm3EqfZdtH8QOIcTbGSQGxtaK/QkTbC0TYIb5uFiEiplGVwKGQ9h4bmELc8EalSik1q560e8q7jkO5hHZtc75YnNqVUOa17P0QwaISj6XsPtyWs+fzFKafFG5m9vA9373rQDc2h+FoSD29ojS8+5C1heKvCRERKpSyDQyEe3tAaf6inm9TOAYei3VOTV3fzrv0wZdQQbntqC+veT1yac/rok+PVSi7suHLWOMaeciJzJg7j9aZ9hNO0csdKDskyLT6UHLRUpSQipdbngkNyFU2sesnrt6+8R82wwfGpKSqi9f+xtoBbL5vBotk1PLjkHG5bvYVlLzfhXGQw24Kzati6a3P8W31yW0FlRSBhPYdsPY2ylRC8JQwRkVLrcxPvJa8Hfcu86axc15KwwpoB504ezh8a9xJ2qRPqVQSMldednVD94/0Wn/w6+fxLX3yX/3x7N865hLWqM+VXJQQRKSZNvJdGXW0VK66dk1JFtGDZa/G2gAEVAebOGM269z+KlxY6HfESRtgd72EUO6b3wZ3pW31s7egX39lDOMeeRiohiEg5KmlwMLPBwIvALc65VT11nnQP85WLz473DLrCUxUUmw4bIt1Hw9Fv+91tBE63drR6GolIb5VTcDCz5cClwG7n3AxP+sXAr4AgcJdz7rYuDvVd4IE885oXb7XNT//qjIR073TY935zDiuvOzvvKp5Y4/LxUc+op5GI9Fq5lhzuAf4V+F0swcyCwB3AhUArsM7MniASKH6W9P5vAJ8G3gJOKCzLufNOV5Fc95+up9C3vzgp7yqe5AFuX64fl1BCySfvaosQEb/kFByccy+Z2fik5FlAo3OuCcDM7gcuc879jEgpI4GZnQcMBqYBR81stXOuWEsopJWpqyhk7ymUj2J2P80W1ERESqGQNoexwHbP61ZgdqadnXM/ADCzrwN7MwUGM1sMLAaoqSlsWcyuuooWeyxBsRqXswU1EZFSKHlvJefcPV1sX2ZmO4F5lZWVdYWe74qZ1SmN0DHl2lNII6RFxG+FBIcdwDjP6+poWsEKmT4jJrlqJraGQ2+gEdIi4rdCgsM6YLKZTSASFK4GFhUjU4VOvAfwyIbWeLfS3lg1U66lGhHpH3JaJtTMVgCvAVPMrNXMrnHOdQDXA08DW4AHnHObi5GpQpcJbWgO8eD67QWvv9zQHOKO5xtpaA51vbOISB+Sa2+lhRnSVwOri5ojCi85vN60L77QT77rL6vHkIj0ZzmVHEqt0JJDrEE3aJHJ8ubn0d6QrseQiEh/UZZzKxVacihGg656DIlIf9bnZmUtJo1SFpHeSLOy9jD1GBKR/qos2xzMbJ6ZLTtw4EDXO4uISNGVZXAotEFaREQKU5bBQURE/FWWwUHVSiIi/irL4KBqJRERf5VlcBAREX+V9TgHM9sDNOf59uHA3iJmp7fpz9ffn68d+vf19+drh+PXX+ucG1HIgco6OBTCzNYXOgikN+vP19+frx369/X352uH4l6/qpVERCSFgoOIiKToy8Fhmd8Z8Fl/vv7+fO3Qv6+/P187FPH6+2ybg4iI5K8vlxxERCRPfTI4mNnFZrbVzBrN7Ea/81NsZjbOzJ43s7fMbLOZ/WM0/VQze9bMtkX/rYqmm5ndHv08/svMZvp7BYUzs6CZ/dHMVkVfTzCztdFrXGlmldH0gdHXjdHt433NeBGY2Slm9pCZvW1mW8zs7P5y783sv0f/5jeZ2QozO6Ev33szW25mu81skyet2/fazL4W3X+bmX0tl3P3ueBgZkHgDmAuMA1YaGbT/M1V0XUA33HOTQPmAN+OXuONwBrn3GRgTfQ1RD6LydGfxcBvSp/lovtHImuXx/wc+Bfn3CQgBFwTTb8GCEXT/yW6X2/3K+A/nHN/BnyayOfQ5++9mY0F/gGod87NAILA1fTte38PcHFSWrfutZmdCvwImA3MAn4UCyhZOef61A9wNvC05/X3gO/5na8evubHgQuBrcDoaNpoYGv09zuBhZ794/v1xh+gOvqf4s+BVUSWCt8LVCT/DQBPA2dHf6+I7md+X0MB1z4UeC/5GvrDvQfGAtuBU6P3chXwF3393gPjgU353mtgIXCnJz1hv0w/fa7kwPE/oJjWaFqfFC0qfwZYC4x0zu2MbvoQGBn9va99Jr8E/hcQjr4eBux3znVEX3uvL37t0e0Hovv3VhOAPcD/i1ar3WVmg+kH9945twP4BdAC7CRyLxvoP/c+prv3Oq+/gb4YHPoNMzsJeBi4wTl30LvNRb4i9LmuaGZ2KbDbOdfgd158UgHMBH7jnPsMcJjj1QpAn773VcBlRALkGGAwqVUu/UpP3uu+GBx2AOM8r6ujaX2KmQ0gEhjudc49Ek3eZWajo9tHA7uj6X3pM/ks8CUzex+4n0jV0q+AU8wstuyt9/ri1x7dPhTYV8oMF1kr0OqcWxt9/RCRYNEf7v0FwHvOuT3OuXbgESJ/D/3l3sd0917n9TfQF4PDOmBytAdDJZEGqyd8zlNRmZkBdwNbnHP/17PpCSDWE+FrRNoiYulfjfZmmAMc8BRLexXn3Pecc9XOufFE7u1/Oue+AjwPXBndLfnaY5/JldH9e+23aufch8B2M5sSTTofeIt+cO+JVCfNMbNB0f8DsWvvF/feo7v3+mngIjOripa+LoqmZed3Y0sPNeBcArwDvAv8wO/89MD1nUukKPlfwJvRn0uI1KeuAbYBzwGnRvc3Ij243gU2Eunt4ft1FOFzOA9YFf19IvAG0Ag8CAyMpp8Qfd0Y3T7R73wX4brPBNZH7/9jQFV/uffA/wbeBjYBvwcG9uV7D6wg0r7STqTUeE0+9xr4RvRzaAT+Npdza4S0iIik6IvVSiIiUiAFBxERSaHgICIiKRQcREQkhYKDiIikUHAQEZEUCg4iIpJCwUFERFL8f8D0j3161Hl/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(T2,'.')\n",
    "plt.axhline(T2_lim,ls='--',c='r');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "df4a2bac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T2: 0.0\n"
     ]
    }
   ],
   "source": [
    "print(f'T2: {(T2>T2_lim).sum()/T2.shape[0]}')"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}