Examples/Exp_Model_Classifier.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# =============================================================================================================\n",
    "# Function: train experimental data, predict model data; or train model data, predict experimental data \n",
    "# The code reproduces Fig.3c,f,g,h.\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 scipy.io import loadmat\n",
    "from scipy import stats\n",
    "from sklearn.utils import resample\n",
    "from sankey import sankey\n",
    "import ecp_func as ef # import self-written sub functions  \n",
    "%matplotlib inline\n",
    "import random\n",
    "\n",
    "plt.rcParams.update({'font.size': 9})\n",
    "sns.set_style(\"whitegrid\", {'axes.grid' : False})\n",
    "apsp_pal={0:[52/255.0,58/255.0,195/255.0],1:[220/255.0,53/255.0,57/255.0]}\n",
    "cmap = sns.cubehelix_palette(light=1, as_cmap=True)\n",
    "\n",
    "colorDict = {\n",
    "    'FS':'steelblue',\n",
    "    'RS':'crimson',\n",
    "    'AP':[52/255.0,58/255.0,195/255.0],#'steelblue',\n",
    "    'SP':[220/255.0,53/255.0,57/255.0],#'crimson',\n",
    "}\n",
    "\n",
    "cre_pal={    \n",
    "    'Pvalb':'lightskyblue',\n",
    "    'Sst':'thistle',\n",
    "    'Nr5a1': 'peachpuff',\n",
    "    'Rorb':  'khaki',\n",
    "    'Scnn1a': 'lightcoral',#'#ffb3e6',\n",
    "    'RS':'crimson',\n",
    "    'FS':'steelblue'\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The setting of parameters \n",
    "random_num = 100  # The number of random_state for seperate the data into test and train data \n",
    "rep_num = 100     # The number of resampling "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "savefigdir = 'Figures/'\n",
    "df_exp = pd.read_csv('data/exp_features.csv',index_col=0) # Load experimental data \n",
    "df_model = pd.read_csv('data/model_features.csv',index_col=0)  # Load Model data\n",
    "#print(df_model.keys())\n",
    "\n",
    "#Experimental features and cell types\n",
    "one_ch_exp =  df_exp[['TPW','REP']].values\n",
    "matrix_exp = df_exp[['PL_below','PL_above']].values\n",
    "label_exp = np.array(df_exp['cell_type0'].values)\n",
    "label2_exp = np.array(df_exp['cell_type'].values)\n",
    "FS = np.where(label2_exp==0)[0]\n",
    "RS = np.where(label2_exp==1)[0]\n",
    "\n",
    "#Model features and cell types\n",
    "one_ch_model =df_model[['TPW','REP']].values\n",
    "matrix_model = df_model[['PL_below','PL_above']].values\n",
    "cell_type = np.array(df_model['cell_type'].values)  # 0-Aspiny; 1-Spiny\n",
    "SP=np.where(df_model['cell_type']==1)[0]\n",
    "AP=np.where(df_model['cell_type']==0)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# zscore plot of experimental data and model data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36\n",
      "1846\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKcAAACNCAYAAADSIWq8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfBklEQVR4nO2deVRU99nHv/fOnQFkcXBcRkZBM6AoiiIaXI8b6XFrk2I00Yi2J27VJn1TqRxTNWrfJLVNjKetPVn0mKh1eW1jbPU18tpYtYmQiAqyiSAMZlirYRsYZrvvH5N7vbPc2Zjl4rmfczwyzF1+wHee3/M8v+f3XIKmaRoiIgKEDPUARET4EMUpIlhEcYoIFlGcIoJFFKeIYBHFKSJYBCnOe/fuhXoIIgJAkOI0mUyhHoKIABCkOEVEAFGcIgJGFKeIYBHFKSJYRHGKCBZRnCKCJeji/Pzzz7F58+Zg31akDxJUcTY0NODUqVMQS0hFPCFo4rRYLPjtb3+L1157LVi3FOnjBE2c77//PrKysjBgwACn7586dQpZWVnIysrCo0ePgjUsEQFDBGubRlZWFiIjI9HT0wONRoNXXnkFq1atcnpseXk5xowZE4xhiQgYKlg3+vTTTwEA3377Ld566y1eYYrYotVq0djYCKlUCqPRCKVSCZVK5fZ4d8f1BYImziedQIhCq9UiPz8fZrOZ/V5NTQ2mTZvm9B7c42trazF16tQ+LdCgp5KGDRuGP//5z8G+bUBhRFFdXY38/HxotVq/XLexsdFGmIA1sCwtLXV7vNlsRmNjo1/GESrEJLwfCJQolEolJBKJw/dbW1udfgC4x0skEiiVSr+MI1SI4vQDgRKFSqXC1KlTnV6vuroahYWFNiJljler1X1+SgeCGK17gxCjdXc+ZaADkeLiYlRUVAAACIJgFzIkEskTIURniOL0AG6g4akYfBWrq/O0Wi2qq6sd3Aa1Wo309HTPfyA/jTXQiNG6BzjzKd2lczyNmrnCAODyPJVK5dSfdeVGeGLxhRrhi+L0AKVSidraWtZyuvMpPRWzvTAGDRrk9jzuWAiCwOjRo6FSqVirClgtKWD1S5ubm2GxWHiF5+0HL5i4FKder8fFixdx69YttLW1ITY2FmlpacjMzERERESwxhhymEDD06nPUzHbC0Ov17P+pP15XAtoPxatVovr16/DYrEAAJqamkAQBPuaub47sQstwuf1Of/yl7/g5s2byMjIgFqthkKhQEdHB8rLy5Gfn4+JEydi9erVARmUEHxOb/wwZ8d6Mp1yLRtJkqBpmg10kpOTkZqayh7ryuctLCxkrSYfrnxlofqcvOKsqKhAcnIy74nu3u8NoRanNwGQr8EScw5Jkhg8eDD0ej1aW1vZY+RyOX7wgx8AcBSffQCk1Wrx1VdfOS1FJAgCQ4YMgVqtFpTwPIE3z8kIz2Kx4Ny5czh06BBu3ryJtrY2m/efRLxJqvuSgOeeY7FYEBkZifDwcJtjuIl2V3lUxur179/f6b0IguiTwgQ8SMK//vrrqK+vR15eHnQ6HbZs2RKMcYUUb5Lq7oRjnygHAKlU6vCaCWK4MEJXqVRQKpWQSqUOrgOzbNra2gqCIACA/R+wit/dlC9U3Iqzvr4e69evR1hYGGbNmoXu7u5gjCuo2IvIm5UWvmNdrbcbjUabaxiNRqhUKpvZiBG6VqtFXl4etFotjEYjtFotzp8/j7y8PNy8edNm7Z2macjlcowePRok+fhP29TUhKtXr/ptzT9YuE0l9evXD1euXIHJZEJhYSEiIyODMa6gwfX/qqur2UBEpVJ5PBU6O9ZVioYvQmYCoPr6ekRHR9sETFx0Oh3vWFpbW1kr29PTA8Aq2sbGRrS0tAgqj+kOtytEjx49wocffoj79+8jPj4e69evx+DBg72+kV6vx5YtW9De3g6ZTIZ9+/bx+knBDIjsgw2CIDB9+nSv/4D2Ea+7QMk++c7UbN67d8+hEsmf9HY1KZi4tZxbt27FwYMHe32js2fPYtKkSXj55Zdx+vRpnDp1CuvXr+/1dXuLUqm0ESdjZbxdcnS2ysKXG+UWEHOtI3fN3B6CICCRSHrV5ExoeUx3uBUnTdP46KOPMHLkSNbRnj9/vtc3WrJkCesHmUwmUJQwFqcYX+/u3btOk9+ewDeFO5vunRUQM7iaxGJjY9lMCZcBAwa43XM1YMAAxMbGCi6P6Q63CklLS0NPTw9bEQP4Jk7GV62ursaJEydw7Ngxr68RKFJTU6FQKHxORHuzyuKsgJiBJEmMGjUKRqMRjY2NNr6lL5v+SJLE0KFDMWPGDPZ7Qk24O8Otz0nTNC5duoSamhrEx8djwYIFPt+svLwc27Ztw7vvvuuQOjl16hROnToFANiyZYvNL1SoOPMbPdnjw5cwl8vlUCgUkEqlqKysdAiEvCUmJgbjx493yCB4s2AQStyKc8eOHYiMjERqaipu374Nk8mEnTt3en2j+vp6bNiwAQcOHEB8fLzLY0O9QuQJvflDc2szuTA+pyvf01uY4pDU1FS3K01Cw+20Xltbi6NHjwIAFi1ahBUrVvh0o0OHDkGn0+HXv/41ACAzMxNr1qzx6VpCwNdqHiZfGRYWxqZ6GBhB+ipMZ6KmaRoVFRVQKBSCLvJwhltxWiwW3L17F6NHj8bdu3d9DmR27NiBHTt2+HSuEOH7Q7srFuZWD9lDkiT7HvdrT3El6tLSUigUCiQlJdlsMRayD+p2Wq+ursa7776LhoYGKJVK5OTkOF1q8yeBmNY1Gg10Oh0iIyORkJDgl2va7ynn5imZ6VShULBTqX1xhz3cyLt/VTViauvQPiIebYn++31zXRCh+6BuzaBOp8OcOXOwfPlyvPHGGzAYDMEYl1/RaDTo6OgAAHR0dKCsrAwWiwURERFOP2jV1dXo7u5GREQEwsPDXQY9Op3O6SoOM5164z9yhTny3AVITCYoSkpRs2ShjUDdCXfE388hRvMA7QnDUfujJTbvmc1m3Llzh/1ZhFpoDHhgOZcuXYp9+/YhISEBDQ0N2LJlC44fPx7QQfnbcjJidIa9QBlhMnR1daG2tpbN0VosFkgkEiQlJfV6NUcmkzn9sA+/dBmDbxexr5snTsCDzLkAbIVrpigb4Q699iUG3SoCZTCAAEADeDQq0UGgABx+Bua1J11FgoVby0kQBDsNDh06tNfpjVAQGRnJWk577AtZ7F8zpWz2VeX19fUOwiRJEnK53OOcpL0wGYtoCpPBTFGsANtHPM5uxNTWQfL9KpHEZEJMbR3aEtUYeu1LDC34BgTnegSAGM0Dp/eOulsJSf43mDJ9KlpGjLBxSYSyl8itOBcvXowVK1YgOTkZVVVVWLRoUTDG5VcSEhJYnxOwFZr9dpOIiAgbgdI0jTFjxsBgMKCmpoa1nHFxcdDpdDYCtVgsDh9epVLp1tcEHC1iU3oaqB6Dw9TdPiIeipJSB+HKq+7bCBOwWs72hOEu72Uq+BrJ/70bFZH9BDfFuxXnT3/6Uzz//POora2FSqXibWEodLhBENenpCgKZWVlbKCkVqttpnamVjMsLAxjxoyBXq+3mfbs85VMXSVN0yBJEmq1Go2NjW7FaW8RqR4DO5VzaUtUo2bJQgefszXxKUQ8fMhO5xaJBK3qkU6n9Nhv69l70Xo9ugq+hnLlizbZB6lUisLCwpBO8W7FeeHCBQBAT08PcnNzsXz5cvzkJz8J9LgCCuNj2gdKGo2GFSgAlJSU2JxH07RN0ppPcPZuPDftxAefRXRGW6LaIRBqmGVdUZNX3Udr4lPsa6eMTwHu3AF6DCDCw9Ev42kM4hSqCGWKd1ts/PHHH2P27Nk4c+YMzp07h3/+85/BGFdQsK+LtH/NrSi3f63VatHc3OxwTfsqdCYynjp1qsP1uDAWsXniBIfo3B39q6ox/NJldA1Von7WdFA9BvSv4q9+r1UMQPWiheicMQ3da3+Kish+0Gq1UKlUSE9Ph9FoFERDMI8y6lVVVYiLi0NXV5fTypi+in2gZF9IbZ8GYsSl0WjQ3t6OoUOH2lSXy+VyAEB7ezvre7a3tyM/Px9KpdJlSsnXvCbXfxx4p8RaWWWxOE1BcWlVj0SreqT1RXW1jYVUKpWsf02SZMhWktxaztWrV+PgwYPYtGkTDh8+/EQ9CSMhIQHR0dEgSRLR0dEOyXmuWMlbt0EdO46yo8fQ0dEBgiAQExNjM921traitbWVrc1kcGd9GIENvl2EkecuuLR69nB9VdJshuT7DwUTyXuKEFsmehStL168GADwyiuvBHxAwYYvUFKr1WyU3/3vL0F9cBCEwQDy6r9h2LAWlrSJIAiCd9uKszVuPvjSQ57A9VUtEglrOZ35rQRBsLsPuNYdsKbBdDodu5+eec9isYQschdGxa8A4Ebo3d3dKC8vh0qlQkJCAlr+569o/z4nSRgMkJSWwZI2ka0gGjNmDDsFms1mVFZWOlzfVX7Ym2DIHvvoHQCve0DTNJtNiIuLQ3h4OKRSKVpbW9Hc3IzGxkY0NzfbfJBCOa17Lc7m5maf9hAJHfvku9lsRl1dHeLj49Ev42m0nf0HCIMBtEwGc8pYAFZLRJKkzf8SiQSjRo1yKlDAmkc1m81sAp7xNZvS01DeMgjFZBLiyHAkosfp+c6wj97dWV2apqHVajFjxgyoVCoUFhbaWEougwcPDlkqya3POWXKFFy9epV9nZOTE9ABBRuNRoOysjLe9+vq6lAzIBaGDWthmjeHndIZ7PeKMwIF4FDB1b+qGoPO/S9i79eyrxlf89vCbnxcPwsFGhX+cSEWVffD/PhTOocpSOHuvSdJkl2qlUgkAS/ycYVby6lWq3H+/HlUVFRg/fr1PtcaGo1G5OTkoKWlBYmJidi9e7fL1Eow4OY53WFJm2gjSneMHTsWsbGxbImczQrQnVLov5+KGV+z1KKG0WwViMlEolYThsSnPLeevUGlUiEpKQn19fWIi4vr1ZYVf+LWcspkMuzduxcmk6lXVjMvLw9JSUk4fvw4ZDIZrl+/7vO1/IWr/d+9gfnQ6XQ6TJs2DUqlEv01jkFP+4h4mL+3rilkNaQSa26RoiwYkRB4YTKpL61Wi8rKSrS3t7PuSHp6uvCXL5kBbtq0CRcuXMC2bdt8ulFRURGeeeYZAMD06dNx48YNTJ8+3adr+QtXBSG+8rfPunGryIi0CVIsfQ4oLrHgfq0K6sQfYmBpGQiDkQ16uMHMgBHxWEK2olYThhEJPUGxmswKl310zkz33jz/KBDwivPmzZuYNGkS3n77bfZ7CxcuxMKFCwEAN27cwOTJkz2+UWdnJ5t2iYiIQFdXl8379hvcAk17e7tDzyJPuVVkRGmZESljpUib8Pgaf/usG+cuWEWlre9BY5MZxSUmGAxAeJgc//VsLlRV/0BbwuNImhvMJCI4omRobm522qJGr9c7bF8OxTImrzjr6upw5MgRpKam2vTnLCsrQ1FREebOneuVOKOiolhBdnV1ISoqyub9F154AS+88AIAaz1nIGlvb8eDBw988p9vFRnx/kEdDAbg2lcGbFwbibQJUtwqMuLyFVth3Sm1ChMA9D0WXH8Yi6fnOxZzhAomh6lWq236hAJwqANwV6kUiO0evOJ87rnnsGTJEly9ehU3btxAW1sb5HI5Jk2ahNWrV3ttdcaPH4+CggJMnjwZ169fx+zZs3s9eF/p7Oz0WpiMtWz5j4UVnMEAlJZZm3IxguViNAJSioDRRCM8jMSMaQpYTI51oMGGu8tTKpVCpVJh2rRp7DTuLA3makNcoPrKu/Q5KYrCvHnzMG/evF7faMGCBcjNzcWLL76IkSNHYubMmb2+pi9UVFR43dKFay0pCqBIC0wWEjLKgpSxUpSWGR2ECQAWCzAtQ464uHBkTJFj9iwFtNpINDY2Qq/Xh6TrG9NoQavVgqZp3Lt3DwqFgu1Ows15Ao/30ruyiIHa7hG0FSKpVIp9+/YF63ZO8UaYXL+SKz6TCZiISgwxP8Q41CDVMhMYm4JrXxlgMABSqdUqmUxAeBiJHz+rxOxZCva63BY1fI9uCSQ0TdtYbncd8FJSUvzWA99b3D6wgCRJyGQyv9ws1HgjTK5f+YP5YZDJrNO4jDQh05CPyXQZYAFMpQOQtmoiNq6NhOZBOEYM14MGUFpmwpzZtsK0hxHqZ599FrSNgzRNQ6/XQyKROBWTtw9n8PUcT+AV5yeffIIzZ86AIAhs3rwZmZmZfrlhKKEoykGgziJvrqU0GIDubhob10aitMyIceFNSL9YBRjALmXeKjLibiWN9EkUOjvBXi91nONzK53x1FNPOe0AEijCw8ORkpLCKyZvepP25hx38Irz4sWLOHPmDLq7u58YcSYnJ9tM7XyRd8pYKTtNy2RghWsV72gYRqyFpLQM5pSxKCRT2Gt8caWbndKvfWXAoEFmxMW5H5dCwW9de0tYWBhommYtM7N1JBBi8je8K0QURYEgCPTr1y/k0aU/4ba2treQTOQNAMmjKKSOo1jBcrGkTYRx1UpY0ibaXMNotAqTuV5xiWdTdSB9TqbQhGHUqFFsQwVn/eqFBK84zWYzOjs70dHRwX7N/OvLcKfPlLFSMO60TAZERBB474+d+POHOhSXmFBR6d5H5V5DKiUgpQj2esPiPPtd2T+62hOLplQq3T7RhCAIhIeH24iT6SsfiOfD+xvepgrZ2dlOu1UQBIEjR44EdFCB6jJn3zABeOxzRkQQyPtnj0NKaP4cGVat6OfyundKaTzQ9kPGFDm02m9RfKeHdQUoiuIVEeNiUBSF6OhoGx+wuLiYbWgLWKdno9HIbk1mconFxcXQaDQgSRLd3d1sGkgulyMlJQUAHFrONDY29oluc7w+J9NZzp6+2FSBgREmees26zOmpU1E2gQpjp3ochAm4286IywsDNHR0VAqlRg37vH3S0oakDrusZj5MgRc39dkMuG7775DREQE+z1nDW2drcIoFAp2pyRTGGz/3CFnkXRf6DbHO61z17c/+ugj9uu+vC2YIAiQt25D9sFBUF/8C7IPDoK8dRuA7fRMUUDqOAo/m9eEp++cZo9hkMlkSEpKcvijajQah3vydeVzJlqLxcJuUQbA7oYErA9WAGDzmhGr/QO3+B6MwHzfm0fZhBJey8nd9nrt2jWsW7cOgO+9I0PFlWsPUfBNKzKmyFFXZ0D5WWCCMRGTUWaz5SJtgpRNF0VEEOiuawH1+f+BMhZD8u+v2CJjmUyGUaNG2UzJycnJTmtDXU3pztJaDNxSPvulQW5/I+Y1X87S1bJiX4jWecXJLQR2tj22L3Dl2kNsf+Mu9D0WnDnbAAAwmUbgCrUSr5qOI11axW65AMBG5dbUUAyu0cvxKmHCZEMZBjU2YdD387f9lFxRUeHg7pAk6TJgsU9rceFumrNfGuT2aDKbzTAajbwJcKF3kXMHrziZCJ2maVgsFvbrvpRWKvimFfoeq2i4GjBAhqLhczD+2TkO1e026SVChjtEEtKlVeiX8TR7jL2gTCYToqOjXe6BdwZXvHz9Q+2XBrk9mhhLyWcF+1onY3t4xUmSJH72s585fM1NeQidjCly/ONcE/Q9FjCun8lkDXTG/DgFlgmOwY5NAp6yYOzYMJjnboBWOQQGu/Y0DBRF2TQL86VBLd/xzpYGPd1GEahlxWDhtj+nPQ8ePMDw4Y6dy/yJP1NJV649xMW8GoxJtn6onBUJ22O/pMnXRxNw7VeK9A5ecV65cgV/+tOfEBUVhd27d2PIkCHYv38/vvjiC1y8eDGgg/J3nrO0tNSnQI7pAmLf0IvLOG4eScSv8E7r+/fvx8GDB9HS0oI9e/agtbUVM2fOxNmzZ3260f79+5Gfnw8A2L59e1D+qJWVlTAYDD4FcfHx8YiJiXGaHmKw7+0p4l94xRkVFQWFQgGFQoH79+/jd7/7nVfbMrjU1NSgoqICJ0+ehEajwZ49e3Do0CGfB+0JjDABsFXfBEF4tIhAkiRiYmLQ2NjIuwGOr5+8iP/wKJWkUql8FiYADBs2DHv37gVgjWx93VjmDfY+Ik3T7HKeq2kaeLwK5kyYjEUVCTy84tRoNHj77bdB0zTq6upsdmG62x58+vRpnD59mn09duxY7Nq1CzqdDtu3b0dubq7DOf7efWkfxDAF0+3t7W7PZVZ1oqOjbR5kFR0dLQoziPAGRF9//TXvSU8//TTve3y0tbVh06ZNyM7Odvv8TH8FRMzUzqzqAMDdu3dhNBpdnsf1h5mpnVlHFwkeXqeSfMFsNmPVqlVYt26dR5vlAvnsS3eR+8CBA0URCgS37Wj8waVLl1BVVYXDhw8jOzvb564h/sB+vzyXiIgI6HQ6lJSUoKSkxGWkLhJ4gmI5vSWQlpMbxXMZOHAgdDqdQ72ns47HIsHhiWgey608crXbkU+YgPVpIfbCBALX7EvEPUGZ1gMJU3l0+m8N2P7GXVy59pD3WFfbb3U6ndOkuicFHCKBoc+Lk1t5pO+xoOCbVt5jXe2/j4yMhFqtthGoOKWHlj4/rXMrj8LDSGRMkfMey7TD5m6TBWBTRSSu+giHJyIg8tTnFOlbPBHiFHky6fM+p8iTiyhOEcEiilNEsIjiFBEsojhFBIsoThHBIopTRLCI4hQRLEEV54MHDzB3rnCewyMibIImTpPJhHfeeeeJefiBSOAJmjgPHDiAlStXIjw8PFi3FOnjBKQqyX73ZVhYGCZNmoSMjAzec7i7L3Nzcz1+xOCjR48wYMCA3g04wIhj5IeiKCQlJTl9LyiFHzt37sS9e/dAURRKSkqwaNEivPnmm365dlZWFj799FO/XCtQiGP0jaDUc+7Zs4f9+tlnn/WbMEWebMRUkohgCbo4fW0ExgfzGGwhI47RNwRZbCwiAojTuoiA6fMb3ELR99MTjEYjcnJy0NLSgsTEROzevVtwD3vQ6/XYsmUL2tvbIZPJsG/fPvTv3z/Uw2Lp05aT2/dz7969eO+990I9JJa8vDwkJSXh+PHjkMlkuH79eqiH5MDZs2cxadIkHD16FAsWLGDzzEKhT1vOUPT99JSioiI888wzAIDp06fjxo0bmD59eohHZcuSJUvY7dHMM5WEhLBG4wZf+n6Gis7OTrZbSEREBLq6ukI8IkeY8VVXV+PEiRM4duxYiEdkS58S57Jly7Bs2TKb7zF9P9esWYOJEyeGZmBOiIqKYgXZ1dXlsrtdKCkvL8e2bdvw3nvvCa4xbp/2Oc1mMzZu3IiXX37ZbUPaYDN+/HgUFBQAAK5fv44JEyaEeESO1NfXY+vWrfjDH/4gyE4nfTrPefHiRWzfvp19DtCwYcNs2oOHEqPRiNzcXNTX12PkyJF46623BBet/+Y3v8Hly5fZh2dlZmZizZo1IR7VY/q0OEWebPr0tC7yZCOKU0SwiOIUESyiOEUEiyhOEcEiipOHN998E9nZ2ZgxYwaysrKQnZ2Nr7/+GtnZ2VixYgWys7Px4osv4uOPP8aXX36JN954gz133rx5OHHiBADr3pyXXnrJ4fq7du2C2WwO2s/D5Y9//COamppCcm+voEVckpubS5eVlbGvV61aRbe1tdE0TdMmk4n+0Y9+ROt0Ovr555+naZqma2tr6U2bNtGvvvoqTdM0fenSJXr//v021/zXv/5Fnzx5Mkg/gSONjY30zp07Q3Z/T+lTy5dCg3kupkQigVQqhV6vR35+PhYtWoTjx4+Dpmncvn3boeDjr3/9K3bt2oX79++zFler1WLevHnYvn07e9wvf/lLNDc3gyRJvPPOO+jXrx9ycnLw6NEjDBw4EPv378fp06dx5swZEASBDRs2IDMzE8uXL0d0dDRmzpyJkSNH4v333wdN01i9ejUWL16MIUOGQKPRoKenB2FhYcH7hXmJKE4f2LBhAyiKAkmS+PnPf85ufS4pKUFBQQFef/11FBYWoqKiAmVlZdi8ebPN+TU1Nezjwo8ePYqGhgb86le/wi9+8Qv2mI6ODtTV1eGTTz5BZWUlOjo68Pnnn2PWrFl46aWXcPr0aRQXF+Pvf/87Tp48CYPBgOXLl2Pu3Ln4z3/+gwMHDkChUGDp0qU4ceIEJBIJsrOzMX/+fISHh0OlUqG8vFxQ9Qj2iOL0gQ8++MChSGLKlCkoLi5mrVpGRgZu3boFiUTi0EhCIpGwX+t0OuTk5GDPnj2Ijo5mvx8dHY21a9fitddeA0VR2Lp1K2pqavDcc88BsBbBFBUVYfTo0aAoChRFYciQIXj06BEiIyMxaNAgPHz4EPX19Vi3bh0Aq+BbWlowfPhwKBQKtLW1Beg35B/EgMhPpKenIy8vD4mJiQCsT1Y+e/Ysxo8f73AsI06LxYKtW7di3bp1DoUXzc3N0Gg0+PDDD7F06VKcPHkSw4cPZ5tNHD58GFqtFpWVlTCbzejq6kJDQwPkcjlboymXy5GQkIDDhw/jyJEj+OEPf4jBgwcDsAo1NjY2YL8PfyBaTj8RFRUFo9GIqVOnAgBiY2Oh1+udPv572LBhePjwIQoKCnDjxg10dnbi0KFDGDVqFHbs2AEAGDRoEKqrq7Fs2TKEh4dj+/btUKlU2Lp1K86fPw+FQoHf//73+O6777BixQqYTCbk5ubaFFxLJBKsXbsWq1atgsFgQGZmJutj1tbWsgUzQkUs/AgBly9fRkNDA1auXBmS+zc1NeHAgQM2zS6EiDith4A5c+agvLw8ZHnOkydPYuPGjSG5tzeIllNEsIiWU0SwiOIUESyiOEUEiyhOEcEiilNEsIjiFBEs/w/BeRzdIWYTMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 158.4x129.6 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Upsample minority class\n",
    "SP_upsampled = resample(SP, \n",
    "                        replace=True,     # sample with replacement\n",
    "                        n_samples=len(AP)-len(SP),    # to match majority class\n",
    "                        random_state=27) # reproducible results\n",
    "\n",
    "# Combine majority class with upsampled minority class\n",
    "model_upsampled = np.concatenate([np.concatenate([SP, AP]),SP_upsampled])\n",
    "print(len(model_upsampled))\n",
    "X_model = stats.zscore(df_model[['TPW','REP']].iloc[model_upsampled].values)\n",
    "Y_model = df_model['cell_type'].iloc[model_upsampled].values\n",
    "\n",
    "\n",
    "# Upsample minority class\n",
    "FS_upsampled = resample(FS, \n",
    "                        replace=True,     # sample with replacement\n",
    "                        n_samples=len(RS)-len(FS),    # to match majority class\n",
    "                        random_state=27) # reproducible results\n",
    "# Combine majority class with upsampled minority class\n",
    "exp_upsampled = np.concatenate([np.concatenate([RS, FS]),FS_upsampled])\n",
    "print(len(exp_upsampled))\n",
    "X_exp = stats.zscore(df_exp[['TPW','REP']].iloc[exp_upsampled].values)\n",
    "Y_exp = df_exp['cell_type'].iloc[exp_upsampled].values\n",
    "\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(2.2,1.8))\n",
    "ax=plt.subplot(111)\n",
    "ax.scatter(X_exp[Y_exp==1, 0], X_exp[Y_exp==1, 1], color=\"darkgrey\",s=10,label='RS')\n",
    "ax.scatter(X_exp[Y_exp==0, 0], X_exp[Y_exp==0, 1], color=\"lightgrey\",s=10,label='FS')\n",
    "ax.scatter(X_model[Y_model==1, 0], X_model[Y_model==1, 1], color=[220/255.0,53/255.0,57/255.0],s=10,label='SP')\n",
    "ax.scatter(X_model[Y_model==0, 0], X_model[Y_model==0, 1], color=[52/255.0,58/255.0,195/255.0],s=10,label='AP')\n",
    "plt.xlabel('TPW (z score)')\n",
    "plt.ylabel('REP (z score)')\n",
    "plt.xlim(-2.5,3.5)\n",
    "plt.ylim(-4,5)\n",
    "ef.simpleaxis(ax)\n",
    "#fig.savefig(savefigdir+'APSP_model_aa_EAP.pdf',bbox_inches = 'tight')\n",
    "#fig.savefig(savefigdir+'APSP_model_aa_EAP.eps',bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model based classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The initilized parameters for model-based classifier\n",
    "accuracy_XY_me = 0\n",
    "conf_matrix_XY_me = np.array([[0, 0],[0, 0]])\n",
    "w_me=[]\n",
    "Coef_name_me=[]\n",
    "y_exp_dictlist = [dict((j,0) for j in range(2)) for x in range(len(df_exp))]\n",
    "random.seed(100)\n",
    "for ii in range(rep_num):\n",
    "    # Upsample minority class\n",
    "    SP_upsampled = resample(SP, \n",
    "                            replace=True,       # sample with replacement\n",
    "                            n_samples=len(AP)-len(SP),    # to match majority class\n",
    "                            random_state=27+ii) # reproducible results\n",
    "    # Combine majority class with upsampled minority class\n",
    "    model_upsampled = np.concatenate([range(len(df_model)),SP_upsampled])\n",
    "    #print(len(model_upsampled))\n",
    "    X_model = stats.zscore(df_model[['TPW','REP']].iloc[model_upsampled].values)\n",
    "    Y_model = df_model['cell_type'].iloc[model_upsampled].values\n",
    "\n",
    "\n",
    "    # Upsample minority class\n",
    "    FS_upsampled = resample(FS, \n",
    "                            replace=True,        # sample with replacement\n",
    "                            n_samples=len(RS)-len(FS),    # to match majority class\n",
    "                            random_state=27+ii)  # reproducible results\n",
    "    # Combine majority class with upsampled minority class\n",
    "    exp_upsampled = np.concatenate([range(len(df_exp)),FS_upsampled])\n",
    "    #print(len(exp_upsampled))\n",
    "    X_exp = stats.zscore(df_exp[['TPW','REP']].iloc[exp_upsampled].values)\n",
    "    Y_exp = df_exp['cell_type'].iloc[exp_upsampled].values\n",
    "    \n",
    "    # Model-based classifier: Training model, predict exp\n",
    "    accuracy_XY1,conf_matrix_XY1,accuracy_newXY1,conf_matrix_newXY1,w1,Coef_name1,dict_y_pred1=ef.classifier(X_model,Y_model,X_exp,Y_exp,'svm',random_num)\n",
    "    accuracy_XY_me = accuracy_XY_me + accuracy_XY1\n",
    "    conf_matrix_XY_me = conf_matrix_XY_me + conf_matrix_XY1\n",
    "    w_me.append(w1)\n",
    "    Coef_name_me.append(Coef_name1)\n",
    "    for k in range(len(df_exp)):\n",
    "        y_exp_dictlist[k][0]+=dict_y_pred1[k][0]\n",
    "        y_exp_dictlist[k][1]+=dict_y_pred1[k][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy of test data set:81.67111111111099%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAABoCAYAAADSHkMnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAANgUlEQVR4nO3df1hUdb7A8fc58wNFdFMv4w9IvCooCqJL4AYm66oxumbYYnpFXNNupjczDTVjs7Uer7/KupfUfcofCVZqz2Jpu6jl6hWLII1AIZUrm2ImyxXUBUmYH/cP5CgwMzDCzJni+3qe7/PMfOfMfD/Ah+/5zplzPiNZrVYrggDIagcgeA6RDIJCJIOgEMkgKEQyCAqRDIJCq3YAQuusXr2aiIgIYmJiSEpKoqysjAEDBrBy5UpMJlOTPkmS7L6WmBl+ompra1m8eDGffvopAIcOHSIwMJD3338fvV5PVlaWzT5HRDL8RNXW1hIfH8/kyZMByMvLY8SIEQBERUVx4sQJm32OuHQ3kfWf21z58k6JTJqudggNaPQdmvQNDYhRbievm8/u3buV+1OnTmXq1KnKfW9vb6Kiojh58iQAlZWVdOrUCYCOHTty8+ZNm32OiDWDB5FljXK78R+/OT4+Psof++bNm/j4+GCxWJr0ORz/HmIWXEQna5XmrNDQULKzswHIysoiLCzMZp8jIhk8iFbWKs1ZRqOR8+fPM23aNKqqqhg5cqTNPkckV35qKdYM9tlaM4wa9Khy+9iZj90ZDiDWDB5FI6k7UYtk8CB6jU7V8UUyeBDNXe8m1CCSwYPoNer+OUQyeJB7eRfRpuOrOrrQgE7sJoR6WpEMQj2dRiSDcJvaM4Pqh6O7DerLwN+NASBgbCRDEicwZMZ4dJ06NtjuvgH3EzrrEUJmTsTb0M1uX2vV1NTw7HOLmDHz96xavYbGB2i///57npjzpHLfarWycNFivj1zptVj6zUapalB1WTQd+5Ej+EDAbivvz8SEgVpf6Uk8xs6dOvSYNv7Rw6j4L0Mzu09Qp+YX9rta62MAwcJCQlhZ+oOqqurKSgsVB77MjubxUuWcu36DQBu3LjB3HnzOXW6oE3G1siy0tRgd9SysjI2bNjAxo0bqaysdMngAWMjuXi07vP4Ln16Yq6pJXhaLL5D+lN5uazBtqd2fIKl1oS+szemW7V2+1rr0UmPMOeJWZjNZsrLy/G5fT4A1B0UeudPm5X71dXVPLvgGX4VGdkmY+s1WqWpwW4yLFu2jEGDBmEwGFi1alWbD+wXHUZZXhGmmz8CoO3ohb6zN9/uOsitG5UYwoIaPsFqxRAWRPDUh6koumi/rw3Isszk+ClUVFzD19dX6Y+IeIAuXe7MWD169CBkyJA2G1cry0pTg91RzWYzEyZMYMqUKVy+fLnNB+4WFEDvB0MJjPs1ne/vQfeBfbn+97pxrhV/j7eha5Pn/CPvHCff2o1fdBiyXmu3r7UkSWLf3nQemxzHO1vd98mrViMrTQ12R7174eSKT7lPbd9H4XsZFH10lH+WlHI+43O69O0FgI+fL9VXr9/ZWJIInhaLJEtYTGasFgtA0z5L6+PctWcP+z/5C1B3qpg79986WaM0Ndj9V7pw4QKrV6/GarUqt+stX768zQMpP3uB+/r5ETJzIqbqHzn30VE6dO2MYdhALh45wf8VFjMkcSJYLfzw5WksNaamfSZzq+MwPvwwS19Yzp/T0+nUyYd5Tz/FaxveIGnxojb4KR2TZfunsbuD3ZNbcnJyGtwvLy/nypUrpKens2/fvha9uDi5xT5bJ7esi1up3F760cvuDAdwMDNERkZiNps5cuQI6enpnDx5ktmzZ7Np0yZ3xteuqLVWqGd39DVr1jB9+nTy8vJYsGABoaGhzJ07F39/f3fG166ovYC0OzOUlJRgMBjo06cPvXr1cnhZltA21HpLWc/u6Bs3buSVV16hpqaGZ555hoKCAtLT07lx44Y742tXdFpZaWpwOGrXrl1JSEhg586d7Nmzh9LSUhITE90VW7ujkSWlqaHFKejv78+8efP4+GP3n8LdXui0GqWpQXyE7UHUPs4gksGDqLVWqCeSwYNoNGJmEG5T+6CTSAYPImYGQaFxYmaoqalh9OjR9OvXD4CkpCS2bdvW4vpNtqh+DqRwh1YjKa05RUVFGI1G0tLSSEtL49KlS07Vb7I5/r0E3VKD4yNc+fJOOb/rE7VDaCBoZnyTPs1d7yZ2797tsIzP2bNnyc/PJyEhgZCQEKxWK+PGjQPu1G+KiopyKiaxm/Agd68Zmivj07t3b5577jmio6N59dVXOXLkCHFxcUDL6jfZIpLBg8hOrBnCwsLQ6eou4R81ahQajcap+k02x3f6GYLLaHWy0prz5ptvcuDAAQBOnDjhdP0mW0QyeBBZlpTWnDlz5vDhhx+SmJjI9evXiY2Ndap+ky1iN+FBNC2YEeoZDAZ27NjRoG/Dhg2tGl8kgwdx5jiDK4hk8CCy+KBKqCcORwsKMTMICpEMgkKj0ulu9UQyeBBJrBmEerJOzAy8seVdhg8ZTNQDw5k0ex59/XsDsHDO7wke0F/Z7mhWDql/3ovFYiV+QiwTx4622eesWpOJ9R/t5lpVFf179mb2GCOz31qPf/e62gxzxo5nQC8/ZfucojPsyvwbsiwz3ziJfj3r4q01m0ja/icWTnxM6XNGu14zmEwm/vhGCqfOFDF8yGCKL5QwZuSDLJk7x+b273ywh63rVqHVapixcAnjR4+y2adxsiZSZuEpAnv5MyU6hv/+JJ2/l/7AyOAQ5sY+YnP7DzL/xuoZT/LPH2+yKWMfL0+dCcB7/3MYSyvKFzjzQZUrqJoMtSYTk8aNoY9f3X/R/353gcJzRTz1wksED+jPwtkzke+65CzllZfo0MELk9mM1WpFlmWbfc76zdDhmC0WzBYL16oqKS79gXOXL/FC6jsM6OXH7LFG5Lsqvr8262k0ssx3ZVfo5FV3NfXXxUV01Ovp17PXPf8+ZJUXkA5/c66u69SxQwcihw1V7vc0+PL0jH/j7TWvYjKbOXTs8wbbd7vvFwBseHsbk8b9BkmSbPbdC1mSeHZLCjeqb9LNpzMzYsaxZua/Y7aYOVZwqsG2GlnmYO5X/HHXDkYEBXO9qoqDuV8RHxVj59VbRtJISlODw2RwdV2nxkKCAgkfGgJAdPhwii82rNNksVhYs+ltdDodCZMn2e27F5IksfGphYwLCyfvQjFD+/4rAOH9g7hYVtpk+9jhEby7YBl7Pj/CscJ8/nGtgj+8t42vzxeR8pe9VNfccjoGWadVmhocJoOr6zo1tnnnBxw+XnfuXm7BtwzoG9Dg8bfe3YmPtzeLnpzlsM9ZGSezOXr6GwC8dHokJI5/exqAgpLv6GvoqWxrtph56f3tmMxmdFoNsiTz8LAHeGPOf7A68Ul+2T+QBb+dTEe9l9NxSFpZaWpwmIKuruvUWOJjj7Li9f9i74FDBPj5MSb6QUou/8Deg5+REPcIu/b/ldCBQTy9fAUALy9a0KTvtT8sa1CuryWig0N5/eM9HPrmBN56L+aPf5TX933Iga+/wq/7vxAdHMLl8qsczP2KJ8YYGTU4lGWpbyNLMpN/9RBeurb50hBZ5XLBDr+jKiYmBqPRiNVq5eDBgxiNRuWxltR1un7uVLPbuEvpl2fVDqEBWyfElmXfWSP5joh2ZzhAMzPD+vXrKSwsJDg4mFGjRnHp0iV0Oh1+fn6OnibcI7XWCsr4jh7MyckhOzub8PBwDhw4QF5eHqWlpRw+fNhd8bUrkkZWmhocpmJubi5bt26lpqaGzz77jGPHjqHX65k+3bMqp/1cyG209rhXLZqXcnJyGDp0KHq9HqirmSy0PUnlmk4OkyEwMJC1a9eSlZXFkiVLqKioYMuWLQwePNhd8bUrkifPDEuXLiUzM5NJkyYRHBzMuXPn6N69OwsXLnRXfO2K2oejHSaDLMvExNw5xBoUFERQUJCDZwitIWnFt9cJt4lkEBSy+JJToZ6YGQSFLL4YXajnzMxQW1tLUlJSq8r2NCauwvYgkkantOYcOnSo1WV7GhPJ4EFknV5pzcnLy2PEiBHAnbI9reXS3cQvgkJd+fJO8aRY7JGkltd0qqyspNPt8zbutWxPY2LN4KGaq+nk4+PT6rI9jYndxE9UW5TtaUwkw0+U0Whsddmexhye9ia0L2JmEBQiGQSFx76bmD9/Pv7+/rz44otkZ2ezZMkSAgICkCQJWZZZu3YtPXr0cGkMZrOZdevWUVRUhNVqJTw8nIiICFVicQePnBnKy8uRZZnc3Fxqa2sBiI2NJS0tjdTUVKZMmcLOnTtdHkdmZiaSJLFt2za2b9/O6dOnKSwsVCUWd/DIZMjIyOChhx4iMjKSo0ePNnm8qqoKLy/nr1hylsFg4IsvviAzM5Nbt26RkpJCcHCwKrG4g0e+m0hISCAlJYUrV66QkpLCrFmzGkzNvXv3Jjk5mc6dO7s8lszMTFJTU8nPz2fMmDGMHz+e5ORkVWJxNY9LhpKSEh5//HHlpNv8/HxWrlxJbm4uycnJbo3l7Nmz+Pr60q1bN27dusWKFSsYNmwYxcXFbo/FHTxuAbl//36ef/554uPrLj/bvHkzpaVNr4J2h+PHj3P16lWWLl2Kl5cXffr0USq5/xx53JohIyODsWPHKvcnTpxIenq6KrEkJiZy7do14uLimDZtGhUVFT/rSws9bjchqMfjZgZBPSIZBIVIBkEhkkFQiGQQFCIZBIVIBkEhkkFQ/D9B/97iLOhGGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 129.6x93.6 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Confusion Matrix\n",
    "df_cm_XY = pd.DataFrame(conf_matrix_XY_me/rep_num, \n",
    "  index = [ 'AP', 'SP'],\n",
    "  columns = ['AP', 'SP'])\n",
    "print('The accuracy of test data set:'+str(accuracy_XY_me/rep_num*100)+'%')\n",
    "fig, ax = plt.subplots(figsize=(1.8,1.3))\n",
    "ax.set_aspect(1)\n",
    "res = sns.heatmap(df_cm_XY, annot=True, vmin=0.0, vmax=100, fmt='.2f', cmap=cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '${\\\\beta}$ Coeff')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGoAAABlCAYAAACld3ZCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJI0lEQVR4nO2dbUiTXx/Hv+rmQ/M5U4skDI16oYItDUGCHsWw7my2jAYRPUippeVTliZZRpYmcSdBPkAorBdSQUEqQhQ5Y2KKqTTyr1n0Ql3o1KabO/eL7v9oqe2aumsdPZ83Y9eunfOdn13HXdf2O8eOEELA+Ouxt3UABjeYKEpgoiiBiaIEJooSmChKoE6USqWydQSbQJ0ovV5v6wg2gTpRNKBQKJCeng6FQrFobfIqqqioCI2NjSbb3rx5A4lEAqlUio6ODj7jWI3q6mq0t7ejurp60drkRZROp0N6ejoaGhpmPFZWVoaqqircu3cPJSUlfMSxOhMTEya3i4Fg0Vr6AzqdDhKJBIGBgSbbNRoNXFxc4ObmBjc3N4yNjUGv10Mg4CUWVfByRK1YsQJRUVEzto+NjUEkEhnvOzk5QavV8hHJYvQGg03btulbVyQSmQwPk5OTcHFxmbGfXC6HXC4HAFy4cIG3fL8isLfHf9/9w2nfEa3eeMvlOWcjAs3uY1NR7u7umJiYgEajgVarhYuLCxwcHGbsJ5VKIZVKAQDd3d18x/wrsImo/v5+yOVyZGZmIi0tDcePH4fBYMDly5dtEWfRETg5mdwuBna0fXHY3d2NTZs22aRvrkPft54OqF6/RHD0HqzeGGp2/79+6FuqrN4YykmQJbArE5TARFECE0UJTBQlMFGUwERRAhNFCUwUJTBRlMBEUYJZUTdv3gQAPH/+3OphGHNj9lqfUqlEYWEhGhsbZ/ymIScnx2rBGKaYFVVRUQGVSoXOzk7s2LGDj0yMWTArqry8HNnZ2ZDJZIiIiOAjE2MW2NBHCRYNfTt37gRl3zMuGcx+6vPw8IBYLMbDhw/R2tqKhoYGjI+PY82aNXzkY/wfzudR2dnZCAgIwIcPH7By5UpkZWVZMxfjNziLGh0dxd69eyEQCBAaGgp7e3auzCec/9o+Pj6oqanB+Pg4nj17Bi8vL2vmYvwGZ1FFRUUwGAwIDQ2FWq02XrFg8ANnUYQQTE9Pw2AwYHp6GkKh0Jq5GL/BWVRmZiYcHBwQHx8PgUCAjIwMa+Zi/Abn3/V9//4dMpkMABAWFjZrCc1c6HQ6XLx4EYODgwgKCkJBQQHs7OwAAKWlpXj9+jVEIhE2bNiAK1euWPgSlgecjyihUIimpiYMDw+jqakJThb8XLe+vh7BwcGora2Fo6MjmpubjY99/PgRFRUVePToEZP0BziLKi4uRnNzM3JycvD27VsUFRVx7qS9vR2RkZEAgKioKCiVSuNj/f39yM3NhUwmWzIVh9bArKje3l4kJCTA0dERubm5SEhIQFtbG378+MG5k1/roFxcXIylNoQQxMXFoaSkBEVFRcjPz5/ny1j6mBVVUFCAa9euwc3NDQCwa9cu5OfnIy8vj3Mnrq6uJuWSrq6uxsdkMhmcnZ2xdu1aCAQCTE1NzXi+XC5HfHw84uPjoVarOfe7lDArymAwYOPGjSbbQkNDYbCgAi8kJAQtLS0AgObmZoSFhQEARkZGkJiYiOnpaQwPD4MQAkdHxxnPl0qlqKurQ11dHby9vTn3u5QwK8rBwQFjY2Mm2zQaDaanpzl3EhMTg0+fPuHw4cMYHx9HQEAAbt26BU9PT2ORWnJyMnJzcy1/BcsEs/VRSqUSxcXF2L17N1avXo0vX76gqakJqamps9blWhsa6qMshUt9lNkjSiwWo7y8HB4eHhgYGICvry/u379vE0nLGU4nvN7e3pBIJNbOYjUUCgUeP36MQ4cOYevWrbaOMy+WRcVhdXU1VCoVJiYmqBVlduibmprCwMAAH1mshjVmUuEbs0fUiRMnsGrVKnz79g0xMTGoqqqCv78/8vPzZ3xs5xO9wQCBlb68tGbb88WsKEII7ty5g8+fPyMuLg6vXr3C1NQUsrOzUVlZyUfGWbH1BB18Y/ZtIxKJ8OLFC/j5+SE1NRWenp7w9fW16BISY+GYPaJu376Nqqoq1NbWYmpqCiqVCh4eHtBoNBgdHYW7uzsfOReENSbo4BuzolxdXZGSkoKUlBQQQtDX14fu7m4IhUKkpaVheHgYT5484SHq/Nm06z/GCTpoxaKP53Z2dggMDERgYCBiY2OtlWnRscYEHXzzd320YcwJE0UJTBQlMFGUwERRAhNFCUwUJTBRlMBEUQITRQlMFCUwUZTARFECE0UJvC1LdO7cORw5cgR5eXkmc1UsxfWjrAEvov5UH7UU14+yBryImqs+6tf1o3x9fY3rRzFmwououeqjaFo/ytbw8kvZueqj5rN+VFZWlnFpou1u3DOo1WrOJTtzLX3ER38CgQDBwcEz9uFF1L/1UWKxGM3Nzdi2bRuA+a0fNV/i4+NRV1e3oDZs2R8vQ99c9VEAjOtHJSUlIT09nY84dEKWCQcOHKC6v2VzwrvQodPW/VG3IttyZdkcUbRDfSHb9evX0dPTg97eXvj5+UEkEuHdu3cIDw+HQCCATqdDTEwMgoODUV9fj4KCAgDA9u3bcfLkSSQmJkKtViMlJQU1NTWc+mxpaUFGRgbWrVsH4OdchufPn0dlZaXJgs8+Pj4oLS2FTCYzbv83z7Fjxyx7oYv6H8+GZGVlka6uLkIIIUePHiUjIyOEEEL0ej3Zt28fGR8fJxKJhBBCSF9fHzlz5gxJTU0lhBDS2NhI7t69y7kvhUJBCgsLjfcHBwfJwYMHTfr9lbnyWMKSH/omJycB/JyGQSgUQqvVQqFQIDY2FkNDQyCE4P379wuaKnxoaMjkCgvXPJZA/dA3F6dPn4ZAIIC9vT2Sk5Ph5OSE8PBwdHZ2oqWlBZcuXUJrayt6enrQ1dWFs2fPWtT+y5cv0d3dja9fv2L9+vW4evUq8vLyjP0CQHR0NE6dOjVnHktYsqIePHgwo3Zry5Yt6OjogFqtho+PDyIjI9HW1gYHBwc4Oztb1P6ePXuQm5sLpVKJGzduwN/ff85+/7SdK0t+6PuVzZs3o76+HkFBQQCAiIgIPH36FCEhIfNuUywWIzo6GmVlZYsVc1aW7BE1G66urtDpdMYpDLy8vKDVahe8lEVSUhL2798PrVZrMvQBP6dOWAzYCS8lLKuhj2aYKEpgoiiBiaIEJooSmChKYKIogYmihP8Bz3WH6Bm5++kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 72x93.6 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Beta Coefficient\n",
    "df=pd.DataFrame(np.mean(np.array(w_me),axis=0),columns=['Beta'])\n",
    "df.insert(0, column=\"Coef_name\",value = Coef_name1)\n",
    "fig, ax = plt.subplots(figsize=(1.0,1.3))\n",
    "ax = sns.barplot(x=\"Coef_name\", y=\"Beta\", data=df,color='skyblue')\n",
    "ax.set_xticklabels(['TPW','REP'])\n",
    "ef.simpleaxis(ax)\n",
    "ax.set_xlabel('')                \n",
    "ax.set_ylabel('${\\\\beta}$ Coeff')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF0AAABfCAYAAACKucvIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIb0lEQVR4nO2cX0xUVx6Av3vvzADDKDilCFosBWKiq4XqSsxKNj5sNlPfSn0RQjAxkU1KbXdN1KZNajUaNaJR3Afd2If6oi+1dktsrdkEq9FqYpBoRGylaqaCCjPDDMjM/bcPY/0TlGFwxuPo+ZJLJmcOv3vOdw8/7rn3nqvYtm0jea6oohvwKiKlC0BKF4CULgApXQBSugCkdAFI6QKQ0gXgSFfgnnIf1vAI2fNnp2sXaWfk8iXUN7N463/HUho3bdKt8DDYNo7pU9O1i/TT0Yk1NJTysDK9CEBKF4CULgApXQBSugCkdAFI6QKQ0gUgpQtASheAlC4AKV0AUroApHQBSOkCkNIFIKULQEoXQNpu172M7NixA9u2URSFy5cvs3LlSnbv3o3T6aSyspJr166xbNkyFi5cOGYcKX2cdHV10dPTQ2trKwAHDhxg+vTpVFdX43a7WbFiBTdu3ODjjz/m66+/HjOWTC/jJC8vjwsXLtDe3o5pmjQ0NDBt2rTH6gwMDJCfn58wlhzp46S4uJjNmzezb98+PvnkE95//31WrVoFwKlTp4hEIoTDYTZs2JAwVkLp58+f54svvqC0tJSysjK6u7upr69n7ty5bN68meLiYgYHB4lEImzZsuXZe/cCU1NTQ01NDX19fXz00Ud88803ACxatIgVK1aMO05C6fPmzWPWrFksXrwYn89HV1cXGzZswOfzMXv2bBoaGgDYunXrxHqSIfj9ftra2li5ciVTp05l4cKF6Lo+oVhJ5/RAIIDX68Xr9dLW1kZXVxcAa9eunVADMgWPx0NnZyc7d+6kpaUFv9/PggULOHfuHKdPn6ajo2Pcscad048ePcqZM2e4fv06O3bsID8/n1AoxOrVqzEMg6amJmprayfSn4wgLy+PPXv2jCr/6quvko417pH+7rvvsn79etxuN7/88guKolBfX09bWxstLS1s3bqV3377LekGvIokffbS1NTErl278Pv9lJaWUlVVxZw5cygvL8c0zXS08aUjofSOjg6uXLlCLBajsrKSt99+m2g0ypEjR5g8eTIVFRWEQiGqq6spLy9/Hm3OeJR0Ld79tfCvYNt4av+WjvDPhUjbjzDdpvznUymNK2ekApDSBSClC0BKF4CULgApXQBSugCkdAFI6QKQ0gUgpQtASheAlC4AKV0A6XsE4/4VY8Pfl7iuojzcVBVFVUDTQFNRNA3FoYFDA01DmWBzLNsmqpvohoVuWhimRaKL2naWG4djYjefxyLtz71EL/06/srKgx88ZtcmfhBVBUVzgFNDcTpRspwoTkf8oPxxcFQNtPgvx3SL8IhOZERnIrcNInkFePMyUHru3/+Ssli2aWLrBugGdkzHjt0Xomn3K9gYukFoWKffVIg5XRgeN46iHIqLpuDJc5OTm01Wtit+gBLw7//c5OfoPP6bsh7EyagnvBRNQ9E0yM4a9V3MMLl+J0JEt1E0lUk5GjMmucjRVBQrht17C363sW2bEduO//UkYF7/HS44SlPej4yS/iRMy+Lm3SGCQ1FQFKZNyaVgcjaaOtHs/5CYGkhBC0eTsdJtIBiJcv1uBLApmuKmMC8HVXl22ekmI6XrhkX370FipkVuloOyqZNxaJlz9ptR0i3b5lZgmDuheyiKQkVRHpNynKKblTQZId20LG6HRugLDYMNBZOzmebNzYhU8iTSLt0cCCWu9MjkKD4xUkFRGTYsboWjDEUNAPLcLt4o8ODMoFTyJNImXcnNwdYNshK9l9GOn38bMZ17kRGGw8Pcvh1Cj9zDEYvh0nWmuV3ke7JRFQP6AxiWDZaF4nKC04nicqA4nfEZbAaM/rRJ/+c//kXMsKh8wztmPduO5+pHcTk0/lxeQHXFVIqn5EBUx4oMY4eHsMJDmMEw5t0Axu+3MW7dxbwTwAqEsKM6qPFLCSg8nMnaNrZlodjEyxUl/iHB8cmyUj8bhTRKH9HjD5P6+8d4maQCqqLg0NT4pio4HSqmZfPT5V5+utw7xh5yYdJb8W1mvEQ1DFzRKM4/Nl1H03UcuoFqmqimiWboqJaNYlkJ+3BdL+FOYMqo8mPHjvHll19y8OBBALZt28bhw4epq6sjGAySlZXFmjVrnho3femF+ECbXTL2SE89kx58su5vEx2vp0/2MXw2Z1T5pUuXCIfD3Lx5k5KSEurq6jhx4gQffvghAHV1dXR3dzNz5swnxs3s/0gCGBwcpKCggNraWo4cOTLqe8uyCIfDeDyep8bIiFPGF4njx4/j8/lQVZXGxkaam5sB6O/vp7W1lVAoRHNz86jljo+S1Oq6iooKAM6dO0dTUxPff/89hYWFXL16FZ/Px5IlS1LUtReX8+fP4/f7AdB1nc7OTrxeL6+99tqD9JKIpFfXAbS3t7N9+3YOHDiAx+Oht7eXH3744Rm6khn09vZSVVXF0qVLASgqKuLbb79l+fLlScVJOqfv27ePd955B6/Xy/79+wkGgxQVFdHY2JhsqIxC13U2bdrEwMDAg7JgMMh3331Hc3MzAwMDHDp0aFyxxrUSY926dQQCAWbMmMHFixfZu3cvIyMjtLa2cvz4cebMmcO6deseW/7i29iGDSz+09Nz24vOjyf7iJz1cPZkTUrjjnukv/fee3z66ads27YNt9tNYWEhGzdupL29nfnz5/PZZ5+ltGEvM0mnl5KSEhwOB+vXrwfA5XKxZMkSDMNIddteWpJaXVdWVvbghP/111/n888/p6CggJ6enjFnYJLHSSi9qqqKw4cPjyr/4IMP0tKgVwE5IxWAlC4AKV0AUroA0nbBy53lYCRmcmfwXrp2kXZUQyU3N/GTYMmStncDSJ6OTC8CkNIFIKULQEoXgJQuACldAFK6AKR0AUjpAvg/aqvn/OaVeFIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 72x108 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sankey plot \n",
    "cell_type_str=[]\n",
    "for i in range(len(label2_exp)):\n",
    "    cell_type_str.append('FS' if label2_exp[i]==0 else 'RS')\n",
    "\n",
    "y_pred_str = []\n",
    "for i in range(len(label2_exp)):\n",
    "    d=y_exp_dictlist[i]\n",
    "    maxoccur = max(d.items(), key=lambda x: x[1])[0]\n",
    "    if maxoccur==0 and sum(d.values())!=0:\n",
    "        y_pred_str.append('AP')\n",
    "    elif maxoccur==1 and sum(d.values())!=0:    \n",
    "        y_pred_str.append('SP')\n",
    "    else:\n",
    "        y_pred_str.append('Unsure')\n",
    "df = pd.DataFrame(cell_type_str[::-1],columns=['cell_type'])\n",
    "df.insert(0, column=\"label\",value=y_pred_str[::-1])\n",
    "fig=sankey(df['cell_type'],df['label'],fontsize=9,colorDict=colorDict,fig_size=(1.0,1.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment-based classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The initilized parameters for experiment-based classifier\n",
    "accuracy_XY_em = 0\n",
    "conf_matrix_XY_em = np.array([[0, 0],[0, 0]])\n",
    "w_em=[]\n",
    "Coef_name_em=[]\n",
    "y_model_dictlist = [dict((j,0) for j in range(2)) for x in range(len(df_model))]\n",
    "\n",
    "random.seed(100)\n",
    "for ii in range(rep_num):\n",
    "    # Upsample minority class\n",
    "    SP_upsampled = resample(SP, \n",
    "                            replace=True,       # sample with replacement\n",
    "                            n_samples=len(AP)-len(SP),    # to match majority class\n",
    "                            random_state=27+ii) # reproducible results\n",
    "    # Combine majority class with upsampled minority class\n",
    "    model_upsampled = np.concatenate([range(len(df_model)),SP_upsampled])\n",
    "    #print(len(model_upsampled))\n",
    "    X_model = stats.zscore(df_model[['TPW','REP']].iloc[model_upsampled].values)\n",
    "    Y_model = df_model['cell_type'].iloc[model_upsampled].values\n",
    "\n",
    "\n",
    "    # Upsample minority class\n",
    "    FS_upsampled = resample(FS, \n",
    "                            replace=True,        # sample with replacement\n",
    "                            n_samples=len(RS)-len(FS),    # to match majority class\n",
    "                            random_state=27+ii)  # reproducible results\n",
    "    # Combine majority class with upsampled minority class\n",
    "    exp_upsampled = np.concatenate([range(len(df_exp)),FS_upsampled])\n",
    "    #print(len(exp_upsampled))\n",
    "    X_exp = stats.zscore(df_exp[['TPW','REP']].iloc[exp_upsampled].values)\n",
    "    Y_exp = df_exp['cell_type'].iloc[exp_upsampled].values\n",
    "    \n",
    "    # Experiment-based classifier: Training exp, predict model\n",
    "    accuracy_XY0,conf_matrix_XY0,accuracy_newXY0,conf_matrix_newXY0,w0,Coef_name0,dict_y_pred0=ef.classifier(X_exp,Y_exp,X_model,Y_model,'svm',random_num)\n",
    "    accuracy_XY_em = accuracy_XY_em + accuracy_XY0\n",
    "    conf_matrix_XY_em = conf_matrix_XY_em + conf_matrix_XY0 \n",
    "    w_em.append(w0)\n",
    "    Coef_name_em.append(Coef_name0)\n",
    "\n",
    "    for k in range(len(df_model)):\n",
    "        y_model_dictlist[k][0]+=dict_y_pred0[k][0]\n",
    "        y_model_dictlist[k][1]+=dict_y_pred0[k][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy of test data set:99.62160173160166%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIUAAABoCAYAAADfADNgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAO0ElEQVR4nO2de0xU177HP/PircLwUhCL4IOrAtVa3x7O6fV6W5vWpqZqbdqm1XRarbFW66tX0V5TR+IjEZuqrY22NcYmVo0aqeec1moLSH2A4BWKgDwERAWGxwjDwNw/uOwjjDLOADOby/okO9mzZq+9voFvfmvtvX/7NwqLxWJBIHgIpasFCOSHMIXACmEKgRXCFAIrhCkEVghTCKwQpuilXL16lZdffpmkpCQALBYLO3bsYM+ePaxfv57KykoAamtr2bBhA3v37kWv19Pc3Gzz3MIUvZTS0lJGjRolfU5JScFgMPD+++8zc+ZMEhMTAfj666+ZNGkSOp0OjUbD2bNnbZ5bmKKX8uKLL7b7nJqaypgxYwCIiYkhOTkZgIsXL7ZrT0lJsXludTdrbccv/7W3J09vF3/brHO1BJvEPBUn7X+asJgjR45In+fNm8e8efMe27eqqoqoqCgAvL29qa6ultq9vLyk9qqqKps6etQUAvtQKlXSvi0TdMTPz4/6+noA6uvr8fX1ldqNRqPU7ufnZ1uHHZoFPYxGqZY2e5k0aRJZWVkAXLt2jSlTpgAwceLEdu2TJ0+2eS4RKWSE2g4zHD9+nJycHEwmE0FBQUyePJmUlBR2795NWVkZK1asAGDRokUkJCRQUFCAyWRi5syZNs+t6MmnpGJNYR9/iZot7Z/PPuEyHSJSyAiVQh6zuTCFjHBTaVwtARCmkBWqh64+XIkwhYxwU8nj3yEPFQLAvquPnkQeKgQAaMT0IeiIWphC0BGNSphC0IE+HynG6V6hxdya8HE/p4ji3zKI+I8JNJua8PDtR/7f02gyNrTr497fm/DnnqGhuhZQUHjuCgC+Q0MIHDUUS0sLVfml3M8pdFjXyZMnuXHjBvX19cyZM4eYmBjpO5PJREJCAoGBgZSVlbF69Wo8PT0pLy9n9+7dhISEYLFYWLJkiUNju8kkUrjsFlplbjHp+0+Svv8kxb9l4BcZitrDjcJfr3L3ej7hz4236hMxcyJll7IpPHcV7yA/+oUGggKGvziFm0kp3ExKIWLGsyjVjv1x6+rqOHDgAKtWreKTTz5h48aN7b4/ceIEQUFB6HQ6oqKiOHToEADbtm3jtddeY/HixeTl5XHt2jWHxlcpldLmSlw2uk+wP0OmxTL038fj4dcPv4hQakvvAlBTUoE2MtSqj19ESLtj/CJD8dQOwNxgwtLcAhZ4UF2Lz6AAhzSlp6cTERHRqs/HB7PZTEVFhfR9amoq0dHRQPuElYsXLzJ69GirdntxU6mlzZV0aoqkpCS+++47APLz81mwYAHz5s3j0qVLXR648NcrFP2Wwe2L/0P0gv9E4+VBc6MZgGZTExpvD6s+Kjc1lpbW53fNjU1ovDxa+5mapGPa2h3h4YQUsE5KeVzCitFoRK1WP7KPPaiVSmlzJZ2OfubMGelR69atW4mLi2Pbtm0cPny4ywPXlt4DwFRnRKFWYmlpQeXe+odVuWloqm+w6tNsMqNQKlqPcdfQZGygydiAyu1fzwza2h3h4YQUsE5KeVzCipeXF2az+ZF97EGtUkqbK+l09OHDhxMcHIzBYCAzM5OFCxcSFhbG0KFDuzSoV4AvA8eOAEChUqLSqLl34xb9QgIB6D84iMq82wCoPd1RalrNUpVf2u6YqrzbPKg0oPZwQ6FSggI8fftRV3bPIV1PP/00eXl5QOv6QqVS4efnx717reebNGkSmZmZQPuElYkTJ3L9+nWrdnvRKFXS5ko6nbwqKiowGo0cPHiQGTNmSCHy/v37XRrU3GgiICoc9/7eePj2I/fU71TeLMF3aAjhfxuH+wAf8s+mARA2NYYHlTWUX8kh/+xFwp97Bu2IMIx3q6i93bq+yD2dzPBZU7C0tJD/jz+kqxp78fHx4Z133kGv12M0Gtm4cSOpqamcPn0avV7P7Nmz2bp1K1988QUVFRWsWbMGgJUrV5KYmMj58+eJiIhod8ViD8r/i4KuptMkm+TkZLZt24ZWq0Wv1xMQEMCmTZsICgrigw8+sHlykWRjHwmvbJL2Vx2Pd5mOTiPFvn37WLZsGePGjUOhUNDS0kJ8vOvE/n/H1WuJNjpVMX78eOLi4jh8+DA6nY6ysjJn6eqT9IqFZhvvvfceo0ePJjS09d7BzZs3e1RUX0Uul6SdTh95eXnSa2ZFRUXS/k8//cT27dt7Xl0fQ6OWx/TRqSkKCgr45ZdfgNZr9Lb9goKCnlfWB1HJ5OqjU1PEx8czduxYq/aMjIweE9SX0Tj4zKa76dQUjzIEQGxsbI+I6evI5T6FyKeQEb1iTSFwLiqViBSCDrj6/kQbwhQyQkQKgRUqOyPF3LlzcXd3ByAuLo6FCxeyc+dOvLy8uH37NsuXL0er1dqtQ5hCRqjtjBTTp09n6dKl0ufk5GQMBgMff/wxFy5cIDEx0aFnVT1qCjk9mTy3fp+rJbTjr//9nlWb6qGrjyNHjtgsb5STk8NXX30lJRl3rHu1efNmh7SJSCEjHl5TPEl5I51OR3R0NBUVFSxatIjY2NhH1r2yF2EKGaG0c03RlkQcFBSEyWRCqVQ+su6V3Toc6iXoEdQapbTZIi8vjx9//BFofR+lsbGRGTNmPLLuld06HOol6BHsuc3t4+PDzz//zJ07dygtLeXTTz9l2rRppKWlWdW9shdhChmheoII0UZwcDC7d++2anfUCA8jTCEj7L1P0VMIU8gIpXggJuiIuM0tsEJECoEVwhQCK1S9IR1P4FwUYk0h6IhSIyIF0Hk5oaNHj5KVlUVAQAB5eXls2bIFd3d3Ll26xA8//MCQIUMoKSlh8+bN0svPjtJvcCBjF80mJeE7mpuaGTZrMg/u1+AVMICCf17CVGtsd7xSo2bwlGgs5mb8IkO5eSaFFnMz4c+Np778Pt4D/bmdkimVXHgS5LKmcKkKW+WEzGYz69atY8mSJbi7u3PhwgUAdu3axbvvvsuHH35ITU0NaWlpXdKhVKsIjh2Oqab1YVLohFE8uF9D8W8ZVGTmMXTGs1Z9nvrrWO6k51L8+zWyj/2KqdaIxsuD8svZFP9+jdspmUTMnGifDpVS2lyJSyPF48oJBQUFAUiPji0WC+Xl5QwePBgArVYrVYsxmUyEhYV1SceQvzxN8W8ZBIx8CgAPv34YisoBaDDU4Rs+yKqPb3gIDZW1qD3cMDeYKLucLZVGAEChoLnJbJcOR2t1dTd2WbKoqIj8/PxuG9xWOSGAW7dusXLlSqKiohg5ciQAS5cuZd++fVJZhIEDBzqsYUD4IBqqamk01EttNSUVUnEUn4H+qNytq+v7DPKnvqKS4t+v0T8sGP+RQ9p9P+iZKAp/uWyXFoVKIW2upFNTrF27loSEBACOHTvG3Llz2bBhA7t27eqWwW2VEwIIDw9n+/bteHh4SPW3PvroI/R6PfHx8ajVao4fP+6wBm1kKG4+ngyZHovKw43QSWOoKamg0VDP4CnRqNw0NNbUW/VrNjVJ64WakgoGPPUvY4Y8+29U5ZXYtZ6A1nVK2+ZKOjWFj48Pq1atAmDv3r188803fP/999TU1HTL4LbKCR07dkw6NiQkhNLSUgDu3bsnJZD4+/vT1NSEoxT88xJFFzIoupBBc4OJ26lZmB80UpF5k5LkTBqqarib1RodHy61ZCgsx2OADwDuA7x5UNn6NwmOHU6zyczd6wVoh9s3rSnUSmlzJZ1acsCAAQBkZ2ejVqulH0d1NKOnI7bKCZWWlrJ9+3Y8PT3Jzs5m7dq1AKxZs4bPPvuM0NBQiouLWbRoUZe1hE2NQeXhxuApMVTeLGHI9Fiqb5Wh0qgpOp8uHdNWaikvKZUh02JprKnHzduTwnNX6BcayLBZk6m/U8mgZ0aiUKmozC1+Yg1KmRRX7bS80bp16xg9ejSnTp3ipZdeYsGCBWRkZHDgwAF27tzpTJ1dpjck7t69+Lu0HzhxqjPltKPTOLVu3TpaWlp48803WbBgAdBahqC7pg9Be+Sypuh09MzMTJKSkvD392fq1Knk5uZy6tQpTCaTs/T1KRS9Icnm0KFDrF+/nqqqKpYvX463tzfLli2TsogF3YtS0wt+WG7EiBHSewRnz56V3ja6c+cOwcHBPa+uj6Fwca2rNmyWN2qrc1VdXS1qXvUwit4QKR6ueeXh4SFqXvUwcrnNLWpeyQhFF5/0dhei5pWM6BWmEDgXpfixWkFHRKQQWKFU9YKrD4FzsTdSdJbK2BWEKWSEwo5I0ZbKePToUerq6njrrbek0gRdRR630AQAKDVu0mYLW7+M2BX6TKR41KNquaFQPHnNq8elMrblt3aFPmOK3oatmldPksroKGL66KU8KpWxO6IE2Mi8EsibkydPcv36dYxGI3PmzOm2O83CFAIrxPQhsEKYQmCFrK4+rly5wqZNmwgPD2fYsGEA/PHHH+h0OpKSkggKCiI3N5fnn3+eWbNmOU1LREQEf/75J2+88QbR0dF8/vnnDBo0iJqaGurq6tDr9T2qxelYZMbq1astZ86ckT6fO3fO8sorr1hqa2stFovFUlZWZjlw4IDTtdy4ccPy+uuvWw4ePGj59ttvpWP0er1TtDgTWU8f+/btY+zYsWi1Wvbv3091dTUDBw7k7bffdrqWqqoqtFotWq2W06dPk52dDcDq1audrqWnkdX00caxY8e4fPkyWVlZzJ8/ny1btpCYmMgLL7zAmDFjWLNmDZGRkU7RcubMGVJTUyksLGTHjh34+vpiMBhYsWIFZrMZnU7Hq6++6hQtTsPVoaojD4fsoqIiS1NTk/RdY2Oj5csvv7TMnz/f6VoWL15sSUtLa/d9ZmamZcKECZaCggKn6HEWsp4+wsLCUKvVUjETNzc3Zs2ahdlsX92H7kCn07Fnzx6OHz9Oeno6AGPGjCEyMpLm5man6+lJZDV9pKenk5OTg8lkIiIighEjRgAQGBhIfHw8AQEBFBQUSG/CO0tLbGwsMTExNDY2cuLECfr378+wYcMwGAxMmDDBaVOZsxB3NAVWyHr6ELgGYQqBFcIUAiuEKQRWCFMIrBCmEFghTCGwQphCYMX/AsJ3gOu8kLzAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 129.6x93.6 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Confusion Matrix\n",
    "df_cm_XY = pd.DataFrame(conf_matrix_XY_em/rep_num, \n",
    "  index = [ 'FS', 'RS'],\n",
    "  columns = ['FS', 'RS'])\n",
    "print('The accuracy of test data set:'+str(accuracy_XY_em/rep_num*100)+'%')\n",
    "fig, ax = plt.subplots(figsize=(1.8,1.3))\n",
    "ax.set_aspect(1)\n",
    "res = sns.heatmap(df_cm_XY, annot=True, vmin=0.0, vmax=100, fmt='.2f',cmap=cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '${\\\\beta}$ Coeff')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGMAAABlCAYAAABZYl2IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAI80lEQVR4nO2dfUhTXRzHv75Upubro6ZmbL6RFIZhqRE133o0FMogi8BMLYWaaP8YJNj+moFJZGEZNktIBhZlMAMrlMyGRFohs9R8mWVW821zzmWe549oTz76tLvprnfufv5yh8P9fuXLOfeeu985syKEELAwAuvlNsDyL2wYDIINg0GwYTAINgwGwYbBIBgZRldX13JbWBYYGcbMzMxyW1gWGBmGMUilUpw5cwZSqXS5rRiN7XIbWCqqqqrQ1dUFtVqNyMjI5bZjFIwfGTOzs2Z13cXA+JFha22Nq629evv9tTMRoz+s8dfOvyn1P7WDuxT2lhTGh0EV702h8N4Uutw2FgXjpylLgg2DQdAWxuvXr7F582aMjIzQJWl20BKGRqNBXV0dvLy86JAzW2gJo6KiApmZmXRImTUmD6O1tRW+vr7w8fExtZTZY/Iwnj9/DoVCgYqKCiiVSlRXV+PDhw/z+onFYqSkpCAlJcVi7ytWdBYkxMTEoLa2Fm5ubn/sJ5PJEBISovtMZRFnKExc9NH2NFVZWQmlUomqqioolUq6ZM0K2lbgmZmZ7E1cD+yij0GwYTAINgwGoTeMc+fO4cePH1AoFHT4sWj0huHj4wMbGxvU1NTMab99+7bJTFkqep+mbGxscOzYMcjlcrS2tgIACCEYGhpCWlqayQ1aEnrDkMvlqKysxKNHj5CUlKRrf/LkiUmNWSJ6pylfX1/Y2tqir69vTvvHjx9N5cli0TsyrK2t2WmKJvSGkZOTg5ycHDQ0NCA+Pl7Xzk5TSw/ldcbvQQBAbGzskpuxdCi/m+rs7ERxcTFcXV2RmJgIe3t77Nq1y5TeLA7KI0MsFkMoFCIkJAR79+7Fs2fPTOnLIqEchpeXF7y9vWFlZQUAWLdunclMWSqUp6nx8XEUFhZieHgYw8PDsLe3N6Uvi4RyGAUFBWhsbERPTw8CAwOxZ88eU/qySCiHIZPJcOXKFQwODmLjxo3w9PSc89Uoy+Ix6AZ+/fp1SKVSXL16dd6Lw5XCcla9Ux4Z3t7ecHd3BwB4eHis2II0qlXvhkKlAIJyGGNjYygqKoKvry8GBwfh4OCwKHMs8zH4Bt7d3Y2YmBjweDwT2rJM9IYxNjaGz58/IzAwEDweDzweDwMDAwaJvHnzBjU1NQgICIBMJkNeXh78/PyMNr1S0XsDv3DhAmpra/H9+3ddW3l5Oerq6iiLfP36FSdOnEBWVhbi4+NRXl5unNsVjt4wHBwcUFhYiLVr1+rahEIh3r59S1kkNjYW/v7+AIDZ2dk512L5F71h/F8ppouLi8FihBBIJBK2mO1/0BvGxMQEpqen57RptVqjSjTLysqQnp6+YEU6W/hM4QZ++PBhZGRkYMuWLXB3d8e3b9/Q2dkJgUBgkNCNGzcQHh6O8PBwNDU1zXudkpqaitTUVAA/V/uWiN4wOBwOKisr8fTpU3z69AlhYWHIz883aN5vaGiASCRCQEAAysvL4eTkxL7bWgBK6ww7Ozvs27fPaJH4+Ph53xSyzIct72QQesPQarWQy+V0eLF49E5TWVlZ8PDwwNDQEBISEiASibB+/XoUFRVh06ZNdHi0GPSGQQjBxYsXMTAwgOTkZDQ1NUGr1eLs2bO4efMmHR4tBkorcIlEAi8vL+Tm5sLFxQWenp6Ympqiw59FoXdklJSUQCQS4c6dO9Bqtejq6oKzszOUSiUmJibg5OREh0+LQG8Yjo6O4PP54PP5IISgr68PMpkMq1atQn5+PhQKBe7fv0+D1ZWPQRssrayswOVyweVyF7XuYFkYdp3BINgwGAQbBoNgw2AQbBgMgg2DQbBhMAg2DAbBhsEg2DAYBBsGg2DDYBBsGAyCtmPxHj58CJlMhsnJSRw8eBChoeZ9iLwpoCUMlUqFqqoq3L17FyqVCmlpabh37x4d0mYFLdNUe3u7rvDZ0dERMzMz+PLlCx3SZgUtYYyOjs7Zquzg4IDR0VE6pM0KWqYpV1dXqNVq3efJyUm4urrO6SMWiyEWiwH83CX1e71tDMX9/yMjI3oPMP7Fn+p5Ta1na2uLoKCg+Z0IDSiVSnLgwAHd3/v37zeJzi8NulhqPVpGhqOjI44fP47i4mKo1WqcP3+eDlmzg7ZH2+TkZCQnJ9MlZ5asqEXfr/0d5qpH668EsPyZFTUyzB3G/05fdnY2/P390dHRAa1Wi7CwMHR0dGB8fBwcDgdcLhc9PT04ffo0mpqacOvWLQiFQuzevRulpaVwcnJCVlYWmpubUVFRgcLCQgQHB+vVffXqFQQCATgcDvz9/fH+/XscPXoUCoUCAoEAhw4d0u3eUqlUyMvLg0AgQFtbG5KSkjA0NITg4GCkp6dT/2eX9NnMBDx+/JgQQohIJCKXL1/WtRUUFJD6+npCCCEtLS2Ez+cTQgjh8XhkamqKEEJIbm4uycjIIIT8fKSurq42SPt3DZlMRo4cOUIIISQ6OpooFApdv8bGRkIIIVKplJw8eZIQQsjs7CyJjo4mSqWSsh7jp6mFDqb8b9vIyAicnZ0BANu2bcPLly8xMTGBsLAw9Pb2QqvVQiqVYseOHUb7GB0dXXCBV1JSsuD+RI1GAwBYvXo1ZQ3GT1N/or6+Hp2dndBoNMjPzwcAREZG4sWLF5iamkJERATa29vR1taG7u5uxMXFGaUhlUrR39+P0tJSXfulS5ewZs0aqFSqOf37+vpQVlaGsbExCIVCywkjMTERCQkJc9qioqJQU1MDOzs7xMXFYefOnWhpaTH6VIZfGqdOnUJ3dze2b98OAMjLy4Obmxt6e+ceh8ThcMDn843SYvw0ZSgbNmyASqWCWq2GlZUVoqKiIJFI4OnpuajrZmdn49q1a/PauVwumpubF3XtX5jFyJDL5WhpaYFGo0Fvby/Gx8fx7t07aLVaBAUFISAgYE7/iIgIcLk/D9vy8/MDIQRRUVEGaba3t+s0tm7ditDQUExPTyMxMVH3Q152dnYAAIVCgfDwcDx48AD9/f2QSCRGbZlgF30MYsVNU+YMGwaDYMNgEGwYDIINg0GwYTAINgwGwYbBIP4BH1rKssei7ccAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 72x93.6 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Beta Coefficient \n",
    "df=pd.DataFrame(np.mean(np.array(w_em),axis=0),columns=['Beta'])\n",
    "df.insert(0, column=\"Coef_name\",value = Coef_name0)\n",
    "fig, ax = plt.subplots(figsize=(1.0,1.3))\n",
    "ax = sns.barplot(x=\"Coef_name\", y=\"Beta\", data=df,color='skyblue')\n",
    "ax.set_xticklabels(['TPW','REP'])\n",
    "ef.simpleaxis(ax)\n",
    "ax.set_xlabel('')                \n",
    "ax.set_ylabel('${\\\\beta}$ Coeff')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF0AAABfCAYAAACKucvIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJNUlEQVR4nO2cXUxcVR7Af+feO3dgmOFjoFBK61Zgu9GtBWlkadZkq8m60/ZlbX3QkqZmG+Wh+JGY1DW20UJalZiayIshcV/6ZLIrttFYtY3BaGIWt6FEt5TaYluxQxcowwzzce/ce/dhKlGrDFCmF+r5JZdkhnP+nPObM/977rn3IBzHcZDcVBS3G/BrREp3ASndBaR0F5DSXUBKdwEp3QWkdBeQ0l1Ay0XQofsfwE4kyFv7+1yEv2kk//1flJiP288dW9C4OZFuT02B46BVVOQi/M3D6ceOxhc8rEwvLiClu4CU7gJSugtI6S4gpbuAlO4CUroL5OTi6Fbk5MmT7N+/n9WrV1NdXc3g4CDNzc3cddddHDx4kMrKSiYnJ4nFYrz88sszxpLSZ0lDQwN33HEHGzduJBQKMTAwQFtbG6FQiDvvvJMdO3YA8Morr2SNJdPLPLl69SrBYJBgMMh7773HwMAAAM8++2zWurMa6YcOHcJxHIQQnD59mscff5zXX38dj8dDXV0d58+f55FHHqGpqenGerIEeP/99/n888+5cOEChw4dori4mEgkwjPPPEM6naalpYWtW7fOGCOr9IGBAYaGhujs7ATg8OHDVFVV0djYiM/nY9euXVy8eJGnn36at99+e2F6tojZtGkToVCI3bt38/XXX3PPPffQ3NxMc3MzX375Jbt27aKhoYHVq1f/Yoys6aWoqIhTp07R09ODZVns2LGDFStW/KjM+Pg4xcXFN9qfJUVLSwtvvPEG77zzDn19fQCsXbuWmpoaLMuasW7WkV5ZWcnBgwfp6uriueeeY9u2bTz55JMAfPbZZ8RiMaLRKG1tbTfek0VMX18fZ86cwTAM6urqWLduHalUiiNHjlBYWEhtbS2RSITGxkZqampmjCXm8ljdyMgITz31FNu2bSMcDk+nl59y7g9/BMfB/8Cf5967RUTs6EcQ1qi58smCxs2aXoaHh+nq6gKgoqKCpqYmTNNc0Eb82siaXvx+P/39/bz22mvYtk04HGbLli20t7ej6zrr16+nvr7+JjT11mFO6WW2yPQyM/LiyAWkdBeQ0l0gpwte8d4vshcSgBAIRQFFQagqaBpC0xCqitC0zGtlfuPDccBM2xgpm5ThYBg2luVg2w7ZTmZJtZJVjM7r785ETqXr1bej+Hy/XMC59sNxcGwbLBvHsiCdBkWAqmas2TZKIIBnZRV6TTWe227Ds3IlarAEIcR1YQ3DZmAwRu8XE5zqn8QwbPCD7YDuUdA0gaIKrq/5Y0Y++g97xT9uRMHPklPpQvei5OffcBzHccA0Mc5/Q+r0AChK5oPw+dBra9B/twa9qoqIp4RPTln0fDqOYdioqiAQ0AgENJJJm2TKxjBsDAOEAmTRrigKin+GQTNPcjplzN+wATXgX+jw0zimiR2Pk4olGJuwiU6a2AhS/iAxvYRRiojrxSSUPAJlAZZVBVhWWYC/OB/dq2Y+vBnQfTrrmlYueLuX9E0MS9H4LqoTDmfGjV6gkDbSqPEpypIR1gQEBYogz+dBSSlwnszxfVrLwqlLE7Tdu51//n3zgrZ7SUp3gIkJk6Fv4qTTDh5dwTQyqaOk2EtFRSEBv8bPpPs5oQz9j2Ry4Zc8lpz0eMLi4qUEkUgaj0cgBJiGTWmpzopKL/l5qttNzMqSkZ5IWnx3OcXYmIEQ4PEITNOhsFDjN7fl48tf/LK/Z1FLtx2Hyck04ZEUk5NpgOlU4jhQW+MjGNSzTv0WGzmRrhQUYBsGTnyKdHwqewUhMvNyoZBOOySSNrGoyVTMIO2o2KpOnq6TNAWmYVNRoVNVmY+mLTXdGXIi/W+e/aSFzZ6/1s5YzrEdkkmTaMRgfCTO2OUojpFA1xLkFyUJlsbxx0fJm7xCYWqMFQU2pWVePKqBiJk4BT6Ex5OLLuSUnEifmrJwHPhX92V078xzYUURqKpA9yh4i0vQSoPEUzZXomkuCPAFVTb+qZSGDSWU6Abpy2HM4WFSA2cwBs9iRSKZb4ptIXQvwutF6Pq8lw1uBjnN6aqWkflLOGRGu2namCYkkza27RAM6tx/Xxn16wqprSlAVb9PI17UQADvmt/iv28jjuNgXb1Kevg7zOFhjKFvML/9FmtsHLiWsnDAssG2MxdDiopQBLOZT4oc/YOQnEq/ciWVtYwiBKom8GgCzaOg6wqjowbHT4xy/MRcFpsqrx0gsMlLT5FvTZGXjuO1EuhWEq+VwGOn0GwTBTur1CG7HEtZ+FlRTqXXrSvKZfgs5ANlWED82jFXTvRepkjMYiIwRxZv4ruFkdJdQEp3ASndBaR0F5DSXUBKdwEp3QUW9dLuYuKHG71qazMLeb29vbS0tHDs2DHKy8s5e/YsoVCIzZtnvr0npc+Sn270Aujp6eHVV1/l8OHD+P1+wuEwH3zwQdZYMr3Mk66uLu6++26CwSBvvvkmExMTLF++nJ07d2atK6XPke7ubg4cOMDHH38MwEsvvcTo6CibNm3iscce49y5c1ljzFr6hx9+yMMPPzz9uqOjgw0bNtDZ2Ul7ezsdHR3z6MLS48EHH+T555+no6MDn89HeXk57e3t9PT0sH79evbu3Zs1xqylf/XVV0SjUS5dugTA9u3bKS0t5YknnmDfvn309fUxODg4/94sMVatWoWmabz44osA6LrO5s2bSafTWevO6kQ6OTlJWVkZW7du5ciRI7S2tv7o97ZtE41G8ftz9zSX2/xwo1d1dTVr1qwBYNmyZbzwwguUlZUxNDTEnj17ssaalfTjx48TCoVQFIWdO3dOSx8bG6Ozs5NIJEJra+t1Wx1vJerr6+nu7r7u/d27d8851qyknzx5kuHhYQBM06S/v59gMDidXiRzI6v0cDhMfX09Dz30EADLly/n6NGjPProo7lu2y3LjCdS0zQ5cOAA4+Pj0+9NTEzw7rvv0trayvj4OG+99VbOG3mrkZNHpRvv/RTHgb/8edlCh76pnOi9TFHTFMf2bVnQuPLiyAWkdBeQ0l1ASncBKd0FpHQXkNJdQEp3ASndBaR0F5DSXUBKdwEp3QWkdBeQ0l0gJ094FRSoJBIWIyPZN3otZjz54PMuvKKc3MSQzIxMLy4gpbuAlO4CUroLSOkuIKW7gJTuAlK6C0jpLvB/Yh1G7xsYEL8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 72x108 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cell_type_str=[]\n",
    "for i in range(len(df_model)):\n",
    "    cell_type_str.append('AP' if Y_model[i]==0 else 'SP')\n",
    "\n",
    "y_pred_str = []\n",
    "for i in range(len(df_model)):\n",
    "    d=y_model_dictlist[i]\n",
    "    maxoccur = max(d.items(), key=lambda x: x[1])[0]\n",
    "    if maxoccur==0:\n",
    "        y_pred_str.append('FS')\n",
    "    else:    \n",
    "        y_pred_str.append('RS')\n",
    "        \n",
    "# Sankey plot mapping AP/SP to FS/RS        \n",
    "df = pd.DataFrame(cell_type_str[::-1],columns=['cell_type'])  # [::-1] : reverse the array for plotting\n",
    "df.insert(0, column=\"label\",value=y_pred_str[::-1])\n",
    "fig=sankey(df['cell_type'],df['label'],fontsize=9,colorDict=colorDict,fig_size=(1.0,1.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI0AAADMCAYAAABOWiBsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlVklEQVR4nO2deZBc1XX/P+e+16+7p3t2aaTRAlqRJSEQKOziZ/jF4SeBEjCEykJRScoBYSNZKWLjkskvbBWBAwg5JFUgij8SkvKvyjb7ZmyHyLLNLiRAyyCBJKQZrbP3TG/v3fv74/W0ZtBIMz3To1GP3qeqa2Z63nt9u/v7zj333HPuFWOMISCgANRoNyCg9AhEE1AwgWgCCiYQTUDBBKIJKJhANAEFE4gmoGAC0QQUTCCagIIJRBNQMIFoAgomEE1AwQSiCSiYQDQBBWOPdgMKpf2hh8DzCJ17buEni4BtI6EQEokgsRgqFkPKyvzf43EkHC5+o8cYJScaMhkATDI5tPONwWiN0Rq0BmNA5Qyu1kg4jKqtxZowATV+PFZtLRKJFKnxY4PSE00OcZzhnd/Pc8YY8Dz0kSN4TU2+ZQKkvBz77LOxJ09GjR+P2CX7sRWFM/vdfwXJdV/YNhKNAjkhZTJkP/2U7KefglLYU6diz5iBNWkSEgqNcqtPPYFoBkBEwHHyls1ojbtvH+7evb6AZswgNHs2qq7OP/YMIBBNgYhSSCwGgPE83M8/x921C4nFCM2fjz1tGqqsbJRbObIMSjRr167FGIOIsH37dh5++GFqa2uH/eIfffQR9957L9/5zndYsmTJsK93qhHLQuJxvwvLZsm8+y6Z997DmjABe9YsrClTULlubiwxoGh27NjB7t27eeKJJwB49tlnSafTRXnxpqYm5s2bV5RrjSa9uzBjDLq5mfShQwComhrs6dOx6utR1dWIZY1ya4fPgKKprKxky5YtbNiwgcWLF3PrrbcCsH79elKpFMlkkrq6Os4991x++MMfcvPNN7Nz50601jz22GO8+eab3H///axYsYKNGzcyffp0vv/97wNw3XXXsXHjxj6v9/bbb/PSSy8xefJkRIQ777xzBN72yCEiEI0i+E606ewk88EH/khMBKuuDmvSJFRtLaqqConFSs4XGlA09fX1rFmzhvXr17N69Wpuuukmli1bxrvvvsszzzwDwCuvvMJFF13EokWLWLBgAbfffjt//Md/TEtLC9dccw3PPvssV1xxBX/2Z3/G17/+9bxo+iMej3P33XdTXV3NihUrSCQSxOPx447TuXhNQSiFKAUip+SLEhGIRPJxHqM1uqUF79AhX0TGgAiqosIXUGWlH2wMh33LZdtgWUNrr1JIRcWIvM9B+TSLFy9m8eLFHDp0iFWrVlFRUcHkyZPz/1+2bFn+9ylTpgBQVVVFd3c3NTU1AEydOhURITTAEDUajbJ+/Xqqq6tpamqira2tX9F4+/YNpul9McaPBjuOP6TuHW/RGizL/6JCIV9cRUaU6iMi8IVkUinc/fth795eBw/vy069vZmO1z9G7BDTP39jWNf6KgOKprGxkVdffZXbb7+dCRMmcOmllxKJRGhsbMwf88ILL3DDDTcAnFDZg1X8Y489xq233srll19OQ0PDCY8LzZw5qOv1pifmolMpjAihadOwZ88G10UnEuiWFnRzM7qtjXzhqYj/JVvWiNy1olSfIX2xMJ5/E+jO7qJeFwYhmng8zscff8zjjz+O1pqDBw9y++23k0wmefzxx3Fdl6uvvprPPvuMhoYGXnvtNS6++OK82BYtWkRjYyOvvfYa5eXldHZ28sorr7Bs2TJeeOEFGhoayGQy1NXVceGFF3LdddfxzDPPsGnTJj7//HNeffVVli9fXpQ3KyIQDmOFw368ZfduvMZGIlddhdPLITdaYzo68Jqb8Zqa8BobMYmEP4IMhXxrUWJ+SDGRUqvlbr//foChTVj2g0mnMek0zkUXEZo/v18x9Di07oEDeLt34x08mO/qTlcBdb3xGzp/vR2TSDHz8G+Keu0zPrgn4TDYNpn33sOkUjiLFh0nAhFBKipwKipgzhxMOo3X2Ei2oQHv0CHfAkUiRe9iTlfOeNGAH6SjvJzsxx+DUjgXXHBS6yHhMPaMGdgzZqC7unC/+ILs1q3ozk5/1BKJjIl4zIkIRJNDlPKFs3kzEovhzJkzqPNULIazYAGh+fPRhw6R3bULd/duP/VCBIlGx5yAAtH0QpSCWIzM229jVVVhTZhQ0LlWfT1WfT3mssvwDhzA3bMHd+9ejOf5PpBtQzg8IsP5U0kgmq8gtg2OQ+q//5vo9dcPafJRbNtPn5g6FXPFFfmAnrdvH97hw74V6jk2FPLjRSM0pB8JSlY0OpEo+ByxLD9wN0ASlYTD6ESC9MaNRK65ZlhfpiiFNW4c1rhxMH8+2vMw7a14rUfQR4+gjzaj29sxXem+mWFG+7EiA0LhA1xVEwFrZERYeqJxHD9HeJA+Rw/GGEwi4d/pyaTvb5xk3kdisfwIyfna1wb5Gi5Gd2N0CqOTaK8brRO555L+8yYNGKhSUAnMAigDLwpJFzIepD1wNeIaP1KtwRQqnGgIXhNoL+y0wVByoqlcvXpY5xtj0EeOkN26FXfPHt/y9JO+ICK+f/Pee9iTJyPxOJgsxqRyokihvS6014H2OjA6gdFpkB5/xYDRgIWIBSgQC5ETCFUBRUwCTB/a7c9tjQAlJxoA7SXAeEM40/h3bI1NaPF5WHPqSf/uQ3T7EYhHc72Dxhidu76HzqZIvPUscnmt///8F25yE44WggViIypeMn7JcCg50bQfeALIYjlThn+xMJj/ZWE+8jBfHoIyC1HCMedCICKYg92oQzWoSeXDf80xQMmJBvyUCKVixbmcA+biODp2BL29GRML5YTjIwImDPqjQ0hdDLFLe7hcDIJPAN9/UfPHo+aPh64sRvf1BcSxIOmid7eNTgNPMwLR5BAR1Nxa1MwqpCvLV+dxTcTGbD2KSQ/FlxpbBKLphYigFk5AJsSgO9v3f7bCuBq9q2WUWnf6EIjmK4gS1CX1SNjGpN2+/4zamIYWTMrt/+QzhEA0/SBhG+vSSZDVffwbsZQf5znDrU0gmhMg48pQc2uRr3RTRGzMztbjrdAZRCCak6Dm1ELc6dMdiaUwnkF/0TZ6DRtlAtGcBLEV1kX14PbtpojYmM9aMK4+8cljmEA0AyC1UdS0Skge66bEVpDV6L0jMBtYAhRVNJs2beL6669n1apV/PjHP+bOO+/k97///aDObW1tZeXKlTz33HPFbFJRUOeOR5T0sSzGsfyRlC6pvPyiUNRphAsvvJC5c+dy1VVXsWTJEnbs2MEDDzzA5ZdfPuC51dXVXH311cVsTtGQiI3MGwefHIG4nzwujoVJZDCHu5CJxxfznW5s2rSJ+++/n2nTpjFjxgw+++wzbrnlFhYsWMCaNWuor6+no6ODRCLBww8/fNJrjejcU2trKzU1NbS0tPD4449z1llnsW/fPu666y527tzJD3/4Q6677jr27NnD+PHjmTt3Lh988AFNTU18+umnfO9732PWrFkj2cRBo2ZW4+1sxWQ9JJTL+bUUpqEFSkA0J7qhlyxZwrx58/I1+j/60Y8GvNaI+DSvv/469913H+vXr+fBBx/k6aef5pJLLuG2227j0ksv5amnnsrXfi9atIh169Zx0003ATBx4kRWrFjBHXfcwdq1a0eieUNCbIUsGI/0nkaIWOij3ZjOIdSVjzI9N3RNTQ2vvvoqO3bsAOAHP/jBgOeOiKVZunQpS5Ys4c4772TXrl3s2rWLa665BvBrvZ9//vn8sdOmTQPga1/7Gtu2bcvXiE+dOpXdu3ePRPOGjJpSjrf1KCbjIY6f02sM6N1tWOfVjXbzBsXrr7/OO++8w969e1m7di1VVVW0t7fz93//97iuy/Lly7nxxhtPeo0RHT0tX76cJ598klmzZrF//34A9u3bx+zZs/PHfDVpqampCYAvv/wyL6jTBbEUMn8ckullbaI25ou2khl+L126lPvuu4+ysjJ27dqFiHDLLbfw6quv8thjj/GjH/2IPXv2nPQaRbU0mzdvztdmn3/++Zx33nmk02m6u7t5++23aWxsZN++fXz/+9/P137/13/9F6tWrSKZTPLWW29hWRb/8i//wtatW/ne975XzOYVBTWlHO/TI8esjaUg6WIOdiFTSidJa/ny5fz4xz+msbGRadOmsXDhQs4991xmzpyJ5518Jr+oolm4cGGfrgfgP//zP/s9tqqqqs+xkUgkv9rW6YxYyh9JfXgQHN8hNrZC72pBncaiOdEN/eKLL1JRUcGsWbNob2/n4osvZuYAK3KUYObe6KOmVuB9fBjjaj/QF7YwR5OYzgxSfnrWcxdyQw9EIJohILZCzqnBbD0Kccd3iAX0njasBYU7xH7CV+9HbiX1/M8htDE6cusbB6IZImp6FWZ7M0YbP6c4YmM+b8PMHZcrUvPAeBi8XCmLAeRYNYPJlbgYL/e7QYwc+8lXEtwLJFw1HhVy8BjidgAnIRDNEDDGQBg4qwz2dEDMzk0zZPH2H0adFUdUFGVVIiqGWHGURBAjSDYDmRSSSSEGMP4AVsTya6ZE5Rd1HA7dX7wDmZEZ0Z3Rojm+WzB5q2DQ/VuIXKWjqCj2OVNwv/wCseKICkHUQ+2JEzn/RpRSvgXJdkOyFbqbwcvkLiegosc28igxSlY02fRQAn/Sq7AthIjj/y52rujNQSzHf16FEVWGUhFEHESFEYkgKgwS9uNL1dBd9zKmvR2JRjERg25PYA42QbkFicM5oRhQIVDOsC3I6UDJisYKTS3wjB5L4mGMh983CJZdSyg6BytUj0jhd76zYAGpt97y16ExLsZNkXn3LaJ/MAeUPWaE0puSFY1Sw2+6MQbPbcHr+C1ileFEz8MOn1VQaa01eTJYgulqQcQgjuAdbkenNCpesh/vSSnNTrVIiAhKRVFWHLRLOvE2yfZf4rltA5+sXUgcRo5uI1QfxyTTfhenbBAhs2v/iLd/tChZ0ZhcN1PY48QxD1EOouJor51k+y/IJLf7CwH0fVFIJ6BlDzRthtY9oD1CZ0/O/TvnJEcc3H2H0N2pkXnzo0wJ2k8HcDFkBzzyOIybi5Xl4iESAnHy3ZGIIFKGMR6Zri142UNEYpcingvJNug6Cl7WD5uoUH5ZEVVuo6rKMYluCDuIUhgg+/l+wgtOj3ygYlJyoqmsXznkc/2VpdJorxMvexQ3sxfttmMERKK5dWSMH0rTIbzuvXQnviSqp6NM2BeK3f/GqM6sKaQ+3H4sHBdxyO45QGjmFFTZ2NoDs+REMxxEBCSCJWEsqxonNB2dbSGb+oJsdjdGu4i2OBaPVWjJkrQ+J8IcLE785VsTasC2MZ7n7wOl/MK6zM59RM6ffcLzSpHSE82uX/tBt/j4ws/tCd1rDz9gBwohbMI4zCSjjpCVoxhRCCEEQRHGkCXFjpxw+k/tFMsiNK2e7M59SMxfWUuiEdwvD6BnTkbFx85uc6XnCBsPP95C4Q8ExAIr7HczVhgsB+wwYscJq+lEZX5OKOn8OndCCFCkaMCj64RNC51dDyLHHOLcAknpbadXBuJwKT3R9JCfpxnK48RxGIsyoswjxMSccPyEJMEGhBQN6BNMAqqyCHb9OEzq2M57Eg3jHWzGax47NVKlK5oRxO+UphJmJqbXSE0IYTAk+QxN/8nkoZmTQZtj1kYEsS3SWz7rs35wKROI5iSEqCHKXEDyIlE4OR9nV94K9UZVlaOqKzDpXhWZYQedSJLd03Sqmj6iBKIZAIsYUeYiWHnhCA6aLtLsPW59XxHBmXMWeH2DiRJxyGzbje4qfn7LqSYQzSBQRIjytbxwBEEI49JMlsPHHW+Nr/ZHS5neq01YgJDe/FnJl/IGohkkx4SjMGRzwnHI8OVxIyoRwZk7DeO6x1kb72g72S9Ke15qSHGad999lzVr1vBHf/RHrFixAoCf/exnvPjii3znO9/hsssu63P8P/7jP+YL3+rq6njsscdOeO2Ojg4effRRDh8+zJNPPjmU5o0YiggRziHJDsDNjagsUuyijPm5v32sCbWoeBSTykLYz9f1t18Ok9m2G6umEqumIpcDbPzYUa+kL20M2tNo99geCYVg6uvQ1sjkCQ9JNJdccglXXHEFL7/8MpdddhmLFi3iT//0Tzlw4MBxggEYP348DzzwwKCuvWXLFq688kp+9rOfDaVpI45FjDAzSbMTPzQYQpMizZeEmY7kJhJECc78GaTe3QqOfWx+y1JgK1LvfkL0ygWoaNiPHTkxCEVxjU1Xa4buRBZjxN90THrnCw+Orp3NkNFDyC4emCF3T9FolEceeYTVq1fT0dGRf/7JJ5/kxhtv5JFHHuG6665j+/btdHV18dRTT/Gv//qvvP322wC0tLSwYsUKnn76ae655x4yuX22r7zySioqKob5tkaWEFU4TMWQwWBy/s1RXPquxWfV1fgTmelMLoHcf0gkisEi9fEBTO08mHwhZvwcOrM1HPhS095poe2YLyQ7irEiGCtc0MNLjZzfNCyf5rzzzuPP//zPuffee/PP3XHHHTQ3N3PXXXfx9NNPM3nyZP7kT/6Ev/3bv+Xb3/4269atY//+/ViWxe23385tt93GOeecw29/+9thv5lTSYgJ2NRgSOf9mzR7jgX+jEZ0hvDcqf5IyopAuAKiNRCpQFXWoju7SP3Pb/DSaZr3NNO6vxWxBbuXZTodGbYj/Dd/8zd0d3fz85//PP/c1KlTsSyLSZMmUVFRwbx587AsC8uymDNnDp988glKKX75y1/y1FNP8eGHH9LW1jbcppxS/PHTNBSR3IjKAgMp/RnG7fatSqwOa87F2LPnY7IGrFCfaLTE43gHD9L23Mt0HW7FClt+QvppzrBbKCI89NBDPP3007S3t+ef682jjz6a/72xsZEpU6bw85//HKUUy5cvZ/HixcNtxqggWESYBcZgdBZlLLTKkipLYCaeD9VnQziOc/HFABj3qyuCChlCmOZmYp++g+ru9I8zBi/r4aZd3Izr/yzwIcZAOoMqL/5E6ZAc4ffff5/333+fF154gRtuuIGamhoefPBBXnjhBV566aX85u3XXnstAEeOHGHdunVks1kuuOACFixYQCQS4YEHHmDt2rV88cUXiAhLlixh7969vPjii+zdu5ef/OQn/MVf/EVR33BRMRrlKRypJ2MdgFAVoiw8fYRMehvhsvMAUPE4zoUXkvngAyTnrxkg053G8wxSFkMyKZz3N5I5+xzciWfhxKPEx8UJx8JYoSFsVXj+VGTlDcV9vzlKbjN3dr7p/4xPHJ3XN35FA9r1u5rYeEysjlTqA7zsAZSKY4zG6ARO/GKciF9MbzyP5EsvoRMJVFkZ2VSWTDKDKIWQm4jPZrHcNKHKCpwLFxI6+2xUrEi7zRSRQDSDRXugs4D4aRXlE32n1vKNtdEputve8EdHKpLbmjBJuPwKQmG/3MZrbib58stoJ0K6O4soOSYYbbDDNk40BNksJu3PlEs4jJSXF2xpXtjcyP/zJmJFwjx39/8p4gdRiklYpwpjfGtivFxFpA3lk6CsGkJlx6VXiIoQKb+SZMevwLh+AZ6KkE78HlhMKDwZq7YWe+FCun77DhIpQ3JbmhqtcaIOdiTkx1UcB3EcP6jneZiOjoK3Pp3mdlKm6mgdgZXVA9FA32J8NP74wEAo5oskXAmh6IBFb1aoFif2B2QS74GK+4nrQDrxW+ByQuGpdNdMwauoxUq0YiJlGG1wyhxC4eOjtyLib788wO6+/TEyYT2f0hVNurPwc3r6gj5F9rmfVhic8tyjzLcmyir4JULhGRi3jWxqZ1/hdP6OTPI8Es0xrAUXoj7YCKkkTkV5v4I5nTn9gwInwujCH9rkovEmv7yHP98jvkBUKFdKa+XLUwpFRHBiC7Gcen8H3Z5SGRUh3fE+TmwXEg6TmvsHWJbCUqWXmFW6liZSWZzr9EwYumnIdIMcIr82TLgcotUQqfAt0SCdURGLSPnlJNvfQnttiMTwMuC5EUKR/XjZTtzaBcSXLSH9q19hRBDn9FxBqz9K19IUCxHfqig7l2we8QWibL+asnUvHPgEDn4CHU2QHVzVpEiIaMXXUaocz+0im8wiojA6irI6KB/3IdSmCf/hH2KyWUyqdKoxA9GcCFF+2N/OVSxoF9r3w8GP4dCnfrWlPvkqmKLCRCuvxsuGQZIgBqMhFK5AWQ7pxHtk49sIf+MSEEEnEictHT5dKL3uSSx/lOOlBz62P/JfihyzMqKO/d3vawqI7VsfYyCbhtbdfi13vM5/2P0X0nmuTceRc4mWf4Ky2hGrDDti+ykUykZ7HaSdj5D/XY1s6kIfbEdFoki4/0rO04HSE82sP4RUGxSc2Z8TS0+xnJcBNwNeyvdntO41usrFZfpzhkV8C9Rzrc5D/qOsGsrrj4vhdBzuAB0ik7wA29lKWWUbGD+Hxq8dj2KMwYS6MBe5mCYLvbUd0w4iNhKO+q9XYJ+gCo7sDJ7SEw1ApKq41+sJ5LkpyCYh3QGpDj8C3DNEzy0h0gdRfvdlTG6JtFY/B6Z8IkSqcF1N4mgCK2ThZTyi5ZcQjh8i0/2xv/CA8q1Jj3hQYKaWoSelkCPdmANJ9JEj0OUWnLk3aYKL3W44QaXNsChJ0TR2abJD2R5bQAnYArYSHOWvHx22wLJC/h0dLve7G5MbUaU7fEGkO/znROVWjOgloJ44jzG+6Jp3gbJIZ2PYWGgiiBIq66uwQrVYoXGkOn+P9hKIKuuzApeIhWXHoD4G9UP/jPZ/8SXu/6hANABrN2fJGpgaLyzimVtusQ/iL8eIMYJjQdyGyrBQE4YKR6hwwsRj45F4nW+JUh3+OnqZzl7r6PUKAOa7rhBae1ipI9REBa0FqRiH5XWBFccK1VFWdS2Z7k/JpnditOTEc/omXvWm5ESTyZnpeFGCqP6XZIwvnm4XOrKG3R2gxHdwFEJtRJhYJoyPVlM7rpqQzkKyGToP+taon+7Ly2g8Hcq5RZqItMPRXFqsE0cilYSdWYScs8ikduBmGnNLnkT8eavTmNO7dacIEcECLMtfMqnXf9DG0JoxHE6CiMEYqI0oziqvo76mjirdgXQ25QKDAiqEEcimsogiN7cURuwe59n4x6Y7QQRlIKKiaGsWWasN1xxB91RuivgVDkNYVzhsBY7wqKFEiFgQyfVCxhg6s/DREcNHQNgqZ3r5HKZXdFOVOYhKteO5HsZoRBSiBDvc62Pu1YXlLggYlKcJZytwiKMljSaJlm50fhGCwkRQF47iDGHubDAEoikQESGcc54BXG1oaIMdEiUkM5gVSzM+sY8aaUGZLHY45KdeniwGxLGotAAWUSyqhtXOluYmsl5PHVVxCUQzTGwlPXuh4mrDtvYQmdQURCZztnQwL95NdbadkNK+nyT2gMudnO4EoikithLCbgYLjTFwMFTDgVQdFoZKk2Cy1c4kaadM0lhKsMWgLNsfcg9xVn00CERTRIw2uFmP3MCLsqidXw0rY8rZacrZ6k3BNhmqJEkFCarooookIfEtkaVAof0AcM4BHopNCodsRiRIQyCaouLm9rY0Br/gTfWqcRIhJBBSAGGShElSxUFj8LRBGZeQzhLSGULGxcHFEQ8bL5cUWhj7iZKWA4yEcALRFAl/HjOLSE40kcF9tCKCbQng4OHgEaMYSRK/6ziKqw8W4UrHE4imSGhP59edsSyFZZWOj1IoY/ednWLcXgsY2ZHSyvktlMDSFAF/btPNj6Kt0MgE1U4XAktTBLxsLwc4HCrlEMygCERTBNxeBWm2M/aNdyCaYWK0wXN9S6MshbLGuJkhEM2wcXtlg4UGOcwudQLRDANjwE0dW2TaCgWiGTSbNm3i+uuvZ9WqVaxbt46VK1fS0NAwqHNbW1tZuXIlzz33HFu3buWb3/wm+/eXxpKp2tPoXGzGdqwx7wD3UBTRXHjhhcydO5elS5fyd3/3d/zlX/4l//Zv/zaoc6urq7n66qsBmD9/PnPmzClGk04JZ5oD3MOIdE8tLS1UVlZy4403sn//fjZs2MAtt9zCkSNH+MEPfsCTTz7JPffcw969e/s9/+WXX2bt2rV897vfpavrxFvljCbG+AG9nqQ6ZY/t2Exviiqa119/nXXr1vHJJ59wxx13kM1mmTx5Mo7j8O1vf5vx48ezbNky7rjjDm699VZ+8pOf9HudCy64gLvuuotzzz2Xn/70p8VsYtHoiQCfKbGZ3hTVpi5dupQlS5bk/16wYAHvvPMO77zzDqtWrQJg165dbNmyhWQyecIVPSdP9nefnTJlCu+++24xm1gUjPFzgHuK686krglGePR000038e///u/EYjGUUmzbto0333yTFStWcPPNN5/wvKYmf4ubffv2MX369JFs4pDQrudPTpozJzbTm6LcIps3b6ahoYFMJsPs2bOZOdNfnHDRokX8wz/8A/fccw8A06dPp6ysjDVr1iAiNDQ08OGHH/LWW28BcM4559DQ0MDGjRvZuHEje/fu5aGHHipGE4tKNnXMAQ6FzywrAyW4UOPDH/lxkXnVo3N3a8+Q7Ejm82bKqspOS3/mdzuPsmPzp6RTad74v9cV9dpBcK9Asrkd44zxrczpKJiRJhBNAWhtcikQvlLsElsrr1gEoimAnmCeMQbLPvMc4B4C0QwSo41fattjZcZ4dt7JCEQzSI6NmAyiBOsMigB/lUA0g0Brk6s0EN8BjpxZEeCvEohmEGSTPekPBpEzLwL8VQLRDIDnaX9iUvVYGeeMtjIQiOakGAPZrkz+DxH6LhtyhhJ8AifBy7h4rodYyl+cKFq6VmbTpk3cf//9TJs2jVmzZgH+Zm/Lly/njTfeoK6ujp07d7JkyZL85m4nIhDNCTDa5DbxkjFhZXoS5a666qp8JsKGDRt49NFHefbZZ4nH4xw8eJBf/OIXA14r6J76wRhId2f8RapEMMYQKmEr0x/r16/nggsuoKamhmeeeYa2tjYmTpzIX/3VXw14biCafnAzLl7G851fbRClxsyI6fnnn+ef/umf8pkFDz30EEePHmXp0qXcdtttfP755wNeIxDNV/A8TaY7c2yNIQPhsrFjZb75zW9yzz338M///M+UlZVRV1fHgw8+yIYNG/KpLAMRiKYXRhvSidyeC+JbGduxxmRt9tSpU7Ftm/vuuw8Ax3G49tprcY/bBvp4xobNLQLGQCqRznVHkt94I1RWOvswnYzeiXIzZszgnHPOAWD8+PHce++9jBs3jt27d3P33XcPeK1ANOQc30QK7en86lXGgFPmoNTY6JcWLlzI888/f9zzd955Z8HXKjnROAJZA4nswMcOhDZ+nq+XyiDaYOUWg8YYrJA1ZpzfYlNyn8pdC0O0pg3ZIWz9aAxoDJ4GV0OiK8ORQwm6jZASm6RYdGtBYXAch0zW4CgIWWCNFU+4CJScaACqw0P/Av1V5jN0HO1A2rupsiyUrYA02tNktSY+YyIpy6ItbTiUhOaUwTUGjKDEELHAGiPd1lAoOdF89pvPMMYwcc7Egs910y6pzpS/CFFutronqcpog3Y1ddNqiVX6zu/UOCzg2NL2LSlDY5f/6HZ9R9mxwFGUzA4qxaDkRNOzGGLbgbaCzhMkt5OywnKsPl+y1hov41EztYZYdez4c0WocPztfKZV+CJqy8CBbs3uDkNbGgzG3zvqDBBQyYmmh2JtgK5djfY0tWfVEh8XH9Q5IkJ1GKrDFvOqoStr2N+l+aLD0JLCFydjtxsrWdEMF2MMXtafKqibWUekov+NSwdDLCTMqbKYUwVJ13A4adiXMDR1GZKewZBzri1/V7tSt0RnnGh6xIKBaGWU6inVRR1aR23h7HLh7HL/tTpyvtCBbl9MiWzPvlF+IbitwBJyy9uXhqBKVjRdzYUvQeLf84KyFVbIItmeJNmeHIHW9cUCpuQeLtCNRTcWXVgkcr8nUHi+53XcXgj9baU4EBPLwmwfftP7pWRFM5TuREQQS0b1braBCFCDBjRwLEqpDWQQXATPCB6gh7QzAuxsbyGbDjbU6MNYnERUAhH8neaGtPVKLxqzxdhhoX+CWe6AgglEE1AwgWgCCiYQTUDBBKIJKJhANAEFE4gmoGAC0QQUTCCagIIJRBNQMKMyjbB27VqMMYgI27dv5+GHH6a2trbPMfv372fHjh184xvfGI0mBpyEUy6aHTt2sHv3bp544gkAnn32WdLp9HHHNTY28qtf/SoQzWnIKe+eKisr2bJlCxs2bMDzPG699VYcx2H16tWsX7+e7373uwC88sorbN++nSeeeOKEeygEjA6n3NLU19ezZs0a1q9fz+rVq7npppuYN28e4XCYb33rW2zZsgWAZcuWkc1mWbly5aluYsAAjIojvHjxYv7jP/6D559/nvfff5/Ozk5qamq4+eab+fWvf43WQyhqCjhlnHLRNDY2sn79egAmTJjApZdeSmtrK3/913/NT3/6UzZv3sy2bduwLMvP+m9r49ChQ6e6mQEn4ZR3T/F4nI8//pjHH38crTUHDx7k5ptvZu3atdTU1DBp0iRmz55Nd3c3TU1NPPLII3zrW9861c0MOAkltwtLw//4G6qW15WPcktObzbvPMwTnzRhjAl2YQkYfQLRBBRMIJqAgglEE1AwgWgCCiYQTUDBBKIJKJhANAEFE4gmoGAC0QQUTCCagIIJRBNQMIFoAgomEE1AwQSiCSiYQDQBBROIJqBgAtEEFEzJLdSoLIX2NG564B3QzmS6XY0xhrIR2OG35HKEA0afoHsKKJhANAEFE4gmoGAC0QQUTCCagIIJRBNQMIFoAgomEE1AwQSiCSiYQDQBBROIJqBgAtEEFEwgmoCCCUQTUDD/H4egt+j2lTGYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 108x252 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sankey plot mapping crelines to FS/RS\n",
    "df = pd.DataFrame(df_model['cre_lines'].values[::-1],columns=['cell_type'])\n",
    "df.insert(0, column=\"label\",value=y_pred_str[::-1])\n",
    "fig=sankey(df['cell_type'],df['label'],fontsize=9,fig_size=(1.5,3.5),colorDict=cre_pal)"
   ]
  }
 ],
 "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
}