view notebooks/test_hubness.ipynb @ 55:98cd5317e504 branch-tests

updated notebooks

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pickle\n",
    "from scipy.stats import pearsonr\n",
    "from scipy.stats import skew\n",
    "import sys\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "sys.path.append('../')\n",
    "import scripts.outliers as outliers\n",
    "import scripts.utils_spatial as utils_spatial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_list, Y, Yaudio = pickle.load(open('../data/lda_data_melodia_8.pickle','rb'))\n",
    "#ddf = outliers.load_metadata(Yaudio, metadata_file='../data/metadata.csv')\n",
    "#w, data_countries = utils_spatial.get_neighbors_for_countries_in_dataset(Y)\n",
    "#w_dict = utils_spatial.from_weights_to_dict(w, data_countries)\n",
    "X = np.concatenate(X_list, axis=1)\n",
    "\n",
    "# global outliers\n",
    "#df_global, threshold, MD = outliers.get_outliers_df(X, Y, chi2thr=0.999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8200, 380)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## distance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "D = pairwise_distances(X, metric='mahalanobis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8200, 8200)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.savetxt('../data/D_mahal.csv', D)\n",
    "D = np.loadtxt('../data/D_mahal.csv')\n",
    "D.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEACAYAAAB78OvLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGWtJREFUeJzt3W2sXdV95/HvDwyUPIwd05F5sCFWa6Q4StrADM50WvUy\nJGCiCniRgKsJeBKrL0pnYPqiE8xIYKtRO54KUaoRRGrzYFCDQdASoqBgB7iaSiNwEiClIYzNqJ7B\nl9pEBpO0o05t8Z8XZ128ubm2r303vvccfz/SkddZe6919jpczu/svfbeJ1WFJEl9OGWuN0CSNDoM\nFUlSbwwVSVJvDBVJUm8MFUlSbwwVSVJvZhQqSRYleSjJj5K8mGRVksVJtiXZkWRrkkWd9dcn2Znk\npSSXd+ovTvJCW3ZXp/6MJA+0+qeTXNBZtra9xo4kN3Tqlyd5prXZkuS02b8dkqTZmOmeyl3AY1X1\nIeCjwEvALcC2qroQeKI9J8lK4DpgJbAauDtJWj/3AOuqagWwIsnqVr8O2Nfq7wQ2tb4WA7cBl7TH\n7UkWtjabgDtamzdaH5KkOXTUUGkf4r9WVV8BqKqDVfUmcBWwua22Gbimla8G7q+qA1W1C3gZWJXk\nHOD9VbW9rXdvp023r4eBy1r5CmBrVe2vqv3ANuDKFlKXAg9N8/qSpDkykz2V5cCPk3w1ybNJ/jTJ\ne4ElVbW3rbMXWNLK5wK7O+13A+dNUz/R6mn/vgKD0ALeTHLWEfpaDOyvqrem6UuSNEdmEioLgIuA\nu6vqIuAfaIe6JtXgXi8n6n4v3ldGkuapBTNYZzewu6q+254/BKwH9iQ5u6r2tENbr7XlE8CyTvul\nrY+JVp5aP9nmfODVJAuAhVW1L8kEMNZpswx4EngdWJTklLa3srT18Q5JDCBJOkZVlaOvNb2j7qlU\n1R7glSQXtqpPAD8EvgmsbXVrgUda+VFgTZLTkywHVgDbWz8/aWeOBbge+EanzWRfn2Yw8Q+wFbi8\nnX32AeCTwONtz+gp4DPTvP7U7R/Jx+233z7n2+D4HJ/jG73HbM1kTwXgPwB/nuR04H8BnwNOBR5M\nsg7YBVzbPsRfTPIg8CJwELixDm3pjcDXgDMZnE327Vb/ZeC+JDuBfcCa1tfrSX4fmNxL2liDCXuA\nLwBbknwReLb1IUmaQzMKlar6AfAvp1n0icOs/wfAH0xT/33gI9PU/z9aKE2z7KvAV6ep/1tg1RE3\nXBw6m3ugj28iknQ4XlE/pMbGxo5h7RN5HkU/jm18w8fxDbdRH99sZJS/uSapUR7fTAz2VCbfg7in\nIumIklDv5kS9JEkzNdOJeo2I7hyLey2S+uaeykln+OZXJA0PQ0WS1BtDRZLUG0NFktQbQ0WS1BvP\n/hpBU6+il6QTxT2VkeVZXpJOPENFktQbQ0WS1BtDRZLUG0NFktQbQ0WS1BtDRZLUG0NFktQbQ0WS\n1BuvqD+J+dsqkvrmnspJzavuJfXLUJEk9cZQkST1xlCRJPXGUJEk9cZQkST1xlCRJPXGUJEk9WZG\noZJkV5K/TvJcku2tbnGSbUl2JNmaZFFn/fVJdiZ5KcnlnfqLk7zQlt3VqT8jyQOt/ukkF3SWrW2v\nsSPJDZ365UmeaW22JDlttm+GJGl2ZrqnUsBYVX2sqi5pdbcA26rqQuCJ9pwkK4HrgJXAauDuHLp0\n+x5gXVWtAFYkWd3q1wH7Wv2dwKbW12LgNuCS9rg9ycLWZhNwR2vzRutDkjSHjuXwV6Y8vwrY3Mqb\ngWta+Wrg/qo6UFW7gJeBVUnOAd5fVdvbevd22nT7ehi4rJWvALZW1f6q2g9sA65sIXUp8NA0ry9J\nmiPHsqfynSTfS/JbrW5JVe1t5b3AklY+F9jdabsbOG+a+olWT/v3FYCqOgi8meSsI/S1GNhfVW9N\n05ckaY7M9IaS/7qq/i7JPwe2JXmpu7CqKsmJuomUN6uSpHlqRqFSVX/X/v1xkr9kML+xN8nZVbWn\nHdp6ra0+ASzrNF/KYA9jopWn1k+2OR94NckCYGFV7UsyAYx12iwDngReBxYlOaXtrSxtffyMDRs2\nvF0eGxtjbGxsutWGXveOw5I0U+Pj44yPj/fWX452y/Mk7wFOraqfJnkvsBXYCHyCweT6piS3AIuq\n6pY2Uf91BsFzHvAd4Bfb3swzwE3AduBbwJ9U1beT3Ah8pKp+O8ka4JqqWtMm6r8HXMRgTuf7wEVV\ntT/Jg8DDVfVAki8Bz1fVl6Zse50st3QfhMrkWI+9fLK8T5KOLAlVddzfUmcSKsuBv2xPFwB/XlV/\n2D7wH2Swh7ELuLZNppPkVuDzwEHg5qp6vNVfDHwNOBN4rKpuavVnAPcBHwP2AWvaJD9JPgfc2l7/\ni1W1ubNdWxjMrzwLfLaqDkzZdkPFUJF0DN71UBlmhsrMyyfL+yTpyGYbKl5RL0nqjaEiSeqNoSJJ\n6s1Mr1PRiOuekuz8iqTj5Z6KmsLrSiXNlqEiSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknqjaEi\nSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknq\njaEiSeqNoSJJ6o2hIknqjaEiSerNgrneAM0/Sd4uV9UcbomkYeOeiqZR7SFJx2ZGoZLk1CTPJflm\ne744ybYkO5JsTbKos+76JDuTvJTk8k79xUleaMvu6tSfkeSBVv90kgs6y9a219iR5IZO/fIkz7Q2\nW5KcNts3QpI0ezPdU7kZeJFDX19vAbZV1YXAE+05SVYC1wErgdXA3Tl0LOUeYF1VrQBWJFnd6tcB\n+1r9ncCm1tdi4Dbgkva4PcnC1mYTcEdr80brQ5I0x44aKkmWAp8C/gyYDIirgM2tvBm4ppWvBu6v\nqgNVtQt4GViV5Bzg/VW1va13b6dNt6+Hgcta+Qpga1Xtr6r9wDbgyhZSlwIPTfP6kqQ5NJM9lTuB\n3wPe6tQtqaq9rbwXWNLK5wK7O+vtBs6bpn6i1dP+fQWgqg4CbyY56wh9LQb2V9Vb0/R1Ukny9kOS\n5oMjnv2V5DeA16rquSRj061TVZXkRM3qHvPrbNiw4e3y2NgYY2NjPW7OfDD5lhgsko7d+Pg44+Pj\nvfV3tFOKfwW4KsmngJ8D/lmS+4C9Sc6uqj3t0NZrbf0JYFmn/VIGexgTrTy1frLN+cCrSRYAC6tq\nX5IJYKzTZhnwJPA6sCjJKW1vZWnrY1rdUJEkvdPUL9sbN26cVX9HPPxVVbdW1bKqWg6sAZ6squuB\nR4G1bbW1wCOt/CiwJsnpSZYDK4DtVbUH+EmSVW1O5HrgG502k319msHEP8BW4PIki5J8APgk8HgN\nLpx4CvjMNK8vSZpDx3rx4+Sxlv8CPJhkHbALuBagql5M8iCDM8UOAjfWoavnbgS+BpwJPFZV3271\nXwbuS7IT2McgvKiq15P8PvDdtt7GNmEP8AVgS5IvAs+2PiRJcyyjfMV0khrx8fHOOZX+y6P8/kn6\nWUmoquOepPWKeklSbwwVSVJvDBVJUm8MFUlSbwwVSVJvDBVJUm8MFUlSbwwVSVJvDBVJUm8MFUlS\nbwwVSVJvDBVJUm8MFUlSb4711vc6yXR/qtg7Fks6GvdUdBTFcfyKs6STlKEiSeqNoSJJ6o2hIknq\njaEiSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknqjaEiSeqNoSJJ6o2hIknqjaEiSerNEUMlyc8l\neSbJ80leTPKHrX5xkm1JdiTZmmRRp836JDuTvJTk8k79xUleaMvu6tSfkeSBVv90kgs6y9a219iR\n5IZO/fK2XTuTbElyWl9viCTp+B0xVKrqH4FLq+qXgY8Clyb5VeAWYFtVXQg80Z6TZCVwHbASWA3c\nnUM/yHEPsK6qVgArkqxu9euAfa3+TmBT62sxcBtwSXvcnmRha7MJuKO1eaP1IUmaY0c9/FVV/7cV\nTwdOZfAhfhWwudVvBq5p5auB+6vqQFXtAl4GViU5B3h/VW1v693badPt62Hgsla+AthaVfuraj+w\nDbiyhdSlwEPTvL4kaQ4dNVSSnJLkeWAv8FRV/RBYUlV72yp7gSWtfC6wu9N8N3DeNPUTrZ727ysA\nVXUQeDPJWUfoazGwv6remqYvSdIcOurPCbcP719uh54eT3LplOWV5ET9NOAxv86GDRveLo+NjTE2\nNtbj5kjScBsfH2d8fLy3/mb8G/VV9WaSbwEXA3uTnF1Ve9qhrdfaahPAsk6zpQz2MCZaeWr9ZJvz\ngVeTLAAWVtW+JBPAWKfNMuBJ4HVgUZJTWuAtbX1MqxsqkqR3mvple+PGjbPq72hnf/385JldSc4E\nPgk8BzwKrG2rrQUeaeVHgTVJTk+yHFgBbK+qPcBPkqxqcyLXA9/otJns69MMJv4BtgKXJ1mU5APt\ntR+vqgKeAj4zzetLkuZQBp/Rh1mYfITBRPgp7XFfVf1ROzPrQQZ7GLuAa9tkOkluBT4PHARurqrH\nW/3FwNeAM4HHquqmVn8GcB/wMWAfsKZN8pPkc8CtbXO+WFWbW/1yYAuD+ZVngc9W1YFptr+ONL5h\ndOhkukmT48sJKHdedcTeV0kDSaiqqR80M28/yh8OoxsqJzJIpi+P2vsqaWC2oeIV9ZKk3hgqkqTe\nGCqSpN4YKpKk3hgqkqTeGCqSpN4YKpKk3hgqkqTeGCqSpN4YKpKk3hgqkqTeGCqSpN4YKpKk3hgq\nkqTeGCqSpN4YKpKk3hgqkqTeLJjrDdBw6v6ssb8CKWmSeyo6TsWhnxeWpAFDRZLUG0NFktQbQ0WS\n1BtDRZLUG0NFktQbQ0WS1BtDRZLUG0NFktQbQ0WS1JujhkqSZUmeSvLDJH+T5KZWvzjJtiQ7kmxN\nsqjTZn2SnUleSnJ5p/7iJC+0ZXd16s9I8kCrfzrJBZ1la9tr7EhyQ6d+eZJnWpstSU7r4w2RJB2/\nmeypHAB+t6o+DHwc+J0kHwJuAbZV1YXAE+05SVYC1wErgdXA3Tl0o6h7gHVVtQJYkWR1q18H7Gv1\ndwKbWl+LgduAS9rj9iQLW5tNwB2tzRutD0nSHDpqqFTVnqp6vpX/HvgRcB5wFbC5rbYZuKaVrwbu\nr6oDVbULeBlYleQc4P1Vtb2td2+nTbevh4HLWvkKYGtV7a+q/cA24MoWUpcCD03z+pKkOXJMcypJ\nPgh8DHgGWFJVe9uivcCSVj4X2N1ptptBCE2tn2j1tH9fAaiqg8CbSc46Ql+Lgf1V9dY0fUmS5siM\nb32f5H0M9iJurqqfTrn1eSU5UbesPabX2bBhw9vlsbExxsbGet4cSRpe4+PjjI+P99bfjEKlTYI/\nDNxXVY+06r1Jzq6qPe3Q1mutfgJY1mm+lMEexkQrT62fbHM+8GqSBcDCqtqXZAIY67RZBjwJvA4s\nSnJK21tZ2vr4Gd1Q0bvD31aRhtfUL9sbN26cVX8zOfsrwJeBF6vqjzuLHgXWtvJa4JFO/ZokpydZ\nDqwAtlfVHuAnSVa1Pq8HvjFNX59mMPEPsBW4PMmiJB8APgk8XoNPrqeAz0zz+jrh/G0VSQM52jfL\nJL8K/Hfgrzn0ybEe2A48yGAPYxdwbZtMJ8mtwOeBgwwOlz3e6i8GvgacCTxWVZOnJ58B3MdgvmYf\nsKZN8pPkc8Ct7XW/WFWbW/1yYAuD+ZVngc9W1YEp216j9s15kMeTY5p/5VF7v6WTTRKqKkdf8zDt\nR/lDYFRCpXt4aWDuw8NQkUbTbEPFK+qHhoeYJM1/hookqTeGiiSpN4aKJKk3hookqTeGiiSpN4aK\nJKk3hookqTeGiiSpNzO+S7E0E95cUjq5uaeinnnlv3QyM1QkSb0xVCRJvTFUJEm9MVQkSb0xVCRJ\nvTFUJEm9MVQkSb0xVCRJvTFUJEm9MVQkSb3x3l9613gfMOnk456K3kXeB0w62RgqkqTeGCqSpN4Y\nKpKk3hgqkqTeGCqSpN4cNVSSfCXJ3iQvdOoWJ9mWZEeSrUkWdZatT7IzyUtJLu/UX5zkhbbsrk79\nGUkeaPVPJ7mgs2xte40dSW7o1C9P8kxrsyXJabN9IyRJszeTPZWvAqun1N0CbKuqC4En2nOSrASu\nA1a2Nnfn0MUK9wDrqmoFsCLJZJ/rgH2t/k5gU+trMXAbcEl73J5kYWuzCbijtXmj9TFSkrz9kKRh\ncdRQqaq/YvDB3XUVsLmVNwPXtPLVwP1VdaCqdgEvA6uSnAO8v6q2t/Xu7bTp9vUwcFkrXwFsrar9\nVbUf2AZc2ULqUuChaV5/xHidh6ThcrxzKkuqam8r7wWWtPK5wO7OeruB86apn2j1tH9fAaiqg8Cb\nSc46Ql+Lgf1V9dY0fWmecs9LOjnMeqK+BvffOFFfp/3aPrTc65JOBsd776+9Sc6uqj3t0NZrrX4C\nWNZZbymDPYyJVp5aP9nmfODVJAuAhVW1L8kEMNZpswx4EngdWJTklLa3srT1Ma0NGza8XR4bG2Ns\nbOxwq0rSSWd8fJzx8fHe+stMbvSX5IPAN6vqI+35f2Uwub4pyS3Aoqq6pU3Uf53BxPp5wHeAX6yq\nSvIMcBOwHfgW8CdV9e0kNwIfqarfTrIGuKaq1rSJ+u8BFwEBvg9cVFX7kzwIPFxVDyT5EvB8VX1p\nmu2uYb2R4eAw0eS2j1Z5WP+bSCeDJFTVcR+nPmqoJLkf+HXg5xnMn9wGfAN4kMEexi7g2jaZTpJb\ngc8DB4Gbq+rxVn8x8DXgTOCxqrqp1Z8B3Ad8DNgHrGmT/CT5HHBr25QvVtXmVr8c2MJgfuVZ4LNV\ndWCabTdU5mF5WP+bSCeDdz1UhpmhMj/Lw/rfRDoZzDZUvKJektQbQ0WS1Bt/+VEnnL8IKY0u91Q0\nB7xmRRpVhookqTeGiiSpN4aKJKk3hookqTee/aU55Zlg0mhxT0VzzDPBpFFiqEiSeuPhr3nEH7CS\nNOwMlXmnexPGk4vzK9Lw8/CX5hHnV6RhZ6hIknpjqEiSeuOciuYl51ek4eSeiuYp51ekYWSoSJJ6\nY6hIknrjnIrmPedXpOHhnoqGgPMr0rAwVCRJvfHw1xzzfl/HxkNh0vzmnsq84OGdmfO9kuYz91Q0\ntNxrkeYfQ0VD7NAdnQ93GNGwkU4sQ0UjovuTAdOHjQEjvfuGek4lyeokLyXZmeQLc709mo8OzcEk\nefsh6d0xtKGS5FTgvwGrgZXAbyb50Nxu1cz08+E23tfmzFPj70Kf8ydgxsfH5+R1TxTHd/Ia2lAB\nLgFerqpdVXUA2AJcPcfbdAxmexbTeE/bMV+Nv8v9z23AjPqHkuM7eQ3znMp5wCud57uBVVNX6h5H\n97CHpnf0Cf8Z9eKcjTTUoTKj/4NPOWWwM7Zr1y4uuOCCd3WDjsRAGxbTT/jPpHys/403btx4/Jt5\nGDMNtiNtq+Go2ciw/gEl+TiwoapWt+frgbeqalNnneEcnCTNoao67m/BwxwqC4D/CVwGvApsB36z\nqn40pxsmSSexoT38VVUHk/x74HHgVODLBookza2h3VORJM0/w3xK8WGN2kWRSZYleSrJD5P8TZKb\nWv3iJNuS7EiyNcmiud7W45Xk1CTPJflmez5KY1uU5KEkP0ryYpJVIza+9e1v84UkX09yxjCPL8lX\nkuxN8kKn7rDjaePf2T5zLp+brZ65w4zvj9rf5w+S/EWShZ1lxzS+kQuVYb4o8ggOAL9bVR8GPg78\nThvTLcC2qroQeKI9H1Y3Ay9y6HSqURrbXcBjVfUh4KPAS4zI+JJ8EPgt4KKq+giDQ9FrGO7xfZXB\n50fXtONJshK4jsFnzWrg7iTz/XN1uvFtBT5cVb8E7ADWw/GNb74P/ngM+UWRP6uq9lTV863898CP\nGFyncxWwua22GbhmbrZwdpIsBT4F/BmD83RhdMa2EPi1qvoKDOYCq+pNRmR8wE8YfOl5Tzt55j0M\nTpwZ2vFV1V8Bb0ypPtx4rgbur6oDVbULeJnBZ9C8Nd34qmpbVb3Vnj4DLG3lYx7fKIbKdBdFnjdH\n29K79s3wYwz+wy+pqr1t0V5gyRxt1mzdCfwe8FanblTGthz4cZKvJnk2yZ8meS8jMr6qeh24A/g/\nDMJkf1VtY0TG13G48ZzL4DNm0ih83nweeKyVj3l8oxgqI3vmQZL3AQ8DN1fVT7vLanDGxdCNPclv\nAK9V1XMc2kt5h2EdW7MAuAi4u6ouAv6BKYeChnl8SX4B+I/ABxl8AL0vyWe76wzz+KYzg/EM7ViT\n/Gfgn6rq60dY7YjjG8VQmQCWdZ4v451JO5SSnMYgUO6rqkda9d4kZ7fl5wCvzdX2zcKvAFcl+Vvg\nfuDfJLmP0RgbDP72dlfVd9vzhxiEzJ4RGd+/AP5HVe2rqoPAXwD/itEZ36TD/T1O/bxZ2uqGTpJ/\nx+Aw9L/tVB/z+EYxVL4HrEjywSSnM5hkenSOt2lWMrinxpeBF6vqjzuLHgXWtvJa4JGpbee7qrq1\nqpZV1XIGE7xPVtX1jMDYYDAfBryS5MJW9Qngh8A3GYHxMTjp4ONJzmx/p59gcMLFqIxv0uH+Hh8F\n1iQ5PclyYAWDC7GHSpLVDA5BX11V/9hZdOzjq6qRewBXMrja/mVg/VxvTw/j+VUG8w3PA8+1x2pg\nMfAdBmdrbAUWzfW2znKcvw482sojMzbgl4DvAj9g8E1+4YiN7z8xCMoXGExinzbM42Owx/wq8E8M\n5mc/d6TxALe2z5qXgCvmevuPY3yfB3YC/7vz+XL38Y7Pix8lSb0ZxcNfkqQ5YqhIknpjqEiSemOo\nSJJ6Y6hIknpjqEiSemOoSJJ6Y6hIknrz/wF0zsvts73EjAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb84d839fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(D.ravel(), bins=100);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## n-occurrence and stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def n_occurrence_from_D(D, k=10, n_items=None):\n",
    "    if n_items is None:\n",
    "        n_items = len(D)\n",
    "    sort_idx = np.argsort(D, axis=1)\n",
    "    D_k = sort_idx[:, 1:(k+1)]  # nearest neighbour is the item itself\n",
    "    N_k = np.bincount(D_k.astype(int).ravel(), minlength=n_items)\n",
    "    return N_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "60\n",
      "10.0869227031\n",
      "0.0\n",
      "60.0\n",
      "389.813680947\n",
      "142\n",
      "7.70305908648\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb823513f90>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFnJJREFUeJzt3X+MXeWd3/H3Jzg0hhhbNJX5KeF2hwWv2AbcwdkmKUPL\nUqfaAH+BkRZZWzdS5N2StlIaO1ILf7EkUrUhqkzlbRJMlHjlzSbINKyXCeFKG7XFCYGNg+MFt+s0\nno2HNNnCrjdWQXz7x318fHc6eMb2eOZeeL+kEc95znPO+T7WcD9znnPvTKoKSZIA3rHUBUiShoeh\nIEnqGAqSpI6hIEnqGAqSpI6hIEnqzBkKSbYleSHJ/iRfTvK3klycZDLJi0meTLJqxviXkhxMcutA\n/7p2jpeSPHSuJiRJOnOnDIUkVwEfAW6oquuA84CNwFZgsqquBp5q2yRZC9wFrAU2ANuTpJ3uYWBz\nVY0BY0k2LPhsJElnZa47hVeB14ALkiwDLgD+HLgN2NnG7ATuaO3bgV1V9VpVHQYOAeuTXAqsqKp9\nbdyjA8dIkobEKUOhqn4G/Afgf9EPg/9TVZPA6qqabsOmgdWtfRlwZOAUR4DLZ+mfav2SpCEy1/LR\n3wP+FXAV/Rf2dyf59cEx1f89Gf6uDEl6C1g2x/5/APzXqvopQJKvAr8CHE1ySVUdbUtDL7fxU8CV\nA8dfQf8OYaq1B/unZrtgEgNGkk5TVWXuUXOb65nCQeB9SZa3B8a3AAeAx4FNbcwm4LHW3gNsTHJ+\nkjXAGLCvqo4CryZZ385zz8Ax/5+qGsmv++67b8lrsP6lr8P6R+9rlGuvWtifo095p1BVf5LkUeA7\nwBvAd4EdwApgd5LNwGHgzjb+QJLd9IPjdWBLnax4C/AIsBx4oqr2LuhMJElnba7lI6rq08CnZ3T/\njP5dw2zjHwAemKX/WeC6M6hRkrRI/ETzApqYmFjqEs6K9S8t6186o1z7QstCr0edrSQ1bDVJ0jBL\nQi3Sg2ZJ0tvInM8UlkKv1wPghhtu4KKLLlraYiTpbWQol49WrryJn//8Wb71rW8yPj6+1CVJ0lB7\nyy8fvfJKj+XLr1nqMiTpbWcoQ0GStDQMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQ\nJHUMBUlSx1CQJHUMBUlSx1CQJHXmDIUkv5jkuYGvV5Lcm+TiJJNJXkzyZJJVA8dsS/JSkoNJbh3o\nX5dkf9v30LmalCTpzMwZClX1p1V1fVVdD6wD/hr4GrAVmKyqq4Gn2jZJ1gJ3AWuBDcD2JCd+z/fD\nwOaqGgPGkmxY6AlJks7c6S4f3QIcqqofAbcBO1v/TuCO1r4d2FVVr1XVYeAQsD7JpcCKqtrXxj06\ncIwkaQicbihsBHa19uqqmm7taWB1a18GHBk45ghw+Sz9U61fkjQk5h0KSc4HPgz8/sx91f+bnsP1\ndz0lSadt2WmM/RDwbFX9pG1PJ7mkqo62paGXW/8UcOXAcVfQv0OYau3B/qnZL3U/x49PsWPHDo4d\nO8bExMRplClJb229Xo9er3dOzp3+D/nzGJj8HvCHVbWzbX8a+GlVfSrJVmBVVW1tD5q/DNxIf3no\nG8AvVFUleQa4F9gHfB34bFXtnXGdgmLlynEmJ7czPj6+QFOVpLemJFRV5h45t3ndKSS5kP5D5o8M\ndD8I7E6yGTgM3AlQVQeS7AYOAK8DW+pk8mwBHgGWA0/MDARJ0tKaVyhU1THgPTP6fkY/KGYb/wDw\nwCz9zwLXnX6ZkqTF4CeaJUkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS\n1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEmdeYVCklVJvpLkB0kOJFmf\n5OIkk0leTPJkklUD47cleSnJwSS3DvSvS7K/7XvoXExIknTm5nun8BDwRFVdC/wycBDYCkxW1dXA\nU22bJGuBu4C1wAZge5K08zwMbK6qMWAsyYYFm4kk6azNGQpJVgIfrKrPA1TV61X1CnAbsLMN2wnc\n0dq3A7uq6rWqOgwcAtYnuRRYUVX72rhHB46RJA2B+dwprAF+kuQLSb6b5HeTXAisrqrpNmYaWN3a\nlwFHBo4/Alw+S/9U65ckDYll8xxzA/BbVfXtJJ+hLRWdUFWVpBaurPs5fnyKHTt2cOzYMSYmJhbu\n1JI04nq9Hr1e75ycO1Wnfi1Pcgnw36pqTdv+ALAN+LvAzVV1tC0NPV1V1yTZClBVD7bxe4H7gB+2\nMde2/ruBm6rqozOuV1CsXDnO5OR2xsfHF3K+kvSWk4Sqytwj5zbn8lFVHQV+lOTq1nUL8ALwOLCp\n9W0CHmvtPcDGJOcnWQOMAfvaeV5t71wKcM/AMZKkITCf5SOAfwl8Kcn5wP8AfgM4D9idZDNwGLgT\noKoOJNkNHABeB7bUyduRLcAjwHL672bau0DzkCQtgDmXjxaby0eSdHoWdflIkvT2YShIkjqGgiSp\nYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhI\nkjqGgiSpYyhIkjqGgiSpM69QSHI4yfeSPJdkX+u7OMlkkheTPJlk1cD4bUleSnIwya0D/euS7G/7\nHlr46UiSzsZ87xQKmKiq66vqxta3FZisqquBp9o2SdYCdwFrgQ3A9iQn/qD0w8DmqhoDxpJsWKB5\nSJIWwOksH2XG9m3AztbeCdzR2rcDu6rqtao6DBwC1ie5FFhRVfvauEcHjpEkDYHTuVP4RpLvJPlI\n61tdVdOtPQ2sbu3LgCMDxx4BLp+lf6r1S5KGxLJ5jnt/Vf04yd8BJpMcHNxZVZWkFq6s+zl+fIod\nO3Zw7NgxJiYmFu7UkjTier0evV7vnJw7Vaf3Wp7kPuCvgI/Qf85wtC0NPV1V1yTZClBVD7bxe4H7\ngB+2Mde2/ruBm6rqozPOX1CsXDnO5OR2xsfHz3KKkvTWloSqmrnEf0bmXD5KckGSFa19IXArsB/Y\nA2xqwzYBj7X2HmBjkvOTrAHGgH1VdRR4Ncn69uD5noFjJElDYD7LR6uBr7U3EC0DvlRVTyb5DrA7\nyWbgMHAnQFUdSLIbOAC8Dmypk7cjW4BHgOXAE1W1dwHnIkk6S6e9fHSuuXwkSadnUZePJElvH4aC\nJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKlj\nKEiSOoaCJKljKEiSOoaCJKljKEiSOvMKhSTnJXkuyeNt++Ikk0leTPJkklUDY7cleSnJwSS3DvSv\nS7K/7Xto4aciSTpb871T+BhwAKi2vRWYrKqrgafaNknWAncBa4ENwPYkJ/6Y9MPA5qoaA8aSbFiY\nKUiSFsqcoZDkCuCfAf8ZOPECfxuws7V3Ane09u3Arqp6raoOA4eA9UkuBVZU1b427tGBYyRJQ2I+\ndwq/A3wceGOgb3VVTbf2NLC6tS8DjgyMOwJcPkv/VOuXJA2RZafameTXgJer6rkkE7ONqapKUrPt\nO3P3c/z4FDt27ODYsWNMTMx6aUl6W+r1evR6vXNy7lS9+et5kgeAe4DXgXcBFwFfBcaBiao62paG\nnq6qa5JsBaiqB9vxe4H7gB+2Mde2/ruBm6rqo7Ncs6BYuXKcycntjI+PL+B0JemtJwlVlblHzu2U\ny0dV9cmqurKq1gAbgW9W1T3AHmBTG7YJeKy19wAbk5yfZA0wBuyrqqPAq0nWtwfP9wwcI0kaEqdc\nPprFiduKB4HdSTYDh4E7AarqQJLd9N+p9DqwpU7eimwBHgGWA09U1d6zK12StNBOuXy0FFw+kqTT\ns2jLR5KktxdDQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLU\nMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSZ1ThkKSdyV5JsnzSQ4k+e3Wf3GSySQvJnky\nyaqBY7YleSnJwSS3DvSvS7K/7Xvo3E1JknSmThkKVXUcuLmq3gv8MnBzkg8AW4HJqroaeKptk2Qt\ncBewFtgAbE9y4o9JPwxsrqoxYCzJhnMxIUnSmZtz+aiq/ro1zwfOA/4CuA3Y2fp3Ane09u3Arqp6\nraoOA4eA9UkuBVZU1b427tGBYyRJQ2LOUEjyjiTPA9PA01X1ArC6qqbbkGlgdWtfBhwZOPwIcPks\n/VOtX5I0RJbNNaCq3gDem2Ql8EdJbp6xv5LUwpZ1P8ePT7Fjxw6OHTvGxMTEwp5ekkZYr9ej1+ud\nk3Onav6v50n+HfBz4F8AE1V1tC0NPV1V1yTZClBVD7bxe4H7gB+2Mde2/ruBm6rqo7Nco6BYuXKc\nycntjI+Pn+UUJemtLQlVlblHzm2udx+958Q7i5IsB34VeA7YA2xqwzYBj7X2HmBjkvOTrAHGgH1V\ndRR4Ncn69uD5noFjJElDYq7lo0uBnUneQT9AvlhVTyV5DtidZDNwGLgToKoOJNkNHABeB7bUyVuR\nLcAjwHLgiarau9CTkSSdndNaPloMLh9J0ulZtOUjSdLbi6EgSeoYCpKkjqEgSeoYCpKkjqEgSeoY\nCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKk\nzpyhkOTKJE8neSHJ95Pc2/ovTjKZ5MUkTyZZNXDMtiQvJTmY5NaB/nVJ9rd9D52bKUmSztR87hRe\nA/51Vf0S8D7gN5NcC2wFJqvqauCptk2StcBdwFpgA7A9yYk/KP0wsLmqxoCxJBsWdDaSpLMyZyhU\n1dGqer61/wr4AXA5cBuwsw3bCdzR2rcDu6rqtao6DBwC1ie5FFhRVfvauEcHjpEkDYHTeqaQ5Crg\neuAZYHVVTbdd08Dq1r4MODJw2BH6ITKzf6r1S5KGxLL5DkzybuAPgI9V1V+eXBGCqqoktXBl3c/x\n41Ps2LGDY8eOMTExsXCnlqQR1+v16PV65+TcqZr7tTzJO4H/AvxhVX2m9R0EJqrqaFsaerqqrkmy\nFaCqHmzj9gL3AT9sY65t/XcDN1XVR2dcq6BYuXKcycntjI+PL9hkJemtKAlVlblHzm0+7z4K8Dng\nwIlAaPYAm1p7E/DYQP/GJOcnWQOMAfuq6ijwapL17Zz3DBwjSRoC81k+ej/w68D3kjzX+rYBDwK7\nk2wGDgN3AlTVgSS7gQPA68CWOnk7sgV4BFgOPFFVexdoHpKkBTBnKFTVt3jzO4pb3uSYB4AHZul/\nFrjudAqUJC0eP9EsSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEg\nSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkzpyhkOTzSaaT7B/ouzjJZJIXkzyZZNXA\nvm1JXkpyMMmtA/3rkuxv+x5a+KlIks7WfO4UvgBsmNG3FZisqquBp9o2SdYCdwFr2zHbk6Qd8zCw\nuarGgLEkM88pSVpic4ZCVf0x8Bczum8Ddrb2TuCO1r4d2FVVr1XVYeAQsD7JpcCKqtrXxj06cIwk\naUic6TOF1VU13drTwOrWvgw4MjDuCHD5LP1TrV+SNETO+kFzVRVQC1CLJGmJLTvD46aTXFJVR9vS\n0Mutfwq4cmDcFfTvEKZae7B/6s1Pfz/Hj0+xY8cOjh07xsTExBmWKUlvPb1ej16vd07Onf4P+nMM\nSq4CHq+q69r2p4GfVtWnkmwFVlXV1vag+cvAjfSXh74B/EJVVZJngHuBfcDXgc9W1d5ZrlVQrFw5\nzuTkdsbHxxdkopL0VpWEqsrcI+c2551Ckl3ATcB7kvwI+PfAg8DuJJuBw8CdAFV1IMlu4ADwOrCl\nTqbOFuARYDnwxGyBIElaWvO6U1hM3ilI0ulZyDsFP9EsSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKk\njqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEg\nSeoseigk2ZDkYJKXknxisa8vSXpzixoKSc4D/iOwAVgL3J3k2jcbf+ONN5Kk+xp2vV5vqUs4K9a/\ntKx/6Yxy7Qttse8UbgQOVdXhqnoN+D3g9lMfUu1r+I36N5b1Ly3rXzqjXPtCW+xQuBz40cD2kdY3\nL6N01yBJo2jZIl9vXj/yX3TRh/n5z188xeHzD4aq/jFnGiQnjh8081yzjZGkUZTFfEFL8j7g/qra\n0La3AW9U1acGxvgKK0mnqaoWZAllsUNhGfCnwD8B/hzYB9xdVT9YtCIkSW9qUZePqur1JL8F/BFw\nHvA5A0GShsei3ilIkobb0HyieVg/1Jbk80mmk+wf6Ls4yWSSF5M8mWTVwL5tbQ4Hk9w60L8uyf62\n76FFqv3KJE8neSHJ95PcO2L1vyvJM0meT3IgyW+PUv0D1z4vyXNJHh+1+pMcTvK9Vv++Uao/yaok\nX0nyg/b9s36Eav/F9m9+4uuVJPcuSv1VteRf9JeSDgFXAe8EngeuXeq6Wm0fBK4H9g/0fRr4t639\nCeDB1l7ban9nm8shTt6N7QNubO0ngA2LUPslwHtb+930n+dcOyr1t2td0P67DPjvwAdGqf52vX8D\nfAnYM0rfP+1afwZcPKNvJOoHdgL/fOD7Z+Wo1D5jHu8AfgxcuRj1L9rE5pj0rwB7B7a3AluXuq6B\neq7ib4bCQWB1a18CHGztbcAnBsbtBd4HXAr8YKB/I/CflmAejwG3jGL9wAXAt4FfGqX6gSuAbwA3\nA4+P2vcP/VD42zP6hr5++gHwP2fpH/raZ6n5VuCPF6v+YVk+OqsPtS2B1VU13drTwOrWvox+7Sec\nmMfM/ikWeX5JrqJ/x/MMI1R/knckeb7V+XRVvcAI1Q/8DvBx4I2BvlGqv4BvJPlOko+0vlGofw3w\nkyRfSPLdJL+b5EJGo/aZNgK7Wvuc1z8soTCyT7urH79DXX+SdwN/AHysqv5ycN+w119Vb1TVe+n/\nxP2Pktw8Y//Q1p/k14CXq+o5YNb3kA9z/c37q+p64EPAbyb54ODOIa5/GXADsL2qbgCO0V+B6Axx\n7Z0k5wMfBn5/5r5zVf+whMIU/fWyE67kb6bbsJlOcglAkkuBl1v/zHlcQX8eU6092D+1CHWS5J30\nA+GLVfVY6x6Z+k+oqleArwPrGJ36/yFwW5I/o/+T3j9O8kVGp36q6sftvz8Bvkb/95eNQv1HgCNV\n9e22/RX6IXF0BGof9CHg2fbvD4vwbz8sofAdYCzJVS0Z7wL2LHFNp7IH2NTam+iv1Z/o35jk/CRr\ngDFgX1UdBV5t734IcM/AMedMu9bngANV9ZkRrP89J95dkWQ58KvAc6NSf1V9sqqurKo19JcAvllV\n94xK/UkuSLKitS+kv7a9fxTqb9f8UZKrW9ctwAvA48Ne+wx3c3Lp6ESd57b+xXxgMsfDlA/Rf3fM\nIWDbUtczUNcu+p++/r/0n3v8BnAx/YeHLwJPAqsGxn+yzeEg8E8H+tfR/x/qEPDZRar9A/TXsp+n\n/2L6HP1fWz4q9V8HfLfV/z3g461/JOqfMZebOPnuo5Gon/66/PPt6/sn/r8cofr/Pv03J/wJ8FX6\nD59HovZ23QuB/w2sGOg75/X74TVJUmdYlo8kSUPAUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLH\nUJAkdf4fQ1xF+tlkFYgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb849123390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFhJJREFUeJzt3X9sXtd93/H3J1bkKokiwcgmyz8Aay29WIO7JOqodE1g\nZnM9tWht/2XLQAyh0wIEWpdkAdJKATb7L88JUDQKDBsYm8RykCjT0tSQV1el4voBOnSI8sNuVTOq\npG1KImaigySNmzLCJOS7P3h09YSj9ZASRdHi+wUQPvd7z73PPcfS8+E996GYqkKSJIDXXekLkCQt\nHYaCJKljKEiSOoaCJKljKEiSOoaCJKkzMBSS7EryUpLDST6f5Nok1yU5mORokrEka2f0P5bkSJK7\n+uqb2jmOJdl9uQYkSbp4FwyFJLcA7wPeUVW3A9cAW4GdwMGquhV4rm2TZCNwP7AR2AI8niTtdE8A\n26tqCBhKsmXBRyNJuiSD7hReAc4Ab0iyAngD8F3gbmBP67MHuLe17wH2VtWZqjoBHAc2J1kPrK6q\nQ63fU33HSJKWiAuGQlX9APg94NtMh8HfVtVBYF1VTbZuk8C61r4BONl3ipPAjbPUJ1pdkrSEDFo+\n+nngQ8AtTL+xvynJe/v71PS/k+G/lSFJV4EVA/b/EvAXVfV9gCRfAn4ZOJXk+qo61ZaGXm79J4Cb\n+46/iek7hInW7q9PzPaCSQwYSZqnqsrgXoMNeqZwBHhnklXtgfGdwDjwDLCt9dkGPN3a+4GtSVYm\n2QAMAYeq6hTwSpLN7TwP9h3z/6mqZf/10EMPXfFrWCpfzoVz4Vxc+GshXfBOoar+MslTwNeAnwLf\nAP4zsBrYl2Q7cAK4r/UfT7KP6eA4C+yo81e8A3gSWAU8W1UHFnQkkqRLNmj5iKr6OPDxGeUfMH3X\nMFv/R4BHZql/Hbj9Iq5RkrRI/InmJWpkZORKX8KS4Vyc51yc51xcHlno9ahLlaSW2jVJ0lKWhFqk\nB82SpGXEUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwF\nSVLHUJAkdQwFSVLHUJAkdQwFSVJnYCgk+cdJXuj7+lGSDyS5LsnBJEeTjCVZ23fMriTHkhxJcldf\nfVOSw23f7ld7zQ996CN8+MMf4dvf/valj1CSNGcDQ6Gq/qaq3l5Vbwc2AVPAHwE7gYNVdSvwXNsm\nyUbgfmAjsAV4PMm5XxP3BLC9qoaAoSRbZnvN3bv/IU888QVOnTp1aaOTJM3LfJeP7gSOV9V3gLuB\nPa2+B7i3te8B9lbVmao6ARwHNidZD6yuqkOt31N9x8zwEa69dv08L02SdKnmGwpbgb2tva6qJlt7\nEljX2jcAJ/uOOQncOEt9otUlSUvEnEMhyUrgN4H/OnNfVRVQC3hdkqQrYMU8+v4a8PWq+l7bnkxy\nfVWdaktDL7f6BHBz33E3MX2HMNHa/fWJ2V/qYU6fnmB0dJSpqSlGRkbmcZmSdHXr9Xr0er3Lcu5M\nf5M/h47JF4A/qao9bfvjwPer6mNJdgJrq2pne9D8eWCY6eWhLwO/UFWV5CvAB4BDwB8Dn6yqAzNe\np6BYs2aYsbHHGB4eXqChStLVKQlVlcE9B5vTnUKSNzL9kPl9feVHgX1JtgMngPsAqmo8yT5gHDgL\n7KjzybMDeBJYBTw7MxAkSVfWnEKhqv4eeMuM2g+YDorZ+j8CPDJL/evA7fO/TEnSYvAnmiVJHUNB\nktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQx\nFCRJHUNBktQxFCRJHUNBktQxFCRJnTmFQpK1Sb6Y5JtJxpNsTnJdkoNJjiYZS7K2r/+uJMeSHEly\nV199U5LDbd/uyzEgSdLFm+udwm7g2aq6DfhF4AiwEzhYVbcCz7VtkmwE7gc2AluAx5OknecJYHtV\nDQFDSbYs2EgkSZdsYCgkWQO8u6o+DVBVZ6vqR8DdwJ7WbQ9wb2vfA+ytqjNVdQI4DmxOsh5YXVWH\nWr+n+o6RJC0Bc7lT2AB8L8lnknwjyWiSNwLrqmqy9ZkE1rX2DcDJvuNPAjfOUp9odUnSErFijn3e\nAfx2VX01ySdoS0XnVFUlqYW7rIc5fXqC0dFRpqamGBkZWbhTS9JrXK/Xo9frXZZzp+rC7+VJrgf+\nR1VtaNvvAnYB/wh4T1WdaktDz1fVW5PsBKiqR1v/A8BDwLdan9ta/QHgjqp6/4zXKyjWrBlmbOwx\nhoeHF3K8knTVSUJVZXDPwQYuH1XVKeA7SW5tpTuBl4BngG2ttg14urX3A1uTrEyyARgCDrXzvNI+\nuRTgwb5jJElLwFyWjwD+HfC5JCuB/wn8FnANsC/JduAEcB9AVY0n2QeMA2eBHXX+dmQH8CSwiulP\nMx1YoHFIkhbAwOWjxebykSTNz6IuH0mSlg9DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLU\nMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSZ05hUKSE0n+\nKskLSQ612nVJDiY5mmQsydq+/ruSHEtyJMldffVNSQ63fbsXfjiSpEsx1zuFAkaq6u1VNdxqO4GD\nVXUr8FzbJslG4H5gI7AFeDzJuV8o/QSwvaqGgKEkWxZoHJKkBTCf5aPM2L4b2NPae4B7W/seYG9V\nnamqE8BxYHOS9cDqqjrU+j3Vd4wkaQmYz53Cl5N8Lcn7Wm1dVU229iSwrrVvAE72HXsSuHGW+kSr\nS5KWiBVz7PcrVfV/kvwD4GCSI/07q6qS1MJd1sOcPj3B6OgoU1NTjIyMLNypJek1rtfr0ev1Lsu5\nUzW/9/IkDwE/Bt7H9HOGU21p6PmqemuSnQBV9WjrfwB4CPhW63Nbqz8A3FFV759x/oJizZphxsYe\nY3h4GEnSq0tCVc1c4r8oA5ePkrwhyerWfiNwF3AY2A9sa922AU+39n5ga5KVSTYAQ8ChqjoFvJJk\nc3vw/GDfMZKkJWAuy0frgD9qHyBaAXyuqsaSfA3Yl2Q7cAK4D6CqxpPsA8aBs8COOn87sgN4ElgF\nPFtVBxZwLJKkSzTv5aPLzeUjSZqfRV0+kiQtH4aCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaC\nJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOnMKhSTX\nJHkhyTNt+7okB5McTTKWZG1f311JjiU5kuSuvvqmJIfbvt0LPxRJ0qWa653CB4FxoNr2TuBgVd0K\nPNe2SbIRuB/YCGwBHk9y7pdJPwFsr6ohYCjJloUZgiRpoQwMhSQ3Ab8O/AFw7g3+bmBPa+8B7m3t\ne4C9VXWmqk4Ax4HNSdYDq6vqUOv3VN8xkqQlYi53Cr8PfAT4aV9tXVVNtvYksK61bwBO9vU7Cdw4\nS32i1SVJS8iKC+1M8hvAy1X1QpKR2fpUVSWp2fZdvIc5fXqC0dFRpqamGBmZ9aUlaVnq9Xr0er3L\ncu5Uvfr7eZJHgAeBs8DPAW8GvgT8M2Ckqk61paHnq+qtSXYCVNWj7fgDwEPAt1qf21r9AeCOqnr/\nLK9ZUKxZM8zY2GMMDw8v4HAl6eqThKrK4J6DXXD5qKo+WlU3V9UGYCvwZ1X1ILAf2Na6bQOebu39\nwNYkK5NsAIaAQ1V1Cnglyeb24PnBvmMkSUvEBZePZnHutuJRYF+S7cAJ4D6AqhpPso/pTyqdBXbU\n+VuRHcCTwCrg2ao6cGmXLklaaBdcProSXD6SpPlZtOUjSdLyYihIkjqGgiSpYyhIkjqGgiSpYyhI\nkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqG\ngiSpc8FQSPJzSb6S5MUk40n+U6tfl+RgkqNJxpKs7TtmV5JjSY4kuauvvinJ4bZv9+UbkiTpYl0w\nFKrqNPCeqnob8IvAe5K8C9gJHKyqW4Hn2jZJNgL3AxuBLcDjSc79MukngO1VNQQMJdlyOQYkSbp4\nA5ePqmqqNVcC1wA/BO4G9rT6HuDe1r4H2FtVZ6rqBHAc2JxkPbC6qg61fk/1HSNJWiIGhkKS1yV5\nEZgEnq+ql4B1VTXZukwC61r7BuBk3+EngRtnqU+0uiRpCVkxqENV/RR4W5I1wJ8mec+M/ZWkFvay\nHub06QlGR0eZmppiZGRkYU8vSa9hvV6PXq93Wc6dqrm/nyf5D8BPgH8DjFTVqbY09HxVvTXJToCq\nerT1PwA8BHyr9bmt1R8A7qiq98/yGgXFmjXDjI09xvDw8CUOUZKubkmoqgzuOdigTx+95dwni5Ks\nAn4VeAHYD2xr3bYBT7f2fmBrkpVJNgBDwKGqOgW8kmRze/D8YN8xkqQlYtDy0XpgT5LXMR0gn62q\n55K8AOxLsh04AdwHUFXjSfYB48BZYEedvxXZATwJrAKeraoDCz0YSdKlmdfy0WJw+UiS5mfRlo8k\nScuLoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSO\noSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqTOwFBIcnOS55O8lOSvk3yg1a9LcjDJ0SRjSdb2HbMr\nybEkR5Lc1VfflORw27f78gxJknSx5nKncAb491X1T4B3Av82yW3ATuBgVd0KPNe2SbIRuB/YCGwB\nHk9y7hdKPwFsr6ohYCjJlgUdjSTpkgwMhao6VVUvtvaPgW8CNwJ3A3tatz3Ava19D7C3qs5U1Qng\nOLA5yXpgdVUdav2e6jtGkrQEzOuZQpJbgLcDXwHWVdVk2zUJrGvtG4CTfYedZDpEZtYnWl2StESs\nmGvHJG8C/hD4YFX93fkVIaiqSlILd1kPc/r0BKOjo0xNTTEyMrJwp5ak17her0ev17ss507V4Pfy\nJK8H/hvwJ1X1iVY7AoxU1am2NPR8Vb01yU6Aqnq09TsAPAR8q/W5rdUfAO6oqvfPeK2CYs2aYcbG\nHmN4eHjBBitJV6MkVFUG9xxsLp8+CvApYPxcIDT7gW2tvQ14uq++NcnKJBuAIeBQVZ0CXkmyuZ3z\nwb5jJElLwFyWj34FeC/wV0leaLVdwKPAviTbgRPAfQBVNZ5kHzAOnAV21PnbkR3Ak8Aq4NmqOrBA\n45AkLYCBoVBV/51Xv6O481WOeQR4ZJb614Hb53OBkqTF4080S5I6hoIkqWMoSJI6hoIkqWMoSJI6\nhoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIk\nqTMwFJJ8OslkksN9teuSHExyNMlYkrV9+3YlOZbkSJK7+uqbkhxu+3Yv/FAkSZdqLncKnwG2zKjt\nBA5W1a3Ac22bJBuB+4GN7ZjHk6Qd8wSwvaqGgKEkM88pSbrCBoZCVf058MMZ5buBPa29B7i3te8B\n9lbVmao6ARwHNidZD6yuqkOt31N9x0iSloiLfaawrqomW3sSWNfaNwAn+/qdBG6cpT7R6pKkJeSS\nHzRXVQG1ANciSbrCVlzkcZNJrq+qU21p6OVWnwBu7ut3E9N3CBOt3V+fePXTP8zp0xOMjo4yNTXF\nyMjIRV6mJF19er0evV7vspw709/oD+iU3AI8U1W3t+2PA9+vqo8l2Qmsraqd7UHz54FhppeHvgz8\nQlVVkq8AHwAOAX8MfLKqDszyWgXFmjXDjI09xvDw8IIMVJKuVkmoqgzuOdjAO4Uke4E7gLck+Q7w\nH4FHgX1JtgMngPsAqmo8yT5gHDgL7KjzqbMDeBJYBTw7WyBIkq6sgaFQVQ+8yq47X6X/I8Ajs9S/\nDtw+r6uTJC0qf6JZktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNB\nktQxFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktQZ+JvXrqTNmzf/zPZcfp+0JOnivQbuFKp9SZIu\nt0UPhSRbkhxJcizJ7y7260uSXt2ihkKSa4DHgC3ARuCBJLfN4/ju62rX6/Wu9CUsGc7Fec7Fec7F\n5bHYdwrDwPGqOlFVZ4AvAPfM/fDls5TkH/jznIvznIvznIvLY7FD4UbgO33bJ1tt3vrvGpbTHYQk\nXU6L/emjOX2b/+Y3/yY/+cnROZ4qP9O+1GDo/4TTxZzrQsf76SlJS10W840qyTuBh6tqS9veBfy0\nqj7W18d3Tkmap6pakKWSxQ6FFcDfAP8S+C5wCHigqr65aBchSXpVi7p8VFVnk/w28KfANcCnDARJ\nWjoW9U5BkrS0LZmfaF5uP9SW5OYkzyd5KclfJ/lAq1+X5GCSo0nGkqztO2ZXm58jSe66cle/8JJc\nk+SFJM+07WU5DwBJ1ib5YpJvJhlPsnk5zkcb10tJDif5fJJrl9M8JPl0kskkh/tq8x5/kk1tDo8l\n2T3whavqin8xvZR0HLgFeD3wInDblb6uyzzm64G3tfabmH7WchvwceB3Wv13gUdbe2Obl9e3eToO\nvO5Kj2MB5+PDwOeA/W17Wc5DG+Me4F+39gpgzXKbjzaW/wVc27b/C7BtOc0D8G7g7cDhvtp8xn9u\nJegQMNzazwJbLvS6S+VO4RJ/qO21p6pOVdWLrf1j4JtM/8zG3Uy/KdD+e29r3wPsraozVXWC6f/p\nw4t60ZdJkpuAXwf+gOnPGMMynAeAJGuAd1fVp2H6OVxV/YjlNx+vAGeAN7QPqLyB6Q+nLJt5qKo/\nB344ozyf8W9Osh5YXVWHWr+n+o6Z1VIJhQX7obbXoiS3MP0dwVeAdVU12XZNAuta+wam5+Wcq2mO\nfh/4CPDTvtpynAeADcD3knwmyTeSjCZ5I8tsPqrqB8DvAd9mOgz+tqoOsszmYRbzHf/M+gQD5mWp\nhMKyfdqd5E3AHwIfrKq/699X0/d7F5qb1/y8JfkN4OWqeoHzdwk/YznMQ58VwDuAx6vqHcDfAzv7\nOyyH+Ujy88CHmF4KuQF4U5L39vdZDvNwIXMY/0VZKqEwAdzct30zP5tuV6Ukr2c6ED5bVU+38mSS\n69v+9cDLrT5zjm5qtde6fw7cneR/A3uBf5Hksyy/eTjnJHCyqr7atr/IdEicWmbz8UvAX1TV96vq\nLPAl4JdZfvMw03z+Xpxs9Ztm1C84L0slFL4GDCW5JclK4H5g/xW+pssq0/8GxqeA8ar6RN+u/Uw/\nUKP99+m++tYkK5NsAIaYfoD0mlZVH62qm6tqA7AV+LOqepBlNg/nVNUp4DtJbm2lO4GXgGdYXvNx\nBHhnklXt78qdwDjLbx5mmtffi/bn6ZX2CbYAD/YdM7sr/YS976n6rzH9CZzjwK4rfT2LMN53Mb2G\n/iLwQvvaAlwHfBk4CowBa/uO+WibnyPAv7rSY7gMc3IH5z99tJzn4Z8CXwX+kunvkNcsx/kAfofp\nQDzM9EPV1y+neWD6zvm7wP9l+pnrb13M+IFNbQ6PA58c9Lr+8JokqbNUlo8kSUuAoSBJ6hgKkqSO\noSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6vw/qYCC2+xLCD4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb823a72dd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGbNJREFUeJzt3X+MXeV95/H3x7/wTzBOyPjnBm81JLhJF+IEJwGWS5c4\nbkWARhEYqZabWFFaNw27aUnsaCXGqpTSpqvGVQVblSSYNPGumzSWUcCxIdxu1G4yCdjBZnBt05i1\nB2YMwcHmh8OM57t/nMeHk+nAjGfO3HNn5vOSrua5z33OOd9zx57vfb7PufcqIjAzMwOYVHUAZmbW\nPJwUzMws56RgZmY5JwUzM8s5KZiZWc5JwczMcoMmBUnvkLSncHtR0mckzZO0W9JBSbskzS1ss1HS\nIUkHJK0s9C+XtC89tnm0TsrMzIZH5/I+BUmTgE7gCuCPgOcj4i8kfR64MCI2SFoGfBN4H7AIeAho\njYiQ1A58OiLaJT0A/HVE7Cz5nMzMbJjOtXx0HXA4Io4CNwBbUv8W4KbUvhHYGhE9EXEEOAyskLQA\nmBMR7WncfYVtzMysCZxrUlgNbE3tlojoTu1uoCW1FwLHCtscI5sx9O/vTP1mZtYkhpwUJE0DPgL8\nQ//HIqtB+fMyzMzGuCnnMPa3gEcj4rl0v1vS/IjoSqWh46m/E1hS2G4x2QyhM7WL/Z39DyLJycXM\n7BxFhMrYz7mUj27l9dIRwA5gbWqvBbYX+ldLmiZpKdAKtEdEF3BS0gpJAtYUtvkVEdFUtzvuuKPy\nGBzT+ImpWeNyTGM3pjINaaYgaRbZIvMnC913AtskrQOOADenP+gdkrYBHUAvsD5ej3o9cC8wA3gg\nfOWRmVlTGVJSiIiXgbf263uBLFEMNP6LwBcH6H8UePe5h2lmZo3gdzQPQa1WqzqEf8cxDU0zxgTN\nGZdjGppmjKlM5/TmtUaQFM0Wk5lZM5NEVLDQbGZm45yTgpmZ5ZwUzMws56RgZlaST30KzpypOoqR\n8UKzmVlJJk2Cnh6YPLmxx/VCs5lZE4rIEsNYNsbDNzNrDmcLHCrl9Xp1nBTMzErQ1zf2ZwngpGBm\nVoq+vrE/SwAnBTOzUoyH9QRwUjAzK4XLR2ZmlnP5yMzMci4fmZlZzuUjMzPLOSmYmVnOawpmZpbz\nmoKZmeVcPjIzs5yTgpmZ5bymYGZmuQm1piBprqRvSXpSUoekFZLmSdot6aCkXZLmFsZvlHRI0gFJ\nKwv9yyXtS49tHo0TMjOrwkQrH20GHoiIS4HfAA4AG4DdEXEJ8HC6j6RlwC3AMmAVcJeUT6ruBtZF\nRCvQKmlVaWdiZlahCZMUJF0AXB0RXwWIiN6IeBG4AdiShm0BbkrtG4GtEdETEUeAw8AKSQuAORHR\nnsbdV9jGzGxMm0hrCkuB5yR9TdJjkv5O0iygJSK605huoCW1FwLHCtsfAxYN0N+Z+s3MxrzxsqYw\nZYhj3gN8OiJ+LOnLpFLRWRERkqKsoNra2vJ2rVajVquVtWszs1HRyPJRvV6nXq+Pyr4V8eZ/yyXN\nB/5vRCxN968CNgL/Ebg2IrpSaeiRiHinpA0AEXFnGr8TuAN4Oo25NPXfClwTEb/f73gxWExmZs3m\n8GFYtSr72WiSiIhSileD5rWI6AKOSrokdV0HPAHcD6xNfWuB7am9A1gtaZqkpUAr0J72czJduSRg\nTWEbM7MxbbysKQylfATwR8A3JE0DngI+DkwGtklaBxwBbgaIiA5J24AOoBdYX3jpvx64F5hBdjXT\nzpLOw8ysUuNlTWHQ8lGjuXxkZmPRk0/CRz+a/Wy0hpaPzMxscBPmfQpmZja48bKm4KRgZlaCvj6Y\nPLnqKEbOScHMrAQuH5mZWe7MGScFMzNLXD4yM7Ocy0dmZpZzUjAzs5zXFMzMLOc1BTMzy7l8ZGZm\nOZePzMws5/KRmZnlXD4yM7Ock4KZmeW8pmBmZjmvKZiZWc7lIzMzy7l8ZGZmOZePzMws5/KRmZnl\nXD4yM7PchCofSToi6XFJeyS1p755knZLOihpl6S5hfEbJR2SdEDSykL/ckn70mObyz8dM7NqTLTy\nUQC1iLg8Iq5IfRuA3RFxCfBwuo+kZcAtwDJgFXCXJKVt7gbWRUQr0CppVUnnYWZWqYmWFADU7/4N\nwJbU3gLclNo3AlsjoicijgCHgRWSFgBzIqI9jbuvsI2Z2Zg20dYUAnhI0k8kfTL1tUREd2p3Ay2p\nvRA4Vtj2GLBogP7O1G9mNuaNlzWFKUMcd2VEPCvpImC3pAPFByMiJEVZQbW1teXtWq1GrVYra9dm\nZqOikeWjer1OvV4flX0PKSlExLPp53OSvgNcAXRLmh8RXak0dDwN7wSWFDZfTDZD6EztYn/nQMcr\nJgUzs7GgkeWj/i+WN23aVNq+Bz0FSTMlzUntWcBKYB+wA1ibhq0Ftqf2DmC1pGmSlgKtQHtEdAEn\nJa1IC89rCtuYmY1pE6l81AJ8J11ANAX4RkTskvQTYJukdcAR4GaAiOiQtA3oAHqB9RFxtrS0HrgX\nmAE8EBE7SzwXM7PKjJeFZr3+97o5SIpmi8nMbDBf+hJ0d8Nf/mXjjy2JiOh/heiwjIO8ZmZWvd5e\nmDq16ihGzknBzKwEPT0wZajXczYxJwUzsxL09HimYGZmSW+vZwpmZpZ4pmBmZjkvNJuZWc4LzWZm\nlnP5yMzMci4fmZlZzuUjMzPLuXxkZmY5v0/BzMxynimYmVnOC81mZpbzQrOZmeVcPjIzs5wXms3M\nLHf6NEyfXnUUI+ekYGZWgldfhRkzqo5i5JwUzMxKcPq0k4KZmSWvvurykZmZJV5TMDOz3IQqH0ma\nLGmPpPvT/XmSdks6KGmXpLmFsRslHZJ0QNLKQv9ySfvSY5vLPxUzs+qcPg3nnVd1FCM31JnCbUAH\nEOn+BmB3RFwCPJzuI2kZcAuwDFgF3CVJaZu7gXUR0Qq0SlpVzimYmVUrAvr6YPLkqiMZuUGTgqTF\nwG8D9wBn/8DfAGxJ7S3ATal9I7A1Inoi4ghwGFghaQEwJyLa07j7CtuYmY1pZ87ApEmQvwQew4Yy\nU/gr4Hagr9DXEhHdqd0NtKT2QuBYYdwxYNEA/Z2p38xszDtzZnzMEgDe9E3Zkq4HjkfEHkm1gcZE\nREiKgR4brra2trxdq9Wo1QY8tJlZU2h0UqjX69Tr9VHZtyLe+O+5pC8Ca4BeYDpwPvCPwPuAWkR0\npdLQIxHxTkkbACLizrT9TuAO4Ok05tLUfytwTUT8/gDHjDeLycys2Zw6BQsWwEsvVXN8SUREKcWr\nNy0fRcQXImJJRCwFVgPfj4g1wA5gbRq2Ftie2juA1ZKmSVoKtALtEdEFnJS0Ii08rylsY2Y2pk2Y\n8tEAzr6EvxPYJmkdcAS4GSAiOiRtI7tSqRdYX3jZvx64F5gBPBARO0cWuplZcxhPSeFNy0dVcPnI\nzMaa556DSy+F55+v5vgNKx+ZmdngxtNMwUnBzGyEnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxy\nTgpmZpZzUjAzs5yTgpmZ5ZwUzMws56RgZmY5JwUzM8s5KZiZWc5JwczMcidOwOzZVUdRDicFM7MR\nevVVOP/8qqMoh5OCmdkI9fbC1KlVR1EOJwUzsxHq6YEp5/rlxk3KScHMbIQ8UzAzs5xnCmZmlvNM\nwczMcp4pmJlZbsLMFCRNl/QjSXsldUj6s9Q/T9JuSQcl7ZI0t7DNRkmHJB2QtLLQv1zSvvTY5tE7\nJTOzxpowM4WIOA1cGxGXAb8BXCvpKmADsDsiLgEeTveRtAy4BVgGrALukqS0u7uBdRHRCrRKWjUa\nJ2Rm1mgTZqYAEBGvpOY0YDJwArgB2JL6twA3pfaNwNaI6ImII8BhYIWkBcCciGhP4+4rbGNmNqZN\nmJkCgKRJkvYC3cAjEfEE0BIR3WlIN9CS2guBY4XNjwGLBujvTP1mZmPeeJopDJrbIqIPuEzSBcD3\nJF3b7/GQFGUG1dbWlrdrtRq1Wq3M3ZuZlarRM4V6vU69Xh+VfQ/5NCLiRUnfBZYD3ZLmR0RXKg0d\nT8M6gSWFzRaTzRA6U7vY3/lGxyomBTOzZtfomUL/F8ubNm0qbd+DXX301rNXFkmaAXwI2APsANam\nYWuB7am9A1gtaZqkpUAr0B4RXcBJSSvSwvOawjZmZmPaeFpTGOw0FgBbJE0iSyBfj4iHJe0Btkla\nBxwBbgaIiA5J24AOoBdYHxFnS0vrgXuBGcADEbGz7JMxM6tCT8/4WVPQ63+zm4OkaLaYzMzezNq1\ncO218Hu/V83xJRERGnzk4PyOZjOzERpPMwUnBTOzEertHT9rCk4KZmYj5JmCmZnlPFMwM7Pca6/B\ntGlVR1EOJwUzsxHq6XFSMDOz5LXXvKZgZmaJZwpmZpbzTMHMzHKeKZiZWe7ll2HGjKqjKIeTgpnZ\nCB0/DhddVHUU5fAH4pmZjcCZM9l6Qm8vTKroZbY/EM/MrEm8/DLMmlVdQijbODkNM7NqvPRSlhTG\nCycFM7MReOUVmD276ijK46RgZjYCr7wCM2dWHUV5nBTMzEbglVfGz+Wo4KRgZjYip087KZiZWfLL\nX8J551UdRXmcFMzMRuD0aZg+veooyuOkYGY2Ap4pmJlZ7vTpCZYUJC2R9IikJyTtl/SZ1D9P0m5J\nByXtkjS3sM1GSYckHZC0stC/XNK+9Njm0TklM7PG+eUvJ175qAf4bxHx68D7gT+UdCmwAdgdEZcA\nD6f7SFoG3AIsA1YBd0k6+5kcdwPrIqIVaJW0qtSzMTNrsAk3U4iIrojYm9ovAU8Ci4AbgC1p2Bbg\nptS+EdgaET0RcQQ4DKyQtACYExHtadx9hW3MzMakiThTyEm6GLgc+BHQEhHd6aFuoCW1FwLHCpsd\nI0si/fs7U7+Z2Zg13mYKU4Y6UNJs4NvAbRFx6vWKEERESCrt867b2trydq1Wo1arlbVrM7NSVXH1\nUb1ep16vj8q+h5QUJE0lSwhfj4jtqbtb0vyI6EqloeOpvxNYUth8MdkMoTO1i/2dAx2vmBTMzJrZ\n6dON/0C8/i+WN23aVNq+h3L1kYCvAB0R8eXCQzuAtam9Fthe6F8taZqkpUAr0B4RXcBJSSvSPtcU\ntjEzG5NOnYI5c6qOojxDmSlcCfwu8LikPalvI3AnsE3SOuAIcDNARHRI2gZ0AL3A+sJXqa0H7gVm\nAA9ExM6SzsPMrBInTsCFF1YdRXn8dZxmZiOwciV89rOwqsIL7P11nGZmTeKFF2DevKqjKI+TgpnZ\nCBw/Dm97W9VRlMflIzOzYYrIvkvhxIlqv1PB5SMzsyZw8iRMm+Yv2TEzM6C7e3yVjsBJwcxs2Lq7\noaVl8HFjiZOCmdkwHT/upGBmZonLR2ZmlnP5yMzMck4KZmaWe/ppWLiw6ijK5aRgZjZMzzwDb397\n1VGUy0nBzGyYnn0WFiyoOopy+WMuzMyGoacHZs7MvmRn8uRqY/HHXJiZVayrCy66qPqEUDYnBTOz\nYRiPpSNwUjAzGxYnBTMzyzkpmJlZzknBzMxyTgpmZpYbj5+QCk4KZmbDcvIkXHBB1VGUz0nBzGwY\nTp6E88+vOoryDZoUJH1VUrekfYW+eZJ2SzooaZekuYXHNko6JOmApJWF/uWS9qXHNpd/KmZmjXPy\nJMyZU3UU5RvKTOFrwKp+fRuA3RFxCfBwuo+kZcAtwLK0zV2Szr71+m5gXUS0Aq2S+u/TzGzM+MUv\nJmj5KCJ+AJzo130DsCW1twA3pfaNwNaI6ImII8BhYIWkBcCciGhP4+4rbGNmNqb09sKJE9nHXIw3\nw11TaImI7tTuBs6uwS8EjhXGHQMWDdDfmfrNzMaco0ezK4+mTKk6kvKN+JQiIiSV+rGmbW1tebtW\nq1Gr1crcvZnZiDz/PMyfX93x6/U69Xp9VPY93KTQLWl+RHSl0tDx1N8JLCmMW0w2Q+hM7WJ/5xvt\nvJgUzMyazdGj1b5xrf+L5U2bNpW27+GWj3YAa1N7LbC90L9a0jRJS4FWoD0iuoCTklakhec1hW3M\nzMaU738fPvCBqqMYHUO5JHUr8C/AOyQdlfRx4E7gQ5IOAr+Z7hMRHcA2oAN4EFhf+Mac9cA9wCHg\ncETsLPtkzMwa4fHH4YMfrDqK0eFvXjMzO0fvehds3QrvfnfVkWT8zWtmZhX6+c/hLW+pOorR4aRg\nZnYOIrKkMG9e1ZGMDicFM7Nz8MILMH16dhuPnBTMzM7BI4/AVVdVHcXocVIwMzsH7e3j93JUcFIw\nMzsn9Tpcc03VUYweX5JqZjZEfX3ZdygcOwZz5w4+vlF8SaqZWQWeeSb7DoVmSghlc1IwMxuivXth\n2bKqoxhdTgpmZkP04IPw4Q9XHcXo8pqCmdkQRGQfl/3QQ83z8RZneU3BzKzBjh6Fl17KPvdoPHNS\nMDMbgv374corQaW8Hm9eTgpmZkOwf//4nyWAk4KZ2ZDU63D55VVHMfq80GxmNog9e+D66+Gpp5rz\ng/C80Gxm1kB/+qdw++3NmRDK5pmCmdmbeOyxbJZw6BDMmlV1NAPzTMHMrAEi4I//GO64o3kTQtmc\nFMzM3sDf/z08/zx84hNVR9I4U6oOwMysGT34YDZLePBBmDq16mgaxzMFM7N+/vZv4dZb4dvfhuXL\nq46msRqeFCStknRA0iFJn2/08c3M3kh7O3zkI/ClL8EPfgBXX111RI3X0KQgaTLwN8AqYBlwq6RL\nGxnDcNTr9apD+Hcc09A0Y0zQnHFN5JgefRT+4A/gd34n+1a1/fvf+EPvmvF5KlOjZwpXAIcj4khE\n9AD/C7ixwTGcs2b8R+CYhqYZY4LmjGuixNTbC/v2wT33wLp18M53wsc+ln1xzk9/Cn/yJ2/+foRm\nfJ7K1OiF5kXA0cL9Y8CKBsdgZmNQBJw5k/1R7+2Fl1+G1157/fbKK3DyJJw6BS++mN1/4QV47jk4\nfhy6urJPOj16FJYsgfe+F666Cj71KXjf+8b/B90NVaOTwpDelXb99YPsZIjvbStr3KFD8MMfNvaY\ng4156in4539u7DEHG/ezn8E//VNjjznYuKefhocfbuwxhzLm6FH43vcae8zBxnV2wne/O/L99fXB\n6dPZ2LO3vr7B2wM9duoU3HlnlgTOnIHJk7PblCkwc2b2in7atOw2Y0b2/cmzZ8OFF2b3586FRYvg\nssuy70JYvBje/vaJ856D4WjoO5olvR9oi4hV6f5GoC8i/rwwxm9nNjM7R2W9o7nRSWEK8K/AfwGe\nAdqBWyPiyYYFYWZmb6ih5aOI6JX0aeB7wGTgK04IZmbNo+k+EM/MzKrTNO9obuSb2iR9VVK3pH2F\nvnmSdks6KGmXpLmFxzamuA5IWlnoXy5pX3ps8whjWiLpEUlPSNov6TNVxyVpuqQfSdorqUPSn1Ud\nU2F/kyXtkXR/E8V0RNLjKa72ZohL0lxJ35L0ZPodrqj439Q70vNz9vaipM80wfO0Mf3f2yfpm5LO\na4KYbkv72i/pttQ3+jFFROU3slLSYeBiYCqwF7h0FI93NXA5sK/Q9xfA51L788Cdqb0sxTM1xXeY\n12dY7cAVqf0AsGoEMc0HLkvt2WRrL5c2QVwz088pwA+Bq6qOKe3js8A3gB3N8PtL+/gZMK9fX9W/\nvy3AJwq/wwuqjqkQ2yTgWWBJlTGl/f4bcF66/7+BtRXH9C5gHzCd7O/jbuDXGhHTiH6pZd2ADwA7\nC/c3ABtG+ZgX86tJ4QDQktrzgQOpvRH4fGHcTuD9wALgyUL/auB/lhjfduC6ZokLmAn8GPj1qmMC\nFgMPAdcC9zfL748sKbylX19lcZElgH8boL/y5yrtZyXwg6pjAuaRvQi7kCxx3g98qOKYPgbcU7j/\n34HPNSKmZikfDfSmtkUNjqElIrpTuxtoSe2FKZ6zzsbWv7+TkmKWdDHZTOZHVcclaZKkvenYj0TE\nE1XHBPwVcDvQV+irOibI3ofzkKSfSPpkE8S1FHhO0tckPSbp7yTNqjimotXA1tSuLKaIeAH4H8D/\nI7sq8hcRsbvKmID9wNWpXDQT+G2yF0OjHlOzJIWmWu2OLKVWEpOk2cC3gdsi4lTVcUVEX0RcRvYP\n8j9LurbKmCRdDxyPiD3AgNdlV/j7uzIiLgd+C/hDSb/ycWoVxDUFeA9wV0S8B3iZbBZeZUwASJoG\nfAT4h/6PVfBv6teA/0pWPVgIzJb0u1XGFBEHgD8HdgEPkpWGzjQipmZJCp1kdcWzlvCr2a0RuiXN\nB5C0ADj+BrEtTrF1pnaxv3MkAUiaSpYQvh4R25slLoCIeBH4LrC84pg+CNwg6WdkrzJ/U9LXK44J\ngIh4Nv18DvgO2Wd9VRnXMeBYRPw43f8WWZLoqvq5Ikucj6bnCqp9nt4L/EtE/DwieoF/JCtpV/o8\nRcRXI+K9EXENcAI4SAOep2ZJCj8BWiVdnF5B3ALsaHAMO8gWl0g/txf6V0uaJmkp0Aq0R0QXcDJd\nzSFgTWGbc5b28RWgIyK+3AxxSXrr2asbJM0gq7PuqTKmiPhCRCyJiKVk5YfvR8SaKmMCkDRT0pzU\nnkVWL99XZVxpX0clXZK6rgOeIKuZV/ZcJbfyeuno7LGriukA8H5JM9K+rgM6qPh5kvS29PM/AB8F\nvkkjnqfhLIKMxo3slcO/kq2abxzlY20lqx2+RraW8XGyxaaHyLLxLmBuYfwXUlwHgA8X+peT/cc/\nDPz1CGO6iqxGvpfsD+8eso8Yrywu4N3AYymmx4HbU3+lz1Vhn9fw+tVHVf/+lqbnaS9ZPXhjk8T1\nn8guEPgp2SvgC5ogplnA88CcQl/VMX2OLGHuI7tia2oTxPR/Ukx7gWsb9Tz5zWtmZpZrlvKRmZk1\nAScFMzPLOSmYmVnOScHMzHJOCmZmlnNSMDOznJOCmZnlnBTMzCz3/wFsJUIB2cVNtAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb846884750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# take k mean of country sample\n",
    "uniq_countries, uniq_counts = np.unique(Y, return_counts=True)\n",
    "k = np.int(np.round(np.mean(uniq_counts)))\n",
    "#N_k = n_occurrence_from_D(D, k=k)\n",
    "\n",
    "print k\n",
    "print skew(N_k)\n",
    "print np.median(N_k)\n",
    "print np.mean(N_k)\n",
    "print np.std(N_k)\n",
    "print len(np.where(N_k>1000)[0])\n",
    "print skew(N_k[N_k<1000])\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(N_k, bins=100);\n",
    "plt.figure()\n",
    "plt.hist(N_k[N_k<1000], bins=100);\n",
    "plt.figure()\n",
    "plt.plot(np.sort(N_k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## what if we use a different distance metric?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.20652171267\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFadJREFUeJzt3W2MXNd93/HvT5IJWXFdVlBBUw8oCVSCzcIAhcZEUOdh\ngiSqFLSS/MaSgwREqwYGlNqGi6QmDbRcJ0CsBLDrF4X9xnLAuA0bIoEFKoFlUaomcF6YgltSkrVm\nJQJmYbrSKm2VVIIqmKr+fTF3yeFqd2d2HnZ35n4/wGDPnvsw5+wl/3Pmf8+9N1WFJGm+XbXVDZAk\nTZ/BXpJawGAvSS1gsJekFjDYS1ILGOwlqQWGCvZJrk5yOsmjze/XJzmZ5IUkjyfZ2bfu4SQvJjmb\n5I5pNVySNLxhR/afAhaB5Un5h4CTVXUb8GTzO0n2AfcB+4A7gS8n8duDJG2xgYE4yc3ALwNfBdJU\n3w0cbcpHgXub8j3Asaq6WFXngXPAgUk2WJK0ccOMuv8d8FvA2311u6pqqSkvAbua8o3Ahb71LgA3\njdtISdJ41g32Sf4J8EpVnebyqP4K1bvfwnr3XPB+DJK0xa4ZsPwfAXcn+WXgWuC9Sb4OLCV5X1W9\nnGQ38Eqz/o+AW/q2v7mpu0ISPwAkaQRVterAe5B1R/ZV9dmquqWq9gL3A/+5qn4NOAEcbFY7CDzS\nlE8A9yfZkWQvcCvw9Br7ntvXkSNHtrwN9s/+ta1vbejfOAaN7N8Ro5ufDwHHkzwAnAc+2gTwxSTH\n6c3ceQt4sMZtoSRpbEMH+6r6C+AvmvL/Bn5xjfV+F/jdibROkjQRzoGfgk6ns9VNmCr7N7vmuW8w\n//0bR7Yiy5LE7I4kbVASahonaCVJ88FgL0ktsNHZOK2RXPlNybSTpFnmyH5dgy4OlqTZYLCXpBYw\n2EtSCxjsJakFDPaS1AIGe0lqAYO9JLWA8+w3yPn3kmaRI/uROP9e0mwx2EtSCxjsJakFDPaS1AIG\ne0lqAYO9JLXAusE+ybVJTiU5k2Qxyeeb+oUkF5Kcbl539W1zOMmLSc4muWPaHZAkDTbwsYRJrquq\nN5JcA/wl8JvALwCvVdUXV6y7D/gj4EPATcATwG1V9faK9bb9Ywl78+mX25hL8+nXqpekaZvqYwmr\n6o2muAO4Gnh1+X1XWf0e4FhVXayq88A54MAoDZMkTc7AYJ/kqiRngCXgqap6vln0iSTPJHk4yc6m\n7kbgQt/mF+iN8CVJW2iYkf3bVbUfuBn42SQd4CvAXmA/8BLwhfV2MYF2SpLGMPS9carqb5L8OfCT\nVdVdrk/yVeDR5tcfAbf0bXZzU/cOCwsLl8qdTodOpzNsUySpFbrdLt1udyL7WvcEbZIbgLeq6q+T\nvBv4FvA54PmqerlZ59PAh6rqV/pO0B7g8gnav7/ybKwnaCVp48Y5QTtoZL8bOJrkKnopn69X1ZNJ\n/jDJfnpR7wfAxwGqajHJcWAReAt4cNtHdUlqgYFTL6fypo7sJWnDpjr1UpI0+wz2ktQCPqlqinyq\nlaTtwpH91PlUK0lbz5H9BDiCl7TdObKfGEfwkrYvg70ktYDBXpJawGAvSS3gCdpN1H8i15O4kjaT\nI/tN5UlcSVvDYC9JLWCwl6QWMNhLUgsY7CWpBQz2ktQCrZp66T1sJLVVC0f2Tn+U1D4tDPaS1D7r\nBvsk1yY5leRMksUkn2/qr09yMskLSR5PsrNvm8NJXkxyNskd0+6AJGmwgQ8cT3JdVb2R5BrgL4Hf\nBO4G/mdV/X6SzwB/p6oOJdkH/BHwIeAm4Angtqp6e8U+t+SB4xt5WPhGHjg+7Lo+qFzSOKb6wPGq\neqMp7gCuBl6lF+yPNvVHgXub8j3Asaq6WFXngXPAgVEaJkmanIHBPslVSc4AS8BTVfU8sKuqlppV\nloBdTflG4ELf5hfojfAlSVto4NTLJgWzP8nfBr6V5OdXLK8k6+UkVl22sLBwqdzpdOh0OsO0d+44\nHVTSWrrdLt1udyL7Gpizv2Ll5N8A/xf4F0Cnql5OspveiP/9SQ4BVNVDzfqPAUeq6tSK/WybnP1a\nwXazcvYbOY8gqd2mlrNPcsPyTJsk7wZ+CTgNnAAONqsdBB5pyieA+5PsSLIXuBV4epSGba7tNfc+\nyRUvSRrXoDTObuBokqvofTB8vaqeTHIaOJ7kAeA88FGAqlpMchxYBN4CHtySIfxc6P9GIEnj2VAa\nZ2Jvuu3SOMOnVzYjjWNqR9Jqpjr1UpI0+wz2ktQCBntJagGDvSS1gMFeklrAYC9JLWCwl6QWMNhL\nUgsY7CWpBQz2ktQCBntJagGDvSS1gMFeklrAYC9JLWCwl6QWMNhLUgsY7CWpBQY9llDbyFoPR5ek\nQRzZz5zt9XB0SbNhYLBPckuSp5I8n+R7ST7Z1C8kuZDkdPO6q2+bw0leTHI2yR3T7MAgSS69NrK+\nJM2TYdI4F4FPV9WZJO8B/kuSk/SGl1+sqi/2r5xkH3AfsA+4CXgiyW1V9faE274B/Q/9nuS6kjQb\nBo7sq+rlqjrTlF8Hvk8viMPqEfEe4FhVXayq88A54MBkmrv9+C1A0izYUM4+yR7gduA7TdUnkjyT\n5OEkO5u6G4ELfZtd4PKHwxwyfy5p+xt6Nk6TwvkT4FNV9XqSrwC/3Sz+HeALwANrbP6OiLiwsHCp\n3Ol06HQ6wzZFklqh2+3S7XYnsq8MM30vybuAPwO+WVVfWmX5HuDRqvpgkkMAVfVQs+wx4EhVnepb\nvzZr2mAvzdKfh79crqp1l69cd+39bWTd9ddZrU3rrSupPZJQVSPljoeZjRPgYWCxP9An2d232keA\n55ryCeD+JDuS7AVuBZ4epXGSpMkYJo3zYeBXgWeTnG7qPgt8LMl+ekPNHwAfB6iqxSTHgUXgLeDB\nTRvGS5JWNVQaZ+JvahrHNI6kDRsnjePtEqbA6ZiSthtvlzAV3tJA0vZisJekFjDYS1ILGOwlqQUM\n9pLUAgZ7SWoBg70ktYDBXpJawGAvSS1gsJekFmh1sPe2BpLaotXB3lsaSGqLlgd7SWoHg70ktYDB\nXpJawGAvSS3gw0u2yLRmAvXv1ydZSVpmsB/S5INz/+MKJ235UYmS1DMwjZPkliRPJXk+yfeSfLKp\nvz7JySQvJHk8yc6+bQ4neTHJ2SR3TLMDm8enT0maXcPk7C8Cn66qfwD8FPAbST4AHAJOVtVtwJPN\n7yTZB9wH7APuBL6cxHMDkrSFBgbhqnq5qs405deB7wM3AXcDR5vVjgL3NuV7gGNVdbGqzgPngAMT\nbrckaQM2NOJOsge4HTgF7KqqpWbRErCrKd8IXOjb7AK9DwdJ0hYZ+gRtkvcAfwp8qqpeWzHro5Ks\nl9B+x7KFhYVL5U6nQ6fTGbYpktQK3W6Xbrc7kX1lmOl5Sd4F/Bnwzar6UlN3FuhU1ctJdgNPVdX7\nkxwCqKqHmvUeA45U1am+/dVmTQvsfSj1z3xZWR60fJjyJPbRK1fVO9q8/Ldaq371/q6+XNLsSkJV\njTTVbpjZOAEeBhaXA33jBHCwKR8EHumrvz/JjiR7gVuBp0dpnCRpMoZJ43wY+FXg2SSnm7rDwEPA\n8SQPAOeBjwJU1WKS48Ai8Bbw4KYN4yVJqxoqjTPxNzWNs2bZNI6ktUw1jSNJmn0Ge0lqAe+NMwd8\nvKKkQRzZzw3v3SNpbXMzsl85uvXkpCRdNmcje0e3krSaOQv2kqTVGOwlqQUM9pLUAnNzgnaeObVS\n0rgc2c8ETzxLGo/BXpJawGAvSS1gsJekFjDYS1ILGOwlqQUM9pLUAgZ7SWoBg70ktcDAYJ/ka0mW\nkjzXV7eQ5EKS083rrr5lh5O8mORskjum1fB55dWykqZhmJH9HwB3rqgr4ItVdXvz+iZAkn3AfcC+\nZpsvJ/Hbw4ZM/krZJFe8JLXPwEBcVd8GXl1l0WpR4x7gWFVdrKrzwDngwFgt1IR4ywWpzcYZdX8i\nyTNJHk6ys6m7EbjQt84F4KYx3mNkjmIl6bJR73r5FeC3m/LvAF8AHlhj3VWHkwsLC5fKnU6HTqcz\nYlPWsvy2BnxJs6nb7dLtdieyrwzzrNYke4BHq+qD6y1Lcgigqh5qlj0GHKmqUyu2qUk/I7Y3iu8P\n8MOWR91u0vsYZX/93rmP5b/xyr+Nz+eVZlMSqmqkEexIaZwku/t+/QiwPFPnBHB/kh1J9gK3Ak+P\n8h4ahnl4ScMZmMZJcgz4OeCGJD8EjgCdJPvpRZofAB8HqKrFJMeBReAt4MGJD+ElSRs2VBpn4m9q\nGmdT2mQaR5ovm57GkSTNFoO9JLWAwV6SWsBgL0ktYLCXpBYY9QpazYGVt5Nwlo40vxzZt54XZklt\nYLCXpBYw2EtSCxjsJakFDPaS1AIGe0lqAYO9JLWAwV6SWsBgL0kt4BW0c2zUB657Za00fxzZz7Vx\nro71ylppnhjsJakFDPaS1AIDg32SryVZSvJcX931SU4meSHJ40l29i07nOTFJGeT3DGthkuShjfM\nyP4PgDtX1B0CTlbVbcCTze8k2QfcB+xrtvlykol+e0hyxUsb599Oap+Bgbiqvg28uqL6buBoUz4K\n3NuU7wGOVdXFqjoPnAMOTKapV7QKTx6Ow7+f1Dajjrp3VdVSU14CdjXlG4ELfetdAG4a8T0kSRMy\n9jz7qqok6w0TV122sLBwqdzpdOh0OuM2RZLmSrfbpdvtTmRfGeaCmSR7gEer6oPN72eBTlW9nGQ3\n8FRVvT/JIYCqeqhZ7zHgSFWdWrG/GvVCnV6ueXnbXLrgZ2X98OVRt5v0Pra+TWv9Ldc7Vl6AJW2e\nJFTVSCfcRk3jnAAONuWDwCN99fcn2ZFkL3Ar8PSI76GZ4TkAabsbmMZJcgz4OeCGJD8E/i3wEHA8\nyQPAeeCjAFW1mOQ4sAi8BTw48hBekjQxQ6VxJv6mpnG2aZv6bSSNM9y6ksazFWkczSXTMdK8mvm7\nXnpx0ObwRKw02+ZgZO9odPP4t5Zm1RwEe0nSIAZ7SWoBg70ktYDBXpJawGCvgZzxJM2+mZ96qc2w\nfJHWYP0fDE7PlLYPR/aaMKdnStuRwV6SWsBgL0ktYM5eU7daHt/bL0iby5G9Nslqwdz8vrRZDPaS\n1AIGe0lqAYO9JLWAwV6SWmCs2ThJzgP/B/h/wMWqOpDkeuCPgb9H83zaqvrrMdupGeRtFqTtY9yR\nfQGdqrq9qg40dYeAk1V1G/Bk87taydk20nYxiTTOyuHb3cDRpnwUuHcC7yFJGsMkRvZPJPlukl9v\n6nZV1VJTXgJ2jfkeQC8lYFpAkkYz7hW0H66ql5L8XeBkkrP9C6uqkkzwe/zwd1+UJF02VrCvqpea\nn3+V5BvAAWApyfuq6uUku4FXVtt2YWHhUrnT6dDpdMZpiiTNnW63S7fbnci+Muo9SZJcB1xdVa8l\n+QngceBzwC8C/6uqfi/JIWBnVR1asW1t9H17KZzlkf3ytpMob5d9zFabrrzHzfD7WGs7740jDZaE\nqhopvTHOyH4X8I0mj34N8B+r6vEk3wWOJ3mAZurlGO8hSZqAkUf2Y72pI/s5aFO/yY3svRumtLZx\nRvZeQasRTXMOvfPzpUnb1vezd6qlJE3GDIzsHeVJ0rhmINhLksa1rdM4mj+m5qSt4chem8y0nLQV\nDPaS1AIGe0lqAXP2mgmDLrbyYixpfQZ7bWtXBvHLV9yufqJ3rSt819+vHwxqA9M42ubWOqFbA5YP\nu2+pHRzZa1uY1pRM0ztSjyN7bRPea0eapm0xsnf0JUnTtS2Cfc+gk2/S8Pw3JF1pGwX7fsPPqpBW\nt/F/Q8N8w/RbqGaVOXvpCsPk9z0HoNmzTUf20nSMkt4ZZhvn7Wu727Jgv7i4uFVvrVYbNUW4+nbv\nvOjL1KO2p6kE+yR3Al8Crga+WlW/t3Kd/ft/kmuv3cPFi69OownS2Ib7FrD2h8da+X3z/toKE8/Z\nJ7ka+PfAncA+4GNJPrByveuu+xCvvbbIm28uTLoJ20B3qxugsXSbn5PIza93BfDmB/lut7vp77mZ\n5r1/45jGCdoDwLmqOl9VF4H/BNwzhffZxrpb3QBtQLJyum93q5oCXG7PynaNU79s3oPhvPdvHNMI\n9jcBP+z7/UJTJ21Tkx1lv/PDY7j6K5dv9BvBoHsI9Xzuc59b9UNA828aOfuh/te8+eZzvPe9/5Qf\n//i/8+abU2iFtGXWyuMPU//Ok7xrBeWNBOu17h46zL6qaqQPhvXOUUzifIbXRWxMJt35JD8FLFTV\nnc3vh4G3+0/SJmnvX1ySxlBVI30lm0awvwb4b8AvAP8DeBr4WFV9f6JvJEka2sTTOFX1VpJ/CXyL\n3tTLhw30krS1Jj6ylyRtP5t+b5wkdyY5m+TFJJ/Z7PeftCTnkzyb5HSSp5u665OcTPJCkseT7Nzq\ndg4rydeSLCV5rq9uzf4kOdwcy7NJ7tiaVg9vjf4tJLnQHMPTSe7qWzZr/bslyVNJnk/yvSSfbOpn\n/hiu07e5OH5Jrk1yKsmZJItJPt/UT+bYVdWmveildc4Be4B3AWeAD2xmG6bQpx8A16+o+33gXzfl\nzwAPbXU7N9CfnwFuB54b1B96F82daY7lnubYXrXVfRihf0eAf7XKurPYv/cB+5vye+idP/vAPBzD\ndfo2T8fvuubnNcB3gJ+e1LHb7JH9vF5wtfLs+N3A0aZ8FLh3c5szuqr6NrDyHhZr9ece4FhVXayq\n8/T+sR3YjHaOao3+weo3tZnF/r1cVWea8uvA9+ld5zLzx3CdvsH8HL83muIOeoPjV5nQsdvsYD+P\nF1wV8ESS7yb59aZuV1UtNeUlYNfWNG1i1urPjfSO4bJZPp6fSPJMkof7vibPdP+S7KH3LeYUc3YM\n+/r2naZqLo5fkquSnKF3jJ6qqueZ0LHb7GA/j2eDP1xVtwN3Ab+R5Gf6F1bv+9bc9HuI/sxiX78C\n7AX2Ay8BX1hn3ZnoX5L3AH8KfKqqXutfNuvHsOnbn9Dr2+vM0fGrqreraj9wM/CzSX5+xfKRj91m\nB/sfAbf0/X4LV34yzZyqeqn5+VfAN+h9jVpK8j6AJLuBV7auhROxVn9WHs+bm7qZUlWvVAP4Kpe/\nCs9k/5K8i16g/3pVPdJUz8Ux7Ovbf1ju27wdP4Cq+hvgz4F/yISO3WYH++8CtybZk2QHcB9wYpPb\nMDFJrkvyt5ryTwB3AM/R69PBZrWDwCOr72FmrNWfE8D9SXYk2QvcSu8iupnS/Ada9hF6xxBmsH/p\n3R/gYWCxqr7Ut2jmj+FafZuX45fkhuUUVJJ3A78EnGZSx24LzjbfRe8s+jng8Fad9Z5QX/bSOxt+\nBvjecn+A64EngBeAx4GdW93WDfTpGL0rn39M7/zKP1uvP8Bnm2N5FvjHW93+Efr3z4E/BJ4Fnmn+\nI+2a4f79NPB282/ydPO6cx6O4Rp9u2tejh/wQeC/Nv17Fvitpn4ix86LqiSpBXzguCS1gMFeklrA\nYC9JLWCwl6QWMNhLUgsY7CWpBQz2ktQCBntJaoH/Dzx1pB+a1db3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7616b7dd90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_cos = pairwise_distances(X, metric='cosine')\n",
    "uniq_countries, uniq_counts = np.unique(Y, return_counts=True)\n",
    "k = np.int(np.round(np.mean(uniq_counts)))\n",
    "N_k = n_occurrence_from_D(D_cos, k=k)\n",
    "print skew(N_k)\n",
    "plt.figure()\n",
    "plt.hist(N_k, bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sort_idx = np.argsort(D, axis=1)\n",
    "k = 100\n",
    "D_k = sort_idx[:, 1:(k+1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4650, 2942, 3520, ..., 3488, 2864, 6361],\n",
       "       [1933, 6143, 6757, ..., 2346, 3441, 6857],\n",
       "       [3170, 2549, 4860, ..., 2260, 2978, 7433],\n",
       "       ..., \n",
       "       [6016, 2243, 1616, ..., 3486, 3441, 4554],\n",
       "       [7027, 4860, 6346, ..., 7312, 6972, 5532],\n",
       "       [5119, 1563, 4035, ..., 6253, 7433, 3232]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 129  186  261  373  432  468  509  515  552  562  627  655  701  710  740\n",
      "  743  861  913  997 1012 1064 1128 1204 1273 1422 1426 1563 1566 1622 1761\n",
      " 1846 1933 1936 1949 1980 1989 2034 2042 2053 2110 2147 2329 2335 2346 2462\n",
      " 2546 2549 2568 2651 2864 2903 3003 3044 3069 3170 3215 3239 3295 3383 3397\n",
      " 3441 3486 3488 3517 3520 3686 3712 3731 3859 3892 3982 4026 4035 4167 4208\n",
      " 4390 4421 4456 4549 4554 4571 4639 4650 4782 4822 4859 4860 4869 4940 5020\n",
      " 5042 5091 5119 5146 5189 5352 5432 5458 5493 5532 5707 5873 5894 5914 5926\n",
      " 5961 6056 6123 6143 6160 6242 6253 6346 6356 6439 6476 6487 6678 6819 6857\n",
      " 6859 6876 7027 7051 7053 7091 7312 7333 7361 7464 7478 7617 7627 7678 7783\n",
      " 7849 7884 7888 7891 7920 8035 8106]\n",
      "['India' 'Kazakhstan' 'Lesotho' 'Somalia' 'Kazakhstan' 'Kazakhstan'\n",
      " 'Kazakhstan' 'Nigeria' 'South Sudan' 'South Sudan' 'Nigeria' 'Kazakhstan'\n",
      " 'Nigeria' 'Kazakhstan' 'Kazakhstan' 'Kazakhstan' 'United Kingdom'\n",
      " 'South Sudan' 'Uganda' 'Solomon Islands' 'United States of America'\n",
      " 'Lesotho' 'Australia' 'Solomon Islands' 'Fiji' 'Solomon Islands'\n",
      " 'Swaziland' 'Kazakhstan' 'Zambia' 'Kiribati' 'Kazakhstan' 'Kazakhstan'\n",
      " 'Kyrgyzstan' 'Kazakhstan' 'Iraq' 'Sierra Leone' 'Kazakhstan' 'Kyrgyzstan'\n",
      " 'Swaziland' 'United Kingdom' 'Kazakhstan' 'Ireland' 'Kazakhstan'\n",
      " 'Sierra Leone' 'Australia' 'Lesotho' 'Swaziland' 'Western Sahara' 'Ghana'\n",
      " 'Kazakhstan' 'South Sudan' 'Ireland' 'Solomon Islands' 'Kazakhstan'\n",
      " 'Belgium' 'Swaziland' 'United Kingdom' 'Iraq' 'Kenya' 'Sierra Leone'\n",
      " 'Solomon Islands' 'Kazakhstan' 'Kenya' 'Democratic Republic of the Congo'\n",
      " 'Swaziland' 'Lesotho' 'Western Sahara' 'Somalia' 'Somalia'\n",
      " 'United Kingdom' 'United Kingdom' 'Nigeria' 'Kyrgyzstan' 'Kazakhstan'\n",
      " 'Kazakhstan' 'United Kingdom' 'Australia' 'Kazakhstan' 'United Kingdom'\n",
      " 'Nigeria' 'Kazakhstan' 'Kazakhstan' 'Solomon Islands' 'Sierra Leone'\n",
      " 'Somalia' 'Kazakhstan' 'United Kingdom' 'Botswana' 'China' 'Swaziland'\n",
      " 'New Zealand' 'Kazakhstan' 'Pakistan' 'Kazakhstan' 'Swaziland'\n",
      " 'Solomon Islands' 'Kyrgyzstan' 'Kyrgyzstan' 'Uganda' 'Solomon Islands'\n",
      " 'Solomon Islands' 'Kazakhstan' 'Australia' 'Kazakhstan' 'Sierra Leone'\n",
      " 'South Sudan' 'Zambia' 'Western Sahara' 'Nigeria' 'Somalia' 'Nepal'\n",
      " 'South Africa' 'Uganda' 'Nigeria' 'Kazakhstan' 'Swaziland' 'Iraq'\n",
      " 'Australia' 'Sierra Leone' 'Pakistan' 'Swaziland' 'Solomon Islands'\n",
      " 'Uganda' 'Ireland' 'Kazakhstan' 'Nigeria' 'Kazakhstan' 'Kazakhstan'\n",
      " 'Kazakhstan' 'Kazakhstan' 'Somalia' 'Ireland' 'South Sudan' 'Nigeria'\n",
      " 'Somalia' 'Solomon Islands' 'Kazakhstan' 'Australia' 'Iraq' 'Nigeria'\n",
      " 'Solomon Islands' 'Ghana']\n",
      "['D:/_Audio/Music-from-India/025A-RKDATX0070XX-0300A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-3200A0.mp3'\n",
      " 'D:/_Audio/Kevin-Volans-South-Africa/025A-C0740X008X01-0700A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0021XX-1000A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0500A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1600A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0700A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043663XX-0100A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0109XX-0001A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0156XX-0001A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043641XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-2500A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043646XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-2200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0084XX-3400A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-2600A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0085XX-0600A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0117XX-0001A0.mp3'\n",
      " 'D:/_Audio/Peter-Cooke-Uganda/025A-C0023X0006XX-0200A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0011XX-1300A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0309XX-2000A0.mp3'\n",
      " 'D:/_Audio/Kevin-Volans-South-Africa/025A-C0740X004X04-0100A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X1025XX-0100A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0011XX-1800A0.mp3'\n",
      " 'D:/_Audio/GK-Roth-Fiji/025A-C0028X0006XX-0200A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0034XX-0100A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-0300A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1900A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0006XX-0400A0.mp3'\n",
      " 'D:/_Audio/Simon-Seligman-Kiribati/025A-C0205X0004XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-2000A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0087XX-2200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0084XX-3800A0.mp3'\n",
      " 'D:/_Audio/Middle-Eastern-music-on-shellac-disc/025A-9CS0025122XX-0002A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043740XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1900A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0087XX-1600A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0004XX-3700A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0338XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1400A0.mp3'\n",
      " 'D:/_Audio/Terry-Yarnell-Collection/025A-C1005X0002XX-1200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1700A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043758XX-0100A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X1076XX-0100A0.mp3'\n",
      " 'D:/_Audio/Kevin-Volans-South-Africa/025A-C0740X004X24-0200A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-2000A0.mp3'\n",
      " 'D:/_Audio/Violeta-Ruano-portraits-of-Saharawi-Music/025A-C1640X0028XX-0001A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043498XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1700A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0145XX-0001A0.mp3'\n",
      " 'D:/_Audio/Terry-Yarnell-Collection/025A-C1005X0006XX-2900A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0012XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-3400A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0097XX-0700A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-1800A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0254XX-2500A0.mp3'\n",
      " 'D:/_Audio/Middle-Eastern-music-on-shellac-disc/025A-9CS0026163XX-0002A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0003XX-1500A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043717XX-0200A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0001XX-0500A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-3100A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0003XX-1700A0.mp3'\n",
      " 'D:/_Audio/Ken-Gourlay-Uganda/025A-C0105X0005XX-0400A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-1100A0.mp3'\n",
      " 'D:/_Audio/Kevin-Volans-South-Africa/025A-C0740X004X06-0500A0.mp3'\n",
      " 'D:/_Audio/Violeta-Ruano-portraits-of-Saharawi-Music/025A-C1640X0025XX-0001A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0012XX-0300A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0024XX-0800A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0264XX-3300A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0323XX-0900A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043639XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0087XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-2400A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1500A0.mp3'\n",
      " 'D:/_Audio/Roy-Palmer-collection/025A-C1023X0006XX-0700A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X1080XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1500A0.mp3'\n",
      " 'D:/_Audio/Reg-Hall-Archive/025A-C0903X0094XX-0200A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043662XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0300A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0030XX-0100A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043712XX-0200A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0017XX-0800A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-3000A0.mp3'\n",
      " 'D:/_Audio/Bob-and-Jacqueline-Patten-Collection/025A-C1033X0281XX-3200A0.mp3'\n",
      " 'D:/_Audio/John-Brearley-Botswana/025A-C0065X0043XX-0900A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0081XX-1500A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-1300A0.mp3'\n",
      " '/import/c4dm-scratch/matthiasm/data/smithsonian/12932/26829/SampleAudio.wav'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0400A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0031XX-0500A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-2400A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-0500A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0032XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0084XX-0800A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0088XX-0400A0.mp3'\n",
      " 'D:/_Audio/Ken-Gourlay-Uganda/025A-C0105X0021XX-1200A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0011XX-0800A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0011XX-1700A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1400A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X1020XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-3000A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043759XX-0200A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0143XX-0001A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0006XX-0600A0.mp3'\n",
      " 'D:/_Audio/Violeta-Ruano-portraits-of-Saharawi-Music/025A-C1640X0047XX-0001A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043665XX-0100A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0020XX-0800A0.mp3'\n",
      " 'D:/_Audio/Carol-Tingey-Nepal-Collection/025A-C1465X0072XX-1400A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0018XX-0200A0.mp3'\n",
      " 'D:/_Audio/Peter-Cooke-Uganda/025A-C0023X0002XX-3300A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043644XX-0100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-0900A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0041XX-0800A0.mp3'\n",
      " 'D:/_Audio/Middle-Eastern-music-on-shellac-disc/025A-9CS0025068XX-0002A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X1071XX-0100A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043769XX-0100A0.mp3'\n",
      " 'D:/_Audio/Nicholas-Pierce-Balochistan/025A-C0140X0012XX-0200A0.mp3'\n",
      " 'D:/_Audio/Rycroft/025A-C0811X0005XX-1500A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0003XX-0200A0.mp3'\n",
      " 'D:/_Audio/Peter-Cooke-Uganda/025A-C0023X0036XX-0700A0.mp3'\n",
      " 'D:/_Audio/Terry-Yarnell-Collection/025A-C1005X0002XX-1000A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-1300A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043583XX-0200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-2500A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1200A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0085XX-2900A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0086XX-1000A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0012XX-1100A0.mp3'\n",
      " 'D:/_Audio/Terry-Yarnell-Collection/025A-C1005X0002XX-0200A0.mp3'\n",
      " 'D:/_Audio/Dinka-songs-from-South-Sudan/025A-C1580X0148XX-0001A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043691XX-0100A0.mp3'\n",
      " 'D:/_Audio/John-Low-East-Africa/025A-C0027X0012XX-0200A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0011XX-1100A0.mp3'\n",
      " 'D:/_Audio/Colin-Huehns-Pakistan/025A-C0485X0084XX-3700A0.mp3'\n",
      " 'D:/_Audio/Ethnographic-wax-cylinders/025A-C0080X0816XX-0100A0.mp3'\n",
      " 'D:/_Audio/Middle-Eastern-music-on-shellac-disc/025A-9CS0025050XX-0001A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043654XX-0100A0.mp3'\n",
      " 'D:/_Audio/Firth-Tikopia/025A-C0333X0018XX-0200A0.mp3'\n",
      " 'D:/_Audio/Decca-West-African-recordings/025A-1CS0043520XX-0100A0.mp3']\n"
     ]
    }
   ],
   "source": [
    "large_hubs_idx = np.where(N_k>1000)[0]\n",
    "print large_hubs_idx\n",
    "print Y[large_hubs_idx]\n",
    "print Yaudio[large_hubs_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7160,)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where(D_k==515)[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## let's get the audio url to listen to tracks identified as large hubs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homes/mp305/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2822: DtypeWarning: Columns (0,1,2,4,5,6,7,8,10,11,12,13,14,15,16,17,19,21,22,23,24,25,26,27,29,31,35,38,39,40,41,44,45,48,55,56,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,93,95,96) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  if self.run_code(code, result):\n"
     ]
    }
   ],
   "source": [
    "ddf = outliers.load_metadata(Yaudio, metadata_file='../data/metadata_BLSM_language_all.csv')\n",
    "#ddf = outliers.load_metadata(Yaudio, metadata_file='/Users/mariapanteli/Documents/'+\n",
    "#                             'QMUL/Code/MyPythonCode/MergeBL-Smith/data/metadata_BLSM_language_all.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8200, 108)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddf.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"songurls_Album\" in ddf.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### first, fix the url for BL tracks (because it was changed recently and the metadata.csv file is not updated) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mariapanteli/anaconda/lib/python2.7/site-packages/pandas/core/indexing.py:121: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  self._setitem_with_indexer(indexer, value)\n"
     ]
    }
   ],
   "source": [
    "bl_inds = np.where(np.isnan(ddf['BuyLinkTrackDownload']))[0]\n",
    "for bl_ind in bl_inds:\n",
    "    ddf['songurls_Album'].iloc[bl_ind] = ('https://sounds.bl.uk/World-and-traditional-music/' + \n",
    "                                        ddf['Folder'].iloc[bl_ind] + '/' + \n",
    "                                        ddf['MetaFile'].iloc[bl_ind].split('.')[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 'https://sounds.bl.uk/World-and-traditional-music/Decca-West-African-recordings/025M-1CS0043663XX-0100V0',\n",
       "       'https://sounds.bl.uk/World-and-traditional-music/Rycroft/025M-C0811X0005XX-2000V0',\n",
       "       'https://sounds.bl.uk/World-and-traditional-music/Colin-Huehns-Pakistan/025M-C0485X0085XX-3100V0',\n",
       "       'https://sounds.bl.uk/World-and-traditional-music/Rycroft/025M-C0811X0005XX-1300V0',\n",
       "       'https://sounds.bl.uk/World-and-traditional-music/Colin-Huehns-Pakistan/025M-C0485X0031XX-0500V0'], dtype=object)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "large_hubs_idx = np.array([515, 2549, 3486, 5020, 5119])\n",
    "ddf['songurls_Album'].iloc[large_hubs_idx].get_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}