Examples/NP_EAP_FeatureAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==================================================================================\n",
    "# Function: Read EAP matrix, extract EAP features and run EAP feature clusters\n",
    "# The code reproduces Fig.2d\n",
    "# @Author: Yina Wei, 2021\n",
    "#==================================================================================\n",
    "import pandas as pd\n",
    "from pandas import Series, DataFrame\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.cluster import KMeans\n",
    "from time import time\n",
    "from scipy import stats\n",
    "import ecp_func as ef # import self-written sub functions  \n",
    "\n",
    "plt.rcParams.update({'font.size': 9})\n",
    "sns.set_style(\"whitegrid\", {'axes.grid' : False})\n",
    "my_pal_FSRS = {0:[5/255.0,23/255.0,182/255.0],1:[23/255.0,163/255.0,95/255.0],2:[213/255.0,9/255.0,21/255.0],3:[254/255.0,200/255.0,46/255.0]}\n",
    "fs_pal=['steelblue','seagreen','mediumpurple']\n",
    "rs_pal=['crimson','darkorange','saddlebrown']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The set of parameters \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load  EAP waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1313, 82, 40)\n",
      "[75. 16. 27. 42. 36. 73. 53. 62. 47. 40. 32. 62. 36. 84. 58. 67. 41. 52.\n",
      " 68. 49. 96. 58. 54. 53. 32.]\n",
      "25\n"
     ]
    }
   ],
   "source": [
    "# File paths\n",
    "path = 'data/wt_EAP_data_info/'\n",
    "savedir='data/wt_EAP_data_info/'\n",
    "savefigdir='Figures/'\n",
    "\n",
    "# session_id \n",
    "session_id=[732592105,739448407,744228101,750332458,750749662,754312389,754829445,757216464,757970808,759883607,761418226,\n",
    " 763673393,766640955,767871931,768515987,771160300,771990200,774875821,778240327,778998620,779839471,781842082,793224716,799864342,847657808]\n",
    "\n",
    "unit_num = np.zeros(len(session_id))\n",
    "for i in range(len(session_id)):\n",
    "    with open(path+str(session_id[i])+'.npy', 'rb') as np_file: # EAP waveform\n",
    "        data=np.load(np_file)\n",
    "        unit_num[i]=data.shape[0]\n",
    "        if i==0:\n",
    "            waveform0=data\n",
    "        else:\n",
    "            waveform0=np.concatenate((waveform0,data))\n",
    "print(waveform0.shape)\n",
    "print(unit_num)\n",
    "[N,T,channels]=waveform0.shape\n",
    "print(len(session_id))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1287\n",
      "1287\n"
     ]
    }
   ],
   "source": [
    "fr0 = []\n",
    "depth0 = []\n",
    "unit_id0 =[]\n",
    "good0=[]\n",
    "for i in range(len(session_id)):\n",
    "    data = pd.read_csv(path+str(session_id[i])+'unitinfo.csv',index_col=0)\n",
    "   # print len(data['fr'])\n",
    "    fr0.append(data['fr_5min'])\n",
    "    depth0.append(data['depth'])\n",
    "    unit_id0.append(data['units'])\n",
    "    good0.append(data['good'])\n",
    "fr0 = np.array(np.concatenate(fr0))\n",
    "depth0 = np.array(np.concatenate(depth0))\n",
    "good0= np.array(np.concatenate(good0))\n",
    "unit_id0= np.array(np.concatenate(unit_id0))\n",
    "\n",
    "bad_unit = np.array([915960371, 915961088, 915960802, 915960838, 951868026, 951868229,951868527, 951869594, 951816073, 951002808, 951814989, 951815653,\n",
    "                     951815724,951815709, 951816276, 951816492, 951877391, 951796229,951768170,951013844,951026044,951027930,951896781,951896911,951800909,950938909,\n",
    "                     950935800,950929267,950927995,950935732,950929283,950930487,950931566,951801768,951762773,951095926,951768557,950928046])\n",
    "\n",
    "for i in range(len(unit_id0)):\n",
    "    if unit_id0[i] in bad_unit:\n",
    "        good0[i]=0\n",
    "\n",
    "good = np.where(good0==1.0)[0]\n",
    "fr0=fr0[good]\n",
    "depth0=depth0[good]\n",
    "unit_id0=unit_id0[good]\n",
    "waveform=waveform0[good,:,:]\n",
    "[N,T,channels]=waveform.shape\n",
    "print(len(fr0))\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "[-100  -80  -60  -40  -20    0   20   40   60   80  100]\n"
     ]
    }
   ],
   "source": [
    "eleRowN = int(channels/2)\n",
    "print(eleRowN)\n",
    "sampling_rate=30000 #sampleing rate is 30kHz\n",
    "fs=sampling_rate\n",
    "dt=1.0/sampling_rate\n",
    "t_STA=np.array(range(0,T))*dt*1000-0.6 \n",
    "Th=5\n",
    "dist=np.array(range(-100,101,20))\n",
    "print(dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate EAP features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "features0=np.zeros((len(dist),len(good),4))     # dist*num_cells*num_features  #dist only includes -100,-80,-60,-40,-20,0,20,40,60,80,100\n",
    "all_features0 = np.zeros((len(good),11))        # each cell has 11 features\n",
    "max_waveform0 = np.zeros((len(good),T))         # save the maximal EAP waveform\n",
    "norm_oneColEAP = {}\n",
    "badch=[]\n",
    "saveEAP = 0 # 0: dont save EAP, 1: save EAP \n",
    "for num in range(len(unit_id0)):\n",
    "\n",
    "    EAP=waveform[num,:,:]  # amplify =0.195, the EAP of the index[k] cell \n",
    "    dist0 = np.zeros((eleRowN,1))\n",
    "    \n",
    "    #Select the one column\n",
    "    oneCol_EAP = np.zeros((eleRowN,T))  # channels * times -- one column\n",
    "    for i in range(0,eleRowN):\n",
    "        temp0 = ef.butter_highpass_filter(np.transpose(EAP[:,i*2]),300,fs)  # even index\n",
    "        oneCol_EAP[i,:]= temp0\n",
    "        dist0[i]=i*20    \n",
    "    \n",
    "    # dist = 0, the channel cloest to soma, which have the largest negative peak\n",
    "    maxCh = 10\n",
    "    maxCh0=np.unravel_index(oneCol_EAP.argmin(),oneCol_EAP.shape)[0]\n",
    "    if maxCh0!=10:\n",
    "        badch.append(num)\n",
    "    dist0 = np.squeeze(dist0)- dist0[maxCh]\n",
    "    indx = np.where(np.logical_and(dist0>=-100,dist0<=100))[0]\n",
    "    \n",
    "    if saveEAP==1:\n",
    "        f=1/abs(min(oneCol_EAP[maxCh,:]))*3\n",
    "\n",
    "        fig, ax = plt.subplots(figsize=(2.5*0.46,3.8*0.5))\n",
    "        ax=plt.subplot(111)\n",
    "        for i in range(maxCh-7,maxCh+8): \n",
    "            if i==maxCh:\n",
    "                plt.plot(t_STA,oneCol_EAP[i,:]*f+i,color='grey',linewidth='2') \n",
    "            elif abs(min(oneCol_EAP[i,:]))>=abs(min(oneCol_EAP[maxCh,:]))*0.12:\n",
    "                plt.plot(t_STA,oneCol_EAP[i,:]*f+i,color='grey',linewidth='1') \n",
    "            else:\n",
    "                plt.plot(t_STA,oneCol_EAP[i,:]*f+i,color='grey',linewidth='1') \n",
    "        plt.yticks([])\n",
    "        plt.xticks([])\n",
    "        simpleaxis_all(ax)\n",
    "        plt.title(str(unit_id0[num]),fontsize=9)\n",
    "        fig.savefig(savedir+'EAP/'+str(unit_id0[num])+'.eps',bbox_inches = 'tight')\n",
    "        fig.savefig(savedir+'EAP/'+str(unit_id0[num])+'.pdf',bbox_inches = 'tight')\n",
    "        plt.close()\n",
    "   \n",
    "    #Extract EAP features\n",
    "    #Extract featrues based on EAP data as the distance changes\n",
    "    A,W,TPW,t_NegPeak=ef.cal_EAP_features(t_STA,oneCol_EAP,0.01)       #oneCol_VeSTA,channels*times    \n",
    "    t_NegPeak = t_NegPeak-t_NegPeak[maxCh]  # set the maximal channel as 0\n",
    "    \n",
    "    norm_oneColEAP[unit_id0[num]]=oneCol_EAP/max(A)\n",
    "    \n",
    "    # Preprocess the data of A,W,TPW,t_NegPeak, deal with unphysiological values\n",
    "    dth = 100  \n",
    "    idx_below = np.where(np.logical_and(np.logical_and(dist0<=0,dist0>=-dth),A>=0.0001))[0]\n",
    "    idx_above = np.where(np.logical_and(np.logical_and(dist0>=0,dist0<=dth),A>=0.0001))[0]\n",
    "        \n",
    "    W[idx_below]=ef.outlier_below(dist0[idx_below],W[idx_below],0.3)\n",
    "    W[idx_above]=ef.outlier_above(dist0[idx_above],W[idx_above],0.3)\n",
    "    TPW[idx_below]=ef.outlier_below(dist0[idx_below],TPW[idx_below],0.3)\n",
    "    TPW[idx_above]=ef.outlier_above(dist0[idx_above],TPW[idx_above],0.3)\n",
    "    t_NegPeak[idx_below]=ef.outlier_below(dist0[idx_below],t_NegPeak[idx_below],0.3)\n",
    "    t_NegPeak[idx_above]=ef.outlier_above(dist0[idx_above],t_NegPeak[idx_above],0.3)     \n",
    "    \n",
    "    features0[:,num,0]=A[indx]\n",
    "    features0[:,num,1]=W[indx]\n",
    "    features0[:,num,2]=TPW[indx]\n",
    "    features0[:,num,3]=t_NegPeak[indx]\n",
    "    \n",
    "    #one-channel-features:'A','W','TPW','PTratio','REP'\n",
    "    max_waveform0[num,:] = oneCol_EAP[maxCh,:]/abs(min(oneCol_EAP[maxCh,:]))\n",
    "    max_waveform0[num,:] = max_waveform0[num,:]-max_waveform0[num,0]\n",
    "    one_channel_feature  = ef.cal_one_ch_features(t_STA,oneCol_EAP[maxCh,:])\n",
    "\n",
    "    # multi-channel-features:'spread_below','spread_above','W_below','W_above','PL_below','PL_above'\n",
    "    dth = 100\n",
    "    spread_below,spread_above,W_below,W_above,PL_below,PL_above,intercept0,intercept1,intercept4,intercept5=ef.cal_multi_ch_features(dist0,A,W,t_NegPeak,maxCh,dth,Th,0.005,3)\n",
    "    all_features0[num,:]=one_channel_feature+[spread_below,spread_above,W_below,W_above,PL_below,PL_above] # one channel features + multiple channel features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1204\n"
     ]
    }
   ],
   "source": [
    "allindex = np.arange(len(all_features0))\n",
    "\n",
    "neg_time=[]\n",
    "for i in range(len(max_waveform0)):\n",
    "    neg_time.append(ef.find_neg_peak_time(max_waveform0[i,:],t_STA))\n",
    "    \n",
    "outlier=np.unique(np.concatenate([ef.outlier_index(all_features0[:,1]),ef.outlier_index(all_features0[:,2]),ef.outlier_index(all_features0[:,4]),ef.outlier_index(all_features0[:,7]),ef.outlier_index(all_features0[:,8]),ef.outlier_index(all_features0[:,9]),ef.outlier_index(all_features0[:,10])]))\n",
    "remove_indices = np.unique(np.concatenate([badch,outlier,np.where(abs(np.array(neg_time))>0.1)[0],np.where(all_features0[:,0]<15)[0],np.where(all_features0[:,0]>300)[0]]))\n",
    "\n",
    "allindex = np.arange(len(all_features0))\n",
    "index0 = np.setxor1d(np.indices(allindex.shape), remove_indices)\n",
    "\n",
    "all_features = all_features0[index0,:]\n",
    "features = features0[:,index0,:]\n",
    "max_waveform = max_waveform0[index0,:]\n",
    "\n",
    "FR = fr0[index0]\n",
    "depth_info = depth0[index0]\n",
    "unit_id = unit_id0[index0]\n",
    "\n",
    "print(len(index0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df0 = pd.DataFrame(all_features,columns=['A','W','TPW','PTratio','REP','spread_below','spread_above','W_below','W_above','PL_below','PL_above'])\n",
    "df0.insert(0, column='Spread',value=df0['spread_above'].values+df0['spread_below'].values)\n",
    "df0.insert(0, column='Diff_spread',value=df0['spread_above'].values-df0['spread_below'].values)\n",
    "df0.insert(0, column='FR',value=FR)\n",
    "PL_dist=ef.sym_dist_fc(df0['PL_below'].values,df0['PL_above'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering between FS and RS based on one-channel features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_r = df0[['TPW','REP']].values   # one-channel features\n",
    "X_r2 = ef.pre_process(X_r)\n",
    "\n",
    "numk=9\n",
    "[K2,wcss2,fk2] = ef.selectionK_kmeans(X_r2,numk) # calculate elbow method and density function\n",
    "elbow_idx2  = ef.findElbow(wcss2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAABnCAYAAACacLX/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPpElEQVR4nO2deVBT1xfHvyEaK1hrmQatHbDKJogrEtDaMjJa6YgjIGUx0lqnWq1xV5yWYVyKKx1RwI46hVaxIqXFCm4joaMIRBCLKEWQUkUHUFaJgmAI5/cHJRrgJWmSF2J/+czkj/feefecfHNzc+99595wiIhgQu+Y9XcA/1VMwrKESViWMAnLEiZhWcIkLEv0q7C7du2CWCzuzxBYo1+ElclkWL9+PTIyMvrDvUHoN2EDAgLg5+en1ra8vFzpePv27WhubmYrNL2hk7C3b98G0CVUcnIyUlNToclAztzcHNOnT2e8npycDH9/f/j7+6O2traXzx9//FGXsA2C1sLGxcUhJiYGcrkc27ZtQ2FhIR4+fIidO3fqHFRQUBBSU1ORmpoKS0tLpWsikQgHDx5EZ2fni5PV1Tr71DcDtL2xsLAQ8fHxeP78OcRiMbKyssDj8bBw4UJ9xteL6dOnw9zcHGKxGB9++GHXyXfeAYxsykPnNjY/Px8TJkwAj8cDADx79kznoFTB4XAgEokQFxfHqh9d4Wg7u7V7925wOBxIJBJs2rQJzs7O+P777/H48WPs2LFDbwHevn0bTk5OSudaW1thY2ODgoICvPvuuwCHY3Q1FqQlcrmcLl26RCUlJUREVFZWRvHx8dTe3q5tkX3SXX5P1q9fT2FhYV0HS5fq1ac+0LopMDMzg5WVFZycnCCTyVBYWIhhw4Zh4MCB+vzcGVmxYgUSEhK6mp4jRwzi899glL0CTbCzs4ObmxuSk5MBV1eD+PxXaFvVlyxZQkRE7e3t5O7urmgCQkJC9PNd+gempoCI6MyZM+Tq6kqk/dtgjVeuV/Ay3t7eaGpqMpi/f4PWwtrb22PPnj349ttv8emnn6KpqQlRUVFwdnbWZ3wq4XK5WLFiBZoGDzaYT43Rtqr3d6+gm4aGBuLz+RQZGUnPnz/Xq29d0KlXIJfL4ejoCADo6OiAjY2NokkwFJaWligXCpGVlQUPDw/cunXLoP6Z0KlXkJ6ejvb2dgDAkCFDcPr0aRw6dEhvwWnKG/v348KFC/jyyy/h5eWFyMhIyGQyg8ehhLZVPTAwsNe5jo4OWrBggU5foZ6oawqISKlXcP/+fZozZw45OTlRREQEXblypV+aCK1rbF9feS6Xi0GDBun0QeuKtbU1zp8/j++++w4ymQxr1qwBn8/H/PnzERsbi4KCAoPUZq3nCnbs2AE7OzsEBQUpzqWkpKCwsFCvg4S+5gp6cf26ykFCbW0tMjMz8fvvv0MikeDevXuYMmUKpk2bhsmTJ8PGxgbW1tZ4++23MWCA1hN+SmgtbHZ2Nq5cuYLs7GyYmZmBiODl5QWRSKTXHzB9CNuT5uZm5OfnQyKR4NatW3jw4AEePHiAuro6DB8+HCNHjgSfz4eVlRX4fL7S66233lK8hg4dyuhDa2G3b9+OkpISvPbaa3B1dYVAIMDkyZPViiqTybBx40bU1dXBzs4O27ZtA4fDYbTXSFg9zW7JZDJUVVWhpqYGdXV1qKurQ21tLerq6lBfX694dV97+vQpc0jaCtvNs2fPUFxcjBs3biArKwtmZmY4evQoo/3Zs2dx9+5diEQiREZGwsvLS+VjGkMKq090alAqKiqQm5uL/Px8tLa2wtHRER4eHirvKSoqwuzZswF0PQ0oKChQKeyritbCvv/++5g0aRJ8fHywfft2vPnmmxrd9/TpU1hYWAAABg8ejNbW1l42ycnJXbNWADZs2KC+0C1bNA/cQGgt7Jo1a5Cfn4+kpCTk5eXBw8MDAoEAw4YNU3nfkCFDFGK2trZiyJAhvWyCgoIUvY28vDzF02AmGmfPhqUaGwBobGzs9XBSF7sBAwbA3t6+74v66AxXVFTQ8ePHKTAwkHx9fVXapqWlUVxcHBERffPNN5SVlaXS3s/PT61/TWzYsFOF1jW2sbERRUVFKCoqQnFxMYgIU6dOhUAgUHmft7c3Nm/ejODgYIwePRozZszQNgSjRmthV65cCYFAADc3Nyxbtgzm5uYa3Tdw4EDs27dPW7evDFoLm5SUpM84GHl5ZKeLDRt2qtC5H2uib0z5sSxhEpYlTMKyhElYljAJyxImYVnC6IXtmSrfJyrmc/sLoxe2o6Ojv0PQCqMXViMOH+7vCHrx3xB22TK1Jm1tbbh586ZGi0/0gX4eSTKwa9cuuLm5YdasWYpz2dnZ2L9/P7hcLsLDwzFhwgTdHal5NFNSUoKQkBBUVVXB1tYWYWFh8PX1BZfLVbIrKyvDr7/+irt37yo9NOTz+XBxcYGNjY3GIbFSY1UtkDtw4AB++OEHxMbGajXLVV5ernGtIyIcOXIEnp6eWLVqFR49eoTNmzcjKioKY8eOxeHDh1FYWIitW7di/PjxmDlzJmpqajB16lS8/vrrqKqqQmZmJmJiYjB16lSMGzcOGzduRGZmpiIDSJVzjWhra9N4krelpYVycnIoJiaGMjIyFOelUimFhoYqjhcsWEAymazX/SdPniQ/Pz/y8/Oj7OxspWvvvfce/fbbb8o39PE2GhoayN/fnyZOnNgrm6azs5OysrJo7ty5NGbMGFq3bh3l5OSQXC5nfE9yuZzy8/Np69at5O7uTkOHDlWpgUph09PT6fPPPycfHx8KDg6mgIAAWrp0KaWnp1NnZ6fKgomol7DV1dW0fPlyxfHChQvpyZMnKsvoKUpaWhq5uLhQR0fHi5M+Pko2VVVVZGNjQ6tXr6Znz56pjVMbGhoaVF5nFHbLli109uxZamlpUTrf2tpKZ86coYiICLXOewrb3NxMn3zyieJ4wYIFygL1QV+1zcPDg3766SfGewICAujrr79WGx+bMArL9ElXV1drXHhPYYm63rRUKqXa2lpatGiR2jL6SorLzMwkOzu7F8luL9XYM2fOkK2tLbW2tmocJxswCisUCnt9TVNSUmjmzJkaF94t7L1792jPnj1ERJSTk0MBAQHk7+9Pf/zxh9oymLINvby86MiRI10H/7SxT58+pVGjRtHFixc1jpEtGIU9d+4cBQUFUVNTE9XW1tKyZctIJBLRo0ePDBkfo7ASiYSsra27vln/CLtp0yYSCoWGDI8RlT9eOTk55OvrS3PnzqULFy4YKiYlVOXHzps3j6Kjo4kAunHjBvH5fHr48KHhglOB2u5WUVERBQUFUU1NjSHi6YUqYW/cuEHDhw+n5uZmcnd3f9E0GAGMDxN9fX3B4XBARGhra4OZmRl4PB44HA5OnTr1rzv22qIuKS4kJAQuubk4b22tSMozCpgUV9VZJiKN+rH6QF2qfFlZGRFAt27dMkg8msL48UZGRiIlJaXXArXGxkYkJSVhi5Ekojk4OAAAXFxc+jkSZRgnYXg8HqysrDB9+nTY2dlhwIAB4HK5eOedd+Dv74+QkBBDxvnKwShsQUEBOjo6YGVlhWnTpinOczgceHp6GiQ4jUlL6+8IesEobHx8PMrLy1FcXIzZs2cbbB5TK4xw9TcrKUaq1hlER0fjypUrsLCwgIODAyIiIlSW9X+ZKs/ExYsXYW9vjwMHDiAyMhISiUSRDn/nzh3Ex8drnAH+qsJKp6+oqAju7u4AXqwz6KayshLh4eEIDQ3FzZs32XBvFLBSY5nWGRAR5s2bh88++wz19fVYtWpVn4ONnmsQ1KbKZ2f/N1Ple7Jjxw66du0aERGJxWKKjY0loq5BxcszZgEBAWqX4b+qqfKsNAXjx49HXl4eAEAikWDixIkAulYEhoSEQC6Xo6GhAURk8GX4hoIVYb29vVFRUYHg4GC0tLTA2toae/fuxbBhwxQrYkQiEcLDw9lwbxzoXOdZ5uTJk3qxYcNOFaZUeZYwkjm2/x4mYVmC1RQjXdBk+X1bWxs2bNgAqVQKHo+Hffv24Y033mAs88KFC0hPT8fBgwcZbXbu3Ini4mKYmZkhOjoafD5fq/iNtsZ2D4tPnDgBHo8HiUTSy+b06dOYMmUKEhMT4e3trRhU9EVNTQ2Sk5NVTiZdvnwZRIQTJ05g5cqVqKys1Dp+oxVW1bC4Gx8fH8VGwB0dHYzbjnR2dmL37t1Yt26dSp/5+fmwsLDAkiVLkJ6erlPCntEKq8nyewsLCwwePBgVFRVISkpCQEBAn2UdOnQI/v7+aoepjx8/xqNHj5CQkICRI0ciJSVF6/iNto3VZPk90DWt+NVXXyE6OppxjxaxWAyJRIL29nZUVlbi+PHjWLRoUS+7oUOHKrYQnDFjhk4PTY22xjINi1+muroaYWFhiImJga2tLWNZqampSExMxL59++Dq6tqnqAAwceJEXL16FUBXUzR69Git4zdaYXsOi/tafh8fH4+WlhbFNKSqvWg0oXtLlaCgIOTm5iIwMFDrskwjL5Yw2hr7qmMSliVMwrKESViWMAnLEkY7QOjm8uXLaGpqwogRI1BRUQGhUKi49vz5c+zduxelpaVob2+Hp6cnRCIRUlNTIZVKsXjxYo18lJaWorOzU6/7ixt9jS0pKYGzszNKS0sxbtw4pWuxsbGwtbXF8ePH8fPPP6O8vFyrf7TLyMhAtZ7/Ycmo+7EikQjXr1+HnZ0d/vrrL7i5uSEmJkZx/aOPPsK5c+cU04ktLS0YNGgQ0tLSIJVK4eTkBLFYjPDwcOTl5UEsFuOLL77A2rVrQUQYMWIEtm7divnz58Pc3BzHjh3DiRMnkJOTAy6Xi4iICDg6OmLOnDmwsrJCYGAgTp8+jba2NlhYWCA6OppxWyyjrrFxcXGYNm0aEhMTIRAIlEQFunZYfnmO1sLCQu3Gujdv3oSDgwOOHTuGWbNmgcPhwM/PD2vXrkVtbS3u3LmDpKQkREVFISoqCgBQX1+P+Ph4ODg4gMfjISEhAYsWLcKTJ08Y/Ri1sMuXL0dBQQFCQkJw/fp1rF69Wul6Z2en0h+mVVZW4s8//+yzrO4vpqenJ4YPH44lS5YgNzdXKQP877//RklJCUJDQxEWFgapVAoAGDlyJHg8HhwdHfHBBx9g+fLlOHXqVK+1uC9j1MKGhoZi9erV2LBhA1asWNGrxnp5eeGXX34BAMjlckRFRSm1lYMGDUJ9fT2AFxtKXLt2Dba2tjh69CjMzc1x9epVxZIAGxsbCAQCJCYmIi4uDt7e3gCgEL+srAxcLhcJCQmYNGkSzp8/zxi7UQtbWlqKsWPH4s6dO4rM7ZdZtWoVCgsLsXDhQnz88cdwdnZWTKQAXVneUqkUQqFQIay9vT0SEhIgFApx//59CAQCjB07FjExMbC0tASfz4dQKMTixYsxZswYJX+jRo1CRkYGgoODcenSJSVfPTHqH69XGaOusa8yJmFZwiQsS5iEZQmTsCxhEpYlTMKyxP8AJPpdQhCc+U4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 46.08x69.12 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(0.64,1.2*0.8))\n",
    "ax=plt.subplot(211)\n",
    "#Plotting the Elbow Curve by k(number of clusters)\n",
    "plt.plot(K2,wcss2/max(wcss2),'k-',markerfacecolor='None',linewidth=1)\n",
    "plt.plot([K2[elbow_idx2],K2[elbow_idx2]],[0,1],'r--',linewidth=1)\n",
    "plt.ylabel('WCSS')\n",
    "#plt.title('Elbow Method',fontsize=9)\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xticks(range(0,8,2))\n",
    "ax=plt.subplot(212)\n",
    "#Plotting the Density Function fk by k(number of clusters)\n",
    "plt.plot(K2,fk2,'k-',markerfacecolor=\"None\",linewidth=1)\n",
    "plt.plot([K2[fk2.argmin()],K2[fk2.argmin()]],[0,1.5],'r--',linewidth=1)\n",
    "plt.xlabel('# Clusters')\n",
    "plt.ylabel('f(K)')\n",
    "#plt.title('Density Function')\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticks(range(0,8,2))\n",
    "plt.ylim(0.5,1.2)\n",
    "\n",
    "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.5,wspace=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 1 ... 1 0 0]\n",
      "281\n",
      "923\n"
     ]
    }
   ],
   "source": [
    "# K-means\n",
    "kmeans2 = KMeans(init='k-means++',n_clusters=2,n_init=1000,random_state=11400)\n",
    "kmeans2.fit(X_r2)\n",
    "label2 =  kmeans2.labels_\n",
    "print(label2)\n",
    "\n",
    "print(len(np.where(label2==0)[0]))\n",
    "print(len(np.where(label2==1)[0]))\n",
    "\n",
    "if 'cell_type' in df0.keys():\n",
    "    df0=df0.drop(columns=['cell_type'])\n",
    "df0.insert(0, column=\"cell_type\",value = label2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGcAAABPCAYAAAD2vtuiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAM1ElEQVR4nO2deXxU1aHHvzf3zp6Zyb6SkEASICyBCC/VYl2qtBb66lbUvoALtLRa7dOq7WuV5+urPp/4aqlPbV3a+tQ+qNZIwbVIUD9hkVURAoGQBULCZCYLs2aWe/vHDYnzSdBoyc3Nx/l+Pnw+mXOX+d1855x7zrmHiaAoikICXZI01gESnJmEHB2TkKNjEnJ0TEKOjknI0TEJOTpGt3KCwRi9pyJjHWNM0aWc/3nsKDnlGymYuYlnnm8d6zhjhqC3GYLDR/3MnP8usjxYtubp2Vy5KHfsQo0Ruqs5D//vUWQZyiZbWXFDIQA33fohXd3hMU6mPbqSI8sKL2/oAGD5kgJ+fX85lRUOAkGZ2366f4zTaY+u5DQ0+uk9FcVmFVlWXYgoCjyxagaCAC+t7+DQYd9YR9QUXcl5f3cPADnZJuzJEgBzZjm55opcZBluufujMUynPbqSs3NPLwB5Oaa48vt+UoYownvbutm7r5fgzj0cW7yUQ0XTaSidzYkf3E648ehYRB5VdCVnR7+c+VWpceWTJlq54boCFAX+v/oRmi5eiHfDG8RcnUTbTtDz7AscqZxP9x+eH4vYo4Zu5CiKQkOjek/5xqVZQ7bf8+MSvibv5PrDj0M0imnWDPKe/R2Ff/0zlvnnQiRK+y134PrFg1pHHz0UneD29CmG7NcU24TXFb8/MmR7X0urstdZqOy3ZCr3Zt2oRL3euO2eJ55S9luzlP2WTMX10CMapR5ddFNzmluDADgdBiwWMW6boii0//DHGMJB3jFW8iBX8+y67rh90r6/nLzfrgag898foHfNS9oEH0V0I+doSwAAp0NCEIS4bb63NuHfuBnBZES48XsgCNy5sn7I3FvKkmvJ+uW9ALR991YCO3ZpE36U0I2c5tZBOR9HURRcK/8TAMuX5rF81eVUzLDj88e48dYPh5wn/fYf4lxyHcRitC5aTKT95OiHHyX0I+eY2qxlpBnjyn1vbKRv3wGEZBu5v/ovJCmJZx+rQBIFNrzp4oWX2uL2FwSBvEdXYZ5biez10rLwKuS+Ps2u42yiHzn995ypZclx5Z5fPQqAZW4lpmlTASifYueBe6cAsOKOfew7cCruGMFopPCl5xAzMwgfbKDtpptR9DW/OyL0I6e/5lTMcAyUhT78iEDdNgSjkaz/+Fnc/j9aUcSiBVmEwwpfvWL7QLN4Gikrk8JX1oDBgLdmPe5Vq0f/Is4yupHTcTIEQMV0+0BZ91N/BMBUPhXL3Mq4/QVB4IXfzaZiup2e3ijnfX0LDY3+uH0sc2aR//vHAei87wFOvbJ+FK/g7KMLOX5/FK8vhpgEk4usAMiBAL1r/wKA8zvfHtKDA7BYRDbWVFE+JRl3V4SqS+t4baMrbh/nVd8i8+d3AXB86Qr8294f5as5e+hCzkm3+qzGZpMGxjjeDW8g+/xIOdmkfe/GMx7rdBioe+1cvlyVij8Q4/LqXSy77UM63YOdgIyf3Ynz2qshGqV14dWEPtg3uhd0ltCHHJf6i7RZxYEa0vunFwEwzShHMBrPeCyoUt+uqeKnP5qEmATP/bmN0nmb+cWqBtyesNqDe/I32C65CCUYoumSbxLcM7QbrjdG5TG17PMR2LGb4Jbt9B06jJiRTtY9dyOmpQ67f82rHVyzbA8lxVYObL2AWHcPhyaWQyxG4bq1JF9y4Yjfe//BUyy9+QP2HVDn6YwGgSsW5bD48lwumpeM+zvVBN6tQzCbKVy3Btv5552NS/5U+vpiNDYHmJBnxmE3jOgY6dN3+XQix44TqNtGYOv7BLbtoG9/PXGLAADvq29Ssvs9kmy2IccP1Bxbf5P2+lsQjSIVTsB28Vc+U5bpUx3s2nQ+b9V2ctfKeuoP+1lb087amnYkSaBq5u3cXgQlzXU0X3YlKQ89SN4Prh/2nvZZUBQFry/KsbYQLceCNLUEaGwOcPCwj8amAC3Hg8gyzCy3s2vT/BGd83PJUWIxejdvoeW59cTe2YT5ZPwKmSgiDWIRe8QyTkhZLIu+Ttqx43Q/t4b07y8bcr6OfjnJ1n45614FwFRWipD0+VreBRdlsuC9TBqO+Hjo0aP8bbObDlcfdXv8bFWWc4fRRnX4LXrvvJsXV77Ci1NW4JyQTm62mbRUA2aziEESEEUBQaD/n0AoFCMYlDnli9LdE6GrO0y7q4/jbSH8gdgZ8wgCZNhhQs7Iag18BjmKolD/1+00P/p/pO/ciD3cgwEwAF4s7JamsFcs5QOxhANiMSHBRFKSWoGssQC38DLda2uGlXOyU+0QZGaakPv68L29GYCUm5aM+ELORFlJMk+vngVAp7uPdW+4qH3PzWsHV9DWNInb3L9nQWALc/buZ3X9Yp43nIcifL4PhCQJ2G0iDruE02EgLdXA7MkiFwv7KD5Ui/xOLRbPOcArIzvfSN/4L5f9hPJ3/0hR/+vWpGy2m2bRmDaD9knzyC5OZ0pJMhcWWykptjGxwILNKnLoiJ+7rmyG9pfp2Xdk2HO3949xCvLNBLfvRAkEETPScSz82kjjjYjMDBPLqwtYXl3QX3I+wYPVHL3yejKbG/ll6Cn+1bqR2rLFbE+pIhSTkGUFBVBkUABJAqtFwumQSEsxkJtjoqjQysxpdkqKrVgsIrGTLnx/q8W74XV8T9eihEKcrlNyIDB8uGEYsRzBbOJkUhr1tqm4S+fxlfuu556L8hHFT26rKyucXHbbxXj/7UHsgS76GpswTS6O2+ekS605kyZa8W2sAcAwsRDBMPIm4PNimVpG+Ud1dD35BzpX3k9GVxPf3vbfXON04PjWImyXXIhl7hwMhQXDNrGx3lOEDzcS2rGZ3t/uoX3r+/QdOBi3j5SXi3FaGak3LcXxz98YcTZNFhUebQmwZcYCKmMNpD7yMLkrlsZtn3xOLcfaQqz/01xKV/4Lod17Sb35u+Q+fP9oR4tDDoXwrH6crieeIebqjN9oMCDlZCMmqx0aORAk5vEg+/xDTyRJGCbkYyydTOqyJdgXfh1BFIfu9ylotuLzmfwrOLe7jrZvLufStQ/EbXMWv0kwKHNwUyWhf6oAYNL2WszTp2kRbViCuz+g67EnCWzbQdTlQvGfoTkySIhOJ2J6GlJONtbzz8V57WKMxYX/cA/wrHSlR4KQkw3d0Hkwforf748SDMpIokBq637aZRkpNxtT+VStog2LpbKC/GceG3gd83oJNzQS6TgJkQhJNhvSxAIMhRMQzeZRyaCZHGNBPtRD2N0VV97pUe83FksSkR071FD5ef/wp+5sI9rtWM6ZjUXD99Rs+sY5OR8AKRTfRnu61EfNVotIYKs6KWkaw+ZMT2gmJ3u6uig9OeKNKz9dc6wmgdCuPQA4rrpcq1i6RjM5hXOKAEiJnYp7Kunu6u9GJ7mQvT6SkpNJvuDLWsXSNZrJySpXB35pSi9NTYO1x91fc6ZEmgAQc7I0Gd+MBzSTk2Qy0WtwIiGz7+36gfLTzVpJSF3rLGUPXe35RUXT5zl+SxoArp2HBso8p5u1oCrHMnuWlpF0jaZywjYnAL6WjoGyTk8YFIXcXlVO8mWXahlJ12gqR7Grcvo6B8c6bk+YLKUbU58fwWLGcl6VlpF0jaZyhFT1Sajc0ztQ5u6KUBo7BoCYkTFqo+3xiKZyDP2PqcWPDUTdnjAlsjqlI2WmaxlH92gqx5Sl/vKNYXUSMRKR6e6JUCqrNUfKzdEyju7RVI41t19OTH24dnoAWqqoNcdaNU/LOLpHUzn2fFWOJaYuve30hBEUmYmxdgBsn2GVzRcBTeWkTlQHmFZZbdbcnjA5ShdmJYxgs2KeUa5lHN2jqRxHf81xKH78/iidnjBFslprxPQ0BEmzJxjjAk3lSOnqDIFD8dPYHMDtDlPc36SJaWlaRhkXaCpHTE0BwKn4OHzES6cnTLF8AhgUl2AQbQehBgNByYqIQuv+E7g9YYr7mzVzxUwto4wLNF/IHjar/znKU38M18fuOYlpm6FoL8euzhIEW9vwdXSTofQSS5Kwzf+S1lF0j+Zy5HS1Oy24OjCcaAEgkpyC6HB80mFfSDSXk5SdDYDF14WpvRkAQ3qK1jHGBZoPLEz56vyZGPSSF1VX3thyEj214dBcjq0oDwBLxE+Kon5FirGw4JMO+cKiuZyUkny8QKbcQ5qiPtcxV87WOsa4QHM56WUT8AL5soscpQsZgeSvXqB1jHGB5nKMeeo9p0BRV/F7TBmYpk3ROsa4QPPempieRlQY/EwEHNm6WxetFzSXIwgCUm72wGtLXqbWEcYNY/I9BBk3XDfwc9mixLTNmRiTByiZP78LwWTkVM0G0pdVj0WEcYHu/pZBgkF08fUqCYYnIUfHJOTomIQcHZOQo2MScnRMQo6OScjRMX8H4nLdnhlIoHUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x86.4 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5,1.2))\n",
    "ax =plt.subplot(122)\n",
    "ef.plot_shaded_errorbar(ax,t_STA,max_waveform[label2==0,:],my_pal_FSRS[0])\n",
    "ef.plot_shaded_errorbar(ax,t_STA,max_waveform[label2==1,:],my_pal_FSRS[2])\n",
    "ax.axis('tight')\n",
    "ef.clean_axis(ax)\n",
    "#plt.xlabel('Time (ms)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.05, 0.6)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAABfCAYAAAAjzQNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS10lEQVR4nO2de1zUVd7H33PhMtxxgGoVy0S7PJttpOlq+SRaUbb2QlHAmFVxn7zmrre8ZqkhaXlpvayVmsrDqk9l2eajGVAiJpr3NS94R1tXZkCHYQaYGTj7xzi/QBjugzLN+/Xy9RrO7/zOxc/ve37n8j2/IxNCCNy4LPI7XQA3zsUtsIvjFtjFcQvs4rgFdnHcArs4rVLgs2fP3ukitBpapcBWq/VOF6HV0CoFdlN/3AK7OG6BXRy3wC6O0hmJWiwWpkyZglarJSIigrlz5yKTyQBYunQpe/bswdfXl86dO/Pmm2+SnZ3NsmXLUCgUzJo1iy5dujijWL9OhBP4+uuvxfLly4UQQsyfP1/s3btXujZ69GhRWFhYJX5sbKwoKioS169fF8OGDasz/ZMnTzZreV2ZOi340KFDHD58GL1eT3BwME888QSRkZG13nPs2DGee+45AHr27MnBgwfp2bMnAJcvX2bWrFkYDAamTp1Khw4dUKlU+Pv74+/vT3FxMVarFaXSKY3Lrw6H/4uZmZl8/vnnRERE0LFjRx599FEMBgO7d+9mzZo1DBw4kH79+tV4b3FxMb6+vgCoVCpMJhMAQgj+8Ic/MGLECHQ6Ha+//jqrVq2S4gJ4eXlRWlqKn59flTS3bNnCli1bAJg8eXLTav0rwqHAMpmMlStXVguPjo4GbA+AI/z8/CRRTSZTFbE0Gg3e3t60a9cOpVKJp6enFBegrKwMlUpVLc24uDji4uIAOHXqVF31cnMLh73oPn36AFBaWkp+fj6FhYWkpaVx/fp1AKKiohwm+thjj7F//34A9u3bx+OPPw6AXq8nISGB8vJyCgoKEEKgVqsxmUwYDAa0Wi0qlQqFQtFsFWwIVl0BuqUrsOoK6h3/3zPf5tLLsZTlnmt0Ok6lrpf0qFGjxJ49e8SkSZPExo0bRVJSUp0vdrPZLCZOnCji4uLE9OnTxcWLF8XChQuFEEKkpqaKQYMGifj4eHH48GEhhBB79+4VsbGxYuDAgVJYbTirk6Vdslz8pAoV2iXLGxT/J1WouNh/UKPTcSZ1CpyYmCjMZrMkbHx8vNMLVRfNJbBFqxPaJcuFRaur8e/63H/huQHiJ1WouDJqgrgUkyBKz5xtcDrOpM6JDm9vbyZOnMjzzz/Pp59+SkBAQEs0LC3CzdRN5M+ax83UTQAoQ9SETBwPUGMTa296y3LPoVu6AoDwv68jLHkO1qs/Y9yZzr+nz6k1nRZvvut6AkpKSkRubq4QQojc3FxRVlbm9KeuLpxlwXYcNbH28Nzf/V66bk/DmPOjuBSTIIw5P4prM94SF54bIFl35XRauvmuc7C5detWtm3bhtlsRgiBTCbjiy++aIlnz2lYdQXcTN1EkCZBsrTK4f79bSOFIE2CdK0s9xzFGd/j3b0rpfsP4tGxA/79o6VWICx5Dm0/Ws7Pr72OcWe6dJ9PVO8q6dh/Vw5zJnUKvGXLFtLS0qqNS1sblUW1iwJUEbhyuD1ekCYBZYiaf0+fgykzC4+OHQCwnL+IYftOfHr2wOOhCHx69uBm6iaMO9ORBwchU6ko/9e1auWo3Hy3BHUK3KNHD/Ly8mjfvr0U1hrFtotn3PMD9777i4iVqWxdtz8E9747j2tmM16dIhD9+lD20ykseVe5um4c1vMXuTQwAe8HO4C/HxU3biIzW0DljSkzi59fe522Hy1HGaJu2UpTD4FLSkqYPXu2NNskk8nYuHGj0wvW3ARpEjDu+QHjznQMz/Ss0YqUIWpJ3MrNtFVXwI31/4v5wiVMmVn4RPWmJHsfJdn7AJCpVIgbekoPHZXSEkaj7ZqvD8ad6dxM3dSilivVqa4IeXl5bN26tSXK4lSUIWrafrRcanZvx96EVxhN6Ba8D1ClJ1y4bJUU17vLbxFmi01gby+UHe7HcvI0eHmBUgm3xCU4kHafrKbkwCEqjCasuoIWt+I6h0lhYWGsXr2a9PR0MjIyyMjIaIlyOQX7+6+m/+TKTXJY8pwqD4F//2hkgbeGh/5+gAyv/3rY9ndpmU1cmQzKysBoRB6iRn7vPSj8/FAEBiL39UG34H1pONaS1GnB7du3x2q1cvr0aSmsb9++Ti3UnSBIk0CF0URF5Xnx3HO2ce09YQh9kS3QUEzhspXIbp8vr7SHr6KgUPo7b9CrBP0xgZCZU/DvH41u6Qqp49YSOBT48OHDREZGMn58ze+NgwcP0rVrV6cVrDnQFZhZte4SIGNs0v2EqD2l8E/W5jLAksUDg5/HsH0n/v2jKTl8VBri6D//EmVoqK0Z9rklpkIO5RUAiJISxxlXFrvwBoXLVuET1duW94L3qTCakPv6tIjQDgXOy8tj48aNdOnShY4dO6JWqzEYDJw8eZJjx47Rp0+fu17gDZuv8s7i8wD4+iiYPO5BKXzWkivssOhZnJOMV8Z2ijO+x5SZhTy8LRQbsZw5hzI0FGXHDljPXwT5L+I2BlNmFj49niIseQ4VRlONwzRnIBPC8f5gq9VKVlYWR44cQa/XExQURGRkJL169cLDw8OpBauNU6dO8cgjj9QZr7IFx8Xcx9ff5DMsvh0AiSN/JDOniOnD2zDhnh8x513h5oefSPfKQ9W0/79Uir/NlDpdjUIugwoBgQF03L0Tr84RWHUFFK5eC0Cb0SOdasW1voOVSiVRUVG1Lg3ezYSoPZkztTMAi1deYMb8MwBMHvcgPXuFkZlThFLdBv74Gh8N/4TnZH608ZehCAzEeuUq2oVLCBw6uGmFqLhlP/oirk2ajk+Pp2gzeiSA9OCEzX6jaXnUwq/GL8ZuuS+/EMbilReIi7kPXx8F4W296dTtO4ym9pg8nmF40Q6sRQbw9sK4Mx3T7uxmyV8WFIgpMwtTZhYAxqy9AFU6dc6gVoFLS0uRy+V4eno6tRAtQYjak2Hx7Rg54Tg7MrSs33SVz9ZHEvVKDkaT7d2aEfTffGAazNOWo8wrW0cwZYiSUuSBgVTo9Y3PXCZD3NTj3b0rfn2fBZAmSeQ+Pk2tWq04HAdv2LCB+Ph44uLiSE9PdxSt1XDmXDF9XslhR4YWpQLOnDPS68Uf6BYZKMX5Z+l9CLmCPV5PMiFkNjcC7gPA8+FO0JQ+x61ujszDg8AhAzHlHCBoVBIhM6dIzbWzcCjwN998wxdffEFaWhppaWkNStRisfDnP/+ZoUOHMmfOHCr349LT0xkyZAixsbHSDNnSpUsZOHAgGo2G+fPnN7IqtfOXmT9x5pxthslabgsrMpTz/9/qaoz/z7L7yB6zBp+o3liv54PF0izluDZpOqbMLMxnbS4+havXOnVt2GETrVQqkclk+Pj4UF5e3qBEd+3aRadOnfjggw9455132Ldvn+Q2u3z5cjZv3oxSqSQmJoZXXnmF3Nxc1q5dS3BwcNNq4wBdgRldoWOBZLIqQ1eUCkgc0paR4x+Gjc+SnznPNgXZhF2N8uAgSrL3ETwqCZmnJx73t5c6WXJfnxqXLZtjnOzQgsvLyykuLsZgMEi/7f/q4tixY3Tv3h34xS/azrp161CpVMhkMoQQyOVyyVdao9Fw/PjxJlXodnQFZkZOOM6xEwaHcW4fKP6PJpyHImwrZkGaBMKS5xD+WSryNo1/AD063A+AQt2G+7f+3dZ7p/p6MVT3NGkKDi1YLpczduxYSYQxY8YA9VtNcuQXDaBW257I5ORkBg0aBFDNV7omh4LG+kVv2HyVHRnaesfv27sNarVXlSGV3boijv7A+V79KL/yc73Ts6PqFkngoFckMduMHulwNqs5nQIcCpyamlpjeEVF3bM5tflFV1RUMHfuXLy9vUlKSkIIUc1X2mw2V+u5N9Yv+uUXwnh32Tn0htpfM2GhHuRrLfy+WzBjk+7H10chDa3sKEPUBLz0Ajc+XFd7ppWmNPHygrIyZMiqpeVoFqs5nQIcNtGVreTjjz+Wfg8fPrzORB35RQO89957+Pv7M2PGDKBmX+nmHJatWnOpVnE9PGT0eDKQfK2FF/uGMjbpAQpumMn6oZCCG2Ypnq7AzOKVF9D7hAIgDwuVrinuvceW1qMPAaC85fWhCG9L0PChAJSdPddszW5DcChwfn6+9HvPnj3S71pmNiWio6M5f/488fHxGI1GwsPDWbRoETqdjtTUVI4ePYpGo0Gj0aBQKCTrHD9+PLNmzWpilaryj2/ya71usQj6PRtKypsPsfavtl2NscMPsyNDyxtv/bKCtmHzVWbMP8Nmz358NmAJvh9/gnf3rigfaM89i+bhG90Pv2d6ARDwUjRhyXN44B+folSrCZk5hfuWvFttGbIlqHXrip3KolYOd4SHhwdLliypEvbGG7bpuBMnTlSLn5iYSGJiYt2lbSC6AjN6g+Oeb9t7vRjxajtpnhpsQp45Z+ShCF8WzX1Yimtvro2mct75Xo3vMyEMG/AS+bPmof/7pxh3pqOaOUUSURmiRrd0BboF7xOWPAevzhF4db6LPDrsPWchBBUVFdLvhg6Z7iTvr7yAodhxecPbegMy/jLjJBl7bGNRu5DD4ttJy4tgmwmbPO5BdAVm6f0cJGzW6N8/GsMzPat1mFrag7ImHK4maTSamm+4C3yy6rua9OLgA5JwlQkKVPJwJ19ARs7Bm7a4fUNZ+9cuVUR1BRrci75y5YrTCtNc6ArMbNh8lbemdeLCZRMX835ZnJfL4au0ruzdf0MaCvV9Ru2S4kItnazdu3czePBgRowYQV5eHmVlZSxcuJA//elPLVm+RmHvEO3df4OYl++tcq2iAnZ9p2NYfDtmT+7I7MkRpK7+nUuKC7VY8LJly1izZg1arZZ58+Zx8+ZNnn76abZt29aS5WsUld+jtgV/UKnklJTYx/CiylqxK+NQYD8/P9RqNWq1mgsXLrBo0aK73kXHjr1DBPB8n1DSPv0XAf4KVColT/zWH5Opgnnvna3ip+Wq1GuY1LZt21Yj7u2kLD0vvYNf7BtK+3Af6d1b2U/LVXEo8OXLl0lJSUEIQV5eHikpKdI1+yxUa2DR3IcpNllB2H6rgz0xmqyArNpUpCvicJh04MABhzc99dRTTitQfajvMMlNLRZ8p0V00zy4P2Xo4rgFdnHcArs4boFdnFbp+F5UVHTHvnZXWFhImzZt7kjeWq2W3r17N+ymFvnUSzMTExPjzrueuJtoF8ctsIvTKgW2e1e6866bWvcHu2n9tEoLdlN/7mqBa9vElp2dTWxsLHFxcc2+3aWuvFtis5ydlJSUars7G1T3Zu7JNyu1He7R0IM8mjPvmg4WaW7s39zu06eP+Pbbb6tca0jd72oLdrSJzWAwSAd5hIWFSQd5tETegFM3y9mxWCzExsYSExNTJbyhdb+rBXa0ia1yOPxykEdL5C1uHSyyZMkSUlJSeOutt5o1Xzs+Pj7SlltH5YK6635XC+xoE5uvr2+9DvJwRt5Q/WARs9nsKJlmp6F1v6sFdrSJLSAgwOkHeTTkYJGW/IZJQ+t+VwvsaBMbwMSJE0lKSmL06NFMmjSpxfIOCgpy6mY5R1y+fLlRdXdPdLg4d7UFu2k6boFdHLfALo5bYBfHLbCL4/ICJycno9Fo6NWrl7RAcODAATQaDQkJCWg0GuLj41m/fj179+6tMjMVFRXFpk22j6YUFhby6quvVkv/7bffbtRXDxYsWFBlwsJpOGmu/K5j2rRpVU5MS0xMFHq9XgghhNVqFQMGDBBGo1HExsYKIYS4dOmSGDt2rJgwYYIQQoj09HSxbNmyKml+//33YvPmzY0qz5EjR8Tf/va3Rt3bEFzegutDWVkZAAqFAg8PD0pLS8nJyeGll15Cp9MhhODo0aPVtvN89tln0iHZCQkJTJs2jUGDBpGamsqYMWPo378/V65cYf/+/QwZMoS4uDjWr18PwOOPP05WVpbT69Yq3Wabi1GjRqFUKpHL5YwfPx4vLy8iIyM5ceIE+/fvZ+bMmRw6dIjTp09z8uRJxo0bV+X+ixcvSl/uy8/PZ/HixRQVFTF16lS++uorNm7cSHZ2NpcvX2b48OFER0fz5ZdfArbtuZ6enmi1WkJDQ28vWrPxqxb4ww8/rHaaardu3Th+/DiFhYWEhITQvXt3jhw5gkKhwNvbu0rcynPAvr6+/OY3v6GiooLw8HBkMhn+/v4UFRXx2muvsWLFCjZt2lTFr1mtVqPX650qsLuJvo0nn3ySXbt2ERERAdh2WW7bto3HHnusWtzKAtf2/bDt27czbNgwUlNT2bFjh/RBV4PB4LQv7Nr5VVtwTfj5+WGxWOjRowcAwcHBlJaW1ridtl27dhQUFEjNtCMeeeQRJkyYQEBAAN26dZOWHo1GY533NhX3YkMT+O6777h27RpDhw5t8L1HjhwhJydH+oqvs3A30U3g2Wef5dSpU40aB2/fvt3hx+aaE7cFuzhuC3Zx3AK7OG6BXRy3wC6OW2AXxy2wi/Mfhz4OdS05sD0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 80.64x69.12 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(1.6*0.7,1.2*0.8))\n",
    "ax.scatter(all_features[label2==1,2],all_features[label2==1,4], color=my_pal_FSRS[2],s=1,label='RS')\n",
    "ax.scatter(all_features[label2==0,2],all_features[label2==0,4], color=my_pal_FSRS[0],s=1,label='FS')\n",
    "plt.xlabel('TPW (ms)')\n",
    "plt.ylabel('REP (ms)')\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim([0,1])\n",
    "plt.ylim([0.05,0.6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Clustering FS into FS1,FS2,FS3, and RS into RS1,RS2,RS3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_r_sub = df0[['PL_below','PL_above']].values   # Multichannel features\n",
    "Multi_X_r=ef.pre_process(X_r_sub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_r_sub0 = Multi_X_r[label2==0,:]\n",
    "idx_sub0 = np.where(label2==0)[0]\n",
    "\n",
    "X_r_sub1 = Multi_X_r[label2==1,:]\n",
    "idx_sub1 = np.where(label2==1)[0]\n",
    "\n",
    "numk=9\n",
    "[K0,wcss0,fk0] = ef.selectionK_kmeans(X_r_sub0,numk)\n",
    "elbow_idx0  = ef.findElbow(wcss0)\n",
    "\n",
    "[K1,wcss1,fk1] = ef.selectionK_kmeans(X_r_sub1,numk)\n",
    "elbow_idx1  = ef.findElbow(wcss1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAABnCAYAAACacLX/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQPElEQVR4nO2deVBTV/vHvyESBMSthVoqOgphq6KUEpdxGahWOlqLQFnFIiOKI66onY61Ci5VMiNTxBYdQ0dxQKo/HKRaW2JVFDEIRSybAlV0XqxoUQOENXl+f/iSCiQhXnIh+OYzkxlyz/bky7lPznnOOTccIiIY0DlGA23Am4pBWJYwCMsSBmFZwiAsSxiEZYkBFfbbb7+FWCweSBNYY0CEbW9vx6ZNm5CdnT0QzfcLAyasn58flixZ0mveysrKLu9jY2Px9OlTtkzTGX0Stry8HMBLodLT05GRkQFtJnJmZmaYOXOm2vT09HT4+PjAx8cHdXV1XdKeP3+OLVu29MXsfoGxsImJiUhISIBcLkdMTAyKiorw999/Y+/evX02KiAgABkZGcjIyMDo0aO7pMXExODixYu4cuWK5kpqa/tsR18YwrRgUVERRCIR2traIBaLkZOTAx6Ph+DgYF3a1wMLCwt89913iIyMxK1bt2BiYqI643vvAQMYBumzj83Pz4eLiwt4PB4AoLm5uc9G9Ya3tzf4fD7i4uJYb4spHKbRrX379oHD4SAvLw9btmyBs7Mzjh49iufPn2PPnj06M7C8vBxOTk49rtfU1MDNzQ15eXng8/k9C3I4A9pjQQyRy+V0+fJlKisrIyKiO3fukEgkotbWVqZVqqSzflUIhUKaN28eKRSKnokRETq143Vh7AqMjIxgZWUFJycntLe3o6ioCCNHjoSxsbEu/+8aWb9+Perq6pCWltYz8ciRfrNDFXo5KtAWY2NjHD58GJs2bUJVVVXXRDe3frNDJUy7enh4OBERtba20rRp05QuICgoSDf30n/R5Ao6OXz4MNnY2FBlZeW/F5l/NJ3AeLjVyUCMCrqzcuVKEBE8PT1x6dIl2Nra9rsN3WEsLJ/Px/79+5WjgmfPnuHo0aNwdnbWpX1as2rVKgCAh4fHS3HffXdA7FDCtKvrw6hAFUlJSWRjY0NVVVU6teN1YdxjjYyMIJfL4eDgAADo6OjAuHHjlC5hoOjsuf83eTL+WrYMPj4+8PDw6NfRCtDHUUFWVhZaW1sBAMOGDUNmZiaSkpJ0ZhxTVq1aha3NzbC1tcU333yDMWPGICwsDBkZGXjy5En/GMG0q/v7+/e41tHRQb6+vn26hbrzuq5AySsf7eHDh5SQkEALFiygESNGEJ/Pp+XLl9PRo0eptLSUOjo6dGTtvzB2BapueS6Xqz4oMoCMHTsWa9euxdq1ayGXy1FaWorc3FxcunQJ+/btQ21tLSZNmgRXV1e4urpi0qRJ4PP5sLS0BIfDYdQmY2EdHR2Rnp6OgIAA5bVTp05h/PjxTKvULQUFKi9zuVy4uLjAxcUFq1evBgA0NDSguLgYRUVFkEgk+PHHH1FZWYmOjg7w+XzY2dlhwoQJsLGxUb7Gjh2Lt99+W23zjIMw165dw9WrV3Ht2jUYGRkpx5FRUVE6/QJTF4TplcLCPs++6uvrUVVVhcrKSty/fx8PHz7s8nr27JnasoyFjY2NRVlZGYYOHQo3NzcIBAK4urr2Kmp7ezs2b96MJ0+ewM7ODjExMRpvN8bCDnB0i7GwnTQ3N6OkpAS3bt1CTk4OjIyMcOzYMbX5z507h3v37iEqKgq7d++Gp6enxmWawSpsn6a01dXVuH79OvLz8yGTyeDg4IDp06drLFNcXIz58+cDAGbOnImCggKNwg5WGAs7e/ZsTJ06FYsWLUJsbCxGjRqlVbnGxkaYm5sDAExNTSGTyXrkSU9PR3p6OgAgOjqamYE7djArpyMYC7t+/Xrk5+cjLS0NEokE06dPh0AgwMiRIzWWGzZsmFJMmUyGYcOG9cgTEBCgHG1IJBLlarA66uvreyw6IiAA6FZOZT5t61PBkCFDVK9eALqJrVVXV9OJEyfI39+fvL29NeY9e/YsJSYmEhHRrl27KCcnR2P+JUuW9Nq+NnnYyKcJxj22vr4excXFKC4uRklJCYgIH374IQQCgcZyXl5e+PLLLxEYGIgJEyZg1qxZTE3QaxgLu2bNGggEAri7u2PlypUwMzPTqpyxsTEOHDjAtNlBA2NhVa4zscCrM7u+5GEjnyb6PI41oBrD/liWMAjLEgZhWYJVYVXt2BaLxfD394efnx8yMjLYbH5g6fNIWAVtbW20ceNG8vDwoOzs7C5pixcvJplMRm1tbbRw4UJWovf6QJ/3Faiic8f2hAkTeqQlJyfD1NQUHR0dICIYGb2Z3oiVT6Vpx/Zbb70FANizZw98fX1VxmJf3dGdk5PDzAiGSyq6gpUeqwmFQoGYmBgMHToU4eHhKvO8GoTpLQCjr/S7sEKhEBYWFti8eTO7DalZ8+ov+sXB1dTUIC4uDk+fPkVKSgpu3bqF0NBQhIaGoqGhoT9M6Hf0fko7WJdm3syvZBUUFBTAzc0NQqEQjx8/Zr09rYXt3Eo0WNm9ezfmzJmD8vJyODo6wtfXF+fPn4dcLmenQU2D3KysLFqxYgUtWrSIAgMDyc/PjyIiIigrK0v1vn8WYLzFaMcO5Z+lpaVkZWVFTU1NRET04sULOnLkCLm7u9PEiRMpNTWV5HK5Dqz9F7XC7tixg86dO6c0phOZTEY///wzbd++XaeGqIOxsK8QFhZGu3fvVpn2+++/k7u7O7m6utKFCxd01mHUCtvc3Kzyem1trdaV7927t8eU9urVq+Tr60v+/v5UXFzcax2MhX33XSIievDgAY0aNYrq6+vVZlUoFHT69Gmyt7cnDw8POnv2LD169IhZu/9FrbAhISHU0NDQ5dqpU6fIw8Oj10o1xQr8/PxIKpXS48eP6Ysvvui1ru7CSiQSSk1NJZFIRImJiSQUCikxMZHa29u7Fvyvl9uwYQNFR0f32g4RUXt7Ox05coQ++ugjGjVqFFlbW9Onn35KO3fupMzMTLp//77WPVrtBCEkJAQrVqxAUlIS2tvb8fXXX4PH4+HkyZO9+m11sYKGhgaYmprCwsICFhYWaGxsREdHB4YM6WqGpn0FmZmZqK6uhqmpqfJ18+ZN5ObmIiUlBVwuV5n3n3/+wbFjx/Dnn3/2/mWDl8vZERERiIiIABHh3r17KCwsRGFhIX744Qfcvn0bTU1NmDx5MlxcXHDo0CH1lWlSPTc3l7y9vWnhwoV04cIFrf5Tr5KQkNClx9bW1lJkZKTyfXBwcI+7ojvauAKZTEaenp60fPnyf7+EPviAYmJilKd7dEVdXR1dvHiR4uPjNebrNWxYXFxMAQEBjHxOd2FfvHhBy5YtU7739fXtNWyorY9tbGykWbNm0erVq0mhUFBjYyNZWlpSeXn5a9utC9S6Am9vb3A4HBARWlpaEB4eDh6PBw6HgzNnzmh1a3Vn+PDhkMlkaGhoQEtLC0xNTbvcun3B3Nwc586dw/z58xEdHY3wGzcwe/ZsODo66qT+10ad4r2N67Rx4p099v79+7R//34ieule/Pz8yMfHh/74449e63jdUUF9fT1NnTqVCCCJRPJaZXWJ2lhBbGwsnJycMG/evC4b3urr6/Hrr7+ivLwcsbGxrP/jmcQKnj59irctLfVzGyePx4OVlRVmzpwJOzs7DBkyBFwuF++99x58fHwQFBTUn3a+Fpq2sPcXaoUtKChAR0cHrKysMGPGDOV1DoeDuXPn9otxfeI//xnQ5tUKKxKJUFlZiZKSEsyfP1+rh+joFYWFgLX1gDXPSjxW0zmDixcv4tChQ+Byudi8eTOmTZumsS5DPPYVfvvtN/D5fKSmpoLH4yEvL0+Z9v3330MkEiEpKQnx8fFsNK8XsCJscXGxsid2njPoxNHREY2NjZDJZMot828irCwmajpnYGNjA39/fwDArl27VJbvHitgtFW+rGzwb5Xvzp49e+jmzZtERCQWi+ngwYNE9HJK+/HHH1NjYyNJpVJavHix2vBkJ4N1qzwrrmDy5MmQSCQAgLy8PEyZMgUAYGJiAjMzMxgbG8PMzEx5NP9NhBVhvby8UF1djcDAQDQ1NcHGxgZxcXEwMTHBsmXLEBwcjKCgIAQGBr65frbPfZ5lTp48qZM8bOTThN7vKxis/M/sK+hvDMKyRL9vitMWbY7ft7S0IDo6GlKpFDweDwcOHMCIESPU1nnhwgVkZWVpXKvau3cvSkpKYGRkhPj4eFhaWjKyX297rKZpcSeZmZn44IMPkJKSAi8vL+WkQhWPHj1Cenq6xmDSlStXQERITU3FmjVrUFNTw9h+vRVW07S4k0WLFikfBKxqtbcThUKBffv2YePGjRrbzM/Ph7m5OcLDw5GVlQUXFxfG9uutsNocvzc3N4epqSmqq6uRlpYGPz8/lXUlJSXBx8en12nq8+fP8fjxYyQnJ8Pa2hqnTp1ibL/e+lhtjt8DL8OKX331FeLj4zF8+HCVecRiMfLy8tDa2oqamhqcOHECS5cu7ZFv+PDhykcIzpo1i/GiKaDHPVbdtPhVamtrsXXrViQkJGh8UGRGRgZSUlJw4MABuLm5qRQVAKZMmYIbN24AeOmKVB1O0Ra9Fbb7tFjV8XuRSISmpiZs27YNoaGhGp9Fow2dj1QJCAjA9evXlVE4JhhmXiyhtz12sGMQliUMwrKEQViWMAjLEno7QejkypUrePbsGcaMGYPq6mqEhIQo09ra2hAXF4eKigq0trZi7ty5iIqKQkZGBqRSKcLCwrRqo6KiAgqFQqfPF9f7HltWVgZnZ2dUVFTg/fff75J28OBB2Nra4sSJE/jpp59QWVnJ6BftsrOzUavjX1PS63FsVFQUCgsLYWdnh6qqKri7uyMhIUGZ/sknn+D8+fPKcGJTUxNMTExw9uxZSKVSODk5QSwWY9u2bZBIJBCLxVi1ahU2bNgAIsKYMWOwc+dOfPbZZzAzM8Px48eRmpqK3NxccLlcbN++HQ4ODliwYAGsrKzg7++PzMxMtLS0wNzcHPHx8Wofi6XXPTYxMREzZsxASkoKBAJBF1GBlw/ZfTVGa25urjbC1cnt27dhb2+P48ePY968eeBwOFiyZAk2bNiAuro63L17F2lpaRAKhRAKhQBebgsViUSwt7cHj8dDcnIyli5dqvEcsF4LGxkZiYKCAgQFBaGwsBDr1q3rkq5QKKBQKJTva2pqUFpaqrKuzhtz7ty5eOeddxAeHo7r1693eRDFX3/9hbKyMoSGhmLr1q2QSqUAAGtra/B4PDg4OGDOnDmIjIzEmTNnNO5G12thQ0NDsW7dOkRHR2P16tU9eqynpydOnz4NAJDL5RAKhV18pYmJifL3FTt/e/HmzZuwtbXFsWPHYGZmhhs3biiPBIwbNw4CgQApKSlITEyEl5cXACjFv3PnDrhcLpKTkzF16lT88ssvam3Xa2ErKirg6OiIu3fvwt7evkf62rVrUVRUhODgYHz++edwdnZWBlIAYNKkSZBKpQgJCVEKy+fzkZycjJCQEDx48AACgQCOjo5ISEjA6NGjYWlpiZCQEISFhWHixIld2hs/fjyys7MRGBiIy5cvd2mrO3r95TWY0eseO5gxCMsSBmFZwiAsSxiEZQmDsCxhEJYl/h+QYy+DBvFnYgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 46.08x69.12 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAABnCAYAAACacLX/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQlElEQVR4nO2deVBT19vHvyEQyiKlWDEuuLK7VhQVtVTUSou1bGU11WFcaKXWlY5DrYBirczIDGLr/GroIAqiDlYRiyW1iiKCUEwHEbRUcUGxjtAIYQt53j98SQ0kISa5EJx8/iL3PPecJ9+ce+45z33OhUVEBAM6x6i/HXhdMQjLEAZhGcIgLEMYhGUIg7AM0a/CfvvttxAIBP3pAmP0i7AdHR3YuHEj8vPz+6P5PqHfhA0MDISfn1+vtrdv35b7HB8fj8ePHzPlms7QStibN28CeCFUVlYWsrOzoc5CztzcHB4eHkrLs7Ky4O/vD39/fzx58kSuTCqVIjw8HJ2dndq4zjgaC5uSkoLk5GR0dnYiLi4O5eXlePz4MXbt2qW1U8HBwcjOzkZ2djZsbGzkyrZt2wYiwo4dO1RXUlentR/aYKzpieXl5eDz+Whvb4dAIEBBQQE4HA7CwsJ06V8P2Gw2MjIyMG3aNMybNw8LFixQbDhiBNCPYRCtx9iSkhJMnjwZHA4HANDS0qK1U73B5XKRnp4OHo+nt+MtS9Po1u7du8FisVBUVIQtW7bA1dUVBw8eRGNjIxISEnTm4M2bN+Hi4qKwLDY2FgUFBcjPzwebzZYvZLH6tceCNKSzs5MuXLhAlZWVRERUXV1NfD6f2traNK1SIV31K0IikZCXlxd98803PQtXrdKpH6+KxmOskZERbG1t4eLigo6ODpSXl8Pa2homJia6/N1VwmazceTIEcyYMQNjx47FihUr/iv83//6zA9F6OWs4FXgcrnIz89HTEwMMjIy/itwc+tTP3qgaVePiIggIqK2tjaaOXOmbAgIDQ3VzbX0/6gaCl6moqKCuFwuHTt27MUBzb+aTtB4KOiiP2YFipgwYQLy8vKwePFicDgcfNwvXvyHxsI6ODjgu+++k80KGhoacPDgQbi6uurSv1diypQpyM3NxYcffojFNjZ4o988wcCeFSjj6tWrxOVyae7cufTTTz9RU1OTTn1SB41vXkZGRujs7ISTkxMAQCKRYNSoUbIhoT+ZOXMmHqxciY0bN+LEiRMYOXIk1qxZg8uXL6Ojo6NvnND0F9m3bx+tX7+exGIxERHV1tZSVFQU/fDDDzr71Yk067FEJHfzevDgASUkJNDUqVPJ0tKSFixYQHFxcfT7779Tc3Ozjjzt1rymJwYFBfU4JpFIKCAgQCuHuqMLYV+moaGBzpw5Q9HR0TRr1iwyMzOjSZMm0fLlyyk5OZkKCwvp+fPnWnj8Ao1vXooueTabDVNTU62uIKaxtraGj48PfHx8AABtbW2oqKhAWVkZ/vjjDxw6dAg3btyAjY0NnJyc4OzsDCcnJ4wbNw52dnYYOXIkbGxswGKxVLajsbDOzs7IyspCcHCw7Njx48cxevRoTavULaWlapmZmprCzc0Nbi8tKKRSKe7du4fq6mpUV1ejqqoKeXl5uH//Ph48eIC2tjaMHDkSt27dUlqvxsJ6enri0qVL8PHxgZGREYgIXl5eiI2N1bRKvcHIyAhjxozBmDFjsHjx4h7lTU1NePDggco6NBb2/PnzqKysxJAhQ+Dm5gZ3d3e88847vc4KOjo6sHnzZvzzzz+wt7dHXFxcr5eVRkyfzlh0y9LSEs7OziptNA4bdtHS0oKKigpcv34dBQUFMDIyQlpamlL73Nxc3LlzB1FRUdi5cye8vLxUPqZRFTZUST+HDbVa0tbU1ODKlSsoKSmBWCyGk5MTZs2apfIcoVCIRYsWAQA8PDxQWlqqUtiBisbCzps3D1OnTsWSJUsQHx+Pt956S63zmpqaYGFhAQAwMzODWCzuYZOVlYWsrCwAwKZNmzRzcPt2zc7TERoL++WXX6KkpASZmZkoLi7GrFmz4O7uDmtra5XnWVpaysQUi8WwtLTsYRMcHCybbRQXF8ueBivj2bNnPR46IjgY6HaeQjt161OAsbExHBwcFBdqPRMmopqaGjp8+DAFBQWRr6+vStvTp09TSkoKERHt2LGDCgoKVNr7+fn12r46NkzYqULjHvvs2TMIhUIIhUJUVFSAiDB9+nS4u7urPM/b2xtfffUVQkJCMHbsWMydO1dTF/QajYVdu3Yt3N3dMWPGDKxevRrm5uZqnWdiYoK9e/dq2uyAQWNhMzMzdemHUl5e2Wljw4SdKrSexxpQjCE/liEMwjKEQViGYFRYRRnbAoEAQUFBCAwMRHZ2NpPN9y9az4QV0N7eThs2bKD58+dTfn6+XNnSpUtJLBZTe3s7+fj4kEQiYcKFfkfrvAJFdGVsjx07tkdZamoqzMzMIJFIQEQwMno9RyNGvpWqjO3BgwcDABISEhAQEKAwFvtyRndBQYFmTjAR430FGOmxqpBKpYiLi8Mbb7yBiIgIhTYvB2F6C8DoK30ubGJiIgYNGoTNmzcz25Caz7yYok8GuNraWuzZswdPnz5Feno6rl+/Dh6PBx6Ph+fPn/eFC31Pf989e6O3vILOzk76+eefSSqVyhf081cb8Lfk2NhY+Pr64tSpU73a3r17Fw8fPkRjYyPzqUbq/gKtra1M/sBKUdVjT548SXZ2dpSWlkaOjo7U3t7+X2G3r7Z//36ytramYcOGkZWVFbHZbDIxMaFJkybRkSNHdD6fVilsTk4OrVy5kpYsWUIhISEUGBhIq1atopycnJ6XHkMoE7ayspLefvttKi4uJqlUSosWLZI9mSAiou3bZX/evXuXBg8eLFeXVCqllpYWOnfuHM2ZM4ccHR0pLS2NOjo6dOK3UmG3b99Oubm5PZLGxGIxnTlzhrZt26YTB3pDkbCNjY3k6OhIfD5fduz69es0dOhQamxslLOVSqXk7e1NO3fuVNqGVCql8+fP03vvvUfjxo2jxMREEgqFWnUepcK2tLQoPF5XV6d25bt27eqxpL106RIFBARQUFAQCYXCXuvoLmxnZyd99NFH9Pnnn/ewXbFiBW3duvXFh2HDiIgoPT2dJk+eLD9MqKCgoIAiIyPJ3t6ebG1tKTQ0lPh8Pl2+fJmuXbtGQqGQqqqq6M6dOyrrUSpseHh4j6y748eP0/z583t1TlWsIDAwkEQiEdXX19Py5ct7rau7sHFxcTRnzhyFCc73798nGxsbqq2tJQKovr6ebG1tqbS0tNd2FHHnzh06ePAghYaG0uzZs8nNzY0mTpxIDg4OZGdnp/JcpcKePXuWgoODqaGhgZ48eUKrV6+mqKgoqq+v79Wh5uZmKiwspOTkZDlhRSIR8Xg82eeAgACFY9rRo0fJz8+P/Pz86PLly3JlO3fupEePHiltOyYm5kUbAAUHB9OWLVt69ZcJVN68CgsLydfXl3x8fCgvL++VK+8ubF1dHUVGRso+h4WF9ZqL+qr5sSKRiLhcLj0cNozs7e0ZSyzuDZXzWA8PD8TFxcHS0hJTpkzRempnYWEhl/nS1tYGMzMzret9mUGDBmH79u0Y8egRfvzxR7WfHusapbECX19fsFgsEBFaW1sREREBDocDFouFkydPatSYlZUVxGIxnj9/jtbWVpiZmfXcA6sDVq5ciY9zczHsvfd0Xre6KBU2OztbZayUiNROv6ytrUVWVhaio6OxYcMGREREQCqV4uuvv351j9XA2NgYw86cYaRudVH6+Ds+Ph4uLi5YuHChXMLbs2fPcO7cOdy8eRPx8fGMO/japXFyOBzY2trCw8MD9vb2MDY2BpvNxogRI+Dv74/Q0NC+9HPAoVTY0tJSSCQS2NraYvbs2bLjLBYLnp6efeKcVjx82K/NKxWWz+fj9u3bqKiowKJFi9R6iY5eUVYGDB/eb80zkmKkap/Bb7/9hv3794PNZmPz5s2YOXOmyroG6hjLSDz2119/hYODAzIyMsDhcFBUVCQr+/7778Hn83HgwAEkJSUx0bxewIiwQqFQ1hO79hl04ezsjKamJojFYlnK/OsIIw8TVe0zsLOzQ1BQEAAofXdW9z0IGqXKV1YO/FT57iQkJNC1a9eIiEggENC+ffuIiOjff/+l999/n5qamkgkEtHSpUuVhie7GKip8owMBZMmTUJxcTEAoKioSBZnMDU1hbm5OUxMTGBubi7bmv86woiw3t7eqKmpQUhICJqbm2FnZ4c9e/bA1NQUn376KcLCwhAaGoqQkJDXd5zVus8zzNGjR3Viw4SdKgyp8gwx4PMK9BWDsAzR50lx6qLO9vvW1lZs2rQJIpEIHA4He/fuxZtvvqm0zry8POTk5GD//v1KbXbt2oWKigoYGRkhKSkJQ4YM0ch/ve2xqpbFXZw6dQrTpk1Deno6vL29ZYsKRTx69AhZWVkqg0kXL14EESEjIwNr165FbW2txv7rrbCqlsVdLFmyRPYiYIlEAmNjxRegVCrF7t27sWHDBpVtlpSUwMLCAhEREcjJycHkyZM19l9vhVVn+72FhQXMzMxQU1ODzMxMBAYGKqzrwIED8Pf373WZ2tjYiPr6eqSmpmL48OE4fvy4xv7r7RirzvZ74EVYcevWrUhKSoKVlZVCG4FAgKKiIrS1taG2thaHDx/GsmXLethZWVnJXiE4d+5cjR+aAnrcY5Uti1+mrq4O0dHRSE5Oxvjx45XWlZ2djfT0dOzduxdubm4KRQVevBvx6tWrAF4MRYo2p6iL3grbfVmsaPs9n89Hc3MzYmJiwOPxVL6LRh26XqkSHByMK1euyKJwmmBYeTGE3vbYgY5BWIYwCMsQBmEZwiAsQ+jtAqGLixcvoqGhAVwuFzU1NQgPD5eVtbe3Y8+ePaiqqkJbWxs8PT0RFRWF7OxsiEQi+f+LoIKqqipIpVKdvl9c73tsZWUlXF1dUVVVhQkTJsiV7du3D+PHj8fhw4dx7Ngx3L59W6P/aJefn486Hf83Jb2ex0ZFRaGsrAz29vb466+/MGPGDCQnJ8vKP/jgA5w9e1YWTmxuboapqSlOnz4NkUgEFxcXCAQCxMTEoLi4GAKBAGvWrMH69etBROByuYiNjcXHH38Mc3NzHDp0CBkZGSgsLASbzca2bdvg5OSExYsXw9bWFkFBQTh16hRaW1thYWGBpKQkpYnNet1jU1JSMHv2bKSnp8Pd3V1OVODFG5ZfjtFaWFgojXB18eeff8LR0RGHDh3CwoULwWKx4Ofnh/Xr1+PJkye4desWMjMzkZiYiMTERADA06dPwefz4ejoCA6Hg9TUVCxbtkzlPmC9FjYyMhKlpaUIDQ1FWVkZ1q1bJ1culUohlUpln2tra3Hjxg2FdXVdmJ6enhg6dCgiIiJw5coVueTqv//+G5WVleDxeIiOjoZIJAIADB8+HBwOB05OTnj33XcRGRmJkydPqsxG12theTwe1q1bh02bNuGzzz7r0WO9vLxw4sQJAEBnZycSExPlxkpTU1M8ffoUwH//e/HatWsYP3480tLSYG5ujqtXr8q2BIwaNQru7u5IT09HSkoKvL29AUAmfnV1NdhsNlJTUzF16lT88ssvSn3Xa2Grqqrg7OyMW7duwdHRsUf5F198gfLycoSFheGTTz6Bq6urLJACABMnToRIJEJ4eLhMWAcHB6SmpiI8PBz37t2Du7s7nJ2dkZycDBsbGwwZMgTh4eFYsWIFxo0bJ9fe6NGjkZ+fj5CQEFy4cEGure7o9c1rIKPXPXYgYxCWIQzCMoRBWIYwCMsQBmEZwiAsQ/wfLwwejXdZveIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 46.08x69.12 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(0.64,1.2*0.8))\n",
    "ax=plt.subplot(211)\n",
    "#Plotting the Elbow Curve by k(number of clusters)\n",
    "plt.plot(K0,wcss0/max(wcss0),'k-',markerfacecolor='None',linewidth=1)\n",
    "plt.plot([K0[elbow_idx0],K0[elbow_idx0]],[0,1],'r--',linewidth=1)\n",
    "plt.ylabel('WCSS')\n",
    "#plt.title('Elbow Method',fontsize=9)\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xticks(range(0,8,2))\n",
    "ax=plt.subplot(212)\n",
    "#Plotting the Density Function fk by k(number of clusters)\n",
    "plt.plot(K0,fk0,'k-',markerfacecolor=\"None\",linewidth=1)\n",
    "plt.plot([K0[fk0.argmin()],K0[fk0.argmin()]],[0,1.5],'r--',linewidth=1)\n",
    "plt.xlabel('# Clusters')\n",
    "plt.ylabel('f(K)')\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticks(range(0,8,2))\n",
    "plt.ylim(0.8,1.2)\n",
    "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.5,wspace=0.5)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(0.64,1.2*0.8))\n",
    "ax=plt.subplot(211)\n",
    "#Plotting the Elbow Curve by k(number of clusters)\n",
    "plt.plot(K1,wcss1/max(wcss1),'k-',markerfacecolor='None',linewidth=1)\n",
    "plt.plot([K1[elbow_idx1],K1[elbow_idx1]],[0,1],'r--',linewidth=1)\n",
    "plt.ylabel('WCSS')\n",
    "#plt.title('Elbow Method',fontsize=9)\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xticks(range(0,8,2))\n",
    "ax=plt.subplot(212)\n",
    "#Plotting the Density Function fk by k(number of clusters)\n",
    "plt.plot(K1,fk1,'k-',markerfacecolor=\"None\",linewidth=1)\n",
    "plt.plot([K1[fk1.argmin()],K1[fk1.argmin()]],[0,1.5],'r--',linewidth=1)\n",
    "plt.xlabel('# Clusters')\n",
    "plt.ylabel('f(K)')\n",
    "#plt.title('Density Function')\n",
    "ef.simpleaxis(ax)\n",
    "plt.xlim(0,8)\n",
    "ax.set_xticks(range(0,8,2))\n",
    "plt.ylim(0.8,1.2)\n",
    "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.5,wspace=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc0=fk0.argmin()+1\n",
    "nc1=fk1.argmin()+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# K-means\n",
    "kmeans_sub0 = KMeans(init='k-means++',n_clusters=nc0,n_init=1000,random_state=110)\n",
    "kmeans_sub0.fit(X_r_sub0)\n",
    "label_sub0 =  kmeans_sub0.labels_\n",
    "#print(label_sub0)\n",
    "\n",
    "kmeans_sub1 = KMeans(init='k-means++',n_clusters=nc1,n_init=1000,random_state=1001)\n",
    "kmeans_sub1.fit(X_r_sub1)\n",
    "label_sub1 =  kmeans_sub1.labels_\n",
    "#print(label_sub1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-4.0, 7.0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAG8AAABuCAYAAAApmU3FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAV+UlEQVR4nO2deVjU1f7HX8O+LwIuKHTZXHFNQXgy3G7pzYuZuefFss300QrLvNflhoqmaf5u/szuL81yKTW5V7BFMnFBQRQRNyAUYhFFRpABZoBh5vz+ICbJZWZghsXm9TzzzHyXc87nO+8553zOOhIhhMBEu8SstQ0w0XRM4rVjTOK1Y0zitWNM4rVjTOK1Y0zitWNM4rVjLHS5qbq6mkOHDpGWlkZ5eTmurq4MHDiQ0aNHY2tra2wbTTwAibYell27dnHu3DmCg4Px8/PDzc2NiooKMjIySE5OZsCAAfztb39rKXtN3IVW8TIzM+nZs2eTr5swHlrFM9F20anOA9iwYQMJCQlYWFgghEAikfCf//zHmLaZ0ILO4iUnJxMbG4tEIjGmPSb0QGfxnnzySX788Ue8vb0150x1Xeuis3gZGRlkZWXh4OCgObd69WqjGGVCN3QWr7q6mq1btxrTFhN6orO3OX/+fJydnfH399fUe8Zo33333XfMnTuXuLg4AgMDDR7/o4TOOW/EiBEt4qz4+PiwcuVKVCqVxqs1cX90Fi8nJ4ejR49ibm4OgEQi4dlnnzWKUTNmzADgnXfeoUuXLrz11lsmEe+H0JFJkyYJtVqt6+1N5sqVK5rPeXl5on///uKll14SNTU1Rk+7vaHzqEJDUyEzM1PzMjbe3t4kJiZy+/Zttm3bZvT02httvqng4OBATEwMAKdOncLZ2Zk+ffoYPd32QIs0FWpraxkxYgS+vr4ALFq0SC9P0sysvoDIz89n/vz5bN++nb/85S9NsuVRQmfx7O3tWbp0aZOaCtnZ2YwZM4alS5c2zcpfmTp1Kt7e3kycOJHo6GhefPHFZsXX3tGrqdBUsrKyuHDhAjNmzCAwMJBFixZpcpO+hIaGkpycjEKhoLq6GjMzM6ysrJpsW7tGm0dz4MABcevWrWZ5RUlJSSIxMVEIIURUVJSIjY1tdP3rr78WEyZMEBMmTNDcpwvr168Xw4cPF1KptFn2tVe09rAkJSWRnJxMcXExrq6uBAcHExQUhJ2dnc4/EIVCgaWlJRYWFhw7dozU1FTefvvt+96bkZFBr169dIpXpVKxePFiYmJiiIuL0zkcgKJSRUaKjF5BTtg6mOscri2htdgMCQkhJCQEgPLycpKSkoiOjkapVPLBBx/olMjGjRvp27cv48aN4+zZswYbjTA3N2ft2rX07t1bb/EyUmScirsNwKCRrgaxp6XRuc774osvcHZ2RiqVcuXKFYKDg3VOZPbs2bzzzjvs2bMHHx8fxowZ0yRjH8SsWbMAOHz4MBkZGcybN09rj0yvIKdG7w9CWw5t1Rysa/k6ZcoUoVarxcyZM4UQQvNuaO7uYdGXnJwc0adPH/Haa6+J2tpag9iT+lOp+PjNbJH8vbT+9Z1UyCvq7rme+lOpQdLTB51dPiEEP/zwAwEBARQUFCCTyYz5m2oSPj4+nDp1isLCQubNm2eQOHsFOTHkaVdu5FRz5lAZZ+LLyEiRNboe+lc3rTnYGOg8JJSSkkJCQgKvvPIK3333HT179mTw4MEGN0gfh+VBqFQqpFIp9vb23Lhxg4CAgGbFd/r725yJL8OxgwV+/ex5fFSHNuHktLnZY4YQr4HDhw8zY8YMdu7cyZ///Ge9wioqVVxIvAMC6pSCtIQ7AAx5ypXgsW4Gsa+5PNKzx0aPHs2+ffuYPHkyi979B2EDXnioY1FWXEviASlPjHcn93IVZw6VAdAtwBZPPxuKrlXDr35QW2hqPPKzx5588klOnTrFF598y6mC2wgheHxUh/vee2x/CYXZClR1JYRN9CA/S46yRk1htoLAUEcsrczoPtARaBtNjTY9JGQofH19eW/5G9j5ZPPemmkUFUg5d6QMRaWq0X3u3aw077mXqyj8WUF1Vf09d6R15GXISTwgRVGpalVHpYE2PyRkKGwdzJkcMYqfjscSEjKUiNGfAIMb5ZrHR3bAzsFCI0h+ppzCbAXWthJs7CU4u1uQlyHnwok7WFqbaYrM1ipC/1Czx35OrSLEKxJvzwDiT0WzOOgw8JtzoqxRY2ldXxjZOpjTxdeGwmwFNQrB1TQ5AC4dLamuUnEmvgxljZrgsW6tVoS2yJBQa1Aml3Hg0jHGB4YBcODSMZ7uO4yBSnvs7Xry9spZ1Koqidl1kD85Pq3xJgFu5dfwxHh36moFHbtZUSatRVkNVrYS7txSYudYn7tu5FZrilDQ3ltjaFpkSKg1OHDpGB8d36U5/uj4LhTKai5X5nDiYhoSaxUjOvVjRdRKvDv8yDNB7+Hibo2jqyV5GXJUqhIKf1Y0irNWUd+q6tDJEksrM/Iy5BzeXczo6Z0YNNIVRaWKc0fKWqz41CrewYMHGTduHBMmTLjnmhCCuLg4wsPDjWJcc2jIcQ3vAAplNSdy0xjqHYhCWY3awYLXF3zEv9dG8/kPb7AhahdlJfUOirunFQgozP5NQDMLUNeBtLiWKrcazC0l5GXIyUiRMWika4sXn1rF69SpE/Pnz8fDw6PR4srLly9z69atNlt0uto5MSvor5rjWUF/pUxe3611/noWW5L2c/56FmkVObz+j+UI6R2uJFehqKmg54BO9BnqzJXTMkoKq6n5Ncep60Bir+bmtRqqc+qwERbY2JvRta8N+05d48m+noTSch6oVvGGDBnCkCFDyMrKIi0tjV9++QUXFxemT5/e7G4nY3F3fedq99sX6WrnhK2lDcn5lxjmMxBft64k51+iLuAGeY4lXC/K5MddHxHR5QMS/90Pr9Jfc4+tGrlEiZ3cmuy6Umpt1UgcBb2qPAh7oSPHs4v47Kf6ptOkkX4t9pw613k9evSgR48exrTFYNxd392d+wDC/AZxtuAKrwx9llO/XCBi8Diybv1CcuElhg0fyEteK9kctZAeI2YSFjgdEPxiXQ6Al3CiwFpGrZkK6uCq9R2K0+4w7YkABvq4c6eqhnJ5Lc52LTMtQ2fx2hP3q+8aOHbtHCdy08grKyL/TjFDvfuSnH+Jod6BeDp5kORylUEvL+WXhFgyHx+FmflvX9E12/ruMic7C2TyOjq72JJytQSAtFwpablSbCzNsbGy4KkBXkYXUW/xDh8+zK1bt5g+fToymQwnp9brYXgQv6/vGiiTy1AoqxnqHUhy/qX6cwoZj3ftiVKlYk96PAAd/qSi3/PvoFRUkvnjdgKeehFLG3tNPDJ5HQAeTjaM7teN4X086e5ZRLm8liOXiigqq28TTgo1bhGq9xSu48ePU1paCsDmzZt1CqNUKlmwYAHTp09n2bJltPRARplcxvaUOL5OO8SWpP10dnTHwqzelc8qySP1eiap1zPo3zmA/p49cbTsgJpqzK1skZhbkPJ/C5GX3tDEZ2EG3bs4czG/Pid6uTsQPuRPnM+VUlQmx8vNnqcGeBn9ufQWz9HREaVSCUBFRYVOYeLj4wkICGD37t1YWVmRlJSkb7JNpkwuY+n3n/zazqthmM9AzhRcpk6twlxS//h2ljYA3KiQkne7mDxFGgqLa9jbWNHrmdfxCh7H2c//Tl1tNQB1asgrqX/2ywWlFEgr+fDAeQpuV+HlZs/yyYNbpN7TW7y5c+fi6urKqlWrCA0N1SlMenq6Zs5LaGgoZ8+e1TfZJtEg3IncNIb5DMTW0poTuWm42znT1ckDd3sXAOTKelFuVZVxp6YMe3NXbOv8UCjr23zeQc8Q/OoGLKxsKL+eDUBNnRobSzPScm/z/t6zpFwtIcjfg/WzQvFyd7ivPYZG7zpv27ZtTJs2DTc33QckKysrsbevrzNsbW2Ry+WNru/Zs4c9e/YAEBkZqa9J9+X3wq0YOweAyzdzfm2o9yU5/yIS4O5CXIIE/44dsba3Ir/4t/PWjq7U1Si4+M06PHoEMeCZl1Eo1TjaWlBwu4ogfw8Wjh/QYp4mNCHnzZgxg7179/Kvf/2LW7du6RTGwcFBI5hcLm80MgEwZcoUYmJiiImJoUOH+4+16cuBS8c0wkUOf4EDl44BsGLsHN56cgY9Oj4G3Cucr1tX0m9kkaM4CYDZr8OXdtbmeHq4EvzKeipu5HB6x/vU1cipUNTdV7hyeS37Tl2jXF5rkOe5H3rnvJqaGvr37096ejoTJ07kxIkTWsP07duX06dPM3jwYJKSkggLu9eFNzR3Nxca2n1nC66wYuwcTW+LraU1cZePc11W7+4LBH27+KFWq+nrPJwUmQrPDnaUlFcjr1HhYmeNpZ0jYa9F41V2lv5hvbCwsGB8sO89OS7+fMFvDXcjeZ1657zY2FjkcjlPP/00CQkJOoUZM2YM165dY+rUqVRVVfHEE0/obaiuNHiWUN9Ad7VzYnxgGMN8BnIiN02TAwFsLW1Y/cw8+nsG4GBVPwP8puw2uWVFXCw7iZpqCm/LqalT4+VmT1GZHGc7S1ZOD+V/1iznhbDufPL3F7mQevoeO54a4MXLo3oa1+vUd65gamqq2LBhg4iNjRVKpdKAsxDrac68TSGE2Jy4V/RbN0VsTtwrhBCitKpcfH46VuRIC8XmxL1ic+Jezbl+66aIz0/H3nPf3G/WiH7rpoh5X30m/h1/RXyZkCnySyrE7P9NEE9FHRRLdp/WpHfo0CHh4eEhtm3b1iy7m4Le4i1cuFDIZDKRnp4u1q5da3CDDC3e3SLd/blBrNKq8kafhRD3HDeQX1Ihluw+LfJLKhqdz8jIEIGBgSIvL69ZtuuL3uJFRUVpVuUsX77c0PY0W7yHCfEgUR4kqj7U1dXPot62bZuQyWTNegZd0Vs8mUwmNm3aJJYvXy4OHjxocIOaK15TuFuw3xen+qBSqcSrr74q+vbtK3Jzcw1v6O/Q22E5fPgwxcXFdO3aleHDhxuhFtadBuekYZyuqdzdF6pQVvN6yMT7dmprw8zMjC1btjB79mxCQ0PJyspqll1a09M3wIULF4iKimL8+PGsW7fOGDbpTEMT4G4PsrnxbUnaj62lTaNxQH2QSCQsWLCAb775Bl9fX8rKygxi2/3Qq523du1aioqKOHbsGP7+/q2+v/TDhn5aO77Q0FCEEDz77LOEhIQQHR3d5KXcD0LvtQpFRUWaSbelpaUsWbLEoAYZcq1CW0AqlfL888/j4uLCzp077+ldag56/xSkUiklJSUMGjSIN954w2CGPKq4u7sTHx9Pt27duHr1qkHj1lu8Xbt2kZmZSU5OjsE6kR91rKys2LRpEwMGDGDJkiUGGxLTW7wuXbrg5OTE1KlTjdrN9agSEhJCeHg4O3fubH5k2toSu3fvFiqVSnNcU1Mj8vPzxYoVK8Snn35q8LZLa7TzWpqLFy8KHx8f8e233zYrHq3irVmzRoSHh4tDhw41KyFd+SOIJ4QQUqlU1NXViYsXL4qKigrtAe6D1mJz0aJFfPrpp5w4cYIZM2a02Cj4o46bmxvm5uZs376dYcOGUVBQoH8k+iidk5MjFixYIObMmSOuXr3apF+LNv4oOa8BtVot1q1bJzw9PfUOq1c77/r161y9epW4uDi+//57Ll++rP+vRQvnz5/H2tpa5/tLS0sNNvre2ulZWFjoNQtdq3iRkZHk5OSgUCjw8vLC398fX19fAgICGDBgQHPtbTbPPfecZj/ORzG9h6G1eywiIgI/Pz/NBCITbQet4n344Yd0796dUaNGERQUpNkg3ETro1W8L7/8kp9//pkjR47w2Wef4eLiwogRIwgLC8PR0bElbHwoU6ZMeaTTexh6d0xLpVISEhI4ceIEI0eONNr2/Ca0o3f3mLu7O5MmTUKtVpuEa2W0Fptbtmy555wQguzsbKMYpC8zZ87UfA4PD2fSpElGSUepVLJw4UJKSkrw9/fn/fffb/UNhbTmvC+++ILOnTvTqVMnzatz585tZl9nlUrFjh072LFjh9GEg9ZdLPMgtOa83r17ExISQqdOnRqdP3funNGM0pXCwkKkUikRERE4OTnxz3/+U681FPqQnp6u2XyuYbGMrgttjIVW8R60cc6KFSsMbow29u3bx759+zTHzs7OvPzyy0yePJlDhw6xceNGo9mlbbFMa9CuljVPmjSpUdGoVCpRqeqXYYWFhfH5558bLW1ti2VaA8POiGlhjh8/TnR0NABnzpyhe/fuRkurYbEM1O94379/f6OlpSttbrNUfVCr1SxevJjCwkLs7OxYvXo17u7uRklLqVSyaNEiioqK8PHxITo6utW9zXYt3h+ddl1s/tExideOMYnXjjGJ144xideOMYnXjjGJ145pE+LNnj2boUOH3rOX2cyZM0lPT7/nv/Y++OADdu3ahbFISUkhNja2SWGTkpL473//2+S0MzMzdX62NiHeqlWrePfddxudy8rKonv37vj5+ZGTk6M5f/PmTRITE406HWHHjh2MGzeuSWGTkpI0/zfYFHr27ElqaqpOHd9tomO6c+fO95w7cuQII0eOxMHBAYVCgVqtxszMjE2bNjF37lwsLIxjemZmJu7u7piZmRETE8PRo0eprKzExsaG3r17c+zYMYKDgxk2bBjr169HIpEwduxYzX/4FRcXc/LkyQeGW7hw4UPjXbhwIcHBwcTHx2udqdAmct79OHPmDEFBQQB069aN69evk5uby7Vr1zR/nqjtnzOb8qcdaWlp+Pv7a44dHR3Ztm0bZmZm+Pv78/XXX3P06FESEhKYNWsWX331lWbP0Tt37uDi4vLQcNriBfDz89NpvLRNildSUoKLiwuWlpYA+Pv7k5OTw8aNG3nrrbeA+tnb165dY/ny5ezfv5+qqirWrl3L0qVLiY+Pp6CgAE9PT2QyGVFRUURGRpKamsr27du5ceMGkZGRXLhwgQ0bNlBb+9v+YOXl5Y1mRDfMYHZwcOCxxx7D3Nwcc3NzXn31Vc6ePUtERAS3b9fv6J6UlKTZ3fBB4bTFC/XrGMrLy7V+T22i2Pw9CQkJjXaa8Pf3JzY2lpqaGk1uzMzMZPr06QwfPpxly5Yhk8kwNzenY8eOlJSUIJFI6NWrF19++SVz5szBysqKrVu34u3tTX5+PtbW1pSUlGBlZdVoSoeLiwuVlZVabfz222+JiIjgscce47nnnmPatGmkpKQQGRlJfHx8s56/oqICV1ftW/u3CfGWLFlCWloatbW1XLp0CTMzM1atWqW57u/vz7Jlyzhw4IDmXGZmJuHh4RQXF2uK1bvXx3/88cdERERw8OBBPDw8SE9Px8/PDxsbG/bu3csLL7zA1q1bWbBgQSNbBg0axO7du7XOh+nVqxfz58/HycmJIUOG4ODgQGVlpUEGaTMyMhg4cKD2Gw212sWQ7N+/X+s9K1asEKtXrxZLly4VFRUVIi4uTixevFisWrVKlJWViZUrVwohhIiLixPLli0Ta9asETU1NeLkyZPivffeEzdv3hTTpk27b9yvv/66Zkej1uDNN98UlZWVWu8zjefdh6SkJIqLi1tlXmpmZibJycka7/VhmMRrx7RJb9OEbpjEa8eYxGvHmMRrx5jEa8eYxGvHmMRrx5jEa8f8PwpbVPi4tKPdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 80.64x69.12 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(1.4*0.8,1.2*0.8))\n",
    "\n",
    "for i in range(nc0):\n",
    "    ax.scatter(all_features[idx_sub0[label_sub0==i],9],all_features[idx_sub0[label_sub0==i],10], color=fs_pal[i],s=1,label='RS')\n",
    "plt.xlabel('$1/V_{below}$ (ms/mm)')\n",
    "plt.ylabel('$1/V_{above}$ (ms/mm)')\n",
    "ef.simpleaxis(ax)\n",
    "plt.plot([-10,10],[10,-10],'--',color='black',linewidth=1)\n",
    "plt.xlim([-7,4])\n",
    "plt.ylim([-4,7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-4.0, 7.0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAG8AAABuCAYAAAApmU3FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYTklEQVR4nO2deViTV973PyHssllQXCo+SHBcsGofN5y6dZy3juNDq9YqWIsdfTrTam0dpdjrFbW1LsXq41WdztjXtSoufbQVrAt1Cq6xuFCqCG5YWVQkIGuAhHDeP24TFsEkkCi0+VxXrjs5uc85v+R7n3N+Z7tvmRBCYKNVYve0DbDRdGzitWJs4rVibOK1YmzitWJs4rVibOK1YmzitWLsTTmpoqKCo0ePkpycTFFREW3btqV///6MHj0aFxcXa9tooxFkxkZYdu7cycWLFxk8eDABAQF4e3tTUlJCWloaZ8+epV+/frzxxhtPyl4btTAqXnp6Oj169Gjy9zash1HxbLRcTGrzANasWUNCQgL29vYIIZDJZHzzzTfWtM2GEUwW7+zZs8TGxiKTyaxpjw0zMFm84cOH8/333+Pn52cIs7V1TxeTxUtLS+Pq1au4ubkZwlasWGEVo2yYhsniVVRUsGnTJmvaYsNMTPY258yZg6enJwqFwtDuWaN/d+jQIWbNmkVcXBxBQUEWT//XhMklb9SoUU/EWfH39+eTTz5Bp9MZvFobDWOyeBkZGSQmJiKXywGQyWS88sorVjFq6tSpAERERNCxY0fmzp3bskVUqyB1C/R+E1x9jIdbCmEikyZNEtXV1aae3mSuXLlieH/79m3Rt29f8Ze//EVUVlZaPe86lOUJkRQtHY1xarEQnyEda5MULYUnRVvDQmHyrIK+q5Cenm54WRs/Pz9OnTpFfn4+mzdvtnp+dUheDyc+kI71Uavg3CrpCKCvFOpXDr3fhOHR0tEKtPiugpubG/v37wfgzJkzeHp60rt3b6vnaxAiKxHOLIF+s6WqT62CI+Fw65D0/cAI6TuHNjUiqVWS6DJq4lmBJ9JV0Gg0jBo1im7dugEQGRlplidpZydVEJmZmcyZM4etW7cyduzYJtnyWPRtVLcQEEDn4ZBzXHo5tJGESl4vCdd5BGjKpDiuPtJ3elK3wNmPpPf6eFbAZPHatGlDVFRUk7oK169fZ8yYMURFRTXNyodMmTIFPz8/Jk6cyPLly3nzzWZWR/UditQtUlWZlSgJ5DdaOs9vtHSOWgV3TkthZTmSQI4PS1ztktb7TUlYGVarMsHMrkJTuXr1Kj///DNTp04lKCiIyMhIQ2kyl6FDh3L27FnKy8upqKjAzs4OR0fHphmmF0tbJpWQbiFSeLcQ6DJSOmbESgKU58O346DwBngppKNTW7j2LRTdhpR/SHH1Je33S5pmkzkY82gOHDgg7t+/3yyvSKlUilOnTgkhhPj4449FbGxsne93794txo8fL8aPH284zxRWr14tRo4cKVQqVdMM03uUem9R7xU25GnuHS2d8xlCJEYIsalHzed1baXj3tF145jqsZrj2dbCqHhnzpwRa9asEZGRkWLlypUiISFBlJWVmZWJWq0WWq1WCCFEYmKiWL16daPn1u4qGKOqqkpERESIgIAAs+I9Qv0/T+/iJ0QIsW+sEKp0SbDPEGLXcOk8VboQX/rXCLhvrBDZyprza6fzuIuiLE+K04QuhdFqMzg4mODgYACKiopQKpUsX74crVbLp59+alLpXrt2LX369GHcuHGcP3/eYrMRcrmc6OhoevXqRVxcHD179rRIunQLkdq9e0mSs6JKA08/GBABigmSt9nvXXDvAm06QKcXYNAHNV5owTXoORV6hNakd26VVD0rHzoyeicmdYsUx3+s2e2jyW3etm3b8PT0RKVSceXKFQYPHmxyJjNmzCAiIoI9e/bg7+/PmDFjzDLSGNOnTwfg2LFjpKWlMXv2bPNGZPRtH0h/4PG/S3+oZzdo4wclt6SXnQPc+EZq77ISoUotxQkcL7WJmlLw6AZFNyRnJvccjNlWk/6QxY/2+/TvmzAKY7J4hw8fZteuXYSHh7N//36zBqXbt2/Ptm3bzDKsKQQEBPD++++TmprKunXrcHBwMC1i7T/wp4ddAae2UJQhCQjg3lX6XJQhfa5Sg4M7+PQB7z6wZxiU50nxOgSDvaOUjt6b1adfX6D63QwzMNnlE0Jw5MgRAgMDycrKori4uEkZWhN/f3/OnDlDdnY2s2fPfuR7dVE+SXvXoS7KbzwRzcPS5BkgHV19oe8sqCiShHPxxdCD15bA3TNwcLIkHHKofAD3lFIp/c8IqcsAkkAW7qybXPLmzZtHQkICs2bN4tChQyxcuNCihlgKDw8PDhw4gEqlorS0lLt37xIYGAjA5aMxnNgotTmDXnu3JlLtLkNeshTWcTBoCuGuEh5cA22hFF6e+2im2ocXsr0LuHgDAjKPgfxh6XO0Tke9xa0eS0tLs5jjcezYMcLCwoh6axIzIj4GJAEDgsdwU3mEoJfCKC8uIPGLBQwZ0IVf0q9CdiL9XxiGq//QGucCADmgq5uBkzdUNlKK/UbD0KVwJgra94eBHzy9ktcaV4+NHj2ape+EMn/FOlLOn+XzPfEMeu1dkvauM5TArJ9Pc+vCce7d8Ka8KB9oQ45TNeMqDuPq4Aba0oep6R7NwCCcDGk8rRZtfycJl3lMerm2s3jp+9WvHnt1xmwepB7n4pWrHF41izHz1xP0UhgAQS+FERA8hsK7v/Ag6wbOHs9QUVxA5pWf2HxNx/jny+nc1pRcagn3UPDqO+ewy0sCQCMLQt5pAnIL/7YWPSVkKg05Ivqw9IT92Bdl89KQIH6IP8zokcOorLZj0Gvv4urpDYCrVzs8O/hRUVyA3NEZgIoqOfsvuqPWmHmxPiypVXnFFJzsQnHWYHLWtiF7wkI0N25b5gc/pMVPCZlC8oGNKHesQluh5vdvRAI1zknw6xEMn7mYgOAx9Dp9CNmJNIYMGcLemK8oTU0k7Yd9FN6R3H+Z3B6dpsKQbmWVnK/PuaPw1dDfrxJXx8bdA02JnLJrzyITlTh01FFy+zlkvxtMeXoG1eUpVGdkk7dgLZ3/938s9rtb/OoxdVE+l4/GEPRSmKGkmELtqlEfL3jKewRPgS+//JL/Dp/Ga50KAHD2aEtF8QOErgpkdiCqDenklTqQV+pAQamMPz9Xjp0dVKntENUgd6mmukKOJt+ZvMO/o+qBay0L7iBL+R5RVIq8a0d0t+/iFBTY/D+kFk9kSqg5NOre16L/yzNxcHYlIHgMSXvXGQTTn1//Anh98kTcrx/mxvnj3ChzZMGHn3L8/y2hVHUHZLJHfA+Aa7nO9DnUgXY9ClH927+eUA0jikpxGT0En6WzKd3/bwQCXX4hcm+vJv8ftXkiU0LNoXYJ0lNfDL1Qtb3I2kLrLwBthdQBz0lNIiflJLh4ceRcJtnvvM3/6VyF3E4G1Q14lYDAjrT7zugumz6LL+/UDp+ls3FUdAWgcNVWKpMu47thMXJvL3T5hRTvOoRH6NgmCWpUvIMHDzJu3DjGjx//yHdCCOLi4ggJCTE7Y1OpXYL0NFQa1UX5aCvUBL8eQUDwGE5/JQ2a9395JkEvhaGtUPPLhUTupp0DwK//CNzbdeb/ep1iS1IuW1MeMP/lgZSpctCoS6mPS6mcrvmephntaA+aKnR38ig7chrH2V0RD4tz+fHzFG3aj8zVmWp1OYWrtgLQdnbYYxJsGKPi+fr6MmfOHNq1a1dnc2Vqair3799/Khsr9aWwdjV5+WgMyh2rGD5zMTeVR1DuWAWAg7OrQWC9cJ4d/NBUqEmNjwFgeh9P0rNc0FWUUVRUTBsXZ6qrNHXyLHfR4aQzsaLSVCFzdsJ1wmg8QqXlGl4zJgIgQ4ZA8OCjf9I24k28IqYj1BVNqk6NWjNw4EAGDhzI1atXSU5O5pdffsHLy4uwsDDDsNOTpqFqsn71qq8ia1e3AJ6d/AkcOpbz//sPQ5i2rIjnewWQfe8+/3OmkLA+7ii8neo4Lub2bkVFJepv/o3u3VDk3l7Ivb3w/mAGALr8QuxcXfAIHUvxrkMUfPRPZK7OZpe+Vj08Zoonmp91ncQNUQwJm0fO5bMGMQ98/CY5l87gG9gXJzcvKksLyb2ewvWCKmIuFfNSgCtDnnU2pONYKePPV8yfh3Qe/p+4jhrUaLvWnHavVd8NQl8C6wtXu9OeuCGKW0nHOBuz2lC9Ary8aAvDZy6mS78XyEw+TtmD+wAEPmPPOwM8uVVUja665rr2KW5jvoFyOfZdOlDw0T8p3nWo4VO8vWg7O8w6Dkt9jh07xv379wkLC6O4uBgPDw+zM7U2tR2akX9dajjqw3+5kEjn3oPo//JMkg9sBMC9XWdKVXfx7d4Xr9Ji2rW5RYWQ8+2VIsYGuqJtyk0vdDpKYxPwnBVqaPssidninThxgvbt2wPwxRdfsGDBAqNxtFot8+fPJy8vD4VCwUcffWTVMdL6HfSJn+wGwMXjGWkgOukYmcnHDecPeHUWWk05cnsHBDJyr6Vg7+SCQ7kauZ2MLy6UsvE/nsOhtz/IZGjv5kJ+I/OZLk5QoQF9a1SiRnP1lsX6drUxu9p0d3dHq9UCUFJSYlKc+Ph4AgMDiYmJwdHREaVSaW62ZtFQdapvH0f+dSkDXp2FX/8RqIvyUe5Yxf2bl0mJ3Uz2JSU5l87g4ulNVWU5cjsZ04b/jv9q78v09NMUlBfil7CZZ2ZPbTzz8soa4dq44jiwNz5LpYlhXX4hD9bHoMsvtMjvNLvkzZo1i71797Js2TKGDh1qUpyUlBT++Mc/AtK6y/Pnz5sc11LU7qjfv3mZzOTjhu5A+4Ag2gcEce/aRYSAnEtKOvcZitzegfYBQZTc/wdvd3anwN+BH48lMCB0LNXqcoS6ApmrM64vDqZgxUbsu/iivZVDdbEa7eXrOPbyp9P2lYZSp/csoWn9uvqYLd7mzZsJDQ3F29v0ccbS0lLatJEafBcXF9RqdZ3v9+zZw549ewBpxt4a1O5C6KtMXZWWzn2GggwcnFzJ/llJ8OsR+PV9gayUU2QmH6e9Iojg1yOoSLqE+95rhB4K4+VpYURHRxu2uwF03rfW8P7Oq3PRAppzqRRt2s8zH/wFwNDuWar9M1u8qVOnsnv3brRaLVOmTDG0f4/Dzc3NIJhara4zMwEwefJkJk+eDEhdBXMwdeBaX5Xqp41yUn8kM/mE9P7SGcPsg94jzb4kVe0OTq78/o1IdH8upLjHIZRjgpny15mEhISwe/du3N3dH8nLZ+X73AmLRHcrxzCyAjWepaUwW7zKykr69u1LSkoKEydO5OTJk0bj9OnThx9//JEBAwagVCoZMWJEk4xtCFMGrmvj6unN79+IRF2UT/KBjYbOPNQ4OPrhNJCG16DuH3/kyBE+//xzHB0dG9y966joSpfD/zL036yGuYuLN2zYIL7//ntx8+ZNwypoY2g0GjF37lwxefJksWDBgsdu0jR35XNZoUr8uOdzUVbYxCXvQogf93wuVv3RW/y453Oz0q6srBRDhgwRJ06caHLezcFs8S5cuCDWrFkjYmNjTRbPHJq1bN0EGhKkMZEaE7U2R48eFe3atRObN2+2ms2NYbZ48+fPF8XFxSIlJUVER1t+u661xTOnlJlaqtPS0kRQUJC4ffu2VWxuDLP7eR4eHmg0Gp577jnKysqsUZNblaCXwgyOSW30bad++AwwtH+Xj8Y8dqFujx49+Omnn/Dz82PLli0m93+bjblqFxcXi/Xr14vFixeLgwcPWvxqsnbJa4zmVJ16dDqdeOutt0SfPn3ErVu3rGRpDWaLt3//fhEVFSW+/PJLUVpaanGDLC1eUx0afTxV5jWz4ldXV4u1a9eKjh07ivT09KaYbDJmi7dkyRIhhBC5ubli8eLFlrbH4uKZU3IsEU/P6dOnhUajEQUFBU2Kbwpm9fOio6O5c+cOx48fR6FQtIr7Sze0Bsaa8fQMHToUIQSvvPIKwcHBLF++vMlbuRvD7MnYO3fuGBbdFhQUWHzDiSX3KrQEVCoVr776Kl5eXuzYseOR0aXmYPaloFKpyMvL4/nnn+edd96xmCG/Vnx8fIiPj+fZZ5/lxo0bFk3bbPF27txJeno6GRkZVhtE/rXh6OjI+vXr6devHwsXLrTYlJjZ4nXs2BEPDw+mTJnCCy+8YBEjfksEBwcTEhLCjh07mp+YMY8mJiZG6HQ6w+fKykqRmZkpli5dKjZs2GBxD+pp9fOeJJcuXRL+/v7iu+++a1Y6RsVbuXKlCAkJEUePHm1WRqbyWxBPCCFUKpWoqqoSly5dEiUlJU1Kw2i1GRkZyYYNGzh58iRTp07l/PnzzS/uNvD29kYul7N161aGDRtGVlaW+YmYo3RGRoZ47733xNtvvy1u3LjRpKvFGL+VkqenurparFq1SnTq1MnsuGb183Jycrhx4wZxcXEcPnyY1NRU868WI/z00084OTmZfH5BQQHPPPOMxe14GvnZ29ubtQrdqHjz5s0jIyOD8vJyunTpgkKhoFu3bgQGBtKvX7/m2ttsJkyYYLgf568xv8dhdHgsPDycgIAAwwIiGy0Ho+J99tlndO/enT/84Q8MGjSozoopG08Xo+J99dVXXLt2jR9++IGNGzfi5eXFqFGjGDFiRIMrp540+lVnv9b8HofZA9MqlYqEhAROnjzJiy++aLXb89swjtnDYz4+PkyaNInq6mqbcE8Zo9Xmv/71r0fChBBcv37dKgaZy7Rp0wzvQ0JCmDRpklXyedKbZUzBaMnbtm0bHTp0wNfX1/Dq0KFD0+/rbGF0Oh3bt29n+/btVhMOnvxmGVMwWvJ69epFcHAwvr6+dcIvXrxoNaNMJTs7G5VKRXh4OB4eHixZssSsPRTm0BI2y9THqHiN3Thn6dKlFjfGGF9//TVff/214bOnpyczZ87ktdde4+jRo6xdu9ZqdhnbLPM0MHuvwtNk0qRJdapGrVaLTifdN2XEiBFs2bLFankb2yzzNGjVe9JPnDjB8uXLATh37hzdu3e3Wl76zTIASqWSvn37Wi0vU2lxd4Mwh+rqaj788EOys7NxdXVlxYoV+PhY57k9Wq2WyMhI7ty5g7+/P8uXL3/q3marFu+3TquuNn/r2MRrxdjEa8XYxGvF2MRrxdjEa8XYxGvFtAjxZsyYwZAhQ/jiiy/qhE+bNo2UlBT+/ve/1wn/9NNP2blzp9XsSUpKIjY2tklxlUol3377bZPzTk9PN/m3tQjxli1bxgcffFAn7OrVq3Tv3p2AgAAyMjIM4ffu3ePUqVNWXY6wfft2xo0b16S4SqXS8LzBptCjRw8uXLhg0sB3ixiY7tChwyNhP/zwAy+++CJubm6Ul5dTXV2NnZ0d69evZ9asWdjbW8f09PR0fHx8sLOzY//+/SQmJlJaWoqzszO9evXi+PHjDB48mGHDhrF69WpkMhl/+tOfDM/wy83N5fTp043Gmz9//mPTnT9/PoMHDyY+Pt7oSoUWUfIa4ty5cwwaNAiAZ599lpycHG7dusXNmzcND0809uTMpjy0Izk5GYVCYfjs7u7O5s2bsbOzQ6FQsHv3bhITE0lISGD69Ons2rXLcM/RwsJCvLy8HhvPWLogPQfQlPnSFileXl4eXl5ehocXKhQKMjIyWLt2LXPnzgWk1ds3b95k8eLF7Nu3j7KyMqKjo4mKiiI+Pp6srCw6depEcXExH3/8MfPmzePChQts3bqVu3fvMm/ePH7++WfWrFmDRlNzM/CioqI6K6L1K5jd3Nzo2rUrcrkcuVzOW2+9xfnz5wkPDyc/X7rNh1KpNDzRs7F4xtIFaR9DUVGR0f+pRVSb9UlISGDkyJGGzwqFgtjYWCorKw2lMT09nbCwMEaOHMmiRYsoLi5GLpfTvn178vLykMlk9OzZk6+++oq3334bR0dHNm3ahJ+fH5mZmTg5OZGXl4ejo2OdJR1eXl6Ulj56a/76fPfdd4SHh9O1a1cmTJhAaGgoSUlJzJs3j/j4+Gb9/pKSEtq2Nf4EqhYh3sKFC0lOTkaj0XD58mXs7OxYtmyZ4XuFQsGiRYs4cOCAISw9PZ2QkBByc3MN1Wrt/fHr1q0jPDycgwcP0q5dO1JSUggICMDZ2Zm9e/fy+uuvs2nTJt577706tjz//PPExMQYXQ/Ts2dP5syZg4eHBwMHDsTNzY3S0lKLTNKmpaXRv39/4ydaareLJdm3b5/Rc5YuXSpWrFghoqKiRElJiYiLixMffvihWLZsmXjw4IH45JNPhBBCxMXFiUWLFomVK1eKyspKcfr0abFgwQJx7949ERoa2mDaf/vb30RVVZVFf5M5vP/++ybd48Y2n9cASqWS3Nzcp7IuNT09nbNnzxq818dhE68V0yK9TRumYROvFWMTrxVjE68VYxOvFWMTrxVjE68VYxOvFfP/AVzfhHPVC9HDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 80.64x69.12 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(1.4*0.8,1.2*0.8))\n",
    "\n",
    "for i in range(nc1):\n",
    "    ax.scatter(all_features[idx_sub1[label_sub1==i],9],all_features[idx_sub1[label_sub1==i],10], color=rs_pal[i],s=1,label='RS')\n",
    "plt.xlabel('$1/V_{below}$ (ms/mm)')\n",
    "plt.ylabel('$1/V_{above}$ (ms/mm)')\n",
    "ef.simpleaxis(ax)\n",
    "plt.plot([-10,10],[10,-10],'--',color='black',linewidth=1)\n",
    "plt.xlim([-7,4])\n",
    "plt.ylim([-4,7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "nf=features[:,:,3]\n",
    "z=nf.copy()\n",
    "z[~np.isnan(nf)]=stats.zscore(nf[~np.isnan(nf)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAABkCAYAAADOmA9MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmW0lEQVR4nO2de3RUVZbwf/dWpapSlUpV5U14BhJ5JEAMAhEVH+jSXtouW0VQB/Eb12i3M04zdjvpbmW0p7vBR7f9ctas7m9aGxn9ZDkGpaWhpRVFhYhImkAekAeQF+SdVFKVet7z/VHWJUUiKUKqKnHq9w/WSXnv3nXvPvucffbZRxJCCOLEiXNRyLEWIE6cyUjccOLEGQNxw4kTZwzEDSdOnDEQN5w4ccZA3HDixBkD2lgLsHnzZpYuXcq1117L97//fTo6OsjNzeXHP/4xPp9vWJskSbEWOU6c2Hkcr9fL448/zp49ewB47733yMvL4/XXX0en03HgwIER2+LEmQjE1HDuvvtuvvWtbwFw5MgRli9fDsCKFSs4dOjQiG1x4kwEYmY4RqORFStWqJ8HBgYwmUwAJCYm4nQ6R2yLE2ciENYcx+Vy8Ze//IXy8nL6+vqw2Wxcfvnl3HjjjSQmJo6LIElJSaphOJ1OkpKSUBRlWNv5bNu2jW3btgGwYcMGVq5cOeq9Dh48SEtLC16vF4CdjQaEgHaXjCxBukEB4NYZLkwmE0ajkdTUVKZNm4bNZovoPOtC+lRXV2MymZgxYwYAe/fuxel04nA42Nq2D4Cznl68wsd0fRoAErA++3qSkpIwGo309vYydepUioqKIqZDOLqMFb/fzwcffEBPTw89lXkAKH4Jv8OEpPGTnC7h8/rweD1kXn4ajUaDz+dDkiRMJhPXXnstOp3ukuUY1XBee+01Dh8+zPLly/nmN79Jamoq/f39VFdX8+STT1JYWMgDDzxwyYIsXLiQzz77jCuuuIIDBw5w7bXX0tvbO6ztfNasWcOaNWuAwIsVDn6/H7/fjyRJCCG4dYaLVofMruZErkhzU5DiU79rNpspLi4elx87HC6kT0JCAj7fOdksFgvd3d0ArMsMvJTPN71Dls7KusyV6PV6FCXQCaSkpNDc3Ex6ejppaWnRUGVMz2Y0XC4Xfr8fAFt+LQA9Vbn4gZTFNWRNs9DR0YEJSEnJoqurC5vNhl6vZ+nSpWi14xMPG/UqS5Ys4f777x/WvnDhQu655x5qamrGRZBbbrmFkpIS1q5dS05ODldffTU+n29Y23jgdruRJEl9qQAquhMAwWXWcy9mQkICaWlpUTOa0dBqtXg8HvWzJEnIsqy+SHafE6/wM0ufAYDJZKKnpwebzYbD4UCSJCRJwmq1xkL8S8Lr9dLW1kZ9fT12u/1cu1OPt8+MztaHRu+hr68PSZLQaDQYjUYcDgdarZYlS5aMm9FAGIYzb948AF588UX27t2LVqtFCIEkSWzfvl39+1h57LHH1P9+8cUXQ/6WkJAwrG08GBwcDDEaRcAZp4Y0vYJuyKxPp9Oh0WjG/f5jRavVhszzvF4vQ5Pb6wfbAJhnzAYC88K+vj6EEPT09JCWlobH48FsNkdX8DGiKArd3d20traqnqO3tzfkO/31MwEw5zSi1WoxGAx4vV5SU1PVd3XJkiXj3vmFbYJlZWXs2LFjUq+jNDc309fXx8DAgPrC7Ww0MOgDBQmvEvh86wwXAAaDYVx7qUtFq9Wq8zK/3093dzc+n0+d37R7+gD4a89RJElivXwder0ej8eDoijodDqSkpImxTNsbm6mvr6exMREpkyZwvz589m3bx9ut/vc3MYn43eakLRe7HU5gMT8VXYcDgdCCE6dOkVBQcG4zcOHEvZbsXLlSvbs2aNOTIFL9jax4MSJE5y/BcmrBF4k/XnOxefzkZmZGS3RLgpZlhkcHAxp86MgIamGIYTA7/fj8XgQQtDf38+SJUtiIe5F09/fz5QpU7jsssuAgC7B+VwQxa0DBBqD68uWgI4QiNL6/X4yMjIiIl/YhlNdXc3x48dDIlubN2+OiFCRYsqUKeh0OrXXhkDk7KMzOursMrfOcCEP6Yzz8/NHjOTFCr/frw4dJUkiMTERr9erBgZ+17oHh+JmXeZKtFotfr9fHc5ptVoKCwux2WyxVCFssrKyOH78uPp5cHBQ7fBs+bUIv0TH54VokxykLAwECQLDNEUNiuh0OnU5Y7wJ23BcLhd/+MMfIiJEtNBoNMO8DYDLLyEhQowG4MyZM0ybNi1K0o2Ooighc66h8zQIeByNFJikBSfIsizj8XjIyMiIWjRtPLBarbjdbhwOByaTiYqKipC/D7algZCwzOpQ23w+H2azmf7+fmbPnk1LS0vE5qhhG47JZGLjxo3k5uaqQ4HxCENHk4GBgWHDGwC3X0IzwrC/paUFvV5PQUEBshz7fNihHgcI8ZwAfqGg+XJNWwihDtUgMPRRFGVC6BEOkiSRmZlJW1sbU6ZMoampKeTvzrMZSBofGbM0dHUF2tLS0ujo6CAnJ4eEhISIzk/DvvL1118fMSGiQXt7O0ePHg3xODsbDQB0u2WEOPc5GBy4/PLLqaysxO12s3jx4piHpX0+n2o4drtdNYpgcGDA70JGUj//n6k3IMsyGo0Gl8vF8ePHmT9/fmyEHwOZmZnU1NSo89KhC56KW4+c4OHkJ2n4FRu2/Fo6OzuRJImCggJqamoiGhEN23BycnKYN28eBkPg5Rq6njAZOHv2rNrjnj/EUQTDhmkQGNpptVp6e3sZGBggJSUlStKOTFJSEnV1dWRlZdHW1jbs7wKBLJ17WYKGFuwszg/lTnSSkpLo7e1VJ/wAQoDfaQQEst6DXznXESqKQmZmJh0dHZw+fZrs7OyIyRa24bzzzjs8++yzeL1eZs+ezfHjx9mxY0fEBBtvCgoKaG9vD1lEvHWGC0XAKyeMTDP5uXGqW/1+amoqFRUVJCUlUVhYSHJycqxEV8nKysLv93Po0CG0Wq06ZF6XuZIBv4tft/yZwqRZ3GRbpGZFqKvsNpuaMDtZ6OnpCTEaW34tvTWzQcgk5ZzGmNWl/s1kMqEoCpIkUV9fj8ViiWhaUdiG8/TTTwPncoX0en3EhIoEsixjNptRFCXEWw54JUAi3eAP+X5/fz85OTkUFBRMqEXQqVOnotFo2L9/f0j6TbMr8BLlGrIA1EVqIQQ6nY7rrrtu0sxvIOA9Dh48GNLmaMnA02NFn9odYjRpaWl4vV6cTie5ubnU1dUxf/78iK5XXfQvqdFouOmmmyad24fAZPr8cGyHK/ATZCWGDt8WL17M4sWLJ5TRBElISBjWcTW4AkO3Kfpz+gWNZ9GiRSQkJERVxkthYGCA3bt3h2RJeOxGHI1T0RhcJOedAgKdodFoxO12oygKSUlJeDweNSk3koTtcb71rW+Rnp7O/PnzSU1NnZAv1Gh4PJ5hcf2zzoAeNn2o4djt9gkXhXK5XDQ1NVFfXx/iUQBaPD3oJS0GOWAgwYRQg8FATk5OLMUOG7vdTn19PQ0NDeoQE8Dv0dBXnQeygjX/BEFHEuxA0tLSMBgM9PT0cOLEiagMScN+K7Zv386zzz7LsmXL8Hg8aLVa7rrrLtauXRtJ+cYNRVEQQoR4yp2NBurtWkCwp8UQ8v3W1lZaWlqiK+QF8Pl87N+/Xw3PJiQkqEaztW0fnV47fiHUiFpwmJKVlTUpUmw6OjooLy+npaVFNZqeyjx6KvPo/tsChCKjMbiw1+ao0bXs7Gw19Qigu7sbk8kUlVy8UT3Ou+++y2233QYEUtOvuuoqrrrqKiDwMk6WAIEsy8yfP5/Kysqwvn/ZZZcxderUCEsVPlqtliuuuIK6ujp6e3uHRQhlJBTOeU2v14ssyzGPBIbL4OAgmZmZZGZmhix2CgHCr0WS/cjac/rp9Xp6enpwuVz4fD7q6+uZOXMmCxYsiIq8oxpOZmYm//zP/0x6ejpz5sxR9+NUVVXR1tY2qRZBs7OzcTgc9PT0AIGo2s5GPR0ujbp2A4F53Jw5c2Il5leSnJxMUVER3d3dVFVV4XA4gEBUrbTjM6oHW7gv4+qQ7w+NSk1kgh3BnDlzVMOx5dcy2G6jv9+Mec5pDGl96vd9Pg12ux2LxYLVamXBggURS68ZiVENZ+nSpSxdupTjx49TXl7OqVOnsFqt3HvvveTl5UVDxnEl+GCCwwFJgvOTcIL5bBN1Qh30/O+//746TJlpSKd6sIUOj50svRWDwUBiYiK9vb3qfGgi4/f7kWWZgYGBkHbnmUyQFPSp54xGp9OhKAqJiYnk5+czZcqUqOsXdnBg7ty5zJ07N5KyRAVZltUESPhykjfEcjQaDQkJCbhcrglrOEIIGhsbGRwcJDs7m9bWVmYYAnlop90dZOmteL1eNYDjcDgmVLLqSATz8IYmdio+DX5nIjpbnxoQCD4/l8vFddddh9FojIm8EydkFCX8fn9I2k2HS0bhXLpNcMjgdru/4gqx5+jRo5w9exabzaam0OzqKgfg077Ai+fz+bDb7SQlJdHX1/eV15ooKIqCw+GgoyOQtNlTmUd3xVxAwu8+14EF90glJibGNAUqbI/jcDjYunUrXV1drFixgry8vAmVORwusixjMpnURVAhhrt4vV6PxWKJtmhhY7FY6O3tJS0tjYaGBnQ6XWBbNIQECIxGI3a7Xd3TMpHJyMigvLw8JAwtvIHt7JImoJNer8fn8yGEIC8vL6ZLImF7nJKSEqZPn05lZSWpqamUlJREUq6IIUkSS5cuVT8nagWJGqEGBywWC8XFxRN2mAYwc+ZM5s6dS1tbG16vl5SUFNZlrkQjyWTpAgugwZoKRUVFMRvOXAw2m43Zs2erC7vWBbUIJLQmJykFgf02wXWpK6+8knnz5sV03ha24djtdm699Va0Wi2LFi2K2sKg1+vlu9/9Lvfddx//9m//NuJ+movFarWq3tKjSOg0gWsaDAauueaaSZFOlJmZSVFRER6PB5fLNezvycnJXH755RN+bjOUmTNnqnuGFLcOFBmdNVCYQ5ZlsrOzWbVq1YTYVxT225+WlsZrr72Gw+Fgx44dUdtJGKkyuEGv41XA8KXHX7VqVUT2p0eKoHc836NIksTy5csnRGLqxRBMD0pMTMTdHRgq61MDSwdLlizhyiuvnDAjgbANZ/PmzSiKwqJFi+ju7ubZZ5+NpFwqkSqDGyz95Bdg1ApuuOGGqK4DjBcGg4Hi4uKQ2gjXXHPNpCwBBYHnctVVV+HusYKkoDW6+MY3vjHh0obCNhxZlrnqqqt4+umnmTdvXtSGM6OVwd22bRt33nknd95557BiDqORkpoGSGgkMWH24o9FH61WS1FRETpJS49vgNra2ghLGR5jfTZGoxHFq0XSKMyenTMhh85hG86//Mu/UF4eCHnW1NTwxBNPREyooYxUGncoa9asobS0lNLS0otKL3E4HJw4XoMsCRxeiV27do3L/OlSGYs+LpeL8vJy5hmnYvcPUtVYyzvvvBPzWttj0cXv9/P++++jMbgRPg0NDSf55JNPOHXqVEjELdaEbTidnZ3cddddADz44IN0dnZGTKihBEvjAhw4cIDFixdf8jWFELz//vsA6GRw+mWcTqd6n8mEz+ejvLwcm81GsTmQyVHWX4vb7Wbnzp1UVVVNiA4hHIQQHDp0KFB90+QEJBRPAr29vfT09LB//37OnDkzIfQJ23CsViuvv/46FRUVbNu2LWrVIG+55Rbq6+tZu3YtDodjXMrgfvHFF2okyqBRcPkCYc3Gxkaqqqou+frRQlEUjhw5gtVqRQiBLSGJZE0iVc5mtVjHsWPH2LNnz4hFSiYSQgiqq6s5c+YMAAlJgTw8nzMRn89HZ2cn+fn5NDU18dlnn9HV1XWhy0WcsA3n5z//OX6/n3feeYfBwUFeeOGFSMqlEiyD+8Ybb7B58+ZLjt23trZy8uRJ9bNRK/AO2YpTWVk5rKLKREQIQWVlJRqNhnnz5tHe3g5AoWkWLsXLWU8gGqXRBJIhd+3aRV1d3YTorUeisbGRpqYmdTimNQY6Nm9/YH7b09NDbW0tS5cuJScnh+rqag4fPjwsty1ahG04g4ODJCQkMGPGDCRJ4n/+538iKVfEOHLkCEIIdR2qxy2jIPHu6UDKjSzLHD58OJYihkVvby9tbW3MmjULSZJUD1r/5U7Q/9exn61t+/D7/UydOhVJkjhy5Ah1dXWxFHtEHA4HDQ0NwLlacX21swDBYFuauv+msbERu91OZmYmK1aswGKxUFZWFpPCMWEbzj/+4z9it9tJTk4mOTl50hTuPh9ZltXs2q9iMqx/BFPp//a3v9HQ0KCGozWSjBYZl+IJqW5jMBjIzs7m5MmT6naEiUJPT4+61TnYoakDiyEpUTqdTtVJkiQcDgfTpk2LSc5a2LlqZrOZhx9+OJKyRAWtVovFYlGTCaea/Gr52yBZWVmxEi9sJEkiOzublJQUKisr1QjausyVHLCf4IPeY9xgLQACo4XU1FSmTp1KU1MTNTU1FBUVTZitBsHDys6ePatmrtvya2kvKyTB7MA6vx4IVB4KJt+2tLTgdDrJz8+Picxhexyz2czGjRvZsmULr776Kq+++mok5YoYQ7cUQKCmGgi1h1MUJeKFHsYTg8FAUVERs2fPVtsKk2YB8Fl/YFjm8/lCjvnweDwj1mWLFX19fVgslmHPBkL3StntdgYHBxkYGKCuro6FCxfGLNEzbI9z/mloE6W3ulgkSQqpO6Cclx09USfPF0KSJKZPn86xY8fo7+8nUdaRqk2i3nVW3cTm8Xg4evSoemzGiRMnSE1NjXkKi9frxeVyqetzw4bQQx6H0+mkoaGBlpYW8vLyYprpEbbHWbZsGQ6Hg/7+frUayWTD4/HQ09MT8nDOOAM/QXA/DhASdZtMBOdmW9v24RcKXuHn/555Xy3gcfr0afr7+3G5XKSlpU2IZ2i32zGbzciyjNfrVZ9NT2UeCAnfgOlc6VtFobe3F6/XG9EqneEQtuE8/vjj+P1+Pv74Y5qamtQoyGRCq9Wi1+tDvGWg9K3EUEczkQ6TuhiG5qclaQxIQJevP8SL9vb24nQ6ycvLo729/aLTlMYbRVHU33uY95NEyH6p4MHFDofjgsGdaBC24ZhMJtavX09aWhpPPfXUpNhVeD7BAnZDDWOBNVDx/7rsczs+g8U8Jht+vx+rNXBw7vqs67jZVoiCYKYhXf2OEIKGhgba2tpYsGABlZWVw049iCYajQa/34/dbsftdmMwGJAkCVt+LVrjIJJ07pDc9PR0rFYrfr+fioqKmA6rwzacxMREDh48iNfrZdeuXVFLuRlvguWEglh0gR/f4T3Xs01WwwkWXAyGZ4uSckjRJvGp/Tj9vnOZA36/n4MHD9LQ0IDZbA7Z5x9tgsepHz16lNzc3JCjTDSJLoRfRnx5Yl5raytWqxWdTkd7e3tMF3TDNpyf/exnGAwGNdlz48aNkZQrIrS3t6uHyQYxf2k4/UMMRwgxKT2qx+NBlmWmT58OBIIGd6ctRyAo7Qytw+zz+WhqaqKtrY2mpiaam5tjITIajYbe3l6Sk5PVYzKDHVuC2QFI+F2B7OiBgQG6urrUIWlHR0fMpgxhG051dTWLFi1i6tSp/OhHP4qpex8LFRUVfPrpp0DoWHp/W6B3PtSpCwkQBL87mfB4PEiSRGpqqrqQuLvnCCZZT7Oni9+37mFr2z41WKAoCgMDA2g0Gj7//HOqqqqiPnfo6+vD7XYzb968kIXpnso8BtsCywK9x3POVfXs7kan0zE4OEhubi5tbW0xCeaMOgv+4IMPKCsr47333uPmm28GAj94WVkZ1113XaTlGxcURaGpqYmlS5dSVlaGyWRSPU+w51DO8/gTPSlyJGw2G83NzUyZMiXkOBObNgmHx02ff5BEzbm9LRaLBbfbjclkIjk5ma6uLnw+X1RX4js7O9HpdOqxJcE5D4Akf/lQlOBTCowKuru7SUlJweFwsGjRIg4fPhz1jW6jGs6SJUtISkqitbWVG2+8Uc3zeuihh6Ih37jg9XoxmUxMmzYNm82Gy+VSC5bfMt3FlloTs8x+VmSey3maDAUuzmfRokX09fVRW1uL2Wymq6tLPVh3a9tHNLu7uS/jahI0WhRFwWw2qxGqefPmcfToUTweT1QNx+/3o9Vq1cBAcnIynZ2d2PJrUbxaOg8twpDeTdKMM0iShFYbKCY/Z84cBgcH1UNyo82oQzWLxcKyZcu44ooraG5u5ujRo2zatIktW7ZEQ75xwel0YjQa1blN8CRmOLe+dv5y7mQ0HEmSMJlMFBQUDNvwV2jKQUFwcrBdXRRtbW3FZDLR3d2N1Wpl1qxZUU8C9fv9JCYmqh7+QhvehBDo9Xo8Hg+pqan09PRE3dCDhD3H2b17N3fccQcfffQRpaWlYRcvnwg4nU6EEGom7dA5TjBOcH4ixGTMIAC+7JW1LFiwIGS9aq4xsGB42HFSHTUEh0Q+nw+n08n06dPp7++P6tlHPp8vxHBCESH/ALjdbvW0DK/Xy8DAwMQ2HCEEu3fvJi8vj6amJux2eyTlGje6u7upqanBbrcza9YsIGBICQkJ7Gw08JfmQEDgpF0TEhyI1T6PSyV4goHBYFADBFvb9rGtYz9aZOoGz/Lq2Y/Y0roXCEzOtVotNTU1atHz4GG10cDn82E0GtUk1b6+PjUQ0FsTyL9zdaaobSaTCVmW1eBNV1dXTNKGwjac733ve1RUVPBP//RPfPTRR5MiHD0wMMAXX3yB2Wxm8eLFNDU1qUXtRisA4Xa7Q9Z7Jgs6nY7GxkbcbvcwHU0aAwLBoOJRO3EhBG63G6fTidPpZMqUKSiKErUkUEVRQqqmhiwDSEEZJbXBbrej1+tJT0/H5/PR3t4ek9pxoxrOe++9BwTC0VlZWfzpT3/C7/dPiqFaUlISM2fOBAJjZ6PRSFFREYmJiQgRqN552wwXIEhPVEK2FgSLe0825s2bR3p6OocOHWL69OlYLBbWZa7kgaxreWTKTSRIGgYVD/dnXIPNZkOn06HX6+nq6qKiooLBwUHmzp3LiRMnotJxJCQkYDabyc3NRQgROG4yvxZbfi3TlrWApJBg9DCj+CzpixsCWQU2mzrq0ev1JCYmRj2MPqrhBAv0JScnk5iYSHJyMhaLZVw2sv31r39l06ZN6uff//73rF69moceekhdvR+p7WLIy8vDaDRSUVGB1+tFq9Uiy7IaqpUkSJDB4Qud5MQ6a3isyLLMzJkzWbFiBQkJCUiShCzLCCFIkDXcYCnAobg52F9LX1+felLdwoUL6e/vZ+/evbS2tmI2m6OyuKjX69U9Np2dnSHzMovFgpzgw+fSodFoSE5OxmAwqF5pwYIF6HQ6Tpw4wccff0xVVRVdXV1RMaJRDSc/P5/vfe97bNu2jdLSUl5++WUOHjzIjTfeeEk3fuWVV3jhhRfUsXRbWxsHDhzgzTffZO3atfzxj38cse1ikSRJnSh3dnaq6wVDF3B1slALdgSJxXbc8USn07FgwQKuvPLKkMnzEvNszBoDH/VV4fJ5VMOqrq4mJyeH9PR0zp49S1dXF7W1tREvijHUcE6fPq0GNyRJwul0ojG4UTyBYz0yMjLweDwMDAyooWur1UpxcTFLly7FaDRSX1/Pvn37qKyspKOjI2IlpUY1nJ/85CfcfvvtvPHGG7zxxhvs2LGDVatW8dRTT13SjWfNmsUzzzyjfj569ChLliwB4Morr+SLL74YsW0syLLMokWL8Hg8NDU10d8fmjFs0AjcSqjh+Hy+SRsgGEpwfnduS7LEN1OuwCcU3us5gkajQa/XI8syPT09apRKlmUsFgv79++npqYmYseeBA0nuF1FlmXVmAcGBtAanSBkBgd8auTNYDDgcrlCjmwxGo3MmjWLZcuWUVxcjNls5tSpU+zbt49jx47R3t4+rkY06iC+vb192Ca2VatWsXXr1rBv8uabb/Lmm2+qnxcsWMAzzzwTUsdsYGBAneQFoywjtY0VjUaDxWKhtbUVn8+HLMv86VSgJx7wSniVc3tygnOdjz/+mJtvvnlCnTw9Fmw2G3q9nt+f+ovappO0HHGcpq/RSaLRyODgII8Zbufqq6+mpaWF+vp6nE6nukh64MAB8vPzSU9Pv8CdLh6DwYDD4aClpYXExEROfprx5QseKG/l9wYSPrsr5+Co1wB5CMDv87FggWfE6J/BYGDGjBnMmDEDt9tNe3s7jY2NVFZWMnPmzJDdsmNlVMP5qq2pFzOOXL16NatXr77gd5KSktRTnp1OJyaTacS289m2bRvbtm0DAkeRVFdXf+U9UlNTQ7ZFn79dPZjKMZRIZQ53dHSwcuXKYe0Xo89oDNUnNzeX53NzR/1/gnlfwUTRIEajkc7OzhGz4sdLF6vVyoq7Ac5/57x0dw+c92w06ppcOL+RyWTCZDLhdrsv+P2v0uV8RjWc06dPs3nz5pC24FF640lBQQGvv/46QggOHDhAYWHhiG3ns2bNGtasWTMuMtx5552UlpaOy7VG48knnxzxAU1Gff436HI+oxrOVxUevNTgwPlkZWVRXFzM2rVr0ev1/OpXvyIlJWVYW5w4E4FRDWfZsmURu/ny5cvVIzwAHn744WElqEZqixMn1kzuWe84M17Diolyr2jp879RF0lM1mzGOHFiSNzjxIkzBiZfMlYE8Hq9fP/736ejo4Pc3Fx+/OMfR6zgosfj4frrr1fXEkpKSigoKBi363+ddIHo6XOxusQ9DpE7oHckamtrueWWW9i6dStbt24d9xft66QLRE+fi9UlbjhE7oDekTh+/DgVFRXcf//96oHE48nXSReInj4Xq0vccBj9gN7xJDs7mw0bNvDaa6/h8/nYuXPnuF7/66QLRE+fi9UlPsdh9AN6x5PFixerWxZWrlw55sTVr+LrpAtET5+L1SXucYjMAb1fxa9+9St2794NwKFDh5g7d+64Xv/rpAtET5+L1SW+jkMgclNSUkJrays5OTls2rQpYpGo9vZ2nnjiCRRFIScnh6effnpcz3j5OukC0dPnYnWJG06cOGMgPlSLE2cMxA0nTpwxEDecOHHGQNxw4sQZA3HDiRNnDMQNZxQ+++wzVq5cybp163jggQd48MEHY3bU+VtvvQUEas2dPn06ZtcYyttvvx3RtJ6R2LRpU0QzIsJCxLkgZWVl4qc//an6+d133xU///nPYyLL7bffPiGuEcTj8YjHHnts3K4XLuXl5eI///M/o37focRTbi4Sh8OBXq+ntLSU0tJSfD4fv/jFL3jqqacC5VttNp577jl2797Nnj17GBgYQFEUfvOb3yDLMk8++SROpxOHw8GLL76IzWbju9/9Lh6PB7PZzPz583nkkUf44Q9/SHd3N319fTzzzDOUl5fT2NjIf/3Xf1FXV8f69etxOBy88MILSJJEcXExGzZs4Ac/+AF6vZ76+nqsViu//e1v1QXDLVu2DLvGli1b0Ol01NbWsmzZMtrb2zly5AiPP/44N954Iy+99BKffvopGo2GjRs3hqyo79u3j0WLFgHwwx/+kObmZjweD3V1dXzyySc89thjuFwuTCYTv/zlL+nr66OkpASn08ncuXP52c9+xi9/+UvKysqQJIl//dd/paioiJtvvpmMjAzWrl2L0+lUveSGDRsoLi5m8eLFPP/883z729+O/gsQJKZmOwkoKysT11xzjfi7v/s7sW7dOlFSUiLsdrt46623xMaNG4UQQvzkJz8Re/bsEUII8cc//lG8/PLL4q233hIlJSVCCCG2b98ufvvb34ojR46IAwcOCCGEeOWVV8Srr74qXn31VfHf//3fQgghNm/eLH7zm9+I06dPi927dwshhNizZ4/YvHmzEOKctygpKRFVVVVi9erVoq2tTQghxKOPPiqOHTsmSkpKxNtvvy2EEOLv//7vRV1dXYg+518j+H1FUURxcbHo7OwU1dXVYsOGDaK6ulr1KK2treKhhx4KudYLL7wg9u7dG9L2xBNPiL1794qamhrxne98R7jdbrFv3z5x9uxZsWnTJvX7v/vd78SRI0fEo48+KoQQoq2tTdx9991CCCGKioqE2+0WXV1d4r777hN+v1/09/eLe+65R73P+vXrRXt7+0U8yfEl7nHC4Oabb+bJJ58c1h48Pu/UqVM88sgjQCBZsLS0lMLCQrUK6cKFC9m3bx933HEHW7ZsYfv27XR2dnLNNddw+vRp7rjjDiBwolp9fT0Wi4UPP/yQ999/H5fLRWZm5ohyeTweMjIy1PueOnUKCNTLhsDx5uGU8s3Ly0OSJDIyMkhNTWVwcBC3201DQwNVVVWsW7cOYFg1z76+vpBaZ//xH//BnDlz1CMuV65cybe//W2sVis/+tGPOHnyJP/wD/8ABIqw/PnPf1Zzz4LlbSGQqazT6WhubqaxsZH169er9wseJJWamkpfX9+4F0gMl3hw4BIIDoFmzpxJRUUFAH/729/UYn7BYobHjh1j1qxZbNmyhVWrVvHcc8+pZ/XMmDFDPfkh+O/27dvJycnh+eefp7Cw8CvPqklISKCjo2PYfS82l+urvj9jxgyWLVvG1q1beemll7jllltC/m61WtUywbt27aKhoUHtQI4fP45Go+Hll1+msLCQXbt2MX36dLUY4LPPPktmZqb6u7W3t6u5YcHKqdnZ2Vx22WVs3bqVP/zhD9x2221qHez+/n5sNttF6TmexA1nHHj44Yd5/fXXue+++zh8+LDaQ1ZWVvLAAw+wY8cO1q9fz7XXXsuvf/1r7rvvPlpbW+no6GD16tV8+OGHrF+/nmPHjqHRaFi+fDlvv/029957L0eOHFGrZ5rNZl566SX1viUlJTz22GOsXr2agoICdb5xIc6/xoUoKCggPT2d+++/nwcffHBY6diioiLV2H/wgx9w9uxZ1q1bx7p167BarezZs4e1a9fy4YcfctNNN/Hwww/zyiuvcP/99+PxeFiyZAmzZ89m7dq1PProo/z7v/97yPXT0tK44YYbuPfee7nnnnvIyspS/+ZwOEKqskabeJJnhCgtLcVut/Pggw9e8HuHDh3CYDBQUFDAc889R15eHnfeeWd0hLxEPB4Pjz/+eNiGOF6Ul5dTVlbGd77znajedyhxjxNjsrOz+elPf8pdd93FyZMn+cY3vhFrkcJGp9Nx/fXX8/nnn0f1vjt37lTnXbEi7nHixBkDcY8TJ84YiBtOnDhjIG44ceKMgbjhxIkzBuKGEyfOGIgbTpw4Y+D/A1/Nb6yJSPRbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 170.496x63.36 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_pal=fs_pal\n",
    "fig, ax = plt.subplots(figsize=(3.7*0.8*0.8,1.1*0.8))\n",
    "for i in range(nc0):    \n",
    "    ax =plt.subplot(1,3,i+1)\n",
    "    dpl0 = np.ma.masked_array(np.squeeze(z[:,idx_sub0[label_sub0==i]]),np.squeeze(features[:,idx_sub0[label_sub0==i],0])<Th)\n",
    "    ef.plot_errorbar_vertical(ax,dist,np.transpose(dpl0),my_pal[i]) \n",
    "    plt.plot(dpl0,dist,'darkgrey',zorder=-10,alpha=0.8,linewidth=1)\n",
    "    plt.ylim([-100,100])\n",
    "    ax.set_xlim([-3,5])\n",
    "    if i==1:\n",
    "        plt.xlabel('Propagation time (zscore)')\n",
    "    if i==0:\n",
    "        plt.ylabel('Distance (${\\mu}$m)')\n",
    "    else:\n",
    "        ax.set_yticklabels([])\n",
    "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=1,\n",
    "                wspace=0.2)       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAABkCAYAAADOmA9MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoQUlEQVR4nO2de3hU1bnwf3vPNTOTTO43CBBCQoAEwjWIclHw0lPaz6oo6lF8jk+11WOLtjb1wkGPFrS0Sqs95+v5jlrK0UdqCUhtQbGACCUgkBDQcL+GhEwuJJlr5ra+P8bZzTCDCXGSDD3z+yeZtfZee62Z/a7Lu971vpIQQhAnTpwrQh7sCsSJczUSF5w4cfpAXHDixOkDccGJE6cPxAUnTpw+EBecOHH6gHqwK7B8+XKmTp3K7Nmz+fGPf0xzczOjRo3i+eefx+v1hqVJkjTYVY4TZ/BGHI/HwxNPPMHmzZsB+OijjygsLOSdd95Bq9Wya9euiGlx4sQCgyo4d9xxB9/5zncAOHDgAOXl5QDMmDGDvXv3RkyLEycWGDTBMRgMzJgxQ/lss9kwGo0AJCQk4HA4IqbFiRML9GqN43K5+PDDD6murqajo4OUlBQmTpzIvHnzSEhIiEpFTCaTIhgOhwOTyYTf7w9Lu5Q1a9awZs0aABYvXsysWbOu6LlHjx6lsbGR9vZ2urq6lPSilZUgBKYTjQiVjD0/++/3LL4NAEmSUKvVDBs2jDFjxmAwGK6s0RH4uu3pLZ999hkXL16kvb0dgDmW5wFQ+V2kek4iUNGqHYWQVMo92zKXhpQhyzJz584lJSUl4jOupC07duxg0qRJWK1Wjh8/zvTp01GpVOzYsQO1Ws3kyZNZv349F/70inJPl+UkALrMkUpa1vzHkWWZoqIivF4vnZ2dZGdnM2TIEBITE6/gG/pqpJ5s1d5++232799PeXk5BQUFpKWlYbVaqauro6qqirKyMu6///4+V+C1115jzJgxOJ1Ozp49y6OPPsqLL77I7NmzaW9vD0ubOXPmZcuqq6tjzJgxV/T84AtktVrx+XwheSn7jjLyzQ85fc8NtF47LuzexMREpk+fftkX5+vSl/b0lpaWFrZs2RKSpvNe5J8u/BBZ+NmY/QoOTWZIviRJBF8XlUrFjBkzyMnJ6dXzemrLzp07GTt2LLW1tUyYMIHk5GQAtmzZQmpqKikpKezevVu53tPZTOMflmIaO4fUGXcq6bIcmEQNHz6coUOHkpmZiUqlItr0OOJMnjyZe++9Nyy9tLSUO++8k8OHD0elIrfccgsVFRUsXLiQ/Px8rrvuOrxeb1haNBFCcPHiRTQaTZjQAOT8eQ9+jYrW6cVhecXFxZSWll6VWj4hBHv27AlJU/nd3Gh5GrVwsyVjqSI0QWGRJAmVSoXf70cIwdChQ3stNL3l+PHj5OTkKEIDgbWwXq/n+PHjIddaD20FILHkhpB0tVqNyWRi6tSpUa3bpfQoOMXFgZfmlVdeYevWrajVauWLXLdunZLfVx577DHl/1deeSUkT6PRhKVFE4fDgcfjITMzk+bm5pA8XdNFEpou0jxjLFzSY82ZM4fMzNDe+Gqhs7OTDz/8kJCJhvAzu/l5EnxtfJbyMK36wMggSRI6nQ6Xy4XZbCY3N5fDhw8jSRJTpkyJar26urrw+XxMmjQpJN3j8eByuWhtbQ1Jd5zch8pgRpOUrqTl5uai1WoHZC3c632cqqoqNmzYcFX2sJejpaUFr9dLW1tbSHrRykr0Da0IIKGhlaKVlcq6BuDixYsAV53wHD9+nP379wN/X9MAJHrOk+DvwC6nMsLxKSMcn7ItcylCCFwuF2q1Gq/Xi9Vqxe/3R2Ut151g2SUlJbS2tuJ2u1GpVKSkpOD1ejl69CgATR+8CoDf04XfZUU2JCtpWfMfp62tTekQ3G43Wq0Wn89HfX09w4cPj2qdey04s2bNYvPmzQwbNkxJ+7qjzWDT3t6OWq0O680QArXDhV+rRmjCv6L6+vqrTmhaW1v54osvwtL1vjYS/B10SUbs6qyI93q9Xmw2GzabDYDs7GxlLRENJElSBPLgwYPK1Lm0tDTic7ydFgBUhuSQdJfLFWiTXq9MK2tra9FqtVGra5BeC05dXR1HjhwJ0WwtX7486hUaSJxOJ7Isc6l+pOmGiSQeb+DsnbNpu2ZsSF5GRgbXXXcdGo1mIKv6tfB6vVRXVysvFgQ0ZCldx5hnWYJdlc7GnF/hlyK/DrIs4/f7kWWZCRMmUFhYGNX6+Xw+tFotdrudpKQkWlpayM3NxW63h1yXNf9x3BcbubD2BRJGTCRj3nfDysrKylIE+9ChQ0iS1C8Kll4Ljsvl4o033oh6BQYTm82G0+kMS8/6636ELHFxSlFY3syZM1GrB91S6YqoqakJm45qvR3MaX4RH2r+mvliRKGRJImUlBTa2towGo34fD7y8/OjXj+73U5ycjKNjY3k5ORw7tw50tLSaGpqwuv1hlzb9unbgERKN01akPz8fGRZJjk5mSNHjtDV1cXEiROjOjoG6fUbYDQaWbJkCaNGjVLWOV9HDT3YdHV1KXP27khuL6ZTF7AW5EacptXV1ZGfnx9xTynWEEJQU1PDyZMnQ9Il4WNu81LUwsXWjH/DpQ5Vp6tUKoxGIy6Xi/b2dkaOHElGRgaHDh3CZrOFaL2igc1mIyUlBavVSktLCxqNhpaWFhoaGkIEx912HrflJIb8SagN5pAyzGYzmZmZnDhxgsTERNrb25k8eXK/qKLhCgTn+uuv75cKDAZCCGpra1GpVCE/TNHKStQddiS/QG0LH4lkWebChQs0NzeTmppKQUFBVDfVoklrayu7du1SNEzdlQFJnnPo/VasqkxKOv/INn3oHpXP58PpdOL1eikuLqawsJDa2lry8/M5evQoU6dOjaqSyGq1otVqyc3Npbq6mtzcXM6cOaNsSAcVAO7WcwB4bG00ffAqWfMfV8pwOBycPXsWIQTNzc1MmTKlX6fTvRac/Px8iouL0ev1QEBrcbXi9Xppbm5WetXuaNttCAl8Rn3YfUajEY/Hg1arpbW1laSkpJgVnN27d0dUyxq8zej9VlxyIk51eoQ7A3i9XpKSkjh9+jRnzpxBo9Fgs9lISkpStiOihdVqxev1kpeXh8vlQqfThVhxAPg9boTHhawzIavCBcLj8WCxWDCZTEyaNAmdThe1+kWi14Lz/vvv89JLL+HxeBg5ciRHjhxhw4YN/Vm3fiO4F3WptuXEd/+JCT/5f7SPH8nJh78Zdl9XVxd5eXkMHTqUtLS0mF3reL3esA5hW+ZSsp3VzGx5iQ71ED7KXhFiTtOdoDLAZrNhMBjIyckhOzubtLS0funF9Xo9bW1tnD9/HgCLxRKSnzX/cSwf/geeVsj+zlOoTeGWGjk5OajVavLy8qJmBvZV9PqXX7o0YKfk8/nYsmVLv0t0fxLcBb+U9E8PIgEXbpoUlqfX67npppuUETeWOX/+fNii2uRp5LqWFXgkA1uyXris0AQZNWoUBQUFJCUl9fveXUJCAm63m87OTmRZxmAwhGwR+Fw2XPWfo8seFVFo5s2bh9Pp5Pz58wwZMqRf6xrkitUNKpWKG2+8UTEOvFpRqVRh6s6MnZ/jTdDhGJEdkp6VlUVaWtpVITQ+n4+6urqQNLXfwVzLswD8NfMFPLIx4r0mkwmTyUR5eTmTJk3CbDYPyIZ3QkICLpdL2bRsaWkJyW/f+ycQguRrFoSkS5JEeXk5JpOJ48ePM2bMmAHboO/1iPOd73yHjIwMxowZQ1paWr9pKwYKlUqF1WpVPmtaO9G1WQMmNt2+fLPZjNPp7Bc1bH9w9OhROjs7lc+S8HODZSlav40d6T/Gqh0acr0kSSQlJTFy5EhOnDiB0Whk6NChlxbbr8iyjNvtRgiB1+sNGS2F14P92C7U5ix0aXkh9+Xk5JCXl0ddXR3Z2dkDut7steCsW7eOtrY26urqqKurQ61Wc/vtt6PRaHj33Xf7s479gs/nC/mBRq9cC0BCfWhvZzabaWxsJCMjY0Dr11c+//xz5f85lucxeCyY/C3Y5HRGW/9CY0Ko8aNeryc9PZ0TJ05gs9mYMWPGgJtVnTx5UtmE7v6bNH3wKj5nJ/i8gBymSRs2bBher5eGhgauueaaAa1zj4LzwQcfMH/+fABSU1O59tprufbaawHw+/1XrYJg/PjxfPLJJ8pnlcuDAIQ29CtpaWkhKysr4kZprNHV1RW6LyX8GPyt+FDjUIcLvtFoxG63c/78edRqNRMmTCApKWkAaxygrKyMpqamsD01CNilAci60AW/SqVi6NChytmbAwcOMGXKlAFbe/coOFlZWfzgBz8gIyMj5DzOF198QVNT01W7CRo8pxE8TuA1JSBUMkcfvz3kOoPBwIgRI7BYLAO28Owrl9rcWfSlZLoPsyNtMQ2G0JFGrVbT1dVFYmIi48aNo6GhgVGjRg1kdRUSExO57rrr2L59e0h61vzHsXz0f3Gd6yDrWz8KGQkTExNDzt74fD727dvHlClT+sU27VJ6FJypU6cydepUjhw5QnV1NadPnyY5OZm777476jZLA0lQsxYUHE2nHVdGStg1sixjNBppa2vD6/XGrAoa4Ny5c8r/svBSbN2AQ5VKQ0L4EQCz2YzD4aCkpIRTp04xfvz4QbV8v9xU2Ge/iKTRh9XN7XZz+vRpRowYAcDIkSPx+/3s37+fyZMn97stYa/fgtGjRzN69Oj+rMuAo8ynfT5klwdXVnJIvhACjUZDR0cHZrOZ1tZWsrIiWxAPNp2dndTX1yufi6wbUIsudic/GqLsABg6dCgXL15ElmXsdjupqalRN6O5Uvx+P5IkYTAYQrSdfpcVWRd+jEGr1XLu3DnF4gCgoKAAn8+nCE9/dnL/Kx0SdnR00NzcrMypR//ij0iA8dSFsGs7OztpaWnBYDCEqUljhZaWFvbu3auMnnOallLS8R4+1BRaN4ZcK8syHo8Hu92OTqfjwoULMTFzCFqidD+GfuFPr+Czd+B32RSzGwjMBFwuF2PGjKGurk7ZGpEkiYKCAgBqa2v7tb69Fkm73c7q1atpbW1lxowZFBYWDrjaMloIIUL2OjSdgR7OG8HMxu/309nZicfjoago3Fo6FnC73Yo5PoDR24yMn3Z1btho4/f7aWpqUgw5i4uLY2IzOyEhAZ1OF3ISN6AYEEja0H2noAr9888/Jz09XbnvwoULtLS0kJSU1O/r0V4LTkVFBd/4xjfYvn073/rWt6ioqODtt9/uz7r1G8nJySHWzbLPjzvJwJGKu0KuS0xMxGw2U1hYSHp6esyefs3MzFR8P2j8dvT+djrVOXyU88uwa1UqFSNGjGD06NExZeEtyzIzZsxg69atSpo2JRtv2zlybq1AlRDYo1GpVJjNZoxGI5mZmdhsNnbv3o3RaCQ7O5uioqIB6Qh6LTidnZ1885vfZM2aNYwfP75fzjhEwuPx9Isb3IKCAk6ePInmfAsaq5PGeROVPEmSyMrKori4mIyMjJgVmCBqtZr09HRaWlqY1PYGMj72pD4aco0kSYwdO5aioqKYPYSXnp5OXl6eYuXsPHsQdVKmIjRqtZqUlBSSkpKw2WxcuHCB7OxsCgoKBtyqo9dvf3p6Om+//TZ2u50NGzb0m0ukS+kvN7hmsxmz2Uz2x/sQQPP1ZUBA/Txz5kxmzpxJZmZmzAtNkNzcXBI8zQxz/o0W7WjadIF1iyRJTJ06lTvuuINx48bFrNAEmTZtWkBV3nQS4XFhGh3YM1SpVJhMJnQ6HSkpKUydOpXy8nKGDx8+KKZQvR5xli9fzh/+8AfGjx9PW1sbL730Un/WS+HAgQPceOONwN/d4Hb3APp1mDhxIq0HXsKdloQn2cSUKVPIz8+/aoSlO6mpqZRf/A8Adn852syePTtmtYCXIzhlW7vkNwCYiq/FaDRSVFRETk5OzEwvez3iyLLMtddey9KlSwd0QdmTG9w1a9Zw2223cdttt4UdD/4qPB4PWz7+GJXTjS0/C71eT15e3qALTV/bc+jQIVLdJ7ioySdlxCTuvPPOQReavrYlJSUFr60NSa1F0iZQWFhIYWFhzAgNXIHgPP7441RXVwNw+PBhnnzyyX6rVHciucbtzl133UVlZSWVlZWkpqb2qkyLxcL7778Psoxfo0bXZsPlcrFu3TqOHTs2qIf0+tIep9PJ4cOH8SMjCBwr2LNnT5gTkoGmL23x+/28//77JAwvRXjdeNqbqKmp4b333uPYsWOD3qYgvRaclpYWbr89YI7ywAMPDNieRmlpqeL6dNeuXUyYMOFrlVdbW8u2bduUPRyfSY+m3abkV1dXs3XrVo4dOxZ2CjEW8fv9/OUvfwFASGpUwgPA6dOn2bJlS9i5nFjn0KFDABgLAlEq7MeqgMAWQnV1NZWVlRw9ejSiXdtA0mvBSU5O5p133qG2tpY1a9YMmAn3LbfcwokTJ1i4cCF2u73PbnC9Xi8ff/xxmMtet9mI2hZ6WrKjo4P6+nr+9re/Kd5SYhEhBFu3blU2Pv2oUeFR8ltbW9m4cWPI8YlY5uzZs4qrW405A0mTgPPMgZBrfD4fNTU1rF+/nrq6ukHrGHqtHPjFL37BunXreP/998nLy2PFihX9WS+FaLjBbW5uZseOHYpvsO69VVeGGePpJvD5QSWH3GMymfB6vezatYvs7GxGjBgRU4fZamtrQww7/ZIKWYS+SE6nk48//pipU6fG9IZ10GFid0HQZY3Edf4LhNeDpP67NlCSJLxeLwcPHqSuro7CwsIB278J0usRx+l0otFoGDZsGJIk8cc//rE/6xVVdu7cidFoRKfThcyRi1ZWknj0PBIw+tW1YffZbDba29u55pprsFqtVFVVxcwc2+FwhDgin2N5Hr2vHb2vPcSjDQQUIfv376ejo2Ogq9lrDh06FNKhNX3wKl5rCwhB4/s/DzG5EUJgNBqRZRmfz8fRo0fZuXPngP42vRacRx99lM7OTpKSkmLau0skcnJymDp1quJpvztquwu/LCHU4Sda9Xo9EydO5Ny5czidTsrKygZd6xbE5XKFRFiQhA8ZLx453CBSo9EwduzYQTlr0xv8fj92uz3sOICkDowgwhuurJFlWXEKL0kSEyZMGNDfptdTtcTERB566KH+rEu/kZKSwsWLF8MWlKfun8f4JauwXFfCubtD/cYZDAbmzJlDXV0dQgjKy8tjwqYryKXq3VPG68lwH6Eq7QdY9KUhecXFxYN21qY3uN3usDArWfMfx+e0cv7tCgwjJ5Iy9daQfKvVqhwLGT58OGlpaQNY4ysUnKvVk2dKSgqnTp0KUzNnbgksPC1zy0LSdTodRUVFVFdXk5WVxahRowbMxKg3uFwuampqQtJG2D/Bj4pm3diw62NZaAAlqFd3XwkAsj6wf+ezRt4D8vl8aDQazGZzxPz+pNeCM3v27JDPsTJl6QmHw8HJkyc5c+ZMWF7avqO4kwx0ZYaaD3V1dXHu3DlKSkoGfRPxUs6ePUtVVVVoovCT5j7ORc2IiG6fduzYEdM+r8+ePYvX6w13RyzJSCoNPsfl12bBMB5ut5vCwsIBMynq9Tc5bdo0tm7dqqwTYvVsSnd2795NfX192DSgaGUlkseLptNBV7IxLP5NkFgL5bFv3z5OnDgRkjbH8jxqnxMVHrR+G3Msz4fF6uzo6GDPnj1RM1WKJk6nM8RZR5CgMkD4/XQ1nw6Jg9Mdj8eDTqdDlmX+9re/MXny5AGxMOj1/OOJJ57A5/Px6aefcu7cuTBH3rGI0+m8rPpY9gTUnj5D5PxYtFnT6/XIsoxGowmdOn5ZTf9lnAwmJydHDNUYC2i12q8MVCVpdODzIPzh9Q/GavJ4PBQXFzN58uQB2y7oteAYjUYWLVpEeno6zz77bEyrNoOo1eqIU62ji2+j4RvTADhzzw0RR5tY1BoGe9Lue1HbMpeyPeNpABoSJkeMDK3RaGLKzqs7KpWKOXPmhL3wWfMfJ2v+42Tc9D0AEoaPDxttgq5yvV4vPp8Pk8k0YNPRXgtOQkICe/bswePxsHHjxpifqgkh6OrquuwXqbEHrAX8+sgeUfbt23dFhokDQTAI1qVrAbeciB+JRE9j2D2SJNHa2hpzo2d3jEYjhYWFEX8rXVYBkjYBW92nYXkulwtJkmhvb2f79u2cPHlywOwMey04P/vZz9Dr9Yqx55IlS/qzXl8bj8eDSqW6rKte1ZeC49NFXkwGw+DFkqvfoLdLj8cTmiFJeCUDRl94Z+bz+fB4PDHf0RUUFKBWq8O0l5IkYRg5GZ+tFU9neBs0Gk3g0JvTyZkzZ/jkk0+oq6vr9wC6VxTKMOgt8emnn2bbtm39Vaeo4HA4cLlcYSrOopWVAGibA1PN/Lc+BJUcNl2z2WxkZGRQU1PDxIkTB0XleSnBoL3dCVoJyMJDkqde+dx9yubz+Whvb6epqSnmtIRBtFotGRkZNDQ0hKQ3ffAq/i83QJv+9Es05oDCJjhtC7rsEkJgNptxuVycPn2aU6dOkZ2dzejRo5VjKdGkR8HZsmULVVVVfPTRR9x8881AoDeuqqpizpw5Ua9QtLDZbIrLoUimGLLHiwCQLz+FaWpqIikpKWaMPC/tBLojJAlZXN5iWJblmHZvBTBkyBCamprCFBmyWguSjN8d2Zuq3+8nLy8Po9HIpEmTsNvtNDQ0YLfbOXbsGGVlZVGva4+CE1TvNTQ0MG/ePIQQyLLMgw8+GPXKRBtJktBqtSEvfnBkGfdvqwIjzSWeO4OYTCbMZjPXXHNNTGx+ejweJepzd4Ijyz81PIYsvGHKgSAmk2nAjrv3FZPJRGpqKhcu/N1NV3BkOb9mCfh8YQqChIQEZXukq6sLk8lERkYGhYWFOBwO9u7d2y917fGNMJvNTJs2jSlTplBfX8/BgwdZtmwZq1at6pcKRQuLxYJWq42sTRICbbsNZ3bkw1Umk4mbb74ZtVpNY2P4gnugCUaN1mq1l13kq4ULr3wZ1fuXWrhYMhmKRDDcRyQXtrLWgN8bOvIHv4uCggIlavXhw4fZuXMnZ86cQa1W4/P5+mXG0OuudNOmTdx666188sknVFZWhnjFjzWEEFgsFnJzcyPOb7VtVmSfH9uo3Ij3FxcXo1KpGDt2LMeOHQuLbjaQ+Hw+Dhw4gNFoxGAwXNYCWCXcEQ08ASVkY6wLjl6vp6urC4PBEDbKyzoDwhuqFAkqBYYPH861117LuHHjlHCT586dY8eOHfh8vrAIb9Gg18oBIQSbNm2isLCQc+fOfeV8e7BpaWnB7/eTlpYW4uAOvgyQ2xnQuJhrT5J47HyYYiB4zDcxMZG8vDwOHz484Na3EPjODx06hFqtZsiQIRGtBoKohQuTpzGiciAYxcDtdvfLQjlaBL9fjUaDTqfD6XQqFgOe9kbwe7nwp1eU67LmP47f76empoYZM2aQkZFBRkYGnZ2dnDlzhubmZnw+HwcPHqSlpYVhw4aRmpoald+x1yPOj370I2pra/nXf/1XPvnkk5hWRzc2NpKQkKBEib4UlTMwdPsvo4pubW1VeiqbzabMnwcaj8eDx+Nh+PDh1NTUXH7E+HIU8kuR+0G3241GowlXY8cYwdHU7XaH78fIX7btEgWI0WjEYrFw+PBhRamQlJREaWkpM2bMIDMzE1mWycjIoK6uTokz+nWRRA+nfz766CNuuukmZU0T1FJJkhRz1tF1dXWMGTNGCUEyevRoduzYEXZd0cq1GE81Uf2rR8LySktL8Xq9tLS0kJiYSFZWFpmZmYMyzdm3bx+5ubm4XC66urrQ6/UcOHAgTIh1vg7+T8NDfJF4K4eS7w7Jk2WZUaNGDcqI2Z3gb9MT58+fZ9++fciyHLIX03nwr7TvXsuoB1aQnBswhyorK8NsNuPxeDh8+DCdnZ2UlJRcduugvr6ejo4Oxo0bFzH/SuhxxAlG8E1KSiIhIYGkpCTMZnNUTFI+/vhjli1bpnz+r//6LxYsWMCDDz6o7FlESusJnU6H3++/bA+r6bDjM4QLQjAuaHJyMjNmzGDKlCnk5eUN2trA5/ORnJxMW1sbqampeL3eiCEkdb7AnpRdHW6UqtFocDgcMW050J0hQ4ZQXFwMELLOUScGQss7O1rIycnBarUq76BGo6G0tJRRo0ZRU1PD8ePHIzrzCAbSigY9Cs64ceP40Y9+xJo1a6isrOTNN99kz549zJs372s9+K233mLFihXK8NzU1MSuXbt47733WLhwIb/73e8ipvUGvV7/lYKjtrnwJEYO6V1SUjKownIpOp2Ojo6OgK+xL22yLkXvawfArgrfowlGc471aVp3gg79uwuOypgMgOuiBYvFgl6vp6mpKeS+rKwspk+fjtVqZc+ePWFOSoKCE40j1j0KzgsvvMC3v/1t3n33Xd599102bNjA3LlzefbZZ7/Wg0eMGMFzzz2nfD548CCTJ08G4JprrmHfvn0R03pDT4KjcrlxJ4ePmD6fj88++yxmLInVajUdHR0YDAZUKtVlBcDkC7xALnVyWJ4QAoPBEHGtF6tIkkRpaWnI2RqVPrCt4LW20NjYyKhRozh27FiYlxudTkdZWRl5eXns27eP06dPK4Ki1WqRZTkq69UetWoWiyXsENvcuXNZvXp1rx/y3nvv8d577ymfx44dy3PPPaf4S4PATn9wz8VgMOBwOCKm9Yag4FxqpFm0shL8AskvSKhvjngOx2KxsGXLFsrLywf9jL5KpaKjowNZltm+fXtIjM/uGjWDJ6Bundz2W8SXCoLuWjWn00ljYyN5eaFRm2MZlUpFUlISp98LTOXFl0oBa912XOcPs/FDNTqtjiFDfqvExAkiSRJDhgwhNTWVzz//nIsXLzJxYsCpvtFoxOFwfO3jBz0KzuXCsl+JQ7gFCxawYMGCr7zGZDIpGg+Hw4HRaIyYdilr1qxhzZo1QCAUSTDuTTDEXYgHyTcCZ/HdBI6w6IDup/ODawkILFKjpYGJRHNzM7NmzQpLv7Q9Wq2WxMTEsDVlK3/s9n84pYS2BwiJCRRNetuWK31+ZmYm8556EwhvSxC32/2V5QbfmeA1RqORpqamsGleT225lB4F58yZMyxfvjwkTQjB2bNneyz8SigpKeGdd95BCMGuXbsoKyuLmHYpd911F3fddVd4gX3gtttuo7KyMipl9cQzzzwT8Qe6Gtvzv6Etl9Kj4FzO8eDXVQ5cSnZ2NtOnT2fhwoXodDpWrlxJampqWFqcOLFAj4Izbdq0fnt4eXk55eXlyueHHnoozAVVpLQ4cQabwTf7jSGiNa2IlWcNVHv+N7alR8uBOHHihBMfceLE6QOx6aFugOmvAL2RcLvdXH/99YwcORIIqGlLSkqiVv4/Ultg4NpzpW2Jjzj0X4DeSBw7doxbbrmF1atXs3r16qi/aP9IbYGBa8+VtiUuOAQC9Aa1e8EAvf3FkSNHqK2t5d5772X58uVRjyz2j9QWGLj2XGlb4oJDzwF6o0lubi6LFy/m7bffxuv18uc//zmq5f8jtQUGrj1X2pb4GoeeA/RGkwkTJijGi7Nmzeq14Wpv+UdqCwxce660LfERh+gH6P0qVq5cyaZNmwDYu3cvo0ePjmr5/0htgYFrz5W2Jb6PQ0BzU1FRQUNDA/n5+SxbtqzfNFEWi4Unn3wSv99Pfn4+S5cuvawhbV/4R2oLDFx7rrQtccGJE6cPxKdqceL0gbjgxInTB+KCEydOH4gLTpw4fSAuOHHi9IG44PTA7t27mTVrFvfddx/3338/DzzwwGXPq/c3a9euBQK+5iJF0R6oMrqzfv36fjXricSyZcv6PXBUj4g4X0lVVZV48cUXlc8ffPCB+MUvfjEodfn2t78dE2UEcbvd4rHHHotaeb2lurpa/Od//ueAP7c7cZObK8Rut6PT6aisrKSyshKv18svf/lLnn32WTweDykpKbz88sts2rSJzZs3KwGufv3rXyPLMs888wwOhwO73c4rr7xCSkoKP/zhD3G73SQmJjJmzBgefvhhnnrqKdra2ujo6OC5556jurqas2fP8t///d8cP36cRYsWYbfbWbFiBZIkMX36dBYvXsxPf/pTdDodJ06cIDk5mddee03ZMFy1alVYGatWrUKr1XLs2DGmTZuGxWLhwIEDPPHEE8ybN4/XX3+dnTt3olKpWLJkSciO+vbt2xk/fjwATz31FPX19bjdbo4fP86OHTt47LHHcLlcGI1GXn31VTo6OqioqMDhcDB69Gh+9rOf8eqrr1JVVYUkSfzkJz9h0qRJ3HzzzWRmZrJw4UIcDocySi5evJjp06czYcIEfv7zn/O9731v4F+AIIMqtlcBVVVVYubMmeKf//mfxX333ScqKipEZ2enWLt2rViyZIkQQogXXnhBbN68WQghxO9+9zvx5ptvirVr14qKigohhBDr1q0Tr732mjhw4IDYtWuXEEKIt956S/z+978Xv//978X//M//CCGEWL58ufj1r38tzpw5IzZt2iSEEGLz5s1i+fLlQoi/jxYVFRXiiy++EAsWLBBNTU1CCCEeeeQRcejQIVFRUSHWr18vhBDiX/7lX8Tx48dD2nNpGcHr/X6/mD59umhpaRF1dXVi8eLFoq6uThlRGhoaxIMPPhhS1ooVK8TWrVtD0p588kmxdetWcfjwYfH9739fdHV1ie3bt4sLFy6IZcuWKdf/9re/FQcOHBCPPPKIEEKIpqYmcccddwghhJg0aZLo6uoSra2t4p577hE+n09YrVZx5513Ks9ZtGiRsFgsV/BLRpf4iNMLbr75Zp555pmw9Pz8fABOnz7Nww8/DASMBSsrKykrK1O8kJaWlrJ9+3ZuvfVWVq1axbp162hpaWHmzJmcOXOGW2+9FYDx48dz4sQJzGYz27Zt469//Ssul+uy4QfdbjeZmZnKc0+fPg0EXMgCZGRk9CoKc2FhIZIkkZmZSVpaGk6nk66uLk6ePMkXX3zBfffdBxDmAbOjoyPE19lvfvMbCgoKlBCXs2bN4nvf+x7Jyck8/fTTnDp1iu9+97tAwAnLX/7yF8X2LDMzU6lrbm4uWq2W+vp6zp49y6JFi5Tnud1utFotaWlpdHR0kJGR0WP7+oO4cuBrEJwCDR8+nNraWgBqamoUj5lHjhwB4NChQ4wYMYJVq1Yxd+5cXn75ZcVh4rBhw5QgXcG/69atIz8/n5///OeUlZVd1texRqNR4v90f+6V2nJd7vphw4Yxbdo0Vq9ezeuvv84tt9wSkp+cnKyEV9y4cSMnT55UOpAjR46gUql48803KSsrY+PGjeTl5SmOAV966SWysrKU781isSi2YUGf0bm5uRQVFbF69WreeOMN5s+fr0Rrs1qtgxqaMS44UeChhx7inXfe4Z577mH//v1KD/n5559z//33s2HDBhYtWsTs2bP51a9+xT333ENDQwPNzc0sWLCAbdu2sWjRIg4dOoRKpaK8vJz169dz9913c+DAASXUemJiIq+//rry3IqKCh577DEWLFhASUmJst74Ki4t46soKSkhIyODe++9lwceeEA5Vhxk0qRJirD/9Kc/5cKFC9x3333cd999JCcns3nzZhYuXMi2bdu48cYbeeihh3jrrbe49957cbvdTJ48mZEjR7Jw4UIeeeQR/v3f/z2k/PT0dG644Qbuvvtu7rzzTrKzs5U8u91OWlpar9rRH8SNPPuJyspKOjs7eeCBB77yur1796LX6ykpKeHll1+msLCQ22677SvviRXcbjdPPPFErwUxWlRXV1NVVcX3v//9AX1ud+IjziCTm5vLiy++yO23386pU6f4xje+MdhV6jVarZbrr7+ezz77bECf++c//1lZdw0W8REnTpw+EB9x4sTpA3HBiROnD8QFJ06cPhAXnDhx+kBccOLE6QNxwYkTpw/8f3dPhdOYGpEVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 170.496x63.36 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_pal=rs_pal\n",
    "fig, ax = plt.subplots(figsize=(3.7*0.8*0.8,1.1*0.8))\n",
    "for i in range(nc1):    \n",
    "    ax =plt.subplot(1,3,i+1)\n",
    "    dpl0 = np.ma.masked_array(np.squeeze(z[:,idx_sub1[label_sub1==i]]),np.squeeze(features[:,idx_sub1[label_sub1==i],0])<Th)\n",
    "    ef.plot_errorbar_vertical(ax,dist,np.transpose(dpl0),my_pal[i]) \n",
    "    plt.plot(dpl0,dist,'darkgrey',zorder=-10,alpha=0.8,linewidth=1)\n",
    "    plt.ylim([-100,100])\n",
    "    ax.set_xlim([-3,5])\n",
    "    if i==1:\n",
    "        plt.xlabel('Propagation time (zscore)')\n",
    "    if i==0:\n",
    "        plt.ylabel('Distance (${\\mu}$m)')\n",
    "    else:\n",
    "        ax.set_yticklabels([])\n",
    "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=1,\n",
    "                wspace=0.2)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 1. 2. 3. 4. 5.]\n"
     ]
    }
   ],
   "source": [
    "label_new = np.zeros(len(FR))\n",
    "for i in range(nc0):\n",
    "    label_new[idx_sub0[label_sub0==i]]=i\n",
    "for i in range(nc1):\n",
    "    label_new[idx_sub1[label_sub1==i]]=nc0+i\n",
    "print(np.unique(label_new))\n",
    "\n",
    "df = df0.copy()#pd.DataFrame(all_features,columns=['A','W','TPW','PTratio','REP','Spread','Diff_spread','W_below','W_above','PL_below','PL_above'])\n",
    "df.insert(0, column=\"cell_type0\",value = label_new)\n",
    "df.insert(0,column='depth',value=depth_info)\n",
    "df.insert(0,column='unit_id',value=unit_id)\n",
    "#df.to_csv('data/exp_features.csv') #79"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}