E11_understand the T square_pca.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bc02f0a0",
   "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 PCA\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "43afb804",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_normal = pd.read_csv('./data/normal_training.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "59cf3f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_faulty = pd.read_csv('./data/faulty_testing.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ea138ad",
   "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": "49a7cee7",
   "metadata": {},
   "source": [
    "### Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "65f46bf4",
   "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": "3ad9ecb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalizing training data\n",
    "X = np.array(((df_train - mu_train)/std_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9d2bff73",
   "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": "7db73799",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(n_components = 30) \n",
    "pcat = pca.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a6c47531",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 30)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c26fa193",
   "metadata": {},
   "outputs": [],
   "source": [
    "# eigenvalues\n",
    "L = pca.explained_variance_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6987f949",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30,)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6fab5a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# eigenvectors\n",
    "P = -(pca.components_.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "05ed2989",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 30)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9b750ac1",
   "metadata": {},
   "outputs": [],
   "source": [
    "conf_T2 = 0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "edefce0c",
   "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": 16,
   "id": "9ea45a1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "T2_lim = ((30*(X.shape[0]**2-1))/(X.shape[0]*(X.shape[0]-30)))*F"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ce45b1f",
   "metadata": {},
   "source": [
    "### Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "488c83ae",
   "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": "93d70d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xt = np.array((df_test - mu_train)/std_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "33e4efa6",
   "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": "e185691a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.53978611e-01,  1.29465955e-01,  2.61852697e-01, ...,\n",
       "         7.81167187e-02, -1.08431991e-01, -3.28164026e-03],\n",
       "       [ 2.17065504e-02, -5.34921036e-02,  3.49511112e-01, ...,\n",
       "         1.09921796e-01,  1.25702213e-01, -4.63941299e-02],\n",
       "       [ 1.42669943e-01,  7.17529283e-03, -2.93703184e-02, ...,\n",
       "        -9.86086909e-02, -4.66941543e-01,  3.74569511e-01],\n",
       "       ...,\n",
       "       [-1.48038512e-01, -3.34824103e-01, -2.93763008e-02, ...,\n",
       "        -9.58489756e-02,  7.41743781e-02, -3.70198785e-02],\n",
       "       [ 5.46128898e-02, -1.27816766e-01,  4.11518118e-01, ...,\n",
       "         2.79922471e-02, -1.27731718e-01, -4.98397042e-02],\n",
       "       [ 7.94874362e-03,  1.85934228e-02,  5.41205692e-02, ...,\n",
       "        -6.47605390e-05,  2.60086376e-02, -9.76170096e-03]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[:,:30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2e371bd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 30)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[:,:30].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8ddff89f",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = Xt@P[:,:30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "2b04ff69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960, 30)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "eb374995",
   "metadata": {},
   "outputs": [],
   "source": [
    "T2 = np.array([T[i,:]@np.linalg.inv(np.diag(L[:52]))@T[i,:].T for i in range (Xt.shape[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b3272271",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960,)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "612b1787",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA07UlEQVR4nO3de3RUVZ4v8O/vVCWBYIgxIK8QIIA0JLQMiQKtjdo+Rrziu0WxHzOtIjP2nfGOM3dsWxmHtp3ue9sZp+9iGtD2uqZH8IUvvNKt0oitLZEUrZKAaIgkhEeAUITIK0mdff84j5xz6pFKqpKqSn0/a7FIvU9VJft39m/v/duilAIREWUnLdUHQEREqcMgQESUxRgEiIiyGIMAEVEWYxAgIspi/lQfAACMGDFCTZw4MdWHQUSUUQKBwBGl1MhEniMtgsDEiRNRU1OT6sMgIsooItKY6HMwHURElMUYBIiIshiDABFRFmMQICLKYgwCRERZjEGAiCiLMQhkuEBjECs21SPQGEz1oRBRBkqLdQLUN2uqm7DstVroSiHXr+HZu+aickJRqg+LiDIIewIZKtAYxMOv1aJLV9AV0NGlY0tDa6oPi4gyDINAhlq1eTdCeveGQLoCNu86xLQQEfUKg0AGCjQG8faOlrDrP9oTxG2rP2QgIKK4MQhkoHXbmhFtU9DOkGJaiIjixiCQgepb2mPe/kUPtxMRWRgEMkygMYite2KnezbuDE8VERFFwiCQYR5+dXtYKkgT9+X2MyH87M2dA3ZMRJS5GAQyyJrqJuw44E71TBk5DC8u/Qby/O6v8rd1Bwfy0IgoQyU9CIjIdBFZKSIvichfJfv5s9mG2gNh1/3g4jJUTijCX35jouv60nPyB+ioiCiTxRUERORpETkkIrWe668WkV0iUi8iDwCAUmqnUmopgFsBXJT8Q85e5WOGuy7fMGssFs8pBQA8cM10zJ86wr7tvS+O4Jp/f4/TRYkopnh7As8AuNp5hYj4AKwAsADADAC3i8gM87brAPw/AG8m7UizXKAxiCffb7Av+zTgu/Mmuu7TdqrTdXnHgXbcuorrBogouriCgFLqPQBHPVdfCKBeKdWglOoA8ByA6837v66UWgDgjmjPKSJLRKRGRGoOHz7ct6PPIsYK4e7LIR1h6wFGDR8S9riQrvDytub+PjwiylCJjAmMA7DXcbkZwDgRuVREfikiqxCjJ6CUWq2UqlJKVY0cOTKBw8gOLcdPuy6LAHPLil3X3XPJZHgmCgFA1IVlRERJHxhWSr2rlPobpdQ9SqkVyX7+bDXP2+B/syysYmjlhCL89MaZYYFgeB6LxRJRZIkEgX0Axjsul5jXUZIFGoN4+oMv7cs+DbiyfHTE+y6eUxoWCFa+14A11U39fJRElIkSCQJbAUwVkUkikgvgNgCv9+YJRGShiKxua2tL4DAGvy0NregMOSqGRhgPcAqe7AhLAT3tGFQmIrLEO0V0LYAPAUwTkWYRuVMp1QXghwB+B2AngBeUUnW9eXGl1Hql1JLCwsLeHndWmVtWDJ+v+9w+x6+FjQd47+9dRdxw5ARnCRFRmHhnB92ulBqjlMpRSpUopX5tXv+mUuo8pdRkpdRP+/dQs9eug+3Qzb0DNAEeWVgecwexyglFePSGma7rdAX8fANLSRCRG8tGpLlAYxAPvbod1v4xugLq9vecPls8pxQXTnQHio/2BDk2QEQuDAJpbt22ZuieBH+8Uz7Pzs8Nu+75rQwCRNQtpUGAA8M980731AS4eXZJXI8dUZAXdt2OA8c5NkBEtpQGAQ4M9+ym2SXINSuE+jTBozfMjDke4HTz7BLk+txhRNfDdx4LNAaxYlM9gwNRFuIqojRXOaEIa++eiy0NrZhbVhx3ALAe+8h1FXj41e2wZ5h6VhoHGoO446ktONOpQxPg8umjcM8lk3v1OkSUuRgEMkDlhKI+N8rBkx1wLDFASDdmG1nPt6WhFWc6dSgAIQW8taMFb+1owTnDcnFrZQkeuGZ6Et4BEaUrBoFBzloz4Bxc/uXGzwEYAaLdU3nUcvREB1a+14AdB47jP++cMxCHSkQpIEqlvrxYVVWVqqmpSfVhDFo/e3MnVr7X9xXDV81giogoHYlIQClVlchzcHZQmgs0BvHgK9vx41e293ng9oFrpmP66II+H8NbO1pwx1NbOHBMNAhxdlAaCzQGcfvqD7GmugnPVjfh9if73hDn+hP7qs906jHrFRFRZuJisTTmLRzX2dX3hnjRBaUJHYsCUBRh8RkRZTYGgTTW28JxsSyeU4ql88sSOp53dx1K6PFElH4YBNKc9QX54igc15P2M1093mdicT4eu3FmxIDxzs4WjgsQDTKcIprGtjS0osuc26lgTOlMRE/zwHya4PFbZ7kCjXNWka6MvY5Xfy+hyQhElEbYE0hjRfm5ruqhiebkI5WRcPrW1851BYAHrpkeVomUvQGiwYVTRNNY8GSHXUBOQ+I9gcoJRVi7ZB7OLwmfjeX3CZZeMjns+qmj3FNLdWVUNiWiwYFTRNPY3LJi5OVo8AmQm9P3QWGnyglFKB8X/nnfWjU+4njDTbNL4Pf0Hl6o2cveANEgwXRQGqucUIRn75qLv7tqGp69a27SVuzePLsEzmUDuX4tannqyglFWH5dheu6rpDCum3NrD5KNAhwYDjNJVI8LtZzPn/PN7BuWzMExtl+rNeojbCT2UcNrXixZi+6Qgo5fg1r705ekCKigcMgkMYCjcE+lZCOR2+CS6Sh5PrDJ+yfO7p0zhoiylBMB6Upq87/42/tSnndnkjjAl4bPzvEtBBRBmIQSFNbGlrR0aVDV8aZdirr9lROKMKiqvERewSWkK6ydtZQpo6NZOpxU3KlNB0kIgsBLJwyZUoqDyMtJXuNQKJuml2C57c2oUuPfp811U346Muj+MFFk7B4jlGrqD9TWukg0BjE7U9uQWeXnlFjI1ZPs6NLR65fS+rEA8osKQ0CSqn1ANZXVVXdncrjSEfWGgGF5KwRSFTlhCJMHHEW6g99FfN+9Ye+woOvbMd/bdmD78ydiGWv16IrpOD3CZ5fMq/PDU2qgkmgMYh125pxpP0MRhbkhQ2iv7ytGR1mZOzo0vHytuaMaEwj9TQz4bgp+TgwnKasNQLWGWYy1ggkqmzEsB6DgGXHgXY8+Mp2+3JXSGHV5t2455LJvW7M4zlr7Y8gEWgMYtGqP7p6Py8Gml1n+95SHMrx2HTuAaVbTzMV0v07GigMAmnKWiOQTr+k91wyGRs/a0EoRkoolnd3HcKmXYcQ0hX8muDbVeN7nJ4KuM9aOyOcta6pbsKy12qhK5XU1Ma6bc1h6S/v6988uwQv1exFZ0hBE2DzrkO46Oe/x6HjpxHSk3s8yZRuPU2n/grozufM1DRef2AQSGP9sUYgUT5NQ0jvbhkFgKYJQnpP5emADsfeCB0hhWerm/DcR034yQ0z7TGESOaWFSPXH7lXFGgMYtlrtXahvWSlNgKNQdTtC18f4X19qxTHys278faOFuw7dtp1/0hBKx1E62km0gAno/GOJ6D39nW8z7ns2nI8v7UpI9N4/YFBgOK2paEVXSF3AMjxCUIJ7FMdUsCPXzXSRs7B5JWbd+PQ8dNYdEEpFs8pxbN3zbUXt3mPyRmANJE+N2jW/Yvyc7H8jTqc6ex+r5oAV0yPvNdy5YQiHDp+2vt0AACfT0NRfi5WbKqPeRx9aUATaXQj9TSttNuZTh0+TbD8+oqYwdmpN70xa5zFu1AxnoDe2wFt73Oe6dTx8Gu1YSctqd9pPXUYBChuzjNyn5nOAYBnq5sSel6lgGWv1WKauQ+yMw//SbMRIKaNLrAHYddta7b/+NtPdULEeA6/2XBVTijqdYrI2sqzM6QgAjtfLgAunjoC911xXswz0rwo23d2dun22IjfJ1gUIQUWKzURraFPpMG2eHuaWxpacaZThwLQpSs85AnO0URqaJ945/Oon9ntqz+0e4XOMZZYAd15jLFSg84AUz620JzR1v2cCggLALk+wfA8P77762osqBjTq8C3ofZArx6TjhgE0li6DVxFOntcU91k55YTEdKVvRbCm4d/fmsTriofbf/xn+nUcf8LH6P0nHy898UR+37Xfn0Mgic78LM3d2L1HxrshjyeFNG6bc12w+Ts2CgA5WOG22fK1nvfdbDdDjJ+TVwNjZPz2i4zBfb81r2uRvvnG3ZGTE3ECg7eBtsKon35PbEazvqWdji/TF0BD726HXX722KO3XgbbwXgg/oj2LrnKJ69ay4A2DOsGo6ccKUFnQ15UX6u6/O66+JJYQ38vmOn4PdpCIUipwa9A/k9mXLuWbjia+fa+2b8wfx96qlRX1PdZAd372PS7e+2JwwCaaq/BjsT5Tx7DDQGsfyNOtcfbl8DggLQfqoTV5aPhk+MNJHl031tmGf2QqyGb0/rSexpPel6jlc/3h/xuSOdUXrVt7RHve3JPzSgtHgYHllfZ/SCfAKlus8onY1aPJyN9q6D7fhoj3uxlvVs0aafBhqD+HjvMU+DrVyBLlYPImyA1HFm7qUro6fnnRXlZI0vdHTqgNmLsoL1ys27sXnXoajP72zIvYPVBUNzXMdtpYH8muC2C0vDAlOkgfyezJl0DuoOHHddt6H2ABbPKY3ZmG+oPRD1MVaPMscnWJvAtOiBwiCQhvprsDPZrK45YPzRXjR1BBZUjMEjr9f22DD6NGBS8TBXDaIn3/8SV5aPxm0XlrpSTEoZDfFPbpiJX278HAePn4n7GEVgp4iiWVPdFNYQO4UU8PQHX9rvtauXjX4kVqPtXQkugF3R1fsq1V8eNU4OzLUXzsfkegZ3I+XNI12/paE1riDW0aXj3mcD+JvLz8O00QX22EnwZAfmlhXbz9V+qtM+q1YANu5sQbQho+FD/HhgwXRXYBMzCvjN92OlXIbm+LrTQCEFBdg9USslE7uwSWQFeX6c7gy5rltQMQZrqpvssQOfJrj74kk4fqbLXi9SPma43QMAjN4i4O5RdoQUVm7ejSfTvKYWVwynoXhyo+nAO2vHygFPG12Af3zpE1cD77XoglII3IXoQrrC95+uxhXTRyHXJ67GKaSAV/7UjCNfxR8AACOAvPqnZjtV4hz8tRqw57f2PKZx8Ngp12XNMW4w4Zx8NB49GeFR3c7Jz8HRk532ZZ9mfKdF+bmuxuQec2/nFZvqUTG2EH6tOz1mLcTzco5ZBBqDeOKdz+0eU0dn9wmEa4GYmbe3Gi8v5/uzP4PjZ/DgK9vhc8wGExgBfXZpEY6e6EDbqU7XY3SFsJ6d5fjpLix7rRZNrSfw1AdfugJbZ5eOn6yvw8fN3TO0NLOVVwBerNmL4Xl+VxpnXNHQqD3R/FwfTnaEwq53bp8KwN5J76FXt9vvP6SrsPvl+ATzp46wv7tnPtyDK8tHhwWid3a04MevbI9rKnSqiEpgZkeyVFVVqZqamlQfRtqwz9g6dWh9HPQbKLHSDla32KcBELH/yHPN/DYALFr9YcQz6xtmjcXWPUddUy4TGXvI8Qn++boKLH+jDqfNWT8CIC9Hw9fHFcbsCXhp5kC0Nw02eeQwfNl6ImwdhV8zBik/cTRo55cU4rUfXhw2E2ra6ALXgO/s0rN7PLYbZo1F6wkjlfKH+iNhZ96P3TgT00YXYN22Zjy/tanHdR4C4Kc3zsQrf2rG1l58LtFYjWWyW5rRw/Pi6hVeMNEYb3hrR0tSX9+akAAYvxP3XzUNRfm5EQN1ri/+dTG9OwYJKKUS6mowCKSpTBtcisT5HgDYszYK8vyoO3AcCyrGAEDEFM85+Tn46kxX3Pn2c4bl4uiJ2Auevuk4c7NYjXfDkRPQlfHHXDWhqFdBwXLBxCI8sGA6Vm7eje3Nx9By/Iyd3x5z9hBXQFs6vwy7j5zA7z87BGUOLl867Vy0HD+NT5vb7AZTA9DHtXnG4wW47cJSvLyt2Q5+PcnxCe68aFLY2W+6Gef5TKO5Y44xfnCbeVISqZeTKJ8Y6amaxmDM585N8jgBgwBlHOesCiD62b0ArqmasQiA2+eU4rnqpqgNpgjwzSkjXLOJvDQBlnyzu3GOZwGc9zV+esNMe42BQvj7EwGuP38s3tx+oNcDyn2hCTBzXKErsPREAIyKcJYtMAJpvOMIyeZNK4n5ryDfj7aTXa77Wr01n2NaLmCkWj/Ze6xPvYKzcn042RkyThZgPHfInFKsEH9guXLGKMwaf3ZSTvCSEQRYSpoGlHdWRbS/m940MT5NUDG2EJUTo/9BKYWYAQAw/ohX/6EBb+9osQcE508dEf+AozLeX0eXbvcAJhTnux+vgC890ySjsRq5ROgK+KQXAQAwPvvDnrEXTYwA9tGeo3jkugr8w59Pw2M3Giu9vXtNCIBxRUMjPrcxhtC3d7Xw/LG4csYojB6eZx+nDuC4JwAARm/u6yWFEABrP2rCHU9tAQDce9kUjCzIc913yshhuGNOKeZPHRHz9b/qCOHiKSPg18R4bV1h2ugCFObn9Kpn8faOFvzv3+3CrSv/iCX/WZPyUt4MAjSgrBRQPAqGRJ+3YDUjPjHmkz+yvs6VvxYYc8B7y/nHrOsKc8qK8dMbZ+L8kkJMNxezRToWa3xhQcUY5Po1+ATIzdGwZP5k5DgaSZ9PIm7XGclZQ3y4IEZg609KGektq722pn2e7tSxadch3HvZFCyeU4qKsYUYVZCHgiF++DSBT4zP4d5Lp0QMYFbjGYvfJ7jQ8dqWVz/ej407W9Di6aF4n00AfNzchk+b29AZUvZ0VWu/i5tmlyDXr0HM15pTVoybZpfgP++cYw8MR/PpvjZ06cbspJACdh5sx9ETnWH3iyfMhRTw1o4W3LrqjykNBJwiSgPKGuB++v2GmLOHAKDtVPgZHmD84Vrz9DVNjAVIjgniVoPcm6qnQHi6we8Tu4SENa3yhlljXesRrFTUuLOH2t17awql87I1HnK4/UzcqYj20yF8tCfoGoDsLz5NoJuNG2A0rNEqi769owX3PfcnnOgI4W3Pe7lqhlFa4+26gxF7H5rAXuzlMxfZOWOC9XhrptP9L3zsWg8Szxm3gjG7SHmue+6jJlSMLcTiOaVYe/dc/GzDTtQ0BrGmugnrtjVj2bXlKOyhmmrhkBwcOxne6HtdOWMURhTk2VNKC/L8eK5mb8THhnTgZxt24sWl3+j5zfUDBgEacIvnGPWA1lQ3YfV7u9HYehIKRiP8remjcLozhANtp6M24KMK8uwBwc6QQsNh9/2+XlKIZQvLARiVSztCKq6ZReeNKsDOg92Lxi6ddi6CJztcZQry8/yuaZs5fg03e2Z8eMsxOC9HmjliGVGQiyPt4YPbSoWPLYwrGopxhUMQaDrW67ELp4nF+bhoygiUjy10re9QCvj9rkN2g617ZkNFW5h3ypxzv/oP4YPKAuDRG2a6gmL52ELU7m8LqyMEGJ/bkvmTY35m0UQaj9EVXIv0nD3H0506lkWoKQQAs0oKcfD4aRz+qgNNPUwFtowoyMNjN84E0D3bL9bAfE1jEIHGYEomgTAdRCmzeE4pHr91FvJyjPRJjl/D0ksm4zd3zsEPLpoU9XFDc3yuy5NGnoVcnxiLpnyCZQvL7YZ37ZJ5+Ic/n4bJcaSGvOWUzy3Is9dCWMcncNcVuqWyd1P+bp5dglwzPSSOf7k+warvVGGpuU4gEmcufX/wFD7d14afXF9h58j7Yvw5+fipmdu3akFZukIKpUVDcfn0UfhalFSYV/GwXNz/wsdhZ+wisF8HMFZDr/2oCcvfqMPNs0vw0xtnRvwcF88pxWM3zsQ5w3LCbotFYJwM+DwtnLVIzzs2JTB6lt4QcMOssXj1hxfj6yVnR7xdE2OW0ogCdw+iYmyh/bNV4iMmhZRtIcueAKVU5YQiLLu23F71aTUEVmOxofYAyscMx9N/3GMvSrti+ih8ecSoDZTjEyy9ZDKWRtmsxvr58bd2uV430jTB80YV4OjJTvt1rDNTZ70kwJjqat3HWt3bG+ePPxs1jUHAPP5vO2av1B04HrHXIgJ862vnYuPOFvus/Eynjue3NoUN4vo0xL3ng3OMJtIWovWHT/SYtgO6Zw55ewgCY2zhH82VwUDvdzWbNroAx09HTg06GYPXAmWWWll0QSl2vl6LkOPT9EdZpPfNqeEzxzQBpo4qwJrqpqgpPF0hbJqqJt0nFIHGIN7ddSji96kcJxN5OanbOIpBIE0NhnUC8bDqD3V06di656jdVbeCwm/unAMAuLJ8tKvMs66MM+NHrusuCRHtc1q1eberwZ8+ugCP3jgT67Y144Mvjtirfd/74giWzi9DwdAc1+furZd00+ySiOmLeN6rt8BZl64w9mxjJo21SCxSKsPv03BuQV5YcTvnAjTAyEUvvWQyVm7ebQcML59mnKlaZbotlROKsOiC0j5VhRVBWKpkYnE+Hr91Vthn1NtdzVZt3h1XqQ5NEyy/rsJeCb6loTW8sJ8YvSnnScaCCqPw4Pv1R1y9PKsUxxPvfN7ja7tfwhhLilXQzqp6e9fFk8J+3wYag0AayqZdj7ylga3NWQB3dUarIV6xqd6egqmUimtHrB2e2TidurKf77u/rnaVfKg7cNwOPF7eon439bIXEKnAmVUSxPocrKmlZw/LsWedWGmnm2aXGLVpzEDhXbGcl2Ok0yonFOHJ71XZJSTed6zWPd8cL4n2+3TT7BK8WLO31+sANBFcXT7atcBsyfzwvReA+HY1c5b3eMdzFn7hxCJMGVWAun1trvUPoZDx+3DvZd1laJxFB437dPc8rLEpq4qqc8Daubp3QcWYsEWGsYR0Zae5YhW005VCwdAc1/GmAmsHpaFM3by8L7z1h7ybs1jVGaPdv6cutFF+2P2cZSOG2T97/8CjTWFNtKjfmuom/NaThwbc5ZKt9yUC17RDv0/swWcrNeVd8PT1CI175YQi3HfFedi656j9ecUKANZj1i6ZZ5d+BoBNuw6hM0ZQcO5nUFo8rMca+z3tn+0sdKeJhC0AnDKqAI/dODOsAmqkXd+szYheCjRHLT9t9cA0AS6PsHGQc0abtbLcogFhARmAK6BHiwPpsrdzSoOAUmo9gPVVVVV3p/I40o33Fyr1a7r7jzfnvutgu72RDBDeKEfa0yCWl7c1h50x33PJZPuyNy0QreFKpKifd5W0k1V4zPm+3qo76ErzzDD3MwC6U1OBxiDe/fxwj417bz8v52tYrLPy9lOdeOr9L+09FC6ddi5GFuS50mLW2XVPzx/rmJy9Q6WUa+quFRCt57ECVrT0nPVebp5dEvH1nPsyhBTw+88OuX4/LIvnlKJuf5trfERgrAVZdm05ave3ob6luyy4royU26M3zsRDr2yHbt5/8rlnYfehr9Jqb2emg9KQc/PyHMcv/WDlbHSs/2M1yt5GKhZvAL1ixqiwx8bTcDlr5ltF/eI9Bu9MFCfnpirWv6L8XFcgXHRB5M9g7d3xNe69+bx6erw1NpNoDjvWMTnHDBSAu79ZhuNnuqJOI+1LYLPMLSu21ysA4fsyOHl/l7y9rxWb6rF1T9DVwN972RTXuhHnZkS5cfRkBwKDQBqyznCyYWA4knga5Xh5A+rSCGd58ejLGbXFm3K6cGIRPm5ui5ieAOLvnSTauPfFQLxmpM1lHrhmer+8VuWEIiy/viKuhtn7u+TtfUVLc1mfmTWmZJUkWXZt7NTcQGEBORr00mGmlXc/2nQ4pnRl5emtxnQgdtWL9/vo6X6xSqsvWvWh3eOwSk8nOijMKqJENCgNtiC5YlM9fvG7XXZKya8Jnr8n8ZLSyQgCTAcRUdpJRaqrPyUyptTfGASIiPpZImNK/Y1BgIhoAKRr74YF5IiIshiDABFRFmMQICLKYgwCaSjQGMSKTfUp33uUiAY/DgynGWfxrNwBWihDRNmLPYE04y2tnKrdhogoOzAIpBnvdobpUGCKiAYvpoPSTDovKiGiwSelPQERWSgiq9va2nq+MxERJR03lUkzHBgmooHEMYE0w4FhIhpIDAJphgPDRDSQODCcZjgwTEQDiUEgDaVrtUEiGnyYDiIiymIMAkREWYxBgIgoizEIEBFlMQYBIqIsxiBARJTFGASIiLIYgwARURZjECAiymIMAmmIewwT0UBh2Yg0w1LSRDSQ2BNIMywlTUQDiUEgzbCUNBENJKaD0gxLSRPRQGJPIM0EGoMMAEQ0YNgTSCMcFCaigcaeQBrhoDARDTQGgTTCQWEiGmhMB6URDgoT0UBjEEgz3F+YiAZS0oOAiNwA4L8BGA7g10qpt5L9GkRElBxxjQmIyNMickhEaj3XXy0iu0SkXkQeAACl1KtKqbsBLAWwKPmHTEREyRLvwPAzAK52XiEiPgArACwAMAPA7SIyw3GXh8zbiYgoTcUVBJRS7wE46rn6QgD1SqkGpVQHgOcAXC+GnwPYoJTaFu05RWSJiNSISM3hw4f7evxERJSARKaIjgOw13G52bzuvwO4AsAtIrI02oOVUquVUlVKqaqRI0cmcBhERNRXSR8YVkr9EsAvk/28RESUfIn0BPYBGO+4XGJeR0REGSKRILAVwFQRmSQiuQBuA/B6b55ARBaKyOq2trYEDoOIiPoq3imiawF8CGCaiDSLyJ1KqS4APwTwOwA7AbyglKrrzYsrpdYrpZYUFhb29riJiCgJ4hoTUErdHuX6NwG8mdQjIiKiAcMCckREWYxBgIgoi6U0CHBgmIgotVIaBDgwTESUWkwHERFlMQYBIqIsxiCQZtZUN+G7v67GmuqmVB8KEWUB7iyWRtZUN+HBV7YDAP7wxREAwOI5pak8JCIa5Dg7KI1sqD0Q8zIRUbJxdlAaWVAxJuZlIqJkYzoojVipnw21B7CgYgxTQUTU7xgE0sziOaVs/IlowHB2EBFRFmMQICLKYpwdRESUxTg7iIgoizEd5BFoDGLFpnoEGoOpPhQion7H2UEOgcYg7nhqCzq6dOT6NTx711xUTigasNfe0tCKuWXFA/aaREQMAg5bGlrR0aVDV0Bnl44tDa0D0iCnMvgQUXZLjyCwaxdw6aXu6269FfjrvwZOngSuuSb8MX/xF8a/I0eAW24Jv/2v/gpYtAjYuxf47nfDb7//fmDhQuO177kHAPC901244MBx6Eph5cW3Yf+xUux48z3M+F/Lwh//2GPAN74B/PGPwIMPht/+xBPArFnAO+8Ajz5qX91+ugvHT3ci+K//BxXfuhBYvx6jH34Uzxw9ad9n9Pp8YP1LwPjxwPPPA7/6Vfjzv/QSMGIE8Mwzxj+vN98E8vOB//gP4IUXwm9/913j/1/8AnjjDfdtQ4cCGzYYP//kJ8DGje7bi4uBdeuMn3/0I+DDD923l5QA//Vfxs/33Qd8/LH79vPOA1avNn5esgT4/HP37bNmGZ8fAHznO0Bzs/v2efOAf/kX4+ebbwZaW923X3458PDDxs8LFgCnTrlvv/Za4O//3vjZ+3sHpOR3z+Whh4ArrjA+t/vuC7+9j797tlWrgGnTgPXrgccfD7/9N7/h7x6Qmb97fcAxAYeCIX5MHzMc5w4fAohg7UdNePCVT9F+uispz99+ugs7DxzH3qMn8Q8vfWKPOwwfkgNNBACgiWD4kJykvB4RUU9EKZXqY0BVVZWqqalJ9WEAMFIzT7zzOT6oPwJdAT4B/u6qabj3sikJP/eKTfV4/K1dEZ+XYwJE1FsiElBKVSXyHClNB4nIQgALp0xJvIFNBis3f6ZThwKgCZDj11CUn4sVm+oTbqDnlhUj16+hs0tHjl/D3LJi+7bKCUVs/IlowKU0CCil1gNYX1VVdXcqj8NiDQwrGHmyi6aMwIKKMVj+Rl1SBm0rJxTh2bvm8oyfiNIGxwQcivJzoYlAEyA3R8N9V5yH4MkOe8bQmU4dT7zzeUJrCConFOHey6ZEDQBcp0BEAyk9Zgf1g97m2AONQSx/ow4hXcGnCZZdW47KCUXYdbAdujlsogC8/8URbN1zNO4eQW+Og1NFiWigDcogEGgM4vbVH6IzpJDjE6xdMq/HxtSZCtKVQvBkBwAgeLIDAiMAwPw/3jUEkRp1AFi3rRkC4KbZJa7nSNU6BSLKXoMyCKzavBsdIaPZ7ggprNvW3GNjWpSfa5/x68q4DBiDuXk5Gjo6dejoHix2DupG423U121rxks1e+1jezHQjLV3d5/txxo4JiLqD4MuCAQag9i4s8V13ZH2Mz0+znnGr5mXAfdgblF+LoInO8JSO9FSPt5GXQB0hrqn5HrP9jlwTEQDbdAFgS0NrfYZveXdzw8j0BiM2ajOLStGjtlg++OYvmk1/EX5uVFnD3kbdQB4fmsTunTjOXw+wf5jp7CmuskVXNj4E9FAGVRBINAYxL5jp+zG3IoFoVD0/LqzMYe1cK6HBXTOXL8mAl0pV8rHeSbvbNQDjUFomgboOgSAriusqW6CAiAwUk0V4wqx6AJuMUlEA2PQLBZzNsx+TXDljFF49/PDCIWi59edjxEAVqamM6Rw/wsfY8n8yVg8pzQs3ePM9evKmE2kmYHjxZq9COkq4uyeLQ2t6AoZ3QAFd6xR5ut/0tyGT5q3AwADARH1u0GzWMzZMId0hREFebilsiRsFo6zQXc+xnVcAPa0nsSDr2xHU+sJPPPhHle6Z25ZMfya2AO8gIJoAl1XCJnXWfl+69jmlhXb6xCUUojd1wA21B5gECCifjdo0kHOQVifJngp0IyukNFw3zS7BIA5dfTJLfZA7SMLy5Hr1+wyEZH8tu5g2LTNey+bgm9XjbdTOboOAN0Nu6C73ISzdwIRhHQFTQBNEztgQBAWiBZUjEn6Z0RE5DVoVgxbg7B/d9U0fLtqPLpC3at8120zysG+vK3ZXgvQ0aWjbn8bnr1rLm6fUwqfJhGf9+ry0cj1a9AAiIg9dfSm2SXIy9HgEyDHJ8jxGz/n+gSL55Ti2bvmulYbd4aU/dpKATPGDMcVM0bB79dcAUATYOn8MvYCiGhADJqeANA9iyfQGMSL5nx8BSNPf/PskrCzfevyuLOH4u6LJ+Gp97+Erowz9fKx3QO0pcXDsOy1WuhKYfkbdZg2uiDizB/v1M5dB9shIoAn/aMD2L6vDXX7jb0LnARAwVCWkiaigTGogoBTafEw1B/6CoBxFv7ytmbcPLsEL9XstVcSV4wtdK3oXX59Rdg6gEBjEBtqD6DLPF0/09k908g7ndO7duCR12sR8uZ5TLrZJRBPB8Tn4yIxIho4gy4IWDN+TnfqrusVjEZ67ZJ5EQeGO7t0BE92uPYNsMpPdDgWeCkA7ac6oy4Qs67fvOuQ63FeAqNI3fypI/H2jhZ7mugtlSVcJ0BEA2bQBQGrYXfK9QluNgeHnSmjfcdOwe/Tok4jXbetOWJD/mFDa9iMIes5vUEjEr9P8K1p52JkQR7KxxbivS8O24PV1nESEQ2EQRcErFlCHZ06RIDLp4/CPZdMdqV31m1rtmcP+TXBbReWhhVzCzQGUbevLeJrjBo+BNv3tblmDAHA8vV1EQPAWXk+nOwIQVfG2f6fjT8bG3e2QFfGoPIj1xlpqKL8XPu52BsgooEwaGYHWSonFGHZteXQNIEC8N4Xh+3brFTRmuom15qCsWcPDQsAdzy1BZ82dwcBTYDzSwqxdH4ZRhTkwe8zZwY5poJ+0hw5aHxnzgR7hpEmQM2eIELKSC11hBRq97dhblkxlr9Rh8ff2oU7ntrC/QSIaEAMup4AYBR/i1TKYf+xUzjjGSuINBDrLCtt8WuCRReU2nWCnD2ILQ2tYc9ruXBiEa4sH23PMArp4QvF6lvaWUaaiFJiUAQB7yBttIVjfk2MRVrmjB0BcMl5I8NSMNbjnYvIQrrChtoDUXsQPk3sGUROW/cYvYr5U0dGDAAAEGg6htmlRdDEqGPKMtJENFAyPgh4awZ9u2o8bppdYs/h33/sFNZ+1ARdAV0hhbKRw/Bl60kopeDXBJs/P4yNO1vCNn2ZP3UkANj1h3yaYGiOL+JAcuWEItx18SSseq/BtWpYmf/OdOrYuLPFvs1n1JCzL+u6wlPvfxm2qxkRUX/L+AJyzjRKR8ioyrluWzOevWsu7r1sij0QbG0KU3/4BPw+wW0XlOKLlnZ8tMfIvZ/p1LFy827X1M5cs7RE7f42vBRoxjs7W6AJMLpwCIbm+rHrYDt2HWzH0+83YPfhE3ajrgkwbVQBdh5sBwBztzLzPQNYdEEpKsYW4mEzPQSB3UtQjl3NiIj6W0oHhpVS65VSSwoLC/v8HFbqxlpz5dz+EegeKD53eJ79mK6QcgUA63G//+yQa9OXji4dtfvbMO7soXYZii4d2HfsNOoPfYUHX9mOB1/ZjnpHAACMBv8zMwBYNE3gEyAvx5gGOm10AcR8lDIHiXuzaxkRUTJk/Owgq3zDFTNGGSWdPQ3pmuomLHutFgePu3cXqz98Iuy5lFnczemlQDOK8nONAnC94M39TyrOx20XltprCrY0tCLkGUu+aMoIbi5PRAMq44NAoDGIlZt34/efHYKuK2jSnVMPNAax7LXaiAO2Q3Pcb91awfuTG2bi/JLunkmn2Rv4dtX4uI/JJwgLJg1HTtiF7ABzJzNf951y/Rruu+I8BgAiGlAZPTAcaYWu7sipG2fbkVfvHmg7Db9mFIqbV1aMgqE59uyiaaML7OdVMHoDjywsxxBzw3lrwNfrhlljkZ/nxws1exEKKWgAys49Cw2Hvwqb+mmVsFi3rTlszwMiooGS0UFgS0OrK4cPAJqInQqyzrY7QgoC4PpZY9F6ogMf1B+xV+9eVT7aVS8IMFJMzv0CQiGjrpBzw/nlb9SF1SdqPdGB/Dw/usxj0gFMGjEMzcGTdlmInvYuJiIaSBkdBJyNPGDM1V9+fYWrYbWaaQXgzdqDeGRhObbuORqxUXa6aXYJXgw0G2sNzAVlzkZ72ugCrNy8G2/vaLEfs6BiDOr2u1cNn1uQ5yo5zUafiNJJRgeBSCkVAFixqd6uEhpy9BSsSqHeRjlaRdBYG89XTijC0ksmAwAOHT9t7z0QaAzawSPH3NWMZ/xElK4yOggA7pSKNUZg7RfwyHUVrp6Cphk7g3kf49xTwDl7p8ucux/SVVgZB+8itbr9bQg0Bo3AdDfP/IkoM2R8EAC6y0Z8sveY3eBbhdnWLpmHlZt3G1U7dYVHXq8FALtqp7MUhHPg1ll6IlLaKNYiNe+Zf9SeBhFRimV8EHCekXtn8guMnsLIgjx7xW5HSHWv1HXwri/wbh/pbbznlhXDr3X3MpyL1KL1GJw9DSKidJDxQcB5Ri4AzC194dME5WMLI+4L4A0AAmOhlneefqRcvnVWX5Sfa7yYOVk02mpfVgclonSW8UHAuYmMDtgT+O3Uj0jYTmOadNfyAYyAEc9CLedZvSZibxKvIXIQsY7Prwk6Q0ZxOJaEIKJ0kvErhq20zUVTR7hW6SoYG8x3egMAjN3G/JpAYOwTcNfFk7BuWzN+/Mr2mJu5OM/qu3QFEaMeUG5OD6t9rd3kvbvKExGlWMb3BAAjECyoGIMPd7dCKWUXY/NrAoigq8voJWhilGdYeslkXDbtXGyoPYDyMcPx9Adf2rn9FwPNWHt35Lx9UX6ua9xBxKgIGmu175aGVnSFdHvRGdNBRJROBkUQCDQGsfyNOrse/10XT7LLQACwc/jBkx32ddYOYR/udpeWiJa3t1/DkUZSEbam9OpplhFnDhFRKg2KIODcDlLXFeoOHHelZ7yN64pN9XZaB8qoHGo17tFWEVuvYZEY93WKNcuIM4eIKNUyflMZwD1dUwfwQf0RbN1zNGqj6j07X3atsXFMrEJu3i0rrR3M4mm0o60Y5swhIkq1lAYBpdR6AOurqqruTvjJHNM1e2pUI52dW2mZaHpaN9AXPaWKiIj626BJB3WFepeqiad0hDdfn+waQP0RWIiIeiPjg0CgMYh9x07ZG8D3NlUDRE7LABiQfD2LyxFRKmV0EPAWcbvtwtjTNaOJlJZhvp6IskFGB4GXtzXjjLnTVyiO6ZrRREvLMF9PRINdxgaBQGMQL9bstbd5tDZ+6StvWob5eiLKBhkbBKx6/4AxEHxLZfL36GW+nogGu4ytHVSUn2sXgVMAKsYWpvR4iIgyUcYGgeDJDruOj2ZeJiKi3snYIDC3rBh5OZpdxZMDt0REvZexYwLWwK21yTwREfVexvYELC9va8baj5pwx1NbYu4FQERE4TI6CERb6UtERPHJ6CBgrfT1Rdnfl4iIYsvYMQGAC7qIiBKV0UEA4IIuIqJEZHQ6iIiIEsMgQESUxRgEiIiyWNYFgUBjECs21XNNARERBsHAcG9E20aSiChbZVVPgIvLiIjcsioIcHEZEZFbVqWDuLiMiMgt6UFARMoA/BhAoVLqlmQ/f6K4uIyIqFtc6SAReVpEDolIref6q0Vkl4jUi8gDAKCUalBK3dkfB0tERMkV75jAMwCudl4hIj4AKwAsADADwO0iMiOpR0dERP0qriCglHoPwFHP1RcCqDfP/DsAPAfg+iQfHxER9aNEZgeNA7DXcbkZwDgRKRaRlQD+TER+FO3BIrJERGpEpObw4cMJHAYREfVV0geGlVKtAJbGcb/VAFYDQFVVlUr2cRARUc8S6QnsAzDecbnEvI6IiDJEIj2BrQCmisgkGI3/bQAW9+WJAoHAERFp7ONxjABwpI+PHQyy+f1n83sHsvv9Z/N7B7rf/4REnyiuICAiawFcCmCEiDQD+Cel1K9F5IcAfgfAB+BppVRdXw5CKTWyL48zj61GKVXV18dnumx+/9n83oHsfv/Z/N6B5L7/uIKAUur2KNe/CeDNZBwIERENvKyqHURERG6DIQisTvUBpFg2v/9sfu9Adr//bH7vQBLfvyjF2ZlERNlqMPQEiIiojxgEiIiyWEYHgUhVTAcTERkvIptEZIeI1InI35rXnyMib4vIF+b/Reb1IiK/ND+PT0VkdmrfQeJExCcifxKRN8zLk0Sk2nyPz4tIrnl9nnm53rx9YkoPPAlE5GwReUlEPhORnSIyL1u+exH5H+bvfK2IrBWRIYP5u49Uqbkv37WIfN+8/xci8v14Xjtjg0CWVDHtAnC/UmoGgLkA7jXf4wMANiqlpgLYaF4GjM9iqvlvCYBfDfwhJ93fAtjpuPxzAP+mlJoCIAjAKlt+J4Cgef2/mffLdP8O4LdKqa8BOB/G5zDov3sRGQfgbwBUKaUqYKxDug2D+7t/Bp5Kzejldy0i5wD4JwBzYBT4/CcrcMSklMrIfwDmAfid4/KPAPwo1cfVz+/5NQBXAtgFYIx53RgAu8yfVwG43XF/+36Z+A9GKZKNAL4F4A0AAmOVpN/7OwBj0eI882e/eT9J9XtI4L0XAvjS+x6y4btHd3HKc8zv8g0Afz7Yv3sAEwHU9vW7BnA7gFWO6133i/YvY3sCiFLFNEXH0u/MLu6fAagGMEopdcC86SCAUebPg+0zeQLA/wSgm5eLARxTSnWZl53vz37v5u1t5v0z1SQAhwH8XzMd9pSIDEMWfPdKqX0AfgGgCcABGN9lANnz3Vt6+1336Xcgk4NA1hCRswCsA3CfUuq48zZlhPxBN89XRK4FcEgpFUj1saSIH8BsAL9SSv0ZgBPoTgcAGNTffRGMvUkmARgLYBjCUyVZpT+/60wOAllRxVREcmAEgGeVUi+bV7eIyBjz9jEADpnXD6bP5CIA14nIHhgbFn0LRo78bBGxyp0435/93s3bCwG0DuQBJ1kzgGalVLV5+SUYQSEbvvsrAHyplDqslOoE8DKM34ds+e4tvf2u+/Q7kMlBwK5ias4SuA3A6yk+pqQSEQHwawA7lVL/6rjpdQDWyP/3YYwVWNd/z5w9MBdAm6M7mVGUUj9SSpUopSbC+G5/r5S6A8AmALeYd/O+d+szucW8f8aeJSulDgLYKyLTzKsuB7ADWfDdw0gDzRWRfPNvwHrvWfHdO/T2u/4dgKtEpMjsTV1lXhdbqgdDEhxIuQbA5wB2A/hxqo+nH97fxTC6gJ8C+Nj8dw2MfOdGAF8AeAfAOeb9BcaMqd0AtsOYXZHy95GEz+FSAG+YP5cB+AhAPYAXAeSZ1w8xL9ebt5el+riT8L5nAagxv/9XARRly3cP4J8BfAagFsBvAOQN5u8ewFoY4x+dMHqBd/bluwbwA/NzqAfwl/G8NstGEBFlsUxOBxERUYIYBIiIshiDABFRFmMQICLKYgwCRERZjEGAiCiLMQgQEWWx/w/58ADcvZOMpwAAAABJRU5ErkJggg==\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": 27,
   "id": "d8baaea2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T2: 0.8291666666666667\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
}