beta-matrix-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deserialize_matrix(path):\n",
    "    with open(path, 'r') as f:\n",
    "        lines = f.readlines()\n",
    "        matrix = []\n",
    "        for l in lines:\n",
    "            l = l.strip()\n",
    "            l = l.replace(' ', '')\n",
    "            tokens = l.split(',')\n",
    "            row = [float(el) for el in tokens if el != \"\" and el != \" \"]\n",
    "            matrix.append(row)\n",
    "    return np.asarray(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deserialize_id_map(path):\n",
    "    with open(path, 'r') as f:\n",
    "        lines = f.readlines()\n",
    "        idmap = dict()\n",
    "        for l in lines:\n",
    "            tokens = l.split(\",\")\n",
    "            col_id = int(tokens[0])\n",
    "            col_name = tokens[1].strip()\n",
    "            idmap[col_id] = col_name\n",
    "    return idmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform_matrix_dataframe(matrix):\n",
    "    idx = [i for i in range(matrix.shape[0])]\n",
    "    matrix_df = pd.DataFrame(data=matrix, index=idx, columns=idx)\n",
    "    return matrix_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def name_pair(pair, idmap):\n",
    "    i, j = pair\n",
    "    return idmap[i] + \" ~ \" + idmap[j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_matrix = \"temp/mit_matrix_n2_k64.dat\"\n",
    "input_idmap = \"temp/mit_matrix_n2_k64.dat_nameMap\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "matrix = deserialize_matrix(input_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "idmap = deserialize_id_map(input_idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1599"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 'Fclt_building_hist_1.csv->Date Occupied'),\n",
       " (1, 'All_olap_descriptor_types.csv->Descriptor Type'),\n",
       " (2, 'Sdo_units_of_measure.csv->Unit Of Meas Type'),\n",
       " (3, 'Sdo_crs_engineering.csv->Data Source'),\n",
       " (4, 'short_fclt_rooms_hist.csv->Fclt Floor Key')]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(idmap.items())[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "aggregated_jc = [np.sum(row) for row in matrix]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "448.19444333"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(aggregated_jc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(aggregated_jc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_zeros = [len(matrix) - np.count_nonzero(row) for row in matrix]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1599"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(num_zeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1133"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(num_zeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ra-mit/development/virtualenvs/fuzzylazo/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "avg_jc_per_column_norm_nonzeros = [aggr/(len(matrix) - zeros) for aggr, zeros in zip(aggregated_jc, num_zeros)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(avg_jc_per_column_norm_nonzeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06719495726666666"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(avg_jc_per_column_norm_nonzeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_jc_per_column = [aggr/len(matrix) for aggr in aggregated_jc if aggr > 0]  # not counting 0 aggr values!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2802967125265791"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(avg_jc_per_column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0006173730456535334"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(avg_jc_per_column)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proposal of computing beta and theta matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_matrix = []\n",
    "for row in matrix:\n",
    "    aggr_jc = sum(row)\n",
    "    total_d = len(row)\n",
    "    if aggr_jc > 0:\n",
    "        beta = math.log((total_d/aggr_jc))\n",
    "        theta_row = []\n",
    "        for jc in row:\n",
    "            theta = (jc * beta) / math.log(total_d)\n",
    "            theta_row.append(theta)\n",
    "    else:\n",
    "        theta_row = [0 for jc in row]\n",
    "    theta_matrix.append(theta_row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_matrix = np.asarray(theta_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_df = transform_matrix_dataframe(matrix)\n",
    "theta_matrix_df = transform_matrix_dataframe(theta_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x111e062b0>]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFd5JREFUeJzt3X+sXOV95/H319cm3DiAk2KacjFrkhhTAlUhV0AWqaFKWBySNSipWrxB2aQUq9mS/UFlCRSaUpoVzaKl2kpsU9ONkqYJlLSR4wpXVpQli5TFXi51C5jUreM02BcUnICtNHbAP777x8z1HcYz5547njMzd+b9kizPnHPmzFePru/Hz/Oc85zITCRJamdRvwuQJA02g0KSVMigkCQVMigkSYUMCklSIYNCklSosqCIiM9HxEsR8Wyb/RERfxQRuyPi6Yi4vKpaJEmdq7JH8QVgTcH+9wOr6n/WA39cYS2SpA5VFhSZ+TjwcsEhNwB/ljXbgGUR8XNV1SNJ6sziPn73BLC34f2++rYXmw+MiPXUeh0sXbr0XRdddFFPCpSkYfHUU0/9MDOXd/LZfgZFaZm5EdgIMDk5mVNTU32uSJIWloj4fqef7edVT9PAiob359W3SZIGSD+DYjPw0frVT1cBBzPzpGEnSVJ/VTb0FBEPAdcAZ0fEPuB3gSUAmfk5YAtwPbAbOAR8vKpaJEmdqywoMnPdHPsT+K2qvl+S1B3emS1JKmRQSJIKGRSSpEIGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQgaFJKmQQSFJKmRQSJIKGRSSpEIGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQgaFJKmQQSFJKmRQSJIKGRSSpEIGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQpUGRUSsiYhdEbE7Iu5osf/8iHgsInZExNMRcX2V9UiS5q+yoIiIMeAB4P3AxcC6iLi46bC7gEcy8zLgJuB/VlWPJKkzVfYorgB2Z+aezHwNeBi4oemYBM6svz4LeKHCeiRJHagyKCaAvQ3v99W3NbobuDki9gFbgE+2OlFErI+IqYiY2r9/fxW1SpLa6Pdk9jrgC5l5HnA98KWIOKmmzNyYmZOZObl8+fKeFylJo6zKoJgGVjS8P6++rdEtwCMAmfkEcDpwdoU1SZLmqcqgeBJYFREXRMRp1CarNzcd8zzwXoCI+HlqQeHYkiQNkMqCIjOPArcBW4HvULu6aWdE3BMRa+uH/TZwa0T8PfAQ8LHMzKpqkiTN3+IqT56ZW6hNUjdu+3TD6+eAq6usQZJ0avo9mS1JGnAGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQgaFJKmQQSFJKmRQSJIKGRSSpEIGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQgaFJKmQQSFJKmRQSJIKGRSSpEIGhSSpkEEhSSpkUEiSChkUkqRCBoUkqZBBIUkqZFBIkgoZFJKkQgaFJKlQpUEREWsiYldE7I6IO9oc86sR8VxE7IyIr1RZjyRp/hZXdeKIGAMeAK4F9gFPRsTmzHyu4ZhVwJ3A1Zn5SkScU1U9kqTOVNmjuALYnZl7MvM14GHghqZjbgUeyMxXADLzpQrrkSR1oMqgmAD2NrzfV9/W6ELgwoj4dkRsi4g1rU4UEesjYioipvbv319RuZKkVvo9mb0YWAVcA6wDHoyIZc0HZebGzJzMzMnly5f3uERJGm1VBsU0sKLh/Xn1bY32AZsz80hmfg/4R2rBIUkaEFUGxZPAqoi4ICJOA24CNjcds4lab4KIOJvaUNSeCmuSJM1TZUGRmUeB24CtwHeARzJzZ0TcExFr64dtBX4UEc8BjwEbMvNHVdUkSZq/yMx+1zAvk5OTOTU11e8yJGlBiYinMnOyk8/2ezJbkjTgDApJUiGDQpJUyKCQJBWac62niDgduAV4J3D6zPbM/PUK65IkDYgyPYovAW8FrgP+D7Ub535cZVGSpMFRJijekZm/A/wkM78IfAC4stqyJEmDokxQHKn/fSAiLgHOAlwOXJJGRJnnUWyMiDcDd1FbguNNwO9UWpUkaWCUCYpv1p8X8TjwNoCIuKDSqiRJA6PM0NNftdj2l90uRJI0mNr2KCLiImqXxJ4VER9q2HUmDZfJSpKGW9HQ02rgg8Ay4N82bP8xtUeYSpJGQNugyMyvA1+PiHdn5hM9rEmSNEDKzFH8KCK+GRHPAkTEL0TEXRXXJUkaEGWC4kHgTur3U2Tm09SeVidJGgFlguKNmfn/mrYdraIYSdLgKRMUP4yItwMJEBG/ArxYaVWSpIFR5oa73wI2AhdFxDTwPeDmSquSJA2MOYMiM/cA74uIpcCizHTlWElaIDbtmOa+rbs47a3veFen5yjzPIo3AB8GVgKLIwKAzLyn0y+VJFXrrk3P8OXtz5N56ucqM/T0deAg8BTw6ql/pSSpSndteoY/3/Z8185XJijOy8w1XftGSVJXzAwrvXDgMOcuG2fDdau58bIJHtq+t6vfUyYo/m9EXJqZz3T1myVJHdu0Y5o7v/YMh48cA2D6wGHu/Frt1/Sxbow3NShaFPBZ4Hj9mI9HxB5qQ08BZGb+QlcrkSS11dx7OPTa0RMhMePwkWPct3UXYxFdDYuiHsUE8Itd+yZJUkda9R7aeeHAYT5y1fk9m6P4XmZ+v2vfJElqq918A8B9W3ed1Hto59xl43zmxksBeGj7Xo5l1noYh3+8v9PaioLinIi4vd3OzLy/0y+VJM0qmm+48bIJXijoQTQaXzLGhutWA/CZGy89ERgA8Qc/6LiLUbSExxi152Of0eaPJKkLWvUYZuYboNZLaGXZ+BImlo0TwMSyce790KUneiHdVNSjeNGb6iSpeu16DDPbN1y3+nU9Dqj1Hu5e+85KgqFZUY8iKv92SVLbHsPM9hsvm+DeD13ak95DK0U9ivf2pAJJGnHtegwz8w1QC4teBUOzokehvtzLQiRpVDVe3dTqqqd+K3NntiSpYv3sMcylzIOLJEkjrNKgiIg1EbErInZHxB0Fx304IjIiJqusR5I0f5UFRUSMAQ8A7wcuBtZFxMUtjjsD+E/A9qpqkSR1rsoexRXA7szck5mvAQ8DN7Q47veBzwI/rbAWSVKHqgyKCaBxUfR99W0nRMTlwIrMfLToRBGxPiKmImJq//6OlyuRJHWgb1c9RcQi4H7gY3Mdm5kbgY0Ak5OT3V1oXdLQKlpoT+VVGRTTwIqG9+fVt804A7gE+Fb9OdxvBTZHxNrMnKqwLkkjYK6F9lRelUNPTwKrIuKCiDgNuAnYPLMzMw9m5tmZuTIzVwLbAENCUlfMtdCeyqssKDLzKHAbsBX4DvBIZu6MiHsiYm1V3ytJMPdCeyqv0jmKzNwCbGna9uk2x15TZS2SFrb5zjecu2y85ZPg2i3Ap/a8M1vSwJuZb5g+cJhkdr5h047ptp/ZcN1qxpeMvW5b80J7KsegkDTwOplv6PfS3MPERQElDbxO5xsGeaG9hcQehaSBN9eDfVQtg0LSwHO+ob8cepI08Ab9wT7DzqCQtCA439A/Dj1JkgrZo5A0bzM3v00fOMxYBMcymXA4aGgZFJJK27Rjmrs37+TA4SMnth3L2oLOLro3vBx6klTKzN3RjSHRzEX3hpNBIamUVndHt+Kie8PHoJBUStkA8Ca44eMchTSC7tr0DA9t38uxTMYiWHflCj5z46WFn2m3Gmsjb4IbTgaFNAIal+geX7KIQ0eOn9h3LJM/3/Y8QGFYbLhu9eueGDdjUcDxxKuehphBIQ255keCNoZEo4e27y0MCu+OHl0GhTQk2j3Yp+wk9MxlrkW8O3o0GRTSArdpxzS/99c7eeXQ7GWrjfc0lJ2EHouopD4tfAaFtAA13hkdQKu+wMw9DWUmoQHWXbmi63VqOBgU0gLQOKy07I1L+JefHuXI8Vo8FA0YvXDgMH/4a7940iT0oqh9LpPSVz1pdBkU0gBqDIazxpfwk9eOcuRYLRIah5jmcu6ycSehdcoMCmlAtBtOKloyo0jjPQ1OQutUGBTSAGi+hHXu64+KLRtfwt1r32k4qCsMCmkAlL2EtZ2ZHog3vakKBoU0AOa7kN6SsWDpaYs5ePiIcw6qnEEhVegjDz7Bt7/78on3V7/9LXz51nefdNxcl7AuWRS86fTFHDhkMKj3DAqpIs0hAfDt777MRx584qSwaLWOksNJGhQGhdShxhVYG830GppDYkar7V7CqkFmUEjzdO393+KfXvpJ2/0zvYb58hJWDSqDQirprk3PnFiOey7tehPSQmRQSC00r8S68mfG5/3L/+q3v6XlZ65++1u6VabUEwaFROuJ5xnTBw6XWlSv2ZdvfXfpq56kQWZQaKTNZzhpPmZ6DYaChoFBoZHTuKZSFew1aNhUGhQRsQb4H8AY8KeZ+QdN+28HfgM4CuwHfj0zv19lTRo9RUt0d8Oqc5byjduv6dr5pEFTWVBExBjwAHAtsA94MiI2Z+ZzDYftACYz81BEfAL4b8CvVVWTRkurJ7/NZ4nuRuNLxrj8/LPYtucVjmX6DAeNlCp7FFcAuzNzD0BEPAzcAJwIisx8rOH4bcDNFdajIdZqMrrdk9/KWrIIjh7Hm9808qoMiglgb8P7fcCVBcffAvxNqx0RsR5YD3D++ed3qz4NiXZXLJ1KSDjPIM0aiMnsiLgZmATe02p/Zm4ENgJMTk52b3BZC07jshkzwz/durnN4SSptSqDYhpofFr7efVtrxMR7wM+BbwnM1+tsB4tcM2Xsh7L7OjSVpfoluanyqB4ElgVERdQC4ibgH/XeEBEXAb8CbAmM1+qsBYtQEU3wXXKJ79J81dZUGTm0Yi4DdhK7fLYz2fmzoi4B5jKzM3AfcCbgK9GBMDzmbm2qpo02BovY128CI4cP/VzulS3dOoqnaPIzC3AlqZtn254/b4qv1+Dq7m3sOqcpex75acnnscwn5AYi+Cqt735pN6H4SB1x0BMZmt0bNoxzYav/t1JQVC0bPdcnICWqmVQqCda3fzWqbEIb3qTesigUNc1L9H9yxct56+emn7dYz475f0NUu8ZFOqqTTumX/fs5+kDh/nytudP6ea3GYaE1B8Ghbrqvq27Tuo5lA2JVecs5dBrx31mtDRgDAqV0jyc1O6X+AsdLt1tb0EaXAaF5tRqOOnOrz0DcFJYnLtsvOVzHlot0OfNb9LCYFCopcYexKL6VUaNDh85xn1bd530S37DdatfFypQW6L7w++a4LF/2O+wkrQAGRSa8yql5pCY0WqYaeaXf5lhKkkLg0Exwq69/1sn3eg2n6uUzl023nL7jZdNGAzSEDEoRsTKOx4tfWyZkBhfMsaG61Z3XpCkBcOgGAHzCYkiYxEcz3Q4SRoxBsWQaXywz6lovkppfMkY937oUsNBGkEGxRC48r9+gx/8+LWunc+rlCQ1MigWoJmrlFrdr3CqFoE9B0mvY1AsIO+481GOVvjE8FXnLOUbt19T3RdIWpAMigHV7eGkIme+YYynf29NT75L0sJjUAyYKnsNrqckqRMGRZ9dcMejXVmCey6GhKROGRR9VFVI/OwZp7H9U9dWcGZJo8ig6IGqJ6FnLA7Yfe8Hqv8iSSPFoKhYlSFhMEjqBYOiYt0KCYeTJPWLQTHg7DVI6jeDYgBNuGSGpAFiUFRscRQPP9ljkDToFvW7gGG3+94PsDha7zMkJC0E9ih6wDCQtJDZo5AkFTIoJEmFDApJUiGDQpJUyKCQJBUyKCRJhQwKSVKhSoMiItZExK6I2B0Rd7TY/4aI+Iv6/u0RsXKucz4zfZCVdzzKyjseraJkSVKTyoIiIsaAB4D3AxcD6yLi4qbDbgFeycx3AH8IfHY+32FYSFL1quxRXAHszsw9mfka8DBwQ9MxNwBfrL/+S+C9EdFmwQtJUj9UuYTHBLC34f0+4Mp2x2Tm0Yg4CPwM8MPGgyJiPbAeYNH4mbz4xf88u++zH3yq65UvHGfT1FYjzLaYZVvMsi1mre70gwtirafM3AhsBIiIqVcPHZzsc0kDISKmMtO2wLZoZFvMsi1mRcRUp5+tcuhpGljR8P68+raWx0TEYuAs4EcV1iRJmqcqg+JJYFVEXBARpwE3AZubjtkM/Pv6618B/ndmVvSEaUlSJyobeqrPOdwGbAXGgM9n5s6IuAeYyszNwP8CvhQRu4GXqYXJXDZWVfMCZFvMsi1m2RazbItZHbdF+B94SVIR78yWJBUyKCRJhQY2KKpY/mOhKtEWt0fEcxHxdER8MyL+VT/q7IW52qLhuA9HREbE0F4aWaYtIuJX6z8bOyPiK72usVdK/Bs5PyIei4gd9X8n1/ejzqpFxOcj4qWIeLbN/oiIP6q309MRcXmpE2fmwP2hNvn9XeBtwGnA3wMXNx3zH4DP1V/fBPxFv+vuY1v8MvDG+utPjHJb1I87A3gc2AZM9rvuPv5crAJ2AG+uvz+n33X3sS02Ap+ov74Y+Od+111RW/wScDnwbJv91wN/AwRwFbC9zHkHtUfh8h+z5myLzHwsMw/V326jds/KMCrzcwHw+9TWDftpL4vrsTJtcSvwQGa+ApCZL/W4xl4p0xYJnFl/fRbwQg/r65nMfJzaFaTt3AD8WdZsA5ZFxM/Ndd5BDYpWy39MtDsmM48CM8t/DJsybdHoFmr/YxhGc7ZFvSu9IjOHfcXIMj8XFwIXRsS3I2JbRKzpWXW9VaYt7gZujoh9wBbgk70pbeDM9/cJsECW8FA5EXEzMAm8p9+19ENELALuBz7W51IGxWJqw0/XUOtlPh4Rl2bmgb5W1R/rgC9k5n+PiHdTu3/rksw83u/CFoJB7VG4/MesMm1BRLwP+BSwNjNf7VFtvTZXW5wBXAJ8KyL+mdoY7OYhndAu83OxD9icmUcy83vAP1ILjmFTpi1uAR4ByMwngNOpLRg4akr9Pmk2qEHh8h+z5myLiLgM+BNqITGs49AwR1tk5sHMPDszV2bmSmrzNWszs+PF0AZYmX8jm6j1JoiIs6kNRe3pZZE9UqYtngfeCxARP08tKPb3tMrBsBn4aP3qp6uAg5n54lwfGsihp6xu+Y8Fp2Rb3Ae8CfhqfT7/+cxc27eiK1KyLUZCybbYCvybiHgOOAZsyMyh63WXbIvfBh6MiP9CbWL7Y8P4H8uIeIjafw7Ors/H/C6wBCAzP0dtfuZ6YDdwCPh4qfMOYVtJkrpoUIeeJEkDwqCQJBUyKCRJhQwKSVIhg0KSVMigkE5BRPxL/e8LI2JLRPxTRPxtRDwSET/b7/qkbhjI+yikhSQiTgceBW7PzL+ub7sGWA78oI+lSV3hfRTSKaj3KP4jcE1mfrTf9UhVcOhJOnWXAE/1uwipKgaFJKmQQSGdup3Au/pdhFQVg0I6dV8B/nVEfGBmQ0T8UkRc0seapK4xKKQO1Z+D8mpmHgY+CHyyfnnsc9Se6T6Ky1hrCHl5rNS5dwLfBcjMfwCG9VGjGnH2KKQORMRvAg8Bd/W7Fqlq3kchSSpkj0KSVMigkCQVMigkSYUMCklSIYNCklTo/wORiBIoS9yJfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.ylim((0, 1))\n",
    "plt.xlim((0, 1))\n",
    "ax.set_xlabel(\"JC\")\n",
    "ax.set_ylabel(\"Theta\")\n",
    "ax.plot(matrix[1], theta_matrix[1], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11389fd68>]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHkxJREFUeJzt3X2UHXWd5/H3J51HAjEEQzaExESIsFmF4PSYYJw5CGaJqBAdFsyJbnaGNTOrOwvCRhLJGcCRFZfdgDPrcQyDY3ZEnrFBYMFMJI5yMGNiB1oeMgHkqQ1JlAQwRkjCd/+o6nC70w91u2/dp/q8zunT99atW/ebOnA/Xb/6PSgiMDOz4hpW6wLMzKy2HARmZgXnIDAzKzgHgZlZwTkIzMwKzkFgZlZwDgIzs4JzEJiZFZyDwMys4IbneXBJnwf+MxBAB/CnwGTgZuAoYBPw6Yh4o7/jvP3tb4/p06fnWaqZWdPZtGnTryNi4kD7Ka8pJiRNAX4CzIqIvZJuBe4DzgLujIibJf0d8EhEfKO/Y7W2tsbGjRtzqdPMrFlJ2hQRrQPtl3fT0HBgjKThwGHANuB04Pb09TXAwpxrMDOzfuQWBBHRCfwv4HmSAHiFpClod0TsT3d7EZjS2/slLZW0UdLGnTt35lWmmVnh5RYEko4EzgFmAMcAY4EFWd8fEasjojUiWidOHLCJy8zMBinPpqEPAb+MiJ0RsQ+4E5gHjE+bigCOBTpzrMHMzAaQZxA8D8yVdJgkAWcAjwMPAuem+ywB7sqxBjMzG0Bu3UcjYoOk24GfA/uBdmA1cC9ws6Qvp9tuyKsGM7NGdNLl9/Pq6wcOPh83qoVHr8zcsl623LqPVpK7j5pZUfQMgS6DCYN66T5qZmZl6C0E+tteCbmOLDYzK6qVbR3ctOEFDkTQIrFozlS+vPA9tS6rVw4CM7NBamvv5Iq7H2P33n0AHHnYCC7/2L9j43Mv852fPn9wvwMRB5/XYxi4acjMbBDa2jtZdtsjB0MAYNfv9rHs9ke4sSQESt204YUBjztuVEtZ2yvBQWBmNgjXPLCFfW8e2tlm34Ggry44BzJ0znn0ygWHfOnn3WvITUNmVjht7Z1c88AWfrV7L8eMH8OyM09g4Sm9znbTp1/t3lv257ZImfbL80u/Nw4CMyuEri//zt17ERz8q71z915W3NkBUFYYHDN+DJ19hMHYkS3seePQXj6L5kwtt+yqcNOQmTW9tvZOVtzZcfCLu2cDzd59B7jmgS1lHXPZmScwYtihf+GPaBFXffw9fGrutINXAC0Sn5o7rS5vFIOvCMyswfXVc6f0r/trHtjC3n3998Mvt6mn6/h9ffbCU6bU7Rd/Tw4CM2sYpc07fenquQNvfVln+ZI/ZvyYsuvp+sJvdA4CM6trK9s6uvXJz2LfgeCaB7Yc/JLurz0fYMyIFpadecKQ6mxkDgIzqyuLr3+Yh55+ecjHKb0KWHbmCay4s6Nb81DXDeMpg+w11EwcBGZWUz3b+CultKmn60t+qF1Gm5WDwMyqqq29k4tv2cybOX7GiBYd0tTTLO35eXAQmFnustzkrZTeeg1Z/xwEZlZx1fzir/eZPRtBbkEg6QTglpJN7wT+Cvi/6fbpwLPAeRGxK686zCx/be2dXPn9x9j1u8q28/dGwOI6HpzViPJcqnILMBtAUgvJIvXfA5YD6yLiaknL0+eX5lWHmeUnrxu9peYdN4EbP3Nqbse36jUNnQE8HRHPSToHOC3dvgZYj4PArCGUTtY2/rAR/Pb3+3udgXOo3M5fXdUKgk8CN6WPJ0XEtvTxS8CkKtVgZoM0f9V6tu7Y021bpZuB3J+/dnIPAkkjgbOBFT1fi4iQ1OufE5KWAksBpk2blmuNZtbdYEbzlstf/PWjGlcEHwZ+HhHb0+fbJU2OiG2SJgM7entTRKwGVgO0trZW/trTzLqpxpe/2/vrUzWCYBFvNQsB3A0sAa5Of99VhRrMrBcr2zq48afP97mi1lCMHzOCK852O38jyDUIJI0F5gN/XrL5auBWSRcAzwHn5VmDmXXX1t7JRbdsrsixRrSIsSOH88refZ62oYHlGgQRsQc4qse235D0IjKzKqnkl3+XFsE1557sL/4m4JHFZk2sUjN59jTz6LGsvfi0ih/XasNBYNZk2to7ueTWzRzIoeG/npdbtMFzEJg1iTyaf7qMG9XCo1cuyOXYVnsOArMGllfTj//yLxYHgVkDmnPVWra/9kbFj3vd+bN987eAHARmDaKtvZNlt21mXwVXdJl0xEg2XDa/cge0huQgMKtjeXz5u9nHenIQmNWhSrf9DwNWudnH+uAgMKsTeazl6/7+loWDwKzG8ljcxc0/Vg4HgVkNVXLStxbB/z7PzT9WPgeBWZWtbOvgxg3PExX49h8xDK75D/7yt6FxEJhVyYmX3cfvKzTvg9v+rZIcBGZVMH35vRU71rNXf6RixzIDB4FZLiq92pdX9rI8OQjMKqhS/f9bJBbNmeqeP1YVDgKzCqjEFcCRh43g8o95aUervryXqhwP/D3wbiCAPwO2ALcA04FngfMiYleedZjlaf6q9WzdsWfQ7z9sxDD+xydOcgBYzeR9RfA14P6IOFfSSOAw4IvAuoi4WtJyYDlwac51mFVUW3sn1zywhc7de4d0HN/4tXqQWxBIehvwx8B/AoiIN4A3JJ0DnJbutgZYj4PAGshQrwC6OASsXuR5RTAD2An8g6STgU3AhcCkiNiW7vMSMKm3N0taCiwFmDZtWo5lmmUz1BvB7vtv9SrPIBgOvBf4y4jYIOlrJM1AB0VESOp1hE1ErAZWA7S2tuaw+qpZNkOZClrAtZ710+pcnkHwIvBiRGxIn99OEgTbJU2OiG2SJgM7cqzBbEja2jtZcWfHoELAff+tUeQWBBHxkqQXJJ0QEVuAM4DH058lwNXp77vyqsGsXJVo/x87soWrPv4eXwVYw8i719BfAjemPYaeAf6UZI2MWyVdADwHnJdzDWaZDCUE3P5vjSzXIIiIzUBrLy+dkefnmpVr8fUPDzoE3ARkjc4ji63wBtsbyAFgzcJBYIU0mC9/AceMH8OyM09w+781FQeBFUpbeyeX3vEor+8vrxuQl360ZuYgsEIYyrrADgFrdg4Ca2pt7Z1c+f3H2PW78gPA9wCsKBwE1rSGsjC8Q8CKxEFgTWdlWwff3fA8b5aZAP7yt6JyEFhTGewCMQ4BKzIHgTWVmza8kHlfrwhmlnAQWFM5EAO3B40fM4IrznYAmHVxEFhTaZH6DAMBi90V1OwQw2pdgFklLZoztdftY0YM49rzZzsEzHrhKwJrKl1f9DdteIEDEbRILJoz1QFg1g9FhjbVWmttbY2NGzfWugwzs4YiaVNE9DYDdDe+IrC6trKtw3/dm+XMQWB1q+eYgAMRB587DMwqJ9ebxZKeldQhabOkjem2CZLWStqa/j4yzxqscfU1JqCcsQJmNrBq9Br6YETMLmmnWg6si4iZwLr0udkh+uoGmmWsgJllV4vuo+cAa9LHa4CFNajBGkCLVNZ2MxucvO8RBPADSQF8MyJWA5MiYlv6+kvApJxrsDrVc42AnlM+LJoztdd5g/oaK2Bmg5N3EHwgIjolHQ2slfRk6YsREWlIHELSUmApwLRp03Iu06ptzlVr2f7aG9227frdPpbd/ggAC0+Z4jEBZlVStXEEkq4Afgt8BjgtIrZJmgysj4gT+nuvxxE0l/mr1rN1x54+X58yfgwPLT+9ihWZNaeKjSOQNBG4FJgFjO7aHhH9/p8qaSwwLCJeSx//e+BLwN3AEuDq9PddA9VgzaGtvZMv3P4Ibxzo/4+PX+3eW6WKzAyyNQ3dCNwCfAT4C5Iv750Z3jcJ+J6SG3vDge9GxP2SfgbcKukC4DngvMEUbo2lt6agvhwzfkzO1ZhZqSxBcFRE3CDpwoj4EfCj9Mu8XxHxDHByL9t/A5xRfqnWiI5fcS/7y2x9XHZmvy2FZlZhWYKga9XvbZI+AvwKmJBfSdYsBhMC846b4HUCzKosSxB8WdLbgEuAvwXGARflWpU1vMXXP1x2CHzKawWY1USWINgVEa8ArwAfBJA0L9eqrGG1tXey7LbN7Hsz+3u8XrBZbWUJgr8F3pthmxXc9OX3lv2ecaNaHAJmNdZnEEg6FXg/MFHSxSUvjQNa8i7MGstgQsBXAmb1ob8rgpHA4ek+R5RsfxU4N8+irHG0tXdy0S2by3rPcMFTX/lIThWZWbn6DIKSrqLfjojnJB0WEb+rYm1W58oZG9DFIWBWf7LcIzhG0v8juTqYJulk4M8j4rP5lmb1bP6q9WWFgJuBzOpXliC4DjiTZGoIIuIRSX+ca1VWtxZf/zAPPf1yWe+57vzZHhtgVscyrUcQET2XhDqQQy1W5wYTAp+aO80hYFbnslwRvCDp/UBIGgFcCDyRb1lWb1a2dZQVAgKu9ZWAWUPIEgR/AXwNmAJ0Aj8APpdnUVY/BtMraNyoFh69ckFOFZlZpQ0YBBHxa2BxFWqxOrOyraPXFcL642kizBpP1vUIPgNML90/Iv4sv7Ks1k66/H5efT37raBnr3aXULNGlaVp6C7gx8A/4ZvEhTBj+b2UM1/cvOM8Ga1ZI8sSBIdFxKW5V2J14fgV5YXA6BZ5fIBZg8vSffQeSWflXonVVFt7J9OXl7d+wKQjRvLkVf5Pw6zR9Tfp3GtAkPQE/KKk10kWqREQETEuywdIagE2Ap0R8VFJM4CbgaOATcCnI6K8eQqsogbTM8iDxMyaR39NQ++OiOcq8Bld4w66guOrwLURcbOkvwMuAL5Rgc+xQSh3vqCZR49l7cWn5VeQmVVdf01D3xvqwSUdS7Lo/d+nzwWcDtye7rIGWDjUz7HBKTcErjt/tkPArAn1FwSqwPGvA74AdK1XdRSwOyL2p89fJBmoZlU2mBBwU5BZc+qvaWiKpL/p68WI+G/9HVjSR4EdEbFJ0mnlFiZpKbAUYNq0aeW+3fpRziIyAn7pMQJmTa2/INhLcjN3sOYBZ6c9jkaT3CP4GjBe0vD0quBYkmkrDhERq4HVAK2trWUug269KfemsNcOMCuG/oLgNxGxZrAHjogVwAqA9Irgv0fEYkm3kaxwdjOwhGTAmuWs3BAY3SJ3DTUriP7uEeTVpfNS4GJJT5HcM7ghp8+xVLkh4PEBZsXS31KVcyv1IRGxHlifPn4GeF+ljm39K3fiOE8aZ1Y8WaaYsAZVzk1hcAiYFZWDoEmVGwLuHmpWXA6CJlRuCHgKabNiy7RmcU+S7ql0IVYZDgEzK9eggoBkoRqrMyddfn9Z+zsEzAwGCAJJLZJu7Lk9IrblV5INxvEr7s28ophwCJjZW/oNgog4ALxD0sgq1WODcOJl92VeR2DSESM9ZYSZdZPlZvEzwEOS7gb2dG2MiFW5VWWZlbOs5LzjJng1MTM7RJYgeDr9GQYckW85Vo5ybgw7BMysLwMGQURcWY1CrDwOATOrlAGDQNL34ZDWh1dIlp/8ZkT8Po/CrG/lhIAHipnZQLJ0H30G+C1wffrzKvAa8K70uVWRQ8DMKi3LPYL3R8Qfljz/vqSfRcQfSnosr8LsUA4BM8tDliA4XNK0iHgeQNI04PD0tbymqrYeygkBjxEws3JkCYJLgJ9IeppkLNIM4LOSxpIsPm85cwiYWZ6y9Bq6T9JM4MR005aSG8TX5VaZAQ4BM8tfll5Dn+ix6ThJrwAdEbEjn7IMksFiWc07bkKOlZhZM8vSNHQBcCrwQ5KmodNIFrWfIelLEfGPvb1J0mjgn4FR6efcHhGXS5pBsl7xUelxPh0RvtfQw/Erso8YBjxOwMwGLUv30eHAv42IcyPiT4BZJOMK5pCsP9yX14HTI+JkYDawQNJc4KvAtRFxPLCLJGisxEmX35957iBwk5CZDU2WIJgaEdtLnu9It70M7OvrTZH4bfp0RPoTwOnA7en2NcDCsqtuYuXMIjrvuAkOATMbsixNQ+vThWhuS5+fC/wo7TW0u783Smohaf45Hvg6yZxFuyNif7rLi0Cvnd0lLQWWAkybNi1DmY1v/qr1ma8EZh491s1BZlYRWa4IPgf8A0nzzmxgTUT8l4jYExEf7O+NEXEgImYDxwLv462eRwOKiNUR0RoRrRMnTsz6toa1sq2DrTv2DLwjyZXA2otPy7cgMyuMAYMgbeK5IyI+HxGfB7ZL+no5HxIRu4EHSW46j5fUdSVyLNBZbtHNZmVbB9/56fOZ9h3dIl8JmFlFZVqqUtIpkv6npGeBLwFPZnjPREnj08djgPnAEySBcG662xLgrkHU3VSyhoCAJ686K99izKxw+rxHIOldwKL059fALYAGag4qMRlYk94nGAbcGhH3SHocuFnSl4F24Iah/AMaXdYBY+NGtfDolQtyrsbMiqi/m8VPAj8GPhoRTwFI+nzWA0fEo8ApvWx/huR+QeGVM2DMIWBmeemvaegTwDbgQUnXSzqDpHXCKmDx9Q9nHjDmLqJmlqc+rwgiog1oS7uJngNcBBwt6RvA9yLiB1WqsenMuWot21/LNpjaIWBmecvSa2hPRHw3Ij5G0sunnf5HFFs/HAJmVm8y9RrqEhG70v79Z+RVULNzCJhZvSkrCGxojl+R7ebw6BbfijGz6nEQVMmJl92XafqI4fJYATOrLgdBFcxftZ7fHxg4BWYePZanvuImITOrriyTztkQlHNz2PMHmVkt+IogR4uvf9g3h82s7jkIcvTQ0y9n2s8hYGa15CDISdY5hNxDyMxqzUGQgzlXrc28r3sImVmtOQgqrK290/cFzKyhOAgq7KJbNmfazyFgZvXCQVBBWe8LOATMrJ44CCrkxMvuy7Tfp+ZOy7kSM7Py5BYEkqZKelDS45Iek3Rhun2CpLWStqa/j8yrhmppa+/MNHIY4MsL35NzNWZm5cnzimA/cElEzALmAp+TNAtYDqyLiJnAuvR5w2pr7/R9ATNraLkFQURsi4ifp49fI1m4fgrJIjdr0t3WAAvzqqEaHAJm1uiqco9A0nSS9Ys3AJMiYlv60kvApGrUkIeVbR2Z9nMImFk9yz0IJB0O3AFcFBGvlr4WEQG9L90raamkjZI27ty5M+8yB+U7P31+wH1mHj22CpWYmQ1erkEgaQRJCNwYEXemm7dLmpy+PhnY0dt705XQWiOideLEiXmWOSgzMnYV9YyiZlbv8uw1JOAG4ImIWFXy0t3AkvTxEuCuvGrIy4zl9/Z+GdODm4TMrBHkuR7BPODTQIekrjuqXwSuBm6VdAHwHHBejjVU3EmX358pBK47f3butZiZVUJuQRARPwH6mlrzjLw+N2+vvn5gwH3GjWph4SlTqlCNmdnQeWRxGbJOIfHolQtyrsTMrHIcBBllnULC9wXMrNE4CDJY2daRaQqJcaNaqlCNmVllOQgyyDJeQLhJyMwak4NgAFlHD//STUJm1qAcBAPIejVgZtaoHAT9OH5Ftl5Cvhows0bmIOjDyrYO9mcYOeZeQmbW6BwEfcjSJOReQmbWDBwEvZhz1dpM+7mXkJk1AwdBD23tnWx/7Y0B93OTkJk1CwdBD1lXHDMzaxYOghJt7Z2Z9vPVgJk1EwdBiUtuHfhqwCuOmVmzcRCkVrZ1MNB0QsIrjplZ83EQpLJ0F/XAMTNrRg4CYPH1Dw+4j+8LmFmzynPN4m9J2iHpFyXbJkhaK2lr+vvIvD6/HA89/XKtSzAzq5k8rwi+DfQccbUcWBcRM4F16fOayrLq2LzjJlShEjOz2sgtCCLin4Gef2qfA6xJH68BFub1+VnMyBACo1vEjZ85tQrVmJnVRrXvEUyKiG3p45eASX3tKGmppI2SNu7cuTOXYgaaU250i3jyqrNy+Wwzs3pRs5vFERH0810cEasjojUiWidOnFjxz8/SJOQQMLMiqHYQbJc0GSD9vaPKn29mZj1UOwjuBpakj5cAd1X584Fss4v6BrGZFUWe3UdvAh4GTpD0oqQLgKuB+ZK2Ah9Kn1ddltlFfYPYzIpieF4HjohFfbx0Rl6fmUWWewPXnT+7CpWYmdUHjyzuQcDCU6bUugwzs6opVBBkmWba8wmZWdEUKgi86IyZ2aEKEwQnXX7/gPv43oCZFVFhguDV1w8MuI/vDZhZERUiCLJMM+1xA2ZWVIUIgizTTHvcgJkVVSGCYCCTjhhZ6xLMzGqm6YMgywCyDZfNr0IlZmb1qemDwMzM+tfUQTB/1foB9/FaxGZWdE0dBFt37Kl1CWZmda+pg2AgqnUBZmZ1oGmDIMtNYs8rZGbWxEFgZmbZFDYIfJPYzCxRkyCQtEDSFklPSVpeixrMzCxR9SCQ1AJ8HfgwMAtYJGlWteswM7NELa4I3gc8FRHPRMQbwM3AOZX8gIFuFLtZyMzsLbUIginACyXPX0y3mZlZDdTtzWJJSyVtlLRx586dFTuurwbMzLqrRRB0AlNLnh+bbusmIlZHRGtEtE6cOLFqxZmZFU0tguBnwExJMySNBD4J3F2DOszMjBoEQUTsB/4r8ADwBHBrRDxWyc/oq/nHzUJmZocaXosPjYj7gPvy/Ax/6ZuZZVO3N4vNzKw6HARmZgXnIDAzKzgHgZlZwTkIzMwKThFR6xoGJGkn8Nwg3/524NcVLKfR+Xx05/PRnc9Hd41+Pt4REQOOyG2IIBgKSRsjorXWddQLn4/ufD668/norijnw01DZmYF5yAwMyu4IgTB6loXUGd8Prrz+ejO56O7QpyPpr9HYGZm/SvCFYGZmfWjqYNA0gJJWyQ9JWl5reupBknfkrRD0i9Ktk2QtFbS1vT3kel2Sfqb9Pw8Kum9tau88iRNlfSgpMclPSbpwnR7Uc/HaEn/IumR9HxcmW6fIWlD+u++JZ0eHkmj0udPpa9Pr2X9eZHUIqld0j3p88Kdj6YNAkktwNeBDwOzgEWSZtW2qqr4NrCgx7blwLqImAmsS59Dcm5mpj9LgW9UqcZq2Q9cEhGzgLnA59L/Bop6Pl4HTo+Ik4HZwAJJc4GvAtdGxPHALuCCdP8LgF3p9mvT/ZrRhSRT4ncp3vmIiKb8AU4FHih5vgJYUeu6qvRvnw78ouT5FmBy+ngysCV9/E1gUW/7NeMPcBcw3+cjAA4Dfg7MIRkwNTzdfvD/G5I1Q05NHw9P91Ota6/weTiW5I+B04F7ABXxfDTtFQEwBXih5PmL6bYimhQR29LHLwGT0seFOUfpZfwpwAYKfD7SZpDNwA5gLfA0sDuSBaOg+7/54PlIX38FOKq6FefuOuALwJvp86Mo4Plo5iCwXkTy50yhuopJOhy4A7goIl4tfa1o5yMiDkTEbJK/hN8HnFjjkmpG0keBHRGxqda11FozB0EnMLXk+bHptiLaLmkyQPp7R7q96c+RpBEkIXBjRNyZbi7s+egSEbuBB0maPsZL6lqtsPTffPB8pK+/DfhNlUvN0zzgbEnPAjeTNA99jQKej2YOgp8BM9MeACOBTwJ317imWrkbWJI+XkLSVt61/T+mvWXmAq+UNJk0PEkCbgCeiIhVJS8V9XxMlDQ+fTyG5H7JEySBcG66W8/z0XWezgV+mF5BNYWIWBERx0bEdJLvhx9GxGKKeD5qfZMizx/gLOBfSdpBL6t1PVX6N98EbAP2kbRvXkDSjrkO2Ar8EzAh3VckPaueBjqA1lrXX+Fz8QGSZp9Hgc3pz1kFPh8nAe3p+fgF8Ffp9ncC/wI8BdwGjEq3j06fP5W+/s5a/xtyPDenAfcU9Xx4ZLGZWcE1c9OQmZll4CAwMys4B4GZWcE5CMzMCs5BYGZWcA4CKyRJCyWFpKqNrJX0rKS3p4//jaSbJT0taZOk+yS9q1q1mJVyEFhRLQJ+kv6uiJLRqAPtJ+B7wPqIOC4i/oBkUsRJ/b/TLB8eR2CFk849tAX4IPD9iDgh3T4M+D8kUw28QDIo71sRcbuks4BVwB7gIZLBRB+VdAVwHMkgpOcjos9gSacyaCUZ2HVFRPxxPv9Cs/L4isCK6Bzg/oj4V+A3kv4g3f4Jkim8ZwGfJpmHB0mjSaao/nD61/vEHsebBXyovxDo4d1A4Sc6s/rhILAiWkQyyRjp764v8A8At0XEmxHxEsmcM5DM0PlMRPwyfX5Tj+PdHRF78yzYLE+Z2jTNmoWkCSRNP++RFEALEJKWDeGwe8rc/zHemtTMrOZ8RWBFcy7wjxHxjoiYHhFTgV8Cf0TS9v8nkoZJmkQyERkk9xPeWbJG7flDrOGHwChJS7s2SDpJ0h8N8bhmg+IgsKJZRNJjp9Qd6fY7SGZsfRz4DslSjq+kzT6fBe6XtAl4jWR1qkNIOkbSfX189nDg9Uh6aHwc+FDaffQx4Cskq6WZVZ17DZmVkHR4RPxW0lEkUw3Pi4iXSrZ3TVW9NSKuLeO4E4HNEdFUS19ac/AVgVl396Rr+v4Y+Ov0pjHAZ9Ltj5GsTPXNrAeUdHZ6vBWVLtasEnxFYGZWcL4iMDMrOAeBmVnBOQjMzArOQWBmVnAOAjOzgnMQmJkV3P8H37gbdgRbS5QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "#plt.ylim((0, 1))\n",
    "#plt.xlim((0, 1))\n",
    "ax.set_xlabel(\"Aggr. JC\")\n",
    "ax.set_ylabel(\"Aggr. Theta\")\n",
    "ax.plot([sum(x) for x in matrix], [sum(x) for x in theta_matrix], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11160df60>]"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHnNJREFUeJzt3X2UVPWd5/H3x+ZBNCIB0T0i2ESIDkZHjh3RmGRNDCNJJsIkzAKDu54dj55s4tnJOOMGAhvRjWvc7NFkJ85OIOaMMSaSEO1tI9FlVs1JPMjQpIk9kDA2SAKdJ4KIT4CC3/2jbmNRdHfd6q5bj5/XOX2ouvW7Vb+fhf3h3t+TIgIzM7PBnFDtCpiZWe1zWJiZWVEOCzMzK8phYWZmRTkszMysKIeFmZkV5bAwM7OiHBZmZlaUw8LMzIoakeWbS5oDfAVoAb4eEV8seH008E3gYmAvsCAidkpqBX4ObEuKPhMRnxzss0477bRobW0ta/3NzBrdpk2b/hARE4uVyywsJLUA9wCzgd3ARkkdEbE1r9h1wL6ImCZpIXAnsCB5bXtEXJT281pbW+ns7CxT7c3MmoOkX6Ypl+VtqEuAnojYERGvAw8CcwvKzAXuSx6vAa6UpAzrZGZmQ5BlWEwCduU9350c67dMRBwG9gMTktemSuqS9CNJ7+vvAyTdIKlTUueePXvKW3szMzuqVju4fwNMiYiZwE3AtyWNLSwUESsjoi0i2iZOLHrLzczMhijLsOgFJuc9Pys51m8ZSSOAU4G9EXEoIvYCRMQmYDvwzgzramZmg8gyLDYC0yVNlTQKWAh0FJTpAK5NHs8HnoiIkDQx6SBH0juA6cCODOtqZmaDyGw0VEQclnQj8Di5obPfiIgtkm4DOiOiA7gXuF9SD/ACuUABeD9wm6Q3gDeBT0bEC1nV1cysHl14y2O8dOjI0edjR7fw7K1zMvksNcpOeW1tbeGhs2bWLAqDok+pgSFpU0S0FStXqx3cZmY2iP6CYrDjw5XpDG4zMxvY8vZuvrNhF0ciaJFYNGsyX5h3QbWr1S+HhZlZhtq7evnS49voffEALRJHIpg0bgytE8bw9Pa3umKPRPCtZ34FUJOB4dtQZmYZae/qZelD3fS+eADIBQJA74sHjgmKfN/ZsKvf44XGjm4p6fhwOSzMzDLypce3ceCN0voQjqQcdPTsrXOOC4YsR0P5NpSZWT/6bh/9+sUDnDluDDdfdS7zZhauWDS4XydXFKVoKWF5vKyCoT8OCzOzRH7/goC+f+P3vniApQ91A5QUGGeOG3P0FlRai2ZNLl6oCnwbysyM4/sXCm8GHXjjCF96fNvxJw7i5qvOZczI/vsQxoxs4fJzxh+9kmiRuObSKTXZuQ2+sjCzJjDQiKT8W0tp+hdKva2U/96DfXY9cFiYWUPJD4b+5I9Iyr+1lCYIzhw3puT6zJs5qa5CYSAOCzOre+1dvdz8vc288WZp5/XdWpo3c1LR/oUxI1u4+apzh1nT+uWwMLO6s7y9m29v+BVvlmFpu74ripuvOpelD3Ufcyuqr5O7Hm8blZvDwszqQq4D+lkOlHr5UETfraX8/oXhDJdtVA4LM6s5WQVDocJbS43Sv5AFh4WZ1Yz2rl5ufWQL+157I7PPqOcRSdXksDCzqmnv6uVzDz3LaxlfQdT6iq71wGFhZhXT3tXLsoe7efX1bPZcyHf5OeN54PrLMv+cZuGwMLNMFZv3UA5jRp7AHR+/0LeUMuSwMLOyyl+Ab9xJI3nl4GHeKMcY1wICFtfw8hiNxmFhZmWxvL376OY9fbLoqHbHdHU4LMxsSNq7evnM6s2Zvb8EEQ6HWuGwMLOSzL7rKZ77/auZvf9JI0/gv7v/oeY4LMxsUO1dvazo2MKLB8p/S8n9DvXDYWFm/Wpd8mjZ3mvcmJHsP/CGl9CoYw4LMzvqvGVrOXikvCOXPN+hMTgszJpce1cvN63eTBZzqGt55zcrjcPCrEll1VHtpTUak8PCrIlk2Vk9dnQLz946p+zva7XBYWHWBGbdvo7fvfx6Wd9z9IgTuPMTHuLaLBwWZg2qvxnV5eB+iObksDBrMItXrefp7S+U9T0dEOawMGsA5V5644xTRrFh2eyyvZ/Vv0zDQtIc4CtAC/D1iPhiweujgW8CFwN7gQURsTPv9SnAVmBFRPzPLOtqVo/KfRUxbsxIVlx9vvsh7DiZhYWkFuAeYDawG9goqSMituYVuw7YFxHTJC0E7gQW5L1+F/DDrOpoVo+yuM3kiXNWTJZXFpcAPRGxA0DSg8BcclcKfeYCK5LHa4CvSlJEhKR5wPNAdiuWmdWR9q5ebv7eZsq1A6kX7LNSZBkWk4Bdec93A7MGKhMRhyXtByZIOgh8ltxVyd8O9AGSbgBuAJgyZUr5am5WQ9q7ern1kS1l2xviywsuckBYyWq1g3sFcHdEvCJpwEIRsRJYCdDW1lb+rbjMqqTce1V7NJMNV5Zh0QtMznt+VnKsvzK7JY0ATiXX0T0LmC/pfwDjgDclHYyIr2ZYX7OqK3d/xM4vfrRs72XNLcuw2AhMlzSVXCgsBP6ioEwHcC2wHpgPPBERAbyvr4CkFcArDgprdA4Kq2WZhUXSB3Ej8Di5obPfiIgtkm4DOiOiA7gXuF9SD/ACuUAxawrlXqfJI5osS8r9Q77+tbW1RWdnZ7WrYZZKuVZ89QqvNlySNkVEW7FytdrBbdZwynElcfKoFm7/sws8mskqzmFhVgHlWI7DI5qsmhwWZhlq7+rlS49vo/fFA8N6H3dWW7U5LMwy0N7Vy39Z8zNeL8N+1g4KqwUOC7MyKsfw1+mnn8y6m64oT4XMysRhYVYGw123ScDdXobDapjDwmyY2rt6WfpQ95CDwh3XVg8cFmYlKtccCU+is3risDArwXCDwgFh9cphYZbS4lXrhxwUEiye5dtNVr8cFmYpDHWUk68krFE4LMwG0N7Vy+ceepbXSuy5njRuDDdfda5HNllDcViYFRjOznQe2WSNymFhlmd5ezcPPPMrSp137dVfrdE5LMwY+tWE+ySsWTgsrOkN9WrCQWHNxGFhTWl5ezff2bCLI0PY/Mv9EtaMHBbWdJa3d/OtZ35V8nmeK2HNzGFhTec7G3alLitgsa8kzBwW1nzS3noaN2YkK64+3/MlzHBYWBNqkQYNDF9NmB3vhGpXwKzSFs2aPOBrk8aN4e4FFzkozAr4ysKaTl8Q9I2G8oQ6s+IUQxg6WIva2tqis7Oz2tUwM6srkjZFRFuxcr6ysIaQP2/CVwpm5eewsLpXOG/iSMTR5w4Ms/JwB7fVvYHmTZQyn8LMBuewsLo30DDYoSzlYWb9c1hY3WuRSjpuZqVzn4XVvPauXlZ0bOHFA7nlw99+0khu+dhbM6sXzZrc71pPg82nMLPSOCysZrV39fLXqzcft3T4vtfe4OY1PwNg3sxJnjdhVgGeZ2E1qb2rl8+s3jxomUnjxvD0kg9WqEZmjals8ywkTQQ+C8wATuw7HhFF/y+VNAf4CtACfD0ivljw+mjgm8DFwF5gQUTslHQJsLKvGLAiIh4u9nnWGGbf9RTP/f7VouV+/eKBCtTGzCBdB/cDwM+BqcCtwE5gY7GTJLUA9wAfJhc0iyTNKCh2HbAvIqYBdwN3Jsf/BWiLiIuAOcDXJPmWWRNoXfJoqqAAOHPcmIxrY2Z90vwCnhAR90r6q4j4EfAjSUXDArgE6ImIHQCSHgTmAlvzyswFViSP1wBflaSIeC2vzIlQ8o6XVmemLX2UwyV8yyNbxM1XnZtdhczsGGmuLPp2sP+NpI9KmgmMT3HeJCB/VtTu5Fi/ZSLiMLAfmAAgaZakLUA38Mnk9WNIukFSp6TOPXv2pKiS1aLWJaUFhYAvzf9j7zNhVkFpriy+IOlU4G+AvwPGAp/JtFZARGwAzpf0R8B9kn4YEQcLyqwk6dtoa2vz1UedSdOJXeiMU0axYdnsjGpkZgNJExb7ImI/uX/1fwBA0uUpzusF8ge6n5Uc66/M7qRP4lRyHd1HRcTPJb0CvAvwcKcGUeo+2NNPP5l1N12RXYXMbFBpbkP9XcpjhTYC0yVNlTQKWAh0FJTpAK5NHs8HnoiISM4ZASDpbOA8ch3r1gBalzxaUlAADgqzKhvwykLSZcB7gImSbsp7aSy5obCDiojDkm4EHk/KfyMitki6DeiMiA7gXuB+ST3AC+QCBeC9wBJJbwBvAp+KiD+U3jyrJUO57XRii/jF7R/JqEZmltZgt6FGAW9LypySd/wlclcBRUXEWmBtwbHP5z0+CPx5P+fdD9yf5jOsPqSdO5Fv5xc/mlFtzKxUA4ZF3jDZf4yIX0o6qWBIq1lRQwkJcFCY1Zo0fRZnStoK/AJA0h9L+vtsq2WNoNSgEHDNpVMcFGY1KM1oqC8DV5F0TkfEzyS9P9NaWd1bvGp9yVcUzzskzGpWqiU0ImKXjt0b4Eg21bFGcOEtj/HSodL+ivhqwqy2pQmLXZLeA4SkkcBfkVsryuwYpc6dAM+fMKsXacLik+RWjp1EbhLd/wU+nWWlrP6ct2wtB4+UNon+mkuneM8JszpRNCyS+Q2LK1AXq1OtSx4tqbxDwqz+pN3P4nqgNb98RPxldtWyelDqbafLzxnPA9dflmGNzCwraW5D/R/gx8A/4Y5tS5Q6LHbs6BYHhVkdSxMWJ0XEZzOvidWNWbev43cvv566/Ikt4tlb52RYIzPLWppJeT+Q5MV5DMh1ZJcSFGNHt3htJ7MGMNhCgi+T26FOwOckHSK3EZKAiIixlami1YoLb3mspBFP7sg2axyD3YZ6V0T8smI1sZpWyoinEYKeOzzJzqyRDHYb6uGK1cJqWilBMXZ0i4PCrAENdmWhQV6zJlBqR7ZnY5s1rsHCYpKk/zXQixHxnzOoj9WIUmdkew6FWWMbLCwOAJsqVRGrHaV2ZDsozBrfYGGxNyLuq1hNrCaUunTHlxdcxLyZkzKqjZnVisHCIv3NamsIpQaFlxU3ax4DjoaKiEsrWRGrrlKCwrvZmTWfNDO4rcGVEhQCT7Qza0IOiyZXalB461Oz5pRqW1VrPItXrefp7S+kLu85FGbNbUhXFpJ+UO6KWOU4KMysVEO9DXV9WWthFbO8vbukoLj8nPEOCjMbPCwktUh6oPB4RPwmuypZVkrd2e6aS6d4sp2ZAUX6LCLiiKSzJY2KCM+7qGPtXb0lBYWHxppZvjQd3DuApyV1AEf30YyIuzKrlZVVqes8OSjMrFCasNie/JwAnJJtdazcpi19lMPpc8JBYWb9KhoWEXFrJSpi5bd41frUQXFii7z9qZkNqGhYSHqE3Paq+fYDncDXIuJgFhWz4Sll1JODwsyKSTN0dgfwCrAq+XkJeBl4Z/J8QJLmSNomqUfSkn5eHy1pdfL6BkmtyfHZkjZJ6k7+/GBpzWpui1etT92Z7aAwszTS9Fm8JyLenff8EUkbI+LdkrYMdJKkFuAeYDawG9goqSMituYVuw7YFxHTJC0E7gQWAH8APhYRv5b0LuBxwOtgp1DKhDvvQ2FmaaW5snibpCl9T5LHb0ueDjac9hKgJyJ2JMNuHwTmFpSZC/TtmbEGuFKSIqIrIn6dHN8CjJE0OkVdm9p5y9amDgqBg8LMUktzZfE3wE8kbSf3O2Yq8ClJJ/PWL/r+TAJ25T3fDcwaqExEHJa0H5hA7sqizyeAn0bEoRR1bVql7m7nBQHNrBRpRkOtlTQdOC85tC2vU/vLmdUMkHQ+uVtTfzLA6zcANwBMmTKlvyJN46VDR1KX9fBYMytVmtFQHy84dE5yBdAdEb8f5NReYHLe87OSY/2V2S1pBHAqsDf53LOAh4H/EBHb+/uAiFgJrARoa2srYTZB4yh1CQ8HhZkNRZrbUNcBlwFPkLsNdQWwCZgq6baIuH+A8zYC0yVNJRcKC4G/KCjTAVwLrAfmA09EREgaBzwKLImIp0trUvNwUJhZpaTp4B4B/FFEzI+ITwAzyM27mAV8dqCTIuIwcCO5kUw/B74bEVsk3Sbp6qTYvcAEST3ATUDf8NobgWnA5yVtTn5OH0L7GpqDwswqRRGD372RtDUiZuQ9F7AlImZI6oqImVlXMo22trbo7OysdjUqppQd7hwUZjYQSZsioq1YuTS3oZ5KNjv6XvJ8PvCjZDTUi8Ooow2Rg8LMKi1NWHwa+Djw3uT5fRGxJnn8gUxqZQOattRBYWaVl2bobADfT36Q9D5J90TEp7OunB3rwlseS70woIPCzMopzZUFkmYCi4B/BzwPPJRlpex4U5c8etxqjgNxUJhZuQ0YFpLeSS4gFpGbUb2aXIe4bz1V2LSlDgozq67Brix+AfwY+NOI6AGQ9NcVqZUdtby9O/WtpzNOGZVtZcysaQ02z+LjwG+AJyWtknQluUl5VkGlzKXYsGx2hjUxs2Y2YFhERHtELCS3JtSTwGeA0yX9b0n9rtVk5TXVQ2TNrEYUncEdEa9GxLcj4mPk1nfqYpCZ21Yere7QNrMakma5j6MiYl9ErIyIK7OqkMGs29elKjdCDgozq4ySwsKyt7y9m9+9PNieUm/pucNBYWaV4bCoMWk7tH1FYWaV5LCoIWnXfHJQmFmlOSxqRNqgGOHBy2ZWBQ6LGpA2KIT7KcysOhwWVXbesrWpyo0d3cLzvv1kZlXisKii5e3dHDySbjbFs7fOybg2ZmYDc1hUSXtXr0c+mVndcFhUyWdWb05V7ssLLsq4JmZmxTksqiBtP8Xl54xn3sxJGdfGzKw4h0WFzbp9Xap+CgEPXH9Z9hUyM0vBYVFBpSzl4ZFPZlZLHBYV5A5tM6tXDosKSTvx7vJzxmdcEzOz0jksKiDtJkYntsj9FGZWkxwWGWvv6k21iZGAX9z+kayrY2Y2JA6LjKWdT+EObTOrZQ6LDKWdT+GJd2ZW6xwWGVm8an2q+RTTTz/ZE+/MrOY5LDLy9PYXipY545RRrLvpiuwrY2Y2TA6LDCxv705VbsOy2RnXxMysPBwWGUgz+e6aS6dUoCZmZuWRaVhImiNpm6QeSUv6eX20pNXJ6xsktSbHJ0h6UtIrkr6aZR3LLc3ku2suncIX5l1QgdqYmZVHZmEhqQW4B/gwMANYJGlGQbHrgH0RMQ24G7gzOX4Q+K/A32ZVvyxceMtjqco5KMys3mR5ZXEJ0BMROyLideBBYG5BmbnAfcnjNcCVkhQRr0bET8iFRl1YvGo9Lx06UrSch8maWT3KMiwmAbvynu9OjvVbJiIOA/uBCRnWKTNpRj+NHd3iYbJmVpfquoNb0g2SOiV17tmzp2r1mH3XU6nKeR9tM6tXWYZFLzA57/lZybF+y0gaAZwK7E37ARGxMiLaIqJt4sSJw6zu0D33+1eLlvHoJzOrZ1mGxUZguqSpkkYBC4GOgjIdwLXJ4/nAExGRZt29mpFm9NMIuVPbzOrbiKzeOCIOS7oReBxoAb4REVsk3QZ0RkQHcC9wv6Qe4AVygQKApJ3AWGCUpHnAn0TE1qzqOxRp137qucOLBJpZfcssLAAiYi2wtuDY5/MeHwT+fIBzW7OsWzmkWfvJt5/MrBHUdQd3Nfn2k5k1E4fFEMy6fV2qcr79ZGaNwmExBL97+fWiZXZ6MyMzayAOixK1dxWO/j3eCFWgImZmFeSwKFGabVJ9+8nMGo3DogSLV60vWubEFl9WmFnjcViUIM36T7+4/SMVqImZWWU5LFJKs/udV5Q1s0blsEhheXt3qt3vvKKsmTUqh0UKaYLCVxVm1sgcFkWkWX58+ukn+6rCzBqaw6KINMuPr7vpiuwrYmZWRQ6LQaTd1MjMrNE5LAaR5qrCy3qYWTNwWAwgzbIenoBnZs3CYTGANMt6eAKemTULh0U/0o6AMjNrFg6LfngElJnZsRwWBdLuq21m1kwcFgXS7KvtEVBm1mwcFnnSLBZ4+TnjK1ATM7Pa4rDIk2YNqAeuv6wCNTEzqy0OixKcccqoalfBzKwqHBaJ1iWPFi2zYdnsCtTEzKz2OCzMzKwohwXprio8AsrMmpnDwszMinJYpODhsmbW7Jo+LNLcgvJwWTNrdk0fFmZmVlxTh0WaPSvcsW1m1uRhkWbPCjMzyzgsJM2RtE1Sj6Ql/bw+WtLq5PUNklrzXluaHN8m6aos6zmQay6dUo2PNTOrOZmFhaQW4B7gw8AMYJGkGQXFrgP2RcQ04G7gzuTcGcBC4HxgDvD3yftV1BfmXVDpjzQzq0lZXllcAvRExI6IeB14EJhbUGYucF/yeA1wpSQlxx+MiEMR8TzQk7yfmZlVQZZhMQnYlfd8d3Ks3zIRcRjYD0xIeW6m3LFtZvaWuu7glnSDpE5JnXv27Kl2dczMGlaWYdELTM57flZyrN8ykkYApwJ7U55LRKyMiLaIaJs4cWIZq25mZvmyDIuNwHRJUyWNItdh3VFQpgO4Nnk8H3giIiI5vjAZLTUVmA78c4Z1NTOzQWQWFkkfxI3A48DPge9GxBZJt0m6Oil2LzBBUg9wE7AkOXcL8F1gK/AY8OmIOFLuOg7UL+H+CjOzYyn3D/n619bWFp2dndWuhplZXZG0KSLaipWr6w5uMzOrDIeFmZkV5bAwM7OiHBZmZlaUw8LMzIpqmNFQkvYAvxzGW5wG/KFM1aklbld9adR2QeO2rd7bdXZEFJ3V3DBhMVySOtMMH6s3bld9adR2QeO2rVHbVci3oczMrCiHhZmZFeWweMvKalcgI25XfWnUdkHjtq1R23UM91mYmVlRvrIwM7OiGj4sJM2RtE1Sj6Ql/bw+WtLq5PUNklrzXluaHN8m6apK1ruYobZLUqukA5I2Jz//UOm6F5Oibe+X9FNJhyXNL3jtWknPJT/XFp5bTcNs15G876xwqf+qStGumyRtlfSspP8n6ey812r2+4Jht61mv7MhiYiG/QFagO3AO4BRwM+AGQVlPgX8Q/J4IbA6eTwjKT8amJq8T0u121SGdrUC/1LtNgyzba3AhcA3gfl5x8cDO5I/3548fnu12zTcdiWvvVLtNgyjXR8ATkoe/6e8v4s1+30Nt221/J0N9afRrywuAXoiYkdEvA48CMwtKDMXuC95vAa4UpKS4w9GxKGIeB7oSd6vFgynXbWuaNsiYmdEPAu8WXDuVcC6iHghIvYB64A5lah0CsNpVy1L064nI+K15Okz5Ha+hNr+vmB4bWs4jR4Wk4Bdec93J8f6LRO5DZv2AxNSnlstw2kXwFRJXZJ+JOl9WVe2RMP5717v39lgTkz2m39G0rzyVm1YSm3XdcAPh3hupQ2nbVC739mQjKh2BazifgNMiYi9ki4G2iWdHxEvVbtiNqizI6JX0juAJyR1R8T2aleqFJKuAdqAf1vtupTbAG2r++8sX6NfWfQCk/Oen5Uc67eMpBHAqcDelOdWy5DbldxW2wsQEZvI3ZN9Z+Y1Tm84/93r/TsbUET0Jn/uAJ4CZpazcsOQql2SPgQsA66OiEOlnFtFw2lbLX9nQ1PtTpMsf8hdOe0g10Hd10F1fkGZT3NsR/B3k8fnc2wH9w5qp4N7OO2a2NcOch13vcD4areplLbllf1Hju/gfp5cZ+nbk8c10bZhtuvtwOjk8WnAcxR0tNZyu8j9ktwOTC84XrPfVxnaVrPf2ZD/e1S7AhX4wj8C/GvyhS5Ljt1G7l8BACcC3yPXgf3PwDvyzl2WnLcN+HC121KOdgGfALYAm4GfAh+rdluG0LZ3k7t//Cq5q8Ateef+ZdLmHuA/Vrst5WgX8B6gO/ll1Q1cV+22lNiufwJ+l/yd2wx01MP3NZy21fp3NpQfz+A2M7OiGr3PwszMysBhYWZmRTkszMysKIeFmZkV5bAwM7OiHBZmg5A0T1JIOq+Cn7lT0mnJ438j6UFJ2yVtkrRWUi1NorQm4bAwG9wi4CfJn2WRzKhPU07Aw8BTEXFORFwMLAXOKFddzNLyPAuzAUh6G7kJmR8AHomIc5PjJwBfBT5IbqG5N4BvRMQaSR8B7iI3se5pcpMh/1TSCuAccrPmfxURA4aPpJ3k1hm6EFgREe/PpoVm6fnKwmxgc4HHIuJfgb6FFwE+Tm7viRnAvwcuA5B0IvA1crP9Lya3tEq+GcCHBguKAu8CNg2rBWZl4rAwG9gicnsYkPzZ90v+vcD3IuLNiPgt8GRy/DxgR+T2PwH4TsH7dUTEgSwrbJYVL1Fu1g9J48ndZrpAUpDbNS0k3TyMt321xPJbgPlFS5lVgK8szPo3H7g/Is6OiNaImExuVdT3keuL+ISkEySdAVyRnLMNeEfePu4LhlmHJ4DRkm7oOyDpwhrcsMqagMPCrH+LyI1Eyvf95Pj3ya0OuxX4FrnVe/cnt5g+BTwmaRPwMrkdCo8j6UxJawf47BHAociNPvkz4EPJ0NktwB3Ab4fVMrMh8GgosyGQ9LaIeEXSBHJLwF8eEb/NOy7gHuC5iLi7hPedCGyOiFraXtTMVxZmQ/QDSZuBHwP/LenoBrg+Ob6F3O6EX0v7hpKuTt5vabkrazZcvrIwM7OifGVhZmZFOSzMzKwoh4WZmRXlsDAzs6IcFmZmVpTDwszMivr/AoT/XhDBQgsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "#plt.ylim((0, 1))\n",
    "#plt.xlim((0, 1))\n",
    "ax.set_xlabel(\"Aggr. JC\")\n",
    "ax.set_ylabel(\"Aggr. Theta\")\n",
    "ax.plot([sum(x)/len(x) for x in matrix], [sum(x)/len(x) for x in theta_matrix], 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### check in which cases theta is larger than JC, or does it always damp it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Theta: 0.991915148124266\n",
      "JC: 0.99065423\n",
      "Theta: 0.9926692844725535\n",
      "JC: 0.9915254\n",
      "Theta: 0.9889061819891763\n",
      "JC: 0.9871795\n",
      "Theta: 0.9889604055336504\n",
      "JC: 0.9872421\n",
      "Theta: 0.9895498742947397\n",
      "JC: 0.98792267\n",
      "Theta: 0.9927915923396584\n",
      "JC: 0.9916667\n",
      "6\n"
     ]
    }
   ],
   "source": [
    "import itertools\n",
    "total_larger = 0\n",
    "for i, j in itertools.product([x for x in range(matrix.shape[0])], [x for x in range(matrix.shape[0])]):\n",
    "    if theta_matrix[i][j] > matrix[i][j]:\n",
    "        total_larger += 1\n",
    "        print(\"Theta: \" + str(theta_matrix[i][j]))\n",
    "        print(\"JC: \" + str(matrix[i][j]))\n",
    "print(total_larger)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plotting the diff between both"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M: 0.8275879767398286\n",
      "m: -0.0017266819891763596\n",
      "avg: 0.3231151660891969\n",
      "med: 0.26179260480984234\n",
      "p95: 0.7507212539236987\n",
      "p90: 0.6753463256387479\n",
      "p5: 0.0757316684008584\n",
      "p75: 0.4858193135266591\n",
      "p85: 0.6122401335965182\n"
     ]
    }
   ],
   "source": [
    "diffs = []\n",
    "for i, j in itertools.product([x for x in range(matrix.shape[0])], [x for x in range(matrix.shape[0])]):\n",
    "    # ONLY COUNT NON-ZERO\n",
    "    if matrix[i][j] > 0:\n",
    "        diff = matrix[i][j] - theta_matrix[i][j]\n",
    "        diffs.append(diff)\n",
    "diffs = np.asarray(diffs)\n",
    "max_diff = np.max(diffs)\n",
    "min_diff = np.min(diffs)\n",
    "med_diff = np.percentile(diffs, 50)\n",
    "p95_diff = np.percentile(diffs, 95)\n",
    "p90_diff = np.percentile(diffs, 90)\n",
    "p5_diff = np.percentile(diffs, 5)\n",
    "p75_diff = np.percentile(diffs, 75)\n",
    "p85_diff = np.percentile(diffs, 85)\n",
    "avg_diff = np.average(diffs)\n",
    "print(\"M: \" + str(max_diff))\n",
    "print(\"m: \" + str(min_diff))\n",
    "print(\"avg: \" + str(avg_diff))\n",
    "print(\"med: \" + str(med_diff))\n",
    "print(\"p95: \" + str(p95_diff))\n",
    "print(\"p90: \" + str(p90_diff))\n",
    "print(\"p5: \" + str(p5_diff))\n",
    "print(\"p75: \" + str(p75_diff))\n",
    "print(\"p85: \" + str(p85_diff))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For a given threshold, which columns make it using JC that do not make it using Theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def candidates(matrix, threshold):\n",
    "    pairs = []\n",
    "    for i, row in enumerate(matrix):\n",
    "        for j, el in enumerate(row):\n",
    "            if el > threshold and i != j:\n",
    "                pairs.append((i, j))\n",
    "    return pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Threshold: 0.5\n",
      "JC: 111394\n",
      "TH: 2072\n",
      "Threshold: 0.6\n",
      "JC: 88787\n",
      "TH: 1021\n",
      "Threshold: 0.7\n",
      "JC: 67744\n",
      "TH: 453\n",
      "Threshold: 0.8\n",
      "JC: 50823\n",
      "TH: 170\n",
      "Threshold: 0.9\n",
      "JC: 33489\n",
      "TH: 14\n"
     ]
    }
   ],
   "source": [
    "for th in [0.5, 0.6, 0.7, 0.8, 0.9]:\n",
    "    jc_pairs = candidates(matrix, th)\n",
    "    th_pairs = candidates(theta_matrix, th)\n",
    "    print(\"Threshold: \" + str(th))\n",
    "    print(\"JC: \" + str(len(jc_pairs)))\n",
    "    print(\"TH: \" + str(len(th_pairs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "JC pairs total: 67744\n",
      "TH pairs total: 453\n"
     ]
    }
   ],
   "source": [
    "jc_pairs = candidates(matrix, 0.7)\n",
    "th_pairs = candidates(theta_matrix, 0.7)\n",
    "print(\"JC pairs total: \" + str(len(jc_pairs)))\n",
    "print(\"TH pairs total: \" + str(len(th_pairs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'All_olap_descriptor_types.csv->Descriptor Type ~ Sdo_crs_engineering.csv->Coord Ref Sys Name'"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# bad\n",
    "name_pair(jc_pairs[0], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'All_olap_descriptor_types.csv->Descriptor Type ~ Sdo_coord_op_methods.csv->Coord Op Method Name'"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# bad\n",
    "name_pair(jc_pairs[1], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'All_olap_descriptor_types.csv->Descriptor Type ~ short_library_reserve_catalog.csv->Catalog Publisher'"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pair(jc_pairs[2], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Fclt_building_list.csv->Building Type ~ Fac_building.csv->Building Type'"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pair(th_pairs[0], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Fclt_building_list.csv->Building Type ~ Fclt_building_hist_1.csv->Building Type'"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pair(th_pairs[1], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Fclt_building_list.csv->Building Type ~ Fclt_building.csv->Building Type'"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pair(th_pairs[2], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Fclt_building_hist_1.csv->Access Level Code ~ Fclt_building.csv->Access Level Code'"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pair(th_pairs[3], idmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "195\n"
     ]
    }
   ],
   "source": [
    "c = 0\n",
    "for i, j in th_pairs:\n",
    "    iname = idmap[i].split(\"->\")[1]\n",
    "    jname = idmap[j].split(\"->\")[1]\n",
    "    if iname != jname:\n",
    "        c += 1\n",
    "        #print(name_pair((i,j), idmap))\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing JC across ngram sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "im2 = \"temp/mit_matrix_n2_k64.dat\"\n",
    "im3 = \"temp/mit_matrix_n3_k64.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 = deserialize_matrix(im2)\n",
    "m3 = deserialize_matrix(im3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hypothesis 1: Non-zero JC values in m3 must always be smaller than non-zero JC values in m2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we need to keep in mind the 0.05 normal error from the lazo method, so we can relax that\n",
    "fail = 0\n",
    "for r2, r3 in zip(m2, m3):\n",
    "    for el2, el3 in zip(r2, r3):\n",
    "        if el2 > 0 and el3 > 0:\n",
    "            if (el3 - 0.05) > el2:  # D value of lazo method\n",
    "#                 print(str(el3) + \"-\" + str(el2))\n",
    "                fail += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_els = 1599*1599\n",
    "perc_fail = (float(fail) / float(total_els))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0015887040094242767"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.5 per thousand fail --> need to understand why"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "im3 = \"temp/mit_matrix_n3_k64.dat\"\n",
    "im4 = \"temp/mit_matrix_n4_k64.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "m3 = deserialize_matrix(im3)\n",
    "m4 = deserialize_matrix(im4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "fail = 0\n",
    "for r3, r4 in zip(m3, m4):\n",
    "    for el3, el4 in zip(r3, r4):\n",
    "        if el3 > 0 and el4 > 0:\n",
    "            if (el4 - 0.05) > el3:  # D value of lazo method\n",
    "#                 print(str(el4) + \"-\" + str(el3))\n",
    "                fail += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_els = 1599*1599\n",
    "perc_fail = (float(fail) / float(total_els))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0006867957263783924"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 0.6 per thousand fail, trend goes in the direction of supporting the hypothesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "im4 = \"temp/mit_matrix_n4_k64.dat\"\n",
    "im5 = \"temp/mit_matrix_n5_k64.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "m4 = deserialize_matrix(im4)\n",
    "m5 = deserialize_matrix(im5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m5.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "fail = 0\n",
    "for r4, r5 in zip(m4, m5):\n",
    "    for el4, el5 in zip(r4, r5):\n",
    "        if el4 > 0 and el5 > 0:\n",
    "            if (el5 - 0.05) > el4:  # D value of lazo method\n",
    "#                 print(str(el5) + \"-\" + str(el4))\n",
    "                fail += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_els = 1599*1599\n",
    "perc_fail = (float(fail) / float(total_els))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0002718240488798307"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 0.2 per thousand fail, trend goes in the direction of supporting the hypothesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "im5 = \"temp/mit_matrix_n5_k64.dat\"\n",
    "im6 = \"temp/mit_matrix_n6_k64.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "m5 = deserialize_matrix(im5)\n",
    "m6 = deserialize_matrix(im6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m5.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 1599)"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m6.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "fail = 0\n",
    "for r5, r6 in zip(m5, m6):\n",
    "    for el5, el6 in zip(r5, r6):\n",
    "        if el5 > 0 and el6 > 0:\n",
    "            if (el6 - 0.05) > el5:  # D value of lazo method\n",
    "#                 print(str(el6) + \"-\" + str(el5))\n",
    "                fail += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_els = 1599*1599\n",
    "perc_fail = (float(fail) / float(total_els))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0006957913423844875"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 0.6 per thousand fail, trend goes in the direction of supporting the hypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TESTING"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_zero_row = [el for el in m2[0] if el > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(non_zero_row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9895834,\n",
       " 0.16666667,\n",
       " 0.48958334,\n",
       " 0.5260416,\n",
       " 0.21875,\n",
       " 0.4635417,\n",
       " 0.47916666,\n",
       " 0.203125,\n",
       " 0.46875,\n",
       " 0.40625,\n",
       " 0.3541667,\n",
       " 0.5260416,\n",
       " 0.16666667,\n",
       " 0.17708333,\n",
       " 0.16666667,\n",
       " 0.5364584,\n",
       " 0.40104166,\n",
       " 0.15104166]"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "non_zero_row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "aggr = np.sum(non_zero_row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.89062506"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aggr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm = [el/aggr for el in non_zero_row]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxnorm = max(norm)\n",
    "minnorm = min(norm)\n",
    "beta_scores = [(el - minnorm) / (maxnorm - minnorm) for el in norm]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0,\n",
       " 0.01863355066856898,\n",
       " 0.4037266886678772,\n",
       " 0.44720485828171175,\n",
       " 0.08074534250376136,\n",
       " 0.372670822563943,\n",
       " 0.391304313605188,\n",
       " 0.06211180376065718,\n",
       " 0.37888196239342836,\n",
       " 0.30434780742101164,\n",
       " 0.24223605136221368,\n",
       " 0.44720485828171175,\n",
       " 0.01863355066856898,\n",
       " 0.03105590188032862,\n",
       " 0.01863355066856898,\n",
       " 0.4596273764499785,\n",
       " 0.29813661988966705,\n",
       " 0.0]"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_scores = [jc * beta for jc, beta in zip(non_zero_row, beta_scores)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.9895834, 0.9895834)\n",
      "(0.0031055918402066655, 0.16666667)\n",
      "(0.19765786068515945, 0.48958334)\n",
      "(0.2352483591782849, 0.5260416)\n",
      "(0.017663043672697797, 0.21875)\n",
      "(0.1727484666316885, 0.4635417)\n",
      "(0.1874999809937905, 0.47916666)\n",
      "(0.012616460138883491, 0.203125)\n",
      "(0.17760091987191953, 0.46875)\n",
      "(0.12364129676478598, 0.40625)\n",
      "(0.08579194293198572, 0.3541667)\n",
      "(0.2352483591782849, 0.5260416)\n",
      "(0.0031055918402066655, 0.16666667)\n",
      "(0.005499482521121854, 0.17708333)\n",
      "(0.0031055918402066655, 0.16666667)\n",
      "(0.24657096696655315, 0.5364584)\n",
      "(0.1195652049473411, 0.40104166)\n",
      "(0.0, 0.15104166)\n"
     ]
    }
   ],
   "source": [
    "for el in zip(final_scores, non_zero_row):\n",
    "    print(el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}