ZMotif.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "7e3c9588-9618-45ab-8351-4d5769dc21a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from keras import initializers\n",
    "from keras.layers import Input, Lambda, Conv1D, maximum, GlobalMaxPooling1D, Dense, GaussianNoise, MaxPooling1D, Flatten\n",
    "from keras import regularizers\n",
    "from keras.constraints import non_neg\n",
    "from keras.models import Model\n",
    "from keras import backend as K\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import logomaker\n",
    "\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "01d864b3-e883-42a0-8ec4-6d6189ca91cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_model(num_kernels=32,\n",
    "                    kernel_width=24,\n",
    "                    seq_len=None,\n",
    "                    dropout_prop=0.0,\n",
    "                    kernel_initializer=initializers.RandomNormal(stddev=0.001, seed=12),\n",
    "                    optimizer='adam',\n",
    "                    activation='linear',\n",
    "                    num_classes=1,\n",
    "                    l1_reg=0.0,\n",
    "                    l2_reg= 0.0,\n",
    "                    gaussian_noise = 0.0,\n",
    "                    spatial_dropout = 0.0,\n",
    "                    rc = True,\n",
    "                    padding=\"same\",\n",
    "                    conv_name=\"shared_conv\"):\n",
    "    if rc:\n",
    "        seq_input = Input(shape=(seq_len,4))\n",
    "        rc_op = Lambda(lambda x: K.reverse(x,axes=(1,2)))\n",
    "        seq_rc = rc_op(seq_input)\n",
    "        if gaussian_noise > 0.0:\n",
    "            noisy_seq = GaussianNoise(gaussian_noise)(seq_input)\n",
    "            noisy_seq_rc = rc_op(noisy_seq)\n",
    "        \n",
    "        shared_conv = Conv1D(num_kernels, kernel_width,\n",
    "                             strides=1, padding=padding, \n",
    "                             activation=activation,\n",
    "                             use_bias=False,\n",
    "                             kernel_initializer=kernel_initializer,\n",
    "                             kernel_regularizer=regularizers.l1_l2(l1=l1_reg,\n",
    "                                                                   l2=l2_reg),\n",
    "                             bias_initializer='zeros',\n",
    "                             name=conv_name)\n",
    "\n",
    "        if gaussian_noise > 0:\n",
    "            conv_for = shared_conv(noisy_seq)\n",
    "            conv_rc = shared_conv(noisy_seq_rc)\n",
    "        else:\n",
    "            conv_for = shared_conv(seq_input)\n",
    "            conv_rc = shared_conv(seq_rc)\n",
    "            \n",
    "\n",
    "        merged = maximum([conv_for, conv_rc])\n",
    "        pooled = GlobalMaxPooling1D()(merged)\n",
    "        if dropout_prop > 0.0:\n",
    "            dropout = Dropout(dropout_prop)(pooled)\n",
    "            output = Dense(1, activation='sigmoid',\n",
    "                       use_bias=False,\n",
    "                       kernel_initializer=initializers.RandomUniform(minval=0.0, maxval=0.001, seed=12), \n",
    "                       kernel_constraint=non_neg(), \n",
    "                       bias_initializer='zeros',\n",
    "                       name=\"dense_1\")(dropout)\n",
    "        else:\n",
    "            output = Dense(1, activation='sigmoid',\n",
    "                           use_bias=False,\n",
    "                           kernel_initializer=initializers.RandomUniform(minval=0.0, maxval=0.001, seed=12), \n",
    "                           kernel_constraint=non_neg(), \n",
    "                           bias_initializer='zeros',\n",
    "                           name=\"dense_1\")(pooled)\n",
    "        model = Model(inputs=seq_input, outputs=output)\n",
    "        model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n",
    "        model.summary()\n",
    "        return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9845218e-5780-4e51-a222-187fef5bf3b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from  keras.callbacks import Callback\n",
    "class SWA(Callback):\n",
    "\n",
    "    def __init__(self, epochs_to_train, prop = 0.2, interval = 1):\n",
    "        super(SWA, self).__init__()\n",
    "        self.epochs_to_train = epochs_to_train\n",
    "        self.prop = prop\n",
    "        self.interval = interval\n",
    "        self.n_models = 0\n",
    "        self.epoch = 0\n",
    "        \n",
    "    def on_train_begin(self, logs=None):\n",
    "        self.nb_epoch = self.params['epochs']\n",
    "        self.models_weights = []\n",
    "    def on_epoch_end(self, epoch, logs=None):\n",
    "        self.epoch += 1\n",
    "        if epoch % self.interval == 0:\n",
    "            self.models_weights.append(self.model.get_weights())\n",
    "            self.n_models += 1\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "    def on_train_end(self, logs=None):\n",
    "        if self.epoch > 10:\n",
    "            num_models_to_average = int(np.ceil(self.prop * self.epoch))\n",
    "            if num_models_to_average < 1:\n",
    "                print(\"Taking last model\")\n",
    "                num_models_to_average = 1\n",
    "#         print(len(self.models_weights))\n",
    "#         print(len(self.models_weights[0]))\n",
    "#         print(self.models_weights[0][0].shape)\n",
    "#         print(self.models_weights[0][1].shape)\n",
    "            avg_conv_weights = np.mean([weights[0] for weights in self.models_weights[-num_models_to_average:]], axis=0)\n",
    "            avg_dense_weights = np.mean([weights[1] for weights in self.models_weights[-num_models_to_average:]], axis=0)\n",
    "        #print(len(avg_conv_weights))\n",
    "#         print(avg_conv_weights.shape)\n",
    "#         print(avg_dense_weights.shape)\n",
    "            self.model.set_weights([avg_conv_weights, avg_dense_weights])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dfae91e6-7813-4be1-bdc8-749b0b6ccb56",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def construct_scan_model(conv_weights):\n",
    "    kernel_width = conv_weights.shape[0]\n",
    "    num_kernels = conv_weights.shape[2]\n",
    "    seq = Input(shape=(None,4))\n",
    "    conv = Conv1D(num_kernels, kernel_width, \n",
    "                  name = 'scan_conv',\n",
    "                  strides=1, \n",
    "                  padding='valid', \n",
    "                  activation='linear', \n",
    "                  use_bias=False, \n",
    "                  kernel_initializer='zeros', \n",
    "                  bias_initializer='zeros',\n",
    "                  trainable=False)\n",
    "    \n",
    "    conv_seq = conv(seq)\n",
    "    \n",
    "    \n",
    "    model = Model(inputs=seq, outputs=conv_seq)\n",
    "    model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n",
    "    model.get_layer('scan_conv').set_weights([conv_weights])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b1b185c7-f2b7-4a82-889d-3a04d636aaaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_score_model(conv_weights):\n",
    "    kernel_width = conv_weights.shape[0]\n",
    "    num_kernels = conv_weights.shape[2]\n",
    "    seq = Input(shape=(None,4))\n",
    "    rc_op = Lambda(lambda x: K.reverse(x,axes=(1,2)))\n",
    "    seq_rc = rc_op(seq)\n",
    "    \n",
    "    conv = Conv1D(num_kernels, kernel_width, \n",
    "                  name = 'score_conv',\n",
    "                  strides=1, \n",
    "                  padding='valid', \n",
    "                  activation='linear', \n",
    "                  use_bias=False, \n",
    "                  kernel_initializer='zeros', \n",
    "                  bias_initializer='zeros',\n",
    "                  trainable=False)\n",
    "    \n",
    "    conv_for = conv(seq)\n",
    "    conv_rc = conv(seq_rc)\n",
    "    \n",
    "    merged = maximum([conv_for, conv_rc])\n",
    "    pooled = GlobalMaxPooling1D()(merged)\n",
    "    # flat = flatten()(pooled)\n",
    "    model = Model(inputs=seq, outputs=pooled)\n",
    "    model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n",
    "    model.get_layer(\"score_conv\").set_weights([conv_weights])\n",
    "    print(model.summary())\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "66527a78-cdf6-4cdb-b4d8-9d866d2cc495",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "DNA_SEQ_DICT = {\n",
    "    'A' : [1, 0, 0, 0],\n",
    "    'C' : [0, 1, 0, 0],\n",
    "    'G' : [0, 0, 1, 0],\n",
    "    'T' : [0, 0, 0, 1],\n",
    "}\n",
    "\n",
    "def encode_sequence(seq, N = [0, 0, 0, 0], seq_dict = None, useN = None):\n",
    "    if seq_dict is None:\n",
    "        seq_dict = DNA_SEQ_DICT\n",
    "    if useN == 'uniform':\n",
    "        N = [(1/len(seq_dict)) for _ in seq_dict]\n",
    "    elif useN == 'zeros':\n",
    "        N = [0 for _ in seq_dict]\n",
    "    d = { **seq_dict, 'N' : N }\n",
    "    return np.array([d[nuc] for nuc in list(seq)]).astype('float32')\n",
    " \n",
    "def decode_sequence(encoded_seq, seq_dict = None):\n",
    "    if seq_dict is None:\n",
    "        seq_dict = DNA_SEQ_DICT\n",
    "    seq_list = encoded_seq.astype('int').tolist()\n",
    "    def decode_base(encoded_base):\n",
    "        for letter,onehot in seq_dict.items():\n",
    "            if np.array_equal(encoded_base, onehot):\n",
    "                return letter\n",
    "        return \"N\"\n",
    "    return \"\".join(decode_base(b) for b in encoded_seq.astype('int'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5e096620-f3d1-4ad8-b365-30e91cf872f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# altschulEriksonDinuclShuffle.py\n",
    "# P. Clote, Oct 2003\n",
    "\n",
    "def computeCountAndLists(s):\n",
    "\n",
    "    #Initialize lists and mono- and dinucleotide dictionaries\n",
    "    List = {} #List is a dictionary of lists\n",
    "    List['A'] = []; List['C'] = [];\n",
    "    List['G'] = []; List['T'] = [];\n",
    "    # FIXME: is this ok?\n",
    "    List['N'] = []\n",
    "    nuclList   = [\"A\",\"C\",\"G\",\"T\",\"N\"]\n",
    "    s       = s.upper()\n",
    "    #s       = s.replace(\"U\",\"T\")\n",
    "    nuclCnt    = {}  #empty dictionary\n",
    "    dinuclCnt  = {}  #empty dictionary\n",
    "    for x in nuclList:\n",
    "        nuclCnt[x]=0\n",
    "        dinuclCnt[x]={}\n",
    "        for y in nuclList:\n",
    "            dinuclCnt[x][y]=0\n",
    "\n",
    "    #Compute count and lists\n",
    "    nuclCnt[s[0]] = 1\n",
    "    nuclTotal     = 1\n",
    "    dinuclTotal   = 0\n",
    "    for i in range(len(s)-1):\n",
    "        x = s[i]; y = s[i+1]\n",
    "        List[x].append( y )\n",
    "        nuclCnt[y] += 1; nuclTotal  += 1\n",
    "        dinuclCnt[x][y] += 1; dinuclTotal += 1\n",
    "    assert (nuclTotal==len(s))\n",
    "    assert (dinuclTotal==len(s)-1)\n",
    "    return nuclCnt,dinuclCnt,List\n",
    "\n",
    "\n",
    "def chooseEdge(x,dinuclCnt):\n",
    "    z = random.random()\n",
    "    denom=dinuclCnt[x]['A']+dinuclCnt[x]['C']+dinuclCnt[x]['G']+dinuclCnt[x]['T']+dinuclCnt[x]['N']\n",
    "    numerator = dinuclCnt[x]['A']\n",
    "    if z < float(numerator)/float(denom):\n",
    "        dinuclCnt[x]['A'] -= 1\n",
    "        return 'A'\n",
    "    numerator += dinuclCnt[x]['C']\n",
    "    if z < float(numerator)/float(denom):\n",
    "        dinuclCnt[x]['C'] -= 1\n",
    "        return 'C'\n",
    "    numerator += dinuclCnt[x]['G']\n",
    "    if z < float(numerator)/float(denom):\n",
    "        dinuclCnt[x]['G'] -= 1\n",
    "        return 'G'\n",
    "    numerator += dinuclCnt[x]['T']\n",
    "    if z < float(numerator)/float(denom):\n",
    "        dinuclCnt[x]['T'] -= 1\n",
    "        return 'T'\n",
    "    dinuclCnt[x]['N'] -= 1\n",
    "    return 'N'\n",
    "\n",
    "def connectedToLast(edgeList,nuclList,lastCh):\n",
    "    D = {}\n",
    "    for x in nuclList: D[x]=0\n",
    "    for edge in edgeList:\n",
    "        a = edge[0]; b = edge[1]\n",
    "        if b==lastCh: D[a]=1\n",
    "    for i in range(3):\n",
    "        for edge in edgeList:\n",
    "            a = edge[0]; b = edge[1]\n",
    "            if D[b]==1: D[a]=1\n",
    "    ok = 0\n",
    "    for x in nuclList:\n",
    "        if x!=lastCh and D[x]==0: return 0\n",
    "    return 1\n",
    "\n",
    "def eulerian(s):\n",
    "    nuclCnt,dinuclCnt,List = computeCountAndLists(s)\n",
    "    #compute nucleotides appearing in s\n",
    "    nuclList = []\n",
    "    for x in [\"A\",\"C\",\"G\",\"T\",\"N\"]:\n",
    "        if x in s: nuclList.append(x)\n",
    "    #create dinucleotide shuffle L\n",
    "    firstCh = s[0]  #start with first letter of s\n",
    "    lastCh  = s[-1]\n",
    "    edgeList = []\n",
    "    for x in nuclList:\n",
    "        if x!= lastCh: edgeList.append( [x,chooseEdge(x,dinuclCnt)] )\n",
    "    ok = connectedToLast(edgeList,nuclList,lastCh)\n",
    "    return ok,edgeList,nuclList,lastCh\n",
    "\n",
    "\n",
    "def shuffleEdgeList(L):\n",
    "    n = len(L); barrier = n\n",
    "    for i in range(n-1):\n",
    "        z = int(random.random() * barrier)\n",
    "        tmp = L[z]\n",
    "        L[z]= L[barrier-1]\n",
    "        L[barrier-1] = tmp\n",
    "        barrier -= 1\n",
    "    return L\n",
    "\n",
    "def dinuclShuffle(s):\n",
    "    ok = 0\n",
    "    while not ok:\n",
    "        ok,edgeList,nuclList,lastCh = eulerian(s)\n",
    "    nuclCnt,dinuclCnt,List = computeCountAndLists(s)\n",
    "\n",
    "    #remove last edges from each vertex list, shuffle, then add back\n",
    "    #the removed edges at end of vertex lists.\n",
    "    for [x,y] in edgeList: List[x].remove(y)\n",
    "    for x in nuclList: shuffleEdgeList(List[x])\n",
    "    for [x,y] in edgeList: List[x].append(y)\n",
    "\n",
    "    #construct the eulerian path\n",
    "    L = [s[0]]; prevCh = s[0]\n",
    "    for i in range(len(s)-2):\n",
    "        ch = List[prevCh][0]\n",
    "        L.append( ch )\n",
    "        del List[prevCh][0]\n",
    "        prevCh = ch\n",
    "    L.append(s[-1])\n",
    "    #t = string.join(L,\"\")\n",
    "    t = \"\".join(L)\n",
    "    return t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2bdcec2a-94c7-4934-9a4e-c38a38f151f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_information_content(x):\n",
    "    ic = x * np.log2((x + .001) / .25)\n",
    "    if ic > 0:\n",
    "        return(ic)\n",
    "    else:\n",
    "        return(0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4bc7980e-3d08-428c-acd6-c1822835e8f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from  keras.callbacks import Callback\n",
    "class SGDRScheduler(Callback):\n",
    "    '''Cosine annealing learning rate scheduler with periodic restarts.\n",
    "\n",
    "    # Usage\n",
    "        ```python\n",
    "            schedule = SGDRScheduler(min_lr=1e-5,\n",
    "                                     max_lr=1e-2,\n",
    "                                     steps_per_epoch=np.ceil(epoch_size/batch_size),\n",
    "                                     lr_decay=0.9,\n",
    "                                     cycle_length=5,\n",
    "                                     mult_factor=1.5)\n",
    "            model.fit(X_train, Y_train, epochs=100, callbacks=[schedule])\n",
    "        ```\n",
    "\n",
    "    # Arguments\n",
    "        min_lr: The lower bound of the learning rate range for the experiment.\n",
    "        max_lr: The upper bound of the learning rate range for the experiment.\n",
    "        steps_per_epoch: Number of mini-batches in the dataset. Calculated as `np.ceil(epoch_size/batch_size)`. \n",
    "        lr_decay: Reduce the max_lr after the completion of each cycle.\n",
    "                  Ex. To reduce the max_lr by 20% after each cycle, set this value to 0.8.\n",
    "        cycle_length: Initial number of epochs in a cycle.\n",
    "        mult_factor: Scale epochs_to_restart after each full cycle completion.\n",
    "\n",
    "    # References\n",
    "        Blog post: jeremyjordan.me/nn-learning-rate\n",
    "        Original paper: http://arxiv.org/abs/1608.03983\n",
    "    '''\n",
    "    def __init__(self,\n",
    "                 min_lr,\n",
    "                 max_lr,\n",
    "                 steps_per_epoch,\n",
    "                 lr_decay=1,\n",
    "                 cycle_length=10,\n",
    "                 mult_factor=2,\n",
    "                 shape=\"cosine\"):\n",
    "\n",
    "        self.min_lr = min_lr\n",
    "        self.max_lr = max_lr\n",
    "        self.lr_decay = lr_decay\n",
    "\n",
    "        self.batch_since_restart = 0\n",
    "        self.next_restart = cycle_length\n",
    "\n",
    "        self.steps_per_epoch = steps_per_epoch\n",
    "\n",
    "        self.cycle_length = cycle_length\n",
    "        self.mult_factor = mult_factor\n",
    "        \n",
    "        self.shape = shape\n",
    "        self.history = {}\n",
    "        self.learning_rates = []\n",
    "\n",
    "    def clr(self):\n",
    "        '''Calculate the learning rate.'''\n",
    "        fraction_to_restart = self.batch_since_restart / (self.steps_per_epoch * self.cycle_length)\n",
    "        #print(fraction_to_restart)\n",
    "        if self.shape == \"cosine\":\n",
    "            lr = self.min_lr + 0.5 * (self.max_lr - self.min_lr) * (1 + np.cos(fraction_to_restart * np.pi))\n",
    "        else:\n",
    "            if fraction_to_restart < 0.5:\n",
    "                lr = fraction_to_restart * (self.max_lr - self.min_lr) / 0.5 + self.min_lr\n",
    "            else:\n",
    "                lr = (1 - fraction_to_restart) * (self.max_lr - self.min_lr) / 0.5 + self.min_lr\n",
    "        self.learning_rates.append(lr)\n",
    "        return lr\n",
    "\n",
    "    def on_train_begin(self, logs={}):\n",
    "        '''Initialize the learning rate to the minimum value at the start of training.'''\n",
    "        logs = logs or {}\n",
    "        if self.shape == \"cosine\":\n",
    "            K.set_value(self.model.optimizer.lr, self.max_lr)\n",
    "        else:\n",
    "            K.set_value(self.model.optimizer.lr, self.min_lr)\n",
    "\n",
    "    def on_batch_end(self, batch, logs={}):\n",
    "        '''Record previous batch statistics and update the learning rate.'''\n",
    "        logs = logs or {}\n",
    "        self.history.setdefault('lr', []).append(K.get_value(self.model.optimizer.lr))\n",
    "        for k, v in logs.items():\n",
    "            self.history.setdefault(k, []).append(v)\n",
    "\n",
    "        self.batch_since_restart += 1\n",
    "        K.set_value(self.model.optimizer.lr, self.clr())\n",
    "\n",
    "    def on_epoch_end(self, epoch, logs={}):\n",
    "        '''Check for end of current cycle, apply restarts when necessary.'''\n",
    "        if epoch + 1 == self.next_restart:\n",
    "            self.batch_since_restart = 0\n",
    "            self.cycle_length = np.ceil(self.cycle_length * self.mult_factor)\n",
    "            self.next_restart += self.cycle_length\n",
    "            self.max_lr *= self.lr_decay\n",
    "            self.best_weights = self.model.get_weights()\n",
    "\n",
    "    def on_train_end(self, logs={}):\n",
    "        '''Set weights to the values from the end of the most recent cycle for best performance.'''\n",
    "        self.model.set_weights(self.best_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ed5cf3db-7fd4-4e94-9c3c-c498c316b3a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyfaidx import Fasta\n",
    "posFasta = Fasta(\"CTCF.pos.fasta\", \n",
    "                 as_raw=True,\n",
    "                 sequence_always_upper=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1d812bca-ef9a-494b-852a-3493f4d65cc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "posSeqs = [seq[:] for seq in posFasta]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "ccced61a-cd93-4bf3-9b11-7a7d02b76c85",
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.utils import Sequence\n",
    "import random\n",
    "class dataGen(Sequence):\n",
    "    def __init__(self, posSeqs, \n",
    "                 negSeqs=None,\n",
    "                 batchSize = 32,\n",
    "                 seqsPerEpoch=10000,\n",
    "                 padBy = 24):\n",
    "        \n",
    "        self.posSeqs = posSeqs\n",
    "        self.nPos = len(self.posSeqs)\n",
    "        \n",
    "        if negSeqs is not None:\n",
    "            self.negSeqs = negSeqs\n",
    "            self.nNeg = len(self.negSeqs)\n",
    "        else:\n",
    "            self.negSeqs = None\n",
    "            \n",
    "        \n",
    "        if self.negSeqs is not None:\n",
    "            self.L = np.max([len(x) for x in self.posSeqs + self.negSeqs])\n",
    "        else:\n",
    "            self.L = np.max([len(x) for x in self.posSeqs])\n",
    "        print(\"Maximum sequence length = {}\".format(self.L))\n",
    "        \n",
    "        self.batchSize = batchSize\n",
    "        self.b2 = self.batchSize // 2\n",
    "        self.seqsPerEpoch = seqsPerEpoch\n",
    "        self.padBy = padBy\n",
    "        \n",
    "        self.labels = np.array([1 for i in range(self.b2)] + [0 for i in range(self.b2)])\n",
    "        \n",
    "    def __len__(self):\n",
    "        return(int(np.floor(self.seqsPerEpoch / self.batchSize)))\n",
    "    \n",
    "    def on_epoch_end(self):\n",
    "        pass\n",
    "    \n",
    "    def __getitem__(self, index):\n",
    "        posSample = random.sample(self.posSeqs, self.b2)\n",
    "        \n",
    "        if self.negSeqs is not None:\n",
    "            negSample = random.sample(self.negSeqs, self.b2)\n",
    "        else:\n",
    "            negSample = [dinuclShuffle(x) for x in posSample]\n",
    "        \n",
    "        X = 0.25 * np.ones((self.batchSize, 2*self.padBy + self.L, 4))\n",
    "        for i,seq in enumerate(posSample + negSample):\n",
    "            l = len(seq)\n",
    "            start = self.padBy + (self.L - l) // 2\n",
    "            stop = start + l\n",
    "            X[i,start:stop,:] = encode_sequence(seq)\n",
    "            \n",
    "        return(X, self.labels)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "b5060042-fa23-40d1-a7ea-3580587dffaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum sequence length = 200\n",
      "Maximum sequence length = 200\n"
     ]
    }
   ],
   "source": [
    "padBy = 24\n",
    "w =24\n",
    "num_kernels = 1\n",
    "holdOutP = 0.1\n",
    "n = len(posSeqs)\n",
    "random.shuffle(posSeqs)\n",
    "trainGen = dataGen(posSeqs[:int((1 - holdOutP)*n)])\n",
    "testGen = dataGen(posSeqs[int((1 - holdOutP)*n):],\n",
    "                  seqsPerEpoch=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "f01d77d9-8473-4b96-8168-43c3ee255989",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_6 (InputLayer)            (None, None, 4)      0                                            \n",
      "__________________________________________________________________________________________________\n",
      "lambda_4 (Lambda)               (None, None, 4)      0           input_6[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "shared_conv (Conv1D)            (None, None, 1)      96          input_6[0][0]                    \n",
      "                                                                 lambda_4[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "maximum_4 (Maximum)             (None, None, 1)      0           shared_conv[0][0]                \n",
      "                                                                 shared_conv[1][0]                \n",
      "__________________________________________________________________________________________________\n",
      "global_max_pooling1d_4 (GlobalM (None, 1)            0           maximum_4[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dense_1 (Dense)                 (None, 1)            1           global_max_pooling1d_4[0][0]     \n",
      "==================================================================================================\n",
      "Total params: 97\n",
      "Trainable params: 97\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "model = construct_model(num_kernels=num_kernels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "d37e992f-f6ef-4775-b70d-cc54baca5b72",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.954992586021436\n"
     ]
    }
   ],
   "source": [
    "epochs = 100\n",
    "min_lr = .001\n",
    "max_lr = .1\n",
    "lr_decay = (min_lr / max_lr) ** (1 / epochs)\n",
    "schedule = SGDRScheduler(min_lr=min_lr,\n",
    "                             max_lr=max_lr,\n",
    "                             steps_per_epoch=trainGen.__len__(),\n",
    "                             lr_decay=lr_decay,\n",
    "                             cycle_length=1,\n",
    "                             mult_factor=1.0, \n",
    "                             shape=\"triangular\")\n",
    "\n",
    "swa = SWA(epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "25809158-195e-44a2-bde0-352b110891b5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.4793 - acc: 0.7287 - val_loss: 0.2898 - val_acc: 0.8880\n",
      "Epoch 2/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.3079 - acc: 0.8771 - val_loss: 0.2815 - val_acc: 0.8871\n",
      "Epoch 3/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2903 - acc: 0.8842 - val_loss: 0.2846 - val_acc: 0.8859\n",
      "Epoch 4/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2964 - acc: 0.8842 - val_loss: 0.2800 - val_acc: 0.8885\n",
      "Epoch 5/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2914 - acc: 0.8842 - val_loss: 0.2880 - val_acc: 0.8861\n",
      "Epoch 6/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2808 - acc: 0.8868 - val_loss: 0.2828 - val_acc: 0.8896\n",
      "Epoch 7/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2834 - acc: 0.8877 - val_loss: 0.2895 - val_acc: 0.8868\n",
      "Epoch 8/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2873 - acc: 0.8858 - val_loss: 0.2938 - val_acc: 0.8866\n",
      "Epoch 9/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2893 - acc: 0.8854 - val_loss: 0.2756 - val_acc: 0.8915\n",
      "Epoch 10/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2835 - acc: 0.8862 - val_loss: 0.2865 - val_acc: 0.8902\n",
      "Epoch 11/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2815 - acc: 0.8888 - val_loss: 0.2827 - val_acc: 0.8900\n",
      "Epoch 12/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2819 - acc: 0.8877 - val_loss: 0.2836 - val_acc: 0.8865\n",
      "Epoch 13/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2785 - acc: 0.8894 - val_loss: 0.2771 - val_acc: 0.8904\n",
      "Epoch 14/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2846 - acc: 0.8875 - val_loss: 0.2705 - val_acc: 0.8935\n",
      "Epoch 15/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2736 - acc: 0.8931 - val_loss: 0.2848 - val_acc: 0.8864\n",
      "Epoch 16/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2845 - acc: 0.8863 - val_loss: 0.2692 - val_acc: 0.8943\n",
      "Epoch 17/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2894 - acc: 0.8855 - val_loss: 0.2806 - val_acc: 0.8898\n",
      "Epoch 18/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2763 - acc: 0.8884 - val_loss: 0.2805 - val_acc: 0.8908\n",
      "Epoch 19/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2805 - acc: 0.8915 - val_loss: 0.2654 - val_acc: 0.8962\n",
      "Epoch 20/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2746 - acc: 0.8925 - val_loss: 0.2830 - val_acc: 0.8874\n",
      "Epoch 21/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2773 - acc: 0.8898 - val_loss: 0.2854 - val_acc: 0.8860\n",
      "Epoch 22/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2702 - acc: 0.8922 - val_loss: 0.2792 - val_acc: 0.8904\n",
      "Epoch 23/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2785 - acc: 0.8889 - val_loss: 0.2834 - val_acc: 0.8907\n",
      "Epoch 24/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2779 - acc: 0.8909 - val_loss: 0.2697 - val_acc: 0.8966\n",
      "Epoch 25/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2797 - acc: 0.8902 - val_loss: 0.2724 - val_acc: 0.8966\n",
      "Epoch 26/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2644 - acc: 0.8930 - val_loss: 0.2736 - val_acc: 0.8894\n",
      "Epoch 27/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2781 - acc: 0.8911 - val_loss: 0.2734 - val_acc: 0.8928\n",
      "Epoch 28/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2719 - acc: 0.8912 - val_loss: 0.2765 - val_acc: 0.8923\n",
      "Epoch 29/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2733 - acc: 0.8915 - val_loss: 0.2867 - val_acc: 0.8916\n",
      "Epoch 30/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2855 - acc: 0.8873 - val_loss: 0.2819 - val_acc: 0.8922\n",
      "Epoch 31/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2742 - acc: 0.8908 - val_loss: 0.2679 - val_acc: 0.8946\n",
      "Epoch 32/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2690 - acc: 0.8950 - val_loss: 0.2738 - val_acc: 0.8962\n",
      "Epoch 33/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2656 - acc: 0.8974 - val_loss: 0.2693 - val_acc: 0.8952\n",
      "Epoch 34/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2731 - acc: 0.8893 - val_loss: 0.2780 - val_acc: 0.8924\n",
      "Epoch 35/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2760 - acc: 0.8954 - val_loss: 0.2764 - val_acc: 0.8909\n",
      "Epoch 36/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2699 - acc: 0.8931 - val_loss: 0.2711 - val_acc: 0.8963\n",
      "Epoch 37/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2767 - acc: 0.8921 - val_loss: 0.2789 - val_acc: 0.8915\n",
      "Epoch 38/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2717 - acc: 0.8930 - val_loss: 0.2789 - val_acc: 0.8938\n",
      "Epoch 39/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2679 - acc: 0.8953 - val_loss: 0.2795 - val_acc: 0.8931\n",
      "Epoch 40/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2673 - acc: 0.8931 - val_loss: 0.2803 - val_acc: 0.8915\n",
      "Epoch 41/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2786 - acc: 0.8883 - val_loss: 0.2739 - val_acc: 0.8945\n",
      "Epoch 42/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2769 - acc: 0.8907 - val_loss: 0.2910 - val_acc: 0.8860\n",
      "Epoch 43/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2753 - acc: 0.8922 - val_loss: 0.2847 - val_acc: 0.8889\n",
      "Epoch 44/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2721 - acc: 0.8911 - val_loss: 0.2736 - val_acc: 0.8938\n",
      "Epoch 45/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2756 - acc: 0.8889 - val_loss: 0.2760 - val_acc: 0.8920\n",
      "Epoch 46/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2666 - acc: 0.8951 - val_loss: 0.2849 - val_acc: 0.8871\n",
      "Epoch 47/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2755 - acc: 0.8901 - val_loss: 0.2746 - val_acc: 0.8963\n",
      "Epoch 48/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2634 - acc: 0.8966 - val_loss: 0.2732 - val_acc: 0.8922\n",
      "Epoch 49/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2814 - acc: 0.8887 - val_loss: 0.2847 - val_acc: 0.8859\n",
      "Epoch 50/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2778 - acc: 0.8894 - val_loss: 0.2647 - val_acc: 0.9014\n",
      "Epoch 51/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2739 - acc: 0.8898 - val_loss: 0.2765 - val_acc: 0.8937\n",
      "Epoch 52/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2819 - acc: 0.8873 - val_loss: 0.2843 - val_acc: 0.8928\n",
      "Epoch 53/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2720 - acc: 0.8934 - val_loss: 0.2713 - val_acc: 0.8969\n",
      "Epoch 54/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2685 - acc: 0.8928 - val_loss: 0.2776 - val_acc: 0.8933\n",
      "Epoch 55/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2636 - acc: 0.8955 - val_loss: 0.2858 - val_acc: 0.8917\n",
      "Epoch 56/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2654 - acc: 0.8941 - val_loss: 0.2760 - val_acc: 0.8932\n",
      "Epoch 57/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2661 - acc: 0.8941 - val_loss: 0.2851 - val_acc: 0.8887\n",
      "Epoch 58/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2587 - acc: 0.8963 - val_loss: 0.2663 - val_acc: 0.8966\n",
      "Epoch 59/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2729 - acc: 0.8902 - val_loss: 0.2716 - val_acc: 0.8927\n",
      "Epoch 60/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2760 - acc: 0.8912 - val_loss: 0.2801 - val_acc: 0.8912\n",
      "Epoch 61/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2725 - acc: 0.8934 - val_loss: 0.2816 - val_acc: 0.8918\n",
      "Epoch 62/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2729 - acc: 0.8925 - val_loss: 0.2765 - val_acc: 0.8925\n",
      "Epoch 63/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2790 - acc: 0.8910 - val_loss: 0.2753 - val_acc: 0.8952\n",
      "Epoch 64/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2654 - acc: 0.8986 - val_loss: 0.2721 - val_acc: 0.8965\n",
      "Epoch 65/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2710 - acc: 0.8930 - val_loss: 0.2746 - val_acc: 0.8930\n",
      "Epoch 66/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2722 - acc: 0.8942 - val_loss: 0.2687 - val_acc: 0.8952\n",
      "Epoch 67/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2728 - acc: 0.8934 - val_loss: 0.2759 - val_acc: 0.8892\n",
      "Epoch 68/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2802 - acc: 0.8906 - val_loss: 0.2791 - val_acc: 0.8927\n",
      "Epoch 69/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2683 - acc: 0.8928 - val_loss: 0.2755 - val_acc: 0.8949\n",
      "Epoch 70/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2734 - acc: 0.8918 - val_loss: 0.2731 - val_acc: 0.8980\n",
      "Epoch 71/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2810 - acc: 0.8867 - val_loss: 0.2776 - val_acc: 0.8972\n",
      "Epoch 72/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2674 - acc: 0.8958 - val_loss: 0.2793 - val_acc: 0.8900\n",
      "Epoch 73/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2638 - acc: 0.8951 - val_loss: 0.2694 - val_acc: 0.8971\n",
      "Epoch 74/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2674 - acc: 0.8953 - val_loss: 0.2690 - val_acc: 0.8957\n",
      "Epoch 75/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2700 - acc: 0.8945 - val_loss: 0.2776 - val_acc: 0.8902\n",
      "Epoch 76/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2709 - acc: 0.8936 - val_loss: 0.2707 - val_acc: 0.8951\n",
      "Epoch 77/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2676 - acc: 0.8951 - val_loss: 0.2764 - val_acc: 0.8920\n",
      "Epoch 78/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2699 - acc: 0.8927 - val_loss: 0.2833 - val_acc: 0.8908\n",
      "Epoch 79/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2683 - acc: 0.8943 - val_loss: 0.2753 - val_acc: 0.8955\n",
      "Epoch 80/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2721 - acc: 0.8904 - val_loss: 0.2808 - val_acc: 0.8901\n",
      "Epoch 81/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2653 - acc: 0.8948 - val_loss: 0.2859 - val_acc: 0.8870\n",
      "Epoch 82/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2781 - acc: 0.8907 - val_loss: 0.2751 - val_acc: 0.8975\n",
      "Epoch 83/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2695 - acc: 0.8911 - val_loss: 0.2706 - val_acc: 0.8956\n",
      "Epoch 84/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2761 - acc: 0.8898 - val_loss: 0.2673 - val_acc: 0.8985\n",
      "Epoch 85/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2765 - acc: 0.8927 - val_loss: 0.2779 - val_acc: 0.8916\n",
      "Epoch 86/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2681 - acc: 0.8957 - val_loss: 0.2701 - val_acc: 0.8987\n",
      "Epoch 87/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2710 - acc: 0.8927 - val_loss: 0.2787 - val_acc: 0.8908\n",
      "Epoch 88/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2608 - acc: 0.8980 - val_loss: 0.2686 - val_acc: 0.8971\n",
      "Epoch 89/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2697 - acc: 0.8944 - val_loss: 0.2761 - val_acc: 0.8909\n",
      "Epoch 90/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2688 - acc: 0.8940 - val_loss: 0.2729 - val_acc: 0.8980\n",
      "Epoch 91/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2715 - acc: 0.8922 - val_loss: 0.2776 - val_acc: 0.8931\n",
      "Epoch 92/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2705 - acc: 0.8919 - val_loss: 0.2820 - val_acc: 0.8934\n",
      "Epoch 93/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2634 - acc: 0.8964 - val_loss: 0.2718 - val_acc: 0.8960\n",
      "Epoch 94/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2760 - acc: 0.8907 - val_loss: 0.2641 - val_acc: 0.9000\n",
      "Epoch 95/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2794 - acc: 0.8862 - val_loss: 0.2726 - val_acc: 0.8956\n",
      "Epoch 96/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2796 - acc: 0.8890 - val_loss: 0.2797 - val_acc: 0.8939\n",
      "Epoch 97/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2690 - acc: 0.8943 - val_loss: 0.2715 - val_acc: 0.8985\n",
      "Epoch 98/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2827 - acc: 0.8907 - val_loss: 0.2816 - val_acc: 0.8937\n",
      "Epoch 99/100\n",
      "312/312 [==============================] - 7s 21ms/step - loss: 0.2646 - acc: 0.8956 - val_loss: 0.2708 - val_acc: 0.8965\n",
      "Epoch 100/100\n",
      "312/312 [==============================] - 6s 21ms/step - loss: 0.2790 - acc: 0.8876 - val_loss: 0.2660 - val_acc: 0.8978\n"
     ]
    }
   ],
   "source": [
    "history = model.fit_generator(trainGen, \n",
    "                              steps_per_epoch = trainGen.__len__(), \n",
    "                              verbose=1, \n",
    "                              epochs=epochs,\n",
    "                              callbacks = [schedule, swa], validation_data = testGen,\n",
    "                              validation_steps = testGen.__len__())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "8b92f57e-2138-44de-8b58-ae8d1d5808dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "conv_weights = model.get_layer(\"shared_conv\").get_weights()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "ec61c8dd-0ecd-49e5-8a9b-aa56a735b14d",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save_weights(\"weights.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "04d4c225-d40b-4649-926c-70e8c1f10503",
   "metadata": {},
   "outputs": [],
   "source": [
    "scan_model = construct_scan_model(conv_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "a35c466e-19c3-4ece-8715-38d6eb245b41",
   "metadata": {},
   "outputs": [],
   "source": [
    "pfms = [np.zeros((w,4)) for i in range(num_kernels)]\n",
    "n_instances = {i : 0 for i in range(num_kernels)}\n",
    "d = []\n",
    "for i, seq in enumerate(posSeqs):\n",
    "    encoded_seq = np.vstack((0.25*np.ones((w,4)), encode_sequence(seq), 0.25*np.ones((w,4))))\n",
    "    encoded_seq_rc = encoded_seq[::-1,::-1]\n",
    "    \n",
    "    conv_for = scan_model.predict(np.expand_dims(encoded_seq, axis = 0))[0]\n",
    "    conv_rc = scan_model.predict(np.expand_dims(encoded_seq_rc, axis = 0))[0]\n",
    "    \n",
    "    for k in range(num_kernels):\n",
    "        matches_for = np.argwhere(conv_for[:,k] > 0)[:,0].tolist()\n",
    "        matches_rc = np.argwhere(conv_rc[:,k] > 0)[:,0].tolist()\n",
    "        for x in matches_for:\n",
    "            motif_start = x - w \n",
    "            motif_end = motif_start + w\n",
    "            score = conv_for[x,k]\n",
    "            pfms[k] += encoded_seq[x:x+w,:]\n",
    "\n",
    "        for x in matches_rc:\n",
    "            motif_end = x + w\n",
    "            motif_start = motif_end - w \n",
    "            score = conv_rc[x,k] \n",
    "            pfms[k] += encoded_seq_rc[x:x+w,:]\n",
    "            n_instances[k] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1bfa073-cdc0-44f9-b7b2-522fb875a03d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "03e9f69b-6275-49b4-ac03-fb05ab140096",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZIAAAI4CAYAAADwNBxtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACDw0lEQVR4nOzdeZxcVZ3//9cnnQTIAgkkYUnCathE1mZRVKIo26i4oAIqroOoOM44zg++M24jzqIzKiIoMoiIM4qjMhIwghsCimASDZCwNgFCCEtCICEJJOn0+f1xK+mq6urq6u6q6ttdr+fj0Q/63nuq6nO6K2/y+NTJuZFSQpIkSZIkSZKk3owa6gIkSZIkSZIkSflmI1mSJEmSJEmSVJWNZEmSJEmSJElSVTaSJUmSJEmSJElV2UiWJEmSJEmSJFVlI1mSJEmSJEmSVJWNZEmSJEmSJElSVTaSlWsR8d8R8URErImIByLiQ4Xz74qItUVf6yMiRcQRhesREV+KiGcKX1+OiCh63kMj4taIWB0RyyLis0XXXhMRd0fEc4XH/l9ETG/+7CVJkiRJkqR8sJGsvPs3YM+U0vbAm4AvRsQRKaX/SSlN2PIFfBRYAvy58LizgTcDhwAHA28APlz0vD8AbgF2BI4DPhIRbypcuwc4MaU0CdgNeBD4VuOmKEmSJEmSJOWbjWTlWkppcUppw5bDwtc+FYa+F7gqpZSKjr+SUlqWUnoc+ArwvqLxewL/k1LanFJ6CPg98NLCaz6VUlpeNHYz8JI6TUmSJEmSJEkadmwkK/ci4psRsR64D3gCmFt2fQ/g1cBVRadfCtxZdHxn4dwWFwJnRcSYiNgPeDnw66Ln3D0ingNeAD4FfLle85EkSZIkSZKGGxvJyr2U0keBicCrgGuADWVDzgJuTSk9XHRuArC66Hg1MKFon+TrgdPIGsX3Ad9JKc0res2lha0tpgCfLoyRJEmSJEmSWpKNZA0LhS0ofg/MAD5Sdvks4Htl59YC2xcdbw+sTSmliNgRuAH4ArAtMBM4MSI+WuF1VxWe+9qIGF2XyUiSJEmSJEnDjI1kDTejKdojOSKOJbsh3k/Kxi0mu9HeFocUzgHsDWxOKV2VUupMKS0DrgZOqfKa0yhtTEuSJEmSJEktw0aycisipkXE6RExISLaIuJE4Azgt0XD3gv8NKX0fNnDrwI+GRHTI2I34O+BKwvXHsiePs6MiFERsQvwTgp7KkfEWyNiv8K1qcBXgb8UVidLkiRJkiRJLcdGsvIskW1jsQx4FvhP4G9TStcCRMS2wDvoua0FwLeB64C7gUXAzwvnSCmtAd4K/F3heRcWxvxL4bHTyba+eL7w+C7gLfWenCRJkiRJkjRcREppqGuQJEmSJEmSJOWYK5IlSZIkSZIkSVXZSJakYSYiroiIpyNiUS/XIyIuioiOiLgrIg5vdo2S1CrMZEnKDzNZkhrLRrIkDT9XAidVuX4yMKvwdTbwrSbUJEmt6krMZEnKiysxkyWpYWwkS9Iwk1K6BVhVZcipwFUpczswKSJ2bU51ktRazGRJyg8zWZIaa/RQvfCUKVPSnnvuOVQvrxEmJXjqKYiAnXce6mrUihYsWLAypTR1qOsomA48VnS8rHDuifKBEXE22WoMxo8ff8T+++/flAIlqdFylMtmsqSWZyZLUn4MJpOHrJG85557Mn/+/KF6eQ3A88/DaafBX/4CF18M73jHUFeU2bwZ3vlO+POfs+MjjoBrroFRrrdXE0XEo0NdQ5GocC5VGphSugy4DKC9vT2Zy5JGihzlspksqeWZyZKUH4PJZFttqsnGjfDWt8IvfwkrVmSN27lzh7qqzBe/CD/9affxtdfCf//30NUj5cAyYGbR8Qxg+RDVIkmtzkyWpPwwkyVpEGwkqyb/9V/w61+XnvuHf4CurqGpZ4sNG7LV0eW+9KWhr63p1i+DDdW2A1MLmQOcVbgr9THA6pRSj3+uJ0lqCjNZkvLDTJakQRiyrS00fKQE3/52z/P33JOtSn7DG5pf0xbXXgsrV/Y8f889cMMNcMopza9pSDx0Ofzpr2HMDvDy78P0Nw51RWqgiPghMBuYEhHLgM8BYwBSSpcCc4FTgA5gPfD+oalUkkY+M1mS8sNMlqTGspGsPt1xB9x9d+VrF100tI3k//qv3q/94Act0kh+6iaYd072/abV8Pu3wQnzYfLBQ1uXGialdEYf1xPwsSaVI0ktzUyWpPwwkyWpsdzaQn36znd6v3bTTbBqiHZTWL0afvOb3q9ff322t/OIt+gCSJu7j7s2wV3/NHT1SJIkSZIkacSxkaw+VWvWdnZmDduhcPvt2bYbvVm9Gn73u6aVMzTWPgxP39Tz/PLrYdWC5tcjSZIkSZKkEclGsqp64gl4+OHqY37+8+bUUu4Pf+h7zOLFja9jSD38vd6vdVzavDokSZIkSZI0otlIVlW1NGtXr258HZXUUtuI99hPq1z7v2ybC0mSJEmSJGmQbCSrqt//fqgrqKyzM7sJYEvbtAZWV1lyvfGZ7EZ8kiRJkiRJ0iDZSFZVeV31u3gxrFs31FUMsWfmAVU2iQZYeVtTSpEkSZIkSdLIZiNZvershLvuGuoqKhvxex/X4pnbaxjUR6NZkiRJkiRJqoGNZPVqyRLYuHGoq6jsnnuGuoIcWFlLI1mSJEmSJEkaPBvJ6lWem7V5rq1pVs0b6gokSZIkSZLUImwkq1f33z/UFfTu4YeHuoIhtvFZePGpoa5CkiRJkiRJLcJGsnr16KNDXUHvli4d6gqG2Jocd/klSZIkSZI04thIVq/y2qxduxZWrRrqKoaYjWRJkiRJkiQ1kY1k9SqvK5Lz2uBuqvX+EFpZRJwUEfdHREdEnF/h+g4RcV1E3BkRiyPi/UNRpyS1AjNZkvLFXJakxrGRrF7ltWHbW11Tp8Jb39rcWobM+mVDXYGGSES0AZcAJwMHAmdExIFlwz4G3JNSOgSYDXwlIsY2tVBJagFmsiTli7ksSY1lI7lFPf88nHEG7LMPfO1r0NVVen31alizpufjXvpSWLQoe+xQqdRIHj0arr0WfvpT+Nd/bX5NTfeCjeQWdhTQkVJaklLaCFwNnFo2JgETIyKACcAqoLO5ZUpSSzCTJSlfzGVJaiAbyS3ouefgta+Fq6+GJUvgk5+Ef/qn0jG9bWvx3/+dNZOvvBKOPLLRlVZWqZH82tfCy1+efX/++XDYYc2tqenWPTbUFWjoTAeK3wDLCueKXQwcACwH7gY+kVIq+7hIklQHZrIk5Yu5LEkN1GcjOSKuiIinI2JRL9dnR8TqiFhY+Pps/ctUPX31qzB/fum5//gPuPfe7uNKzdqXvAQOOST7fuxY+PSnG1djNY9V6KH+1V91fx8BH/5w8+oZEpVWJO90DBz8LzBqTPPrUTNFhXOp7PhEYCGwG3AocHFEbF/xySLOjoj5ETF/xYoV9axTklqBmSxJ+VK3XDaTJamnWlYkXwmc1MeYW1NKhxa+vjD4stQoq1fDRRf1PL95c7aSd4unnuo55o1vzJq0W5x0EkyaVPcS+/Tccz3PFTeSIdsredRIXW/fuQ42Plt6bpspcNx18NJ/hKOvHJKy1DTLgJlFxzPIVlMUez9wTcp0AA8D+1d6spTSZSml9pRS+9SpUxtSsCSNYGayJOVL3XLZTJaknvpstaWUbiHbM0gjwDe/mTWTK7nuuu4tLSrtj/yGN5Qejx0Lb3lLfeurRXlt++6b7fVcbOrUbLsLIOt+9+cr7yrdaG/Gm7NmMsCeZ8K045pakppqHjArIvYq3BTkdGBO2ZilwPEAEbEzsB+wpKlVSlJrMJMlKV/MZUlqoHqt2Xx5RNwZEb+IiJf2Nsh/GjL0fvrT3q+llO19DD2btW1tcOyxPR/z1rfWrbSaldd2+OGVx514YuNrGRIvPt3z3G5vLD1+yTnNqUVNl1LqBM4FbgTuBf43pbQ4Is6JiC2/+AuAV0TE3cBvgPNSSiuHpmJJGrnMZEnKF3NZkhprdB2e48/AHimltRFxCvAzYFalgSmly4DLANrb28v3KVKDPfUULFhQfcx118HnPtezWbvPPrDNNj3Hv/KV2SrnZnr++dLj/Sv+41B41avgttsaX0/TdZb9ANq2hV1eV3pu+qkwemLzalJTpZTmAnPLzl1a9P1y4IRm1yVJrchMlqR8MZclqXEGvSI5pbQmpbS28P1cYExETBl0Zaq7X/6y7zGdndl/yxvJvTVrJ02Cgw8eVFn9Vmtthx8O48Y1vp6m21TWSN7pGBhdNtHR28HMIVguLkmSJEmSpBFp0I3kiNglIttYNiKOKjznM4N9XtXfDTfUPrbWVb+QrfxtpvJG8n77VR43Zgwcc0zj62m68hXJOxxYedyufd0jU5IkSZIkSapNn1tbRMQPgdnAlIhYBnwOGANb/3nIacBHIqITeAE4PaXkthU5dPPNtY8tb9buvXfvYyvtndwomzbBCy+Untt9997HH3RQY+sZEuUrkrfvpcu/82thzX2Nr0eSJEmSJEkjXp+N5JTSGX1cvxi4uG4VqSGefBIef7z28eWN5OnTex87adKAShqQ8pXS22wDkyf3Pr6trbH1DInyFckTe2kkbzsNJr2s8fVIkiRJkiRpxBv01hYaHvq6yV658kbybrvVr5bBKK9r110h21ilhWxaW3o8Yc/ex06b3chKJEmSJEmS1CJsJLeIe+/t3/hKDds8yGtdTVW+InnbKj+EbXZqbC2SJEmSJElqCTaSW8SDD/ZvfHHDdtQomDatvvUMVPnWFrvsMjR1DKniPZJHT4QxE4auFkmSJEmSJLUEG8ktoj+N5JRKG8mTJ+dnr+HyFck77DA0dQyp4hXJ27XikmxJkiRJkiQ1m43kFtFbI/mww2D8+NJzL7wAmzd3H0/I0YLX8kZynmprmuIVydsOcEl26oJFX4SfTIK5B8PyuXUpTZIkSZIkSSOTjeQWsH49LFvW8/wnPpHdhO+PfyzdIiLPzdo819Y0xSuSxw5wSfaCT8Ddn4FNq2H13XDzX8Fj19SnPkmSJEmSJI04NpJbwEMP9Ty37bbwhS9ABLzsZXDJJd3X8tysLa+tfDV1Syhekdw2gB/Ac3fDg5f0PP+nD8H6Cp84SJIkSZIkqeXZSG4BS5f2PHfqqbD99t3Hb34z7L139n35De3y3EjOU21NU7wiefQAGsl3/j8g9Ty/8Vm4+3MDLkuSJEmSJEkjl43k4S4l6Pgv+PmB8NvXwdM39xjy7LM9H/bud5cejxoF55yTfZ/nVb95rq1pNg2ikfzCk9X3Q374KlhX4ZMHSZIkSZIktTQbycPdnf8P5p0Na+6Fp34Dv3kNLP1JyZDyRnIEzJ7d86ne9a7sWp5X/eZ5tXRTpASda7uP+9tIfvw6Kq5G3vr8nd54bxiIiJMi4v6I6IiI83sZMzsiFkbE4ojo+QmTJKkuzGRJyhdzWZIaZ/RQF6BBePYuuO8/yk4muOMDMPkQmDgrG1bWSN5338oN2N12g8MOy3cjeUC1pQqN04jq1/Oqcx0ljeDR/fzlPP6zGgYNo59HC4qINuAS4PXAMmBeRMxJKd1TNGYS8E3gpJTS0oiYNiTFStIIZyZLUr6Yy5LUWK5IHs4W/gOkrp7nO5+Hhf/f1sPyRvJhh/X+lKecAmvXlp7L0/YRL7xQepyn2pqis+yX058VyZtfhCd/Xd96NBSOAjpSSktSShuBq4FTy8acCVyTUloKkFJ6usk1SlKrMJMlKV/MZUlqIBvJw9W6x+DJX/Z+fdnPYNUCoGcj+cADe3/Y618PnZ2l5/K0Innz5tLjPNXWFJ1le3v0p5G8agF0baxvPRoK04HHio6XFc4V2xeYHBG/i4gFEXFWb08WEWdHxPyImL9ixYoGlCtJI5qZLEn5UrdcNpMlqScbycPVYz/pe8yK3wM9G8l77dX7Q3bYId/N2vLaWm5F8uYXS4/b+vEDWHl7fWvRUIkK58r3IxkNHAH8FXAi8JmI2LfSk6WULksptaeU2qdOnVrfSiVp5DOTJSlf6pbLZrIk9eQeycPVYz+ueWh5I3mPPaqPH06N5DzV1hSp7AfQnxXJqxfVtxYNlWXAzKLjGcDyCmNWppTWAesi4hbgEOCB5pQoSS3DTJakfDGXJamBXJE8HG1YBSv/WPPw8kbynntWH1/erB03ruaXarjy2rbddmjqGDLle2L352Z7a+6rby0aKvOAWRGxV0SMBU4H5pSNuRZ4VUSMjohxwNHAvU2uU5JagZksSfliLktSA7kieThaeVu/hpc3knfZpfr48mbt6By9S7rK+qhtbUNTx5AZ6IrklGBNb383Cnr+ay/lVUqpMyLOBW4E2oArUkqLI+KcwvVLU0r3RsQNwF1AF3B5Sskl6ZJUZ2ayJOWLuSxJjZWjFqFq9swd/Rpe3EgeNw7GjKk+vryRPCpH69bzXNuARKUtvKpYUbYSfdTY2h734lOwaXXpuTGT4LjrYPsD4Y73w+PlH9Qrr1JKc4G5ZecuLTv+D+A/mlmXJLUiM1mS8sVclqTG6bMNFxFXRMTTEVHxE7rIXBQRHRFxV0QcXv8yVWL14pqHbtwI69d3H2+/fd+PKW/W5mnVb55ra4ryrS2i6I9w57psy5OSr8IN9iqtRt7/kzD1lbDNjvDKn8CO7Y2rW5IkSZIkScNaLSuSrwQuBq7q5frJwKzC19HAtwr/VaP0o5Fcvq3FDjv0/Zg8r/rNc21NUb61RXEjee0S+NUrSq+PGgPv3AgvlN1fYtRYeMmHS8cd+h/w29fUt15JkiRJkiSNCH224VJKtwCrqgw5FbgqZW4HJkXErvUqUGU2vwhrO3qef8mHof0S2HbnktPljeSBrEgeimbtunWlK6m3cEVyeSO5xh/AhmdKj6e/CbadVnpu59mw82sHXJokSZIkSZJGrnq0CKcDjxUdLyuc6yEizo6I+RExf8WKFXV46Ra05v6e2xvscQYceSnM+iic+GfYrruPX49GcjObtTffDCecABMmwPjx0N4ON97Ye22ttyK5ytYW1Wws+yxo6qsqjytepSxJkiRJkiQV1KMNV+luYanSwJTSZSml9pRS+9SpU+vw0i1o7ZKe5w44r/v7cbvBoV/ZejicViRffTW87nXwq191n1uwAE46Cb72tey4q6yP2nqN5LJfTq1/hMtXJE95eeVx098EYyf3uyxJkiRJkiSNbLXskdyXZcDMouMZwPJexmqwNqwsPZ42GyYfUnpuj3fC/V8Fem4PUUsjOSp9NNBgd94J7353zyb2Fp/8JOy9d8/aUsWPLEawanskV7OxrJG8w0GVx7VtC7u9of91SZIkSZIkaUSrx3rOOcBZkTkGWJ1SeqIOz6tKyhuCu7y+55gYBfv9LQCdnaWXarnZXvkq3/JVwPXW1QXnnNN7E3mLP/6x5zYbja4tf8o75zV2/TcUbW0xbgaM3q73sWMm9LsqSZIkSZIkjWx9rkiOiB8Cs4EpEbEM+BwwBiCldCkwFzgF6ADWA+9vVLGi54rk3va6nfk2WPrjHs3ZWlYklzdr+2rwDtacOXD77bWNbXZtudNjBXKNS7KLP4CYMKtu5UiSJEmSJKk19NlITimd0cf1BHysbhWpuvJG8uRDK49r2xZmvLnHiuTx4/t+iWav+v32t2sf2/IrkqPsB1B+873eFN9sb/zu9atHkiRJkiRJLaHVblU2/BXfNG273WDMxN7HjpnYY8Xu6Bp2xW7mqt/nnoPf/Kb28S2/IrnHH9kaG8nF7xtvpidJkiRJkqR+spE83BSvSJ64b5/Dy1cklzdiK2nmqt+f/xw2bSo9N25cdn7Vqmzbi332GZracmnUAFYkd3XCptXdx2Mm1bUkSZIkSZIkjXw1rE9VMy1bBjfeCE88AdOmwQknwJ57Fg0obiSPm9Hn85Wv2M1bI/l3vys9Hj0arrkGTjwxO37jG+Goo+DYY7Pj8hsB5npF8qbn4aHvwLKfAqNgp6Ng1sdgwp7dY1LZHsdRdvO88usrfl92vYYfwMZnS4/HTur7MZIkSZIkSVIRVyTnxKZN8KlPwcyZ8KEPwWc+Ax/+MOy1F7zvfdkWEEDpTdPG7tjn89ZjRXL5iuFa3XorvPOdMGUK7LADvOIV8J3vlDamH3ig9DHvfnd3E3mLnXeGH/84azKX17Zx48Bqa7jO9XDTCfCXv8uavytugfv+E+a+FB7+/iCeuOyPbC0rkov3RwZXJI8AEXFSRNwfER0RcX6VcUdGxOaIOK2Z9UlSKzGTJSlfzGVJahxXJOfAmjXwhjdkjddKvvc9+NOf4J5FnaWrS2vY67YeK5LXr+/7MeW+/W342MdKX/+Pf8y+vv99+L//g8mT4cEHSx/3oQ9Vfr7DDoOddsp+DsXWret/bQ3X1Qm3nQ7P3N7z2ub1cPtZMGoM7HF6/5+7x832aliRXLw/MrgieZiLiDbgEuD1wDJgXkTMSSndU2Hcl4Abm1+lJLUGM1mS8sVclqTGckVyDvzd3/XeRN7i3nupsEVB/xvJ5VtDVFLeSF67tu/HFPvKV+Ccc3rfduLmm+H887PnfeKJ7vM77ABHH9378+6+e8/actlIXvJdePy66mMevXpgz13eSO7a0PdjOst+SDaSh7ujgI6U0pKU0kbgauDUCuM+DvwUeLqZxUlSizGTJSlfzGVJaiAbyUPs5pvhiitqHFy8PzLU1EgeyH7Hg2nWPvAA/NM/9T1uwwbo6Cg9d/TR2fYV/amtv03uhuvaDPd+uXHPH2V/ZMubxJWksv1N3NpiuJsOPFZ0vKxwbquImA68Bbi0ryeLiLMjYn5EzF+xYkVdC5WkFmAmS1K+1C2XzWRJ6slG8hD7z//seW7HHbO9gt/6Vth++6IL5XvdDqCRXMvN6cqbuf1p1v7TP2VN4lqU74+81159P2bMmNLj3DWSn/wlrO3oe9xAla9IHlAjeYf61aOhEBXOld2VkQuB81Lqe++TlNJlKaX2lFL71KlT61GfJLUSM1mS8qVuuWwmS1JP7pE8hB5/HK6/vvTckUfCT36SbeMA2f7Jf/M32T7JdJXdWa6GRnJ5U7iWRvLEiaXHtTZrly/P9j4ud+ih0N6eNY5vuaX7fPn+yHvs0fdrlDTWyeHWFo9fW/l8jKrtxnh9GVX2Cy1pJAe0bZu9TvF7pav8jovbDr4ODaVlwMyi4xnA8rIx7cDVEQEwBTglIjpTSj9rSoWS1DrMZEnKF3NZkhrIRvIQKr9x3OTJMHcuTJnSfW777eG734UJE+jZEGzQiuTyZm2tjeQf/KDn859zDlx0UfdK4j/8AU4r3BP30UdLx+65Z+Nqa4qUeu6NPHoiHHcdTDkWVs2HP/8tPHPHwF+jbVzpcXEjedJB8I4XYM198PMDiuoq+6UUr2p+6ib4/dtKr283HU65e+A1qtHmAbMiYi/gceB04MziASmlrev7I+JK4Hr/YixJDWEmS1K+mMuS1EBubTGE5s8vPf7IR0qbyFtEwNe/Ts8tCtrG9/ka9ViRXOuq3/L5nHgifPObpdtRHHtsti/05Mnwwgul47eswq4m1yuSX3gcXij6sHvUGJh9A0w7LltJPOUYOP53sMuJA3+NMWW/nIFsbVHcSO7qzG7iWPy1afXA61PDpZQ6gXPJ7jB9L/C/KaXFEXFORJwztNVJUmsxkyUpX8xlSWosVyQPofLG69veVnkcFFYWlzcEi7c5WH0PLPpC6fUxO9DW9u2SU51lT1HJQFf9LlxYevxP/5Q1wcvtuy/88z/DBz9Yen7SpMbV1hTPl+2NvPvpMPUVpefatoVjfwR3f25grzG6rJG8eX3fjylfyV6+z7KGnZTSXGBu2bmKNwtJKb2vGTVJUqsykyUpX8xlSWocG8mNsn45LLkcnvotxBjY8XDY+4Ow/b5AtgtCcSO5rQ0OPLCP5+zRECz69b34NCz9Uen1bXdm9OjSRnItK3gH0qxdt6705nm77pqtPq72Gi++WHpum20aU1vTrH2o9HjPd1ceN3YHOPiCgb3GqLHZ733LhwqDXZEsSZIkSZIk1cBGciOsvgdueh288ET3uad+Dfd9FQ76LLz00yxbFqxa1X151izYtq97oA2gIVi+R/KaNX0+ZEDN2rvvzprjWxx2GIzqY+OUDRtKj4u3wOjNQLfdaIriRvKobWDqq3ofW75FRa0issdufDY7rqWRLEmSJEmSJA2SjeR6W3M//OY42LCy57XUCXd/FmIUz27+p5JLBxzQc3gPA9iioHyP5IE0kmtp1pZva7H//n0/prxx3MgbATZFcSN5wj4wervGvM7ookby5hp+OVH2Jii/+Z4kSZIkSZLUBxvJ9fbnv63cRC72/EOsS6WnatkfGFLfQ8qMG1d6XEsjedttswb0lv2Un3++78csXVp6XEtjvHwFdvkK5UrKG8m11NY0xXskj9+jca8zekL397WsSB5lI1mSJEmSJEmD08fmA+qXp2+GJ26oaegLL5Qe17I/8EAagpMnlx6vXt33y0SUNmyffbbvm/StL7vn2+679/065Y3kgezf/PTTfT+maYpXJI+v4QcwUMXbYtTSSHZFsiRJkiRJkgbJRnI9PXRFzUPLG8ebNtXwoAE0BMsbybWsSIbShm1KfTdsyxvJ5SuhKylvJC9b1r+6AJ54ovK4pktdsKmoSz+ukSuSixrJW7a4qKbHBxB9fCogSZIkSZIklampkRwRJ0XE/RHRERHnV7g+OyJWR8TCwtdn61/q0NnctZnOrs6tXxV1vgDLrik9N35POGEevH09HH9Lyc3Xxo8vHVrTXr9D1EgGWL68+viB3Dhv551Ljx95pO/H9LeupunaWHq83S6lxxuf6/nVNcCVwWOKtrZ48cm+x5e/bzbXsIeIJEmSJEmSVKTPPZIjog24BHg9sAyYFxFzUkr3lA29NaX0hgbUOOT2v2R/OlZl+99OGDuB5857jrZRZTe6W303dBZ1g8fNhNf/sbuhOO1V8NrfwvyPQNdmJk4sffiDD9ZQyABWljarkbxd2X3latnveN99S48ffbT/da1Zk22JUd6Yb7rNL5Yejypbcv7THemxx/Up98AOtdxlsUzxiuQXaliSXX5Txk017G8iSZIkSZIkFallRfJRQEdKaUlKaSNwNXBqY8vKj1UvrNraRAZYu3Et9628r+fA1feWHh/4jz1XpY4aDe2Xwm6nsMcepc3Xe+6BzX0tUC1fWdq5vvK4ImPHlm4zUcseydCzYfv4472PXbGi51YWtaywnjWr9LiWFcnbbNNztXO1Jnf5TQAbpnyV76ixjXutkj2S18KmPn7Y5e+bTc/VvSRJkiRJkiSNbLU0kqcDjxUdLyucK/fyiLgzIn4RES+tS3U5MH/5/B7n5i2f13PgmqIF2jEa9nhn5Scc1Qa7n8bo0XDYYd2nX3wRHnqo8kNKnrdYLfvjUroqef362vZjLm8kL1nS+9jf/77niuBq47cobyQvXJjtx1xN+Y0AoffV3CnBrbf2XUefNq2Bjm/Dr46Fn+0GP90Jrt8Pfv92WHJlTQ39uhpdtpz9xT5WJY8p+4FtLPo0YefXwGmr4dXX16c2SZIkSZIkjUi1NJKjwrnydt+fgT1SSocA3wB+VvGJIs6OiPkRMX/FihX9KnSo/OnxP9V0jjVFK5InHwJjJ/ccU+bII0uP587t4wFtZatci1eWbjsN9jgTdnl9j4eVb2/xZA3b6pY3a++9t/I4yJq15U3hauO3mDoVdtih+/iRR6o307fc0K98W5D7KiwQh6zB3NdNAvvU8W34v11h3jmw8rZsK4mNq+D5B+Cxn8Ad74ebT4G2sq0syvdMrqfRE0qP1/bRtd9mx9Lj4vfNqNFZo3n0UO8NIkmSJEmSpDyrpZG8DJhZdDwDKNlMIKW0JqW0tvD9XGBMREwpf6KU0mUppfaUUvvUqVMHUXbzVFp9XHFFcvHWFhNm9bxeQXt76fEll/S+t/DllwNjyxqCxSuSdzgQXvE/cNDnezy2vJFcyxYStTZrV67MmsaHHlrb+GIRPRvQ3/9+5bErVsA3v5l9X2uTe9CrkZdcmTWQNxetOG7bDsbvUbrKd/OL0LZt6WPL90yupzFlv5w191cfX/6+cY9kSZIkSZIk9VMtjeR5wKyI2CsixgKnA3OKB0TELhERhe+PKjzvM/UuttlSSltXH++43Y4cNO0gAO588k42dBZ1fDdvgHVFq0InvqSm5y9fkdzRAeeck21z0V1D1kA95xxg7E6lD9j4XE2vU95IHshN7R5+GF54oee4W27Jatx//9K9ixcsqO2Ge+WN5G9+Ex57rPRcZyd88INZ07pSbYsXV37um27q+/V7tWkNLDi3+3jMJDjyMnjrSnjTI/C25+Dku2Dfv8lurFe+J/KLNSz7HqjyrS2ef6D6+DGTKPmHBTW+b5RfEXFSRNwfER0RcX6F6++KiLsKX7dFxCFDUacktQIzWZLyxVyWpMbps5GcUuoEzgVuBO4F/jeltDgizomIcwrDTgMWRcSdwEXA6Sn1tdtt/j3+/OM8uTZrCB62y2EcvuvhAGzq2sSdT93ZPXDTGkhd3ccT9qnp+ffdF444ovTclVfCscfCN74BX/4yvOIV8LGPFW7EV75FwQD2SIasKdyb557L/lverO3qgptv7jn+pz/N/jt2LLy0aGfs1avh5z/v/XW27KFc3kxfuRJOOCFbgQywahW8611w3XXdY8pr+9Ofuuve4sUXSx/Tb49fB53ruo9fdQ285K9hdOGughEw6WVwxNfh1ddCjCrdzmR92V3+TvgTnDCv5z7XAzGmbGuL1Yuqjx/VBmMndR8PcEXy5s3wk5/Apz8N554LF1yQrfoe/n/Sh5eIaAMuAU4GDgTOiIgDy4Y9DByXUjoYuAC4rLlVSlJrMJMlKV/MZUlqrJq6WoXtKuaWnbu06PuLgYvrW9rQm/d49xYWh+96ONMnTueqO6/aeu2o6UdlFzeXLdXdpmhXj80vwr3/0fPJ9/s7YswE/v7v4cwzSy/9+c/ZVw+jxsCYHbobgcV73VZR3ki+557K4wB+9St4+9t7NmsBrr8eTjqp+/iFF2DOnKzxDdn2FgsXdl///Ofhr/4KtinbPviGG7KG5OWXw1veAp/8ZOn1++7LViq/5CXZthXry+5lV17b5s3wy1/CO97Rfe6Xv4Q1a3qfZ5+W/rj7+x1elt2UrjdbmrQT9oFVhZszri9bVr1TYR+TGNVzh/H+Kl+RvPI22PR8zy0vSmrcsfuDhxo/gNiiqyv7wODzn6/83jnooGxblle/ul9Pq4E7CuhIKS0BiIirgVOBrb+dlNJtReNvJ9uSSJJUf2ayJOWLuSxJDVTL1hYtq/imeoftchiH7XpY97XlRTfcK98Pt3i/3M0vwt2f7flVWO36jnf0swFX3KQe4Irkv/yl97FbVhHvskvPa9dfX7r6dM4cWLu2+3hLQ3mLu++G44/vvuHdhg1w0UXwpjdl21UA7LknHHYYPaxenW2PUd5EBthtt8q1Fbvqqp5j+mXlH7q/3/XE0mtdm3p+pa7SlejrylYk19OYHUqPuzbBk7+u/pjifZLX1bC3yZan7oL3vS97n/b2AcSiRfCpT9X8lFXNe3wen73ps1u/7nrqrvo88cgyHSj+pGJZ4VxvPgj8oqEVSVLrMpMlKV/MZUlqoDr8O/uRq/imeofveji7Tty1+1rRamVSZ+kDo63m12hry7azeOUrYfny3sdtvfndNjvB2oey79cvq+k1yhvJDzyQNYAnlO2Q8Pjj3SuK99ij5/M8+ijcfju8/OVZg/eCC0qvv/Od8Hd/V9pc/sMfYJ99stXFDz+cNYjLvfe91Zvb5Xbfvee566/PtsOYOjXbbuGaa2p/vh5SKm3Sjyv7gHrOHvDCE6XnXvvb0kby8w9WvglfPWxXoZP+xFyY+Zbs+3WP9bxe3Eje8v6pwX/+Z+83QGyEi+ddvHXVP8CzLzzLN075RvMKGB6iwrmK69wj4jVkfzl+Za9PFnE2cDbA7pX+cEmSqjGTJSlf6pbLZrIk9eSK5F50pS7mL8+2KZgwdgKzdprF9ttsz6wds7vD3bfyPtZsKOydUN4s3FzDXeaK7LVXtpXF7NmVrx90ENy25R/fFN9w7/kHa3r+8kZySvC73/Uc99//3f39zJmVn+td78oee9ppPW9yN3Fitnq13Nq1WYO6UhMZ4AMfgF13rXytkkr/D3/22aym7343W/E8qH17I7JtRLbo2lTb44obyV0bYMVtvY8djEqN5Ed/lL0fujbDvf/e83rx/trrl8LmjX2+zPLl8IUvlJ4bMyb7sOBHP4JvfQve8IZ+1l7F5q7NzH2wZAcdrn/wekbAduv1tgwo/hM6A+jxMVREHAxcDpyaUur15qcppctSSu0ppfapU6fWvVhJGuHMZEnKl7rlspksST3ZSO7Fg888yOoNWefzkJ0PYVRkP6ot21skEguWL8gGjypvJFfYj6EPO++c7U98zTVw9tnZ3sIf+AD8+MfZat2DDioMLN7aYv2jNTUEd9yx57nipjFkWxh8+9vdx9tvD5Mm9Xzcww/Da14Dv+jlH/98/vPZquD+mDgRLryw73HTpmX/7e3D4FtuyX5m5TfeG5DiZu1zZdsr7P1B2PdvYHTZku6JZTdZfPR/6lBIBW1jYdudS891Pg+/mQ03HgFP/bbnY4o/gEhdsObePl/mhz+EdUX3G5w2De66C7761Wyri3POyW5oePPNlbdC6a87Hr+DletXlpx75LlHWLxicS+PaFnzgFkRsVdEjAVOB+YUD4iI3YFrgPeklB4YgholqVWYyZKUL+ayJDWQjeRelG9rsfX7XQ7vOWb0uNIHF++PG6Oylarj9+r5IhElX6PHBG95a/Dty4Lrfx5854rgtLdn54nCv9DZpqwhuHpRn3OptLr4Zz8rXSH8f/+XNYmLVdreoi877QTf+Q6MHVt9XPkK5Le/Hf76r3sf/8Y3wic+MfC6+m23v+r+/okboato+5KDL4Ajvl66XQSUrkgGeOR/YOUfG1Nf+XYbAC8sh+furDx+m7Jan7m9z5f4ddm2yxddBPvv33Pcq18NP/hBn0/Xp+sf6N7o+u0Hvr3ieUFKqRM4F7gRuBf435TS4og4JyLOKQz7LLAT8M2IWBgR84eoXEka0cxkScoXc1mSGstGci+Kb7T3vTu/x4yvzmDGV2fwpT98qeeYMTuU3gCteA/aMdvDGzvglLvrU1jximQovSlcL/bZp+e5DRvgM5/JtoC480742Md6jhnoNlBvfCPcdFPllckR8MlPwhe/2PP8t78Nn/tc6d7NY8bA3/89/PSn2X7SkDWrt9tuYLXVbGZ3I5MNT8Pif+n7MdtNh3FFXfuuDfC7U+Cp38GmNbDke9DV9wrymmzXzxsLlze9n7616vANG7IV3ltMnpxtHdKb8v22B+K6B64DYPSo0fzb8f/W47y6pZTmppT2TSntk1L6l8K5S1NKlxa+/1BKaXJK6dDCV/vQVixJI5eZLEn5Yi5LUuOM/JvtRaW99qso7MdavCJ5zYY13fshF9k6JgK2P6B7lWc/bmbWb8VbFEC2jcG+H6/6kHHjYMYMWFZ2b75vfCPblqCjI7t5XrnBrPzd99CVXP7zR/jVtVOY/4cdGMt4Dj14LO9/Pxx6aGFQ2e8mgM8Df88EfsPxBImjNv2JXb/yJHylMCglIrIm9/33D7y+Pk09FnZ4GawufACw6POw4vfwkr/OGsbrH4NNz5U+JgJ2ewN0fKv73Kbn4LevqX9946rdeLiC8kbysmuyGwqOnVxx+B13lL4njjmmu5HfCI889wiLns5W1x8781j22XEfDtvlMP7y5F/442N/ZOX6lUwZN6WPZ5EkSZIkSVKjuCK5go2bN/KXJ/7S57ilq5fy1NqnsoPtD+i+sHpRbXd7S6n0q5br5SuSH7+udCuNXsyaVfn8XXdVbiLD4BrJF95+IafOPZKLx+zF7bN35DWf+1e+/vWiJnIVE1nLm7mWU5nDrjxZ99pqEqPg2KtLf95P/Rr+8E749SvhtjOyVcblpr+xwYUV9HdF8viyH9jmF2DJld3HG1eVXH6o7LOQo4/u38v1V/H2FSfuc2LJfxOpx034JEmSJEmS1Fw2kitY9PQiNmzeAEAQbDt625KvLTfeg6JVyTsUNZLXPwbP1Wkri3Ll+/CmzXD/17uP13ZUfFhvjeRq9tuv/4/Z4hcdv6h6PFiV9uqtux0OhNffBru/A6KXxftjdoD9PwWTs5swsssJMHEAP+z+mlBhz+1qJlb4ZS7+Ijx1Ezy7EBb+fyWX1pT1yOtxM71qihvJJ73kJABOfMmJFa9LkiRJkiSp+Yb11hYn/8/J3PJo90auN733Jo6aflTpoPKVvuVbXVRYCTzv8e5tLc4+4mwufcOlJde/eMsX+cxNn9k69g37vqF0RTLAgxfDUZfVOJN+2H4/sk0giuq+/2vZDf+2mQJ3fabiwwbSSD7wwAFVyJNrn+TPT/y55Ny8x+exYt0Kpo4v2ji50irs4t9PlVXdA62t3ybOgmN/xNy7rmDCU79g4qZniYCD9ziBUTseAVNfBW1FdxYc1ZY1lud9uPrzbjMVtp02iLr27d/4bafBmEml23FsXAW/fW3F4RvLtnIe3cCkeH7D89z0yE0AjG0by70r7+X+Z+5n4+aNBEEicUPHDWzcvJGxbX3cxVGSJEmSJEkNMWwbyQ8/+zA3dNxQcu67f/luz0byABTfaO+IXY/ocf3wXQ/vHru8MHaHg0oHPXQ5TD4cZhVuDNvVOei6gKxhPGEvWLuk6GTKVpdWMZDG6557wrbbwosv9u9xxb+X3XfYnaWrl5JI3PjQjbz74Hf3v5AKmtZIBtZuXMvbr/846zd17wFy84Gf5dW7vLryA/b+ICz/OTw+p/L1tm3hVT+DbXaqfL0W/W0kR8D2+3fv493X008sPV67tn8v1x+/XvJrNm7OOtcbN2/kXde8q8eY5zc+z62P3srxex/fuEIkSZIkSZLUq2G7tcX37vze1u+3rFL84aIf8sKmFwb93MU32ituGlc6N+/xeaSUYMKeMOmQolEJ5n8Efv92WHhetq9uvezw0n4/5PCe0+hTWxscdlj/H1e8jcUXZn+h4vnBOvRQGNWkd++c++eUNJEBfnD3D3p/wKg2eMUPsg8Syo2bAcffClNfMbiixkyAcbv37zHb174fyE5lPe5Fi/r3Uv1x3QPX1XWcJEmSJEmS6m9YrkjuSl1cufBKACZtO4mPH/VxLrjlAlZvWM3/3fd/nPmyMwf83J1dnUwbP43Ze84mCA6adhB0bYbnH8i2Akib2WX0RN70kpNY05kt1V25fmW2ZcMeZ8Bzd5Y+4WM/GXAtvdr+wOwme/2wyy4wfTo8/nj/XurYY+GPf6x9fGdXJ7986JcATB03lXcf/G4+9atPsXL9Sm7ouIHNXZtpG9XWvyIqmDgRDj4YFi4c9FP1qbhpPHXcVFasX8GP7/kxF518Ue9bLYweD8ffzPrF36Vz6c/Z0DWJteNeTdce72PmxHHUZYOGnY6C9X3caHHUNt3fl6+ar+LlLy897s97oNiaNdnvaNWqbLuMCRNg5sxs/+2xY7M/yz9/8Odbx3/31O+yTVt3zc9vfJ4PX59tE3LdA9fxtRO/RpRvTyNJkiRJkqSGG7pG8oIFPfcrrqZov9ybHr6JR1c/CsBpB5zG+w59HxfccgEA31343UE1kkePGs2vz/p1thXFI9+H378ZVt4Gm0rvPnbtKGD7neElH4Et+/7u8yG45197jK27nQa2fcdxx8EPqiykreTYY+E//7P28bcvu53nXnwOgNft/TraRrVx/F7H86PFP2LVC6uYt3wex8w4pn9FVKmt0Y3kletXcuNDNwJw1PSjOG6P4/iP2/6DVS+s4pcP/TLbH7tISvC738GVV8Ltt0/ggQc+Dny8ZMzYsXDEEfCTn8Buuw2iuJ2O6vuDip2LtoKYcnTNTz1jRnZDw/vuy44XL4ZHHsm2O+nL44/D174GN96YPa7SVtdjx8Lf/A2c9ol5PL3uaQDad2vnfYe+r8fYr/7xq9z/zP0seXYJ9628jwOmHtBjjCRJkiRJkhpr6La2OOKIrMNU/FWsyrUrFl6x9fszX3Yme0/ee2tz8jdLfsOjzz06uNo618Fvj4c7PgBP3ACda2H6m+Cgz0P7JXDIl7Km8dgd4Zk7uh+3zU5w4P/r+/kn7ANjJvY9rjdTBrYtwskn9/8xxx7b95jibRB+8WD39hWv2/t1Jf8tvz5YtdQ2bRD3swP4yT0/obOwv/XbDngbbzvgbVuvlW9v0dkJ73oXvPa1cNVV8MADWaP4Yx+DL38ZvvEN+PSn4YQT4K67YPXqwdVW0wcKM07t/n7y4TCq9rXQr3td6fF551VuCm/eDN8r7DTzy19m+1d/5SvZdhijR8OZZ2aN5R/+EL77XfjMZ7IVzx0dMPfBuYyKUYyKUZz8krI3aOoC4OSXnLx1TPHqZUmSJEmSJDXPsNva4rkXn+Oae68BYPrE6Ry353EAvOtl7+L2ZbeTSHzvzu/x2eM+O/AXufc/YMUt2fdjJ8Prbu19X+INK0uP9/8HWPMAPPzdyuO3PwBm35jdNG+gttsFJh0Kzy2sPi5Kt5A44YS+n3rMmNLjqVNh332zpmhv3vSm7u/ndszd+n2lRvLcjrn882v+ue9CatBXI7mtbWDN82LFzeK3HfA29p68NzO3n8ljax7j2vuvZe3GtUwYOwGACy/MmqVb/O3fZg3k8p8pwLp1ddjjefIREKO2Nlx7iFGwW9GK6dHjYOfXZh+OVLPNFADe8Aa4+OLu0//7v9l/v/KVbMVyVxfcdhv8v/8HGzbAqafCe9+bbWcBsPPO2fW99678MuvXw7hx/5y9H569E5b+CH53Cqy+G158Cro2QYzia2Mm87UjX5b9GTziAzX8YKrbtHkTb/3ft3Lfyvu2nvvC7C9wxsvOGPRzS5IkSZIkjVTDrpF89aKrebGwN/ELnS9w3JVZI/n5Dc9vHfPdhd/l06/+NKNiAJ261AX3Fe3lMOvc0iby6gr/Vr9tPIzeLvt+VBscfTmMnQQd34LNWa2MmQQv+TAc9NnBNZG3mPnWvhvJu7+j5HDatGwh+IIFvT/kjW/see71r++9kTx2bHez9onnn2Dhk1lNo0eN5lO//NTWcW3Rxua0mfnL5/P0uqeZNn6QS4WB3XfP9tq9//7K12fPhh13HPjzL129lFuX3grAobscyj477gPAWw94K1+/4+us37SeOffP4cyXnUlnJ3z+892P3WUX+Ld/q9xEBhg/fuB1bTVmQvaBwrN/rnx9yith26ml56a/qe9G8qyPAtkHD7NnZ1t1bPG//wvXXJM1kp99tntV9ZFHwne+A08+2T32//v/em8iA4wbR/ZnacHH4cFLCnPaAfY4E6a+EradBoyCjSuzD2eeuSP71wLbDOKXCnz6t5/m+geuB+CgaQex6OlF/PV1f81hux7G/lNqvyGhJEmSJElSK8lNI/nyP19Ox6/O23r8tsfnceT0I3uM++7C7pW+q15Yxe+X/r7HmEeee4SbH7mZ1+z1mv4X0rk+a1ZtMfElpddvOAK6NpSeO2khTD6EjRuzPXsXLhzFPfd8lY3rPsvUtoWsWDOZx9fsx9jttmXKlOxmY3/7t4VG2kDNfBvcXWXVdYyG/T7Z4/Rpp/XeSI6A972v5/kPfQguuaTyY048EbbfPvv+ho7uBmVnVyc/vufHFR9zY8eNvOeQ9/Reez986EPwD/9Q+drb3lb5fK2uXnT11u/Xb1rPe3/2XiB7f23xg7t/wJkvO5Pnn89WGW+x116w7baDe/2a7HFG743kfT/W89xubwA+2vvzbTMF9ngXkL0fLr8824ZixYruIZ2d2X7J5YrHQPUm8lYPf6+7iTxqG3j9H2GHxu2BfP0D1/Pl274MwEfbP8pnjvsML/vWy1i5fiVv//HbueNDdzBuTB0+6JEkSZIkSRphctFI/uItX+QzN30GgCN3O5J5y+dx6fxL+fVZv6Z9t/at4xY/vZg/Pf4nAGbtOIt/nl26RcLCJxdubRJ9d+F3tzaS16+HP/8ZHnwQHvqnxJo12bkXXwTena0a3X77bBuHs86awO4T94PnC8tcn7gR9jqr+0WO/i7QBXd8sKSh/D//A5/4BDzzTHZ82mnwlrdM4oADZrP99tlesc8/D088kd3ArKuX3Qhqtv0BsOORsGpe5ev7fAjGz+xx+iMfgS99CZ57rudD3vSmbJVvuUMPhaOOgj/9qee1T3yi+/vibS2qmdsxt26N5Pe+F/7xH2HTptLzkyfDewb5Ev9z9/9s/f6BZx7ggWd6Lsu+8aEbWbl+JTtNmsIBB8C992bnFyyARx+FPfYYXA192uNMuPO8nttbjNsdZry15/jxM2GXE+HJGys/3wHnd6+uB/bZJ7tp3l//dfWV7AcfnN3M8Utf6j734x/Dm9/cR/1P39z9/aSDSpvId30aHvxm6fg93wNHfL2PJ61s6eqlWz8MmLztZN524Nt45LlHOPvws/nX3/8ri55exMfnfpzvnPqdAT2/JEmSJEnSSFZTIzkiTgK+DrQBl6eU/r3sehSunwKsB96XUuplmWSpL//hy1ubyF8/6et8+IgPc8oPTuG3D/+WE75/Ar997285dJdDgdLVyB847AM99jR9035v4pJ5l7Bu0zp+cs9P+MbJ3+DSr+/AP/8zvPBCNubTn4azzoKXvAQmTMgavBs3Zv9M/7HHYLvtgEO/BLe+OXvAoz+E7feH/T4BY7aHPQuvOe+crY3kx5eP5t3v7q7jmGOyLQAies73oIOyrSIGLQIO/TL8tsKq63G7Z3OoYIcd4G/+Br7whdLzbW3w7/9e8SEAnHNOz0by4YdnN5aDbAXyrx76VfZc0caCsxeUrOxcv2k9h192OF2pi18+9Es2d22mbVTpHs4DMXUqvPWt8KMflZ7/6Eez3+9ALXp6EXc9dRcAu0zYheP2OK7k+r0r7+Wup+7KVl4v/jEfOfIjfPvb8OpXZ9c3bsz2DP7KV7KfUfF74bnn4Be/yM7vvPPAawRg3G6w8+t7NoYP+BSM6uWP9yH/UrmRPPnw7H1e5rDDYN48mDMHvvhFuPPOrHE/YUI2349+FE45JZvjm98MP/tZ9rgf/CCb36c+ld10cIuurmzlfkcHvOMVp8LDV2YXnl0IK++AKUdnx3t/CHb7K3jwW/DI97Nzm9fX9nMps2nzJt75k3ey6oVVtEUbU8dP5SM//8jW67tN3I3lzy/nioVXcNyex3HWIWfBoz+CJ38F6x7JttzY5fXZDTVHjYW1DxW2rUmQNkPbdr2+diM0MpMlSf1jJktSvpjLktQ4fTaSI6INuAR4PbAMmBcRc1JK9xQNOxmYVfg6GvhW4b9Vfe2PX+O8X2fbWWw3ejuue+A6rnvgOtZuXAvAsy8+y+uueh03vfcmXrbzy9h3Qjsf3fMbrFoF6/5wGl9fANOnZw2tUaNg1Kjx/MM+/83TLy4jJbj34ee49dYdtjaRx47NVk0efnjpjc7Gjs0aXtOmwebNQNub4IhvwMLzssbV3Z+Fe/4Vxu2R7Tm74RnoXLv18bvu0sn553c3YufNy/bLPe002H//7n1yu7rg6aezFcnHHFOHrQ92np1tQ/Bo98pZ2raDV/wga3r34u/+Dm64obQxfN55Wa29OessmDsXfvKTwsu0wde+1t0g/eNjf2T1hmzD3KOmH8UhuxzS4zmO3O1I7nj8Dla9sIo/Pf4nXj7z5bXOtKqvfz37mS9Zkh3PnJk1Lwfjh3d33zVvyxYIxW559Jat+3P/YNEP+MiRH+FVr4Jrr4Vzz80+lLjzTnjd62CnnbKV3tttB089BQ8/nL0X7rmnDo1kgMO/Br86BjYV7nI3bTa85CO9j9/xCNjnr+Gh/+o+N3YyvPyqXpvPEXDqm4NTgS6C9c93MX58zw9Lvv99+OQn4Yorsj9LX/ta9vuZPj3bN/rFF7P5r12b7cf9jne8GQ7+1+zPWOqE37w6ax5PObZw078ELzzeo55Hn3uU3zz8G37z8G/4w9I/sP+U/Tlo2kHc8ugtrN24ltfs+RqO3/t4XrPna5i83WQuvP1CHn3uUXadsCt/d8zf8Q/Hlu6Hsm7jOmZ/bzaPr3mcC265gNl7zmb37feHF55gc9s4nn9xFU89cz9p7BTGjhnHi+ufYPzobdh1/QOMjS5oa952GI3MZElS/5jJkpQv5rIkNVak8hvHlQ+IeDnw+ZTSiYXj/weQUvq3ojHfBn6XUvph4fh+YHZK6Ynennefl+6TXvhA1uF9y/5v6dGoe37D85z+09N54vknmDZ+GtefeT13/GoGv/51tj/rc8/BW96S7cPa3UjO7t2VUtbEetnLsgbWrbfC739f2NrioayJt+uuWTN0zJhsfGdntor0c5+DQ7b0QF98GpbPhadvgVXzYdNzEKNg9PYwfneYdAjs/BqYdhyMGsMjj8Bvfwt/+UvWJFy9OmuYpZStfI7ItlyYMQO+/e3uvYWLftilx338bgDo2gyLvwCP/AAmHwIH/wtsv1+fD1u7tnuV8bnnZl/FzfVKNmyAD34wm9+FF5aurP7yH77MhbdfCMDfHP03nP/K83s8/l9v/Vcu/tPFAPz9y/+ev3/F3/d8keKfQS3zL+joyLa52GYb+N73smZyVVVeJ6XEq777KpY8m3Wmf3PWbzhgaum+vZu7NnPotw/lmfXPEBEsOHsBu0zYBcjeR7/9Ldx+O9xxR9bgXr8+ew9MnQqzZmUfJHzgAxVuujeQ9wDA8htg4aeyBuwRF0HbNtXHpy5Y/EV4+Puw05FwyL9n7+lq+vG7Wbq0+2ewYAGsWpX9XMaPz343hxySrV5+5SsLD1i/HJZdk61Kfu4u2LACujZl8xi7E4ybCTscCNPfxB3L92LO75fw5BPQ+fyOHDR9b/bZO5g4MfsznWIzj6y7h41dL5K6RvGq/V7Kyw6osGK4j/k89+JzPLTqIZauXsr4seOZ/OoTmbYO2hI8uy2snncrT619ip3G7cQ+k/dh90m7L0gptfd4ojprVCYDtLe3p/nz5zesdklqpohoeC6byZJUm2ZkcuF1GpLLZrKkkWQwmVxLI/k04KSU0ocKx+8Bjk4pnVs05nrg31NKvy8c/wY4L6U0v+y5zgbOLhzuB9xf4SWnACsHMpkRpNV/Bq0+f/BnMBznv0dKaWqjX6SemVy4VpzLBwGLGjyFvBmO77V6cN6tpVXnvV9KaWIjX8BMbohWfL+24pzBebeahmcyNLR/YSa3FufdOlpxzjCITK5lj+QKO/1S3n2uZQwppcuAy6q+WMT8ZnxSmWet/jNo9fmDP4NWn38f6pbJUJrLrfhzb8U5g/Me6jqarZXn3YyXqXDOTB6EVpx3K84ZnPdQ19FsTcpkaFD/opV/b867dbTivFtxzjC4TO5jMwMg21eoeKOAGcDyAYyRJA2emSxJ+WEmS1K+mMuS1EC1NJLnAbMiYq+IGAucDswpGzMHOCsyxwCr+9r3TZI0IGayJOWHmSxJ+WIuS1ID9bm1RUqpMyLOBW4E2oArUkqLI+KcwvVLgbnAKUAHsB54/yBqqrr1RYto9Z9Bq88f/Bm0+vx71eBMbsWfeyvOGZx3q3HeDWImN0QrzrsV5wzOu9U0Zd4NzGV/b63FebeOVpwzDGLefd5sT5IkSZIkSZLU2mrZ2kKSJEmSJEmS1MJsJEuSJEmSJEmSqspVIzkiToqI+yOiIyLOH+p6mi0iHomIuyNiYUTMH+p6miEiroiIpyNiUdG5HSPiVxHxYOG/k4eyxkbr5Wfw+Yh4vPBeWBgRpwxljY0UETMj4qaIuDciFkfEJwrnW+p90Cx95WzhpiMXFa7fFRGHD0Wd9VbDvN9VmO9dEXFbRBwyFHXWW63/X42IIyNic0Sc1sz6GqWWeUfE7EK+Lo6Im5tdY73V8B7fISKui4g7C3MezP0scqPS/0PLruc608xkM7mXcWbyCNCKuTzcMxnM5VbKZTPZTC67bibXmmcppVx8kW2E/xCwNzAWuBM4cKjravLP4BFgylDX0eQ5vxo4HFhUdO7LwPmF788HvjTUdQ7Bz+DzwKeGurYmzX9X4PDC9xOBB4ADW+190KSfdZ85S3bjkV8AARwD3DHUdTdp3q8AJhe+P7lV5l007rdkN545bajrbtLvexJwD7B74XjaUNfdhDn/45YcBaYCq4CxQ117Hebe4/+hZddzm2lmsplsJm8dM6IyuR/zHnG5PJwzuR+/t1zPoYHzHlG5bCabyWbywPMsTyuSjwI6UkpLUkobgauBU4e4JjVYSukWsj+gxU4Fvlf4/nvAm5tZU7P18jNoGSmlJ1JKfy58/zxwLzCdFnsfNEktOXsqcFXK3A5Miohdm11onfU575TSbSmlZwuHtwMzmlxjI9T6/9WPAz8Fnm5mcQ1Uy7zPBK5JKS0FSCkN97nXMucETIyIACaQ/X+ns7ll1l8N/w/Nc6aZyWaymZwZaZkMLZrLwzyTwVxupVw2k81kM3mAeZanRvJ04LGi42WFc60kAb+MiAURcfZQFzOEdk4pPQFZkxGYNsT1DJVzC/+84IpokW0dImJP4DDgDnwfNEItOTsSs7i/c/og2Sezw12f846I6cBbgEubWFej1fL73heYHBG/K/w/96ymVdcYtcz5YuAAYDlwN/CJlFJXc8obUnnONDM5YyYXmMkjJpPBXO5N3jPNXM60Qi6byRkzuZuZXGOejW5YOf0XFc6lplcxtI5NKS2PiGnAryLivsInCGo93wIuIPszcAHwFeADQ1pRg0XEBLJPe/82pbQm+yBQdVZLzo7ELK55ThHxGrK/HL+yoRU1Ry3zvhA4L6W0eQT9matl3qOBI4Djge2AP0bE7SmlBxpdXIPUMucTgYXAa4F9yP6ecWtKaU2Daxtqec40M7mbmZy5EDN5JGQymMu9yXummcvdRnoum8ndzOSMmZzpM8/y1EheBswsOp5B9klAy0gpLS/89+mI+D+y5fet2Eh+KiJ2TSk9UVhWPxL+KUW/pJSe2vJ9RPwXcP0QltNwETGGrIn8PymlawqnW/590AC15OxIzOKa5hQRBwOXAyenlJ5pUm2NVMu824GrC385ngKcEhGdKaWfNaXCxqj1fb4ypbQOWBcRtwCHkO3RPhzVMuf3A/+esg3ROiLiYWB/4E/NKXHI5DnTzOSMmdzNTB4ZmQzmcm/ynmnmcqYVctlMzpjJ3czkGvMsT1tbzANmRcReETEWOB2YM8Q1NU1EjI+IiVu+B04AKt5ZsQXMAd5b+P69wLVDWMuQKNuX5i2M4PdCYQ+i7wD3ppS+WnSp5d8HDVBLzs4BzircwfUYYPWWLUaGsT7nHRG7A9cA7xnmn7YX63PeKaW9Ukp7ppT2BH4CfHSY/+UYanufXwu8KiJGR8Q44Giy/dmHq1rmvJRsZQkRsTOwH7CkqVUOjTxnmplsJpvJmZGWyWAu9ybvmWYut04um8lmspk8wDzLzYrklFJnRJwL3Eh2R8UrUkqLh7isZtoZ+L/Cp12jgR+klG4Y2pIaLyJ+CMwGpkTEMuBzwL8D/xsRHyT7w/z2oauw8Xr5GcyOiEPJ/lnBI8CHh6q+JjgWeA9wd0QsLJz7R1rsfdAMveVsRJxTuH4p2R2JTwE6gPVkn8wOazXO+7PATsA3CzncmVJqH6qa66HGeY84tcw7pXRvRNwA3AV0AZenlIbtB3Y1/q4vAK6MiLvJ/hnbeSmllUNWdJ308v/QMZD/TDOTzWQzeWRmMrRuLg/nTAZzmRbKZTPZTDaTB55nka3aliRJkiRJkiSpsjxtbSFJkiRJkiRJyiEbyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSaqqz0ZyRFwREU9HxKJerkdEXBQRHRFxV0QcXv8yJUlbmMuSlB9msiTlh5ksSY1Vy4rkK4GTqlw/GZhV+Dob+Nbgy5IkVXEl5rIk5cWVmMmSlBdXYiZLUsP02UhOKd0CrKoy5FTgqpS5HZgUEbvWq0BJUilzWZLyw0yWpPwwkyWpsUbX4TmmA48VHS8rnHuifGBEnE32qR/jx48/Yv/996/Dy0vS0FuwYMHKlNLUoa6jwFyW1PJylMtmsqSWZyZLUn4MJpPr0UiOCudSpYEppcuAywDa29vT/Pnz6/DykjT0IuLRoa6hiLksqeXlKJfNZEktz0yWpPwYTCbXskdyX5YBM4uOZwDL6/C8kqSBMZclKT/MZEnKDzNZkgahHo3kOcBZhbufHgOsTin1+GchkqSmMZclKT/MZEnKDzNZkgahz60tIuKHwGxgSkQsAz4HjAFIKV0KzAVOATqA9cD7G1WsJMlclqQ8MZMlKT/MZElqrD4bySmlM/q4noCP1a0iSVJV5rIk5YeZLEn5YSZLUmPVY2sLSZIkSZIkSdIIZiNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdXUSI6IkyLi/ojoiIjzK1zfISKui4g7I2JxRLy//qVKksBMlqQ8MZMlKV/MZUlqnD4byRHRBlwCnAwcCJwREQeWDfsYcE9K6RBgNvCViBhb51olqeWZyZKUH2ayJOWLuSxJjVXLiuSjgI6U0pKU0kbgauDUsjEJmBgRAUwAVgGdda1UkgRmsiTliZksSfliLktSA9XSSJ4OPFZ0vKxwrtjFwAHAcuBu4BMppa7yJ4qIsyNifkTMX7FixQBLlqSWVrdMBnNZkgbJTJakfLF/IUkNVEsjOSqcS2XHJwILgd2AQ4GLI2L7Hg9K6bKUUntKqX3q1Kn9LFWSRB0zGcxlSRokM1mS8sX+hSQ1UC2N5GXAzKLjGWSf3BV7P3BNynQADwP716dESVIRM1mS8sNMlqR8MZclqYFqaSTPA2ZFxF6FDehPB+aUjVkKHA8QETsD+wFL6lmoJAkwkyUpT8xkScoXc1mSGmh0XwNSSp0RcS5wI9AGXJFSWhwR5xSuXwpcAFwZEXeT/VOS81JKKxtYtyS1JDNZkvLDTJakfDGXJamx+mwkA6SU5gJzy85dWvT9cuCE+pYmSarETJak/DCTJSlfzGVJapxatraQJEmSJEmSJLUwG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqqqZGckScFBH3R0RHRJzfy5jZEbEwIhZHxM31LVOStIWZLEn5YSZLUr6Yy5LUOKP7GhARbcAlwOuBZcC8iJiTUrqnaMwk4JvASSmlpRExrUH1SlJLM5MlKT/MZEnKF3NZkhqrlhXJRwEdKaUlKaWNwNXAqWVjzgSuSSktBUgpPV3fMiVJBWayJOWHmSxJ+WIuS1ID1dJIng48VnS8rHCu2L7A5Ij4XUQsiIizKj1RRJwdEfMjYv6KFSsGVrEktba6ZTKYy5I0SGayJOWL/QtJaqBaGslR4VwqOx4NHAH8FXAi8JmI2LfHg1K6LKXUnlJqnzp1ar+LlSTVL5PBXJakQTKTJSlf7F9IUgP1uUcy2Sd4M4uOZwDLK4xZmVJaB6yLiFuAQ4AH6lKlJGkLM1mS8sNMlqR8MZclqYFqWZE8D5gVEXtFxFjgdGBO2ZhrgVdFxOiIGAccDdxb31IlSZjJkpQnZrIk5Yu5LEkN1OeK5JRSZ0ScC9wItAFXpJQWR8Q5heuXppTujYgbgLuALuDylNKiRhYuSa3ITJak/DCTJSlfzGVJaqxIqXy7oOZob29P8+fPH5LXlqR6i4gFKaX2oa5jMMxlSSPJcM9lM1nSSGImS1J+DCaTa9naQpIkSZIkSZLUwmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqaqaGskRcVJE3B8RHRFxfpVxR0bE5og4rX4lSpKKmcmSlB9msiTli7ksSY3TZyM5ItqAS4CTgQOBMyLiwF7GfQm4sd5FSpIyZrIk5YeZLEn5Yi5LUmPVsiL5KKAjpbQkpbQRuBo4tcK4jwM/BZ6uY32SpFJmsiTlh5ksSfliLktSA9XSSJ4OPFZ0vKxwbquImA68Bbi0fqVJkiowkyUpP8xkScoXc1mSGqiWRnJUOJfKji8Ezkspba76RBFnR8T8iJi/YsWKGkuUJBWpWyaDuSxJg2QmS1K+2L+QpAYaXcOYZcDMouMZwPKyMe3A1REBMAU4JSI6U0o/Kx6UUroMuAygvb29PMwlSX2rWyaDuSxJg2QmS1K+2L+QpAaqpZE8D5gVEXsBjwOnA2cWD0gp7bXl+4i4Eri+0l+OJUmDZiZLUn6YyZKUL+ayJDVQn43klFJnRJxLdjfTNuCKlNLiiDincN19hSSpScxkScoPM1mS8sVclqTGqmVFMimlucDcsnMVAzil9L7BlyVJ6o2ZLEn5YSZLUr6Yy5LUOLXcbE+SJEmSJEmS1MJsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpqaiRHxEkRcX9EdETE+RWuvysi7ip83RYRh9S/VEkSmMmSlCdmsiTli7ksSY3TZyM5ItqAS4CTgQOBMyLiwLJhDwPHpZQOBi4ALqt3oZIkM1mS8sRMlqR8MZclqbFqWZF8FNCRUlqSUtoIXA2cWjwgpXRbSunZwuHtwIz6lilJKjCTJSk/zGRJyhdzWZIaqJZG8nTgsaLjZYVzvfkg8ItKFyLi7IiYHxHzV6xYUXuVkqQt6pbJYC5L0iCZyZKUL/YvJKmBamkkR4VzqeLAiNeQBfF5la6nlC5LKbWnlNqnTp1ae5WSpC3qlslgLkvSIJnJkpQv9i8kqYFG1zBmGTCz6HgGsLx8UEQcDFwOnJxSeqY+5UmSypjJkpQfZrIk5Yu5LEkNVMuK5HnArIjYKyLGAqcDc4oHRMTuwDXAe1JKD9S/TElSgZksSflhJktSvpjLktRAfa5ITil1RsS5wI1AG3BFSmlxRJxTuH4p8FlgJ+CbEQHQmVJqb1zZktSazGRJyg8zWZLyxVyWpMaKlCpuF9Rw7e3taf78+UPy2pJUbxGxYLj/BdRcljSSDPdcNpMljSRmsiTlx2AyuZatLSRJkiRJkiRJLcxGsiRJkiRJkiSpKhvJkiRJkiRJkqSqbCRLkiRJkiRJkqqykSxJkiRJkiRJqspGsiRJkiRJkiSpKhvJkiRJkiRJkqSqbCRLkiRJkiRJkqqykSxJkiRJkiRJqspGsiRJkiRJkiSpKhvJkiRJkiRJkqSqbCRLkiRJkiRJkqqqqZEcESdFxP0R0RER51e4HhFxUeH6XRFxeP1LlSSBmSxJeWImS1K+mMuS1Dh9NpIjog24BDgZOBA4IyIOLBt2MjCr8HU28K061ylJwkyWpDwxkyUpX8xlSWqsWlYkHwV0pJSWpJQ2AlcDp5aNORW4KmVuByZFxK51rlWSZCZLUp6YyZKUL+ayJDXQ6BrGTAceKzpeBhxdw5jpwBPFgyLibLJP/AA2RMSiflU7MkwBVg51EUPAebeWVpz3fk16nbplMpjLtOZ7FZx3q2nVeTcjl83k+mvF92srzhmcd6sZdn9XNpOB1n2/Ou/W0YpzhkFkci2N5KhwLg1gDCmly4DLACJifkqpvYbXH1Gcd2tx3q0jIuY366UqnBtQJoO53IpzBuc91HU0WyvPuxkvU+GcmTwIrTjvVpwzOO+hrqPZhuPflc1k5z3UdTRbK867FecMg8vkWra2WAbMLDqeASwfwBhJ0uCZyZKUH2ayJOWLuSxJDVRLI3keMCsi9oqIscDpwJyyMXOAswp3Pz0GWJ1S6vHP9SRJg2YmS1J+mMmSlC/msiQ1UJ9bW6SUOiPiXOBGoA24IqW0OCLOKVy/FJgLnAJ0AOuB99fw2pcNuOrhzXm3FufdOpoy5wZmMvh7ayXOu7U47wYxkxuiFefdinMG591qhvvflf29tRbn3Tpacc4wiHlHShW3aJMkSZIkSZIkCahtawtJkiRJkiRJUguzkSxJkiRJkiRJqqrhjeSIOCki7o+Ijog4v8L1iIiLCtfviojDG11TM9Qw73cV5ntXRNwWEYcMRZ311Neci8YdGRGbI+K0ZtbXKLXMOyJmR8TCiFgcETc3u8ZGqOE9vkNEXBcRdxbmXeuekLkVEVdExNMRsaiX67nPMzO5dTIZzOVWyuVWzGQY/rlsJpvJvYwzk0eAVszl4Z7JYC63Ui6byWZy2XUzudY8Syk17Itsc/uHgL2BscCdwIFlY04BfgEEcAxwRyNrasZXjfN+BTC58P3Jw33etcy5aNxvyW5wcNpQ192k3/Uk4B5g98LxtKGuu0nz/kfgS4XvpwKrgLFDXfsg5/1q4HBgUS/Xc51nZnLrZHKt8y4aZy7noPYGz3nEZXJhLsM2l81kM9lM3jpmRGVyP+Y94nJ5OGdyP35vuZ5DA+c9onLZTDaTzeSB51mjVyQfBXSklJaklDYCVwOnlo05FbgqZW4HJkXErg2uq9H6nHdK6baU0rOFw9uBGU2usd5q+V0DfBz4KfB0M4troFrmfSZwTUppKUBKaSTMvZZ5J2BiRAQwgSyIO5tbZn2llG4hm0dv8p5nZnLrZDKYy62Uyy2ZyTDsc9lMNpPN5MxIy2Ro0Vwe5pkM5nIr5bKZbCabyQPMs0Y3kqcDjxUdLyuc6++Y4aa/c/og2acAw1mfc46I6cBbgEubWFej1fK73heYHBG/i4gFEXFW06prnFrmfTFwALAcuBv4REqpqznlDZm855mZnGmFTAZzeYtWyGUzuXd5zjQzOWMmF5jJIyaTwVzuTd4zzVzOtEIum8kZM7mbmVxjno1uWDmZqHAuDWDMcFPznCLiNWRB/MqGVtR4tcz5QuC8lNLm7EOeEaGWeY8GjgCOB7YD/hgRt6eUHmh0cQ1Uy7xPBBYCrwX2AX4VEbemlNY0uLahlPc8M5O7jfRMBnO52EjPZTO5d3nONDO5m5mcuRAzeSRkMpjLvcl7ppnL3UZ6LpvJ3czkjJmc6TPPGt1IXgbMLDqeQdbd7++Y4aamOUXEwcDlwMkppWeaVFuj1DLnduDqQghPAU6JiM6U0s+aUmFj1PoeX5lSWgesi4hbgEOA4RzEtcz7/cC/p5QS0BERDwP7A39qTolDIu95ZiZnWiGTwVzeohVy2UzuXZ4zzUzOmMndzOSRkclgLvcm75lmLmdaIZfN5IyZ3M1MrjHPGr21xTxgVkTsFRFjgdOBOWVj5gBnFe4WeAywOqX0RIPrarQ+5x0RuwPXAO8Z5p/sbNHnnFNKe6WU9kwp7Qn8BPjoMA9hqO09fi3wqogYHRHjgKOBe5tcZ73VMu+lZJ9iEhE7A/sBS5paZfPlPc/M5NbJZDCXWymXzeTe5TnTzGQz2UzOjLRMBnO5N3nPNHO5dXLZTDaTzeQB5llDVySnlDoj4lzgRrK7JF6RUlocEecUrl9KdvfLU4AOYD3ZpwDDWo3z/iywE/DNwidcnSml9qGqebBqnPOIU8u8U0r3RsQNwF1AF3B5SmnR0FU9eDX+vi8AroyIu8n+ycR5KaWVQ1Z0HUTED4HZwJSIWAZ8DhgDwyPPzOTWyWQwl2mhXG7VTIbhnctmsplsJo/MTIbWzeXhnMlgLtNCuWwmm8lm8sDzLLJV25IkSZIkSZIkVdborS0kSZIkSZIkScOcjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRV1WcjOSKuiIinI2JRL9cjIi6KiI6IuCsiDq9/mZKkLcxlScoPM1mS8sNMlqTGqmVF8pXASVWunwzMKnydDXxr8GVJkqq4EnNZkvLiSsxkScqLKzGTJalh+mwkp5RuAVZVGXIqcFXK3A5Miohd61WgJKmUuSxJ+WEmS1J+mMmS1Fij6/Ac04HHio6XFc49UT4wIs4m+9SP8ePHH7H//vvX4eUlaegtWLBgZUpp6lDXUWAuS2p5OcplM1lSyzOTJSk/BpPJ9WgkR4VzqdLAlNJlwGUA7e3taf78+XV4eUkaehHx6FDXUMRcltTycpTLZrKklmcmS1J+DCaTa9kjuS/LgJlFxzOA5XV4XknSwJjLkpQfZrIk5YeZLEmDUI9G8hzgrMLdT48BVqeUevyzEElS05jLkpQfZrIk5YeZLEmD0OfWFhHxQ2A2MCUilgGfA8YApJQuBeYCpwAdwHrg/Y0qVpJkLktSnpjJkpQfZrIkNVafjeSU0hl9XE/Ax+pWkSSpKnNZkvLDTJak/DCTJamx6rG1hSRJkiRJkiRpBLORLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpqaiRHxEkRcX9EdETE+RWu7xAR10XEnRGxOCLeX/9SJUlgJktSnpjJkpQv5rIkNU6fjeSIaAMuAU4GDgTOiIgDy4Z9DLgnpXQIMBv4SkSMrXOtktTyzGRJyg8zWZLyxVyWpMaqZUXyUUBHSmlJSmkjcDVwatmYBEyMiAAmAKuAzrpWKkkCM1mS8sRMlqR8MZclqYFqaSRPBx4rOl5WOFfsYuAAYDlwN/CJlFJX+RNFxNkRMT8i5q9YsWKAJUtSS6tbJoO5LEmDZCZLUr7Yv5CkBqqlkRwVzqWy4xOBhcBuwKHAxRGxfY8HpXRZSqk9pdQ+derUfpYqSaKOmQzmsiQNkpksSfli/0KSGqiWRvIyYGbR8QyyT+6KvR+4JmU6gIeB/etToiSpiJksSflhJktSvpjLktRAtTSS5wGzImKvwgb0pwNzysYsBY4HiIidgf2AJfUsVJIEmMmSlCdmsiTli7ksSQ00uq8BKaXOiDgXuBFoA65IKS2OiHMK1y8FLgCujIi7yf4pyXkppZUNrFuSWpKZLEn5YSZLUr6Yy5LUWH02kgFSSnOBuWXnLi36fjlwQn1LkyRVYiZLUn6YyZKUL+ayJDVOLVtbSJIkSZIkSZJamI1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVWUjWZIkSZIkSZJUlY1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVWUjWZIkSZIkSZJUlY1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVVVTIzkiToqI+yOiIyLO72XM7IhYGBGLI+Lm+pYpSdrCTJak/DCTJSlfzGVJapzRfQ2IiDbgEuD1wDJgXkTMSSndUzRmEvBN4KSU0tKImNageiWppZnJkpQfZrIk5Yu5LEmNVcuK5KOAjpTSkpTSRuBq4NSyMWcC16SUlgKklJ6ub5mSpAIzWZLyw0yWpHwxlyWpgWppJE8HHis6XlY4V2xfYHJE/C4iFkTEWfUqUJJUwkyWpPwwkyUpX8xlSWqgPre2AKLCuVTheY4Ajge2A/4YEbenlB4oeaKIs4GzAXbffff+VytJqlsmg7ksSYNkJktSvti/kKQGqmVF8jJgZtHxDGB5hTE3pJTWpZRWArcAh5Q/UUrpspRSe0qpferUqQOtWZJaWd0yGcxlSRokM1mS8sX+hSQ1UC2N5HnArIjYKyLGAqcDc8rGXAu8KiJGR8Q44Gjg3vqWKknCTJakPDGTJSlfzGVJaqA+t7ZIKXVGxLnAjUAbcEVKaXFEnFO4fmlK6d6IuAG4C+gCLk8pLWpk4ZLUisxkScoPM1mS8sVclqTGipTKtwtqjvb29jR//vwheW1JqreIWJBSah/qOgbDXJY0kgz3XDaTJY0kZrIk5cdgMrmWrS0kSZIkSZIkSS3MRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSaqqpkZyRJwUEfdHREdEnF9l3JERsTkiTqtfiZKkYmayJOWHmSxJ+WIuS1Lj9NlIjog24BLgZOBA4IyIOLCXcV8Cbqx3kZKkjJksSflhJktSvpjLktRYtaxIPgroSCktSSltBK4GTq0w7uPAT4Gn61ifJKmUmSxJ+WEmS1K+mMuS1EC1NJKnA48VHS8rnNsqIqYDbwEurfZEEXF2RMyPiPkrVqzob62SpDpmcmGsuSxJA2cmS1K+2L+QpAaqpZEcFc6lsuMLgfNSSpurPVFK6bKUUntKqX3q1Kk1lihJKlK3TAZzWZIGyUyWpHyxfyFJDTS6hjHLgJlFxzOA5WVj2oGrIwJgCnBKRHSmlH5WjyIlSVuZyZKUH2ayJOWLuSxJDVRLI3keMCsi9gIeB04HziwekFLaa8v3EXElcL0hLEkNYSZLUn6YyZKUL+ayJDVQn43klFJnRJxLdjfTNuCKlNLiiDincL3P/d4kSfVhJktSfpjJkpQv5rIkNVYtK5JJKc0F5padqxjAKaX3Db4sSVJvzGRJyg8zWZLyxVyWpMap5WZ7kiRJkiRJkqQWZiNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdXUSI6IkyLi/ojoiIjzK1x/V0TcVfi6LSIOqX+pkiQwkyUpT8xkScoXc1mSGqfPRnJEtAGXACcDBwJnRMSBZcMeBo5LKR0MXABcVu9CJUlmsiTliZksSfliLktSY9WyIvkooCOltCSltBG4Gji1eEBK6baU0rOFw9uBGfUtU5JUYCZLUn6YyZKUL+ayJDVQLY3k6cBjRcfLCud680HgF5UuRMTZETE/IuavWLGi9iolSVvULZPBXJakQTKTJSlf7F9IUgPV0kiOCudSxYERryEL4vMqXU8pXZZSak8ptU+dOrX2KiVJW9Qtk8FclqRBMpMlKV/sX0hSA42uYcwyYGbR8QxgefmgiDgYuBw4OaX0TH3KkySVMZMlKT/MZEnKF3NZkhqolhXJ84BZEbFXRIwFTgfmFA+IiN2Ba4D3pJQeqH+ZkqQCM1mS8sNMlqR8MZclqYH6XJGcUuqMiHOBG4E24IqU0uKIOKdw/VLgs8BOwDcjAqAzpdTeuLIlqTWZyZKUH2ayJOWLuSxJjRUpVdwuqOHa29vT/Pnzh+S1JaneImLBcP8LqLksaSQZ7rlsJksaScxkScqPwWRyLVtbSJIkSZIkSZJamI1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVWUjWZIkSZIkSZJUlY1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVWUjWZIkSZIkSZJUlY1kSZIkSZIkSVJVNpIlSZIkSZIkSVXZSJYkSZIkSZIkVVVTIzkiToqI+yOiIyLOr3A9IuKiwvW7IuLw+pcqSQIzWZLyxEyWpHwxlyWpcfpsJEdEG3AJcDJwIHBGRBxYNuxkYFbh62zgW3WuU5KEmSxJeWImS1K+mMuS1Fi1rEg+CuhIKS1JKW0ErgZOLRtzKnBVytwOTIqIXetcqyTJTJakPDGTJSlfzGVJaqBaGsnTgceKjpcVzvV3jCRp8MxkScoPM1mS8sVclqQGGl3DmKhwLg1gDBFxNtk/HQHYEBGLanj9kWYKsHKoixgCzru1tOK892vS69Qtk8FcpjXfq+C8W02rzrsZuWwm118rvl9bcc7gvFvNsPu7spkMtO771Xm3jlacMwwik2tpJC8DZhYdzwCWD2AMKaXLgMsAImJ+Sqm9X9WOAM67tTjv1hER85v0UnXLZDCXW3HO4LyHuo5ma+V5N+FlzOQ6a8V5t+KcwXkPdR3NNhz/rmwmO++hrqPZWnHerThnGFwm17K1xTxgVkTsFRFjgdOBOWVj5gBnFe5+egywOqX0xECLkiT1ykyWpPwwkyUpX8xlSWqgPlckp5Q6I+Jc4EagDbgipbQ4Is4pXL8UmAucAnQA64H3N65kSWpdZrIk5YeZLEn5Yi5LUmPVsrUFKaW5ZGFbfO7Sou8T8LF+vvZl/Rw/Ujjv1uK8W0fT5tygTAZ/b63EebcW591AZnLdteK8W3HO4LxbzXD/u7K/t9bivFtHK84ZBjHvyDJUkiRJkiRJkqTKatkjWZIkSZIkSZLUwhreSI6IkyLi/ojoiIjzK1yPiLiocP2uiDi80TU1Qw3zfldhvndFxG0RcchQ1FlPfc25aNyREbE5Ik5rZn2NUsu8I2J2RCyMiMURcXOza2yEGt7jO0TEdRFxZ2Hew37vsYi4IiKejohFvVzPfZ6Zya2TyWAut1Iut2Imw/DPZTPZTO5lnJk8ArRiLg/3TAZzuZVy2Uw2k8uum8m15llKqWFfZJvbPwTsDYwF7gQOLBtzCvALIIBjgDsaWVMzvmqc9yuAyYXvTx7u865lzkXjfku2Z9VpQ113k37Xk4B7gN0Lx9OGuu4mzfsfgS8Vvp8KrALGDnXtg5z3q4HDgUW9XM91npnJrZPJtc67aJy5nIPaGzznEZfJhbkM21w2k81kM3nrmBGVyf2Y94jL5eGcyf34veV6Dg2c94jKZTPZTDaTB55njV6RfBTQkVJaklLaCFwNnFo25lTgqpS5HZgUEbs2uK5G63PeKaXbUkrPFg5vB2Y0ucZ6q+V3DfBx4KfA080sroFqmfeZwDUppaUAKaWRMPda5p2AiRERwASyIO5sbpn1lVK6hWwevcl7npnJrZPJYC63Ui63ZCbDsM9lM9lMNpMzIy2ToUVzeZhnMpjLrZTLZrKZbCYPMM8a3UieDjxWdLyscK6/Y4ab/s7pg2SfAgxnfc45IqYDbwEuZeSo5Xe9LzA5In4XEQsi4qymVdc4tcz7YuAAYDlwN/CJlFJXc8obMnnPMzM50wqZDObyFq2Qy2Zy7/KcaWZyxkwuMJNHTCaDudybvGeauZxphVw2kzNmcjczucY8G92wcjJR4VwawJjhpuY5RcRryIL4lQ2tqPFqmfOFwHkppc3ZhzwjQi3zHg0cARwPbAf8MSJuTyk90OjiGqiWeZ8ILAReC+wD/Coibk0prWlwbUMp73lmJncb6ZkM5nKxkZ7LZnLv8pxpZnI3MzlzIWbySMhkMJd7k/dMM5e7jfRcNpO7mckZMznTZ541upG8DJhZdDyDrLvf3zHDTU1zioiDgcuBk1NKzzSptkapZc7twNWFEJ4CnBIRnSmlnzWlwsao9T2+MqW0DlgXEbcAhwDDOYhrmff7gX9PKSWgIyIeBvYH/tScEodE3vPMTM60QiaDubxFK+Symdy7PGeamZwxk7uZySMjk8Fc7k3eM81czrRCLpvJGTO5m5lcY541emuLecCsiNgrIsYCpwNzysbMAc4q3C3wGGB1SumJBtfVaH3OOyJ2B64B3jPMP9nZos85p5T2SintmVLaE/gJ8NFhHsJQ23v8WuBVETE6IsYBRwP3NrnOeqtl3kvJPsUkInYG9gOWNLXK5st7npnJrZPJYC63Ui6byb3Lc6aZyWaymZwZaZkM5nJv8p5p5nLr5LKZbCabyQPMs4auSE4pdUbEucCNZHdJvCKltDgizilcv5Ts7penAB3AerJPAYa1Guf9WWAn4JuFT7g6U0rtQ1XzYNU45xGnlnmnlO6NiBuAu4Au4PKU0qKhq3rwavx9XwBcGRF3k/2TifNSSiuHrOg6iIgfArOBKRGxDPgcMAaGR56Zya2TyWAu00K53KqZDMM7l81kM9lMHpmZDK2by8M5k8FcpoVy2Uw2k83kgedZZKu2JUmSJEmSJEmqrNFbW0iSJEmSJEmShjkbyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSarKRrIkSZIkSZIkqSobyZIkSZIkSZKkqmwkS5IkSZIkSZKqspEsSZIkSZIkSaqqz0ZyRFwREU9HxKJerkdEXBQRHRFxV0QcXv8yJUlbmMuSlB9msiTlh5ksSY1Vy4rkK4GTqlw/GZhV+Dob+Nbgy5IkVXEl5rIk5cWVmMmSlBdXYiZLUsP02UhOKd0CrKoy5FTgqpS5HZgUEbvWq0BJUilzWZLyw0yWpPwwkyWpsUbX4TmmA48VHS8rnHuifGBEnE32qR/jx48/Yv/996/Dy0vS0FuwYMHKlNLUoa6jwFyW1PJylMtmsqSWZyZLUn4MJpPr0UiOCudSpYEppcuAywDa29vT/Pnz6/DykjT0IuLRoa6hiLksqeXlKJfNZEktz0yWpPwYTCbXskdyX5YBM4uOZwDL6/C8kqSBMZclKT/MZEnKDzNZkgahHo3kOcBZhbufHgOsTin1+GchkqSmMZclKT/MZEnKDzNZkgahz60tIuKHwGxgSkQsAz4HjAFIKV0KzAVOATqA9cD7G1WsJMlclqQ8MZMlKT/MZElqrD4bySmlM/q4noCP1a0iSVJV5rIk5YeZLEn5YSZLUmPVY2sLSZIkSZIkSdIIZiNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdlIliRJkiRJkiRVZSNZkiRJkiRJklSVjWRJkiRJkiRJUlU2kiVJkiRJkiRJVdXUSI6IkyLi/ojoiIjzK1zfISKui4g7I2JxRLy//qVKksBMlqQ8MZMlKV/MZUlqnD4byRHRBlwCnAwcCJwREQeWDfsYcE9K6RBgNvCViBhb51olqeWZyZKUH2ayJOWLuSxJjVXLiuSjgI6U0pKU0kbgauDUsjEJmBgRAUwAVgGdda1UkgRmsiTliZksSfliLktSA9XSSJ4OPFZ0vKxwrtjFwAHAcuBu4BMppa66VChJKmYmS1J+mMmSlC/msiQ1UC2N5KhwLpUdnwgsBHYDDgUujojtezxRxNkRMT8i5q9YsaKfpUqSqGMmg7ksSYNkJktSvti/kKQGqqWRvAyYWXQ8g+yTu2LvB65JmQ7gYWD/8idKKV2WUmpPKbVPnTp1oDVLUiurWyaDuSxJg2QmS1K+2L+QpAaqpZE8D5gVEXsVNqA/HZhTNmYpcDxAROwM7AcsqWehkiTATJakPDGTJSlfzGVJaqDRfQ1IKXVGxLnAjUAbcEVKaXFEnFO4filwAXBlRNxN9k9JzksprWxg3ZLUksxkScoPM1mS8sVclqTG6rORDJBSmgvMLTt3adH3y4ET6luaJKkSM1mS8sNMlqR8MZclqXFq2dpCkiRJkiRJktTCbCRLkiRJkiRJkqqykSxJkiRJkiRJqspGsiRJkiRJkiSpKhvJkiRJkiRJkqSqbCRLkiRJkiRJkqqykSxJkiRJkiRJqspGsiRJkiRJkiSpKhvJkiRJkiRJkqSqbCRLkiRJkiRJkqqykSxJkiRJkiRJqspGsiRJkiRJkiSpKhvJkiRJkiRJkqSqamokR8RJEXF/RHRExPm9jJkdEQsjYnFE3FzfMiVJW5jJkpQfZrIk5Yu5LEmNM7qvARHRBlwCvB5YBsyLiDkppXuKxkwCvgmclFJaGhHTGlSvJLU0M1mS8sNMlqR8MZclqbFqWZF8FNCRUlqSUtoIXA2cWjbmTOCalNJSgJTS0/UtU5JUYCZLUn6YyZKUL+ayJDVQLY3k6cBjRcfLCueK7QtMjojfRcSCiDir0hNFxNkRMT8i5q9YsWJgFUtSa6tbJoO5LEmDZCZLUr7Yv5CkBqqlkRwVzqWy49HAEcBfAScCn4mIfXs8KKXLUkrtKaX2qVOn9rtYSVL9MhnMZUkaJDNZkvLF/oUkNVCfeySTfYI3s+h4BrC8wpiVKaV1wLqIuAU4BHigLlVKkrYwkyUpP8xkScoXc1mSGqiWFcnzgFkRsVdEjAVOB+aUjbkWeFVEjI6IccDRwL31LVWShJksSXliJktSvpjLktRAfa5ITil1RsS5wI1AG3BFSmlxRJxTuH5pSuneiLgBuAvoAi5PKS1qZOGS1IrMZEnKDzNZkvLFXJakxoqUyrcLao729vY0f/78IXltSaq3iFiQUmof6joGw1yWNJIM91w2kyWNJGayJOXHYDK5lq0tJEmSJEmSJEktzEayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqrKRLEmSJEmSJEmqykayJEmSJEmSJKkqG8mSJEmSJEmSpKpsJEuSJEmSJEmSqqqpkRwRJ0XE/RHRERHnVxl3ZERsjojT6leiJKmYmSxJ+WEmS1K+mMuS1Dh9NpIjog24BDgZOBA4IyIO7GXcl4Ab612kJCljJktSfpjJkpQv5rIkNVYtK5KPAjpSSktSShuBq4FTK4z7OPBT4Ok61idJKmUmS1J+mMmSlC/msiQ1UC2N5OnAY0XHywrntoqI6cBbgEurPVFEnB0R8yNi/ooVK/pbqySpjplcGGsuS9LAmcmSlC/2LySpgWppJEeFc6ns+ELgvJTS5mpPlFK6LKXUnlJqnzp1ao0lSpKK1C2TwVyWpEEykyUpX+xfSFIDja5hzDJgZtHxDGB52Zh24OqIAJgCnBIRnSml/7+9OwqV9D7rOP57yBpQFFOaKGU3S1dZq3uRQFzTIFVbvTDJzSL0Im1psBRCkIiXDV7Ui97oXSm2hiWE0htzoaFGSQ2CaIUYzRbaNLGkrCk0awoltVRowbDt48VMuof1nDnvnnPeOe/M+/nAgTPz/tnzf5jDl8N/Zmc+fxSbBODHNBlgOjQZYFp0GWBEQw6SX0hytqrOJPmvJA8k+eDOBd195q3vq+qzSf5OhAFGockA06HJANOiywAj2vcgubuvVtUjWXya6U1Jnujul6vq4eX1fd/vDYCjockA06HJANOiywDjGvKK5HT3M0meue6+XQPc3b9/+G0BsBdNBpgOTQaYFl0GGM+QD9sDAAAAAGDGHCQDAAAAALCSg2QAAAAAAFZykAwAAAAAwEoOkgEAAAAAWMlBMgAAAAAAKzlIBgAAAABgJQfJAAAAAACs5CAZAAAAAICVHCQDAAAAALCSg2QAAAAAAFZykAwAAAAAwEqDDpKr6t6qeqWqLlfVo7tc/1BVvbj8eq6q7jz6rQKQaDLAlGgywLToMsB49j1Irqqbknw6yX1JziX5QFWdu27ZN5L8VnffkeQTSS4e9UYB0GSAKdFkgGnRZYBxDXlF8t1JLnf3q939ZpInk1zYuaC7n+vu7y5vPp/k1NFuE4AlTQaYDk0GmBZdBhjRkIPkk0le23H7yvK+vXw0yRcOsykA9qTJANOhyQDTossAIzoxYE3tcl/vurDqfVmE+D17XH8oyUNJcvr06YFbBGCHI2vyco0uAxycJgNMi/MLgBENeUXylSS377h9Ksnr1y+qqjuSPJ7kQnd/Z7d/qLsvdvf57j5/2223HWS/AHN3ZE1OdBngkDQZYFqcXwCMaMhB8gtJzlbVmaq6OckDSZ7euaCqTid5KsmHu/vrR79NAJY0GWA6NBlgWnQZYET7vrVFd1+tqkeSPJvkpiRPdPfLVfXw8vpjST6e5O1JPlNVSXK1u8+Pt22AedJkgOnQZIBp0WWAcVX3rm8XNLrz58/3pUuXjuVnAxy1qvrSpv8BqsvANtn0LmsysE00GWA6DtPkIW9tAQAAAADAjDlIBgAAAABgJQfJAAAAAACs5CAZAAAAAICVHCQDAAAAALCSg2QAAAAAAFZykAwAAAAAwEoOkgEAAAAAWMlBMgAAAAAAKzlIBgAAAABgJQfJAAAAAACs5CAZAAAAAICVHCQDAAAAALDSoIPkqrq3ql6pqstV9egu16uqPrW8/mJV3XX0WwUg0WSAKdFkgGnRZYDx7HuQXFU3Jfl0kvuSnEvygao6d92y+5KcXX49lOQvjnifAESTAaZEkwGmRZcBxjXkFcl3J7nc3a9295tJnkxy4bo1F5J8rheeT3JLVb3jiPcKgCYDTIkmA0yLLgOM6MSANSeTvLbj9pUk7x6w5mSSb+1cVFUPZfGMX5L8b1W9dEO73Q63JnnjuDdxDMw9L3Oc+11r+jlH1uRElzPP39XE3HMz17nX0WVNPnpz/H2d48yJuedm4/5W1uQk8/19Nfd8zHHm5BBNHnKQXLvc1wdYk+6+mORiklTVpe4+P+DnbxVzz4u556OqLq3rR+1y34GanOjyHGdOzH3c+1i3Oc+9jh+zy32afAhznHuOMyfmPu59rNsm/q2syeY+7n2s2xznnuPMyeGaPOStLa4kuX3H7VNJXj/AGgAOT5MBpkOTAaZFlwFGNOQg+YUkZ6vqTFXdnOSBJE9ft+bpJA8uP/30niTf6+7/99/1ADg0TQaYDk0GmBZdBhjRvm9t0d1Xq+qRJM8muSnJE939clU9vLz+WJJnktyf5HKSHyT5yICfffHAu95s5p4Xc8/HWmYescmJx21OzD0v5h6JJo9ijnPPcebE3HOz6X8re9zmxdzzMceZk0PMXd27vkUbAAAAAAAkGfbWFgAAAAAAzJiDZAAAAAAAVhr9ILmq7q2qV6rqclU9usv1qqpPLa+/WFV3jb2ndRgw94eW875YVc9V1Z3Hsc+jtN/MO9b9WlX9sKrev879jWXI3FX13qr6clW9XFX/vO49jmHA7/jPVtXfVtVXlnMPfU/IyaqqJ6rq21X10h7XJ98zTZ5PkxNdnlOX59jkZPO7rMmavMc6Td4Cc+zypjc50eU5dVmTNfm665o8tGfdPdpXFm9u/59JfiHJzUm+kuTcdWvuT/KFJJXkniT/Nuae1vE1cO5fT/K25ff3bfrcQ2bese4fs/iAg/cf977X9FjfkuQ/kpxe3v654973mub+4yR/tvz+tiT/neTm4977Ief+zSR3JXlpj+uT7pkmz6fJQ+fesU6XJ7D3kWfeuiYvZ9nYLmuyJmvyj9dsVZNvYO6t6/ImN/kGHrdJzzDi3FvVZU3WZE0+eM/GfkXy3Ukud/er3f1mkieTXLhuzYUkn+uF55PcUlXvGHlfY9t37u5+rru/u7z5fJJTa97jURvyWCfJHyb56yTfXufmRjRk7g8meaq7v5kk3b0Nsw+Zu5P8TFVVkp/OIsRX17vNo9XdX8xijr1MvWeaPJ8mJ7o8py7PssnJxndZkzVZkxe2rcnJTLu84U1OdHlOXdZkTdbkA/Zs7IPkk0le23H7yvK+G12zaW50po9m8SzAJtt35qo6meT3kjy2xn2Nbchj/UtJ3lZV/1RVX6qqB9e2u/EMmfvPk/xKkteTfDXJH3X3j9azvWMz9Z5p8sIcmpzo8lvm0GVN3tuUm6bJC5q8pMlb0+REl/cy9abp8sIcuqzJC5p8jSYP7NmJ0bazULvc1wdYs2kGz1RV78sixO8ZdUfjGzLzJ5N8rLt/uHiSZysMmftEkl9N8jtJfjLJv1bV89399bE3N6Ihc/9uki8n+e0kv5jkH6rqX7r7f0be23Gaes80+Zptb3Kiyztte5c1eW9TbpomX6PJC5+MJm9DkxNd3svUm6bL12x7lzX5Gk1e0OSFfXs29kHylSS377h9KovT/Rtds2kGzVRVdyR5PMl93f2dNe1tLENmPp/kyWWEb01yf1Vd7e7Pr2WH4xj6O/5Gd38/yfer6otJ7kyyySEeMvdHkvxpd3eSy1X1jSS/nOTf17PFYzH1nmnywhyanOjyW+bQZU3e25SbpskLmnyNJm9HkxNd3svUm6bLC3PosiYvaPI1mjywZ2O/tcULSc5W1ZmqujnJA0mevm7N00keXH5a4D1Jvtfd3xp5X2Pbd+6qOp3kqSQf3vBndt6y78zdfaa739nd70zyV0n+YMMjnAz7Hf+bJL9RVSeq6qeSvDvJ19a8z6M2ZO5vZvEsZqrq55O8K8mra93l+k29Z5o8nyYnujynLmvy3qbcNE3WZE1e2LYmJ7q8l6k3TZfn02VN1mRNPmDPRn1FcndfrapHkjybxackPtHdL1fVw8vrj2Xx6Zf3J7mc5AdZPAuw0QbO/fEkb0/ymeUzXFe7+/xx7fmwBs68dYbM3d1fq6q/T/Jikh8leby7Xzq+XR/ewMf7E0k+W1VfzeK/THysu984tk0fgar6yyTvTXJrVV1J8idJfiLZjJ5p8nyanOhyZtTluTY52ewua7Ima/J2NjmZb5c3ucmJLmdGXdZkTdbkg/esFq/aBgAAAACA3Y391hYAAAAAAGw4B8kAAAAAAKzkIBkAAAAAgJUcJAMAAAAAsJKDZAAAAAAAVnKQDAAAAADASg6SAQAAAABY6f8AolirgEtg5dQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 16 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(4, 4, figsize=(20,8), tight_layout=True)\n",
    "for i, (pfm, ax) in enumerate(zip(pfms, axes.ravel())):\n",
    "    ppm = pfm/pfm.sum(axis=1, keepdims=True)\n",
    "    ppm = pd.DataFrame(ppm, columns=[\"A\", \"C\", \"G\", \"T\"])\n",
    "    ic = ppm.applymap(get_information_content)\n",
    "    logomaker.Logo(ic, ax=ax)\n",
    "    ax.set_title(Ns[i])"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}