class_04_hits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HITS & ECI"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![hits.png](hits.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:17.465741Z",
     "start_time": "2019-07-04T14:30:17.456668Z"
    }
   },
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import numpy as np\n",
    "import pickle \n",
    "import pandas as pd\n",
    "\n",
    "fr = open('MatrixD.pkl','rb')  \n",
    "dat = pickle.load(fr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:18.115860Z",
     "start_time": "2019-07-04T14:30:18.110681Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(102, 242)"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = dat[1995]\n",
    "len(df[0]), len(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MATRIX TO EDGELIST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:19.427176Z",
     "start_time": "2019-07-04T14:30:19.377452Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24684"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat = []\n",
    "for i in range(len(df)):\n",
    "    for j in range(len(df[i])):\n",
    "        #if df[i][j]:\n",
    "        dat.append(['row'+str(i), 'col' + str(j), df[i][j]])\n",
    "\n",
    "len(dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:21.474200Z",
     "start_time": "2019-07-04T14:30:21.399558Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Name: \\nType: DiGraph\\nNumber of nodes: 344\\nNumber of edges: 24684\\nAverage in degree:  71.7558\\nAverage out degree:  71.7558'"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = nx.DiGraph()\n",
    "for i in dat:\n",
    "    G.add_edge( i[1], i[0], weight = i[2])\n",
    "\n",
    "nx.info(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:26.863726Z",
     "start_time": "2019-07-04T14:30:25.386439Z"
    }
   },
   "outputs": [],
   "source": [
    "h, a = nx.hits(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:27.741296Z",
     "start_time": "2019-07-04T14:30:27.727046Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "102"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "country_ids = [i for i in list(G.nodes) if 'col' in i]\n",
    "df1 = [(i, h[i]) for i in country_ids]\n",
    "df1 = pd.DataFrame(df1, columns = ['id', 'hub'])\n",
    "df1['ids'] = [int(i.split('col')[1]) for i in df1['id'].tolist()]\n",
    "df1 = df1.sort_values(by = 'ids', ascending = True)\n",
    "len(df1['hub'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:30:28.677689Z",
     "start_time": "2019-07-04T14:30:28.664020Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "242"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "country_ids = [i for i in list(G.nodes) if 'row' in i]\n",
    "df2 = [(i, a[i]) for i in country_ids]\n",
    "df2 = pd.DataFrame(df2, columns = ['id', 'aut'])\n",
    "df2['ids'] = [int(i.split('row')[1]) for i in df2['id'].tolist()]\n",
    "df2 = df2.sort_values(by = 'ids', ascending = True)\n",
    "len(df2['aut'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ECI & PCI\n",
    "\n",
    "https://nbviewer.jupyter.org/github/hashc/Skills/blob/master/ecomplexity.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:15:29.149291Z",
     "start_time": "2019-07-04T14:15:28.921966Z"
    }
   },
   "outputs": [],
   "source": [
    "def RCA(Xcp):\n",
    "    \"\"\"\n",
    "    Xcp is a numpy.matrix type\n",
    "    --|steel|oil|chip|bean|clothes\n",
    "    --|--|--|--|--|--\n",
    "    USA|1|1|1|1|0\n",
    "    China|1|0|0|1|1\n",
    "    Viet Nam|0|0|0|1|1\n",
    "    \"\"\"\n",
    "    if type(Xcp) is not np.matrix:\n",
    "        if type(Xcp) is np.array or list:\n",
    "            Xcp = np.matrix(Xcp)\n",
    "        else:\n",
    "            raise ValueError('Xcp must be matrix')\n",
    "    else:\n",
    "        pass\n",
    "    B = Xcp.sum(1)*Xcp.sum(0) \n",
    "    Sum = Xcp.sum()\n",
    "    return Xcp*Sum/B\n",
    "    \n",
    "def fliterRCA(R):\n",
    "    M = R>1\n",
    "    return M.astype(float)\n",
    "\n",
    "\n",
    "def Get_eci_pci(M):\n",
    "    d = M.sum(1).T.tolist()[0]\n",
    "    u = M.sum(0).tolist()[0]\n",
    "    D = np.diag([1.0/i if i>0 else 0. for i in d])\n",
    "    U = np.diag([1.0/i if i>0 else 0. for i in u])\n",
    "    mcp1 = D * M\n",
    "    mcp2 = M * U\n",
    "    \n",
    "    Mcc = mcp1 * mcp2.T\n",
    "    Mpp = mcp2.T * mcp1\n",
    "    eigvals, eigvecs = np.linalg.eig(Mpp)\n",
    "    eigvecs = np.real(eigvecs)\n",
    "    # Get eigenvector corresponding to second largest eigenvalue\n",
    "    eig_index = eigvals.argsort()[-2]\n",
    "    kp = eigvecs[:, eig_index]\n",
    "    kc = mcp1 @ kp\n",
    "    s1 = np.sign(np.corrcoef(M.sum(1).reshape(-1), kc.reshape(-1))[0, 1])\n",
    "    eci = s1 * kc\n",
    "    pci = s1 * kp  \n",
    "    return eci.T.tolist()[0],pci.T.tolist()[0]\n",
    "\n",
    "def Get_eci_pci_sparse(M):\n",
    "    d = M.sum(1).T.tolist()[0]\n",
    "    u = M.sum(0).tolist()[0]\n",
    "    nd,nu = len(d),len(u)\n",
    "    D1 = scipy.sparse.csc_matrix(([1.0/i if i>0 else 0. for i in d], (range(nd), range(nd))),shape=(nd,nd))\n",
    "    U1 = scipy.sparse.csc_matrix(([1.0/i if i>0 else 0. for i in u], (range(nu), range(nu))),shape=(nu,nu))\n",
    "    M = sparse.csc_matrix(M)\n",
    "    mcp1 = D1 * M\n",
    "    mcp2 = M * U1\n",
    "    Mcc = mcp1 @ mcp2.T\n",
    "    Mpp = mcp2.T @ mcp1\n",
    "    A=scipy.sparse.csc_matrix(Mpp)\n",
    "    eigvals, eigvecs = scipy.sparse.linalg.eigs(A, k=2)\n",
    "    eigvecs = np.real(eigvecs)\n",
    "    # Get eigenvector corresponding to second largest eigenvalue\n",
    "    eig_index = eigvals.argsort()[-2]\n",
    "    kp = eigvecs[:, eig_index]\n",
    "    kc = mcp1 @ kp\n",
    "    s1 = np.sign(np.corrcoef(M.sum(1).reshape(-1), kc.reshape(-1))[0, 1])\n",
    "    eci = s1 * kc\n",
    "    pci = s1 * kp\n",
    "    return eci,pci\n",
    "\n",
    "def Get_z_score(x):\n",
    "    x=np.array(x)\n",
    "    std=np.std(x)\n",
    "    mean=np.mean(x)\n",
    "    return (x-mean)/std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:16:54.648823Z",
     "start_time": "2019-07-04T14:16:54.619472Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/datalab/Applications/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in true_divide\n",
      "/Users/datalab/Applications/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:22: RuntimeWarning: invalid value encountered in greater\n"
     ]
    }
   ],
   "source": [
    "c = dat[1995] # MatrixD[1995]\n",
    "R = RCA(c)\n",
    "M = fliterRCA(R)\n",
    "eci,pci = Get_eci_pci(np.matrix(M))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:25:40.707438Z",
     "start_time": "2019-07-04T14:25:40.703158Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(242, 102, 102)"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(eci), len(pci), len(hub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:31:05.960583Z",
     "start_time": "2019-07-04T14:31:05.955896Z"
    }
   },
   "outputs": [],
   "source": [
    "pci_z = Get_z_score(pci)\n",
    "eci_z = Get_z_score(eci)\n",
    "\n",
    "hub = df1['hub'].tolist()\n",
    "aut = df2['aut'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:20:53.306082Z",
     "start_time": "2019-07-04T14:20:53.300928Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:29:47.918047Z",
     "start_time": "2019-07-04T14:29:47.270811Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAF9ZJREFUeJzt3X+MLWddx/HP9+69q2xbIN171QjdXQjEiJUIbDTGgEWIkoZQRVHjtpZIvHIbjSb+Ac0m8ofZSEI0KRDarFA03COGIApCCT8SCQjRuJcAKRYIEPbaQGDbGihcQ2n79Y/Zwz27d2bOMzPP/HrO+5VMds+cc2aeZ8/ud579Pj/G3F0AgHSc6LsAAIC4COwAkBgCOwAkhsAOAIkhsANAYgjsAJAYAjsAJIbADgCJIbADQGII7ACQmJN9nPT06dO+sbHRx6kBYLQuXLjwgLufmfe6XgL7xsaG9vb2+jg1AIyWme2HvI5UDAAkhsAOAIkhsANAYgjsAJAYAjsAdGEykTY2pBMnsq+TSWun6mVUDAAslMlEOntWunQpe7y/nz2WpK2t6KejxQ4AbdvevhzUpy5dyva3gMAOAG27eLHa/oYI7ADQtrW1avsbIrADQNt2dqSVlaP7Vlay/S0gsANA27a2pN1daX1dMsu+7u620nEqMSoGALqxtdVaID+OFjsAJIbADgCJIbADQGII7ACQGAI7ACSGwA4AiWkc2M3sOjP7NzO7z8w+b2Z/GqNgAIB6Yoxjf1TSn7v7p83sGkkXzOwj7v7fEY4NAKiocYvd3b/h7p8+/P5hSfdJekrT4wIA6omaYzezDUnPkfSfMY8LAAgXLbCb2dWS/knSn7n7d3KeP2tme2a2d3BwEOu0AIBjogR2MzulLKhP3P09ea9x911333T3zTNnzsQ4LQAgR4xRMSbpbZLuc/e/aV4kAEATMVrsvyTpFkm/YmafOdxujHBcAEANjYc7uvu/S7IIZQEARMDMUwDpm0ykjQ3pxIns62TSd4laxY02AKRtMpHOnpUuXcoe7+9nj6XObnzRNVrsAOoZSyt4e/tyUJ+6dCnbnygCO4Dqpq3g/X3J/XIreIjB/eLFavurGOjFjcAOoLqyVvDQgt3aWrX9oQZ8cSOwA6iuqLU7DW5DCnY7O9LKytF9KyvZ/iYGnOIhsAOorqi1u7Q0vGC3tSXt7krr65JZ9nV3t3nHaZspnoYI7ACqK2oFP/ZY/uv399svU5mtLelrX5Mefzz7GmM0TFspnggI7ACqK2oFLy3lv75o/5i1leKJwNy985Nubm763t5e5+cF0DIrmYTeQ6xp3WSSpZkuXsxa6js7rY6NN7ML7r4573W02IFF0GSkyvH33nZb8bHW1/OPUbS/bW2P0GkjxRODu3e+Pe95z3MAHTl/3n1lxT1rM2fbykq2v857j2+zx2pyrqp1Wl93N8u+5h2/q7J0SNKeB8RYAjsQIiSQDNX6en5AXl+v/96yY7X9swoN2E3qHauckX8OoYGdHDswz/G1RqSskyzGkLkunDiRn982y1IIdd5b51ixbGzkj7JZX8/SIVNN6t1US78z5NiBWAY8ESVIk2F5oUP3uhziFzp+vM/hiD3/zhDYgXkGPBElSJNheXnvPa7rIX6hAbvP4Yg9/84Q2IF5BjwRpdDsaJDtbenWW+vNvMwbr37uXPxZnFWEBuy2ZpyG6Pt3JiQRH3uj8xSjMrbRFWMrbx1D78xu6TMQo2KAiIYWSMrK0/doEGQYFQMg2LwRF32MBul4BuaiYlQMkKp5Iy66zu82WZd8aGu3J4LADozNvBEXXY8GqTu0b8A3qhg7AjswNvNa5F2PBqk7tG/s8wMGjMAOjE1Zi3ya2rjllmz/O97R/uJUdVM/Y58fMGAEdmAsZoP2E54gra4ebZFL/aQ26qZ+yi4I5N6bCRk6M2+TdLekb0m6N+T1DHcEKgoZFx1rmGOdYXp135NXp3Pn0h+HX5O6HMcu6QWSnktgB1oSErTN8l9jFn6eric35V0QGIdfKDSwRxvHbmYbkt7v7tfPey3j2IGKQsamh656WCbGMZrqc1XGgWMcO5CSkA7KGMMch9Ch2fc6KwnoLLCb2Vkz2zOzvYODg65OC6QhJGjHGOY4hKA64JtEj0ZIviZkk7QhcuxAe7pYr6Zqjr2tMg1tbZ6BUNeLgBHYgUSEBtWm91INPQcB/odCA3uUzlMze6ekGySdlvRNSa9z97cVvZ7OUyABdTtaQ28bN/ZbErYgtPOU1R0B1GNW/FxZXAm9IAxhhM7AMCoGQLuWlqrtnwodeTOEETojRWAHFkXsafqPPVZt/1ToyJshjNAZKQI7kLrJRDp9Wrr55rjryKyvV9s/FTqckWGPtRHYgZRNOyAffPDK55oskTuZSN/97pX7QwJv6Hj7Pm9GPXJ0ngIpK+qAnKozTT9vtIqUrTZ5xx0E3hbReQpgfkdjnSVy826QIUlXX91vUGep3x862XcBALRoba24xb6yIt1449HW9zT3LhUH6SGOVjn+X0RIPRJGix1IWV4HpJSlTXZ3pXvuqX57uiGOVuE2e0cQ2IGU5XVAnj8vPfBA9lyd1vcQR6sM8b+IHhHYgdRtbWUzNR9//Mr7n157bf57ivZPjze00SpD/C+iRwR2IFVNOhO///3y95ZdLELKctttcTs6h/hfRJ9CVgqLvbG6IxZGX6sThq68mHcLurytye3x8soS8/iz50l8JUh1fWu8KhjHjoXQ5+qEIQtoTSbSLbeUL9hV9N4YZYl1/AXC6o5A3/pcnbDJPVKL1L3naFFZYh1/gTBBCehbnyM1QjoTq5Zj+t6qufvQDswF7ehsA4EdaEufIzVCOhOLyrG6WvzeaXqpymJiRWPpy8qGZkIS8bE3Ok+xEJrcOi7W+cs6E8vKV/Te9fX8zs/19WplOXcu+Y7ONojOU2AAJpNs9uPFi1kLeWdnWFPcq5YvJHeP1tB5CiA+blfXKzpPAcQ3lIlArORYisAO4KiyoDmE5QTqdOAuGAI70IaxtijnBc0h9BmwkuNcBHYgtpAWZZeBv8q5yoLmUFrKrOQ4X8jQmdgbwx2RtHlDArscBln1XGb5ZZ8OSyxa66XLIYt1h1wmQIHDHWmxA7HNa1F2mUqoeq6ySVVlLeIuW+9D6cAdsCiB3cxeYmZfNLMvm9lrYxwTqK3vNMe8GaddphKqnqssaM6bMdtVnnsIHbhDF9KsL9skLUn6iqSnS1qW9FlJzyp7D6kYtKatNEfeTMy8c5m5v+hF5WXoMpVQ51xFs05Dlt81i18H/JACUzExAvsvSvrQzOPbJd1e9h4CO1rTRtAsulisrhYHt7Ip80POsYccb16uHa3pMrD/lqS3zjy+RdKby95DYEdryjr/6ioLZHUDXJc3hWjjXHUuGAtwI4y2dRnYX5ET2N+U87qzkvYk7a2trXXwI8BCaqPFXnSxGHpKou1AWuX4fS+IlojQwB6j8/R+SdfNPH6qpK/n5PJ33X3T3TfPnDkT4bRAjjZGTJQtb2sW9p7YHbrzjlc25jxWWY7f91QqPi6TiroVEv3LNkknJX1V0tN0ufP0Z8reQyoGrYrdUi1rbZ47d2WLfvq4rJM19j1Ejx+v6D+X1dX2OpfLjttGimwBqatUTHYu3SjpS8pGx2zPez2BHaNTdrGY7VA8HsDKOlnrpodC0k1V00dNOz3nlWmBJxXFFBrYWbYXiKXve4jOHu/0aenBB9svS2iZ+ryxd0JYthfoWt17iFY1bwLUZCI9/PCVz586lfULxCxLaJmYVNQpAjsQIqTDsc49ROuY10G8vS098siV73viE6U77mhnOn5Ip/XxztZYQX2sK2m2KSRfE3sjx45RCe38rHMP0SZlKjrevI7KGGXJO8a5c+5LS9m5lpayx7HOV1aOBRpGqS47T6tuBHaMSpWOvyFMwmm7ozIvmJ465b68XHxBaSvwLlinbGhgp/MUmGdsN3Buu6Oyaidxnlj3SB3bZ9MQnadALPM6Boem7Y7KGKtQxlrJcmyfTUcI7MA8Q1v/O6SzsEpHZdXOxxhBM1bgHdpnMxQh+ZrYGzl2jM4QcufTcnQ9izXkPadO5ee687bYnZtD+Ww6IHLsQIKK8tt1c9Z1j5d3U+tPflK6667inLf75XH0Dz3U382wRyw0x05gB8Ykdmdh7ONNA/7+vrS0JD32WHaRmKZGmH3aCIEdSNFQWux1dHmuRDEqBkhR7M7CnR1pefnovuXlsONV7XTt8l6vC47ADoxJG0MZH300/3FZ4C5b770IQxM7QyoGWGRFq0BedVUWsIvy4XXSKqzw2BipGADzFS3t+73vld/xqE5ahRUeO0NgB1Jz223SyZNZ8Dx5MnscyzRwF6VPrr22/P1trfCIIwjsWExjXOo15D6nV18t3XlnNsxQyr7eeWdxcC9an/1EQWiYBvSdnWx99+MefngcP8vUhcxiir0x8xS9GuNSr/PKnPf87La0VHzc46syLi9nS+7O+xnFvuUf5hLL9gIFxrjUa1GZl5ayqfTTddDLtiJFU/LnTdXnBtWdCw3sjIrB4hnjUq9FZQ61tHTlsMammHDUOUbFAEXGOJ66adnOng1/bWj/AysrDhaBHYtnaAEpJJDmlTnEiRPSuXPSW94SXpbQiUcMXxyukHxN7I0cO3o3lKVeQzpFp+VcXc22spz6NOdet05j7H9YICLHDoxAWZ56Z6d4pqbUzizOslz+tEy0yHtDjh0Yg7IZnNvbxbM/20qDlOXyQ9aDyZOXahrjPIIxCWnWF22SXiHp85Iel7QZ+j5SMcChstRHrOGEVdJOeePam6Rl8lJNy8tX3nFp6PMIBkKBqZimLfZ7Jb1c0scbHgdYTGUduTFG71RdhXFrS7rmmvJjVllmN++/jkcekX7wg6P7ZtehQWONAru73+fuX4xVGGDhlKVUYozeKUvnFHnoofJjzl5Y5qVUqlwE9vezn8Hp06Rmmgpp1s/bJH1MpGKA+JqO3qmTzilKD4UsY3A8pVJ2rLJteZnUTA7FSsWY2UfN7N6c7aYqFxAzO2tme2a2d3BwUOsiBCycpqshzkvn5LW4i8bMX3XV0Q7akP8G8o61vJy/gNisRx4hNdPA3MDu7i929+tztvdWOZG777r7prtvnjlzpn6JgXn6HnHR9/lnlaVzivLvknTrrVlaZNbxYZAha7LnpZruvlt6+9uz78twy7z6Qpr18zaRisFQ9L1yY9/nLypTXjqnbERO0XOrq5ePG2MyU1mqpq1JUUOZnFaDuljdUdJvSLpf0vclfVPSh0LeR2BHa/qeOdn3+asoy78XPSdVy7HPc/78lUMf28yxD/HCW0Engb3uRmBHa/peSrbv81dRp8V+/CIVo/V7/vzRtd1XV9sLtGO68OYIDewsKYC09L2UbN/nr6Ls5tKSdPPN+e8b8vLG84xxyeYZLCmAxdT3yo19n7+KsjH0W1vFt80b8vLG84xxyeYaCOxIS99LyfZ9/qrKhlPeccd4LlKhxnThbYDAjvQ0Hfs99vPHUnSRkoYznLOqsV14ayLHDiBcWV4+seA4ROTYAcRXZ+0ZdI7ADiBcyGxT9I7ADiDcEEeVDGkJh4EgsAMIN7RRJVXXmw895sgvFAR2AOGGNqokds6/jQtFDwjsQOpit0CHNJwzds4/kc5hAjuQskRaoIVi5/wT6RwmsAMpS6QFWih2zn+IncM1ENiBsaiTUkmkBVoods5/aJ3DNRHYgTGom1JJpAVaKmbOP5FlFAjswBjUTakk0gLt1PELhTS6fgoCOzAGdVMqQxueOEYj7Kc42XcBAARYW8u/gUdISmW6vjrqGWE/BS12YAxIqfRnhP0UBHZgDEip9GeEF1VSMcBYkFLpx/Rnvr2dpV/W1rKgPuDPgsAOAPOM7KJKKgYAEkNgB4DEENgBIDGNAruZvcHMvmBmnzOzfzazJ8cqGACgnqYt9o9Iut7dny3pS5Jub14kAEATjQK7u3/Y3R89fPgfkp7avEgAgCZi5tj/QNIHi540s7NmtmdmewcHBxFPCwCYNXccu5l9VNJP5Dy17e7vPXzNtqRHJRUud+buu5J2JWlzc9NrlRYAMNfcwO7uLy573sxulfRSSS9ydwI2APSs0cxTM3uJpNdI+mV3vzTv9QCA9jXNsb9Z0jWSPmJmnzGzuyKUCQDQQKMWu7s/I1ZBAABxMPMUABJDYAeAxBDYASAxBHYASAyBHQASQ2AHgMQQ2AEgMQR2AEgMgR0AEkNgB4DEENgBIDEEdgBIDIEdABJDYAeAxBDYASAxBHYASAyBHQASQ2AHgMQQ2AEgMQR2AEgMgR0AEkNgB4DEENgBIDEEdgBITKPAbmZ/aWafM7PPmNmHzewnYxUMAFBP0xb7G9z92e7+c5LeL+kvIpQJANBAo8Du7t+ZeXiVJG9WHABAUyebHsDMdiT9vqRvS3ph4xIBABqZ22I3s4+a2b05202S5O7b7n6dpImkPy45zlkz2zOzvYODg3g1AAAcYe5xsidmti7pA+5+/bzXbm5u+t7eXpTzAsCiMLML7r4573VNR8U8c+bhyyR9ocnxAADNNc2xv97MfkrS45L2Jb26eZEAAE00Cuzu/puxCgIAiIOZpwCQGAI7ACSGwA4AiSGwA0BiCOwAkBgCOwAkZjyBfTKRNjakEyeyr5NJ3yUCgEFqvAhYJyYT6exZ6dKl7PH+fvZYkra2+isXAAzQOFrs29uXg/rUpUvZfgDAEeMI7BcvVtsPAAtsHIF9ba3afgBYYOMI7Ds70srK0X0rK9l+AMAR4wjsW1vS7q60vi6ZZV93d+k4BYAc4xgVI2VBnEAOAHONo8UOAAhGYAeAxBDYASAxBHYASAyBHQASY+7e/UnNDpTd/FqSniTp2zNPzz7O+35232lJD9QowvFzVnlN3v66dahb/rLyhbxmXh2K6pP3mjbrUPZ82c/8+ON53/dRhxi/R7Pfj/1vQRpmHULq0+Xf85Pd/czco7l7r5uk3aLHed8f27cX45xVXpO3v24d6pa/7ToU1aegLq3Voez5sp95yGfQdx1i/B7FqMNQ/haGWoeQ+vT995y3DSEV868lj/O+P/76GOes8pq8/anVoag+Za+pY94xyp4v+5kffxzyfV116xDj9yjk/PPwt1C+L6Q+fdfhCr2kYmIxsz133+y7HHWNvfwSdRgK6tC/IZV/CC32Jnb7LkBDYy+/RB2Ggjr0bzDlH3WLHQBwpbG32AEAxxDYASAxBHYASEySgd3Mnm9md5nZW83sU32Xpw4zO2FmO2b2JjO7te/y1GFmN5jZJw4/ixv6Lk9dZnaVmV0ws5f2XZaqzOynD3/+7zazc32Xpw4z+3Uz+1sze6+Z/Wrf5anDzJ5uZm8zs3d3cb7BBXYzu9vMvmVm9x7b/xIz+6KZfdnMXlt2DHf/hLu/WtL7Jf19m+XNE6MOkm6S9BRJP5B0f1tlLRKpDi7pu5J+VOOtgyS9RtK72illsUh/C/cd/i38tqTOh+JFqsO/uPsfSnqlpN9psbi5ItXhq+7+qnZLevSEg9okvUDScyXdO7NvSdJXJD1d0rKkz0p6lqSfVRa8Z7cfm3nfuyQ9cYx1kPRaSX90+N53j7QOJw7f9+OSJiOtw4sl/a6yoPLSsZX/8D0vk/QpSb83xs9g5n1/Lem5I69DJ3/Lg7uDkrt/3Mw2ju3+eUlfdvevSpKZ/aOkm9z9ryTl/ntsZmuSvu3u32mxuLli1MHM7pf0yOHDx9orbb5Yn8Oh/5X0I22Us0ykz+GFkq5S9kf7f2Z2j7s/3mrBD8X6DNz9fZLeZ2YfkPQP7ZU499wxPgOT9HpJH3T3T7db4itF/lvoxOACe4GnSPqfmcf3S/qFOe95laS3t1ai6qrW4T2S3mRmz5f08TYLVkGlOpjZyyX9mqQnS3pzu0ULVqkO7r4tSWb2SkkPdBXUS1T9DG6Q9HJlF9Z7Wi1ZuKp/C3+i7D+nJ5nZM9z9rjYLF6jq57AqaUfSc8zs9sMLQGvGEtgtZ1/pzCp3f11LZamrUh3c/ZKyi9OQVK3De5RdoIak8u+SJLn738UvSi1VP4OPSfpYW4WpqWod3ijpje0Vp5aqdXhQ0qvbK85Rg+s8LXC/pOtmHj9V0td7Kktd1GEYxl6HsZdfog6tG0tg/y9JzzSzp5nZsrLOrPf1XKaqqMMwjL0OYy+/RB3a13UPc0AP9DslfUOXh/m96nD/jZK+pKwnervvclIH6kD5qcNQ68AiYACQmLGkYgAAgQjsAJAYAjsAJIbADgCJIbADQGII7ACQGAI7ACSGwA4AiSGwA0Bi/h+BBWHaMeKVVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hub, pci_z, 'ro')\n",
    "plt.xscale('log');\n",
    "#plt.yscale('log');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:31:15.534518Z",
     "start_time": "2019-07-04T14:31:14.824009Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAH3lJREFUeJztnW2MJllVx/+ne3uEHt6yvaPZsEw3BCUgEl46GkJAwhLdrASQICINGXTjuENi1g9E2fQH9UPHGCMhwIdN86Ij3UEIorwqL4ENEATtQSCLCwTI9LiRuA0T1E0TFnaOH6ofu6am7q17q+6tuvc+/19S6X7qqafq1K26/3vq3FP3iqqCEEJIOSxMbQAhhJCwUNgJIaQwKOyEEFIYFHZCCCkMCjshhBQGhZ0QQgqDwk4IIYVBYSeEkMKgsBNCSGFQ2AkhpDCum+KgN9xwg66trU1xaEIIyZYLFy58T1VPdW03ibCvra1hb29vikMTQki2iMi+y3YMxRBCSGFQ2AkhpDAo7IQQUhgUdkIIKYxgwi4iiyLybyLy4VD7JIWwuwusrQELC9Xf3d2pLSKkaEJmxdwB4F4Ajwq4T5I7u7vA2bPA4WH1eX+/+gwAGxvT2UVIwQTx2EXkJgC/BuDtIfZHCmJz81jUZxweVusJIVEIFYp5E4A/BHAl0P5IKVy65LeeEDKYwcIuIi8CcL+qXujY7qyI7InI3sHBwdDDklw4fdpvPSFkMCE89ucAeLGIXATwtwBeICI7zY1UdVtV11V1/dSpzjdiSSlsbQHLy1evW16u1hNCojBY2FX1TlW9SVXXALwSwKdU9dWDLSNlsLEBbG8Dq6uASPV3e5sdp4REZJKxYsicsbFBISdkRIIKu6reDeDukPskhBDiB988JYSQwqCwE0JIYVDYCSGkMCjshBBSGBR2QggpDAo7IYQUBoWdEEIKg8JOCCGFQWEnhJDCoLATQkhhUNgJIaQwKOyEEFIYFHZCCCkMCjshhBQGhZ0QQgqDwk4IIYVBYSeEkNDs7gJra8DCQvV3d3fUw3NqPEIICcnuLnD2LHB4WH3e368+A6NNEUmPnRBCQrK5eSzqMw4Pq/UjQWEnhJCQXLrktz4CFHZCCAnJ6dN+6yNAYSeEkJBsbQHLy1evW16u1o/EYGEXkYeJyL+IyFdE5Gsi8qchDCOEkCzZ2AC2t4HVVUCk+ru9PVrHKRAmK+ZHAF6gqg+IyBKAz4nIP6rqFwLsmxBC8mNjY1QhbzJY2FVVATxw9HHpaNGh+yWEENKPIDF2EVkUkS8DuB/AJ1T1iyH2SwghxJ8gwq6qD6nq0wHcBOAXReSpzW1E5KyI7InI3sHBQYjDEkIIaSFoVoyq/gDA3QBuafluW1XXVXX91KlTIQ9LCCGkRoismFMi8pij/x8O4IUAvj50v4QQQvoRIivmRgDnRWQRVUPxXlX9cID9EkII6UGIrJivAnhGAFsIIYQEgG+eEkJIYVDYCSGkMCjshBBSGBR2QggpDAo7IYQUBoWdEEIKg8JOCCGFQWEnhJDCoLATQkhhUNinZncXWFsDFhaqv7u7U1tESFmMWccSqc8hxoohfdndBc6eBQ4Pq8/7+9VnYNLZVwgphjHrWEL1WaoJkMZlfX1d9/b2Rj9ucqytVRe/yeoqcPHi2NYQUh5j1rERjiUiF1R1vWs7hmKm5NIlv/WEED/GrGMJ1WcK+5ScPu23nhDix5h1LKH6TGGfkq0tYHn56nXLy9V6QshwxqxjCdVnCvuUbGwA29tVDE6k+ru9zY5TQkIxZh1LqD6z85QQQjKBnaeEkLxJJCc8R5jHTghJj4RywnOEHjshJD02N49FfcbhYbWedEJhjwUfIwnpT0I54TlCYY/B7DFyfx9QPX6MpLgT4kZCOeE5QmGPAR8jCRmGa044n4xbGSzsIvI4Efm0iNwrIl8TkTtCGJY1fIwkZBguOeF8MjYyOI9dRG4EcKOqfklEHgngAoCXquq/m35TfB47B/ciJD5zWM9Gy2NX1e+q6peO/v9fAPcCeOzQ/WZNQq8WE1Ispifg/f25D80EjbGLyBqAZwD4Yst3Z0VkT0T2Dg4OQh42PRJ6tZiQYrn+evN3cx6aCSbsIvIIAH8H4A9U9X+a36vqtqquq+r6qVOnQh02XTY2qsfBK1eqv6mJOjudyDwwp0kLQd48FZElVKK+q6rvD7FPEhG+1UdK4PJlt+3mMGkhRFaMAHgHgHtV9Y3DTSLRYTpmmczbU5hrTvsc5r6HCMU8B8BrALxARL58tNwaYL8kFkzHLI95TP1rS1JoMqdJCyGyYj6nqqKqT1PVpx8tHw1hHIkE3+orj3l8CmtLUjh3jkkL4OiO88nW1tUxdmBuPZtimNensI2NuRTuLjikwDzCdMzy4FMYqUFhn1dST8ckfvClOFKDwk5ICfApjNRgjJ2QUmC8mRxBj51My7zlXhMyAhR2Mh3zmHtN0qYQR4PCnisl3IAl5F6XcB1IRUGOBoW9Ti6VtJQbMPfc61Kuw7zQVb9LcDRmqOroy7Oe9SxNjp0d1eVl1aqKVsvycrU+NVZXr7ZztqyuTm2ZH7mfR+72zxMu9Vuk/XqKTGd3AwB76qCx9Nhn5NRa5+7pzsg997qU69AklydXH1zqd0EveVHYZ+RUSUu5AXPPvS7lOtSJGV6assFwqd9tjoZIVQa5NXAubn3oJclQTE6P1TmFjUqmxOsQqx5MXVau57Wzc7xtMzSTwLWFYyiGwj5j6hvPl9kNKFL9TdXO0intOsSKM9uEdYwy9K3fMRq4AOdJYe/DmJXU91ilCQhJk5CCVr9n2/ZZF9gxHKpmHTp3zlynQjdwgRxHCnvK+F7kqZ4m2JjMH6Hutbb9tC2Li+bv+txztnu2/t3KiuqJE+bzDO2xB9ofhT1lfC/yFPH/3EJTZBhN0VtZGdagm+5Zm6du2sb1+LZ71rWhmdWp0Pd/oCcACnvK+F7kKfJrc+pMJsOI0Yjbwi/1BsOlAXC952z3rMtxmnUq5BMrPfY5IAePPYOXNUggYtxfPlkoXZ606z1nu2e74vyx69TIMXbmsU+B74s5txrmBjetD0GJOdqknRjvcLje4/V3GUy43nO2e9ZlH232hcq9H/udDRf1D73Mvceu6veYl1KM3ZZJQPIkZu66b+bXEK/WN8a+tGTvS0iwnwkMxRTEGGGRtkrYlh421Y3ODJ14pCRgQ6+za1aMy74T7GcaVdgBvBPA/QDucdmewu5J7BvMtWJPdaOnJDylwobzWhLsZ3IVdqm2HYaIPA/AAwD+RlWf2rX9+vq67u3tDT7u3DAbv6M+iNHycrgY3dpaNR5Gk9XVaqLrGQsL1a3dRKSaFDsWrvYREpIE7zsRuaCq613bBek8VdXPALgcYl+khVnHy8rK8bqHPzzc/l07z6bqUM1pgDZSDhmPPsqsmJz44Q+P///+98ONuucq2LFu9K7MA2bokCnIefRRl3iNywJgDZYYO4CzAPYA7J0+fTpqHKpIYsa3fWLYoWOxLsdmjJ2odt97c9BPgLGzYrqEvb6w87QHsTtypqoUvsOpplhpfQaXIv3oatznpPGnsJdGgqlXQUgw88ALlzcnCxQYKzEa4a77v9T60cBV2IPE2EXk3QD+GcCTROQ+EbktxH5JjYw7cqzkHj9vm3KtSapTLMYg1gxMXR3opu/bslrmgFBZMb+lqjeq6pKq3qSq7wixX1Ij544cG7k3WK6ZOfOSwRNr7uAuB8D0vUheU9oFglkxObGxUeXPXrlS/c1d1IH8G6yh45jMKGUC6Viec5cDsLVV3T9NVOfnaamOS7wm9MIYe4OUOwaJnRAx9pI6/kyxbpEwGVS2emIq/1z6axwAx4rJhJIq9bwyNCsmlUG4QrCzY+4Qj92R6VKOmTtRFPZcmJPefFKjKS4xPM0pHYapPOc5SImksOdC7ul+JjL3jKLRJi4xPNwpHYYpj2277wpwoijsuVDAzXYNBXhG0bDFoF3Ly6XRnNJhSPX6F+BEuQp7PlkxpWQNNMk93a+NWClvIZj6PjJljahemxkEXGura574lO8HDMl0CnF9TPvI/Z0JH1zUP/Ti7bGn6gHMiDk5QI74ekZjnX8K99GQuUCXl6sZf4b8PuV7K4TNvrMopV4mDVBUKCblcEUBN0twfK6Xb/kNaQRSuI9cpxw0CbhpaWs062W1smKfBm7oOYVomENcn6597OxcXbYrK1nV1bKEPeXYWApikRo+Yh2zEWiSyn3kMuWg72K732I6HyH3HeL6dO0jc0esLGFPWTxTEYvUcM3t9im/ofeBaxhjbGwpjy5LmzDVy39xMd55h6ybY3jsKWuJA2UJe8qtbOY3yiiY3s6chQdcy29I7H5lRfW666797dLS9PeR6bxclrY3Ol3ehg3lfNiuiW+IJnaMvcveDChL2FXT7WBMudFJBZtHeuJEJa4u5Tc0bNO2rKxEPnkHTI3byZNmb7u+NMvK9Qkgpse+smLO17fV3xD1vOBc9vKEPWWYFWOnyyNdWXE7/xCx+z7efuxrYhL2WWenr0C7/CZmjF1E9RGPGOf4IezNyBGjsOdC5jeaE10i6/MY7Cq4ruGNGJ20vtjCAy4NVLP8un4TuqE6d65fOGkqLzljR4rCnguZPxo60RUWqaeihapwLoIYIuQTAtvxXEJKbTnsIRrRofb7NkikEwp7Ltg8nQw9CiPN/OGmuIb2kl0E0dfbjyVELoNXzcTTdeiBkyfHa5z6dv4uLpZxb48IhT0XXB+1z52b2tIwmLzyGF7yzk6/VL8pnqJcn1ZcttvZqTqlm/bHygAakq5ZWtgxMhT2XPBJTSu5AsTykvs8CeTe72HLVImB6z1sWkoKO0aGwp4T9Uftea0AMb3kPrH7HDrYTDZOkattCxe5LMQJCntuuHg9JXc25e4l+xAqV9tUXlN2yJuOvbiourBg/o44QWHPjXn32FWvFqVZbHxMbzmml27qPO7bgNkE9Ny5a1/6GusNW9vTAj32wYwq7ABuAfANAN8C8Iau7SnsLXTd+C6VP/SLUr5zd4ZgKs891nFtgt4UZJ9ytt0vS0vtwyeM0QFve1qYh9TeyIwm7AAWAXwbwBMAnADwFQBPsf2Gwq7Xiqit8o8xzoZLKKjPuB2+DcNUlT9WVk6fTkWXcu6TiTJGB7ztPpyncFskxhT2ZwP4WO3znQDutP1m7oW97Qb3GTOljaHCFHp8kb6VeKpBmkIOZjVjSBpg10tbfRuNMbxjW3nl0CmdMGMK+8sBvL32+TUA3mr7zdwLuy0dre9NP1QQXTMZXPfXt6FJzWNvG8zKtcENMWpjMx/9xImrxd1lkLAxG0gSFVdhDzHnqbTNuHfNRiJnRWRPRPYODg4CHDZjTPNeXr4MXLwIXLlS/XWZI3LG0PkcQ29nOkfT+hl95oANMU+m6bhA//lbXcpK2qrP0W/vuAN48MGr1z/4YLUeqO6P8+fN++hrUxdTzxtLunFRf9sChmL8GSueO2WMfcg5+jyuh4zbth13iOdrK9PZlGw2+23Xok6fcEyzfHzefDWNrT+0s55hmU4wYijmOgDfAfB4HHee/rztN3Mv7DEzMFLJimkb8a9tbk/TK/Gux+1qQIaUyc6OOZzi09fQ93xdhb1PLL85Fk2I4ZCHOhLsSO1kNGGvjoVbAXwTVXbMZtf2cy/sqmV4Kz4deyKqN9/cXZl9K3xXp+cQ8TCJWD27JOZ1tI3TXsdU3i6eu+08Z+dTPz/XfXbB1MdejCrsvguFvQD6vPnoMiCXb4WPmTdtE8euMgjBzo77i0a+AjwT/67z9G0whnbWs3PXCoWdxMUmmr7ZIPXK7DLLfDNc5BujDpXZ4+Pp9vXwY6Zadp2HqSG2XV9bo1k/l5gTbBcMhZ3Epc+sP6bKXE/ztFV4k4fcFrcPFR/vMzFym6c7s3OsuHJXZ7hLjN3WKNjG1u9jT8yyKAgKO+mPi5fY5a26CtvSUvvY4W0V3ie04hIfdzlnW1mYjmEb7GpML7Vu+8rK8RyqrlkxITumbY19zv1MI0NhJ/1w9Yq7vE9XoTR1ELZVeJ+4bFd83OWcXcIkzRh4nyW1uHJd1F1nbOqCMfUgUNh9KSFLJQQ2T7dL7GNX+BAe+9Btm7gM8DWVx94HW4bNkHrBLJggUNh9YE7tMT4dnz6V0tRw+lR4n+vUtW3dHts5dt0DruXlG2P3CaOEJJYAs44FgcLuQy7exBhPFT4vu7g+RtsqtW+F9ykDnzx7myD3iR2bxv1xfWHJteMzNDFDJnwqHgyF3Ycc4n9jeTw+L7uEehllzAq/s+M/cFZXCl+fvgYbIV8Eato6pFOcTA6F3YccbuYxbfTJFXchRsPZRzB9PHUfO32eDJrldu7ccUMzm/3IJbzjU3Y7O+7piQyZJA2FvYlNCHK4mad+qhjiVYdulPpeL58wU31ZXOx3L3SlhJ482f69aX2fsutqzEx9GQyZJAmFvY6LEKR+M+fwVGFiSDph24tHfbJLdnb6iXqfht6WA17fn+m7hYVwMfYuO1IKN5JOKOx1chbFGTk8VdjwbTjbzrdtlikXkdrZ6X4JymVxuV9cwj0uMf4hWTGuGT+51QFCYb+KqcMYoUj9qSIkfcImJpHyyTXv6926eOldnnrdY+97nUNm/Axhnu7VEaGw1ynBY583fAcSq08Z1ySEqAPtKYymjknT/ebSADSHJAgZeqmfS59UURdyf7pMGAp7Hd5o3aTmYfl67M0xyuv0EfFm6KZtTJuu0JDJiTB51QsL5o5TVyfE1iAOCYONObE6MUJhb5KacKVEig2fKcZuEy0TfUIx9d/MYtx9GghTWZrux6FhwxCiOnQfpYQ+E4TCTtxJ1cNqy6d38Yjb9uPzUlLbuDh9RX2Wm+7K0GsRopFOoXEhrVDYiTs5eFg2UTcNxVvH1eM2lYXv26p9hTWEMA99Ok2hcSGtUNiJO6l7WDZRny0zfEMcdUF3yfmuf/ZJoTSNHdPG1GHDFBoX0gqFnbiTuofV5S3Xx5wxnYfNY3eZwakp7jOxch0CYIgXPwUU5iShsBM/Uq7IrgJpe/KwCbvvZBnNJ5l62fmEbFJ5IiLZ4CrsCyAEADY2gIsXgStXqr8bG+PbsLsLrK0BCwvV393dav3iovk329vHtl661L7NpUvA5cvmffz4x9eusx2zeZx62Z09C4iYf2vbDyGBGCTsIvIbIvI1EbkiIuuhjCJzyO5uJYr7+5U/u79ffZ6tb+Pmm4HNzeOG4Prr27c7fbpafLhyBVhdNe/PdA7nz1f2zxABTp70289sXzfcUP1epPp/1tANxdSAknJwcetNC4AnA3gSgLsBrLv+jqEYcg1dHbjN4W1vvtltLBnbhB5dYRLfvgfTOays+E8m0tYxu7Q0PESWen8KsYIxY+wUdjIYn8mnVe0iahue2UXU6+mTPn0PtrRRn/3YsnOGxuVTz4AiVlyFXapthyEidwN4varuuWy/vr6ue3tOm5J5YW2tCr80EQHe9a5rY/4LC1eHPOrbX7nif5wmfeqFad+rq1UM3hXTuQHd59d330P3S0ZBRC6oamfYuzPGLiKfFJF7WpaXeBp0VkT2RGTv4ODA56dkHtjaau90VK3i6E1M8emuWPrWFrC8bN/GFFvvom3fy8vVeh9s5+DbV+D6+6H7JWnh4tZ3LWAohoTAFhppMiRW3BWSGRJvDpE2yhg7MQDG2El2+MZ/bSLaJbC2GH0KNF+qMg2z26chSfmdBWJlFGEH8OsA7gPwIwD/BeBjLr+jsJNWQnmTrlMh5u65lnAOxItRPXbfhcJeOEMnaRg6poqr5x/bc429f2a4zB0UdjINY3mRpvHau0ZxHMubNZXDuXPhxD6HUTlJUCjsZBrG8iL7zIk6ZqjCZF/bWO997UnlyYSMhquwc6yY3HF9PXys18ht47WMcZwuDg/b0ydDY7JPGznkQ+xxSa+0DdVAysVF/UMv9NgD4Rr2GLOTLXWPfaxQhY99Q+zpm/3DOHyWgKGYOcC10o5ZuaeMsbsuY4ham32mmHhMexiHLwpXYWcoJmdcwx6xwyP1MM/mJnDmTPX2pkj1tz60big2Nqr9zo6zsgIsLXX/rh6qiBmeatq3ugrcfnuYN1N94Jum84mL+ode6LEHIgWPPaVc6q5Jq+uhiqnsHrsjM6XrQwYDhmLmgBRi7KnFcNtexz9xYr5jz8yKKQZXYWcoJmfaHvfbwh6u2/VhrCwYH1TtnwE/u3OfmCKF2bHIqAQZttcXDttbEKGGqg2Fqz2u283SBQ8Pj9ctL8fpNyCkg2DD9hJixXeo2tjer6sn7mr35ubVog7EyYXP/amApIVLvCb0whh7YbjGcG0pgKFivz6x86bdba/7D0kXrO9/ZaVaTCNRsoOTOAB2npLk6HppJ4SY9RVJ0+9MY890dbJ25dnXbZqnjlwyCFdhZyiGjEdXh6pPiMMUuujbUWwKuQD9cs/b9tfc9+xcU+yAJllDYSfj4fJSjEtWyuteZx//pE8WiElEL1+uXrhaXKw+Ly5Wn7v26SLKs234EhEJDIWdjIfLfKNNMWsbxOquu8J3aJpE9PrrgfPngYceqj4/9FD1uatz00WUZ9uEmiuVkCMo7GQ86mES4NrJq12zUkwpukNCFyZxBfo1Il2NWP1cY75nQOYTl0B86IWdp0RV3bJpTFkpMTob2+wZIyuGEEfArBgyOSFeZR9jwgqbrcxYIQnhKuwMxZA4hJrgwRQiuf32cKELm62Mf5MccVH/0As99jkgpKc79aTQHESLJAIcPXaOFUPisLDQ3skpUqUhpoTJVqB6Grh0qcpg2dpihyaZlFHGihGRvxCRr4vIV0Xk70XkMUP2Rwoip9xsk00inCuUZMnQGPsnADxVVZ8G4JsA7hxuEimCnGLTbbaKXOvFjzURNiEDGSTsqvpxVf3J0ccvALhpuEmkCFLLzbaNnthmqyk0w9f8SQYEi7GLyIcAvEdVd7q2ZYydjEqfMdVTG2eeEASMsYvIJ0XknpblJbVtNgH8BIAxACkiZ0VkT0T2Dg4OXM+DkOGYBvg6c8YcM88plERIg8Eeu4icAXA7gJtV1TKc3TH02Mmo2LJebJ777m7VKDArhiSCq8c+SNhF5BYAbwTwy6rq7IZT2MmomMIqMxheIZkw1tR4bwXwSACfEJEvi8hdA/dHSHi6BuRihygpjOuG/FhVnxjKEEKiMQufnDlzPPxunRRz6wkZAMeKIfPBxkY1jjo7RMkcQGEncbDljU9Farn1hESCwk7CE2pkxxh0TZuXYoNEiCcUdhIeU9546q/jp9wgEeIBhZ2Ex5Rlknr2Sa4NEiENKOwkPDmN7Fgn1waJkAYUdhKeXF/Hz7VBIqQBhZ2EJ9fsk1wbJEIaDHpBiRAjGxvpC3mTmb0cH4ZkDoWdkDo5NkiENGAohhBCCoPCTgghhUFhJ4SQwqCwE0JIYVDYCSGkMIJNZu11UJEDAD8A8N9Hqx7d8n9z3RKA73keqr4P1++b60yfbTbfENhW03c2W13KtI+tMcu0zdbUy7TN5nkt0z62xirTLltzLdNVVT3VuUdVnWQBsG37v7kOwN6QY7h+31xn+myzObStpu9strqUaR9bY5apoSyTLtMQtpZSpn1sjVWmfa5/TmXatUwZivlQx/+m7/sew/X75jrT5y6bfbH91vSdzdYcy7T+fy5lWv9/3svU9P0UZdr129zL1MokoZg+iMieOkzimgK0NTy52AnkY2sudgL52JqKnTl1nm5PbYAHtDU8udgJ5GNrLnYC+diahJ3ZeOyEEELcyMljJ4QQ4gCFnRBCCoPCTgghhVGEsIvIaRH5oIi8U0TeMLU9NkTkuSJyl4i8XUQ+P7U9JkRkQUS2ROQtInJmantsiMjzReSzR+X6/KntsSEiJ0Xkgoi8aGpbbIjIk4/K830icm5qe2yIyEtF5G0i8gER+ZWp7TEhIk8QkXeIyPtiH2tyYT8S4/tF5J7G+ltE5Bsi8i0Hsf45AB9R1d8B8JSUbVXVz6rq7QA+DOB8qnYCeAmAxwL4MYD7YtgZ0FYF8ACAh8WyNZCdAPBHAN4bw8aaTSHu03uP7tNXAIiWvhfI1n9Q1d8F8FoAv5mwnd9R1dti2Nd2sEkXAM8D8EwA99TWLQL4NoAnADgB4CuoBPsXUAlifflpACsAPg3gUwB+O2Vba797L4BHpWongDcA+L2j374v5TIFsHD0u58BsJuwnS8E8EpUAvSilMv06DcvBvB5AK9K3daj3/0lgGdmYGe0+vT/x4h9AMdCW2sU2LMBfKz2+U4Ad1p+/3oAzxuj0IbaerTNaQBvS9lOAK8G8Iqj/9+Tsq217U7EvP4BynQLwJsAfBzAB3DUIKVoa2NfH0n5+gMQAH8O4IUp21nbLrqwpzo13mMB/Eft830Afsmy/T8B+BMReRWAixHtasPXVgC4DcBfRbOoHV873w/gLSLyXACfiWlYC162isjLAPwqgMcAeGtc067Cy05V3QQAEXktgO+p6pWo1l2Nb5k+H8DLAPwUgI9GtexafO/V30f1NPRoEXmiqt4V07gavmW6gqpxf4aI3KmqfxbLsFSFXVrWGd+kUtV7ALw8njlWvGwFAFX940i22PAt00NUDdAU+Nr6flQN0dh4X3sAUNW/Dm9KJ75lejeAu2MZ04GvrW8G8OZ45hjxtfP7AG6PZ84xk3eeGrgPwONqn28C8J8T2dJFLrbmYieQj6252AnQ1hgka2eqwv6vAH5WRB4vIidQdTh9cGKbTORiay52AvnYmoudAG2NQbp2xg7iO3QkvBvAd3GcVnfb0fpbAXwTVa/z5tR25mRrLnbmZGsudtLW+bZztnAQMEIIKYxUQzGEEEJ6QmEnhJDCoLATQkhhUNgJIaQwKOyEEFIYFHZCCCkMCjshhBQGhZ0QQgqDwk4IIYXxf5pfjC5ZurbIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(aut, eci_z, 'ro')\n",
    "plt.xscale('log');\n",
    "#plt.yscale('log');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# iteration methods of eci\n",
    "\n",
    "https://github.com/hashc/Skills/blob/master/iter_eci.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:34:26.796560Z",
     "start_time": "2019-07-04T14:34:26.787331Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2.75, 2.75, 2.75, 2.75]), array([2.25, 2.25, 2.25, 2.25]))"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as  np\n",
    "import copy \n",
    "c=np.array([[1.,1.,1.,1.],[0.,1.,0.,0.],[0.,0.,1.,1.],[0.,0.,0.,1.]])\n",
    "def eci(c):\n",
    "\tc.astype(np.float64)\n",
    "\ta=c.sum(1)\n",
    "\tb = c.sum(0)\n",
    "\ta0 = copy.deepcopy(a)\n",
    "\tb0 = copy.deepcopy(b)\n",
    "\tfor _ in range(100):\n",
    "\t\ta1 = copy.deepcopy(a)\n",
    "\t\tfor i in range(len(a)):\n",
    "\t\t\ta1[i] = 1.0*sum([b[j] for j in range(len(a)) if c[i,j]])/a0[i]\n",
    "\t\tb1 = copy.deepcopy(b)\n",
    "\t\tfor j in range(len(b)):\n",
    "\t\t\tb1[j] =  1.0*sum([a[i] for i in range(len(b)) if c[i,j]])/b0[j]\n",
    "\t\t#print(a1,b1)\n",
    "\t\ta=a1\n",
    "\t\tb=b1\n",
    "\treturn a,b\n",
    "\n",
    "eci(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# END"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:23:52.105008Z",
     "start_time": "2019-07-04T14:23:52.093165Z"
    }
   },
   "outputs": [],
   "source": [
    "import pickle \n",
    "\n",
    "fr = open('cdata_res.pkl','rb')  \n",
    "dat = pickle.load(fr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:23:52.596741Z",
     "start_time": "2019-07-04T14:23:52.592556Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['ECI', 'density', 'cog', 'PCI'])"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-04T14:23:53.074770Z",
     "start_time": "2019-07-04T14:23:53.070748Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.11142756821655171"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat['ECI']['~019343c06c6a15aa47']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:anaconda]",
   "language": "python",
   "name": "conda-env-anaconda-py"
  },
  "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.5.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}