train_rand.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential\n",
    "from keras.layers import Activation, Dense, Dropout\n",
    "from keras.layers.core import Flatten, Dense\n",
    "from keras.layers import Conv2D, MaxPooling2D\n",
    "from keras.metrics import categorical_crossentropy\n",
    "from keras.optimizers import Adam\n",
    "from keras.preprocessing.image import ImageDataGenerator\n",
    "from keras.layers.normalization import BatchNormalization\n",
    "from keras.layers.convolutional import *\n",
    "from keras.models import load_model\n",
    "\n",
    "import keras\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import shutil\n",
    "import os\n",
    "import random\n",
    "\n",
    "from PIL import Image\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn import datasets\n",
    "\n",
    "import time\n",
    "from PIL import Image\n",
    "from numba import jit\n",
    "import pickle\n",
    "from pathlib import Path\n",
    "import skimage.measure\n",
    "%matplotlib inline\n",
    "\n",
    "from keras.optimizers import Adam\n",
    "\n",
    "from keras.models import load_model\n",
    "from sklearn import metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the nth ('n') of an integer 'num'\n",
    "def check_nth_bit(num, n):\n",
    "    return (num>>n)&1\n",
    "\n",
    "# flip the bit of the integer 'int' at a position 'p'  (staring at 0)\n",
    "def bitflip(int,p):\n",
    "    mask = 1 << p\n",
    "    return (int^ mask)\n",
    "\n",
    "# count the numbers of 1 on a binarie representetion of an interger n\n",
    "def CountBits(n):\n",
    "    n = (n & 0x5555555555555555) + ((n & 0xAAAAAAAAAAAAAAAA) >> 1)\n",
    "    n = (n & 0x3333333333333333) + ((n & 0xCCCCCCCCCCCCCCCC) >> 2)\n",
    "    n = (n & 0x0F0F0F0F0F0F0F0F) + ((n & 0xF0F0F0F0F0F0F0F0) >> 4)\n",
    "    n = (n & 0x00FF00FF00FF00FF) + ((n & 0xFF00FF00FF00FF00) >> 8)\n",
    "    n = (n & 0x0000FFFF0000FFFF) + ((n & 0xFFFF0000FFFF0000) >> 16)\n",
    "    n = (n & 0x00000000FFFFFFFF) + ((n & 0xFFFFFFFF00000000) >> 32) \n",
    "    return n\n",
    "\n",
    "# defining n!\n",
    "@jit\n",
    "def factorial(n):\n",
    "    if n == 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return n * factorial(n-1)\n",
    "\n",
    "# finding the states that have the maginetization \"z\"\n",
    "@jit\n",
    "def statep(y,z):\n",
    "    a = 0\n",
    "    for x in range(0,2**l):\n",
    "        if (2* CountBits(x)-l) == z:\n",
    "            y[a] = x\n",
    "            a += 1\n",
    "    return (y)\n",
    "\n",
    "# numbers of states with magnatization x\n",
    "def num_mag(x):\n",
    "    y = p = int(factorial(l)//((factorial((x+l)/2))*(factorial(l -((x+l)/2)))))\n",
    "    return y\n",
    "            \n",
    "# find the position in the list of states \"s\" that are iqual to \"b\"\n",
    "@jit\n",
    "def find_state(s,b):\n",
    "    for z in range(0,p+1):\n",
    "        if s[z] == b:\n",
    "            return(z)\n",
    "\n",
    "# power method\n",
    "@jit\n",
    "def power_iteration(A):\n",
    "    # Choose a random vector\n",
    "    y_0 = np.random.rand(A.shape[0])\n",
    "    #erro precision\n",
    "    e= 0.000000001\n",
    "    #convergence constant\n",
    "    er= 1.0\n",
    "    it = 0\n",
    "    while (er >= e):\n",
    "        # calculate the matrix-by-vector product Ab\n",
    "        z_1 = np.dot(A, y_0)\n",
    "        # calculate the norm\n",
    "        z_1_norm = np.linalg.norm(z_1)\n",
    "        # renormalize the vector\n",
    "        y_1 = z_1 / z_1_norm\n",
    "        z_2 = np.dot(A, y_1)\n",
    "        z_2_norm = np.linalg.norm(z_2)\n",
    "        lamb_1 = np.dot(z_2,y_1)/np.dot(y_1,y_1)\n",
    "        y_2 = z_2/z_2_norm\n",
    "        z_3 = np.dot(A, y_2)\n",
    "        z_3_norm = np.linalg.norm(z_3)\n",
    "        lamb_2 = np.dot(z_3,y_2)/np.dot(y_2,y_2)\n",
    "        er = abs(lamb_1 - lamb_2)/abs(lamb_2)\n",
    "        y_0 = z_3/z_3_norm\n",
    "        it += 1\n",
    "        #if( it%10 == 0):print('# of interations:',it,';','Precision:', er)\n",
    "    return (y_0, lamb_2)\n",
    "\n",
    "@jit\n",
    "def POWER_METHOD(A,dim):\n",
    "    for x in range(0,dim):\n",
    "        A[x,x] -= 2*l\n",
    "        \n",
    "    C = np.asarray(power_iteration(A))\n",
    "    C[1] += 2*l\n",
    "    return (C)\n",
    "\n",
    "def chart_regression(pred,y,sort=True):\n",
    "    t = pd.DataFrame({'pred' : pred, 'y' : y.flatten()})\n",
    "    if sort:\n",
    "        t.sort_values(by=['y'],inplace=True)\n",
    "    a = plt.plot(t['y'].tolist(),label='expected')\n",
    "    b = plt.plot(t['pred'].tolist(),label='prediction')\n",
    "    plt.ylabel('output')\n",
    "    plt.legend()\n",
    "    plt.title('RMSE ='+ str(np.round(100000*score)/100000))\n",
    "    plt.savefig('predct_pool.png')\n",
    "    plt.show()\n",
    "    \n",
    "#plt.switch_backend('agg')\n",
    "\n",
    "def plot_history(history):\n",
    "    loss_list = [s for s in history.history.keys() if 'loss' in s and 'val' not in s]\n",
    "    val_loss_list = [s for s in history.history.keys() if 'loss' in s and 'val' in s]\n",
    "    acc_list = [s for s in history.history.keys() if 'mean_squared_error' in s and 'val' not in s]\n",
    "    val_acc_list = [s for s in history.history.keys() if 'val_mean_squared_error' in s and 'val' in s]\n",
    "\n",
    "    if len(loss_list) == 0:\n",
    "        print('Loss is missing in history')\n",
    "        return\n",
    "\n",
    "    ## As loss always exists\n",
    "    epochs = range(1,len(history.history[loss_list[0]]) + 1)\n",
    "\n",
    "    ## Loss\n",
    "    plt.figure(1)\n",
    "    for l in loss_list:\n",
    "        plt.plot(epochs, history.history[l], 'b', label='Training loss (' + str(str(format(history.history[l][-1],'.5f'))+')'))\n",
    "    for l in val_loss_list:\n",
    "        plt.plot(epochs, history.history[l], 'g', label='Validation loss (' + str(str(format(history.history[l][-1],'.5f'))+')'))\n",
    "\n",
    "    plt.title('Loss')\n",
    "    plt.xlabel('Epochs')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.savefig('loss_retrain_pool.png')\n",
    "    plt.show()\n",
    "    \n",
    "\n",
    "    ## Accuracy\n",
    "    #plt.figure(2)\n",
    "    #for l in acc_list:\n",
    "        #plt.plot(epochs, history.history[l], 'b', label='Training Mean Square Error (' + str(format(history.history[l][-1],'.5f'))+')')\n",
    "    #for l in val_acc_list:\n",
    "        #plt.plot(epochs, history.history[l], 'g', label='Validation Mean Square Error (' + str(format(history.history[l][-1],'.5f'))+')')\n",
    "\n",
    "    #plt.title('Mean Square Error')\n",
    "    #plt.xlabel('Epochs')\n",
    "    #plt.ylabel('Mean Square Error')\n",
    "    #plt.legend()\n",
    "    \n",
    "    \n",
    "# ------------------------------------------------------    \n",
    "# ------------------PLOT CLECIO-------------------------\n",
    "# ------------------------------------------------------\n",
    "\n",
    "def plot_regression_results(\n",
    "        results,results_residuals,\n",
    "        truth_values,individual_truth_values,\n",
    "        individual_medians,plot_std=True,\n",
    "        bin_edges=10,real_data=None,\n",
    "        plot_dropout_err=True, plot_name='regression.png',\n",
    "        xlabel='$\\\\theta_{E}$',ylabel='$\\\\theta_{E}$',\n",
    "        xmin=None, xmax=None,\n",
    "        ymin=None, ymax=None,\n",
    "        color_std='blue', color_err='green',\n",
    "        color_pop='magenta',plot_label=None,\n",
    "        fill_label_std=None, fill_label_err=None,\n",
    "        figsize=(6, 8),alpha_hist=0.4, \n",
    "        pop_error=None,pop_error_residuals=None,\n",
    "        plot_real_residuals=False):\n",
    "    \"\"\" \n",
    "    Funtion to create Plots of true values vs predicted values.\n",
    "    It also plots residuals and the histogram of true and \n",
    "    predicted.\n",
    "    INPUT:\n",
    "        results: <list> - a list with the format [bin centers, \n",
    "        bin medians, bin standard deviations,bin 1 sigma \n",
    "        intervals, i.e., [low,high] limit ,bin 2 sigma intervals \n",
    "        ,bin 3 sigma intervals].\n",
    "        results_residuals: <list> - a list with the format [bin \n",
    "        centers, bin median residuals, bin standard deviation\n",
    "        residuals,bin 1 sigma interval of residuals,bin 2 sigma \n",
    "        intervals of residuals, bin 3 sigma intervals of \n",
    "        residuals].\n",
    "        truth_values: <array> - same as np.unique\n",
    "        (individual_truth_values).\n",
    "        individual_truth_value: <array> - The non bined truth \n",
    "        values for histogram.\n",
    "        individual_medians: <array> - The non bined predicted \n",
    "        medians for histogram.\n",
    "        plot_std <bool> - If True will plot a standard deviation\n",
    "        in results and results_residuals\n",
    "        bin_edges: <array> - Bin edges, used to define a scale \n",
    "        only, optional.\n",
    "        real_data: data points, instead of binned curves to add\n",
    "        in the plot. if plot_real_residuals=True it must follow \n",
    "        the format [true_valules,median,lowlim,highlim,\n",
    "        median_residual,lowlim_residual,highlim_residual] if \n",
    "        plot_real_residuals=False the format is [true_valules,\n",
    "        median,lowlim,highlim]\n",
    "        plot_dropout_err: <bool> -  if True will plot the \n",
    "        1-sigma, 2-sigma and 3-sigma error in results_catalog \n",
    "        and results_residuals.\n",
    " \n",
    "        plot_name: <str> : output plot name.\n",
    "    OUTPUT:\n",
    "       Saves a plot with a a path defined in the input plot_name.\n",
    "    \"\"\"\n",
    "\n",
    "    pop_err=pop_error\n",
    "    pop_err_r=pop_error_residuals\n",
    "    x_bin,y_bin,y_bin_std,y_bin_percentile1sig,y_bin_percentile2sig,y_bin_percentile3sig=results\n",
    "    x_bin,yr_bin,yr_bin_std,yr_bin_percentile1sig,yr_bin_percentile2sig,yr_bin_percentile3sig=results_residuals\n",
    "\n",
    "\n",
    "    if type(xmin)==type(None):\n",
    "        if type(bin_edges)!=type(10):\n",
    "            xmin=min(bin_edges)\n",
    "            xmax=max(bin_edges)\n",
    "        else:\n",
    "            xmin=min(x_bin)\n",
    "            xmax=max(x_bin)\n",
    "\n",
    "\n",
    "\n",
    "    imed=individual_medians\n",
    "    itruth=individual_truth_values \n",
    "\n",
    "    ty_bin_percentile1sig=np.transpose(y_bin_percentile1sig)\n",
    "    ty_bin_percentile2sig=np.transpose(y_bin_percentile2sig)\n",
    "    ty_bin_percentile3sig=np.transpose(y_bin_percentile3sig)\n",
    "\n",
    "    tyr_bin_percentile1sig=np.transpose(yr_bin_percentile1sig)\n",
    "    tyr_bin_percentile2sig=np.transpose(yr_bin_percentile2sig)\n",
    "    tyr_bin_percentile3sig=np.transpose(yr_bin_percentile3sig)\n",
    "    \n",
    "    # ================================\n",
    "    # Initialize figure canvas\n",
    "    # ================================\n",
    "    ymax = xmax\n",
    "    ymin = xmin\n",
    "    # set up figure\n",
    "    fig = plt.figure(figsize=figsize)         \n",
    "    \n",
    "    # set up grid of panels in figure\n",
    "    grid = plt.GridSpec(3, 2, hspace=0.0, wspace=0.0, \n",
    "                        left=0.1, right=1.0, bottom=0.1, top=1.0,\n",
    "                        height_ratios=[1, 4, 1], width_ratios=[4,1])\n",
    "         \n",
    "    # main axis: predicted statistic vs. true statistic\n",
    "    ax_main = fig.add_subplot(grid[1,0])\n",
    "    \n",
    "    # hist axis [x]: histogram of data on x axis on top\n",
    "    ax_hist_x = fig.add_subplot(grid[0, 0], sharex=ax_main)\n",
    "    \n",
    "    # hist axis [y]: histogram of data on y axis on right\n",
    "    ax_hist_y = fig.add_subplot(grid[1, 1], sharey=ax_main)\n",
    "    \n",
    "    # res axis: residual of statistic vs. true stratistic\n",
    "    ax_diff = fig.add_subplot(grid[2,0], sharex=ax_main)\n",
    "    \n",
    "\n",
    " \n",
    "    \n",
    "    # calculate upper and lower bounds [for systematic error bars]\n",
    "    if plot_std==True:\n",
    "        lo_err, up_err = measure_uplow_std_limits(y_bin, y_bin_std, nsigma=1)#(ys_bin, ys_err_bin)\n",
    "\n",
    "    # calculate upper and lower bounds on diff [for systematic error bars]\n",
    "    if plot_std==True:\n",
    "        lo_err_diff, up_err_diff = measure_uplow_std_limits(yr_bin, yr_bin_std, nsigma=1)\n",
    "        \n",
    "\n",
    "    # ================================\n",
    "    # Plot data [sim]\n",
    "    # ================================\n",
    "    # plot: dropout medians [main axis]\n",
    "    ax_main.plot(x_bin,y_bin, color=color_err, alpha=0.9, label=plot_label)\n",
    "    \n",
    "    # plot: dropouts percentile errors [main axis]\n",
    "    if plot_dropout_err==True:\n",
    "        ax_main.fill_between(x_bin, ty_bin_percentile1sig[0], ty_bin_percentile1sig[1], color=color_err, alpha=0.5, label=fill_label_err)\n",
    "        ax_main.fill_between(x_bin, ty_bin_percentile2sig[0], ty_bin_percentile2sig[1], color=color_err, alpha=0.4)\n",
    "        ax_main.fill_between(x_bin, ty_bin_percentile3sig[0], ty_bin_percentile3sig[1], color=color_err, alpha=0.3)\n",
    "    \n",
    "    # plot: one-to-one line [main axis]\n",
    "    ax_main.plot([xmin, xmax], [ymin, ymax], color='black', ls='dashed', alpha=0.9)\n",
    "    \n",
    "    # plot: errors from standard deviation [main axis]\n",
    "    if plot_std==True:\n",
    "        ax_main.fill_between(x_bin, lo_err, up_err, color=color_std, alpha=0.4, label=fill_label_std)\n",
    "    \n",
    "    # plot: deltas [diff axis]\n",
    "    ax_diff.plot(x_bin,yr_bin, color=color_err, alpha=0.9)\n",
    "    \n",
    "    # plot: statistical errors [diff axis]\n",
    "    if plot_dropout_err==True:\n",
    "        ax_diff.fill_between(x_bin, tyr_bin_percentile1sig[0], tyr_bin_percentile1sig[1], color=color_err, alpha=0.5, label=fill_label_err)\n",
    "        ax_diff.fill_between(x_bin, tyr_bin_percentile2sig[0], tyr_bin_percentile2sig[1], color=color_err, alpha=0.4)\n",
    "        ax_diff.fill_between(x_bin, tyr_bin_percentile3sig[0], tyr_bin_percentile3sig[1], color=color_err, alpha=0.3)\n",
    "    \n",
    "    # plot: systematic errors [main axis]\n",
    "    if plot_std==True:\n",
    "        ax_diff.fill_between(x_bin, lo_err_diff, up_err_diff, color=color_err, alpha=0.9, label=fill_label_std)\n",
    "    \n",
    "    # plot: errors from medians\n",
    "    if pop_error!=None:\n",
    "        xp_bin,yp_bin,yp_bin_std,yp_bin_percentile1sig,yp_bin_percentile2sig,yp_bin_percentile3sig=pop_err\n",
    "        xp_bin,ypr_bin,ypr_bin_std,ypr_bin_percentile1sig,ypr_bin_percentile2sig,ypr_bin_percentile3sig=results_residuals=pop_err_r\n",
    "        typ_bin_percentile1sig=np.transpose(yp_bin_percentile1sig)\n",
    "        typr_bin_percentile1sig=np.transpose(ypr_bin_percentile1sig)\n",
    "        ax_main.fill_between(xp_bin, typ_bin_percentile1sig[0], typ_bin_percentile1sig[1], color=color_pop, alpha=0.4)\n",
    "        ax_diff.fill_between(xp_bin, typr_bin_percentile1sig[0], typr_bin_percentile1sig[1], color=color_pop, alpha=0.4)\n",
    "    # plot: zero line [res axis]\n",
    "    xmin_line, xmax_line = ax_diff.get_xlim()\n",
    "    ax_diff.axhline(y=0.2, ls='dashed', color='magenta', lw=0.5)\n",
    "    ax_diff.axhline(y=0.1, ls='dashed', color='red', lw=0.5)\n",
    "    ax_diff.axhline(y=0, ls='dashed', color='black', lw=0.5)\n",
    "    ax_diff.axhline(y=-0.1, ls='dashed', color='red', lw=0.5)\n",
    "    ax_diff.axhline(y=-0.2, ls='dashed', color='magenta', lw=0.5)\n",
    "\n",
    "    #ax_diff.plot([xmin, xmax], np.zeros(2), color='black', ls='dashed', alpha=0.4) \n",
    "\n",
    "    \n",
    "    # ================================\n",
    "    # Plot data [simulation]\n",
    "    # ================================\n",
    "    # plot: histogram on upper x axis\n",
    "    ax_hist_x.hist(itruth, bin_edges, \n",
    "                   color=color_err, alpha=alpha_hist, \n",
    "                   histtype='stepfilled')\n",
    "    \n",
    "    # plot: histogram on right y axis\n",
    "    blank = ax_hist_y.hist(imed, bin_edges, \n",
    "                   color=color_err, alpha=alpha_hist,\n",
    "                   histtype='stepfilled',\n",
    "                   orientation='horizontal')\n",
    "        \n",
    "    # plot: sky data main_axis\n",
    "    if type(real_data)!=type(None):\n",
    "        print('real data')\n",
    "        print(np.array(real_data).shape)\n",
    "        if plot_real_residuals==False:\n",
    "            real_tr,real_me,real_lo,real_hi=np.transpose(real_data)\n",
    "        else:\n",
    "            real_tr,real_me,real_lo,real_hi,realr_me,realr_lo,realr_hi=np.transpose(real_data)\n",
    "        yerr_lo=real_me - real_lo\n",
    "        yerr_hi= real_hi - real_me\n",
    "\n",
    "        yerr_lor=realr_me - realr_lo\n",
    "        yerr_hir= realr_hi - realr_me\n",
    "\n",
    "        #res_me=np.divide(real_me-real_tr,real_tr)\n",
    "        ax_main.errorbar(real_tr, real_me, yerr=[yerr_lo,yerr_hi], fmt='ko', markersize=5,label='DES')\n",
    "        ax_diff.errorbar(real_tr, realr_me, yerr=[yerr_lor,yerr_hir], fmt='ko', markersize=5,label='DES')\n",
    "\n",
    "\n",
    "    #if xr is not None and yr is not None:\n",
    "    #    if yr_err is not None:\n",
    "    #        ax_main.errorbar(xr, yr, yerr=yr_err, color='black', fmt='o', label='DES')\n",
    "    #    if yr_diff is not None:\n",
    "    #        ax_diff.scatter(xr, yr_diff, color='black')\n",
    "    #    if yr_diff is not None and yr_err is not None:\n",
    "    #        ax_diff.errorbar(xr, yr_diff, yerr=yr_err, color='black', fmt='o')\n",
    "\n",
    "    # ================================\n",
    "    # Set plot options\n",
    "    # ================================\n",
    "    # set plot options [main axis]\n",
    "    plt.setp(ax_main.get_xticklabels(), visible=False)\n",
    "    plt.setp(ax_main.get_yticklabels(), visible=True)\n",
    "    plt.setp(ax_diff.get_yticklabels(), visible=False)\n",
    "    plt.setp(ax_hist_x.get_xticklabels(), visible=False)\n",
    "    plt.setp(ax_hist_x.get_yticklabels(), visible=False)\n",
    "    plt.setp(ax_hist_y.get_yticklabels(), visible=False)\n",
    "    plt.setp(ax_hist_y.get_xticklabels(), visible=False)\n",
    "    ax_main.set_ylabel(ylabel)\n",
    "    #ax_main.autoscale(enable=True,axis='both')\n",
    "    ax_main.axis('equal')\n",
    "    ax_main.set_xlim((xmin, xmax))   \n",
    "    ax_main.set_ylim((xmin, xmax))\n",
    "    ax_main.set_ylabel(ylabel)   \n",
    "    ax_main.grid(True, alpha=0.2) \n",
    "    ax_hist_x.grid(True,alpha=0.2)\n",
    "    ax_hist_y.grid(True,alpha=0.2)\n",
    "    ax_main.legend()\n",
    "    #ax_main.set_aspect('equal')\n",
    "    \n",
    "    # set plot options [residual axis]\n",
    "    ax_diff.set_ylim((-0.8,0.8))\n",
    "    ax_diff.set_xlabel(xlabel)\n",
    "    ylabel_plain=ylabel.replace('$','')\n",
    "    ax_diff.set_ylabel('$\\\\frac{\\\\delta '+ylabel_plain+'}{'+ylabel_plain+'}$')\n",
    "    ax_diff.grid(True, alpha=0.2)\n",
    "\n",
    "\n",
    "    plt.savefig(plot_name)\n",
    "    \n",
    "    \n",
    "def make_bin(col1_dat,col2_dat,bin_edges,return_vals=False,method=\"mean\"):\n",
    "    \"\"\"\n",
    "    col1= array like of the x values\n",
    "    col2= array like of the y values\n",
    "    Funtion that bins in x of a sample of (x,y) pairs.\n",
    "    USe the return_vals if you want a list of y values in each bin\n",
    "    \"\"\"\n",
    "    print (type(col2_dat))\n",
    "    print (np.array(col2_dat).shape)\n",
    "    #col1_dat=fits_catalog_col(cat,col1)\n",
    "    #col2_dat=fits_catalog_col(cat,col2)\n",
    "    #xy=zip(col1_dat,col2_dat)\n",
    "    bin_vals = [[] for x in range(0,len(bin_edges)-1)]\n",
    "    bin_valsx = [[] for x in range(0,len(bin_edges)-1)]\n",
    "    bin_median=numpy.zeros(len(bin_edges)-1)\n",
    "    bin_mean=numpy.zeros(len(bin_edges)-1)\n",
    "    bin_std=numpy.zeros(len(bin_edges)-1)\n",
    "    bin_center=numpy.zeros(len(bin_edges)-1)\n",
    "\n",
    "    bin_percentile1sig=[[] for x in range(0,len(bin_edges)-1)]\n",
    "    bin_percentile2sig=[[] for x in range(0,len(bin_edges)-1)]\n",
    "    bin_percentile3sig=[[] for x in range(0,len(bin_edges)-1)]\n",
    "\n",
    "    for i in range(0,len(col1_dat)):\n",
    "        for j in range(0,len(bin_edges)-1):\n",
    "            if bin_edges[j] <= col1_dat[i] <= bin_edges[j+1]:\n",
    "                #print \"data size\"\n",
    "                #try:\n",
    "                #    print col1_dat.size\n",
    "                #    print col2_dat.size\n",
    "                #except:\n",
    "                #    print len(col1_dat)\n",
    "                #    print len(col2_dat)\n",
    "                bin_vals[j].append(col2_dat[i])\n",
    "                bin_valsx[j].append(col1_dat[i])            \n",
    "\n",
    "    for i in range(0,len(bin_edges)-1):\n",
    "        try:\n",
    "            bin_mean[i]=numpy.mean(bin_vals[i])\n",
    "        except:\n",
    "            print(type(bin_vals[i]))\n",
    "            print(\"Error!\")\n",
    "        bin_median[i]=numpy.median(bin_vals[i])\n",
    "        bin_std[i]=numpy.std(bin_vals[i])\n",
    "        if method==\"percentile\":\n",
    "            if len(bin_vals[i])!=0:\n",
    "                #median=np.percentile(measurements,50.0)\n",
    "                lowlim=np.percentile(bin_vals[i],15.87)\n",
    "                highlim=np.percentile(bin_vals[i],84.13)\n",
    "                lowlim2sig=np.percentile(bin_vals[i],2.25)\n",
    "                highlim2sig=np.percentile(bin_vals[i],97.6)\n",
    "                lowlim3sig=np.percentile(bin_vals[i],0.15)\n",
    "                highlim3sig=np.percentile(bin_vals[i],99.7)\n",
    "            else:\n",
    "                lowlim=0\n",
    "                highlim=0\n",
    "                lowlim2sig=0\n",
    "                highlim2sig=0\n",
    "                lowlim3sig=0\n",
    "                highlim3sig=0\n",
    "            bin_percentile1sig[i]=[lowlim,highlim]\n",
    "            bin_percentile2sig[i]=[lowlim2sig,highlim2sig]\n",
    "            bin_percentile3sig[i]=[lowlim3sig,highlim3sig]\n",
    "            \n",
    "        if bin_std[i]<0.001:\n",
    "            print(\"zero bin std data:\")\n",
    "            print(bin_vals[i])\n",
    "        bin_center[i]= (bin_edges[i]+bin_edges[i+1])/2\n",
    "\n",
    "    if return_vals==False:\n",
    "        if method=='median':\n",
    "            return bin_center,bin_median,bin_std\n",
    "        elif method=='percentile':\n",
    "            return bin_center,bin_median,bin_std,bin_percentile1sig,bin_percentile2sig,bin_percentile3sig,       \n",
    "        else:\n",
    "            return bin_center,bin_mean,bin_std\n",
    "    else:\n",
    "        # if you want the a list of values inside each bin,\n",
    "        # so you can perform your own bin metrics.\n",
    "        return bin_center,bin_valsx,bin_vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_folder = 'matrix_del_rand/'\n",
    "ex = '.npy'\n",
    "\n",
    "if not os.path.exists(path_folder):\n",
    "    os.mkdir(path_folder)\n",
    "\n",
    "c_i = 8.\n",
    "\n",
    "l = 8\n",
    "j = 1\n",
    "delt = -1\n",
    "\n",
    "row = []\n",
    "for i in range(0,1000):\n",
    "    \n",
    "    h = c_i*np.random.randn(l)\n",
    "    \n",
    "    Hxxz = np.zeros((2**l , 2**l ))\n",
    "\n",
    "    # constructing the full xxz hamiltonian (2^L, 2^L) \n",
    "    for x in range(0,2**l):\n",
    "        for y in range(0,l):\n",
    "            k = (y+1)%l\n",
    "            if check_nth_bit(x,y) == 0:\n",
    "                Hxxz[x,x] += -j*h[y]/2\n",
    "            else:\n",
    "                Hxxz[x,x] += j*h[y]/2\n",
    "            #if y != l-1:    \n",
    "            if check_nth_bit(x,y) == check_nth_bit(x,k):\n",
    "                Hxxz[x,x] += -(j*delt)/2\n",
    "            else:\n",
    "                Hxxz[x,x] += (j*delt)/2\n",
    "                b = bitflip(x,y)\n",
    "                c = bitflip(b,k)\n",
    "                Hxxz[x,c] += -j\n",
    "    #print(Hxxz)\n",
    "       \n",
    "    \n",
    "    mag_array = np.zeros(l+1)\n",
    "    E_array = np.zeros(l+1)\n",
    "    a = 0\n",
    "    for mag in range(-l,l+1,2):\n",
    "    \n",
    "        p = num_mag(mag)\n",
    "        state = np.array(range(p))\n",
    "        s = statep(state,mag)\n",
    "        \n",
    "        Hxxz_p = np.zeros((p,p))\n",
    "        for x in range(0,p):\n",
    "            for y in range(0,l):\n",
    "                if check_nth_bit(s[x],y) == 0:\n",
    "                    Hxxz_p[x,x] += -(j*h[y])/2\n",
    "                else:\n",
    "                    Hxxz_p[x,x] += (j*h[y])/2\n",
    "                k = (y+1)%l\n",
    "                #if y != l-1:\n",
    "                if check_nth_bit(s[x],y) == check_nth_bit(s[x],k):\n",
    "                    Hxxz_p[x,x] += -(j*delt)/2\n",
    "                else:\n",
    "                    Hxxz_p[x,x] += (j*delt)/2\n",
    "                    b = bitflip(s[x],y)\n",
    "                    c = bitflip(b,k)\n",
    "                    D = np.where(s == c)\n",
    "                    Hxxz_p[x,D] = -j\n",
    "                                            \n",
    "        #print(Hxxz_p)\n",
    "        mag_array[a] = mag\n",
    "        A = POWER_METHOD(Hxxz_p,p)\n",
    "        E_array[a] = A[1] \n",
    "        psi_s = -mag_array[np.where(E_array == np.min(E_array))]\n",
    "        a += 1\n",
    "   \n",
    "    row.append([str(i+1) + ex, np.min(E_array)/l])\n",
    "    np.save(os.path.abspath(path_folder+str(i+1)), Hxxz)   \n",
    "    #print('h:',h)\n",
    "    #print('MAG_GS:', psi_s[0])\n",
    "    #print('rhoE_array:', np.min(E_array)/l)\n",
    "    \n",
    "    \n",
    "row = np.asarray(row)\n",
    "\n",
    "data_set = pd.DataFrame(row)\n",
    "data_set.to_csv(\"rand.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = pd.read_csv('rand.csv')\n",
    "\n",
    "path_folder = 'matrix_del_rand/'\n",
    "\n",
    "#writing the columns of the csv in an list\n",
    "X_1 = dataset.iloc[:,1].values\n",
    "Y_1 = dataset.iloc[:,2].values\n",
    "\n",
    "W = np.random.permutation(np.c_[X_1.reshape(len(X_1), -1), Y_1.reshape(len(Y_1), -1)])\n",
    "\n",
    "X = W[:, :X_1.size//len(X_1)].reshape(X_1.shape)\n",
    "Y = W[:, X_1.size//len(X_1):].reshape(Y_1.shape)\n",
    "\n",
    "x_train = np.array([X[i] for i in range((len(X)//10),len(X))])\n",
    "y_train = np.array([Y[i] for i in range((len(Y)//10),len(Y))])\n",
    "\n",
    "x_test = np.array([X[i] for i in range((len(X)//10))])\n",
    "y_test = np.array([Y[i] for i in range((len(Y)//10))])\n",
    "\n",
    "train_data = []\n",
    "\n",
    "load_train_data = []\n",
    "for file in x_train:\n",
    "        path = os.path.abspath(path_folder+ file)\n",
    "        H = np.load(path)\n",
    "        load_train_data.append(H)\n",
    "load_train_data = np.array(load_train_data)\n",
    "train_data = load_train_data[:,:,:,np.newaxis]\n",
    "train_data -= np.mean(train_data)\n",
    "train_data /= np.std(train_data)\n",
    "\n",
    "train_label = []\n",
    "for file in y_train:\n",
    "    train_label.append(file)\n",
    "train_label = np.array(train_label)\n",
    "\n",
    "test_data = []\n",
    "\n",
    "load_test_data = []\n",
    "for file in x_test:\n",
    "        path = os.path.abspath(path_folder+ file)\n",
    "        H = np.load(path)\n",
    "        load_test_data.append(H)\n",
    "load_test_data = np.array(load_test_data)\n",
    "test_data = load_test_data[:,:,:,np.newaxis]\n",
    "test_data -= np.mean(test_data)\n",
    "test_data /= np.std(test_data)\n",
    "\n",
    "test_label = []\n",
    "for file in y_test:\n",
    "    test_label.append(file)\n",
    "test_label = np.array(test_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Sequential()\n",
    "model.add(Conv2D(64,(3,3), activation='relu', input_shape=(256,256,1)))\n",
    "model.add(Conv2D(64,(3,3),activation='relu'))\n",
    "model.add(Conv2D(64,(3,3),activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "model.add(Conv2D(128,(3,3),activation='relu'))\n",
    "model.add(Conv2D(128,(3,3),activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "model.add(Conv2D(256,(3,3),activation='relu'))\n",
    "model.add(Conv2D(256,(3,3),activation='relu'))\n",
    "model.add(Conv2D(256,(3,3),activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(Conv2D(512,(3,3),activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "model.add(Flatten())\n",
    "model.add(Dense(4096, activation='relu'))\n",
    "model.add(Dropout(0.5))\n",
    "model.add(Dense(4096, activation='relu'))\n",
    "model.add(Dropout(0.5))\n",
    "\n",
    "model.add(Dense(1, activation='linear'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 810 samples, validate on 90 samples\n",
      "Epoch 1/15\n",
      "810/810 [==============================] - 369s 456ms/step - loss: 6.8917 - val_loss: 3.6280\n",
      "Epoch 2/15\n",
      "810/810 [==============================] - 287s 354ms/step - loss: 3.8130 - val_loss: 1.0445\n",
      "Epoch 3/15\n",
      "810/810 [==============================] - 287s 354ms/step - loss: 1.3473 - val_loss: 1.2583\n",
      "Epoch 4/15\n",
      "810/810 [==============================] - 287s 354ms/step - loss: 0.7850 - val_loss: 0.2348\n",
      "Epoch 5/15\n",
      "810/810 [==============================] - 287s 355ms/step - loss: 0.5184 - val_loss: 0.2316\n",
      "Epoch 6/15\n",
      "810/810 [==============================] - 288s 355ms/step - loss: 0.3455 - val_loss: 0.3061\n",
      "Epoch 7/15\n",
      "810/810 [==============================] - 289s 356ms/step - loss: 0.2986 - val_loss: 0.1903\n",
      "Epoch 8/15\n",
      "810/810 [==============================] - 284s 351ms/step - loss: 0.2562 - val_loss: 0.1819\n",
      "Epoch 9/15\n",
      "810/810 [==============================] - 287s 355ms/step - loss: 0.2178 - val_loss: 0.1705\n",
      "Epoch 10/15\n",
      "810/810 [==============================] - 287s 354ms/step - loss: 0.2163 - val_loss: 0.1762\n",
      "Epoch 11/15\n",
      "810/810 [==============================] - 286s 353ms/step - loss: 0.2081 - val_loss: 0.1747\n",
      "Epoch 12/15\n",
      "810/810 [==============================] - 288s 355ms/step - loss: 0.1988 - val_loss: 0.1486\n",
      "Epoch 13/15\n",
      "810/810 [==============================] - 288s 356ms/step - loss: 0.1929 - val_loss: 0.1508\n",
      "Epoch 14/15\n",
      "810/810 [==============================] - 288s 356ms/step - loss: 0.1788 - val_loss: 0.1384\n",
      "Epoch 15/15\n",
      "810/810 [==============================] - 286s 353ms/step - loss: 0.1644 - val_loss: 0.1323\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEWCAYAAABliCz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4FFW6x/Hvm+4AgQTCGpagYQlL2AKJCAKBoKiJBhD1jiguqIM67o6OuD7g1RnckdE7rrgyZhwRVxYdDNugIEGICMgaMCwSkS2QQJZz/6hOm5CFDulKd6ffz/PU093V1ad+CeGt6lNVp8QYg1JKqfovxNcBlFJK1Q0t+EopFSS04CulVJDQgq+UUkFCC75SSgUJLfhKKRUktOArpVSQ0IKvgpKIZIvIeb7OoVRd0oKvlFJBQgu+UmWIyB9FZIuI/CYin4pIe9d8EZHnRWSfiBwSkSwR6e16L1VE1ovIERHZJSL3+vanUKpyWvCVchGRkcDfgP8B2gE7gHTX2+cDSUA3IBL4A7Df9d4bwE3GmAigN/B1HcZWymNOXwdQyo9cBcw0xqwGEJEHgAMiEgMUAhFAD2ClMWZDmc8VAnEistYYcwA4UKeplfKQ7uEr9bv2WHv1ABhj8rD24jsYY74GXgReAn4RkVdFpKlr0UuBVGCHiCwWkcF1nFspj2jBV+p3u4EzS1+ISBOgJbALwBgzwxiTAPTC6tq5zzX/O2PMGKAN8DHwQR3nVsojWvBVMAsVkUalE1ahnigi8SLSEPgrsMIYky0iZ4nI2SISChwFCoBiEWkgIleJSDNjTCFwGCj22U+kVDW04KtgNhfILzMNAx4BZgN7gC7AFa5lmwKvYfXP78Dq6nnG9d7VQLaIHAZuBibUUX6lakT0BihKKRUcdA9fKaWChBZ8pZQKElrwlVIqSGjBV0qpIOFXV9q2atXKxMTE+DpGOUePHqVJkya+juERzWqfQMobSFkhsPL6Y9bMzMxfjTGtPVnWrwp+TEwMq1at8nWMchYtWsSIESN8HcMjmtU+gZQ3kLJCYOX1x6wisuPUS1m0S0cppYKEFnyllAoSthV8EekuImvKTIdF5C671qeUUqp6tvXhG2N+AuIBRMSBNQDVHLvWp1RdKCwsJCcnh4KCAl9HqVKzZs3YsGHDqRf0E4GU15dZGzVqRHR0NKGhoafdRl0dtD0X2GqM8fjgglL+KCcnh4iICGJiYhARX8ep1JEjR4iIiPB1DI8FUl5fZTXGsH//fnJycujUqdNpt1MnY+mIyExgtTHmxUremwRMAoiKikpIT08/eRGfysvLIzw83NcxPKJZ7VOat1mzZnTp0sVviz1AcXExDofD1zE8Fkh5fZnVGMPWrVs5dOhQufnJycmZxphET9qwveCLSAOsccZ7GWN+qW7ZxMREo6dlnj7Nap/SvBs2bKBnz56+jlOtQNpjhsDK6+uslf39iYjHBb8uztJJwdq7r7bYn67jx+Hpp+Grr+xoXSml6o+6KPjjgfftarxBA3jqKZg1y641KOUf9u/fT3x8PPHx8bRt25YOHTq4X584ccKjNiZOnMhPP/1U7TIvvfQSs7z0H2ro0KGsWbPGK22dbNWqVdx0002A1d3xpz/9ia5du9KvX78q1zl58mSio6OJjIys8N77779PXFwcvXr14pprrin33qFDh2jXrh33339/hc+lpqYSHx9fYf60adMQEQ4ePAjAxx9/zP/+7//W+Of0JlsP2opIY2AUcJN964Bhw2DJErvWoJR/aNmypbuQTZkyhfDwcO69995yyxhjKCkpqbKNN99885TrufXWW2sXtI488cQTPP744wB89tln/Pzzz2zZsoVly5Zx66238t///rfCZ8aMGcNtt91G7969y83fuHEjzzzzDMuXLycyMpJ9+/aVe//BBx8kOTm5QnsffPABkZGR7N69u9z87OxsFi9eTIcOHcqte8qUKdx33300atTotH/u2rB1D98Yc8wY09IYc+jUS5++pCTYvh1+/tnOtSjln7Zs2ULv3r25+eabGTBgAHv37mXSpEkkJibSq1cvHnvsMfeypXvcRUVFREZGMnnyZPr168fgwYPdRe7hhx9m+vTp7uUnT57MwIED6d69O8uXLwesMWUuvfRS+vXrx/jx40lMTDzlnvx7771Hnz596N27Nw8++CAARUVF/PGPf3TPnzFjBgDPP/88cXFx9OvXjwkTKt5A7NChQ2zcuJFevXoB8Mknn7j3yocOHcrevXvJzc2t8LnBgwfTtm3bCvNfffVVbr/9dveef5s2bdzvrVy5koMHDzJy5Mhynzl8+DAzZszggQceqNDe3XffzdNPP11unogwbNgw5s6dW/UvyWZ+NZbO6UpKsh6XLoUrr/RtFhU87roLvN1bER8PrlpbI+vXr+fNN9/k5Zdf5siRI0ybNo0WLVpQVFREcnIyl112GXFxceU+c+jQIYYPH860adO45557mDlzJpMnT67QtjGGlStX8umnn/LYY48xf/58/v73v9O2bVtmz57N2rVrGTBgQLX5cnJyePjhh1m1ahXNmjXjvPPO4/PPP6d169bs37+fH374AcDd/fHUU0+xY8cOGjRo4J5X1sqVK+nbt6/79a5du+jYsaP7dXR0NLt27aJ1a4/GFGPTpk2EhoYyZMgQSkpKmDp1Kueffz7FxcXcd999vP/++xUK9UMPPcT9999PWFhYufmzZ8+mc+fOFb5FACQmJrJ06VLGjRvnUS5vqxdDK/TrBxER2q2jgleXLl0466yz3K/ff/99BgwYwIABA9iwYQPr16+v8JmwsDBSUlIASEhIIDs7u9K2S4tT2WWWLVvGFVdYt/vt16+fe0+7KitWrGDkyJG0atWK0NBQrrzySpYsWULXrl3ZsmULd955JwsWLKBZs2YA9OrViwkTJjBr1qxKLzTas2dPuWJe2dmGNTl1tqioiG3btrF48WLee+89rr/+eg4fPszf//53xowZQ/v27cstn5mZSU5ODmlpaeXm5+Xl8dRTTzFlypRK19OmTZsK3T91qV7s4TscMGSItYevVF05nT1xu5QdsnfLli288MILrFy5ksjISCZMmFDplcENGjRwP3c4HBQVFVXadsOGDSssU9PTuatavmXLlixfvpxly5YxY8YMZs+ezauvvsqCBQtYvHgxn3zyCY8//jjr1q0rd/57WFhYuZ8pOjqan3/+mUGDBgHWN4qTi3R1oqOjGTFiBE6nky5dutClSxe2bt3Kt99+y/Lly5kxYwZ5eXmcOHGC5s2b065dO1asWEFMTAxFRUXs27ePc889l2effZbt27fTp08fAPbu3Uvfvn3JzMykdevWFBQUVPhGUJfqxR4+WN0669dDJd12SgWV0nPFmzZtyp49e1iwYIHX1zF06FA++OADAH744YdKv0GUNWjQIDIyMti/fz9FRUWkp6czfPhwcnNzMcZw+eWXM3XqVFavXk1xcTE5OTmMHDmSp59+mtzcXI4dO1auvZ49e7Jlyxb369GjR/POO+8A1rePqKgoj7tzAMaOHUtGRgYA+/btY+vWrXTq1In09HR27txJdnY206ZNY8KECTzxxBPcdttt7N69m+zsbBYtWkRcXBwLFy4kPj6effv2kZ2dTXZ2Nm3btiUrK8udZdOmTZV29dSVerGHD7/34y9bBpdc4tssSvlSfHw8cXFx9O7dm86dOzNkyBCvr+P222/nmmuuoW/fvgwYMIDevXu7u2MqEx0dzWOPPcaIESMwxpCWlsZFF13E6tWrmThxIiKCiPDkk09SVFTElVdeyZEjRygpKeH++++vcLFTr169yM3Ndd+QJC0tjXnz5tGlSxeaNGnC22+/DVhXxp599tnu+2zcc889fPDBBxw+fJjo6GhuvvlmHn74YS666CK++uor4uLicDqdPP/885WeullbGRkZPPfcc15v12PGGL+ZEhISzOkqKDCmUSNj7rrrtJuoVEZGhncbtJFmtU9p3vXr1/s2iAcOHz5s+zoKCwtNfn6+McaYTZs2mZiYGFNYWHhabZ1u3qeeesq8+eabp/XZ01Wb3+2uXbvMqFGjarX+yv7+gFXGwxpbb/bwGzaEQYP0wK1SdSEvL49zzz2XoqIijDG88sorOJ11W05uu+02PvroozpdZ238/PPPPPPMMz7NUG8KPljdOo8/DocPQ9Omvk6jVP0VGRlJZmamTzOEhYVx1VVX+TRDTZx99tm+jlB/DtqCdcVtSQm4rg1RSilVRr0q+IMHg9Op3TpKKVWZelXwmzSBhAQt+EopVZl6VfDB6sdfuRLy832dRCml/Eu9LPiFhVbRV0op9bt6V/CHDLGGTNZuHVUfjRgxosKVs9OnT+dPf/pTtZ8rvZ3k7t27ueyyy6ps+1R3nJs+fXq5q15TU1MrHdyspqZMmWLbKYv5+fkMHz6c4uJiAN5++21iY2OJjY11X6B1sn//+9/06tWLkJCQcr+TVatWue9B0K9fP+bMmQNYp1wmJyfTs2dPevXqxQsvvOD+zCOPPELfvn2Jj4/n/PPPd4+lM2vWLPr27Uvfvn0555xzWLt2LQAnTpwgKSmpyqEuaqPeFfzmzaFPHy34qn4aP348J9/3OT09nfHjx3v0+fbt2/Phhx+e9vpPLvhz58615YpUb5o5cybjxo3D4XDw22+/MXXqVFasWMHKlSuZOnUqBw4cqPCZ3r1789FHH5FUegm/S1xcHKtWrWLNmjXMnz+fm266iaKiIpxOJ88++ywbNmzg22+/5aWXXnIPN3HfffeRlZXFmjVruPjii93DVXfq1InFixeTlZXFI488wqRJkwBrjKNzzz2Xf/3rX17/XdSr8/BLJSXBzJlW104lA+0p5RV3zb+LNXu9Oz5yfNt4pl9Y9ahsl112GQ8//DDHjx+nYcOGZGdns3v3boYOHUpeXh5jxoxh//79FBcX8/jjjzNmzJhyn8/Ozubiiy9m3bp15OfnM3HiRNavX0/Pnj3JL3Pg65ZbbuG7774jPz+fyy67jKlTpzJjxgx2795NcnIyrVq1IiMjg5iYGFatWkWrVq147rnnmDlzJgA33ngjd911F9nZ2aSkpDB06FCWL19Ohw4d+OSTT6odQGzNmjXcfPPNHDt2jC5dujBz5kyaN2/OjBkzePnll3E6ncTFxZGens7ixYu58847AWt0zCVLllQYhmHWrFn885//BGDBggWMGjWKFi1aADBq1Cjmz59fYYNZ1X2LGzdu7L7ArKCgwD0iZ7t27WjXrh0AERER9OzZk127dhEXF0fTMhcFHT161P2Zc845xz1/0KBB5OTkuF+PHTuWBx54wOvXGdS7PXywCv6xY7B6ta+TKOVdLVu2ZODAgcyfPx+w9u7/8Ic/ICI0atSIOXPmsHTpUjIyMvjzn/9c7aiW//jHP2jcuDFZWVk89NBD5S6keuKJJ1i1ahVZWVnuvdA77riD9u3bk5GR4R5orFRmZiZvvvkmK1as4Ntvv+W1117j+++/B2Dz5s3ceuut/Pjjj0RGRjJ79uxqf8ZrrrmGJ598kqysLPr06cPUqVMB65aB33//PVlZWbz88ssAPPPMM7z00kusWbOGpUuXVtiQnDhxgm3bthETEwNUPW5+TaxYsYJevXrRp08f9waorOzsbL7//vtyF1o99NBDdOzYkVmzZpW7IU2pN954wz1UNVjfML777rsa5fJEvdzDHzbMelyyBPzg4jZVT1W3J26n0m6dMWPGkJ6e7t6rNsbw4IMPsmjRIpxOJ7t27eKXX36p9A5PAEuWLOGOO+4AcPcll/rggw949dVXKSoqYs+ePaxfv77c+ydbtmwZl1xyiXuY5nHjxrF06VJGjx5Np06d3Pd8rW7cfbBuynLw4EGGDx8OwLXXXsvll1/uznjVVVcxduxYxo4dC8CQIUO45557uOqqqxg3bhzR0dHl2vv111/LdTlVtgGsybj5YF0x++OPP7JhwwauvfZaUlJS3LcszMvL49JLL2X69Onl9uyfeOIJnnjiCf72t7/x4osvujdiYA2o9sYbb7Bs2TL3PIfDQYMGDdwjn3pLvdzDb9sWunXT8fFV/TR27FgWLlzI6tWryc/Pd99tatasWeTm5rJkyRLWrFlDVFRUpePgl1VZsdu+fTvPPPMMCxcuJCsri4suuuiU7VT3TaJ0PH2oftz9U/niiy+49dZbyczMJCEhgaKiIiZPnszrr79Ofn4+gwYNYuPGjeU+U9W4+aVqOm5+WT179qRJkyasW7cOgMLCQi699FL3xqcyV155ZblvOFlZWdx444188skntGzZstyyx48f9/q9b+tlwQdrL3/pUmuoBaXqk/DwcEaMGMH1119fru/50KFDtGnThtDQUDIyMtixY0e17SQlJTFr1iwA1q1bR1ZWFmDdq7VJkyY0a9aMX375hXnz5rk/ExERwZEjRypt6+OPP+bYsWMcPXqUOXPmMKz0q3YNNGvWjObNm7PUtbf27rvvMnz4cEpKStxnwjz11FMcPHiQvLw8tm7dSp8+fbj//vtJTEysUPCbN29OcXGxu+hfcMEFfPnllxw4cIADBw7w5ZdfcsEFF3icLzs7273B2rFjBz/99BMxMTEYY7jhhhvo2bMn99xzT7nPbN682f38008/pUePHgDs3LmTcePG8e6779KtW7dyn9m/fz+tW7eu9G5ftWFrl46IRAKvA70BA1xvjPnGznWWSkqCN96Adeugmm+iSgWk8ePHM27cuHJn7Fx11VWkpaUxfPhwBgwY4C4sVbnllluYOHGi+5TBgQMHAtYtC/v370+vXr0qjKc/adIkUlJSaNeuXbl+/AEDBnDddde527jxxhvp379/td03VXn77bfdB207d+7Mm2++SXFxMRMmTODQoUMYY7j77ruJjIzkkUceISMjA4fDQVxcXLl+8FLnn38+y5Yt47zzzqNFixY88sgj7ttBPvroo+4DuDfeeCM333wziYmJzJkzh9tvv53c3Fwuuugi4uPjWbBgAd988w1XXHEFoaGhhISE8H//93+0atWKZcuW8e6779KnTx9399Vf//pXUlNTmTx5Mj/99BMhISGceeaZ7uMPjz32GPv373efUut0Ot2ngGZkZJCamlrj390peTqO8ulMwNvAja7nDYDI6pavzXj4J9u+3Rgw5u9/r107gTRuu2a1j46Hbx+7865evdpMmDDBK23V1e/2kksuMRs3bqwwv7bj4dvWpSMiTYEk4A3XhuWEMab2V2h46MwzoWNHPR9fqWDXv39/kpOT3Rde+bsTJ04wduxYunfv7vW2xdTwZsQeNywSD7wKrAf6AZnAncaYoyctNwmYBBAVFZVw8kUltfHEEz1ZvTqSDz/8hhoeiHfLy8tzX6Xo7zSrfUrzNmvWjC5dutT4zI66VFxcXO6G3/4ukPL6Mqsxhq1bt3Lo0KFy85OTkzONMYkeN2LHBCQCRcDZrtcvAP9b3We82aVjjDGvvGJ162zadPptBFLXg2a1T2nebdu2mdzcXFNSUuLbQNXQLh37+CprSUmJyc3NNdu2bavwHn5yi8McIMcYs8L1+kNgso3rq6Ds+fixsXW5ZlVfRUdHk5OTQ25urq+jVKmgoMDrp/PZKZDy+jJro0aNKlxnUFO2FXxjzF4R+VlEuhtjfgLOxereqTM9ekCrVlbBv+GGulyzqq9CQ0Pp1KmTr2NUa9GiRfTv39/XMTwWSHkDKWtl7L7S9nZglog0ALYBE21eXzki1umZeuBWKaVsvvDKGLPGGJNojOlrjBlrjKk4LJ3NkpIgOxvKXFynlFJBqd5eaVuqdHRTHWZBKRXs6n3B79sXmjbVbh2llKr3Bd/hsO6CpQVfKRXs6n3BB6tbZ8MG2LfP10mUUsp3gqbgA5QZbloppYJOUBT8xERo1EgP3CqlgltQFPwGDWDwYO3HV0oFt6Ao+GANs7BmDZw07pBSSgWNoCn4SUnW3a+WL/d1EqWU8o2gKfiDBoHTqd06SqngFTQFv0kT6+CtHrhVSgWroCn4YHXrrFwJ+fm+TqKUUnUvqAr+sGFQWAgrVpx6WaWUqm+CquAPGWINmaz9+EqpYBRUBb95c2swNS34SqlgFFQFH6x+/G++sbp2lFIqmARlwT92DFav9nUSpZSqW0FX8Mve2FwppYJJ0BX8qCjo1k0LvlIq+ARdwQerW2fpUigu9nUSpZSqO7YWfBHJFpEfRGSNiKyyc101kZRkDaK2bp2vkyilVN2piz38ZGNMvDEmsQ7W5RG9sblSKhgFZZfOmWfCGWdoP75SKrjYXfAN8KWIZIrIJJvXVSPDhlkF3xhfJ1FKqbohxsaKJyLtjTG7RaQN8BVwuzFmyUnLTAImAURFRSWkp6fblqeszz5rx3PPdeedd1bQsWPVo6nl5eURHh5eJ5lqS7PaJ5DyBlJWCKy8/pg1OTk50+Muc2NMnUzAFODe6pZJSEgwdWXDBmPAmNdeq365jIyMOsnjDZrVPoGUN5CyGhNYef0xK7DKeFiHbevSEZEmIhJR+hw4H/Cb82K6d4fWrfXArVIqeDhtbDsKmCMipev5pzFmvo3rqxER62wdPXCrlAoWthV8Y8w2oJ9d7XtDUhLMng07d1pn7SilVH0WlKdlliodV0e7dZRSwSCoC37fvtC0qXbrKKWCQ1AXfIcDhg7VPXylVHAI6oIPVj/+hg2wb5+vkyillL204LvG1Vm2zLc5lFLKbkFf8BMSICxM+/GVUvVf0Bf8Bg1g0CAt+Eqp+i/oCz5Y3Tpr1lhj5CulVH2lBR+r4BsDy5f7OolSStlHCz5Wl47Tqd06Sqn6TQs+0LgxnHWWFnylVP2mBd9l2DD47js4dszXSZRSyh5a8F2SkqCwEFas8HUSpZSyhxZ8lyFDrCGTdZgFpVR9pQXfJTIS+vXTfnylVP2lBb+MpCTr1MwTJ3ydRCmlvE8LfhnDhkF+Pqxe7eskSinlfVrwyyi9IYp26yil6iMt+GVERVk3N9cDt0qp+kgL/kmSkqyCX1zs6yRKKeVdWvBPkpRkDaK2bp2vkyillHfZXvBFxCEi34vI53avyxtKb4ii/fhKqfqmLvbw7wQ21MF6vOKMM6xJC75Sqr6xteCLSDRwEfC6nevxtqQkq+Ab4+skSinlPWJsrGoi8iHwNyACuNcYc3Ely0wCJgFERUUlpKen25bHU59/3o5nn+3OO++soHnzXMLDw30dySN5eXma1SaBlDeQskJg5fXHrMnJyZnGmESPFjbG2DIBFwP/53o+Avj8VJ9JSEgw/mDjRmPAmNdeMyYjI8PXcTymWe0TSHkDKasxgZXXH7MCq4yHddnOLp0hwGgRyQbSgZEi8p6N6/Oabt2gTRvtx1dK1S+2FXxjzAPGmGhjTAxwBfC1MWaCXevzJhHrqlst+Eqp+kTPw69CUhLs2AG//NLQ11GUUsor6qTgG2MWmUoO2PqzESOsx8zM5j7NoZRS3uJRwReRLiLS0PV8hIjcISKR9kbzrT59oGNH+Oablr6OopRSXuHpHv5soFhEugJvAJ2Af9qWyg+IQFoarFrVgoICX6dRSqna87TglxhjioBLgOnGmLuBdvbF8g9paVBQ4ODrr32dRCmlas/Tgl8oIuOBa4HSMXFC7YlUc0dPHOW3/N+83u6IEdCoUTGffeb1ppVSqs55WvAnAoOBJ4wx20WkE+AX59QfKzxG1DNRPPfNc15vu1EjOOus3/j8cx1mQSkV+Dwq+MaY9caYO4wx74tIcyDCGDPN5mweaRzamAHtBjB381xb2h88eD85ObBmjS3NK6VUnfH0LJ1FItJURFoAa4E3RcT7u9SnKTU2le/3fs+eI3u83vagQfsRQbt1lFIBz9MunWbGmMPAOOBNY0wCcJ59sWompWsKAPO3zPd6282bFzJoEHz6qdebVkqpOuVpwXeKSDvgf/j9oK3f6BvVl/YR7Zm7xZ5unbQ0yMyE3bttaV4ppeqEpwX/MWABsNUY852IdAY22xerZkSE1K6pfLn1SwqLC73eflqa9fi5323qlFLKc54etP23MaavMeYW1+ttxphL7Y1WMymxKRw+fphvcr7xetu9ekFMjPbjK6UCm6cHbaNFZI6I7BORX0RktutuVn7jvM7n4Qxx2nK2TulVt//5Dxw75vXmlVKqTnjapfMm8CnQHugAfOaa5zeaNmzKsDOG2XZ65ujRUFBgFX2llApEnhb81saYN40xRa7pLaC1jblOS0rXFH7Y9wM5h3O83nZSEjRtqt06SqnA5WnB/1VEJoiIwzVNAPbbGex0pMamAjBv8zyvt92gAVxwgXXgtqTE680rpZTtPC3412OdkrkX2ANchjXcgl+Jax3HGc3OsPX0zL17rVM0lVIq0Hh6ls5OY8xoY0xrY0wbY8xYrIuw/IqIkNI1hf9s+w8nik94vf3UVAgJ0W4dpVRgqs0dr+7xWgovSo1NJe9EHst2LvN62y1bwjnn6FW3SqnAVJuCL15L4UUjO42kgaOBrWfrrF0LO3fa0rxSStmmNgXfLwcMDm8QTtKZSczb4v0Dt6BX3SqlAle1BV9EjojI4UqmI1jn5Ff32UYislJE1orIjyIy1avJq5HaNZX1uevJPpjt9ba7d4euXbUfXykVeKot+MaYCGNM00qmCGOM8xRtHwdGGmP6AfHAhSIyyFvBq2Pn6ZmlV91+/TXk5Xm9eaWUsk1tunSqZSylJTHUNdVJN1C3lt3oFNnJ1m6dEyfgq69saV4ppWwhxsZ794mIA8gEugIvGWPur2SZScAkgKioqIT09HSvrPuFzS8wf+98PhnyCQ1CGpx2O3l5eYSHh5ebV1QkXHLJOQwd+iv33/9TbaN6TWVZ/VUgZYXAyhtIWSGw8vpj1uTk5ExjTKJHCxtjbJ+ASCAD6F3dcgkJCcZbvtj0hWEKZsGWBbVqJyMjo9L548cb07q1MUVFtWreq6rK6o8CKasxgZU3kLIaE1h5/TErsMp4WItt69I5aaNyEFgEXFgX6wMYETOCho6GtvTjg9Wtk5sLK1fa0rxSSnmdbQVfRFqLSKTreRjWLRE32rW+kzUObUxyp2Tbhlm48EJwOPRsHaVU4LBzD78dkCEiWcB3wFfGmDo9ez3jPQRuAAAZFUlEQVS1ayqb9m9i629bvd528+YwbJgWfKVU4LDzLJ0sY0x/Y90pq7cx5jG71lWVlFjr5uZ2nq2zbh1s325L80op5VV10ofvK11bdCW2RaxtwyyUXnWre/lKqUBQrws+WBdhZWRnkF+Y7/W2Y2OhRw8t+EqpwFDvC35K1xQKigpYlL3IlvbT0mDxYjh82JbmlVLKa+p9wR8eM5wwZ5it3TqFhbBggS3NK6WU19T7gt/I2YhzO5/L3C1zSy8C86rBg6FFCx0jXynl/+p9wQerW2fbgW1s/m2z19t2OuGii2DuXCgq8nrzSinlNUFT8AFbu3V++w2++caW5pVSyiuCouB3at6Jnq162lbwL7gAQkP1bB2llH8LioIP1l7+4h2LOXriqNfbbtoUhg/Xgq+U8m9BU/BTY1M5UXyCr7d/bUv7aWmwcSNs9v5hAqWU8oqgKfhDzxhKeINwvepWKRW0gqbgN3Q25NxO5zJvyzxbTs/s1Al699aCr5TyX0FT8MHq1tlxaAcbft1gS/tpabB0KRw4YEvzSilVK0FV8Ovi9MziYpg/35bmlVKqVoKq4Hds1pHebXrbNlzywIHQurV26yil/FNQFXywboqydMdSDh/3/mhnDsfvV90WFnq9eaWUqpXgK/ixqRSWFLJw20Jb2h89Gg4dgmXLbGleKaVOW9AV/HM6nkPThk1t69YZNQoaNNBuHaWU/wm6gh/qCGVU51HM3WzP6Jnh4TBypFXwbWheKaVOW9AVfLC6dXYd2cUP+36wpf20NNiyBX76yZbmlVLqtNhW8EWko4hkiMgGEflRRO60a101dWHXCwGYt9mebp2LL7YedYx8pZQ/sXMPvwj4szGmJzAIuFVE4mxcn8faR7Qnvm08c7fYcz7+GWdAv37aj6+U8i+2FXxjzB5jzGrX8yPABqCDXeurqdSuqfx35385WHDQlvZHj4bly2H/fluaV0qpGhM7DlxWWIlIDLAE6G2MOXzSe5OASQBRUVEJ6enptucB+OHQD9yx5g6mxE1heOvhVS6Xl5dHeHh4jdvfuDGCW25J4IEHNnD++b/UJqrHTjerLwRSVgisvIGUFQIrrz9mTU5OzjTGJHq0sDHG1gkIBzKBcadaNiEhwdSVwuJCEzkt0kz8eGK1y2VkZJxW+8XFxrRta8zll5/Wx0/L6Wb1hUDKakxg5Q2krMYEVl5/zAqsMh7WY1vP0hGRUGA2MMsY85Gd66opZ4iTC7pcwLwt8ygxJV5vPyTEOng7fz6cOOH15pVSqsbsPEtHgDeADcaY5+xaT22kdE1hb95e1u5da0v7aWlw5AgsXmxL80opVSN27uEPAa4GRorIGteUauP6aqz09Ey7Rs887zxo1EjP1lFK+Qc7z9JZZowRY0xfY0y8a7Knsp6mqPAoEtsn2nZ6ZuPGVtHXq26VUv4gKK+0LSulawrf5nzLb/m/2dJ+WhpkZ8OPP9rSvFJKeSzoC35qbColpoQvt35pS/ulV91qt45SyteCvuCf1f4sWoa1tK0fv317SEjQYRaUUr4X9AXfEeLggq4XMH/LfFtOzwSrW2fFCti3z5bmlVLKI0Ff8MEaZiH3WC6ZuzNtaX/0aOug7Rdf2NK8Ukp5RAs+cEHXCxDEtm6d+HiIjtZ+fKWUb2nBB1o1bsXADgNtuwuWiHXw9ssvoaDAllUopdQpacF3SY1NZeWuleQezbWl/bQ0OHoUMjJsaV4ppU5JC75LamwqBsOCrQtsaX/kSOtCLO3WUUr5ihZ8lwHtBtCmSRvbunUaNYLzz4fPP9erbpVSvqEF3yVEQriw64XM3zKf4pJiW9aRlgY//wxr7RmrTSmlqqUFv4zUrqn8lv8bK3ettKX9iy6yDuBqt45Syhe04JcxqssoQiTEtm6dqCgYONC66la7dZRSdU0LfhktwlowOHqwbefjA/zP/8CqVXDdddZZO0opVVe04J8kNTaVzD2Z7M3ba0v7d94JU6bAu+9ae/vr19uyGqWUqkAL/klSuqYAsGCLd0/PzC/MZ97medy14Ha29L2aLxYU8OuvcNZZ8N57Xl2VUkpVyunrAP4mvm087cLbMXfLXK6Nv7ZWbW39bStzN89l3pZ5ZGRnUFBUQJgzjPyifJoPbM73389g/Hi4+mpYsgReeAHCwrz0gyil1Em04J9EREjpmsJHGz+iqKSoRp8tKCpgyY4l7iK/af8mAGJbxDJpwCRSY1NJOjOJBxc+yPQV0xnVeRQLF6bx6KPwt7/BypXwwQfQrZsdP5lSKthpwa9ESmwKM9fM5Nucb0+5bPbBbOZtnsfcLXP5evvXHCs8RkNHQ5I7JXPrWbeS0jWF2Jax5T4z7bxpLNqxiImfTGTtzWv56187MGwYTJgAiYnw+uvWwV2llPImLfiVGNV5FA5xMHfzXM53nF/uvRPFJ1i6Yynztsxj7ua5bPh1AwCdIjsxMX4iqbGpjIgZQePQxlW239DZkPRL0xnw6gCunnM1X139FSkpDtasgT/8wZqWLIFnn4WGDW39UZVSQcS2gi8iM4GLgX3GmN52rccOzRo1Y+gZQ62C3+N8cg7nuPfi/7PtP+SdyKOBowHDzxzOHwf8kdTYVLq17IaIeLyO7q2682LKi1z/6fU8+d8neXDYg3TsCIsXwwMPWMX+22/h3/+GTp1s/GGVUkHDzj38t4AXgXdsXIdtUrqmMHnhZK7Pu57ti7cDcEazM5jQZwIpsSmM7DSS8AbhtVrHdfHX8eW2L3k041GSY5IZ3HEwoaHwzDMwbJh1rn7//vD22zBmjBd+KKVUULPttExjzBLgN7vat9tlcZcR1SSKyNBInh71NOtuWUf2ndn84+J/MLr76FoXe7AOEL980ct0bNaRKz+6kkMFh9zvjRkDq1dDbCyMHQt//jMUFtZ6lUqpIKbn4VehS4su7L13L8/1e457z7mXXm161ajLxlPNGjXj/Uvf5+dDP3PT5zdhyoy50KkTLFsGt98Ozz0HSUmwc6fXIyilgoQYGwd1EZEY4PPq+vBFZBIwCSAqKiohPT3dtjynIy8vj/Dw2u/Nn8qsnbN4ffvr/KX7X0hpm1Lh/UWLWvP0091xOg0PPLCBQYMqfnmqq6zeEEhZIbDyBlJWCKy8/pg1OTk50xiT6NHCxhjbJiAGWOfp8gkJCcbfZGRk1Ml6ioqLzMi3R5rGTzQ2G3M3VrrMpk3GxMcbA8ZMnmxMYWH59+sqqzcEUlZjAitvIGU1JrDy+mNWYJXxsMZql46fcIQ4ePeSdwlzhnHF7Cs4XnS8wjKxsbB8OUyaBNOmwbnnwu7dPgirlApIthV8EXkf+AboLiI5InKDXeuqL9pHtOetsW+xZu8aJv9ncqXLhIXBK69Y4+9kZkJ8PPznP3UcVCkVkOw8S2e8MaadMSbUGBNtjHnDrnXVJxd3u5g7Bt7B9BXT+WLTF1Uud9VV8N130KaNdevEKVOg2J4bdSml6gnt0vFDT456kn5R/bjuk+vYc2RPlcv17AkrVsA118DUqXDTTYm88AL88ksdhlVKBQwt+H6okbMR6Zelc6zwGFfPuZoSU1Llsk2awFtvWV08DofhrrugfXu48EJrXl5e3eVWSvk3Lfh+qkerHrxw4Qss3L6Qp//79CmXv+oqeOWVTNavt4Zm2LjRGnY5KsoalG3ePCiq2eCfSql6Rgu+H7uh/w1cHnc5D2c8zIqcFR59pmdPePxx2LYNli61iv7cuZCaCh06WHfc+u47vaeuUsFIC74fExFeTXuVDhEdGD97fLmhF04lJASGDoWXX4Y9e+Djj60rdV95xbq1Yvfu8NhjsHWrjT+AUsqvaMH3c5GNIvnnpf9k56Gd3PLFLeWGXvBUw4bW2Dz//jfs3WuNt9+hg3VmT9euMHgwvPQS5OZ6P79Syn9owQ8A53Q8hykjpvD+uvd5Z23tBh+NjIQbboCMDNixA5580jqwe9tt1sHetDT417/g2DEvhVdK+Q0t+AHigaEPMPzM4dw691b3rRNrq2NH+Mtf4IcfYO1auPtu+P57uOIK62DvdddZF3XpKJ1K1Q9a8AOEI8TBe+Peo6GzIVd8WPnQC7XRty889ZS11//119YtFufMgVGjrG8FI0fCI4/AggVw+LBXV62UqiNa8ANIdNNoZo6eyfd7v+fBhQ/asg6HA5KT4Y03rP7+2bPhxhvh0CH461+t8/ubN7eGdLjtNnj/ffj5Z1uiKKW8TO9pG2DG9BjDrWfdynPfPsd5nc8jJbbiUMreEhYG48ZZE8CRI9aVvf/9rzVO/9tvWwd7weoeGjoUhgyxHnv3tjYeSin/oQU/AD096mmW7FjCtR9fS9YtWbQNb1sn642IgPPOsyawLuTKyvp9A7B4sbXHD9C0qXX2T+lGYOBA66pgpZTvaMEPQGGhYaRflk7iq4lcM+ca5k+YT4jUfe+c0wkDBljT7bdbF3Pt2GEV/9KNwKOPWvOdTuv+vKUbgJKSBhQX67cApeqSFvwAFdc6jukXTuemz2/i2eXPct+Q+3wdCRGIibGmCROseQcOwDff/L4B+Mc/4PnnAc4BrG8NkZG/T82bl39d3fymTa0LzJRSntGCH8D+OOCPfLn1Sx78+kFGxIzwdRy340XH2Zu3l91HdrMnbw+7W+3GjNxNp7N20+DwHrbl7uaXw/toEdKVpsd7EHa0B46DPTie24MdP3ciK8vJgQPWgeLqiFhFv+zGoFmz36emTcs/VjYvIkI3Gip4aMEPYCLCa2mvsXLXSsbPHs8LcS/Yur4TxSd+L+RH9rD7yG5ryiv/en/+/gqfdYiDdhHtaB/Rnj4duhDbpAMFjQr4af9c9jacCS2AzhAaEkrXFl3p36oH3Vr0ICa8B22dPWhJd4qPNuPgQetbw8GDlU/btlmnjR4+bG0wSqoeaNQtIuLUG4hffunAtm3WshER1vzS56VTgwbe/50r5U1a8ANc87Dm/PPSfzL8reFMXDWRNpvaICKESAiCICKVPoZISJXvlf08wP78/ew5sofcYxXHXnCIg7bhbWkf0Z7OzTszpOMQ2ke0p31Ee3eBbx/RnlaNW5U7zrBo0SJGjBgBwMGCg/z0609s/HWjNe3fyIZfN/DZps8oKvl9iM924e3o0aoH3Vt3p0fPHpzTqgc9WvWgY7OOlR7DMAaOHv29+Jc+ln1e2XsHDkB29u/zrKuOY0/5b9GwYdUbhLKvyz53uv4Hinjv8YcfWlBQYB0fcTisbzClz2s6r7L5Tufv80vXqwKDFvx6YOgZQ3n3knd5ZckrtGrVyrphMabcY4kpqTCvuvdKTIn7dUxkDIOjB/9eyMPLF3JHSO2OvEY2iuTs6LM5O/rscvMLiwvZfnD77xsC1/SvH//FgYID7uXCnGF0b9WdHq16cEbTMwh1hBIiITjEgSPEUf6xqYOQZtZ7jUMcRIiDM8osU9nnTImDNas30rVzP04UODme76Tg2O+P+cec5Oc5yT/q5FieNR09Yk17DjjZvMNJ3mEnRw45OXrEAXhSJQ1IyWlOkcBOVzOl65Lyz09+z/26mveMAwrDoLiBe37pxqB0I1DZhqGq16Xzjh7tT2Sk1ZY3J6cTQkPLT5XNq8n7mzaF07Jl9cuUtuOPG0Mt+PXElX2upP3+9u695vog1BFKt5bd6NayG6O7j3bPN8aQeyy3wreClbtW8tGGjyguKabY2HC/xy2neD/MNbWufjGHOHCGOAnB2lCWmBJKKMGYEvdzfyYITsJwmjCcNMJpwnAQhrMkDIdxTSVhhJSEEVJc/lGKwpDiMExRGMXFYRQXhlFwOJ/DTcJdOxoGU4L16HptPbf+3UuKXPNLwGDcy5WdV/q8uMhBSaGT4kInRSdCKS50YoqdUBwKJc4y00mvy71/8rLxUGLwZKNdduNxqg1Jmzbw+ef2/9tpwVcBR0Ro06QNbZq0YdiZw6pczhhDsSl2bwCqeiwxJadc5rtV39E3vi9FJUUUlRRRbIrdz2szlXaveXP66aef6NGjh3tk1dIiWPq89Hdz8uvq3gMoKikivyif/ML88o8V5v1KQVEB+YX5HC3zXkFRgbf/FHxKCCGEEKTMROlz47BemxCKCKHIhHDc9bri5GB7SWtgie2ZteCrektEcIoTZ0jt/8yPbT5W7cbFnyw6vIgR/Uf4OkYFJaaE40XHK2wgVny3gsTExHLHkYAKx5Y8mQe455eYEopKiigsKXRvYAuLC8ttcD19r/T9zds2c+aZZ1rfxKqYSnciPJlKl23aoGmd/BvYWvBF5ELgBcABvG6MmWbn+pRS/itEQggLDSMsNMzq+nLZH76f+LbxvgtWA4uKFwV0t6ltZyCLiAN4CUgB4oDxIhJn1/qUUkpVz85LTgYCW4wx24wxJ4B0YIyN61NKKVUNOZ1b5nnUsMhlwIXGmBtdr68GzjbG3HbScpOASQBRUVEJ6enptuQ5XXl5eYSHh/s6hkc0q30CKW8gZYXAyuuPWZOTkzONMYmeLGtnH35l5y1V2LoYY14FXgVITEw0/tY/VvYCIX+nWe0TSHkDKSsEVt5AyloZO7t0coCOZV5HA7ttXJ9SSqlq2FnwvwNiRaSTiDQArgA+tXF9SimlqmFbl44xpkhEbgMWYJ2WOdMY86Nd61NKKVU9W8/DN8bMBebauQ6llFKese0sndMhIrnADl/nOEkr4Fdfh/CQZrVPIOUNpKwQWHn9MeuZxphTjOBk8auC749EZJWnpzz5mma1TyDlDaSsEFh5AylrZfReP0opFSS04CulVJDQgn9qr/o6QA1oVvsEUt5AygqBlTeQslagffhKKRUkdA9fKaWChBZ8pZQKElrwKyEiHUUkQ0Q2iMiPInKnrzOdiog4ROR7EamDO2PWjohEisiHIrLR9Tse7OtMVRGRu11/A+tE5H0RaeTrTGWJyEwR2Sci68rMayEiX4nIZtdjc19mLKuKvE+7/hayRGSOiET6MmOpyrKWee9eETEi0soX2U6XFvzKFQF/Nsb0BAYBtwbAzVvuBDb4OoSHXgDmG2N6AP3w09wi0gG4A0g0xvTGGiLkCt+mquAt4MKT5k0GFhpjYoGFrtf+4i0q5v0K6G2M6QtsAh6o61BVeIuKWRGRjsAoYGddB6otLfiVMMbsMcasdj0/glWQOvg2VdVEJBq4CHjd11lORUSaAknAGwDGmBPGmIO+TVUtJxAmIk6gMX424qsxZgnw20mzxwBvu56/DYyt01DVqCyvMeZLY0yR6+W3WCPr+lwVv1uA54G/UMlw7/5OC/4piEgM0B9Y4dsk1ZqO9QdY4usgHugM5AJvurqgXheRJr4OVRljzC7gGaw9uT3AIWPMl75N5ZEoY8wesHZegDY+zlMT1wPzfB2iKiIyGthljFnr6yynQwt+NUQkHJgN3GWMOezrPJURkYuBfcaYTF9n8ZATGAD8wxjTHziKf3U5uLn6vscAnYD2QBMRmeDbVPWXiDyE1Z06y9dZKiMijYGHgEd9neV0acGvgoiEYhX7WcaYj3ydpxpDgNEiko113+CRIvKebyNVKwfIMcaUfmP6EGsD4I/OA7YbY3KNMYXAR8A5Ps7kiV9EpB2A63Gfj/OckohcC1wMXGX89+KgLlgb/7Wu/2/RwGoRaevTVDWgBb8SIiJYfcwbjDHP+TpPdYwxDxhjoo0xMVgHFL82xvjtXqgxZi/ws4h0d806F1jvw0jV2QkMEpHGrr+Jc/HTA8wn+RS41vX8WuATH2Y5JRG5ELgfGG2MOebrPFUxxvxgjGljjIlx/X/LAQa4/qYDghb8yg0BrsbaW17jmlJ9HaoeuR2YJSJZQDzwVx/nqZTrW8iHwGrgB6z/L351ab2IvA98A3QXkRwRuQGYBowSkc1YZ5NM82XGsqrI+yIQAXzl+r/2sk9DulSRNaDp0ApKKRUkdA9fKaWChBZ8pZQKElrwlVIqSGjBV0qpIKEFXymlgoQWfFXviUhxmdNr14iI167sFZGYykZTVMofOX0dQKk6kG+Mifd1CKV8TffwVdASkWwReVJEVrqmrq75Z4rIQtf47AtF5AzX/CjXeO1rXVPpMAsOEXnNNW7+lyIS5lr+DhFZ72on3Uc/plJuWvBVMAg7qUvnD2XeO2yMGYh1ted017wXgXdc47PPAma45s8AFhtj+mGN//Oja34s8JIxphdwELjUNX8y0N/Vzs12/XBKeUqvtFX1nojkGWPCK5mfDYw0xmxzDZa31xjTUkR+BdoZYwpd8/cYY1qJSC4QbYw5XqaNGOAr181GEJH7gVBjzOMiMh/IAz4GPjbG5Nn8oypVLd3DV8HOVPG8qmUqc7zM82J+PzZ2EfASkABkum6iopTPaMFXwe4PZR6/cT1fzu+3MrwKWOZ6vhC4Bdz3EG5aVaMiEgJ0NMZkYN2cJhKo8C1DqbqkexwqGISJyJoyr+cbY0pPzWwoIiuwdn7Gu+bdAcwUkfuw7s410TX/TuBV16iJxVjFf08V63QA74lIM0CA5/38Vo4qCGgfvgparj78RGPMr77OolRd0C4dpZQKErqHr5RSQUL38JVSKkhowVdKqSChBV8ppYKEFnyllAoSWvCVUipI/D85YbUIe7p2vgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Score (RMSE): 0.5194892048729426\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4HNXZt++j3otV3GRZcrfcccMNY4pNDZjQqyH0kBASEnhDSAipX+gQmgm99+YAoYONce+9Spas3nvbPd8fZ0e7Wq2klbWSLOu5r0vXaGd2Z84KPL95utJaIwiCIAh+Pb0AQRAE4ehABEEQBEEARBAEQRAEByIIgiAIAiCCIAiCIDgQQRAEQRAAEQRBEATBgQiCcNSglEpXStUopSqVUrlKqReUUhEux19QSmml1E/cPvewY/8Sx+sgpdQDSqksx7kOKqUeauU61s+/fbD+2xzrLlNKPaeUCm7lfSmO9bpe/26X4xcqpVYqpaqVUt96+PzZSqltjs+tVEqltXKdrx3XCXDZN1sptUYpVaGU2qKUmtvZ7y0cO4ggCEcbZ2utI4DJwBTg/9yO7wGusl44bnYXAPtd3vN/wDRgBhAJLAA2erqOy88tnVm0UmoRcCdwMpACDAP+3M7HYlyu/xeX/cXAw8A/PVxnJPAqcCMQA3wMfOR603e87zLAfV8/4CPgPsdn/wV8rJSK9fJrCsc4IgjCUYnWOhf4H0YYXPkYmONyEzsN2ALkurxnOvC+1jpbG9K11i918ZKvAp7VWm/XWpcAfwGWHMmJtNZfaq3fArI9HF4ELNdar9BaNwL/DxgMzLfeoJSKBv4E/M7ts7OBPK3121prm9b6FaAAOO9I1ikce4ggCEclSqkk4HRgn9uhWsxT7sWO11cC7jf7VcCvlVI3K6UmKKVUJ9ZxqVKqtI2fZMdbxwGbXT66GeivlIpr4/QZDrfW80qpeG+X5Phxfz3eZd/fgSdpLpKePmvtG48gIIIgHH18oJSqADKBfMyTrjsvAVc6noTnAx+4Hf8H5sn5MmAdcFgpdZXbez5wu7Ff52kxWuvXtNYxbfwccrw1Aihz+aj1e6SH0xZirJihwFTHe171dH0PfAHMV0qdqJQKAn4PBAFhAEqpacAc4DEPn10JDFJKXaKUCnT8TYZbnxUEEQThaONcrXUkcCIwBmjx5Ky1XgEkAH8Almmta9yO27TWj2ut52B85X8DnlNKjXW7juuN/ZlOrrsSiHJ5bf1e4WH9lVrrdVrrRq11HnALsFApFeX+Xg+f3YVxT/0byMH8fXYAWUopP+AJ4FaHO8n9s0XAOcCvgTyMu+1LIMvrbykc04ggCEclWuvvgBeA+1t5yyvAb2jpLnI/T43W+nGgBPCYjdMWSqnL3LKB3H8sl9F2YJLLRydh/PVFXlzGajnslWtLa/2O1nq81joOY0ENBdZiRGga8KZSKtexD4xYzHN89jut9XStdT/gCmA0sMab6wrHPgHtv0UQeoyHgXSl1GSt9Sa3Y48Cy4Hv3T+klPoVsAlYDTRgXEeRtMw0ahet9at45855CXhBKfUq5sn9DxhBa4FSaiZQCuwFYjHf5VutdZnjuD8QiPn36aeUCgFsWusGx/Gpju/XD2MpfKy13uWIlQxyudQQzM1+KiZ4jFJqCrANCAXuBbK01v/z6o8hHPOIhSActWitCzA32rs9HCvWWn+lPQ/0qAEewARVC4GfAz/VWh9wec/Hbk/673dyrZ9h0ji/ATIcP03xD6XUdkcqKJiU1M8w7qRtQB1wicvprnB8hyeBeY7fXV1aj2AEZbdje51jDVprnWv94BABjKVS7/j9d5i/SSYwEFjcme8tHFsoGZAjCIIggFgIgiAIggMRBEEQBAEQQRAEQRAciCAIgiAIQC9LO42Pj9cpKSk9vQxBEIRexfr16wu11gntva9XCUJKSgrr1q3r6WUIgiD0KpRSGd68T1xGgiAIAiCCIAiCIDgQQRAEQRCAXhZD8ERDQwNZWVnU1tb29FKOGUJCQkhKSiIwMLCnlyIIQjfS6wUhKyuLyMhIUlJS6MQcFMGB1pqioiKysrJITU3t6eUIgtCN9HqXUW1tLXFxcSIGPkIpRVxcnFhcgtAH6fWCAIgY+Bj5ewpC3+SYEARBEIRjlbKaBv7+yU4OFFR2+bVEEHoh6enpvPbaax3+3JIlS3jnnXe6YEWCIHQVqw4UsfT7AxRW1rf/5k4igtALOVJBEAShl9BQA6uXwppnqFjzKmcGbWRyv4Yuv6wIgg945ZVXmDFjBpMnT+aGG24gIyODkSNHUlhYiN1uZ968eXz++eekp6czZswYrrrqKiZOnMj5559PdXU1AOvXr2f+/PlMnTqVRYsWkZOTA8C+ffs45ZRTmDRpEscddxz79+/nzjvvZPny5UyePJmHHnoIm83Gb3/7W6ZPn87EiRN5+umnAZMxdMstt5CWlsaZZ55Jfn5+j/2NBEHoAHs/h09/C5/czvkZ9/K4330EFWzt8sv2+rRTV/788XZ2ZJf79Jxpg6L409njWj2+c+dO3nzzTX744QcCAwO5+eab+e6777jjjju48cYbmTlzJmlpaSxcuJD09HR2797Ns88+y5w5c7jmmmt44oknuPXWW/nFL37Bhx9+SEJCAm+++SZ33XUXzz33HJdddhl33nknixcvpra2Frvdzj//+U/uv/9+li1bBsDSpUuJjo5m7dq11NXVMWfOHBYuXMjGjRvZvXs3W7duJS8vj7S0NK655hqf/n0EQegCig8CkLdkFRc+vZJfzR3A4qTpXX7ZY0oQeoKvvvqK9evXM326+Y9VU1NDYmIi99xzD2+//TZPPfUUmzY558MPGTKEOXPmAHD55Zfz6KOPctppp7Ft2zZOPfVUAGw2GwMHDqSiooLDhw+zeLEZexsSEuJxDZ9//jlbtmxpig+UlZWxd+9evv/+ey655BL8/f0ZNGgQJ510Upf9HQRB8CGlGRAay/eFEWToAYw5bh6ERHX5ZY8pQWjrSb6r0Fpz1VVX8Y9//KPZ/urqarKysgCorKwkMjISaJnSqZRCa824ceP48ccfmx0rL/fO2tFa89hjj7Fo0aJm+z/55BNJIRWE3khJBsQMZeX+IuIjghjdP7JbLisxhE5y8skn88477zT554uLi8nIyOCOO+7gsssu49577+W6665rev+hQ4eabvyvv/46c+fOZfTo0RQUFDTtb2hoYPv27URFRZGUlMQHH3wAQF1dHdXV1URGRlJRUdF0zkWLFvHkk0/S0GCCTnv27KGqqooTTjiBN954A5vNRk5ODt988023/E0EQfCS3Z/CS+eAzS1gXJqBjh3Kin2FzB4ej59f9zzYiSB0krS0NP7617+ycOFCJk6cyKmnnkp6ejpr165tEoWgoCCef/55AMaOHcuLL77IxIkTKS4u5qabbiIoKIh33nmHO+64g0mTJjF58mRWrlwJwMsvv8yjjz7KxIkTmT17Nrm5uUycOJGAgAAmTZrEQw89xLXXXktaWhrHHXcc48eP54YbbqCxsZHFixczcuRIJkyYwE033cT8+fN78k8lCIIrlQXwwc1w4Fso2ufcb7dD6SFKggZRUFHH3BHx3bYkpbXutot1lmnTpmn3ATk7d+5k7NixPbSijpGens5ZZ53Ftm3benop7dKb/q6C0Ct56yrYYax/LngBxplYIWWH4aE0fhhzF5dtGseKOxaQFBvWqUsppdZrrae19z6xEARBELqb7R8YMTjht4CCgj3OY6VmuNnqkghS4sI6LQYdQQShG0lJSekV1oEgCF1IVRH89zcwcDLMvxNihkDhbufxEiMIX+SEMqcb3UUggiAIgtA92G2w8RV4ai7UlsG5T4B/ACSMaW4hlKSjUeyvj+3W+AH0kCAope5TSu1SSm1RSr2vlIrpiXUIgiB0C7lb4al58OHPIWoQLPkv9HekycePgsI9RjAAW3E6hSqO+OhI5o9O6NZl9pSF8AUwXms9EdgD/F8PrUMQBKHr+fIeqMw1weNrv4Tkmc5jCaPBVtcUO8jN2M0BWzx/OXc8YUHdWyrWI4Kgtf5ca93oeLkKSOqJdQiCIHQ5dZVw8HuYeLHJJHIvFo0fbbYFe9hfUIlfWQY6Zignj+3f7Us9GmII1wCftnZQKXW9UmqdUmpdQUFBNy6r54iIiAAgOzub888/v833Pvzww00N8gDOOOMMSktLu3R9giC4UJ4N1cWtHz/wDdjqYfRpno8njAJAF+zmj++up78qYfy4iV2w0PbpMkFQSn2plNrm4eccl/fcBTQCr7Z2Hq31Uq31NK31tISE7vWn+RKbzdbhzwwaNKjd+QXugvDJJ58QEyMhGUHoNl45Hz64qfXjuz+D4GhInuX5eGgsRPRn/451ZGfsxQ9NRP/hXbPWdugyQdBan6K1Hu/h50MApdRVwFnAZbo3Vcd5oLW21ikpKdx7773MnTuXt99+m/3793PaaacxdepU5s2bx65duwA4ePAgs2bNYvr06dx9993Nzjt+/HjACMrtt9/OhAkTmDhxIo899hiPPvoo2dnZLFiwgAULFgAmtbWwsBCABx98kPHjxzN+/HgefvjhpnOOHTuW6667jnHjxrFw4UJqamq6888lCMcOVUWQv924hBo9DLCx22Hv/2DkKeAf2OppKiOHUZm1g7OG1JkdsUO7aMFt0yPN7ZRSpwF3APO11tXtvd9rPr3TRPN9yYAJcPo/232bp7bWYDqUrlixAjB9j5566ilGjhzJ6tWrufnmm/n666+59dZbuemmm7jyyit5/PHHPZ5/6dKlHDx4kI0bNxIQEEBxcTH9+vXjwQcf5JtvviE+vnl62vr163n++edZvXo1WmtmzpzJ/PnziY2NZe/evbz++us888wzXHjhhbz77rtcfvnlnfxDCUIfJHO12TZUQ/YGSD6++fHD66GqAEad3uopyqob+Co/moV+mxk5KdCk3MSmdNmS26KnYgj/BiKBL5RSm5RST/XQOnyGe1trSwQuuugiwHQ8XblyJRdccEHTIB1rCM4PP/zAJZdcAsAVV1zh8fxffvklN954IwEBRsP79evX5npWrFjB4sWLCQ8PJyIigvPOO4/ly5cDkJqayuTJkwGYOnUq6enpnfjmgtCHyVwNfo7n6oPLWx7f8ykof2MheEBrze3vbGZL/UAiqCY8bx34B0PEgC5cdOv0iIWgtR7RJSf24km+q/DU1hogPDwcALvdTkxMTLPZCG193h2tdYdaWbflhQsODm763d/fX1xGgnCkZK6GQVOgoRYOfgfzf9v8+O7PTOwgNLZpl92u2ZhZwje7CvhqVz47c8p5cvYs2PA87P3CVC779cyz+tGQZXRM4KmttStRUVGkpqby9ttvA+aGvXnzZgDmzJnDG2+8AcCrr3qOry9cuJCnnnqKxkaTrVtcbLIa3FthW5xwwgl88MEHVFdXU1VVxfvvv8+8efN88E0FQQCgsQ4Ob4AhMyF1HmSuMcJgUXrIxBdcsouKq+q58rk1/PTJH3ni231EBPtzz9lpnHbiCeYNNcU95i4CEQSf4amttTuvvvoqzz77LJMmTWLcuHF8+OGHADzyyCM8/vjjTJ8+nbKyMo/nv/baa0lOTmbixIlMmjSJ1157DYDrr7+e008/vSmobHHcccexZMkSZsyYwcyZM7n22muZMmWKj7+1IPRhcjabgrLk4yFlnvk9a63z+O7PzNYRP9h2uIyzH1vBmoPF/OnsNDbevZC3b5zNkjmpqMgBJhMJIKZnAsog7a99Qm9qa+0tR8PfVRC6nbpKeP40OP0+GNpKmqjFD4/CF3fD7XvBPwj+lWq6ly74Pdga4el5YG+EW9by6dYcfvXmJvqFB/HU5VOZNMRDavh/ToWsNXDqvTDnVp9+LWl/LQiC0FEKdptMRWtOQVtkrobYVIhIhNAYGDjJpJ8CrFkK+Tvg5D+yYm8hv3xjI+MGRfHxL+Z6FgNoKlDrSQtBBMEHSFtrQThGcPQTakonbQ2tzXtc00xT5kHWOijaD9/8HUacyrbIedzw8jqGJ0Tw/NUziI8Ibv2cVgsLiSF0jt7k9uoNyN9T6LNYgpCzBeqrWn9f8QFTXzDEpUld6nywN8BrF4KtnsOz7mHJC2uJCQvixWtmEB3aemEaYPoczbje2QW1B+j1ghASEkJRUZHcxHyE1pqioiJCQkJ6eimCcGTYbVBbfmSfdQynQdsge2Pr77MsCFcLIfl4U5NQtI+KaT/n4nfzabRrXrxmBv2jvPj3FDMEzrivzYrmrqZH6hB8SVJSEllZWfSVxnfdQUhICElJ0oC2T9JYZwKhQeE9vZKOoTXkbIItb8O2d83T/Q3fQVwHewKVHjIum5J0c9NPmev5fYdWQUi0080DEBwBQ47HVnqIC7YdT2lVA69cO5MRiRFH+q26nV4vCIGBgaSmpvb0MgTh2OCDm6HkIFz3dU+vpGOsfMxk/PgFwsiFkL7CfJerPwE/f+/PU5phgsP+waauwJ38nbDjQ9j5ESTNaFFAVnLmUm54eS2HKjQvXTOj9QDyUUqvFwRBEHxEbTns/Njk01cXQ1jb7VGOKrLWQkwyXP+dWffmN+D9G+DHx2HOL707h91uLITRZ0BwFOxaZvb5+Zk00pfPhfTlgIKhs+GkuwDjZt12uJx3N2Tx0eZsquqCeX7JNKal9KK/nwMRhN5IfRVseNkEoHqoxF04BtnzmREDMC6RMWf07Ho6QmmGcd9YIjbxIiNuX/8FRp4KiV7U1FTmmbkFMclmrOXGl6Fon0kH3fGBEYP5d8C0n/Hc5mp++LyQ/IoV5JTVUlhZR1CAH6em9efaualMSY5t/3pHIXI36Y3s/QI+uwNyt/T0SoRjie3vQ+RAU2R16MeeXk3HKElvnq6pFJz1EARHwvs3Ns0rbhMrwyg2xZk9lLnaxCeWPwgJY2D+nby5q457l+0gvaiKuIggTh6TyN8Wj2ft70/h8UuP67ViAGIh9E7qHL2Laj23uRCEDlNbDvu+hOnXweF1vUsQakrMvwX3/P2IRDjpblj2K+P7HzC+7fOUHjLbmKEQN8I0pMtcDeEJpifR4qfZml3B3R9uZ+6IeF68Zgb+ft43nOwNiCD0Rqz86LojTK0TBHd2f2rcJePONWmPPz4O9dUQFNbTK2ufknSz9TRUJtXRNO7w+vYFwUo5tbqNJs0wglCwG6KTKR12Njc9sZr48CAeuXjyMScGIC6j3kl9pdnWtexyKvRStrwNjx9vgpc9wfb3ISoJBk8z7ZrtDeYm2htoEoSUlsf6DTPpod58l9J0iOgPgaHm9ZAZULgHstawNukKrn55E/nldTxx+VTi2qo47sWIIPRGLAvhSItvhKOP9OVQsNPcgLqb2jLY/5WxDvz8IHkmoDruNtIaig92yRLbxBIETz2AlILBU800s/YoPUR95BAWP/EDM//+Jdd+Y9JVC3QUl68fyaGiav750wlM7mWppB1BXEa9EXEZHXtYAc3cLdA/rXuvbbmL0s41r0NjITGt44Kw4SX4+Jdw5UcwbL7v19kaJRkQFgchUZ6PDzoOVjzUrgvMXpLBj9Up7K2r5MwJAwlXc6naGUPWqBv5dMGppMaHd2hIVW9EBKE3YrmMJKh87GD5r3M2w6SLu/faOz4y7qIkl+7IycfDljeNC8vfi9tEXSV8/Vfz+w8Pty4IjXXw+sUw/06HJeID3DOM3Bk81bSiyN3ScuaxA21rQJdmsq1xMg9fOplT0vqbA7b9TPHm+x8jiMuoNyIxhGMLuw3KsszvOT2QSnx4vQm+uj79Dp1t/j/L2+rdOVY+BlX5kHYO7P+69e+Rv9Mc3/Vx59dtUZLedsvowceZ7eHW3Uavf7kKf+yMTZvgFAPwTgyPIUQQeiPiMup6bI3w1b1QVdT11yrPNkHcwHDzFGu3d/01LSrzoTK3ZQZOsmM4zKFV7Z+jIhdWPmpcTmc/CkGR5rUnCnabbf7OI1+zK7ZGKMts20KIHGAsINfAstaQt4Oaeht//2QnH31nvueCme3OkDmmEUHojUhQuevJ2QzLH4CdH3b9taz4waiFRuRL07v+mha5DgtgwITm+6MHQ3QyZKxs/xzf/A1sDXDKn8ygmGlLYNt7TjeYKwW7zDZvR6eW3UT5YdOMr70ZAoOnNBME+9Z34MlZ/P3+f7L0+wOcP8xkdylPqat9CBGE3kiTy0gEocsoP2y2Rfu7/lpWlszYs822O91GliD095Cjnzyz/XTN/F2w8RWYfq1J8QSYeZNxP616ouX7LUGoyDYFZZ3FQ8qp1pqdOeV8uzuft9dl8tR3+/m2MhlKDnL7S99w+kPfcuDdPwFwVePbvHHdTM4frkH5GUuiD9O3HGTHCk0uI4khdBnl2WbbHWmUJRmAMl06/QKMdTLu3LY/U1MCH99q+lm11qLZG3K3QvQQz43sYoaaVtJ2W+sdQ7e+ZdZ+wu3OfdGDYcKFJuto/h3Nz12wC0JioLbUWAkpc4587dC83YSD//fZbp76rrmQz/WP5cRA8MveyFmRmhHqMNmJ8xmR/x0jbOvNf4PIQRAQ1Ln19HLEQuiNiMuo67EshOJusBBKMyA6yfTdSRjjXY+qPf8zbZhfPg92LjP7tDZP6w+Mce5rj9ytLd1FFpEDQNuhuo04yr6vTAFXeHzz/TOvh4Zqs06LhhojsGk/Ma/zfeA2KkkH5Q9RgwHYX1DJf5Yf4MyJA3n3ptl8/9sFbP/zIl6++3pA8a/jG/i53/sQN4JB171lxPD7+81/gz7uLgIRhN6JBJW7HlcLoauDvCUZziyZgZOMhdDeBMCDy00F7oAJ8NYVJsvn1fPhw59DRQ4c+Lb969ZXQ9He1gUhItFsK3I9H68sMENphp/c8tiAiRAcDZkuQenCvYCGYQvM2n0lCDFDmrKB/vbfnYQE+nPP2eOYOjSW5LgwwoMDUCHRpoPp2mdM5tS830BgCMy5FbLWmFnIPTjc/mhBBKG3YbebGILyM1tvujgKHccSBFud01rwBetfgP/e3nyf69PpgIlmVm9rN2GL9OVmqPuVH8KwE+HzP0DGj3D6fTBoimnb3B75O40F0KogDDDbynzPxw98Y7YjTmp5zM8fhkyHQy7D6q0Mo8SxpvDNF4FllxqEb3fn8/WufH558ggSIj20lhg81fxtY4bChAvMvimXQ3iiyfISC0EEodfRUG22EY5caYkjdA0V2c6/cfEB3513+/uw7jnnf7eGWvNE72ohgLESWqP0kBGRlHlmbOMlb8KZD8BNPxhXTfwo7wTBck21ZyFUtiJO+74yFcIDp3g+nny8acdRXWxeF+w0MZJ+w40g5O9s3xJqD4cgNNjs/GXZDlLiwlgyu5UJilY9wtzbnHOLA0Nh1s/N7zHJnVvLMYAIQm/DchdFDjRbcRt5x74v4bnTvWsep7WxEKxgrS/jCCXppmrWGtJutVxushDGA6rtOMLB5WZrrS8gyJHl47gRxo00ufn11W2vJXermQzWmqsk0rIQ8loes9tNgdmwBa0PaRriqArOWmu2BbuNGAQEmfYcdWXOgrwjoa4CqotoiEzm3o93sL+girvOTCMooJX1TLgAFv4NJl/WfP+M64xIjDrtyNdyjCCC0NuwUk4tQZDAsndkrIRDK6G6sP33VheZ3j5J081sXV9ZCLZGKM00v6f/YLZWlox1Uw6ONOmbbVkI6SsgtJ95yvaENVi+vXVbAeXW+vMEhpo4QIUHQcjbaiqTR3iIH1gMnmosAqsnUv5OSBxjfk8c59jXQbeR1k6rwlHn8NeV1by8KoMrZw3llLGJrX82NAZm39IykygoHE65p3eNDO0iJO20t2FZCFFiIXSIqkLn1nrybQ0rZhCdZJ66i3wkCOVZxjoAZ8GXp17+AyeZIKcntHbED+a2/mQeP9Jsi/a2PgPAboe87XDcFW2vOSLRs4Ww7yuzHe4hfmARFGa+y6HVxjVWchAmnG+OWcKQvwNGLWp7DZ/eYfoqNdRCYy16wHg2Tv0XP6xZwy+APP8BvH7d8cwaHtf2eYR2EUHobbRwGUkMwSus1ElvLAQroBw1yLg4fGUhWDf/wVNNwVdDjbEQ/IOdAVyAgRNh+3um1iA0tuU5yjJNdkxrWAVibcURSg5CQ1Xr8QOLyAGeBWH/19B/QvvimjwL1jxjJo5pOySMNvtDY02qaHuBZbsdNr+OjhlKQcLxHC5rICXzPcYuO4dGNQqAR25eTHCkiIEvEJdRb6PJQhhktl3lMlr3HKz8d9ecuz3WPOO7XjcWTYLgRW8iy0KIGmwshBIfpZ5agjD5UpPVkrXWkXKa3PxpP8ExEL5wb8tzpLvFDzwRFG4qbgvbEIT2AsoWniyEugrT48hTdpE7Q2aaTK3Nb5rXCS7D7hPT2nUZVefugdoy/pQzmxlrT2TxnlP5Rcy/KUuYygy2QXA0wREiBr5CLITeRlMMwfFkVtdFLbA3vW6epmff0jXnbw1bI3xyO8y6BRb9zXfntYTAm2Z15dnG9x2eYPzxjbUm6yi6k20NStLBLxDGnWdSTzNWei6Islw+hXtM0Zcr6SvMuhLGtH2tuOFtWwi5W813bO88EQNaxhAOLjeC5qn+wB2r3fSWN0wBmRXfABNYPvCt6YNkZf24sCu3nPdeeoPfAzEjZ/LYhCnMTO1HYlQI2M+CNU+btOtjfEZBdyIWQm+jSRAcFkJXuYyqCkwAtLvrHKwZD76OjVgxBK8shGzjkvPzd3G/+CDTqCTdWANh/cyTefoKz62bY4Ya4XC3ELQ2N+OUue3fBONHmhhCa2mdOVuMGAS0MwoyItG4llz/Pzu83tzcW5kt0OLz/YaZ/65xw5tfL3GcERYPwvXu+izOffwHhtXvxhYQxq8vPYezJw0yYgDGojr+pu5/YDnGEUHobVguo/B484TXVS6jqkLzj9Xyp3cXVsMzXw7/sduc5/UqhnDY6ZLr52XGjje4DnJJmWvcLrVlLS0E/wBz83QXhOIDxlJJmdf+teJGmnNXefi+hXvNk/lQL/oIRXooTivaZ75He2JiYaWfWvEDi0SH+8jNbfTW2kx+8/ZmpgyJ5af98/AffFzrvZQEnyKC0NuwLISgCJND3hVZRg01UO94IrT6mTHMAAAgAElEQVT83t1FjaOIyZeCUFMCOJ6UPd0g3SnPdgpC1GBH6qmPLATr5j90thFc8FwHED+y5XzlzDVma80qaIu4EWbr/vStNXx2p0kpPeH2lp9zx1P7iqJ9zvN7g2VJuMYPwAhEUAR8/bem4PJHm7O5470tnDAqgReumkRg/lZnQZnQ5Ygg9Dbqq4y5HhBscta7wkKoKnD+3u2C0EkLQeuWbi5XN1F7LiOrKM3RLA0/PxNY7mzX05pS890sCyF5tvOYp17+8aNMMNvW4NyXswkCw1o+aXsi3hIENytjz2emSO/EO503+7aIcCtOs9uN+6wjgpB6gnGBJU1vvj8gGC57G+or0f85mXXLlnLbm5uYntKPpy+fSnDhTlMPMniq99cSOkWPBJWVUn8BzgHsQD6wRGvdzb6JXkp9lWlXoJQZKt4VMYTeLAjf/gP2fg7Xf+vcZ1kFQZHtC0JtqWkPYqX1gnEbdTaG4N6mOTzOPDEX7PTcQydupBn8UpLuDDJnbzK9jrxxn0QPMZaNq4XQUGusg/jRpm22N7hXK1dkQ2ONU3C8oV8q3L4HQmOx2zXrMkrIKauhrKaBosp4DkY/wlVZ9zBt3W/5Xb+buGzJXwkN8nfOYhBB6DZ6KsvoPq313QBKqV8CfwRu7KG19C7qK42ZDaaKtCtcRq5ulVIPU6+6ks4KQu42k0Hj2sPfEoGEUe23SnCtQbDolwr7vzJPx60Vg4FxqwSGGaF2x8MgF0acbP77udcagLEQwLiN4kea75O7BY67su31W1gBcdfU0x//bdZxxfses3o8Ehprnu4tQbDiGh2xEIBCezhvfbef19ccIrO4ptmxEYkRfDT5aYZm38H1lW+h1B+BADMDOTyx89ldgtf0iCBorV3vYuE0OXiFdqmvMnnmYG48VisEX2JZCLGpPWshaN3xlMLKPPNkXVXgfLq1Asnxo52tpVs7b5MgDHbu8yb11G6H/5xigsWLn2p53JMgnHQ3zPmV5/NZT+CFe4AzzbahGgZO9vx+T7gGposPmL7/Y85qu7rYHaVMkz8r9dSyOOJGtvkxm12zYl8hK/cVsupAEVsPl2HXMGtYHL9dNIa0gVHEhAUSFRLo7D2U8Sd4/nQz02Hm9cZCGDxV0kq7kR6rQ1BK/Q24EigDFrTxvuuB6wGSk6UbYTNBCI7smjoESxCSpnnXV9+XWIJgbzQ3QOu7ekuVIxum/LCLIDgshHiHG6a21PNTufU5cLMQHKmnxQdaF4SsNaaC+MC3ngWnJN1cMyTauS8wxPx4IiTa+O+tJ/zsTWY7qAOCED/SDKixNcBHvzRWwRn3ef95C9fitKJ9EBjeaoWy3a753/ZcHvxiD3vzKwn0V0weEsPPF4zgnMmDGZEY0fp1kmeZQraVj5kWF4V7nG2qhW6hywRBKfUl4On/mru01h9qre8C7lJK/R9wC/AnT+fRWi8FlgJMmzZNLIk6V5dRV8UQCs0/+sSxsPVtc83gNv4h+xLXObu1ZR0TBK3N0BaA8hywHvKrikz8wLrJVxe3IQjZgGp+w7Oehnd/ZgKkntjxodlW5JgOpu5xAdeUU29xzTSyAsqWK8kb4kaYTKav7jUVzmc/0lzovCVygLMra9E+Y3m4Cd6homq+3JnHO+uz2JFTzvCEcB65eDIL0waYeIA3KAVzfw2vXwT/uwvQkmHUzXSZIGitT/Hyra8B/6UVQRDcqK90/qMOiTJZRkfiWmmLqgJT52DdwEoPmarS7sBdEDpyA6urMAFPaF4/UV1kgrhhjjGPVYXNK2ZdKT9sXCSuPvbowXDcVbDqcdOsbdJFzT9jtxtB6DfMWBGZqz0LgjXrwFviR8K298x/3+yN3geULSwhW/kopM433+FIiEh0prwW7WsW5P1yRx73/W83u/PMg8mYAZE8eOEkzpk8GH+/I/h/cuRC09Ji82vmtQhCt9IjaadKKVcH5E+AXT2xjl5JM5dRlOme2dBO3/uOUlVg2iNYgtCdcYSaElNwBx1PqXXNjnKdclZdaAa5WO2N28o0cq1BcOWM+01B2Ee3mIIyV7I3mOvNu938N7HaPVvYbQ6rIaVDX4f4Uca9VZlvAuUdcReBMzspMAx+8uiRPzREDDB/s/oq8z3iRlDbYOOej7Zz7Uvr0GjuPiuN7357Ip/96gTOOy7pyMQATNDeiqvEjWjdkhO6hJ6qQ/inUmqbUmoLsBBoo3Wj0Az3GAL4vhbBEoSYFPO6uwUheoj5vaOZRq5N2CpynL9XFxnrwBoE31a1cnmOZ0EICIILXzJre+PS5n+THR+YTJwxZ5pce3fBKD9sYhdH4jIC2LWs4wFlMAI4+kwzTa2j13YlIhHQxvLRdg7ogZz7+A+8sDKda+ak8vEv5vKzuakMjetgvKc1xp9nxKA195zQZfSIIGitf6q1Hq+1nqi1Pltr7cOhtcc49VXOGIIVoPR1HKGq0Nw8w/oZ33t3C4J18+qwIDgCysFRzV1GVUUOC8HFZdQarkVp7oT1g0vfMk/8r15o1qq1cRcNX2AGsCTPMq0YXF1fnjKMvMFy+Wx922w7aiEAXPKa6a7aGRzxlIwNnwNw6xeVFFbW8dySafzx7DSCA3zcVsI/EG5YDqf/y7fnFdpFKpV7E1o76hBcXEbg21oErZ0WglLmJtZdtQh2u6nobRKE0o593nIZDZzk5jJyxBCCwiAgtHWXUV2FydpqK24RPwIuesXECt68wrSwLj0EaeeY41abhsy1zs8cqSBED4GAEOOC6mhA2Qd8vj2XS59ZxQ3vm9qNvK1fA3DhohP5/ncLOGlM/667eFCY97USgs8QQehNNNQA2sVCcAiCL/v+1JYa90Z4gnkdO7RzFkJNCXz9V8fa26GuDNDOgOyRuIyUn+kkWp7jENAqE2gOc/TMD49vXRDKXOYgtEXqPPjJYyZz5/WLTcxj9BnmmPvYSDB/P+VvZhR0BD8/p5XQ0YByJ/l+TwE3v7qB7NIaBgw2/z2mBhxAR/TnihMnEBYknfOPRUQQehNWp1P3GIIvLQQrbbNJEFIcg+GPMON3xcPw/X2w/5v232u5WSIHmifjI3EZhcWbWoHGGnM+6+ZvuYvC+nl2GVUWwAc3OgSllbGTrky+BObfac6feoIzYN00NtIljlCwG2KGmC6mHcWKIxyJu+gI2ZJVyo2vrDcVxL+Yy58vNWVC/vYGVAcrlIXehch8b8LqQOpahwC+jSFYbhcrABubYqp0K/PaH5foTk0prH3W/J63Dcac0c77HYJgFXAdiSBEJDpdPhU50FhnfrcshDAPFkLRfnjlp6b1xMWvOdsyt8eJdxohGDq7+f7kWbD2P+baOz82QWFvewe5YwlCRwPKXlJTbyO9qIrS6ga01lTX27jzvS3EhgXx4jUziApxuG1CY81/HxGEYxoRBF/x6gUmz3vsWV13DXcLocll5EMLwRIEqxOma+ppRwVh7X+MiAVHO0c2tkVnBaHKIQjW8CCryAycAhce37wDaOkhePZU8/uSZaY621uUgpk3tNyffLzpG7RmKXz1F9PZdOFfO/ZdLAZNMe4m98lpbXCgoJLduRUcLq0hu7SWmoZGGm0am11T22ijqs5GdX0j2aW1ZJfVtDD+YsMCeelnM+gf5VJFHTFABKEPIILgCxrrTYfNuBHdKwhBXeAyqvLgMgIz+9ebCVkWDTWw6kkYcapxo+Rsbv8zNY4g8hFbCAXG525ZCOWHjesJXCyEOFOpbLHrE2Mx3LQS+o/r2PVawxoI8/kfTD+oi17xfpiMO6NOg9u2eV2g9+LKdO75eHvTTT400J+IkAAC/RT+/oqQAH/CggMIC/RnWkosw+KHMCwhnLiIIPyUQgHDEyOIj3Bbb0Si6cwqgnBMI4LgC6yhNTUdzIo50utYLiM/PyMKPrUQCgEFoQ6fuFUT0NHA8sZXTL7/3Nvg0EqTmllb7rkTqIW7heB6424PrY1bKyLRYckoE1i2rucqCPWVphV0YIipAI4Y4DsxAIhIMMJUlW/6/Yd3Ygi8Ul6Jgd2u+cenO3lm+UFOTevPr04ZSVJMGFGhAShfVLFb1mH8yLbfJ/RqRBB8geXD72iaZEdxtxDA9zMRqgqMX9wKgAaGGBdMRwTB1gA/PGoalQ2d7RSyvO0wtI1pX5YghMQYQejIUJq6crDVGUHwDzTb8sNmwIpfgLNmwxKG6kITfM7e0DXtES543ly3C26gtQ021qWXsOFQCbUNNmxaszOngu/3FLBkdgp3n5V25JXCrWGlwHqa7iYcM4gg+IJusxA8CEJwlG87nlo1CK5YmUbesuczKDsEZ/zLPOH2d2Tt5G1rXxCCo4wYhXRw1kNTdpQj9hE1yMQQlJ8RAeUWS6guMtcq3Ns1HTUHTOj0KSpqG/jvlhyW7y2krtGO1pqK2kY2ZZVS32gHIMBP4eenCAnw4w9njuVnc1N9YxG4M/sW4w4NCPL9uYWjBhEEX9DtFoJL51GrwZ2vqCpsKQj9hpkYiTtF++Hg9zDt6ub7M340T5MjHP0NowYZF1R7geWaElPtC+Zm3ZGZCFbbCisYHjnIjKAMDHWmnELzauXackCbwO1RxM6ccp5ZfoBPt+ZS02BjcEwo0aGB+PlBcIA/Vx4/lDkj4pme2o+I4G76JxwaK32F+gAiCL6grrtjCK4WghdjITtCVUHLPPzEsbDpFWdLC4tVT8LaZ2DUouZ+7qy15iZrVZoqZc6Zu63ta9eUOG86IdHG3dNYa27q7a7b0bYiwsVCyFhhzmPVCICLy6jY2e/oKBGEQ0XVPPjFbj7cnE14UADnThnMBdOSmDIkpmue+gXBDREEX1DfTRZCXaVxgbjeIIOjfNtryJPLyAq45m2HYfOd+3O3mm3GSjPQBEzGVc5mM/HKlQETTRqqrbH1Ai13QQBjJVjf9/VLjE/+1HtbftbqYxThaKcQNch8tjQThrgMd3dtcJe9EaKTm4tcN1PXaOPb3QV8tDmbz7fn4qcUN5wwnBvnDyMmTNwzQvciguALLAuhodrcELvKz2o1tnN9WjwSl1HeDmcPnrJMGP9T85TfWG9EzRtBsNtNTAAg4wenIORtNcHdwW75/AMmmKf94v2QMNrsq3P0ZbK+T02Jy6wHF0GIHGBcR+krmjeNc6Uy3+TrW9lRTamnWRB2mvN9ITFGVKsKHQHlrrEOtNY02jWNNk2D3U5BRR378ivZl1/J4dIaiivrKa6uZ1dOOeW1jfQLD+KymUO56cThzfP/BaEbEUHwBZYrB8wN1XJbdMV13CeIBUd1LPhqt8Fzi8xnlL9x62RvMoNJLNeT+xNzRKLxvedvd+4rOWjWo/wg/Qfn/qx1Zps0vfk5rMBy7lYjCNkb4dlFcP5zztoN10lmIY5YglWLUFNi1lzWSmPcqnyzbj9HNxZXF5ZrDMHPz4hG0T5jWU1d4vl8XtBos1PTYCOjqJpduRXszi3nYGE1mcXVZJZUU11v8/i5uPAg+oUHERsexKJxAzhz4kDmjIgn0F86yQg9iwiCL6hzEYSarhSEKs+C0FjrvWVSuNfcWM+4H6ZeDRtfgmW3NQ/4ulsIYKyEPBdBsNxFY86CnR+ZLJ+IBGN5RA4yU8ZciR8F/kHmc+MWw7JfG0siY6URBK1bdxmBs+NqRbYRNfdGb1bbCotIV0FwqwMIj4cDjt5Kg7xLOW2w2Vl1oIhPt+Xy9c58CivraLQ3L/ENCvAjNS6cIf3CmDMintiwQAL8/Qj0V8SEBTEiMYLhCeFEhkgXT+HoRATBF9S71AF0ZRzBkyCEuPQzCvCiAMqqGE6Za3z5aefCp3fA5jdhxEnmWLgHQes/DtY977wZ520zFsb0a40gHFppWkBnrYOkqS0/HxBkLIPcrbDhReOuCQxzClFdhZn+1poglDgEwd5obv5RA5ufvzK/+bpdj7sXhoXFQYFjSF8rYy1tds2u3HJWHSjmx/1FrD5YREVtI2FB/iwYnUhKfBjBAf6EBPoxOCaM0QMiSYkLI0Ce8oVejAiCL3C3EDpDY71xxbhmxli4DsexaGpwV9b8xme3wTd/h+OubD7fN2ezmQlgtVUO62fcRdvecc5N9hRk7T/OdBAtSTfziHO3mqf+5FnmfBkrYegc40pyT0O1GDARdn8Ch9ebcZRxw2H7+07rADwIguPv6TqTofxwS0GoKnDGJsAIZ0iM+by7hWC97jccQmNosNk5WFjFrtwKduaUs+lQKZuzSptcPkPjwjhzwkBOGpPICaMSCAnsvjbUgtCdiCD4gvpKTBM13XkLYfn9sPwBM5933m+au4HqK1u6o1prcJe1zpzL3gin/tm5P2ezSQF1zfSZeJHpyLntXfPak8so0SEWeducgjB0jlnfkBkmjjDMtEluET+wGDABNr1qKnjPfMAEide/YG72rn2MoHULAaAsq3kTOte2Fa5EDYLaUnIbI1i5IYv9BZXY7LCoyJ8pwI91Q7nr/m/JLKmmwWbcPwF+irRBUVwwNYnJyTHMTI1jUIwXaa+CcAzglSAopW7VWj/S3r4+S12lSXeszO28hbD/a+Nr/+6fxhVzzr/N0BVoxWXkCL5aaZcWB7832/Tlzn12uxGESRc1f+/IheYGvP9r8A92zllwJWEMoEyG0tC55indqsYdOge+/Qfs+9K4kVpr1TxgotnOusU8zVsFfTlbnNe0BCEwxKzFNYYQM9Rsy90Cy7VlpmbBxWVks2uyGmMYCvzkuV3kk4e/n8JfKSL9YYofbLEPY/TgSBaNH8Do/pGMHhDJsIRw34+EFIRegrcWwlWA+81/iYd9fZP6ChNErcztnIXQWOfI4b/B3HSX3QYvnAW/3mkqeD1lGQ2aYqqCD3wDoxY69x/8zmyzNzmbypUcNGt195sHhphYwoYXnaMz3QkKM5ZB3jaTWgpOQUiZA2jz9N9/nHmvJ4bONoPqRznSQBPTTJZS7hanBeJaDRsS7bR8SjLMuqsKWmYaudUgHCqq5tdvbeKn+cEMDYBbz57JccP6M6p/pOnxsyodPnufGy4+H4Z6iHcIQh+lzQiYUuoSpdTHQKpS6iOXn28AH5bH9nLqKs2TelBESwth+QOw5W3vzpOz2TzpJs2A0aeZBmkN1ebJHTzHEILCzMSu3Z86p5o11EDmGnMD1TbnOMecTWbrKZA66WKzbatIKzHNDJDPdROEwVONVdNQ3bq7CIzQpJ3jbAUd5JgTnLOlZQwBnC2w7XZTMxE71Iy3LM9qeovdrjmQYZrgvbennjvf3cLpj3zP7rwKUmefhx7/Uy6bM5KxA6OcDd9GLoTJl3dNUztB6MW0ZyGsBHKAeOABl/0VgBcTT/oI9ZUQk+wMYrqy+mnTQ36iFw3UMteYrTUMJWm6uUHu/cKkatZXthQEMDe4vZ+b3Pr4keY8tjoTg3j3WuM+GrXICI5/ECR4mAg25HjjkoluY+5v//FmAtihVSat0xKPwFBTiHZoZduC4IkBE00swfrOlgsMnIJQmWe+T8xQY4mVHWb53gLe22C2M6tX8HgQPLW+guLwPGYNj+fP54xjcEwocGXLa8YNh3Mf79g6BaEP0KYgaK0zgAygjRaVAnWVEBxh3DquFoLNkSJp91yg1IKsNaaVgtV73s8fhp8M+74wT9/a3tJlBOZm/8ntsOd/RhAOfmd8+cNPMi2orXhCzmbzlO+pXsHPD6762AhGa/RPA7QRn2EnNj82dLZDEDowcQxg4ETY+pYRs8Aw476ysATByjCKTaE2dAC1h77kimfXEBsWyAmjErjaPxy2w0d3nkdIdBfVgAhCH8CrpGmlVIVSqtzxU6uUsimlfNhis5dTX2kG1bhbCFX5gDZ9c7wZ9pK5tnnfHTA3+6oCZzWwJwshJhkSx5m202AEYPBUE6hNmWdcPNXFRhBaybsHHC6Zga0ft/z8tvqW7Z2PvwnOeaLjE7WsQPOB71p20wwxHU8rcvcB8NoeeHF7A1ENRfzmpGGs+v3JPHLxFCb3qwflT0hkz/UkEoRjAa8EQWsdqbWOcvyEAD8F/t21S+slaG0EwZOFYHXTBCjY3fZ5yrJMFW6S2+zc4ScDyuTrg2cLAUxA+dCP5jyHNzh7DqXOAzRsedP46dsShPaITTVP8dBSEMLjYcpl3rWqxgx3X/r9fu7b4rBIKrLJqg3mhpfX8bMX1nL+kyv5eHc1hUX5LP3QVBX/eXklRCXhpzS/mOGSDVSZb4LhflIUJgid4Yj+BWmtPwBO8vFaeidNrpyIlhZCRa7z98I9bZ+nKX7gZiGExxm//K7/mtetCsJppubg67+aQHLqCWb/4KmmcGzVE+Z1aymh3uDnZ1phg/PJ/gjYl1/BOY+v4O+f7OLFjeVkY+oeCmxhpBdWk1teS1CAH8GRsUSrGs5KbqAuNJEVd53BDWc7vld5tvOEJektW2UIgtBhvK1DOM/lpR8wDdCtvL1vYVUpt2chtCcIWWvNjdvTjXbkQhNfAM8uI3AGoDe/YfL3LUsjIBiSZ8KBb01cwapGPlIGTjb+/tjUDn9Ua817Gw7zhw+2ERbkz8s/m8G8kQnwxkzYtYwpo4bxv4tOcH5g+Q/wVT2j/bMhPpWEyGCocdz4yxyZRrZGU/k85YrOfS9BELy2EM52+VmEyTI6p6sW1atoGlrjiCE0VJmZwmAsBOVnfO/tuYwy1zQfKuOKa31BaxaCnz+MOBXQRgBcg7Mp88w2YYx3w2ba4qQ/wDX/67B7ZsOhEi5auorfvL2ZiUnRfHLrPCMG4BTBFjEER7Vy3nbnLN8ohyBYxWl524yVNsTN1SYIQofxykLQWrfSnEZoqra1LARwdDxNMBZCRH/jZrHaQnuiodYEfGfd7Pn4gIkQMcAUvrUmCGAC0FvfcrqLLFLnA3/pXPzAIqyf5z5LLtjtmu3Z5RwqruZwaTVrDpbw5c484iOCuPeccVw6I7l5E7iBrQmC4+/ZUO3sxxQSZfo3WcVpWWvNdsjMTn4xQRC8dRkNw1QlH49xFf0I3Ka1PtCFa+sdNFkIES49/C1ByDUppPGjYNt7pmDM0xN6zmawN7QMKFsoBSNPgY2vtC0Io0+HKZfDpEua7x802bSXGHt2x79fB9mSVcqfPtrOxkNO11lsWCC/PnUUP5ubSrinGcCWheBeFGdZCOC0EMBRnOYQhMzVjnbbbdRPCILgFd62rngNeBxY7Hh9MfA6II9lrjEEW7353YojVOSalND4UYA2swgGeogRZK4227bcHlOvNu0botoIngaFwzkeCq78A+HqT9r9Kh1Fa01ZTQM5ZbXkltXy+Y5c3libSVx4MH9bPJ4pQ2IZHGsGxLdJ9GC45I2WT/muguDasTV6sDOGkLnaBOJl5rAgdBpvBUFprV92ef2KUuqWrlhQr8NyGQVFmhYL4Mw0qsgxN3mrLXPhnpaCoDXs+NC0o25rsE7SNFiyzLdr95KskmrWpheTU1ZLTmkt2aU1ZJXUkFVSTZXLVDB/P8XP5qTyy1NGEtXRITCjT2+5ry0LIWczlOeYlhYzb+zgNxIEwRPeCsI3Sqk7gTcwLqOLgP8qpfoBaK29qLo6Rql3iSE0zQYuNY3qqosgcqDpu6/8PGcaHfwODq+DMx/svjV7yY7scp7+fj/LtuRgc0wHiw4NZGB0CMlxYcweEcfgmFAGxYTSPyqEoXFhxEcE+24BliAo/+aWUXSSo1hvhXkt8QNB8AneCoLVL/kGt/3XYARimM9W1Nuoc4kh+DmeimtLTf8dMDGEwBDzhOsp02j5AyZgPPmy7llvG9jsmq2Hy/hudwHf7clnw6FSwoP8uWZOCudPHcKQfqGEBXXjCA1LEKKTms9vsOYl7/jApNh2oiZCEAQn3v7rHqu1rnXdoZQKcd/XJ3ENKlsB45pSZ1FapKMVRMLolhZC5lrTZmLhX5unifYAe/Mq+PlrG9iTV4lSMGFwNL87bTSXzRhKdFgPzQAOCDG9lVzjB+C0FvZ+4Wj/7cUsaUEQ2sVbQVgJuPcK9rTv2ENryN5oUjbdB7uDsRACwx15+X7m99pSZ1Ga1agufpRpY+06IH75/SbVcmrPZvW+tyGLu97fRniwP/dfMImTxiTSL/wouMkq5czScsXKKLLVSf2BIPiQNgVBKTUAGAyEKqWmYOZEAkQBrUxBOcbI2QzPLDA37bMeapnNUl9h4gcWVrVyk4XgcG8kjDZZSK4zifd8Bgvuav75LqTBZievvJacsloOFVWTUVTF1sNlfLO7gBmp/fj3JVNIjOpZS6UFV37Ysj7BNZ4g8QNB8BntWQiLMJPRkgDXqGcF8PvOXlwpdTtwH5CgtS7s7Pm6hCLTaZP1zxvf9fzfNT9e5zajwOpnVJFjYgpWEZf1lFu4x7xn2W0mM2nGdT5bqtaag4VV5JTVUlBRR155LRnF1RwsqCK9qIrc8tqmGToAfgoGx4byy5NG8MuTRzYvFjta6OchPBUUZkSipkQsBEHwIe3NQ3gReFEp9VOt9bu+vLBSaghwKnDIl+f1OWWZZjtuMXzzN1N5PPUq53Gr06mFq4UQOdBpUViCsGsZfHanSZk8b2nLp98O0mCzszu3gk+35fDfLTmkF1U3Ox4bFkhKfDizhsWRFBvKwJhQkyXUL4yk2DCCAo5CEfCGqCQTdG4rVVcQhA7hbQxhvFJqnPtOrfW9nbj2Q8DvgA87cY6upzTT3LTPe8bUHCz7FQwY7xx8X+eYhWAREm1y48v9nfEDMEIR0d9UG0f0h6s/haQjm+ebWVzNA5/vZmdOBQcKK2mwafz9FLOGxXHtvGGMSIwgITKYhMjgjtcD9Bbm/qqnVyAIxxzeCkKly+8hwFnAziO9qFLqJ8BhrfVm1U6FqVLqeuB6gOTk5CO95JFTlgnRQ0y177lPwf0jIONHpyDUVzT3aYfEQM0W0+DOKkizGH6ycUFd8MIRt2vOK6/lsv+spriqnpmp/VgwJpHRAyI4YWQCcb6sATjamXB+T69AEI45vH7HoFsAABB9SURBVG1u5zpPGaXU/cBHbX1GKfUlMMDDobsw8YeFHo55uvZSYCnAtGnTur/ldmmm048dHm/iBZYbCRzjM10shFBHDKGuouWYycVPdmopxVX1XP6f1RRV1vHqdcczeUhM+x8SBEHwkiOtMgqjnWI0rfUpnvYrpSYAqYBlHSQBG5RSM7TWuZ4+02NobW7+1vQxpUxvolIXQXAffB8S46xNiPSkhx1dgqaoqp69eZX8/ZOdHCqu5oWrZ4gYCILgc7ztdroV50AcPyAR+MuRXFBrvdXxeevc6cC0ozLLqLbU3Nyjhzj3RQ+BMpc4eJ2HoLJFZBvziT1QWFnHBxsPs+ZgMaU1DZRVN5BXUUtptZmvEOTvx9NXTGXW8Lgj+TaCIAht4q2FcBYQC8wDYoBPtNbru2xVRwuWJRDjIggxQ5zdSW2N0FjjFlR2FQTvLIS16cU8/d0Bvt2dT6NdMywhnISIYFLiw5iaEsuIhAhG9o9g7MAo3/YKEgRBcMFbQTgHuA54D1Oc9rxS6hmt9WOdXYDWOqWz5+gyrFiBu4VQWwq15WaWMhyxhVBcVc8/PtnJ2+uziI8I5mdzUzl/ahIj+0e2+TlBEISuwFtBuBY4XmtdBaCU+n+YITmdFoSjmlIPgmBZC2WZzuZr7jEEizYshM+35/K7d7dQWdvIzScO5xcnjSQ0yENrDEEQhG7C63kIgM3ltQ1nG4tjl7JMM/jedZJXtCP1tTQTYh03cE8WQmBY837+LuzLr+CXb2xkZGIkD1w4iVFiEQiCcBTgrSA8D6xWSr3veH0u8GzXLOkooizTNFJzrZVwtRAsofAUQ4gc4HGKV22DjVte20hYUADPXjXt6OsdJAhCn8XbOoQHlVLfAnMxlsHVWuuNXbmwbsFuM0NsWmt/UJrZPKAMEJ5oWjKXHoK4EWafJwuhlfjBPz/dxa7cCp5fMl3EQBCEowqv6xC01huADV24lu5n48vw6Z1w27aWA97BWAEDxjff5+dnrIayzOazECwCgo2byUP84NOtObywMp2r56SwYIz04BEE4eiiG8dfHYUcWmXSRjNWQtpPmh9rqDFjGqM9tMuwitOsaWnu7asnnA+p85tebjhUwqNf7eXb3QWMGxTFnaeP8fEXEQRB6Dx9WxByt5mtJ0EoyzJbd5cRmKyjvZ+b9hTQPIYAvJN0Jyt3FVK4bg15ZbXszqugX3gQvzttNFfOSiE4QLKJBEE4+ui7gtBYDwW7zO8ZP7Q8XuqoRo72IAgxyWZmcnWRee1iIdQ32vn9e1sJC/YnuV8Yg2NDOX9qEpfOTCY8uO/+uQVBOPrpu3eogl1gb4D40WZ6WU1p86Ky9iwEgIKdoPzN7F8He/IqqLfZefDcSZw1cVAXfgFBEATf0kuno/iA3K1mO/MGQDvbUViUZYLy85wtZIlE/k5jHbikl249XAaYIfWCIAi9ib4tCAGhMPEiM+rS3W1UmmnmIft7GDBjWQhF+1vED7ZklREVEkByv74xcloQhGOHvisIedug/zjzhD94qgksu1LmoQbBImqQsR60rUWG0bbDZUxIiqa9wT+CIAhHG31TELSG3C0wYIJ5PXQ2ZG+E+irne0ozPQeUwVgNkY74gEsNQl2jjV255YwXd5EgCL2QvikIZZlQW+YiCHPA3ghZa81rWyOUH27dQgCTaQTNLIQ9uWa+scQPBEHojfRNQbACygMmmu2QGcYFlO6II1TmGndQdFLr57DEwsVCsALKEwfLNDNBEHoffVQQtgEK+qeZ1yFRMHCSiSNUF8O658x+T1XKFpY7yWWe8tbDpUSHBjKkX2jXrFsQBKEL6Zt1CLlbIG44BIU79w2dA6uehAdGg63evE4+vvVzxHgShDImDJaAsiAIvZM+aiFsdcYPLNLONS6iqVfDjT/A1Z+07FHkSnRzl1Fdo43duRUSUBYEodfS9yyEmlIozYCpVzXfP2Q6/GqL9+dxCyrvzq2QgLIgCL2avmch5G03WyugfKTEDDVupaTpgClIA5iYJIIgCELvpO9ZCNmOuT7uLqOOEhBk3EoOth0uIzo0kKRYCSgLgtA76VsWQlURrHgIBh0HEf19emoJKAuC0NvpW4Lw6e9MQdo5j3ucd3wkaK15bsVBduSUc9zQWJ+cUxAEoSfoOy6jnctg2zuw4C5n/UEnabDZ+eOH23l9zSEWpvXnxvnDfHJeQRCEnqBvCEJ1MSy7zcQN5t7W6dM12Ox8v6eAJ7/dz7qMEn6+YDi/OXU0fn7iLhIEoffSNwThi7uhphguf9dzO2svySyu5uVVGby34TCFlXXEhQfx8EWTOXfKYB8uVhAEoWfoG4Iw5zaTIjrwyFJN12eU8J/lB/jf9lz8lOKkMYlcMG0IJ45OINC/b4VhBEE4dukbghA/wvwcAS/8cJB7Pt5BdGggN8wfzpWzhjIwWlJLBUE49ugbgnCEvLU2k3s+3sHCtP48fPFkwoLkzyUIwrGL3OFaYdmWbO58bwvzRsbz2KVTCA7w7+klCYIgdCniAPfAjuxyfvXGJqYN7cfSK6aJGAiC0CcQQfDAR5uzAXj6iqmEBokYCILQNxBB8MAXO3I5flgcseFBPb0UQRCEbkMEwY0DBZXsL6ji1DTf9joSBEE42hFBcOOLHXkAnCKCIAhCH0MEwY0vd+aRNjCKwTFSayAIQt9CBMGFoso61meUiLtIEIQ+SY8IglLqHqXUYaXUJsfPGT2xDne+2pWPXSOCIAhCn6QnC9Me0lrf34PXb8EXO/IYFB3CuEFRPb0UQRCEbkdcRg5q6m0s31vAKWn9ZeqZIAh9kp4UhFuUUluUUs8ppVodNaaUul4ptU4pta6goKDLFrPqQBG1DXZOGSvuIkEQ+iZdJghKqS+VUts8/JwDPAkMByYDOcADrZ1Ha71Uaz1Naz0tISGhq5ZLelEVgLiLBEHos3RZDEFrfYo371NKPQMs66p1eEtBRR3+forYMKlOFgShb9JTWUYDXV4uBrb1xDpcKaysIz4iSMZgCoLQZ+mpLKN/KaUmAxpIB27ooXU0UVBRR0JkcE8vQxAEocfoEUHQWl/RE9dti4LKOuIjRBAEQei7SNqpg4KKOhJEEARB6MOIIAB2u6aosl5cRoIg9GlEEIDSmgYa7VoEQRCEPo0IAsZdBEgMQRCEPo0IAiblFBALQRCEPo0IAk4LQQRBEIS+jAgCIgiCIAggggCYGoSgAD8ig3uyG7ggCELPIoIAFDpqEKTttSAIfRkRBIyFIO4iQRD6OiIISB8jQRAEEEEAjCBIDYIgCH2dPi8IjTY7xdXStkIQBKHPC0JxVT1aS8qpIAhCnxeEfKsGIUImpQmC0Lfp84IgbSsEQRAMfV4QmqqUI0J6eCWCIAg9iwiCw0KIjxSXkSAIfRsRhIo6woP8CQuSthWCIPRt+rwgFMqkNEEQBEAEgYKKWhEEQRAERBCkbYUgCIIDEQRpWyEIggD0cUGoa7RRXttIggiCIAhC3xaEwsp6QIrSBEEQoI8LgozOFARBcNKnBaHQIQgSQxAEQYA+V42VX17L+xsPk1NWy8bMUkAsBEEQBOiDgvDIV3t5dfUhIoMDGBAdwrmTB9E/SvoYCYIg9ClB0Frzza58Fqb1Z+mV03p6OYIgCEcVfSqGsCevkuyyWk4ak9jTSxEEQTjq6FOC8M3ufABOHC2CIAiC4E6fEoSvd+UzduD/b+/uQqSqwziOf3+sSlbk9mJZulaC9kKU2m7Yu5oXWZFdlBVFIoUgURZFWF1EF0FF9EYRiG4ZhBUmJV0ElWJ1kZUJWZkplrlmrlLaK760Txfnv7RsM5nrHE+e8/vAMvv/z+yc5+GZnWfmf86cOYIhg7zPwMyst8o0hB1/7GbFhp+YeOrgokMxM/tfqkxD+GDtNv7sCiZ4ucjMrKbKNIQlX3UyaGB/Rrc0Fx2Kmdn/UmENQdJtktZI+kLSo3luq6srWPZ1JxePGky/psr0QDOzfVLI5xAkTQCmAGdGxE5Jua7jrNq0g22/7mKC9x+YmdVV1MvlmcDDEbETICI689zY0jWdSHDRSDcEM7N6imoIo4ALJS2XtExSW54bO2HQQKae3cLRPomdmVlduS0ZSXoHGFLjqvvTdo8ExgFtwKuSRkRE1LifGcAMgOHDh/cplqltLUxta+nT35qZVUVuDSEiJtW7TtJMYFFqAB9J6gKOAbbWuJ85wByA1tbWfzQMMzNrjKKWjF4HJgJIGgUMALYVFIuZmVHc2U7bgXZJnwO7gGm1lovMzOzAKaQhRMQu4MYitm1mZrX5U1pmZga4IZiZWeKGYGZmgBuCmZklOpgO7pG0FdjQxz8/hmoe2lrFvKuYM1Qz7yrmDPue94kRsddz9xxUDWF/SPokIlqLjuNAq2LeVcwZqpl3FXOG/PL2kpGZmQFuCGZmllSpIcwpOoCCVDHvKuYM1cy7ijlDTnlXZh+CmZn9uyq9QzAzs3/hhmBmZkBFGoKkSyWtkbRO0uyi48mDpBZJSyWtlvSFpFlp/ihJb0tamy6PLDrWRpPUJGmlpDfT+OT0bXxrJb0iaUDRMTaapGZJCyV9lWp+btlrLenO9Nj+XNICSYeUsdaS2iV1prNBd8/VrK0yT6fnts8kjd2fbZe+IUhqAp4FJgOnA9dLOr3YqHKxB7grIk4j+ya6W1Oes4F3I2Ik8G4al80sYHWP8SPAEynnn4CbC4kqX08Bb0XEqcBZZPmXttaShgK3A60RcQbQBFxHOWv9AnBpr7l6tZ0MjEw/M4Dn9mfDpW8IwDnAuohYn067/TIwpeCYGi4iNkfEp+n3X8ieIIaS5To/3Ww+cFUxEeZD0jDgcmBuGovsy5cWppuUMecjgIuAeZCdTj4itlPyWpOdrn+gpH7AocBmSljriHgP+LHXdL3aTgFejMyHQLOk4/u67So0hKHAxh7jjjRXWpJOAsYAy4HjImIzZE0DOLa4yHLxJHAP0JXGRwPbI2JPGpex3iPIvm72+bRUNlfSYZS41hGxCXgM+I6sEewAVlD+WnerV9uGPr9VoSGoxlxpj7WVdDjwGnBHRPxcdDx5knQF0BkRK3pO17hp2erdDxgLPBcRY4DfKNHyUC1pzXwKcDJwAnAY2XJJb2Wr9d409PFehYbQAbT0GA8Dvi8ollxJ6k/WDF6KiEVpekv3W8h02VlUfDk4H7hS0rdkS4ETyd4xNKdlBShnvTuAjohYnsYLyRpEmWs9CfgmIrZGxG5gEXAe5a91t3q1bejzWxUawsfAyHQ0wgCyHVGLC46p4dLa+TxgdUQ83uOqxcC09Ps04I0DHVteIuLeiBgWESeR1XVJRNwALAWuTjcrVc4AEfEDsFHSKWnqEuBLSlxrsqWicZIOTY/17pxLXese6tV2MXBTOtpoHLCje2mpLyrxSWVJl5G9cmwC2iPioYJDajhJFwDvA6v4ez39PrL9CK8Cw8n+qa6JiN47rA56ksYDd0fEFZJGkL1jOApYCdwYETuLjK/RJI0m25E+AFgPTCd7gVfaWkt6ELiW7Ii6lcAtZOvlpaq1pAXAeLJTXG8BHgBep0ZtU3N8huyopN+B6RHxSZ+3XYWGYGZme1eFJSMzM/sP3BDMzAxwQzAzs8QNwczMADcEMzNL3BDMzAxwQzAzs+QvDrWNZMxdajEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.compile(Adam(lr=.0001),loss='mean_squared_error')\n",
    "\n",
    "\n",
    "history = model.fit(train_data, train_label, epochs=15, batch_size=128, validation_split=0.1, verbose=1)\n",
    "\n",
    "model.save('rand_train.h5')\n",
    "plot_history(history)\n",
    "\n",
    "\n",
    "\n",
    "pred = model.predict(test_data, batch_size=128)\n",
    "score = np.sqrt(metrics.mean_squared_error(pred,test_label))\n",
    "print(\"Score (RMSE): {}\".format(score))\n",
    "\n",
    "chart_regression(pred.flatten(),test_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}