Create_Synthetic_dataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "funky-trinidad",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import time\n",
    "import sys\n",
    "sys.path.append('./BASS/')\n",
    "sys.path.append('./Utils')\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "from GMM import GMM_model, GMM_synthetic\n",
    "import bass as md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "perceived-string",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 7 50 5 10000 0.0 0.0 0.15 4 0.0 1\n",
    "\n",
    "Sigma = 7\n",
    "dict_size = 50\n",
    "wordlength = 5\n",
    "L = 5000 \n",
    "std = 1.0/4\n",
    "\n",
    "eps_true = 0.0\n",
    "p_d_true = 0.5\n",
    "p_ins_true = 0.2\n",
    "p_b = 0.0\n",
    "\n",
    "params_true = np.array([eps_true,p_d_true,p_ins_true,0,0,Sigma])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "alone-ceiling",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(121)\n",
    "alphfreqs = np.random.dirichlet(5*np.ones(Sigma))\n",
    "model_true = GMM_synthetic(params_true, std)\n",
    "w_dict_true = md.generate_w_dict(alphfreqs,dict_size,wordlength)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "greater-version",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([4, 6, 6, 6, 5]), array([3, 3, 6, 6, 5, 0, 4]), array([0, 2]), array([3, 6, 6, 4]), array([6, 0, 5, 5, 5, 6, 6, 2, 6, 3]), array([3, 5, 6, 6]), array([3, 1, 2, 2, 5]), array([4, 3, 5, 5]), array([6, 2]), array([5, 2, 6, 4, 4, 1, 3, 1, 5]), array([3, 2, 0, 5, 0, 0, 5]), array([4, 5, 1, 3, 6, 3]), array([3, 5, 0, 5, 4, 6, 4]), array([5, 6, 3, 1]), array([6, 4, 4]), array([2, 6, 4, 0]), array([0, 1]), array([4, 5, 4, 1, 5, 0]), array([5, 6, 5, 6, 3, 2, 3]), array([2, 3, 4]), array([0, 4, 6]), array([5, 6, 6, 6, 6]), array([3, 1, 5, 6, 2]), array([6, 5, 6, 6, 0, 2]), array([6, 3, 1, 3, 4, 4]), array([3, 5]), array([6, 6, 2, 6, 0]), array([3, 3, 2, 5, 1]), array([3, 1, 6, 5, 0, 0]), array([0, 6, 6, 2]), array([1, 4, 3, 6]), array([1, 6]), array([4, 0, 1]), array([5, 2, 6, 6, 1, 6, 6]), array([4, 6, 6, 2, 6]), array([6, 1, 6, 2]), array([6, 4, 3, 2]), array([4, 6, 0]), array([3, 6, 4, 5, 5, 1, 3]), array([5, 5, 3]), array([6, 5, 1, 5, 0, 4]), array([0, 6]), array([2, 1, 1, 1, 3, 0, 0]), array([6, 1, 6, 1, 3]), array([5, 6, 5, 3, 6, 3]), array([5, 1, 3, 3]), array([4, 3, 4, 1, 2, 4, 5, 1, 6, 3]), array([3, 6, 6]), array([6, 4, 6, 6]), array([3, 5, 3, 6]), array([0]), array([1]), array([2]), array([3]), array([4]), array([5]), array([6])]\n"
     ]
    }
   ],
   "source": [
    "print(w_dict_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "complicated-reform",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(121)\n",
    "P_w_true = np.zeros(len(w_dict_true))\n",
    "P_w_true[:-Sigma] = np.random.dirichlet(1.0*np.ones(len(w_dict_true) - Sigma))*(1-p_b)\n",
    "P_w_true[-Sigma:] = alphfreqs*p_b\n",
    "Y,words_true,Y_data = md.generate_Y(L, P_w_true, w_dict_true, params_true,model_true)\n",
    "lengths_Y = np.array([len(Y)],dtype = int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "opponent-pasta",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_dict_sorted = []\n",
    "w_dict_sorted = []\n",
    "num_instances = []\n",
    "idx = np.argsort(-P_w_true)\n",
    "motif_lengths = [len(w) for w in w_dict_true]\n",
    "lmean = np.mean(motif_lengths)\n",
    "\n",
    "for i in idx:\n",
    "    P_dict_sorted.append(P_w_true[i])\n",
    "    w_dict_sorted.append(w_dict_true[i])\n",
    "    num_instances.append(int(P_w_true[i]*len(Y)/lmean))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "absent-binding",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Probability  Number of instances                       Sequences\n",
      "0      0.075831                   87     [5, 2, 6, 4, 4, 1, 3, 1, 5]\n",
      "1      0.066811                   76           [3, 2, 0, 5, 0, 0, 5]\n",
      "2      0.064301                   73                 [6, 6, 2, 6, 0]\n",
      "3      0.058106                   66              [4, 5, 4, 1, 5, 0]\n",
      "4      0.042684                   49                    [5, 6, 3, 1]\n",
      "5      0.042446                   48                       [5, 5, 3]\n",
      "6      0.038826                   44  [6, 0, 5, 5, 5, 6, 6, 2, 6, 3]\n",
      "7      0.038707                   44              [6, 3, 1, 3, 4, 4]\n",
      "8      0.037651                   43           [5, 2, 6, 6, 1, 6, 6]\n",
      "9      0.033304                   38           [3, 6, 4, 5, 5, 1, 3]\n",
      "10     0.032497                   37                          [3, 5]\n",
      "11     0.029969                   34                    [4, 3, 5, 5]\n",
      "12     0.029939                   34                    [0, 6, 6, 2]\n",
      "13     0.029192                   33                 [3, 3, 2, 5, 1]\n",
      "14     0.027206                   31              [6, 5, 1, 5, 0, 4]\n",
      "15     0.024497                   28              [3, 1, 6, 5, 0, 0]\n",
      "16     0.023352                   26                 [6, 1, 6, 1, 3]\n",
      "17     0.022614                   26                       [4, 0, 1]\n",
      "18     0.021332                   24                    [2, 6, 4, 0]\n",
      "19     0.019752                   22           [2, 1, 1, 1, 3, 0, 0]\n",
      "20     0.019048                   21              [6, 5, 6, 6, 0, 2]\n",
      "21     0.018345                   21                       [2, 3, 4]\n",
      "22     0.017770                   20                 [3, 1, 2, 2, 5]\n",
      "23     0.015694                   18                    [1, 4, 3, 6]\n",
      "24     0.015007                   17                    [6, 4, 6, 6]\n",
      "25     0.014739                   16              [4, 5, 1, 3, 6, 3]\n",
      "26     0.013870                   15                          [1, 6]\n",
      "27     0.011942                   13                       [4, 6, 0]\n",
      "28     0.011477                   13                    [3, 5, 6, 6]\n",
      "29     0.011392                   13           [5, 6, 5, 6, 3, 2, 3]\n",
      "30     0.011367                   13                       [6, 4, 4]\n",
      "31     0.010752                   12                    [3, 5, 3, 6]\n",
      "32     0.010472                   12                    [6, 4, 3, 2]\n",
      "33     0.007666                    8                          [0, 1]\n",
      "34     0.006256                    7                          [6, 2]\n",
      "35     0.006176                    7                       [0, 4, 6]\n",
      "36     0.006154                    7                       [3, 6, 6]\n",
      "37     0.005808                    6                          [0, 2]\n",
      "38     0.005603                    6                 [5, 6, 6, 6, 6]\n",
      "39     0.005183                    5           [3, 3, 6, 6, 5, 0, 4]\n",
      "40     0.003609                    4                    [3, 6, 6, 4]\n",
      "41     0.003443                    3                 [4, 6, 6, 2, 6]\n",
      "42     0.002585                    2                 [4, 6, 6, 6, 5]\n",
      "43     0.001669                    1              [5, 6, 5, 3, 6, 3]\n",
      "44     0.001652                    1                    [5, 1, 3, 3]\n",
      "45     0.001250                    1                          [0, 6]\n",
      "46     0.000855                    0                 [3, 1, 5, 6, 2]\n",
      "47     0.000404                    0           [3, 5, 0, 5, 4, 6, 4]\n",
      "48     0.000399                    0                    [6, 1, 6, 2]\n",
      "49     0.000394                    0  [4, 3, 4, 1, 2, 4, 5, 1, 6, 3]\n",
      "50     0.000000                    0                             [5]\n",
      "51     0.000000                    0                             [0]\n",
      "52     0.000000                    0                             [1]\n",
      "53     0.000000                    0                             [2]\n",
      "54     0.000000                    0                             [3]\n",
      "55     0.000000                    0                             [4]\n",
      "56     0.000000                    0                             [6]\n"
     ]
    }
   ],
   "source": [
    "full_dict = pd.DataFrame({'Probability':P_dict_sorted,'Number of instances':num_instances,'Sequences':w_dict_sorted})\n",
    "print(full_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "valid-language",
   "metadata": {},
   "outputs": [],
   "source": [
    "full_dict.to_csv('./Data/Synthetic_validation/Synthetic_truedict_condition1.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "golden-cannon",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.save('./Data/synth_dataset_condition1',Y_data)\n",
    "# np.save('./Data/synth_lengths_condition1',lengths_Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "decreased-quarter",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5003]\n"
     ]
    }
   ],
   "source": [
    "print(lengths_Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "seventh-separate",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x13b298e80>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9hklEQVR4nO2df5Ac5Xnnv8/MtqRZ+aJZmY2NxggBcUkxh7Vr9kBGqZSFXcggA2t+ycRU7CunSO7iukLm9m5dpowgrmNzKlvkYt85nOOKcxAsQGYjInIijpTylRwRr9hdy7KlBBASDCTIllYxuyNpdva5P2beUU9Pv91vT/fM9Ew/nyqVdmd6p9/p6Xne932e7/M8xMwQBEEQup9UuwcgCIIgtAYx+IIgCAlBDL4gCEJCEIMvCIKQEMTgC4IgJISedg/Ai4suuohXrVrV7mEIgiB0DAcPHvw5M/e7PRdrg79q1SpMTEy0exiCIAgdAxEd1z0nLh1BEISEIAZfEAQhIYjBFwRBSAhi8AVBEBKCGHxBEISEEGuVjiAoxifz2LbnKN6cKWBFNoORjasxPJhr97AEoaMQgy/EnvHJPL74vUMoFEsAgPxMAV/83iEA6DqjLxOb0EzE4AsN0yrjtG3P0aqxVxSKJWzbc7SrjGGSJjahPYgPX2gIZZzyMwUwLhin8cl85Od6c6YQ6PFOxWtiE4QoEIMvNEQrjdOKbCbQ451KUiY2oX2IwRcaIgrjND6Zx/qxvbhsdDfWj+3V7g5GNq5GxkrXPJax0hjZuNp8wB1AUiY2oX1EYvCJ6NtE9DYR/UTz/EeI6AwRTVX+fTmK8wrtI9trBXrcSRCX0PBgDo/cdhVy2QwIQC6bwSO3XRXIr206ubSTbpnYOuFaJ5WogrZ/BuDrAP7c45j/x8yfiOh8QpvRtUI2bZGscwnd/9Q0tuyYqgsCDw/mXA28SeC4U4KhaiydrNLplGudVCIx+Mz8AyJaFcVrCZ3BmUIx0ONOdK6fUmXGMDEUpsal2SqfKNVKuokt7qhrkHf5XLtRUdWptFKW+WEimgbwJoD/zMyH3Q4ionsB3AsAK1eubOHwhCCsyGZcv9xu/mY3g6j7ezt+hsLUkDczGNrsFW2zpK9Rvq7zGrghged40Kqg7UsALmXmtQD+GMC47kBmfoyZh5h5qL/ftYa/oKGVvlNTf7POV79hTX/d37vhZShMDXkzg6HNVCs1S/oa9eu6XQMnEniOBy0x+Mz8r8z8TuXn5wFYRHRRK84dd6Iy0q3UxSsW91y4ffp6LddAqs4g7jtysiYQmyZyPYeXoTA15M0MhjZz99CsySTq1/XbqVkpwtz5eQnixoCWGHwiei9R+RtNRNdUzvuLVpw7zkRppFupi1fjnrH5688WF1yP9TKIw4M57B+9Hts3D+BXMvXeRT+jbGrIw6p8vCblZu4edNcuP1MItUiIepLSTdYAkM1YAAGn54otW4gIeiLx4RPRkwA+AuAiInoDwIMALABg5m8CuAPAfyCieQAFAJ9iNtVzdC9RBhObsdLU+XmDjNvP16/z/2YzFj6x9mJs23PUVbUD1Kpa8jMFpIlqJjnnsY34qP189CMbV9eNP6rdg+7aES6squ3jAcwUPkHiLyaUPL7KSxf31CwMAAnitpNIVvjMfDczX8zMFjO/j5n/lJm/WTH2YOavM/OVzLyWmdcx8w+jOG+nE6WRjnql6bX7CDJuv1W4zv9LBOw8mPfd/Sijm7HSdQqfKFaRfjunILuHoO47t2tHAJzmtVAsYeuuw8a7Rd3rqp1D0OuW09xjuWxGsodjhhRPayNRrrSiXmnqDN1Dzx1Gish1VbcsU5905act1/l/T8/Vyzt1K8Nmyi5NDJbJ7sFUzePcVd1+dQ77jpys/q67Xs5VNKC/Bs6dkX0ScY7LS81jl2I6JyJ17+mkmhLEbQ9i8NtIlEY6iHvDBJ2hczPEitnz8xifzLsaGN3505rJQ0d+plB3jma6s3Qjs7ukTNwoJpOS26Sw82C+Zsewfmyvb5DUju4aqM/E7fXs941uknI+x7iw+8g5rkOzXF5CcMTgt5GoMyvdvmCN6sJNdPJOiiUOvKoOYuwVzvcTtU/aT1dupQkjG1cH0uB7BWAvG92NFdkMZs/N+04KIxtXY+TpaRQXLlw3K0V415Ie18nY7xp4TZZ+7iznc8rY7x+9vvpYN2QPdxNi8NtM1JmVUbk3Rjauxsgz0yiWghlkZUBMV765BiYWNyPYbHeWnaWLeqqr47DBawBVn7sO+zXduutwjbEHgAUAmz54MXYezAe+Bl6TZSM7J7fnOjV7uBuRapkdhEnQLyr3xvBgDksXBV8PrMhmAslNN6zph17Up8fpQw9bXE332m6o8hFBg9eNvE+g9pq6+epLC4xnX8rj9qtzVYlkmggfWrkM2/Yc9bxfvILqXkIA3XOmxfOE9iAGv0MwNaJB1Dp+E4hpXRyFPVBnkhMwPpnHzoN5rZ/cC/v7sQcPU0TIV9wRD4wfakiv7ucGUc8HudbDg7mG3qfXNbUze76Exw+cqLrISszY/8opI5WTbrJ0mwysNGH23Lx2R/LO2XnR2McYMfgdgqkRDVvywP5l1a3W0kQglLXyfb1WnaEwXfmapOS7YX8/9vcB1BZfe/zAibr3ZzIJuF1Dt3MHzeDVyRft9PVarsY3Chmj835RE/59O6bwz2fO1k1Izsmgr9cC2F0RpCguMLbuci2TJcQA8eF3CKZG1DRI5jeBPPTcYdcgoJUmbLtjrafOXLeSTRHVKGwaMWIE4ParczXv03TSKBRLeOLACa0EUeGmeCox16lPggYk3WINdjJWGg/efGWgZKmgqCDxsoyF2fPz1RiNrkqp3f++fmyvp0pLMVMo1nzO0pg9PlCcE16HhoZ4YmKi3cOIBTo5nlMVYcplo7u1hjljpbVGSUnvnEYQqJff6V5brVqDSgwV9vfs9T4aeb1mYzd+yzIWiICZuaKvITSpSNkM+nqt6iQU5Fqra+o2bvs9IEQPER1k5iG352SF3yFErUTRrRhTVC+3s6O+8M4V4RIrZWSM7CoWvxWvDvu4s72W0arT9PVMaHTF6vZ3wIUdwrY9RzFx/FRNotWGNf3V35dlLCyxUqHfbxBOzxUx8sw0gGC7DJW1G6Y+vuwMokdW+B1Es2uYW2kKLMNslJyLwTM9c5oIrzxyE8Yn8/jCjim4l20zhwBs3zzgeS29skp1CUfOv6+73ikCCIGuuZUm9KQIBU2xumaRy2aw6t0Z7H/llNHxbiUgnM8fG9ukfV52Bo0jK/wuIUo9s5v/+dTsuZYZfLUzeOS2q6ruFFMXj9pdbNtz1NPYr79iuZGB4sprmbpTnFfILyagxurcyTj19CYUS9yyz8hOfqYQaCfkN0I/JVSzu5QlFVHpJBhVnvjY2CaMbFztuWr005BnM5ZRQxM7TtXIyMbVsNL+avVcRZfuZ4BeOnGmrCwxwCuAHDQwrN6TXfYaRcC1GVhtsAAmrkgputYcxOALAMqqHB3ZjIXtmwd8DfrtV+eq0kPTJKO6L7DB4nVm7jxGnp72Pa5QLOFssWQ0EXmtOIMamTcr9X7sste4Mt/iwZkmxTWzz0CSEYMvAPAuirb1litrNNlAObhrZ6ZQxM6DeYxsXI3XxjZh++YBo/Pai5Dd/9S0kZtj9nzJ2B1SKC5UJyKVO5B2DN5eG8dNox/UyKzIZhrOMYgCIuDRzQPlGIEPrQrhZaw0Ht08gP2j1xu5ZJrZpSzJiA9fMMapyXarsnjfjils23PU+Is5snE1Hhg/VKOPj5p9R05W4wTjk3mMPD0Nuykulhhf/N6PMb9wwT+enyngvh1TeOi5w9j0wYsDja/d7puqEW+0lkNEpImwwNyQwECKrjUHUekkHLv6RIdSn9i/fCZ9TBsJSjaLRysqnEa0/3F7LyY0UpSuGbzmocQRmoOodARXTJN5lPrEXgbYT3bXagPpZ5SVeqaRoF+nGftFaYqFsSegrRm3ouOvR1b4XYzfDW+62g3apCSuqPhDHIxhMzHZkeSyGbx5ptAyH77bAiFjpes6ekVllJOs4/da4YvB70BMVi4mN/yq0d2+50oBoROb4obf7qSTyViplidlhUGXxKarX2RK1KVIOommu3SI6NsAPgHgbWb+ty7PE4A/AnATgDkAn2Xml6I4d9Iw7bLkVxzNWWXTCQGVcgmdYzxMicLYx3XS6LTPS5fEpivmZoro+N2JSpb5ZwA+7vH8jQDeX/l3L4D/FdF5E4dpmWSvlnr2csI6jo1twvlWi7RbTMZK1Zb+NSRF8TT23Yrb/e2H6PjdicTgM/MPAHjlsN8K4M+5zAEAWSK6OIpzJw3TlYvuxlbNzb3IZsrGrxv89l4UigsY2bgax8Y2odeju5eaEAjlgGiHxXBjSxDVaNC4i+j43WlV4lUOwOu239+oPFYHEd1LRBNENHHy5MmWDK6TMF256G54EyNObdZvtxK1cvTa6jOAs8UFXHfFcpxvQx2bbqSv18J1Vyw3NvpK8aPw69YWddvLbiF2mbbM/BgzDzHzUH9/f7uHEzvcesC6rVx0N7xJ16WZFpbfbTfK0Ptt9QvFknGlSB2N1BvqRlSN/ZdOnDF2jakCd4B5u097rSjTDN9up1U6/DyAS2y/v6/ymBAAtx6wzg5QdnTVNf2098syFtaP7Y1gxPFHvVeT3AJT1IRsf62MlQb59BroRpzXVHX1aqT0hJqcpZJm47Rqhb8LwG9TmXUAzjDzWy06d9fgdqMzyqUDTFErf+Wnd2KlCLPn9U2quw37e2VEU42AUa6vb7/GhWKppY1L4gABuO6K5TWB8cU9KWzZMdXQ/aV2YV6ChKAN65NGVLLMJwF8BMBFRPQGgAcBWADAzN8E8DzKksyXUZZl/vsozps0opSanZuvl+8ppUqSDJOztnwUK/xcNoOJ46c8m30nAUa5RPUjt12FieOnQtVLIqDqtvQq7WF38QDBpJxJIBKDz8x3+zzPAH4/inMlGd2NHlRqpttO9y7qSbxOOSyEcpzliQMn2j2UWKAK6png5VK77orlVeNt0hpTXDzuSC2dGOGXQet2o1tpwuy5eVw2uruuB6ouC9drpxCkb6lQj3KxiZbHHAKq96qukN9rv7jwmLOSpu5ay+KlHjH4McEkg9Z5o2d7Lbxzdr7qOsjPFPC4bWWp29p67RQabSwulCESQxMEZ6mDLZrdgPOa+pXqBiTJyo3YyTKTimkGrV1q1ruox7dIlttr+CWlLO6R26JRmIHeRSK9NMUpJ24kQ9btflauNaEWWeHHhEYCsqYrSbfVEVDfXALwl2wmHStVLurlNc/OnpfrZ8p9O6Zw/1PT1UJpG9b0Y+fBfF3RP68M2eHBXF1QmAHsPJjH0KXLxY9vQwx+TGgkIGvqb3d7DTeN/vqxvWLsPSAAm6+5BMdOvhM6CUu4gL1Q2s6D+YZKJrvFTSRwW48Y/Jjg5jv3W9mY+Nut1IV+rSogZi89u2FNP3b/+K1ESTEbhQHs/vFbicpEbjWFYqmmJaUpUh3TDDH4McHuZlFG2e5/12XS2v/GjUU9KWzddbhGE25fUT0u8sFANDoxqv6uqS5pJtNMGjHSUUmWux2JzsWI4cFcNQDlrAeuyxxUQVxdhujs+VLiE4DiQKnSzPvuay+BlU5QdboGaMRIS3VMM8TgxwxTtY6TbIB67qZkMxZSYpsiIz9TwI4fvQ7LcVEzVqpafjndxaVKey1/c6NckEGR6phmiEsnZjTiixyfzOOds/ORjkNV1ZTdQbQUS1xTziGdIpwtLqBQXOia3sE65gy6cb1rSU/DRlpXLFBHEpuci8GPGaa+SPvNGrVfWK2ydEkwnYxbJct2UrLpO7vZ2Jtyeq6I9WN7m2Z87eIFeymHpNTfEZdOzDDxRTrrgTdqKPp6LdyzbmVNe79sxsK2O9dieDCHjMEWvNPYvnkAx8Y2BWppKLQWv7hVo9i/N0D9pN9IK8VOQ1b4MUOXFKUeH5/MVxNVgmKXY9pf8yvDV7keX3CpqNnpqPf84M1X4gtPTUm7wpjSDA29SQ3+bpdxisFvIzofos4X+cD4oYZKzGasNB65rWzU1fm85J6KbvMw5LKZmi29EG+iNr4mr9ftMk4x+G3CpFia83g/Y6+03ssyFojKrQp1ZRPs5wPcdxTdFETMWGlsWNMvpSM6iKiNr19mehJknGLw20TQNm3b9hz1NPZqFT9x/BSefPF1lJiRJsKGNf0YHsy5lk0oFEvYuuswzs0vuE4E6y7v64oSAsqF1UhbPbfXAiA7hCZjpRuTZyrcds9umekqcOt0c3YrYvDbRFD5pdd2NE1UNfb2zNkSc/V33d+7yS7jHLy6Z91K4+zge9atrMYnxifzoY20WgE6r3MjEMr6eyVV7Ou1sOmDF0vmc4WlixqXZ45P5jHy9HS1kmx+poCRp6ex7c61eOS2q3ylmN0s1xSD3yb85JfOmy7ba7mm9ROAr95VVtXc/9S067mefPH1wI1N4raCJQCfrhhwE6NIAJ44cAL7jpzEhjX9+IuQhlRNqsODucgmw5/+wY3Vn9Xn3ckQyqWho6gWeiZE/sfWXYfryoYXFxhbdx3G1IM3eBrvoK7WTqP7dHcdgpf80im7zM8U8M7Z+bqUfGUE1Y2o87eXmLXn85InxinnU5W7HZ/MY6lBvXnGhf6mjx84gTB6o4yVxlfvWgtA32wjKCmiasPtB8YP1cgFO5FcNoPtmwcwF1Fp6DD+e12yoEkSYaOZ7p2CrPDbhJf80s3fXlxgZDNl46xu3GyvhaFLl1eP0QVZCRduZKc0E9DXwI9buFb1R21FfkCKUJVsLu5JYeL4qbo67U6cPVkJwK/96lK8/PZs3bW010oK09w7LuRnCrhvxxRSFF7dZRI8bZbbpdurbkZi8Ino4wD+CEAawLeYeczx/GcBbAOgMim+zszfiuLcXkRxUzTTn6eTX3r52+2r9NNzRYw8M42tuw7jTKGInhRQcvmyqZUuUDY0KiA2PJjD+GQeS6xURylXCgYp+mGxewRmCkVfo+zWgJsB/NPbs77n6nRjbydsXoNJ8NTP7dKncX+aJNt1e9XN0EslIkoD+AaAGwF8AMDdRPQBl0N3MPNA5V9LjL3TLRI0ey+K12gE3c2lSibbKZYYM4UiGICpHSyWGA89d7j6/jq1Fn7GSlV3PUC5n2yz8LJjuWwmtkbbSlNb/LZposAuQdUPeNueo3XfsfHJPNaP7cVlo7tx/1PTnm6XB2++ss79aaUJD958pe8Yur3qZhT3wjUAXmbmV5n5PIDvArg1gtcNRRS+uHb583Q3XZSa+NNzxUhkiu2kUFzA7LkLRePakTLQ12th/+j1sSzVkM1Y2HbHWnxt80DLy2SUmHFsbFOgv2GG68LKtJSI2hkPD+aw7Y61NZUzt92x1mhn3u1VN6Nw6eQAvG77/Q0A17ocdzsR/SaAfwSwhZlfdzkGRHQvgHsBYOXKlQ0PKgpfXLv8eU7/vkqkito4d4Nf0q+Je7NhLhukMzHcJS1dfEHaODyYw6rR3S09//qxvQ3/baFYwpanpvDQc4eNd6D2nbHzO+SXWe7lulW7i26QabYqaPscgCeZ+RwR/S6A7wBw7WHGzI8BeAwAhoaGGv42R+GLa6c/T/n3nf7KKLHrwIPi5rNOImcKRWzddTiUCqhZOCf0XEBpbljCnovZvMOYs46+swxJfqaALTumcN+Oqbo4gVdMANBnqHei0Y9in5cHcInt9/fhQnAWAMDMv2Dmc5VfvwXg6gjO60kUvjjTypXKt7h+bG9k/n31uvftmNIa+75eq66ZRhDmigsNNzjpNZBGhiUustA0UU2swM6KbCa2PQMYqLknN6zpb++AmontZtGVIbEb//t2TGHgoReqK3ud67bbZJrEIR2fRNSDspvmoygb+h8B+C1mPmw75mJmfqvy8ycB/FdmXuf32kNDQzwxMdHw2Jqt0nFbfbsVKgua0WeyqicAx8Y21bxGo59kxkpVlS9LI0qc6RasNGHbHWUNvi4tP+5krDRuvzrnKys1IZfN4PTsucA7w1bUZYq6b7BX7wT1/YsjRHSQmYfcngvt0mHmeSL6PIA9KMsyv83Mh4noYQATzLwLwH8iolsAzAM4BeCzYc9rQtAOOEFfQzf7e9WnMdlGmgRTVUtD9XrORuVBsMscTYx9qwxdNmO1f/VceaPqOtt9yp1g7IHyPRlVyYYNa/rxRMDXShPhlUduAlDbgCRqlJGPamJZ4TG5dapMM/QKv5mEXeE3m8tGdwf+0iv/oc5A6zTETlIAvrZ5AIA+cSos7V7tK+VLu2Wj6jMLEkDsVhotHvdaZTdqr3HTKPakuGaRsdL40MplrsUDUwR87a6B2PrwvVb4UlrBBVO/fCOzfH6mgJFnprUrV1ODsgBofYxRcba4UNcRq5XMzBVdNdWtRu2+km7sgcr9u3F14M9kfDLvWuOmEZpp7O1SzB9qKsUyOjNgC4jBryNIstWGNf11gUW/+jQAappYhyE/U2iq6qLEjJ0H8zjbJq3+imzGVVN9z7rG5bqN4Jbw1g00Mo2mVXZbwFt4256jLXfPqeQv03sml83g2Ngm7B8tCwh1bzHGThFfpJaOA9M69eOTeew8mK+rnXL71TkMXbq8axpttPM9nJ49h4GHXsCZQrmRy/bNF7bRrSwj3GlNYNIEvHeZtwSzUYlmiRnb9hwNvFJvR2G4u6+9pKZ9574jJ7XjSAHVwoV+MYZ0M1O6m4ys8B2YJlu5TQyM8k1lz9YTGmeuuFAtG6F01A+MH/L9uySTThG+etcA9o9ej0c3D7jKih/dXH6+0fuzU6p6PvkPr9fszN1k1op0mjBx/JRR1dK7r73E8/k4Iwbfgc4vn+21avz6upsiP1PA+rG92LJjCkC5CYfuJutm0kRGZYyDwCiv7Mcn8y0vFdApLCwwtuyYqma5epUJGNm4OlQeR9wpLZRrRimGB3O4/Wp333uxxHjyxdd9d7QE1FSo7TREpePATQNvpQng2jR+nTTR+biK9v/wlVMdI+OLO9lKqQkJonqTAoCKooUIyPSU8y3sfY5HnpmOLKYUV16r6OWjylq30oSli3qqrsa4lVrwUumIwXfBmRA1e27eNeDkVv/c7Wp2UzPwuNApSU9Rk81YWLq4J7RbhYBKWew4FoWIFntPYz/ffCPfU5VsGRejn0iDH2Udey+9fS6bqZ6jU3yb3UBSJ1HpfdsYfgsEK03Y/O8uaTgbOZfNVNU97SZxOvyo69jr/PrqQz42tgkjG1fHpvZLEkiisQfKbqydB5vbj6Eb8btblFvLGfMwzUHplMqzXWnwoy54ZFJEbdueo4l0MUSBlSI8unkAj24eqH7ZshlLJlAN3SD3bRZhJJOqZIRaxO0fvR4P3nylkeiiU0otdKXBj7qOvUlTBHHnNM67lpTTQey7MiXHFIQgLFT6NTcCo1wrya7GA8qrfq+JpJM6YnWlD18nm2ymn+2KLz6fWDeDnaQGU+NICohlnf5mogK091Vk0WGxV791U/hkMxa23nJlbAK2QAJ9+O3oSynGvsx1V3SuRrnbUMX1koKVpqo4I6o0DXuWvXOX/+jmAUw9eEMgY9+s/hmmdGVpBWd7syi0svaUa6UQsXfOaXU3objy96+6F5wSWkuuUocoTNnsjqOy5hqfzCNKtam9V25YG9Lu7lld6dKJGq+EDfuW7wtPTTW9bKsgmHDPupWetWO6maglu71WCj/9gxtDv06rXM2Jc+lEjVcJYvuWT2y9EBceP3AikcYeiN69OldciMT1ErWYpBG60qUTNX4fiPpied1nfb0WzhZLichsNKGv18LMXL0SJwWAUoSSbJWEGOGsluuFW9InAG3rxVZKOmWFb4DJBzI+mYeXBPj0XFGMPS7UYD+jkV1mFqWx7rK+Vg5JEHzJzxQw+PALvit9t6TPkaenMfLMtKuxb7WkUwy+AV5lVRUPPXe4oxsjNBOy/a8ukW4BP3u+5NpWThDazem5IkaemfY0+m7u3+ICuxaoSxO1vAZPYg1+EHmUkmR5IZUb9bDjfyF+SLVpM4olxtZdh7XPB/HHLzC3XL8fycdMRB8noqNE9DIRjbo8v5iIdlSef5GIVkVx3kYYn8xj8OEXcN+OqUC1dpT0UhC6EfE2mjNTKEbS57od5RhCG3wiSgP4BoAbAXwAwN1E9AHHYZ8DcJqZfw3AdgB/GPa8jaD8a26rcZNaO52SPi0Iii7ub9JW7LbC7i2YOz9f11TGSlFd0/d2lWOIYoV/DYCXmflVZj4P4LsAbnUccyuA71R+fgbAR4la3xjSS14J+G/HhgdzkXdxEoRmImKn5qBshTNIe3quCNCF4n+5bAbb7lyLbXes9azF1SqikGXmALxu+/0NANfqjmHmeSI6A+DdAH7ufDEiuhfAvQCwcqV/p/kg+Bl0ky3WgkRmBSHxKFvhGqQtMZYu7sHUgzfUPG5i4KPs4+FG7EI1zPwYMw8x81B/f3+kr+1l0E23WCKtbD6q0bbETIS4MntuHuOTeaNkKlOBSNR9PNyIYoWfB2Bv4/6+ymNux7xBRD0AlgH4RQTnDsTIxtWeFe+Acvqzml03rOnHviMna34XgqPkmLlsBqdmz3lOmn29Ft45W4ys2qEgNIOZQhFf/N4hZHst15igWlwGqZ/j1ccjqlV+FCv8HwF4PxFdRkSLAHwKwC7HMbsAfKby8x0A9nIbivh4VbybOH4KWxzKHZWern5/QtrKBaav18L2zQN4rdIVzMvYp1OE03NFUYy0mKTFddMRRbILxRKY4VmZN0gzplaUXgi9wq/45D8PYA+ANIBvM/NhInoYwAQz7wLwpwD+DxG9DOAUypNCU/DzgblVvBufzOOJAyd8deLivQ/OzFx5tf7Qc4fxztl5z2OlnEJ7SNpVj/I+O1MoYvvmAa3NCWLEdX2xo5RvRlJLh5mfB/C847Ev234+C+DOKM7lhcn2yW1CkPaEzUNdV0lMEzoFArDEShnF61ZUylDrXC5BjLibyzlq+WbsgrZh8Ns+6YIiSa0qKAhCPQxgiZWu09O74RfXC9KMyaSVali6qlqm3/ZJNyFEXT9bEITO5vRcEVaakM1YOFMoa+vdTMS+Iyc9XydoM6awTVb86CqD77d90k0IJWZkrHTNZEAot+t77ReF6gclOwFBSA7FEuOXZ+fx6XUr8bhGsJGfKeCy0d2ehrzZRjwIXWXw/XxgOqOds/nyvWbhwYdfEF90xBABK5bJZCrEkxKzrzrP7h4GWteusBG6yofv5wPz8qcND+awf/R6HBvbhP2j17t+aA/efGVdTQzBG7+rxVz2g8pVFeKKqbPXpB5Xu0lcT9uwqcv2ZuaCP/esW4mhS5d7JlLlKkltJtJYQYg7uWymaaURTPDqaZs4gx8l45N5bNkx5WqkJBBcZumiNObOlzwNOQE4NrYJD4wf0vpKhfZjb2DTDHqtFOY6POvOeY3sWeatMv7SxNxGkMYnfgwP5rB984Crm+juay9BNmOFHW7HM+tj7AFgWeU6+SkehPbCQM09vXSRmXTRhGzGwk//4EY8unkAmSZ3Ywk6YtPj3SZE9Xsz6uI0QlcFbf3QJWZNHD9VUzPHbya2u3XUSl79r9wTOw/mPUsxCxeYPe9diEpojKWL0iicLyGqNXNfr4WzthX47PlSJDGtFKFay6pM8yI6yiib7la8jrNShHct6cHMXNFIxRd1XZxGSJRLZ/3YXtcPxfmhZqy0NuHBOWnA5e/Exx8cVRlTrlu0pFMUSSkBK01YuqgHM4V6lVqKwtXdp8oXcEU2g9lz867niIIgLqk0Ed67bIn2fnRz0ejsi3MMx8Y2GY6iMcSlU0G3gnTeBF7Rdq8mKurvZKUanDdnChjZuDoyF4FQJrK6MQytIQ57CuYL0sZmGXsgWPzh7msv0X6PCXBV8rmpAJ20o62hnUQYfOW3D/KB6+pZ+83gyi0kBEPVJNl251r0Skft2FHs0sJ2hLKSLF1pwJcmwj3rVuIrw1dpv8fZXvfYnF0Wrl7bTrvaGtrpepeOlwsG0G/zctkM9o9e7/v3TlRt/SB/I1wgar+zIHihvudujE/mMfLMNIqlWgthpQjb7lzr64tvdvcqHV4una4P2nq5YHQBVr961l4QXci0u/+p6Y6SZqaAlhlanWx19rxMkkI9zZA5e624lbF2GnugvNsxCb7GqaSCousNvp8fDgCGLl0euJ61jplK6QX195200l9AeYfSTD9q9VwVRZMEaQU/nHWuwkKA54rbZFffqXG6rjf4JvWoG6ln7XU+++sCFyrlLWuRMXVDZf/5rZHOzbdmcpJidIIXaSIsMNf0rIjqfvFzr5js6js1Ttf10bEg9ahN/16HlSLMnZ+vSepSNXq2bx7A0sXm82uUSVuE8vswuUlb1aR9w5p+UJMEOSlKwI3dxVhpwlfvWltT1yrI99APvyQov9W73X5EmcjZCrp+hR+0HrXu7x967nBdpUylTT5TKGJZxsLs+fnqMfbqeUC9a8dPE7x0cQ+WLu6JZFWTsVKYOH4Kp2bPhX6tqGhmCQWVuSuVTTuTpYt66r6f9u9xFN8JZxKUPcCa8ogX2PX3QRqUx4WuV+lEiVfUXZd04ZVQ5BWIIgDbNw8EigGkCPjw5ctx+M1fNuQ6slLl8bRbgZexUli+dLHvl0/oTtS97/Zdi7J4oUqCMvHZuyVjen3ndcqfVpBolU6UePn6G+k472XIlC4dMF/VLDDw0okzoAZLXG27cy2enjiB/a+caujvo+KR2z5Ys/LqpMC3EJ5sr6UtgRJlyRK1E9T57J1xBOd3v5HvfLsJZfCJaDmAHQBWAXgNwF3MfNrluBIA5d84wcy3hDlvHPEKDs/ZXD0mKJ87UDvJXDa629eUN/plUIknL50409DfR0Vfr1XzxRoezGHi+Cn8xYsn2r7z6CbuWbcytuWozxZLdbGkQrGEJ198PdLdnooh6Qz0ArNnGYQgDcrjQtjY1iiAv2Xm9wP428rvbhSYeaDyr+uMPeAdHA56j153xXLXnUQzb6QSM+7bMdX2lfTZYgmDD79QDYI9MH4IO370uhj7iHn8wAn0xDSyrRMOBDH2JnoAJaHWfa/8vm9hBSHtIOxHfiuA71R+/g6A4ZCv17F4dds6E9Cf/tKJM67R/iiVCnGlUFzA6blitbbK4wdOuCa/mCKVefQUF5qnZrLS1FCZ4yg+r4yVxqfXrfRVuimD3qjh9uuwF0fC+vDfw8xvVX7+ZwDv0Ry3hIgmAMwDGGPmcd0LEtG9AO4FgJUrV4YcXmvR+fiDas4LxRLuf2oaW3ZM1fgPo1Qq9PVaYI+CWFEQhyYwsjHwJkoRbsZK4WxxoeaeXTW62/jvrRSFrtnTa6Xw3267CgCw86BeImk36Op7tXXX4er3YYnHZNWukglR4Gvwiej7AN7r8tSX7L8wMxOR7tO6lJnzRHQ5gL1EdIiZX3E7kJkfA/AYUFbp+I2vE3Brru6HMpROqZf65xbMJJTdQT985ZSroevrtTD55Ruqv18W4MvYCO029kJrOV9ibN88UGP8TLOpVQ2qIIsZ94Yj5T2CX0kVNyN9bv7C9Hd6rugqsexEKaYdX4PPzB/TPUdE/0JEFzPzW0R0MYC3Na+Rr/z/KhH9HYBBAK4GvxtRgcdGg2RujRO88gseGD9Ud66MlcaDN5ebTKgViphjIUpKlRozQG12uZUmX7fc0sUXtPemiyO3V/QrUW4vqWLHbYJw+96ZHhdXwrp0dgH4DICxyv9/6TyAiPoAzDHzOSK6CMB6AP895Hk7jn1HToYysG43sM6F9JXhqzB06fK6LerE8VOuCWReKLdMs/uZCt1BfqZQ07De1GWYnyngstHdWJHN4Parc6ES89QCKIiCRjdB5GcK1Yx5r+PiLMW0E9bgjwF4iog+B+A4gLsAgIiGAPweM/8OgF8H8CdEtIByjGiMmX8a8rwdh+kNoTOsuht1fDJfY8TV1hgotw5UnJ4rBv4S2YtWBWkLlyQyVqpl5Si6HRWo33kwH6qhudrtOncKXoFYrzib3WXTiVJMO6GC9Mz8C2b+KDO/n5k/xsynKo9PVIw9mPmHzHwVM6+t/P+nUQy80zC9IRio6xOqu1FVvW77in2mUMTI09P40rOHQqlb0kR1W9ckGftsxRXhRZqopserEA2FYgmF+cauq/quBFXQeCng7B3wOlGKaUcybVtEkMDt0kXlOjpeKoDxyby23n5xgVEMWVc+6QHXT6y9uFo2W7fyW3d5H177RUGqfjYB09vP2Uh8w5p+bNtztE7h5oc6xu6OsqN26GFrc7UbMfgtZImVMjL4ZwpFTD14QzW4umXHFLbtOVpXtCnpRrmZ7P7xWxi6dDlmz81rj9n/yin09VotbRyTdLIZS7sYCqugGR7MaSd403LqcUcMfgsIWg9mRTajvXknjp+KPMVcqOf0XFG72nMeJ7SGjJXG1luuDFTHPqiCJqjvv9NIrMFvZfJE0DaJquGD280b1/ongtAM/AqY2YlCQdPpLhs/EmnwW508EeSGU8XDtmhWl40ae1HYCJ2IXwEzO1EpaDrZZeNHTMsnNRevrV/UjE/mkdK0dnI+ak+OikrmlU4RHt08gO2bB9DXG10XLUFoBUG+B52uoGkFiVzhtyp5wiu4mrHSuP3qHPYdOem6dXTzJTaySlfZj6pVnGlnHxPUdlt2DkIzSBECGWsvd4y9cYpKJtSVWOhmEmnwW5U84dVYwa+qntvNu2FNf10DCJNJQE1kzrhFWDmhBI6FZvIrS6zAxtjNHeN04erqVCWBRBr8VkXivRormMrEnMcpbXgQo70im8Gn//ff13Syys8UxK8vxJqgZcV1eIkmOqkOThQk0uC3KhLfjJ2EcxLQ9dVUZKw0Vr0749q20LRcgkwM0WJaQTLpZCOKOfm5ajulDk4UJDJoC5QN5/7R63FsbFPVvx01UQWRxifzWD+2t9oFyt4cxe0cKhis0sn//lV9j1rGhUbrXscI0fDo5gEJIhpyeq6IgYdecG0GFAS/BVan1MGJgkSu8FtFFDsJPwmpyTm8Eohy2Qz2j16P8ck8vrBjSjJGm8y2PUcxd16fvSvUompDAY372b3KmiRNxUMc48Db0NAQT0xMtHsYbUXnslGG2o/xybynwX+00rDCzzUkCO1EKcKWZSwQoVo7x3QBlSSVDhEdZOYht+dkhR9zGpWQ2m/wMOcRhDiglDX2+vqq9v5Dzx3GgzfrSy4A3Z1MFQQx+DGnkcBvkNo9yj20LGM1tb+tIDQLXTtCoZ7EBm07BV1QdsOafu3fBKndo2RpmmRgQegImpUp322IwY85w4M53H51rqYMAwPYeTCvVS8E9cXnZwpS9VGIJUHWIeKW9EdcOjHGyw+vSxgJK2EThDjwWqVgWhD3ZJLklY0iBj+mmNzobiuabtnWpghYiK+ATGgiGeuC48EpO16WsTB3fh7nXdp3rnp3BuvH9nZlWeOoEIMfU0z88G4rmm7Z1oqxTy5LHDEr0+zyH75yqpokmMQ6OSaE8uET0Z1EdJiIFojIVfdZOe7jRHSUiF4motEw50wKfoZblzAi21qh05nxiSfpvhvONYIEcusJG7T9CYDbAPxAdwARpQF8A8CNAD4A4G4i+kDI83Y9XoZblUxwW7kkKWswKFaq7Bu+Z93KQMFAobUsy3jX0AmyqOmWHW9UhDL4zPwzZvabQq8B8DIzv8rM5wF8F8CtYc6bBHR1eB7dPOBZ+2d4MGfc6MRKE+5Zt9K3lk63UFwAHhg/hH1HTkp9oBgze36+Tnxgryc1d34eVqp2ytZN4LLjraUVsswcgNdtv79RecwVIrqXiCaIaOLkyZNNH1xcGR7M4ZHbrkIumwHBe1XvRHXN8qNYYuw7chIb1vTXTS7dyhMHTkgJiZhTLHGNK0YJGPIzBTAqjeMJyGas6nfj0+tWSrcrA3yDtkT0fQDvdXnqS8z8l1EPiJkfA/AYUK6lE/XrdxJ+6eBejdhNSxrnZwqJaoyelPfZ6dib9tz/1HRds51iibF0cQ+mHryh+pizV4SodOrxNfjM/LGQ58gDuMT2+/sqjwkhcKuied+OKfyXZ6bBHMywiREU4saKbMazRShQ75+Xejn+tEKW+SMA7yeiy1A29J8C8FstOG/X4LaS18k23fTJgtBJKFeMnzRZ/PPBCSvL/CQRvQHgwwB2E9GeyuMriOh5AGDmeQCfB7AHwM8APMXMh8MNOzk4/ZdKX9yIHzpjpZH1UUAIZqSoNkFICI8zVuXXyU3888EJtcJn5mcBPOvy+JsAbrL9/jyA58OcK6m4rXJMC6M5KRRLWGKlkLHSgRuhC7UsMFAoLsi1iwi3/g6qbr0bpgIGoRZZosScqHXEM3PFOvWPTuGgJJvqOBPtepoIhGBFrzoZMfbh0a3WdcYekOzZRpHSCjFHVw8/zOu5Bbf8FA5+nbMUC8zYvnnA6FhB8Oo6pWv2npS8kWYgBj/mePXjDIqX39M+Cagg8ZYdUzVBYhNWZDOSzi6EZnwyj9lz9b1/xXcfDjH4McdeLTDMSt+0f6euabrJhEMIXou/mxB/fnCcRc7GJ/N46LnDrv0Z+not31aGgjdi8DsAtfoOUhtcQQC2VxqVOzGVe5qeL+nGLunvv1HsRc687u/eRT1i7EMiQdsOwl5uwQQC8Ol1K7XGPiq5Z6eRtBpCncCbMwVf3b0UQguPrPA7DLXa19UETxNhgdk3tVy3kveSwnUDTreA7jo2QoqAxT3pSOItSWNFNuNr0CXRKjyywu9QdNU0v3rXWhwb2+RZURPQr5ZKzK6vm/LQWbYzmctrXE5y2Qwmv3xDzXUxWTWanmKBEWgHJpQhlO9nL4MuwdpoEIPfoYSppgnoV0vqdZyv+1vXrnQ9/p51K7H1liuNq21Grc837YzlZjDGJ/NIkfeIMlYan7blI6Q9js9VJK/7R6+PhdE3LZPdbnoXle8dt0UMUM5olkSraBCXTgcTpFiUM0C7YU0/dh7M17gflFF0e131+5Mvvo4SM9JEuPvaS/CV4auqx2zZMeUZuCQA112xvKYVnRe9VgpzxQWj9+dFrvJ+7VJT9f693FdpojpDc9nobu3x9gllZONqjDw9jWIbezW6KV3iyOz5Er74vUN45LarcPvVOZfqrUlJ42s+xDH21w4NDfHExES7h9HxuKl7MlYat1+dw74jJyMpJ+uXmKUCyPuOnDTymatUezVRNepnVyol5/v3k1BmrLTrqlLn889mrJpSvQAw+PALHWN044DaFemSrZylFwR3iOggM7u2nJUVfpejqydeKJaw78jJyL5EXslWaSJ89a61GB7Mea6QFVaaqqtl+25j4KEXMFMIZkBVIpgzkOpl7NNEuP1q993ThjX9dSvQjJXG1lvqm8749WbtdghAT6rcacwEr3iKKHSiQXz4XUzQeuJh8HqtEjPuf2oaq0Z3+/rM+3otbLtjrauxPRPQ2CsXVdD3WWLGzoN51zZ7Ow/ma4w9AdrJQRcnWbooHasOY71WyihGYYq9RlMQd8yKbEZ7zUShEw1i8LuYVtYT93stNel4+czdVDR2lniUI85lM3h084BrELuR92lPBlLodgp/Nf1Wtd/q+rG91YliZONqWOl6gzd7vlSVwAJld1AQtVHUFBcYIxtX49jYJnz1rrWhJqNcNlNVie07ctI1hpGpVGytfaw8OevUZ6LQiQYx+F2M18o26i/RyMbVdY2lg5KfKdQZTcX4ZB4FD9+Aij/sH72+TpaqMyJLF3kbNuf1013PmUKxLoFtfDKP4cEcli7Se02VBHbrLVfia3cNNHXVr6qfun1E9h6ySv3ViMLHShNmz81XP0Nd3OVscUGrMAurPhO8ER9+F6OrtOmmPvHqj2vC8GBOWwMlCHajqV4X8I4R2I/zes75/rb4VPR07gxMK5eq3cHwYM7XDaWOVbGUsDWT7GQzFs4UijWf5xMHTrgea5/M7KU81Hj8gtx9vRbeOTtfjbF4/Y2uYqvz/EL0iMHvYtwqbbqpT3QF04BgdceDBCn9DIjdaALeuxUTzbubEfEyrm47oCCVS9V4TSYJdWx1R/LMNIohW1XqVC268bi5vZwVVN2C/8AFL73TfcOo/5zFPdNexKXTxZhuj3VlFoKWOTb1lTuTmXTYjbzutVWWZiPoEn2yGcv1OrldT53rQ41Xdw63Y4HyZxHE2Gczlq/Pe3wyX40xzJ6br4srmBjh4cEcFjTxF4Ze88+AuGdihKzwuxyT7bFu9WyibrG7gpYZlFhwK9Os8/faDaHb6tqrOJwJOlePn3vIa3cE1BpQv/LWdgkqEEw5ZZeD6t6Dc3wzhSKsFKGv18LMXDGQ+66RZjyin48XYvCFQNt8O27GRAcBODa2yfU5nevJbggbMc4mhPUXm4zL7hOvy751LJr9jGquUmTMeZ4gRfKKC4zeRT2Y/PINrn+jQ8U93Nb52YyFc/MLnp+h0H5CGXwiuhPAVgC/DuAaZnZNiyWi1wD8EkAJwLwuC0xoDyYG1w0/2acdr8nD1Ji3M5jnFdQ2Hde2PUfr/NzFBa6JVYxsXK3NWO7rtQKvlnWTRyOB4eHBHCaOn/JMPIt6QhaiJewK/ycAbgPwJwbHbmDmn4c8n9AEGl09m7ofTH3EcTUOUQW1TVxnyqg+7lDTWGnCgzfXZ/Oq8ek+O12560YTrL4yfJVn/+O4foZCmVAGn5l/BgAUQXae0F4aMbg690Nfr4XeRT1ds9LzCmoHeV+mrjMvo+pXBM85GekS3cL0PAg7OYeVAAuN0yofPgN4gYgYwJ8w82O6A4noXgD3AsDKle4leYV4oHMFdVvf0TBBbTtBXGduRtVtp1FfWbJ2MsppJpl2lW+OarckNIavLJOIvk9EP3H5d2uA8/wGM38IwI0Afp+IflN3IDM/xsxDzDzU398f4BRCq0lKVmRU9V3cslgX95gro4MUgVOTUdxKFUQlARYaw3eFz8wfC3sSZs5X/n+biJ4FcA2AH4R9XaH9xNn3HhWNBrV1nLWViJgpFI1XuEF2FGoyMo3PtMrNEtVuSWiMprt0iGgpgBQz/7Ly8w0AHm72eQUhKqKUhIaJB+hiAH7ZrH6TcivdLI1KgIVoCCvL/CSAPwbQD2A3EU0x80YiWgHgW8x8E4D3AHi2EtjtAfAXzPx/Q45b6EA6OVgX1U4mzApXt9PQNbIxvd5RBaVNiHq3JAQjrErnWQDPujz+JoCbKj+/CmBtmPMInY8E68qEWeEG2WkEud6tdLM0K4FOMEMybYWW0MpVZJwJu8INkuRler1b7WZJQtwnrkjxNKElSLCuTKuUTUGud9yUPELzkBW+0BIkWHeBVqxwg5ZBBlDTzyCIXFToHORTFVpCJ6wi7WWE3bpudRKNXG83uWgnXwOhHjH4QkuIe5KWCnK6tSrsRIJeb0mISgbi0hFaRpyDdd0YVA5yvSXGkgxkhS8IEIMXVfkIId6IwRcEiMHrhBiLEB4x+IIAMXhxj7EI0SA+fEFAd2SAhi1dEecYixANYvAFoUInGzwpXSGYIC4dQegCRFYpmCAGXxC6gKSrjAQzxOALQheQdJWRYIYYfEHoApKuMhLMkKCtIHQB3aAyEpqPGHxB6BI6WWUktAZx6QiCICQEMfiCIAgJQQy+IAhCQhCDLwiCkBDE4AuCICQEYuZ2j0ELEZ0EcLzd47BxEYCft3sQAZDxNpdOGm8njRWQ8YbhUmbud3si1gY/bhDRBDMPtXscpsh4m0snjbeTxgrIeJuFuHQEQRASghh8QRCEhCAGPxiPtXsAAZHxNpdOGm8njRWQ8TYF8eELgiAkBFnhC4IgJAQx+IIgCAlBDL4HRHQnER0mogUi0kquiOg1IjpERFNENNHKMTrGYTrejxPRUSJ6mYhGWzlGxziWE9HfENE/Vf7v0xxXqlzbKSLa1eIxel4rIlpMRDsqz79IRKtaOT6X8fiN97NEdNJ2PX+nHeOsjOXbRPQ2Ef1E8zwR0f+ovJcfE9GHWj1Gx3j8xvsRIjpju7ZfbvUYfWFm+af5B+DXAawG8HcAhjyOew3ARZ0wXgBpAK8AuBzAIgDTAD7QpvH+dwCjlZ9HAfyh5rh32jQ+32sF4D8C+Gbl508B2NHGz99kvJ8F8PV2jdExlt8E8CEAP9E8fxOAvwZAANYBeDHm4/0IgL9q93X1+icrfA+Y+WfM3DFdoA3Hew2Al5n5VWY+D+C7AG5t/uhcuRXAdyo/fwfAcJvGocPkWtnfwzMAPkpE1MIx2onTZ+sLM/8AwCmPQ24F8Odc5gCALBFd3JrR1WMw3tgjBj8aGMALRHSQiO5t92B8yAF43fb7G5XH2sF7mPmtys//DOA9muOWENEEER0gouHWDA2A2bWqHsPM8wDOAHh3S0ZXj+lne3vFRfIMEV3SmqE1RJzuVVM+TETTRPTXRHRluwfjJPEdr4jo+wDe6/LUl5j5Lw1f5jeYOU9Evwrgb4joSGU1EDkRjbdleI3X/gszMxHpNMKXVq7v5QD2EtEhZn4l6rEmhOcAPMnM54jod1HenVzf5jF1Cy+hfK++Q0Q3ARgH8P72DqmWxBt8Zv5YBK+Rr/z/NhE9i/LWuikGP4Lx5gHYV3XvqzzWFLzGS0T/QkQXM/Nbla3625rXUNf3VSL6OwCDKPuqm43JtVLHvEFEPQCWAfhFC8bmhu94mdk+tm+hHEeJKy29V8PCzP9q+/l5IvqfRHQRM8elqJq4dMJCREuJ6N+onwHcAMA1ih8TfgTg/UR0GREtQjnQ2FLli41dAD5T+fkzAOp2KETUR0SLKz9fBGA9gJ+2aHwm18r+Hu4AsJcrEbw24Dtehw/8FgA/a+H4grILwG9X1DrrAJyxuQBjBxG9V8VviOgalO1ruyZ/d9odNY7zPwCfRNlveA7AvwDYU3l8BYDnKz9fjrIaYhrAYZRdK7Edb+X3mwD8I8qr5HaO990A/hbAPwH4PoDllceHAHyr8vN1AA5Vru8hAJ9r8RjrrhWAhwHcUvl5CYCnAbwM4B8AXN7me9ZvvI9U7tNpAPsArGnjWJ8E8BaAYuW+/RyA3wPwe5XnCcA3Ku/lEDyUcjEZ7+dt1/YAgOvaOV63f1JaQRAEISGIS0cQBCEhiMEXBEFICGLwBUEQEoIYfEEQhIQgBl8QBCEhiMEXBEFICGLwBUEQEsL/B0xahgv6SIV6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(Y_data[:,0],Y_data[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "matched-journal",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "Dictionary length 33\n",
      "[]\n",
      "Pruned length 17\n",
      "1 iter, w_dict length = 17, Train -logL = 1.686\n",
      "Dictionary length 50\n",
      "[35 36 45 46 48]\n",
      "Pruned length 30\n",
      "2 iter, w_dict length = 30, Train -logL = 1.571\n",
      "Dictionary length 64\n",
      "[49 53 54]\n",
      "Pruned length 36\n",
      "3 iter, w_dict length = 36, Train -logL = 1.365\n",
      "Dictionary length 64\n",
      "[52 53]\n",
      "Pruned length 32\n",
      "4 iter, w_dict length = 32, Train -logL = 1.133\n",
      "Dictionary length 48\n",
      "[]\n",
      "Pruned length 38\n",
      "5 iter, w_dict length = 38, Train -logL = 1.026\n",
      "Dictionary length 55\n",
      "[]\n",
      "Pruned length 42\n",
      "6 iter, w_dict length = 42, Train -logL = 0.910\n",
      "Dictionary length 54\n",
      "[50 51]\n",
      "Pruned length 41\n",
      "7 iter, w_dict length = 41, Train -logL = 0.851\n",
      "Dictionary length 47\n",
      "[]\n",
      "Pruned length 44\n",
      "8 iter, w_dict length = 44, Train -logL = 0.814\n",
      "Dictionary length 48\n",
      "[47]\n",
      "Pruned length 45\n",
      "9 iter, w_dict length = 45, Train -logL = 0.808\n",
      "Dictionary length 46\n",
      "[]\n",
      "Pruned length 44\n",
      "10 iter, w_dict length = 44, Train -logL = 0.801\n",
      "Dictionary length 44\n",
      "[]\n",
      "Pruned length 44\n",
      "11 iter, w_dict length = 44, Train -logL = 0.801\n",
      "Dictionary length 44\n",
      "[]\n",
      "Pruned length 44\n",
      "12 iter, w_dict length = 44, Train -logL = 0.801\n",
      "[]\n",
      "Final length of w_dict = 44\n",
      "Done, w_dict length = 44\n"
     ]
    }
   ],
   "source": [
    "entropy = md.get_entropy(Y)\n",
    "w_thr = 1e-5*np.exp(-entropy)\n",
    "eps = 0.0\n",
    "p_d = 0.5 \n",
    "p_ins = 0.2 #Unused\n",
    "\n",
    "Sigma = 7 \n",
    "Jthr = 0.3\n",
    "params = np.array([eps,p_d,p_ins, w_thr, Jthr, Sigma], dtype=float)\n",
    "P_w, w_dict = md.solve_dictionary(Y,lengths_Y,params,model_true,15) #Solves for the dictionary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "absolute-beginning",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Dictionary\n",
      "[3 2 6 4 0 4]\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "True Dictionary\n",
      "[3 2 6 4 0 4]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n\")\n",
    "#Output\n",
    "print(\"Dictionary\")\n",
    "print(w_dict[22])\n",
    "\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"\\n\")\n",
    "#Output\n",
    "print(\"True Dictionary\")\n",
    "print(w_dict_true[13])\n",
    "print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "hidden-marking",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[13 30 36 32  0 34 29 11 44 15 18  1 16 28 35  3 42 14 38 31  8 26 10 41\n",
      "  5 20 22 39 37 45  6 33 49  9  2 40 46 17  4 21 19 43 12 25 47 23  7 48\n",
      " 24 27 54 50 55 51 53 52 56]\n",
      "0.06784963962068925\n"
     ]
    }
   ],
   "source": [
    "idx = np.argsort(P_w_true)[::-1]\n",
    "print(idx)\n",
    "print(P_w_true[13])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "existing-sugar",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[22 24 16 19 20 17  7 33 13  8 29 30 15 27 11  2  9  1 31 38 14 18 32 34\n",
      " 25 40 26 37 12  3 41  5 36 28 35  4 21  6 10  0 39 43 42 23]\n",
      "[0.0633771  0.05454514 0.05163899 0.05162916 0.05051882 0.04265023\n",
      " 0.03948246 0.0377097  0.03746888 0.03735276 0.03573001 0.03572898\n",
      " 0.03474793 0.03143759 0.02676955 0.0266027  0.02495036 0.02474071\n",
      " 0.02084391 0.02084376 0.01947906 0.01687981 0.0158901  0.0156344\n",
      " 0.01498817 0.0132585  0.01246848 0.01089545 0.01079265 0.01022463\n",
      " 0.00991927 0.00968817 0.00909676 0.00890605 0.00850647 0.00790321\n",
      " 0.00790096 0.00786443 0.00773506 0.00771853 0.00759289 0.00694526\n",
      " 0.00592468 0.00501831]\n"
     ]
    }
   ],
   "source": [
    "idx = np.argsort(P_w)[::-1]\n",
    "print(idx)\n",
    "print(P_w[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "quiet-commerce",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "only integer scalar arrays can be converted to a scalar index",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-aa413daa3dd6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0midx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: only integer scalar arrays can be converted to a scalar index"
     ]
    }
   ],
   "source": [
    "print(w_dict[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "moved-panel",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}