VNSGM_function_v4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HSfbgraphcore.txt       VN_via_SGM_2.ipynb      VN_via_SGM_6.ipynb\r\n",
      "HSfrgraphcore.txt       VN_via_SGM_3.ipynb      ego.py\r\n",
      "VNSGM_function_v3.ipynb VN_via_SGM_4.ipynb      ego_testing.ipynb\r\n",
      "VN_via_SGM.ipynb        VN_via_SGM_5.ipynb      fig6_data.npy\r\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from graspologic.plot import heatmap\n",
    "import igraph\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from graspologic.match import GraphMatch as GMP\n",
    "from graspologic.utils import pass_to_ranks\n",
    "from scipy.stats import rankdata\n",
    "import pyreadr\n",
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helper Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "def permute_adj_matrix(adj_mat):\n",
    "    assert adj_mat.shape[0] == adj_mat.shape[1]\n",
    "    indexs = np.random.permutation(np.arange(adj_mat.shape[0]))\n",
    "    ress = np.zeros(adj_mat.shape)\n",
    "    for ii, ind in enumerate(indexs):\n",
    "        ress[ii] = adj_mat[ind][indexs]\n",
    "    return ress, indexs\n",
    "\n",
    "def ego(graph_adj_matrix, order, node, mindist=1):\n",
    "    #Note all nodes are zero based in this implementation, i.e the first node is 0\n",
    "    dists = [[node]]\n",
    "    for ii in range(1, order+1):\n",
    "        clst = []\n",
    "        for nn in dists[-1]:\n",
    "            clst.extend(list(np.where(graph_adj_matrix[nn]==1)[0]))\n",
    "        clst = list(set(clst))\n",
    "        \n",
    "        #Remove all the ones that are closer (i.e thtat have already been included)\n",
    "        dists_conglom = []\n",
    "        for dd in dists:\n",
    "            dists_conglom.extend(dd)\n",
    "        dists_conglom = list(set(dists_conglom))\n",
    "        \n",
    "        cn_proc = []\n",
    "        for cn in clst:\n",
    "            if cn not in dists_conglom:\n",
    "                cn_proc.append(cn)\n",
    "        \n",
    "        dists.append(cn_proc)\n",
    "    ress = []\n",
    "        \n",
    "    for ii in range(mindist, order+1):\n",
    "        ress.extend(dists[ii])\n",
    "        \n",
    "    return list(set(ress))  \n",
    "\n",
    "def ego_list(graph_adj_matrix, order, node, mindist=1):\n",
    "    print(type(node))\n",
    "    if type(node) == list:\n",
    "        total_res = []\n",
    "        for nn in node:\n",
    "            ego_res = ego(graph_adj_matrix, order, nn, mindist=1)\n",
    "            total_res.extend(ego_res)\n",
    "        return list(set(total_res))\n",
    "    else:\n",
    "        return ego(graph_adj_matrix, order, node)\n",
    "\n",
    "def shuffle_adj_matrix(adj_mat, shuffle_lst):\n",
    "    assert adj_mat.shape[0] == adj_mat.shape[1]\n",
    "    ress = np.zeros((len(shuffle_lst), len(shuffle_lst)))\n",
    "    for ii, ind in enumerate(shuffle_lst):\n",
    "        ress[ii] = adj_mat[ind][shuffle_lst]\n",
    "    return ress\n",
    "\n",
    "def pass2ranksuplus(matt):\n",
    "    out_matt = np.zeros(matt.shape)\n",
    "    for rr in range(matt.shape[0]):\n",
    "        out_matt[rr] = rankdata(matt[rr])\n",
    "    return out_matt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation of VN via SGM algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vnsgm(x, seeds, A, B, h, ell, R, g, pad=0, sim=False, verb=False, plotF=False):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    x: int or ndarray\n",
    "        verticies of interest (voi)\n",
    "    \n",
    "    seeds: ndarray\n",
    "        seeds vector\n",
    "        \n",
    "    A: ndarray\n",
    "        Adjacency matrix of `G_1`, the graph where voi is known\n",
    "        \n",
    "    B: ndarray\n",
    "        Adjacency matrix of `G_2`, the graph where voi is not known\n",
    "        \n",
    "    h: int\n",
    "        distance between voi used to create induced subgraph on `G_1`\n",
    "        \n",
    "    ell: int \n",
    "        distance from seeds to other verts to create induced subgraphs on `G_1`\n",
    "        \n",
    "    R: int \n",
    "        number of restarts\n",
    "    \n",
    "    g: float\n",
    "        gamma to be used, max tol for alpha, tollerable dist from barycenter\n",
    "        \n",
    "    \n",
    "    \n",
    "    '''\n",
    "    \n",
    "    nv1 = A.shape[0]\n",
    "    nv2 = A.shape[0]\n",
    "    \n",
    "    nv = max(nv1, nv2)\n",
    "    \n",
    "    nsx1 = set(range(nv1)).difference(set([seeds[0], x]))\n",
    "    \n",
    "    vec = [seeds[0], x]; vec.extend(nsx1)\n",
    "    \n",
    "    AA = shuffle_adj_matrix(A, vec)\n",
    "#     print(AA[:, -10:])\n",
    "\n",
    "\n",
    "    ns2 = set(range(nv2)).difference(set([seeds[1]])) # <- setdiff(1:nv2,seeds[,2])\n",
    "    vec2 = [seeds[1]]; vec2.extend(ns2) # <- c(seeds[,2],ns2)\n",
    "    \n",
    "#     print(B[:, -10:])\n",
    "\n",
    "    BB = shuffle_adj_matrix(B, vec2)\n",
    "#     print(BB[:, -10:])\n",
    "\n",
    "    nrow_seeds = seeds.shape[0]\n",
    "\n",
    "    S = list(range(nrow_seeds))\n",
    "    \n",
    "    voi = list(range(nrow_seeds, nrow_seeds+2))\n",
    "    \n",
    "    #print(\"inside 1 w/ shapes \", AA.shape, BB.shape)\n",
    "    \n",
    "    P = vnsgm_ordered(voi,S,AA,BB,h,ell,R,g,pad=pad,sim=sim,verb=verb,plotF=plotF)\n",
    "    \n",
    "    P['x'] = x\n",
    "    P['S'] = seeds\n",
    "#     print(seeds, P['Sx'],P['Sxp'])\n",
    "#     P['Sx']  = seeds[P['Sx'],0]\n",
    "#     P['Sxp'] = seeds[P['Sxp'],1]\n",
    "#     print('vec2 = ', vec2)\n",
    "#     print('vec = ', vec)\n",
    "#     print('Cxp = ', P['Cxp'])\n",
    "#     print('labelsGx = ', P['labelsGx'])\n",
    "    \n",
    "    if P['Cxp'] is not None:   \n",
    "        foo = []\n",
    "        for kk in P['Cxp']:\n",
    "            foo.append(vec2[kk])\n",
    "        P['Cxp'] = foo\n",
    "        \n",
    "    if P['labelsGx'] is not None:\n",
    "        foo = []\n",
    "        for kk in P['labelsGx']:\n",
    "            foo.append(vec[kk])\n",
    "        P['labelsGx'] = foo\n",
    "    \n",
    "    if P['labelsGxp'] is not None:\n",
    "        foo = []\n",
    "        for kk in P['labelsGxp']:\n",
    "            foo.append(vec2[kk])\n",
    "        P['labelsGxp'] = foo \n",
    "    return P\n",
    "\n",
    "def vnsgm_ordered(x,S,A,B,h,ell,R,g,pad=0,sim=True,verb=False,plotF=False):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    x: list \n",
    "        Verticies of intrerest\n",
    "        \n",
    "    S: list\n",
    "        Vector of seeds\n",
    "        \n",
    "    A: ndarray\n",
    "        `G_1` where voi are known\n",
    "    \n",
    "    B: ndarray\n",
    "        `G_2` where voi unknown\n",
    "        \n",
    "    rest same as before\n",
    "    '''\n",
    "    s = len(S)\n",
    "    Nh = ego_list(A, h, x, mindist=1) # TODO: make able to take a list in <- unlist(ego(g1,h,nodes=x,mindist=1))\n",
    "    \n",
    "    Sx1 = None; Sx2 = None\n",
    "    \n",
    "    Sx1 = list(set(Nh).intersection(set(S))); Sx1.sort(); sx = len(Sx1)\n",
    "    \n",
    "    Sx2 = list(set(list(range(B.shape[0]))).intersection(set(Sx1))); Sx2.sort(); sx2 = len(Sx2)\n",
    "    #print(\"SX1, Sx2, \", Sx1, Sx2)\n",
    "               \n",
    "    case = 'possible' if sx2>0 else 'impossible'\n",
    "    \n",
    "    if case == 'possible':\n",
    "        Cx2 = ego_list(B, ell, Sx2, mindist=1)\n",
    "        Cx2 = set(Cx2).difference(set(Sx2))\n",
    "        print(Cx2)\n",
    "        if sim:\n",
    "            case = np.array(['possible' if _x in Cx2 else 'impossible' for _x in x])\n",
    "            \n",
    "        \n",
    "        if len(np.where(case == 'possible')[0])>0:\n",
    "            Nx1 = ego_list(A, ell, Sx1, mindist=0)\n",
    "            Nx2 = ego_list(B, ell, Sx2, mindist=0)\n",
    "            \n",
    "            #print(\"NX1, Nx2, \",len(Nx1), len(Nx2), Nx1, Nx2)\n",
    "            \n",
    "            if sim:\n",
    "                #print(\"in opt 1\")\n",
    "                wxp = np.where(case == 'possible')[0]\n",
    "                print(wxp)\n",
    "                xp = x#x[wxp]\n",
    "                \n",
    "                print(Sx1, x, Nx1, np.setdiff1d(Nx1, np.concatenate((Sx1, x))))\n",
    "                ind1 = np.concatenate((Sx1, x, np.setdiff1d(Nx1, np.concatenate((Sx1, x)))))\n",
    "                print(Sx2, xp, Nx2)#, np.setdiff1d(Nx2, np.concatenate((Sx2, xp))))\n",
    "\n",
    "\n",
    "                ind2 = np.concatenate((Sx2, xp, np.setdiff1d(Nx2, np.concatenate((Sx2, xp)))))\n",
    "            else:\n",
    "                #print(\"in opt 2\")\n",
    "                ind1 = Sx1 + x + list(set(Nx1).difference(set(Sx1+x)))\n",
    "                \n",
    "                ind2 = Sx2 + list(set(Nx2).difference(set(Sx2)))\n",
    "                \n",
    "            if verb:\n",
    "                print(\"seed = \", Sx1, \", matching \", ind1, \" and \", ind2)\n",
    "                \n",
    "            #print(\"inside 2 w/ params Sx1 = \", Sx1, \", x = \", x, \", Nx1 = \", Nx1)\n",
    "            #print(\"inside 2 w/ indlens \", len(ind1), len(ind2), ind1, ind2)\n",
    "            n_iters = 20\n",
    "            As = shuffle_adj_matrix(A, ind1)\n",
    "            Bs = shuffle_adj_matrix(B, ind2)\n",
    "            \n",
    "            sgm = GMP(n_init=n_iters)\n",
    "            \n",
    "            # Note we must padd because the inputs are note guarenteed to be the same\n",
    "            totv1 = As.shape[1]\n",
    "            totv2 = Bs.shape[1]\n",
    "\n",
    "            if (totv1>totv2):\n",
    "                #print(\"in pad 2\")\n",
    "                diff = totv1-totv2\n",
    "                Bs = np.concatenate((Bs, np.full((Bs.shape[0], diff), pad)), axis=1)\n",
    "                Bs = np.concatenate((Bs, np.full((diff, Bs.shape[1]), pad)), axis=0)\n",
    "            elif (totv1<totv2):\n",
    "                #print(\"in pad 2\")\n",
    "                diff = totv2-totv1\n",
    "                As = np.concatenate((As, np.full((As.shape[0], diff), pad)), axis=1)\n",
    "                As = np.concatenate((As, np.full((diff, As.shape[1]), pad)), axis=0)\n",
    "        \n",
    "            #print(\"inside 2 w/ shapes \", As.shape, Bs.shape)\n",
    "            \n",
    "            corr = sgm.fit_predict(As, Bs, seeds_A=np.arange(len(Sx1)), seeds_B=np.arange(len(Sx1))) #multistart_sgm(A, B, R, len(Sx1), g, pad=pad, n_iters=n_iters)\n",
    "#             P_outp = sgm.P_final\n",
    "            Cx2 = list(Cx2)\n",
    "        else:\n",
    "            ind1 = ind2 = P_outp = As = Bs = corr = None\n",
    "    else:\n",
    "        ind1 = ind2 = P_outp = Cx2 = As = Bs = corr = None\n",
    "                \n",
    "    out_dic = {'case':case, 'x':x, 'S':S, 'Sx':np.array(Sx1), 'Sxp':np.array(Sx2), \n",
    "               'Cxp':Cx2, 'labelsGx':ind1, \n",
    "               'labelsGxp':ind2, 'As':As,\n",
    "                'Bs':Bs, 'corr':corr} #'P':np.array(P_outp),\n",
    "#     print(\"*\"*50)\n",
    "#     print(out_dic['case'])\n",
    "#     if P is not None:\n",
    "#         print(P)\n",
    "#         print(sgm.score_)\n",
    "#         print(sgm.P_final)\n",
    "    return out_dic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4])"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.concatenate(([1,2,3], [4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial Reproduction\n",
    "https://github.com/dtborders/graspologic.git\n",
    "\n",
    "## Goal: Compare outputs to R output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 4\n",
    "h = 1\n",
    "ell = 1\n",
    "R = 100\n",
    "gamma = 0.1\n",
    "mc = 2\n",
    "\n",
    "np.random.seed(seed=1234+mc)\n",
    "\n",
    "HSfbgraphfull = igraph.load(\"HSfbgraphfull.txt\", format=\"edges\")\n",
    "HSfbgraphfull.to_undirected()\n",
    "\n",
    "HSfriendsgraphfull = igraph.load(\"HSfriendsgraphfull.txt\", format=\"edges\")\n",
    "HSfriendsgraphfull.to_undirected()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "coremap = pyreadr.read_r('coremap.rds')\n",
    "coremap = coremap[None]\n",
    "coremap = np.array(coremap).astype(int)\n",
    "coremap = np.stack((coremap[0:82, 0],coremap[82:, 0]), axis=1)\n",
    "coremap = coremap-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(82, 2) [[ -1  -1]\n",
      " [  0   0]\n",
      " [  4   4]\n",
      " [  5   5]\n",
      " [  6   9]\n",
      " [  8  11]\n",
      " [ 10  12]\n",
      " [ 11  15]\n",
      " [ 12  16]\n",
      " [ 13  18]\n",
      " [ 17  19]\n",
      " [ 18  20]\n",
      " [ 19  22]\n",
      " [ 20  23]\n",
      " [ 21  24]\n",
      " [ 22  25]\n",
      " [ 24  27]\n",
      " [ 26  28]\n",
      " [ 27  29]\n",
      " [ 28  30]\n",
      " [ 29  32]\n",
      " [ 33  34]\n",
      " [ 34  35]\n",
      " [ 36  36]\n",
      " [ 37  38]\n",
      " [ 38  39]\n",
      " [ 39  40]\n",
      " [ 40  42]\n",
      " [ 41  43]\n",
      " [ 42  44]\n",
      " [ 43  46]\n",
      " [ 44  47]\n",
      " [ 45  49]\n",
      " [ 47  54]\n",
      " [ 48  56]\n",
      " [ 50  57]\n",
      " [ 53  61]\n",
      " [ 54  63]\n",
      " [ 56  64]\n",
      " [ 58  65]\n",
      " [ 59  66]\n",
      " [ 63  70]\n",
      " [ 65  71]\n",
      " [ 66  72]\n",
      " [ 67  73]\n",
      " [ 68  74]\n",
      " [ 70  75]\n",
      " [ 73  78]\n",
      " [ 76  81]\n",
      " [ 77  82]\n",
      " [ 78  84]\n",
      " [ 79  85]\n",
      " [ 80  87]\n",
      " [ 83  88]\n",
      " [ 86  90]\n",
      " [ 87  91]\n",
      " [ 88  92]\n",
      " [ 89  93]\n",
      " [ 91  94]\n",
      " [ 92  95]\n",
      " [ 93  97]\n",
      " [ 96 100]\n",
      " [ 97 101]\n",
      " [100 102]\n",
      " [105 104]\n",
      " [106 105]\n",
      " [108 106]\n",
      " [110 108]\n",
      " [112 110]\n",
      " [113 111]\n",
      " [117 114]\n",
      " [119 115]\n",
      " [120 117]\n",
      " [122 118]\n",
      " [123 119]\n",
      " [130 122]\n",
      " [131 124]\n",
      " [140 126]\n",
      " [145 127]\n",
      " [148 128]\n",
      " [151 131]\n",
      " [153 132]]\n"
     ]
    }
   ],
   "source": [
    "print(coremap.shape, coremap-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(156,) [  0   1   5   6   7   9  11  12  13  14  18  19  20  21  22  23  25  27\n",
      "  28  29  30  34  35  37  38  39  40  41  42  43  44  45  46  48  49  51\n",
      "  54  55  57  59  60  64  66  67  68  69  71  74  77  78  79  80  81  84\n",
      "  87  88  89  90  92  93  94  97  98 101 106 107 109 111 113 114 118 120\n",
      " 121 123 124 131 132 141 146 149 152 154   2   3   4   8  10  15  16  17\n",
      "  24  26  31  32  33  36  47  50  52  53  56  58  61  62  63  65  70  72\n",
      "  73  75  76  82  83  85  86  91  95  96  99 100 102 103 104 105 108 110\n",
      " 112 115 116 117 119 122 125 126 127 128 129 130 133 134 135 136 137 138\n",
      " 139 140 142 143 144 145 147 148 150 151 153 155]\n",
      "(134,) [  0   1   5   6  10  12  13  16  17  19  20  21  23  24  25  26  28  29\n",
      "  30  31  33  35  36  37  39  40  41  43  44  45  47  48  50  55  57  58\n",
      "  62  64  65  66  67  71  72  73  74  75  76  79  82  83  85  86  88  89\n",
      "  91  92  93  94  95  96  98 101 102 103 105 106 107 109 111 112 115 116\n",
      " 118 119 120 123 125 127 128 129 132 133   2   3   4   7   8   9  11  14\n",
      "  15  18  22  27  32  34  38  42  46  49  51  52  53  54  56  59  60  61\n",
      "  63  68  69  70  77  78  80  81  84  87  90  97  99 100 104 108 110 113\n",
      " 114 117 121 122 124 126 130 131]\n"
     ]
    }
   ],
   "source": [
    "nnodes_fb = g1.shape[0]\n",
    "nnodes_fr = g2.shape[0]\n",
    "perm_fb = np.concatenate((coremap[:, 0], np.setdiff1d(np.arange(nnodes_fb), coremap[:, 0])), axis=0)\n",
    "perm_fr = np.concatenate((coremap[:, 1], np.setdiff1d(np.arange(nnodes_fr), coremap[:, 1])), axis=0)\n",
    "print(perm_fb.shape, perm_fb)\n",
    "print(perm_fr.shape, perm_fr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "HSfbgraphfull = HSfbgraphfull.permute_vertices(list(perm_fb))\n",
    "HSfriendsgraphfull = HSfriendsgraphfull.permute_vertices(list(perm_fr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "adjfb = np.array(HSfbgraphfull.get_adjacency().data)\n",
    "adjfr = np.array(HSfriendsgraphfull.get_adjacency().data)\n",
    "# adjfr[np.where(adjfr>=1)]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f98ef0c2f50>"
      ]
     },
     "execution_count": 189,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAHrCAYAAADyu/bLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3defwkVXnv8e8zoyJIFBmDuL5GDCooGjRxiwh6URMNiOCShKvgeseICUbjEvEG1BDXKIhXohI1iWvUQVHk6lVGXHCLiSiD4MIoRATCKjCoyLl/dP+w/U11V53qU+c5VfV5v17z+kF3dZ1TS3effp6zWAhBAAAAHtZ4VwAAAIwXDREAAOCGhggAAHBDQwQAALihIQIAANzQEAEAAG5oiABTZva7ZvZZM7vCzIKZHW1m+03/+3Dv+qXS5pjM7M/N7Ltm9vPpa9dHlrl+5ZyuejyY2btj9gVgWGiIoBdmvjxftGCbYGafmPn/tWb2VDP7opn91MyuN7MLzex0M3ulmW03s+3NJH1E0u6SXiHpqZI+2uEh9YaZPULSWyV9V9IGTc7Npa6VAjAYN/OuANCh90l6sqQvSXqjpCsk3VXSAyX9jaTjJf18uu1u038vDCGcsLIDM9svY31L9ajp32eEEC53rQmAwaEhgkEyswdo0gjZGEI4uOL520u6auahXad/Xb9ozWytpO1CCNd51mOVXSWJRgiALpCawVDtPv37uaonQwgXhxB+KUlmtknS56dPvWua4tmmH4SZPd/MzpumeM4zs+dX7dvMDjGzb023+7GZ/a2Z7b+6X4aZHT59bH8ze4WZ/UDS9Zo0oGRmjzazD5rZD81sq5ldaWafNrN9K8rcZGZbzGw3M/uYmV1lZleb2UYz223eSTKzp5vZ2dO+Hz8ysxfPPLfezIKkp0//f+W8bJots2Kflf1BAKAKERH0zQ5mdrsG2/1g+vdJZvbeEMIVC7b9O03SN38j6e2SvjB9/FJJ66f//XxNIgP/KOlnkv5U0vFmtnMI4ZiVHZnZUyS9f1r+MZJukHSYpAMWlP8GSTeX9A5JV0s6d/r44ZJ2lvTPki6UdCdJz5L0WTN7RAjhC6v2cytJp0v6mqSXadIY+3NJDzazvUMIP121/QZJt5d0kqQrJf1PSa81swtDCO+bHv9TJT1H0j7T/5akixccCwBEoSGCvjlm+m+hEMLXzewUTRoAF5rZlyV9dfrvs7OpjxDCZ8zsl5o0RM4MIfzrynNmtvKf95C0Rwjhwunjb5X0RUlHmdlJIYQLpx1e/0GTL/AHrjR+zOxtks5aUN3tJe1dkY55dgjh2tkHzOxESWdr0tBY3RC5naTjQghHzmx/hiadbo/WpOEx666S9gwhXDnd9p8k/UiTRtf7pmX/q5ntL2mf2fMCAKmQmkHfvF2TzpNV/1Y7RJMv1e9I2k/SyyV9XNJPzeyFkeW+d6URIkkhhF9IepMmjfmVaMcDJN1R0rtnIzAhhGsknbhg32+r6hMy2wgxsx3NbJ2kX2nSmHrQnH29ZtU+NmoSYTmoYtt3rTRCptteJ+kr+nVaCwA6R0QEffO9EML/q3piJnohSZr2ATlB0glmtr0mDYXHatI4eYOZ/SSE8P6G5Z5T8djm6d+VPhh3m/49t2LbqsdWnFf1oJndXZO00WMk7bTq6VDxkisr0i/SpO4HmdmtVkVYflix7WWS1i2oKwAkRUMEoxBC2KpJKuWLZna6pE9LeqYm/Tka7aLiMav5/6a2iYaY2Y6SztCk38ebJX1bk74pN2qSlnlkwzouqtevomvarDw+VwA0xgcGxugr0793injNnhWP7TH9uxJZOH/6954V21Y9tsj/0CTN84wQwrtmnzCzV895zW3NbNeKqMi9JF2yur9JApdrEmVabe4oHQBYjT4iGCQz293MfmfO0yv9JTbPeb7KoWZ255n930LSCzSJKqzM5voNSRdJOtzMbjuz7Y7atqNonZVoxW9EM8zs0ZrfP0SSXrpq+ydo0gg6ObL8Js6T9Ftm9sCZ8tZocl4AoBEiIhiq+0n6oJl9XtImTYa/3kqTL/Ena5LmeGXE/s6T9NXpqJWfSfozSb8v6VUhhAskKYRww3QK+vdK+pqZnaTJ8N3DNel7cTfNT2es9kVJP5X0xul8JhdK+l1NhtB+W9JeFa/5b0kHm9kdp8e8Mnz3Yk1GzaT2dkkvlLTRzI6T9AtJTxSfKwAi8IGBoTpD0l9rMprmGZrMl2GSLpD0LkmvDyF8P2J/b5F0a006ut5V0o8lHRlCOG52oxDC+8zsBklHaTLM+GJN5uk4S5NhtFubFBZCuNLMHiPpddMybybp3zXpbPtMVTdErtWk78ibNBk9Y5JO02Ta+osijrWREML5ZnaQpGMlvUqTxta/SPonTdalAYBaFkLTH2gA2poOF36DpIeEEL5St32L/W+StD6EsD71vgGgS/QRARIys1tM14uZfWxHSc/TJGLwTZeKAUChSM0Aae0m6VNm9gFNRtHcQZMp3u8m6bnTidAAAFM0RIC0LtVkePChknbRpLPqtyW9NITwIc+KAUCJ6CMCAADc0EcEAAC4oSECAADc0BABAABuaIgAAAA3NEQAAIAbGiIAAMANDREAAOCGhggAAHBDQwQAALihIQIAANzQEAEAAG4WLnp3/datNy1Ec+QOe2zz/JuvO6eDKm1rtux5Za5sk6pOKY63ah+z+2lyXDHlzO5jmX2Xfq3bvi7F+U51zer2PStFOTH1TnXvVL22y+vU5Pylvn5d3g9V5aT+fMv1vu5SF9fglttvb0l2hEaIiAAAADc0RAAAgJuFqZmqkNe88GeXPMKHKcqs20eq40p9fvp4rXPdI12WMy+FkHrfKbdN+dpl95EjNbK6nCGkNvqOa9B/REQAAIAbGiIAAMDNwtTMrBy9rOeFP5v0iu5Tj/XV++1i38vwSMksy/seSaGU8+59X3Y5cqpq+9LvC6m7Onpfa0AiIgIAABxZCGHuk7PziKSwTMSjZB7173KuhS71/VrPGtKxNOF9n8/Tx3Pf5RxCufTp/m/y3TPrxLCFeUQyco+I9OlmRjpca6Qw5vuolFReX3H+yuHeEAEAAOPVeWrGM+KxTBoiZprz0qeN7nJ/s3J1EMxxfkqZbr3N9su+riSeyw2kkOtap95HW128b2KWsEhVPlO850VEBABQJNIn40BDBAAAuMk6aqYt7xBzn+YbKJH39eujXOk5rk0/9P0zKPcqxcuWQ2omLyIiAADADQ0RAADgJllqpvQQb+n1m1V6XUuvX2pjO15g7EjN5EVEBAAAuKEhAgAA3CyVmun75Dt91ve1ZoZwrbl3u1c3j0Rfz/sQPju9y09p9X3GWjN5EREBAABuaIgAAAA3vZjQLEbq5bVn95NrUqguyyklnOqxBkdbTcqOSXl5r/2TO93Rp8nUSqpTKXUZYz0YNZMXEREAAOCGhggAAHCTddRMjtCaRxixT2X2ddRMl/XOfS5LCXX3wRhGzcxqe//1adRMTKoz17o0qzFqJi8iIgAAwA0NEQAARs7M7mxmx5nZF83sGjMLZrZfxOsfYGafNbNrzewKM/uAmd2p0WtjUzN9DQfG7iPHyIZU56/LfbetRym962d5pnFi91f3+BBSEh5h97bnsi5F1HQ/bfX9uvfNGEfNTBsdH5L0TUk/l3SgpEeEEDY1eO0ekr4m6euSXi/pVpL+TpNgx94hhGsWvf5my1QcAAAMwhkhhF0kycwO0qQh0tQxkn4m6YAQwrXTfXxH0tmSnifptYtePLh5RLpUyi/+ZXjP6TDUzqp96mTYV6VE/lCuVO+tMUZEZk0bIhvVICJiZjeXdJWkd4YQ/mLVc1+SdPMQwgMX7YOICAAAA2FmV9ZtE0LYKWGRu0naXtJ3Kp47S9JhdTugIQIAgKMNtj5l9uGqhPtqYt307+UVz10uaXsz2z6EsHXeDhY2RGKmOa97Xde67MxVte8U58HjXM7bx7xjS1F2zDE02TZ1XdseW0wnx9lt20693rTMun13Oc9J23s+9rpXvW4ej+UYUqOzajOcn+TRjqiiWz7H8F0AANDaZdO/6yqe21nS1hDC9Yt2QEMEAAC09UNJWyXdp+K5vVTdd+Q3FDVqZiyjCQizoi3unWa8lzJAv+UeNfO8Nen6iLz1xuWnp48ZNTPd/t8kPVTS7iGE66aP3UPSOZJeHkJ4zaLX01kVAADIzJ44/c/fn/7d18xuJ+naEMKnpttskaQQwvqZl/6tJhOafdzM3qBfT2i2RdJb68qlIQIAACTp31b9/9HTvz+StH7ei0IIm83sEZpMXPYRSb+U9GlJLwwh/Kyu0KJSMyUZS5qoNG1XdOYalaGU900p9UA/5U7N/MXauyX7rj3+V+f3bjI2OqsCAAA3NEQAAICbxhOarejrqq7L9KJPEf6vm3gp9XlNMSlVn0LaMautrn68FKlXjO1Sk4nEcvCoR4nXiVFC/ba2d8mUtIiIAAAANzREAACAG7dRM6WkfWIss05GKeXkXhsjVTldhp7b1rXt2jAlrXeS435IdT/n+MwoKX1XSl1i1sTq8v7L+b7JPWrmRTffLdl37Rt++cPeJXqIiAAAADc0RAAAgJtkqZlSwohDUFLoHul1eT241vE4Z+kM5bMrd2rmJbdIl5p57S9IzQAAADTWeURkCJ0j+/Irqa/nehkxneRSdEr1Pt7cvN/XKXhed4/z1KdrUyoiInkNbtE73ngAgD5Za71rOyRFagYAALhJFhGZF4nwCIWWGBXpcg6E3NNJl3h+Z6UOTXvMGZNq+2WleF97aDLHS937JsV19zhPpV8bYDVSM0APcZ8Dw8FaMwAAAE7cpngfmxxzR/ArGfDDaJXhyD1q5hW3vHuy79pXXf+D3sVXBpeaAQCgTxg1AwAA4ITUDDpXYsi6bqIu73qWeM6Ascidmjl6+99J9l179Nbv9y68QmoGAABHY09NjP34AQCAo4URkZjwcNuJmZaZNChmkrB5Exy1rWsTTSZViqlHTDkp9tf2Gqx+bZeT2tXdR/O0nagr18q5KfbR5YRrpaxBE3vO6tJwy1yDquMdwlozpAnRNVIzAAA4YtQMAACAk+iISJMwXZNtqh6P3V9MSHWZtEqX65Wk0FU6oe11XLSfNvVoUpfY+7Itj9B0zL0d+x6qO2+p15pZ5jq1LTPF8c6T+nMnhWVSqm33l0KKz6lZdfcXaaZykJoBAMARa80AAAA4cZvQrC481mU6Ycy8z0/bNE3qcP28fefYn/c16CtC6khp0fsw94Rmr7nV7sm+a1967fd6F18hNQMAgKOxj5rpxRTvQ/0Fmarjb5f4FbqcLufi6aoeXZaTev4R7/dHbswRsq0ujiF3ROT1O94j2XftX19zXu9aNUREMFfq0T1j1tcPeQDdo7MqAACAk6VSMymmTY99XRN9Tyd0mbLpstNnTPmp5mvoKjzd5XwXQwiH51LiPB0YhpI6q/7Db6VLzfzVz0jNuCOdAADok7F3ViU1AwAA3PRi1AyGjVRFP6RYGbl03IvbGuo5KSk1c9yt75nsu/Yvrz63d+GVwaVmAADoE0bNAAAAOFmYmtlg62960jMk5xEaTD2qp+3IjJh9D2HU0ZgnthrCNPa5r0OTzumeE7ihn3KnZk64TbrUzBFXkZoBAAARGDUDAADgpKiIyLwwZ+oe+rEh8LYh17rwdeqQtneYuJQUWpfheo8UR9X2sWWnvh6lpCFizlmTfczTdn4i7wkEgT4oqiECAMDYMGoGAADAycKISInhw9Sh8S4nZopJvZR0rmOWb89Rduz2XU7zzwiMiRJTDk3SNE226apOAKqRmgEAwBGpGQAAACcLIyJdTsJVJfWS97HlxIhdJr6kdEadqlRUkzI8Jt6KmXiuxJWZ56UKSk/TeN/zVeXNk2J0Wu6UTi6l32cYB1IzAAA4GvuEZqy+m0mXvzz4VdNfXLvulTK9Pde6P3JP8f7udXsk+649/LJzeteqoY/IgPDhBgDoG1IzAAA4GvuomWSpmdhwZZcrvKZeOTemvBKjEiXN+dD2uudembbLzsi50nSzupwvZ0i838ve5VeJqVMX79O2S2Iscy5zp2b+5XbpUjNP/W9SM+6G1KMdAMaMz/NxIDUDAIAjRs30YNRM38OVfTWGYyxJ7vSTh7YjTmL3NzRdpZvHcv5i5U7NvO+390z2Xftnl27uXatmcKkZAADQH6RmAABwxKgZp9RMVz30CTWWYWyTO3mOTEoxLX7T17aROh3TJe86pb42KY4nZoTKMqPhYkbkdH1tcqdmPrhLutTMUy4hNQMAANAYqRkAABwxamaJ1MwyIeEcq/nm0qdj8J7cLEVKrpTzXWpYuU5JdalTN49E6fWfJ1fKpEt9H+m06N46MWzJ2jL48O3vnSw188SLz+5dq4bUDAAAI2dmO5rZ8WZ2kZltNbNvmNmBDV97iJl92cyumP4708ye3LRsGiIAADhaa+n+LWGjpEMlHSXpcZI2S9poZo9d9CIzO0zShyX9RNKfTf/9l6QPmtkzmhS8MDWzwdbf9GSJvf9TlDMrdY/+mLDymEO1KfbhcU09JpRiPZjf1GREzqy2613FTjXel5FiJdapBLlHzWzcNV1q5gk/jU/NTBsbn5R0cAhh4/Qxk/QFSetCCHPfAGa2SdJ6SbuFEG6cPrZG0g8lbQkh7FdXPhERAADG7QmSrpL0sZUHwiRK8R5J9zKzPRe89peSrllphExfe6OkayT9vEnhjJoBAMBRylEzZnZl3TYhhJ1WPXQfSZtnGxNTZ80+P2d3J0j6qJm9XNI/SjJJz5F0T0kvalLnhQ2RIYXqPEKQTULFVdvGKiVcn6L8tvuYd65TX3fvc1wK79FXKcrr8n6oek/mSjfHlMP9jKl1ks6rePzymecrhRA+Nu3U+q+SXj19+FpJTwohnNakcCIiAAAMREW0o/FL2zxnZo+S9D5J75f0EUlrNen0+n4ze2II4ZN1BdMQAQDA0Rr/Cc0uU3XUY+fp38srnlvp0PoeSZ8LIWyYeeo0M7uzpLdo0gl2ocajZqqUFNYrJT0xJH3oUd925EOpx7OiT3WN4X1cVSmlPtw7qT/f+LxcLPeomU/cca9ko2b++CffbjNq5p2SDtFkhMyNM48/S9I7JN07hLBNHxEz21XSRZKOCiH83arnXq9JH5HtQwjXLyqfUTMAAIzbRkk7STpg1eNPk3RuVSNk6gpJ10t6YMVzD5Z0WV0jREq4+q73L4YqqeqU+9eDx7nsskzvuTdi9lHKL8W282CU0ok0lRzvhRI/u7qQYkmCXPP8zNsmppxlVgrOHRH55F3umywi8rgLzmoTETFJn5V0X0kvlnS+pMM0aYg8PoRwynS7TZL2DSHYzGvfJOlISSdpMrHZ2unrnqyKSEkV+ogAPRQ7wRbQR0NuGJYkhBDM7CBJx07/7aTJcN2DVxohC7xI0nc1GbL7REk3ajIC56mS3tukfBoiAACMXAjhaklHTP/N22a/isd+pcn8If/Ytmy31XerpA4B5koxeKQeupyCPsd8B6mm1fZM9bStX19TAd4pwyp9On+zcqUhUhtCvatWfl8t9+q7n1p/v2SpmT/a8i33ITix6KwKAADc0BABAABu3PqI1IXHUkz1nZrH1PCp9Skt0KfznWs6/xSYDjydPr2fUoi5z5vsI5fSr42tHXdMYNxHDwAAXNEQAQAAbpJNaJba2EKeXSppkquhTsnepboe/2M+N0NQ4rUsZWRLF/socUKz/3uP+yf7rn3Med9k1AwAAEBTNEQAAICbxqvveoaEU4Xb2pZZV3bMPmb3U2JI1luKSfK6nASuS23Xl4ndd9t1RlKXk2IfsZPK9XUiwNRSfL6VoovrkTs185k9HpAsNfOoc/6d1AwAAEBTNEQAAICbxhOalb7aZ+qwYt3EPU1CvHV16nLkSOlpiHli69pk3QhPqVOGpV/L3MvLt63HbJnLrHuUOz3d5edcXw3hGGzNuGMC4z56AADgioYIAABwU+yEZkPWVZgVALC83KNmPnffByb7rn3kWV9j1AwAAEBTyVbfzdUJra1S6lFC+YhT0r0DAEOTrCECAADi2dreZVOSIjUDAADcdB4RKSWU7TFVNYZhzNe/pPdBjk7eqecLKXVJDM/9AauRmgEAwJGtHXdyYtxHDwAAXBU1j0hJIcAuw5sp9zu7b+9zNgbeYW+udTOpVzKOXeUX/ZZ7HpFNv/eQZN+1+33jzN71fCU1AwCAozWMmgEAAPBRVGomhZLSO6kN6diGdCyS7+qsbdM7q3VV19InO5zlXacSR+002UfblGFM/XJem9ypmTMe9NBk37UP/+qXexdeITUDAIAjW9O7tkNSpGYAAICbTlIz3iMLUpfT95EK3uHmFHLdU6Vc6yFcsxQ4DxOl3JdD0ORc5k7NfOEhf5AsNbPPmV/qXXiF1AwAAI7WMKEZAACAj6VSM7E9vKvCrH3oUV83GZlHKqptmX2dmKnLkR6pz2XVPjxH1fRZ3f1a6mdGnVyjWbqU4n2Tut5tR+GsdmLYkjW98aWH7ZMsNfMHX/wCqRkAANCcMaEZAACAD7cJzfrYC7yvoVAPfaprCkxo1qy80lMpJdWplPRNlxOapa5HKrlHzXx534cn+6596OfP6F14hdQMAACOxp6aGdwU795K+kUFpDSGlYJ5/25rqOdk0XHljoic+ch9k33XPuRzn+9dq4Y+IgAAwA2pGQAAHI19QrOFqZkNtv6mJz1Dck06u81KPbdE2/3W1S9VyDN1CHyoIfXSQ8xN7udSjqGUDqi55sUp5XiRR+7UzFcf9YhkqZkHfeZ0UjPe+DDAGHCfAxgKUjMAADhi1AyjZkaFEDNQFt575cmdmvnaHz4y2XftA0/7XO9aNYNLzQAAgP4gNQMAgKM1a3oXxEjKrSESs6LtvG3q9pt6FMm8enisvptarqmZU5RXSii7xNVHhyZm9d2+KvEYukzhxoyCHMqU8ViM1AwAAHBDagYAAEc28gnNkjVEYsNcVZN61W3rJaauucScb+96N0m51Sk9jJprcr0Sjx3xSllZdx6P+yzF5wT6adzNMAAA4IrUDAAAjtYwoVn5E5ql6GUdGzrPMdKjy9DqMvtmrRkfMWvNeNe/y1RUinrM8jhXua9T6utRyvX1kntCs/846NHJvmv3PvnTvWvVkJoBAABuSM0AAOBo7GvNLGyI5E4nxIrpZR1b19QjPbpUFQbu62iN2DKrjn3ePvoaok9xH6W+ljkms1qmnNjyU+8v973WNv06bx/en+0xSvpOQjukZgAAgJuiOqum6GDZ5LUeHRu7nL4499TIJcr1q8i7HH79zRfbcTX1XDxcj+HI3Vn1W0/+o2Tftff70Kd6l+cZXESEDwOMAfd5OpxLwNfgGiIAAKA/BjNqJuZXzdB+AZU4BX1ufV2luG1nxyN32KPI+9gzdRTbwbzLTrue84iUeF9gsbFPaEZEBOghvmwADAUNEQAA4KbxPCIruvwllmKuhVRTvMfsu61U4dSYuTRi69V2H22lmEekT2JHbZWyUmspU4CXkoJpUo73iCuUy9b4p2bMbEdJx0p6kqSdJJ0t6ZUhhI83eK1Jerak/yVpD0m/kPRdSX8VQvhy3esH00cEAAC0tlHS/SW9WNL5kg6XtNHMDgghnFrz2ndKOkTS6yR9WdKtJD1g+rcWDREAABytWevbS8LMHitpf0kHhxA2Th87XdJukt4oaW5DxMwO0aTR8rAQwpkzT32ycfklTWjmbaghTe/j8i4/Rimpnj6dM6AUqd43uSc02/y0A5J91+75z6dE193M3qFJSmbnEMKNM48/W9LbJd07hLB5zms3SQohhEe0qzGdVQEAGLv7SNo82wiZOmvm+W2Y2c0lPVjSt83sWDO72MxuMLOzzeywpoWTmgEAwFHK1XfN7Mq6bUIIO616aJ2k8yo2vXzm+SrrJG0n6TBJF0o6QtKVkp4p6d1mdosQwjtq65wqNdM2JJZ6TY0uQ9q51v/wOIZcSkl91Gk7eqfJ9oyeaKcv984y+n7NUo+mWuYzd5m65E7NnPP0A9OlZt59ylV126xuiJjZeZLODSEcsOrx3TVpoDw3hHDi6v2Y2R0l/Zcmo2TuEUL40fRxk/RVSXcIIdylrj5ERAAAGIiKaEcTl6k66rHz9O/lFc9J0hWSgqTvrjRCpnUIZnaapFeY2S4hhEsWFU5DBAAAR+Y8akaTOUMOMbM1q/qJ7DX9+52qF4UQtprZ9+fscyWqtLrfybYbMmrm1/oeFp3H+7i8y69SYp1mlV4/YMhyp2a++6yDkn3X3uudJ7cZNfM4SZ+QdFAI4WMzj58haZcQwr0WvPY1kv5Kk9TMluljJunrkm4bQrh7XflERAAAGLdTJZ0u6SQzW6fJhGaHSXqYpMevbDQdqrtvCGG2sfN6SYdKOs3MjtGvO6s+QNKfNCmchggAAI5sjW9qZtqn4yBNpng/VpMp3jdrMsHZKTWvvczM9tGkQfJ/JG0v6duSnhBCOLlJ+e6jZoagT6NmsK3c55vri6Eaymdh7tTMuc85OFlq5p5v/6j/wjWR3HvIAACA8eo8IpJr7o06/ApNh3OJIRvDfCVYLHdE5HvPfWKyiMjub/swEZFYvNkBlGLeJFgAuuPeEAEAAOO1MDWzwdbf9OSQ0iex0wOX3umqlFByKSmbUurRRJdTyY9BbASjKj0ce96b7Bv9ljs18/0jnpwsNfM7J3yI1AyA7vGlB2AoaIgAAAA3S42aSRFW7tPqu7nK7DJ1lHp1zFht00je9a5T4uq7Q1CXEvEcibeMFHX1SDGn2Hcp12nRvXVi2JI1vfGDv/yTZKmZux/3AVIzAAAATdEQAQAAbpZaayY2rDYvjJpi30OSK8zqOTIp9lhi7p1YqUPMpU/WlzvNtcx5qrt36raNlevaNXmP123vnXpuOwKp7Xuly2vj/X3jvdaMt3EfPQAAcEVDBAAAuFkqNZOCd0gsRpMQZZ949F6vGi3lre2xtw0x9/V+KcWQz1+XKcgSeV7Lku4jW7vWuwquiIgAAAA3NEQAAICbpSY0y6WUCXBmlVin1Pp0jH2q6zxDOIY6uSYw9D6Xpaz/VMf7PKXQxdbLnbUAABpPSURBVDHkXmtmy0sOS/Zdu/6172FCMwAAgKZoiAAAADe9SM1gW30NqfYlZA10KddaM2gnd2rmxy97erLv2rv+/btIzQAAADRFQwQAALhJNqFZivUuUvein7fvtnL10G+7BkeqfZco93WNXWI9dZ3aKrFOy8ixJpXnxH7L1KXEa9rk/kuxvkyK9+eie+vEsGVhvVKzteOOCYz76AEAgKuiOqvm+mXi8QuoqwhQk33nUmK0xSNaNTYlnpOYX+bLrFg8b991+yjlPKFa7s6qFxz1zGTftXd59Um966zqvtYMAABjRmoGAADAycLUzAZbf9OTMZ34UovpmJSqTilSKU06WrXdR5c8Uz2pOz17a3tf1t07qfaXWu5O5ilSJk2UMpU8qcZtDWGK9wv/9tnJUjN3PuYdpGYAAEBztmbcyYlxHz0AAHDVOCLSJARatW2KlTRziTnGFKHxmPJilTSfRCnh69QjKZro0/2foh6ljxxJkcrr8n1bV2aXqecur80yqezU5aA8pGYAAHDEqBkAAAAnySY08w7rdanvK8bmCovm0nZK9iEZ2jWdJ8V04JioOz+pp2ePlSs10+TzPPeomZ+8+rnJRs3c8ai3MWoGAAA0R2oGAADASe/Wmkkx0VPbSbOWCUV2uYpo6km9PCfNSjGh2TylTwBV0kintnJPuBU7giXFfZlifyhb7tTMRX//vGTftXd42VtJzQAAgObWkJoBAADwUVRqplSMmoGnMV+bMR97ajGjzTzW1hnzqJmLX/f8ZN+1t3/xW0jNAACA5lhrBgAAwIlbaqardMeQQ7l9OraY0Q6lH0sTfRpVkbuuHmH+tkqqUyl1SZXuaLvvNtsuK3dq5pI3/GWy79pdXnQcqRkAANDc2Cc0G1xEpKqMNuV4/hrJPRcD8urT9Si9c2IqJUZoSzo/njwGC+SOiFz6phck+6797Re8qXcRkXE3wwCMXuykaMiHazMOpGYAAHA09tRMsoZIbBgxR5htmTLGHAodm9JD4J7167K8LpZj6KoeufZd+r1YJ3VH6Ddfd07lPpkbaVjG3QwDABSL1Mw4kJoBAMDR2Cc0WzhqZoOtv+nJmNDW0EZ9VPXazhUa9Dj21Kv55jCEUG0pc5F4jJRpW2aTX8ylz12C8uQeNXPZCX+dbNTMuiNez6gZAACApkjNAADgaM3atd5VcDXoCc2GJnXo3nOyqqGFyD3TKqVPmsUU783LibmPhpAeZor3iSve9tJk37W3fe5rSM0AAAA0RWoGAABHTGjWgSGt0zK7z1LCx1KaupR0PMinxJTI2Azt/dv2M7KktX3gZ9zNMAAA4IrUDAAAjsaemmHUTIEILwIoxRg/j3KPmrnqpKOSfdfe5pmvZtQMAABAU6RmAABwxFozkWvNeK99MpbQoKc+rjUzNENbrylHPYaw3tA8qVPZpMYXy52aufpd/ztZaubWT38lqRkAANAvZrajmR1vZheZ2VYz+4aZHRi5DzOzz5lZMLM3N36dV2fV3GJ/CfX1l1NunKeycD0mUq/mW+L06LkM4Rhi5Y6I/Oyfj072XftbTzu6Vd3N7DOS7i/pxZLOl3S4pEMlHRBCOLXhPp4j6WhJd5B0XAjhyCavo48IAAAjZmaPlbS/pINDCBunj50uaTdJb5RU2xAxsztJep2kZ0r6cEz5pGYAABi3J0i6StLHVh4Ik3TJeyTdy8z2bLCPt0k6I4TwkdjCl0rNjDFkV6f0c1Ji/UqsE5ZHh0iUouozZtHnTu7UzDXvfWWy1MyOh/7v6Lqb2ZmatD0euurxB0n6iqSnhBA+tOD1fyrp7ZL2DCFcYGZBpGYAABgfM7uybpsQwk6rHlon6byKTS+feX5eebeTdJykl4cQLmhaz1mkZgAAwKKozKLnjtekc+sJbQteKiJSUsi1Lgyca9TM7Lbz9tFVyLpJnUu6ZiuanLN5Ys4lKaBupB5pkkKuOuV+jy9TpzGrOg8lnZuUE5pVRDuauEzVUY+dp38vr3hOZvYoSU+R9EhJtzb7jazQdma2k6RrQgg3LCqciAgAAON2tqQ9zGx1m2Cv6d/vzHndvTVpR2ySdMXMP0naMP3v/esKp48IAADjtlGTYbcHaGbkjKSnSTo3hLB5zus+LOk/Kx4/XdJHNEnXnFVX+GgaIrFhuBRhu5JCf30xtnPW1zB6n+pdl65dJq1X4rGXWCcsZmvWelfhVE0aDyeZ2TpN+nwcJulhkh6/spGZbZK0bwjBJCmEcKGkC1fvbJqiuTCEsKlJ4aNpiAAAgG2FEIKZHSTp2Om/nSRt1mSCs1O6Lp+GCAAAIxdCuFrSEdN/87bZr+G+ouYyWdgQKSX8Wtczvesy68puMlKnaj+pz+8QVuaMPSd125RyD8/T9p4q8Vg8LPMZ0HY0XBfbY+T8UzOuGDUDAADc0BABAABullprBtsqPRVQIs5ZP4zhOo3hGGMN9ZyUtNbMdR95Q7Lv2h0OeVHWuqdARAQAALihIQIAANy4Dd/1XG+lS12W6X1sMerq2qdjacJzvZVl1lGaVffaob1XS5fi/OQ6x20/z2PqV+oowxRsLaNmAAAAXCzsrLrB1i/sQJN7JcsuymyiqrWf+5dGbDmp6teneStKrKv39WuzP4/3W4rzNKvt3D5NyunT50Cu/Q1N7s6qW09+U7LOqtsf9ILedVZlZlXM1eWEcamVXlc+7MvVp2uTuq59OvZBY0IzAAAAH27ziHQZRi8xRJ9Cn8KpdFadyN1ZtUnZuevap2tdUl1L6aza5dIDpaYMs6dmPn58utTMgX9BasZb6SF6AAB+A6kZAAAAH72Y4r2kcOmKvvaW70vZsUqta0zIusuUSSnnh9Ef7Qx13qXUUh1P9tTMJ96aLjXzx88jNQMAAJqzNeNOToz76AEAgKuFEZGqMJdHKM87ZFgVFh3ysQ911JGHUs7hUOsRO7nY0FIR4JoOAakZAAA8MWoGAADAx1KjZpaZSMkz1ROrLlVR+tocTeZWyTXyp215nhOGzZNrHZm661fSfdmlFPdRKUqZrGwZXd7/Oa71ovfViWFL1pEn13/6pGSjZm756Gf2btQMEREAAOCGhggAAHBTbGfVUpZPn92Pdyi0rXkjC1Jock66LD+F0q9r2zqVeCyx6tK5Q9D2/isxJRf7eTDPEO7dGMwjAgAA4ISGCAAAcJMsNTOGUNoYjjHWEM4Ja570w1An2ks9asbz/MSOAhvatWyNeUQAAAB80BABAABulprQLJeYCXBShfrqJtRJsa5FqrUxch97l1JM9tXlaIK29Vtm4rKq+26ZCeFyTV4Xc8+nvjazPI+39PKGkCbp4hhuuf32WScF+/mm9yb7rt1uv0OZ0AwAAKCpXqy+G+PIHfbIUq9cv0iHNuY+R4fDUn79pyq/6pd+l/e59wrbbZcvaDJfTZeRlxy6jHoCXoqd0KytPn0pA21xnwPDYWsZNQMAAOAiOiLiHepOEVL0Ti91Oc36ihTH1UX6yXPK6dRSdFaNXRk5R0i9SZqkbT26TDXGdoRN0ZG4bvtSO3kTUUNJBpeaAQCgV1hrBgAAwMfCiEgp4bsU8w0sE9JOsfppinlJcil9iui67T3mk5gn9QipHPWIfV2J90tJ76cVJY3mGpIxHOPQkZoBAMATa80AAAD4WDjF+wZbv3Da2VJ7hHdV/jJl554cbpl9e07x3pb3PZKa51Te3vdlivRrSe+nEvXxPZ5T7inef3HmR5JN8X6LhxzSuyneSc0AAODISM0AAAD4iB410+WESvPCn30KD8eUn3rdiFT7q0odxdQjZV2ayjUCqa+rn3Y5eif1MXQ1Sm2Zfc/yXOU39XscKAGpGQAAPDGhGQAAgI+Fo2au37o1WU9eD0Pr6d4W56EsfUodATH6mrpcLfeomV9+/ePJvmtv/vsHMmoGAAA0x6gZAAAAJ26pma6WyU4VyvOcVCl2OXNPbeta0nowdWJHWc3KcTyx90VMXbucGK/E+3mWRwot9X3U5edR1TZdjrjKmfbJnZq54d8/mey79mYPeBypGQAAEGHkqZlBd1adtcyvxhJ/rZWC8zRsY7u+YzveJsZ4TrJHRP7jtHQRkb3/sHcREfqIAAAAN6RmAADwNPIJzRqnZnKs1pgrBOgxxXuX569uKvZcHd1ydYyL2cesEuc1aDKNfukdPEtZ2XdeOXVSdMJts582ZZZyzYcue2rmW59Ol5q536NJzQDoHl9IAIaC1AwAAI5s7bhHzTRuiJT+CywmdFnKsXQZ4k01h0pupVwbD03mkyhRlyv7xlgmzRWz75h9pHqPj/l9geEjNQMAANyQmgEAwNPIJzRb2BAppYe+dz1STF9cxyMdU3dePc91bPlNQuDe91Gd2NSCJ88lEObtL5eYYyjlegElIzUDAADckJoBAMDTyFMzo1lrZp5SwvUpVjklDAwgRimff6XJPaHZr777hWTftWvvtU+rupvZjpKOlfQkSTtJOlvSK0MIH6953bMkHSjpfpJ2kXShpE9JelUI4dImZZOaAQAAGyUdKukoSY+TtFnSRjN7bM3rjpF0taSXSfpDSf8g6cmSvm5mOzUpmNQMAACOzHmtmWljY39JB4cQNk4fO13SbpLeKOnUBS/fO4Rwycz/f97MNkvaJOmpkt5SW/5QUjOkJ8YZZq277h7rglSVP7Tr4X1e65Rev9TG+N7vUu7UzI3nfSnZd+2ae/xBdN3N7B2apGR2DiHcOPP4syW9XdK9QwibI/a3g6RrJb0uhPCS2jrHVhgAAAzKfSRtnm2ETJ0183yMR07/fqfJxqRmAADwlHDUjJldWbdNCGF13411ks6r2PTymeeblr+zpOMlfU/Shxq9JlVqpvTQoMeS9l3urxRDPa6+4r7tHudk+LKnZr7/lXSjZnZ/yFV126xuiJjZeZLODSEcsOrx3TVpoDw3hHBi3X6nKZnTJN1X0sNDCGfVvEQSEREAAAajItrRxGWqjnrsPP17ecVzv8HMtpf0cUl7S3pM00aIlLAhUvovg9JXu+zrPCKlX/dZJf6STV2n1MdVynlaRonXHfgN5t5d82xJh5jZmlX9RPaa/l3Y18PMbinpY5IeIumxIYQvxxTufvQAAMDVRk0mMTtg1eNP0yRlM3fEjJltJ+lkSftIenwI4fOxhZOaAQBg3E6VdLqkk8xsnaTzJR0m6WGSHr+ykZltkrRvCGG2D82HJT1G0islXWNmD5557tIQwg/qCl+qs2qKDqCETdvpcpXaLqVII5Vyz8TUo5Q6S2XVpU7d6rqp69/Xc5P7/bua5wrMbS2qR/bOqj/8Rrp5RHb7vbZTvN9akynen6hJdGSzJlO8nzyzzSataoiY2aK6vyeEcHhd2UREAAAYuRDC1ZKOmP6bt81+FY8t3WijjwgAAHDjPsV7KWG6sagLdUvlXAfuDcwa6nT5fTXk92f21XfP/2a6eUTudv+sdU+BiAgAAHBDQwQAALgZzRTvfTLvXPbpHI+triWOoGlyH80q/Trl0Kf7toncx5NriYGujyt7ambLf6ZLzaz/XVIzAAAATdEQAQAAbtzWmumqB/wQQqt9rXeMkq5T7nuwy/K6LKetIaQax6zJdWr7eR7zXhn0/WK9y6YkRUQEAAC4oSECAADcLBw1s8HWb/Nk2zVlmiCEu63U5zJH2V3tZ0w8z1mXZZc0kicmnTC0e7jqfPf1uLq4NtlHzVzw7XSjZu6yV+/yPEREAACAGxoiAADAjftaM33lPSkVhodrHY9zthzOX7XcqZkbLjw72Xftze58b1IzAAAATS0VEYltTdd1VKNF3lyKTqyzcp37FNe6lF9xbe9/7/u8lPPXxFA/M/p0Debx7Eif2urPxRPDFiIiGSWb0KytJsvSA33HfZ4O5xKDY+NOToz76AEAgKvBdFYtMdyX2xjPQV2I3jsVBR/e173vqaOxy95Z9b/OSZeaudMepGa88cYHAPQKqRkAAAAfySIipacFUo3wKZlHPUs/T02mEUc8j+sek+7wvu4lvheAUg0uNQMAQK+QmgEAAPAxmFEzY1B6GmRWXV37dCzzNAn5l3hsuUeU9GlVbe86dXltYo4tth5VabPYlZbrymtS/1TXL/uomYu+l27UzB12Z9QMAABoLpCaAQAA8EFqZgS8w81YTtX145qiiRTpmHna3ncxKZsuUjNNRl/lTs388qc/SPZde/Nd705qBgAARCA1AwAA4GNhamaDrV8YLhpLSLiP60YsE7qvC5GWaGipCs/j6bJs73339T7p42dQLl1c0+ypmYvPT5eauf3dSM0AAIAI1ru2Q1KkZgAAgJuFEZG6iajGLEU4sMte4ClClISBlxcTUo95b3WZYujyuje5z2PqUtLnUY60T19TS13iPPQfqRkAADyNfNQM84g0wK8Q5JZ7GnYMX+o5RWKmc19m6ve6/cVqch6yd1a99MfpOqv+9l171+Fk3M2wSHwJAMB8bdeUybW/efuGL1IzAAA4GvtaMzRERmAI8yvEhH6HoMkvwaEeO9JJ0QHas0PwkTvsEXV/l/45hmrjboYBPVXSaBGgKzQmxoGICAAAntaMOybg1hDpKqycKjTXdj9dzi/Sdt9Nts1d19hRIaX/MupylEvVPpa5L+qiKSnuuXnlpNp3V0oarVT6NP8ponIekT1WsC7PuJthAADAFakZAAA8jXzUDBOajUBJ4ea2hnAMbVWFkAkr95vH9UudLm2beosZnRObymu7JMZquSc0+8UVP032XXuL2+7KhGYAAABNkZoBAMATqZn2qZlleusTYl5OilE9s3KHh4dwrbl3uxczwqdPUrwPc91/qT8zUqVP2u67atvVTgxb8qZmrrwkXWpmp11IzQAAADRFagYAAE+kZoY1ambI4fI+HVtMD/3Sj6UJz5RXqknHuqpr6gn6ulRSnbqcHLFu2yZi9td27ZrRjJq56r/TpWZucztSMwAAAE2RmgEAwFEgNZMmNVNSSHNFiXVqoq/17qvco7k8RjsM4T4qfX2qFEqqy7I81t1KlSLNnZr5+dWXJ0vNbHfrnXuXmiEiAgCAp5FHRMZ99AAAwNWgIyJ9D23CD/dOebq6JiVda8+00zKjZurWiWlSp7ryY+tXt5+SrvvYDbohAgBA8ax33TqS6mQekSF1uBoCrkf/cM2AeKneN9k7q15zVbrOqjvepnetGvqIAAAAN6RmAADwNPJRM500REpcQXJoUqy+63G+U6yqOcvznunyGvBemBjqit25pm+Ped+k6Ay6zBT+dR1eU133Jsd5YtjSev+IN+5mGAAAcEVDBAAAR8HWJPvXlpntaGbHm9lFZrbVzL5hZgc2fO3dzexkM7vKzH5mZqea2Z6Ny2b13f7o07H1qa4peB7vMiFwz9WBS7xHSqpT7tV3l9lHjrk5cl6b3KNmrr/u2mTftbfc4Vat6m5mn5F0f0kvlnS+pMMlHSrpgBDCqQtet4uk/5R0iaSjJd0g6ShJd5e0dwjhwrqy6awKAMCImdljJe0v6eAQwsbpY6dL2k3SGyXNbYhIepGk20r6vRDCT6avPVOTxszLJT23rnxSMwAAeLI16f618wRJV0n62MoDYZIueY+ke9WkWZ4g6TMrjZDpay+TdIqkg5sUPrjUTBf6PiVwbO/6Eo+z9PqlVjfaYWznA+10OWqmKs0Wu0JuTKon1Yq7TT7Ps6dmEn7Xtqn7NIIRQggPXfX4gyR9RdJTQggfqnjd9pKulfSaEMLfrHruJZJeI+n2IYRLFpVPagYAgIEwsyvrtgkh7LTqoXWSzqvY9PKZ56vcVpLNbDfvtTREAAAoVeIITG1DZI5FUZm6iM0yr/VviPQhxFxqvVLq0zH24Z5ZVt1xDfW4VxvDte6rVKvhVu0v5lrHrODbpGyp3xOaVUQ7mrhM1VGPnad/qyIeknSFJg2NNq+9CZ1VAQAYt7Ml7WG2TW/XvaZ/v1P1ohDCVkk/lHSfiqf3knRpXf8QiYYIAABjt1HSTpIOWPX40ySdG0LYXPPaR5nZrisPmNnO0319tEnhyUbNlBhCLbFOTfS13n2Vez2TXPue1df7KMd7oaT3W0l1aSPX/df1SMDco2a8mZlJ+qyk++rXE5odpklD5PEhhFOm222StG8IwWZee3tJ35L0E0nH6NcTmt1DkwnNflxXPhERAABGbDpnyEGSPiDpWEmf0qRRcvBKI2TBay+WtI+kCyT9i6QPatJh9uFNGiFSAZ1VAQCArxDC1ZKOmP6bt81+cx7/nqTHty3bbUKzriYJ63toc5E+HVtdXft0LE14pkRSLLfe9LVtsNZMO6nXiZlVyoRmMWU0kWpdnLGlZryRmgEAAG4G3Vm1rziXeZWyGm1fy0E6fb9mXUTbYqLnqconIpIXEREAQJFST5qGMtEQAQAAbhamZjbY+pue7GOYMJWYjlhD0MfVhku9Hh5h5b6LuZZNfjEP4fylfk+W+B4v6T1MaiYvIiIAAMANDREAAOCmk1Ezs1KMf/cO05UYxowROx1y3bYp6uJ9Lku6v6rUvZ9S1b+U81BKPYam7edyrk6idfOLNJl/JHaK97pUuySdGLaQmsmIiAgAAHBDQwQAALhhQrMCcS7zYkIzII1cKfqulw1g1ExeREQAAIAbGiIAAMDNUqkZQr/9wzXb1lBTM6n1td5ALFIzeRERAQAAbmiIAAAAN8lGzQxZKZNwtRU7oVnq4ySkv5yq88c5Ray6e6bJiJe2k4TFPh7zunnHwKiZ/iAiAgAA3NAQAQAAbtxHzaQOMXukGDyWb297nN5LzceEcKtet1op66PMilnHo09pldh6p0hp1q150tc0Yp/ugdSTh5V47Kw144uICAAAcHOznIVVtYTffN05SVd67LKFPW/fpbTq+6rJ+Sv9HJdevxRij7GP5yRXnUuP/s7yuI45oi0lRmbGyj0ikmu5aQBAv9BAGAf3hggAABgvGiIAAMDNwlEzAAAAXSIiAgAA3NAQAQAAbmiIAAAANzREAACAGxoiAADADQ0RAADg5v8DI1ihSS2FuZkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAHrCAYAAADscsZ4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debwkZXno8d8z48ZihBkENCQZFg2K3LjGEDVA1BCN4ChGgxsYk7jvCnEhDqAYFW/QxIsXQ4JblGgyKIrixrhrNG4XUFA2QURkgAnggAvP/aPrDD1Nd586faq73u7+fT+f/vSZ6tq6urrrmfd53rciM5EkSSrZirZ3QJIkaTEGLJIkqXgGLJIkqXgGLJIkqXgGLJIkqXgGLJIkqXgGLJo6EXHfiPhMRFwbERkR62osc2pEtN6Hv5T9GJeIuCQiNixh/iV/ln3WsSEiLumZNtPHWZpHBixqVEQcUF14Xj5knoyIj3b9e2VEPC0ivhgRV0bETRFxeUScHRHHRsQdu+a9HfAfwD2Ao4GnAf85xrekMfGzlLQUt2t7ByTg34AnAl8C3gJcC/w28PvAq4C3ATdX8+5RPV6Wmf+0hG38NfDspnZYjRj1s5Q0hwxY1KqIeACdYGV9Zj6+z+u7AJu6Ju1aPV9TY90BbJeZN2TmL4FfNrDLUysi7pyZ17e9H11qf5aSZEpIbbtH9fzZfi9m5k+rYIOqNuJz1Uv/WqWWMiLWdKWijoiI50XEecBNwMurZfvWNETE3SLipIj4UUT8IiKuiIiTI2LnnvnWVev/3Yg4vkpZ3RwR34mIR/dZ750i4s3V+jZHxH9FxJ/0e48RsU9EfDAiflyt88oqHfZnPfOtiYj/iIj/iYhNEfHhiNi9X91Ita+nRsTDq1TbDcAZ1Wt3j4i3RMS3q9qRmyLivIg4KiJW9qzniGpdj6iOwaXVPn43Iv6i3/uplts7Ij4WEddX+/qhiNi16/UNDP4sF7Z5QJ/13qZeRdJ8sIVF47JtROxUY74Lq+c/j4j3Zea1Q+Z9PZ200auAk4EvVNN/Bqyp/n4xsBp4J3AlcNmglUXEbwNfAe4AnFLty17Ac4ADI+KBmbmpZ7F30WmpOaFa7sXA6RFxz8y8pGu+9wNr6QQJZwF70qnPuLhnH1Zza7D2DuBSYCfggcCDgY91zfcFYJdqvu8BDwPOBrYb8BYfCBxaHYt3dU3/X8DjgfXVe7498Cjg7+mkaJ7VZ11vrLZzEpDAM4D3R8SdMvPUnnl/E9hQrf8VwO9V6/wNYCFoG/ZZStJtZaYPH409gAPoXNAWe3y0a5mPVNNuBD4FvA44GNh2yPqPGDD9GmDnPsud2jndt5r2YeAqYLee6Q8EfgWs65q2bmG/geia/qBq+hu6pv1JNe3UnvWuXXj/XdMOqaY9cZHj+qZqvqcMmL6hZ/rCcX5En3Vt0/0euqa/B/g1cLeuaUdU67kUuEvX9LtU064Btumafkm/9wO8vZq+d43PcmGbB/TZxw3AJTWm3ebz9uHDx3Q/TAlpXE4GHjng0etQ4AXAOXQuYq+mE8RcGREvW+J2352ZVy02U0TcBXhMtZ2bImKnhQedi+4PubU1oNtbM3NLaikzvw5cz62pLegEJgBv7l4wM08Hzu9Z30ILzqMi4jeG7PLBwE/otNx0O2HIMt/JzE/3TszMzQvvISLuEBGrqvd9Fp008QP7rOuk7Gptqv5+B7Ajnc+s2xWZ+e890xZakfYasr+SNJApIY3LD/pdLAE6tbC3yk6Nyj8B/xQR2wAPAB5NJ4g5ISKuyMzeC/UgF9Sc73fpXJyfWT36uajmtGvopKEW7AHcMmBfvldtG4DM/FxEvJtOq8JTIuLrwKeB0zLzvK7ldgf+KzNv6V5ZZl4VEdcN2P++x6LqTvy3wNPpBBDRM8uOA/a718L+7dEzvd8x2lg9r+7zmiQtyoBFRcnMzcAXgS9GxNnAJ+kEFHUDlp/XnG/hIv1etq7v6La5z7RfL7K+3r+HzQdAZh4eEW+mE6Q9FHgZ8OqIeHEur7vvoGPxv+kEg6fRqSW5ik5dzv3p1Kr0a3ntNwjboPc56BgNW2axbS3wN0uaU375VbKvVs+/OYZ1/5DOhfEOg1qCluFCOumkewLn9ry2d78FMvMcOimxN0XEDsDXgL+PiLdX6ZtLgL0iYkV3K0vVm2mHJe7f04DPZ+ZWvXwiYli65t500mfd7lU992tRWY6Fbs6r+ry2O3PePV2aV9awqFURcY8hF8qFWpDzBrw+sszcCJwJPD4i/qDPfkVE3HXE1X+4en5FzzrX0pUOqqatioitvoeZeR2d3kTbAneqJp8B3A04rGdbA0cUHuLX9LR0RMR2wEuGLPOcqu5nYf670BmI7zpu7Z7clIVU1iN69vEw4O4Nb0vSlLCFRW37PeC0iPgcnd4el9PpPvtgOgPKXQ8cO6ZtP4dO+unzVR3Jt+gE8XsAjwXeTad30JJk5lkRcQZweESsAj5Bp1vzs+i0otyna/anAy+JiPV0Wn1+CewPHAT8e5Uig06q5sl0xiz5feD7dNJHDwGuZngapdeHgGdFxGl06mV2Af6SW+tM+rka+FpE/AudYOcZdEYj/qvMrJuGqyUzz4+IT1f7GMC3gfsCj6NzjG7f5PYkTQcDFrXt83RaIh5J56K5C50L4mXAvwJvzswfjmPDmXlZdEbaPYpOgPJUOoPNXUanRaO3p8tSPIlO9+yn0Hlv59DpDXUYWwcsG4D70emxdDc6rR8X02k52VK/kplXR8RD6dy64C/pBChnAwcCX6d/vc0gL6UTCD6Rzvu+jE6vroWC336OojPuy/PpfEY/oNPF+t+WsN2leBrwj3SO39PojNNyIJ1xYNaMaZuSChZdPTQlTZlqQLmrgf+bmY3fKykijqATOB6YmRuaXr8k1WUNizQlqi7fvY6qnj81yX2RpEkzJSRNj49HxKXAN4CVwMPppJK+DJze5o5J0rgZsEjT4ww6Rbpr6QyvfzmdmpZjMnPY2CeSNPWsYZEkScWzhkWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBXPgEWSJBVv6M0Pb9q8ecuNhja+6cVbpq8+8sRlb/jF295ry98n/vx7y17fOHXvK5S/v01/VsvVvT9Q3j6VsD+j6j03u5V+npag7u9Q3fNlmn7XNNidttkm2t4H3ZYtLJIkqXgGLJIkqXiRmQNf/PExz9ryYhPN5rPSDK+lKT2lVvr+aXya+E2axTTQLL6npTAlVCZbWCRJUvEMWCRJUvEMWCRJUvGG1rB0d2uua1bqVGblfQwz73lqaRTWPM0+a1jKZAuLJEkqngGLJEkqXuMpoRJMW5PtNKVmmk6Vlb6+eeAxk7ZmSqhMtrBIkqTiGbBIkqTijZQSGvVmdoNSH9OWwpkV40wFTFOaS0s3z5/vNP9elfa5lXhjVDAlVCpbWCRJUvEMWCRJUvEMWCRJUvFmsltzW6Y5t61blZbn1+yxK3nZrGEpky0skiSpeAYskiSpeLVTQnWbMKe5OX2a971kpXZdLJnn4mwb5fOd5DlRd1uzOlSFKaEy2cIiSZKKZ8AiSZKKNzQl9OxYs+XFaWvSm2e9zbEL/Axn27Q3w4/LsOOy3NTHqPOpbKaEymQLiyRJKp4BiyRJKp4BiyRJKl7jI92WnsMtIRddt9Zg2roDWzsjla303+dSWMNSJltYJEmaYxHx8Ig4NSLOj4ifR8TlEfGfEbFvzeX3jIjTI2JTRFwfEWdGxL0HzPvCiLggIm6OiAsj4siIqBWLGLBIkjTfng38NvAPwKOAl1b//npE/MGwBSNiZ+ALwBrgcOAwYBXwuYjYrWfe11Tb+ABwEHAK8Hrg+Do76c0PNXXsvls20w6advOWEoqInTPzqp5pOwAXA5/NzEOHLPsm4AXAnpl5RTVtdbXs+zLzOV3TLgdOzswXdS3/euBIYPfMvHzYftrCIknSHOsNVqpp1wE/AHa77RJbeRzwqYVgpVp2I3AG8Piu+f4UuBPwrp7lTwVuBxyy2H4asEiSNCMi4rrFHjXXc1fgPsA5Q+bZBthzwDzfBXauUkZU60rg3O6ZMvMHwObq9aFuV2fHpZIMSzPUvUmnFjdqLzXTQEszKylOU4Gj6x5VvgGblruCiAjgZDqNGicMmXVHIIBr+ry2MG01cFX1/PPMvLnPvNdWrw9lwCJJhfBCr+XKzB0aWM2bgbXAMzKzzkk5LODKAX8vZR2AKSFJklSpimBfBrwoM09dZPZr6QQa/VpHVlXPCy0tG4HtIuKOfebdkf6tNFsxYJEkSUTEscCrgCMz822LzZ+Zm4GL6F9/si/ws66C3nPppI/26dnmXsA2DKmV2TLvNHdrtl5hOk1rrnvctQbzdj43fR40ffym9TzV8k26W/PzVjRXw/L2Wy4Zad8j4rXAOuDozHzdEpZ7M/B8Ot2Sr6ymraLTrfn9mfnsatpq4MfASZn5kq7ljwP+FtgjMy8bti1rWCRJmmMR8TI6wcpHgU/3DBZ3c2Z+q5pvA7B/ZnYHRScATwPOjIhjgF8Br6metwwIl5kbI+INwNERsQk4G9gPOAo4cbFgBQxYJEmadwdXz4+pHt0upTOKbV+Z+dOIeBidwOU9dEpNvgD8UWb+qGf2Y+n0Ynoe8ErgCuC1wBvr7ORUBSzTdjNAzRbTAmXz90DTamW0O7BuZh6wnPmqsVQeW2P5BE6sHktm0a0kSSqeAYskSSreVKWEbPKdDUevO6jWfPPWa2belJ5iK33/NDtWztWtFkdnC4skSSqeAYskSSreVKWEJEmaNW33EpoWjYx064iQyzPPx2+e33sprBWStjbpkW5ffvs9Ghvp9oRfXjSz0Y8pIUmSVDxTQpIktcheQvVM9c0PpUka980PJdXTdCq597v9jhztBoKjOuoOzaWE3vgLU0KSJEmtMSUkSVKL7CVUTxEpoRKb2sfZe6XE99ut9P1r2rjfrz2htBTjPF/sEVbPpHsJveqOezZ2rT3+5gtnNvoxJSRJhegNniXdypSQJEktspdQPbawSJKk4hVRw6JmWTNRLj+byRh2nKf1M7D+ZHImXcNy9J2aq2E57qbZrWExJSRJUovsJVSPKSFJklS8iaaEprUpVpo0vytq07yff5NOCa3bZq/GrrXrNv9wZptrTAlJktQiUx31eJwkSVLxRkoJdVerw/Ir1ic5smoJTZ3zNpKsxseeI1LzJp0SOrbBlNDfmRKSJEnjYC+hekwJSZKk4hmwSJKk4g2tYXl2rNnyonUWmiUl1DI1YdjN8tp6X+Osq2n6c6u7vqbnU9kmXcPy99vdo7Ealr+98Qczm1+yhUWSJBXPgEWSJBXPmx9K6qvp4QukaTHplNCbt79nY9faV9xwgSkhSZKktjgOiyRJLVo5s20izTJgmQGOnNuOWe8R0lYKyFSUpH5MCUmSpOIZsEiS1KKVEY09RhURu0XEWyPiixFxQ0RkRBxQc9kc8vhE13xrhsz3p4ttx5SQJEnaCzgM+CbwGeCQJSy7X59pfwIcA5ze57UTgdN6pi2aVx9rt2ZrK2bPPNwdeNZrUyQNN+luzW/9jd9trFvzi/7n/JH2PSJWZOYt1d9rgfXAgZm5YcT1fQR4BHC3zNxUTVsDXAy8JDOXfAGxhUWSpBaV0EtoIVhpQkTsAjwK+MBCsNIEa1gkSVKTDqfTIPIvA15/dUT8IiJujIizI+LhdVY61hYWm9PHp620RXcaaFZTfrPyPkpgem36jfMznNXfkDZFxHWLzZOZO4x5N54BXARs6Jl+M/BO4JPAlcDuwMuAT0XEoZm5fthKTQlJktSi5fTuKU1E/CGwN3B09hTJZuZPgL/pmvTFiPgP4NvAm+nUzQxkwCJJ0oyYQOvJYv4SuAU4tc7MmfnziPgQ8MqIuGtm/mzQvAYsU6qEptO6+1D6yKU2SzfLNND0m9Rn6PkxWyJiO+CJwCcz8/IlLLpQTzu08NeARZKkFpXQS6ghfw7cmcHFtrcREdsChwI/zMyNw+Y1YJEkSUTEE6o/H1Q97x8ROwE3ZubHq3kuAcjMNX1W8ZfARuDDA9b/FjqtKV8GfgasAV4C7AGsXWz/DFgkSRLAB3v+va56vpROcDFQROwJPAx4W2b+YsBs5wLPotPt+c7AJuArwHMz80uL7dxYR7qta5I1BHVzs/Och5+3mo5xv995Ppe0NOM+F+dhpOomTHqk21NW7d3YtfaZ13x/dhJMPRw4TpIkFc+ARZIkFa+RlFBpTd7zltKYZqWdO/PINIG0tUmnhE5dfa/GUkJHbPyeKSFJkqS2GLBIkqTiFdFLSO0YlI4xpaalKH0k46b5/SjLONLKk04JvWen5lJCT7valJAkSVJrDFgkSVLxHOlWkqQWrYyZzeI0qogalnnLgc+qUXLJs9itubfG4eh1B235u61zu/c71s3vWzn8LbxVm93tJ13D8m93vXdj19on/+y8mY1+TAlJkqTimRKSJKlFK2e2TaRZRaSEND0cFVW97OarNs1Ct+bTdm4uJfSkq0wJSZIktcaUkCRJLbKXUD2mhMaotB4ww5rul9vDZynLaX6U9h2YNnWP3zhTteP8DEv9DZl0SuhDu+zT2LX2CT89d2ajH1NCkiSpeKaEJElqkb2E6rGFRZIkFc8almUyRz+//Oyl2TTpGpb1uzZXw/K4K2e3hsWUkCRJLbKXUD2mhCRJUvGGtrBMa1e5UfZh1P0oPRVQwnGeVaMcz3nsItqkUbr5HrfurIHz1b0xZd3ttvV9a3q7pY9oXfr+aTxMCUmS1KIVpoRqMSUkSZKK11ovIVMV4+OxnU/dzeRgU/mCJkZ4LmHEWU3OpHsJffTu+zZ2rX3MFf9vZptrTAlJktSicOS4WkwJSZKk4rXawmK6YnkWmqn7Hcfl9GBZyrKmn8pSYhpilPOqacO2XXe/6s5X9zNYSB+V+Jn1GvYZLvd9+BuiuloLWDwxl6c3L99tud1tx7mMxqfEC18J50jJwco0GPYZLvd9lHB+lGCFKaFaTAlJkqTiGbBIkqTijdSteVg3Qbv1LY/5XKls8zCi8LybdLfmT+x5v8a6Nf/phd+a2fySLSySJKl4BiySJKl4rY10q+UZJXXUdFO2I6uqTfN2/rWVbp/HG3ZOOiV01j3v39i19qALvjnSvkfEbsArgAcA9wW2Aw7MzA01lj0VOLzPS1/LzD/omff2wGuAI4BdgR8C/zszT1lsO450K0mS9gIOA74JfAY4ZInL3wA8smfa9X3mOwl4MvBq4FvAY4B/jojbZ+Y7hm3AgEWSJH0+M3cGiIi1LD1g+XVmfnXYDBGxD/BM4KWZ+Q/V5A0RcTfg+Ig4NTNvGrj8vKWESmmCHNTMWsr+jWLYIFKz3lw/KnvVaZJmNY3W9Pdo0imhT93rAY1dax/5vf9e9r5XAct6lpYSWpuZOywy36uB44CdM/PqrumPBD4J/Flmnjlo+bkuui0hWBlmmoKVYWblR7Fp0zTaqWbPrHwvDfq3FhHXLfYY06a3j4ifRsSvI+LSiHhLRGzfM899gCu7g5XKd7teH8iUkCRJWo7vAN8GzgFW0qlleQHwsIh4SGb+sppvNXBNn+Wv6Xp9IAMWSZJaFCuaS3YslpYZh656lAVnRcT5wMnAk4D3ds/ebxVDXtti7mpYJEmzOap2UzWAk65h+fQ+D2rsWvuIc78+8RqWAevYhk7PoXdk5vOqae8H9s/Mu/fMuwtwJXBUZr5p0DrnuoZFkiSNxULgdEvXtHOBXSOiN/Wzb/V8zrAVGrBIktSiFSujsUdBnkonxuju6nw6nUDmqT3zHg5cB5w9bIWN1LBYpd0x7i6DHuf+BvW6mpVm7ra6orbVxX5Wu95OKgVT93diVr4fs5jaaktEPKH680HV8/4RsRNwY2Z+vJrnEoDMXFP9+3eA9wDvBy6kU3T7COD5wFeA0xbWn5nnVF2g3xARwa0Dxz0VeH5mbh62fxbdStKY1R3KQGrZB3v+va56vhRYM2CZ/wGuBo4CdqHTgnIR8PfA32fmr3rmfxZwOfDSav6LgL/JzHcutnMGLJIktSgKSeVk5qI7stCy0vXva4HHL2EbvwCOrh5LYi+hMSrxpmGD9mlWm+HVLpvrp9+4U4MlniOT7iV09v0e3Ni19sBvfa2M6GcMLLqdAaPerXkQgxVJ/Yz7P17SMKaEJElqUay07aAOj5IkSSpeETUs03yH4nGyW+lsKK07ut+3co37u1divUiJJl3DsuGB+zV2rT3gG1+Z2RoWU0KSJLWosAHfimVKSJIkFa/xFhabHJvj8ZsNJaSBunlelau0c0UqiSkhSZJaFCtMCdVhSkiSJBVvaC+hHx/zrC0vHrfurK1eG6VZubTeEpNkz5tyDes1M8/n7Dxo4vOd1Dkyyd+QeU/tT7qX0Bf2e0hjvYQe9pUvzWxzzUgpoXk8gZvkhU/ScvgbMltWOHBcLR4lSZJUPAMWSZJUvCJGuq1r2I2yjl530Ja/bS4V3Dbn312HtdwarGFKr4UoRdN1Ek0fv3mv42hbm9+HSdewfHn/P2rsWvuHn/v8zNaw2MIiSZKKZ8AiSZKK10hKqLSm7Gm7udugpudpex9anra+R3VTH6V9z1WWWTo/Jp0S+sof799YSmi/z37OlJAkSVJbDFgkSVLxpqqXkFSSWWoCl3SrSaeEvvbIAxu71j74U2ebEpIkSWqLAYskSSreSPcSkiRJzYiVM5vFadTMBCzWE2gxTZ8jnmeSNDmmhCRJUvFmpoVFkqRptGKFKaE67NYsSTQ/srQjVU+vSXdr/sajH97YtfaBZ35mZqMfU0KSNGYGK9LymRKSJKlFsdK2gzqKDFjq3oxNkprS9G+Nv11SswzrJElS8YpsYZEkaV6scOC4WmxhkSRJxSuyhcXcb7lKry8qff9G0Vb32N7tHr3uoC1/O8rv0jU50vKwc2JWvgOOXj5ZEbEb8ArgAcB9ge2AAzNzwyLLrQReDBwE7APsCFwKfBB4c2Ze3zXvGuDiAat6VGZ+Yti2igxYJGmWdF98pV6F3EtoL+Aw4JvAZ4BDai63DfBa4P3AycDVwIOAo4FHRcR+mfmrnmVOBE7rmbZodG3AIkmSPp+ZOwNExFrqByybgd0zc2PXtA0RcRVwKvAo4IyeZS7NzK8udQcbCVjG2XRns6AWMytN4PNs1r/nq488cevz9Mjlra87PbeU16RBMvOWEZf7NbCxz0tfr553G3mnelh0K0lj1ltzInWLlSsaexTkj6vnc/q89uqI+EVE3BgRZ0fEw+us0JSQJEkzIiKuW2yezNxhzPuwB3As8LnM/ELXSzcD7wQ+CVwJ7A68DPhURByameuHrtebH2rS2mr+X+52vZldswYVoo56TjSdGpz1NNUkjfLZ9J4fk/wMJn3zw+888VGNXWvv+8FPbFpsnsUClqqGZT01egn1WfauwOeB7YH9MvPyRebfFvg2sCIz9xo2ry0skqSizXrA2OTAceNuPRkmIlYDnwZ2APZfLFgByMyfR8SHgFdGxF0z82eD5jVgkSRJyxIRq+h0h74bnZaZC5aw+ELxzdDCXwMWSZI0sojYkU7Lym/RCVbOXcKy2wKHAj/s6Rp9GwYsmrjj1p215e9h3TuXW0PQZg58FPPWPXur86CB99t0d97Sz5c2jPqd6v58R/leT9t3ealiRREDxxERT6j+fFD1vH9E7ATcmJkfr+a5BCAz11T/3gY4C/g94AXAthHxB12rvXwhNRQRb6HTmvJl4GfAGuAlwB7A2sX2z4BFkiRBZzj9buuq50vpBBf97MKtAc7b+7x+TNd6zgWeBRwO3BnYBHwFeG5mfmmxnTNgkSSpRSsKGT8lMxdt6lloWen69yVArSaizPwX4F9G2TdYpFvzj4951pYXZ60JbinszqpJq9tsPmxAMs/TjmHphFGO87DjOm9pvVk16W7N5z394Ma6Nd/73WeUkV8agzLCOkmSpCFMCUmS1KJC7tZcvNoj3c5KU2db72Oem+5NqUnTZVZ+70c16ZTQ955xSGMpoXv960dmNvoxJSRJkopnSkiSpBYVdpflYnmUJElS8cZ6t+Z5q12Yt/crSXVM252vJ13D8v2/WttYDcve/3z6zNawmBKSJKlFscJkRx0eJUmSVLyxtrDMW0pk3t5vXbN+4zJpqUrvNrzcEYBNj2scTAlJktSiUu4lVDqPkiRJKt5Yewk1bdRmxtKbX+fNpHoMTNvowss9LtPWDN/9fo9bd9aWv0fd76a/572pzAWmNJdulM+mzfN50r2EfvCcJzR2rb3HSR+yl5AkSW0oPfheLgeOq8ejJEmSimfAIkmSijdVNSwlKqE+ZpprF7rrAez+LPA8UPsmXcNy4Yv+orFr7Z5v/cDM1rDYwiJJkopnwCJJkopnL6FlKiH9UsI+NGHem/5LSC9OUtPdmqVp5b2E6vEoSZKk4hmwSJKk4tlLqGDT1vtH5Zq2UX8nZd7ScKpn0r2ELn750xq71u5+wnvsJSRJktQWAxZJklQ8ewlJktQi7yVUz9wFLCXkrOvWphy97qBJ7E5jJnUXZi2d9Rn9lXZcxl235ndU08ywTpIkFW/uWlgkSSrJCgeOq2WuuzXParfhEtJe0jiMcm733kxxFkfV9TvfrEl3a/7RK5/R2LX2t9/wr3ZrliRJaospIUmSWmQvoXrmOiWk/nqb0O1NUJbSenrMampVs2uxc3bSKaHLXvPMxq61v/W6U0wJSZIktcWARZKkFsXKFY09Rt6HiN0i4q0R8cWIuCEiMiIOWMLyD4iIz0TEjRFxbUR8ICJ+s898t4+IYyLi0oi4OSLOjYhn1tmGAYskSdoLOAy4AfjMUhaMiHsBG4AAngD8NXA/YENEbN8z+0nAK4ATgYOAjwP/HBHPXnQ701TD0ltbMUgJeX0NNuhzHPfnttyun8Nqe+xWujxN1OWUVtuj6TXpGpbLX/vXjV1rdzvmnSPte0SsyMxbqr/XAuuBAzNzQ41l/x14KHCPzLyxmrY3cC7wqsx8YzVtH+Ac4KWZ+Q9dy78PeBRw98y8adB2bGGRJKlFsWJFY49RLQQrS3EA5LUAABjASURBVN73iNsDjwE+tBCsVOv7PvBV4NCu2dcCCbynZzWnAjsCfzxsWwYskiRpVHsA29BpOen1XeA+Xf++D3BlZl7dZz565r2NiY7Dstxm89JTBqqnreb67s90lPTBsPk8X5aniXNinOfVvKWb/C2cXhFx3WLzZOYODW5ydfV8TZ/XrgG2iYhtMnNzNe+g+brX1ZcDx0mS1KIZGThuWB1ODvi7d9rQWh4DFkmSZkTDrSd1bKye+7WOrAI2dxXSbqR/2mdYK80WEw1YSm9aLH3/SlBi03jvqJULhn2epex7ySY54nGJ51W3Jvap6TTLoPU18bmV8FvoiNtT4yJgM/0DkX3ZurblXOBJEbE6Mzf2zAf962C2mIl2KEkq2aCgetzrm5WL/Ky8j0FKGDhuVJn5S+BjwKERse2W9xRxT2A/4D+7Zj+dzlgtT+1ZzeHAdcDZw7ZlSkiSJBERT6j+fFD1vH9E7ATcmJkfr+a5BCAz13Qt+lrgv4CPRMQJwHbA64FLgLcvzJSZ50TEqcAbIiKAb9HpEv1U4PlVYe5ABiySJAnggz3/Xlc9XwqsGbRQZp4XEQcCbwT+A/gl8EngZZl5fc/szwIuB14K7EInpfQ3mfnOxXZuqka6lWaJOXqpHaXdrfmnb3pBY9faXY78R+/WLEmS1BYDFkmSVLyhNSzzPNphd3P9cevO2uq1eTsWs6iEc7vpFNBizdyTMs4uypPqGjzu+dSu0j6b5dwDaJ54lCRJUvEMWCRJUvFa69ZcetPpVk3Z65od9Gncuo/t0esO2vK3vVCWp5SUS7cSv0fjPM9GeY8lfm5Stxm5l9DYeZQkSVLxDFgkSVLxHOlWkqQWmRKqx5FuJS1L3RoRa0k0DuOo45r0SLdXv+1ljV1rd3rhWxzpVpIkqS2mhCRJapEDx9VjwFIwb443nUpPfTS9f3WXL+04aDZ4Xs0PwzpJklQ8W1gkSWrRipUr296FqWDAosaUOOqqOvxs+hvnjRolNcuUkCRJKp4tLJIktciB4+rxKEmSpOI1MtJtaXngaesOXNrxmyRrKySVZtIj3V538qsaG+l2h785fmZHujUlJElSi0wJ1eNRkiRJxfPmh5KK0jsSb7d5SxuaMm3HpFNCm055TWPX2rs883WmhCRJUvO8l1A9HiVJklS84ltYSuxBU+I+TUqJTdSDUgil7N9yTVuvt+Walc+t16DvzrDPd9ixKO27OG/nqSav+IBFkqbdsLqcbnUv8nXX1xaDlaWxl1A9HiVJklQ8AxZJklQ8uzWrVU3UA/Xmzpe7Pi1unuu4Sldabcs0mnS35hved2xj19rtn/J3M9ut2RYWSZJUPAMWSZJUvMZTQjZHSrPJNFCH3Xdn36RTQje+/3WNpYS2O+w1poQkSZLaYsAiSdKci4jtI+JtEfGTiNgcEd+IiENqLHdJROSAx/d75h0037Nr7aO9hMoyrNnddNts6R38q+nPtLTzZdzvV/NpHOf5pFNCPz/tDY1da7d90itH2veI+BRwf+BI4GLgCOApwMGZeeaQ5e4H3LFn8r7AycAbM/Nvu+ZN4DSgN496UWZetdg+OtKtJElzLCIeDTwCeHxmrq+mnQ3sAbwFGBiwZOa3+qzvydWf/9JnkSsz86uj7KcpIUmS5tvjgE3AhxcmZCf98i5g74i4d90VRcQdgCcDX8zMC5rcSQMWSZLatGJlc4/R3Ac4LzNv6Zn+3a7X61oLrKZ/6wrA06samZsi4msR8cS6Ky6yhqW03HsT3RhLe0+aDXY11jwa9+/pxGtYPvim5ro1P/GoTYvNk5k7dP87Ii4ALsjMx/RMvwdwAfDczDypzvYj4hPAQ4BdM/PGntfeSye9dBlwN+C5wP7AizPzrYut2xoWSZI0LGiqFVBFxG7AI4F/7Q1WADLzqT3zfwjYALwuIk7OzM3D1m/AIklSm1Y0V53R23pS00Y6aZxeq6rna2qu5wg6pSaD0kFbycxbqlaXh9FJO3192PxFBiylpUxGbWofZ7PlKOtuulup3VTbZxpoOk1TKm9Sv2NLWb+/NY07Fzg0Ilb01LHsWz2fs9gKIiLoBCzfz8wvL2HbC9Fab/3MwBklSdJ8Wg/sABzcM/3pwPmZeV6NdewP7EnN1hWAiFhBZ6yX6+kETUMV2cIiSdK8iJUj9+5pypnA2cApEbGazsBxhwMPBR67MFNEbAD2z8x+Rcl/CfwKeHe/DUTEy4HfBT4L/ATYFXhOtY3nZeZNi+1kkb2EJC3fpHqmmRpUm2ZhpNvNp/9DY9fabda+ZNSRbn8DOB54Ap3WlvOAYzPz9K55NtAnYImIOwNXAp/OzMfSR0QcTGcU3b2r9d8I/DdwYmaeUWcfbWGRJGnOZeb/AM+vHoPmOWDA9OuB7RZZ/xlArcBkEAMWSZLaNPqAb3PFoltJklS8qaph6c2VdzNvPj0mVVvRO0LxcevOWtZ2h414PE3dVGdV05+Bn2m72jz+E69h+cjbmqthOeSFE933STIlJElSm0wJ1WJKSJIkFW+qUkLSYnrTNt2W26w8rPtuE2muQU3gTdx8cx6UfoPRJvevxK7kpR//pZh4Suijb28uJfSY55kSkiSNZlj9nRQN3ktolnmUJElS8UwJLdO89SSYt/c7zCwei1lq1pdGNemU0E1nntTYtfZOj36OKSFJkjQG9hKqxZSQJEkqni0skiS1yRaWWqxhkbSFNSxl8fNox8RrWD55SnM1LH/yzJmtYTElJEmSimdKSJKkFjkOSz0GLDWUOKrkMDYjl6vpz6bp9c3i+TLs+1u3a3pbXdjrfh5+5zUPDOskSVLxbGGRJKlN9hKqxV5CAmZz1FZpKaYt9avxmXQvoZs/++7GrrV3/OOn20tIkiSpLaaEJElqkymhWmxhkSRJxZuqGpbuOguw1mJaTapeprcm4eh1By1ru+OucZi3OqLu93vcurO2/G3tyOwZpdt1m9+HidewbHhfczUsBzxlZmtYTAlJktSiWGlKqA5TQpIkqXhTlRJq2qx2Y5y31IK0FLM4Kuwsvqc2TTol9IsvfKCxa+0dHvYXpoQkSdIYeC+hWjxKkiSpeK2lhGzCHJ/eVNcCj/Ot2kqbjXO7s5ri7NbW70bd7fq7NhsmnhL60r83lxJ6yBNNCUmSpDFw4LhaTAlJkqTiGbBIkqTi2a25y6zknM2jl8vPZjLm7TjP2/sdt0nXsPzya6c3dq29/YPXzmwNiy0skiSpeAYskiSpeENTQs+ONVteXO6N43pNUxPmtKWOSrupXIk3rZyV7sWDurAP226Jn0ddbXVHH/Z71Xs8F0zTcW3atP1m9pp4SujrH2kuJfSgQ0ba94jYHjge+HNgB+Bc4NjM/Mgiy60DXtvnpZ9m5q595n8h8Hzgd4DLgf8LnJCZtyy2j3ZrliRJ64H7A0cCFwNHAOsj4uDMPLPG8o8Ebuj69y96Z4iI1wDHAK8HPgv8YfX3KuBvF9uAAYskSXMsIh4NPAJ4fGaur6adDewBvAWoE7B8IzOvG7KN1cCrgX/KzL+rJm+IiO2AIyPinzLz8qH7Oc+9hDQ9pimFqMkoPe1Q+v5psEmnhH713x9r7Fp7uwf82ZL3PSLeSScVtKo7NRMRfw2cDOyTmecNWHYdnZTQjosELE8B3gs8IDO/2TX9HsAFwPMy8/8M20+LbiVpzAxWVLj7AOf1qSP5btfri/leRPw6In4SEe+MiJ37bCPp1MZskZk/ADbX2YYpIUmSZkREDGzlWJCZO/RMWk2nlaPXNV2vD3Ih8CrgW3TqVh5Cpw7m4RHxgMy8tmsdP8/Mm/us49pFtgEYsEiS1K4y7iU0LC018LXMfE/PpM9GxFeBTwLPA1633G0smKqAZZQunCrPKPUoo3y+TXffnebuwCVabhfl0r/zpe+fZlOf1pM6NtK/hWNV9XxNn9eG7cOnIuInwH4929guIu7Yp5VlxzrbsIZFkqT5di5wr4jojQn2rZ7PGWGdK4DumphzgQD26Z4pIvYCtqmzDQMWSZLatGJFc4/RrKczWNzBPdOfDpw/qIfQIBHxJ8AuwFe7Jn8cuBl4Ws/shwO/As5YbL1DU0KTGg10nGmBJkxb98RBqbNS9ntS+zHsnC29m3Td7940p0lH+U0Z5Tdp2Pe36fOg9POqBG2NVqyhzgTOBk6pxku5mE4g8VDgsQszRcQGYP/MjK5p3wLeDZwP/JLOYHAvB34IvH1hvszcGBFvAI6OiE3V9vYDjgJOzMzLFtvJqaphkSRJzcrMjIi1dIbmP55Oa8t5dAaSW6zl4/vAc4G7A7cHLgP+GTiuz7gsxwKb6BTjvhK4gs4YLm+ss58GLJIktShWtt9LKDP/h849fp4/ZJ4D+kw7bAnbSODE6rFkI410O8km1nkzyePnZ6VR2ayvWTbpkW5/fc5nGhvpduV9Hj7RfZ8ki24lSVLxTAlJktSmMgaOK54tLJIkqXi1a1isd9CoHCFWpbIWp795Py4Tr2E5b0NzNSz3PmBma1hMCUmS1CZTQrWYEpIkScUb2sIyzjTQNKWYpnmk26PXHbTl77aadktsUp6m82+YaR7pdlKa+P7OyvlSV4nfWcmUkCRJLYrR7wE0VzxKkiSpeCONdCtN2rw1yWvypi31q/GZdC+hWy74UmPX2hX3fIi9hCRJ0hjYS6gWU0KSJKl4BiySJKl41rBoSUqsJekelfO4dWdt+buU/VuutmorJrndEs+rkllvM14Tr2G58L+aq2HZ8/dntobFFhZJklQ8AxZJklS8saaEem96N4vN9dPEmxAunamK5Zn3m+g1aZLHsu55P2i+aU9ZTTwldNE3mksJ7fFAU0KSJEltMWCRJEnFK76XkE3yUjtMIWopZum3etIpoV9f/M3GrrUrd7+/KSFJkqS2GLBIkqTieS8hSZLaFLYd1FF8DUvTSsiz1u3yN201BHZh7a+Ec04Cv6N1TbyG5ZJvN1fDsua+1rBIkiS1xZSQJEltipltFGnU3KWEZoVpBo1DiedVifs0T+YxjTTxlNCl32kuJfQ7vzez0Y8pIUmSVDxTQpIktWmFbQd1TDQlVFrT7rT1whmncTb7epylrZWeZuntybig93e79PcxqomnhC77f82lhH5rX1NCkiRJbTElJElSi9KB42rxKEmSpOLZrVlTrzvffvS6g7b8PSs59bojI0/zdkurbxuHeXiPs2LSNSy/uvzcxq61t9ttH2tYJEmjGVTEKgGdewk19Rh1FyK2j4i3RcRPImJzRHwjIg6psdxfRcRHIuLSarkfVOu5a595c8Dj2XX20RoWSZK0Hrg/cCRwMXAEsD4iDs7MM4csdwxwNvBK4MfAvYHXAodExH0z87qe+U8Depu/L6qzg6aEptSw7oQ2PZfLz2Z56h6/YfP5GWgxE08J/fh7zaWEfvNeS973iHg08DHg8Zm5vpoWwBeA1Zk5sIkwInbOzKt6pu0PbABemJn/2DU9gbdm5tZjXdRkSkiSpDa1nxJ6HLAJ+PDChOy0ZrwL2Dsi7j1owd5gpfL16nm3UXeoHwMWSZJmRERct9ijz2L3Ac7LzFt6pn+36/Wl+OPq+Zw+rz29qnW5KSK+FhFPrLtSU0JakhJHtpz1XkKljBRsKkXzYuIpoSvOb+xae/vf3HvTYvNk5g7d/46IC4ALMvMxPdPvAVwAPDczT6qz/YhYBXwD+BWwb2be3PXae4EzgcuAuwHPBfYHXpyZb11s3RbdSpLUpgYHjusNRpay6IivbRER2wKnA6uAP+oOVqp9e2rP/B+iU+vyuog4OTM3D1u/KSFJkubbRmB1n+mrqudrFltBRGwDfAS4H/DozPzuIotQpaDeC2xPjbSTAYskSfPtXOBeEbdp6tm3eu5Xi7JFRNyJTsHufsBjMvPLS9j2wjZ762duux1rWCRpvrU1mnKpJl3D8ssrL2yuhmXXPUfp1vxnwEeBtZn54a7pnwd2zsy9hyx7RzrByv7AwZn56SVsdwXwOeD3qu3cNGx+a1gkSZpvZ9IZ/O2UiFhNZ+C4w4GHAo9dmCkiNgD7Z2Z3UPQh4CDgWOCGiPiDrtd+lpkXVsu+HPhd4LPAT4BdgedU23jeYsEKGLBIkjTXMjMjYi1wfPXYATiPzkByZyyy+ELPor+rHt3eRWfEXIDz6QQ/a6v13wj8N3BIjW0AY0gJjdL1cZLNkSV0yy2lm+ogo+7foPulzHvzcgnnnKaH50v7Jp4S+unFzaWEdtndmx9Kksar9z8Lkm5lwCJJkopnLyEtyhFOJU2DpsoLJp4SuuqS5lJCO68xJSRJktQWAxZJklQ8uzVLktSmBu8lNMusYZGkGTVtI9iWUi838RqWn/2ouRqWu/62NSySJEltMSUkSVKL0pRQLbVTQk001ZXS3DfIOPdvku+9tJEyS2+WHrZ/bR3L0r8rutUo58hyRwQv5ZwYtE+lf+cXM+mU0C+uvryxlNAddtrNlJAkSVJbTAlJktSmFbYd1GEvoTEqLcU0rJm29JtWajqVmMaYJt3ppuPWnbXVayUfz7o3UC31N2TiKaFrrmguJbTq7qaEJEllKOXCLk2SKSFJktpkL6FaPEqSJKl41rBobpTW3VuaR9NQLzfxGpZrr2yuhmXHXWe2hsWUkCRJbTIlVItHSZIkFW+iKSG7OLarbldD3cpzVuMwynk1yZTmcs/7Ursr1zXxlNB1VzWXEtphZ1NCkiRpDEwJ1eJRkiRJxSu+l5BN8rPNnjvlmmQKcZ6/59OePplFE08Jbbq6uZTQXXYyJSRJGi+DlfmUpoRq8ShJkqTi2cIiSVKbbGGpxW7N0hCes7Otrc+39POqe/+OXnfQlr9765hmtfvzpGtYbr7+usautXe88w4zW8NiWCdJkopnSkiSpDbFzDaKNKr4bs3SqEZtbh7U1brU5mt1lJ5mUUfdz6n3+9Zt1m5+ePMNm5pLCW1/l5mNfkwJSZKk4pkSkiSpTfYSqsWUkOaGo+pqVJ47HaOmRaetN9HEU0I3Xt9cSmi7O5sSkiRJaosBiyRJLcpY0dhjVBGxfUS8LSJ+EhGbI+IbEXFIzWX3jIjTI2JTRFwfEWdGxL0HzPvCiLggIm6OiAsj4siIejtuwCJJktYDTwFeA/wZcB6wPiIePWyhiNgZ+AKwBjgcOAxYBXwuInbrmfc1wD8AHwAOAk4BXg8cX2cHG69hMdd7q3F2sxxl3cPywE2vTwK7Gs+CefxtmHQNy00/v7GxGpY7bbvdkve9Cko+Bjw+M9dX04JOILI6Mwf2MY+INwEvAPbMzCuqaauBi4H3ZeZzuqZdDpycmS/qWv71wJHA7pl5+bD9tIVFkqQ2xYrmHqN5HLAJ+PDChOy0ZrwL2HtQeqdr2U8tBCvVshuBM4DHd833p8CdqnV2O5VOj+VF008GLJIkzbf7AOdl5i0907/b9fptRMQ2wJ7AOX1e/i6wc5UyWlhHAud2z5SZPwA2D9rGVtuzW7MkzbdpT+E0beIpoQavtdtsu+2mxebJzB26/x0RFwAXZOZjeqbfA7gAeG5mntS7noi4O/Bj4BWZeULPa38NnAzcOzO/FxEnA0/OzO37rOfHwBcz80nD9tuB4yRJalHDAdJ1Iy43LGhaLKCqu+xytmHAIknSrOhtPalpI7C6z/RV1fM1A5a7lk6gUWfZjcB2EXHHzLy5Z94dh2xjCwOWgnX3uILZ6XXVdE+y3uNUxyjbbavZ3POgnfXNk3Gfy/YWK965wKERsaKnjmXf6rlfjQqZuTkiLqJ//cm+wM8y86qubQSwD/DNhZkiYi9gm0Hb6GbRrSRJ8209sANwcM/0pwPnZ+Z5iyz7yIjYdWFCRKyq1vWfXfN9HLgZeFrP8ocDv6LTq2goW1gkSZpvZwJnA6d0jaFyOPBQ4LELM0XEBmD/zOyuuTmBThByZkQcQyf4eE31vGVAuMzcGBFvAI6OiE3V9vYDjgJOzMzLFttJAxZJkuZYZmZErKUTYBxPp7XlPDoDyQ1t+cjMn0bEw+gELu+hk7n5AvBHmfmjntmPpTPey/OAVwJXAK8F3lhnP+3WrLlhjYO0PE3crXmU5Uc16nYn3a1Z9VjDIkmSimfAIkmSijcXKaF561Jn6kNtmrfvm6Zfb+roHXmJKaEC2cIiSZKKZ8AiSZKKNxcpoWkyrDl9Uk3tw0ZWtbl/ekzqs6rbc8Qb7Gla2EuoTLawSJKk4hmwSJKk4hmwSJKk4lnDokVZtyLNHr/Xg1nDUiZbWCRJUvEMWCRJUvFMCUmSitZEl/ilpMBMCZXJFhZJklQ8AxZJklS827W9A5IkjYM9oWaLLSySJKl4BiySJKl4BiySJKl4dmuWWjJqV03z8tJ42a25TLawSJKk4hmwSJKk4s1FSsgmdEmaHeP+TTclVCZbWCRJUvEMWCRJUvEMWCRJUvGG1rBIkiSVwBYWSZJUPAMWSZJUPAMWSZJUPAMWSZJUPAMWSZJUPAMWSZJUvP8PG+Iq9Cn6AQUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "heatmap(adjfb, title=\"HSfbgraphfull\")\n",
    "heatmap(adjfr, title=\"HSfriendsgraphfull\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "len W =  82\n",
      "<class 'list'>\n",
      "<class 'list'>\n",
      "{1, 97, 33, 44, 55}\n",
      "<class 'list'>\n",
      "<class 'list'>\n",
      "[0]\n",
      "[0] [1] [1, 2, 131, 12, 15, 149, 22, 23, 28, 29, 31, 37, 39, 44, 45, 55, 60, 63, 69, 77, 78, 85, 89, 94, 97, 100, 105, 106, 108, 115, 124] [  2  12  15  22  23  28  29  31  37  39  44  45  55  60  63  69  77  78\n",
      "  85  89  94  97 100 105 106 108 115 124 131 149]\n",
      "[0] [1] [1, 97, 33, 44, 55]\n",
      "[ 0  2  3 26 22  7  1 19 28 23 18 13 10 29  8 12 14 24  6  5 25 31 30 27\n",
      "  9  4 17 16 15 11 21 20]\n"
     ]
    }
   ],
   "source": [
    "W = np.arange(coremap.shape[0]); print(\"len W = \", len(W))\n",
    "x = 1\n",
    "W = np.setdiff1d(W, [22])\n",
    "S = [0, 7, 17, 24]\n",
    "\n",
    "NBDS = vnsgm_ordered([x], S, adjfb, adjfr, h, ell, R, gamma) \n",
    "\n",
    "print(NBDS['corr'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAHMCAYAAADoEsNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAY5UlEQVR4nO3dfYxld3kf8O8zWwIGq9p4U7dV+WPlRK1N7SqEtiKVCQYZpTX4ndBSFExDE7mqK0UqiVRB01KllkpDlTpp66ihitUXmhZYjMuGNG1sqARRcaTKsTe8tLA0JDSkXmzzshCIf/1j7tJhfGfm7t3Z+8yZ+/lIV9c+5577O2fundmvfs8556kxRgAAOmx07wAAsL4EEQCgjSACALQRRACANoIIANDmD+228qtnzy51Sc2PPveq5fYGAJrdO05X9z6sEzMiAEAbQQQAaCOIAABtBBEAoI0gAgC0EUQAgDaCCADQRhABANoIIgBAG0EEAGgjiAAAbQQRAKCNIAIAtNm1++6yXXR/+iu/udR2Uxhv2bGYbwqf+YWMtyz7eTAc5uObyu/esqbwGbDJjAgA0EYQAQDaCCIAQBtBBABoI4gAAG0EEQCgjSACALQRRACANoIIANBGEAEA2ggiAEAbQQQAaCOIAABtdu2+u6ypdHVcZXfGKexjcvj3c9X8PPeXn+czTaXL7FS69q6rqnp+kh9L8qIk353keUleNsZ4aMHtX5TkbUlenOT3k/xykr89xvjtvbY1IwIAfFeS1yb5UpL/ej4bVtVVSR5KUkleneSHk7wwyUNVdele21+UGREAYFI+NMa4PEmq6pYkN53Htm9N8sUkN44xvjx7j0eTPJbkbyb5R7ttbEYEANbcGOPpZbarqmcleVWSd50LIbP3+1iSX0ty+17vYUYEAA6Jqnpir9eMMY7u45BXJLkkyaNz1j2S5I693kAQAYBGd9bxsY9v9+Q+vtcijs2ez8xZdybJJVV1yRjj7E5vIIgAwCGxz7Md5zX0kuucIwIALO3x2fOxOesuS3J2jPHV3d5AEAEAlvWpJGeTXD1n3TWZf+7ItxBEAKDRkdq/x6qNMb6e5P1Jbq+q555bXlV/Msn3JnnPXu/hHBEAIFX16tl//rnZ80ur6juSfHmM8Uuz15xOkjHG8S2b/r0k/z3J+6rqp7J5V9Z/mOR0kn+217iCCACQJP9x2////dnzZ5Ic32mjMcapqnpZNm9c9u4kX0/yn7N5i/cv7jWoIAIAjY5UQ01ljjHGnjuybSZk6/KPJnn5MuM6RwQAaHNRZkSm0kV3mfGm0sVTR879pQPo/jrsn98y3+tl99HPkqlTmgGARh1XuxwkSjMAQBtBBABoozQDAI0OylUzXcyIAABtBBEAoI3SDAA0ctUMAEATQQQAaKM0AwCNXDUDANBEEAEA2ijNAECjdb9q5qIEkSl00V12vKl0npxKt1g/z/0dbyqdSnWx7jeFv7esB6UZAKCN0gwANHLVDABAE0EEAGijNAMAjdZ9RmDdjx8AaCSIAABtlGYAoJGrZgAAmggiAEAbpRkAaLTuvWbMiAAAbQQRAKDNRSnNrLrL4iq7SK66Y+Vh76I7lW64y/J9AfbiqhkAgCZOVgWARk5WBQBoIogAAG2UZgCgkZNVAQCaCCIAQBulGQBo5KoZAIAmgggA0EZpBgAauWoGAKCJIAIAtDlQpZkpdA6dSpfZqYy3aqv+jq16vKl0M57K92wKf5OWNZVj69jPe8fplY7nqhkAgCaCCADQ5kCVZgBg3SjNAAA0EUQAgDZKMwDQyA3NAACaCCIAQBulGQBo5KoZAIAmgggA0EZpBgAauWoGAKDJgZoRmULH2GXHmkoX3WVNpZvqVH6eU+lmfNit8nOYyndsKvvJdByoIAIA68ZVMwAATQQRAKCN0gwANHLVDABAE0EEANZcVV1aVfdU1eeq6mxVPVxVNy247e1V9eGq+sLs8ZGqes2iYwsiANDoSO3f4wKcSPK6JG9J8sokp5KcqKobdtuoqu5I8q4kv5Pkr84ev53kF6vqhxYZ2DkiALDGZmHj+iS3jTFOzJY9mOSKJG9PcnKXzf9aks8kec0Y4+nZtr+c5FNJXp/kX+01vhkRAFhvtyZ5Msn95xaMMUaS+5JcWVUv2GXbryf50rkQMtv26SRfSvK1RQY3IwIAjfbzqpmqemKv14wxjm5bdHWSU1vDxMwjW9fv8HY/m+Q9VfXmJD+XpJL8SJI/leRNi+yzIAIA6+1Ykk/MWX5my/q5xhj3z05q/TdJfnK2+MtJfmCM8YFFBhdEAOCQmDPbsfCmy6yrqlck+XdJ3pnk3UmOZPOk13dW1avHGO/fa2BBBAAabfTf0OzxzJ/1uGz2fGbOulRVZfM8kl8dY9y5ZdUHqur5SX4myXoEkSl0kZxCZ+ELseqOnIfdVDqcTuX7yTP5nWWLx5LcXlUb284TuWb2/OgO2/3RJH88ycNz1j2c5Lqqes4Y46u7De6qGQBYbyeSHE1y47blr0/y8THGTieqfiHJV5P8+TnrXpzk8b1CSHJIZkQAYKrqAu9Etg9OJnkwyTuq6liSTye5I8m1SW4+96KqeijJS8cYlSRjjK9V1b1JfrSqfj6bNzY7ks0Ac202b462J0EEANbYGGNU1S1J7p49jmbzct3bxhgP7LH5m5J8LJuX7L46ydPZvALnB5P820XGF0QAYM2NMZ5KctfssdNrrpuz7A+yef+Qn1t2bEEEABpt9JdmWjlZFQBoI4gAAG2UZgCgUR1Z7zmB9T56AKCVIAIAtFGaAYBGB+CGZq3MiAAAbQ7UjMhhbsJ0mI8t0cTsoJhCA0jm++mv/KbPgbV0oIIIwLoSQtaXG5oBADQRRACANkozANCoNtZ7TmC9jx4AaCWIAABtlGYAoJGrZgAAmggiAEAbpRkAaKTXDABAE0EEAGijNAMAjerIes8JHKggMoUOrss2plp1V9RV7+eqTWU/V+2wf+7LmsLxrfpvhE7NHBTrHcMAgFYHakYEANaNG5oBADQRRACANkozANCoNpRmAABaCCIAQBulGQBotLHmNzRb76MHAFoJIgBAG6UZAGhUbmgGANBDEAEA2lyU0syquzOusovkqrt4TqFr6JRMpTut/TwYpnB8uvbON5XvWKI0Y0YEAGgjiAAAbVw1AwCN3NAMAKCJIAIAtFGaAYBGrpoBAGgiiAAAbZRmAKDRxobSDABAC0EEAGijNAMAjcoNzQAAetQYY8eVd9bxnVfuYirdGdk/U+p0uYypfDf97u0vP8/9M6Wf5b3j9ErPHv31G1+x1L+187zogV+Z3JmvSjMA0GjDDc0AAHoIIgBAG6UZAGik1wwAQBNBBABoozQDAI3c0AwAoIkgAgC0UZoBgEZuaAYA0EQQAQDaCCIA0Kg2at8eS+9D1aVVdU9Vfa6qzlbVw1V104LbVlX9SFX9elV9paqeqKpfq6q/sMj2F+Uckal0nlxlx9hlfya62s636p/LlDqH8kyr/vyWGW/VvwtTGW9Zh/1v50VwIsn3JPnxJJ9O8oYkJ6rqxjHGyT22/fkktyd5W5IPJ3lekhfNnvfkZFUAaLTRfB+RqrohyfVJbhtjnJgtezDJFUnenmTHIFJVt2cztFw7xvjIllXvX3R8pRkAWG+3Jnkyyf3nFowxRpL7klxZVS/YZdu/leRD20LIeRFEAGC9XZ3k1Bjj6W3LH9my/hmq6llJXpzkN6rq7qr63ar6RlU9VlV3LDq40gwANNrP7rtV9cRerxljHN226FiST8x56Zkt6+c5luTZSe5I8tkkdyV5Iskbk/xCVX3bGONf7rU/gggAMJZYd66q8pwkN4wxPpMkVfVfsnl+yU8kEUQAYF3Mme1YxOOZP+tx2ez5zJx1SfKFbIaUj50LIbN9GFX1gSR/t6ouH2N8frfBBREAaHQAuu8+luT2qtrYdp7INbPnR+dtNMY4W1X/c4f3PFdv2n7eyTO0Hz0A0OpEkqNJbty2/PVJPj7GOLXLtu9JclVVHT+3oKoqyV9K8qkxxv/da3AzIgCw3k4meTDJO6rqWDZvaHZHkmuT3HzuRVX1UJKXjjG2nl37j5O8LskHquqt+f8nq74oyV9ZZHBBBAAa1UZvcWJ2TsctSe6ePY4mOZXNG5w9sMe2j1fVS7IZSP55kkuS/EaSW8cY711kfEEEANbcGOOpbF5+e9cur7luh+Wnk/zAsmM7RwQAaGNGBAAadfea6XZRgsiqux4u29VxCp1Rp7CPF2IqnUPZX1P5GzGF78tUuuhO4WdJj/WOYQBAK6UZAGh0AG5o1mq9jx4AaCWIAABtlGYAoJHSDABAE0EEAGijNAMAjbp7zXRb76MHAFoJIgBAG6UZAGhUR45070IrMyIAQBtBBABoc6BKM1PoNLvqjpXsr1V/DjqOzrfqDq5+b59pKsfWsZ/3jtMrHc8NzQAAmggiAECbA1WaAYB1s+GGZgAAPQQRAKCN0gwANHLVDABAE0EEAGijNAMAjZRmAACaCCIAQBulGQBoVG5oBgDQ46LMiKy6s+ayltnPqXRv1VmT83HYv2fLWuXfJH9b5tPB+vBTmgGARq6aAQBoIogAAG2UZgCgkdIMAEATQQQAaKM0AwCNNpRmAAB6CCIAQBulGQBopNcMAEATQQQAaKM0AwCN1v2GZhcliEylW+Iy+zmVzsLLWnVHzqn8XJjP5wdcqPWOYQBAK6UZAGi07qWZ9T56AKCVIAIAtFGaAYBGbmgGANBEEAEA2ijNAECjjSNHunehlRkRAKCNIAIAtFGaAYBGbmgGANBEEAEA2hyo0syyHVynYNXHdti76B7m70oyna62vtf7Ox7750I+g3vH6f3bkQUozQAANBFEAIA2B6o0AwDrRq8ZAIAmgggArLmqurSq7qmqz1XV2ap6uKpuOs/3qKr61aoaVfXTi26nNAMAjQ7IVTMnknxPkh9P8ukkb0hyoqpuHGOcXPA9fjjJlec7sCACAGusqm5Icn2S28YYJ2bLHkxyRZK3J9kziFTVn0jytiRvTPKu8xn/QMQwAKDNrUmeTHL/uQVjjJHkviRXVtULFniPf5HkQ2OMd5/v4GZEAKDRASjNXJ3k1Bjj6W3LH9m6fqeNq+q1SV6WZJHA8gyCCAAcElX1xF6vGWMc3bboWJJPzHnpmS3rdxrvO5L80yRvHmP81qL7uVV7DAMA2o0l192TzZNbf3bZgc2IAECj/byh2ZzZjkU8nvmzHpfNns/MWZeqekWSv5zk5Un+cFVtXf3sqjqa5EtjjG/sNrgZEQBYb48luaqqtmeCa2bPj+6w3Z/OZo54KMkXtjyS5M7Zf1+/1+C1eWLsfHfW8d2mY+CbdCk9GHSnnW8q+8nBcO84XXu/av98+Z0/uW//1j7vtW85732vqlcm+U9Jbhlj3L9l+YeSXD7GmHtvkKp6fpLvmrPqwSTvzma55pExxtwZlXOUZgCgUW0c6d6Fk9kMD++oqmPZPOfjjiTXJrn53Iuq6qEkLx1jVJKMMT6b5LPb32xWovnsGOOhRQYXRABgjY0xRlXdkuTu2eNoNi/XvW2M8cDFHl8QAYA1N8Z4Kslds8dOr7luwfc6r/KQIAIAnfpLM61cNQMAtBFEAIA2SjMA0Gkfb2g2Ret99ABAK0EEAGijNAMAjeqIq2YAAFoIIgBAG6UZAOi05jc0uyhBZNlOl8tadefQZUyl2+iqTaUr6lT2c9UOe9feZS1zfIf52GA3SjMAQBulGQDotOalGTMiAEAbQQQAaKM0AwCNSq8ZAIAegggA0EZpBgA6uWoGAKCHGREA6GRGBACghyACALRRmgGARut+H5EDFUSm0n1yGTpWrqepfO6r7oY7hY7ZUzGVTsZT6ZjN6q13DAMAWh2oGREAWDuumgEA6CGIAABtlGYAoJPSDABAD0EEAGijNAMAjeqI0gwAQAtBBABoozQDAJ3WvNfMeh89ANBKEAEA2hyo0sxUOoDyTFPpdLnq79hUfi6H/fiWNYXjm8pnN4WfZRs3NAMA6CGIAABtDlRpBgDWTSnNAAD0EEQAgDZKMwDQyQ3NAAB6CCIAQBulGQBo5KoZAIAmgggA0EZpBgA6Kc0AAPS4KDMiU+lqu8pukFPpkLlqq/65HPbxlrXq41u1w/x7dNg/u2Ud5s/8sFGaAYBObmgGANBDEAEA2ijNAECjOuKqGQCAFoIIANBGaQYAOrmhGQBAD0EEAGijNAMAnZRmAAB6CCIAsOaq6tKquqeqPldVZ6vq4aq6aYHt/npVva+qPjPb7pOz9/kji44tiABAo9rY2LfHBTiR5HVJ3pLklUlOJTlRVTfssd1bkzyV5O8k+YtJ/kmS1yT5aFUdXWTgi3KOyGHvVMr+8V05GFbdidXnMF1T+Z29kO/0veP00ttO0SxsXJ/ktjHGidmyB5NckeTtSU7usvkLxxif3/L/H6yqU0keSvKDSX5mr/HNiADAers1yZNJ7j+3YIwxktyX5MqqesFOG24LIed8dPb8/EUGd9UMAHTax6tmquqJvV4zxtheMrk6yakxxtPblj+ydf157MbLZ8+PLvJiMyIAsN6OJTkzZ/mZLesXUlWXJbknySeT/IdFtjEjAgCHxJzZjoU3XXLdN1XVc5O8N8llSb5vjPG1RbYTRACgU7UXJx7P/FmPy2bP82ZLvkVVXZLkfUlemOT7xxiP7LHJN7UfPQDQ6rEkV1U9IxFdM3ve9VyPqnpONk90/d4krxpjfPh8BhdEAGC9nUhyNMmN25a/PsnHxxg7nqhaVc/OZjnmJUluHmN88HwHV5oBgE79pZmTSR5M8o6qOpbk00nuSHJtkpvPvaiqHkry0jFGbdn2XUm+P8k/SPKlqnrxlnW/N8b4X3sNLogAwBobY4yquiXJ3bPH0WxernvbGOOBPTZ/1ez5J2aPre5L8oa9xhdEAGDNjTGeSnLX7LHTa66bs6zmvPS8CCIA0Gj0l2ZarffRAwCtBBEAoE1t9rWZ7846vtDd1Lbr6JZIr1V3U/Vd4bA57H83p3R8947TF3zew/n4g9P/Y6l/a+c5cvy7V7rv+8GMCADQRhABANq4agYAOtXkqin7yowIANBGEAEA2ijNAECnjfWeE1jvowcAWgkiAEAbpRkAaKTXDABAE0EEAGijNAMAnZRmAAB6HKgZkSl1Z+RbLfsZrLprr+/YtK3681vl93MK+5is/nfB7+zhd6CCCACsHaUZAIAegggA0EZpBgA6Kc0AAPQQRACANkozANBIrxkAgCaCCADQRmkGADopzQAA9BBEAIA2SjMA0Kmqew9a1Rhjx5V31vGdVx4Ch7mzpg6Z++swd329ED6/flP5riyr43foOZdcstJk8PXPn963f2ufdfnxyaUapRkAoI3SDAB0ctUMAEAPQQQAaKM0AwCN9JoBAGgiiAAAbZRmAKDTxnrPCaz30QMArQQRAKCN0gwAdHLVDABAD0EEAGhzUbrvrrob5BQ6ZDLfVLra+o7tr6l0jF3l92wqHbPXwaq77/7+E5/ft+6733b0ct13AQAWJYgAAG1cNQMAnVw1AwDQQxABANoozQBAo6E0AwDQw4wIAHQyIwIA0EMQAQDaKM0AQKea3F3Z95UZEQCgzUWZEVl1Q7JlafrUbyrN63xX5lt1o7apfH5T2c8p6Gg4ee84vfIx15nSDMABIISsMVfNAAD0EEQAgDZKMwDQyC3eAYC1VlWXVtU9VfW5qjpbVQ9X1U0LbvudVfXeqnqyqr5YVSer6gWLji2IAAAnkrwuyVuSvDLJqSQnquqG3TaqqsuT/Lckx5PckeS1SS5L8sGqev4iAyvNAECn5tLMLGxcn+S2McaJ2bIHk1yR5O1JTu6y+ZuSfHuSPzvG+J3Zth9J8ukkb07yN/Ya34wIAKy3W5M8meT+cwvGGCPJfUmu3KPMcmuSXzkXQmbbPp7kgSS3LTK4IAIAh0RVPbHXY85mVyc5NcZ4etvyR7asnzfWJUm+M8mjc1Y/kuTyWelmV0ozANDoOZdcsp/NZuYFjb0cS/KJOcvPbFk/z7cnqS2v22nbz+82uCACAIfEGOPospsuue5Ct1WaAYA193jmz3pcNnueN+ORJF/IZtBYZttvEkQAYL09luSqqmdcvnPN7HneOSAZY5xN8qnMP4fkmiS/N8bYtSyTHLDSzGHuWHmYj61DR0fOZfjc19Myn99h70Q9lfHW1Ikkb0xyY7ZcOZPk9Uk+PsY4tce2d1XVHxtj/J8kqarLZu/1zkUGNyMCAOvtZJIHk7yjqn6oql5WVb+Q5NokP3buRVX1UFVtP+fjp7J56e/Jqrq5ql6Z5P1JvpHk7kUGF0QAYI3N7hlyS5J/n83w8EtJ/kw2b3D2wB7b/m6SlyT5rST/OskvZvPKne8bY/zvRcY/UKUZAGD1xhhPJblr9tjpNdftsPyTSW5edmwzIgBAG0EEAGgjiAAAbQQRAKCNIAIAtBFEAIA2gggA0EYQAQDaCCIAQBtBBABoc1Fu8a5zKIuaSifPw24qnUpX/flNofPrYf9dOOzHhxkRAKCRIAIAtBFEAIA2gggA0EYQAQDaCCIAQBtBBABoI4gAAG0EEQCgjSACALQRRACANoIIANBGEAEA2ggiAECbGmN07wMAsKbMiAAAbQQRAKCNIAIAtBFEAIA2gggA0EYQAQDa/D+rNCaGTmJUMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAHMCAYAAAD22/P2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfN0lEQVR4nO3dfbRnVX3f8ffnTq2MYNcwo2hakoxAG1FYrUldDYkGn58q4yA2XSYIRlsfgtFEDTZRghDFGGWJ1Ky0GO0YzUqImkFtoNHojEqURBNSyoyikYdKxKfBoQKX8YFv//idiz9//O79nblzmbO59/1a66wz95y9zzmcf/is/d1n/1JVSJIktWxu6AeQJEmaxcAiSZKaZ2CRJEnNM7BIkqTmGVgkSVLz/slSJ++Yn1/WJ0R7fvdXl/c0q9imMy8Y+hEkST0csn59hn4G3Z0jLJIkqXkGFkmS1DwDiyRJap6BRZIkNc/AIkmSmmdgkSRJzTOwSJKk5hlYJElS8wwskiSpeQYWSZLUPAOLJElqnoFFkiQ1z8AiSZKat+SvNS/3V5cP9i8T/+r9jl1Wvwtu/9wKP4kkSbonOMIiSZKaZ2CRJEnNM7BIkqTmGVgkSVLzDCySJKl5BhZJktQ8A4skSWqegUWSJDXPwCJJkppnYJEkSc0zsEiSpOYZWCRJUvMMLJIkqXlL/lrzvcXB/NVlfxlakrSaJHk88BzgBOBHgZuBvwHOrqr/06P/0cD5wGMZDYR8EnhlVe2e0valwEuAHwduBP478OaqunPWfRxhkSRpbXsR8GPAW4CnAi/v/v5Mkp9eqmOSIxgFlM3A6cCzgY3Ax5McOdH2Nd09/gR4MvAO4PXAeX0eclWMsEiSpGU7o6q+Pn4gyYeB64BfB05Zou8rgcOBf1tVX+n6frrr+2rgxd2xTd3fb6uq3+r67kxyKHBmkrdV1Y1LPaQjLJIkrWGTYaU7thf4InDk3Xv8kJOBjyyEla7vHuBDwDPH2j0FOAR410T/bYwGT7bMek4DiyRJq0SSvbO2ntd5IHAccPUSbdYDRy/S5irgiK5kRHetAnaNN6qqLwLz3fklWRKSJGlAL8rmWsHL3XKgF0gS4CJGgxpvXqLp4UAYTdKdtHBsE/D1bn97Ve2b0vZb3fklGVgkSVolqmrDClzmTcBW4Jeqqs8nrksFrlrk3/tzDcCSkCRJ6iR5PfAK4GVVtW1G828xChrTRkc2dvuFkZY9wKFJ7jul7eFMH6X5IQYWSZJEknOB3wTOrKoLZ7WvqnngWqbPPzke+MbYhN5djMpHD5+45zHAepaYK7PAwCJJ0oDWZeW25UpyNnAWcFZVvWk/um4HnpjkwWPX2gicBPzZWLvLgH2MFqgbdzrwPUZfFS3JOSySJK1hSV4BvBb4n8BfTiwWt6+qruza7QROrKrxaPRmRiHk0iTnMAofr+n2dy0IV1V7krwBOCvJLcAORivrvgq4oKq+POs5DSySJK1tJ3X7p3fbuBsYrWI7VVV9LcmjGQWXd/ODpfl/rqr+70Tzcxl9xXQG8BvAV4CzgTf2eUgDiyRJA1qXA6jlrICqesyBtOvWUnlGj/4FXNBt+805LJIkqXn3yAjLav5F43vDM0qStNpYEpIkaUAH8nXPWmJJSJIkNc/AIkmSmmdJSJKkAQ39ldC9hSMskiSpeQYWSZLUPEtCkiQNyK+E+nGERZIkNc/AIkmSmmdJSJKkAfmVUD+OsEiSpOYZWCRJUvMsCUmSNCC/EurnHgksq/kXjVfzL1FLktQqS0KSJKl5loQkSRqQXwn14wiLJElqnoFFkiQ1z5KQJEkDcuSgH9+TJElqnoFFkiQ1z5KQJEkD8iuhfhxhkSRJzTOwSJKk5lkSkiRpQP6WUD+OsEiSpOYZWCRJUvOaKgkd7F9CXs79/NVlSdJK8iuhfhxhkSRJzWtqhEWSpLXGSbf9OMIiSZKaZ2CRJEnNsyQkSdKAnHTbjyMskiStcUmOTPLWJJcnuTVJJXlMz761xPa/xtptXqLdU2bdxxEWSZJ0DPBs4O+AjwJb9qPvCVOOPQk4B7hkyrkLgIsnjs1cM8TAIknSgBr5SugTVXUEQJKt7EdgqaorJo8l+U1gHvjjKV1umNZnFktCkiStcVV150pdK8mDgKcC76+qW1bqugYWSZJWiSR7Z20H4TFOZ1TBeeci51+d5DtJbkuyI8nj+1zUkpAkSQNahV8J/RJwLbBz4vg+4O3Ah4GvAg8BXgF8JMkpVbV9qYsaWCRJWiWqasOQ90/yM8BDgbOqqsbPVdVNwAvGDl2e5P3A3wNvApYMLJaEJEnSSnkecCewrU/jqrodeB9wdJIHLtV2yRGWTWde0PP5VsbB/iVkf3lZkjS0Rr4SOmBJDgV+HvhwVd24H10XBk+WnPjrCIskSVoJ/wG4P4tPtr2bJPcDTgH+oar2LNXWOSySJIkkz+r++chuf2KSBwC3VdVlXZvrAapq85RLPA/YA3xgkeufz2ig5FPAN4DNwK8BRwFbZz2fgUWSpAE1VBJ678Tfr+32NzAKF4tKcjTwaODCqvrOIs12AS9k9Nnz/YFbgE8Dv1xVfzXr4QwskiSJqpoZnRYZWaGqvgQs2b+q3sl+lIsmOYdFkiQ1zxEWSZIGtAoXjrtHOMIiSZKaZ2CRJEnNsyQkSdKAGvpKqGmOsEiSpOYZWCRJUvMsCUmSNCC/EurHERZJktQ8A4skSWqeJSFJkgbkV0L9OMIiSZKaZ2CRJEnNsyQkSdKA/EqoH0dYJElS8wwskiSpeZaEJEkakF8J9eMIiyRJap6BRZIkNc+SkCRJA/IroX4cYZEkSc0zsEiSpOZZEpIkaUBzloR6cYRFkiQ1z8AiSZKaZ0lIkqQBxZXjenGERZIkNc/AIkmSmmdJSJKkAc1ZEurFERZJktQ8A4skSWqeJSFJkgaUdY4d9OFbkiRJzTOwSJKk5hlYJEkaUNZlxbZlP0NyZJK3Jrk8ya1JKsljevbd1rWf3K6Y0vY+Sc5JckOSfUl2JXl+n/s4h0WSJB0DPBv4O+CjwJb97H8r8MSJY9+e0u73gV8AXg1cCTwd+IMk96mq/7bUDQwskiTpE1V1BECSrex/YPl+Vd1tRGVckocDzwdeXlVv6Q7vTPIjwHlJtlXVHYv1tyQkSdKA5tZlxbblqqo7V/A/aTFbgQLePXF8G3A48LilOhtYJElaJZLsnbXdQ7c+LMnXkny/m59yfpLDJtocB3y1qr45cfyqsfOLsiQkSZIOxP8G/h64GljHaC7LrwCPTvKzVfXdrt0m4OYp/W8eO78oA4skSQPK3MoVO6pqw4pdrP893zJx6C+SXANcBPxH4D3jzaddYolzd7EkJEmSVtp7gDuBE8aO7WH6KMrCsWmjL3cxsEiSpJW2MAN4fDLvLuDBSSZDy/Hd/uqlLmhgkSRpQC18JXQPOJVRxhj/1PkSRkHm1Im2pwN7gR1LXdA5LJIkiSTP6v75yG5/YpIHALdV1WVdm+sBqmpz9/ePM/pM+Y+BLzGadPsE4CXAp4GLF65fVVcn2Qa8IUn4wcJxpwIvqar5pZ7PwCJJkgDeO/H3a7v9DcDmRfr8P+CbwKuABzEaQbkW+B3gd6rqexPtXwjcCLy8a38t8IKqevush0vV4pNy75ifX3LGriRJq80h69cf1NrKjkf8uxX7f+1jr/zrpupCK8k5LJIkqXkGFkmS1DznsEiSNKCsc+ygD9+SJElqnoFFkiQ1z5KQJEkDamzBt2Y5wiJJkppnYJEkSc2zJCRJ0oAyZ0moD0dYJElS8wwskiSpeZaEJEka0JwLx/XiW5IkSc0zsEiSpOZZEpIkaUBx4bheHGGRJEnNM7BIkqTmWRKSJGlAloT6cYRFkiQ1z8AiSZKaZ0lIkqQBuXBcP74lSZLUPAOLJElqniUhSZIG5FdC/TjCIkmSmmdgkSRJzbMkJEnSgObmLAn14QiLJElqnoFFkiQ1z5KQJEkDigvH9eJbkiRJzTOwSJKk5lkSkiRpQHMuHNeLIyySJKl5BhZJktQ8A4skSQPKuqzYtuxnSI5M8tYklye5NUkleUyPfuuSvCLJh5P8Y5Lbk3wuyblJ7j/RdnN33WnbU2bdyzkskiTpGODZwN8BHwW29Oy3Hjgb+GPgIuCbwCOBs4CnJjmhqr430ecC4OKJY5+bdSMDiyRJ+kRVHQGQZCv9A8s88JCq2jN2bGeSrwPbgKcCH5roc0NVXbG/D2hgkSRpQC0sHFdVdy6z3/eBPVNOfabbH7nsh5ow/FuSJEkrIsneWdtBepTHdfurp5x7dZLvJLktyY4kj+9zQQOLJElaMUmOAs4FPl5Vnxw7tQ94O/BiRoHmRcDhwEeSnDzrupaEJEka0EouHFdVG1bsYsuQ5IHAnwO3AaeOn6uqm4AXjB26PMn7gb8H3gRsX+rajrBIkqQDlmQT8JfABuDxVXXjrD5VdTvwPuDoLuwsyhEWSZJ0QJJsZPQ59I8Aj62qL+xH94XBkyUn/hpYJEkaUObu3b8llORwRiMrP8oorOzaj773A04B/mHi0+i7MbBIkiSSPKv75yO7/YlJHgDcVlWXdW2uB6iqzd3f64G/AP418CvA/ZL89Nhlb1woDSU5n9FoyqeAbwCbgV8DjgK2zno+A4skSQOaa2Adls57J/5+bbe/gVG4mOZB/CDg/N6U8+eMXWcX8ELgdOD+wC3Ap4Ffrqq/mvVwBhZJkkRVzaxNLYysjP19PdCrplVV7wTeuZxnA78SkiRJ9wKOsEiSNKAD+ZXltcQRFkmS1DwDiyRJap4lIUmSBtTCrzXfG/iWJElS8wwskiSpeZaEJEkaUOYcO+jDtyRJkppnYJEkSc2zJCRJ0oAa+i2hpvmWJElS8wwskiSpeZaEJEkakAvH9eNbkiRJzTOwSJKk5lkSkiRpQJaE+vEtSZKk5hlYJElS8ywJSZI0IH9LqB/fkiRJap6BRZIkNc+SkCRJA8q6dUM/wr2CIyySJKl5BhZJktQ8S0KSJA3IheP68S1JkqTmGVgkSVLzLAlJkjSgOReO68W3JEmSmmdgkSRJzbMkJEnSgPxKqB/fkiRJap6BRZIkNc/AIknSgLJubsW2ZT9DcmSStya5PMmtSSrJY/aj/08l+WiS25J8K8mfJPkXU9rdJ8k5SW5Isi/JriTP73MPA4skSToGeDZwK/DR/emY5FhgJxDgWcB/Bh4B7Exy2ETz3wd+HbgAeDJwGfAHSV406z5OupUkSZ+oqiMAkmwFtuxH33OAbwMnVdVt3TWuBnYBZwBv7I49HHg+8PKqekvXd2eSHwHOS7Ktqu5Y7CaOsEiSNKDMza3YtlxVdeeynj25D/B04H0LYaW73ueBK4BTxppvBQp498RltgGHA49b6l4GFkmStFxHAeuBq6ecuwo4buzv44CvVtU3p7Rjou3dWBKSJGmVSLJ3Vpuq2rCCt9zU7W+ecu5mYH2S9VU137VdrN34taYysEiSNKBVsnBc9Tw3rV0tce4uBhZJklaJFR496WNPt582OrIRmB+bSLuH6WWfpUZp7rIqYp0kSRrEtcA804PI8fzw3JZdwIOTTIab47v9tHkwdzGwSJI0oBYWjluuqvou8OfAKUnud9d/U/KvgBOAPxtrfgmjtVpOnbjM6cBeYMdS97IkJEmSSPKs7p+P7PYnJnkAcFtVXda1uR6gqjaPdT0b+Bvgg0neDBwKvB64Hvi9hUZVdXWSbcAbkgS4ktEn0acCL+km5i7KwCJJkgDeO/H3a7v9DcDmxTpV1e4kj2W0QNz7ge8CHwZeUVXfnmj+QuBG4OXAgxiVlF5QVW+f9XAGFkmSBjTXyFdCVZUebTYvcvwzzFj4rWv3HeCsbtsvbbwlSZKkJRhYJElS8ywJSZI0oAP5DaC1xLckSZKaZ2CRJEnNsyQkSdKAVslvCd3jfEuSJKl5BhZJktQ8S0KSJA3IklA/viVJktQ8A4skSWqeJSFJkgbkwnH9+JYkSVLzDCySJKl5loQkSRrQ3Lp1Qz/CvYIjLJIkqXkGFkmS1DxLQpIkDciF4/rxLUmSpOYZWCRJUvMsCUmSNCBLQv34liRJUvMMLJIkqXmWhCRJGpC/JdSPb0mSJDXPwCJJkppnSUiSpAH5lVA/viVJktQ8A4skSWqeJSFJkgZkSagf35IkSWqegUWSJDXPkpAkSQNy4bh+fEuSJKl5BhZJktQ8A4skSQPK3LoV25b9DMlhSS5MclOS+SSfTbKlR7/rk9Qi2+cn2i7W7kV9ntE5LJIkaTvwk8CZwHXAc4HtSU6qqkuX6HcycN+JY8cDFwGXTGl/MXDBxLFr+zyggUWSpDUsydOAJwDPrKrt3bEdwFHA+cCigaWqrpxyvV/o/vnOKV2+WlVXLOc5LQlJkjSkuXUrty3PycAtwAcWDlRVAe8CHprkYX0vlOSfAr8AXF5VX1juA01jYJEkaZVIsnfWNqXbccDuqrpz4vhVY+f72gpsYvroCsBp3RyZO5L8dZKf73thA4skSWvbJuDmKcdvHjvf1/OAW4E/nXLuj4CXAE8CTgPmgYuTvKzPhZ3DIknSkFZw4biq2rDcrss8d5ckRwJPBP5HVd025dlOnWj/PmAn8LokF1XV/FLXd4RFkqS1bQ/TR1E2dvtpoy/TPJdRrlisHPRDuhLUe4DD6FF2MrBIkrS27QKOTTKZCY7v9lfPukCSMAosn6+qT+3HvRfuOTl/ZtGGkiRpAFm3bsW2ZdoObABOmjh+GnBNVe3ucY0TgaPpOboC0AWkXwS+zSg0Lck5LJIkrW2XAjuAdyTZxGjhuNOBRwHPWGiUZCdwYlVlyjWeB3wP+MNpN0jySuAngI8BNwEPBl7c3eOMqrpj1kMaWCRJWsOqqpJsBc7rtg3AbkYLyX1oVv8k9wdOAS6tqq8t0uwaRuFna3f924C/Bbb0uQdARmvDTHfH/HyvmcGSJK0Wh6xfP20E4R4z/8ELV+z/teu3vPSgPvvB5BwWSZLUPAOLJElqnnNYJEka0vJ/A2hNcYRFkiQ1z8AiSZKaZ0lIkqQBZQV/S2g18y1JkqTmGVgkSVLzLAlJkjQkvxLqxREWSZLUPEdYJEkakiMsvTjCIkmSmmdgkSRJzbMkJEnSgFyHpR/fkiRJap6BRZIkNc+SkCRJQ/IroV4cYZEkSc0zsEiSpOZZEpIkaUiWhHpxhEWSJDXPwCJJkppnSUiSpAFlnSWhPhxhkSRJzTOwSJKk5lkSkiRpSP6WUC++JUmS1DwDiyRJap4lIUmShuTCcb04wiJJkppnYJEkSc2zJCRJ0oBiSagXR1gkSVLzDCySJKl5loQkSRqSC8f14luSJGmNS3JYkguT3JRkPslnk2zp0e+1SWrK9tVF2r80yReS7EvypSRnJumVRRxhkSRJ24GfBM4ErgOeC2xPclJVXdqj/xOBW8f+/s5kgySvAc4BXg98DPiZ7t8bgf8y6wYGFkmSBjT0V0JJngY8AXhmVW3vju0AjgLOB/oEls9W1d4l7rEJeDXwtqr6re7wziSHAmcmeVtV3bjUDSwJSZK0tp0M3AJ8YOFAVRXwLuChSR62Avd4CnBId81x2xgNnswsPxlYJElaJZLsnbVN6XYcsLuq7pw4ftXY+Vk+l+T73RyYtyc5Yso9Ctg1frCqvgjM97mHJSFJkoY0/MJxm4AvTDl+89j5xXwJ+E3gSkbzVn6W0TyYxyf5qar61tg1bq+qfVOu8a0Z9wAMLJIkrRpVtWG5XZdzrqrePXHoY0muAD4MnAG87kDvscCSkCRJa9sepo9wbOz2N085t6iq+ghwE3DCxD0OTXLfKV0O73MPA4skSUOam1u5bXl2AcdOWQ/l+G5/9XL+q4DxOTG7gAAPH2+U5BhgfZ97GFgkSVrbtgMbgJMmjp8GXFNVu/fnYkmeBDwIuGLs8GXAPuA5E81PB74HfGjWdZ3DIknS2nYpsAN4R7deynWMgsSjgGcsNEqyEzixqjJ27ErgD4FrgO8yWgzulcA/AL+30K6q9iR5A3BWklu6+50AvAq4oKq+POshDSySJA0o64b9SqiqKslW4Lxu2wDsZrSQ3KyRj88Dvwz8c+A+wJeBPwB+e8pCcucyWu/lDOA3gK8AZwNv7POcGa0NM90d8/MzZ+1KkrSaHLJ+fWa3Wjnfv/qjK/b/2nXHPf6gPvvB5BwWSZLUPEtCkiQNafiF4+4VHGGRJEnNM7BIkqTmWRKSJGlIloR6cYRFkiQ1z8AiSZKaZ0lIkqQBZfm/AbSm+JYkSVLzDCySJKl5loQkSRqSXwn14giLJElqnoFFkiQ1z5KQJElDimMHffiWJElS8wwskiSpeZaEJEkakiWhXnxLkiSpeQYWSZLUPEtCkiQNqCwJ9eJbkiRJzTOwSJKk5lkSkiRpSJaEevEtSZKk5hlYJElS8ywJSZI0pGToJ7hXcIRFkiQ1z8AiSZKaZ0lIkqQhzTl20IdvSZIkNc/AIkmSmmdJSJKkAflbQv34liRJUvMMLJIkqXkGFkmShpS5lduW+wjJYUkuTHJTkvkkn02ypUe//5Tkg0lu6Pp9sbvOA6e0rUW2F/V5RuewSJKk7cBPAmcC1wHPBbYnOamqLl2i3znADuA3gH8EHgacDWxJ8m+qau9E+4uBCyaOXdvnAQ0skiStYUmeBjwBeGZVbe+O7QCOAs4Hlgosj6iqr4/9/fEku4GdwHOA/zrR/qtVdcVyntOSkCRJQxq+JHQycAvwgYUDVVXAu4CHJnnYYh0nwsqCz3T7I5f7QNMYWCRJWiWS7J21Tel2HLC7qu6cOH7V2Pn98bhuf/WUc6d1c13uSPLXSX6+70UNLJIkrW2bgJunHL957HwvSTYCFwJfBP504vQfAS8BngScBswDFyd5WZ9rO4dFkqQhreDCcVW1Ybldl3nuLknuB1wCbAR+rqr2TTzbqRPt38dorsvrklxUVfNLXd8RFkmS1rY9TB9F2djtp42+/JAk64EPAo8AnlZVV83oQleCeg9wGD3KTgYWSZLWtl3AscndhnqO7/bT5qLcJckhjCbsngA8vao+tR/3Xrjn5PyZRRtKkqQBVOZWbFum7cAG4KSJ46cB11TV7sU6JrkvozLQo4FnVNXH+960C0i/CHybUWhaknNYJEla2y5ltPjbO5JsYrRw3OnAo4BnLDRKshM4saoy1vd9wJOBc4Fbk/z02LlvVNWXur6vBH4C+BhwE/Bg4MXdPc6oqjtmPaSBRZKkNayqKslW4Lxu2wDsZrSQ3IdmdH96t/+tbhv3LkYr5gJcwyj8bO2ufxvwt8CWHvcAIKO1Yaa7Y36+18xgSZJWi0PWr8/sVivnu1+7bsX+X3ufBz3koD77weQcFkmS1DwDiyRJap5zWCRJGlJWbRVnRTnCIkmSmmdgkSRJzbMkJEnSkFbwt4RWM9+SJElqnoFFkiQ1z5KQJEkDOoDfAFpTfEuSJKl5BhZJktQ8S0KSJA1pzrGDPnxLkiSpeQYWSZLUPEtCkiQNya+EevEtSZKk5hlYJElS8ywJSZI0JEtCvfiWJElS8wwskiSpeZaEJEkakiWhXnxLkiSpeQYWSZLUPEtCkiQNqCwJ9eJbkiRJzXOERZKkITnC0otvSZIkNc/AIkmSmmdJSJKkISVDP8G9giMskiSpeQYWSZLUPEtCkiQNya+EevEtSZKk5hlYJElS8ywJSZI0IJfm78e3JEnSGpfksCQXJrkpyXySzybZ0rPv0UkuSXJLkm8nuTTJwxZp+9IkX0iyL8mXkpyZ9EtsBhZJkrQd+EXgNcC/B3YD25M8balOSY4APglsBk4Hng1sBD6e5MiJtq8B3gL8CfBk4B3A64Hz+jxgqmrRk3fMzy9+UpKkVeiQ9esP6kpuK/n/2uU8exdK/hx4ZlVt746FURDZVFXHLtH3d4FfAY6uqq90xzYB1wF/VFUvHjt2I3BRVb1srP/rgTOBh1TVjUs9pyMskiStbScDtwAfWDhQo9GMdwEPXay8M9b3Iwthpeu7B/gQ8Myxdk8BDumuOW4bo/m0M8tPBhZJklaJJHtnbVO6HQfsrqo7J45fNXZ+2r3WA0cDV085fRVwRFcyWrhGAbvGG1XVF4H5xe4xzq+EJEka0AqXoKYFklk2AV+YcvzmsfPTHA5krN1ifb/e7W+vqn1T2n5riXvcxcAiSdIqUVUbltt1mef2p++B3MOSkCRJa9wepo9wbOz200ZQYDQyUj377gEOTXLfKW0PX+IedzGwSJK0tu0Cjp2yHsrx3X7aHBWqah64lunzT44HvlFVXx+7R4CHjzdKcgywfrF7jDOwSJK0tm0HNgAnTRw/DbimqnbP6PvEJA9eOJBkY3etPxtrdxmwD3jORP/Tge8x+qpoSc5hkSRpbbsU2AG8Y2wNldOBRwHPWGiUZCdwYlWNTxJ+M6MQcmmScxiFj9d0+7sWhKuqPUneAJyV5JbuficArwIuqKovz3pIF46TJGnMwV44rgVJ/hmjgPEsRqMtu4Fzq+qSsTY7uXtgIcm/ZBRcHsuocvNJ4JVVtWuiXYCXAWcAPwZ8BbgIeOOUT6rv/owGFkmSfmAtBpZ7A+ewSJKk5hlYJElS8wwskiSpeQYWSZLUPAOLJElqnoFFkiQ1z8AiSZKaZ2CRJEnNM7BIkqTmGVgkSVLzDCySJKl5BhZJktQ8A4skSWqegUWSJDXPwCJJkppnYJEkSc0zsEiSpOYZWCRJUvMMLJIkqXkGFkmS1DwDiyRJap6BRZIkNc/AIkmSmmdgkSRJzTOwSJKk5qWqhn4GSZKkJTnCIkmSmmdgkSRJzTOwSJKk5hlYJElS8wwskiSpeQYWSZLUvP8PDM++KYmP0OcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "heatmap(NBDS['As'])\n",
    "heatmap(NBDS['Bs'])\n",
    "# heatmap(NBDS['As']-shuffle_adj_matrix(NBDS['Bs'], list(NBDS['corr'])))\n",
    "print(NBDS['corr'][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 1. 1. 1. 1. 1.]\n",
      " [1. 0. 0. 2. 0. 0.]\n",
      " [1. 0. 0. 0. 0. 2.]\n",
      " [1. 2. 0. 0. 0. 0.]\n",
      " [1. 0. 0. 0. 0. 2.]\n",
      " [1. 0. 2. 0. 2. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print(NBDS['Bs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 1., ..., 0., 0., 0.],\n",
       "       [0., 1., 0., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 232,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = NBDS['P']\n",
    "# P[np.where(P<1.8e-25)]=0\n",
    "# np.transpose(P)\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cxp =  [ 2 34 45 56 98]\n",
      "labelsGx =  [  1   2   3  13  16  23  24  29  30  32  38  40  45  46  56  61  64  70\n",
      "  78  79  86  90  95  98 101 106 107 109 116 125 132 150]\n",
      "labelsGxp =  [ 1  2 34 45 56 98]\n"
     ]
    }
   ],
   "source": [
    "Cxp = NBDS['Cxp']\n",
    "Cxp.sort()\n",
    "print(\"Cxp = \", np.array(Cxp)+1)\n",
    "\n",
    "labelsGx = NBDS['labelsGx']\n",
    "labelsGx.sort()\n",
    "print(\"labelsGx = \", np.array(labelsGx)+1)\n",
    "\n",
    "labelsGxp = NBDS['labelsGxp']\n",
    "labelsGxp.sort()\n",
    "print(\"labelsGxp = \", np.array(labelsGxp)+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 4 reproduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(82, 82) <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "hhop = 2\n",
    "ell = 2\n",
    "g = 0.1\n",
    "R = 50\n",
    "n_iter = 20\n",
    "voi = 27\n",
    "\n",
    "HSfrgraphcore = igraph.load(\"HSfrgraphcore.txt\", format=\"edges\")\n",
    "HSfrgraphcore.to_undirected()\n",
    "\n",
    "adj2a = np.array(HSfrgraphcore.get_adjacency().data)\n",
    "print(adj2a.shape, type(adj2a))\n",
    "\n",
    "adj2, perm2  = permute_adj_matrix(adj2a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mmp = np.transpose(np.arange(adj2.shape[0]))\n",
    "mmp = np.transpose(np.stack((mmp, perm2)))\n",
    "\n",
    "mmp2 = np.transpose(np.arange(adj2.shape[0]))\n",
    "mmp2 = np.transpose(np.stack((perm2, mmp2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp = [(ii, jj) for ii, jj in zip(perm2, list(range(adj2.shape[0])))]\n",
    "sp.sort(key = lambda x: x[0])\n",
    "sp = np.array(sp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "map2sort = mmp2[sp[:,1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = np.diag(list(range(1, adj2.shape[0]+1)))[perm2]-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproduce figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error2\n",
      "error2\n",
      "error2\n",
      "error 4\n",
      "error 4\n",
      "error 4\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n",
      "Error\n",
      "(66, 66) [17, 0, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 42, 43, 46, 47, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 81] 81\n",
      "error 3\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-a2afc87a6d39>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    133\u001b[0m         \u001b[0mnrvois_h2l2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msi2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrankvoi_h2l2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    134\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 135\u001b[0;31m         \u001b[0mL_h2l3\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvnsgm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvoi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0madj2a\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0madj2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mhhop\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mell\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    136\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mL_h2l3\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'case'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'possible'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    137\u001b[0m \u001b[0;31m#             print(L)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-3-55ef50255cda>\u001b[0m in \u001b[0;36mvnsgm\u001b[0;34m(x, seeds, A, B, h, ell, R, g, pad, sim, verb, plotF)\u001b[0m\n\u001b[1;32m     61\u001b[0m     \u001b[0;31m#print(\"inside 1 w/ shapes \", AA.shape, BB.shape)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     62\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 63\u001b[0;31m     \u001b[0mP\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvnsgm_ordered\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvoi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mAA\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mBB\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mell\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpad\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpad\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mverb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mplotF\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mplotF\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     64\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     65\u001b[0m     \u001b[0mP\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-3-55ef50255cda>\u001b[0m in \u001b[0;36mvnsgm_ordered\u001b[0;34m(x, S, A, B, h, ell, R, g, pad, sim, verb, plotF)\u001b[0m\n\u001b[1;32m    179\u001b[0m             \u001b[0;31m#print(\"inside 2 w/ shapes \", As.shape, Bs.shape)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 181\u001b[0;31m             \u001b[0mcorr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msgm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_predict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mAs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mBs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseeds_A\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mSx1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseeds_B\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mSx1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m#multistart_sgm(A, B, R, len(Sx1), g, pad=pad, n_iters=n_iters)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    182\u001b[0m             \u001b[0mP\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msgm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mP_final\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    183\u001b[0m             \u001b[0mCx2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mCx2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-12-4e12f9c4b5c5>\u001b[0m in \u001b[0;36mfit_predict\u001b[0;34m(self, A, B, seeds_A, seeds_B)\u001b[0m\n\u001b[1;32m    339\u001b[0m             \u001b[0mThe\u001b[0m \u001b[0moptimal\u001b[0m \u001b[0mpermutation\u001b[0m \u001b[0mindices\u001b[0m \u001b[0mto\u001b[0m \u001b[0mminimize\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mobjective\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    340\u001b[0m         \"\"\"\n\u001b[0;32m--> 341\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mB\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseeds_A\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseeds_B\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    342\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mperm_inds_\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    343\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-12-4e12f9c4b5c5>\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, A, B, seeds_A, seeds_B)\u001b[0m\n\u001b[1;32m    275\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    276\u001b[0m                 alpha = minimize_scalar(\n\u001b[0;32m--> 277\u001b[0;31m                     \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbounds\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"bounded\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    278\u001b[0m                 ).x  # computing the step size\n\u001b[1;32m    279\u001b[0m                 \u001b[0mP_i1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mP\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mQ\u001b[0m  \u001b[0;31m# Update P\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.7/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize_scalar\u001b[0;34m(fun, bracket, bounds, args, method, tol, options)\u001b[0m\n\u001b[1;32m    788\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdisp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbool\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    789\u001b[0m             \u001b[0moptions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'disp'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdisp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 790\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_scalar_bounded\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbounds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    791\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'golden'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    792\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_scalar_golden\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbracket\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m_minimize_scalar_bounded\u001b[0;34m(func, bounds, args, xatol, maxiter, disp, **unknown_options)\u001b[0m\n\u001b[1;32m   1844\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1845\u001b[0m         \u001b[0msi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msign\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrat\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrat\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1846\u001b[0;31m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxf\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0msi\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtol1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1847\u001b[0m         \u001b[0mfu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1848\u001b[0m         \u001b[0mnum\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "outdic = {'nrankxp0':[], 'nrankxpBad':[], 'imp2':[], 'potseeds':[], \n",
    "         'numcand':[], 'Nh1':[], 'Nh2':[], 'nrvois':[], 'nrvois_h1l2':[],\n",
    "         'imp2_h1l2':[], 'nrvois_h2l2':[], 'imp2_h2l2':[], 'nrvois_h2l3':[], \n",
    "         'imp2_h2l3':[]}\n",
    "\n",
    "for vvi in range(1, adj2.shape[0]+1):\n",
    "    voi = vvi - 1\n",
    "    \n",
    "    fb_voi = map2sort[voi,0]\n",
    "    fr_voi = map2sort[voi,1]\n",
    "    \n",
    "    Nhtmp = ego(adj2a,hhop,voi,mindist=1)\n",
    "    \n",
    "    \n",
    "    imp2 = 0\n",
    "    nrankxp0 = 0\n",
    "    nrankxpBad = 0 \n",
    "    \n",
    "    nrvois = np.full(len(Nhtmp), None)\n",
    "    \n",
    "    imp2_h1l2 = 0\n",
    "    nrankxp0_h1l2 = 0\n",
    "    nrankxpBad_h1l2 = 0 \n",
    "    \n",
    "    nrvois_h1l2 = np.full(len(Nhtmp), None)\n",
    "    \n",
    "    \n",
    "    for si, seed in enumerate(Nhtmp):\n",
    "        S = map2sort[seed]\n",
    "        frtmp = fr_voi + (1 if S[1]>fr_voi else 0)\n",
    "        #print(S)\n",
    "#         print(\"iter \", ii, \" w shapes \", adj2a.shape, adj2.shape)\n",
    "        L = vnsgm(fb_voi, S, adj2a, adj2, hhop, ell, R, g) # L <- localnbd(fb.voi,S,HSfbgraphcore,g.frscrambled,hhop,ell,R,g)\n",
    "#         print(L)\n",
    "        if L['case'] == 'possible':\n",
    "#             print(L)\n",
    "            if(frtmp in L['labelsGxp']):\n",
    "                try:\n",
    "                    Ps = L['P'][0:len(L['labelsGx']), 0:len(L['labelsGxp'])]\n",
    "                    rankPs = pass2ranksuplus(Ps) \n",
    "                    rankvoi = rankPs[1, np.where(L['labelsGxp']==frtmp)]\n",
    "                    rankvoi = (rankvoi - 1)/(len(L['Cxp'])-1)\n",
    "                except:\n",
    "                    rankvoi = None\n",
    "                    print(\"error2\")\n",
    "            else:\n",
    "                rankvoi = None\n",
    "                imp2 = imp2 + 1\n",
    "            \n",
    "            if rankvoi is not None:\n",
    "                if rankvoi == 0:\n",
    "                    nrankxp0 = nrankxp0 + 1\n",
    "                else:\n",
    "                    if rankvoi >= 0.5:\n",
    "                        nrankxpBad = nrankxpBad + 1\n",
    "        else:\n",
    "            rankvoi = None\n",
    "            imp2 = imp2 + 1\n",
    "            \n",
    "        nrvois[si] = rankvoi\n",
    "        \n",
    "        L_h1l2 = vnsgm(fb_voi,S,adj2a,adj2,hhop,(ell+1),R,g)\n",
    "        if L_h1l2['case'] == 'possible':\n",
    "#             print(L)\n",
    "            if(frtmp in L_h1l2['labelsGxp']):\n",
    "                try:\n",
    "                    Ps = L_h1l2['P'][0:len(L_h1l2['labelsGx']), 0:len(L_h1l2['labelsGxp'])]\n",
    "                    rankPs_h1l2 = pass2ranksuplus(Ps) \n",
    "                    rankvoi_h1l2 = rankPs_h1l2[1, np.where(L_h1l2['labelsGxp']==frtmp)]\n",
    "                    rankvoi_h1l2 = (rankvoi_h1l2 - 1)/(len(L_h1l2['Cxp'])-1)\n",
    "                except:\n",
    "                    print(\"error 4\")\n",
    "                    rankvoi_h1l2 = None\n",
    "            else:\n",
    "                rankvoi_h1l2 = None\n",
    "                imp2_h1l2 = imp2_h1l2 + 1\n",
    "            \n",
    "            if rankvoi_h1l2 is not None:\n",
    "                if rankvoi_h1l2 == 0:\n",
    "                    nrankxp0_h1l2 = nrankxp0_h1l2 + 1\n",
    "                else:\n",
    "                    if rankvoi_h1l2 >= 0.5:\n",
    "                        nrankxpBad_h1l2 = nrankxpBad_h1l2 + 1\n",
    "        else:\n",
    "            rankvoi_h1l2 = None\n",
    "            imp2_h1l2 = imp2_h1l2 + 1\n",
    "            \n",
    "        nrvois_h1l2[si] = rankvoi_h1l2\n",
    "                \n",
    "    Nhtmp2 = ego(adj2a,hhop+1,voi,mindist=1)\n",
    "    \n",
    "    imp2_h2l2 = 0\n",
    "    nrankxp0_h2l2 = 0\n",
    "    nrankxpBad_h2l2 = 0 \n",
    "    nrvois_h2l2 = np.full(len(Nhtmp2), None)\n",
    "    imp2_h2l3 = 0\n",
    "    nrankxp0_h2l3 = 0\n",
    "    nrankxpBad_h2l3 = 0 \n",
    "    nrvois_h2l3 = np.full(len(Nhtmp2), None)\n",
    "    \n",
    "    for si2, seed in enumerate(Nhtmp2):\n",
    "        S2 = map2sort[seed]\n",
    "        frtmp = fr_voi + (1 if S2[1]>fr_voi else 0)\n",
    "        #print(S)\n",
    "#         print(\"iter \", ii, \" w shapes \", adj2a.shape, adj2.shape)\n",
    "        L_h2l2 = vnsgm(fb_voi, S, adj2a, adj2, hhop+1, ell+1, R, g) # L <- localnbd(fb.voi,S,HSfbgraphcore,g.frscrambled,hhop,ell,R,g)\n",
    "#         print(L)\n",
    "        if L_h2l2['case'] == 'possible':\n",
    "#             print(L)\n",
    "            if(frtmp in L_h2l2['labelsGxp']):\n",
    "                try:\n",
    "                    Ps = L_h2l2['P'][0:len(L_h2l2['labelsGx']), 0:len(L_h2l2['labelsGxp'])]\n",
    "                    rankPs_h2l2 = pass2ranksuplus(Ps) \n",
    "                    rankvoi_h2l2 = rankPs_h2l2[1, np.where(L_h2l2['labelsGxp']==frtmp)]\n",
    "                    rankvoi_h2l2 = (rankvoi_h2l2 - 1)/(len(L_h2l2['Cxp'])-1)\n",
    "                except:\n",
    "                    print(\"error 3\")\n",
    "                    rankvoi_h2l2 = None\n",
    "            else:\n",
    "                rankvoi_h2l2 = None\n",
    "                imp2_h2l2 = imp2_h2l2 + 1\n",
    "            \n",
    "            if rankvoi_h2l2 is not None:\n",
    "                if rankvoi_h2l2 == 0:\n",
    "                    nrankxp0_h2l2 = nrankxp0_h2l2 + 1\n",
    "                else:\n",
    "                    if rankvoi_h2l2 >= 0.5:\n",
    "                        nrankxpBad_h2l2 = nrankxpBad_h2l2 + 1\n",
    "        else:\n",
    "            rankvoi_h2l2 = None\n",
    "            imp2_h2l2 = imp2_h2l2 + 1\n",
    "            \n",
    "        nrvois_h2l2[si2] = rankvoi_h2l2\n",
    "        \n",
    "        L_h2l3 = vnsgm(voi,S,adj2a,adj2,hhop+1,ell+2,R,g)\n",
    "        if L_h2l3['case'] == 'possible':\n",
    "#             print(L)\n",
    "            if(frtmp in L_h2l3['labelsGxp']):\n",
    "                try:\n",
    "                    Ps = L_h2l3['P'][0:len(L_h2l3['labelsGx']), 0:len(L_h2l3['labelsGxp'])]\n",
    "                    rankPs_h2l3 = pass2ranksuplus(Ps) \n",
    "                    rankvoi_h2l3 = rankPs_h2l3[1, np.where(L_h2l3['labelsGxp']==frtmp)]\n",
    "                    rankvoi_h2l3 = (rankvoi_h2l3 - 1)/(len(L_h2l3['Cxp'])-1)\n",
    "                except:\n",
    "                    rankvoi_h2l3 = None\n",
    "                    print(\"Error\")\n",
    "                    print(rankPs_h2l3.shape, L_h2l3['labelsGxp'], frtmp)\n",
    "\n",
    "\n",
    "            else:\n",
    "                rankvoi_h2l3 = None\n",
    "                imp2_h2l3 = imp2_h2l3 + 1\n",
    "            \n",
    "            if rankvoi_h2l3 is not None:\n",
    "                if rankvoi_h2l3 == 0:\n",
    "                    nrankxp0_h2l3 = nrankxp0_h2l3 + 1\n",
    "                else:\n",
    "                    if rankvoi_h2l3 >= 0.5:\n",
    "                        nrankxpBad_h2l3 = nrankxpBad_h2l3 + 1\n",
    "        else:\n",
    "            rankvoi_h2l3 = None\n",
    "            imp2_h2l3 = imp2_h2l3 + 1\n",
    "            \n",
    "        nrvois_h2l3[si2] = rankvoi_h2l3\n",
    "        \n",
    "    outdic['nrankxp0'].append(nrvois_h2l3)\n",
    "    outdic['nrankxpBad'].append(nrankxpBad)\n",
    "    outdic['imp2'].append(imp2)\n",
    "    outdic['potseeds'].append(len(Nhtmp)) \n",
    "#     outdic['numcand'].append(len(L['Cxp']))\n",
    "    outdic['Nh1'].append(Nhtmp)\n",
    "    outdic['Nh2'].append(Nhtmp2)\n",
    "    outdic['nrvois'].append(nrvois)\n",
    "    outdic['nrvois_h1l2'].append(nrvois_h1l2)\n",
    "    outdic['imp2_h1l2'].append(imp2_h1l2)\n",
    "    outdic['nrvois_h2l2'].append(nrvois_h2l2)\n",
    "    outdic['imp2_h2l2'].append(imp2_h2l2)\n",
    "    outdic['nrvois_h2l3'].append(nrvois_h2l3)\n",
    "    outdic['imp2_h2l3'].append(imp2_h2l3)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "na =  np.array(outdic['imp2_h2l2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrvois_h2l2_nonone = []\n",
    "for ii in outdic['nrvois_h2l2']:\n",
    "    midres = []\n",
    "    for jj in ii:\n",
    "        if jj is not None:\n",
    "            midres.append(jj[0][0])\n",
    "    nrvois_h2l2_nonone.append(np.array(midres))\n",
    "# print(nrvois_h2l2_nonone)\n",
    "zero =  np.array(list(map(lambda x: len(np.where(x==0)[0]), nrvois_h2l2_nonone)))\n",
    "chance =  np.array(list(map(lambda x: len(np.where(np.array(x)>0.5)[0]), nrvois_h2l2_nonone)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "ssizehnbdx = np.array(list(map(lambda x: len(x), outdic['nrvois_h2l2'])))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "between = ssizehnbdx - chance - zero - na"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[31  0 14 26  9 14 17  1 20  0  5 17 13  2 24 17  0  0 15  0  0  0  0  0\n",
      "  0  0  0 36  0 29  0  4  0 26  0  0  0 23 17 28 32 29  2  9  1 19 47  2\n",
      "  0 32  0 23 37 29 14 15 29 15  0  1 18 28  0  0  0  0 35  0 13 24  0  0\n",
      "  0  6 17  0  2  0  0  0  1  1]\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0]\n",
      "[ 0  1  2  0  4  0  0  0  0  1  2  0  0  0  0  0  3 11  0  1  3  0  2  3\n",
      "  8  1  3  0  1  0  0  0  2  0  4  0  6  0  0  0  0  0  1  2  3  2  1  0\n",
      "  0  0  3  0  0  0  0  0  0  0  4  0  2  1  7  7 10  3  0  4  4  0  2  5\n",
      " 13  8  0  0  0  2  1  2  0  0]\n",
      "[12 29 17  7 28  0  0  0 29 22 41  0  7  0  0  0 33 39  0 22 38 27 45 27\n",
      " 42 33 50  0 46  0  0 13 41  9 47  0 23  0  0  0  0  0 25 23 32  8  3  0\n",
      "  0  0 39  0  6  0  0  0  0  0 38  0 19  9 40 43 48 23  0 34 24  0 18 42\n",
      " 37 39 13 29  0 31 43 41  0  0]\n"
     ]
    }
   ],
   "source": [
    "print(na)\n",
    "print(zero)\n",
    "print(chance)\n",
    "print(between)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fdaf2681290>"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAcDUlEQVR4nO3df3QU9b3/8ecbSAn+oPy2aOQmfMVeREOEoCCKP1DQ1tYfxYqiDRdbsOVa7I/Tq/V7KlTsqUevv/1ev5yrggpYwSpe+21FAUXFCgSDDQSLtmijCBgFoYgk+P7+sUOykIQsIbOzu/N6nLMnO5OZnffuTt75zHtm32vujoiIxEe7qAMQEZH0UuIXEYkZJX4RkZhR4hcRiRklfhGRmOkQdQCp6NGjhxcWFkYdhoikSXl5w/3Bg6OLI9uVl5d/7O4995+fFYm/sLCQlStXRh2GiKSJWcN9/em3npm919R8lXpERGJGiV9EJGaU+EVEYiYravxNqa2tpbq6ml27dkUdSlbJz8+noKCAvLy8qEMRkYhkbeKvrq7myCOPpLCwEEs+EyTNcndqamqorq6mqKgo6nBEJCJZW+rZtWsX3bt3V9I/CGZG9+7ddZQkEnNZm/gBJf1W0GsmIlmd+EVE5OBlbY2/kbYeyep7CkQkR2nEn2Z79uyJOgQRiTkl/lZ68MEHKSkpoaSkhKKiIs4++2wWLlzIsGHDGDRoEJdddhk7duwAEi0nfv3rX3P66aczb948KioqGDp0KMXFxVxyySV8+umnET8bEYkTJf5Wuvbaa6moqGDFihUUFBQwYcIEpk+fzosvvsiqVasoLS3lzjvvrF8+Pz+fV199lbFjx/K9732P2267jbfeeouTTjqJadOmRfhMRCRucqfGH5EpU6Zwzjnn0LVrV9auXcvw4cMB2L17N8OGDatf7vLLLwdg27ZtbN26lTPPPBOAsrIyLrvssvQHLiKxpcR/CGbOnMl7773H/fffzx/+8AfOO+885s6d2+Syhx9+eJqjExFpmko9rVReXs4dd9zB448/Trt27Rg6dCivvfYa77zzDgA7d+7kr3/9a6P1vvrVr9K1a1deeeUVAB577LH60b+ISDrkzog/zZdf3n///XzyySecffbZAJSWljJz5kyuuOIKvvjiCwCmT5/O8ccf32jdWbNmce2117Jz50769u3LI488ktbYRSTecifxp1lzyXrFihWN5m3YsGGf6ZKSEv785z+HEZaISIuU+CV0yZ+t0+fiJKvMSdp5r8ydnVc1fhGRmFHiFxGJGSV+EZGYUeIXEYkZJX4RkZjJmat6bFrbtmX2m1s+g79hwwYuvPBCKisrU3rMmTNnMmrUKI4++uhDDU9EpNU04k+jmTNn8uGHH0YdhojEXKiJ38w2mNlfzKzCzFYG87qZ2Qtmtj742TXMGMJWV1dHWVkZxcXFjBkzhp07d1JeXs6ZZ57J4MGDGT16NBs3bmT+/PmsXLmScePGUVJSwssvv8yll14KwIIFC+jUqRO7d+9m165d9O3bF4B3332X888/n8GDB3PGGWewbt06ALZs2cJ3vvMdhgwZwpAhQ3jttdcAmDp1KhMmTOCss86ib9++3HvvvdG8KBnCrOGWFeZYw00kROkY8Z/t7iXuXhpM3wAscvd+wKJgOmu9/fbbTJw4kbfeeovOnTvzwAMPcN111zF//nzKy8uZMGECN910E2PGjKG0tJTZs2dTUVHB8OHDefPNNwF45ZVXOPHEE1mxYgVvvPEGp556KgATJ07kvvvuq+8L9KMf/QhIdAT9yU9+wooVK3jqqaf4/ve/Xx/PunXreP7551m+fDnTpk2jtrY2/S+KiGS0KGr8FwFnBfdnAS8B/xFBHG3i2GOPrW/FfNVVV/Gb3/yGyspKzjvvPCDxjVu9e/dutF6HDh047rjjqKqqYvny5fz0pz9l6dKl7NmzhzPOOIMdO3awbNmyfVo27+0B9OKLL7J27dr6+Z999hnbt28H4Jvf/CYdO3akY8eO9OrVi02bNlFQUBDa8xeR7BN24ndgoZk58H/dfQZwlLtvBHD3jWbWq6kVzWwiMBGgT58+IYfZerZfHeHII49kwIABvP766y2ue8YZZ/DHP/6RvLw8zj33XMaPH8+ePXu44447+PLLL+nSpQsVFRWN1vvyyy95/fXX6dSpU6PfdezYsf7+7t3tWbWqjo8+gtLSRou2mRZbMkxNfo1y52Pv0nYypa1HpsQRtrBLPcPdfRBwATDZzEakuqK7z3D3Uncv7dmzZ3gRHqL333+/PsnPnTuXoUOHsmXLlvp5tbW1rFmzBkj8U9g7MgcYMWIEd999N8OGDaNnz57U1NSwbt06BgwYQOfOnSkqKmLevHkAuDurV68GYNSoUdx///31j9PUPwcRkeaEOuJ39w+Dn5vN7GngFGCTmfUORvu9gc1tsq0ULr8MQ//+/Zk1axaTJk2iX79+XHfddYwePZof//jHbNu2jbq6Oq6//noGDBjA+PHjufbaa+nUqROvv/46p556Kps2bWLEiMT/w+LiYnr16lV/FDF79mx++MMfMn36dGpraxk7diwDBw7k3nvvZfLkyRQXF1NXV8eIESN48MEHI3n+IpJ9zEM6njGzw4F27r49uP8C8GtgJFDj7r81sxuAbu7+iwM9Vmlpqa9cuXKfeVVVVfTv3z+U2HNF8kuWXOpp69eupcPj5M9YpPMfdNYdtudoJ8hU7P9eRfXeNdpulr8nZlaedGFNvTBH/EcBTwej1w7AHHf/k5mtAJ40s2uA94Hc/MLZmqSs2z3EAruIyEEKLfG7+9+AgU3MryEx6hcRkQjok7siIjGjxC8iEjM506Qtq6j+L3Jg+uxHqDTiFxGJmZxJ/MkNudri1hrjx49n/vz5bfvERETaWM4kfhERSY0S/yF49NFHKS4uZuDAgVx99dUALF26lNNOO42+gy9i/rOLANixYwcjR45k0KBBnHTSSSz4fy8DsOH9D+nfvz8/+MEPGDBgAKNGjeLzzz8H4J133uHcc89l4MCBDBo0iHfffReA22+/nSFDhlBcXMzNN9+ctucaVovjUFsnT7WGm2ScTGib3VIMNs7rb7lEib+V1qxZw6233srixYtZvXo199xzDwAbN27k1Vdf5bk5d3HDLYl+Ovn5+Tz99NOsWrWKJUuW8LOb72bvJ6bXr1/P5MmTWbNmDV26dOGpp54CYNy4cUyePJnVq1ezbNkyevfuzcKFC1m/fj3Lly+noqKC8vJyli5dGs0LICJZS1f1tNLixYsZM2YMPXr0AKBbt24AXHzxxbRr144Tvt6XTZs/ARIN1n75y1+ydOlS2rVrxwcbt7Bpcw0ARUVFlJSUADB48GA2bNjA9u3b+eCDD7jkkkuAxD8OgIULF7Jw4UJOPvlkIHEksX79+vpePyIiqVDibyV3b9SSGfZti+zBZWizZ89my5YtlJeXk5eXR2Gfo9n1xe5Gy7dv357PP/+c5vonuTs33ngjkyZNasunEi1dtieSdir1tNLIkSN58sknqalJjNw/+eSTZpfdtm0bvXr1Ii8vjyVLlvDePzYe8LE7d+5MQUEBzzzzDJD4ApadO3cyevRoHn74YXbs2AHABx98wObNbdLcVERiJGdG/OnuvjhgwABuuukmzjzzTNq3b19ffmnKuHHj+Na3vkVpaSklJSX8a7/CFh//scceY9KkSfzqV78iLy+PefPmMWrUKKqqqhg2bBgARxxxBI8//ji9ejX5XTYiIk3KmcQfhbKyMsrKypr9/Y73Eidee/Tose83ciV9creysrL+/s9//vP6+/369WPx4sWNHnPKlClMmTLlUMIWkZhT4pfUqBYvkjNU4xcRiRklfhGRmFHiFxGJGdX4RURClmnf/6wRv4hIzOTOiH9OG3d6ujID/i2LiIRAI34RkZjJnRF/BDZs2MAFF1zA6aefzrJlyzjmmGNYsGABjz/+ODP+z13srq3juKICHvvdcxx22GH16638e8PXLZZ2jyLy3JBpddOMlnxEnIlHs7n+OZH9nl/U+65G/IeoqbbKl156KStefJTVL8+h//FFPPTQQ1GHKSJSTyP+Q9RUW+XKykr+9w3Xs3Xbdnb883NGX5AXcZQiIg2U+A9RU22Vx48fzzMzb2Xgicczc+7/8NLK9yKMMA1y/TA9C0RVOoi6ZHHIYrrvqtQTgu3bt9P7qB7U1tYxe/6fog5HRGQfuTPib+6EVVInTLqXNr1MG7vllls4dfS/8S8FX+OkE45je21aNisikpLcSfwRKCwsbLat8g+/O6RhwTT9wxERSYUS/4GEdbRwdNLjkhn/FLKyVhvT+qy0INMvXYXI993Qa/xm1t7M3jSz54LpIjN7w8zWm9nvzOwrYccgIiIN0nFydwpQlTR9G3CXu/cDPgWuae0DN/el5NI8vWYiEmriN7MC4JvAfwfTBpwDzA8WmQVc3JrHzs/Pp6amRonsILg7NTU15OfnRx2KiEQo7Br/3cAvgCOD6e7AVnevC6argWOaWtHMJgITAfr06dPo9wUFBVRXV7Nly5YDR/DPjxvub65qfrk0rvvx1obfVW07yMc9CB8nhVAVbCY/P5+CgoLQtimSsqlt3FhRUhZa4jezC4HN7l5uZmftnd3Eok0O2d19BjADoLS0tNEyeXl5FBUVtRzInBMa7h/siZ6Q1j1hWsPv/ObwjlhOSApBB0YisleYI/7hwLfN7BtAPtCZxBFAFzPrEIz6C4APQ4xBRET2E1qN391vdPcCdy8ExgKL3X0csAQYEyxWBiwIKwYREWksiuv4/wN4wsymA28Cal0ZFl3nHh+Z+F5nw/X0LcnE17UNpCXxu/tLwEvB/b8Bp6RjuyIi0piatImIxIxaNkhCjh7SSsgOZr8JqfSTle1GIqYRv4hIzCjxi4jEjBK/iEjMqMYvIrJXTM51acQvIhIzSvwiIjGjxC8iEjNK/DnGrOGWsebYvtd0x8ne5x7X5y8ZQYlfRCRmlPhFRGJGiV9EJGZ0Hb/kjixoA2zjGuLyKyMMJNvNicf19mHRiF9EJGaU+EVEYkalHmlzapObJAvKT1kvJm0W2pJG/CIiMaPELyISM0r8IiIxE+8a//7111yox07NgFYAqrmK7Ns2ZWpUUTRNI34RkZhR4hcRiRklfhGRmEmpxm9mi9x9ZEvzco0+Xi9NOpRzQW11/iMXzkdJZA6Y+M0sHzgM6GFmXYG9e1tn4OiQYxMRkRC0NOKfBFxPIsmX05D4PwMeCDEuEREJyQETv7vfA9xjZte5+31piklEREKUUo3f3e8zs9OAwuR13P3R5tYJykRLgY7BOvPd/WYzKwKeALoBq4Cr3X13q5+BiEgmyoTP1DQjpat6zOwx4A7gdGBIcCttYbUvgHPcfSBQApxvZkOB24C73L0f8ClwTStjFxGRVkj1k7ulwAnuqfdaDJbdEUzmBTcHzgH2XiMzi8Rn2v4r1ccVEZFDk+p1/JXA1w72wc2svZlVAJuBF4B3ga3uXhcsUg0c08y6E81spZmt3LJly8Fu+gAxNdyk9fQ6SuzMsYZblkt1xN8DWGtmy0mUcABw928faCV33wOUmFkX4Gmgf1OLNbPuDGAGQGlpqS5UFhFpI6km/qmHshF332pmLwFDgS5m1iEY9RcAHx7KY4uIyMFJ9aqelw/2gc2sJ1AbJP1OwLkkTuwuAcaQuLKnDFhwsI8tIiKtl2rLhu00lGS+QuJE7T/dvfMBVusNzDKz9iTOJTzp7s+Z2VrgCTObDrwJPNTq6FsjTZdY2fqG+6pTHYQ0tSJo06+HbG3MUbVdOJjtqjVEs2xaw2vjN2fXa5PqiP/I5Gkzuxg4pYV13gJObmL+31paV0REwtOq7pzu/gyJyzJFRCTLpFrquTRpsh2J6/qz69hGRESA1K/q+VbS/TpgA3BRm0eTafQVggKqiWcZtVNvWao1/n8LOxAREUmPVHv1FJjZ02a22cw2mdlTZlYQdnAiItL2Uj25+wjwLIm+/McA/xPMExGRLJNqjb+nuycn+plmdn0YAUlmaPEa5QxuOSsiB5bqiP9jM7sqaLrW3syuAmrCDExERMKRauKfAHwX+AjYSKLlgk74iohkoVRLPbcAZe7+KYCZdSPxxSwTwgpMJBvpUsJwtGmLDUl5xF+8N+kDuPsnNNGOQUREMl+qib+dmXXdOxGM+FM9WhARkQySavL+T2CZmc0n8RHW7wK3hhaViIiEJtVP7j5qZitJNGYz4FJ3XxtqZCLZaL82H6r5SyZKuVwTJHolexGRLNeqtswiIpK9lPhFRGJGV+ZkG7X9FTmw/dup62+mEY34RURiRolfRCRmlPhFRGJGNf40Se41wtSoothXJsaUNumqA0f19Z25VtfW16C2KY34RURiRolfRCRmVOpJl0z8xqpMjElEQqcRv4hIzCjxi4jEjBK/iEjMqMYvIqFQS+rMFdqI38yONbMlZlZlZmvMbEowv5uZvWBm64OfXVt6LBERaTthlnrqgJ+5e39gKDDZzE4AbgAWuXs/YFEwLSIiaRJa4nf3je6+Kri/HagCjgEuAmYFi80CLg4rBhERaSwtNX4zKwROBt4AjnL3jZD452BmvZpZZyIwEaBPnz7pCDNz5NrH7UUyyN5WJR7jP63Qr+oxsyOAp4Dr3f2zVNdz9xnuXurupT179gwvQBGRmAk18ZtZHomkP9vdfx/M3mRmvYPf9wY2hxmDiIjsK8yregx4CKhy9zuTfvUsUBbcLwMWhBWDiIg0FmaNfzhwNfAXM6sI5v0S+C3wpJldA7wPXBZiDAek64xFQqRWyhkrtMTv7q8CzXUBGxnWdkVE5MDUskFEJGbUskFE4qW+BBXf8pNG/CIiMaPELyISM0r8IiIxE+8a/6FcbtaWbRUO5bHi1N4hXc81Tq+pxJJG/CIiMaPELyISM0r8IiIxo8QvIhIzSvwiIjGjxC8iEjNK/CIiMRPv6/jjRtenSyrmNNdUF+1DOUIjfhGRmFHiFxGJGZV6DsaBDoFFspXKN7GjEb+ISMwo8YuIxIwSv4hIzKjGL1nD1jfcVyVa0i2X9j+N+EVEYkaJX0QkZpT4RURiJudr/LlUlxMRaQsa8YuIxIwSv4hIzCjxi4jETGiJ38weNrPNZlaZNK+bmb1gZuuDn13D2r6IiDQtzBH/TOD8/ebdACxy937AomBaRETSKLTE7+5LgU/2m30RMCu4Pwu4OKzti4hI09J9OedR7r4RwN03mlmv5hY0s4nARIA+ffqkKTyR9LFxDRcY+5URBpImcXu+mSxjT+66+wx3L3X30p49e0YdjohIzkh34t9kZr0Bgp+b07x9EZHYS3fifxYoC+6XAQvSvH0RkdgLrcZvZnOBs4AeZlYN3Az8FnjSzK4B3gcuC2v7Ej9qz5HhpiZ/daneoSiFlvjd/YpmfjUyrG2KiEjLMvbkroiIhEOJX0QkZnK+LbNIxlLNWyKiEb+ISMwo8YuIxIwSv4hIzKjGnyRTrgPPlDhEJDdpxC8iEjNK/CIiMaNSj0gM7FM+7BddHJIZNOIXEYkZJX4RkZhR4hcRiRklfhGRmFHiFxGJGSV+EZGYUeIXEYkZXccvkiFsXEODDr8ywkAk52nELyISM0r8IiIxo8QvIhIzqvGLZAp9FaOkiUb8IiIxo8QvIhIzKvVkujlJh/9X6vA/VvTeZ5cser804hcRiRklfhGRmFHiFxGJGdX420oW1fdEQpP8dyCpa+51CymXRDLiN7PzzextM3vHzG6IIgYRkbhKe+I3s/bAA8AFwAnAFWZ2QrrjEBGJqyhG/KcA77j739x9N/AEcFEEcYiIxJK5p7cebWZjgPPd/fvB9NXAqe7+7/stNxGYGEx+HXj7EDfdA/j4EB+jrSmm1Cim1GViXIopNWHE9C/u3nP/mVGc3G3qLEaj/z7uPgOY0WYbNVvp7qVt9XhtQTGlRjGlLhPjUkypSWdMUZR6qoFjk6YLgA8jiENEJJaiSPwrgH5mVmRmXwHGAs9GEIeISCylvdTj7nVm9u/A80B74GF3X5OGTbdZ2agNKabUKKbUZWJciik1aYsp7Sd3RUQkWmrZICISM0r8IiIxk/OJP1PaQ5jZw2a22cwqk+Z1M7MXzGx98LNrGuM51syWmFmVma0xsylRxxRsP9/MlpvZ6iCuacH8IjN7I4jrd8GFAWllZu3N7E0zey4TYjKzDWb2FzOrMLOVwbyo378uZjbfzNYF+9awiPfzrwevz97bZ2Z2fdSvUxDbT4J9vNLM5gb7flr2qZxO/BnWHmImcP5+824AFrl7P2BRMJ0udcDP3L0/MBSYHLw2UcYE8AVwjrsPBEqA881sKHAbcFcQ16fANWmOC2AKUJU0nQkxne3uJUnXf0f9/t0D/Mnd/xUYSOL1iiwmd387eH1KgMHATuDpKGMCMLNjgB8Dpe5+IokLXcaSrn3K3XP2BgwDnk+avhG4McJ4CoHKpOm3gd7B/d7A2xHGtgA4L8NiOgxYBZxK4hONHZp6X9MUSwGJBHEO8ByJDyJGHdMGoMd+8yJ7/4DOwN8JLhrJhJj2i2MU8FomxAQcA/wD6Ebi6srngNHp2qdyesRPw4u7V3UwL1Mc5e4bAYKfvaIIwswKgZOBNzIhpqCkUgFsBl4A3gW2untdsEgU7+PdwC+AL4Pp7hkQkwMLzaw8aHEC0b5/fYEtwCNBSey/zezwiGNKNhaYG9yPNCZ3/wC4A3gf2AhsA8pJ0z6V64k/pfYQcWZmRwBPAde7+2dRxwPg7ns8cWheQKKpX/+mFktXPGZ2IbDZ3cuTZzexaLr3reHuPohEKXOymY1I8/b31wEYBPyXu58M/JP0l5qaFNTKvw3MizoWgOCcwkVAEXA0cDiJ93F/oexTuZ74M709xCYz6w0Q/Nyczo2bWR6JpD/b3X+fCTElc/etwEskzkF0MbO9HzhM9/s4HPi2mW0g0U32HBJHAFHGhLt/GPzcTKJufQrRvn/VQLW7vxFMzyfxjyAT9qkLgFXuvimYjjqmc4G/u/sWd68Ffg+cRpr2qVxP/JneHuJZoCy4X0aizp4WZmbAQ0CVu9+ZCTEFcfU0sy7B/U4k/kCqgCXAmCjicvcb3b3A3QtJ7EOL3X1clDGZ2eFmduTe+yTq15VE+P65+0fAP8zs68GskcDaKGNKcgUNZR6IPqb3gaFmdljwt7j3tUrPPhXFSZY0n0T5BvBXEnXimyKMYy6JWl4tiZHRNSTqxIuA9cHPbmmM53QSh5FvARXB7RtRxhTEVQy8GcRVCfwqmN8XWA68Q+JwvWNE7+NZwHNRxxRse3VwW7N3386A968EWBm8f88AXTMgpsOAGuCrSfMijSmIYRqwLtjPHwM6pmufUssGEZGYyfVSj4iI7EeJX0QkZpT4RURiRolfRCRmlPhFRGJGiV+kDZnZ0WY2P+o4RA5El3OKiMSMRvwiLTCz28zsR0nTU83sZ2Z2e9BL/S9mdnnwu0JL+s4FkUykxC/SsieAy5Omv0uifW4JiZ7z5wK37+39IpLplPhFWuDubwK9gvr9QBJfkFECzPVEJ9FNwMvAkCjjFElVh5YXERESnSbHAF8jcQTwv6INR6T1NOIXSc0TJDpzjiHxT2ApcHnwpTE9gREkmmuJZDyN+EVS4O5rgjbIH7j7RjN7msRX460m0eX0F+7+UfBtZiIZTZdziojEjEo9IiIxo8QvIhIzSvwiIjGjxC8iEjNK/CIiMaPELyISM0r8IiIx8/8BXjjxX6G+8aoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure()\n",
    "plt.xlabel('voi')\n",
    "plt.ylabel('count')\n",
    "plt.bar(np.arange(len(zero)), zero, color='r')\n",
    "plt.bar(np.arange(len(between)), between, bottom=zero, color='g')\n",
    "plt.bar(np.arange(len(chance)), chance, bottom=zero+between, color='b')\n",
    "plt.bar(np.arange(len(na)), na, bottom=zero+between+chance, color='orange')\n",
    "plt.legend(['zero', 'between', 'chance', 'na'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modify Graph Match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from scipy.optimize import linear_sum_assignment\n",
    "from scipy.optimize import minimize_scalar\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.utils import check_array\n",
    "from sklearn.utils import column_or_1d\n",
    "from graspologic.match import SinkhornKnopp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GraphMatch(BaseEstimator):\n",
    "    \"\"\"\n",
    "    This class solves the Graph Matching Problem and the Quadratic Assignment Problem\n",
    "    (QAP) through an implementation of the Fast Approximate QAP Algorithm (FAQ) (these\n",
    "    two problems are the same up to a sign change) [1].\n",
    "\n",
    "    This algorithm can be thought of as finding an alignment of the vertices of two\n",
    "    graphs which minimizes the number of induced edge disagreements, or, in the case\n",
    "    of weighted graphs, the sum of squared differences of edge weight disagreements.\n",
    "    The option to add seeds (known vertex correspondence between some nodes) is also\n",
    "    available [2].\n",
    "\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    n_init : int, positive (default = 1)\n",
    "        Number of random initializations of the starting permutation matrix that\n",
    "        the FAQ algorithm will undergo. n_init automatically set to 1 if\n",
    "        init_method = 'barycenter'\n",
    "\n",
    "    init_method : string (default = 'barycenter')\n",
    "        The initial position chosen\n",
    "\n",
    "        \"barycenter\" : the non-informative “flat doubly stochastic matrix,”\n",
    "        :math:`J=1*1^T /n` , i.e the barycenter of the feasible region\n",
    "\n",
    "        \"rand\" : some random point near :math:`J, (J+K)/2`, where K is some random doubly\n",
    "        stochastic matrix\n",
    "\n",
    "    max_iter : int, positive (default = 30)\n",
    "        Integer specifying the max number of Franke-Wolfe iterations.\n",
    "        FAQ typically converges with modest number of iterations.\n",
    "\n",
    "    shuffle_input : bool (default = True)\n",
    "        Gives users the option to shuffle the nodes of A matrix to avoid results\n",
    "        from inputs that were already matched.\n",
    "\n",
    "    eps : float (default = 0.1)\n",
    "        A positive, threshold stopping criteria such that FW continues to iterate\n",
    "        while Frobenius norm of :math:`(P_{i}-P_{i+1}) > eps`\n",
    "\n",
    "    gmp : bool (default = True)\n",
    "        Gives users the option to solve QAP rather than the Graph Matching Problem\n",
    "        (GMP). This is accomplished through trivial negation of the objective function.\n",
    "\n",
    "    padding : string (default = 'adopted')\n",
    "        Allows user to specify padding scheme if `A` and `B` are not of equal size.\n",
    "        Say that `A` and `B` have :math:`n_1` and :math:`n_2` nodes, respectively, and\n",
    "        :math:`n_1 < n_2`.\n",
    "\n",
    "        \"adopted\" : matches `A` to the best fitting induced subgraph of `B`. Reduces the\n",
    "        affinity between isolated vertices added to `A` through padding and low-density\n",
    "        subgraphs of `B`.\n",
    "\n",
    "        \"naive\" : matches `A` to the best fitting subgraph of `B`.\n",
    "\n",
    "    Attributes\n",
    "    ----------\n",
    "\n",
    "    perm_inds_ : array, size (n,) where n is the number of vertices in the fitted graphs.\n",
    "        The indices of the optimal permutation (with the fixed seeds given) on the nodes of B,\n",
    "        to best minimize the objective function :math:`f(P) = trace(A^T PBP^T )`.\n",
    "\n",
    "\n",
    "    score_ : float\n",
    "        The objective function value of for the optimal permutation found.\n",
    "\n",
    "    n_iter_ : int\n",
    "        Number of Frank-Wolfe iterations run. If `n_init` > 1, `n_iter_` reflects the number of\n",
    "        iterations performed at the initialization returned.\n",
    "\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] J.T. Vogelstein, J.M. Conroy, V. Lyzinski, L.J. Podrazik, S.G. Kratzer,\n",
    "        E.T. Harley, D.E. Fishkind, R.J. Vogelstein, and C.E. Priebe, “Fast\n",
    "        approximate quadratic programming for graph matching,” PLOS one, vol. 10,\n",
    "        no. 4, p. e0121002, 2015.\n",
    "\n",
    "    .. [2] D. Fishkind, S. Adali, H. Patsolic, L. Meng, D. Singh, V. Lyzinski, C. Priebe,\n",
    "        Seeded graph matching, Pattern Recognit. 87 (2019) 203–215\n",
    "\n",
    "\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        n_init=1,\n",
    "        init_method=\"barycenter\",\n",
    "        max_iter=30,\n",
    "        shuffle_input=True,\n",
    "        eps=0.1,\n",
    "        gmp=True,\n",
    "        padding=\"adopted\",\n",
    "    ):\n",
    "\n",
    "        if type(n_init) is int and n_init > 0:\n",
    "            self.n_init = n_init\n",
    "        else:\n",
    "            msg = '\"n_init\" must be a positive integer'\n",
    "            raise TypeError(msg)\n",
    "        if init_method == \"rand\":\n",
    "            self.init_method = \"rand\"\n",
    "        elif init_method == \"barycenter\":\n",
    "            self.init_method = \"barycenter\"\n",
    "            self.n_init = 1\n",
    "        else:\n",
    "            msg = 'Invalid \"init_method\" parameter string'\n",
    "            raise ValueError(msg)\n",
    "        if max_iter > 0 and type(max_iter) is int:\n",
    "            self.max_iter = max_iter\n",
    "        else:\n",
    "            msg = '\"max_iter\" must be a positive integer'\n",
    "            raise TypeError(msg)\n",
    "        if type(shuffle_input) is bool:\n",
    "            self.shuffle_input = shuffle_input\n",
    "        else:\n",
    "            msg = '\"shuffle_input\" must be a boolean'\n",
    "            raise TypeError(msg)\n",
    "        if eps > 0 and type(eps) is float:\n",
    "            self.eps = eps\n",
    "        else:\n",
    "            msg = '\"eps\" must be a positive float'\n",
    "            raise TypeError(msg)\n",
    "        if type(gmp) is bool:\n",
    "            self.gmp = gmp\n",
    "        else:\n",
    "            msg = '\"gmp\" must be a boolean'\n",
    "            raise TypeError(msg)\n",
    "        if isinstance(padding, str) and padding in {\"adopted\", \"naive\"}:\n",
    "            self.padding = padding\n",
    "        elif isinstance(padding, str):\n",
    "            msg = 'Invalid \"padding\" parameter string'\n",
    "            raise ValueError(msg)\n",
    "        else:\n",
    "            msg = '\"padding\" parameter must be of type string'\n",
    "            raise TypeError(msg)\n",
    "\n",
    "    def fit(self, A, B, seeds_A=[], seeds_B=[]):\n",
    "        \"\"\"\n",
    "        Fits the model with two assigned adjacency matrices\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        A : 2d-array, square, positive\n",
    "            A square adjacency matrix\n",
    "\n",
    "        B : 2d-array, square, positive\n",
    "            A square adjacency matrix\n",
    "\n",
    "        seeds_A : 1d-array, shape (m , 1) where m <= number of nodes (default = [])\n",
    "            An array where each entry is an index of a node in `A`.\n",
    "\n",
    "        seeds_B : 1d-array, shape (m , 1) where m <= number of nodes (default = [])\n",
    "            An array where each entry is an index of a node in `B` The elements of\n",
    "            `seeds_A` and `seeds_B` are vertices which are known to be matched, that is,\n",
    "            `seeds_A[i]` is matched to vertex `seeds_B[i]`.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        self : returns an instance of self\n",
    "        \"\"\"\n",
    "        A = check_array(A, copy=True, ensure_2d=True)\n",
    "        B = check_array(B, copy=True, ensure_2d=True)\n",
    "        seeds_A = column_or_1d(seeds_A)\n",
    "        seeds_B = column_or_1d(seeds_B)\n",
    "\n",
    "        # pads A and B according to section 2.5 of [2]\n",
    "        if A.shape[0] != B.shape[0]:\n",
    "            A, B = _adj_pad(A, B, self.padding)\n",
    "\n",
    "        if A.shape[0] != A.shape[1] or B.shape[0] != B.shape[1]:\n",
    "            msg = \"Adjacency matrix entries must be square\"\n",
    "            raise ValueError(msg)\n",
    "        elif seeds_A.shape[0] != seeds_B.shape[0]:\n",
    "            msg = \"Seed arrays must be of equal size\"\n",
    "            raise ValueError(msg)\n",
    "        elif seeds_A.shape[0] > A.shape[0]:\n",
    "            msg = \"There cannot be more seeds than there are nodes\"\n",
    "            raise ValueError(msg)\n",
    "        elif not (seeds_A >= 0).all() or not (seeds_B >= 0).all():\n",
    "            msg = \"Seed array entries must be greater than or equal to zero\"\n",
    "            raise ValueError(msg)\n",
    "        elif (\n",
    "            not (seeds_A <= (A.shape[0] - 1)).all()\n",
    "            or not (seeds_B <= (A.shape[0] - 1)).all()\n",
    "        ):\n",
    "            msg = \"Seed array entries must be less than or equal to n-1\"\n",
    "            raise ValueError(msg)\n",
    "\n",
    "        n = A.shape[0]  # number of vertices in graphs\n",
    "#         P_tot = np.zeros((n, n))\n",
    "        \n",
    "        n_seeds = seeds_A.shape[0]  # number of seeds\n",
    "        n_unseed = n - n_seeds\n",
    "\n",
    "        score = math.inf\n",
    "        perm_inds = np.zeros(n)\n",
    "\n",
    "        obj_func_scalar = 1\n",
    "        if self.gmp:\n",
    "            obj_func_scalar = -1\n",
    "            score = 0\n",
    "\n",
    "        seeds_B_c = np.setdiff1d(range(n), seeds_B)\n",
    "        if self.shuffle_input:\n",
    "            seeds_B_c = np.random.permutation(seeds_B_c)\n",
    "            # shuffle_input to avoid results from inputs that were already matched\n",
    "\n",
    "        seeds_A_c = np.setdiff1d(range(n), seeds_A)\n",
    "        permutation_A = np.concatenate([seeds_A, seeds_A_c], axis=None).astype(int)\n",
    "        permutation_B = np.concatenate([seeds_B, seeds_B_c], axis=None).astype(int)\n",
    "        A = A[np.ix_(permutation_A, permutation_A)]\n",
    "        B = B[np.ix_(permutation_B, permutation_B)]\n",
    "\n",
    "        # definitions according to Seeded Graph Matching [2].\n",
    "        A11 = A[:n_seeds, :n_seeds]\n",
    "        A12 = A[:n_seeds, n_seeds:]\n",
    "        A21 = A[n_seeds:, :n_seeds]\n",
    "        A22 = A[n_seeds:, n_seeds:]\n",
    "        B11 = B[:n_seeds, :n_seeds]\n",
    "        B12 = B[:n_seeds, n_seeds:]\n",
    "        B21 = B[n_seeds:, :n_seeds]\n",
    "        B22 = B[n_seeds:, n_seeds:]\n",
    "        A11T = np.transpose(A11)\n",
    "        A12T = np.transpose(A12)\n",
    "        A22T = np.transpose(A22)\n",
    "        B21T = np.transpose(B21)\n",
    "        B22T = np.transpose(B22)\n",
    "\n",
    "        for i in range(self.n_init):\n",
    "            # setting initialization matrix\n",
    "            if self.init_method == \"rand\":\n",
    "                sk = SinkhornKnopp()\n",
    "                K = np.random.rand(\n",
    "                    n_unseed, n_unseed\n",
    "                )  # generate a nxn matrix where each entry is a random integer [0,1]\n",
    "                for i in range(10):  # perform 10 iterations of Sinkhorn balancing\n",
    "                    K = sk.fit(K)\n",
    "                J = np.ones((n_unseed, n_unseed)) / float(\n",
    "                    n_unseed\n",
    "                )  # initialize J, a doubly stochastic barycenter\n",
    "                P = (K + J) / 2\n",
    "            elif self.init_method == \"barycenter\":\n",
    "                P = np.ones((n_unseed, n_unseed)) / float(n_unseed)\n",
    "\n",
    "            const_sum = A21 @ np.transpose(B21) + np.transpose(A12) @ B12\n",
    "            grad_P = math.inf  # gradient of P\n",
    "            n_iter = 0  # number of FW iterations\n",
    "\n",
    "            # OPTIMIZATION WHILE LOOP BEGINS\n",
    "            while grad_P > self.eps and n_iter < self.max_iter:\n",
    "\n",
    "                delta_f = (\n",
    "                    const_sum + A22 @ P @ B22T + A22T @ P @ B22\n",
    "                )  # computing the gradient of f(P) = -tr(APB^tP^t)\n",
    "                rows, cols = linear_sum_assignment(\n",
    "                    obj_func_scalar * delta_f\n",
    "                )  # run hungarian algorithm on gradient(f(P))\n",
    "                Q = np.zeros((n_unseed, n_unseed))\n",
    "                Q[rows, cols] = 1  # initialize search direction matrix Q\n",
    "\n",
    "                def f(x):  # computing the original optimization function\n",
    "                    return obj_func_scalar * (\n",
    "                        np.trace(A11T @ B11)\n",
    "                        + np.trace(np.transpose(x * P + (1 - x) * Q) @ A21 @ B21T)\n",
    "                        + np.trace(np.transpose(x * P + (1 - x) * Q) @ A12T @ B12)\n",
    "                        + np.trace(\n",
    "                            A22T\n",
    "                            @ (x * P + (1 - x) * Q)\n",
    "                            @ B22\n",
    "                            @ np.transpose(x * P + (1 - x) * Q)\n",
    "                        )\n",
    "                    )\n",
    "\n",
    "                alpha = minimize_scalar(\n",
    "                    f, bounds=(0, 1), method=\"bounded\"\n",
    "                ).x  # computing the step size\n",
    "                P_i1 = alpha * P + (1 - alpha) * Q  # Update P\n",
    "                grad_P = np.linalg.norm(P - P_i1)\n",
    "                P = P_i1\n",
    "                n_iter += 1\n",
    "            # end of FW optimization loop\n",
    "\n",
    "            row, col = linear_sum_assignment(\n",
    "                -P\n",
    "            )  # Project onto the set of permutation matrices\n",
    "            perm_inds_new = np.concatenate(\n",
    "                (np.arange(n_seeds), np.array([x + n_seeds for x in col]))\n",
    "            )\n",
    "\n",
    "            score_new = np.trace(\n",
    "                np.transpose(A) @ B[np.ix_(perm_inds_new, perm_inds_new)]\n",
    "            )  # computing objective function value\n",
    "            \n",
    "#             m_foo = len(seeds_A)-1\n",
    "#             n_foo = n - m_foo\n",
    "#             prob_mat = np.diag(np.ones((n_foo)))\n",
    "#             prob_mat = np.concatenate((np.concatenate((np.diag(np.ones((m_foo))), np.zeros((m_foo, n_foo))), axis=1), np.concatenate((np.zeros((n_foo, m_foo)), prob_mat[perm_inds_new]), axis=1)), axis=0)\n",
    "#             P_tot+=prob_mat\n",
    "\n",
    "            if obj_func_scalar * score_new < obj_func_scalar * score:  # minimizing\n",
    "                score = score_new\n",
    "                perm_inds = np.zeros(n, dtype=int)\n",
    "                perm_inds[permutation_A] = permutation_B[perm_inds_new]\n",
    "                best_n_iter = n_iter\n",
    "\n",
    "        permutation_A_unshuffle = _unshuffle(permutation_A, n)\n",
    "        A = A[np.ix_(permutation_A_unshuffle, permutation_A_unshuffle)]\n",
    "        permutation_B_unshuffle = _unshuffle(permutation_B, n)\n",
    "        B = B[np.ix_(permutation_B_unshuffle, permutation_B_unshuffle)]\n",
    "        score = np.trace(np.transpose(A) @ B[np.ix_(perm_inds, perm_inds)])\n",
    "\n",
    "#         self.P_final = P_tot/self.n_init\n",
    "        self.perm_inds_ = perm_inds  # permutation indices\n",
    "        self.score_ = score  # objective function value\n",
    "        self.n_iter_ = best_n_iter\n",
    "        return self\n",
    "\n",
    "\n",
    "    def fit_predict(self, A, B, seeds_A=[], seeds_B=[]):\n",
    "        \"\"\"\n",
    "        Fits the model with two assigned adjacency matrices, returning optimal\n",
    "        permutation indices\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        A : 2d-array, square, positive\n",
    "            A square, positive adjacency matrix\n",
    "\n",
    "        B : 2d-array, square, positive\n",
    "            A square, positive adjacency matrix\n",
    "\n",
    "        seeds_A : 1d-array, shape (m , 1) where m <= number of nodes (default = [])\n",
    "            An array where each entry is an index of a node in `A`.\n",
    "\n",
    "        seeds_B : 1d-array, shape (m , 1) where m <= number of nodes (default = [])\n",
    "            An array where each entry is an index of a node in `B` The elements of\n",
    "            `seeds_A` and `seeds_B` are vertices which are known to be matched, that is,\n",
    "            `seeds_A[i]` is matched to vertex `seeds_B[i]`.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        perm_inds_ : 1-d array, some shuffling of [0, n_vert)\n",
    "            The optimal permutation indices to minimize the objective function\n",
    "        \"\"\"\n",
    "        self.fit(A, B, seeds_A, seeds_B)\n",
    "        return self.perm_inds_\n",
    "\n",
    "\n",
    "\n",
    "def _adj_pad(A, B, method):\n",
    "    def pad(X, n):\n",
    "        X_pad = np.zeros((n[1], n[1]))\n",
    "        X_pad[: n[0], : n[0]] = X\n",
    "        return X_pad\n",
    "\n",
    "    A_n = A.shape[0]\n",
    "    B_n = B.shape[0]\n",
    "    n = np.sort([A_n, B_n])\n",
    "    if method == \"adopted\":\n",
    "        A = 2 * A - np.ones((A_n, A_n))\n",
    "        B = 2 * B - np.ones((B_n, B_n))\n",
    "\n",
    "    if A.shape[0] == n[0]:\n",
    "        A = pad(A, n)\n",
    "    else:\n",
    "        B = pad(B, n)\n",
    "\n",
    "    return A, B\n",
    "\n",
    "\n",
    "def _unshuffle(array, n):\n",
    "    unshuffle = np.array(range(n))\n",
    "    unshuffle[array] = np.array(range(n))\n",
    "    return unshuffle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1,2,3],[4,5,6],[7,8,9]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4, 5, 6],\n",
       "       [1, 2, 3],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[[1, 0, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}