JaimieMurdock/LDA Model Comparison.ipynb

## LDA Model Comparison.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LDA Model Comparison\n",
    "This notebook presents a framework for comparing LDA models.\n",
    "\n",
    "This follows a 4-step process:\n",
    "1. Load in the two LDA models under investigation.\n",
    "2. Perform heirarchical clustering among all topics in the two models.\n",
    "3. Plot the topic alignment.\n",
    "4. Plot the topic alignment agreement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the two models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running from notebook, using serial load function.\n",
      "Loading LDA data from /tmp/ap/models/ap-nltk-en-freq5-N2000-LDA-K20-document-20.npz\n",
      "Running from notebook, using serial load function.\n",
      "Loading LDA data from /tmp/ap/models/ap-nltk-en-freq5-N2000-LDA-K20-document-20.npz\n"
     ]
    }
   ],
   "source": [
    "from vsm import *\n",
    "from configparser import ConfigParser\n",
    "\n",
    "def load_from_config(config_file, k):\n",
    "    \n",
    "    config = ConfigParser()\n",
    "    config.read(config_file)\n",
    "\n",
    "    # path variables\n",
    "    path = config.get('main', 'path')\n",
    "    context_type = config.get('main', 'context_type')\n",
    "    corpus_file = config.get('main', 'corpus_file')\n",
    "    model_pattern = config.get('main', 'model_pattern') \n",
    "\n",
    "    c = Corpus.load(corpus_file)\n",
    "    m = LdaCgsMulti.load(model_pattern.format(k))\n",
    "    v = LdaCgsViewer(c, m)\n",
    "    return v\n",
    "\n",
    "v1 = load_from_config('ap.ini', k=20)\n",
    "v2 = load_from_config('ap.ini', k=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells are more or less automatically importing everything you need to do the alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from vsm import *\n",
    "import numpy as np\n",
    "import itertools\n",
    "import copy\n",
    "from scipy.stats import spearmanr, rankdata\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "from random import randrange\n",
    "import random\n",
    "\n",
    "\n",
    "import os.path\n",
    "from configparser import ConfigParser\n",
    "\n",
    "def is_valid_filepath(parser, arg):\n",
    "    if not os.path.exists(arg):\n",
    "        parser.error(\"The file %s does not exist!\" % arg)\n",
    "    else:\n",
    "        return arg\n",
    "    \n",
    "\n",
    "def deep_subcorpus(labels):\n",
    "    # resolve labels to indexes\n",
    "    docs_labels = [v._res_doc_type(d) for d in labels]\n",
    "    docs, labels = zip(*docs_labels)\n",
    "    \n",
    "    # get lengths of all contexts\n",
    "    lens = np.array([len(ctx) for ctx in v.corpus.view_contexts('book')])\n",
    "    \n",
    "    # get the context_type index for use with context_data\n",
    "    ctx_idx = v.corpus.context_types.index(v.model.context_type)\n",
    "    \n",
    "    # get original slices\n",
    "    slice_idxs = [range(s.start,s.stop) for i, s in enumerate(v.corpus.view_contexts('book',as_slices=True)) \n",
    "                      if i in docs]\n",
    "    \n",
    "    new_corpus = copy.deepcopy(v.corpus)\n",
    "    # reduce corpus to subcorpus    \n",
    "    new_corpus.corpus = new_corpus.corpus[list(itertools.chain(*slice_idxs))]\n",
    "    \n",
    "    # reinitialize index fields\n",
    "    for i,d in enumerate(docs):\n",
    "        new_corpus.context_data[ctx_idx]['idx'][d] = lens[list(docs[:i+1])].sum()\n",
    "    \n",
    "    # reduce metadata to only the new subcorpus\n",
    "    new_corpus.context_data[ctx_idx] = new_corpus.context_data[ctx_idx][list(docs)]\n",
    "    \n",
    "    \n",
    "    \n",
    "    return new_corpus\n",
    "\n",
    "\n",
    "from vsm.spatial import *\n",
    "\n",
    "import numpy as np\n",
    "import scipy.cluster.hierarchy as sch\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "__all__ = ['model_dist','avg_log_likelihood','perplexity','model_stats','plot_topic_similarity','compare_models']\n",
    "\n",
    "def topic_overlap(v1,v2):\n",
    "    \"\"\"\n",
    "    Calculates the overlap of two corpora and recalculates normalized topic matricies\n",
    "    including only the overlapping words, ordered by the overlap order.\n",
    "    \n",
    "    Returns the joint vocabulary, v1.topics() and v2.topics() filtered and renormed.\n",
    "    \"\"\"\n",
    "    vocab = set(v1.corpus.words)\n",
    "    t1 = np.array(v1.topics())['value']\n",
    "    t2 = np.array(v2.topics())['value']\n",
    "    \n",
    "    if v1.corpus.words_int != v2.corpus.words_int:\n",
    "        print(\"corpus.words_int different, aligning words\")\n",
    "        vocab = vocab.intersection(v2.corpus.words)\n",
    "        print(\"preserving {}% of words in v1; \".format(100 * len(vocab) / float(len(v1.corpus.words))))\n",
    "        print(\"{}% of words in v2 \".format(100 * len(vocab) / float(len(v2.corpus.words))))\n",
    "        \n",
    "        t1 = t1[:,np.array([v1.corpus.words_int[word] for word in vocab])]\n",
    "        t1 = (t1.T / t1.sum(axis=1)).T\n",
    "\n",
    "        t2 = t2[:,np.array([v2.corpus.words_int[word] for word in vocab])]\n",
    "        t2 = (t2.T / t2.sum(axis=1)).T\n",
    "\n",
    "    return (vocab, t1, t2)\n",
    "\n",
    "def model_dist(v1,v2, dist_fn=JS_dist):\n",
    "    \"\"\"\n",
    "    Takes two LdaViewer objects and a distance metric and calculates the topic-topic distance.\n",
    "    \"\"\"\n",
    "    vocab, t1, t2 = topic_overlap(v1,v2)\n",
    "    #t1 = np.array(v1.topics())['value']\n",
    "    #t2 = np.array(v2.topics())['value']\n",
    "    combined = np.concatenate((t1,t2))\n",
    "    \n",
    "    # NOTE: Doing this by row to reduce memory requirements and time requirements\n",
    "    D = np.column_stack(np.lib.pad(dist_fn(combined[i:,:],row.T), \n",
    "                                   (i,0), 'constant', constant_values=0)\n",
    "                                if i + 1 < len(combined) else np.zeros(len(combined))\n",
    "                            for i,row in enumerate(combined))\n",
    "    return D + D.T - np.diag(D.diagonal())\n",
    "    # Old simple version:\n",
    "    # return np.column_stack(dist_fn(combined,row.T) for row in combined)\n",
    "\n",
    "def doc_overlap(v1,v2, context_type, norm=True):\n",
    "    context_label = context_type + '_label'\n",
    "    ids = np.intersect1d(v1.corpus.view_metadata(context_type)[context_label], \n",
    "                         v2.corpus.view_metadata(context_type)[context_label])\n",
    "    d1 = v1.doc_topic_matrix(ids)\n",
    "    d2 = v2.doc_topic_matrix(ids)\n",
    "\n",
    "    # renormalize so that each topic is now a document probability\n",
    "    if norm:\n",
    "        d1 = (d1 / d1.sum(axis=0))\n",
    "        d2 = (d2 / d2.sum(axis=0))\n",
    "    \n",
    "    # original d1 and d2 are doc_topic, switch to topic_doc\n",
    "    return (ids, d1.T, d2.T)\n",
    "\n",
    "def model_doc_dist(v1, v2, context_type, dist_fn=JS_dist):\n",
    "    ids, d1, d2 = doc_overlap(v1,v2, context_type)\n",
    "    combined = np.concatenate((d1,d2))\n",
    "\n",
    "    # NOTE: Doing this by row to reduce memory requirements and time requirements\n",
    "    D = np.column_stack(np.lib.pad(dist_fn(combined[i:,:],row.T),\n",
    "                                   (i,0), 'constant', constant_values=0)\n",
    "                                    if i + 1 < len(combined) else np.zeros(len(combined))\n",
    "                                for i,row in enumerate(combined))\n",
    "    return D + D.T - np.diag(D.diagonal())\n",
    "\n",
    "\n",
    "def pearson(v1, v2, context_type):\n",
    "    context_label = context_type + '_label'\n",
    "    ids, d1, d2 = doc_overlap(v1, v2, context_type)\n",
    "\n",
    "    r_all = []\n",
    "    for id in ids:\n",
    "        sim1 = v1.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim1['doc'], ids).reshape(sim1['doc'].shape)\n",
    "        sim1 = sim1[np.where(ix)]\n",
    "        sim1 = sim1[sim1['doc'].argsort()]\n",
    "        \n",
    "        sim2 = v2.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim2['doc'], ids).reshape(sim2['doc'].shape)\n",
    "        sim2 = sim2[np.where(ix)]\n",
    "        sim2 = sim2[sim2['doc'].argsort()]\n",
    "        \n",
    "        r, pval = pearsonr(sim1['value'], sim2['value'])\n",
    "        r_all.append(r)\n",
    "    \n",
    "    return sum(r_all)/len(r_all)\n",
    "\n",
    "def spearman(v1, v2, context_type):\n",
    "    context_label = context_type + '_label'\n",
    "    ids, d1, d2 = doc_overlap(v1, v2, context_type)\n",
    "\n",
    "    r_all = []\n",
    "    for id in ids:\n",
    "        sim1 = v1.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim1['doc'], ids).reshape(sim1['doc'].shape)\n",
    "        sim1 = sim1[np.where(ix)]\n",
    "        sim1 = sim1[sim1['doc'].argsort()]\n",
    "        \n",
    "        sim2 = v2.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim2['doc'], ids).reshape(sim2['doc'].shape)\n",
    "        sim2 = sim2[np.where(ix)]\n",
    "        sim2 = sim2[sim2['doc'].argsort()]\n",
    "\n",
    "        r, pval = spearmanr(rankdata(sim1['value']), rankdata(sim2['value']))\n",
    "        r_all.append(r)\n",
    "\n",
    "    return sum(r_all)/len(r_all)\n",
    "\n",
    "def recall(v1,v2, context_type,N=10):\n",
    "    context_label = context_type + '_label'\n",
    "    ids, d1, d2 = doc_overlap(v1, v2, context_type)\n",
    "\n",
    "    r_all = []\n",
    "    for id in ids:\n",
    "        sim1 = v1.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim1['doc'], ids).reshape(sim1['doc'].shape)\n",
    "        sim1 = sim1[np.where(ix)]\n",
    "\n",
    "        sim2 = v2.dist_doc_doc(id)\n",
    "        ix = np.in1d(sim2['doc'], ids).reshape(sim2['doc'].shape)\n",
    "        sim2 = sim2[np.where(ix)]\n",
    "\n",
    "        sim1 = np.array(sim1[1:N+2])['doc']\n",
    "        sim2 = np.array(sim2[1:N+2])['doc']\n",
    "        r = np.where(np.in1d(sim1, sim2))[0].size / float(N)\n",
    "        r_all.append(r)\n",
    "    \n",
    "    return sum(r_all)/len(r_all)\n",
    "\n",
    "def avg_log_likelihood(viewer):\n",
    "    \"\"\" Calculates the average log likelihood per token. \"\"\"\n",
    "    return viewer.model.log_probs[-1][1] / len(viewer.corpus.corpus)\n",
    "\n",
    "def perplexity(viewer):\n",
    "    \"\"\" Calculates the perplexity. \"\"\"\n",
    "    return np.exp(-1*avg_log_likelihood(viewer))\n",
    "\n",
    "def model_stats(*viewers):\n",
    "    \"\"\"\n",
    "    Prints a table of avg log likelihood and perplexity for each viewer.\n",
    "    \"\"\"\n",
    "    print(\"model\", \"k\", \"tokens\", \"types\", \"avg-log-likelihood\", \"perplexity\")\n",
    "    for i,v in enumerate(viewers):\n",
    "        print(\"M{}\".format(i), v.model.K, len(v.corpus.corpus), len(v.corpus.words), avg_log_likelihood(v), perplexity(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Cluster Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jaimie/anaconda3/lib/python3.6/site-packages/IPython/core/magics/pylab.py:160: UserWarning: pylab import has clobbered these variables: ['copy', 'angle', 'random', '__version__']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "import scipy\n",
    "import scipy.cluster.hierarchy as sch\n",
    "\n",
    "def create_dendogram(D, xdim=None, method='ward', show=True):\n",
    "    # create_dendogram\n",
    "    fig = pylab.figure(figsize=(3,9))\n",
    "    Y = sch.linkage(D, method=method)\n",
    "    Z = sch.dendrogram(Y, orientation='right')\n",
    "    \n",
    "    # create plot\n",
    "    if show:\n",
    "        ax = fig.gca()\n",
    "        ax.set_xticks([])\n",
    "        if xdim is not None:\n",
    "            ax.set_yticklabels([\"M2 \" + str(idx - xdim) if idx >= xdim else \"M1 \" + str(idx)\n",
    "                                for idx in np.array(Z['leaves'])])\n",
    "        fig.show()\n",
    "    else:\n",
    "        pylab.close()\n",
    "    return Z\n",
    "\n",
    "import pylab\n",
    "def plot_dendogram(D, Z, xdim,ydim, filter_axis=False, show_self=False):\n",
    "    # Compute and plot dendrogram.\n",
    "    fig = pylab.figure(figsize=(12,10))\n",
    "    \n",
    "    D = np.copy(D)\n",
    "    \n",
    "    # Now make plot better\n",
    "    if not show_self:\n",
    "        D[:xdim,:xdim] = 0\n",
    "        D[xdim:,xdim:] = 0\n",
    "\n",
    "    # Plot distance matrix.\n",
    "    axmatrix = fig.add_axes([0.3,0.1,0.6,0.8])\n",
    "    index = np.array(Z['leaves'])\n",
    "    if filter_axis:\n",
    "        D = D[index[index < xdim],:]\n",
    "        D = D[:,index[index >= xdim]]\n",
    "    else:\n",
    "        D = D[index,:]\n",
    "        D = D[:,index]\n",
    "    im = axmatrix.imshow(D, interpolation='none', cmap=mpl.cm.Greys_r, vmax=0.25, vmin=0.0)#aspect='auto', origin='lower')\n",
    "    if filter_axis:\n",
    "        axmatrix.set_xticks(arange(ydim))\n",
    "        axmatrix.set_xticklabels([str(i -xdim) for i in index[index >= xdim]])\n",
    "        axmatrix.set_yticks(arange(xdim))\n",
    "        axmatrix.set_yticklabels([str(i) for i in index[index < xdim]])\n",
    "    else:\n",
    "        axmatrix.set_xticks(arange(xdim+ydim))\n",
    "        axmatrix.set_xticklabels(index)\n",
    "        axmatrix.set_yticks(arange(xdim+ydim))\n",
    "        axmatrix.set_yticklabels(index)\n",
    "    \n",
    "    # create an axes on the right side of ax. The width of cax will be 5%\n",
    "    # of ax and the padding between cax and ax will be fixed at 0.05 inch.\n",
    "    divider = make_axes_locatable(axmatrix)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.15)\n",
    "\n",
    "    fig.colorbar(im, cax=cax)\n",
    "\n",
    "    # Display and save figure.\n",
    "    fig.show()\n",
    "\n",
    "\n",
    "def plot_dendogram(D, Z, xdim,ydim, dist=None, filter_axis=False, show_self=False, alignment=None):\n",
    "    # Compute and plot dendrogram.\n",
    "    fig = figure(figsize=(12,10))\n",
    "    \n",
    "    D = np.copy(D)\n",
    "    \n",
    "    # Now make plot better\n",
    "    if not show_self:\n",
    "        D[:xdim,:xdim] = 0\n",
    "        D[xdim:,xdim:] = 0\n",
    "\n",
    "    # Plot distance matrix.\n",
    "    axmatrix = fig.add_axes([0.3,0.1,0.6,0.8])\n",
    "    index = np.array(Z['leaves'])\n",
    "    if filter_axis:\n",
    "        D = D[index[index < xdim],:]\n",
    "        D = D[:,index[index >= xdim]]\n",
    "    else:\n",
    "        D = D[index,:]\n",
    "        D = D[:,index]\n",
    "        \n",
    "    # generate palette\n",
    "    palette=cm.Blues_r\n",
    "    palette.set_bad(alpha=0.0)\n",
    "    MAX = np.sort(D.flatten())[int(.1*D.shape[0]*D.shape[1])]\n",
    "    MAX = np.max(np.diagonal(dist[:,dist.argsort(axis=0)[1]]))\n",
    "    MAX *= 1.1\n",
    "    im = axmatrix.imshow(D, interpolation='none', cmap=palette, vmax=MAX, vmin=0.0)#aspect='auto', origin='lower')\n",
    "    if filter_axis:\n",
    "        axmatrix.set_xticks(arange(ydim))\n",
    "        axmatrix.set_xticklabels([str(i - xdim) for i in index[index >= xdim]])\n",
    "        axmatrix.set_yticks(arange(xdim))\n",
    "        axmatrix.set_yticklabels([str(i) for i in index[index < xdim]])\n",
    "    else:\n",
    "        axmatrix.set_xticks(arange(xdim+ydim))\n",
    "        axmatrix.set_xticklabels(index)\n",
    "        axmatrix.set_yticks(arange(xdim+ydim))\n",
    "        axmatrix.set_yticklabels(index)\n",
    "    \n",
    "    if alignment is not None:\n",
    "        axmatrix.autoscale(False)\n",
    "        ys,xs = zip(*alignment)\n",
    "        xindex = index[index >= xdim] - xdim\n",
    "        yindex = index[index < xdim]\n",
    "        \n",
    "        xs = np.array([np.squeeze(np.where(xindex == x)) for x in xs])\n",
    "        ys = np.array([np.squeeze(np.where(yindex == y)) for y in ys])\n",
    "        \n",
    "        axmatrix.scatter(xs, ys, marker='o', s=125, color='w', lw=4, edgecolor='k')\n",
    "    \n",
    "    title(\"Jensen-Shannon Distance from topic to topic\")\n",
    "    ylabel(\"Model 1\")\n",
    "    xlabel(\"Model 2\")\n",
    "    # Plot colorbar.\n",
    "    #axcolor = fig.add_axes([0.91,0.1,0.02,0.8])\n",
    "    divider = make_axes_locatable(axmatrix)\n",
    "    cax = divider.append_axes(\"right\", size=\"3%\", pad=0.1)\n",
    "    colorbar(im, cax=cax, extend='max')\n",
    "\n",
    "def dendogram_demo():\n",
    "    # Generate features and distance matrix.\n",
    "    x = scipy.rand(40)\n",
    "    D = scipy.zeros([40,40])\n",
    "    for i in range(40):\n",
    "        for j in range(40):\n",
    "            D[i,j] = abs(x[i] - x[j])\n",
    "    \n",
    "\n",
    "    gram = create_dendogram(D, show=False)\n",
    "    #dendo = plot_dendogram(D,gram,10,30, show_self=True)\n",
    "    dendo = plot_dendogram(D,gram,10,30, D, show_self=False)\n",
    "    #dendo = plot_dendogram(D,gram,10,30, show_self=True, filter_axis=True)\n",
    "    dendo = plot_dendogram(D,gram,10,30, D, show_self=False, filter_axis=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Model comparison code\n",
    "\n",
    "def plot_topic_similarity(v1, v2, dist=None, dist_fn=JS_dist, sorted=True, alignment=None):\n",
    "    # Calculate distance between all topics\n",
    "    if dist is None:\n",
    "        dist = model_dist(v1, v2, dist_fn)\n",
    "    # print dist\n",
    "\n",
    "    if sorted:\n",
    "        dendo = create_dendogram(dist, xdim=v1.model.K, method='ward') \n",
    "        plot_dendogram(dist, dendo, v1.model.K, v2.model.K, dist=dist, filter_axis=True, alignment=alignment)\n",
    "    else:\n",
    "        # plot distance heatmap\n",
    "        figure(figsize=(12,10))\n",
    "        \n",
    "        # filtering axis\n",
    "        xdim = v1.model.K\n",
    "        index = np.arange(0,len(dist))\n",
    "        D = np.copy(dist)\n",
    "        D = D[index[index < xdim],:]\n",
    "        D = D[:,index[index >= xdim]]\n",
    "        \n",
    "        # generate palette\n",
    "        palette=cm.Blues_r\n",
    "        palette.set_bad(alpha=0.0)\n",
    "        MAX = np.sort(D.flatten())[int(.1*D.shape[0]*D.shape[1])]\n",
    "        MAX = np.max(np.diagonal(dist[:,dist.argsort(axis=0)[1]]))\n",
    "        MAX *= 1.1\n",
    "        \n",
    "        # plot data\n",
    "        im = imshow(D, cmap=palette, vmax=MAX, vmin=0.0, interpolation='none')\n",
    "        if alignment is not None:\n",
    "            ax = gca()\n",
    "            ax.autoscale(False)\n",
    "            ys,xs = zip(*alignment)\n",
    "            scatter(xs,ys, marker='x', s=200, color='w', lw=5)\n",
    "        \n",
    "        # label data\n",
    "        xticks(np.arange(D.shape[1]))\n",
    "        yticks(np.arange(D.shape[0]))\n",
    "        title(\"Jensen-Shannon Distance from topic to topic\")\n",
    "        xlabel(\"Model 1\")\n",
    "        ylabel(\"Model 2\")\n",
    "        \n",
    "        # create heatmap\n",
    "        divider = make_axes_locatable(gca())\n",
    "        cax = divider.append_axes(\"right\", size=\"3%\", pad=0.1)\n",
    "        colorbar(im, cax=cax, extend='max')\n",
    "\n",
    "def compare_models(v1, v2, context_type='document', dist_fn=JS_dist, sorted=True):\n",
    "    model_stats(v1,v2)\n",
    "    plot_topic_similarity(v1,v2,dist_fn=dist_fn,sorted=sorted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Create asymmetric model comparison matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def alignment_fitness(topic_pairs, v1, v2, dist=None, dist_fn=JS_dist):\n",
    "    \"\"\"\n",
    "    Takes a list of topic pair tuples and returns the sum of the JS_dist between them\n",
    "    \"\"\"    \n",
    "    if dist is None:\n",
    "        dist = model_dist(v1, v2, dist_fn)\n",
    "    if dist.shape[0] == (v1.model.K + v2.model.K):\n",
    "        dist = filter_dist(v1, v2, dist)\n",
    "    \n",
    "    return sum([dist[t[0]][t[1]] for t in topic_pairs])\n",
    "\n",
    "def filter_dist(v1,v2,dist):\n",
    "    xdim = v1.model.K\n",
    "    index = np.arange(0,len(dist))\n",
    "    D = np.copy(dist)\n",
    "    D = D[index[index < xdim],:]\n",
    "    D = D[:,index[index >= xdim]]\n",
    "    return D\n",
    "\n",
    "def plot_alignment(v1, v2, dist, alignment=None, fn_name=None):\n",
    "    # Calculate distance between all topics\n",
    "    if dist is None:\n",
    "        dist = model_dist(v1, v2, dist_fn)\n",
    "    if dist.shape[0] == (v1.model.K + v2.model.K):\n",
    "        dist = filter_dist(v1, v2, dist)\n",
    "    \n",
    "    Xs = cm.jet_r(dist)\n",
    "    \n",
    "    if alignment is None:\n",
    "        alpha = 1.0\n",
    "    else:\n",
    "        alpha = np.zeros(dist.shape)\n",
    "        alpha[zip(*alignment)] = 1\n",
    "    \n",
    "    Xs[:,:,3] = alpha\n",
    "    \n",
    "    # plot distance heatmap\n",
    "    figure(figsize=(12,10))\n",
    "    imshow(Xs, interpolation='nearest', cmap='jet_r', vmin=0, vmax=1.0)\n",
    "    colorbar()\n",
    "    #imshow(alpha, interpolation=None, cmap=get_cmap('binary'), vmin=0, vmax=1.0, alpha=0.5)\n",
    "    xticks(np.arange(dist.shape[1]))\n",
    "    yticks(np.arange(dist.shape[0]))\n",
    "    \n",
    "    if fn_name is not None:\n",
    "        title(\"Topic Alignment %s() using Jensen-Shannon Distance\" % fn_name)\n",
    "    else:\n",
    "        title(\"Topic Alignment\")\n",
    "    xlabel(\"Model 1\")\n",
    "    ylabel(\"Model 2\")\n",
    "    show()\n",
    "\n",
    "\n",
    "# In[46]:\n",
    "\n",
    "def basic_alignment(v1, v2, dist=None, dist_fn=JS_dist, debug=False):\n",
    "    \"\"\"\n",
    "    Simply aligns to the closest topic, allowing for multiple assignment. \n",
    "    Properties:\n",
    "        non-surjective, non-injective\n",
    "    \"\"\"\n",
    "    if dist is None:\n",
    "        dist = model_dist(v1, v2, dist_fn)\n",
    "    if dist.shape[0] == (v1.model.K + v2.model.K):\n",
    "        dist = filter_dist(v1, v2, dist)\n",
    "        \n",
    "    alignment = []\n",
    "    for i, topic in enumerate(dist):\n",
    "        # topic = a[i]\n",
    "        #s = topic[topic.argsort()]\n",
    "        #sim = topic.argsort()[s < 0.05]\n",
    "        closest = topic.argsort()[0]\n",
    "        alignment.append((i, closest))\n",
    "        if debug:\n",
    "            print(i, closest, topic[closest])\n",
    "    \n",
    "    return alignment\n",
    "\n",
    "\n",
    "# In[47]:\n",
    "\n",
    "def naive_alignment(v1, v2, dist=None, dist_fn=JS_dist, debug=False):\n",
    "    \"\"\"\n",
    "    First naive overlap detector just goes to next closest element if the first topic has already been assigned\n",
    "    \n",
    "    Properties: \n",
    "        k1 < k2: injective, non-surjective\n",
    "        k1 == k2: bijective\n",
    "    \"\"\"\n",
    "    if v1.model.K > v2.model.K:\n",
    "        raise ValueError(\"Models must have k1 <= k2\")\n",
    "    if dist is None:\n",
    "        dist = model_dist(v1, v2, dist_fn)\n",
    "    if dist.shape[0] == (v1.model.K + v2.model.K):\n",
    "        dist = filter_dist(v1, v2, dist)\n",
    "    \n",
    "    alignment = []\n",
    "    aligned = []\n",
    "    for i, topic in enumerate(dist):\n",
    "        # topic = a[i]\n",
    "        #s = topic[topic.argsort()]\n",
    "        #sim = topic.argsort()[s < 0.05]\n",
    "        topic_idx = 0\n",
    "        closest = topic.argsort()[topic_idx]\n",
    "        if debug:\n",
    "            print(i, closest, topic[closest])\n",
    "        \n",
    "        while closest in aligned:\n",
    "            topic_idx += 1\n",
    "            closest = topic.argsort()[topic_idx]\n",
    "            if debug:\n",
    "                print(i, closest, topic[closest])\n",
    "        \n",
    "        \n",
    "        aligned.append(closest)\n",
    "        alignment.append((i, closest))\n",
    "    \n",
    "    return alignment\n",
    "\n",
    "def compare(sample_v, v, context_type='document'):\n",
    "    sample_size = len(sample_v.labels)\n",
    "\n",
    "    print(\"{k}\\t{N}\\t{seed}\\t{LL}\\t{corpus_size}\\t\".format(k=sample_v.model.K, \n",
    "        N=sample_size, seed=sample_v.model.seed, \n",
    "        LL=sample_v.model.log_probs[-1][1],\n",
    "        corpus_size=len(sample_v.corpus)))\n",
    "\n",
    "    # compute similarity on topic-word matrix - given a topic, what is its\n",
    "    # distribution over words?\n",
    "    dist = model_dist(sample_v, v)\n",
    "    basic = basic_alignment(sample_v, v, dist=dist)\n",
    "    naive = naive_alignment(sample_v, v, dist=dist)\n",
    "    m1, m2 = zip(*basic)\n",
    "    \n",
    "    print(\"{fitness}\\t{naive_fitness}\\t{overlap}\\t\".format(\n",
    "        fitness=alignment_fitness(basic, sample_v, v, dist=dist),\n",
    "        naive_fitness=alignment_fitness(naive, sample_v, v, dist=dist),\n",
    "        overlap=len(set(m2))))\n",
    "\n",
    "    # Compute similarity on topic-document matrix - given a topic, what is its\n",
    "    # distribution over documents?\n",
    "    dist = model_doc_dist(sample_v, v, context_type)\n",
    "    basic = basic_alignment(sample_v, v, dist=dist)\n",
    "    naive = naive_alignment(sample_v, v, dist=dist)\n",
    "    m1, m2 = zip(*basic)\n",
    "    \n",
    "    print(\"{fitness}\\t{naive_fitness}\\t{overlap}\\t\".format(\n",
    "        fitness=alignment_fitness(basic, sample_v, v, dist=dist),\n",
    "        naive_fitness=alignment_fitness(naive, sample_v, v, dist=dist),\n",
    "        overlap=len(set(m2))))\n",
    "\n",
    "    # Calculate Spearman, Pearson, top-10 recall, and top-10-percent recall\n",
    "    # for each document - more of an IR-related search\n",
    "    \"\"\"\n",
    "    print(\"{spearman}\\t{pearson}\\t{recall}\\t{recall10p}\".format(\n",
    "        spearman=spearman(sample_v, v,context_type),\n",
    "        pearson=pearson(sample_v, v, context_type),\n",
    "        recall=recall(sample_v, v, context_type, N=10),\n",
    "        recall10p=recall(sample_v,v,context_type, N=int(np.floor(0.1*sample_size)))))\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Model Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running from notebook, using serial load function.\n",
      "Loading LDA data from /tmp/ap/models/ap-nltk-en-freq5-N2000-LDA-K20-document-20.npz\n",
      "Running from notebook, using serial load function.\n",
      "Loading LDA data from /tmp/ap/models/ap-nltk-en-freq5-N2000-LDA-K20-document-20.npz\n",
      "model k tokens types avg-log-likelihood perplexity\n",
      "M0 20 460795 10602 -8.471835631897047 4778.278584687178\n",
      "M1 20 460795 10602 -8.471835631897047 4778.278584687178\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jaimie/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: ClusterWarning: scipy.cluster: The symmetric non-negative hollow observation matrix looks suspiciously like an uncondensed distance matrix\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/home/jaimie/anaconda3/lib/python3.6/site-packages/matplotlib/figure.py:459: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure\n",
      "  \"matplotlib is currently using a non-GUI backend, \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\t2250\t1263529745\t-3903779.5\t460795\t\n",
      "0.0\t0.0\t20\t\n",
      "7.290601886635126e-07\t7.290601886635126e-07\t20\t\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANgAAAH+CAYAAAD+osByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHrRJREFUeJzt3X+wZGV95/H3RxAKYQgg6oQwxbCBsO7OUNcwRkEsR92VWFYUqtQkl9rKkKEmG9kYxJBVSUAZUWLKqiRmjTiyUKQyIdZGTLQI6AhjORTCMnpX4g8oLXGDYVwi6gwyRJTv/nFOm753uvv2+fX0Oac/r6pb07fvOf08XTXferrP+Z7PUURgZs14xqwnYNZnLjCzBrnAzBrkAjNrkAvMrEEuMLMGucDMGuQCM2uQC8ysQYfPegKTnHjiibF+/fpZT8PsEHv37v2XiHjOatu1usDWr1/PfffdN+tpmB1C0rem2c4fEc0aNLHAJG2W9GNJz81/f6GkkLRe0q9I+rykuyW9dcS+V0j6Z0nvHnruZkm7832W6n87Zu0yzQq2BLwuf3wBMPjM9n+AlwDnAK+V9DMr9vsIcOHwExHxaxGxGXgf8MmSczbrjGkK7A7glfnj/wh8GSAi/m9E/CSy611+Ajw9vFNEfAcYdy3MBcDHSs3YrEOmKbAfAU9KejHw1ZV/lPRq4OsRcWCaASUdDmyMiC+suvEDD0zzkmatNe1BjluBD7Fi1ZH074DfB95SYMyXA7sLbG/WWUUKbC/wvwdPSFoD3AhsjYgfFhjzAuCWAtubddZUBRYRj0fE1lieL/DfgFOB/5kfGTx1eB9JW4H3AxdK+h/5cwLOBvbUMnuzllObMzk2rVkT9x2Y6qudWVKS9kbEptW2a3UnBwcPwubN9b7m4iJs21bva5qNMV+dHEtLsHPnrGdhc6TdK9hRR8Hu3fW9Xt2rodkq5msFM0vMBWbWoCrNvi+QdL+kh8bs+0FJj0q6eOi53fnPFyV9vNZ3YtZC03wHGzT77mB5s+/XgRcDt4/Zbztw7/AYeaMvkt4C+Pi79V6VZt8Dkzo4IuKRCa/5WuDvVh35jDOmmJ5Ze1Vu9i0q/7gZEfFo1dcya7tKzb4lvY5pVi+zHijd7FvB+YAPcNhcKN3sK2mdpF3ABkm7JK0f3kfSFcDlwGWSrsyfOxY4LiKmCgwx67qJRxEjYjcrrt2KiC1Dv/6nCfteA1yz4rn9ZDEDZnOh3a1STVhacgOxJeNOjqrcQGwTzN8KtrDgBmJLxiuYWYNSB49ukLRH0l2Szqz/7Zi1S9LgUbL+xF8H3pg/Nuu11MGjJ0TEP0XEt4GVBWnWO6mDR58x5vFoDh61jksdPPr0mMdmvTTtYfpbgfMYHTy6pUDw6GOSTiYrrh8UmKdZJ01VYBHxOLAVIMsOBZYHjwJcFBHfHPwxDx59E3CCpOMj4hLgKuBmQMAlNb0Hs9aar+DRwUnhJk401/ma1noOHh1laSnr5DBLZL46ORYWssZcs0TavYLVHTxqlth8rWBmiaXORXxF3rt4Z3643qzXqvQiDnIRHx6z33ayyIBhfwi8Cngb8PZCMzXroGm+gw16EXewIhcRlp0XWyYiHhn+m6RnAQfz/e6RdG2lmbdJE1dJz5qv0q5FylzE44H9Q78ftuoeDh6dDV+lXZsirVIfAraRdWeU8T3g2KHf+9OLWPdV0rPWt9V4hpLlIkbEE8BRko6R9EvAV8q+lllXlO5FlLQOuIE8FxG4OCIeGuyT5yIuZg91UkRcTRbj9mngSeA36nsbZu2UOhdxF7Cr4BzNOssnms0a5AIza5ALzKxBLjCzBqXORfxTSZ+VdI8k3wTCei91LuLvRcTLyHIR31FqxmYdkjQXMSKeyh8eQ1agZr2WOhcRSbcAn2Ka82HORbSOS52LSERcQHaZy3um3cesq0r3Ig7lIm6dNhdR0pH5wwPAtFmKZp2VOhfxb/KDIYfjCy5tDlTpRXxv/jNu3+uB61c8d36JOZp1VrtTperORWyCsxZtAndyVOWsRZug3SuYcxGt49pdYDYzm5f+BDbPehazU1fmT+pcxI/mvYh7JDnRxlqpzsyfaVawQS/iDkbnIt4+Zr/twL0rxrgwIp6S9DLgzfgWRq21e+HSuf14XudxtSq9iAcmnWCOiEdGPDfci/ilYlM1655pVrCVvYhryw4m6Qiygj2JbDWczLmI1nGVehGLiogfRcS5wBuAq6u8llkXJMtFVOaZ+a/7gYNlX8usK5LlIgLvA26TFGTXifkAh/Ve0lxE5vrMis0jt0qZNcgFZtYgF5hZg1xgZg1ygZk1qHCB1R1GatZnZVewOsNIzXqrbIHVFkZq1mdlC6zWMNKxHDxqHVflIEdtYaRmfVW1wCqHkZr1WelMjhrDSM16q3CB1R1GatZn7U6V6kLwaB85TLU27uSwQzlMtTbtXsEcPGod5xXMrEF19yIWDiM167O6exEHYaQPj9lvO3B5yTHNOqfuXsTCYaRmfdZIL2JtHDxqHVd7L6KZ/ZtaexHNbLnSBRYRj0fE1vzaLyALI81DSDdI2iVp/fA+eRjp5cBlkq4sO7ZZV9TdiwjFw0jNessnms0a5AIza5ALzKxBLjCzBlW5CXrhDERJJ+Q3Qr8jP6Jo1mvTrGB1ZiBeBVwZEa/Ijyia9VqVm6CXyUDcALxD0p2Szq40c7MOqHwT9IIZiOcAvwg8BvwtcO6kjR/4rnMRrdsq3QS9RAbigxHx1Xx1e3rVrc06rvRN0EtmID4o6WclHU3b4wrMalD6JuiUy0C8Cvhr4CjgXTW9B7PW0lCvbuusOXVNHPhmtXh7s6IGSYGT8pYk7Y2ITau9Vqs/ph186iCbb9xcaJ/FjYtsO2tbMxMyK6hXnRxL+5bYef/OWU/D7KdavYId9cyj2L1l99TbF13tzJrWqxXMrG2q9CIWzkCUdIqkT+adHFtrfSdmLTTNR8RBL+IORmcg3j5mv+3AvSvGuIbscP6jpWZr1jFVehELZSBKeiZwCnCdpNsl/UK5KZt1R+VexAJOBM4ETgOeC7wPOH/SDmc827mI1m2VehEL+j7wlYh4NCK+DDy7wmuZdULpXsSiIuIg8LikZ0n6OWB/2dcy64rSvYiS1gE3kGcgAhdHxEODffIrlhezhzopIq4G3k12UORw4Hfqextm7TSxwOrOQIyIzwIvLThHs85qdSdHGUv7lqbq6HDPoqUwl50c7lm0VHq3gi2sXVi1f9E9i5bKXK5gZqm4wMwalDR4NH/+KEn7JI09AmnWF6mDRwG2Af9YfKpm3ZM0eFTSEcCLgD3TTM65iNZ10xTYxBueFwwevQj4y2JTNOuuZMGjkg4HzouIfyg6SbOumvY82K3AeYwOHt0yZfDo84B1km4ju2TlNXn01feKTdmsO1IHj74w/9s7gT0uLuu7Ks2+781/xu17PXD9mL+9c/opmnVXq1uligaPLu1bYmHtQnMTMiuoV50cC2sXWNy4OOtpmP1Uq1ewosGjZm3TqxXMrG2SBo/mz58k6UlJp9X2Lsxaqkov4iB49OEx+20HLh/x/KXA5wvM0ayzpvkONuhF3MGK4FFYdl5smYh4ZOXfJD0HWAM8VHbCdZk2WmDAEQNWRuVexIIuBf582o3bEjzqiAErq0ir1IfILjV5U5mBJB0HrIuIL49b9VKaJlpgwBEDVlbpXsQSzgBOz3sRNwInMyH2zawPUgePnp3/7UayEFKzXksaPDrmNcx6yyeazRrkAjNrkAvMrEEuMLMGJc1FlHSdpLsk7ZF0Zv1vx6xdUuciXhsRLyFLl7qq1IzNOiRpLuJQZsdT+T5mvZY6F3HgvcCfrbaRg0et65LlIg7tcynZzdCnSvc167KUuYhIehXZd7ZfLTZNs25KnYv4AWA/cKekByLit+p6I2ZtlDQXMSLacYGXWSKtTpUqmos4LecnWipz2cnh/ERLpdUrmHMRrevmcgUzSyVpLqKkT0j6nKTPSDq51ndi1kKpcxHfHBEvBa6lwMlps66q0ot4YNIJ5oh4ZMRzg/NkP8a9iDYHUuciIukw4ArgutW2bUsuollZlXoRS3o/cFNEfKOG1zJrtSIFtpdquYiD9qmIiJuqvI5ZV0xVYBHxeERsza/9ArJcxDwPcYOkXZLWD++T5yJeDlwm6cr86Q8CmyTtlvSuWt6BWYslzUWMiCOLTtCsy3yi2axBLjCzBrnAzBrkAjNrkAvMrEGpg0ffIOleSfdIet3Kfcz6JnXw6FuAzfnPZcWna9YtSYNHgQeAo4FjyMJvJnIuonXdNFc0r2z2XTv8x4LBox8DvkBW2BcVnKtZ56QOHt1Otgo+H7hylW3NOq90s+9Q8OjWaYNHgX8FngB+CBwx/TTNuil18OhfAHflm3y4jjdg1mapg0dvJFv1zOZCq2PbigaPLm5cZNtZ25qbkFlBvenkWNq3xM77d856GmbLtHoFKxI82kTEtllVvVnBzNpo4gomaTOwCzgpIv6fpBcC95IdPTweuAlYExHrR+z7QeANwNsj4iP5c78LLJJ1fVwWEXfX91aatbRvqbZV0t8V50fq4NEtwNnA68lOUM8df1ecL9N8Bxv0Iu5gRfAoLDsvtkxEPDLib18HjgSOA75basYzsrB2oZYbUfi74nyp3ItY0GeAr+Xjvnq1jR08al2XLHhU0rHAbwKnAy8iy6c367WUwaNPA09ExI+AH5BdtmLWa6V7ESWtA24gDx4FLo6Ihwb75MGji9lDnRQRV0v6lKS7gcOAq+t8I2ZtlDp49D3AewrO0ayzenWieXCu6sN73ahv7dCrAgOfZ7J2aXUvYlELaxdmPQWzZXq3gpm1SeECq5iVeJ2kuyTtkXRmHW/ArM3KrmBlsxKvjYiXkCVKXVVybLPOKFtgZbMSB5kdT+GboNscKFtgE2+MPkVW4nuBP1ttEAePWtdVOchRKitR0qXAVyJiT4WxzTqhymH6W4HzGJ2VuGVUVqKkV5F9P/vVCuOadUbpFWzUjdFZnpW4W9KpK3b7QP73OyVdV3Zss64ovIJVzEr0BV42V1rdyVEkF3Fp35I7Oax1etPJsbB2gcWNi7OehtkyrV7BiuQiDrjR19qkNyuYWRvV3Yv4Akn3S3pozL6fkPQ5SZ+RdHLFuZu1Xt29iKtlJb45Il5KFngz7U37zDqr7HewslmJg17EH9NgL+K4FF4n6lpqjfQiTiLpMOAKYNUTzXXmIvpKZ5uFqq1SHwK2kd3JclrvB26KiG9UGHuiUSm8TtS1Waja7FsoKzG/rWxExE0VxjXrjFp7ESWtyzMSN0jaJWn9it0+CGzK+xTfVXZss66ouxcRJmclHll0PLMu84lmswa5wMwa5AIza5ALzKxBqXMRT5D0UUl35HdfMeu1sieaB72IOxidi/g0sFvSRyLiB0P7XQVcGRFfKzmuWackzUUENgDvkHSnpLNLjm3WGWVXsIn3bZ6Qi3gO8IvAY8DfAudOGsS5iNZ1qXMRH4yIr0bEdzh0dTPrnVp7EYdyEbeOykUEHpT0s5KOpuVxBWZ1KP2ffNR9m1meiwhw0dA1YJAd5Phr4CjAvYjWe6lzEb8CbC46pllXtfpjWpFcxAHnI1qb9K6Tw/mI1iatXsHK5CKatUnvVjCzNmn1Cla3cWlT0+zn73VWRurg0VMkfTJvldpace7J+HudlVV3s+8gePT2MftdQ3Zu7NGS41YyKm3KrEl1N/seGNPBgaRnAqcA10m6XdIvlBzbrDMaafYd40TgTOA04LnA+4DzJ+1QZ/Co2SzU3uw7wffJbn7+aER8GXh2hbHNOiFZ8GhEHAQel/QsST8H7K8wtlknpA4efTfZAZD/Rdb4a9ZrqYNHPwu8tOiYZl3lTg6zBrnAzBrkAjNrkAvMrEGpg0dvlHRPfvsiN/dZ79V9E/RB8Og5wGsl/cyIfS+MiM0R4fu5Wu+lDh4N4CZJn5B0SsmxzTqjkZugTwgefWtEnAP8Edm9midy8Kh1XdLg0Yh4LP93D9M1CJt1WtLgUUnH5v+eQdb8a9ZrqYNH/0rS8WTfxX677NhmXZE6ePRXio5n1mWtDr0pEzw6joNrbBbmppPDwTU2C61ewRw8al03NyuY2SykzkXcLemz+b+vqDh3s9aruxdxkIv48IR9X5n3It5RcmyzzkiWi5h7Gtgl6WZJJ5Qc26wzUuYiArw+Ih7LL1X5A+CySRv3NRexbEZ+GYsbF9l21rYkY9mhUuYi/rQXEbgF2FBhbJvC0r4ldt7vq4Jmqcph+luB85gyFxGyXsSI2E92zdg3Kozdaaky8lOtkjZerb2IktYBN5DnIgIXR8RDQ7vdIekg8CSwpezYZl2ROhdxU9HxzLrMJ5rNGuQCM2uQC8ysQS4wswYlzUXMtz9K0j5JYw+GmPXFLHIRtwH/WHJcs05Jmoso6QjgRcCekuOadUrqXMSLgL+cdhDnIlrXJctFlHQ4cF5E/EOFMc06pdZexKFcxC0jLlt5HrBO0m3AacBrJO2NiO9VmINZqyXLRYyIbwMvzLd/J7DHxWV9lzQXcWj7dxYd16yLWp0qVTYX0RcZWlv0rpPDFxlam7R6BSuTi+iLDK1NereCmbWJC8ysQUmDR/PtT5L0pKTTKszbrBNmETx6KfD5kuOadUrS4FFJzwHWAA9NM0hfcxFtfjTS7DvBpcCflxzTrHOSBY9KOg5YFxFfrjCmWaekDB49Azg9b/bdCJzMhIg3sz4ovYJFxOMRsTW/uBLIgkfzwNENknZJWj+0/T0RcXZE/DLwaeC/Vpi3WSckDR4ds71Zb7W6VaqsSXcvcSOwpTRXnRxuBLbUermCjbt7iRuBLbW5WsHMUksaPCrpT/OboN8j6SV1vAGzNksdPPp7EfEy4I3AO0qObdYZSYNHI+Kp/OExZMVo1mupg0eRdAvwKWDXaoM4eNS6Llnw6EBEXEB2Sct7Koxt1glVC2wvo4NHt466bEXSkfnDA8DYy1rM+iJZ8Gjub/IDH4cDby87tllXJA0ejYjzi45n1mWt7uQoEzy6tG+JhbULzUzIrKDedXIsrF1gcePirKdhBrR8BSsTPGrWJr1bwczaJGkuoqSP5r2IeyQ5Msp6L3Uu4oV5L+IVwJtLjm3WGWW/gw16EXewIhcRlp0XW2ZFL+KXSo5dybirnX2lszWhbIGt7EVcO81Oko4gK86TyFa+iVIFjy7tWwJwgVntqsa2fQjYBrxpmh0i4kfAuZLOAq7m3z5mJjPqamdf6WxNqbUXcRJlnpn/uh84WGFss06otRdR0jrgBvJcRODiiHgo3+VI4DZJAQRwSflpm3VDslzEiHgS2Fx0PLMu84lmswa5wMwa5AIza5ALzKxBqXMR/3P+9zsl/fs63oBZm6XORbySrMVqEXhXybHNOiNpLmK+zQ8j4hHg50uObdYZjfQirpKL+DzgeOD5qw3iXETrutS5iL8P3Ay8DbirwthmnVDrPZqHchG3jMpFjIi7gZdLOp0s4s2s15LmIkq6gqyV6rvAb5Ud26wrUuciXgNcU3RMs65qdapUmVzESZyZaKnNVSeHMxMttVavYM5FtK6bqxXMLDUXmFmDUgeP/pe8Efg2SVMlUZl1WbLgUUmHk6VPnUsWPPrfS45t1hl1N/seGNXBkXs28HBE/ISs6/7Fqw2SKhfRrCmN3AR9jH8BTpV0NPBy4ISSY5t1RrLg0Yj4iaSr8/2+CDxYYezajYvUbmIcn+yeH8mCRwEi4u/zmz98HPhchbE7yye750vK4FEkfYDsO9u3mDJuO5VRkdpmVSULHs23+52i45l1mU80mzXIBWbWIBeYWYNcYGYNmlhgFUNGr5D0z5LePfScg0dtrkyzgpUNGf0IcOGK5xw8anNlmgIrFTIaEd8hu9EeK5538KjNjWnOg5UOGR3FwaM2T6Y9yFEmZHQUB4/aXJm2k6NwyOgoDh61eTNVgZUMGd1K1m94gqTjI+ISB4/avJlYYBVDRq8Hrl/xnINHba60OratzuDRxY2LbDtrWy2vZTatuejkWNq3xM77d856GjaHWr2C1RU8muJKZbNR5mIFM5uVKr2Iq2UgflDSo5IuHnrOuYg2V6b5iDjoRdzB6AzE28fstx24dzDGilzEBbJcxGlPUFtJS9f+CZtvnPUsumVpCRZqyiWq0os4KQORvN9wWOFcRLNZWFiAxZpyiSr3IhYwnIt4DlPkIjp4tLqFt13qMJ8ZqtSLWES+cg1yEV9Dy3IRzZpQpMAKZSCO4lxEmzelexGnyEC8guzCSkk6KSKubnMuolkTqvQiwuQMxEP6Dp2LaPOm1Z0cdZqUPe+8eGuKOzlwXrw1Z25WMGfP2yx4BTNrUOpcxEOeM+uz1LmIo54z662kuYjjshLN+mqaApt4P+aiuYhFOBfRui51LqLZXCndiziUi7h12lxEs3mTOhfxkOdqeydmLZQ6F/GQ58z6rNWdHHXlIrrX0GZlLjo53Gtos9LqFayuXESzWZmLFcxsVlLnIh7ynFmfVelFHOQiPjxmv+3A5VM8Z9ZbKXMRRz5n1meVexGb5FxE67pkuYhm8yhpLqLZvEmdi3jIc3W+GbO2SZ2L6Hs021zxiWazBrnAzBrkAjNrkAvMrEEuMLMGpQ4evVHSPZJ2S/IFWtZ7qYNHAS6MiM0RsbPEfM06JWnwaP77TZI+IemUSjM364DUwaNvjYhzgD8C3r/axg4eta5LGjwaEY/l/+4B1k4/TbNuSho8KunY/N8zgO8XmqlZByUNHgX+StLxZN/FfruuN2HWVsqOUbTTmlPXxIFv1n5PibkxyJR0Mlf9JO2NiE2rbdfq2LYiwaOLGxfZdta2ZidkVlAvOjmW9i2x836fVrP2afUKNm3waB3x2mZN6MUKZtZWqYNHd0v6bP7vK2p9J2YtNM1HxEEv4g5GB4/ePma/7cC9I8Z4ZUT8uPhUzbonafAoWb/iLkk3Szqh6GTNuiZ18OjrI2Iz8PfAH6y2sYNHreuSBo8OehGBW4ANVV7LrAuSBo8OehHJriP7RpXXMuuCpMGjwB2SDgJPAlvqextm7ZQ6eHTV3i2zPml1J0cRS/uWRnZ0uEfRZqnXnRzuUbRZ680KtrB24ZC+Rfco2qz1egUzm7XCBVYxK/HmvA/xbklLdbwBszYru4KVykqMiF/LOzneB3yy5NhmnVG2wEplJQ65AN+O1uZA2QIrnZUo6XBgY0R8YbVBnItoXVflIEfZrMSXs+LktVlfVS2wMlmJF5A1+5r1XukCi4jHI2JrLM99G85K3C3p1OF9lDUyng3sKTuuWZcUPtG8Sn/ie/OfcfsG8IKiY5p1Vas7OabNRVzat8TC2oXmJ2RWUC86ORbWLrC40ffzs/Zp9Qo2bS6iWVv1YgUza6u6exFXy0rcnf98UdLHK87drPXq7kUcZCU+PGqn/N7Mm4GbcC+izYG6exEnZiUOeS3wdyXHLmRwpfOH9344xXBmyzTSizhJ/tEyIuLR1batKxfRVzbbrNTeiziF15Fo9YLsEL7Pkdms1NqLOKXzAR/gsLlQay+ipHV5RuIGSbskrR/eJw8ePS4ivlV2XLMuqbsXESZnJe4nu+LZbC74RLNZg1xgZg1ygZk1yAVm1iAXmFmDUgePbpC0R9Jdks6s4w2YtVnS4FGyG6P/OvDG/LFZr6UOHj0hIv4pIr4NrCy+QzgX0boudfDoM8Y8Nuul1MGjT495bNZLVTI5bgXOY3Tw6JYx14U9JulksuL6QYWxzTqhdIGNujE6y4NHAS6KiG8O7XYVcDMg4JKyY5t1Rerg0S8B5xYd06yrWh3bNm3w6DgOJLVZ6/WRPAeS2qy1egWrK3jUeRw2K60usDoN0qXmiT8iz17q4NFX5H2Kd+aH661B/og8e2VXsEEv4g5GB4/ePma/PwReBfwH4O0kPFS/sHbBOfeWXLLgUUnPAg7m29xDVmRmvZYyePR4YP/Q74ettkNdwaNms5IyePR7wLFDv7sX0XovWfBoRDwBHCXpGEm/BHylwthmnVBrL6KkdcAN5MGjwMUR8dDQbtcAnwaeBH6j7NhmXZE6eHQXsKvomGZd1etWKbNZc4GZNcgFZtYgF5hZg1LnIr5B0r2S7pH0ukNf3axfUucivgXYnP9cVnJss85InYv4AHA0cAzL26bMeqnsieaVvYhrh/84IRfxY8AXyAr7otUGcfCodV3qXMTtZCve84ErK4xt1gm19iIO5SJuHXPZyr8CTwA/BI6oMLZZJ6TORfwL4K788YfLjm3WFalzEW8kW+HM5kKrQ2+q5iIOOPzFZmUuOjkc/mKzouyUVTutOXVNHPjmyiP9ZrMnaW9EbFptu7lYwcxmJXUu4u/mfYh3Szq74tzNWq/uXsRBLuLDY/bbApwNvJ7sZLRZryXLRcx9HTgSOA74bsmxzTqjkV7ECT4DfC0f99WrbexcROu6ZLmIko4FfhM4HXgRcG2Fsc06IVkuItmlK09ExI/I7s98dIWxzTqhdIFFxOMRsTWGTqRJWpfnIW6QtEvS+uHtgU9Jupvso+Ifl5+2WTe0+kTzpk2b4r777lt9Q7PEfKLZrAVavYJJehT41qznYTbCKRHxnNU2anWBmXWdPyKaNcgFZtYgF5hZg1xgZg1ygZk1yAVm1iAXmFmDXGBmDXKBmTXo/wNZ64djZvQjCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 216x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAItCAYAAACaWBZfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xm4ZFV57/Hvj54YZOxuBplVHHBCaUFjYogMojLlRm9E48WI6atPnBJNxDglqIkm3jhEM7TSgkZxwAkUg6ghxFwltkoUBAMiSANKN80kQ4/v/aP28RbHc7pP1+ld+3T199PPfk7VHtZau6pO9XvetfbaqSokSZI0fNt13QBJkqRtlYGYJElSRwzEJEmSOmIgJkmS1BEDMUmSpI4YiEmSJHXEQEwaUJJK8rCu29G2JP+Y5E1dt2NQSV6W5OdJfpFkftftmY4kf5bkQ123YzqS/EaSH3XdDmmmMBDTUCW5JMlLum7HVCTZLcnSJD9LcneS/07yuq7btSUluT7Jfc353ZHk/yZ5aZJffjdU1Uur6q1TLOuYdlu8eZLMAf4WOK6qHlRVtw25/oOagH32liivqv6yqjb792e6v3db8r2tqn+vqkdsibKkUWAgJk3u3cCDgEcBuwInAT/utEXtOLGqdgYOBN4BvA44q9smbTF7AdsDV060cUsFSJI0KAMxdSbJCUku78vEPK5v2/VJXpvk+0nuTPLJJNs32xYk+WJz3Kok/z6WwUny4CSfSbIiyU+SvLKvzD9P8qkkH2kyQFcmWbSRJj4J+HhV3V5VG6rq6qo6b9w+xyS5JsntST6QJE1dD03y9SS3JVmZ5GNJdpvi+R2VZHmS1yS5NcktSX6/79hdm3NYkeSGJG/sO/8XJflGknc1bfpJkmdO5f2oqjur6nzgd4HTkjymKfPsJG/b2Guf5KPAAcAFTRfgnzb7f7rJKN6Z5NIkj+47j7Ob1+xLzftxWZKH9m1/dJKLm3p+nuTPmvXbJTkjyY+b1/dTSfYYfz5JHg6MdYHdkeTrzfpK8odJrgGuadb9WpJvN+38dpJf6yvnkiRvaz6jv0hyQZL5zXt6V7P/QZO8rJf21f+LJE9p2v/G5r27tXkvd23qGsugLU5yc/Pev6avLX+e5J/7nv960647ktyY5EUTvA5vB34DeH/Thvdv6pzHHT/Ze3tS8zt0R/MaParvmOuTvD7JD5vP4YfHf7779t0/yWebz/NtY+2TthlV5eIytAW4BHgJ8ETgVuBIYBZwGnA9MK/Z73rgP4EHA3sAVwEvbbb9FfCPwJxm+Q0g9P6w+A7wZmAu8BDgOuAZzXF/DtwPPKup86+Ab22krR+il0n5feCQCbYX8EVgN3r/Ua0Ajm+2PQw4FpgHLKT3H/J7+o7d2PkdBawDzmzO71nAvcDuzfaPAF8AdgYOAv4bOL3Z9iJgLfAHzTm+DLgZyCTneD1wzATrfwq8rHl8NvC2jb32k5UFvLhp5zzgPcDlfdvOBlYBRwCzgY8Bn2i27QzcAryGXkZrZ+DIZturgW8B+zXl/hNw7iTnd1DzPs0e975d3LzuOzQ/bwde2LTj1Ob5/L7P7LXAQ+llRn/YvObHNPt/BPjwZtT/4qa8h9DLuH4W+Oi4/c8FdgIeS+9zdUzfZ/ifm8cHAHc37Z0DzAcO29jvXd/zjZ7zpj4nwMOBe+h9xucAf9qc09y+/a8A9m/q+g/+/2foKGB583gW8F/0ss87Ne/1r3f9PeXiMszFjJi68gfAP1XVZVW1vqrOAVYDT+7b531VdXNVrQIuAA5r1q8F9gEOrKq11RtzUvQyWAur6syqWlNV1wEfBJ7XV+Y3qurCqloPfBR4/Eba+Ap6wcHLgR8muXaC7NI7quqOqvop8K9jbayqa6vq4qpaXVUr6I1T+s1xx052fmPneGZzfhcCvwAekWQWvYzV66vq7qq6Hvg/9P5DHXNDVX2wOcdzmtdqr42c50Rupvcf6HiTvfYTqqqlTTtX0wsiHj+W/Wl8tqr+s6rW0Xutx16DE4CfVdX/qar7mzIua7b9b+ANVbW8r9znZPO6Gf+qqlZV1X3As4FrquqjVbWuqs4FrgZO7Nv/w1X146q6E/gy8OOq+mrT7k8DT9iMul8A/G1VXVdVvwBeDzxvXPv/oqruqaofAB+mFyhNVM5Xq+rc5r24raoun2IbpnLOG/O7wJeaz/ha4F30gtr+rNr7q+rG5vP99knO4Qh6f4z8SXO+91fVN6bYBmkkGIipKwcCr2m6Ne5Icge9v54f3LfPz/oe30svewDwN/T++v5KkuuSnNFX5oPHlflnPDAIGV/m9klmJ3lB0+3yiyRfBqiq+6o3OPpwetmGTwGfHtcNNmEbk+yZ5BNJbkpyF/DPwIJxr8Fk5wdwW/Of/PjtC+hl+27o23YDsO9E5VbVvc3D/rKnYl962arxJnvtf0WSWUne0XQh3kUvSwIPfB0mew32Z/LxeAcCn+t7j68C1rN5weaNfY8fzANfT/jV1/TnfY/vm+D55ry+4+u7gV5Wqr/9N47b3v97MWZjr9HmtmGsnn0n2HeTx1fVBnpt7j9+qudww7jPurRNMRBTV24E3l5Vu/UtOzZ/mW9Ukx15TVU9hN5f8H+c5OimzJ+MK3PnqnrWFMr8WPWuqntQVf3KmKqqugv4S3rdJwdP4fz+il4X0+Oqahfg9+h1n07XSnpZqQP71h0A3LQFygYgyZPo/Yf6K5mJjbz20Dvffs8HTqbXhbcrvW43mNrrcCO9rsDJtj1z3Pu8fVVtzmvQ39abeeDrCVvuNZ0oWzi+vgPodUX3B3f7j9t+8wTlbOw12lQ7NvecN3p8ktBrc//xUz2HAzYzmymNFAMxdeWDwEuTHJmenZI8O8nOmzowvUH+D2u+/O+ilw1ZT2/M1V1JXpdkhyYj85gmsNhsSd6U5ElJ5jYDjV8F3MH/HwC+MTvT6068I8m+wJ8M0obxmu7GTwFvT7JzkgOBP6aXcZuWJLskOQH4BL1xSD+YYJ/JXnvoBRIP6dt9Z3rdzbcBO9ILZKfqi8DeSV6dZF5zrkc22/6R3vkf2LRpYZKTN6Ps8S4EHp7k+U129HeBQ5s2TNcKYAMPfF3OBf4oycFJHkTvdfnkuKzQm5LsmN7FDb8PfHKCsj9G72KR/9m0e36SwybYD371vdnccx5//KeAZyc5Or0pQl5D773+v337/GGS/ZoM8p9Ncg7/SW8s4Dua74Dtkzx1kjZII8lATF2oqlpGb5zY++kNEr6W3kDzqTgE+Cq9QOebwN9X1SVNkHIivXFGP6GXPfoQvWzMQO2kNz5nJb2/5o8Fnt2M69mUv6B3QcKdwJfoDcjeUl5Bb6D0dfSyVh8Hlk6jvAuS3E0vO/EGeuPZfn+SfSd87ZttfwW8sekyfC29Qew30MuS/JDeAPspqaq76b3eJ9LrvrwG+K1m83uB8+l1j97dlHvkROVMsa7b6I1Jew29oPFPgROqauWgZfaVfS+98VH/0bwuT6b3Xn2U3gUcP6F3Ackrxh36b/R+J74GvKuqvjJB2T+ldyHHa+h1I1/O5GMe30tvHN3tSd43wDk/4L2tqh/Ry/L+Hb3fjxPpTYOypu+YjwNfofc5vQ542wTnMPY7+zB6F4gspzf+TNpmjF3tJA1Fku/SG4T++a7bIs006U2D8RNgztY8birJ9fSu0vxq122RZjozYhqappvlUcD3um6LJEkzgYGYhiLJO+l1U7yuqsZfrSVJ0jbJrklJkqSOmBGTJEnqyFYxd8sOu+xeO+851XkGB7Pb9u2/FLOzJaaR2oSWq9hhzqx2K5A0qbXr2+/BmLVd+99TG4bQEzN7COehTbvhhutZuXLlZr8ZS8/+aB1w4IEc81tPG/k3cqsIxHbec1/+5998utU6TnnkwlbLB9ht+zmt17Fdy8Heofvt0mr5kib3szvub72O3XZs/3vqntXtXxA6f+d5rdehTXvqkYsGOu5v3/XXPPZxj+eY33raFm7RzGPXpCRJmjGSPOKggw/myit+QJK5XbenbVtFRkySJG0b3vqX77h6rz334vLLv8df/fW7VtP6oJtumRGTJEkzxoVfvIBnPvsETjr5FM7//Oe6bk7rDMQkSdKMkOTB8+bNY8GCBTz113+D//sf3yDJSMcqdk1KkqQZ4b1/9/c3rV27FoDZs2dz2BOeyBVX/OBJwGXdtqw9BmKSJGlGuOD8z/P+f1jyy+cnnfLb7LX33t9ihMeJjXS6T5IkbR2S7HbXXXdx4IEH/nLdMccex9e/enGHrWqfGTFJktS5pef88+03XP+TB6zbYYcd2P/AA0nyqKq6qqOmtaqTQCzJHwEvAQr4AfD7VdX+TIWSJGlGOv8Ln+ONb/rzX1l/8sm/zeGHL/ohI9o9OfSuyST7Aq8EFlXVY4BZwPOG3Q5JkjQzJNnhJz/+MYc++tG/su2Zzz6BC7/0xQ5aNRxdjRGbDeyQZDawI3BzR+2QJEkdO+9z59/79GOOJRPcpm/33Xdnp512Isl+HTStdUMPxKrqJuBdwE+BW4A7q+or4/dLsjjJsiTL7rtr1bCbKUmShuT8z3+OE086ZdLtJ550Cn/7nr+7cYhNGpouuiZ3B04GDgYeDOyU5PfG71dVS6pqUVUt2mGXPYbdTEmSNARJZn/3O8s48slPnnSfE046mS9e8IUhtmp4uuiaPAb4SVWtqKq1wGeBX+ugHZIkqXtPffJTfo3ttps8JNlvv/249957STJymZkuArGfAk9OsmN6ncFHAyN5SaokSdq4P3zFqy458eTJuyXHPOvZJ/ChpefcNoQmDVUXY8QuA84Dvktv6ortgCUbPUiSJI2cJLn0kn/lqN96+ib3PfGkUzj//M8PoVXD1ck8YlX1FuAtXdQtSZJmjMMe9ehHM3fu3E3u+MhHPYobb7iBJDtW1b1DaNtQeIsjSZLUide/4U3f3djVkuMdfexxfPK8z93TYpOGzkBMkiR14isX/QvPOP6ZU97/xJNO4YIvjFb3pIGYJEkauiQP3XPPPdl5552nfMyiJz2Jy7/3XZoJ4TdV/vFJfpTk2iRnTLD9aUm+m2RdkueM23Zakmua5bS+9Ycn+UFT5vsy0Qy0m8lATJIkDd1fvvNvrt2cbkmA7bbbjl976q8DPG1j+yWZBXwAeCZwKHBqkkPH7fZT4EXAx8cduwe9cexHAkcAb2nmQAX4B2AxcEizHL9ZJzABAzFJkjR0F37xAp594kmbfdxJp/w2wKYiuCOAa6vquqpaA3yC3mTyv1RV11fV94EN4459BnBxVa2qqtuBi4Hjk+wD7FJV36yqAj4yhXZskoGYJEkaul/cfTd77rnnZh+36ElHADxj7DaIzbJ43G77Av23RFrerJuKyY7dt3k8SJmT6mT6is212/azOeWRC1ut4/NXr2i1fIAXP2Ha79cmrVk3PrDXRH5+5/2t1zH/QZu+HHu6br9nbet1LNxlXqvl9/6wbNfqIfxebD9nVut1rF3f/nncdf+69uu4t/3P7fyd2/3cavo2NpP+FI67paqO2shuE43dmuqXzWTHTqfMSZkRkyRJQ1dA1WDLFCwH9u97vh9w8xSbNtmxy5vHg5Q5KQMxSZI0ar4NHJLk4CRzgecB50/x2IuA45Ls3gzSPw64qKpuAe5O8uTmasn/BUz7TuQGYpIkqRM14L9Nllu1Dng5vaDqKuBTVXVlkjOTnASQ5ElJlgPPBf4pyZXNsauAt9IL5r4NnNmsA3gZ8CHgWuDHwJen+xpsFWPEJEmSNkdVXQhcOG7dm/sef5sHdjX277cUWDrB+mXAY7ZkOw3EJElSNwYZ6t7+9T1DZdekJElSR8yISZKkTpgQMyMmSZLUmU4CsU3diFOSJI2+FucR22oMPRCb4o04JUmSRl4XGbFN3ohTkiSNvrbmEduadBGITedGnJIkSSOji0BsSjfNTLJ47K7qd666bQjNkiRJQ1UDLiOki0BsSjfirKolVbWoqhbtusf8oTVOkiS1b9AYbMTisE4CsenciFOSJGlkDH1C16pal2TsRpyzgKVVdeWw2yFJkjo04FQUozZ9RScz6090I05JkqRtjbc4kiRJnRhkKgqnr5AkSdIWYUZMkiR1wjFiZsQkSZI6YyAmSZLUEQMxSZKkjjhGTJIkDV3hGDHYSgKx2Qm7bT+n1Tpe/IT27zu+9Hs3tV7HcQ/bvfU6RsGuO7b7eQKYPav9hPMwzqNtG4bxpTqEOu5fu771Otatb/9Edpw30e2At6x5c+yMEfQmonD6Cn8bJEmSOrJVZMQkSdLoGaibcbQSYmbEJEmSumJGTJIkdcKEmBkxSZKkzpgRkyRJ3TAlZkZMkiSpK2bEJElSJ5xHrIOMWJJHJLm8b7kryauH3Q5JkqSuDT0jVlU/Ag4DSDILuAn43LDbIUmSulPlLY6g+zFiRwM/rqobOm6HJEnS0HU9Rux5wLkdt0GSJHXAiyY7zIglmQucBHx6ku2LkyxLsuz2VbcNt3GSJKl9NeAyQrrsmnwm8N2q+vlEG6tqSVUtqqpFu+8xf8hNkyRJal+XXZOnYrekJEnbLKev6CgjlmRH4Fjgs13UL0mSNBN0khGrqnsB+xslSdqGOX1F99NXSJIkbbO6nr5CkiRto5y+woyYJElSZ8yISZKkoSscIwZmxCRJkjpjRkySJHXEUWIGYpIkafjKrknYWgKxwHZJq1WsWbeh1fIBjnvY7q3X8ZVrb2+1/OMetXer5QOs39D+b9nt96xtvY6FO7f7mQW4Z/W61uuYO3tuq+XXEL5VZ23X/nuxfgjnMXd2+6NJVq9d33oda9eP2P+k0jRsHYGYJEkaOXZMOlhfkiSpM2bEJElSJxwjZkZMkiSpMwZikiSpEzXgv6lIcnySHyW5NskZE2yfl+STzfbLkhzUrH9Bksv7lg1JDmu2XdKUObZtz+m+BgZikiRppCSZBXwAeCZwKHBqkkPH7XY6cHtVPQx4N/BOgKr6WFUdVlWHAS8Erq+qy/uOe8HY9qq6dbptNRCTJEndqAGWqTkCuLaqrquqNcAngJPH7XMycE7z+Dzg6ORX5so6FTh36ie0+QzEJEnS1mafJMv6lsXjtu8L3Nj3fHmzbsJ9qmodcCcwf9w+v8uvBmIfbrol3zRB4LbZWgvEkixNcmuSK/rW/XmSm/r6Vp/VVv2SJGnmGiQZ1pcUu6WqFvUtS8YVP1GAND6fttF9khwJ3FtVV/Rtf0FVPRb4jWZ54SZPdBPazIidDRw/wfp39/WtXthi/ZIkadu0HNi/7/l+wM2T7ZNkNrArsKpv+/MYlw2rqpuan3cDH6fXBTotrQViVXUpDzwhSZKkX6oabJmCbwOHJDk4yVx6QdX54/Y5Hzitefwc4OvV3HMtyXbAc+mNLaNZNzvJgubxHOAE4AqmqYsxYi9P8v2m63LSmy8mWTzW93v7bbcNs32SJKl1g05eselIrBnz9XLgIuAq4FNVdWWSM5Oc1Ox2FjA/ybXAHwP9U1w8DVheVdf1rZsHXJTk+8DlwE3AB6f7Kgx7Zv1/AN5Krw/2rcD/AV480Y5Nf+8SgEMf94QRm0dXkiS1qRn+dOG4dW/ue3w/vazXRMdeAjx53Lp7gMO3dDuHmhGrqp9X1fqq2kAvipx236okSdpKTWO0/qgYaiCWZJ++p7/NFuhblSRJ2lq11jWZ5FzgKGBBkuXAW4CjmtsEFHA98L/bql+SJM1sgyS3Riwh1l4gVlWnTrD6rLbqkyRJ2toMe7C+JEkSTH0qigceNmIpMW9xJEmS1BEzYpIkaeh6F0BufnprkGNmMjNikiRJHTEjJkmSuuFlkwZikiSpG8ZhW0kgtsOcWRy63y5dN2OrcNyj9m61/NdecFWr5QMsXrRf63Xsvdv2rddx9c13t17HTvPa/xWeN7vdEQxr17f/tbrrjnNar2P9mvWt17HXrvNar+Oan/2i9Tp222lu63WMgvuH8JlacffqVstfs25Dq+WPgq0iEJMkSaNnoKkoRiwl5mB9SZKkjpgRkyRJnXD6CjNikiRJnTEjJkmSuuEYMTNikiRJXTEjJkmSOmFCzIyYJElSZ8yISZKkoSsGm0dsoLnHZrDWMmJJlia5NckVfeuem+TKJBuSLGqrbkmSpK1Bm12TZwPHj1t3BfA/gEtbrFeSJG0FasB/o6S1rsmqujTJQePWXQWQpK1qJUnS1qBwtD4zeLB+ksVJliVZtmLliq6bI0mStMXN2ECsqpZU1aKqWrRwwcKumyNJkrawGnAZJTM2EJMkSRp1Tl8hSZI64fQV7U5fcS7wTeARSZYnOT3JbydZDjwF+FKSi9qqX5IkaaZr86rJUyfZ9Lm26pQkSVuPQaaiGLXpKxwjJkmS1BHHiEmSpG44j5gZMUmSpK6YEZMkSUPnxPo9BmKSJKkTTl9h16QkSVJnzIiNmPUb2v1TYfGi/VotH2DJsuWt1/G633xI63XMmdX+3zkbhvCn4c/uXN16HW1beffWfw4AC3eZ13oda9ZtaL2Ou+9b23od++y2fet1tG3u7Pa/Q3ac124YsN122cjWcvoKzIhJkiR1xoyYJEnqhqP1zYhJkiR1xYyYJEnqhAkxM2KSJEmdMSMmSZKGr5xHDMyISZIkdcaMmCRJGrreLY6cR8yMmCRJUkdaC8SSLE1ya5Ir+tY9Psk3k/wgyQVJdmmrfkmSNMPVgMsIaTMjdjZw/Lh1HwLOqKrHAp8D/qTF+iVJ0gxmHNZiIFZVlwKrxq1+BHBp8/hi4Hfaql+SJGmmG/YYsSuAk5rHzwX2n2zHJIuTLEuybMXKFUNpnCRJGp6qwZZRMuxA7MXAHyb5DrAzsGayHatqSVUtqqpFCxcsHFoDJUmShmWogVhVXV1Vx1XV4cC5wI+HWb8kSZo5asB/U5Hk+CQ/SnJtkjMm2D4vySeb7ZclOahZf1CS+5Jc3iz/2HfM4c0Fh9cmeV+STPc1GGoglmTP5ud2wBuBf9z4EZIkSZsnySzgA8AzgUOBU5McOm6304Hbq+phwLuBd/Zt+3FVHdYsL+1b/w/AYuCQZhl/UeJma3P6inOBbwKPSLI8yen0Xoj/Bq4GbgY+3Fb9kiRphmvvsskjgGur6rqqWgN8Ajh53D4nA+c0j88Djt5YhivJPsAuVfXNqirgI8ApU2rNRrQ2s35VnTrJpve2VackSdom7JNkWd/zJVW1pO/5vsCNfc+XA0eOK+OX+1TVuiR3AvObbQcn+R5wF/DGqvr3Zv/l48rcd7on4i2OJElSJwa5ALI55paqOmoju02U2Rpf3WT73AIcUFW3JTkc+HySR0+xzM3mLY4kSdKoWc4Dp8jaj96QqAn3STIb2BVYVVWrq+o2gKr6Dr0LCx/e7L/fJsrcbAZikiRp6ArYULXZS01tIrFvA4ckOTjJXOB5wPnj9jkfOK15/Bzg61VVSRY2g/1J8hB6g/Kvq6pbgLuTPLkZS/a/gC9M93Wwa1KSJHViGl2TG9+nN+br5cBFwCxgaVVdmeRMYFlVnQ+cBXw0ybX07gT0vObwpwFnJlkHrAdeWlVjdwp6Gb1bOO4AfLlZpsVATJIkjZyquhC4cNy6N/c9vp/eXX7GH/cZ4DOTlLkMeMyWbKeBmCRJGr4Bb1c0arc4MhBr/PzO+1uvY9cd57Rex+33rG21/L13277V8gFe95sPab2Od/7bda3XcfoTp31V8ybdu3pD63U89oBdW69jFExx3Mq0bIFJvDdp3z12aL2Oe1evb72OUbDddu2/3/MfNLfV8mcP4Ry2dgZikiSpA1O/XdH4o0aJV01KkiR1xIyYJEkaut70FZt/3CDHzGRmxCRJkjpiRkySJHXCMWJmxCRJkjpjRkySJHVioHnEtnwzOmVGTJIkqSNmxCRJUicGGyM2WjrJiCV5VZIrklyZ5NVdtEGSJHVnbPqKQZZRMvRALMljgD8AjgAeD5yQ5JBht0OSJKlrXWTEHgV8q6rurap1wL8Bv91BOyRJUodqwH+jpItA7ArgaUnmJ9kReBaw//idkixOsizJshUrVwy9kZIkSW0beiBWVVcB7wQuBv4F+C9g3QT7LamqRVW1aOGChUNupSRJalvVYMso6WSwflWdVVVPrKqnAauAa7pohyRJUpc6mb4iyZ5VdWuSA4D/ATyli3ZIkqSODJjdGrWMWFfziH0myXxgLfCHVXV7R+2QJEnqTCeBWFX9Rhf1SpKkmaGADQNcATnIMTOZtziSJEnqiLc4kiRJnXCMmIGYJEnqxGCTszqhqyRJkrYIM2KSJKkTdk2aEZMkSeqMGbHG/AfNbb2O2bPaj3sX7pxWy7/65rtbLR9gzhBep9OfuG/rdZz13Ztar+OUR7Z/+6816za0Wn4N4c/bpN3fC4D71qxvvY6d5s1qvY677/uVO85tcWvXt/uZGhV33LOm9Tra/v1eu37y3+9Bp69wjJgkSZK2CDNikiSpE44RMyMmSZLUGTNikiSpE4Mkt0YsIWZGTJIkqStmxCRJUicGuWp6GFdaD5MZMUmSpI6YEZMkSUPXm0ds843aLHRDD8SSbA9cCsxr6j+vqt4y7HZIkqQOlV2T0E1GbDXw9Kr6RZI5wDeSfLmqvtVBWyRJkjoz9ECseqHsL5qnc5pltMJbSZK0SU5f0dFg/SSzklwO3ApcXFWXddEOSZKkLnUSiFXV+qo6DNgPOCLJY8bvk2RxkmVJlq1YuWL4jZQkSa2qqoGWUdLp9BVVdQdwCXD8BNuWVNWiqlq0cMHCobdNkiSpbUMPxJIsTLJb83gH4Bjg6mG3Q5IkdWds+opBllHSxVWT+wDnJJlFLxD8VFV9sYN2SJIkdaqLqya/Dzxh2PVKkqSZZLDxXo4RkyRJ0hbhLY4kSVInBklujVhCzEBMkiR1wwld7ZqUJEnqjBkxSZI0dAVsGKCfcZBjZjIzYpIkSR0xEJMkSZ2oAZepSHJ8kh8luTbJGRNsn5fkk832y5Ic1Kw/Nsl3kvyg+fn0vmMuacq8vFn2HPTcx9g12bj9nrWt17HrjnNar+Oe1etaLX+nee1/ZIaRdr53dftzM5/yyPZvzfX5q9u/D+uTH7JHq+Xfv67993v9hvbf72F8bu9Zvb71OtZtaP887lvT/nkQBJd9AAAgAElEQVSMgt12mtt6Hfe3/F7M2i6tlj+ZZtL4DwDHAsuBbyc5v6p+2Lfb6cDtVfWwJM8D3gn8LrASOLGqbm7uhX0RsG/fcS+oqmVbqq1mxCRJ0vBVqzf9PgK4tqquq6o1wCeAk8ftczJwTvP4PODoJKmq71XVzc36K4Htk8zbAmc8IQMxSZK0tdknybK+ZfG47fsCN/Y9X84Ds1oP2Keq1gF3AvPH7fM7wPeqanXfug833ZJvSjLtlJ9dk5IkaejGbvq9uZpjbqmqozay20QB0vhU2kb3SfJoet2Vx/Vtf0FV3ZRkZ+AzwAuBj2yy0RthRkySJI2a5cD+fc/3A26ebJ8ks4FdgVXN8/2AzwH/q6p+PHZAVd3U/Lwb+Di9LtBpMRCTJEmdqNr8ZYqXTX4bOCTJwUnmAs8Dzh+3z/nAac3j5wBfr6pKshvwJeD1VfUfYzsnmZ1kQfN4DnACcMU0Th8wEJMkSSOmGfP1cnpXPF4FfKqqrkxyZpKTmt3OAuYnuRb4Y2BsiouXAw8D3jRumop5wEVJvg9cDtwEfHC6bXWMmCRJ6kSbM+tX1YXAhePWvbnv8f3Acyc47m3A2yYp9vApN3SKDMQkSdLQFU1X4+YeN1p3OOquazLJrCTfS/LFrtogSZLUpS4zYq+i12+7S4dtkCRJHfGm3x1lxJrLQp8NfKiL+iVJkmaCrjJi7wH+FNh5sh2aWXIXA+x/wAFDapYkSRqKgkFubTqE26EO1dAzYklOAG6tqu9sbL+qWlJVi6pq0cIF7d88WZIkadi6yIg9FTgpybOA7YFdkvxzVf1eB22RJEkdGPiqyS3ekm4NPSNWVa+vqv2q6iB6M91+3SBMkiRti5xHTJIkdaDYMEB+a5BjZrJOA7GqugS4pMs2SJIkdcWMmCRJ6oQz6xuISZKkDhROXwEd3uJIkiRpW2dGTJIkdcJbHJkRkyRJ6owZMUmSNHzlYH0wIyZJktQZM2KNhbvM67oJW8Tc2XNbLX/e7PZj95/dubr1Oh57wK6t17Fm3YbW63jyQ/ZovY4zLvxRq+Wf/sR9Wy0fhvMXdNJ+HfN3bv97au4Qfsd3mjca37eaHq+a7DEjJkmS1BEzYpIkqRPlIDEzYpIkSV0xIyZJkjrhGDEzYpIkSZ0xIyZJkobOqyZ7BsqIJXnzlm6IJEnattSA/0bJoF2TL9mirZAkSdoGTdo1meSuyTYBO0yn0iTXA3cD64F1VbVoOuVJkqStTNk1CRsfI3YH8KSq+vn4DUlu3AJ1/1ZVrdwC5UiSJG2VNhaIfQQ4EPiVQAz4eDvNkSRJ2wrnc93IGLGqemNV/eck2143zXoL+EqS7yRZPNEOSRYnWZZk2YqVK6ZZnSRJ0szT1fQVT62qm5PsCVyc5OqqurR/h6paAiwBOPzwRSMW/0qStG3rTV+x+f+9D3LMTNbJhK5VdXPz81bgc8ARXbRDkiSpS0MPxJLslGTnscfAccAVw26HJEnq1oYabBklG5u+Yo+NHVhVqwascy/gc0nG6v94Vf3LgGVJkiRttTY2Ruw79LpwM8G2Ah4ySIVVdR3w+EGOlSRJo8OrJjcSiFXVwcNsiCRJ2nZUlYP1mcIYsfT8XpI3Nc8PSOLgekmSpGmaymD9vweeAjy/eX438IHWWiRJkrYJVYMto2Qq84gdWVVPTPI9gKq6PcncltslSZI08qYSiK1NMoveAH2SLAQ2tNoqSZI08gYJJkYtAJlK1+T76E26umeStwPfAP6y1VZJkiRtAzaZEauqjyX5DnA0vaksTqmqq1pv2ZDVEDqdhzEJXdvnsXb9iHXOt2gYn6n717Vfx+lP3LfV8s/67k2tlg/w/Mfs3Xods7abaKafLeue+9e1Xsd9a9a3XsdO23d1dz2Nt/3cWa2Wv7Ffi0FvcTSM79ZhmuqErrcC5/Zvm8aErpIkSWLqE7oeANzePN4N+CngPGOSJGlgTui6kTFiVXVwVT0EuAg4saoWVNV84ATgs8NqoCRJ0qiaymD9J1XVhWNPqurLwG+21yRJkrQt8KbfU5u+YmWSNwL/TK+r8veA21ptlSRJ0jZgKhmxU4GF9Kaw+DywZ7NOkiRpIL1Z8muApeuWb1lTmb5iFfCqJLsAG6rqF+03S5IkjbpBuhm3uQldkzy2ub3RD4Ark3wnyWPab5okSdJom0rX5D8Bf1xVB1bVgcBrgCXtNkuSJI26NgfrJzk+yY+SXJvkjAm2z0vyyWb7ZUkO6tv2+mb9j5I8Y6plDmIqgdhOVfWvY0+q6hJgp+lUmmS3JOcluTrJVUmeMp3yJEmSxjT3yP4A8EzgUODUJIeO2+104PaqehjwbuCdzbGHAs8DHg0cD/x9kllTLHOzTSUQuy7Jm5Ic1CxvBH4yzXrfC/xLVT0SeDwwcrdMkiRJkysGHaw/pZTYEcC1VXVdVa0BPgGcPG6fk4FzmsfnAUcnSbP+E1W1uqp+AlzblDeVMjfbVAKxF9O7avKz9K6cXAj8/qAVNoP+nwacBVBVa6rqjkHLkyRJ25x9kizrWxaP274vcGPf8+XNugn3qap1wJ3A/I0cO5UyN9tUrpq8HXjldCvq8xBgBfDhJI+ndyulV1XVPf07NS/qYoD9DzhgC1YvSZJmgkGugGyOuaWqjtrIbhPdbnx8Km2yfSZbP1HyatqTaWzspt/nb+zAqjppGnU+EXhFVV2W5L3AGcCbxpW/hOaigMMPXzRis4ZIkqQWLQf273u+H3DzJPssTzIb2BVYtYljN1XmZttYRuwp9FJw5wKXMXGEOIjlwPKquqx5fh69QEySJG1Dpjjea5Bjvg0ckuRg4CZ6g++fP26f84HTgG8CzwG+XlXVJKI+nuRvgQcDhwD/SS8O2lSZm21jgdjewLH0ZtF/PvAl4NyqunI6FVbVz5LcmOQRVfUj4Gjgh9MpU5IkaUxVrUvycuAiYBawtKquTHImsKyqzqc3Vv2jSa6llwl7XnPslUk+RS82WQf8YVWtB5iozOm2ddJArKn0X4B/STKPXkB2SZIzq+rvplnvK4CPJZkLXMc0Bv9LkqStUDHQ7YqmekxVXQhcOG7dm/se3w88d5Jj3w68fSplTtdGB+s3Adiz6QVhBwHvo3f15LRU1eXAoumWI0mStk5FsWGASGyQY2ayjQ3WPwd4DPBl4C+q6oqhtUqSJGkbsLGM2AuBe4CHA6/szXEG9AarVVXt0nLbJEnSCGuza3JrsbExYlOZ7FWSJEkD2uSErpIkSW1ocfqKrYZZL0mSpI6YEZMkSZ1wjJiB2C+tXjfIHa820xA+PLO221I3QJjYrjvOabV8gJV3r269jmHou8ClNes3tP+5bftL7/mP2bvdCoCPX/Gz1us4+ZELWq9jp3ntf2XPmd1+R8m69UP4vpW2EgZikiRp6KoGmxPMMWKSJEnaIsyISZKkTgyS2xqtfJgZMUmSpM6YEZMkSZ1wHjEDMUmS1IECNgwQUw1yzExm16QkSVJHzIhJkqRO2DVpRkySJKkzQw/Ekuyf5F+TXJXkyiSvGnYbJElS96oGW0ZJF12T64DXVNV3k+wMfCfJxVX1ww7aIkmS1JmhB2JVdQtwS/P47iRXAfsCBmKSJG0rqhwjRsdjxJIcBDwBuGyCbYuTLEuybMXKFcNumiRJUus6C8SSPAj4DPDqqrpr/PaqWlJVi6pq0cIFC4ffQEmS1JqxecQGWUZJJ4FYkjn0grCPVdVnu2iDJElS14Y+RixJgLOAq6rqb4ddvyRJmhkcI9ZNRuypwAuBpye5vFme1UE7JEmSOtXFVZPfADLseiVJ0swySG5rtPJh3uJIkiR1ZMMA3YyDHDOTeYsjSZKkjpgRkyRJQzfo7YpGLCFmRkySJKkrZsQkSVInnL7CjJgkSVJnzIg1tp8zq/U67l+7vvU61rf8l8L6Ne2fwzAM4y+q+4bwWg3j6qG0PNnMrO3an83m5EcuaL2OL1y9svU6Xv3Ug1qvYxif27T9odJWY6AxYlu+GZ0yIyZJktQRM2KSJGnoejf9dh4xM2KSJEkdMSMmSZI6MVBya7QSYgZikiSpA1VOX4Fdk5IkSZ0xIyZJkjqxYYDk1iDHzGRmxCRJkjrSWiCWZGmSW5Nc0bfusCTfSnJ5kmVJjmirfkmSNHMVUAP+GyVtZsTOBo4ft+6vgb+oqsOANzfPJUmStkmtjRGrqkuTHDR+NbBL83hX4Oa26pckSTPbQLc4Gq2E2NAH678auCjJu+hl435tyPVLkiTNGMMerP8y4I+qan/gj4CzJtsxyeJmHNmyFStXDK2BkiRpOKqZS2xzl1Ey7EDsNOCzzeNPA5MO1q+qJVW1qKoWLVywcCiNkyRJGqZhB2I3A7/ZPH46cM2Q65ckSTPEhhpsGSWtjRFLci5wFLAgyXLgLcAfAO9NMhu4H1jcVv2SJEkzXZtXTZ46yabD26pTkiRtHaoGu2+kY8QkSZK2gF4wtvnLdCTZI8nFSa5pfu4+yX6nNftck+S0Zt2OSb6U5OokVyZ5R9/+L0qyopm0/vIkL5lKewzEJEnStuQM4GtVdQjwteb5AyTZg96QqiPpXVj4lr6A7V1V9UjgCcBTkzyz79BPVtVhzfKhqTTGQEySJHViQ9VAyzSdDJzTPD4HOGWCfZ4BXFxVq6rqduBi4Piqureq/hWgqtYA3wX2m05jDMQkSdK2ZK+qugWg+bnnBPvsC9zY93x5s+6XkuwGnEgvqzbmd5J8P8l5SfafSmOGPbO+JEkSMK1bHO2TZFnf6iVVtWTsSZKvAntPcPgbplhNJqq6r/zZwLnA+6rqumb1BcC5VbU6yUvpZduevqmKDMQkSdLW5paqOmqyjVV1zGTbkvw8yT5VdUuSfYBbJ9htOb0puMbsB1zS93wJcE1Vvaevztv6tn8QeOfGTmCMXZOSJGnoisFub7QFpq84n96dfmh+fmGCfS4CjkuyezNI/7hmHUneBuxK7/7Zv9QEdWNOAq6aSmO2iozY2vXFz+64v+U6NrRaPsC69e3PfTJ3drux9V67zmu1fICFu7RfRzJR1nnL2mnerNbruGf1+tbrmL9zu+/HPfeva7V8gJ3mtf9V9+qnHtR6He/5j+tbr+NFhz249TpGbWZ0bXXeAXwqyenAT4HnAiRZBLy0ql5SVauSvBX4dnPMmc26/eh1b14NfLf5v+T9zRWSr0xyErAOWAW8aCqN2SoCMUmSNGIGnBNsugmxpgvx6AnWLwNe0vd8KbB03D7LmXj8GFX1euD1m9seuyYlSZI6YkZMkiR1YrBbHLXQkA4ZiEmSpE500TU509g1KUmS1BEzYpIkqRODdU2OVkrMjJgkSVJHzIhJkqShK8yIQYsZsSRLk9ya5Iq+dX+T5Ormhpifa26YKUmStE1qs2vybOD4cesuBh5TVY8D/psBJj6TJEkjoJnQdXMXRish1l4gVlWX0pviv3/dV6pq7H4m36J3E01JkqRtUpdjxF4MfHKyjUkWA4sBHrzf/sNqkyRJGpKBxoiNWEqsk6smk7yB3k0xPzbZPlW1pKoWVdWi+fMXDq9xkiRJQzL0jFiS04ATgKNr1C59kCRJU+bM+kMOxJIcD7wO+M2quneYdUuSJM00rQViSc4FjgIWJFkOvIXeVZLzgIuTAHyrql7aVhskSdJMVc4jRouBWFWdOsHqs9qqT5IkbT16E7oOcNxoxWHe4kiSJKkr3uJIkiQNX3mLIzAjJkmS1BkzYpIkqROOETMjJkmS1BkzYpIkqROOETMjJkmS1JmtIiM2a7uw245zWq3jrvvXtVo+wI7z0nodq9eub7X8a372i1bLB1izbkPrdey7xw6t13H3fe1/ptZtaP8vw7mz2/177b417X5mAea0fA4wnPN40WEPbr2Osy+/ufU6nnvoXq3Xoa2DY8TMiEmSJHVmq8iISZKk0dKbWd8xYgZikiRp+MquSbBrUpIkqTNmxCRJUifsmjQjJkmS1BkzYpIkqROOETMjJkmS1JnWArEkS5PcmuSKCba9NkklWdBW/ZIkaSYrqgZbRkmbGbGzgePHr0yyP3As8NMW65YkSZrxWgvEqupSYNUEm94N/Cm9udwkSdI2aGxCVzNiQ5TkJOCmqvqvKey7OMmyJMtWrlwxhNZJkiQN19ACsSQ7Am8A3jyV/atqSVUtqqpFCxYsbLdxkiRpuJqZ9QdZRskwM2IPBQ4G/ivJ9cB+wHeT7D3ENkiSJM0YQ5tHrKp+AOw59rwJxhZV1cphtUGSJM0czqzf7vQV5wLfBB6RZHmS09uqS5IkbX3smmwxI1ZVp25i+0Ft1S1JkrQ18BZHkiRp6ArYsGGArskBjpnJvMWRJElSR8yISZKk4RtwvNeojREzIyZJktQRM2KSJKkTA01fMWJ3SDQjJkmS1BEzYpIkqROOEdtKArENVdyzel2rddx179pWyweYN6f9BOTa9e1+QnfbaW6r5QPcfV/778W9q9e3Xsfa9Rtar+O+Ne2fx07z5rVb/vbtfw2tG8J7kaT1OoZx1f5zD92r9To+/cOft17HkQ/do/U6pC1hqwjEJEnSqClvcYRjxCRJUgeKbm5xlGSPJBcnuab5ufsk+53W7HNNktP61l+S5EdJLm+WPZv185J8Msm1SS5LctBU2mMgJkmStiVnAF+rqkOArzXPHyDJHsBbgCOBI4C3jAvYXlBVhzXLrc2604Hbq+phwLuBd06lMQZikiRp+KrXzTjIMk0nA+c0j88BTplgn2cAF1fVqqq6HbgYOH4zyj0PODpTGDxqICZJkrYle1XVLQDNzz0n2Gdf4Ma+58ubdWM+3HRLvqkv2PrlMVW1DrgTmL+pxjhYX5IkdWIag/X3SbKsb/WSqloy9iTJV4G9Jzj8DVOsZqJM1lhjX1BVNyXZGfgM8ELgI5s4ZlIGYpIkaWtzS1UdNdnGqjpmsm1Jfp5kn6q6Jck+wK0T7LYc6C9/P+CSpuybmp93J/k4vTFkH2mO2R9YnmQ2sCuwalMnYtekJEnqRg24TM/5wNhVkKcBX5hgn4uA45Ls3gzSPw64KMnsJAsAkswBTgCumKDc5wBfrymk/FoLxJIsTXJrkivGrX9Fc9nnlUn+uq36JUmSJvAO4Ngk1wDHNs9JsijJhwCqahXwVuDbzXJms24evYDs+8DlwE3AB5tyzwLmJ7kW+GMmuBpzIm12TZ4NvJ9eug6AJL9F76qCx1XV6rG5NyRJ0raniwldq+o24OgJ1i8DXtL3fCmwdNw+9wCHT1Lu/cBzN7c9rWXEqupSfrVv9GXAO6pqdbPPRP2ykiRJ24RhjxF7OPAbzYyz/5bkSZPtmGRxkmVJlt22csUQmyhJkoaho3nEZpRhB2Kzgd2BJwN/AnxqssnOqmpJVS2qqkXzFywcZhslSZKGYtjTVywHPttcRfCfSTYACwBTXpIkbUMGzW6ZEZuezwNPB0jycGAusHLIbZAkSTPAYF2TXbd6y2otI5bkXHqToS1IspzezTOXAkubKS3WAKdNZY4NSZKkUdRaIFZVp06y6ffaqlOSJG1FBknFjFj6xpn1JUmSOuK9JiVJUiccrG9GTJIkqTNmxCRJUifMiJkRkyRJ6owZMUmSNHxlRgzMiEmSJHVmq8iIzd4uzN95Xqt1tF2+pm6f3bbvugmSWnTkQ/dovY7XXnBVq+W/68RHtVr+NsN5xLaOQEySJI2WwntNgl2TkiRJnTEjJkmSOmFGzIyYJElSZ8yISZKkTpgRMyMmSZLUGTNikiRp+Aad0HXE5q8wIyZJktSRTgKxJEuT3Jrkii7qlyRJM0ANuIyQrjJiZwPHd1S3JEnSjNDJGLGqujTJQV3ULUmSZgavmpzBY8SSLE6yLMmyFStXdN0cSZKkLW7GBmJVtaSqFlXVooULFnbdHEmStAWN3WtycxdGLCPm9BWSJKkTdk3O4IyYJEnSqOtq+opzgW8Cj0iyPMnpXbRDkiR1ZNCpK0YrIdbZVZOndlGvJEnSTOIYMUmS1AnHiDlGTJIkqTNmxCRJUifMiJkRkyRJ6owZMUmS1AkzYmbEJEmSOmNGTJI0ct514qNaLf+1F1zVavnQ/jkA3Hnv2lbLX7dh8uzV2C2ONteoZcQMxCRJ0vANOjnraMVhdk1KkiR1xYyYJEnqhF2TZsQkSZI6Y0ZMkiR1woyYGTFJkqTOmBGTJEmdMCNmRkySJKkzrQViSZYmuTXJFX3r3prk+0kuT/KVJA9uq35JkjST9SZ0HWQZJW1mxM4Gjh+37m+q6nFVdRjwReDNLdYvSZI0o7UWiFXVpcCqcevu6nu6EyM3P64kSZqSmsYyQoY+RizJ25PcCLwAM2KSJGmIkuyR5OIk1zQ/d59kv9Oafa5JclqzbudmeNXYsjLJe5ptL0qyom/bS6bSnqEHYlX1hqraH/gY8PLJ9kuyOMmyJMtWrFwxvAZKkqTWFXQ1RuwM4GtVdQjwteb5AyTZA3gLcCRwBPCWJLtX1d1VddjYAtwAfLbv0E/2bf/QVBrT5VWTHwd+Z7KNVbWkqhZV1aKFCxYOsVmSJGkYBgvEpl3tycA5zeNzgFMm2OcZwMVVtaqqbgcuZty49ySHAHsC/z6dxgw1EGsaPeYk4Oph1i9JkrZ5e1XVLQDNzz0n2Gdf4Ma+58ubdf1OpZcB6w8Nf6eZHeK8JPtPpTGtTeia5FzgKGBBkuX0UnzPSvIIYAO9dN5L26pfkiTNcIOkt3rH7JNkWd/aJVW1ZOxJkq8Ce09w9BumWEsmqnnc8+cBL+x7fgFwblWtTvJSetm2p2+qotYCsao6dYLVZ7VVnyRJ2mbcUlVHTbaxqo6ZbFuSnyfZp6puSbIPcOsEuy2nl0wasx9wSV8ZjwdmV9V3+uq8rW//DwLv3MQ5AM6sL0mSulAFtWGwZXrOB05rHp8GfGGCfS4Cjkuye3NV5XHNujGnAuf2H9AEdWNOAq6aSmO816QkSdqWvAP4VJLTgZ8CzwVIsgh4aVW9pKpWJXkr8O3mmDOrqn9u1P8JPGtcua9MchKwjt48qi+aSmMMxCRJUjcGugRyepdNNl2IR0+wfhnwkr7nS4Glk5TxkAnWvR54/ea2x65JSZKkjpgRkyRJ3RhkvNf0x4jNKGbEJEmSOmJGTJIkdWPwecRGhoGYNsv9a9a3Xsfc2e0narfbbqK5+rasO+5Z03odu+00t/U6JP2qd534qNbreO0FU5r9YFpe+ZQDWy1//YaNBU1l1yR2TUqSJHXGjJgkSRq+wowYZsQkSZI6Y0ZMkiR1w8H6ZsQkSZK6YkZMkiR1wKsmwYyYJElSZ8yISZKkbjhGrL2MWJKlSW5NckXfuj2SXJzkmubn7m3VL0mSNNO12TV5NnD8uHVnAF+rqkOArzXPJUnStqg2DLaMkNYCsaq6FFg1bvXJwDnN43OAU9qqX5IkaaYb9hixvarqFoCquiXJnkOuX5IkzQRVjhFjBg/WT7IYWAyw/wEHdNwaSZK0xQ00fcVoBWLDnr7i50n2AWh+3jrZjlW1pKoWVdWihQsWDq2BkiRJwzLsQOx84LTm8WnAF4ZcvyRJminGuic3dxkhbU5fcS7wTeARSZYnOR14B3BskmuAY5vnkiRJ26TWxohV1amTbDq6rTolSdLWwlscgbc4kiRJ6syMvWpSkiSNOKevMCMmSZLUFTNikiRp+ArHiGFGTJIkqTNmxCRJUjc2DDDea5BjZjADMUmS1AGnrwADMW2mFXevbr2OHee1/7Gc/6C5rdexZl37Xxb3r1nfeh1t237urK6bIG22O+9d23odr3zKga3X8b5v3tBq+bfes6bV8keBgZgkSeqGGTEH60uSJHXFjJgkSRq+wgldMSMmSZLUGTNikvT/2rv3YKvKMo7j318iCVYgKEWComVk4xgpOZZKBmpihllqmjU02jg1XbxMF42mMGtGU7OaZixTy8rIMjPTNAg1ayZRQUAQCk3S4wWkvFTOiJenP9736Pa4D5y997vO3mfz+8zsOev6PGvttfY+737XWu9rZm3Q5FOT+B4xMzMzMyvANWJmZmbWHr5HzDViZmZmZu1SWY2YpEuBw4H1EbFHnnYFMDkvMhp4PCKmVLUNZmZm1sHcjlillyZ/DHwP+EnvhIj4UO+wpPOBJyrMb2ZmZtbRKiuIRcQtkibVmydJwDHA9Krym5mZWScL3yNG+27WPwBYFxFr+ltA0knASQATd9ppsLbLzMzMBkPQ5KXJ7iqItetm/eOAeZtaICIuioipETF1h+13GKTNMjMzMxs8g14jJmkY8AFg78HObWZmZh3ElybbUiN2ELA6InrakNvMzMysY1RWEJM0D/grMFlSj6QT86xj2cxlSTMzM+t2uYujZl5dpMqnJo/rZ/rHqsppZmZmNpS4iyMzMzNrD98j5i6OzMzMzNrFNWJmZmbWHu7iyDViZmZmZu3igpiZmZkNviDd79XMqwWSxkhaIGlN/rtdP8vdIOlxSdf2mb6LpEV5/SskDc/TX5nH78nzJw1ke1wQMzMzszZoW/MVpwMLI2I3YGEer+dc4KN1pp8DXJDXfwzobZ7rROCxiHgjcEFebrN8j5iZmZkNvuee5pkHbmp4tXj+eYBnWsh8BHBgHr4MuBn44svyRCyUdGDtNEkCpgMfrll/LnBhjjs3T78S+J4kRWy6Cm9IFMSWLFm8YcTW+mcDq2wPbKhqe5yj43J0wz44R2fl6IZ9cI7Oib8l59i5vxmx8cldn/vX3WOa3I53S7qjZvyiiLhogOu+NiIeBoiIhyWNayDvWODxiHg2j/cAO+bhHYEHctxnJT2Rl9/k+zUkCmIR0VCv35LuiIipVW2Pc3RWjm7YB+forBzdsA/O0TnxnaO+iLgPuK/J1RcD5/U3U9IfgdfVmTWnyXwvhK4zLQYwr19DoiBmZmZmNlARcVB/8yStkzQ+14aNB9Y3EHoDMFrSsFwrNgF4KM/rASYCPZKGAaOAf28uoG/WN7nMprYAAAmBSURBVDMzsy3JNcDsPDwb+O1AV8z3e90EHFVn/dq4RwE3bu7+MOjegthArxM7R3fk6IZ9cI7OytEN++AcnRPfOTrL2cDBktYAB+dxJE2VdHHvQpL+DPwKmCGpR9J78qwvAqdJuod0D9glefolwNg8/TT6fxrzJTSAwpqZmZmZVaBba8TMzMzMOp4LYmZmZmZt0nUFMUmnSlopaYWkeZK2qSDHoZL+lrsxGNA14AbjT5a0tOb1pKRTCsS9VNJ6SStqps2V9GBNrsMKxz86H4/nJbX8yHM/Od4q6a+S7pL0O0mvaTVPn5wn5/NpZYnjUCf+NpJuk7Qs5zizdI6cZytJd/btrqNg/LX5GCzt075PyRyjJV0pabWkVZLeUTj+REk35dgrJZ1cKG6983aKpFt73y9J+1SQ49z8Xi2X9BtJo0vnqJn3OUkhafvS8SV9Jn/nrpT0zWbjDzRfFXElnZWPw1JJ8yW9voIcA+q6p8UcV9T8v1graWkrOQyIiK55kRpTuw8Ykcd/CXyscI6tgHuBXYHhwDLgLRXu01bAI8DOBWJNA/YCVtRMmwt8rtC21ou/OzCZ1HLx1Ipy3A68Kw+fAJxV8P3fA1gBjCQ19/JHYLfCx1jAq/Lw1sAiYN8KzqXTgJ8D15aOneOvBbavInZNjsuAj+fh4cDowvHHA3vl4VcDfy/x+e7nvJ0PzMzDhwE3V5DjEGBYHj4HOKd0jjx9IvAH4J+tnAP97MO78+fulXl8XMHjXXd/Kjrer6kZ/izw/QpyfBM4PQ+fXtXxrpl/PvCVku/dlvjquhox0j/LEUpteIzkxfY9StkHuCci/hERG4FfkLo1qMoM4N6IaKRngboi4hYG0KZJyfgRsSoi/lZlDlJB75Y8vAD4YKl8pILkrRHxVKQ2Y/4EHFkwPpH8N49unV9Fn6KRNAF4L3Dx5pbtVLmmcxr5CaWI2BgRj5fMEREPR8SSPPwfYBUvtprdStx6520AvbW3o2jxu6qfz9/8eLEF8FtJbR4VzZFdAHyBFs/bfuJ/Ejg7Ip7OyzTS5lMz+SqJGxFP1oxuSzXv1RGkHyvkv++vIAfwQlc/xwDzWslhXXZpMiIeJLW0ez/wMPBERMwvnOaFLgyy2u4NqnAs1Z/on85V5pe2WpXdJiuAWXn4aNKv85Kxp0kaK2kkqeaiZHzghcuGS0kNCy6IiEWFU3yb9I+y5d5yNyGA+ZIWSzqpgvi7Ao8CP8qXWC+WtG0FeQCQNAl4G6mGsgqnAOdKeoD0vXVGRXl6nQBcXzqopFnAgxGxrHTs7E3AAZIWSfqTpLdXlKdykr6Rj/fxwFcqSPGSrnuARrruadQBwLqIWFNhji1CVxXEciHiCGAX4PXAtpI+UjpNnWmVtAEiaTipgPGrKuJnFwJvAKaQCq/nV5irKicAn5K0mHQ5aWOpwBGxinRJZwFwA+lS9LObXKm5PM9FxBRSjcU+kvYoFVvS4cD6iFhcKmY/9ouIvYCZpOMxrXD8YaTLJBdGxNuA/zHAdnoaJelVwK+BU/rUZJT0SeDUiJgInMqLbREVJ2kO6by9vHDckaQuY6ooVPQaBmwH7At8Hvhlro0ZciJiTj7elwOfbvf2tOg4XBtWRFcVxICDgPsi4tGIeAa4Cnhn4Ry9XRj0qu3eoLSZwJKIWFdRfCJiXS4EPA/8kHTpdUiJiNURcUhE7E36Yri3cPxLImKviJhGqqav7BdgvtR2M3BowbD7AbMkrSVdSp8u6WcF4wMQEQ/lv+uB31D+XOoBempqC68kFcyKkrQ1qRB2eURcVTp+jdmk7yhIP7Yq+exJmg0cDhwfEaV/NL6B9MN3WT6/JgBLJNXr469ZPcBV+RL+baRa3aYfCOgQP6fsLRS91il12YMa77pnwPKtPx8Arqgi/pam2wpi9wP7ShqZfzHNIN3jUdLtwG6Sdsk1VseSujWoQuW/OHo/tNmRpEtxQ4qkcfnvK4AvA9+vKP5OpC+fosdE0g69T7NJGkH6QbG6VPyIOCMiJkTEJNL5emNEFK0plrStpFf3DpNuEi96LkXEI8ADkibnSTOAu0vmyN8blwCrIuJbJWPX8RDwrjw8nQoK+JIOJbUCPisiniodPyLuiohxETEpn189pIcdHimY5mrS+4OkN5Ee0thQMP6gkLRbzegsCn7GazTddU+DDgJWR0RPRfG3LO1+WqD0CziTdIKvAH5KftKmcI7DSE9T3QvMqWg/RgL/AkYVjDmPdPnxGdIX5on5PboLWE76EI8vHP/IPPw0sA74QwX7cHI+Hn8ndVWhwsfiz6R/+MuAGRUc6z2BO/MxWEGFTyEBB1LBU5Ok+7eW5dfKCj8XU4A78nt1NbBd4fj7k241WA4sza/DCsStd97uDyzO79kiYO8KctxDuqe1d19afVLvZTn6zF9La09N1tuH4cDP8mdjCTC94PHe5P4U3o9f531YDvwO2LGCHGOBhaRC/UJgTBXvD/Bj4BOljsOW/nIXR2ZmZmZt0m2XJs3MzMyGDBfEzMzMzNrEBTEzMzOzNnFBzMzMzKxNXBAzMzMzaxMXxMysYZJC0k9rxodJelTStQ3GWStpk41z1lsmtxV4naTVklZKOruxPTAz6wwuiJlZM/4H7JEboAU4GHhwkLfhvIh4M6k/yP0kzRzk/GZmLXNBzMyadT3w3jz8kl4gJI2RdHXuTP5WSXvm6WMlzc+ddv+Amr5bJX1E0m2Slkr6gaSt+kscEU9FxE15eCOpoc8J5XfRzKxaLoiZWbN+ARwraRtS7wCLauadCdwZEXsCXwJ+kqd/FfhLpE67rwF2ApC0O/AhUsfhU4DngOMHshG5e6j3kVoSNzMbUoa1ewPMbGiKiOWSJpFqw37fZ/b+5E6NI+LGXBM2CphG6q+TiLhO0mN5+RnA3sDtqbtHRjCADotz58PzgO9GxD9a3Sczs8HmgpiZteIa4DxSH5Zja6arzrLR528tAZdFxBkN5r8IWBMR325wPTOzjuBLk2bWikuBr0XEXX2m30K+tCjpQGBDRDzZZ/pMYLu8/ELgKEnj8rwxknbeVGJJXwdGAaeU2RUzs8HngpiZNS0ieiLiO3VmzQWmSloOnA3MztPPBKZJWgIcAtyf49wNfBmYn9dZAIzvL6+kCcAc4C3AknyD/8fL7JWZ2eBRRL2rBGZmZmZWNdeImZmZmbWJC2JmZmZmbeKCmJmZmVmbuCBmZmZm1iYuiJmZmZm1iQtiZmZmZm3igpiZmZlZm/wfiYdadzGwF14AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v1 = load_from_config('ap.ini', k=20)\n",
    "v2 = load_from_config('ap.ini', k=20)\n",
    "\n",
    "# prints statistical comparisons\n",
    "compare_models(v1,v2)\n",
    "\n",
    "# renders chart comparing two models.\n",
    "compare(v1,v2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Interpretaiton Notes\n",
    "A perfect alignment between two models would show a dark diagonal all the way down the alignment heatmap.\n",
    "\n",
    "If multiple cells for a particular row or column are dark, this indicates that the topic is captured by multiple topics in the other model.\n",
    "\n",
    "Note that often missing data is as illuminating as the presence of data in a model alignment: if a particular row or column is entirely white, then that indicates that it is not covered in the other model at all. When looking at data aggregated by year from a particular discipline, this may indicate the introduction or removal of a particular topic from discussion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# viewing particular topics in each model\n",
    "print v1.topics([26])\n",
    "print v2.topics([2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}