Supplementaries_GetParameters.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this notebook reads the experimental and simulation data from a defined directory\n",
    "#then a fitting on all data is performed using the best regularization type and strength\n",
    "#ultimately the correction parameters also published in this study are obtained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bussilab \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import scipy\n",
    "from scipy.optimize import minimize\n",
    "import cudamat as cm\n",
    "import concurrent.futures\n",
    "kBT=0.6 #kBT in kcal/mol\n",
    "np.random.seed(1995)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def concatenate_simulation_data(files):\n",
    "    result=None\n",
    "    for t in files:\n",
    "        if result is None:\n",
    "            result=np.load(t)\n",
    "            output=+result\n",
    "        else:\n",
    "            result=np.load(t)\n",
    "            output=np.concatenate((output,result))\n",
    "\n",
    "    return output\n",
    "\n",
    "def read_TLs(files):\n",
    "    corr=[]\n",
    "    for file in files:\n",
    "        with open(file, \"r\") as f:\n",
    "            for line in f:\n",
    "                nums=line.split()\n",
    "                if nums:\n",
    "                    corr.append(np.array([float(i) for i in nums]))\n",
    "    return np.array(corr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_TL_bias(Sequence,key,prune,skip,trajGAGA,FFprefactorsGAGA,trajUUCG,FFprefactorsUUCG):\n",
    "    \n",
    "    directories=[\"nh-n_0.5_nh-o_0.5_oh-bo-nbO_-0.5_dumping\", \n",
    "                 \"nh-n_1.0_nh-o_0.0_oh-bo-nbO_-0.5_dumping_bias-from-half-half\",\n",
    "                 \"nh-n_0.5_nh-o_0.5_oh-bo-nbO_-0.5_dumping_bias-from-one-zero\",\n",
    "                 \"nh-n_1.0_nh-o_0.0_oh-bo-nbO_-0.5_dumping\"\n",
    "                 ]\n",
    "\n",
    "    collection_weights=[]\n",
    "    for index,d in enumerate(directories):\n",
    "        Metadweight=[]\n",
    "        with open(\"./data_loaded/%s/%s/weights.rep0\" %(Sequence,d) ,\"r\") as fp:\n",
    "            for line in fp:\n",
    "                Metadweight.append(float(line))\n",
    "        collection_weights.append(np.array(Metadweight))\n",
    "\n",
    "    MetadPot1=np.concatenate((collection_weights[0],collection_weights[1]))\n",
    "    MetadPot2=np.concatenate((collection_weights[2],collection_weights[3]))\n",
    "\n",
    "    MetadPot=np.c_[ MetadPot1,MetadPot2 ] \n",
    "    MetadPot=kBT*np.log(MetadPot)\n",
    "    \n",
    "    if prune==True:\n",
    "        MetadPot=MetadPot[::skip,:]\n",
    "    \n",
    "    if Sequence=='GAGA':\n",
    "        if key == 'reference':\n",
    "            bias=np.matmul(trajGAGA,np.array(FFprefactorsGAGA).T)+MetadPot\n",
    "        if key == 'proposed':\n",
    "            bias=np.matmul(trajGAGA,(FFprefactorsGAGA-FFprefactorsGAGA[1]).T)+MetadPot\n",
    "        \n",
    "    if Sequence=='UUCG':\n",
    "        if key == 'reference':\n",
    "            bias=np.matmul(trajUUCG,np.array(FFprefactorsUUCG).T)+MetadPot\n",
    "        if key == 'proposed':\n",
    "            bias=np.matmul(trajUUCG,(FFprefactorsUUCG-FFprefactorsUUCG[1]).T)+MetadPot\n",
    "        \n",
    "\n",
    "    del MetadPot1\n",
    "    del MetadPot2\n",
    "    del MetadPot\n",
    "        \n",
    "    return bias   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def activate_cuda(traj,weights, USE_CUDA):\n",
    "\n",
    "    if USE_CUDA:\n",
    "        cm.init(max_ones=traj.shape[0])\n",
    "        cm_traj=cm.CUDAMatrix(traj.T)\n",
    "        cm_weights=cm.CUDAMatrix(np.reshape(weights,(1,-1)))\n",
    "        return cm_traj, cm_weights\n",
    "    else:\n",
    "        return traj,weights\n",
    "\n",
    "def compute_newweights(par,cm_traj,cm_weights,USE_CUDA):\n",
    "    if USE_CUDA:\n",
    "        cm_par=cm.CUDAMatrix(np.reshape(par,(1,-1)))\n",
    "        cm_correction=cm.dot(cm_par,cm_traj)\n",
    "        \n",
    "        m=float(cm_correction.min(axis=1).asarray()[0,0])\n",
    "        cm_correction.subtract(m)\n",
    "        cm_correction.mult(-1.0)\n",
    "        newweights=cm.exp(cm_correction)\n",
    "        newweights.mult_by_row(cm_weights)\n",
    "        newweights.mult(1.0/(cm.sum(newweights,axis=1).asarray()[0,0]))\n",
    "        return newweights \n",
    "    else:\n",
    "        correction=np.matmul(cm_traj,par)\n",
    "        correction-=np.min(correction)\n",
    "        newweights=np.exp(-correction)*cm_weights\n",
    "        return newweights/np.sum(newweights)\n",
    "    \n",
    "def function(par,obs,newweights,traj,cm_traj,USE_CUDA):\n",
    "    gradient=True\n",
    "    if USE_CUDA:\n",
    "        if(isinstance(obs,cm.CUDAMatrix)):\n",
    "            cm_obs=obs\n",
    "        else:\n",
    "            cm_obs=cm.CUDAMatrix(obs)\n",
    "        if(isinstance(newweights,cm.CUDAMatrix)):\n",
    "            cm_newweights=newweights\n",
    "        else:\n",
    "            cm_newweights=cm.CUDAMatrix(np.reshape(newweights,(1,-1)))\n",
    "        av_obs=np.reshape(np.array(cm.dot(cm_newweights,cm_obs).asarray()),(-1))\n",
    "    else:\n",
    "        av_obs=np.reshape(np.array(np.matmul(newweights,obs)),(-1))\n",
    "    \n",
    "    if not gradient:\n",
    "        return av_obs\n",
    "    \n",
    "    dav=np.zeros((len(par),obs.shape[1]))\n",
    "    \n",
    "    if USE_CUDA:\n",
    "        cm_traj_weighted=cm.empty((traj.shape[1],traj.shape[0]))\n",
    "        cm_traj_weighted.assign(cm_traj)\n",
    "        cm_traj_weighted.mult_by_row(cm_newweights)\n",
    "        av_traj=np.reshape(cm_traj_weighted.sum(axis=1).asarray(),(-1))\n",
    "        dav=-cm.dot(cm_traj_weighted,cm_obs).asarray()\n",
    "    else:\n",
    "        weighted_traj=np.copy(traj)\n",
    "        for ipar in range(len(par)):\n",
    "            weighted_traj[:,ipar]=newweights*traj[:,ipar]\n",
    "        av_traj=np.matmul(newweights,traj)\n",
    "        for ipar in range(len(par)):\n",
    "            for iobs in range(obs.shape[1]):\n",
    "                dav[ipar,iobs]=-np.sum(weighted_traj[:,ipar]*obs.T[iobs,:])\n",
    "\n",
    "    for ipar in range(len(par)):\n",
    "        for iobs in range(obs.shape[1]):\n",
    "            dav[ipar,iobs]+=av_obs[iobs]*av_traj[ipar]\n",
    "            \n",
    "    return av_obs,dav"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _logsum(ns):\n",
    "    max = np.max(ns)\n",
    "    ds = ns - max\n",
    "    sumOfExp = np.exp(ds).sum()\n",
    "    return max + np.log(sumOfExp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _softmax(x):\n",
    "    e_x = np.exp(x - np.max(x))\n",
    "    return e_x / e_x.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_and_grad(par,noe_exp,unoe_exp,sugar_exp,backbone1_exp,backbone2_exp,trajGACC,backbone1,backbone2,sugar,noe,unoe,trajGAGA,trajUUCG,weightsGACC,weightsGAGA,weightsUUCG,eps,prefactor,TLweight,obs_weight,populationGAGA,populationUUCG,regularization,alpha,gradient=True,regularize=True,components=False):\n",
    "    func=0.0\n",
    "    chi2_native_GAGA=0.0\n",
    "    chi2_native_UUCG=0.0\n",
    "    chi2_noe=0.0\n",
    "    chi2_unoe=0.0\n",
    "    chi2_backbone1=0.0\n",
    "    chi2_backbone2=0.0\n",
    "    chi2_sugar=0.0\n",
    "    \n",
    "    if gradient:\n",
    "        grad=np.zeros(len(par))\n",
    "        grad_backbone1=np.zeros(len(par))\n",
    "        grad_backbone2=np.zeros(len(par))\n",
    "        grad_sugar=np.zeros(len(par))\n",
    "        grad_noe=np.zeros(len(par))\n",
    "        grad_unoe=np.zeros(len(par))\n",
    "        grad_native_GAGA=np.zeros(len(par))\n",
    "        grad_native_UUCG=np.zeros(len(par))\n",
    "    \n",
    "    \n",
    "    \n",
    "    USE_CUDA=True\n",
    "    TLweight=(noe_exp.shape[0]+unoe_exp.shape[0] + backbone1_exp.shape[0]+ backbone2_exp.shape[0]+sugar_exp.shape[0])\n",
    "    \n",
    "    funcGAGA=0.0\n",
    "    gradGAGA=np.zeros(len(par))\n",
    "    if obs_weight[0] != 0:\n",
    "        chi2_native_GAGA=0.0\n",
    "        cm_traj, cm_weights,=activate_cuda(trajGAGA,weightsGAGA, USE_CUDA)\n",
    "        w=compute_newweights(par,cm_traj,cm_weights,USE_CUDA)\n",
    "\n",
    "        f=function(par,populationGAGA,w, trajGAGA,cm_traj, USE_CUDA)\n",
    "        target_population=0.5\n",
    "        diff=np.log(f[0]+eps)-np.log(target_population)\n",
    "        if(diff<.0):\n",
    "                chi2_native_GAGA+=diff**2\n",
    "                if gradient:\n",
    "                    grad_native_GAGA+=np.matmul(f[1], 2*(np.log(f[0]+eps)-np.log(target_population))/(f[0]+eps))\n",
    "        funcGAGA=obs_weight[0]*TLweight*chi2_native_GAGA\n",
    "        if gradient:\n",
    "            gradGAGA=obs_weight[0]*TLweight*grad_native_GAGA\n",
    "        \n",
    "    \n",
    "    funcUUCG=0.0\n",
    "    gradUUCG=np.zeros(len(par))\n",
    "    if obs_weight[1] != 0:\n",
    "        chi2_native_UUCG=0.0\n",
    "        cm_traj, cm_weights,=activate_cuda(trajUUCG,weightsUUCG, USE_CUDA)\n",
    "        w=compute_newweights(par,cm_traj,cm_weights,USE_CUDA)\n",
    "\n",
    "        f=function(par,populationUUCG,w, trajUUCG,cm_traj, USE_CUDA)\n",
    "        target_population=0.5\n",
    "        diff=np.log(f[0]+eps)-np.log(target_population)\n",
    "        if(diff<.0):\n",
    "                chi2_native_UUCG+=diff**2\n",
    "                if gradient:\n",
    "                    grad_native_UUCG+=np.matmul(f[1], 2*(np.log(f[0]+eps)-np.log(target_population))/(f[0]+eps))\n",
    "        funcUUCG=obs_weight[1]*TLweight*chi2_native_UUCG\n",
    "        if gradient:\n",
    "            gradUUCG=obs_weight[1]*TLweight*grad_native_UUCG\n",
    "    \n",
    "    func+=prefactor*(funcGAGA+funcUUCG)\n",
    "    if gradient:\n",
    "        grad+=prefactor*(gradGAGA+gradUUCG)\n",
    "    \n",
    "    USE_CUDA=True\n",
    "    cm_traj, cm_weights,=activate_cuda(trajGACC,weightsGACC, USE_CUDA)\n",
    "    w=compute_newweights(par,cm_traj,cm_weights,USE_CUDA)\n",
    "    \n",
    "    if obs_weight[2] != 0:\n",
    "        f=function(par,noe,w, trajGACC,cm_traj, USE_CUDA)\n",
    "        chi2_noe=np.sum(((f[0]**(-1/6)-noe_exp[:,0])/noe_exp[:,1])**2)/noe_exp.shape[0]\n",
    "        func+=obs_weight[2]*noe_exp.shape[0]*chi2_noe\n",
    "        if gradient:\n",
    "            grad_noe+=np.matmul(f[1],(2/noe_exp.shape[0] * ((f[0]**(-1/6)-noe_exp[:,0])/noe_exp[:,1]) * (-1/6) * f[0]**(-7/6))/noe_exp[:,1])\n",
    "            grad+=obs_weight[2]*noe_exp.shape[0]*grad_noe\n",
    "    \n",
    "    if obs_weight[3] != 0:\n",
    "        f=function(par,unoe,w, trajGACC,cm_traj, USE_CUDA)\n",
    "        chi2_unoe=0.0\n",
    "        for i in range(len(f[0])):\n",
    "            diff=f[0][i]**(-1/6)-unoe_exp[i,0]\n",
    "            if(diff<.0):\n",
    "                chi2_unoe+=(diff/unoe_exp[i,1])**2/unoe_exp.shape[0]\n",
    "                if gradient:\n",
    "                    grad_unoe+=f[1][:,i]*2/unoe_exp.shape[0]* diff/unoe_exp[i,1]* (-1/6) * f[0][i]**(-7/6) /unoe_exp[i,1]\n",
    "        func+=obs_weight[3]*unoe_exp.shape[0]*chi2_unoe\n",
    "        if gradient:\n",
    "            grad+=obs_weight[3]*unoe_exp.shape[0]*grad_unoe\n",
    "    \n",
    "    if obs_weight[4] != 0:\n",
    "        f=function(par,backbone1,w, trajGACC,cm_traj, USE_CUDA)\n",
    "        chi2_backbone1=np.sum(((f[0]-backbone1_exp[:,0])/backbone1_exp[:,1])**2)/backbone1_exp.shape[0]\n",
    "        func+=obs_weight[4]*backbone1_exp.shape[0]*chi2_backbone1\n",
    "        if gradient:\n",
    "            grad_backbone1+=np.matmul(f[1],2/backbone1_exp.shape[0] * (f[0]-backbone1_exp[:,0])/backbone1_exp[:,1]/backbone1_exp[:,1])\n",
    "            grad+=obs_weight[4]*backbone1_exp.shape[0]*grad_backbone1\n",
    "    \n",
    "    if obs_weight[5] != 0:\n",
    "        f=function(par,backbone2,w, trajGACC,cm_traj, USE_CUDA)\n",
    "        chi2_backbone2=np.sum(((f[0]-backbone2_exp[:,0])/backbone2_exp[:,1])**2)/backbone2_exp.shape[0]\n",
    "        func+=obs_weight[5]*backbone2_exp.shape[0]*chi2_backbone2\n",
    "        if gradient:\n",
    "            grad_backbone2+=np.matmul(f[1],2/backbone2_exp.shape[0] * (f[0]-backbone2_exp[:,0])/backbone2_exp[:,1]/backbone2_exp[:,1])\n",
    "            grad+=obs_weight[5]*backbone2_exp.shape[0]*grad_backbone2\n",
    "    \n",
    "    if obs_weight[6] != 0:\n",
    "        f=function(par,sugar,w, trajGACC,cm_traj, USE_CUDA)\n",
    "        chi2_sugar=np.sum(((f[0]-sugar_exp[:,0])/sugar_exp[:,1])**2)/sugar_exp.shape[0]\n",
    "        func+=obs_weight[6]*sugar_exp.shape[0]*chi2_sugar\n",
    "        if gradient:\n",
    "            grad_sugar+=np.matmul(f[1],2/sugar_exp.shape[0] * (f[0]-sugar_exp[:,0])/sugar_exp[:,1]/sugar_exp[:,1])\n",
    "            grad+=obs_weight[6]*sugar_exp.shape[0]*grad_sugar\n",
    "    \n",
    "    denominator=(obs_weight[0]*TLweight+obs_weight[1]*TLweight+obs_weight[2]*noe_exp.shape[0]+obs_weight[3]*unoe_exp.shape[0] + obs_weight[4]*backbone1_exp.shape[0]+ obs_weight[5]*backbone2_exp.shape[0]+obs_weight[6]*sugar_exp.shape[0])\n",
    "    func/=denominator\n",
    "    \n",
    "    \n",
    "    if gradient:\n",
    "        grad/=denominator\n",
    "    \n",
    "    if regularize:\n",
    "        if regularization == 'L1':\n",
    "            regL1=alpha * np.sum(np.sqrt(eps+par**2))\n",
    "            func+=regL1\n",
    "            if gradient:\n",
    "                grad_regL1=(alpha*par)/(np.sqrt(eps+par**2))\n",
    "                grad+=grad_regL1\n",
    "\n",
    "\n",
    "        if regularization == 'L2':\n",
    "            regL2=alpha * np.sum(par**2)\n",
    "            func+=regL2\n",
    "            if gradient:\n",
    "                grad_regL2=2*alpha*par\n",
    "                grad+=grad_regL2\n",
    "\n",
    "        if regularization == 'KishSize':\n",
    "            USE_CUDA=True\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGACC,weightsGACC, USE_CUDA)\n",
    "            wGACC=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGAGA,weightsGAGA, USE_CUDA)\n",
    "            wGAGA=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajUUCG,weightsUUCG, USE_CUDA)\n",
    "            wUUCG=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            t=[trajGACC, trajGAGA, trajUUCG] \n",
    "            w_list=[weightsGACC, weightsGAGA, weightsUUCG]\n",
    "            new_weights=[wGACC,wGAGA,wUUCG]\n",
    "            \n",
    "            \n",
    "            kish=np.sum((new_weights[0]**2))\n",
    "            kish_GAGA=np.sum((new_weights[1]**2))\n",
    "            kish_UUCG=np.sum((new_weights[2]**2))\n",
    "            \n",
    "            kish_size=_logsum(np.array([kish,kish_GAGA,kish_UUCG]))\n",
    "            chi2_kish=alpha*kish_size\n",
    "            func+=chi2_kish\n",
    "\n",
    "            if gradient:\n",
    "                grad_array=[0.0,0.0,0.0]\n",
    "\n",
    "                for index,elem in enumerate(new_weights):\n",
    "                    w_i=(elem**2)\n",
    "                    c_a_i=t[index]\n",
    "                    c_ave=np.matmul(new_weights[index],c_a_i)\n",
    "                    grad_array[index]=-2*alpha*np.matmul(w_i,(c_a_i - c_ave)) \n",
    "                    grad_array[index]*=_softmax(np.array([kish,kish_GAGA,kish_UUCG]))[index]\n",
    "                    \n",
    "                grad+=(np.sum(np.array(grad_array),axis=0)).reshape(12,)\n",
    "                \n",
    "        if regularization == 'RelativeKishSize':\n",
    "            USE_CUDA=True\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGACC,weightsGACC, USE_CUDA)\n",
    "            wGACC=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGAGA,weightsGAGA, USE_CUDA)\n",
    "            wGAGA=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajUUCG,weightsUUCG, USE_CUDA)\n",
    "            wUUCG=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            t=[trajGACC, trajGAGA, trajUUCG] \n",
    "            w_list=[weightsGACC, weightsGAGA, weightsUUCG]\n",
    "            new_weights=[wGACC,wGAGA,wUUCG]\n",
    "\n",
    "            Nf_list=[len(trajGACC),len(trajGAGA),len(trajUUCG)]\n",
    "            kish=np.sum((eps+(new_weights[0]**2))/(eps+w_list[0]))/Nf_list[0]\n",
    "            kish_GAGA=np.sum((eps+(new_weights[1]**2))/(eps+w_list[1]))/Nf_list[1]\n",
    "            kish_UUCG=np.sum((eps+(new_weights[2]**2))/(eps+w_list[2]))/Nf_list[2]\n",
    "\n",
    "            kish_size=_logsum(np.array([kish,kish_GAGA,kish_UUCG]))\n",
    "\n",
    "            chi2_kish=alpha*kish_size\n",
    "            func+=chi2_kish\n",
    "\n",
    "            if gradient:\n",
    "                grad_array=[0.0,0.0,0.0]\n",
    "\n",
    "                for index,elem in enumerate(new_weights):\n",
    "                    w_i=eps+(elem**2)\n",
    "                    w_i/=(eps+w_list[index])\n",
    "                    c_a_i=t[index]\n",
    "\n",
    "                    c_ave=np.matmul(new_weights[index],c_a_i)\n",
    "                    grad_array[index]=np.matmul(w_i,(c_a_i - c_ave))\n",
    "                    grad_array[index]*=-(2*alpha)\n",
    "                    grad_array[index]/=Nf_list[index]\n",
    "\n",
    "                    grad_array[index]*=_softmax(np.array([kish,kish_GAGA,kish_UUCG]))[index]\n",
    "\n",
    "                grad+=(np.sum(np.array(grad_array),axis=0)).reshape(12,)\n",
    "\n",
    "        if regularization == 'RelativeEntropy':\n",
    "            USE_CUDA=True\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGACC,weightsGACC, USE_CUDA)\n",
    "            wGACC=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajGAGA,weightsGAGA, USE_CUDA)\n",
    "            wGAGA=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            cm_traj, cm_weights,=activate_cuda(trajUUCG,weightsUUCG, USE_CUDA)\n",
    "            wUUCG=compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray()\n",
    "            t=[trajGACC, trajGAGA, trajUUCG] \n",
    "            w_list=[weightsGACC, weightsGAGA, weightsUUCG]\n",
    "            new_weights=[wGACC,wGAGA,wUUCG]\n",
    "\n",
    "            entropy=np.exp(np.sum(new_weights[0]*np.log((new_weights[0])/(w_list[0])+eps)))\n",
    "            entropy_GAGA=np.exp(np.sum(new_weights[1]*np.log((new_weights[1])/(w_list[1])+eps)))\n",
    "            entropy_UUCG=np.exp(np.sum(new_weights[2]*np.log((new_weights[2])/(w_list[2])+eps)))\n",
    "\n",
    "            relative_entropy=_logsum(np.array([entropy,entropy_GAGA,entropy_UUCG]))\n",
    "\n",
    "            chi2_entropy=alpha*relative_entropy\n",
    "            func+=chi2_entropy\n",
    "\n",
    "            if gradient:\n",
    "                grad_array=[0.0,0.0,0.0]\n",
    "\n",
    "                for index,elem in enumerate(new_weights):\n",
    "                    w_i=(elem) \n",
    "                    w_i_0=w_list[index]\n",
    "                    c_a_i=t[index]\n",
    "                    c_ave=np.matmul(eps+new_weights[index],c_a_i)\n",
    "                    grad_array[index]=np.array([entropy,entropy_GAGA,entropy_UUCG])[index]\n",
    "\n",
    "                    part1=-np.matmul(w_i*np.log((w_i/w_i_0)+eps),(c_a_i - c_ave))\n",
    "                    part2=alpha*part1-alpha*np.matmul(w_i,(c_a_i - c_ave))\n",
    "                    grad_array[index]*=part2\n",
    "\n",
    "                    grad_array[index]*=_softmax(np.array([entropy,entropy_GAGA,entropy_UUCG]))[index]\n",
    "\n",
    "                grad+=(np.sum(np.array(grad_array),axis=0)).reshape(12,)\n",
    "    \n",
    "    \n",
    "    if gradient:\n",
    "        return float(func), np.array(grad)        \n",
    "    \n",
    "    if components:\n",
    "        chi2_native_GAGA=float(obs_weight[0]*TLweight*chi2_native_GAGA)\n",
    "        chi2_native_UUCG=float(obs_weight[1]*TLweight*chi2_native_UUCG)\n",
    "        chi2_noe=float(obs_weight[2]*noe_exp.shape[0]*chi2_noe)\n",
    "        chi2_unoe=float(obs_weight[3]*unoe_exp.shape[0]*chi2_unoe)\n",
    "        chi2_backbone1=float(obs_weight[4]*backbone1_exp.shape[0]*chi2_backbone1)\n",
    "        chi2_backbone2=float(obs_weight[4]*backbone2_exp.shape[0]*chi2_backbone2)\n",
    "        chi2_sugar=float(obs_weight[5]*sugar_exp.shape[0]*chi2_sugar)\n",
    "        \n",
    "        return float(func), chi2_backbone1/denominator,chi2_backbone2/denominator,chi2_sugar/denominator,chi2_noe/denominator,chi2_unoe/denominator, chi2_native_GAGA/denominator,chi2_native_UUCG/denominator  \n",
    "    \n",
    "    return float(func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_single_min(data):\n",
    "    \n",
    "    chi2_col=data[1]\n",
    "    chi2_ref_col=data[2]\n",
    "    chi2_native_GAGA_col=data[3]\n",
    "    chi2_native_UUCG_col=data[4]\n",
    "    chi2_noes_col=data[5]\n",
    "    chi2_unoes_col=data[6]\n",
    "    chi2_sugar_col=data[7]\n",
    "    chi2_backbone1_col=data[8]\n",
    "    chi2_backbone2_col=data[9]\n",
    "    populationGAGA_col=data[10]\n",
    "    populationUUCG_col=data[11]\n",
    "    \n",
    "    temp=np.array([chi2_col,\n",
    "chi2_ref_col,\n",
    "chi2_native_GAGA_col,\n",
    "chi2_native_UUCG_col,\n",
    "chi2_noes_col,\n",
    "chi2_unoes_col,\n",
    "chi2_sugar_col,\n",
    "chi2_backbone1_col,\n",
    "chi2_backbone2_col,\n",
    "populationGAGA_col,\n",
    "populationUUCG_col],dtype='object')\n",
    "    \n",
    "    #comment out the following section for the CVSystem execution run\n",
    "    np.save(\"./data_produced/chi2.npy\",temp)\n",
    "    lambdas = np.asarray(data[0])\n",
    "    np.save(\"./data_produced/parameters.npy\", lambdas)\n",
    "    \n",
    "    #uncomment the following section for the CVSystem execution run without bounds\n",
    "    #np.save(\"./data_produced/SingleMinimization_errors_None.npy\",temp)\n",
    "    #lambdas = np.asarray(data[0])\n",
    "    #np.save(\"./data_produced/SingleMinimization_Lambdas_None.npy\", lambdas)\n",
    "    \n",
    "    #uncomment the following section for the CVSystem execution run with bounds\n",
    "    #np.save(\"./data_produced/SingleMinimization_bounds_errors_None.npy\",temp)\n",
    "    #lambdas = np.asarray(data[0])\n",
    "    #np.save(\"./data_produced/SingleMinimization_bounds_Lambdas_None.npy\", lambdas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reorder(input_array):\n",
    "    \n",
    "    temp=input_array\n",
    "    \n",
    "    input_array=np.array([temp[1],temp[4],temp[3],temp[5],temp[0],temp[2],temp[7],temp[10],temp[9],temp[11],temp[6],temp[8]])\n",
    "    \n",
    "    return input_array\n",
    "\n",
    "def opt_func_plumed(rdist):\n",
    "    \n",
    "        r0=1\n",
    "        d0=0.2\n",
    "        dmax=0.3\n",
    "        c=0.8\n",
    "        rdist=(rdist-d0)*1/r0\n",
    "\n",
    "        dmax2=dmax-d0\n",
    "\n",
    "        A = (-c*dmax2*dmax2)/((1-c)*dmax2*dmax2)\n",
    "        B = (2*dmax2)/((1-c)*dmax2*dmax2)   \n",
    "        C = -1/((1-c)*dmax2*dmax2)\n",
    "        D = 1/(c*dmax2*dmax2)\n",
    "\n",
    "        scale=1\n",
    "\n",
    "        if (rdist > dmax2):\n",
    "            result=0.0\n",
    "            return result\n",
    "        if (rdist <= 0): \n",
    "            result=-1.0\n",
    "        else:\n",
    "            result=-1.0\n",
    "            if rdist > c*(dmax2): \n",
    "                result+=(A + B*rdist + C*rdist*rdist)\n",
    "            elif rdist > 0 : \n",
    "                result+=D*(rdist*rdist) \n",
    "\n",
    "\n",
    "        result*=scale\n",
    "\n",
    "\n",
    "        return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main():\n",
    "    \n",
    "    print(\"Load experimental data\")\n",
    "    noe_exp=np.loadtxt(\"./data_loaded/GACC/exp_noes\",usecols=(2,3))\n",
    "    unoe_exp=np.loadtxt(\"./data_loaded/GACC/exp_unobs-noes\",usecols=(2,3))\n",
    "    sugar_exp=np.loadtxt(\"./data_loaded/GACC/exp_jcoupl-sugar\",usecols=(1,2))\n",
    "    backbone1_exp=np.loadtxt(\"./data_loaded/GACC/exp_jcoupl-backbone-1\",usecols=(1,2))\n",
    "    backbone2_exp=np.loadtxt(\"./data_loaded/GACC/exp_jcoupl-backbone-2\",usecols=(1,2))\n",
    "    print(noe_exp.shape)\n",
    "    print(unoe_exp.shape)\n",
    "    print(sugar_exp.shape)\n",
    "    print(backbone1_exp.shape)\n",
    "    print(backbone2_exp.shape)\n",
    "    print('There are no experimental population data for the TLs. The population is enforced to be 0.5 in the cost function.')\n",
    "    \n",
    "    #get the names of the different simulations and there forcefield prefactors\n",
    "    trajectoryNamesGACC=[]\n",
    "    trajectoryNamesGAGA=[]\n",
    "    trajectoryNamesUUCG=[]\n",
    "    FFprefactorsGACC=[]\n",
    "    FFprefactorsGAGA=[]\n",
    "    FFprefactorsUUCG=[]\n",
    "\n",
    "    with open(\"./data_loaded/coefficients.dat\",\"r\") as fp:\n",
    "        for line in fp:\n",
    "            if(line[0]=='#'):\n",
    "                continue\n",
    "            l=line.split()\n",
    "            file=l[0]+\".skip.npy\"\n",
    "            trajectoryNamesGACC.append(file)\n",
    "            FFprefactorsGACC.append(np.array(np.array(l)[1:],dtype='float'))\n",
    "\n",
    "    with open(\"./data_loaded/coefficients_TLs.dat\",\"r\") as fp:\n",
    "        for line in fp:\n",
    "            if(line[0]=='#'):\n",
    "                continue\n",
    "            if line.strip():\n",
    "                l=line.split()\n",
    "                file=l[0]+\".npy\"\n",
    "                if 'gaga' in line:\n",
    "                    FFprefactorsGAGA.append(np.array(l[1:],dtype='float'))\n",
    "                if 'uucg' in line:\n",
    "                    FFprefactorsUUCG.append(np.array(l[1:],dtype='float'))\n",
    "\n",
    "    firstFF='nh-n_0.5_nh-o_0.5_oh-bo-nbO_-0.5_dumping/'\n",
    "    secondFF='nh-n_1.0_nh-o_0.0_oh-bo-nbO_-0.5_dumping/'\n",
    "    trajectoryNamesGAGA.append('./data_loaded/GAGA/'+firstFF+'gHBfix-parameters_state.rep0')\n",
    "    trajectoryNamesGAGA.append('./data_loaded/GAGA/'+secondFF+'gHBfix-parameters_state.rep0')\n",
    "    trajectoryNamesUUCG.append('./data_loaded/UUCG/'+firstFF+'gHBfix-parameters_state.rep0')\n",
    "    trajectoryNamesUUCG.append('./data_loaded/UUCG/'+secondFF+'gHBfix-parameters_state.rep0')\n",
    "\n",
    "    #-------------------------------------------------------------------------------\n",
    "    #load observables\n",
    "    print('There are 14 simulations for GACC.')\n",
    "\n",
    "    trajGACC=concatenate_simulation_data(['./data_loaded/GACC/' + x for x in trajectoryNamesGACC])\n",
    "    backbone1=concatenate_simulation_data(['./data_loaded/GACC/' + re.sub(\"HBfix-energy\",\"jcop_backbone-1\",x) for x in trajectoryNamesGACC])\n",
    "    backbone2=concatenate_simulation_data(['./data_loaded/GACC/' + re.sub(\"HBfix-energy\",\"jcop_backbone-2\",x) for x in trajectoryNamesGACC])\n",
    "    sugar=concatenate_simulation_data(['./data_loaded/GACC/' + re.sub(\"HBfix-energy\",\"jcop_sugar\",x) for x in trajectoryNamesGACC])\n",
    "    noe=concatenate_simulation_data(['./data_loaded/GACC/' + re.sub(\"HBfix-energy\",\"noe\",x) for x in trajectoryNamesGACC])\n",
    "    unoe=concatenate_simulation_data(['./data_loaded/GACC/' + re.sub(\"HBfix-energy\",\"unoe\",x) for x in trajectoryNamesGACC])\n",
    "\n",
    "    #apply forward models\n",
    "    sugar=9.67*np.cos(sugar*np.pi/180)**2 - 2.03*np.cos(sugar*np.pi/180)\n",
    "    backbone1=9.7*np.cos(backbone1*np.pi/180)**2 - 1.8*np.cos(backbone1*np.pi/180)\n",
    "    backbone2=15.3*np.cos(backbone2*np.pi/180)**2 - 6.1*np.cos(backbone2*np.pi/180)+1.6\n",
    "    noe=noe**-6\n",
    "    unoe=unoe**-6\n",
    "    print('There are 2 simulations per TL.')\n",
    "    trajGAGA=read_TLs([x for x in trajectoryNamesGAGA])\n",
    "    trajUUCG=read_TLs([x for x in trajectoryNamesUUCG])\n",
    "    populationGAGA=trajGAGA[:,-1].reshape(-1,1)\n",
    "    populationUUCG=trajUUCG[:,-1].reshape(-1,1)\n",
    "    trajGAGA=trajGAGA[:,:12]\n",
    "    trajUUCG=trajUUCG[:,:12]\n",
    "\n",
    "    #to obtain the same parameters as published set prune=False\n",
    "    prune=False\n",
    "    if prune==True:\n",
    "        skip=5000\n",
    "        print(\"Pruning data, skip\")\n",
    "        trajGACC=trajGACC[::skip,:]\n",
    "        backbone1=backbone1[::skip,:]\n",
    "        backbone2=backbone2[::skip,:]\n",
    "        sugar=sugar[::skip,:]\n",
    "        noe=noe[::skip,:]\n",
    "        unoe=unoe[::skip,:]\n",
    "        populationGAGA=populationGAGA[::skip,:]\n",
    "        populationUUCG=populationUUCG[::skip,:]\n",
    "        trajGAGA=trajGAGA[::skip,:]\n",
    "        trajUUCG=trajUUCG[::skip,:]\n",
    "    else:\n",
    "        skip=0\n",
    "\n",
    "    print(\"Check shape of arrays:\")\n",
    "    print(trajGACC.shape)\n",
    "    print(backbone1.shape)\n",
    "    print(backbone2.shape)\n",
    "    print(sugar.shape)\n",
    "    print(noe.shape)\n",
    "    print(unoe.shape)\n",
    "    print(populationGAGA.shape)\n",
    "    print(populationUUCG.shape)\n",
    "    print(trajGAGA.shape)\n",
    "    print(trajUUCG.shape)\n",
    "    \n",
    "    \n",
    "    def wham(bias_ff,frame_weight=None,traj_weight=None,T=1.0,maxit=1000,threshold=1e-40,verbose=False):\n",
    "        \"\"\"WHAM tool\n",
    "\n",
    "           bias_ff[i,j] is the potential applied on frame i by Hamiltonian j\n",
    "           Returns an array with the weight of each frame\n",
    "\n",
    "        \"\"\"\n",
    "        import numpy as np\n",
    "        # allow lists\n",
    "        if not isinstance(bias_ff, (np.ndarray, np.generic) ):\n",
    "            bias_ff=np.array(bias_ff)\n",
    "\n",
    "        nframes=bias_ff.shape[0]\n",
    "        ntraj=bias_ff.shape[1]\n",
    "\n",
    "        # default values\n",
    "        if frame_weight is None:\n",
    "            frame_weight=np.ones(nframes)\n",
    "        if traj_weight is None:\n",
    "            traj_weight=np.ones(ntraj)\n",
    "\n",
    "        assert(len(traj_weight)==ntraj)\n",
    "        assert(len(frame_weight)==nframes)\n",
    "\n",
    "        # do exponentials only once\n",
    "        expv=np.exp((-bias_ff+np.min(bias_ff))/T)\n",
    "\n",
    "        Z=np.ones(ntraj)\n",
    "        Zold=np.copy(Z)\n",
    "\n",
    "        if(verbose):\n",
    "            print(\"WHAM: start\")\n",
    "        for it in range(maxit):\n",
    "            weight=1.0/np.matmul(expv,traj_weight/Z)\n",
    "            weight/=np.sum(weight*frame_weight)\n",
    "            Z=np.matmul(frame_weight*weight,expv)\n",
    "            Z/=np.sum(Z*traj_weight)\n",
    "            eps=np.sum(np.log(Z/Zold)**2)\n",
    "            Zold=Z\n",
    "            if verbose:\n",
    "                print(\"WHAM: iteration\",it,\"eps\",eps)\n",
    "            if(eps<threshold):\n",
    "                break\n",
    "        if(verbose):\n",
    "            print(\"WHAM: end\")\n",
    "        return weight*frame_weight\n",
    "    \n",
    "    weightsGACC=wham(np.matmul(trajGACC,np.transpose(FFprefactorsGACC-FFprefactorsGAGA[1])),verbose=False,threshold=1e-20,T=kBT)\n",
    "    weightsGAGA=wham(calculate_TL_bias('GAGA','proposed',prune,skip,trajGAGA,FFprefactorsGAGA,trajUUCG,FFprefactorsUUCG),verbose=False,threshold=1e-20,T=kBT)\n",
    "    weightsUUCG=wham(calculate_TL_bias('UUCG','proposed',prune,skip,trajGAGA,FFprefactorsGAGA,trajUUCG,FFprefactorsUUCG),verbose=False,threshold=1e-20,T=kBT)\n",
    "\n",
    "    bounds=[]\n",
    "    for i in range(12):\n",
    "        bounds.append((-1.666666666666666,+1.666666666666666)) #bounds in kbT units set to -1 to 1 in kcal/mol\n",
    "\n",
    "    eps=1e-30\n",
    "\n",
    "    prefactor=1.0\n",
    "    TLweight=(noe_exp.shape[0]+unoe_exp.shape[0] + backbone1_exp.shape[0]+ backbone2_exp.shape[0]+sugar_exp.shape[0])\n",
    "\n",
    "    \n",
    "\n",
    "    #regularization_list=['None','L1','L2','KishSize','RelativeKishSize','RelativeEntropy']\n",
    "    regularization_list=['RelativeKishSize']\n",
    "    regularization='RelativeKishSize'\n",
    "    alpha=0.0\n",
    "    \n",
    "    #uncomment the following line for the CVSystem execution run with and without bounds\n",
    "    #regularization_list=['None']\n",
    "\n",
    "    alpha_star_list=[18.679135990207847] #alpha* for RelativeKishSize regularization identified in this study\n",
    "\n",
    "    obs_weight=[1.0,1.0,1.0,1.0,1.0,1.0,1.0]\n",
    "    ref_chi2=func_and_grad(np.zeros(12),noe_exp,unoe_exp,sugar_exp,backbone1_exp,backbone2_exp,trajGACC,backbone1,backbone2,sugar,noe,unoe,trajGAGA,trajUUCG,weightsGACC,weightsGAGA,weightsUUCG,eps,prefactor,TLweight,obs_weight,populationGAGA,populationUUCG,regularization,alpha,gradient=False,regularize=False,components=False)\n",
    "    print(ref_chi2)\n",
    "\n",
    "    chi2_col=[]\n",
    "    chi2_ref_col=[]\n",
    "    chi2_native_GAGA_col=[]\n",
    "    chi2_native_UUCG_col=[]\n",
    "    chi2_noes_col=[]\n",
    "    chi2_unoes_col=[]\n",
    "    chi2_sugar_col=[]\n",
    "    chi2_backbone1_col=[]\n",
    "    chi2_backbone2_col=[]\n",
    "    populationGAGA_col=[]\n",
    "    populationUUCG_col=[]\n",
    "    collection=[]\n",
    "    for r,reg in enumerate(regularization_list):\n",
    "        regularization=reg\n",
    "        print(regularization)\n",
    "\n",
    "        alpha=alpha_star_list[r]\n",
    "        print(alpha)\n",
    "\n",
    "        #comment the 'bounds=bounds' in following line for the CVSystem execution run without bounds\n",
    "        m=minimize(func_and_grad,x0=np.zeros(12), args=(noe_exp,\n",
    "unoe_exp,\n",
    "sugar_exp,\n",
    "backbone1_exp,\n",
    "backbone2_exp,\n",
    "trajGACC,\n",
    "backbone1,\n",
    "backbone2,\n",
    "sugar,\n",
    "noe,\n",
    "unoe,\n",
    "trajGAGA,\n",
    "trajUUCG,\n",
    "weightsGACC,\n",
    "weightsGAGA,\n",
    "weightsUUCG,\n",
    "eps,\n",
    "prefactor,\n",
    "TLweight,\n",
    "obs_weight,\n",
    "populationGAGA,\n",
    "populationUUCG,\n",
    "regularization,\n",
    "alpha), method='L-BFGS-B', jac=True,bounds=bounds)\n",
    "        print(m.message, m.success)\n",
    "        print(m.x)\n",
    "\n",
    "        collection.append(list(m.x))\n",
    "\n",
    "        t_chi2, t_chi2_backbone1,  t_chi2_backbone2, t_chi2_sugar, t_chi2_noes,  t_chi2_unoes, t_chi2_native_GAGA,t_chi2_native_UUCG  =func_and_grad(m.x,noe_exp,\n",
    "unoe_exp,\n",
    "sugar_exp,\n",
    "backbone1_exp,\n",
    "backbone2_exp,\n",
    "trajGACC,\n",
    "backbone1,\n",
    "backbone2,\n",
    "sugar,\n",
    "noe,\n",
    "unoe,\n",
    "trajGAGA,\n",
    "trajUUCG,\n",
    "weightsGACC,\n",
    "weightsGAGA,\n",
    "weightsUUCG,\n",
    "eps,\n",
    "prefactor,\n",
    "TLweight,\n",
    "obs_weight,\n",
    "populationGAGA,\n",
    "populationUUCG,\n",
    "regularization,\n",
    "alpha,gradient=False,regularize=False,components=True)\n",
    "\n",
    "        chi2_col.append(t_chi2)\n",
    "        chi2_ref_col.append(ref_chi2)\n",
    "\n",
    "        chi2_native_GAGA_col.append(t_chi2_native_GAGA)\n",
    "        chi2_native_UUCG_col.append(t_chi2_native_UUCG)\n",
    "        chi2_noes_col.append(t_chi2_noes)\n",
    "        chi2_unoes_col.append(t_chi2_unoes)\n",
    "        chi2_sugar_col.append(t_chi2_sugar)\n",
    "        chi2_backbone1_col.append(t_chi2_backbone1)\n",
    "        chi2_backbone2_col.append(t_chi2_backbone2)\n",
    "\n",
    "        par=m.x\n",
    "        USE_CUDA=True\n",
    "\n",
    "        cm_traj, cm_weights,=activate_cuda(trajGAGA,weightsGAGA, USE_CUDA)\n",
    "        new_weights=np.reshape(compute_newweights(par,cm_traj,cm_weights,USE_CUDA).asarray(),(-1,))  \n",
    "\n",
    "        cm_traj, cm_weights,=activate_cuda(trajUUCG,weightsUUCG, USE_CUDA)\n",
    "        new_weights=np.reshape(compute_newweights(par,cm_traj,cm_weights,True).asarray(),(-1,)) \n",
    "\n",
    "        populationGAGA_col.append(np.matmul(new_weights,populationGAGA))\n",
    "        populationUUCG_col.append(np.matmul(new_weights,populationUUCG))\n",
    "\n",
    "    save_list=[collection,\n",
    "    chi2_col,\n",
    "    chi2_ref_col,\n",
    "    chi2_native_GAGA_col,\n",
    "    chi2_native_UUCG_col,\n",
    "    chi2_noes_col,\n",
    "    chi2_unoes_col,\n",
    "    chi2_sugar_col,\n",
    "    chi2_backbone1_col,\n",
    "    chi2_backbone2_col,\n",
    "    populationGAGA_col,\n",
    "    populationUUCG_col] \n",
    "    save_single_min(save_list)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Load experimental data\n",
      "(20, 2)\n",
      "(285, 2)\n",
      "(12, 2)\n",
      "(8, 2)\n",
      "(9, 2)\n",
      "There are no experimental population data for the TLs. The population is enforced to be 0.5 in the cost function.\n",
      "There are 14 simulations for GACC.\n",
      "There are 2 simulations per TL.\n",
      "Check shape of arrays:\n",
      "(1400000, 12)\n",
      "(1400000, 8)\n",
      "(1400000, 9)\n",
      "(1400000, 12)\n",
      "(1400000, 20)\n",
      "(1400000, 285)\n",
      "(1000000, 1)\n",
      "(1000000, 1)\n",
      "(1000000, 12)\n",
      "(1000000, 12)\n",
      "21.803205490112305\n",
      "RelativeKishSize\n",
      "18.679135990207847\n",
      "b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH' True\n",
      "[-1.66666667 -1.08691868  0.20668562 -1.66666667  1.29619341 -1.66666667\n",
      " -1.66666667  1.41575784 -1.66666667 -1.66666667  1.4333629  -1.60218027]\n"
     ]
    }
   ],
   "source": [
    "executor=concurrent.futures.ProcessPoolExecutor(max_workers=1)\n",
    "future = executor.submit(main)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Future at 0x7f409435ffa0 state=running>\n"
     ]
    }
   ],
   "source": [
    "print(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAADtCAYAAAB6WSfQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe40lEQVR4nO3de3hU9YH/8U8yIUAIICTcLyIUQkBAYUVBLl4BtYDA2ha5BS8s6G6FVp9dFS2LsLvaolbXC2gFFsU+v/4QIYJoseVaChiBtIQIqyQEEsIlCCEJJAyzf4gnMyEzmWRm4PtN3q/n4XnOmfmefL/mQ+JnDjPnRHk8HgEAAACmir7aCwAAAAACobACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMFhPoyaiokVzzymJud+rVXgIAAEBQoqMV5fe5K7kQAAAAoLoorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxmRWFNTGyiESP66bnnfqZVq2YrN3epPJ5U58+UKXdekXW0aNFUv/jF/dq8+UUdObJEJSUrlJ39ntavn6dp00aoceOGV2QdNtqxY4cef/wx9erVU82bX6OmTRurR4/uSkmZovXr10d8/mPHjunllxdoyJBBat++reLiGqhTp466++47tXDhQp05cybia7AZ+dmL7OxGfvYiu/CK8ng8/p+MGun/ySugVatr9Ne//kadOrUKOC4l5VUtXfpFRNcyduxALVr0uBISmvgdk519TJMmvazNm/dGdC3BcrtTr/YSVFxcrCef/KUWLnw74LgxY8Zq0aJ31Lx587CvYcWKFZo+fZoKCgr8junYsaOWLl2mIUOGhH1+m5GfvcjObuRnL7KruehoRfl7zujCeu21LZWV9bsqx0W6sP7sZ0P04YdP+TyWmXlYR4+eUqdOLX0K9fnzZbrrrtnasiUjYusJ1tUurG63W/ffP1pr165xHouLi1OPHj0UExOjjIwMn1d4/fr104YNmxQXFxe2NXz44YeaOPFBn8eSkpLUunVrZWVlKTs723k8NjZWn3++XoMHDw7b/DYjP3uRnd3Iz15kF5pAhdWKtwRI0rFj3+nTT9P0wgu/1+jR867YvMnJHfTeez939jMzD6tfv5lKTp6h229/Rtdd94juvvs5HT16SpJUv349rVz5rBIT/Z+JrSvmzv13nx/aRx+dpkOHDmv79p3aunWbDh/O1bPPznaeT0tL02OPzQjb/BkZGXrkkYec/aSkJO3Y8aUyMjL1pz9t0LffZmndus/VqtX3LzhKS0s1btwYHT9+PGxrsBn52Yvs7EZ+9iK7yDH6DGvjxg01bNiN2rnzgA4d8v1mejzlZw8jeYb1o4+e0ZgxAyRJx4+f1vXX/7OOHfvusnE9enRUWtoratAgVpL08ssf65e/rPrscCRdzTOseXl56tbtRyouLpYkTZw4SUuX/k+lY59//jnNn//9i5CoqCilpe1Snz59Ql7D2LFjtGrVx5KkxMREpaf/3fkh9bZ3717ddFM/nT9/XpI0c+YsLVjwcsjz24z87EV2diM/e5Fd6Kw9w1pYWKIVK/5yWVm9UpKS2jtlVZJmz36/0rIqSRkZh/Tqq6ud/Rkz7lGTJuE7xW+b1177rfNDGxcXp1deedXv2Nmzn1OHDh0kSR6PRy+99GLI82dmZjo/tJI0d+68Sn9oJalnz5564omZzv7bb7+l06dPh7wGm5GfvcjObuRnL7KLLKML69U2dmx5WS0sLNYHH2wIOH7RonXOdsOG9XXvvf8QqaUZb+XKj5ztBx74ScA3lcfGxiolZaqzv2bNJyotLQ3b/PHx8ZowYULA8Y8+Os3ZPnfunNauXRvS/LYjP3uRnd3Iz15kF1kU1gB+/OObnO0tW/apqOhcwPEHD+YrM/NwpcfXJfv379eBAwec/eHDR1R5zIgR9zjbhYWF2rhxY0hrWLPmE2f71lsHKT4+PuD4zp07KykpqdLj6xrysxfZ2Y387EV2kUdhDaB3707O9rZtmUEd4z2uT5/rwr0kK+zZs8dnf8CAAX5Gluvbt69iY2Od/fT0PQFGB+bxeJSenl6t+SXpllvKx4Uyv+3Iz15kZzfysxfZRR6F1Y8OHVooPr78RgDffJMX1HHe47p2bavo6Lr3Lc7M3Odsx8bGOu/TCaTiuH379gUYHVhOTo6Kioqc/c6duwR1XJcu5eMOHDggt9td4zXYjPzsRXZ2Iz97kV3k1b02FaROnVr67Af7wS/vcfXr11PbtuG/ILDpsrKynO327dsrKsrvh/58dOzY0dnOzs7yP7Aa81f8uoF06FA+rrS0VLm5uTVeg83Iz15kZzfysxfZRR6F1Y+Kn/A/fbo4qOPOnCnx2a+Lt2v1vihy06ZNgz6ucePya9cWFhaGZf7qrKFJE99r54ayBpuRn73Izm7kZy+yizwKqx+NGtX32T93LrhP75WUnPfZj49vELY12cL7nyUaNAj+v79hw/Jyf/bs2bDMX501eM8f6hpsRn72Iju7kZ+9yC7yKKx+1KsX47N/4UJw7+uoOC42tl7Y1mSLCxfKnO2YmJgAI315jw3l8h5lZWU++8GuoeK4UC8xYivysxfZ2Y387EV2kUdh9aO42PdM6Q93sKpKxXFnz5b4GVl7NWxY/naKc+cCXwrMm/fYqi7HEUjFezIHu4aK40JZg83Iz15kZzfysxfZRR6F1Y+zZ31DbNgwuMIaF+f7VoKKX6cu8P4LX1ISfGH/4Q4hFb9GKPNXZw3e84e6BpuRn73Izm7kZy+yizwKqx8nTvi+gblNm+A+7d+6dTOf/ZMnz/gZWXslJiY623l5wV0OTJLy8486282bJ4Rl/uqswXt+SUpIqPkabEZ+9iI7u5Gfvcgu8iisfuzff0QXL1509jt2bBHUcR06lP+lyc8/pe++KwowunbyvnPGyZMnL3sF509OTo6z3b179xrP361bN59LiuTkHKr2/C1btlSzZs0CjK69yM9eZGc38rMX2UUehdWPoqJzysk54ezfcENwd6268cbOzva+fYcDjKy9kpN7+Ozv3r27ymOOHDmi48fLr2HbvXtyjeePj4/3uRhzMPNL0q5du5zt5OSaz2878rMX2dmN/OxFdpFHYQ1g06a9zvagQT0CjPxeTIxLN99c/irL+/i6pH///j6X1Ni6dUuVx2zevNlnf8iQISGtYfDg8uODmb+srEw7dmyv9Pi6hvzsRXZ2Iz97kV3kUVgDWLXqr852jx4ddcMNnQOMlkaNutnnhgMff/zXAKNrr0aNGumOO+509pcv/6DKY7zH9OrVS507B/5eV2XUqNHO9r59+3xeRVZm9erVPhdMHj36/pDmtxn52Yvs7EZ+9iK7yKOwBrB2bZry8085+7Nn/9Tv2OjoaP3rv45z9vfsOahdu76J6PpMlpIy1dlOT09Xamqq37FfffWV1q37tNJja+ree+9Vy5blt9edP3+e37Fut1u//vWLzn7v3r3Vt2/fkNdgM/KzF9nZjfzsRXaRVScLq8eT6vxZvHim33ElJec1f/4fnP1x4wbq5z8fWenY3/zmIfXv383Znz17WdjWa6OxY8f6/OWfMeOflJmZedm43NxcTZ48UW739zdcaNeunaZPn+H367pcUc6fqVNT/I6Li4vTM8886+yvXPmRXnvtt5eN83g8euqpJ7Vz507nsblz/f+Q1xXkZy+ysxv52YvsIivK4/H4fzJqpP8nr5BFi/5Zkybdftnj3hfoLyu7ILf74mVjkpKm69Ch45c97vGUv+pZsuQLTZ36qt/5Y2JcWr9+noYOvd55bPXq7Vq+fKOOHj2lTp1a6ZFHhvm8x/WDDzZo4sQFVf63RZrb7f/V3ZWwa9cuDRkyyPm0ZJMmTTRjxmMaNGiwXC6Xvvxyp95447+Vn58vSXK5XFq1KlX33HOP36/pcpV/CnLy5ClavHiJ37FlZWUaNuwubdq0yXnsxz8eqfHjH1Tr1q2VlZWl9957V1u3bnWeHz/+Qb3/ftX/lFMXkJ+9yM5u5GcvsgtNdLSi/D1nfGFdvHimUlLurHpgJTp1eljZ2ccue7w6hVWSmjdvrPXrX9CNN3apcs7PPvtKo0a9oNLSC9Veb7hd7cIqSampqZowYfxl9zmuKCYmRq+//oamTZsWcFx1fnCl7y8vMmzYXUF9YvLuu4dp1arVql+/fpVj6wrysxfZ2Y387EV2NReosNbJtwRUV0FBoW6++UnNn///VFBQWOmY3NyTmjXrXY0Y8SsjyqopRo4cqbS0XRo+fISioyv/6zZw4EBt2fKXKn9oayIhIUHbtm3X008/4/f6cm3atNGCBS/r00/XGfNDawrysxfZ2Y387EV2kWH8GVbT1KsXo9tu66Xrrmul5s3jdfz4GR04kKstWzJ8bjRgAhPOsHrLzc3Vli1blJt7RG63W23btlP//v3VpUvVZ67DobS0VBs2bNDBgwd16lSBEhNbqGvXrho0aJBcLtcVWYPNyM9eZGc38rMX2VWP1W8JQM2ZVlgBAAD84S0BAAAAsBaFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGi/J4PH6fvHhR/p8EAKAWmuuKutpLAOqkOR6P3x8+zrACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARou52gsI1o4dO7R06RJt2rRRR44ckdvtVrt27dS//82aOHGS7rrrrojOf+zYMb3//jJ9/PFKffvttyooKFDLli3VtWtX/eM//kTjx49XkyZNIroGm5Gf3cjPXmRX+0RFRSlp1Cj1njhRrfr0UZN27VRWXKwzhw/rm88/1+4lS3R8376wz3vNtddqZlZWtY9b/2//pi0vvhj29diK/GomyuPx+H3y4kX5f/IKKS4u1pNP/lILF74dcNyYMWO1aNE7at68edjXsGLFCk2fPk0FBQV+x3Ts2FFLly7TkCFDwj6/zcjPbuRnL7KrubmuqKu9BL+aduyoscuW6doA36+LFy5o07x52vjCC/JcvBi2uW0vPCYgv8DmeDx+f/iMLqxut1v33z9aa9eucR6Li4tTjx49FBMTo4yMDJ05c8Z5rl+/ftqwYZPi4uLCtoYPP/xQEyc+6PNYUlKSWrduraysLGVnZzuPx8bG6vPP12vw4MFhm99m5Gc38rMX2YXG1MLaqGVLPbp9u67p1Ml5rOj4cZ3IzFT9xo3VokcPuWJjned2vvmm1jz+eNjmr1h4/nfduqCO2/Xee9r7hz+EbR22Ir+qWVtYf/Wr5zVv3gvO/qOPTtN//ud/qVmzZpKkoqIivfjif2n+/HnOmEmTJmvJkqVhmT8jI0M33dRP586dk/T9L9tlyz5Qv379nDF//OMfNWXKJOXn50uSEhIStHfvPrVo0SIsa7AZ+dmN/OxFdqExtbCmbNigTkOHSpLcpaVaN2uW0hYt0sULFyRJTdq1031vvqmkUaOcYz6eOlW7lywJy/wVC8+cKDO/T6Yiv6pZWVjz8vLUrduPVFxcLEmaOHGSli79n0rHPv/8c84v3qioKKWl7VKfPn1CXsPYsWO0atXHkqTExESlp/9drVq1umzc3r17ddNN/XT+/HlJ0syZs7Rgwcshz28z8rMb+dmL7EJnYmFNGjVK41etcvZXPfSQdi1efNm4qOhoTfniC3W67TZJ0pkjR/Raly66cOl7HAobCo+pyC84gQqrsVcJeO213zq/cOPi4vTKK6/6HTt79nPq0KGDJMnj8eill0J/r0VmZqbzC1eS5s6dV+kvXEnq2bOnnnhiprP/9ttv6fTp0yGvwWbkZzfysxfZ1U6Dn37a2c7Ztq3SsiNJnosX9cmMGc5+k3bt1Gfy5IivD4GRX+iMLawrV37kbD/wwE8CfiAgNjZWKSlTnf01az5RaWlp2OaPj4/XhAkTAo5/9NFpzva5c+e0du3akOa3HfnZjfzsRXa1T+O2bdWuf39n/8u3A3+Q7kRmprI2bnT2k8eNi9jaUDXyCw8jC+v+/ft14MABZ3/48BFVHjNixD3OdmFhoTZ6hV0Ta9Z84mzfeusgxcfHBxzfuXNnJSUlVXp8XUN+diM/e5Fd7dTtvvsUFV3+v+tgPizzv59+6mxfd/vtqhfGD9ShesgvPIwsrHv27PHZHzBgQJXH9O3bV7Fen65LT98TYHRgHo9H6enp1Zpfkm65pXxcKPPbjvzsRn72IrvaqZXX+4pPffutio4dq/KYnG3bnG1XbKxaJCdHZG2oGvmFh5E3DsjMLL9gbmxsrPMeq0B+GPfNN99IkvaFcNHdnJwcFRUVOfudO3cJ6rguXcrHHThwQG63Wy6Xq8brsBX52Y387EV2tZN3WSm4lFNVTlUYl5icrNy0tLCua8zSpWp/yy1q3LatoqKjVVJQoJMHDih740btWbZMp779Nqzz2Yr8wsPIM6xZXp9ia9++vaKC/CRbx44dne3s7Cz/A6sxf8WvG0iHDuXjSktLlZubW+M12Iz87EZ+9iK72sn7up2nDx0K6pjC3Fy5y8oq/Rrh0mfyZCV066bY+HjVi4tTk/btdd3tt+u2OXP0L19/rfsXL1ZsFW8JqQvILzyMPMPqfUHrpk2bBn1c48blt/crLCwMy/zVWUPF2wuGsgabkZ/dyM9eZFc71ff6/pwP8ioKHo9HZUVFcl1zzfdfo3HjsK+rMC9P3x08qLLiYjW45holJicrtlEjSVJ0TIxuSElR+wEDtGToUJ29dL3duoj8wsPIwur9T0oNGjQI+riGDRs622fPng3L/NVZg/f8oa7BZuRnN/KzF9nVTvUulQhJunDpZgzBKCspUYNLhSccZ8o8Ho8Ob9+uXb/7nfavWaPCCmfCo2Ni1PWee3TH/Plq1auXJCkxKUnjU1P1u4EDnQvk1zXkFx5GviXgwoXy0+AxMcF3au+xoVyapczrNHx11lBxXKiXh7EV+dmN/OxFdrWTq149Z7s6pcF7rPctP2vq9KFDeveWW5T2zjuXlZ0f5vs6NVXv9O+vA16fcm93003q+/DDIc9vK/ILDyMLa8OG5ZdvOFeNVyPeY6u6lEogFe+nHewaKo4LZQ02Iz+7kZ+9yK52Krt0IwhJiqnGmXPvsaVX8Kz1hXPn9P9/+lOff0a++Yknrtj8piG/8DDyLQHev6xKSkqCPq7Y6y9FKL/wKh4b7Bq85w91DTYjP7uRn73IznytevXS3S+9VOW49+8pvz5u6dmzzvsgYyq8fSIQ72t3XsnCI0nnCwv15Vtv6bY5cyR9/0n5ph066HROzhVdR7iR39XLz8jCmpiY6Gzn5eUFfVx+/lFnu3nzhLDMX501eM8vSQkJNV+DzcjPbuRnL7IzX8PmzfWjEVXf0MFb8YkTaty2rSSpcZs2QR0TGx/vfIBGkopPnqzWnOFw8M9/dgqPJCV062Z9YSW/q5efkW8J8L7rycmTJy979e1Pjtc3snv37jWev1u3bj6Xg8nJCe4yFN7zt2zZUs2aNavxGmxGfnYjP3uRXe104uuvne2mQV4qrGmFa/CeyMwM65qCcbbCC5a4Ci9o6gryCw8jz7AmJ/fw2d+9e7cGDhwY8JgjR47o+PHjzn737jW/K0R8fLw6dOigQ5eul7Z79+6gjtu1a5eznVyH70pBfnYjP3uRnfmyNm7UnCCvj/uD4xkZznaLnj0V7XLpotsd8JjWN97os38ihBtC1FTF24mWVeNtKqYiv6vHyDOs/fv397kcytatW6o8ZvPmzT77Q4YMCWkNgweXHx/M/GVlZdqxY3ulx9c15Gc38rMX2dVO2Zs2OduxjRqpTd++VR5z7eDBzvZ32dn6Ljs7ImsLpEUP3xdQRQZcy/NqIL/wMLKwNmrUSHfccaezv3z5B1Ue4z2mV69e6ty5c0hrGDVqtLO9b98+nzMAlVm9erXPxa5Hj74/pPltRn52Iz97kV3tlL1xo0oKCpz9XhMmBBzvqldPyePGOftfr1oVsbUFcv348c52WXGx8qr4u1BbkV94GFlYJSklZaqznZ6ertTUVL9jv/rqK61bV37NMO9ja+ree+9Vy5Ytnf358+f5Het2u/XrX7/o7Pfu3Vt9g3gFVZuRn93Iz15kV/tcdLuV/kH5C4sbp051PsRTmX+YPl2NWrRw9ncvXRrR9VUmadQoJY0c6ewfWLtW7jp6fV3yCw9jC+vYsWN9fnHNmPFPyqzkTce5ubmaPHmi3JfeD9KuXTtNnz7D79d1uaKcP1OnpvgdFxcXp2eeedbZX7nyI7322m8vG+fxePTUU09q586dzmNz5/r/BV1XkJ/dyM9eZFc7bZ4/X6WX7iRWv0kTjVu+vNK7H3W89Vbd+R//4exnrFihvK++8vt1Ow0dqjkej/PnhilTKh3XokcPjXr3XSV061blWntPmKBxy5c7+xcvXNCff/WrKo+rzcgvdEZ+6EqSoqKitGjRuxoyZJCKi4uVl5enAQNu1owZj2nQoMFyuVz68sudeuON/1b+pfdVuFwuLVz4TrVuSRjI9Okz9NFHK7Tp0vtPZs2aqS+++ELjxz+o1q1bKysrS++99662bt3qHDN+/IMa6fWqpK4iP7uRn73IrnY6m5+vz2bN0shFiyR9X1Sm796tHa+/rvy//U2x8fH60fDhuvGhh5wLzp89elSfzZoVlvld9eqp78MPq+/DDys3LU0H//Qn5aenqyg/X2XFxarftKla9e6tng88cNl7NNfNnOnzwaO6iPxCF+XxePw+efGi/D95haSmpmrChPGX3aO6opiYGL3++huaNm1awHEuV/mn+yZPnqLFi5cEHH/y5EkNG3ZXUJ92vfvuYVq1arXq169f5di6gvzsRn72Iruam+uq3qfAr6Tb587V0Oeeq3Jc0fHjWn7ffTridQa7Mp2GDlXKhg3O/scpKZX+E3TrPn00PcirPvygrLhYn/3iF/py4cJqHVebkV9gczwevz98xr4l4AcjR45UWtouDR8+QtHRlS934MCB2rLlL1X+wq2JhIQEbdu2XU8//YzfawO2adNGCxa8rE8/XWfML1xTkJ/dyM9eZFc7/fn55/X+iBHKT0+v9PkL58/r77//vd7q1avKslMdhXl52r1kiU7u31/l2PNnzmjnm2/qrd69KasVkF/NGX+G1Vtubq62bNmi3Nwjcrvdatu2nfr3768uXbpckflLS0u1YcMGHTx4UKdOFSgxsYW6du2qQYMGyeVyXZE12Iz87EZ+9iK76jH5DKu3ltdfr9Y33KDGbdvqQkmJzhw+rKwNG1Ry6lRE541LSFCrPn3UqEULxSUmqn7TpiorLlZJQYGO/e1vyk9Pr/I6oyC/ygQ6w2pVYQUAINJsKaxAbWP1WwIAAABQt1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMRmEFAACA0SisAAAAMBqFFQAAAEajsAIAAMBoFFYAAAAYjcIKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARqOwAgAAwGgUVgAAABiNwgoAAACjUVgBAABgNAorAAAAjEZhBQAAgNEorAAAADAahRUAAABGo7ACAADAaBRWAAAAGI3CCgAAAKNRWAEAAGA0CisAAACMFuXxeK72GgAAAAC/OMMKAAAAo1FYAQAAYDQKKwAAAIxGYQUAAIDRKKwAAAAwGoUVAAAARvs/p065kpGdv/4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEYCAYAAAANjbKIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABSXUlEQVR4nO2dd3gc1dWH36PeJcuWeweCMd2YmkAg9BIgEJOAE3CAUEJCAgHCBwQcIEAg1NC7SWihJJTEdFNDsymmdwOuGFfZsixLOt8fd8Y7Ws1Ku6tZ7Uo+7/PMs7Nz79y5O1rNb++9p4iqYhiGYRi5Sl62O2AYhmEYHWFCZRiGYeQ0JlSGYRhGTmNCZRiGYeQ0JlSGYRhGTlOQ7Q4YmaNfv346cuTIbHfDMAxjLTNmzPhWVetSOceEqhczcuRIpk+fnu1uGIZhrEVEvkz1HJv6SwERKRORvUXkLBF5UES+FBH1tskRXWOAiFwqIh+JyCoRWSwiL4jI0SIiUVzDMAyjJ2EjqtTYBvhvphoXka2Ax4G+3qEVQCXwPW+bICL7q+rqTPXBMAwj17ARVeosAZ4GLgEOBeZH0aiIVAOP4kTqQ2BrVa0EyoFfA2uAPYDLo7ieYRhGT8FGVKnxgqrWBg+IyEURtX0KMBBYBeyjql8AqGoTcI2IVAEXAMeIyBWq+nFE1zUMw8hpbESVAqraksHmD/de7/FFKo6/4aYC84GJGeyHYRhGTmEjqhxARDYEhntvp4bVUdUVIvICsDduCvCcbuqeYazTrF69hscem8Hzz7/H/PlLaWpag8Xybsttt51IZWVZxto3ocoNNgnsv9tBvXdxQjU2s90xDENVueWWJznnnLuYO3dxtruT09x44wkZbd+EKjcYHNif00E9v6xKRCpUdUUG+2QY6yz19Q1MnHgZjzzyWra7YmBClStUBvYbOqgXLKvErVm1QUSOAY4BGD58eHyxYRidsHTpCvbcczKvvWb2SrmCGVP0MlT1RlUdr6rj6+pSilJiGOs8TU1rOPjgi0ykcoxuG1GJSIuq5nfX9XoY9YH9MmB5gnrB1cr6BHUMw0iT00+/g2eemRladvDBO7DPPltRUVFKXp4FiQlSUVGS0fa7c+rP/rKJmRvYH0JioRrivS639SnDiJann36byy9/qN3xkpIi7rzzZA46aIcs9MqAJIRKRB4A3gJmAm+pasoBBT3MoDMxQUu/TYAPEtTzrQPfz2x3DGPdYvnyBiZNurLd8fz8PP71r/9jr722ykKvDJ9kRlTbAj/y9lVElgNve9tb3uu7XgQFIw1U9SMR+QrnS7UXcF98HREpB3b03j7Rjd0zjF7PhRfex+zZ37Y7ftFFR5hI5QCdCpWqDhWRWmALYPPA63FAIW6k1CIiHxEQL1V9MlGbEY7SehN3AGcBPxWR81R1Vlz5CUAF0ALc2c19M4xey6xZC7j88ofbHd955005+eQDAGhtbeWrr5Ywe/ZSWlpsciieHXYYRWFh5kwQklqjUtXFwDPeBoCIFOIcT4PitRdwGE68Oup1jx2liUgf2n4233KyTET6BY43BteRvDQgfjSJUSFC9FfgaFy8v/+IyOGqOkNEioCjgPO8ejdanD/DiI4zzvg7q1evaXOsuLiQW2/9DfPn1/PXvz7D3Xe/wfz5iZaOjUWLLqC2tjxj7adtTKGqa4iJyx3+cREZihOtjs6NfJTWjbwJjAg5fqq3+UwBJiXbqKouE5H9cGk+xgLTRaQeKMHdE3BTfiel0WfDMEL46KPZ3HPPC+2On3TSAUyd+gmnnvowDQ0593t5nSNyqz9VnQ3MTqJeolHaxsTEazOSH6X1eLwR1MbAH4D9gGHASpyxxRTgVlVtzWIXDaNXcfHFD6Jxgfvq6qqYNauFiy66P0u9MuKR+D9SUieJjAf2AYqBT3AjnndVtbmDc9L2o/JHaar6n3TOX1cZP368Wip6wwhn9uxvGT36GNasafvY2myzTZg5c2WWetUzSWXqT0RmqOr4VNpPeUQlIgfgrNIKgGZio5w1IvIe8Abwhqpel2K7uwCnA9/BjSLeAv4NPJjsKM0wDCNZrr76P+1EqqSkmJkz6wkL2lNUlM9mmw2mrKyom3rYcygoyOxkV8ojKhF5FagFfgCsxmW4vRK3rrI7nr9U/OipoxGViOwNPIKLXTcT5xy8kXed94D9E+RoMjrARlSGEU5zcwvDhx/FvHnxUdFrgD7t6h977A6cf/6+9OtX0R3d69WkM6JKJ9bfGOAmVf0aZyoN8Iiq7okzhngLSKkTwLk4o4mRqrqTqu4I1AEHAlXA857xhWEYRpd54ok3Q0QK3OOmLdddN4Hrr/+JiVQWSUeoWoBF3r5v01kBLiAq8DVwcIptboIzFFjqH1DHw8D2OKu3s9Poq2EYRjtuv/3pkKNlxNtrnXfePhx33Pe6pU9GYtIRqlnEstHWA03AgED5k8ChKba5nAQWfao6F7gB2D/FNg3DMNqxeHE9Dz30akhJ2xHTj3+8BWeeuUf3dMrokHTM05/CmYyfo6oqIm9472/yyiuAQSm2+TzwcxG5OoH59de0TS5oGIaRFnff/TxNTfEGynkEkxMMGFDJ9dcfgojw9tvf8sADnzF9+jfMm7eS1laLTBHP888fRHV1ccbaT0eoLgdmiUipqq7CjXZuE5ErgI+Bk0k9aOoFwCvAP0XkOFWND7q1C85owzAMo0vcdlvYtF8FwQQP119/CPPnN3LooU/x5JNfd1vfeiqZDiuVslCp6jzg2sD7KSKyNXCid2gp8LsU23xTRA7HObV+JSKP4iKIrwa+B+wJTE61r4ZhGEHeeWcWM2Z8GlISS7K9++4bMnu2csgh97JmjfnX5wLp+FFtiBsjv6Webbuq/lpELsFFUnhHVZel2q6q3isi7wJnAj8EfuwVrcJZBf451TYNwzCChBtRFHkbiAiDBo3gN79pH1bJyB7pTP1dBGyKc8xdO97zIqB3KQq6qr4HHCYi+cBor/0vVLWl4zMNwzA6Zs2aZv7xj2dDSmKjqY02Gs4dd5jLZq6RjlBtB9wQZvTgxep7CrheVe9Ot1OeMH2S7vmGYRjxTJ06g2++CZvscaF/ROD9BKvr661Xze9+tzlbb92fkpLuTIzeM6iqymy0jnTueA0JRk6qukZEPgB+BSQtVF56jPNw6T/WAB/hopS/iZtiXJpGPw3DMNYSPu1Xju8Zo1qBPwUY5Je/HMtVV+1kApVF0vGjmo2blkvETFz4o1S4FTgWF+RWcWtUlwJPA4tE5HMv2aJhGEbKLFy4jEceeT2kJOg71T500tlnb80NN+xiIpVl0rn7U4HjROQKVV0UUp6PE5xU2Bm4TVWP8g+IyEhcqo8tgC2BcWn01TAMgzvvfJbm5vil7nyg1Nsv8bYYRx89lsmTt0FEMLJLOkJ1MfAz4HER+amqxtt6Hgh8lmKbS4HXgge8DLizcBHUDcMw0kJVufXWp0JKgr5T1W1Ktt12ANde+/21IjV7diNvv72cBQuazOE3hJ/9bDAlJVlORR9EVWeLyL7Aw8B7IvIY4I+p9wR2wDn9psK9OKfeG1Ltj2EYRke88cZnvPNO2LK6P+2XR3AKsLy8kL//fXdUheuv/4rrrvvaS/1hJOKggwbkllABqOrLXibaPwETcGtK4ALW3gJcnWKT5wAvelEprk+nT4ZhGGGEj6aKiRlOVBFcrr/wwu2YP7+VAw54iQ8+sASKuUA6Dr9bqOpbqvoNcLyInIBz9K0CvkrH2RcXgeIB4BovQsUDuFHaW6q6PI32DMMwWLVqNXfd9VxISWVgPzbtt+mmfWlsrOD733+NNJKfGxkinRHVgyKylaouAfD8qbrk6AtcBxyFs/jbEuerpQAi8gXOTP0NVb2wi9cxDGMd4t//foWlS+NHRYLvO+WC7MRM0tdffzinnWYunLlGOkK1ArhHRPbyQygFEZFi4BZV/VkKbU7ApQf5saquEJFROMHytx2AgwATKsMwkiZ82q+c2FRfbDQ1YkQt//pX4rWoLbesYpNNKigqSserp3dTXJzZe5KOUB2Ms9C7CPhDsEBE6nBGFuNwloHJ0gw8oKorALy0818ADwba7p9GXw3DWEeZNWsBTz89M6TEn/YrIDaygi+/LAmpC9tuW80114xlq62qQ8uNzJOyDKrqJ8Ak4BQR8QPHIiKb4ATsO8DeKTb7GLBZJ9f9JsU2DcNYh7n99qdpP+lTQMzNsxrfPN1NBJUSz4knjuDFF7c1kcoyaY3XVPUh4C/ArSIyVkT2Af6HM4rYVlWfSbHJM4B9RGT7dPpjGIYRpLm5hZtvfjKkpJKY71TV2qOqfQjmowI466z1uPLKjSgosKm+bNPp1J+IPAi8hQuN9LY3LQdwFrA1LsxRP2AaMCFNq78PcU6/z4rIDbgpv+n+VKBhGEYqPPLIa8yZExY4pyLw6j/+CglOAQL88pdDOffc9TPXQSMlklmj2gYXbQJARWQF8A5OvJ4DdgJuA45LkEY+Ge4BNgdqgV8DJ3jX+oxYcNo3VfWJNNs3DGMd4rrrpoYcLSP2yAtO5bUdTW29dTVXXz0WEWHBgmbuvLOeZ59dxcyZTSxZ0kpLi9JiiYfaMGfOSGprs+jwq6pDRaQWF3Nv88DrL3E/RcCtWW0vIm/jBOxtVQ0bdye6xlEAXh6qjYjF+NsC2BU4BGeunrk7YRhGr+CTT+by5JNvhZT4U31FxNajCgj6VFVW5vPPf27OokWtnHHGQv7xj3qamzPaXSMJkrL6U9XFwDPeBqzNPTWWtuK1F3AYaYqKl4fqXW/7R+BaQ732DcMwOuTGGx8POVpALOhszIjCZS2KjaYuvnhDnn22mRNPnEd9vXn85gppx65X1TXA2952h3+8A1FJOwSxqs7GpRcxDMNISGNjE7fdFuY7VYV7BOURG1kF9+G7363h7beLuf56MzDONSI3Z1HV2ar6n5DjvcZ0RkQqRWSyiLwjIitEZJmIvC4ivxeRtFJdeu1pEput8BpGAu677yUWLYp32hViRhTVxB57NQQfgUVFdVx/vQWfzUUsG1iKiMgI4FlgpHeoAeeYMd7bJorIrn6IqTRYAyzuoNxmzA0jAdde+9+Qo7EsvjEjCiFoUDFkSA3TpoVP9dXW5nHccdUccEA5G2xQSGGhkJfnUtcbjpKSzN4ME6oU8Iw9HsGJ1DzgcFV9SkTycGGgbsKFfLoT2CfNy/xPVXfuem8NY93ilVc+5JVXPgop8Y0lKojZf1Xji1deHsyZUxVyHvz855Vcfnk/+vY1O65s0mum47qJScCm3v7BqvoUuMC8qnovcKxXtreI7JqF/hnGOssVVzwScrSIWCQKP9W84Kb9HK2tVQQD04IbLV11VT+mTOlvIpUDmFClxhHe6zRVfTmk/B5cjEKAw7unS4ZhfPXVQu6//6WQEt+IIphqvorYZJLg3DdjiMDNN/fnN7+psTT0OYJN/SWJiJQB3/XehnkToqrqZTw+Htiju/pmGOs611zzH1pa4uMN5BMzougTOF4Tt9/2MXjxxX058sjYVOCCBcqrr8KHH0J9vdLUFFWvew/nnCOUlWVO1E2okmcjYiPQdzuo55cNFJFazwctFTYWkXeB9XAZk+cAzwPXquqbKbZlGL2eFStWJfCd8uP6BUMkVRJbp8qjrWjB4YdX8vvf19DUpNx1F9xyi/Lii5npd2/iD3+AsrLMtW9Tf8kzOLA/p4N6wbLBCWslph9OFH1rwu8ARwMzROT8NNozjF7NHXdMS5Ac0R8VBUMkBUdWfQjGJRg6tICrr67joYdg442VX/zCRCpXMKFKnmDu6oYO6gXLKhPWas8nwGnAhkCJqvbF/QzcE5iB+087U0R+31EjInKMiEwXkekLFy5M4fKG0fNoaWnhyivDjCh8k/QCYoJVQcxoooC28f7g0kv7ccIJ8KMfKZ9+mqEOG2lhU385gqreGXKsCXhCRJ7HTf9tDUwWkZsTRalX1RuBGwHGjx9vMWCMXs2///0qH38cNsHhi1ANsdFU0GiiluDv9H32KWfy5DI++CATvTS6SuRCJSLr4YLIbob7NhR2UF1VtaeYcQdd1juajQ2WReLmrqqNInIG8CTuZ+GuBLIfG8a6iKpywQX3hZSU4EZOecQEq5LYaKqItoFo85g+vS/fdBA5acwY2GknGDxYKC42Z994StvnnIyUSIVKRM7B5anKI7nYfj3pF//cwP4QXH6uMIYkOKerBM3hR0fYrmH0SJ588i3eeOOzkJKawKs/agqOpvrS9vFUyzffhD8Kd90Vzj9f2HZbzFQ9i0QmVCIyETjHezsXeNx77S0hfz4AWnHf/E1IYKLulQHMT8PizzCMJAkfTRXhRlRBi74qYhM7JQQnPYqKSqivb7+UXF4O11wjHH64CVQuEOWI6gTv9WHgEG99pdegqg0i8hKwIy6dySXxdcR9o/f03kad5HG7wP4XCWsZxjrA//73Ac89F+YlUkMs8oRv0Re09OuHP5oSEZqaYu99Bg+Gxx4TNt3UBCpXiFKoNsFN5f2qt4lUgCk4odpFRLZV1VfjyicQm5a7gyQREVHVhNOgIlIM/Nl7uxJ4OvkuG0bv48IL7w85WogbLQVHU9XERlOVxKJTgGoN8aGThg+H558XRoxoL1KrVsHKlaAa2wxHXZ2LmZgpohQqBZarapTrMrnGFOC3uHh/D4jIEar6tBeU9mBcUFqAqaraRkxEZDKxqdFRqjorULyTiPwRuB141su/5Sen3Am4EGfxB3Cuqi6N+HMZRo/hzTc/49FHXw8p8RMi+v5RecTWpgS3NuVTSLyzb//+8NRTMZH66it46il4/XX46CNYujTKT9G7mDYNqqs7r5cuUQrVh8AWIlKsqqsjbDdnUNVmEdkfmIaLoP6UiDTg/iP8n2pvAhNTbFpwlny7AojIKtzIKfhzsBW4SFUv7spnMIyeztln3xVytABnEBu09As69NbS9nFXR3DKLy8P7r1X2GAD4bXX4JZbnEAZuUGUQnUzcANu+usfndTtsajqLBHZDDgFOAgYhcsh9R5wN/C3NKY+3/Ha2x43WuuH+7nXALwPvADcqKrvRPEZDKOn8sorHyY5mgo69MaPnqoITgEC/OlPTqROOAFeDgs3bWSVyIRKVW8SkX2Bq0TkK1V9Pqq2cw1VrcdN453TWd3AOZOByQnKFgGXRtE3w+jNnHVWO7943GOsEidQNd6xvsRM04MGEwXER0vfYw/nJ3XIIW4dysg9ojRPPxt4G2dsMM2zkHuVTpxeVfXcqPpgGEbvZdq0mTz99NshJTXE0nXk4UJk+lHTy4kFpAU35Rdb9R88GMaPFy64oHMLv/x8tw7jGw2ImOOvT6bvQ5RTf5OJOfAK8D1iaTE6woTKMIwOUVX++MdEoyk/c68/1eePoMTb96kEYiEU8vNh3Dhh6tTET9kRI+CAA2CHHWD0aCiwoHNZIcrb/jw9K9KEYRg9hKlTZ/DSS2GB+PzI6L44VRATo2AEt3zaWv3BVlsJc+aEi9TAgfDb38Luu2fW7NpIjijXqHaOqi3DMAyf5uYWTjnltpASP89UifcaHEEV0daAou2U34gR0NQUPmW1//5wyilQUdG+zMgONpA1DCOnuemmx/ngg69DSvzRlG9qHjRB70/MgKKKYNik6mqoqpJ2oZHy851AHXKIrT3lGiZUhmHkLMuWrUzgN1WME58K3IgqaIJeQ8z8vJCglV9BAQwcKBQUtFWiggK45BL4/vdjx1pbYc4cmD0b5s6FFStg9Wpo7i3RSyPkyCOhpKTzeuliQmUYRs5y4YX38+23y0NKfPPzOu+9P6pqK0xuZBWb8hs0SCgraytSInDeeU6kVOGTT+DVV+Htt2F52KWNdrS0ZLb9tIRKRG71duep6plxx1JBVfWodPpgGEbv5pNP5nL55Q+FlJTjRlR9iflQ+VN7QWGq9ep572qhX7/2c3pnnul8qaZPh8cfhy+/jO4zGNGQ7ohqEs7C7yPgzLhjyeahEu/VhMowjDaoKr/+9Q00NcXPs/lrUb6xRD4xA4oaYhZ/pQRTzZeWwrBh7R9NJ50E48bBX/4Cn38e5ScwoiRdoboDJzLzQo4ZhmF0ifvvf4knnngzpKQK99jyp/rqcGIVnPLLJ2hMUVAAo0cLeXltheroo6GqCi680CKh5zppCZWqTkrmmGEYRqrU1zfwu9/dHFJSgBs1VeCm+sqJRaAYSGzKbwB+MFoRWG89oaiorUhNmODM0//zn477UlAAG2wAw4a56OolJe6YWQW2pbi48zpdwYwpDMPIKf74xzuZOzcsObZvfu77RPmGFP2IrUXVEgw4O2JEe+OJvfZyRhJfdJB+dORIl4Z+880za81mJIcJlWEYOcMLL7zHVVc9GlLij6DqcI+t/t5rGbG1qOC+s/Dr06etSG27rcsrlcjEfNAg50c1dqyNmnIJEyrDMHKChobVHHnkVbRPdu0bUJThLPwqcVN++bhpPsEZV8TWpWprYcCAtq1suKETqDABEoF994V99rF4frmI/UkMw8gJzjjjDj79dF5ISR/c1F5/nNGEP+U3kFgm3wH4a1TV1TB8eNvIE0OHQnl5uEiVlcFRR8Gmm4b3q7nZOfs2NTkn4NbW9D5fb2bgwJ6Tit4wDCMtpk2bmWDKrxhn6VeHEylfkPoRM0UfgB98tqICRo5sK1IDBrgtTKSGDoXjjnOGEuBEaP58F4li4UL49lvLUZUMhx9ukSkMw+jFfPPNUiZOvCzBlF8/YtN9tThxqiS2FlWHL1ilpe3N0OvqnBiFidS228LPfgZFRbBgAXzwAcya5UZORm5hQmUYRtZobW1l0qQrmTcvzMqvD25dqj/OkMKfAvR9qPrgRMv9ml9vPSE/v+1IKkyk8vKcwcT3v+/CJb37LiwOu7yRM5hQGYaRNS677CGmTp0RUlKCm/Ib6O0PwD2ufH+pKvwgtCUlsP76QmFh2zWpsOm+mhr45S+dWN1/PyxbFvUnMjKBCZVhGFnh6aff5vTTp4SU+D5SdbgR1UDco2oQsRxUfQFpJ1Iizgeqb9/2rY4dC3vuCTNnwqJFyfdTxBlclJa6VCCWSLE9mb4nJlSGYXQ7n346lwkT/kJLS5gJXR1utFSDG0kVExtZleJP/cWLVH4+rLeeC4sURAR22cUJzQsvdN63wkIYMsRFo6irc6buJk7ZxYTKMIxuZdmyley///ksWbIipNSf0huAM6SoILZGVYpv9Vde7gwn/LxSxcWw/vpOjNq0VgWbbeYiUXSWsmPAANhoIxg92nypco2upvnoKpbmwzDWIRobmzj44Av54IPZIaUlOHEaRGxE1RcnXjGRqqqCUaNi1n1VVeHiMnQo9OsHjY0d92nIEBg/vr2DsJE7dDXNR7pBRizNh2GsYzQ3t3DooZfw9NMzQ0oLcEI0BGfNV4cTqT4ERaq2tq0z78CBTmiCRhOFhW7arraWDunfH7be2p1v5DZdTfNhGIbRKa2trRx99N/4979fDSkVnBANxvlKDfBe++Cm/tya1ODBQv/+ICLk5zujiT592rZUXQ0jRjjfqET06eMEasSItgKn6nyoVq1yr2vWuK25ORaRwqJShLPhhm6NMFNElubDMAwjjObmFiZNupI773w2QY3+OJHyp/364USqBqglL08YOVKornaqUlHhpvqCYuSPovr0SRxMtqLCTfGtv74zjmhpcaGR/G3VKhOidMl0Pi9bMjQMI2M0NjZxyCEX88gjryWoUQcMwwnUYJxo1eCm/aopLnbrUaWlTn0GD3YRzoNiVFfnpu8SGUBUVLh0HRtsAKtXuxBJ9fWwcmU0n9HIPCZUhmFkhPnzl7DnnpOZOTNR4qe+wHBgKM78fBAuNFJ/oIzaWhg61EWbKC11U3UVFbGzy8vdKCp4LEh1NYwZ4wwqVq6E996zTL49FROqNBCRSuD3wMHAKKAF+Bi4B/ibqqYdLUxEBgCnAfvh/otXAe8BU4BbtH1ANMPIOV5++UP22+98Fi9OZBPeB/evMxwnUENw1n0DyMsrZNgwobZWEHEjqGB07uJiN4JKNM1XXe0s/kpL3fTeggUZ+IBGtxK5UInI5sAJwPdwP5XKO6iuqtqjxFJERgDPAiO9Qw04j8Tx3jZRRHZV1SVptL0V8DjupybAClwws+952wQR2V9VV3flMxhGpmhubuHMM//BJZc8GBJk1qcvsAEwAjeS8o0o+lJRkcfw4UJxsVBR4UZRvm9UQYETrP792zvgijjh6t8/cTqPZMjLc9crKXHrXkVF7jUvz20illAxjEwaUkDEQiUivwYuwyWJ6XV/ThHJBx7BidQ84HBVfUpE8oAJwE3AlsCdwD4ptl0NPIr7L/4Q+LmqTheRIuCXwOXAHt7rryL5QIYRIdOmvcNPf3oZ33zTUXyiOmBD3GhqqPd+AHl55QwZIvTtC8XFwrBhLi6fSEyg6uraPxALC93UXl1dx5Z+iSgshMpKt1VUuDZMiHKPyIRKRLYFrvTeXgv8B/gvsBg4BPfTaTfgMGA5cCLuYd+TmAT46dUOVtWXAVS1FbjXE6y7gL29UdXTKbR9Cu4erQL2UdUvvLabgGtEpAq4ADhGRK5Q1Y8j+USG0UVeeulDJk26lk8/ndVBrTzc9N4Y3O+8YTgz9L7U1hYweLAzmAiOmIqKnBNuv35tBUrEiVi/fs7ZNxVhKShoK0zFxSZMPYEoR1Qn4kZRV6jqyYDvlNekqs94de4Skatw01vnAeMivH53cIT3Os0XqTjuAf6M+7l4OJCKUB3ut+GLVBx/A87AOZZMBM5JoW3DiJTPP1/AH/7wIA8//BJNTZ3EJqIY9y+xEbAebqpvIOXlZQwdKlRWxgQqP98JyIABsRGVT0WFc+KtrU0+xFFBgTvPF6aSEhOmnkiUQvVdnBPwlXHH23wtVPUtEfkN7qF+KnB2hH3IGCJShvuMAFPD6qiqishjwPG4abpk294Qt6rcUdsrROQFYG+v7ciE6gcb3cYLn3wa1rMEZ3RkzxFepu3eJGpDU/AkT8WuJFUblAT1/cMSViXZa7Sv1/GZUX/Ozuok+BuqAvW4pdNkl0n7AGNxI6nRwBBKS6sYNCifujoYMEDo18+NbGpr3SiprCx2ti9ONTXJTe0VF7s1qvJyE6beRJRCNQBYrapfBo614gJ4xfMvYA1wED1EqHA/B/0l3Hc7qOeXDRSRWlVNJiXbJiHnJ2p7b9x/fmS89NlXNLe8E2WTvR+zveyEIpyxxBa4f53hVFT0YeDAAgYPdgJVU+MEyBeivDw3Aqqqclt1tVtDSkRhoROi8nInbmVl7UdaqmaS3l1kMsJ8lELVENJePVAlIsVBSzVVXSMiDbhvck9hcGB/Tgf1gmWDcWt0UbddJSIVqhoWftowskg+7jfr1sBmiIyiX78ahg7NZ+hQoa7OGT7U1DghKipyQlNZ6cSprKz9CEjEjZRKS9tuFuF83SHKP/UcYGMRqVJVf9L6M5wV3NbAi35FERmM8+xriPD6maYysN9Rv4NllQlrRdN2O6ESkWOAYwCGDx8eX2wYGaIQZwu0HbAV1dUjGD68lFGjhKFDXYw+f/SUaGqusNAJkr+VlLhXs8QzohSqmcDGONvT171jz+IMJs72/H8aPXPrq7xym2+KGFW9EbgRYPz48TbpYWSQItzvzTHk53+fQYPGMmZMGWPGuNh8dXXOtylozFBY6EZCvn9SUVFs38TISESUQvUozvT8J8SE6hqc8++uwGwR+Qj4Ds67T4GrI7x+pqkP7JclrNW2rD5hrY7bTmRGlU7bnTK4zwBmfdM/jTNTfbJE+SRKt60o+pBKG9m4R5m7pkg5RUVDqKvbmO2225pddilno43yGDTICZG/TpSfH3u17LhGV4lSqP4L/AmY7x9Q1S9E5DDgNpw4be8VtQKXqOqdEV4/08wN7A/BjSDDCGa3mZugTmdtJxIqv+3lUa5PfbHguKiaMtYBVi1dyszzzqPxth8z+KEV9K+poayujvzBg8kbOBD69nULUNXVTr18T9rgcCq4ZWoolZcXWwCrrAxfADN6BJEJlbcu9aeQ4/8SkedwkRqGAcuAJ1Q1zB46l/kAJ7B5OCu9UDNyYhZ885O0+IO2ln6beNfqqO33k2zXMCKntKaGbS+9FC69lObGRl4+91xm//WvDH//fYYBtUVFlFZWkj9woAspUV0d884tK3PCVV4es4bwF6eKitzCVHDrynCstdWFSa/3Jh/y8mLCacLVoxCLcZo8IvI8sCPwjKruGlIuwKc4h5E7VPWI+DodtP0lzpfqVlVtl/VYRMpxo9UK4FxV7dSPavz48Tp9+vRku2AYadPa3MyL55/PJ+edx9DW1rUJO6pEKKmuRmprnRdvnz7OqqK83AmRLxpVVeG26L54lZU5U7+ysuico3zvYn/EZU5X3YKIzFDV8SmdY0KVPCJyFHAzbn1te1V9Na78EOBe7+1uqYRQEpHzgLNwln0bq+qsuPLTgL/gIrWPTSaEkgmV0d2sWLCAhw4+mDUvvcQQnKF6H28rKyggr29fN7ryPXn79o1FnS0tdaOvqionHIlEQyQmWr4JYWlp10UmPz9m+x6MShtcbDMh6zImVBlGRAqAN3Dx/uYAR6jq016Mv4NxIlYFTFXVfeLOnUwsmsSoECGqxgWjHYib2jtcVWd4VpJHAVfgzKyuU9WkgtKaUBnZQFV5/ZpreOW3v2VgaysDcDl7a3D/HOUi5NXUuBFWaakTmn79nHD5Qf3y82P27JWVnU8BBtej+vd3ApOpdL0mVu3ZdNOkHdvSEaq01qhE5FZvd56qnhl3LBU0bJorV1HVZhHZH5iGi6z5lOe4nEcsAsebuFh8qba9TET2w8VBHAtMF5F6r11/TuQJ4KQufQjDyDAiwja//jVDttmG+/bck4alS2kAVuKmCypVqVqyhNIlS5CaGreOtXIlzJ7tvIF9oVm0yG2+aPXt60ZbYQTXo+Z6tkkbbAAjR7qyFStccqoosB/33U5aIyoRacVNf32kqmPjjiXzc8Ovp6qa4Uwm0eMlTjwFFwJqFM7I4mPgbhIkTuxsRBWoNwD4Ay5x4jCgEWdsMQW3fpX0z0QbURnZZvns2dy5667Ixx+vHVlV46YCS3CjrEI/HPrgwTErwD59XMZEf1rQp6jIjb5cPpDkOjF8OGyxhRM5X8xWrMjciGtdZLPNMjqiSleobseJzTxVPSPuWEqo6i9S7oCRFCZURi6wavFi/rHrrjS99Rb9ia1ZVeOyqpbjwqxIXp4bTflh1MEJ1uDB7QULnPD4IS+SYdgw2GYbJ3Kq0NDQVrhspJQ+uShURs/AhMrIFVbX1/OP3Xdn5auvUocTqlqcSNXgFl9r8EZXRUUul3xwmq+21h0LC6FeXOwEKz5xVSLWXx/Gj2/bviqsXg2rVrmtqQnWrHFbc7MbfdkILDEmVEa6mFAZucSqJUuYstNOrHn3XfrihKkWKPVeC4gZWwBuNDVkSOwBmJfnDDAGDgwXpLw8t8Y1eHDngiUCG20E48a1zSvSERaKPTEiSRuZZFWoRGQ40KKqHUX/DtYfDBSo6leRdMBohwmVkWusmD+fW3fYAfniC/rgxKovbr2qDy7FYgluWjDPz0M/bJgzW/cpLHSjq759wy8i4kZYI0d2bkCRnw8bb+zWsErCMhIZUZNtoWrFrVkN6bSyq/8FMExVLVh/hjChMnKRJZ9/zk3jx1O2ZMnatapanEj1wY2wCrz9Av9Xer9+bqQUNFOvqnKGEokExh819evn1qM6orAQvvMdV7+2tmsf0OiQdIQq6nCR2YxSahhGD6DP6NEc8uCDLMvPZxkusOUiXM7gJTgz9mbgW2C1/0P622/hk0+gsTHW0PLl8N57zhw9bP1IFd5/H6ZPd6OrMIMMnzVrXFv33w//+he88QYsXBidSbvRJaIeUc1X1cGdVnb15wG1qpqkjamRKjaiMnKZ6ddfz+PHH08tbvRUQWxkVUksSVs1UOaPrPLy3Cgq3tKvrAxGjepYjPbYw9WZOdMJUzL4Plx+iKW8PHfMnH7bMn587jn8dhVvfaoO+CYb1zcMI/uMP+445r3xBu/edBOCyw28GLdm5eewqcRFsW5RpVLEjZxmzYoZVfiC0dDgRk9DhriyMCF54gnnBPzzn8OXX8K773Y+YmppiTkeG4nZYouMplxOu2UR2QnYOe5whYic3dFpuPXTfbz9VzuoaxhGL2evK67g6//9jxXvvUceMbHqR1uxWgG0qlLtC9CCBc6MfMSImIWfqotusWyZGzmFmbJ/8glccgkccwwceqib4vvoI5viy3HSnvoTkXNwkRb8BoTkHX4FF3FhZ1V9La0OGJ1iU39GT2Dh++9z41ZbUdHYSDlOpIq813ycyXqFV7cUNxUovmCVlMDo0e1FKT/fiVUiZ+C8PDj4YNhtN+cz9ckn8OGHsDjZzDxGGw4/PGmrye6e+psFPBd4/31gDfByB+e04tZO3wWm9MCcVIZhREzd2LHsffXVPHr00RTiDCr6EpsGXI6z+ioDVnnnVKs6sWpsdCIzenTb9amWFvj0Uzc9OHRo+4u2tsJ998EXX7iH7CabuG3ZMjctOGeOM6YIGm8YWSNrxhRG5rERldFTUFXuOeAAvnjkEWpx031VOJ8q31i8lljk53Yjq7w8N4KqrKQdFRUuunciA4rBg+H4493aVttOuWC5K1a4rakpFqHCHH/bsvHGPSMyhYgcAaxS1X9G0qDRZUyojJ5E/bx5XLvxxhQsWUIFMZ8q3wJQcCMsf5KvFKgJGk2IuDWrsOm+ggLYccdYtt94SkrgyCOdUYCRUbLqR6WqU0ykDMNIl8pBg9j7qqtYATQBS3H+VPW4BW3FTQc2e/VXAcuDP7RVnUXgkiXtG29uhmefdVaBYTQ2wrXXwk032TpVDhK1w69hGEbabDpxIhvuvz/LcAvaS3ECtQQnUK04sfLde1cCK+Jnhb78MtycXBUeesitRYVZBAK8/jr88Y/wj3/AvHkRfCIjCiIPSisi6wGHAJvhppULO6iuqrprpB0w1mJTf0ZPZNnXX3PNRhtRuHIlVcSm/gpxloBC27UriHMK9hkyxAWpDePkk+GDD5zBREcMGeLWt4YNc/EDS0oy6i/UY+nTp/MszB5Zd/j1TNbPwo3Ukk2gaBiGsZbqYcPYefJknjz1VEpwU39++JrlOFFq9Pb9RB3LgDxVSoJiNWeOe3iGBa+94gr405+c1d/MmYk7M2eO24yOufxyKC/PWPORCZWITCSWwXYuLqX6XGJTyoZhGEmx7W9/y1u3387i996jH27qrz9uqs+PsL4CZ1jhWwIuBfqqupxWPl9/7Yws4gPNtrY6obrsMmdA8a9/JTa0MLJOlCOqE7zXh4FDwtKxG4ZhJEN+YSH7Xncdt++0EytxDr/1uBHUUlz8tXycgNXhHmS+sUU/VfKDYvXVV25kFW8N2NwMp54KV18N550Hjz3mDC7MdyrniNKYYhPcd+VXJlKGYXSVETvuyKaHHcYK3LSMbw3oG1lAzNDCX0No9d+HGVgsX97+IqtXw+9+56wFDzoI/vIXOOQQZ+Zu5AxR+lEtA1pVtU8kDRpdxowpjJ7Osq++4uoNNyS/sdHlp8KNoAS3VuWvipThgoj6lAB94o0rRGDMGJe6Pp6aGmeavt56sWOLF7tAt/46le/0m2zk9XWJc85JOlNyth1+XwW2AKpUdXUkjRpdwoTK6A08feaZvHjBBdTgBCjoAOxP+0HMQdinEqiIF6u8PNh88/DIEjU1LmDtVltF+wGMNmQ7ceLNOAvSCRG2aRjGOs73Tj+d8v79qcdN8dXjpgKV2BQgeOlAAu/rCSRe9GltdVZ+Yb/+ly51oZTuvddCJOUYUUamuAlnSHGVlwLEMAyjyxRXVrLzuefSgrP6AydK4Nas/GPBtSufpbhcVm1oaYHXXnO+P/E0N7t1ql/9Cj7/PIruGxEQ5dTf2ThDnF/jpotfwuWb6tDmU1XPjaQDRjts6s/oLbSsWcO1Y8ey5NNP16b/8IPUCs503ctK1SYtCLhpnr4EAtj61Na6ab5ETr95ebD77nDggTBuHBR2FLvASJZsr1G10jY3FSTh0Kuq+Z3VMdLDhMroTbxz1108OHHi2sjp+TiBCotUUUfbkDhlEEu6GGTLLWGzzeDttzu+eEkJjBzpxM2PwJCXZynpfc4/P2mH32xHpngeizRhGEaG2PgnP+GFCy5g4XvvrbX2W4Ezmmj0tqDzrx9uCaABKI6PXAHw5puwzTYuieIDDyS+eGOjS6xohNOc2bgOkQmVqu4cVVuGYRjx5OXns8u55/LPgw+mHmfltwJn6VeAW7cqwi28r4G1jsI+S3HOwAXxYnXDDXDPPfC3v8Gf/wzz52f8sxipYdHTDcPoMYz50Y8YNG4cq3GGFL4VIDiLvxWBur51oI9vJRg67fPLX7qgsw8+6ByAwwwtjKyRMaESRz8RGZ6paxiGsW4hIuxy/vlATKBW4UQLnFD57rjx5ut4ZcvD1uXr62HCBBep4vDDYepUuOgi2HVXqK6O9kMYKRN5vHoRGYeLoL4bznFcg9cRkT7ARd7xE3tauCURqQR+DxwMjML9kPsYuAf4W7qfR0QmEwvq2xEbqOqn6VzDMHoD6++1F4O22op5M2awGhekdjluTQpv34+X7puvB5f5G4AiVUrjpwBnzoT99nMx/8rLYY893NbaCgsWuPxUK1c6Hyt/MxxJRqVIl6jTfPycmONvKKq6RERGAbsCjwD/ibIPmURERgDPAiO9Qw24/5Px3jZRRHZV1ZAUo0mzBhdbMxEWjd5YpxERdjzjDP558MGswP0DNhEzplhNW8OK5V6d4MNuGVCYl0dBaytteOklJ04PPxxLD5KXB4MGuc3ICpFN/YnIRsBNOJG6Cvfg/jZB9TtwBjkHRHX9TCMi+ThhHQnMA3ZX1XKc5etPcTMRWwJ3dvFS/1PVgR1ss7rYvmH0eMYceCD9xoxhDU6UwAmSP8ZZFthXYg7CBI4taWlB80O8Y15+GbbfHmbMiLzfRnpEuUZ1Ms7o5hpV/Z2qvkHbiCZBnvFet4/w+plmErCpt3+wqj4FoKqtqnovcKxXtreIWNZiw8ggkpfHd08/HYgZUDTj1quANlEswI2yVtGWZmBxczMa5sj76aewww7whz+40EpGVonS4fcz3GhjhKrO9o7NA/qHOfWKyAqgRVV7xEqliDwP7AhMU9UfhJQL8Blu3eoOVT0ixfYn49aonovK1N8cfo3eTMuaNfxt/fVZ9tVXawPWBp2A4yNW5Hnv43+dlw4YQPXSpUhTguXlsjL40Y/ge99zqelLSszRN56ddko6cke2HX4HAyt9kUqCVTgH85xHRMqA73pvp4bVUVUVkceA44E9uqtvhrGukl9YyA6nnsrU3/yGFTihasEtHPtWXPXE0n+04qYHa+LaWbVgAYU//Snljz0Gy+InCYGGBrjzTrcZ4Xz7bfssyhES5dTfaqBI2gXUao+IlOK+LyHfipxkI2L36t0O6vllA0Uk3b/axiLyroisEpEVIvKRiNwkIlum2Z5h9Fq2POooyurqaMY9gIC1UdbBiVbQ+qghUC/I8vvuY82tt7qQSkbOEaVQzcIZUmyQRN19cCPy9yO8fiYZHNif00G9YNnghLU6ph9OGH2Lwu8ARwMzROT8zk4WkWNEZLqITF+YKNimYfQSCktLGX/88UBsraqVtutT8Xl9lxLi9NvSwtLzz0effx7OPNNN7xk5Q5RC9RhuWvi3HVUSkb7AxbjvSk8xTa8M7Dd0UC9YVpmwVjifAKcBGwIlqtoXN4OxJzADd2/PFJHfd9SIqt6oquNVdXxdXV2KXTCMnsfWxx9PflERa2jr+OsbnjcGjoObHgxL6dD8xhusvPJKOO88+PhjOOEEl0zRyDpRGlMMwD1sy4FzgctwjrD9VTXfm+77EfBnYATOdH19VY3/wRNVfyYBt3Whib1V9TGvrcOImZ0ndLgVkd2BJ7y3O6jqy124frDdElzQ361x/4NDVbXTaVMzpjDWFf49aRJvT5lCMS4GILRN91FEzCHYJz7COgAi9Hn4YUr228+9b2yEadPglVeceC1f7qJXGG158EGoqkqqalaNKVR1gfdAvx84G/g/vBGbiLwHjMZ9XwQ3TXxopkQqAwR/gHXkgh0s6zAPVyqoaqOInAE8ifvf2xV4MKr2DaOns93vfsfbU6awGucxX4j7RVeOe+AEHYJ9ltI2wjoAqiw99FBqH3+coh12cFOAe+/tNiNrRBrrT1UfBXbCTVUV4YRQcGsuxd7+m8BOqvp0lNcO4W7cj6Z0t2D/5gb2h3RwzWDZ3IS10iM4OhsdcduG0aMZuMUWjNx5Z6Btxt/gXHz8L8c1tA1i66MrVrB4zz1Z/dxzUXfTSJPIY/2p6mvANiKyGfA9nFFBPjAfeElVu2UuSlVXE27gkw4f4L73ecAmJDBR98oA5qtqR2GQDMOImO1OOolZzz5LI24dKh8nRGW4X8hrcD4xpYFz6oGiggKK4/Ip6YoVLN5tN6ouuYSy3/wGCYtgYXQbGYuerqozVfVaVT1LVf9PVa/sLpGKGlVtAF7y3u4VVsczy9/Te/tEWJ0usl1g/4sMtG8YPZoN9t2XPuutB8RGUi20jUgRDLPks7i5GQ0zPGpuZvlJJ7Fou+1Y9cADaGNj+zpGtxDZiMpL59Giqh2ZbwfrDwYKVPWrqPqQYabgIlPsIiLbquqrceUTiE3J3ZFKwyIi2oFVi4gU44xQwM1sZHra1DB6HHn5+Wz9q1/xxO9/v9bpNw83airFjaqCDsE+CnyzbBmDx4yhJSSL75rp01n64x9Dfj75w4a5CAx5eUhexn7n9zj6/u9/5GXQQjLKqb9ZuGCtHa3hBHkJGBZxHzLJFJzp/abAAyJyhKo+LSJ5uJQfN3n1poatv8Wl8RgVF1x2JxH5I3A78GwgBFUhbs3vQpzFH8C5qro0ws9lGL2GLSZN4pkzz6S5sZFVOEFqwRlS+FN+vnAFZaalqYkVY8dS1bcva156iVBaWmiZNStjfe/RxEehj5iofxKkGgCrxwTMUtVmYH+cIA8BnhKRlbgRzj9x1rBvAhPTaF5wlnx/B74WkQYRWei1/RROpFqBC1T14i5+FMPotZTW1rLxT34CJDakiHcI9lny4IMUnHYaZSeckMEeGumQzbFrCT0st5I3CtoM5yf2Lm7WYA3OyvEUYLs0c1G9453/AM73bBUuxNQq4G3gamALVT2za5/AMHo/4487DoiNpMA9aIIrTEGH4CBfn3giFRdcQO1zz1Gw+eYZ7aeRPFE6/LbirN06DR3krU99BXyTTH0jPczh11gXUVVuHDeO+W+9RRHgB92Md/otJeYcHKTPhAmMuvdeUKXpmWdouO02mp57jtY5SS2/r5MMWLSIvCSD0narw6+I7ATsHHe4QkTO7ug03EhhH28/3iDBMAyjS4gI448/nkePPZYmYg7ATd5W5NVbhTNdL447f8l991H+3e8y4Le/pXi33SjebTcAWhctonXJEmhpQVsSpdpbN5HqzCbCSHtEJSLn4IwD/AaEkFiPiU7HjcR39vyujAxgIypjXaVpxQouHTyYpvp6SonlEyoG+gbq5eO8+9utgYgw6p57qD3kkMx3dh2ju0MozQKCrtvfx/146Si+nZ8S5l1gSqKYeYZhGF2hqKKCzQ8/nNevuYZGXIToPGgTYgncOtZSYtODa1Hli4kTaW1ooN+kSd3TaSMhWVmjMroHG1EZ6zLz3nyTG8eNA5xJrh+IM2xtqpq2vlVB+h17LEMuvJCCPmErWkaqZDvD7y9o6wRuGIaRNQZtuSUDt9iC+W+9RQMxoVqFG2EFH37LgPIhQyDEYOLbG25gyT330GfCBEo23pi84vhVLaPvL35BXgZzeEU2ojJyDxtRGes6r/7tbzx24omAm97zDSnKaJ+SXvLyGD1+PKtes2XzVNl80SIKMmj1lxE/KhEZJyKXiMg0EXnPS60+TUQutpTqhmF0F5sedhj5RU6egg7ADbj1qSDa2srymhr6TEzHZ9/IJJEKlYiUi8hdwOvAyTgDi42Asd7+74HpInKniCSaEjYMw4iEsr59GXPggYAzMw46+YZFp1j4xBOUTpjAsKuuIq+iIqSGkQ0iEyov5t1DwE9w5ufzcVlx/+Jtd+JiAQrwU+DfXsRxwzCMjLHFkUeu3Q+OqlYSHhrn3RNOoOZnP2PjDz6g39FHI6WlIbWM7iRKY4rDgR/grD9/D1yrqm2ilHhidhxwuVf356QYadwwDCMVRu+2G1VDh7J89uy1gWp9p8/ltDdNb5wzh3eOO45x99zDiJtuYuill1I/bRoNb71F8zffoBkOwNoTkQwbmEQpVD/D/e1PVdWrwyp4wnWtiOQDV+LEzYTKMIyMkZefz+aTJvHC+efTgotO4T9WG3G+VfGP2bn//Cd9dtiB0b/9LflVVdQccAA1BxzQnd02AkS5RrU5bn3yps4qAjfjRt1bRHh9wzCMULYIOO02xJUtJTxA7XsnncTcBx7IXKeMpIlSqCqBelXt1JfKq1MP2GqlYRgZp3a99Ri+446AG0EFLf5acH5U7VBlxiGHMOuGGzLfQaNDohSqb4FqEenfWUWvTg2wKMLrG4ZhJGSLX/xi7X78qGoV0JCf3/6k1lbeOe44XjvwQBa//DKta9ZktI9GOFGGULoPOAi4QVV/1Und64BjgftV1aI+Zghz+DWMGKvr67l04EDWNDQghAejHbneejR99lmH7UiYoOUSWTCm3mPBAop6iMPvNThjmmNF5O8isn58BRFZX0T+gRMp9c4xDMPIOMWVlYydMAFwD5/4URXAnDlz6LPTTh22o16aj5zdmpu7fcs0kQmVqj4LXIETq8OAj0Rkloi8JCIvisiXwEfAod4pl6vqc6GNGYZhZIB4o4r4+aQ1jY2sGT2aYYFpQiP7RBqZQlVPBk4EluAEaziwPbADMMw7thj4jaqeEuW1DcMwOmPETjtRM2oU4Cz9wkZVH0+ZwuBf/Ypxd99N8aBB3do/I5zIY/15PlRDgQOA84AbvO0879gwVbUpP8Mwuh3Jy2PzI45Y+34FIabpqjx35JEMPOggdv38cza/+Wb677svRf07tRMzMoRFT+/FmDGFYbRn6axZXOmNqsBFUq8KqbflGWewzZ//3OaYtrSQ08/MLPVNCgpINiJe1vJReaGRxuD+3otV9eMo2jUMw4iampEjGbnLLsyaNg1w03+lxLL++rx5wQX0GzeO0QcfvPaY5OdjAUq7ny5N/YlIoYj8Bbfu9A7wEvCBiCwUkTMt6KxhGLnIFnHGEstob1gBMO3nP2fu8893S5+MxHR1jerfwCm4kZQEtr7AucBtXWzfMAwjcjY66CCKKivXvm/GrVfF07xqFf/dc0/ev/FGC0abRdJeoxKRCcC93ttPgfuA2cBIYCIwGPcj5Qdmhp4dbI3KMBLz8NFH8+Ytt7Q5VldWRn5DmC0g1IwZw+Cddybf0n60Y+vzz6ewrCypuumsUXVFqB4Cfgg8ARygqqsDZVXAM8CWwC2qekxaFzG6hAmVYSRm/ttvc8MWW7Q5JsCoESNo+PLLrPSpp3LEokWU5GhkinG4EdNJQZECUNXlwB9wf3dLPW8YRs4xcPPN2WCffdocU2B5URF141N6jhoZpitC1Q9oVNUPEpRPD9QzDMPIOb57+untjn37yScMnjiRDQOZgY3s0hWhKiZBdHwAVV0WqGcYhpFzjNhxR0aExPZ74c9/ZrvLLmO/Z55hUCex/4zME2WG316NiJQB3we2wk17boULEQXwJ1WdHNF1BgCnAft57a8C3gOm4Nb7ctjb0DB6Hrtfcgk3b7ttm2MN337L4yefzAG33MKQXXZhxezZfPPqq6ycPRttaUnQ0rpLQYYNTEyokmcb4L+ZvICIbAU8jjPvB2cxWwl8z9smiMj+8WuChmGkz5BttmHTiRN558472xx/69Zb2fCHP2TMgQdSMXQoFUOHZqmHRlf9qAaISEuiDbc22VGdzMeHj5YlwNPAJbgo8POjalhEqoFHcSL1IbC1qlYC5cCvgTXAHsDlUV3TMAzHrhdeSGF5ebvjD/3iF3z74YdZ6JERpKtCJRFsPYUXVLVWVXdT1dNU9R5cVuuoOAUYiJvq20dVpwOoapMXxPccr94xIvKdCK9rGOs81cOGsdcVV7Q73rh0KXfsthsL33+/+ztlrKUrU39/iqwXPQBVzfTE9OHe6z2q+kVI+d+AM4AKnEP1OSF1DMNIky2POoqPHn6Yjx95pM3x+jlzuGnrrdnxrLP4zr77UlhejuRFnniiR1M9fDh5Gcx8bNHTu4CIzAJG0EVjChHZEDfdB3CIqt6XoN5/gb2BV1R1+87aNYdfw0iNhkWLuGW77Vj86afZ7kqP4rRFiyjNUYdfIzo2Cey/20E9v2xsBvtiGOssZX37MnHqVCosYWJOYUKVGwwO7M/poJ5fViUiFWEVROQYEZkuItMXLlwYWQcNY12hdv31OfLFF6ndYINsd8XwMKHKDSoD++ERMduXVYZVUNUbVXW8qo6vq6uLpHOGsa7RZ/Rojn3jDbY69lhbj8oBeq0flYhMomtpRvZW1cci6o5hGD2MoooK9rv+enY49VTenjKFL59/npULFtC8enXWMunmKpkW814rVD2M+sB+GbA8Qb1gHP36BHUMw4iQ2vXWY5dzz812N9ZperNQ3Y1zoE2XhHEMM8DcwP4QEgvVEO91uaqG5XkzDMPodfRaofLCDPWUUENBS79NgEQR6X3rQPM+NAxjncFWCXMAVf0I+Mp7u1dYHREpB3b03j7RHf0yDMPIBUyococ7vNefisjIkPITcFEpWoA7Q8oNwzB6JSZUKSAifUSkn78Ru39lweNhPk4iMllE1NtGhjT/V1yQ2zLgP14kdUSkSESOB87z6t2oqh9H/uEMwzByFBOq1HgTWBjYhnnHT407fnWqDXuJJvcDFuEiT0wXkeW4VB/XAkW4Kb+TuvYRDMMweha91piiJ6KqM0RkY+APONEaBqzEGVtMAW5V1dZk25sxY8a3IvJlit3oB3yb4jm9FbsXMexexLB7ESOdezEi1YtYUFqjDSIyPdWAkb0Vuxcx7F7EsHsRo7vuhU39GYZhGDmNCZVhGIaR05hQGfHcmO0O5BB2L2LYvYhh9yJGt9wLW6MyDMMwchobURmGYRg5jQmVYRiGkdOYUBmGYRg5jQlVL0REKr2QTe+IyAoRWSYir4vI70WkqIttDxCRS0XkIxFZJSKLReQFETlaRCSqzxAVmbgXIlIjIgeIyLki8qiIzAuEx5oU8UeIjAzdiyEi8isRuU9EPvW+E6tE5AsRuVtEfhD15+gqGboP3xeRP4vI4yLyiYgsEZE1IvKNiEwTkRNFpDTqz9JVMvmsCLnW9YH/k1kpnayqtvWiDef1/QWg3rYSaAy8fwPok2bbW+G80P226oE1gfePA8XZvgeZvhfApEAb8dukbH/u7roXuMgprXGffyXQEHfsFiA/2/cgw9+JR+M+8wpvCx77HPhOtu9Bpu9FgmvtHPddmZXK+Tai6kWISD7wCDASmAfsrqrluEC3P8UJy5akEX1dRKpx/4x9gQ+BrVW1EigHfo0TrD2Ay7v8QSIgk/fCYz4wFfgzcHBX+5tJMngv8gEBngaOAIZ47VYAGwMPefWOBCZ36UNEQIa/E08BJwLjgCpVrVDVClyIoROBVcAo4F8ikvXnbjf8fwSvVQbcDDQD09NqJNuqblt0G3AUsV8s24eUHxoo3zXFts/zzmsARoWU/59X3kwO/GrM8L0oCDmWsyOqTN0LoBoY10G54MTcH32X9Mb7kOS1jwm0/d3e+p1IcK3LvXbOB27HRlTrPEd4r9NU9eWQ8ntwQ32Aw1Ns269/j6p+EVL+N9xURz4wMcW2M0HG7oWqNnelY1kgI/dCVZep6hsdlCtwq/e2Atgo2bYzRCb/PzrjlcD+0IjbToduuRcish1uRPkxTqjSwoSql+ANr7/rvZ0aVsd7cDzmvd0jhbY3BIZ30vYK4IVU284EmbwXPY0cuBeNgf38iNtOmhy4DzsG9j+LuO2U6K57ISLFuB8qAhyrqo2dnJIQE6rew0bE/p7vdlDPLxsoIrVJtr1JyPkdtT02yXYzRSbvRU8j2/diZ++1CferOlt0+30QkVIR2UBEzgAu9Q4/r6rprdNER3fdi7O9a92iqs+mcf5aLB9V72FwYH9OB/WCZYOBxRlou0pEKrxRVjbI5L3oaWTtXojIKOA47+29qrq8q212gW65DyIyEGecEMYjOIvRbJPxeyEiWwKnAQu81y5hI6reQ2Vgv6GDesGyyoS1uq/tTNDT+ptJsnIvPJ+h+3BWZItwxjbZpLvuQwvu4byAttOe9wGnqWou/BjK6L0QkQLclF8BcKKqLkmte+0xoTIMI1K8B9VdOL+7NcBhqtrRL/deg6ouVNWBqjoQJ9LDcC4MPwRmisgxWe1g93A6sAXwqKr+M4oGTah6D/WB/bIO6gXL6hPW6r62M0FP628m6dZ74fnn/AM4EOeqcJiqPpFuexHS7d8JdcxW1bNwlrCFwHUisnlX2o2AjN0LERkL/BFnAfyr1LsWjglV72FuYH9IB/WCZXMT1upa28uzuD4Fmb0XPY1uuxcBkfoJbgrsZ6p6fzptZYCsfidU9UHgS9wz96io2k2TTN6La4Ai3ChyiYhUBDdidhESOF7YWaMmVL2HD3AhSqCtlV48ftn8FObLg5ZBybT9fpLtZopM3oueRrfcC0+k7sRFNfBF6t5U28kgufCd8B/260fcbqpk8l6M8l4vxI3C4jffx3J44NgJnTVqQtVLUNUG4CXv7V5hdbygsXt6b5OejlHVj4CvOmm7nJivSFanejJ5L3oa3XEvAiIVHEndk3pvM0e2vxNe2/5DPKvTzNm+F+lgQtW7mOK97iIi24aUTwBGe/t3pNi2X/+nIjIypPwEXPSBFiKIDxYBmbwXPY2M3QtPpO7CiVQzMDHXRCpARu6DZzzSGb8ABnr7zybbdgbJyL1Q1ZGqKom2wHW/DBy/IpmGbeslG27+dyYultZsvBhduB8kE4BlXtl/Q86dTCy218iQ8mqcf4gC7wFbeceLgOOB1V7Ztdm+D5m+F16dfnGbX//XccfLeuu9wEWauMsrWwNMyPZnzdJ92Bl4Hvg5MDSubAPgImJZBj4FSnvrvUjiureTRqy/rH95bIt2w0VD/iLwRVqJi9zsvw8N3Z/kwzk+zcdyXMQB/32upfnI5L3QJLfJ2b4PmboXwE6BsiZcRPmOtp/00vuwc9zffBWwkPbpTt5K9H3qLfciiWvejgWlNVR1FrAZcC7OCML/tTsDOAXYTtN0wFPVGbj0DZcDn+DMbVcCLwK/BPZW1dVd/AiRkcl70dPI0L0IPj8KgQGdbFlPHJih+zADF7j1VuBt3GikBmew8BnO2fenuFmIWV39DFHRk/4/xFM5wzAMw8hJbERlGIZh5DQmVIZhGEZOY0JlGIZh5DQmVIZhGEZOY0JlGIZh5DQmVIZhGEZOY0JlGIZh5DQmVIZhGEZOY0JlrNOIyM4ioiJinu9Zxv87iMjO2e6LkVuYUBk9DhGZHHio+VuriCwXkdki8j8RuUZEfiwiRd3ct5Fe/yZ353UNozeTTHh6w8hlFgT2S4HBuMyk2+NSYS8SkT8C12t4vLAG4KMI+zMSOMfbnxxhu+sC/t+hIau9MHIOEyqjR6OqA4PvvfxIY4HdcSk3RgHXAt8TkZ/Fi5WqvgaM6abuGh2gqvZ3MEKxqT+jV6GqLar6jqpehkul7SfxOww4PXs9MwwjXUyojF6LupTbRwBveodOF5HaYJ3OjClEZIyI3CgiH4tIg4isEpGvReQVEblARMYE6s4CpgXex6+j3R4oKxSR3UXkKhGZLiLzRKRJRL4RkcdF5FAvHXhYn9r0WUTWF5FbvX6t9tbpbhKRIR3dHxEpEpGjReQxEVngnTtPRF4WkbNFZFSC86pF5EwReVVElnjnfS0id4vIdh1ds5P+hBpTeOt+ftlIERkgIleKyBci0uj1/Z7g3yLF63bpforIJO/8Wd77HUXkEe9vuVJE3hSRo+LO2VdEnhSRhd736nUR+Uk6/V8nyHbyLttsS3UjkLgtyfo/Jpbo7ci4sp0TtYWbPmwMnNsELAm8b5MYEXgdWBwoi08ceGXYdb2tEaiPO/ZPIC+kX8Fzdwmct5xYJlkF5gBDEtyTUcA7gbqt3mdrDhy7IuS8bb3P4tdp9q4bbOf/0vy7+m3sHHd8ZKBsX9y6pJ/oL/j3WQZsnsZ1u3Q/gUle+SzgaKDFuw9L4/6eF3r1/+S9bwmpc1y2/79ycct6B2yzLdWN1IWqIvAAnhJXtnOitnDJIRWXuXiTwPES3LTiOcAvkm0vrt62wJ3APrikgn5uuFrgRGKpwE8MOTf4YF0MPASM8cqKgEMC4nFHyPlVwMeB838JVHtlhcB3gJOBk+LOG0lMqO8DxgEFXll/XAI+/8F+YBp/12SEajEuUed4r6wA2A2Y65U/n8Z1u3o/JxETztXAlUBd4O95OzFhOs37Lp4ZuOeDgKlenRX+cdsC9zjbHbDNtlQ3UhQq7xz/wfxi3PFQYfEevP7Da1AK1wltL43P6I8CP+3oGsAzhI+6fuOVN/hiEig7j9gobssU+nRfood1oM5JXp230vjMyQjVB0BpyLk/DNQZmuJ1u3o/JwXOvynk3Hzg80CdM0PqVHkipcDPuvLd6Y2brVEZ6wqLvdfaDmvFqMdN34D7xdvd/Md7XU9EOrr+BaraGnL8Ie+1FNggruxI7/VmVX0zmc54a3sHeW8v6qDqHd7r5iIyIJm2U+RSVV0VcnwqbmoWYNMutJ/O/QzS7t6oagvwtPe2EbgipM5y4GXv7WbJdnZdwczTjXWFUMOERKjqKhF5GrdO9ZiIXI8TjzdVtanjs5PskEglcBywH7ARUIObeotnCDAvQTOvJjg+N7C/VpxFZATO1wzgkRS6uz0x46tnEth5xDOCtn5uURD6eVW1WUQW4u5Vsj9Gkm6fBPczjsWq+lmCMv8+vK+qKzup06eD/q2TmFAZ6wr+P/+iFM45GngY2Bz4o7c1icjruF/Yt6jq4g7OT4iIfAf3K3to4HADbnHd/0Xvj0jKE7WjqvUJjjcHxCQofkG/sy+T7/FacQv2qzPKUmg/WUI/r0ez9xom9kmRxv0MkkzfMtr/3opN/Rm9HhGpAEZ7bxP94m2Hqn6FMxjYC7gKmIH7n/kucDHwqYj8IM1u3YYTqVnABKCvqparan91TsxBU+iURoMpoCnUzfdeV6mqJLk9m4E+G+sgJlTGusBexB60z6Zyoqq2qurjqvpbVR2Pm/aZCHyFG6XdlWo8QREZBuzgvT1UVe8PGZkNJDMEpxBHpnDefO+1VETWj647htE5JlRGr8YTkTO8t8uAf3elPVWtV9W7AN+BcwBtF+/XLsQnctgFhgX2Exkz7JZ2JzvAGyXO9t7+MIVT/0dsBPbTSDtlGJ1gQmX0WkSkFOfDsqV36EJVXZrkuZ2NkoKWZy2B/eWB/ZoE5y4L7G8ecu1K4KxOrt8VbvVejxaRLTus6aGq3xCzfDvVW2NLSHwEEMPoCiZURq9CRPJEZBMRORl4DzjUK/o7bl0pWXYQkZkicpKIbCQieV77IiI7ANd59WbjIjz4fEzMTProBKOq93FThwC3ishWgf5vj5uezKTl119xzszFwNMi8ksRqfKuXygi3/FCKJ0Sd97vccYoVcCLInKkiFQH+t5PRA4SkQeBuzPYf2Mdw4TK6NGIyPzAtgQXGeEd4FJcmKBvcWFpDlfPszIFNgUuwwlLo4h8ixOhl7yy5cBhnp8MsDa+4N+9txcDK0TkSxGZJSJ/9eoocALOymtjYLoXE24lboptDJCxuG+eZdte3ufqA9wILBGRxbiR4ke4MD9D4877HGeuPwuoA27xzxORemAh8ADwI+zZYkSImacbPR3fVNoPYTMfN1p5E2f+/Uiafk+v40Ln7AJsgzPPrsM5bH4KPIGL3Tc35NwTgK9x0SVGA8O94/38Cqr6qIjshAul812cKfd8r89/UdWPkvRVSgtV/dyb9jsK9zk3xY2UFuDM1v9LTHCD570pImNxTsMH4qYu++AE/BPcfXvYO98wIkFS/5FpGIZhGN2HDc8NwzCMnMaEyjAMw8hpTKgMwzCMnMaEyjAMw8hpTKgMwzCMnMaEyjAMw8hpTKgMwzCMnMaEyjAMw8hpTKgMwzCMnMaEyjAMw8hp/h/DE6DSCXAJRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAADtCAYAAAB6WSfQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtzElEQVR4nO3deXxU1f3/8XcmCdkTkhDIwhLEgGExbFGMaLB1KQhYUb4VKYJF6wL9ga18xfoVrFptbfVrxUJREYsoFhWVRdCKCEKRVUC2IEsgLCJbAknIRvL7Y77cmSHJZCDLnElez8eDh/dmzp175BNm3nPn3HP8KioqBAAAAJjK5u0OAAAAAO4QWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYLcPegn99i5rzyaVO93QFcosla6O0uAIDPsf2B2OLLJk2SX3WPcYUVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMFqAtzvgqfT0KI0alaTMzBglJQXL399Phw4Vac2aXL399mEtXXqizs/ZoUOo+vaNVo8ekerSJVzJySFq2TJIoaE2FRaW6/jxEm3efFqff35C7757WKdPl9V5HxqL9PQUjRp1ozIzuyopKUb+/jYdOnRSa9Zk6e23l2np0s31dm6bzaYBA3rpzjuvVXp6ihISYhQeHqz8/CIdOXJS69fv1ocf/kcLF65TeXl5vfWjMYlu316JvXsrMT3d/t9evRQUGWk9/nJysnL372+QvsSlpqr7qFHqcPPNimzdWoGhoTp96JB+2LRJ373zjrLmz1dFRUWD9MUXUDvfRv18R2KilJYmtWsnRUZKfn7SmTPSoUPSli3Svn31c15/f6lVK/v5ExLs/42Lk2z/d4kyO1uaNat+zl2f/Nz9Mvn5Lfb6b1pIiE0vvpiqhx5q67bdhx/+oPvv36pTp0rr5LwZGc21atU1HrfPyyvV44/v0rRpB+rk/HVjqrc7oJCQIL344q/00EMD3Lb78MP/6P77p+jUqfw6PX9aWnvNmvWIrryyfY1tt27dr5Ej/1cbN+6p0z5cisla6O0uVCllwADdPmuWQmNj3bZriDdNP5tNmZMm6fonnpAtoPrP3tnLl+ujESOUl5NTr/0xHbXzbdTPM7Y/eD22KCBAuvlmqXdv9+127JAWLJCKiuru3IMGSVdeaQ+t1TE5sE6aJL/qHjN6SIDNJs2d28MlrBYUlGndulytXn1KeXmOcHrHHfH6/PN0hYTUzf+Sn5/r31lZWbm+/75AK1ee1BdfHNe6dbk6c8ZxRTUqKlBTp3bRX/7SqU7O3xjYbDbNnfuYS1gtKCjSunXfa/XqncrLK7B+fscdGfr886cVEhJUZ+fv2bODli9/3iWsnj1brPXrv9cXX2zShg27VVRUYj3WtWs7ffXVc7rqqo511ofGJiwursY3zIYy4NVX1W/yZOsN81xJiY58+60OrFypwuPHrXbJmZkatXy5wuLivNVVI1A730b9fIOfn3Tnna5htaTEflX14EHXcJqaKv3yl/aAW1eio92HVV9m9JCAyZNTNHBgS2t/+vQDmjgxS7m59qAYGuqviRMv05NPXi5J6t07StOmddGoUd/V+twVFRVavz5PH398VMuWndC6dXkqLXX95Obv76fbbmupF1+8QsnJoZKkRx+9TF99dVKLFh2rdR983eTJwzRwYLq1P336Ek2c+JZyc+1BNTQ0SBMn3qknn7xLktS7d4qmTXtIo0a9XOtz+/vbNGvWbxUVFSZJKi0t0+TJ7+qVVxaooMDxihEeHqLx4wdr8uRhCgjwV0REqGbNekRdu45VWdm5WvejMTu1d68Or1+vw+vXS5JueuGFBjt3j3vvVfpDD1n7Oz/+WIvGjNGZw4clSf6Bger1wAO65aWX5B8YqOj27TV07ly9dcMNDdZHk1E730b9zJWZKXV0uuaxYYO0dKkjqAYGStdeK11/vX0/MVG69Vbpk0/qth9lZdLRo9Lhw9KRI1LnztLll9ftORqasUMC4uODtHv39QoLs2fqWbMOaeTILVW2ffrpFCu0lpdXqEePVdqy5UyD9TUhIUibNl2rli3tVwe/+uqEbrhhbYOdv3reGxIQHx+t3btfU1hYsCRp1qwvNXLk/1bZ9umnh1uhtby8XD16jNOWLdm1Ov/NN/fQZ589be2PG/eaXnllQbXtH3nkNr300n3W/oABT2nx4g216kNtmDokILZjRzVPTtbhdet09tQp6+fJmZka9dVX1n59fi0ZEBys/7d7tyKTkiRJ+5Yt06wbb1RFFeOPe/zqV7ptxgxrf87gwcpaUP3vQWNG7Xwb9fOMN4cEhIdLY8dKzZrZ9zdvrj6I9uvnCK0VFdJrr9kDZm116CAVFtqfy7ksgwdL3bvbtxkSUMfGjWtnhdWCgjKNH7+j2rbPPLNbBw6clSTZbH567LHLGqSP5x05UqwpUxwvENddFyN//2r/zpuEceMGW2G1oKBI48e/Xm3bZ575lw4csF+RttlseuyxO2t9/uuv72pt5+bm6+9/X+S2/ZQpC5Wb6xg/e+21nWvdh8boxK5d2vP55y5vmA2t+8iR1htmRXm5Fj38cJVvmJL07Ztv6uA331j7fSdObJA+moja+TbqZ76rr3aE1ZIS6bPPqm+7YoWUl2ff9vOTMjLqpg979tivqDbG+4eNDaxDhsRb23Pn/uD2ZqrS0grNnHnQ2h84sKUCAxs2MG7YcNra9vf3U1xcswY9v2mGDHHcsDZ37kq3N1OVlpZp5swvrP2BA9MVGFi70SpxcY67ZnfsOKhz59z/6y0rO6edOx2/Q87HwyypQ4ZY29nLl+v4zp1u26+fPt3abt2njyISEuqtb3CP2vk26ufeFVc4trdvd38zVXm5tGmTY79jR8dd/KiakX89KSmh6tgxzNpfsqTm8aCLFzsGekdGBigzM6Ze+ladgADXgOx8Q1ZTk5KSqI4dk6z9JUtq/mrd+ev3yMhQZWZ2ddO6Zvn5jleKoKBAj44JDnZ8yKjr2QpQNwJDQ5Xcr5+1v3vJkhqP2b14sbXtZ7Mp5dZb66NrqAG1823Uz72YGMn5nrjdu2s+xrlNUJCUnFzn3WpUjAysaWmuV7dWr86t8ZiNG/NUXOy4inbhc9Q354C8efNpFRQ03Rt20tJcp5BavTqrxmM2btyj4mLHVfQLn+NirVnjOGfXrm0VExPhtn1sbKQ6d25j7a9cub1W50f9iOvcWf7NHB8sDq5eXeMx+UeP6pTThIfxaWn10je4R+18G/Vzr1Ur1/2DB6tu5+zIEfvNUdU9B1wZGVhTU8Ot7eLicuXk1DxJWWlphXJyzjo9R5ib1nWrV69Il6m3/vrXepoN2EekpjqCX3FxqXJyar5CXlpa5tIuNbV1rfrwySdrdPCg/ap7s2aBevnl+6tt6+fnpylTfq1mzexXYrdtO6BPP/XeDVeoXlxqqsv+yT2ezZl7yqldiwueAw2D2vk26uee88xdZWXS6dPVtz2vvNy1XYsWdd+vxsTIwJqcHGJtHzzo+Yy6Bw442jo/R32IiAjQVVdF6YUXOunrr/soNNQ+8dlrrx3Q7NmH6/XcpktOdkxFdj40euL8jVf256jdR83i4lING/YX60aqESNu0LJlz2nAgN6KiYmQzWZTbGykBg26SsuXP69hwzIlSdnZR3XHHc+z4pWhmjt9Z3autFT5R454dFzeAceCHs353s0rqJ1vo37uRUU5tj0Jq+edv/FKkpo3r7PuNEpGzsMaGenolvPiADVxXho1IqJu/9emTeuiBx+sfrWt7OxCPf30Hpebv5qqyMhQazsvr9Dj406fdlwhj4io/QeOlSu3KyPjv/XGG79RRkaq+vXrpn79ulXZtqioRO+/v0oTJrypo0dza31u1A/nJShLzpzxeNnHYqd3kKAI98NDUD+onW+jfu4FOa15U1zs+XHObZs17Xu1a2TkFdawMMcyDUVFnl/pOnvWMW40PLzhsvgPP9intZo374cGO6fJzk9nJcllJamanD3r+JcbHh7spqXnduzI0W23Pat//GOx23Yff/yNpkxZQFg1XGCYY6hP2UWsZ1h61vFhqFl4uJuWqC/UzrdRP/ecw2bZRdxz7dyWwOqekYHVeUqqsjLPJwF2btusWd1Oa7VtW76WLDmmJUuO6d//ti/Nmp9v/02Ljw/Siy+mat++fho6NL6GZ2r8AgMdHzguZrWosjLHh5NmzWr/gcNms+mJJ/5L2dkz9OCD/SXZhwps2LBbS5du1oYNu60bve6663qtXfuS3nnnUZfADbP4BzpmfCi/iHcF57b+vCt4BbXzbdTPPecpqS5mRJlz28a6pGpdMXJIQGGhI+QEB3ueqZ3b5ufX7V36r766X6++6rp6iM0m3XJLnP70p4668spIRUcH6r33usvff7Pee8+z8T2NUWGh40qp81RRNQkOdrwgOk9Ldaneemu8RoywLwdYVFSiJ5+cralTP3XpX1hYsMaMuVVPPz1cQUGBuvvuTCUmxujGG/+nxrlb0fBKCx1DTAKCPf9g4dy2JJ8py7yB2vk26udeqdPoxYCLSFbObUs8/0KySTIysDqHzZAQzz9ynL/xyf4c9T8Panm5tHjxMX355Ql9+mlv/eQnsbLZ/PSPf3TRZ58dd7vYQWPmHDZDQjwPrKGhjkFAtQ2s9957oxVWy8vLNWTIc1UutVpQUKQXXvhQ27cf0IIFkyRJ/fp10/jxt+nFFz+qVR9M1qpbN4/WH5/dv38D9MZzzm94ASGej3MODHWMq/b1N01q59uoX+PkHDYvJrA6XbgmsNbAyMB6/LijagkJQW5auoqPd7Q9caLhwmJxcblGjdqiPXsyFRhoU1RUoIYPT6x0RbapOH7cMcg+IcHzBRzi46Ot7RMnLuI2yypMnOhY3nX+/LVVhlVnCxeu0/z5azR48NWSpHHjBjXqwBoSE6PLf/Yzb3fjohUed8w6ERQRoWZhYSopKKjxuPB4x1CdwhMn6qVvDYXa+Tbq1zg5XYDWxdxb5jQ0WE7DfVEFI8ewZmU5/hG0aNFMISGedbNNG8dXDzt31vwPqS7l5BTpm29yrf1rrmneoOc3SVaWY6aEFi0iFRLi2YeONm0ck9A5L5N6sdq0iXNZaWv+/DUeHffJJ452bdrE1XpqLdS941mui1BEta1+5g6Xdm0ccwPXtJwk6ge1823Uzz3nLB4a6vlVVufpsI57Pgtkk2TkFdbt212/NujePbLG1a4SE4PUsqUjGO3Y0fBfPTgvcNCiReMdXF6T7dtzXPa7d2+v1avdv1AlJsaoZcvm1v6OHZceWJOSXK/q5uR49ipwYbv4+ObKzj56yf0wWfby5XrKr25vTGwIx7a7rkAW3727ju3Y4fYYW0CA4rp0sfaP19DedNTOt1G/xunCsBkfX/NqVxERrldYCazuGXmFde3aPJcpqvr2jXbT2u6661xDyooVJ+u8XzWJinLk/9zcpjl+VZLWrt3lMkVV376dazzmuuu6uOyvWLH1ks9fXOw6ftnTcbTOY2gl15vHYIbc7GyXicjb9u1b4zGJvXqpmdO7wv4VK+qlb3CP2vk26ufeoUOuN155cgH6wjb7m+YoQo8ZGVgLC89p6VLH9fXhwxNrPMa5zZYtp7VvX8MOBvH391OfPs2t/T17PJ8wv7EpLCzW0qVbrP3hw/vVeIxzmy1b9mnfvku/snn4sOs4qd69Uzw6rlevDhc8T8N/6EHNsubPt7Y7Dx3qMt1OVboNH25tF5440ajfNE1H7Xwb9ateaam0z2lV9q5daz6mm9M6NkePSrm5dd6tRsXIwCrJZcWotLRIDRzYstq2PXpEqn9/x/jHmTMP1WvfqvLII8mKjXVcyVu06Jib1o3fzJlfWNtpae01cGB6tW179Oig/v17OR27tFbnPno0V7t2OX4HRo++qdLV0wuFh4fovvtutva3bTvgcvMYzPHtzJnWdlhcnHo98EC1bSOTktR95Ehrf8vs2So/V7dT3sFz1M63UT/3Nm1ybMfHSx07Vt82Pl66/PKqj0XVjA2s8+Yd1YYNjkV2p0/vok6dwiq1S0gI0uzZaQoIsP+vHDxYpGnTDlRqd15FRX/rz8yZVS/TKUnz5/fS1Vc3r7GfNps9rP7pT52sn61YcVKrVp2q8djGbN68/2jDht3W/vTpY9SpU+tK7RISYjR79m8VEGCfkuzgweOaNu3Tap+3omKB9WfmzPHVtnMOzElJsXr//YkKD696KpaIiBB98MFElxkNZsz4vNrnRt1LzszUUxUV1h/nN7oLHdm4UTvmzbP2f/rcc2qTkVGpXVBEhO6YM8daUrKkoEArn3++7jvfxFE730b96s7OndLhw479W2+VYmMrtwsPl26/3bHYwOnT0vr11T/vpEmOP4MH122ffYmRN12dd999W7Vy5dUKCwtQYmKw1qy5RlOnHtDXX5/SuXMVSk+P0tix7azprMrKynX//d+puLj2E74PGtRSgwa11J49hfrss2PatOmMDh4s0unTZWrWzE+tWgWpZ89I3XlnvNq3d8wz9+OPxbrvvu9qff7G4L77pmjlyj8rLCxYiYmxWrPmr5o69VN9/fU2nTtXrvT0FI0dO9Cazqqs7Jzuv/9Va/Wp2vjb3xbonnt+otRU+x2qAwb0VlbWNM2Y8W99802W8vIKFR0dpj59rtDo0Te5TKm1efM+/f3vi2rdh8Zq7M6dat6uncvP/Gyun33HZmVJF6w1nrt/v1694oo66cOS8ePVJiND4fHxCoqI0MilS7XxjTe05/PPVVJQoPi0NF01dqyiL7vMOuazRx5R/tHGeROdp6idb6N+5lu4UBo1yr7MakSENHq0PYweOGAvS2KilJ5uD62SfT73hQulurr43Lat9MtfVv65869Ju3bS739fuc3XX9v/mMrowLpp02nddddmzZmTpvDwAEVFBerxxzvo8ccrty0tLdeYMdu1ZEnd3mbXoUOoHn64Xc0NJW3deka/+MUmff990x2/6mzTpr26664XNGfOBIWHhygqKkyPPz5Ujz8+tFLb0tIyjRnzDy1Z4n6+VE+dPVusW26ZrMWLn1KXLvaR7YmJsXryybvcHvftt3t0661Pq6Sk/hee8FUBQUE1rnQTEFR5CEZVP7tUeTk5mjN4sO5etEhhcXEKCA7WVWPH6qqxYyu1rSgv14pnn9WG11+vs/P7Kmrn26if+X74QZo3TxoyxB5ag4Ol6u5PO3dOWrxY2r276scvhZ9fzVNqVdfGZux37naGd09auPBH9eixSosXH9O5cxVVtlm16pQyMr7R66/nVPn4pZgwYaeWLTvhMltBdTZuzNPDD29Tjx6rKk3J1dQtXLhOPXqM0+LFG3Sumo+Qq1ZtV0bGf+v11z+r03Pn5BxTr17j9fvfz9KBA+7HFGdnH9WECW/q6qsf1ZEj3GzlCw6tW6dp3bpp63vvqay46hkdfti8WbP799eyyZMbuHdwh9r5Nurn3q5d0muv2YNoeTVf+ObkSDNnShs3NmzffJlfRUXVIVCS/PwWV/+gFyQkBKlv32glJQXL399Phw4Vae3aPO3dW39XNAMC/NS1a4Q6dgxVQkKwwsP9VVJSrry8Mu3ff1YbN57WsWOmrqc21dsdcJGQEKO+fTsrKSlW/v42HTp0QmvX7tLevT80yPlTU9uoZ88OiouLUlhYkPLzi/Tjj3nauHGPy2IHJpishd7ugk8JiY5Wcr9+imzdWgEhITpz+LB+2LRJP2699OnR0DConW8zrX62PxgVWxQebv+aPiLCfgXz9Gn7ONdTTfs2l2pNmqRqJyn2qcCKi2VWYIXnCKwAcPFMC6y4OO4Cq/FDAgAAANC0EVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYLQA9w9PbZheAAAA1FL5ZD9vdwG1Mami2oe4wgoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYLcDbHfBUenqKRo26UZmZXZWUFCN/f5sOHTqpNWuy9Pbby7R06eZ6O7fNZtOAAb10553XKj09RQkJMQoPD1Z+fpGOHDmp9et368MP/6OFC9epvLy83vrhy7xZP0m6/fZrdPvt16hPn05q1aq5/P1tOnz4pL77Llvvvfe1Pv74G5WWltVrHxqT6Pbtldi7txLT0+3/7dVLQZGR1uMvJycrd//+BulLXGqquo8apQ4336zI1q0VGBqq04cO6YdNm/TdO+8oa/58VVRUNEhffAG1823erF/zdu00Pjv7oo/7YuJErfzzn+u+Qz6I+l064wNrSEiQXnzxV3rooQGVHrviilBdcUVrjRz5U3344X90//1TdOpUfp2ePy2tvWbNekRXXtm+0mPR0eGKjg5X585tdc89P9HWrfs1cuT/auPGPXXaB1/m7fqlprbRnDkTlJZWuX4pKYlKSUnUkCEZ2rRpr0aMeElbtzbMG7WvShkwQLfPmqXQ2Fhvd0V+NpsyJ03S9U88IVuA60tZbEqKYlNS1GXoUGUvX66PRoxQXk6Ol3pqBmrn20yqHy4e9as9owOrzWbT3LmPaeDAdOtnBQVF2r49R2Vl59S5cxtFRYVJku64I0Pt2sXp+usf19mzxXVy/p49O+jLL/9onUOSzp4t1rZtB5SbW6Do6HB16dJWwcHNJEldu7bTV189pxtvfFJr1+6qkz74Mm/XLz09RV988awiI0Otn+Xnn9V33+1XcXGpUlISlZRkf/Ho3v0yff31n9S372Patu1AnZy/MQqLizPmBXfAq68q/aGHrP1zJSX6cds2lRYUqMUVVyi0RQtJUnJmpkYtX643rr5aBceOeau7XkftfJtJ9Ttv95IlHrU7tXdvPffEfNSv9owOrJMnD3MJO9OnL9HEiW8pN7dAkhQaGqSJE+/Uk0/eJUnq3TtF06Y9pFGjXq71uf39bZo167dWoCotLdPkye/qlVcWqKCgyGoXHh6i8eMHa/LkYQoI8FdERKhmzXpEXbuOVVnZuVr3w5d5s37R0eH66KMnrLBaUlKqxx77p6ZOXaSSEsdX/z/9aZpee22sLrssXs2bh2vRosnq0mWMS41RtVN79+rw+vU6vH69JOmmF15osHP3uPdel8Cz8+OPtWjMGJ05fFiS5B8YqF4PPKBbXnpJ/oGBim7fXkPnztVbN9zQYH00GbXzbd6sn7PZ/ft75by+jvpdGmMDa3x8tH73u59b+7NmfakHH/y7S5vCwmJNmvSOJFmhZ8SIG/TSSx9ry5bsWp3/pz9NU5cuba39Rx99U6+8sqBSu/z8s3r22X+poKBIL710nySpU6fWuumm7lq8eEOt+uDLvF2/Rx+93bp6KknDh7+oDz5YVand0qWbdd11j+nbb/+mli2bq127lvrd736up59+r1bnb6xyVq/W27fcosPr1unsqVPWz5MzMxusDwHBwbrhmWes/X3Llulfd9yhCqfx4+dKS7X21VdVWlio22bMsPexXz91GjRIWQsq/ztuCqidbzOhfrh01K/2jJ0lYNy4wQoLC5Zk/xp5/PjXq237zDP/0oED9q+LbDabHnvszlqf//rru1rbubn5+vvfF7ltP2XKQuXmOsZfXntt51r3wZd5u36jR99kbS9ZsqHKsHre4cMnNXnyu9b+7353uzXMA65O7NqlPZ9/7vKC29C6jxypyKQkSVJFebkWPfywS+Bx9u2bb+rgN99Y+30nTmyQPpqI2vk2E+qHS0f9as/YwDpkyDXW9ty5K93ejFNaWqaZM7+w9gcOTFdgYO0uHsfFOe7a27HjoM6dc3/3f1nZOe3cebDK45sib9avS5e2atUq2tqfM2dFjce8994Ka4aHyMhQ/exnPS/5/KhfqUOGWNvZy5fr+M6dbtuvnz7d2m7dp48iEhLqrW9wj9oBuFRGBtaUlER17Jhk7S9ZUvNX685fv0dGhiozs6ub1jXLz3eMYQwKCvToGOercnV9t7sv8Xb9kpNbuex7MrwgN7dABw+esPZvu63PJZ8f9ScwNFTJ/fpZ+57cNLB78WJr289mU8qtt9ZH11ADagegNowMrBdOQbR6dVaNx2zcuEfFxaXVPsfFWrPGcc6uXdsqJibCbfvY2Eh17tzG2l+5cnutzu/LvF2/qKhQl/0zZ856dNzp04XWdo8el13y+VF/4jp3ln8zxwfDg6tX13hM/tGjOrVvn7Ufn5ZWL32De9QOQG0YedNVaqoj+BUXlyonp+bpTEpLy5STc0yXX574f8/RulZ9+OSTNTp48Lhat26hZs0C9fLL9+uee16qsq2fn5+mTPm1mjWzX4ndtu2APv206d5w5e36XRhQIyJCPDrOefqrTp2SZLPZWAjCMHGpqS77J/d4NufxqT17FN3e/iGoxQXPgYZB7Rqf2//5T/tQjcRE+dlsOnvypE58/732L1+uzW+/bcx0SKiar9XPyCusycktre2DB497fNz5G3fsz9HKTcuaFReXatiwv1g3Uo0YcYOWLXtOAwb0VkxMhGw2m2JjIzVo0FVavvx5DRtmv9MvO/uo7rjj+SYddLxdP+ev9iWpW7d2NR4TFRWm1q0dswoEBzdTQkK0myPgDc2Tk63tc6Wlyj9yxKPj8g445tZ1fg40HGrX+KTdc49iO3ZUs/BwBYaGKrJ1a7W/4Qb1e+op/SYrSz+fOVPNwsO93U1Uw9fqZ+QVVucrXXl5hW5aujp92nFlzdOrau6sXLldGRn/rTfe+I0yMlLVr1839evXrcq2RUUlev/9VZow4U0dPZpb63P7Mm/Xb/PmfcrLK7Dm0L3rruv19tvL3B7zi19cJ5vN9fNbXfwOoW45L2FYcuaMx8t2Fp8+7XiOCPfDe1A/qF3jc+bIEeXu26fSwkIFN2+uFqmpahZmf921BQSo+6hRan3NNXorM1P5R496ube4kK/Vz8jAen46JMkeBD3lvEJSeHiwm5ae27EjR7fd9qyeeeaXevDB6ifZ/fjjbzRlyoImH1Yl79evvLxcc+assOo1YEBv3X77Nfroo6rHzMXHR+sPf7i70s8JrOYJDHOsOldW5PniDqVnHR+GTLpi0JRQO99XUVGhg2vW6NsZM7Rr0SJrsYfzbAEBSunfXz/54x/Vqpv94k6LTp00bMECzcjIUHlZWVVPiwbi6/UzckhAYKC/tX0xq0WVlTm+hm/WrPZZ3Gaz6Ykn/kvZ2TOs8FNcXKoNG3Zr6dLN2rBht3Wj0F13Xa+1a1/SO+886hLYmiIT6venP32gM2ccV3fnzJmg3/xmUKXn7devm1au/LPi4yt//V/bqdFQ9/wDHTN2XMyLp3Nb5xt/0HCone/LO3BAb/Tpow2vv14p7Ej2WmUtWKDXr7pK3zvN8JCUnq6eo0c3ZFdRBV+vn5GBtbDQcaXtYiZwDw52vCA6T0t1qd56a7yefXaEwsKCVVRUogkT3lRMzDD17v2Ibrzxf9S79yOKjb1bjz32lhVc7747UwsXTpK/v5F/tQ3ChPrt3/+jRo36mxWYg4IC9corv9bx4+9o1aoX9OWXf1ROzkwtW/acOnRIUHl5eaXFIfLyCmrVB9S90kLHh5CAYM8/GDq3LclvulPOeRO1azrKior0wS9+4fI18tXjxnmxR7gYptbPyEtIzmElJMTzwBMaGlTlc1yKe++9USNG2NeuLi8v15Ahz1W51GpBQZFeeOFDbd9+QAsWTJJkv2o3fvxtevHFj2rVB19lQv0kad68/6h//6c0a9YjSkiIkSRFRIQqI8P1TuO8vAKNGfMPlZaWacwYxzyPubmNN7C26tbNo/WrTVtr2jmwBIR4PmQjMNQxrtrXQw+1822+Wr+LVXzmjNZPm6Z+Tz0lyT5LRFSbNsrLyfFux2qJ+nmvfkYG1uPHHYPszwcNTzh/rXvixGk3LWs2caJjedD589dWGVadLVy4TvPnr9HgwVdLksaNG9RkA6sJ9Tvviy826bLL7tfIkT/Rrbemq3v39mrRIlLFxaXKzv5R8+ev1WuvLdGhQydcloQtKipp1OORQ2JidPnPfubtbly0wuOOWSeCIiLULCxMJQU1f7AIj493PMeJE25amo/a+TZfrd+l2LdsmRV4JCm2Y0efD6zUj8DqIivLscRpixaRCgkJcrkhpzpt2rSwtp2XSb1YbdrEuazUNH/+Go+O++QTR2Bt0yZOycmtlJ3t/TvrGpq363ehoqISTZ++RNOnu19ZJz09xdretGnvRY2/RcM4nuW6CEVU27Y6tmNHjcdFtXHMDVzTcqCoH9Su6blw6rLQFi2qaQkTmVY/IwPr9u2uCb579/Zavdr9C1ViYoxatmxu7e/YcemBJynJ9apgTo5nc4le2C4+vnmTDKzert+lsNlsuvZax1CBNWt2Nej5G1r28uV6ys/P2924aMe2u64gF9+9e42hxxYQoLguXaz94x6EJJNRO9/mq/W7FM7DOSTXGR98FfXzHiPvDFq7dpfLFbm+fTvXeMx113Vx2V+xYusln7+42PUOVk/HYTqPwZRcbz5qSrxdv0tx003dXYYkvPvu8gY9PzyTm53tMpF82759azwmsVcva25BSdq/YkW99A3uUbumJ66z62t/gQFzecJzptXPyMBaWFispUu3WPvDh/er8RjnNlu27NO+fZf+F3v4sOs4qd69U6pp6apXrw4XPM/JS+6DL/N2/S7F738/1OX8a9c27iusvixr/nxru/PQoS7TJVWl2/Dh1nbhiROEHi+idk1L12HDrO3SwkId+fZbL/YGF8u0+hkZWCVp5swvrO20tPYaODC92rY9enRQ//69nI5dWqtzHz2aq127Dln7o0ffVOnq6YXCw0N03303W/vbth1wufmoqfFm/S7Www8P0PXXd7X2J078Z4OeHxfn25kzre2wuDj1euCBattGJiWp+8iR1v6W2bNVfo6xyd5C7ZqOToMHq9OgQdb+959+qnMlni8kA+8ysX7GBtZ58/6jDRt2W/vTp49Rp06tK7VLSIjR7Nm/VUCAfbL6gwePa9q0T6t93oqKBdafmTPHV9vOOXAlJcXq/fcnKjy86qlYIiJC9MEHE13uiJ8x4/Nqn7sp8Hb9JGnw4KvVokVktY/7+flp/Pjb9Morv7Z+9s47X9U4IwTqXnJmpp6qqLD+OAeVCx3ZuFE75s2z9n/63HNqk5FRqV1QRITumDPHWhK0pKBAK59/vu4738RRO9/maf3iOnfW4DfeUGzHjjU+55XDh+uOd9+19svLyrRs8uQ66zMcmlL9jLzp6rz77puilSv/rLCwYCUmxmrNmr9q6tRP9fXX23TuXLnS01M0duxAa+xhWdk53X//q9Yk/rXxt78t0D33/ESpqfY7VAcM6K2srGmaMePf+uabLOXlFSo6Okx9+lyh0aNvchn/uHnzvkqT0DdF3qyfJD30UH/NnfuYli7drGXLtmj79hzl5RUqJiZcV16ZrGHDMtWlS1ur/YoVW/XrX79aJ+duzMbu3Knm7dq5/MzP5vrZd2xWlnTBWvG5+/fr1SuuqJM+LBk/Xm0yMhQeH6+giAiNXLpUG994Q3s+/1wlBQWKT0vTVWPHKvqyy6xjPnvkESPWw/YmaufbvFk//8BA9Rw9Wj1Hj9bhDRu078svdXTLFhUcParSwkIFRUWp1ZVXqsvQoUro2dPl2CXjx1e66a4pon61Y3Rg3bRpr+666wXNmTNB4eEhiooK0+OPD9Xjjw+t1NY+6fs/tGRJ3VwdO3u2WLfcMlmLFz9lhZrExFg9+eRdbo/79ts9uvXWp1VSwprJ3qzfeUFBgRowoLcGDOjttt1HH63WiBEvNdkb5S5GQFBQjSsVBQRVHkJT1c8uVV5OjuYMHqy7Fy1SWFycAoKDddXYsbpq7NhKbSvKy7Xi2We14fXX6+z8vora+TYT6ifZb4ZL7NWrxnalhYX67Le/1frp0+v0/L6K+tWOsUMCzlu4cJ169BinxYs36Fw145dWrdqujIz/1uuvf1an587JOaZevcbr97+fpQMHjrltm519VBMmvKmrr35UR440zZutquLN+i1dull79hyp9vHy8nKtXr1TP//5HzVkyHMqKKj96lpoOIfWrdO0bt209b33VFZc9QeNHzZv1uz+/Y34OgsO1M73nDlyRJveeksndtV8Q2rx6dNaN3Wqpl15pTFhp6lrDPXzq7jg0rPLg36Dqn/QCxISYtS3b2clJcXK39+mQ4dOaO3aXdq794cGOX9qahv17NlBcXFRCgsLUn5+kX78MU8bN+5xmSwfVfNW/dq2jVNaWnu1atVcLVpE6syZszpy5JTWrMnSoUNmrpwzWQu93QWfEhIdreR+/RTZurUCQkJ05vBh/bBpk37c2rDTo+HiUTvfExobq1ZpaQqLi1NoixYKiopSaWGhzp48qR+/+05Ht2zhBjmDmVy/pyoqqp3k1qcCK9BUEFgBAE2Nu8Bq/JAAAAAANG0EVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAo/lVVFR4uw8AAABAtbjCCgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEb7/0eljsnPjbU8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#in kcal/mol\n",
    "\n",
    "gHBfix=[]\n",
    "for i in range(12):\n",
    "    gHBfix.append(0.0)\n",
    "\n",
    "gHBfix[1]=1.0\n",
    "gHBfix[6]=-0.5\n",
    "gHBfix[8]=-0.5\n",
    "gHBfix=reorder(gHBfix)  \n",
    "data=[gHBfix]\n",
    "cmap = plt.cm.get_cmap('gist_gray')\n",
    "\n",
    "for l,lambdas in enumerate(data):\n",
    "    newlambdas = data[l].reshape(-1,12)\n",
    "    l=newlambdas[0].reshape(2,6)\n",
    "    image = plt.matshow(l,cmap='seismic_r',vmin=-0.5,vmax=0.5)\n",
    "    for (i, j), z in np.ndenumerate(l):\n",
    "        if z > 0:\n",
    "            c='white'\n",
    "        elif z == 0:\n",
    "            c='black'\n",
    "        else:\n",
    "            c='white'\n",
    "        plt.text(j, i, '{:0.1f}'.format(z), ha='center', va='center', fontsize=35, color=c)\n",
    "    plt.axis('off')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "cmap = plt.cm.get_cmap('seismic_r')\n",
    "#cmap = plt.cm.get_cmap('seismic')\n",
    "plt.figure()\n",
    "\n",
    "count=0\n",
    "for eta in range(-10,11):\n",
    "    c=cmap(count/20)\n",
    "    eta/=10\n",
    "    x=np.linspace(0,0.4,500)\n",
    "    y=[]\n",
    "    for e in x:\n",
    "        y.append(((opt_func_plumed(e))+1)*eta)\n",
    "    plt.plot(x,y,linewidth=5,color=c)\n",
    "\n",
    "    count+=1\n",
    "#plt.xlabel('Distance in $\\AA$', fontsize=16)\n",
    "plt.xlabel('Distance in nm', fontsize=25)\n",
    "plt.locator_params(nbins=7)\n",
    "#plt.ylabel('gHBfix potential in kcal/mol', fontsize=16)\n",
    "plt.ylabel(r'Potential in $\\frac{kcal}{mol}$', fontsize=25)\n",
    "plt.xticks(fontsize=25)\n",
    "plt.yticks(fontsize=25)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "error= np.load('./data_produced/chi2.npy',allow_pickle=True)\n",
    "lambdas = np.load('./data_produced/parameters.npy',allow_pickle=True)\n",
    "\n",
    "gHBfix=[]\n",
    "for i in range(12):\n",
    "    gHBfix.append(0.0)\n",
    "\n",
    "gHBfix[1]=1.0\n",
    "gHBfix[6]=-0.5\n",
    "gHBfix[8]=-0.5\n",
    "gHBfix=np.array(gHBfix)   \n",
    "data=lambdas /1.666666666666666\n",
    "data=data+gHBfix\n",
    "cmap = plt.cm.get_cmap('gist_gray')\n",
    "\n",
    "for l,lambdas in enumerate(data):\n",
    "    newlambdas = reorder(data[l]).reshape(-1,12)\n",
    "    l=newlambdas[0].reshape(2,6)\n",
    "    plt.matshow(l,cmap='seismic_r',vmin=-0.5,vmax=0.5)\n",
    "    for (i, j), z in np.ndenumerate(l):\n",
    "        c='white' \n",
    "        plt.text(j, i, '{:0.1f}'.format(z), ha='center', va='center', fontsize=35, color=c)\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAADtCAYAAAB6WSfQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtzElEQVR4nO3deXxU1f3/8XcmCdkTkhDIwhLEgGExbFGMaLB1KQhYUb4VKYJF6wL9ga18xfoVrFptbfVrxUJREYsoFhWVRdCKCEKRVUC2IEsgLCJbAknIRvL7Y77cmSHJZCDLnElez8eDh/dmzp175BNm3nPn3HP8KioqBAAAAJjK5u0OAAAAAO4QWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYLcPegn99i5rzyaVO93QFcosla6O0uAIDPsf2B2OLLJk2SX3WPcYUVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMFqAtzvgqfT0KI0alaTMzBglJQXL399Phw4Vac2aXL399mEtXXqizs/ZoUOo+vaNVo8ekerSJVzJySFq2TJIoaE2FRaW6/jxEm3efFqff35C7757WKdPl9V5HxqL9PQUjRp1ozIzuyopKUb+/jYdOnRSa9Zk6e23l2np0s31dm6bzaYBA3rpzjuvVXp6ihISYhQeHqz8/CIdOXJS69fv1ocf/kcLF65TeXl5vfWjMYlu316JvXsrMT3d/t9evRQUGWk9/nJysnL372+QvsSlpqr7qFHqcPPNimzdWoGhoTp96JB+2LRJ373zjrLmz1dFRUWD9MUXUDvfRv18R2KilJYmtWsnRUZKfn7SmTPSoUPSli3Svn31c15/f6lVK/v5ExLs/42Lk2z/d4kyO1uaNat+zl2f/Nz9Mvn5Lfb6b1pIiE0vvpiqhx5q67bdhx/+oPvv36pTp0rr5LwZGc21atU1HrfPyyvV44/v0rRpB+rk/HVjqrc7oJCQIL344q/00EMD3Lb78MP/6P77p+jUqfw6PX9aWnvNmvWIrryyfY1tt27dr5Ej/1cbN+6p0z5cisla6O0uVCllwADdPmuWQmNj3bZriDdNP5tNmZMm6fonnpAtoPrP3tnLl+ujESOUl5NTr/0xHbXzbdTPM7Y/eD22KCBAuvlmqXdv9+127JAWLJCKiuru3IMGSVdeaQ+t1TE5sE6aJL/qHjN6SIDNJs2d28MlrBYUlGndulytXn1KeXmOcHrHHfH6/PN0hYTUzf+Sn5/r31lZWbm+/75AK1ee1BdfHNe6dbk6c8ZxRTUqKlBTp3bRX/7SqU7O3xjYbDbNnfuYS1gtKCjSunXfa/XqncrLK7B+fscdGfr886cVEhJUZ+fv2bODli9/3iWsnj1brPXrv9cXX2zShg27VVRUYj3WtWs7ffXVc7rqqo511ofGJiwursY3zIYy4NVX1W/yZOsN81xJiY58+60OrFypwuPHrXbJmZkatXy5wuLivNVVI1A730b9fIOfn3Tnna5htaTEflX14EHXcJqaKv3yl/aAW1eio92HVV9m9JCAyZNTNHBgS2t/+vQDmjgxS7m59qAYGuqviRMv05NPXi5J6t07StOmddGoUd/V+twVFRVavz5PH398VMuWndC6dXkqLXX95Obv76fbbmupF1+8QsnJoZKkRx+9TF99dVKLFh2rdR983eTJwzRwYLq1P336Ek2c+JZyc+1BNTQ0SBMn3qknn7xLktS7d4qmTXtIo0a9XOtz+/vbNGvWbxUVFSZJKi0t0+TJ7+qVVxaooMDxihEeHqLx4wdr8uRhCgjwV0REqGbNekRdu45VWdm5WvejMTu1d68Or1+vw+vXS5JueuGFBjt3j3vvVfpDD1n7Oz/+WIvGjNGZw4clSf6Bger1wAO65aWX5B8YqOj27TV07ly9dcMNDdZHk1E730b9zJWZKXV0uuaxYYO0dKkjqAYGStdeK11/vX0/MVG69Vbpk0/qth9lZdLRo9Lhw9KRI1LnztLll9ftORqasUMC4uODtHv39QoLs2fqWbMOaeTILVW2ffrpFCu0lpdXqEePVdqy5UyD9TUhIUibNl2rli3tVwe/+uqEbrhhbYOdv3reGxIQHx+t3btfU1hYsCRp1qwvNXLk/1bZ9umnh1uhtby8XD16jNOWLdm1Ov/NN/fQZ589be2PG/eaXnllQbXtH3nkNr300n3W/oABT2nx4g216kNtmDokILZjRzVPTtbhdet09tQp6+fJmZka9dVX1n59fi0ZEBys/7d7tyKTkiRJ+5Yt06wbb1RFFeOPe/zqV7ptxgxrf87gwcpaUP3vQWNG7Xwb9fOMN4cEhIdLY8dKzZrZ9zdvrj6I9uvnCK0VFdJrr9kDZm116CAVFtqfy7ksgwdL3bvbtxkSUMfGjWtnhdWCgjKNH7+j2rbPPLNbBw6clSTZbH567LHLGqSP5x05UqwpUxwvENddFyN//2r/zpuEceMGW2G1oKBI48e/Xm3bZ575lw4csF+RttlseuyxO2t9/uuv72pt5+bm6+9/X+S2/ZQpC5Wb6xg/e+21nWvdh8boxK5d2vP55y5vmA2t+8iR1htmRXm5Fj38cJVvmJL07Ztv6uA331j7fSdObJA+moja+TbqZ76rr3aE1ZIS6bPPqm+7YoWUl2ff9vOTMjLqpg979tivqDbG+4eNDaxDhsRb23Pn/uD2ZqrS0grNnHnQ2h84sKUCAxs2MG7YcNra9vf3U1xcswY9v2mGDHHcsDZ37kq3N1OVlpZp5swvrP2BA9MVGFi70SpxcY67ZnfsOKhz59z/6y0rO6edOx2/Q87HwyypQ4ZY29nLl+v4zp1u26+fPt3abt2njyISEuqtb3CP2vk26ufeFVc4trdvd38zVXm5tGmTY79jR8dd/KiakX89KSmh6tgxzNpfsqTm8aCLFzsGekdGBigzM6Ze+ladgADXgOx8Q1ZTk5KSqI4dk6z9JUtq/mrd+ev3yMhQZWZ2ddO6Zvn5jleKoKBAj44JDnZ8yKjr2QpQNwJDQ5Xcr5+1v3vJkhqP2b14sbXtZ7Mp5dZb66NrqAG1823Uz72YGMn5nrjdu2s+xrlNUJCUnFzn3WpUjAysaWmuV7dWr86t8ZiNG/NUXOy4inbhc9Q354C8efNpFRQ03Rt20tJcp5BavTqrxmM2btyj4mLHVfQLn+NirVnjOGfXrm0VExPhtn1sbKQ6d25j7a9cub1W50f9iOvcWf7NHB8sDq5eXeMx+UeP6pTThIfxaWn10je4R+18G/Vzr1Ur1/2DB6tu5+zIEfvNUdU9B1wZGVhTU8Ot7eLicuXk1DxJWWlphXJyzjo9R5ib1nWrV69Il6m3/vrXepoN2EekpjqCX3FxqXJyar5CXlpa5tIuNbV1rfrwySdrdPCg/ap7s2aBevnl+6tt6+fnpylTfq1mzexXYrdtO6BPP/XeDVeoXlxqqsv+yT2ezZl7yqldiwueAw2D2vk26uee88xdZWXS6dPVtz2vvNy1XYsWdd+vxsTIwJqcHGJtHzzo+Yy6Bw442jo/R32IiAjQVVdF6YUXOunrr/soNNQ+8dlrrx3Q7NmH6/XcpktOdkxFdj40euL8jVf256jdR83i4lING/YX60aqESNu0LJlz2nAgN6KiYmQzWZTbGykBg26SsuXP69hwzIlSdnZR3XHHc+z4pWhmjt9Z3autFT5R454dFzeAceCHs353s0rqJ1vo37uRUU5tj0Jq+edv/FKkpo3r7PuNEpGzsMaGenolvPiADVxXho1IqJu/9emTeuiBx+sfrWt7OxCPf30Hpebv5qqyMhQazsvr9Dj406fdlwhj4io/QeOlSu3KyPjv/XGG79RRkaq+vXrpn79ulXZtqioRO+/v0oTJrypo0dza31u1A/nJShLzpzxeNnHYqd3kKAI98NDUD+onW+jfu4FOa15U1zs+XHObZs17Xu1a2TkFdawMMcyDUVFnl/pOnvWMW40PLzhsvgPP9intZo374cGO6fJzk9nJcllJamanD3r+JcbHh7spqXnduzI0W23Pat//GOx23Yff/yNpkxZQFg1XGCYY6hP2UWsZ1h61vFhqFl4uJuWqC/UzrdRP/ecw2bZRdxz7dyWwOqekYHVeUqqsjLPJwF2btusWd1Oa7VtW76WLDmmJUuO6d//ti/Nmp9v/02Ljw/Siy+mat++fho6NL6GZ2r8AgMdHzguZrWosjLHh5NmzWr/gcNms+mJJ/5L2dkz9OCD/SXZhwps2LBbS5du1oYNu60bve6663qtXfuS3nnnUZfADbP4BzpmfCi/iHcF57b+vCt4BbXzbdTPPecpqS5mRJlz28a6pGpdMXJIQGGhI+QEB3ueqZ3b5ufX7V36r766X6++6rp6iM0m3XJLnP70p4668spIRUcH6r33usvff7Pee8+z8T2NUWGh40qp81RRNQkOdrwgOk9Ldaneemu8RoywLwdYVFSiJ5+cralTP3XpX1hYsMaMuVVPPz1cQUGBuvvuTCUmxujGG/+nxrlb0fBKCx1DTAKCPf9g4dy2JJ8py7yB2vk26udeqdPoxYCLSFbObUs8/0KySTIysDqHzZAQzz9ynL/xyf4c9T8Panm5tHjxMX355Ql9+mlv/eQnsbLZ/PSPf3TRZ58dd7vYQWPmHDZDQjwPrKGhjkFAtQ2s9957oxVWy8vLNWTIc1UutVpQUKQXXvhQ27cf0IIFkyRJ/fp10/jxt+nFFz+qVR9M1qpbN4/WH5/dv38D9MZzzm94ASGej3MODHWMq/b1N01q59uoX+PkHDYvJrA6XbgmsNbAyMB6/LijagkJQW5auoqPd7Q9caLhwmJxcblGjdqiPXsyFRhoU1RUoIYPT6x0RbapOH7cMcg+IcHzBRzi46Ot7RMnLuI2yypMnOhY3nX+/LVVhlVnCxeu0/z5azR48NWSpHHjBjXqwBoSE6PLf/Yzb3fjohUed8w6ERQRoWZhYSopKKjxuPB4x1CdwhMn6qVvDYXa+Tbq1zg5XYDWxdxb5jQ0WE7DfVEFI8ewZmU5/hG0aNFMISGedbNNG8dXDzt31vwPqS7l5BTpm29yrf1rrmneoOc3SVaWY6aEFi0iFRLi2YeONm0ck9A5L5N6sdq0iXNZaWv+/DUeHffJJ452bdrE1XpqLdS941mui1BEta1+5g6Xdm0ccwPXtJwk6ge1823Uzz3nLB4a6vlVVufpsI57Pgtkk2TkFdbt212/NujePbLG1a4SE4PUsqUjGO3Y0fBfPTgvcNCiReMdXF6T7dtzXPa7d2+v1avdv1AlJsaoZcvm1v6OHZceWJOSXK/q5uR49ipwYbv4+ObKzj56yf0wWfby5XrKr25vTGwIx7a7rkAW3727ju3Y4fYYW0CA4rp0sfaP19DedNTOt1G/xunCsBkfX/NqVxERrldYCazuGXmFde3aPJcpqvr2jXbT2u6661xDyooVJ+u8XzWJinLk/9zcpjl+VZLWrt3lMkVV376dazzmuuu6uOyvWLH1ks9fXOw6ftnTcbTOY2gl15vHYIbc7GyXicjb9u1b4zGJvXqpmdO7wv4VK+qlb3CP2vk26ufeoUOuN155cgH6wjb7m+YoQo8ZGVgLC89p6VLH9fXhwxNrPMa5zZYtp7VvX8MOBvH391OfPs2t/T17PJ8wv7EpLCzW0qVbrP3hw/vVeIxzmy1b9mnfvku/snn4sOs4qd69Uzw6rlevDhc8T8N/6EHNsubPt7Y7Dx3qMt1OVboNH25tF5440ajfNE1H7Xwb9ateaam0z2lV9q5daz6mm9M6NkePSrm5dd6tRsXIwCrJZcWotLRIDRzYstq2PXpEqn9/x/jHmTMP1WvfqvLII8mKjXVcyVu06Jib1o3fzJlfWNtpae01cGB6tW179Oig/v17OR27tFbnPno0V7t2OX4HRo++qdLV0wuFh4fovvtutva3bTvgcvMYzPHtzJnWdlhcnHo98EC1bSOTktR95Ehrf8vs2So/V7dT3sFz1M63UT/3Nm1ybMfHSx07Vt82Pl66/PKqj0XVjA2s8+Yd1YYNjkV2p0/vok6dwiq1S0gI0uzZaQoIsP+vHDxYpGnTDlRqd15FRX/rz8yZVS/TKUnz5/fS1Vc3r7GfNps9rP7pT52sn61YcVKrVp2q8djGbN68/2jDht3W/vTpY9SpU+tK7RISYjR79m8VEGCfkuzgweOaNu3Tap+3omKB9WfmzPHVtnMOzElJsXr//YkKD696KpaIiBB98MFElxkNZsz4vNrnRt1LzszUUxUV1h/nN7oLHdm4UTvmzbP2f/rcc2qTkVGpXVBEhO6YM8daUrKkoEArn3++7jvfxFE730b96s7OndLhw479W2+VYmMrtwsPl26/3bHYwOnT0vr11T/vpEmOP4MH122ffYmRN12dd999W7Vy5dUKCwtQYmKw1qy5RlOnHtDXX5/SuXMVSk+P0tix7azprMrKynX//d+puLj2E74PGtRSgwa11J49hfrss2PatOmMDh4s0unTZWrWzE+tWgWpZ89I3XlnvNq3d8wz9+OPxbrvvu9qff7G4L77pmjlyj8rLCxYiYmxWrPmr5o69VN9/fU2nTtXrvT0FI0dO9Cazqqs7Jzuv/9Va/Wp2vjb3xbonnt+otRU+x2qAwb0VlbWNM2Y8W99802W8vIKFR0dpj59rtDo0Te5TKm1efM+/f3vi2rdh8Zq7M6dat6uncvP/Gyun33HZmVJF6w1nrt/v1694oo66cOS8ePVJiND4fHxCoqI0MilS7XxjTe05/PPVVJQoPi0NF01dqyiL7vMOuazRx5R/tHGeROdp6idb6N+5lu4UBo1yr7MakSENHq0PYweOGAvS2KilJ5uD62SfT73hQulurr43Lat9MtfVv65869Ju3bS739fuc3XX9v/mMrowLpp02nddddmzZmTpvDwAEVFBerxxzvo8ccrty0tLdeYMdu1ZEnd3mbXoUOoHn64Xc0NJW3deka/+MUmff990x2/6mzTpr26664XNGfOBIWHhygqKkyPPz5Ujz8+tFLb0tIyjRnzDy1Z4n6+VE+dPVusW26ZrMWLn1KXLvaR7YmJsXryybvcHvftt3t0661Pq6Sk/hee8FUBQUE1rnQTEFR5CEZVP7tUeTk5mjN4sO5etEhhcXEKCA7WVWPH6qqxYyu1rSgv14pnn9WG11+vs/P7Kmrn26if+X74QZo3TxoyxB5ag4Ol6u5PO3dOWrxY2r276scvhZ9fzVNqVdfGZux37naGd09auPBH9eixSosXH9O5cxVVtlm16pQyMr7R66/nVPn4pZgwYaeWLTvhMltBdTZuzNPDD29Tjx6rKk3J1dQtXLhOPXqM0+LFG3Sumo+Qq1ZtV0bGf+v11z+r03Pn5BxTr17j9fvfz9KBA+7HFGdnH9WECW/q6qsf1ZEj3GzlCw6tW6dp3bpp63vvqay46hkdfti8WbP799eyyZMbuHdwh9r5Nurn3q5d0muv2YNoeTVf+ObkSDNnShs3NmzffJlfRUXVIVCS/PwWV/+gFyQkBKlv32glJQXL399Phw4Vae3aPO3dW39XNAMC/NS1a4Q6dgxVQkKwwsP9VVJSrry8Mu3ff1YbN57WsWOmrqc21dsdcJGQEKO+fTsrKSlW/v42HTp0QmvX7tLevT80yPlTU9uoZ88OiouLUlhYkPLzi/Tjj3nauHGPy2IHJpishd7ugk8JiY5Wcr9+imzdWgEhITpz+LB+2LRJP2699OnR0DConW8zrX62PxgVWxQebv+aPiLCfgXz9Gn7ONdTTfs2l2pNmqRqJyn2qcCKi2VWYIXnCKwAcPFMC6y4OO4Cq/FDAgAAANC0EVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYLQA9w9PbZheAAAA1FL5ZD9vdwG1Mami2oe4wgoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYLcDbHfBUenqKRo26UZmZXZWUFCN/f5sOHTqpNWuy9Pbby7R06eZ6O7fNZtOAAb10553XKj09RQkJMQoPD1Z+fpGOHDmp9et368MP/6OFC9epvLy83vrhy7xZP0m6/fZrdPvt16hPn05q1aq5/P1tOnz4pL77Llvvvfe1Pv74G5WWltVrHxqT6Pbtldi7txLT0+3/7dVLQZGR1uMvJycrd//+BulLXGqquo8apQ4336zI1q0VGBqq04cO6YdNm/TdO+8oa/58VVRUNEhffAG1823erF/zdu00Pjv7oo/7YuJErfzzn+u+Qz6I+l064wNrSEiQXnzxV3rooQGVHrviilBdcUVrjRz5U3344X90//1TdOpUfp2ePy2tvWbNekRXXtm+0mPR0eGKjg5X585tdc89P9HWrfs1cuT/auPGPXXaB1/m7fqlprbRnDkTlJZWuX4pKYlKSUnUkCEZ2rRpr0aMeElbtzbMG7WvShkwQLfPmqXQ2Fhvd0V+NpsyJ03S9U88IVuA60tZbEqKYlNS1GXoUGUvX66PRoxQXk6Ol3pqBmrn20yqHy4e9as9owOrzWbT3LmPaeDAdOtnBQVF2r49R2Vl59S5cxtFRYVJku64I0Pt2sXp+usf19mzxXVy/p49O+jLL/9onUOSzp4t1rZtB5SbW6Do6HB16dJWwcHNJEldu7bTV189pxtvfFJr1+6qkz74Mm/XLz09RV988awiI0Otn+Xnn9V33+1XcXGpUlISlZRkf/Ho3v0yff31n9S372Patu1AnZy/MQqLizPmBXfAq68q/aGHrP1zJSX6cds2lRYUqMUVVyi0RQtJUnJmpkYtX643rr5aBceOeau7XkftfJtJ9Ttv95IlHrU7tXdvPffEfNSv9owOrJMnD3MJO9OnL9HEiW8pN7dAkhQaGqSJE+/Uk0/eJUnq3TtF06Y9pFGjXq71uf39bZo167dWoCotLdPkye/qlVcWqKCgyGoXHh6i8eMHa/LkYQoI8FdERKhmzXpEXbuOVVnZuVr3w5d5s37R0eH66KMnrLBaUlKqxx77p6ZOXaSSEsdX/z/9aZpee22sLrssXs2bh2vRosnq0mWMS41RtVN79+rw+vU6vH69JOmmF15osHP3uPdel8Cz8+OPtWjMGJ05fFiS5B8YqF4PPKBbXnpJ/oGBim7fXkPnztVbN9zQYH00GbXzbd6sn7PZ/ft75by+jvpdGmMDa3x8tH73u59b+7NmfakHH/y7S5vCwmJNmvSOJFmhZ8SIG/TSSx9ry5bsWp3/pz9NU5cuba39Rx99U6+8sqBSu/z8s3r22X+poKBIL710nySpU6fWuumm7lq8eEOt+uDLvF2/Rx+93bp6KknDh7+oDz5YVand0qWbdd11j+nbb/+mli2bq127lvrd736up59+r1bnb6xyVq/W27fcosPr1unsqVPWz5MzMxusDwHBwbrhmWes/X3Llulfd9yhCqfx4+dKS7X21VdVWlio22bMsPexXz91GjRIWQsq/ztuCqidbzOhfrh01K/2jJ0lYNy4wQoLC5Zk/xp5/PjXq237zDP/0oED9q+LbDabHnvszlqf//rru1rbubn5+vvfF7ltP2XKQuXmOsZfXntt51r3wZd5u36jR99kbS9ZsqHKsHre4cMnNXnyu9b+7353uzXMA65O7NqlPZ9/7vKC29C6jxypyKQkSVJFebkWPfywS+Bx9u2bb+rgN99Y+30nTmyQPpqI2vk2E+qHS0f9as/YwDpkyDXW9ty5K93ejFNaWqaZM7+w9gcOTFdgYO0uHsfFOe7a27HjoM6dc3/3f1nZOe3cebDK45sib9avS5e2atUq2tqfM2dFjce8994Ka4aHyMhQ/exnPS/5/KhfqUOGWNvZy5fr+M6dbtuvnz7d2m7dp48iEhLqrW9wj9oBuFRGBtaUlER17Jhk7S9ZUvNX685fv0dGhiozs6ub1jXLz3eMYQwKCvToGOercnV9t7sv8Xb9kpNbuex7MrwgN7dABw+esPZvu63PJZ8f9ScwNFTJ/fpZ+57cNLB78WJr289mU8qtt9ZH11ADagegNowMrBdOQbR6dVaNx2zcuEfFxaXVPsfFWrPGcc6uXdsqJibCbfvY2Eh17tzG2l+5cnutzu/LvF2/qKhQl/0zZ856dNzp04XWdo8el13y+VF/4jp3ln8zxwfDg6tX13hM/tGjOrVvn7Ufn5ZWL32De9QOQG0YedNVaqoj+BUXlyonp+bpTEpLy5STc0yXX574f8/RulZ9+OSTNTp48Lhat26hZs0C9fLL9+uee16qsq2fn5+mTPm1mjWzX4ndtu2APv206d5w5e36XRhQIyJCPDrOefqrTp2SZLPZWAjCMHGpqS77J/d4NufxqT17FN3e/iGoxQXPgYZB7Rqf2//5T/tQjcRE+dlsOnvypE58/732L1+uzW+/bcx0SKiar9XPyCusycktre2DB497fNz5G3fsz9HKTcuaFReXatiwv1g3Uo0YcYOWLXtOAwb0VkxMhGw2m2JjIzVo0FVavvx5DRtmv9MvO/uo7rjj+SYddLxdP+ev9iWpW7d2NR4TFRWm1q0dswoEBzdTQkK0myPgDc2Tk63tc6Wlyj9yxKPj8g445tZ1fg40HGrX+KTdc49iO3ZUs/BwBYaGKrJ1a7W/4Qb1e+op/SYrSz+fOVPNwsO93U1Uw9fqZ+QVVucrXXl5hW5aujp92nFlzdOrau6sXLldGRn/rTfe+I0yMlLVr1839evXrcq2RUUlev/9VZow4U0dPZpb63P7Mm/Xb/PmfcrLK7Dm0L3rruv19tvL3B7zi19cJ5vN9fNbXfwOoW45L2FYcuaMx8t2Fp8+7XiOCPfDe1A/qF3jc+bIEeXu26fSwkIFN2+uFqmpahZmf921BQSo+6hRan3NNXorM1P5R496ube4kK/Vz8jAen46JMkeBD3lvEJSeHiwm5ae27EjR7fd9qyeeeaXevDB6ifZ/fjjbzRlyoImH1Yl79evvLxcc+assOo1YEBv3X77Nfroo6rHzMXHR+sPf7i70s8JrOYJDHOsOldW5PniDqVnHR+GTLpi0JRQO99XUVGhg2vW6NsZM7Rr0SJrsYfzbAEBSunfXz/54x/Vqpv94k6LTp00bMECzcjIUHlZWVVPiwbi6/UzckhAYKC/tX0xq0WVlTm+hm/WrPZZ3Gaz6Ykn/kvZ2TOs8FNcXKoNG3Zr6dLN2rBht3Wj0F13Xa+1a1/SO+886hLYmiIT6venP32gM2ccV3fnzJmg3/xmUKXn7devm1au/LPi4yt//V/bqdFQ9/wDHTN2XMyLp3Nb5xt/0HCone/LO3BAb/Tpow2vv14p7Ej2WmUtWKDXr7pK3zvN8JCUnq6eo0c3ZFdRBV+vn5GBtbDQcaXtYiZwDw52vCA6T0t1qd56a7yefXaEwsKCVVRUogkT3lRMzDD17v2Ibrzxf9S79yOKjb1bjz32lhVc7747UwsXTpK/v5F/tQ3ChPrt3/+jRo36mxWYg4IC9corv9bx4+9o1aoX9OWXf1ROzkwtW/acOnRIUHl5eaXFIfLyCmrVB9S90kLHh5CAYM8/GDq3LclvulPOeRO1azrKior0wS9+4fI18tXjxnmxR7gYptbPyEtIzmElJMTzwBMaGlTlc1yKe++9USNG2NeuLi8v15Ahz1W51GpBQZFeeOFDbd9+QAsWTJJkv2o3fvxtevHFj2rVB19lQv0kad68/6h//6c0a9YjSkiIkSRFRIQqI8P1TuO8vAKNGfMPlZaWacwYxzyPubmNN7C26tbNo/WrTVtr2jmwBIR4PmQjMNQxrtrXQw+1822+Wr+LVXzmjNZPm6Z+Tz0lyT5LRFSbNsrLyfFux2qJ+nmvfkYG1uPHHYPszwcNTzh/rXvixGk3LWs2caJjedD589dWGVadLVy4TvPnr9HgwVdLksaNG9RkA6sJ9Tvviy826bLL7tfIkT/Rrbemq3v39mrRIlLFxaXKzv5R8+ev1WuvLdGhQydcloQtKipp1OORQ2JidPnPfubtbly0wuOOWSeCIiLULCxMJQU1f7AIj493PMeJE25amo/a+TZfrd+l2LdsmRV4JCm2Y0efD6zUj8DqIivLscRpixaRCgkJcrkhpzpt2rSwtp2XSb1YbdrEuazUNH/+Go+O++QTR2Bt0yZOycmtlJ3t/TvrGpq363ehoqISTZ++RNOnu19ZJz09xdretGnvRY2/RcM4nuW6CEVU27Y6tmNHjcdFtXHMDVzTcqCoH9Su6blw6rLQFi2qaQkTmVY/IwPr9u2uCb579/Zavdr9C1ViYoxatmxu7e/YcemBJynJ9apgTo5nc4le2C4+vnmTDKzert+lsNlsuvZax1CBNWt2Nej5G1r28uV6ys/P2924aMe2u64gF9+9e42hxxYQoLguXaz94x6EJJNRO9/mq/W7FM7DOSTXGR98FfXzHiPvDFq7dpfLFbm+fTvXeMx113Vx2V+xYusln7+42PUOVk/HYTqPwZRcbz5qSrxdv0tx003dXYYkvPvu8gY9PzyTm53tMpF82759azwmsVcva25BSdq/YkW99A3uUbumJ66z62t/gQFzecJzptXPyMBaWFispUu3WPvDh/er8RjnNlu27NO+fZf+F3v4sOs4qd69U6pp6apXrw4XPM/JS+6DL/N2/S7F738/1OX8a9c27iusvixr/nxru/PQoS7TJVWl2/Dh1nbhiROEHi+idk1L12HDrO3SwkId+fZbL/YGF8u0+hkZWCVp5swvrO20tPYaODC92rY9enRQ//69nI5dWqtzHz2aq127Dln7o0ffVOnq6YXCw0N03303W/vbth1wufmoqfFm/S7Www8P0PXXd7X2J078Z4OeHxfn25kzre2wuDj1euCBattGJiWp+8iR1v6W2bNVfo6xyd5C7ZqOToMHq9OgQdb+959+qnMlni8kA+8ysX7GBtZ58/6jDRt2W/vTp49Rp06tK7VLSIjR7Nm/VUCAfbL6gwePa9q0T6t93oqKBdafmTPHV9vOOXAlJcXq/fcnKjy86qlYIiJC9MEHE13uiJ8x4/Nqn7sp8Hb9JGnw4KvVokVktY/7+flp/Pjb9Morv7Z+9s47X9U4IwTqXnJmpp6qqLD+OAeVCx3ZuFE75s2z9n/63HNqk5FRqV1QRITumDPHWhK0pKBAK59/vu4738RRO9/maf3iOnfW4DfeUGzHjjU+55XDh+uOd9+19svLyrRs8uQ66zMcmlL9jLzp6rz77puilSv/rLCwYCUmxmrNmr9q6tRP9fXX23TuXLnS01M0duxAa+xhWdk53X//q9Yk/rXxt78t0D33/ESpqfY7VAcM6K2srGmaMePf+uabLOXlFSo6Okx9+lyh0aNvchn/uHnzvkqT0DdF3qyfJD30UH/NnfuYli7drGXLtmj79hzl5RUqJiZcV16ZrGHDMtWlS1ur/YoVW/XrX79aJ+duzMbu3Knm7dq5/MzP5vrZd2xWlnTBWvG5+/fr1SuuqJM+LBk/Xm0yMhQeH6+giAiNXLpUG994Q3s+/1wlBQWKT0vTVWPHKvqyy6xjPnvkESPWw/YmaufbvFk//8BA9Rw9Wj1Hj9bhDRu078svdXTLFhUcParSwkIFRUWp1ZVXqsvQoUro2dPl2CXjx1e66a4pon61Y3Rg3bRpr+666wXNmTNB4eEhiooK0+OPD9Xjjw+t1NY+6fs/tGRJ3VwdO3u2WLfcMlmLFz9lhZrExFg9+eRdbo/79ts9uvXWp1VSwprJ3qzfeUFBgRowoLcGDOjttt1HH63WiBEvNdkb5S5GQFBQjSsVBQRVHkJT1c8uVV5OjuYMHqy7Fy1SWFycAoKDddXYsbpq7NhKbSvKy7Xi2We14fXX6+z8vora+TYT6ifZb4ZL7NWrxnalhYX67Le/1frp0+v0/L6K+tWOsUMCzlu4cJ169BinxYs36Fw145dWrdqujIz/1uuvf1an587JOaZevcbr97+fpQMHjrltm519VBMmvKmrr35UR440zZutquLN+i1dull79hyp9vHy8nKtXr1TP//5HzVkyHMqKKj96lpoOIfWrdO0bt209b33VFZc9QeNHzZv1uz+/Y34OgsO1M73nDlyRJveeksndtV8Q2rx6dNaN3Wqpl15pTFhp6lrDPXzq7jg0rPLg36Dqn/QCxISYtS3b2clJcXK39+mQ4dOaO3aXdq794cGOX9qahv17NlBcXFRCgsLUn5+kX78MU8bN+5xmSwfVfNW/dq2jVNaWnu1atVcLVpE6syZszpy5JTWrMnSoUNmrpwzWQu93QWfEhIdreR+/RTZurUCQkJ05vBh/bBpk37c2rDTo+HiUTvfExobq1ZpaQqLi1NoixYKiopSaWGhzp48qR+/+05Ht2zhBjmDmVy/pyoqqp3k1qcCK9BUEFgBAE2Nu8Bq/JAAAAAANG0EVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAoxFYAQAAYDQCKwAAAIxGYAUAAIDRCKwAAAAwGoEVAAAARiOwAgAAwGgEVgAAABiNwAoAAACjEVgBAABgNAIrAAAAjEZgBQAAgNEIrAAAADAagRUAAABGI7ACAADAaARWAAAAGI3ACgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEYjsAIAAMBoBFYAAAAYjcAKAAAAo/lVVFR4uw8AAABAtbjCCgAAAKMRWAEAAGA0AisAAACMRmAFAACA0QisAAAAMBqBFQAAAEb7/0eljsnPjbU8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "error= np.load('./data_produced/Prepared_chi2.npy',allow_pickle=True)\n",
    "\n",
    "lambdas = np.load('./data_produced/Prepared_parameters.npy',allow_pickle=True)\n",
    "\n",
    "gHBfix=[]\n",
    "for i in range(12):\n",
    "    gHBfix.append(0.0)\n",
    "\n",
    "gHBfix[1]=1.0\n",
    "gHBfix[6]=-0.5\n",
    "gHBfix[8]=-0.5\n",
    "gHBfix=np.array(gHBfix)   \n",
    "data=lambdas /1.666666666666666\n",
    "data=data+gHBfix\n",
    "cmap = plt.cm.get_cmap('gist_gray')\n",
    "\n",
    "for l,lambdas in enumerate(data):\n",
    "    newlambdas = reorder(data[l]).reshape(-1,12)\n",
    "    l=newlambdas[0].reshape(2,6)\n",
    "    plt.matshow(l,cmap='seismic_r',vmin=-0.5,vmax=0.5)\n",
    "    for (i, j), z in np.ndenumerate(l):\n",
    "        c='white' \n",
    "        plt.text(j, i, '{:0.1f}'.format(z), ha='center', va='center', fontsize=35, color=c)\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}