train_big.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "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": [
    "# number of spins\n",
    "l = 12\n",
    "\n",
    "#dimension constant\n",
    "j = 1\n",
    "\n",
    "#anisotropy constant\n",
    "delt = -1\n",
    "\n",
    "path_folder = 'matrix_pool/'\n",
    "\n",
    "if not os.path.exists(path_folder):\n",
    "    os.mkdir(path_folder)\n",
    "    \n",
    "h = -8.\n",
    "row = []\n",
    "while h <= 8:\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/2\n",
    "            else:\n",
    "                Hxxz[x,x] += j*h/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",
    "\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)/2\n",
    "                else:\n",
    "                    Hxxz_p[x,x] += (j*h)/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",
    "                        #if mag == 0:\n",
    "                            #print('[',x,D,']','k:',k,'s[x]:',s[x],'y:',y,'b:',b,'c:',c,'s[5]',s[5])\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",
    "    ex = '.npy'\n",
    "    Hxxz_pool = skimage.measure.block_reduce(Hxxz,(16,16), np.average)\n",
    "    np.save(os.path.abspath(path_folder + str(h)), Hxxz_pool)\n",
    "    row.append([str(h) + ex, psi_s[0], np.min(E_array)/l, h])\n",
    "            \n",
    "        \n",
    "    #print('h:',h)\n",
    "    #print('MAG_GS:', psi_s[0])\n",
    "    #print('rhoE_array:', np.min(E_array)/l)\n",
    "    h += .01\n",
    "    \n",
    "row = np.asarray(row)\n",
    "\n",
    "a_h = np.array(row[:,3]).astype(np.float) #h\n",
    "b_e = np.array(row[:,2]).astype(np.float) #E_GS/l\n",
    "c_m = np.array(row[:,1]).astype(np.float) #mag_GS\n",
    "\n",
    "plt.plot(a_h,b_e,'-')\n",
    "plt.title('GS-Energy density distribution (L='+str(l)+')')\n",
    "plt.xlabel(r'$h$')\n",
    "plt\n",
    "plt.ylabel(r'$E_{GS}/L$')\n",
    "plt.grid()\n",
    "plt.savefig('GS-Energy_big.png')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.plot(a_h,c_m)\n",
    "plt.title('GS-Magnetization distribution')\n",
    "plt.xlabel(r'$h$')\n",
    "plt.ylabel('GS-Mag')\n",
    "plt.grid()\n",
    "plt.savefig('GS-Mag_big.png')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "data_set = pd.DataFrame(row)\n",
    "data_set.to_csv(\"pool.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = pd.read_csv('pool.csv')\n",
    "\n",
    "path_folder = 'matrix_pool/'\n",
    "\n",
    "#writing the columns of the csv in an list\n",
    "X_1 = dataset.iloc[:,1].values\n",
    "Y_1 = dataset.iloc[:,3].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": 5,
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 1296 samples, validate on 145 samples\n",
      "Epoch 1/15\n",
      "1296/1296 [==============================] - 502s 387ms/step - loss: 3.1394 - val_loss: 3.2251\n",
      "Epoch 2/15\n",
      "1296/1296 [==============================] - 496s 383ms/step - loss: 2.9417 - val_loss: 0.2352\n",
      "Epoch 3/15\n",
      "1296/1296 [==============================] - 497s 383ms/step - loss: 0.4030 - val_loss: 0.0775\n",
      "Epoch 4/15\n",
      "1296/1296 [==============================] - 496s 383ms/step - loss: 0.1447 - val_loss: 0.1105\n",
      "Epoch 5/15\n",
      "1296/1296 [==============================] - 495s 382ms/step - loss: 0.0737 - val_loss: 0.0718\n",
      "Epoch 6/15\n",
      "1296/1296 [==============================] - 495s 382ms/step - loss: 0.0554 - val_loss: 0.0574\n",
      "Epoch 7/15\n",
      "1296/1296 [==============================] - 496s 383ms/step - loss: 0.0478 - val_loss: 0.0294\n",
      "Epoch 8/15\n",
      "1296/1296 [==============================] - 496s 382ms/step - loss: 0.0367 - val_loss: 0.0216\n",
      "Epoch 9/15\n",
      "1296/1296 [==============================] - 567s 437ms/step - loss: 0.0346 - val_loss: 0.0161\n",
      "Epoch 10/15\n",
      "1296/1296 [==============================] - 499s 385ms/step - loss: 0.0353 - val_loss: 0.0169\n",
      "Epoch 11/15\n",
      "1296/1296 [==============================] - 498s 385ms/step - loss: 0.0278 - val_loss: 0.0228\n",
      "Epoch 12/15\n",
      "1296/1296 [==============================] - 497s 384ms/step - loss: 0.0296 - val_loss: 0.0216\n",
      "Epoch 13/15\n",
      "1296/1296 [==============================] - 498s 384ms/step - loss: 0.0278 - val_loss: 0.0114\n",
      "Epoch 14/15\n",
      "1296/1296 [==============================] - 498s 384ms/step - loss: 0.0252 - val_loss: 0.0142\n",
      "Epoch 15/15\n",
      "1296/1296 [==============================] - 498s 385ms/step - loss: 0.0221 - val_loss: 0.0119\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOXVwPHfyR4SSIBAgIQETFCULRBAKFYRpS604Poqxb1KUfvi2kqttUWtIrWKVKp1geJSqa9YpWoRq0RcKgrIjgoiSGQLYQnZIMt5/5jJmD0Tcm9mxpyvn/uZO3d57kkMc+a591lEVTHGGGMAwgIdgDHGmOBhScEYY4yPJQVjjDE+lhSMMcb4WFIwxhjjY0nBGGOMjyUFY4wxPpYUjGmAiGwTkTMDHYcxrcmSgjHGGB9LCsY0k4hcJyJbRGS/iCwSkR7e7SIij4jIXhE5JCJrRaS/d9+5IrJRRA6LyLcicntgfwpj6mdJwZhmEJExwAPA/wDdge3AAu/uHwGnAscDicAlQL533zPAz1W1PdAfeLcVwzbGbxGBDsCYEDMJmKuqqwBE5NfAARHpBZQB7YG+wCequqnaeWXASSKyRlUPAAdaNWpj/GQ1BWOapwee2gEAqlqIpzaQoqrvAo8Bc4A9IvKkiHTwHnohcC6wXUTeE5GRrRy3MX6xpGBM8+wE0qveiEgc0Bn4FkBVZ6tqNtAPz22kX3q3f6qqE4CuwKvAS60ctzF+saRgTOMiRSSmasHzYX61iGSJSDRwP7BcVbeJyDAROVlEIoEioBSoEJEoEZkkIgmqWgYUABUB+4mMaYQlBWMa9yZQUm35IfBbYCGwC8gALvUe2wF4Cs/zgu14bis95N13ObBNRAqAKcBlrRS/Mc0iNsmOMcaYKlZTMMYY42NJwRhjjI8lBWOMMT6WFIwxxviEXI/mpKQk7dWrV6DDqKGoqIi4uLhAh+G3UIo3lGKF0Io3lGKF0Io3GGNduXLlPlXt0tRxIZcUevXqxYoVKwIdRg05OTmMHj060GH4LZTiDaVYIbTiDaVYIbTiDcZYRWR700fZ7SNjjDHVWFIwxhjjY0nBGGOMT8g9UzAm2JWVlZGbm0tpaWmgQ2lUQkICmzZtavrAIBFK8QYy1piYGFJTU4mMjDym8y0pGOOw3Nxc2rdvT69evRCRQIfToMOHD9O+fftAh+G3UIo3ULGqKvn5+eTm5tK7d+9jKsNuHxnjsNLSUjp37hzUCcF8P4kInTt3blEt1ZKCMS6whGACpaV/e20mKazbs44737mTAyU2C6IxxjSkzSSFrQe28sAHD/DVga8CHYoxrsrPzycrK4usrCy6detGSkqK7/3Ro0f9KuPqq6/miy++aPSYOXPm8MILLzgRMqeccgqrV692pKzaVqxYwc9//nPAc8/9hhtuIDMzk0GDBjV4zU8//ZT+/fuTmZnJLbfc4tt+6623csIJJzBw4EAuvPBCDh06BMDixYsZMmQIAwYMIDs7m/fff993zrRp00hNTSUxMbHGNWbOnMmJJ57IoEGDGDt2LDt27ABg9+7djBs3ztHfQXO0maSQlpAGwDeHvglwJMa4q3PnzqxevZrVq1czZcoUbrnlFt/7qKgowPPhWFlZ2WAZ8+bN44QTTmj0OjfeeCOTJk1yNHY3/OEPf2Dq1KkA/Otf/2LHjh1s2bKFOXPmcOONN9Z7zpQpU5g3bx6bN29mw4YNvP322wCcddZZbNiwgbVr19KrVy9mzpwJQNeuXXnjjTdYt24dc+fOZfLkyb6yJkyYwMcff1znGkOHDmXVqlWsWbOG8ePHM23aNAC6detGp06dWL58uaO/B3+1maSQnuiZVteSgmmrtmzZQv/+/ZkyZQpDhgxh9+7dTJ48maFDh9KvXz/uuece37FV39zLy8tJTExk2rRpDBo0iJEjR7J3714A7rrrLmbNmuU7ftq0aQwfPpwTTjiBjz76CPCMAXThhRcyaNAgJk6cyNChQ5usETz//PMMGDCA/v37c+eddwJQXl7Odddd59s+e/ZsAB555BFOOukkBg0axGWX1Z3M7tChQ3z++ef069cPgNdee40rrrjCF/Pu3bvJy8urcc6OHTsoLS1l2LBhiAiXX345r776KuBJChERnkabI0aMIDc3F4AhQ4bQvXt3AAYMGEBRURFlZWUAjBw5km7dutWJbcyYMcTGxtYpC+C8885zrBbWXG2mSWrHmI7ERcZZUjCt6uabwem7IllZ4P0sbraNGzcyb948nnjiCQ4fPsyMGTPo1KkT5eXlnH766Vx00UWcdNJJNc45dOgQp512GjNmzODWW29l7ty5vm+11akqn3zyCYsWLeKee+5h8eLF/PnPf6Zbt24sXLiQNWvWMGTIkEbjy83N5a677mLFihUkJCRw5pln8vrrr9OlSxfy8/NZt24dAAcPHgQ8t2C2b99OVFSUb1t1n3zyCQMHDvS9//bbb+nZs6fvfWpqKt9++y1dunRp8pja5s6dy5VXXlln+0svvUR2dnaz+gk888wznHPOOb73Q4cO5b777vP7fCe1mZpCWZnQKTyN7Qf9GhPKmO+ljIwMhg0b5nv/4osvMmTIEIYMGcKmTZvYuHFjnXNiY2N9H1jZ2dls27at3rIvuOCCOsd88MEHXHqpZwrrQYMG+b6xN2T58uWMGTOGpKQkIiMj+elPf8qyZcvIzMxky5Yt3HTTTbz11lskJCQA0K9fPy677DJeeOGFej+Ed+3aVeMDv77ph2u31vHnmOnTpxMfH+/72aqsW7eOu+66i0ceeaTRn7O6+fPns27dOm699Vbftq5du7Jz506/y3BSm6kpPPss7FifxsG93/DZ8TB4cKAjMm3BsX6jd0v14Zy3bNnCo48+yieffEJiYiKXXXZZve3bq55DAISHh1NeXl5v2dHR0XWOae4c8A0d37lzZz766CM++OADZs+ezcKFC3nyySd56623eO+993jttde47777WL9+PeHh4b7zYmNja/xMqamp7NixgxEjRgCemkmPHj1qXKvqmCq1j3nmmWdYsmQJ77zzTo3zvvnmGy644AKef/55/B3ef/HixcycOZP33nuvxu+5tLTUd2uptbWZmsKVV8Ko/ukURnxDdjZcdRXUUyM0ps2o6nXboUMHdu3axVtvveX4NU455RReeuklwPMtur6aSHUjRoxg6dKl5OfnU15ezoIFCzjttNPIy8tDVbn44ouZPn06q1atoqKigtzcXMaMGcMf//hH8vLyKC4urlHeiSeeyJYtW3zvx48fz7PPPgt4ajHJyck1ahIAPXv2JDo6mk8//RRV5bnnnmPChAkAvPHGGzz88MMsWrSImJgY3zkHDhxg3LhxPPTQQ76E05QVK1Zw4403smjRIpKSkmrs+/LLL+nfv79f5TitzSSFyEg4Z2Qa2m4vN/+yhBdfhD594O67obAw0NEZ0/qysrI46aST6N+/P9dddx2jRo1y/Br/+7//y7fffsvAgQP505/+RP/+/X23fuqTmprKPffcw+jRo8nKymLEiBGMGzeOHTt2cM4555CVlcV1113H/fffT3l5OT/96U8ZOHAgQ4YM4Y477qgztES/fv3Iy8ujqKgIgJ/85CekpKSQkZHBDTfcwJw5cwCoqKhg6NChvvMef/xxrrrqKjIzMznxxBMZO3Ys4GlxdfjwYc444wyysrJ8rZceffRRvv76a373u9+RlZXFqFGjyM/PBzzNWHv16kVBQQGpqam+ZwW3336770F8VlYW559/vu/6S5cuDVyzVFUNqSU7O1uP1bOrn1V+j36x7wvdulX1kktUQbVbN9Wnn1YtLz+2cpcuXXrMMQVCKMUbSrGqeuLduHFjoMPwS0FBgevXKCsr05KSElVV/fLLL7VXr15aVlZ2TGUda7wzZ87UefPmHdO5x6olv9vKykodNWqUHjx48JjLqO9vEFihfnzGtpmaAtTsq9C7NyxYAP/9L/TuDddeC0OGgLc5sjHGAYWFhYwaNYpBgwZx4YUX8te//tXXpLO1/OIXvzjmEUMDYe/evfzqV79qtEblpjbzoBnq76swYgR8+CG8/DLccQf86Edw7rnwxz9CrZZ5xphmSkxMZOXKlQGNITY2NiQ62VVJTk5m/PjxAbu+azUFEYkRkU9EZI2IbBCR6fUcEy0i/xCRLSKyXER6uRUPQEr7FASp01dBBC6+GDZt8iSDDz+EgQPh+uvB20/HGGPaBDdvHx0BxqjqICALOFtEaj+W/xlwQFUzgUeAB12Mh8jwSHq078H2Q/X3VYiOhttvhy1b4IYb4OmnITMTHngASkrcjMwYY4KDa0nB+2yjql1PpHep3Qh5AjDfu/4ycIa4POZwWkJak72ak5Jg9mzYsAHGjIE774S+feHvf4dGhosxxpiQJ9rMziXNKlwkHFgJZAJzVPWOWvvXA2eraq73/VfAyaq6r9Zxk4HJAMnJydkLFiw45pju3XgvXxR+wfPDn/f7nNWrE/nLXzLYvLk9ffsWcMMNXzFgwCHf/sLCQuLj4485ptYWSvGGUqzgiTclJYXMzMxAh9KkioqKGh29gl0oxRvoWLds2eIbwbXK6aefvlJVhzZwynf8aaLU0gVIBJYC/Wtt3wCkVnv/FdC5sbJa0iRVVfVXS36lUfdGaUVlRbPOq6hQnT9fNSXF04z1ggtUN2/27AvFZpOhIpRiVbUmqW4KpXgDHWvQN0lV1YNADnB2rV25QE8AEYkAEoD9bsaSlpDG0Yqj7C1q3hPksDC44gr48ku491546y1P66Rbb4WCgjbViMsEudGjR9fpnTxr1ixuuOGGRs+rqpHt3LmTiy66qMGyV6xY0Wg5s2bNqtGz+Nxzz613sLrmuv/++3nooYdaXE59SkpKOO2006ioqAA84xH16dOHPn36MH/+/HrP2b9/P2PHjqVPnz6MHTuWAwc8E3h9/vnnnHHGGURHR9eJ95prrqFr1651eiuvWbOGkSNHMmDAAH7yk59QUFAAwNtvv012drZvnoZ3333Xd86ZZ57pu6aT3Gx91EVEEr3rscCZwOe1DlsEVA0zeBHwrjejuaal8yq0awd33eV5GH3VVfDoo3D77YMcjNCYlpk4cSK1b7EuWLCAiRMn+nV+jx49ePnll4/5+rWTwptvvllngplgM3fuXC644ALCw8PZv38/06dPZ/ny5XzyySdMnz693g/fGTNmcMYZZ7B582bOOOMMZsyYAUCnTp2YOXMmt99+e51zrrrqKhYvXlxn+7XXXsuMGTNYt24d559/Pn/84x8BSEpK4l//+hfr1q1j/vz5XH755b5zLr/8cv7yl7849SvwcfMrbndgvve5Qhjwkqq+LiL34KnGLAKeAZ4TkS14agiXNlycM6r3VRieMvyYy+nWDZ58Erp3h3vuaU9JCQRo/CoTxG5efDOrdzs7dnZWtyxmnd3wSHsXXXQRd911F0eOHCE6Oppt27axc+dOTjnlFAoLC5kwYQIHDhzgyJEj3H///b5xfaps27aNH//4x6xfv56SkhKuvvpqNm7cyIknnkhJtWZ4119/PZ9++iklJSVcdNFFTJ8+ndmzZ7Nz505OP/10kpKSWLp0Kb169WLFihUkJSXx8MMPM3fuXMDzQXjzzTezbds2zjnnHE455RQ++ugjUlJSeO211xodEK5qAqHi4mIyMjKYO3cuHTt2ZPbs2TzxxBNERERw0kknsWDBAt577z1uuukmwDPa6bJly+oMh/HCCy/w97//HYC33nqLsWPH0qlTJwDGjh3L4sWL6yTV1157jZycHACuvPJKRo8ezYMPPkjXrl2JjY317avu1FNPrXeU2S+++IJTTz3Vd72zzjqLe++9l8HVRu7s168fpaWlvv+v48eP54c//CG/+c1vGvw9HQs3Wx+tVdXBqjpQVfur6j3e7Xd7EwKqWqqqF6tqpqoOV9WtbsVTpaqm4NQQ2iee6Hnd6nrkxvinc+fODB8+3PeNdMGCBVxyySWICDExMfzzn/9k1apVvPHGG9x2222NjmT6+OOP065dO9auXctvfvObGh3R/vCHP7BixQrWrl3Le++9x9q1a5k6dSo9evRg6dKlLF26tEZZK1euZN68eSxfvpyPP/6Yp556is8++wyAzZs3c+ONN7JhwwYSExNZuHBhoz/jFVdcwYMPPsjatWsZMGAA06d7ukHNmDGDzz77jLVr1/LEE08A8NBDDzFnzhxWr17N+++/XyfZHD16lK1bt/pGNvV3PoU9e/b4Jtbp3r27b/KhY9G/f38WLVoEwP/93//VGKW1ysKFCxk8eLBvNNqOHTty5MgR3xhLTmlzN8MTohNoH9Xescl2MjI8r1u2QBNDxZs2qLFv9G6quoU0YcIEFixY4Pt2rqrceeedLFu2DPB8AO7Zs6femcEAli1b5pvKcuDAgTUmrHnppZd48sknKS8vZ9euXWzcuLHG/to++OADzj//fN/w3RdccAHvv/8+48ePp3fv3mRlZQGNz9kAnkl/Dh48yGmnnQZ4vqVffPHFvhgnTZrEeeedx3nnnQfAqFGjuPXWW5k0aRIXXHABqampNcrbt29fjdtb9SVJl1vKM3fuXKZOnco999zD+PHjawyjDbBhwwbuuOMOlixZUmN71bwLnTt3diyWNjX2EXj+56YnpvNNgTNJoarlYbXReY0JuPPOO4933nmHVatWUVJS4pvx7IUXXiAvL4+VK1fy4YcfkpycXO8cCtXV94H49ddf89BDD/HOO++wdu1axo0b12Q5jdVIqr79QuNzNjTljTfe4MYbb2TlypVkZ2dTXl7OtGnTePrppykpKWHEiBF8/nnNR5sNzblQpb45F8AzHMWuXbsAz2Q+Xbt2PaaYAfr27cuSJUtYuXIlEydOJKPq26b3+ueffz7PPvtsje3gzrwLbS4pgH8d2PzVsSN06FDGV185UpwxjoiPj2f06NFcc801Ne6FHzp0iK5duxIZGcmyZcvYvr3x26innnqqb67g9evXs3btWgAKCgqIi4sjISGBPXv28O9//9t3Tvv27Tl8+HC9Zb366qsUFxdTVFTEP//5T374wx82+2dLSEigY8eOvP/++wA899xznHbaaVRWVrJjxw5OP/10Zs6cycGDByksLOSrr75iwIAB3HHHHQwdOrROUujYsSMVFRW+xHDWWWexZMkSDhw4wIEDB1iyZAlnnXVWnTjGjx/va5k0f/78Os9mmqPq1lNlZSX33XcfU6ZMATzTjo4bN44HHnigztDmqsru3bv9ntDHX20zKXRwdlrOHj1KrKZggs7EiRNZs2ZNjSkjJ02axIoVKxg6dCgvvfQSffv2bbSM66+/nsLCQgYOHMjMmTMZPtzTOGPQoEEMHjyYfv36cc0119T4wJo8eTLnnHMOp59+eo2yhgwZwlVXXcXw4cM5+eSTufbaa2s8SG2O+fPn88tf/pKBAweyevVq7r77bioqKrjssssYMGAAgwcP5pZbbiExMZFZs2bRv39/Bg0aVGNq0ep+9KMf8cEHHwCe1kO//e1vGTZsGMOGDePuu+/2PXS+9tprfU1yp02bxttvv02fPn14++23ffNW7969m759+/Lwww9z3333kZqa6mtiOnHiREaOHMkXX3xBamoqzzzzDOCZFvX444+nb9++9OjRg6uvvhqAxx57jC1btnDvvfeSlZVFVlaWL4GsXLmSESNGOD/qrD+dGYJpaWnnNVXV+5fdr/weLTxS2OKyVFXHjNmtvXs7UlSrCKUOYaEUq6p1XnOTm/GuWrVKL7vsMsfKa43f7dSpU/U///lPvfuCvvNasKlqlrqjoO4T/mORklLC9u1w9KgjxRljWtngwYM5/fTTfZ3XQkH//v0544wzHC+3TSaFlnZgq61HjxIqK6GJ27OmDVF3+2AaF1xzzTUhM7YSwHXXXVfv9pb+7bXppODUc4WUFE+HHnuuYABiYmLIz8+3xGBanaqSn59PTEzMMZfR5vopAPRo34NwCXesppCS4mm1YC2QDHiaNObm5pKXlxfoUBpVWlraog+P1hZK8QYy1piYmDp9MZqjTSaFiLAIUjqkONZXoWPHo8TFWU3BeERGRtK7d+9Ah9GknJycY279EwihFG8oxVpbm7x9BM72VRDxdGKzmoIxJtS16aTgZF+FjAyrKRhjQl/bTQod0sgtyKWi0pkmaJmZnkHxQqhFmzHG1NFmk0J6YjpllWXsKdrjSHmZmZ5+CvUMpmiMMSGjzSYFp/sqVB8t1RhjQlWbTwpOPVew0VKNMd8HbT4pONdXAaKirAWSMSa0tdmk0CG6A4kxiY4lhfBwOO44qykYY0Jbm00K4G2Wesi5ZqmZmZYUjDGhrc0nBadqCuB52PzVV2BD3hhjQlWbTgrpCemOJoXMTCgqgj3OtHI1xphW16aTQlpCGgdKD3D4SN2pA4+FtUAyxoQ615KCiPQUkaUisklENojITfUcM1pEDonIau9yt1vx1MetvgrWAskYE6rcHCW1HLhNVVeJSHtgpYi8raobax33vqr+2MU4GlQ9KfTr2q/F5aWne1ohWU3BGBOqXKspqOouVV3lXT8MbAJS3LresUhP8EzL6VRNISrKkxispmCMCVWtMp+CiPQCBgPL69k9UkTWADuB21V1Qz3nTwYmAyQnJ5OTk+NIXBVaQbiE8/669zmh8IRjLqewsNAXU6dOA/nsswhyclY5EqMbqscb7EIpVgiteEMpVgiteEMp1jpU1dUFiAdWAhfUs68DEO9dPxfY3FR52dnZ6qRes3rppIWTWlTG0qVLfevXX6/asWMLg3JZ9XiDXSjFqhpa8YZSrKqhFW8wxgqsUD8+s11tfSQikcBC4AVVfaWehFSgqoXe9TeBSBFJcjOm2tzoq3DgAOzf71iRxhjTatxsfSTAM8AmVX24gWO6eY9DRIZ748l3K6b6uNFXAey5gjEmNLn5TGEUcDmwTkRWe7fdCaQBqOoTwEXA9SJSDpQAl3qrOa0mLeG7yXbCw8JbXF71vgrDhrW4OGOMaVWuJQVV/QCQJo55DHjMrRj8kZaQRoVWsPPwTnom9Gxxeccd53m1moIxJhS16R7N4HwHtthYzzDa1lfBGBOK2nxScLqvAthoqcaY0NXmk0LVLSM3Rks1xphQ0+aTQnxUPJ1iOzk+r8Lu3VBY6FiRxhjTKtp8UgDn+ypUtUDautWxIo0xplVYUsD5vgpVo6XacwVjTKixpIA7vZrBkoIxJvRYUsCTFA4dOcSh0kOOlJeQAF262MNmY0zosaSAO81SMzKspmCMCT2WFHC+AxtYXwVjTGiypMB3ScHJZqkZGbBjBxw54liRxhjjOksKQHJ8MlHhUY7XFFTh668dK9IYY1xnSQEIkzB6dujpSl8Fu4VkjAkllhS83GqWai2QjDGhxJKCV1pCmqPPFJKSoEMHqykYY0KLJQWvtIQ0dh7eSVlFmSPliXhuIVlNwRgTSiwpeKUnpFOplew8vNOxMq2vgjEm1FhS8HKrr8LXX0N5uWNFGmOMqywpeLnRVyEz05MQduxwrEhjjHGVJQUvtybbAbuFZIwJHZYUvNpFtqNLuy7WV8EY06a5lhREpKeILBWRTSKyQURuqucYEZHZIrJFRNaKyBC34vGH030VuneHmBhrgWSMCR1u1hTKgdtU9URgBHCjiJxU65hzgD7eZTLwuIvxNMnpvgphYdYCyRgTWlxLCqq6S1VXedcPA5uAlFqHTQCeVY+PgUQR6e5WTE2pqimoqmNl2mipxphQ0irPFESkFzAYWF5rVwpQvW1OLnUTR6tJT0in8GghB0sPOlZmRoZnrubKSseKNMYY10S4fQERiQcWAjerakHt3fWcUudruohMxnN7ieTkZHJycpwOE4CCPE94C99ZSGZ8pt/nFRYWNhhTZWUPSkqOZ+HCj+jS5agTYbZYY/EGm1CKFUIr3lCKFUIr3lCKtQ5VdW0BIoG3gFsb2P9XYGK1918A3RsrMzs7W93ySe4nyu/R1z5/rVnnLV26tMF9S5aogmpOTguDc1Bj8QabUIpVNbTiDaVYVUMr3mCMFVihfnxuu9n6SIBngE2q+nADhy0CrvC2QhoBHFLVXW7F1BQ3ejXbaKnGmFDi5u2jUcDlwDoRWe3ddieQBqCqTwBvAucCW4Bi4GoX42lS17iuRIdHO5oU0tIgIsIeNhtjQoNrSUFVP6D+ZwbVj1HgRrdiaC4RcbyvQkQE9O5tNQVjTGiwHs21ON1XAayvgjEmdFhSqCU9Id3RmgJ811fBwe4PxhjjCksKtaQlpLHr8C6OVjjXfDQjAwoKID/fsSKNMcYVlhRqSUtIQ1FyC3IdK9MGxjPGhApLCrW4NdkOWFIwxgQ/Swq1pCemA84mhd69PXM2WwskY0yws6RQS2qHVMDZpBAdDT17Wk3BGBP8LCnUEhMRQ3JcMtsPOtss1UZLNcaEAksK9UhLSOObAmebpWZk2O0jY0zws6RQj/REd/oq5OV5mqYaY0yw8ispiEiGiER710eLyFQRSXQ3tMBJ6+DOZDtgtQVjTHDzt6awEKgQkUw8I5/2Bv7uWlQBlpaQRnFZMfklzvU2qxot1Z4rGGOCmb9JoVJVy4HzgVmqegsQsGkz3ebmENqWFIwxwczfpFAmIhOBK4HXvdsi3Qkp8NzoqxAfD8nJdvvIGBPc/E0KVwMjgT+o6tci0ht43r2wAsuNmgJYs1RjTPDzaz4FVd0ITAUQkY5Ae1Wd4WZggdQ5tjOxEbGu9FX4z38cLdIYYxzlb+ujHBHpICKdgDXAPBFpaIrNkOebbMeFvgrffgslJY4Wa4wxjvH39lGCqhYAFwDzVDUbONO9sALPrb4KAFu3OlqsMcY4xt+kECEi3YH/4bsHzd9rVX0VnGSjpRpjgp2/SeEe4C3gK1X9VESOAza7F1bgpSWksbtwN6XlpY6VWdUs1VogGWOClV9JQVX/T1UHqur13vdbVfVCd0MLrKpmqU5OttOpE3TsaDUFY0zw8vdBc6qI/FNE9orIHhFZKCKpbgcXSNYs1RjTFvl7+2gesAjoAaQA//Jua5CIzPUmkfUN7B8tIodEZLV3ubs5gbvNraRgo6UaY4KZv0mhi6rOU9Vy7/I3oEsOIHReAAAYkklEQVQT5/wNOLuJY95X1Szvco+fsbSK1A6pCOJKX4Xt26GszNFijTHGEf4mhX0icpmIhHuXy4BGR4tT1WXA/hZHGCBR4VF0b9/dlZpCRYUnMRhjTLDxq0czcA3wGPAIoMBHeIa+aKmRIrIG2Ancrqob6jtIRCYDkwGSk5PJyclx4NJNSySRNdvWNHm9wsJCv2M6fDgBGMwrr6xl+PDA5MzmxBtooRQrhFa8oRQrhFa8oRRrHap6TAtwsx/H9ALWN7CvAxDvXT8X2OzPdbOzs7W1/M///Y/2md2nyeOWLl3qd5m7dqmC6p//3ILAWqg58QZaKMWqGlrxhlKsqqEVbzDGCqxQPz5jWzLz2q0tTEYFqlroXX8TiBSRpJaU6TQ3JttJToa4OHvYbIwJTi1JCtKSC4tINxER7/pwbyzOzWrjgPTEdI5UHCGvOM+xMkU8zxWsWaoxJhj5+0yhPo1+fRaRF4HRQJKI5AK/wzsHg6o+AVwEXC8i5UAJcKk6+ZXcAdWbpXaN6+pYuZmZsHGjY8UZY4xjGk0KInKY+j/8BYht7FxVndjE/sfwPLwOWlVJYfvB7QztMdSxcjMy4PXXPa2QwsMdK9YYY1qs0aSgqu1bK5Bg5Gav5qNHPcNop6U5WrQxxrRIS54pfO91jOlIfFS8jZZqjGkzLCk0ws3JdsBaIBljgo8lhSakJaQ5PtRFaipERVlNwRgTfCwpNMGNyXbCw+G446ymYIwJPpYUmpCemE5ecR4lZc5OrGx9FYwxwciSQhOqWiDtKNjhaLlV8yoEV88MY0xbZ0mhCdX7KjgpIwOKimDvXkeLNcaYFrGk0IT0BM+0nNYs1RjTFlhSaEKP9j0IkzBLCsaYNsGSQhMiwyPp0b6H430V0tMhLMxaIBljgoslBT+40VchKsqTGKymYIwJJpYU/JCekO747SP4rgWSMcYEC0sKfkhLSGNHwQ4qtdLRcjMy7PaRMSa4WFLwQ1pCGkcrjrK3yNn2o5mZsH8/HDjgaLHGGHPMLCn4wa2+ClUtkKy2YIwJFpYU/OBWX4Wq0VLtuYIxJlhYUvCDW5PtHHec59WSgjEmWFhS8ENCTAIdojuw/ZCzt4/atYMePez2kTEmeFhS8FNagvNDaIM1SzXGBBdLCn6yvgrGmLbAkoKf3KopZGTA7t2eEVONMSbQXEsKIjJXRPaKyPoG9ouIzBaRLSKyVkSGuBWLE9IS0sgvyafoqLOf3tYs1RgTTNysKfwNOLuR/ecAfbzLZOBxF2NpMbdaINloqcaYYOJaUlDVZcD+Rg6ZADyrHh8DiSLS3a14WsrtvgpWUzDGBIOIAF47Bag+x2Wud9uu2geKyGQ8tQmSk5PJyclpjfhq2FvqGeLi7U/fJjo3usa+wsLCFsWUkPADli3bx7BhX7YkRL+1NN7WFEqxQmjFG0qxQmjFG0qx1hbIpCD1bKt3xmJVfRJ4EmDo0KE6evRoF8OqX3llOeGfhBOTHEPt6+fk5NTZ1hwnngjFxT0YPbpHy4L0U0vjbU2hFCuEVryhFCuEVryhFGttgWx9lAv0rPY+FdgZoFiaFBEWQUqHFNdaINntI2NMMAhkUlgEXOFthTQCOKSqdW4dBRM3+yp88w0cOeJ40cYY0yyu3T4SkReB0UCSiOQCvwMiAVT1CeBN4FxgC1AMXO1WLE5JS0jjox0fOV5uRgaowrZtcMIJjhdvjDF+cy0pqOrEJvYrcKNb13dDWkIa/9jwDyoqKwgPC3es3OrNUi0pGGMCyXo0N0N6QjrlleXsLtztaLnWV8EYEywsKTSDWx3YkpKgfXt72GyMCTxLCs3gVlIQsYHxjDHBwZJCM/im5XR4XgWwpGCMCQ6WFJqhfXR7OsZ0dK2vwrZtUF7ueNHGGOM3SwrN5OZkO2VlsGNH08caY4xbLCk0U1pCmmu3j8BuIRljAsuSQjO5OdkOWAskY0xgWVJopvSEdA6WHqTgSIGj5fboATExVlMwxgSWJYVmqmqBtOOQszf/w8LguOOspmCMCSxLCs1kzVKNMd9nlhSaya0ObOBJCl995RkczxhjAsGSQjN1b9+dyLBI1x42l5TArqAeQNwY831mSaGZwiSM1A6prtUUwG4hGWMCx5LCMbC+CsaY7ytLCsfArb4KaWkQEWEtkIwxgWNJ4RikJ6TzbcG3lFc6O1BRRAT06mU1BWNM4FhSOAZpCWlUaAW7Djv/RNiapRpjAsmSwjFws69CRoY1SzXGBI4lhWOQnpgOuNdX4dAhyM93vGhjjGmSJYVj0LNDT8CdpGAD4xljAsmSwjGIi4qjc2xn66tgjPnecTUpiMjZIvKFiGwRkWn17L9KRPJEZLV3udbNeJzkVl+F3r09czZbUjDGBEKEWwWLSDgwBxgL5AKfisgiVd1Y69B/qOov3IrDLemJ6WzZ7/wnd0wMpKba7SNjTGC4WVMYDmxR1a2qehRYAExw8XqtKq2DOx3YwJqlGmMCx7WaApACVJ90IBc4uZ7jLhSRU4EvgVtUtc5EBSIyGZgMkJycTE5OjvPRNlPZvjIKjhTw+n9eh1Icjaldu+P57LMkcnI+cqzM6goLC4Pid+iPUIoVQiveUIoVQiveUIq1DlV1ZQEuBp6u9v5y4M+1jukMRHvXpwDvNlVudna2BoOX1r+k/B5ds3uNLl261NGyZ8xQBdVDhxwt1sfpeN0USrGqhla8oRSramjFG4yxAivUj89uN28f5QI9q71PBXbWSkj5qnrE+/YpINvFeBzldl8FsOcKxpjW52ZS+BToIyK9RSQKuBRYVP0AEele7e14YJOL8TjK7cl2wJ4rGGNan2tJQVXLgV8Ab+H5sH9JVTeIyD0iMt572FQR2SAia4CpwFVuxeO0rnFdiQqPYvtB55ulHnec59VqCsaY1ubmg2ZU9U3gzVrb7q62/mvg127G4JYwCaNnh558U/CN58mIg9q3h+RkeP11mDIFEhOdLd8YYxpiPZpbID0x3bVmqdOnw/LlMHgw/Pe/rlzCGGPqsKTQAm5NtgPw85/D++971n/4Q3jgAaisdOVSxhjjY0mhBdI6pLHz8E7HJ9upMmIEfPYZXHgh3Hkn/OhHsMv5KRyMMcbHkkILpCWkUamV7Du6z7VrJCbCggXw1FPw0UcwaBD8+9+uXc4Y08ZZUmiBqr4Ke0r3uHodEbj2WlixArp1g3PPhdtug6NHXb2sMaYNsqTQAlV9FfYe2dsq1zvpJM/D5xtugIcfhlGjrC+DMcZZlhRaoGqyHbdrCtXFxsKcOfDKK55+DIMHwwsvtNrljTHfc5YUWiA2MpaucV3Zc6T1kkKV88+H1ashKwsuuwyuugoKC1s9DGPM94wlhRZKS0hjb2nr3D6qc+00WLoU7r4bnn0WsrM9rZWMMeZYWVJoobSENLYXb2fX4cC0FY2I8HR0e/ddKCryNGN99FHwDDxrjDHNY0mhhcYfP569R/aSPiudq169ijW71wQkjtGjPbeTzjoLbr4Zxo+Hfe61lDXGfE9ZUmihK7Ou5LnhzzFl6BRe3vgyWX/N4sxnz+TNzW9Sqa3bBTkpCV57zVNTWLLE06chVOf5MMYEhiUFB6TEpjD7nNnsuGUHD575IJ/v+5xxfx9H/7/056mVT1FSVtJqsYjA1Knw8ccQHw9jxsBvfwvl7nS6NsZ8z1hScFDH2I78atSv+Pqmr3n+/OeJjYxl8uuTSZuVxu+W/o49ha3XSmnwYFi5Eq68Eu67z3N76Rt3hmkyxnyPWFJwQWR4JJMGTmLFdSvIuTKHH/T8Afcuu5e0WWn87LWfsX7v+laJIz4e5s3z9GNYu9ZzO+mVV1rl0saYEOXqfAptnYhwWq/TOK3XaXyZ/yWPfvwo81bPY+7quZyVcRa3jryVsceNRURcjeOnP4WTT4ZLL/UMrpecPIIf/MDThHXoUM9rUpKrIRhjQoTVFFrJ8Z2PZ864Oey4ZQd/GPMH1uxZw1nPn8WAxwcw97O5lJaXunr9jAz48EN47DHo1+8Q69fDXXfB2WdDly7Qq5cnYdx/v+chdX6+q+EYY4KU1RRaWed2nbnzh3dy28jb+MeGf/Cn//6Jny36Gb9+59f8YtgvmDJ0Cl3iujS73EqtJL84n92Fu+sse4r2+NbzSvKIPiua47sdz9CYVMIKUynZk0reV6l8ujaVVxanQnFnQEhP/64mUfXaqZPzvxNjTPCwpBAg0RHRXDHoCi4feDlLty3lT//9E3fn3M39H9zPFQOv4OYRN9M3qS8FRwrq/XCvvewt2kuFVtS5TmxELN3bd6dbfDeO73w8P+j5A7785kuKyop4P/8ddh7e6Wk62xPPMg6iwmJorykUF6by1p5UFr6bCq+mQkEq3eNSye6TysiBXRk+LIwhQyxRGPN9YkkhwESEMb3HMKb3GDblbWLWx7N4du2zPLnqSWIiYuq9rRQRFkFyXDLd4rvRo30PhnQfQrf4bnSL7+bbXrXER8XXeWaRk5PD6NGjASivLGdP4R5yC3JrLoe9rykfciTzW8oqywDYBbwOvF4SCW/0gBdTiatIJTmuO+2j4omLiiM+uh3tY+LoENuOxHZxJMa1o1P7ODq1b0fnDu3omhhHl0TPa1REpMu/YWNMc1hSCCIndjmRv/7kr9w35j6e+ewZ8ovza3zAd4vvRnJ8Mp1iOxEmzjwOigiLIKVDCikdUjiZk+s9plIrySvKq5E0Nu/NZcM337J1Xy57S1exLWwXlRFFIN7xNcq8S0ETAVREIOVxhFe2I7wyjkhtR6S0I1riiAlvh5aGk/zua8SHJ9E+PImEiC4kRCaRGJVEp+guJEZ3IiY6nMhIiIrC91p9vb7X6GjPiLORkZ6+HcYYD1eTgoicDTwKhANPq+qMWvujgWeBbCAfuERVt7kZUyjoEteFaadMC3QYPmESRnJ8MsnxyWT3yG7wOFWltLyUoqPF7CsoIu9gMXmHithfUMz+wmIOFBZxsLiYguIiCkqLOVxaRNHRYorLiygpK6aksoijlcUc1SKKKaY8bB8VkYfYXv4fCC+ECjzLkeoXFSjpCMVdoDjJsxRVW69v+9F4wJMJRCAmpuVLbKwn0Wze3I2tWyEszFO2U6+e3+93Y1r589rUMevXd6a0tP6k2VhCjYz0xGW+n1xLCiISDswBxgK5wKciskhVN1Y77GfAAVXNFJFLgQeBS9yKybhLRIiNjCU2MpakuM707d7yMqtudRUdKWV3QT57Du9jd0Eeewv3sbdwH3nFeeQX7yO/dB/7S/dx4MhWDh79hENl+yjXsnrLjCCKKOIJI5IwjSJMo6AykiMaxZGKSAoqo6AiCq2IRMuj0HLPa2VZFBVHIqk8HEXF0SjKj3q2UxEJFZ4yqIiE99+CygjnFg33Rq7f1cTwvorWWm/GcVIJYSshrNyzhJd9tx5W1uh2iSgnPLKM8MhywiLKCY8qIyyinLDwck/5KiigKp7EreJb/+617n7fvkrP+b7jK8PQigjCKzchldFIRTRSGYNURhNWEUNYped9WGU0YVXbK2MI12jvejThYWGI1E26tderLw1tb2wfouzf349uyTW3h4fXv96c96ecAmee2eJ/Vo1ys6YwHNiiqlsBRGQBMAGonhQmAL/3rr8MPCYiompjfJqa4qJjyOiSQkaXFL+OV1UKjhSwr3ifb8krzvO8FuVRXFZMWWUZRyuO+l6PVhylrKLaemUxRysO1rO97vFK2/qTVaDcu1QRjUA0HPHWwtSXnNT7+6m+TrXE5b+WjNYSppGEqydxhGsMYRqNaKQ3pso6C9Lw+6r1+l5r/FwqoGGA51Wqrfv21bNNNQwqBMprHjc653rOPPNXLfgtNM3NpJAC7Kj2Phfq3LT2HaOq5SJyCOgM2PiepkVEhISYBBJiEsjolOHqtVSVCq3gaMVRct7LYeSokZRXlju6VDUWEKTGetXPWn3d3+M2bNjA4IGDiQyLJCIsgoiwCCLDq603Y3u4hB9zJ0xVT6Ko/grU2Faplbz73rsMGzGM0vJSjlQc8byWH/G9r75ee1+9x1aUUlZRRpiEOb5s27aN9PR0KrXSF7+q5/WYtqFUVlZyTmavY/odN4ebSaG+v5DaXw38OQYRmQxMBkhOTiYnyIb+LCwsDLqYGhNK8YZSrACVpZWsWd684dMFIdL7X2vKjs0mfld8jW3l3v+CUilsWrGp3l0R3v/iiGv4/AharWlNYVIh8RJf/ydcS+Tj+r8HN39FuXhavldJBXY2cEyuiEQACcD+2gWp6pPAkwBDhw7VquaUwaJ6E89QEErxhlKsEFrxhlKsEFrxhlKstbnZhuBToI+I9BaRKOBSYFGtYxYBV3rXLwLetecJxhgTOK7VFLzPCH4BvIWnSepcVd0gIvcAK1R1EfAM8JyIbMFTQ7jUrXiMMcY0zdU7bKr6JvBmrW13V1svBS52MwZjjDH+sy4oxhhjfCwpGGOM8bGkYIwxxseSgjHGGB8JtRagIpIHbA90HLUkEVq9sEMp3lCKFUIr3lCKFUIr3mCMNV1Vm5zBK+SSQjASkRWqOjTQcfgrlOINpVghtOINpVghtOINpVhrs9tHxhhjfCwpGGOM8bGk4IwnAx1AM4VSvKEUK4RWvKEUK4RWvKEUaw32TMEYY4yP1RSMMcb4WFIwxhjjY0mhBUSkp4gsFZFNIrJBRG4KdExNEZFwEflMRF4PdCxNEZFEEXlZRD73/o5HBjqmhojILd6/gfUi8qKIxAQ6pupEZK6I7BWR9dW2dRKRt0Vks/e1YyBjrNJArH/0/h2sFZF/ikhiIGOsrr54q+27XURURJICEduxsKTQMuXAbap6IjACuFFETgpwTE25Cah/+qrg8yiwWFX7AoMI0rhFJAWYCgxV1f54hooPtmHg/wacXWvbNOAdVe0DvON9Hwz+Rt1Y3wb6q+pA4Evg160dVCP+Rt14EZGewFjgm9YOqCUsKbSAqu5S1VXe9cN4PrT8m1k+AEQkFRgHPB3oWJoiIh2AU/HMuYGqHlXVg4GNqlERQKx3BsF21J1lMKBUdRl1ZzWcAMz3rs8HzmvVoBpQX6yqukRVq+YJ/RjPTI5BoYHfLcAjwK+oZ4rhYGZJwSEi0gsYDCwPbCSNmoXnj7Qy0IH44TggD5jnvd31tIg0MgFv4Kjqt8BDeL4R7gIOqeqSwEbll2RV3QWeLzhA1wDH469rgH8HOojGiMh44FtVbd6E3UHAkoIDRCQeWAjcrKoFgY6nPiLyY2Cvqq4MdCx+igCGAI+r6mCgiOC5vVGD9178BKA30AOIE5HLAhvV95OI/AbPbdsXAh1LQ0SkHfAb4O6mjg1GlhRaSEQi8SSEF1T1lUDH04hRwHgR2QYsAMaIyPOBDalRuUCuqlbVvF7GkySC0ZnA16qap6plwCvADwIckz/2iEh3AO/r3gDH0ygRuRL4MTApyOdyz8DzBWGN999bKrBKRLoFNCo/WVJoARERPPe8N6nqw4GOpzGq+mtVTVXVXngegr6rqkH7bVZVdwM7ROQE76YzgI0BDKkx3wAjRKSd92/iDIL0oXgti4ArvetXAq8FMJZGicjZwB3AeFUtDnQ8jVHVdaraVVV7ef+95QJDvH/TQc+SQsuMAi7H8617tXc5N9BBfY/8L/CCiKwFsoD7AxxPvby1mZeBVcA6PP+ugmqYAxF5EfgvcIKI5IrIz4AZwFgR2YynlcyMQMZYpYFYHwPaA297/509EdAgq2kg3pBlw1wYY4zxsZqCMcYYH0sKxhhjfCwpGGOM8bGkYIwxxseSgjHGGB9LCsZ4iUhFtabFq0XEsR7UItKrvlE0jQk2EYEOwJggUqKqWYEOwphAspqCMU0QkW0i8qCIfOJdMr3b00XkHe8Y/++ISJp3e7J3zP813qVqyItwEXnKO+/CEhGJ9R4/VUQ2estZEKAf0xjAkoIx1cXWun10SbV9Bao6HE/P2lnebY8Bz3rH+H8BmO3dPht4T1UH4RmvaYN3ex9gjqr2Aw4CF3q3TwMGe8uZ4tYPZ4w/rEezMV4iUqiq8fVs3waMUdWt3gEQd6tqZxHZB3RX1TLv9l2qmiQieUCqqh6pVkYv4G3vhDaIyB1ApKreJyKLgULgVeBVVS10+Uc1pkFWUzDGP9rAekPH1OdItfUKvnumNw6YA2QDK70T9RgTEJYUjPHPJdVe/+td/4jvpt2cBHzgXX8HuB58c2J3aKhQEQkDeqrqUjwTICUCdWorxrQW+0ZizHdiRWR1tfeLVbWqWWq0iCzH80VqonfbVGCuiPwSzyxxV3u33wQ86R0tswJPgtjVwDXDgedFJAEQ4JEgn3bUfM/ZMwVjmuB9pjBUVfcFOhZj3Ga3j4wxxvhYTcEYY4yP1RSMMcb4WFIwxhjjY0nBGGOMjyUFY4wxPpYUjDHG+Pw/G74R8ED1kjIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Score (RMSE): 0.13183554133416536\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4FVX6wPHvm04gAUJCDSX0XkMTlSIKKIIgCCoKIopd13UF265lf6uuuy52QAVRmtIFKSqgSBGk904gIUACISG93fP7Yy4a4g2kT8r7eZ775M7MmblvJrn3vXPOmXPEGINSSimVW252B6CUUqp00cShlFIqTzRxKKWUyhNNHEoppfJEE4dSSqk80cShlFIqTzRxKKWUyhNNHKrUEZEwEUkWkQQROSsiX4hIpSzbvxARIyKDsu03ybl+jHPZS0T+KyIRzmOdEJH/5fA6lx8fFkL8f3HGHSci00TEO4dyXiIy3xmHEZFe2bY/IyLHReSSiESKyP9ExCPL9jdEZI+IZIjIqy6O/6Tzd74kIltF5PocYjgoIhEF/b1V2aGJQ5VWtxtjKgHtgQ7AC9m2HwZGX15wfqAOB45lKfMCEAp0AfyA3sAOV6+T5fFEQYIWkX7AROAmoAHQEHjtKrusB0YBZ11sWwp0NMb4A62BdsBTWbYfBZ4HvnMRR1fgLWAYUBn4HFgkIu7Ziv4NiLrW76XKF00cqlQzxpwFVmElkKyWAj1EpKpzuT+wmys/gDsDi4wxkcYSZoz5sohDHg18bozZZ4y5CLwBjHFV0BiTZoyZZIxZD2S62H7MGBPrXBTAATTOsn2GMWYFEO/i8A2AfcaYbcYaPuJLIBCofrmAiIRgJa038/xbqjJNE4cq1UQkGBiA9e06qxTgW2Ckc/l+rA/HrH4FnhWRx0SkjYhIAeK4R0Rir/Ko5yzaCtiVZdddQA0RqVaA170EnMe64piSy11XAO4i0tV5lTEW2MmVifUD4EUgOT+xqbJLE4cqrRaLSDwQjlWV8g8XZb4E7heRykBPYHG27W8CbwP3AluB0yIyOluZxdkSwEOugjHGzDbGVLnK45SzaCUgLsuul5/75fL3dvW6/kBTYDJwLpe7xgMLsKrCUrHO38POqw9EZAjgYYxZlJ+4VNmmiUOVVncYY/yAXkBzrGqWKzireIKAl4FlxpjkbNszjTEfGWN6AFWA/wOmiUiLbK+TNQF8WsC4EwD/LMuXn7uqTso1Y8wRYB/wcS53GYd1ldEK8MKqklomIrVFpCLwb+DJgsSkyi5NHKpUM8b8DHwB/CeHIjOBv/Lnaqrsx0k2xnwEXARa5jUOEbk3W++r7I/LVVX7sKqULmsHnDPGXMjra7rgATTKZdl2wFJjzGFjjMMYsxI4A1wHNMFqA/lFRM4CC4Fazp5gDQohTlXKaeJQZcEk4GYRyd5ADvA+cDOwLvsGZ3fWXiJSQUQ8nNVUfvy5Z9U1GWNmZet9lf1xuarqS+BBEWnpbLh/GSvxuSQi3iLi41z0EhGfy20xIjJORKo7n7fE6iW2Osu+ns593QAP576Xe039BtwmIg3FcjNWddde56MuVoeD9lhXJ+ecz8Pzem5U2aOJQ5V6xphorA/kV1xsizHGrL5cd59NMvBfrAbh88DjwJ3GmONZyizNduVQoDp/5zf7fwNrgZPOx+/tMyKyT0TuzbLLIWecdbB6jyUD9Z3begB7RCQRWO58vJhl30+d5e8GXnI+v8+57UtgLvATcAkrwY43xhw0xmQYY85efgAxgMO5/KfeXar8EZ3ISSmlVF7oFYdSSqk80cShlFIqTzRxKKWUyhNNHEoppfLE49pFSp/AwEDToEEDu8NQSqlSY9u2beeNMUG5KVsmE0eDBg3YunWr3WEopVSpISInc1tWq6qUUkrliSYOpZRSeaKJQymlVJ6UyTYOV9LT04mIiCAlJcXuUMoMHx8fgoOD8fT0tDsUpVQxKjeJIyIiAj8/Pxo0aEAB5utRTsYYLly4QEREBCEhIXaHo5QqRuWmqiolJYVq1app0igkIkK1atX0Ck6pcqjcJA5Ak0Yh0/OpVPlUbqqqlFKq2GSkwtHV0Pgm8PCGqINwYCm4e0KlGlD/OqgcDInR4OYBvtXgWl/EkmLg+FqIPwcY6Hg/eOdrxuEC08RRhoWFhbFx40buueeePO03ZswYBg4cyLBhw4ooMqXKuJUvwNbPoVpjaHYrbJ4MmWnZCgngnNbCyw9qtrESStsRENT0j2IxJ2D1a3DwuyuPsWUqDP0U6naxli+ehMjt0GpIUf5mgCaOMi0sLIzZs2fnOXEopQrg4HIrabQcDGd2wcb3ocUguO2/4FURLoZB2AZIjLKuPhwZEHMcIrbC+v/BhknQ9REI7gzRh6x1bh4Q+iC0GQbVGkHUAVg0Hj6/BTrcCzVaw5p/gocPNLnFep0ipImjGM2cOZP333+ftLQ0unbtyosvvkjfvn3ZtGkTAQEB9OzZk1deeYWmTZvSv39/unbtyo4dO2jatClffvklvr6+bNu2jWeffZaEhAQCAwP54osvqFWrFkePHuWRRx4hOjoad3d35s2bx8SJEzlw4ADt27dn9OjRPPXUU0ycOJGffvqJ1NRUHn/8ccaPH48xhieffJI1a9YQEhKCTu6l1FWcP2J9yMeFW9VQvtWsdWd2giMTzu62rh6GfgrGYX3I1+7wR1VUjVbWw5WEaFjzOmz6iN+vRpr0g4H/g8p1/ihX/zp4ZAP8/DZsngKOmdCwNwz6oMiTBpTRGQBDQ0NN9rGqDhw4QIsWLQB4bek+9kdeKtTXbFnbn3/cnsM/g/P1n3/+eRYuXIinpyePPfYY3bp1Iy0tjZUrV9K1a1eOHj3KlClTCAsLIyQkhPXr19OjRw/Gjh1Ly5Ytefrpp+nZsydLliwhKCiIr7/+mlWrVjFt2jS6du3KxIkTGTJkCCkpKTgcDrZs2cJ//vMfli1bBsDUqVOJiori5ZdfJjU1lR49ejBv3jx27NjBJ598wsqVKzl37hwtW7bks88+y1VVVdbzqlSZlp4MP70FGz+A7DPounlCzdbg6Wt96x/wbwhsnP/XijkBaYngGwB+ta7e/hFz3Loyadr/2u0kVyEi24wxobkpq1ccxWT16tVs27aNzp07A5CcnEz16tV59dVXmTdvHpMnT2bnzp2/l69bty49evQAYNSoUbz//vv079+fvXv3cvPNNwOQmZlJrVq1iI+P5/Tp0wwZYtVt+vj4uIzh+++/Z/fu3cyfPx+AuLg4jhw5wrp167j77rtxd3endu3a9OnTp8jOg1KlQlqSVc10/hBcOgNn98DJ9ZASBx3ugx5PQ5V6VptDQhT41wFP1++7fAnIw71RAQ2tRzEql4njalcGRcUYw+jRo3nzzTevWJ+UlERERAQACQkJ+PlZvSSyd3UVEYwxtGrVik2bNl2x7dKl3F09GWP44IMP6Nev3xXrly9frl1rlYo+DJs+hNPbrOqlrFcVVUOsNos2d0HIDX+s9/C2rWeTncrVfRx2uummm5g/fz5RUVEAxMTEcPLkSSZMmMC9997L66+/zkMPPfR7+VOnTv2eIObMmcP1119Ps2bNiI6O/n19eno6+/btw9/fn+DgYBYvXgxAamoqSUlJ+Pn5ER8f//sx+/XrxyeffEJ6ejoAhw8fJjExkRtvvJG5c+eSmZnJmTNnWLt2bbGcE6VKjGNr4LO+sHeh1WB9w7Nw91x4Zg+8HA1P77TaD7ImjXKsXF5x2KFly5b885//5JZbbsHhcODp6cm7777Lb7/9xoYNG3B3d2fBggVMnz6d3r1706JFC2bMmMH48eNp0qQJjz76KF5eXsyfP5+nnnqKuLg4MjIyeOaZZ2jVqhVfffUV48eP5+9//zuenp7MmzePtm3b4uHhQbt27RgzZgxPP/00YWFhdOzYEWMMQUFBLF68mCFDhrBmzRratGlD06ZN6dmzp92nS6mic3obrHoZMOBZwbovIvogBDWHe76GKnXtjrDEK5eN4yVdWFgYAwcOZO/evXaHck2l6byqcswYSL4Ih1fBsmegQoDVrTU9ybrCCGwKN/wVfPztjtQ22jiulFKXXYqEGbfDhaPWct1uMHIWVAy0N65STBNHCdSgQYNScbWhVImXngJz74X4s3DzGxDYBBrdBB5edkdWqmniUEqVfunJEPEbJMdaVVIpsdbziN+sYThGzIIWA+2OsszQxKGUKt2MgVnDIeyXK9e7eUCFqtDvX5o0CpkmDqVU6XbkBytp9JwALW4HnypWwvCqWKA7qXMSn5LOyQtJnE9IJT3T4DCG+JQMLiamEZOURmxSOqnpmaRmOkjLcJCR6Sj0GHLiX8GT90Z2KPLX0cShlCq9HA5Y/TpUbQA3PFckbRfGGGIS0zgSlcCSnadZuP00qRmuk4GHm1DF1xMfT3e8PNzwcnfD092tKPKXS5nF1ElWE0cpValSJRISEoiMjOSpp576fRgRVyZNmsTDDz+Mr68vALfeeiuzZ8+mSpUqxRWuUkVj/yI4t8caULAAScMYw9aTF1mx5ywRF5M4F59KfHI6l1LSuZSSQZozUXh7uDG0Yx16Ng0iyM8bbw93APx9PKlS0RM/b49yMQqDJo4SJDMzE3d39zztU7t27asmDbASx6hRo35PHMuXL893jEqVGMbAL+9aN+61zv3cMUej4lmx5ywpGZlkOiAhNZ19kZfYcSoWH0836lb1pWZlH+pWrYB/BU/8fTwJ8vOmUVBF2tetQhVf7ZGliaOYhIWFuRwqvWXLlowdO5bvv/+eJ554gs6dO/P4448THR2Nr68vn376Kc2bN+fEiRPcc889ZGRk0L9//yuOe/lmwczMTCZMmMCqVasQER566CGMMURGRtK7d28CAwNZu3YtDRo0YOvWrQQGBvLuu+8ybdo0AMaNG8czzzxDWFgYAwYM4Prrr2fjxo3UqVOHJUuWUKFCBbtOn1J/dnwtnNsLgz8Gt5xHTzodm8y09ScIj0ni7KUUdkfEAVa1kpsIfj4eBPl589qgVtwVWpcKXnn78lYe2ZI4RGQ48CrQAuhijNmaQ7n+wHuAO/CZMeatQglgxURrtMvCVLMNDLh6eIcOHeLzzz//faj0jz/+GLBGs12/fj1gjWk1efJkmjRpwubNm3nsscdYs2YNTz/9NI8++ij3338/H330kcvjT506lRMnTrBjxw48PDyIiYkhICCAd999l7Vr1xIYeOUNT9u2bWP69Ols3rwZYwxdu3alZ8+eVK1alSNHjjBnzhw+/fRT7rrrLhYsWMCoUaMK4UQpVUg2fgCVakKbYeyJiGPDsfOkpjtIzcgkNcP6GZuUzvf7zgHQMKgi/hU8mdC/OcNDgwms5G3zL1B62XXFsRcYCkzJqYCIuAMfATcDEcBvIvKtMWZ/8YRY+FwNlQ4wYsQIwBodd+PGjQwfPvz3fVJTUwHYsGEDCxYsAOC+++5jwoQJfzr+jz/+yCOPPIKHh/VnDQgIuGo869evZ8iQIVSsaE38MnToUH755RcGDRpESEgI7du3B6BTp06EhYXl99dWqvA4Mq0vfdEH4dgaHH3+weT14fz3+8NkOqyWYU93wdvDHW8PN7w93BjSoQ5P9W1CnSp6xVxYbEkcxpgD8Oehw7PpAhw1xhx3lp0LDAYKnjiucWVQVFwNlQ78/sHtcDioUqXKFfNyXG3/7IwxeWqYu9o4Zd7ef3wbc3d3Jzk5OdfHVapQpCVZw4TERVjDlyecs9o0LhyxtvsG8sbZrkzffojb2tbijcGtqVLBEze3st84bbeSPKx6HSA8y3KEc51LIvKwiGwVka3R0dFFHlx+uBoqPSt/f39CQkKYN28eYH2w79q1C4AePXowd+5cAGbNmuXy+LfccguTJ08mIyMDsIZuB/40vPplN954I4sXLyYpKYnExEQWLVrEDTfosNHKRg6HNW3q5/3gzWCYcgPMvRtmDoXFj4K7F9wxGcatZnmvpUzfHsv4ng358O4OBFT00qRRTIoscYjIjyKy18VjcG4P4WJdjl+RjTFTjTGhxpjQoKCg/AVdxC4Pld62bVtiYmJ49NFH/1Rm1qxZfP7557Rr145WrVqxZMkSAN577z0++ugjOnfuTFxcnMvjjxs3jnr16tG2bVvatWvH7NmzAXj44YcZMGAAvXv3vqJ8x44dGTNmDF26dKFr166MGzeODh2K/uYhpXK07h1Y9SJkJMP1z8DwL+ChNTB2FYxdxbl7f2SFey9mRgTxt2Wn6NygKn+7pVm56AJbktg6rLqI/AQ856pxXES6A68aY/o5l18AMMa8mb1sdiVxWPXSNFR6Xth9XlUZcmgFzBlJbJM7+TbkFQ5HJXAxMZ341AyMMVxMSmPv6T9mu6zu583ix3tQW9suCkVZGVb9N6CJiIQAp4GRwD32hqSUKjSOTLhwzBqEcN8izJEfiK7Ughv23E7qnv34O7vJVvL2wM1N8PP25G/9mnFDk0ACK3lTrZLX7zfgqeJlV3fcIcAHQBDwnYjsNMb0E5HaWN1ubzXGZIjIE8AqrO6404wx++yItzDoUOlKOWWmw87ZmHXvIHFWM2ayT3V+qTKcF870pH/7Bkzo35xalX20CqqEsqtX1SJgkYv1kcCtWZaXA4V2m3Neex2pqyuLs0eqIhJ/lrRNU8g4/APeF4/gnpnCfmnCF+kPs9cRwqGUunglejD+pkY807eJvk9LuJJcVVWofHx8uHDhAtWqVdN/ykJgjOHChQv4+PjYHYoqoYwxHN+/lcwNHxASuRwPk8FWRwv2mT5sNK1xNOzLDU2DuKmqL42rVyIksCLu2iuqVCg3iSM4OJiIiAhKalfd0sjHx4fg4GC7w1AlUHxKOsunvsKImE9INl4s9byZM83H0Kh5O9r6ejI4qCLV/fRLR2lVbhKHp6cnISEhdoehVJkUHZ/Kvsg4jkYlEJecTtjONfwvcQrHA3tSYdgnDK2V4y1YqhQqN4lDKVW4Mh2GhdsjmLc1gi1hMXiQQUs5SbDbBf7hNZMMv2AaPvQV+FS2O1RVyDRxKKXy7PC5eP42fze7wmO5sVosy+uvoMnFX/BMi7UKeFSEu+dp0iijNHEopfIk7Hwigz/cgJ9HJqtbrqBh2Bwk3QeaD4Sm/aBaY6haX5NGGaaJQymVa8YYXlq8B1+3dNbVn47P8dXQ6QHo/SJUqm53eKqYaOJQSuXagu2n2XQ0mvXBH+MTtglufx86jbY7LFXMNHEopa7KGMOszadYuD2CneGxPFFjL7XPb4Rb/6NJo5wqycOqK6Vs5nAYXlu6n5cX7yUl3cHjvRrylMdia57v0AftDk/ZRK84lFIuGWOYsGA387ZF8OD1Ibx0awvcDiyGTQfhzs+vOs+3Kts0cSilXHp75SHmbYvgqT6N+cv11ZHt0+GX/0FgU2g1xO7wlI00cSilrpCSnsmHa44y+edj3NetPn/pUgH5uBsknIVqTWDgJHDT4czLM00cSinAqppasfcs//fdAU7HJjO0Qx1e7VcPmd4f0pPggZVQrxvoIKHlniYOpco5YwwHz8Yz6cfDrNp3jha1/HlnWFuua1QNvh4F0Qfh3nlQv7vdoaoSQhOHUuXY0l2R/PO7/Zy7lIq3hxsTBzRn3PUheLi7wbYZcHAZ9PsXNL7J7lBVCaKJQ6lyatbmk7y8eC/tgqvw11ua0atZ0B9DncecgFUvQsiN0PVRewNVJY4mDqXKkfMJqaw7HM2C7RFsOHqBPs2r8/G9HfHxdDZ2Oxyw+2tY8waIGwz+WLvdqj/RxKFUGZaSnskHa45w4Ew8YRcSOR6dCEC9AF+eu6Up43s2wtPdmRiSY2HeaDj+E9TuAHd9CVXq2he8KrE0cShVRiWnZfLQl1vZcOw8LWr60zioEsM71aV7o2q0C6585RTKF47B3HvgwlEY+D/oOEavNFSONHEoVcZkZDr45ch5Plx7lO2nLvLOsHYM65TDFL+pCbDu3/DrJ+BRAUYthIY9izdgVepo4lCqlDLGEJ2QSsTFZBJSMsg0ho1Hz7N4ZyTR8alU9fVk0oj2DG6fw7StMSesq4yo/dD+XujzCvjXKt5fQpVKmjiUKmUyHYZluyP5cM1RjkQlXLHNw03o07w6d3YKpnez6nh5uKhuMgb2LoDlz1nP71sEjfoUU/SqLNDEoVQpseVEDFN+PsaWEzHEp2bQtEYlXhnYkpBAX/x9PBGBkMBKBFT0+vPOjky4GAYnN8CuudbPWu1h+HQIaFjsv4sq3TRxKFXCGWP46teTvL50P4GVvBnYrja9mwXRt0UN3NyuMfxHzHGYNwbO7QdHurXOr7Y13lTH+3XMKZUvmjiUKsEcDsPry/bzxcYwbmpenf+NbI+/j2fudo4NhxmDIS0BrnsCqoZYY00FNtXxplSBaOJQqoTKyHTwt/m7WbTjNGN7hPDybS2ufYUBkHgedn8Dmz6C1EsweinUbl/0Aatyw5bEISLDgVeBFkAXY8zWHMqFAfFAJpBhjAktrhiVstvHPx1j0Y7TPHdLUx7v3fjK+y4AMjOsdouLYRC1DyJ3QOROuHjC2l4nFO6aoUlDFTq7rjj2AkOBKbko29sYc76I41GqRDlxPpEP1x5lYNtaPNGnyZUbwzbAjplweAUkX/xjfZV6VoN3x/ugyS1Qs03xBq3KDVsShzHmAPDnb1BKKYwxvLx4D97ubvx9YMs/NoRvgR9fg5PrwacyNB1g3awX0NBqt/ANsC9oVa6U9DYOA3wvIgaYYoyZmlNBEXkYeBigXr16xRSeUoVv5uZTbDh6gTcGt6K6vw8kxcDKidbgg5VqQv+3odNo8Kxgd6iqnCqyxCEiPwI1XWx6yRizJJeH6WGMiRSR6sAPInLQGLPOVUFnUpkKEBoaavIVtFI22xUeyxtL99OrWRD3dq0PKXHw1R0QdQBu+Ctc/yx4V7I7TFXOFVniMMb0LYRjRDp/RonIIqAL4DJxKFWaJaZmsHLvWf7z/SGC/LyZNKI9bhlJMOsu6x6Mu+dCkwK/pZQqFCW2qkpEKgJuxph45/NbgNdtDkupQheTmEa/SeuIjk+lfjVfPrqnI1V8vWDJXyB8s9UzSpOGKkHs6o47BPgACAK+E5Gdxph+IlIb+MwYcytQA1jkbED3AGYbY1baEa9SRWnmryeJjk9l+gOd6dU0yOo0sneh1XPqhueg5WC7Q1TqCnb1qloELHKxPhK41fn8ONCumENTqlilpGfy5aYwejcLonez6tbK80dg6TMQ3Bl6TbQ1PqVc0ZlalLLRkp2nOZ+QxkM3OAcajDoI028FDy8Y+im453J4EaWKUYlt41CqrEvLcPDpLydoUcuf7u4H4LslsHc+uHvB6GUQEGJ3iEq5pIlDKZv8a/kBjkYlMG+AQWbcC56+ENITbvknBDa2OzylcqSJQ6lilpKeyZKdp/liYxgPXlePzgceBf9gePxX8PazOzylrkkTh1LFxBjDX7/ZxeKdp3EY6FS/Ki/W3ALb98CwaZo0VKmhiUOpYrJ8z1kW7jjN8E7B9GlenV5Vz+M+83Wo3wNaDbU7PKVyTROHUsUgMTWDN5btp0eNDN4Kjcdd4mG2c7ypwR/qxEqqVNHEoVQxeG/1EdwvhTODv+M+I85aWamGNcmSzvmtShlNHEoVsbWHovhi3SFWV52Mh8PAXV9BejI0uB4q17E7PKXyTBOHUkUo4mISf5m7g/f8Z1E3+QCMmAktbrc7LKUKRBOHUkUgMTWDb7aG896Ph/mr4wsGpH1vDYuuSUOVAZo4lCoEcUnpzNpykn2nL3HgzCVOXEjEy6TxSbVv6JO4HLo+Cn1esTtMpQqFJg6lCigpLYPHPl9N87NL6VchirHeybg3qE3TpG34xh2FHk9D39e055QqMzRxKFUAqRmZvDpjBa9H/41GnmfAqxpUqAYxu615we9bBI362B2mUoVKE4dS+bRy50m2L/+U51K/orJXJoz6zuopdZkxepWhyiRNHErlw9IfVtNl/Vj6SywJVZvjfc8XUL3FlYU0aagyShOHUnm05/hpWq5/Ah93IfPuRVRq3FuThCpXNHEolRsOB4nH1rNz/yEydn1DSzlL0vCFuDfpbXdkShU7TRxK5ST5IpzZRVrkXhI2fkZA0gl6ODed7TyBmi00aajySROHUi6Yc/tI+exWKqTH4gWccdRncY2JdO/Ri+b1a1Ozan27Q1TKNpo4lMrKkUnG6Z0kzxhGUjq86f86Feq2pU+nNoxtFGh3dEqVCJo4lAKIOUHMor/hH/ETHiadFOPPio6f89qgmxFt+FbqCrlKHCLytDHmvWutU6rUSIqB7TMgbAOOjFQyT27GyyHM4WYSqjQluNNAxtzY2e4olSqRcnvFMRrIniTGuFinVMkVGw4bP4BzezGntyMZyUT5NuZMsicH0rtzos1TPDOkFxW83O2OVKkS7aqJQ0TuBu4BQkTk2yyb/IALRRmYUoUiPRncvSBqP8wchiMllhMejViXciNzMvtwJLUu1zcO5P7uDRjZsobd0SpVKlzrimMjcAYIBP6bZX08sLuoglKqUBxbC3PuBkcGRoREj6rclfwGJ9zqcf919ZnYMIA2daoQ5Odtd6RKlSpXTRzGmJPASaB7Yb6oiLwD3A6kAceAB4wxsS7K9ceqDnMHPjPGvFWYcagyLCEKFo0n3a8OYYG9CDsTxSvnb6ZD61ZMu70VNSv72B2hUqVWbhvH4wHjXPQCPIFEY4x/Pl/3B+AFY0yGiLwNvABMyPaa7sBHwM1ABPCbiHxrjNmfz9dU5YUjExY9QmZyLIMuPsuBM3Wp4OnOi3e0YFTXetpLSqkCylXiMMb4ZV0WkTuALvl9UWPM91kWfwWGuSjWBThqjDnufM25wGBAE4fKmcMBS5+GY6v5wOdREr2bsXxUJxpXr4SXh5vd0SlVJuTrPg5jzGIRmVhIMYwFvnaxvg4QnmU5Auia00FE5GHgYYB69eoVUmiq1AjbAJHbIXIn7J3PtvrjmHToBqaNaUnL2vm9MFZKuZLbqqqhWRbdgFD+qLrKaZ8fgZouNr1kjFniLPMSkAHMcnUIF+tyfE1jzFRgKkBoaOhVY1NlTPRhmHE7mEwAzrUZz327etO3RSB9mmtfpGqAAAAacklEQVRPKaUKW26vOG7P8jwDCMOqNsqRMabv1baLyGhgIHCTMcbVB30EUDfLcjAQmZtgVTmz+jXw9IXHNvLLGTfGz9lLkJ83b9zR2u7IlCqTctvG8UBhvqizt9QEoKcxJimHYr8BTUQkBDgNjMS6p0SpP5z6FQ4u43Crp5kwO5wdp2JpXtOPLx/sQnU/7TmlVFHIVWuhiDQUkaUiEi0iUSKyREQaFuB1P8S6ifAHEdkpIpOdr1NbRJYDGGMygCeAVcAB4BtjzL4CvKYqa2JPwbK/kFmxBkO2tyMuKZ0XBjTnm0e6a9JQqgjltqpqNlbX2CHO5ZHAHK7SWH01xpjGOayPBG7NsrwcWJ6f11Bl3NHVMH8sGAdz671G6sUKzBzXldpVKtgdmVJlXm77J4ox5itjTIbzMZNrNI4rVWQuRcK8MeBfh/gxq3nzcB1ua1tLk4ZSxSS3Vxxrnd1v52IljBHAdyISAGCMiSmi+JS6kjGw7FnITIeRM5m7FxJSM3johoLUnCql8iK3iWOE8+f4bOvHYiUSfdeqouNwwIZJcGwN+FSGwytI6f0akzanMW3DCbo3rEbrOpXtjlKpciO3iaOFMSYl6woR8cm+TqlCc+EY/PoJVK0PJ9bBke+hekuI2o+pfz137+nIztPHGNKhDhMHNLc7WqXKldwmjo1Ax1ysU6rgjIElT8CpTYABN0+47b8Q+iCIsHRXJDvm7ODfd7blrs51r3k4pVThutZ8HDWxhv6oICId+ONubn/At4hjU+XVnvlwaiPc/h60GATGARWt+b7TMhz8Z9Uhmtf0485OwTYHqlT5dK0rjn5YM/0FA+9mWR8PvFhEManyLCUOfngFarWHDveBmzUbX3xKOpuOXeCH/ec4FZPE9Ac64+6mo9wqZYdrzccxA5ghIncaYxYUU0yqvIk5Dukp1s8VEyDhHNz1Fbi5k5iawYLtEUz68QgxiWmIwOD2tenVNMjuqJUqt3LbxtFaRFplX2mMeb2Q41HlzZndMOVGfr8tKLAZjF3FDtOYv/73J45HJwLQrWEAT93UhPZ1q+Drla9BnZVShSS378CELM99sAYnPFD44ahyZ9t08PCBwR9ac4M37Ydx9+L1TzaSkJLBszc3JbR+Vbo3qqYTMClVQuR2kMOs840jIv8Bvi2SiFT5kZoAu+dBqzugzR9zeW04cp4dp2L5vyGtubdrfRsDVEq5kt9rfl/0pj+VX5dH0d+3CNLioePoKza/v+YItSr7MEx7TSlVIuV2Iqc9/DE2lRtQHXijqIJSZci5fbB5Mog7tL7Tujdj04fgVwsyUqw2jXrdMMbwW9hFvtkazpYTMbx6e0u8Pdztjl4p5UJurzgGAlWBG4AqwHJjzLYii0qVfonnYflz1lWFZ0XAWO0ZAE37Q3IshP+KY8A7TF13nNmbT3EqJglfL3fu6VqPkV10+l+lSqrcJo7BwEPAQqybAKeLyKfGmA+KLDJVeoVtsEavTYmFnhOg6yPg7mkNG1I1BOo4BxxIimH1iTTeWrSNriEBPNO3Cf1b19ReU0qVcLl9h44DuhljEgFE5G1gE6CJQ10pMx0WPwreleC+RVAzy/Stre+8sqxvAN9s20pgJW9mjuuKp3tuR/lXStkp1/NxAJlZljP5Y/gRpf6wZx7EnoR+/7oyabgQFZ/CmoNR3NmpjiYNpUqR3F5xTAc2i8gi5/IdwOdFE5IqtRyZsO4/ULON1Y5xDQu3nybTYbgrVAcqVKo0ye19HO+KyE/A9VhXGg8YY3YUZWCqFNq7EGKOwYiZkOVmvYTUDKIupeAwEH4xiUNn40lKy2TBtghC61elUVAlG4NWSuVVrlshjTHbge1FGIsqzRyZsO4dqN4Kmt0GQEp6Jl9sDOOjNUeJT8340y4ebsKLt7Yo7kiVUgWk3VdU4di/BM4fguFfcCY+lSk/H2fJztNcTErnpubVGdiuFm4i1PT3oXktfypX8LQ7YqVUPmniUAXncFhXG4HNoMVgnp/+G5uPx3BzqxqM6lqf7o2q2R2hUqoQaeJQBbdnHkTth6GfcSEpnQ1Hz/Nor0b8rZ9O6apUWaR9IFXBhG2Ab5+E4M7Qeigr953FYeC2NrXtjkwpVUQ0caj8O7sH5oyEqg3gnm/AzZ3vdp+hYVBFWtTyszs6pVQR0cSh8ifmBMy8E7z94L6F4BtAdHwqvx6/wMA2tXTuDKXKME0cKu9Ob4evhkBmGoxaCJWt4c9X7D1jVVO11WoqpcoyWxrHReQd4HYgDTiGdUNhrItyYUA81hAnGcaY0OKMUzkZA4eWw+ltcGoznFwPFarCPfOgutUAnpqRyZSfj9OmTmWa1tAb+pQqy+zqVfUD8IIxJsM5YOILwIQcyvY2xpwvvtDUn/z6Cax6wZpTo1pjuPkNCH3AqqZymvnrKU7HJvPWnW20mkqpMs6WxGGM+T7L4q/AsJzKKpvEhkN6EsRFwPcvQ/OBMGw6eHj9qeillHQ+XHOE6xsHckOTIBuCVUoVp5JwH8dY4OscthngexExwBRjzNScDiIiDwMPA9Srp5MAFcjeBbBwPDjSreXAZjBkssukAfDvlQe5mJTOhP5634ZS5UGRJQ4R+RGo6WLTS8aYJc4yLwEZwKwcDtPDGBMpItWBH0TkoDFmnauCzqQyFSA0NNS4KqNyYdsMWPo01OsOnUZD/BloM/yKaqmsluw8zcxfT/HQDSG0Ca5czMEqpexQZInDGNP3attFZDTWlLQ3GWNcftAbYyKdP6OcQ7p3AVwmDlUIzu6B756FxjfBXV+Bl+9Vix8+F8/EBXvo3KAqz+vVhlLlhi3dcUWkP1Zj+CBjTFIOZSqKiN/l58AtwN7ii7Kcycyw7gCvUBWGfnrNpJGQmsEjM7dR0duDD+/pqBMxKVWO2NXG8SHgjVX9BPCrMeYREakNfGaMuRWoASxybvcAZhtjVtoUb9mVlgQnN8DB7yByh9UA7htw1V2MMUxYsJuw84nMGteNGv4+xRSsUqoksKtXVeMc1kcCtzqfHwfaFWdc5dI398HRH63nbUdAqyE5FjXGsC/yEv9afoCNxy7wfP9mOvKtUuVQSehVpexy4ZiVNK57Enr8BSq6TgKXUtKZMH83W07EcCExjaq+nrw+uBWjutYv5oCVUiWBJo7ybPsM66a+7k/kmDQA/rlsP6v2neXOjsG0rlOZOzrU0YmYlCrHNHGUVxlpsGMWNBsAfq56TVvWHozim60RPNarkfacUkoBmjjKh6SYKxu8HQ7YOROSzkOnMS53OR2bzLyt4czYGEazGn483bdJ8cSqlCrxNHGUdT+9DT/9yxoypPVQOPGLNWBhwjkIbAqN+vxpl90RsYyc+ivJ6Zlc16garw1qhbeHuw3BK6VKIk0cZdnB76ykUberlTAOLgPPitCkr5VImvYHtysTwonziTww/TcCKnoxe1w36lW7+v0cSqnyRxNHWXV6Oyx8GGp3hPu/hYxkiDpgLXu6vu8iPdPBQ19uxQBfju2iSUMp5ZImjrLoxDqYczdUCIARM61E4ekD9a+76m5f/xbO0agEPr0/lIZBOqeGUso1HSeirAnbADOHWbPyPbgKKtfJ1W6JqRlM+vEIXRoE0LdF9SIOUilVmukVR2lmjFX9dPwnq6HbvzbMvRuq1ocxy696b8ZlaRkOdobHMm9rOOcTUplyXyediEkpdVWaOEqrsA2w4nk4l23cx4rV4d55uUoaAE/N2cHKfWcBGNYpmE71qxZ2pEqpMkYTR2n0w99hw3tQuR7c9l9ocguc2QWHVkDX8VC1Qa4OE3Y+kZX7zjK6e30e692Y6n7eRRu3UqpM0MRR2kTutJJGu3uspHF5+PMq9aDF7Xk61KzNJ/FwEx7v3ZjqOsKtUiqXtHG8tFnzT2vOjAFvX3POjKtJTsvkm60R9GtVU5OGUipP9IqjNEiJg4QouHgSjv4AfV8DH//8Hy49kynrjhGXnM6objrCrVIqbzRxlGSnfoWf3oKwX8CRYa2rWB26PJSnw8QmpZGW4SD8YhILtp9m6a5I4lMy6NIggG4Nrz5pk1JKZaeJo6Q6t8+6H8PHH7o9BjVaWVcedULBq2KuDzN3yyleWLSHy7O6+3i6MaB1Le7sGEz3RtW0661SKs80cZQUiedhy1QIWw8128CBpeBdCR78Idc38WV36kISry3dT+cGAQxuXxs/H0/6NK9OJW/9syul8k8/QUqCoz/C3FHWeFI128C2GeDmAQ8sz3fScDgMz83fhYebMGlEe2pXqVDIQSulyitNHHaLDYcF4yCgIQyfDkHNIDMd0pPAp3K+D/v2yoNsORHDv4e11aShlCpUmjjsFH8W5j8AmRkw4iuo1sha7+4J7vlPGp+vP8GUdce5r1t9hncKLqRglVLKoonDDnERsGKCdae3ccDwL/5IGgWQkengfz8e5qO1x+jfqiavDmqljd9KqUKniaO4HV1tVU1lpsF1T0KHURCY/2lZ0zIcbDh6nh3hsaw9GMWe03GMCK3La4Nb4e6mSUMpVfg0cRSnY2th1nBrJNsRXxUoYaSkZ/Le6iN8/Vs4MYlpuAk0CqrEeyPbM7h9/hrUlVIqNzRxFJeY4zBvjJU0Hvy+QHd+H49O4PHZOzhw5hIDWtdkeGgw3RsGUsFL5wVXShU9TRzFISPV6m4LcPfsAiWN+JR07v70V9IyHEwbE0qf5jUKKUillMod2wY5FJE3RGS3iOwUke9FpHYO5UaLyBHnY3Rxx1koNrwPUftg6FSr220BvL/6CFHxqUwb01mThlLKFnaOjvuOMaatMaY9sAz4e/YCIhIA/APoCnQB/iEipWumoQvHYN070PIOaNqvQIc6fC6e6RvCGBFalw71StdpUEqVHbZVVRljLmVZrAgYF8X6AT8YY2IAROQHoD8wp+gjzCeHAw6vgE0fQ/wZ60Y+dy/o/2a+D2mMYc3BKN5Ytp+K3h483795IQaslFJ5Y2sbh4j8H3A/EAf0dlGkDhCeZTnCuc7VsR4GHgaoV69e4QaaWwlR8PUoCN9sTaxUJxSSY6DDfdZ84HkUn5LOoh2nmb8tgt0RcTQMrMjkUZ0IqOhVBMErpVTuFGniEJEfgZouNr1kjFlijHkJeElEXgCewKqWuuIQLvZ1dWWCMWYqMBUgNDTUZZlCd+EYbPoQzu2HOp3g4FJIiIZBH0K7u8E9/6d3x6mLPDlnBxEXk2lRy5837mjNyM518XTXubeUUvYq0sRhjOmby6Kzge/4c+KIAHplWQ4GfipwYIVh3Tuw5v+saqiabeC3z6yZ+R74zkoi+WCMYUd4LPO2RjBvazg1/H34Znx3uoTonBlKqZLDtqoqEWlijDniXBwEHHRRbBXwrywN4rcALxRHfFf109vw07+g9Z3Q703wq2F1uRX3fF9l7I6I5Z/fHWDLiRh8PN24s2MwL97agsq+noUcvFJKFYydbRxviUgzwAGcBB4BEJFQ4BFjzDhjTIyIvAH85tzn9csN5cUu8QJ8/xKc3AixJ62qqMEfgZvzpjsP73wfevLPx3hrxUGqVfTitUGtGNqxDn4+mjCUUiWTnb2q7sxh/VZgXJblacC04oorR/sXw6450Ow2a4yp0LF/JI18Msbw1sqDTPn5OLe3q82/hrTWhKGUKvH0zvHcCt9szfc9chYU0oiz//3+MFN+Ps6obvV4fVBr3HRQQqVUKaCJI7dO/Qr1uhZa0vhqUxgfrj3KyM51eWNwax3+XClVamjfztyIP2u1a9TtViiHm7PlFH//dh99W1Tnn3do0lBKlS56xZEb4Zutn3W7FvhQlxvCezUL4oO7O+Kh92UopUoZTRy5cWozePhArXYFOswP+8/x1oqD3N6uNv8d3g4vD00aSqnSRz+5rub0NkhLgvBfoXZH8Mj/UB8XE9N4cdEemtf006ShlCrV9IojJ5unwIrnwa8WJEZbXXAL4B/f7uNiYhpfPNBZk4ZSqlTTxJFdegocXgkrJkCjPpCWaI1y27BXvg/59W+n+HZXJH+9uSmtalcutFCVUsoOmjguMwbeqgepztHe63SCEbPAswLEnoKq9fN12L2n43hlyT5uaBLIY70bF2LASillD00cl4lAt8fA3RN8q0GrO8DL19qWj6SR6TAs2BbBv1cdJMDXi0kj2uOuN/gppcoATRxZ9S6c8RN/ORLN/313gINn4+lQrwpvDm1DtUr5H8tKKaVKEk0chSQhNYPle84wf2sEW8JiqBtQgQ/v6cBtbWrpDX5KqTJFE0cB7YmIY9qGE6zce5bk9ExCAivy8m0tuK97fbw9CjYIolJKlUSaOPLpdGwy76w8yOKdkfj5eDCkYx3u7BhMx3pV9ApDKVWmaeLIo+S0TD5Yc4TP158A4LFejXi0VyMdDl0pVW5o4siDuKR0xs74jW0nLzKkQx2e69eMOlUq2B2WUkoVK00cuRSXnM6IqZs4Hp3IJ/d2ZECbWnaHpJRSttDEkUsr957h4Nl4po0JpU/zGnaHo5RSttFBk3JpZ3gslSt40rtZdbtDUUopW2niyKUdp2JpV1d7TCmllCaOXEhMzeDwuXja161idyhKKWU7TRy5sDsiDoeBDpo4lFJKE0du7AyPBaCdJg6llNLEkRs7wy/SoJovARXzPwOgUkqVFZo4rsEYw45Tsdq+oZRSTnofRw7OxCXz3o9HOJ+QRlR8qiYOpZRy0sThwpqD5/jrN7tISXdQL8CXdsGVuamF3vSnlFJgU+IQkTeAwYADiALGGGMiXZTLBPY4F08ZYwYVZVxpGQ7eWXWQT385QYta/nx0TwcaBlUqypdUSqlSx64rjneMMa8AiMhTwN+BR1yUSzbGtC+OgOKS0rl/+hZ2hcdyX7f6vHRbC3w8dT4NpZTKzpbEYYy5lGWxImDsiCMrPx8PGlTzZfyNDblVBzBUSqkc2dbGISL/B9wPxAG9cyjmIyJbgQzgLWPM4qsc72HgYYB69erlOR43N+G9kR3yvJ9SSpU3YkzRfNkXkR+Bmi42vWSMWZKl3AuAjzHmHy6OUdsYEykiDYE1wE3GmGPXeu3Q0FCzdevWAkSvlFLli4hsM8aE5qZskV1xGGP65rLobOA74E+J43KDuTHmuIj8BHQArpk4lFJKFR1bbgAUkSZZFgcBB12UqSoi3s7ngUAPYH/xRKiUUiondrVxvCUizbC6457E2aNKREKBR4wx44AWwBQRcWAluLeMMZo4lFLKZnb1qrozh/VbgXHO5xuBNsUZl1JKqWvTsaqUUkrliSYOpZRSeaKJQymlVJ4U2X0cdhKRaKxG9/wIBM4XYjiFRePKu5Iam8aVdyU1trIUV31jTFBuCpbJxFEQIrI1tzfBFCeNK+9KamwaV96V1NjKa1xaVaWUUipPNHEopZTKE00cfzbV7gByoHHlXUmNTePKu5IaW7mMS9s4lFJK5YlecSillMoTTRxKKaXyRBOHk4j0F5FDInJURCbaHEtdEVkrIgdEZJ+IPO1cHyAiP4jIEefPqjbF5y4iO0RkmXM5REQ2O+P6WkS8bIipiojMF5GDzvPWvSScLxH5i/NvuFdE5oiIj13nS0SmiUiUiOzNss7lORLL+873w24R6VjMcb3j/FvuFpFFIlIly7YXnHEdEpF+xRlXlm3PiYhxjtxdrOfrarGJyJPO87JPRP6dZX3hnjNjTLl/AO5Y83w0BLyAXUBLG+OpBXR0PvcDDgMtgX8DE53rJwJv2xTfs1jzqCxzLn8DjHQ+nww8akNMM4BxzudeQBW7zxdQBzgBVMhynsbYdb6AG4GOwN4s61yeI+BWYAUgQDdgczHHdQvg4Xz+dpa4Wjrfn95AiPN9615ccTnX1wVWYd1kHFjc5+sq56w38CPg7VyuXlTnrMj/WUvDA+gOrMqy/ALwgt1xZYlnCXAzcAio5VxXCzhkQyzBwGqgD7DM+UY5n+VNfsW5LKaY/J0f0JJtva3ny5k4woEArJGolwH97DxfQINsHzYuzxEwBbjbVbniiCvbtiHALOfzK96bzg/w7sUZFzAfaAeEZUkcxXq+cvhbfgP0dVGu0M+ZVlVZLr/BL4twrrOdiDTAmvlwM1DDGHMGwPmzug0hTQKex5pLBaAaEGuMyXAu23HuGgLRwHRnFdpnIlIRm8+XMeY08B/gFHAGiAO2Yf/5yiqnc1SS3hNjsb7Ng81xicgg4LQxZle2TSXhfDUFbnBWg/4sIp2LKjZNHBZxsc72fsoiUglYADxjjLlUAuIZCEQZY7ZlXe2iaHGfOw+sy/ZPjDEdgESsahdbOdsLBmNVD9QGKgIDXBS1/X/NhZLwd0VEXgIygFmXV7koVixxiYgv8BLwd1ebXayz431QFauq7G/ANyIiFEFsmjgsEVj1lpcFA5E2xQKAiHhiJY1ZxpiFztXnRKSWc3stIKqYw+oBDBKRMGAuVnXVJKCKiFyeFMyOcxcBRBhjNjuX52MlErvPV1/ghDEm2hiTDiwErsP+85VVTufI9veEiIwGBgL3Gmcdi81xNcL6ErDL+R4IBraLSE2b47osAlhoLFuwagUCiyI2TRyW34Amzt4uXsBI4Fu7gnF+S/gcOGCMeTfLpm+B0c7no7HaPoqNMeYFY0ywMaYB1jlaY4y5F1gLDLMxrrNAuFjTEQPchDU/va3nC6uKqpuI+Dr/ppfjsvV8ZZPTOfoWuN/ZW6gbEHe5Sqs4iEh/YAIwyBiTlC3ekSLiLSIhQBNgS3HEZIzZY4ypboxp4HwPRGB1YjmLzefLaTHWlzlEpClWJ5HzFMU5K8rGm9L0wOoVcRirx8FLNsdyPdal5G5gp/NxK1Z7wmrgiPNngI0x9uKPXlUNnf+IR4F5OHt1FHM87YGtznO2GOuS3fbzBbwGHAT2Al9h9Wyx5XwBc7DaWtKxPvQezOkcYVVvfOR8P+wBQos5rqNY9fKX//8nZyn/kjOuQ8CA4owr2/Yw/mgcL7bzdZVz5gXMdP6vbQf6FNU50yFHlFJK5YlWVSmllMoTTRxKKaXyRBOHUkqpPNHEoZRSKk80cSillMoTTRxKKaXyRBOHUkqpPPl/TokuioOBrScAAAAASUVORK5CYII=\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('big_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
}