cross entropy beamforming-partial connected.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Load necessary module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io as io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import bernoulli\n",
    "from scipy.stats import multinomial\n",
    "from collections import OrderedDict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Load model that defined in channel realization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "mat = io.loadmat('tempNs4.mat')\n",
    "model = {}\n",
    "model['H'] = mat['H']\n",
    "model['Wopt'] = mat['Wopt']\n",
    "model['Fopt'] = mat['Fopt']\n",
    "model['Ns'] = mat['Ns'][0,0] ## need to change nb of stream in channel realization, this case 1\n",
    "# model['Ns'] = 2\n",
    "model['Nt'] = mat['Nt'][0,0]\n",
    "model['Nr'] = mat['Nr'][0,0]\n",
    "model['Nc'] = mat['Nc'][0,0]\n",
    "model['Nray'] = mat['Nray'][0,0]\n",
    "model['realization'] = mat['realization'][0,0]\n",
    "model['At'] = mat['At']\n",
    "model['Ar'] = mat['Ar']\n",
    "model['Nrf'] = 4\n",
    "archi = 'FC'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H\n",
      "Wopt\n",
      "Fopt\n",
      "Ns\n",
      "Nt\n",
      "Nr\n",
      "Nc\n",
      "Nray\n",
      "realization\n",
      "At\n",
      "Ar\n",
      "Nrf\n"
     ]
    }
   ],
   "source": [
    "for key in model:\n",
    "    print(key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model['Ns']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load values of phase shifter included:\n",
    "- nb_bit: number of bit\n",
    "- nb_ps: number of phase shifter in the system\n",
    "- nb_state: number of state in each phase shifter (2^nb_bit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross Entropy Beamforming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Helper function for Cross entropy Beamforming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phase_shifter(n_bit,n_ant):\n",
    "    \"\"\"\n",
    "    Define phase shifter with N bits resolution\n",
    "    \n",
    "    Arguments:\n",
    "    N -- number of bits resolutions\n",
    "    \n",
    "    Returns:\n",
    "    y -- python dictionnary of phase shifter: index of state and value of phase\n",
    "    \"\"\"\n",
    "    nb_state = 2**n_bit\n",
    "    w = np.exp(1j*2*np.pi/nb_state)\n",
    "    ps = np.zeros(nb_state,dtype=complex)\n",
    "    for i in range(nb_state):\n",
    "        ps[i] = w**i/np.sqrt(n_ant)\n",
    "    return np.around(ps, decimals=3)\n",
    "\n",
    "def generate_candidat(S,model,phase_value,ps_system_proba,archi):\n",
    "    \"\"\"\n",
    "    generate S candidates from the previous pmf\n",
    "    \n",
    "    Arguments: \n",
    "    S -- Number of candidates to generate\n",
    "    model -- channel model\n",
    "    phase_value -- possible phase shifted value\n",
    "    ps_system_proba -- pmf of phase shifters\n",
    "    \n",
    "    Returns\n",
    "    F_rf -- S analog precoder system\n",
    "    \"\"\"\n",
    "#     print(ps_system_proba)\n",
    "    F_rf = np.zeros((S,model['Nt'],model['Nrf']),dtype=complex)\n",
    "    if (archi == 'FC' or archi == 'FCOMP'):\n",
    "        for i in range(S):\n",
    "            for j in range(model['Nt']):\n",
    "                for k in range(model['Nrf']):\n",
    "                    F_rf[i,j,k] = np.random.choice(phase_value, 1, p=ps_system_proba[j,k,:])\n",
    "    if archi == 'PC':\n",
    "        only_PS = np.zeros((S,model['Nt']),dtype=complex)\n",
    "        m = int(model['Nt']/model['Nrf'])\n",
    "        for i in range(S):\n",
    "            for j in range(model['Nt']):\n",
    "                only_PS[i,j] = np.random.choice(phase_value, 1, p=ps_system_proba[j,:])\n",
    "            for k in range(model['Nrf']):\n",
    "                F_rf[i,k*m:(k+1)*m,k] = only_PS[i,k*m:(k+1)*m]\n",
    "    return F_rf\n",
    "\n",
    "def compute_Fbb(F_rf,model,channel):\n",
    "    \"\"\"\n",
    "    Compute baseband precoder base on the effective channel (Heq = H*Frf)\n",
    "    \n",
    "    Arguments:\n",
    "    F_rf -- analog precoder\n",
    "    model -- channel model\n",
    "    \n",
    "    Returns:\n",
    "    F_bb -- S baseband precoder\n",
    "    \"\"\"\n",
    "    S = F_rf.shape[0]\n",
    "    channel_eq = np.zeros((S,model['Nr'],model['Nrf']),dtype=complex)\n",
    "    F_bb = np.zeros((S,model['Nrf'],model['Ns']),dtype=complex)\n",
    "#     channel_index = 0 ### need to adjust \n",
    "    for i in range(S):\n",
    "        channel_eq[i,:,:] =  channel@F_rf[i,:,:]\n",
    "        U, s, V = np.linalg.svd(channel_eq[i,:,:], full_matrices=True)\n",
    "        V = V.conj().T\n",
    "        F_bb[i,:,:] = V[0:model['Nt'],0:model['Ns']]\n",
    "        F_bb[i,:,:] = np.sqrt(model['Ns'])*F_bb[i,:,:]/np.linalg.norm(F_rf[i,:,:]@F_bb[i,:,:],'fro')\n",
    "    return F_bb\n",
    "\n",
    "def optimal_beamformer(model, channel):\n",
    "    \"\"\"\n",
    "    Compute the optimal beamformer by SVD\n",
    "    \n",
    "    Arguments:\n",
    "    model -- channel model\n",
    "    channel_index -- index in 1000 realization channel\n",
    "    \n",
    "    Returns:\n",
    "    Fopt -- optimal precoder\n",
    "    Wopt -- optimal combiner\n",
    "    \"\"\"\n",
    "    Wopt = np.zeros((model['Nr'],model['Ns']),dtype=complex)\n",
    "    Fopt = np.zeros((model['Nt'],model['Ns']),dtype=complex)\n",
    "    U, s, V = np.linalg.svd(channel, full_matrices=True)\n",
    "    V = V.conj().T\n",
    "    Fopt = V[0:model['Nt'],0:model['Ns']]\n",
    "    Wopt = U[0:model['Nr'],0:model['Ns']]\n",
    "    return Fopt, Wopt\n",
    "\n",
    "\n",
    "def update_pmf(F_elites,phase_value,nb_elite,R_elite,ps_system_proba,archi):\n",
    "    \"\"\"\n",
    "    Update new phase shifters's pmf base on previous candidates's results\n",
    "    \n",
    "    Arguments:\n",
    "    F_elites -- elites analog precoder\n",
    "    phase_value -- phase shift value\n",
    "    nb_elite -- number of elites \n",
    "    \n",
    "    Returns:\n",
    "    ps_system_proba -- pmf of phase shifter\n",
    "    \"\"\"\n",
    "    alpha = 0.9\n",
    "    nb_state = phase_value.shape[0]\n",
    "    stat_count = np.zeros(nb_state)\n",
    "    assign = np.zeros(nb_elite)\n",
    "    weight = R_elite/np.mean(R_elite)\n",
    "#     print(weight.shape)\n",
    "#     print(assign.shape)\n",
    "#     print(F_elites.shape)\n",
    "    if (archi == 'FC' or archi == 'FCOMP'):\n",
    "        for i in range(model['Nt']):\n",
    "            for j in range(model['Nrf']):\n",
    "                for k in range(nb_state):\n",
    "                    assign = F_elites[:,i,j] == phase_value[k]\n",
    "#                     print(assign.shape)\n",
    "#                     print(weight.shape)\n",
    "                    proba_est = np.sum(assign*weight)/nb_elite\n",
    "#                     proba_est = np.sum(assign)/nb_elite\n",
    "                    ps_system_proba[i,j,k] = alpha*proba_est + (1-alpha)*ps_system_proba[i,j,k] \n",
    "    if archi == 'PC':\n",
    "        # Convert to PC form\n",
    "        F_RF_conv = np.zeros((F_elites.shape[0],model['Nt']),dtype=complex)\n",
    "        \n",
    "        m = int(model['Nt']/model['Nrf'])\n",
    "        for i in range(F_elites.shape[0]):\n",
    "            f_rf = F_elites[i,:]\n",
    "            f_rf_conv = f_rf[0:m,0]\n",
    "            for j in range(model['Nrf']-1):\n",
    "                f_rf_conv = np.concatenate((f_rf_conv,f_rf[(j+1)*m:(j+2)*m,j+1]))\n",
    "            F_RF_conv[i,:] = f_rf_conv\n",
    "    \n",
    "    \n",
    "        for i in range(model['Nt']):\n",
    "            for k in range(nb_state):\n",
    "                assign = F_RF_conv[:,i] == phase_value[k]\n",
    "                proba_est = np.sum(assign*weight)/nb_elite\n",
    "#                 proba_est = np.sum(assign)/nb_elite\n",
    "                ps_system_proba[i,k] = alpha*proba_est + (1-alpha)*ps_system_proba[i,k]   \n",
    "    return ps_system_proba\n",
    "\n",
    "def data_rate(Ns,SNR,W,H,F):\n",
    "    \"\"\"\n",
    "    Compute the data rate of the beamforming system\n",
    "    \n",
    "    Arguments:\n",
    "    Ns -- number of stream\n",
    "    SNR -- Signal to noise ratio\n",
    "    W -- combiner\n",
    "    H -- channel\n",
    "    F -- precoder\n",
    "    \n",
    "    Returns:\n",
    "    R -- data rate\n",
    "    \"\"\"\n",
    "    R = np.log2(np.linalg.det(np.eye(Ns)+SNR/Ns*np.linalg.pinv(W)@H@F@F.conj().T@H.conj().T@W))\n",
    "    return R.real"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Phase shiter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_bit = 2\n",
    "nb_ps = model['Nt']*model['Nrf']\n",
    "nb_ps\n",
    "nb_state = 2**nb_bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.125+0.   j,  0.   +0.125j, -0.125+0.   j, -0.   -0.125j])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phase_value = phase_shifter(nb_bit,model['Nt'])\n",
    "phase_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Phase shifter system dimension is: NtxNrfxN_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross entropy beamforming Algorithm\n",
    " Step 1: Randomly generate S candidate analog beamformers\n",
    "\n",
    " Step 2: Compute S corresponding digital precoders Fbb based on the effective channel Heq = HFrf\n",
    "\n",
    " Step 3: Calculate the achievable sum-rate\n",
    "\n",
    " Step 4: Sort R in a descend order\n",
    "\n",
    " Step 5: Select elites as Frf_1, Frf_2, ..., Frf_elite\n",
    "\n",
    " Step 6: Update new pmf for each phase shifter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ps_system_proba = np.ones((model['Nt'],nb_state))*(1/nb_state)\n",
    "# F_rf = generate_candidat(S,model,phase_value,ps_system_proba,'PC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# F_RF_conv[0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# F_RF_conv = np.zeros((F_rf.shape[0],model['Nt']),dtype=complex)\n",
    "# m = int(model['Nt']/model['Nrf'])\n",
    "# for i in range(F_rf.shape[0]):\n",
    "#     f_rf = F_rf[i,:]\n",
    "#     f_rf_conv = f_rf[0:m,0]\n",
    "#     for j in range(model['Nrf']-1):\n",
    "#         f_rf_conv = np.concatenate((f_rf_conv,f_rf[(j+1)*m:(j+2)*m,j+1]))\n",
    "#     F_RF_conv[i,:] = f_rf_conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test = F_rf[0,:]\n",
    "# m = int(model['Nt']/model['Nrf'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test2 = test[0:m,0]\n",
    "# print(test2)\n",
    "# for i in range(model['Nrf']-1):\n",
    "#     print(i+1)\n",
    "#     test2 = np.concatenate((test2,test[(i+1)*m:(i+2)*m,i+1]))\n",
    "# #     print(test[(i+1)*m:(i+2)*m,i+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def OMP(_Fopt,_NRF,_At):\n",
    "    \"\"\"\n",
    "    Orthogonal matching pursuit for hybrid beamforming\n",
    "    \n",
    "    \"\"\"\n",
    "    _FRF = np.empty((_At.shape[0],0))\n",
    "    _Fres = _Fopt\n",
    "    for k in range(_NRF):\n",
    "        PU = _At.conj().T@_Fres\n",
    "        ind = np.argmax(np.sum((np.absolute(PU))**2,axis=1))\n",
    "        _FRF = np.append(_FRF,_At[:,[ind]],axis=1)\n",
    "        _FBB = (np.linalg.pinv(_FRF))@_Fopt\n",
    "        _Fres = (_Fopt-_FRF@_FBB)/np.linalg.norm((_Fopt-_FRF@_FBB),'fro')\n",
    "#         _Fres = (_Fopt-_FRF@_FBB)\n",
    "    _FBB = np.sqrt(model['Ns'])*_FBB/np.linalg.norm(_FRF@_FBB,'fro')\n",
    "    return _FRF,_FBB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantize_Frf(model,n_bit,F_rf,channel_index):\n",
    "    \"\"\"\n",
    "    Quantize precoder for low resolutioin phase shifter\n",
    "    \n",
    "    \"\"\"\n",
    "    q = 2*np.pi/2**n_bit\n",
    "    amplitude = 1/np.sqrt(model['Nt'])\n",
    "    Frf_angle = np.angle(F_rf)\n",
    "    angle_quantized = q * np.round(Frf_angle/q)\n",
    "    Frf_quantized = amplitude*np.exp(1j*angle_quantized)\n",
    "    Frf_quantized = np.around(Frf_quantized, decimals=3)\n",
    "    \n",
    "    channel_eq =  model['H'][:,:,channel_index]@Frf_quantized\n",
    "    U, s, V = np.linalg.svd(channel_eq, full_matrices=True)\n",
    "    V = V.conj().T\n",
    "    F_bb = V[0:model['Nt'],0:model['Ns']]\n",
    "    F_bb = np.sqrt(model['Ns'])*F_bb/np.linalg.norm(Frf_quantized@F_bb,'fro')\n",
    "    \n",
    "    return Frf_quantized, F_bb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantize_Frf(model,n_bit,F_rf,channel_index):\n",
    "    \"\"\"\n",
    "    Quantize precoder for low resolutioin phase shifter\n",
    "    \n",
    "    \"\"\"\n",
    "    q = 2*np.pi/2**n_bit\n",
    "    amplitude = 1/np.sqrt(model['Nt'])\n",
    "    Frf_angle = np.angle(F_rf)\n",
    "    angle_quantized = q * np.round(Frf_angle/q)\n",
    "    Frf_quantized = amplitude*np.exp(1j*angle_quantized)\n",
    "    Frf_quantized = np.around(Frf_quantized, decimals=3)\n",
    "    \n",
    "    channel_eq =  model['H'][:,:,channel_index]@Frf_quantized\n",
    "    U, s, V = np.linalg.svd(channel_eq, full_matrices=True)\n",
    "    V = V.conj().T\n",
    "    F_bb = V[0:model['Nt'],0:model['Ns']]\n",
    "    F_bb = np.sqrt(model['Ns'])*F_bb/np.linalg.norm(Frf_quantized@F_bb,'fro')\n",
    "    \n",
    "    return Frf_quantized, F_bb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_pmf_omp(model,nb_bit,Frf_omp_quantize,phase_value,ratio):\n",
    "    \"\"\"\n",
    "    Initialize pmf of phase shifter with solution getting by OMP\n",
    "    \n",
    "    \"\"\"\n",
    "    low_proba = 1/(2**nb_bit-1+ratio)\n",
    "    high_proba = ratio*low_proba\n",
    "    proba_init = np.zeros((model['Nt'],model['Nrf'],2**nb_bit))\n",
    "    for i in range(phase_value.shape[0]):\n",
    "        assign = Frf_omp_quantize == phase_value[i]\n",
    "        proba_init[:,:,i] = assign*high_proba + np.invert(assign)*low_proba\n",
    "    return proba_init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'channel_index' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-7a92592e2cb4>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mmodel\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Fopt'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mchannel_index\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'channel_index' is not defined"
     ]
    }
   ],
   "source": [
    "model['Fopt'][:,:,channel_index].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model['Nrf']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel_index = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timeit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.21632180438771087\n"
     ]
    }
   ],
   "source": [
    "start = timeit.default_timer()\n",
    "\n",
    "FRF_omp, FBB_omp = OMP( model['Fopt'][:,:,channel_index], model['Nrf'], model['At'][:,:,channel_index] )# ps_system_proba_omp = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "Frf_quantize, Fbb_quantize = quantize_Frf(model,nb_bit,FRF_omp,channel_index)\n",
    "ps_system_proba_omp = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "stop = timeit.default_timer()\n",
    "print(stop - start) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.zeros((S))\n",
    "Fopt, Wopt = model['Fopt'][:,:,channel_index],model['Wopt'][:,:,channel_index]\n",
    "FRF_omp, FBB_omp = OMP( Fopt, model['Nrf'], model['At'][:,:,channel_index] )# ps_system_proba_omp = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "Frf_quantize, Fbb_quantize = quantize_Frf(model,nb_bit,FRF_omp,channel_index)\n",
    "ps_system_proba = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "relative_std = 1000 # initialize high value of relative standard deviation\n",
    "i = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_rf = generate_candidat(S,model,phase_value,ps_system_proba,archi)\n",
    "F_bb = compute_Fbb(F_rf,model,model['H'][:,:,channel_index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNR = 10**(10/10)\n",
    "for j in range(S):\n",
    "    R[j] = data_rate(model['Ns'],SNR,Wopt,model['H'][:,:,channel_index],F_rf[j,:,:]@F_bb[j,:,:])\n",
    "R_sort_index = np.argsort(R)[::-1][:nb_elite]\n",
    "if i==0:\n",
    "    R_max = R[R_sort_index[0]]\n",
    "else:\n",
    "    R_max = np.append(R_max,R[R_sort_index[0]])\n",
    "F_elites = F_rf[R_sort_index,:,:]\n",
    "Frf_max = F_rf[R_sort_index[0],:,:]\n",
    "Fbb_max = F_bb[R_sort_index[0],:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test['nb_ite']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ps_system_proba = update_pmf(F_elites,phase_value,nb_elite,R[R_sort_index],ps_system_proba,archi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.84063809776006"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "26.636025176473936/18/0.21632180438771087"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26.636025176473936\n"
     ]
    }
   ],
   "source": [
    "start = timeit.default_timer()\n",
    "\n",
    "test = CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,10**(10/10),'FCOMP')\n",
    "stop = timeit.default_timer()\n",
    "print(stop - start) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,SNR,archi):\n",
    "    \"\"\"\n",
    "    Hybrid beamforming by Cross-entropy method\n",
    "    \n",
    "    \"\"\"\n",
    "    phase_value = phase_shifter(nb_bit,model['Nt'])\n",
    "    Fopt, Wopt = model['Fopt'][:,:,channel_index],model['Wopt'][:,:,channel_index]\n",
    "    if archi == 'FC':\n",
    "        ps_system_proba = np.ones((model['Nt'],model['Nrf'],nb_state))*(1/nb_state)\n",
    "#         proba_evo = np.copy(ps_system_proba[0,0,:])\n",
    "    elif archi == 'PC':\n",
    "        ps_system_proba = np.ones((model['Nt'],nb_state))*(1/nb_state)\n",
    "#         proba_evo = np.copy(ps_system_proba[0,:])\n",
    "#     print(ps_system_proba)\n",
    "    elif archi == 'FCOMP':\n",
    "        FRF_omp, FBB_omp = OMP( Fopt, model['Nrf'], model['At'][:,:,channel_index] )# ps_system_proba_omp = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "        Frf_quantize, Fbb_quantize = quantize_Frf(model,nb_bit,FRF_omp,channel_index)\n",
    "        ps_system_proba = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "\n",
    "    relative_std = 1000 # initialize high value of relative standard deviation\n",
    "    i = 0\n",
    "    while ((i<nb_iterate) and (relative_std>3)):\n",
    "#         print(ps_system_proba)\n",
    "        F_rf = generate_candidat(S,model,phase_value,ps_system_proba,archi)\n",
    "        F_bb = compute_Fbb(F_rf,model,model['H'][:,:,channel_index])\n",
    "#         print(F_bb)\n",
    "        R = np.zeros((S))\n",
    "        for j in range(S):\n",
    "            R[j] = data_rate(model['Ns'],SNR,Wopt,model['H'][:,:,channel_index],F_rf[j,:,:]@F_bb[j,:,:])\n",
    "#         print(R)\n",
    "        R_sort_index = np.argsort(R)[::-1][:nb_elite]\n",
    "        if i==0:\n",
    "            R_max = R[R_sort_index[0]]\n",
    "        else:\n",
    "            R_max = np.append(R_max,R[R_sort_index[0]])\n",
    "        F_elites = F_rf[R_sort_index,:,:]\n",
    "        Frf_max = F_rf[R_sort_index[0],:,:]\n",
    "        Fbb_max = F_bb[R_sort_index[0],:,:]\n",
    "        ps_system_proba = update_pmf(F_elites,phase_value,nb_elite,R[R_sort_index],ps_system_proba,archi)\n",
    "#         if archi == 'FC':\n",
    "#             proba_evo = np.vstack((proba_evo,ps_system_proba[0,0,:]))\n",
    "#         elif archi == 'PC':\n",
    "#             proba_evo = np.vstack((proba_evo,ps_system_proba[0,:]))\n",
    "#         print(R[R_sort_index[0]])\n",
    "        \n",
    "        if i>=5:\n",
    "            R_max_normalize = R_max-np.min(R_max) ## Need to review\n",
    "            relative_std = 100*np.std(R_max_normalize[i-5:i])/np.mean(R_max_normalize[i-5:i])\n",
    "#             print(relative_std)\n",
    "        i +=1\n",
    "    output = {}\n",
    "    output['R_max'] = R_max\n",
    "    output['Frf_max'] = Frf_max\n",
    "    output['Fbb_max'] = Fbb_max\n",
    "    output['nb_ite'] = i\n",
    "#     output['proba_evo'] = proba_evo\n",
    "    output['ps_system_proba'] = ps_system_proba\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_elite = 40\n",
    "S = 200\n",
    "channel_index = 0\n",
    "nb_iterate = 30\n",
    "# if archi == 'FC':\n",
    "#     ps_system_proba = np.ones((model['Nt'],model['Nrf'],nb_state))*(1/nb_state)\n",
    "# elif archi == 'PC':\n",
    "#     ps_system_proba = np.ones((model['Nt'],nb_state))*(1/nb_state)\n",
    "SNR_dB = np.arange(-30,15,5)\n",
    "# SNR = 10**(SNR_dB/10)\n",
    "# SNR = 10**(10/10)\n",
    "# smax = SNR.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Evolution of pmf of 1 phase shifter\n",
    "# np.around(CE_result['ps_system_proba'], decimals=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "rate_CE_PC = np.zeros((SNR_dB.shape[0]))\n",
    "sum_rate_CE_PC = np.zeros((SNR_dB.shape[0]))\n",
    "CE_nb_ite = 0\n",
    "sum_nb_ite = 0\n",
    "for channel_index in range(5):\n",
    "    print('index: ' + str(channel_index))\n",
    "    CE_result = CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,10**(10/10),'PC')\n",
    "    CE_nb_ite = CE_result['nb_ite']\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        rate_CE_PC[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],CE_result['Frf_max']@CE_result['Fbb_max'])\n",
    "    sum_rate_CE_PC = sum_rate_CE_PC + rate_CE_PC\n",
    "    sum_nb_ite = CE_nb_ite + sum_nb_ite\n",
    "    print(sum_nb_ite)\n",
    "    print(rate_CE_PC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_CE_FC = np.zeros((SNR_dB.shape[0]))\n",
    "sum_rate_CE_FC = np.zeros((SNR_dB.shape[0]))\n",
    "CE_nb_ite = 0\n",
    "sum_nb_ite = 0\n",
    "for channel_index in range(2):\n",
    "    print('index: ' + str(channel_index))\n",
    "    CE_result = CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,10**(10/10),'FC')\n",
    "    CE_nb_ite = CE_result['nb_ite']\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        rate_CE_FC[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],CE_result['Frf_max']@CE_result['Fbb_max'])\n",
    "    sum_rate_CE_FC = sum_rate_CE_FC + rate_CE_FC\n",
    "    sum_nb_ite = CE_nb_ite + sum_nb_ite\n",
    "    print(rate_CE_FC)\n",
    "    print(sum_nb_ite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_rate_CE_FC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print('Data rate evolution: ' + str(CE_result['R_max']))\n",
    "# plt.plot(CE_result['R_max'])\n",
    "# plt.title('data rate against iteration - 2 bits PS')\n",
    "# plt.ylabel('data rate (bps/Hz)')\n",
    "# plt.xlabel('iteration')\n",
    "# print('Number of iteration: ' + str(CE_result['nb_ite']))\n",
    "# print('Phase shifter pmf' + str(np.around(CE_result['proba_evo'], decimals=3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test:\n",
    "\n",
    "2 bits -- 17.56 bps/Hz - 16 iterate\n",
    "\n",
    "3 bits -- 17.68 bps/Hz - 18 iterate\n",
    "\n",
    "4 bits -- 17.71 bps/Hz - 18 iterate\n",
    "\n",
    "6 bits -- 17.68 bps/Hz - 19 iterate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Digital Beamforming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model['Ns']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_o = data_rate(model['Ns'],10**(10/10),model['Wopt'][:,:,0],model['H'][:,:,0],model['Fopt'][:,:,0])\n",
    "print('Data rate by digital beamforming: ' + str(R_o))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize condition\n",
    "x = np.linalg.norm(CE_result['Frf_max']@CE_result['Fbb_max'],'fro')\n",
    "print(\"Ns\", model['Ns'] , \"   Frobenius norm FRF*FBB=\", x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digital_rate = np.zeros((SNR_dB.shape[0]))\n",
    "sum_digital_rate = np.zeros((SNR_dB.shape[0]))\n",
    "for channel_index in range(5):\n",
    "    print('index: ' + str(channel_index))\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        digital_rate[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],model['Fopt'][:,:,channel_index])\n",
    "    sum_digital_rate = sum_digital_rate + digital_rate\n",
    "    print(digital_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_CE_OMP = np.zeros((SNR_dB.shape[0]))\n",
    "sum_rate_CE_OMP = np.zeros((SNR_dB.shape[0]))\n",
    "CE_nb_ite = 0\n",
    "sum_nb_ite = 0\n",
    "for channel_index in range(5):\n",
    "    print('index: ' + str(channel_index))\n",
    "    CE_result = CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,10**(10/10),'FCOMP')\n",
    "    CE_nb_ite = CE_result['nb_ite']\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        rate_CE_OMP[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],CE_result['Frf_max']@CE_result['Fbb_max'])\n",
    "    sum_rate_CE_OMP = sum_rate_CE_OMP + rate_CE_OMP\n",
    "    sum_nb_ite = CE_nb_ite + sum_nb_ite\n",
    "    print(rate_CE_OMP)\n",
    "    print(sum_nb_ite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "OMP( Fopt, model['Nrf'], model['At'][:,:,channel_index] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model['At'][:,:,channel_index] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "array([  19.54922669,   58.2922556 ,  158.79208973,  369.28672871,\n",
    "        706.87896155, 1135.99094508, 1609.32176865, 2099.33372726,\n",
    "       2594.96375472])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNR_dB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_rate_CE_OMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "OMP = np.array([ 0.22317937,  0.65307923,  1.7242831 ,  3.89130729,  7.38210255,\n",
    "       12.0446525 , 17.52735967, 23.49692577, 29.73026376])\n",
    "quan_OMP = np.array([ 0.15360772,  0.46233855,  1.28936884,  3.13826287,  6.41532018,\n",
    "       11.08009538, 16.72319897, 22.90968806, 29.3595054 ])\n",
    "sum_rate_CE_PC = np.array([  0.52366061,   1.60358201,   4.64395702,  12.00925971,\n",
    "        26.37698906,  48.61650672,  76.83223108, 108.17246538,\n",
    "       140.75704393])\n",
    "sum_rate_CE_FC = np.array([ 0.32917173,  1.00278803,  2.86268021,  7.16043021, 14.83647874,\n",
    "       25.44730884, 37.70663827, 50.64386273, 63.81787168])\n",
    "sum_digital_rate = np.array([  1.59151382,   4.67524429,  12.36836586,  27.68981093,\n",
    "        51.23774157,  80.39861587, 112.17409253, 144.91729675,\n",
    "       177.98393834])\n",
    "sum_rate_CE_OMP = np.array([  1.05497188,   3.16202751,   8.718125  ,  20.75158987,\n",
    "        41.07679319,  68.25812854,  99.18205758, 131.62253347,\n",
    "       164.58955784])\n",
    "sum_rate_coordiante_descent = np.array([ 0.16609992,  0.50635606,  1.44755711,  3.62588512,  7.50916101,\n",
    "       12.8508677 , 18.9981428 , 25.4733022 , 32.06249374])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXd8zdf/x58nQ5YIEmJEbDGyREJtoURr02r91KxNqbaKttTsQNEaLTGCtkGt2rMJolaQEjNGRIIgJG52cu/5/XGT+01k3UhC8Hk+HveRfM7nc855n09yz/vzOeP1FlJKFBQUFBTeXAxetgEKCgoKCi8XxREoKCgovOEojkBBQUHhDUdxBAoKCgpvOIojUFBQUHjDURyBgoKCwhuO4ggUFBQU3nAUR6CgoKDwhqM4AgUFBYU3HKOXbYA+2NjYyGrVqj1X3ri4OCwsLArXoEJAsSt/KHblD8Wu/FFc7YKC2XbmzJlHUspyeV4opSz2n0aNGsnnxc/P77nzFiWKXflDsSt/KHblj+Jql5QFsw0IlHr0scrQkIKCgsIbjuIIFBQUFN5wFEegoKCg8IbzSkwWZ0dKSgrh4eEkJibmep2VlRWXL19+QVbpj2JX/lDs0g9TU1Ps7OxethkKrxivrCMIDw/H0tKSatWqIYTI8TqVSoWlpeULtEw/FLvyh2JX3kgpiYqKIjw8/GWbovCK8co6gsTExDydgILCm4QQAmtrax4+fPiyTVEoIG02tCEqMep/CWu0P6xNrfH/wL/Q63ul5wgUJ6CgkBnlO/F6kMkJ6JFeUF5pR6CgoKCgUHCKzBEIIUyFEKeEEP8JIS4KIaanpfsIIW4JIYLSPq5FZcPLJDU1FRsbGyZPnpwpPSUlhUmTJuHq6oqjoyONGzdmz549AFSrVg0nJydcXV1xdXVl7Nix+aqzTZs2ODg44OLiQvPmzbl69Wq+8m/bto0rV67ojqdOncrBgwdzzTNw4EA2bdqkqz8wMFDv+mJjYxk+fDg1a9akQYMGtGrVipMnT+bL5sLGx8eHu3fv5itPaGgojo6OOZ5fsGABpqamxMTE5NuenO5pYGBgvv8/FBRyoijnCJKAtlLKWCGEMRAghNiTdm6ClHJTEdadhW3nIpi77yp3oxOoVNqMCV4OdG9Yucjq279/Pw4ODmzcuJHvvvtO98o+ZcoU7t27x8mTJ7GxsSEyMpLDhw/r8vn5+WFjY/Pc9f7xxx+4u7uzfPlyJkyYwPbt2/XKl5qayrZt22jXrh0eHh4AzJgx47nt0IchQ4ZQvXp1QkJCMDAw4ObNmy99BY6Pjw+Ojo5UqlSp0Mr09fXFw8ODrVu3MnDgwCznU1NTMTLK31fR3d0dd3f3QrJQoTgRk5T/B4aCUmRvBGk7nGPTDo3TPrKo6suNbecimLzlAhHRCUggIjqByVsusO1cRIHKnTlzJnXr1qV9+/b06dOHefPm6c75+voybtw47O3tOXHiBADx8fF4e3uzaNEiTExMALC1taV3794FsiM7WrVqxfXr1wFth+7h4YGjoyPDhg1Du/Nc+7T51Vdf0bp1a3788Ue2b9/OlClTcHV15caNG5me9nMqIztWrlzJ+PHjdcfe3t589tlnma65ceMGJ0+eZNasWRgYaP8Na9SoQadOnQCYP38+jo6OODo6snDhQkD75F2vXj2GDh1KgwYN6NChAwkJCbq2TJw4kcaNG1OnTh2OHj0KgFqtZsKECXh4eODs7MyyZct0NsyZMwcnJydcXFyYNGkSmzZtIjAwkL59++Lq6kpCQgJnzpyhdevWNGrUCC8vL+7duwfAmTNncHFxoV27dixZsiTHe3Hjxg1iY2OZNWsWvr6+unQfHx/ef/99unTpQocOHbK1J52//vorS7v8/f3p3LlzjvUqvJo8SXzC0P1DX3i9RbpqSAhhCJwBagFLpJQnhRAjgdlCiKnAIWCSlDKpIPVM33GRS3efZntOrVZzPkJFslqTKT0hRc2Xm87jeyos23z1K5Xi2y4NcqwzMDCQzZs3c+7cOVJTU3Fzc6NRo0bashMSOHToEMuWLSM6OhpfX1+aNm3K9evXsbe3p1SpUqhUqmzL9fT0xNDQEIABAwZk6lDzw44dO3BycgJgzJgxTJ06FYB+/fqxc+dOunTpAkB0dLTujSQkJIR27drRr1+/LOXlVsazfPjhhzg7OzNnzhyMjY1ZvXp1pg4Y4OLFi7i6uurampEzZ86wevVqTp48iZSSJk2a4O7ujp2dHSEhIfj6+uLt7U3v3r3ZvHkzH330EaB9sj516hS7d+9m+vTpHDx4kJUrV2JlZcXp06dJSkqiefPmdOjQgStXrrBt2zZOnjyJubk5jx8/pmzZsixevJh58+bh7u5OSkoKn3zyCX///TflypVjw4YNfP3116xatYpBgwaxaNEi3Nzccn1z8vX1pU+fPrRs2ZKrV6/y4MEDypcvD8Dx48c5f/48ZcuWZc+ePVnsSSe7dim8fjxKeMTQ/UO5o7pDqRKleJqctU+zNrUukrqL1BFIKdWAqxCiNLBVCOEITAbuAyWA5cBEIMs3SQgxDBgG2qdmf3//TOetrKx0nWlKcgpqtTonG7I4gXSS1Zoc86Ukp+TYWQMcPHiQjh07kpqaCoCXlxdJSUmoVCq2bt1KixYtUKvVdOjQgRkzZjBjxgzi4uLQaDSoVCrUanWW8qWU7NixA2vr//2xc7PhWdRqNX369MHMzAx7e3vmzp2LSqVi9+7dLFy4kISEBJ48eUKtWrVo06YNarWaLl26/O8+pqTo7Es/TkhIyLWMjNeo1WpdG1u2bMlff/2Fg4ODbqlvxrYkJCSQmpqabfsOHjzIu+++i0aj/bt16tSJY8eO0alTJ6pWrUrNmjVRqVQ4Ojpy9epVXd0dO3ZEpVLh4ODAzZs3dXYHBwezceNGAJ4+fcp///3HoUOH6NOnj+7vYGxsnKkNKpWKS5cuERwcTLt27XT319bWlvDwcJ48eYKbmxtqtZqePXuya9eubNvy559/8scffxAXF0fnzp1Zt24dw4YNIzExkTZt2ujq3b17d472ZNeu+Pj4HO9fYmIisbGxWb4zxQHFruyJTo1mUeQiotXRDCs3DAczh0y2lSxZUndcFHa+kH0EUspoIYQ/0FFKmT5+kiSEWA18kUOe5WgdBe7u7rJNmzaZzl++fFm3kWdWr5znm1UqFR2XnCYiOiHLucqlzdg0qkV+mwOAiYkJJiYmOhtKlCihO962bRvHjh3TPZE/fvyYwMBAmjVrptvsY2homGUjkhCCkiVL5rhBSa1W6946unbtmuVJ1NDQEF9f30xjx4mJiXz++ecEBgZSpUoVpk2bhpQSS0tLDA0NKVeunK4+Y2NjDAwMMh2bmZlhbGycYxnp16SXZ2FhgaWlJSNHjuS7776jbt26DBkyJEub3N3dmTx5MhYWFrqhoZzurYmJie7epNcFYG5uTmxsrK7uMmXKYGlpSVJSEhqNRpe+ZMkSvLy8MtVx5MiRTGVlvIfpbTA3N6dBgwYcP3480zXR0dG6+6RSqXRteLas8+fPc+PGDXr06AFAcnIyNWrU4PPPP8fU1JTSpUtnudfZ2ZNdu8zNzTEyMsr2f8XU1JSSJUvy7HemOODv76/Y9Qx3Y+/y8b6PiRNxLPdaTiPbRi/ctqJcNVQu7U0AIYQZ8DZwRQhRMS1NAN2B4KKyIZ0JXg6YGWcegjAzNmSCl0MOOfKmRYsW7NixQ/f0tWvXLkD7xBkQEEBYWBihoaGEhoayZMkSfH19MTc35+OPP2bs2LEkJycDcO/ePX7//Xe96jQ0NCQoKIigoCC9J3LTJThsbGyIjY3Vjflnh6WlJbGxsVnS81NGOk2aNOHOnTv8+eef9OnTJ8v5mjVr4u7uzrfffqubbwgJCeHvv/+mVatWbNu2jfj4eOLi4ti6dSvNmjXTq73P4uXlxa+//kpKSgoA165dIy4ujg4dOrBq1Sri4+MBdEMx6Z07gIODAw8fPtQ5gpSUFC5evEjp0qWxsrIiICAA0E7QZ4evry/Tpk3T/R/cvXuXiIgIbt++neXanOxReL258/QOA/cOJCYphuXtszqBF0VRvhFUBNakzRMYABullDuFEP8IIcoBAggCRhShDQC61UGFuWrIw8ODrl274uLiQtWqVXF3d8fKyootW7bQtm1b3WQwQLdu3fjyyy9JSkpi1qxZfPPNN3h4eGBubo6FhUWmTj3jHIGzszNr1659bhsBSpcuzdChQ3FycqJatWq6FUHZ8eGHH/Lxxx+zfPnyTJ19fsrISO/evQkKCqJMmTLZnl+xYgWff/45tWrVwtzcHGtra+bOnYubmxsDBw6kcePGgHZ1kYuLC1FR+d9MM2TIEEJDQ3Fzc0NKSbly5di2bRsdO3YkKCgId3d3SpQowbvvvst3333HwIEDGTFiBGZmZhw/fpxNmzYxduxYYmJiSE1N5dNPP6VBgwasXr2awYMHY2Jiwrvvvptt3evXr9ctDU6nR48erF+/Hltb20zpOdmj8PpyK+YWQ/YNIVmTzAqvFdS3rv/yjNEnaMHL/mQXmObSpUt6BWZ4+vSpXtc9DyqVSkopZVxcnGzUqJE8c+aM3nmL0q6CUJh2derUSR48eLBQynoT7ldhcenSpWIbaEWxS8u1x9dk6/WtZav1reTVx1dzvVYJTFPMGTZsGK6urri5udGrVy/c3NxetknFgujoaOrUqYOZmZluolVBQUHL5ajLDN43GENhyOqOq6lTps7LNunVFZ0rDvz5558v24RiSenSpbl27drLNkNBodhx4eEFhh8cTknjkqzssJIqpaq8bJMARWtIQUFB4YVwNvIsQw8MxaqEFT4dfYqNEwDFESgoKCgUOSfvnWTEwRGUMyuHT0cfKpXMW8JEqtVErVxFuc+/IGrVKmQOe54KA8URKCgoKBQhAREBjD40msolK7O642psLWzzzJMcGsqtXu/xcPFiDOLieLhoMbfee5/k0NAisVFxBAoKCgpFhF+YH2P/GUt1q+qs8lqFjZl+gpKh/9eXpGvXkGlaWjIhgaSrVwn9v75FYqfiCIqIlylDnZ4/r41fioS0IiGtUHTsC93HZ/6f4VDGgRUdVlDGNPv9NNlRompV0DwjjaPRYFK7diFbqeXNWDX0Wwu4fyFregUnGBFQJFW+bBnqokaRkNaSrjWVHxQJ6defHTd28M2xb3Ap58LSdkspWaJk3pnSSIl8QHI2caeFuTmle/UsTDN1vBlvBHaNwbBE5jTDEtr0AlCcZajTefaJd968eUybNi3H69euXVssJaSbNm1a6BLSTZs2LbCEdPfu3bO1Jx1FQvrNY0vIFr4O+Bp3W3d+e/u3/DmBu3e53a8ftyzdeWLrlOncE6vahFC3sM0FXpc3gj2Tsn/iB8zUqSA1oHnmyU2Tqs2zulP2ZVZwgnd+yLHK4ipD3bdvX8zMzAA4dOhQvvIC9OrVi+bNmxeqhHTr1q0pU6ZMgSSkHz16RMeOHQtVQrp169ZMmDAhx3uhj4S0sbGxIiGtoGP9lfXMPjmb5pWas9BzIaZGpnrnTQ4LI2zgINQqFbXfvodfdB/aN7uHvcl5wpMc+Tf6A7yuzgTWFbrdr4cjyAujEmBRHmIj0cbGEdrjZ98S8kFAQADdunXTdboZtfl37tyJp6cn5ubm9OrVi5kzZ7JgwQK9yi3soaH8yFgDWFhY0LZtW3bu3Em9evVISUnRqajqQ0BAAD169MDCwgKAnj17cvToUbp27Ur16tVxddUqxTZq1IjQDCsgevbsmSV9//79nD9/nk2bNunksUNCQjh48CCDBg3C3NwcgLJly2ax4+rVqwQHB9O+fXtA+3ZRsWJFYmJiiI6OpnXr1oA2tsKzekDprF+/nq1bt2JgYEDPnj3566+/GD16NADt27enbNmyqFSqXO3Jrl0KrydrLq5hXuA82ti14ac2P1EiH/1L0s2bhA0chExKwt5nNWa3VpJwdD07n3yLg+k/hCZ54GX9M3b1ika94PVwBLk8uSeoVFqpXtV9+NkFUhPByASGHwHLvJdx5YTMJUKXr68vx44do1q1agBERUXh5+dHs2bNCAsLy3fnnE5eMtTZYWRkpNP1h/8piebGkCFDdBLSgwYNynK+QYMG/Pfff2g0miwS0rndl4xCfIaGhrqhoYznDA0NdePuUkoWLVqEl5cXqvS/I7B3717dnEtOSClzlJDOKy9oJaRDQkJ0jiRdQjrdEaQ7uvS6ciozu3YpvH54n/fml3O/0L5qe35s+SPGhsZ65028eo2wwYMBsF+7FlOHOtxnDMd3BWNEAlcS38bdYiN2ppeh9YYisf/NmCMAsKwArn1BGGh/FsAJwKsjQ21ra8uDBw+IiooiKSmJnTt35pmnsCWkW7ZsqZetz6JISCsUd6SULD63mF/O/UKnGp2Y02pOvpxAwsWLhPXvjzA0pOo6rRO4GxLN9t9uYSiSMRAa3C02EpzQkXC7sQXut3LizXEEAK2/BPu3oPXEAheVUYa6Z8+eecpQb9++XSdDXa5cOV383+7du1OuXDndtZ6enrrln/379y+wncbGxkydOpUmTZrQuXNn6tbVb7Kpd+/eNG/ePFcJ6fv371OrVi2cnJwYOnQolSpVyiQh3aRJE4YMGULDhg2fy/YhQ4ZQv3593NzcaNKkCcOHDyc1NZWOHTvStWtX3N3dcXV11U3Sp0tIu7q6olar2bRpExMnTsTFxQVXV1f+/fdfAFavXs3o0aNp2rSpbmjvWdavX68LKJNOuoT0s+Rkj8LrjZSSBWcWsOz8MnrU6sHs5rMxMtB/kCUhKIiwgYMQFuZU/X0dJjVqcOfKY3b8fAYT9SPUwpSO1gtoYumLV9mf2XehBeFXnxRdY4r7R5GhfnGk21WYEtKFQXG/X8UJRYY6/+TXLo1GI7878Z109HGUM4/PlGqNOl/5406dklcausmQ9h1kckSElFLKW+cfyl9HHZS+o1fKE9NnyzsXIqTcMV5qvrWScsd4eefKY3lmb2i+6kFPGerXY47gJTFs2DAuXbpEYmIiAwYMeC1kqKOjo2nUqBEuLi6KhLSCQjZopIaZJ2ay6dom+tXvxwT3CXrNO6UTe+wY4aPHYFypEvarV2NsW56b5x6yz/s81gY36Ormh+lHK6GEOVT9kpiQE5RuPRE7yzLYOei/KS0/KI6gALyOMtSKhLSCQs6oNWqm/juV7Te2M9RpKJ80/CRfTkDl50fEuE8pUb069qtWYmRtzbXT9zm46iK2Rlfp3PgkJn18tAtaACwrENTwO9oU0dxAOoojUFBQUNCDFE0KXx39ir2hexntOpoRLvmLsvt0/34iPv8C0zp1sF+5AsPSpbn8713+WXuZSsYX6dT0P0p8sALyMdlcWCiOQEFBQSEPktXJTDg8gX/u/MP4RuMZ7Dg4X/ljduzk7qRJmDk5UcV7OYaWlgT7h3N4/TWqlAjinZY3MO65DAxfTpeca61CCDvgQ6AlUAlIAIKBXcAeKaUml7ymwBHAJK2eTVLKb4UQ1YH1QFngLNBPSplcCG1RUFBQKHSS1EmM9xvP0YijTGo8ib718qcAGr15M/e+mYK5uzt2v/6KYUkLgg6GcWzTdaqZnMKrzX2Mui0Gg5e3iDPHmoUQq4FVQDLwI9AHGAUcBDoCAUKIVrmUnQS0lVK6AK5ARyHEW2llLZBS1gaeAB8XRkMUFBQUCpv4lHjGHBpDQEQAU5tOzbcTePznn9z7+hssmjWjyvJlGJa0IHDXLY5tuk5Nk3/p2CEao+4LXqoTgNz3EfwkpewgpfxFSvmvlPK6lDJYSrlFSvkJ0AbIUfM3bfVSbNqhcdpHAm2BdH3kNUD3ArfiJWFoaIirqysNGjTAxcWF+fPn63bx6is13KxZM0AryqbP5HNOssmhoaGYmZnh6upK/fr1GTFiRKYdxUVNtWrVePTo0QurT0GhqIlLiWPkwZGcun+Kmc1n8n6d9/OVP2q1D5EzZlLS0xO7X5ciTE05+fd1Tu64RR1Tfzq8m4Lhu99DPiabi4och4aklMEAQoh+wDYppU4XQQjRWUq5E7ieW+FCCEPgDFALWALcAKKllOl77cOByjnkHQYMA+3uWH9//0znrays9JJqUKvVtF7fmsdJWXd7ljUpy85Oee+0zQkzMzOdouTDhw/5+OOPefDgAV9//TUODg7Mnj07RxvVajUqlYp9+/ahUqm4dOkSa9euzaRZlB2xsbE63Z1n06tXr87Ro0dJTU2lc+fO+Pr60rVr10x1ZicUl51d+UVKSWxsbKaNdHnVk5cthWFXUVMc7Urf7f7sd6Y48KrYFa+J59fIXwlLDqO/TX+swq3wD/fPMf+zWOzeQ8nt20l0cyOyV0+uHztG5DkNUdcE9c0OUNU5miMlekMGCXp9bSsS8tpoAEQDQUC9DGln9dmkkOH60oAf2rmG6xnSqwAX8spf0A1ljj6OOX4KgoWFRabjGzduyLJly0qNRiP9/Pxkp06dpJRSPnjwQL799tuyYcOGctiwYdLe3l7eunUrUxlNmjSRpUqVki4uLnL+/Pny1q1bskWLFrJhw4ayYcOG8tixY1JKKW/duiUbNGiQxZZn0ydOnCh//PFH6efnJ9u0aSP79Okj69WrJ6WUct26ddLDw0O6uLjIYcOGydTUVCmllHv27JEuLi7S2dlZtm3bVkopZVRUlOzWrZt0cnKSTZo0kf/995+UUspHjx7J9u3bS1dXV12bHj58mGv5FhYWcsqUKbJx48by6NGj+brXxXHjlpTF0y5lQ1n+yWjXk4Qn8v3t70vXta7yYGj+NlVqNBoZuWCBvORQV4ZPmCA1KSlSo9ZI/98vysXDD8nD4z+TmiMLntu2/EIhbii7hXYcf5MQYpqU8i8gX+8yUspoIYQ/8BZQWghhJLVvBXbkMrykLz+e+pErj69ke06dR8DnQXuziqoB1C1bl4mN8ydFUaNGDTQaDQ8ePMiUPn36dNq2bcvkyZPZu3cvy5cvz5L3hx9+YN68eTotoPj4eA4cOICpqSkhISH06dNH7whi8fHxHDp0SKdHdOrUKYKDg6levTqXL19mw4YNHDt2DGNjY0aNGsUff/zBO++8w9ChQ9m9ezdOTk46vZxvv/2Whg0bsm3bNv755x/69+9PUFAQ06dPp0WLFkydOpVdu3bp2pRT+f379ycuLg5HR0e9dZIUFF40UQlRDD0wlNsxt/nZ82da2eU2DZoZKSUPfpzDYx8fSr//HhWmT0ci8F8bzOUTD2losZWm7zsi3hpehC14PvRxBFJKeVYI0RrwFUI0AfJ8pxdClANS0pyAGfA22oliP+A9tCuHBgB/P7f1xRCZjfpmQEAAW7duBbS6NDnp92QkJSWFMWPGEBQUhKGhoV6bvG7cuIGrqytCCLp168Y777yDv78/jRs3pnr16oA2RsGZM2fw8PAAtLETypcvz4kTJ2jVqpVOMTVdSjkgIIDNmzcD0LZtW6KiooiJieHIkSNs2bIFgE6dOunalFP5oJ1T6dWrV57tUFB4EbTZ0IaoxKj/Jaz536/L2i+jWaVmepclNRoiZ83iyZ++lOnbF9uvv0JKOLjqPCFnovCw2IBHn6YI94GF14BCRB9HcA9ASvlICOGFtjPPOcjr/6gIrEmbJzAANkopdwohLgHrhRCzgHPAyucz/X/k9uSuUqlotiXnP+jqjqsLWr2OmzdvYmhoSPny5TOFbczOOeTFggULsLW11ck9m5rmHeCiZs2aBAUFZUl/VjJ5wIABfP/995mu2b59e7Y7JLOzPf26nK7PrnwAU1PTfM0LKCgUJZmcwDPkywmo1dybOpWYzVso+/Fgyn/xBRq1ZP/yIG6ej+Yty99p9FEHcPmwMMwuEvJcsySl7JThd42UcoKUUp9856WUDaWUzlJKRynljLT0m1LKxlLKWlLK96WUSQVrQvHg4cOHjBgxgjFjxmTpIFu0aMHGjRsBbbCVJ0+yKghmlFAGiImJoWLFihgYGLBu3bo8h7j0pV27dmzatEk3fPX48WNu375N06ZNOXz4sC54SvrQUKtWrXRSzf7+/tjY2FCqVKlM6Xv27NG1KafyFRReR2RqKncnTiJm8xZsRo2i/BdfoE7VsGfpWW6ej6ZFqdU0GtilWDsByOWNQAixA+1yz2yRUnbN6Vxxw9rUOlvvb21qXaByExIScHV1JSUlBSMjI/r165clvi9ox9n79OnDhg0baN26NRUrVqRkycxxTJ2dnTEyMsLFxYWBAwcyatQoevXqxV9//YWnp2emp/qCUL9+fWbNmkWHDh3QaDQYGxuzZMkS3nrrLZYvX64LHVm+fHkOHDjAtGnTGDRoEM7Ozpibm7NmzZpMbXJzc6N169bY29vnWn7VqlULxX4FheKCTE4m4osJqPbvp9xnn2EzbCgpSWp2Lz5DeEgsrUt74ziwH9R992Wbmicip2GLtDkB0E4MewNDMp6XUua97qmQcHd3l89OlF6+fJl69erlmTdjZKuXRVJSEoaGhhgZGXH8+HFGjhzJ0aNHX7pd2VEc7ld2KHbpz+XLl4mMjKRNmzYv25Qs+Pv7Fxu7nNbkHIL1woDsY6Cno0lKImLsOGIPH8Z28iTKDhhAcmIqO38O5P6tWNqWWUbdj4dDrbcLbGdB7pkQ4oyU0j2v63LbR6Dr6IUQsS+y43/dCAsLo3fv3mg0GkqUKIG3t/fLNklB4Y0m+FHwc+fVJCQQPnoMcf/+S4Vp31Lmww9Jik9hx4LTPLgTR3vrpdQe8hlUf77IfOlsOxfB3H1XiYhOoPKJf5jg5UD3htluuyow+ioc5X+2U0FH7dq1OXfuXKa04rYJSUHhTeFs5FlGHRqFAQZoyLr7PrchY3VsHOEjRhB/9iwVv/+e0j26kxCbzPb5p3h8N56O5ZZQY+hXYN+kQDZuOxfB5C0XSEjRzg1GRCcweYv2LaUonEFucwRlMxwaCiHKkGH/gJRSCcyqoKDwSnH87nHG+Y3D1twW727eVLCoAOg3/KJ++pQ7Q4eREBxM5XlzKfXuu8Q/Tebvn04S8yCOd20XU3X4TKj0fKFZMzJ331WdE0gnIUXN3H1XX6wjQCsNkZGzGX6XQI009XdTAAAgAElEQVRCt0ZBQUGhiDgSfoTxfuOxL2WPdwdvbMxs9M6b+uQJdz4eQmJICJUXLqBU+/bEPknk73kniX0cT+cKi7Eb/iNU0Gdlfd5ERCdkm343h/SCkpsjqCOlTCmSWhUUFBReIPtD9zPxyETqlK3DsreXUdq0tN55Ux89ImzQYJJv36bK4kWUbN2ap48S+HveSRJi4uhaeQkVhy+EcnUKxdYVR2/meK5SabNCqeNZctsPcFwIsU0IMUIIUa1IaldQUFAoYnbc2MGEIxNwtHFkRYcV+XICKZGR3O7Xn+TwcKos+42SrVsTHRnP1jknSHr6lG5VllBx1OJCcQJSShYdCmHWrss425XC1Dhz92xmbMgEL4cC15MdOTqCtCVH49IOFwohTgshFgghOggh9JOYfINJTU3FxsaGyZMnZ0pPSUlh0qRJuLq64ujoSOPGjdmzZw+glXJ2cnLC1dUVV1dXvWSsM5KcnMynn35KzZo1qV27Nt26dSM8PFx3XghBv379MtlYrlw5OnfuDICPjw/Vq1fXSVlnt7rJ399fd31u5FdeW0GhKPjr2l98HfA1HrYeLGu/DMsS+i/1TYmI4PZH/Uh98AD7Fd5YNG3K47txbJ1zgtTYGLpX/Q3bUb9B2YKPkksp+XHvVX46cI2eDSuzZWRzfujpTOW0N4DKpc34vqdTka0aynWHsJTytpTyNylld6AZsAOtZtBRIcSuIrGoiJBqNVErV3HtraZErVqFLKSdujmxf/9+HBwc2LhxYyaZhilTpnDv3j1OnjxJcHAwO3bsyLSCyM/Pj6CgIIKCgvjll1/yVedXX32FSqXi2rVrhISE0L17d3r27Kmr38LCguDgYBIStOOMBw4coHLlzP9YPXv2JCgoCH9/f7766isiIyOfq/3//vsvoDgChZfHukvrmHF8Bi0qt2Bxu8WYG5vrnTf59m1C+/VDHROD/epVmDdqxKNwFVvnnoSEJ/So4Y3NyJVQ2r7Admo0kmnbL/Lb4Rv0bWLPvPddMDI0oHvDyhyb1BafjhYcm9S2yJwA6CExkU7afME54HcpZWPSYgW8CiSHhnKr13s8XLwYdXQ0Dxct5tZ775OcJqfwvMycOZO6devSvn17+vTpw7x583TnfH19GTduHPb29pw4cQLQqoJ6e3uzaNEinW6/ra0tvXv3LpAd6WWvXr2aBQsW6PR8Bg0ahImJCf/884/uunfeeYddu3bpbOzTp0+25ZUvX56aNWvmKg8xbdo0Bg8eTJs2bahRo0Ymx5W+c3rSpEkcPXoUV1dXFixYUOB2Kijow/Lzy5lzeg5v27/Nz54/Y2qUt1ZXOkk3bnD7o37IhESqrvHBzNmZyNCnbJt7CqPkR/So40PZET5QqmKB7VRrJBM3n2fN8dsMbVmdWd0dMTB48YFq8txHkCYf3TXt2iDgoRDisJQyq5bCS+L+d9+RdDl7GepUtZrk8+chNVWXJhMSSLp8mRudu2DeMPulXib16lLhq69yrDMwMJDNmzdz7tw5UlNTcXNzo1GjRoBWeuLQoUMsW7aM6OhofH19adq0KdevX8fe3p5SpUrluI/A09NT15EPGDCA8ePH63UPMpadEXd3dy5evEi7du0A+PDDD5kxYwadO3fm/PnzDB48WBdcJyM3b97k5s2b1KpVK9d6r1y5gp+fHyqVCgcHB0aOHImxsbHu/LPy2goKRYmUkkXnFuF9wZtONToxq/ksjAz0DwifePUqYYMGg6EBVdeuwaR2be5dj2bHz2cw0zygm4MvpYasAfOyeReWBylqDeM3BLHz/D3GtavNp2/XzlbI8UWgzx2yklI+FUIMAVZLbQD680VtWGFiYGaGJpuO18Ds+WfgAwIC6NatG2ZpZWSMLLZz5048PT0xNzenV69ezJw5U++nYT8/P2xs9F/Wlo6UMkc10Izpzs7OhIaG4uvry7vvZtVA2bJlC6dOncLExIRly5bp5KhzolOnTpiYmGBiYkL58uWJjIzEzs4u3/YrKBQUKSVzTs/h98u/06t2L6a8NQVDg9zVbqVazWOfNZRbupT7vXoRs20bBubm2K9ehUn16oRffcKuRWcpKe/SrcEWSg5eB6ZWBbY1MUXNmD/PcvDyAya/U5fhrWsWuMyCoI8jMBJCVAR6A18XsT3PRW5P7iqVCo2fH/emT0fGxevShbk5FaZ8g1XX59POy01a2tfXl2PHjum0/aOiovDz86NZs2aEhYU9965itVqte+vo2rVrpgAvtWrV4vbt21m0b86ePZsl/GXXrl354osv8Pf3Jyoqsxhfz549sw2ckxMZQ1MaGhqSmuHNS0HhRaGRGmaemMmma5voW68vEz0m5vl0nRwaSvin40m+fRuDhASerF0LRkZU/uVnTKpX5/bFKPYsDcJKhNHNZQ/mA9aBSclcy9SH+ORUhq09Q8D1R8zs1oB+TasVuMyCos8cwQxgH9oQk6eFEDWAkKI1q3Ap6emJeObJQBgaUtLT87nLbNGiBTt27NDFh00fd3/69CkBAQGEhYURGhpKaGgoS5YswdfXF3Nzcz7++GPGjh1LcnIyAPfu3eP333/Xq05DQ0PdRPKzUb4sLCwYMGAAn332mU6yeu3atcTHx9O2bdtM1w4ePJipU6fi5JSz6FZh8ay8toJCYZOqSeXrgK/ZdG0TQ5yG6OUEAEL/ry9J164hEzJs0tJouPvZ59wMesjuJUGUETfp7rYf84E+heIEniam0H/lKf698Yh577sUCycAuUtM9AH2p4Wm/Cs9XUp5E3ilwkwZWlricPpUoZbp4eFB165dcXFxoWrVqri7u2NlZcWWLVto27Ztpiflbt268eWXX5KUlMSsWbP45ptv8PDwwNzcHAsLi0ydesY5AmdnZ9auXau3Td9//z1ffPEFderUwcDAgLp167J169YsXwo7OzvGjRuXQymFy7Py2vrOeSgo6EOKOoWJRydy4PYBPmn4CcOc9V/DYlKrFvGnnukXNBqi6rYnaNl5yhldo0vjE5j08QGjgq+YfxKXTP9Vp7h87ymL+rjRybngk82FRW4y1JOADoAxcAjYA5ySzxNuq4AUVxnq2NhYSpYsSXx8PK1atWL58uW4ubnplbc4yheDYld+KY52vSky1EnqJD7z/4wj4UeY4D6B/g365yv//VmzefLM2/h9+5ZcrtGbisaX6dQ0iBIf/AaGxjmUoD8PVIn0W3GKW1Fx/PaRG23r2uqd92XLUP8A/CCEsES7d2Aw8JsQ4jKwF9gnpXy+ReavCcOGDePSpUskJiYyYMAAvZ2AgoJCwYhPiWes31hO3TvFlLem0Nshf0uwo7duI+jwfUqVqU2ZJ9qR7oiKzbla40NKGdylc6vLGPdcBob6rzjKiYjoBPp6n+CBKgmfgR40q5X/xSBFTZ6tlFKqgK1pH4QQ9YF3gLWAV5FaV8xRNkopKLx4VMkqRh8azX8P/2N2i9l0qdkl70wZePzHH0TOnEV5R3fO1RqMV+m5RKVW46rq/zAghVZV9mDcax0Y6L3NKkdCH8XRd8VJniaksO7jxjSqWvBlp0VBni0VQjQXQlik/f4R2khlm6SUuToBIUQVIYSfEOKyEOKiEGJcWvo0IUSEECIo7VP847gpKCgUC2KSYhi6fygXHl5gTqs5+XYCj5Z7EzlzFiXbtcN12Wy8rH9m15OvCVB9jAEpdCrzPVXH/FQoTiAkUkXvZceJT07Fd9hbxdYJgH7LR38FXIQQLsCXwEq0bwOtc80FqcDnUsqzacNLZ4QQB9LOLZBSzsslr4KCgkImHiU8YtiBYdyOuc1Cz4W0rpJXF/Q/pJQ8XPgzUcuWUapzZyp9/x0YGXHH6kNSH2l3HTe0+Bv75m5gWaHAtgZHxNBv5UmMDA3YMLwpdWyL1zzSs+jjCFKllFII0Q34WUq5UggxIK9MUsp7wL2031VpcwtFJ5ahoKDw2nI/7j5D9w8lMj6Sxe0W07RSU73zSo2GyO9/4Mm6dZR+/30qTPsWhAFHfS9z4UZtDEmmocU2ghM6YlfZmYJuhzxz+wkDV5+ilKkxfwxpQjUbiwKWWPTkuGpId4EQh9FODg8CWgEPgSAppd6L0NNkrI8AjsBnwEDgKRCI9q3hSTZ5hpGmZ2Rra9to/fr1mc5bWVnlKX8A2k1Y6csxixOKXflDsUt/rl+/TkREhE7vqTiRvtIuPzxKecTiB4uJU8cxovwIaprmYxeuRkOp3//A7N9/iXu7HbG9eiEl3D2pJvq2AYYk0abyWhzUe/jPchAnwrpQpZnAwvb5pB4uRan5+WwipU0EX3qYYm1W8CGm57ln6Xh6euq1aggpZa4foALazrtl2rE90D+vfBnyl0Qb7axn2rEtYIh2fmI2sCqvMho1aiSf5dKlS1nSsuPp06d6XfeiUezKH4pd+nPp0iXp5+f3ss3IlvzadSv6lmy3sZ1s9mczeeHhhXzl1SQny/Dx4+Ulh7rywS+LpEajkakparln6Vm5ePgh+fcnP8k729dL+fSefDK/qZRP78s7Vx7LM3tD81VPOocu35e1v94t28/3l5ExCc9VRnYU5G8JBEo9+ml93NV4KeV8KeXRNMcRBjTQxxsJIYyBzcAfUsotafkjpZRqKaUG8AYa61NWQTi77zbhVzO/dIRffcLZfTkra+rD/fv3+fDDD6lZsyb169fn3Xff5dq1awUqMzumTZumUzadOnUqBw8efK5ygoKC2L17d77yhIaGYmZmRsOGDalXrx6NGzdmzZo1z1V/YREdHc3SpUtzPP/LL79Qr149+vbtW6R2DBw4ECcnJ3777Tdd2saNG6lfvz4NGjTg//7v/wC4ceMGrq6uxfIJvTgT8iSEgXsHkqJJYZXXKhxt9A8DqUlKIvyTsTzdvYfyX35JuU/GkJqiYffiM9z47wktSq2m67Cq2HX5ACwrENTwO7C0xc6hDG5eVfNt6+4L9xi+7gx1bEuyflhTypfSX+20OKDPHEF7YOIzae9kk5YJod3OuhK4LKWcnyG9otTOHwD0AIL1N/f5KF+tFPu8g/Ea6oidQxnCrz7RHT8vUkp69OjBgAEDSB+2CgoKIjIykjp1nj9aUV5DDc9KS+SHoKAgAgMDsxWby42aNWty7tw5QKtK2rNnTzQaDYMGDXpuWwpCuiMYNWpUtueXLl3Knj17qF69eqb01NRUjIwKvi48IzNnztQF+wkJCeH777/n2LFjlClThgcPHgDa+xcUFKQ4gnxwMeoiww8Mx8TABG8vb2pY6R/8RRMXx53RY4g/eZIK06ZR5sMPSE5IZdeiQO7ejMWztDf1Bw0Ch46FYuvmM+FM2PQfDe3LsHqQB6VMC74B7UWTm8TESGAUUOMZtVFL4F89ym4O9AMuCCGC0tK+AvoIIVwBCYQCw5/D7kwc3XiNR3disz2X3rFaWJVgxy9BmFuVID4mmTIVzDm98xand97KNp9NlZK07J1zh+7n54exsTEjRozQpbm6ugJaJ/Hll1+yZ88ehBB88803fPDBB5nSpZRMnTqVDz74AH9/f6ZPn07FihUJCgri0qVLzJ49m7Vr11KlShXKlSunE5sbOHAgnTt35r333mPGjBns2LGDhIQEmjVrxrJlyxBC0KZNG5o0aYKfnx/R0dGsXLmSJk2aMHXqVBISEggICGDy5Ml07tyZTz75hAsXLpCamsq0adOy6BI9S40aNZg/fz6ff/45gwYNIi4uLksZ3bp14+LFiwwaNIjk5GQ0Gg2bN2+mdu3arF27lnnz5iGEwNnZmXXr1vHw4UNGjBhBWFgYAAsXLqR58+ZMmzaNsLAwbt68ye3btxk/fjxjx45l0qRJuqfs9u3bM3fuXJ19I0aM4ObNm3Tt2pXBgwcTExPD3bt3CQ0NxcbGhlWrVjFy5EgCAwMxMjJi/vz5eHp64uPjw7Zt21Cr1QQHB/P555+TnJzMunXrMDExYffu3XkqsXp7ezN69GjKlCkDaGM6KOSfoAdBjDw4klIlSrGiwwqqlKqid17106fcGTachAsXqPTjD1h17UpiXAo7FpziUXg8Hax/pfbH46CG/iuOcmPdidtM2RZM81rWePd3x7xE4T5ovChys/pPtLIS3wOTMqSrpJSP8ypYShkAZDfjkr+xiULCxNwYc6sSxD5OomRZE0zMC+a1g4ODdZ3zs2zZsoWgoCD+++8/Hj16hIeHB61ateLff//VpYeGhuLp6UmrVq0AOHXqFMHBwVSvXp0zZ86wfv36bGMdZGTMmDFMnToVgH79+rFz506d0mhqaiqnTp1i9+7dTJ8+nYMHDzJjxgwCAwNZvHgxoI1o1rZtW1atWkV0dDSNGzfmyJEjeUomuLm5ceWKNv7D7Nmzs5Tx9ttv89tvvzFu3Dj69u1LcnIyarWaixcvMnv2bI4dO4aNjQ2PH2v/jcaNG8f48eNp0aIFYWFheHl5cfnyZeB/8Q7u3btHo0aNGDlyJD/88APBwcEEBQVlse23335j7969OjnvadOmcebMGQICAjAzM+Onn34C4MKFC1y5coUOHTrohvOCg4M5d+4ciYmJ1KpVix9//JFz584xfvx41q5dy6effprrfUkvp3nz5qjVaqZNm0bHjoXz1PmmcPLeST755xPKm5dnRYcVVLDQfyln6uPHhH08hKTr16m8cAGl2rcnLiaJ7fNPEfMgno7ll1B96BSo4lEoti4/coPvdl+hbd3yLO3rhqlx8Vo0kB9ycwRSShkqhBj97AkhRFl9nMGLIrcn93QtmPThIPd3qxF8JAKPztWxcyhTJPYEBATQp08fDA0NsbW1pXXr1pw+fTpTevny5XXppUqVonHjxrqhjKNHj9KjRw/MzbWh9brmIJXt5+fHnDlziI+P5/HjxzRo0EDnCHr27AlAo0aNCM0hEtv+/fvZvn27bv4hMTGR8PBwKlTI/csnM6w0y66MsLAwmjZtyuzZswkPD6dnz57Url2bf/75h/fee08XbyH9CfvgwYNcunRJV+bTp091iqXp8Q6sra118Q7yS9euXXVxIwICAvjkk08AqFu3LlWrVtV14J6enlhaWmJpaYmVlZXuXjo5OXH+fN4hOFJTUwkJCcHf35/w8HBatmxJcHAwpUvrHyz9TeZo+FHG+4+nimUVvDt4Y2OmvxRDSmQkYYM/JiUigipLl1KyZQtUjxP5e95J4p7E0anCEqoM/wEqFFxxV0rJwoMh/HwohE5OFVnwgSsljAq+OuhlktcbQWe0K34kmZ/uJVDwiM0viIxzAnYOZajsUCbT8fPQoEEDNm3alO25jB2lPumglZHOSF4yuomJiYwaNYrAwECqVKnCtGnTSExM1J1PVz/NLUaAlJLNmzfj4OCgS9NHMvrcuXM6wb/sygCoV68eTZo0YdeuXXh5ebFixYocg+doNBqOHz+u66wzUhjxDjLe29z+BhnrMjAw0B0bGBjoVa+dnR1vvfUWxsbGVK9eHQcHB0JCQvDwKJwn0NeZg7cPMuHIBGqXrs2y9ssoY6r/9zI5PJywgYNQP3mCvfdyzD08iI6M5++fTpKsiqOb3a9UGL4QbGoX2E4pJd/vucLyIzd5r5EdP/ZyxvAlhJYsbHJ0Y1LKzmk/q0spa6T9TP+8Mk4A4EHo00ydvp1DGbyGOvIg9Olzl9m2bVuSkpLw9vbWpZ0+fZrDhw/TqlUrNmzYgFqt5uHDhxw5coTGjRtnSn/06JEu/VlatWrF1q1bSUhIQKVSsWPHjizXpHf6NjY2xMbG5uiUMvJsbAAvLy8WLVqk6xzTJ4RzIzQ0lC+++EL3VJ1TGTdv3qRGjRqMHTuWrl27cv78edq1a8fGjRt1wXDSh4Y6dOigG64Csh3yya0d+aFVq1b88ccfgHYoJywsLIsTe166d++On58fAI8ePeLatWvUqPFKfVVeCrtu7uKLw1/QwLoBK7xW5MsJJN28ye2+H6FRqbD3WY25hwdREbFsmXMCdWw03astpcKopYXiBDQayTfbgll+5Cb93qrKnNfECYCeweuFED2FEPOFED8JIboXtVGFjZtX1SxP/s+7TCwdIQRbt27lwIED1KxZkwYNGjBt2jQqVapEjx49cHZ2xsXFhbZt2zJnzhwqVKiQKb1z58669Cz2urnxwQcf4OrqSq9evWjZsmWWa0qXLs3QoUNxcnKie/fuej11enp6cunSJVxdXdmwYQNTpkwhJSUFZ2dnHB0dmTJlSrb5bty4oVs+2rt3bz755BPdiqGcytiwYQOOjo64urpy5coV+vfvT4MGDfj6669p3bo1Li4ufPaZNuz1L7/8QmBgIM7OztSvXz/TcszssLa2pnnz5jg6OjJhwoQ8252RUaNGoVarcXJy4oMPPsDHxyfTm0BB8PLywtramvr16+Pp6cncuXOxtrYulLJfVzZf28zko5Nxs3VjefvllCpRKu9MaSRevqwNMq9WY792LWZOTkSGPmXr3JMYJEbRo+Yqyo1aBWWe/3ueTqpawxd//ccfJ8MY3roGM7o1eClB5ouMvDYaAEuB/Wh3Fg9Cu8t4iT6bFArro2woe3EodunPgAED5Nq1a/W+3sLCogit+R+vyoay3y/9Lh19HOXwA8NlQkr+NmDFnzsnr3g0ltfaeMqkW7eklFJGXHssl405KNeO+lPG/NxVythHz2XXsySlqOXI3wNl1Yk75c8Hr0mNRpMvWwvKi9hQps9ap9aAY1qhCCHWABeKxCspKLxCWFlZMWvWLOLi4jItI36WGzdu0KtXL2xt9Q9G8rqz4sIKfj77M+3s2zGn1RxKGJbQO2/ciRPcGTUao3I2VF21CuPKlXXxhUuJcLo67aDkwLWFFmR+5O9n8Lv6kG861WNIy9dzqE8fR3AVraxE+jbcKkDeSyheADKHyUcFhRfBzz//rFeEsvQNZS8C+eIDCOZJmw1tiEqM+l9C2sb0EgYlmNt6LsYG+i/lVvn7EzF2HCWq2lNl5UqMy5fnxtkH7F9xgbIGN+nq9g9mH62FEuYFtjsuKZUhawI5cSuK2T0c6duk4ENMxZXcNpTtQLs6yAq4LIQ4lXbcBP02lBUppqamREVFYW1trTgDBQW0TiAqKgpT0+Ilb5DJCWQgWZOcLyfwdO9eIr6YgGndulTxXo5RmTJc+fcu/6y7jK3RFTo3OYPJhz6FEl84JiGFQatPEXQnmvm9XejRsKCapMWb3N4IinW8ADs7O8LDw3n48GGu1yUmJha7LwYoduUXxS79MDU1xc7Ojtu3C6ajVdyI3ryFe1OmYNawIVWW/YZhyZJc8A/nyPpr2JU4z7stQjDutbxQQktGxSbRf9UprkWqWPJ/brzjVHyCzBcVud21IzKP90whhMjrmqIifa12Xvj7+9OwYcMXYFH+UOzKH4pdby6Pf/+DyFmzsGjeHLtFv2Bgbs7ZvaEc33aTaian8Gr7CKMuiwslqljk00Q+WnGSsMfxLO/vjqfDmyETkpsj8BNCbAb+llrFUQCEECWAFsAAwA/wKVILFRQUXlniU+ILlP/RsuU8XLCAkm+3o/L8+QhjY05sDeHMvjvUNj1Cu3fAsMNceM7h4W3nIpi77yoR0QnY/nuQVLWGhBQNPoMa07Tmm7P0NzdH0BEYDPgKIaoD0YAp2lgC+9GGm3wxM2AKCgqvHDFJMYw+lEWhRi+klDxcsJCo5csp1aULlb6bDYZGBGy4wnn/e9Q320/rbjYYtPniue3bdi6CyVsukJCiBiDyaRIAn75d+41yApCLI5BSJqLdQ7A0La6ADZAgpYx+UcYpKCi8mqTHFw6NCcXS2BJVStad4Nam2Xe2UqMh8rvvefL775Tu3ZsK075FIvBbE8yVkw9xMf+b5u/XQzTNecmuPszdd1XnBDLyV2A4n779/FLyryJ6zaxIKVNIiz+soKCgkBvhqnCGHRjGo4RHLGm3JFN8YX9/f9q0aZNjXqlWc2/KVGK2bKHsoEGU/3ICGrXkwIogbgQ9waPkBjz+ryXC7aMC23k3OiFf6a8zr6Z4toKCQrHk+pPrDDswjCR1Eis6rMC5nLPeeWVyMhETJ6LasxebMWOwGT0KdYqGvb+e5fZlFc1KraHhgO7QoEeh2FrWogRRcclZ0iuVzip++LqjOAIFBYVC4fzD84w6NIoSBiXw6ehD7TL6C71pEhOJGPcpsYcPU/7LL7EePIjkxFR2Lwok4kYsbUqvpMHgQVCnQ6HYuvlMOI/jkhECMq57NDM2ZIJX4YgQvkrkud5KCDFGCFE0wv0KCgqvBSfunWDI/iFYGluy5p01+XMCcXHcGT6C2CNHqDB9OtaDB5EYl8Lf805y94aKt62X0WD4qEJzAiuO3uTzv/6jaU1rvu/hROW0N4DKpc34vqcT3RtWLpR6XiX0eSOoAJwWQpwFVgH7XtbeAQUFheLHoduHmHBkAlVLVWV5++WUMy+nd151TIw2tGRwMJXm/IhVly7EP01m+08nefIgno7lf6XGsG/Azr3AdkopmbPvKr/63+Adxwos/NAVEyNDPmxsn+fcxetOnm8EUspvgNpoA9EPBEKEEN8JIWrmlk8IUUUI4SeEuCyEuCiEGJeWXlYIcUAIEZL2U3nbUFB4RdkaspXPDn9GPet6+HT0yZcTSH38mNsDB5F46RKVFy7AqksXVI8T2frjcWIexNK5wmJqjJpVKE4gVa1h8pYL/Op/gz6N7Vn8f26YGL26oSULG7224qW9AdxP+6QCZYBNQog5uWRLBT6XUtYD3gJGCyHqo41/fEhKWRs4ROZ4yAoKCq8Iay+uZeq/U2lSoQne7b2xMtFf7TMlMpLbH/Uj+dYt7H79lVLt2xP9IJ6tPx4n/rGKrnZLqDJ6AVRwLLCdiSlqRv95lvWn7zDGsxbf9XB8bQLKFBZ5Dg0JIcai3UX8CFgBTJBSpgghDIAQ4Mvs8kkp75G25FRKqRJCXAYqA92ANmmXrQH8gYkFaoWCgsILQ0rJ4qDFLD+/nECHT4EAACAASURBVPZV2/NDyx/ylJGWajWPfdZQbulSHvTtS8yuXWiio7Ff4Y25uztREbFsn38KTYKK7tW8KTd8KZS2L7CtqsQUhq09w/GbUUztXJ/BLfKWpXkT0WeOwAboKaXMpGIlpdQIITrrU4kQohrQEDgJ2KY5CaSU94QQb4aYh4LCa4BGavj+5Pesv7qenrV7MvWtqRga5D7EkhwaSvin40m+fRuDhASili8HAwMqL1yIubs7D24/ZfuC0xgmP6FHrbWUHboCLLNG7ssvj2KTGLj6FFfuqVjwweuvIFoQRF7zvkKIt4CLUkpV2rElUF9KeVKvCoQoCRwGZksptwghoqWUpTOcfyKlzDJPIIQYBgwDsLW1bbR+/Xp925SJ2NhYSpYs+Vx5ixLFrvyh2JU/isIutVTz+6PfCYwPpF2pdnQr3U0vCXibLyZgEBeHyNDXSCHQWFhw+8s53DmcgjlRtKuymluNxpBqrH+4ypx4GK9hXmAiTxIloxua4FIu92fe4vp3hILZ5unpeUZKmfckS14hzIBzpDmMtGMD4Kw+4c8AY2Af8FmGtKtAxbTfKwJX8yonu1CV+vIqhOwrTih25Y83xa6ElAQ56uAo6ejjKL3Pe+crXGNov/7ykkPdLJ/TAyfJ30YdkL+PWiNVv74vZUJModh65d5T2Xj2Aen07V55+laUXnmK699RyhcTqlKfyeJMUtNSSg36zS0ItCuNLksp52c4tR3tnANpP//WwwYFBYWXhCpZxfADwzn6/+3dd3gUVdvA4d/ZTe89kAaE3kOvShB8QcWCXVF8FQQVFdTXimLHguUTUBFBARtWVIoiIFV6SQgdQhKSkN7b9vP9sYuEFjbZNMi5r2uv7M7MmXkyhHl2Zs6cJ20jL/V/ifFdx1erGJRrmzbnTMsJ78sOt6vwFSmM7vUXXg8sADfHzwR2peRz+2dbkBJ+eGgAvVsGOLzOpsCeRHBcCPG4EMLZ9poMHLej3SDgXuAqIUSc7XUt8DZwtRDiKHC17bOiKI1QXkUe41aOY2/OXt698l1ub397tdrnf/U1cRuzKaj0gFlmaB8S2tyDB/nc1H8zHmO/qJXSkmsPZzNm3jb8PZz5+eGBdGjmeGJpKuy5WfwQMBN4EWupyjXYrt1XRUq5CbjQ14Zh9gaoKErDyCjNYMKqCWSWZTLzqplcEXGF3W1lpWGkQ3v0Y3ebBxjhN4NCczgHiseiwcKQiN9xu2thrVQV+3VPOv/7MZ72zbxZ+EBfgrwcL1fZlFz0X0BKmQ3cWQ+xKIrSSCQVJTFh1QTKDGV8dvVn9AztaXdbaTKR8fLLFP38i3UY6ScnEPzePSzPnYoJNzQYudb/bVo89m2tJIEv/0ni1aUH6B8dwOdje+PtZn8dZMXKnmv9wcCDQMvKy0spH6i7sBRFaSgH8g7w0KqHEELwxcgv6BDQwe62looK0p98itK1awmaNImgRychJRz3HIMp11rbuYfnb7QY1AN8HKsFLKXkg1VHmPX3MUZ0DuWjO3vg5qyeFq4Je9Lxb8BGYDVwbhUHRVEuGzsyd/DY34/h4+LD3Kvn0tK3pd1tzYWFpD78CBVxcTR7eRr+d92FyWhm9bx4ElOi0aKnh+dv7KsYSUR4Nxzp1W+2SF76bR/fbjvBHb0jeXN0F5y0jtcsbqrsSQQeUkr15K+iXObWp67nqfVPEeYVxtyr59LM0/6HuowZGZx48EGMKScI//BDfEaOQF9hHUb65PFynISO6zouIzx/MeFdI1n5TSAjPPyJaF/9ocb0JjNPfB/HioRMHo5tzTMj2lerF5NyLntS6DJbbx9FUS5Ty44vY/LaybTxa8PCkQurlQT0iYkk33U3pswsIufNw2fkCMoK9Sx5ewuZx4tp67GZ6253I+L+aRT5diLilvGMeLAL2cnF1Y6zVG/igQU7WJGQyYvXdeTZkR1UEqgF9pwRTAZeEEIYAAPWnkBSSqn6ZinKZeC7Q98xfdt0+jTrw8yhM/Fysf8p1vI9e0h76GFwcabF11/h1qEDBZll/P7BdvQlFYxqPofIB1+DsBgA4npMJ9Y7lIj2VPtsIK9Uz/0LdrD/ZDHv39adW3qpISNqiz29hrzrIxBFUeqXlJK5e+cyO242sZGxvDfkPVy19ne7LFm7lvQnnsQpNISo+fNxiYgg83gRy2fuRBiKGd1yPsEPzoIAxwd6Sy+s4N7520gvqOCze3oxvFOow+tUTrP3CeExQCsp5etCiEisQ0Rsr/PoFEWpExZpYcaOGXx98Guuj76e1wa9hpPG/q6chT//Qsa0abh17EjkZ3NwCgwkOSGXlXPi8CCLGzr+jO8DC8AzyOFYj2aVcO/87ZQZTHw9vh991NPCtc6ef/lPAAtwFfA6UAp8DPSpw7gURakjJouJVza/wm+Jv3F3h7t5tu+zaIR9PW6klOR9Po+cDz7Ac+BAwmfOROvlyYF/TrLu64MEaY8zqudGPMYsAlfHB3HbfaKABxbswFmr4YeJA+jYXF2Rrgv2JIJ+UsqeQog9AFLKAiFE1YOPK4rSKOnNep5Z/wx/p/7NI90f4aHuD9l9s1VaLGS9/TYFi77C57rrCHtrOjg7s3P5cbYtTSbSJY6RVyThcvMXoHX8oa71R3J46KtdhPi48tUD/YgKdHwYCuX87EkERiGEFuvwEqceMLPUaVSKotS6MmMZk/+ezLbMbTzX9znGdBxjd1tpMHDy+RcoXr6cgPvGEvLss0gEm749QMLGLNq5reeq65zQDp8FtdCL5/f4kzz1QxxtQrxZ+EAfQrzdHF6ncmH2JIKZwBIgRAjxJnAr1nGHFEW5RBTqCnl49cMczD/I9MHTub719Xa3NZeWkf74Y5Rt3kLI/54iYNw4zCYLq+ftITG+iBiPXxl4exdE/4sOQWaXRVuSefn3/fRpGcC8+3rjo4aMqHP29Br6RgixC+tAcQK4SUp5sM4jUxSlVmSVZTFx1URSS1L5MPZDhkYNtbutKS+P1AkT0R06RPPp0/G7eTT6cqP1QbGkCgb5LCLmvpug800Oxyml5P9WH+WjNUcZ3jGU2XerISPqywUTgRDCR0pZLIQIALKB7yrNC5BS5tdHgIqi1NyJ4hM8+NeDFOoL+XT4p/Rt3tfutobUVE6MH48pK5uIj2fjHRtLWaGepR9uoyBLx9VB82g37jFoZf+opBditkhe+X0/X21N4bZeEbx1c1c1ZEQ9quqM4FtgFLAL2/0BG2H7HF2HcSmKUk2x38eSp8s7PWGh9YdA8N1139E5qLPd69IdPMiJByeA0UjUl1/g0aOH7UGxbdYHxcI+J3LimxBq/zovxGCy8OQPcSzbm8HEK6N57hr1tHB9u2AikFKOsv10/GkQRVHq3BlJoBKJrFYSKNu6jbRJk9D4+BC1cAGurVuTebyIZTN3ojEUMTp6IcEPzga/SIdjLtObeOjrXWw8msvz13Rg4pDWDq9Tqb6LnnsJIUYLIXwrffYTQjh+QVBRlEan+M+VpD74IE7Nm9Hyu29xbd2a5L25/Pb+DlyNmdzSeTHBk76slSSQX2bg7nnb2JyYx7u3dlNJoAHZcxHuZSll0akPUspC4OW6C0lRlIaQ/+23pD/xBG5du9Ly669xbtaMA/+cZMWn8fiL49zSdzW+ExaCh+NP9p4srOC2OZs5mFHMnHt6cXtvxxOLUnP2dB89X7JwvKyQoii1xmQx1bitlJLcWbPI/eRTvIYOJfyD9xFubuxcnsi2pSnWB8WGpuNyw+c1rij26550Zqw8THphBSH/rMZgMmO2wFcP9KVfdGCNY1dqhz1nBDuFEB8IIVoLIaKFEB9ivYFcJSHEF0KIbCHEvkrTXhFCpJ9VzF5RFAeUGkp57O/HatRWmkxkTnuZ3E8+xffWW4iYNRNc3dj47QG2LU2hndt6rrvRhMtN7zuUBJ7/JYH0wgoAskv0FFaYmBgbrZJAI2FPIngM6/DT3wM/Ajpgkh3tFgAjzzP9QylljO21wt5AFUU5V0ZpBmP/HMuWk1vwdPY87zKBbuc/2Fp0OtImT6Hwxx8JfGgizV9/HbMU/DVnNwkbs4jx/I3h97RCO+xZh54WnrHyMBXGc4sbfrcttcbrVGqXPQ+UlQHPVXfFUsoNQoiWNYhJURQ77M/dz6N/P4rOpOOTYZ8wMHzgv/PWrVtHbGzsBduai4pIfWQSFbt3Ezp1KgH33mN7UGwHJ5N0DPL9ipj7b4MOjp+0n7SdCdg7Xal/VT1Q9n9SyilCiKWc+RwBAFLKG2q4zUeFEGOBncBTUsqCGq5HUZqsNSlreG7jcwS4BfD5NZ/Txr+N3W2NWVmkjn8QfXIy4R+8j88111gfFPtgKwXZeuuDYhOegqh+DsepM5pxc9ae94wgzM/d4fUrtUNIec4x3jpDiJ5Syt1CiCHnmy+lXH/RlVvPCJZJKbvYPocCuVgTy+tY6xo8cIG2E4AJAKGhob0WL1580V/mfEpLS/Hycnw43Nqm4qoeFZeVlJK/i//mt8LfaOHSggdDHsRHe+7QzBeKS5uZif/MWYiyMooefghDhw7oiyWpf+ux6M0MC55DXr/bKfd0vBdPod7CrN16EossaAWYKx1qXDTw3y4uDAxrHOMINda/L3AstqFDh+6SUva+6IJSyvO+gDW2n+9caJmLvYCWwL7qzjv71atXL1lTa9eurXHbuqTiqh4Vl5QGs0G+svkV2WVBF/nk2idlhbGiWnGVx8XJw/36y8MDB8nyffuklFJmJBbKzx9fJec//LPMnnGzlEXptRJrQlqh7D99tezw4h9yxd6TcsnuNDnwrTWyxbPL5MC31sglu9NqZTu1pbH+fUnpWGzATmnHMbaqewTNbWcDNwghFmMdWqJyAtld3ewkhGgupcywfRwN7KtqeUVRrIoNxTy17im2ZmxlfNfxPNbjMbuLyQCUbthA2uQpOAUFETV/Hi5RUSTvzWXlZ3F4ksn1XVfg+9/54O7ncKzL92bw1I9xBHi48NPDA+gcZn0e9aYe4Re9d6E0jKoSwTSsN4kjgPc5MxFIrBXLLkgI8R0QCwQJIdKwPoQWK4SIsbVPBibWNHBFaSrSStJ4dM2jpBSn8NrA1xjddnS12hf++isZL76Ea7u2RM2di1NQEAc2pbPuG1tFsf478bjjC3B2bMx/i0Xy0Rrr6KG9Wvgz555eBHvbXwNZaThVJYIMKeU1QohpUsrXqrtiKeVd55k8v7rrUZSmLC47jslrJ2O0GPns6s8uOnqoNJvJX7CQ4E8+IW/SI0izhZz338ejf38iZs9C4+nJzmWJbFuWQqRLPCOH5+Ay6lPQODbcc7nBxFM/xPPHvkxu6xXBG6O74OqkhpC+VFSVCGYCvYCbgGonAkVRHPNn8p9M3TiVUM9QPh72Ma18qx7/0ZCcTNqUJzCkpKCpqCD7gw/BZMLzyiuImD0bnJzZ+M0+Ejbl0M5tHVfd5I92yFsOVxRLL6zgwYU7OZRZzIvXdWTc4FZq9NBLTFWJwCiE+BIIF0LMPHumlPLxugtLUZouKSXzEuYxc89MeoT04KOhH+Hv5n/Rdsl3j8FcWAgWWyVZkwmEQLdvPxahZfWnu0hMKKaH568MuKc/osfdDse6KyWfiV/tQm+0MP+/fRjaPsThdSr1r6pEMAoYjvVewEWHlFAUxXFGs5FXt7zKb4m/cW2ra3l90Ou4aF3sauvapg3l27efOVFKNG07sfSDrZxM0jPI7ytixt0Nba92ONafdqXxwi8JNPdzY/GE3rQJ8XZ4nUrDqKoeQS6wWAhxUEoZX48xKUqTVKQvYsraKezM2snD3R/m4e4PV+sSi0ef3uckAr1vM3YH3kRxUhlXh3xBuwlPQ0Qvh+I0WyTv/HmIuRuOM6hNIB/f3RM/D/uSldI4XbD/mRDiBwApZbwQ4p2z5v1V14EpSlNyovgE96y4h/iceN664i0eiXmkWkmg4Psf2LXsGAX+7f6dVuYRyvbOUygud2JU5FzaTX7T4SRQojMyfuEO5m44ztgBLVhwf1+VBC4DVV0aalvp/dXAs5U+B9dNOIrS9OzK2sXktZMRCD7/z+f0CrX/YG0xGMh6/Q0Kf/yRkHYd2RPzX0b4vYeTMLAx/38YcSPWbz6Rj38M3qEOxZmcW8b4RTtJzi3jzdFdGNOvhUPrUxqPqhLB+ceeuPg8RVHstDRxKS9vfplwr3A+HvYxUT5Rdrc1ZmWR/vhkKuLjCZwwgeA2JwjZ/hEr8p7HJF0AwXCfj2g/sKXDSWDzsVwe+db6DOlX4/oxoLUaPvpyUlUi8BBC9MB6+cjd9l7YXmq0KEVxgJSST+M/5dP4T+nTrA8fxn6Ir6vvxRvalO/aRdrkKVjKywn/v//DZ+QILIUZpGyYgVF6ANDd4zfa+2yHofMcivWrLcm8svQA0UGezLuvNy0Czz/ctXLpqvKBMuAD2/vMSu9PfVYUpQb0Zj3T/pnGiqQV3Nj6Rl4e8DLOWvsGX5NSUvDtt2S99TYu4eG0+PILXNu2pbzYwMrPjnOydBRa9MR4LmV/xX9oGRFCRA3PBoxmC68u3c/XW08wrEMI/3dnDN5ujWOQOKV2VdVraGh9BqIoTUG+Lp8pa6ewJ3sPk3tOZlyXcXbfFLbo9WS+8ipFS5bgNWQIYTPeRevjQ0ZiESs/2UlFmRFnYebawHeJcIojwv0gKxNeYcThAiLaX/w5hMoKygw88s1uthzPY+KQaJ4Z0QGtRj0kdrlStYcVpZ4cLzrOpNWTyC7PZsaQGYxseb4CfudnPHmStMcno9u3j6BHHiHo0UkgBAlrT7DphyN4abLpEpREq5tuJiKzF3JnPBH9ezCibTeyk4urlQiOZpUwbuFOMot0fHB7d27uGVGTX1e5hKhEoCj1YHvGdqasm4KzxpkvRn5B9+Dudrct27ad9ClTkAYDER/PxnvYMIwGM+sWxnFkVyEtXXcx7IpM3G58yzpwXMkzFB3dit+QZ4nw9q9WEvj7UBaPfxeHm7OWxRP70zOqemcSyqVJJQJFqWNLji7htS2v0cKnBbOHzSbC275v2FJKChYtIuvdGbi0aEHE7Fm4RkdTmF3On7O3k5dtop/3j/S6cxCiV6Xe3d7NiOsxndhq3BuQUjJ3w3He/vMQnZr78PnY3qqCWBNSVanKnlU1rEk9AkVpSizSwqw9s5iXMI8BzQfwXux7+LicW03svG0rKsh4aRrFy5bhNXwYYW+/jdbLi6T4HFbPj0MYKxgV8RUtHpgKYTEOxakzmnlhSQK/7E7nuq7NmXFbNzxc1HfEpqSqf+33q5h30XoEitKU6Uw6pm6ayl8pf3Fru1t5od8LOGvs63FjSEsj7bHH0R86RPCUyQROmIBEsG3JYXauTCfYKZmRvbfgc9dn4BHgUJzZJTomfrWLPScKeWJ4Ox4f1kaNHNoEqV5DilLLcitymfz3ZBJyE/hf7/8xttNYuw+upf/8w8knn0JaLETO+RSvIUPQlRpZ9dkOThzV0cF9DUNG+eE0dB5o7K9Qdj770ot4cNFOCsuNfDqmJ9d0be7Q+pRLl13nf0KILkAn4N8SRlLKRXUVlKJcqo4VHGPSmknk6/L5cOiHDIsaZlc7KSX58+eT/cGHuLZuTcTsWbi0aEF2SjF/zt5OWYmZ2MBv6PTfsYi2wx2O80LlJJWm6aKJQAjxMtaSk52AFcA1wCZAJQJFqWTzyc08te4p3JzcWDByAZ2DOtvVzlJWxskXX6Tkjz/xHjmSsDffQOPpyYGNaWz47hDu5HFz+98Ivf9t8LN/CIrzbqtSOcmeUX58dm9vVU5SseuM4FagO7BHSnm/ECIUcOyZdUW5hMV+H0ueLu/0hIWn37bzb8fHwz6mmWczu9ZlSEkh7dHH0CcmEvK/pwgYNw6zycL6L/dwYFsBES4J/OfKk7jf6HhN4crlJG/pGcH0m1U5ScXKnkRQIaW0CCFMQggfIBuIvlgjIcQXWIvbZEspu9imBQDfAy2xFq+/XUpZUMPYFaVBnJEEzrLomkV4Ots3Fk/phg2k/+9phBBEfj4Xr0GDKMnX8eesrWRnWOjp/Sv97uyHptcTDsd8qpzkwcxipl7bkfFXqHKSymn23G3aKYTwAz7HWqlsN7C96iYALADOfnTyOWCNlLItsMb2WVEuG/YkAWmxkDtnDqkTH8I5LIyWP/+E16BBpB7I44dXN1KYWco1EfMZ8OR4NL3GOBzTrpR8bpy9idT8cr64rw8PXhmtkoByhirPCIT1r+UtKWUhMEcI8SfgI6Xce7EVSyk3CCFanjX5Rqz3G8B6Qr2OM+scKMplzVxaRsbzz1GyajU+o0bR/PXXEK5u7Fp2hG3LTuDvlM41vTbjN2Y2uFf/qd5f96QzY+Vh0gsrCN/6N1e2C+LnXek093Pjuwd70zZUlZNUzlVlIpBSSiHEr0Av2+dkB7cXKqXMsK0rQwihKl0rlxSj2VjjtvrjSaQ99hiG5GRCnnuWgPvuw6Azs2bWVpIOVtDGbTNDr/fEZeicGnUN/XVPOs//kkCF0QxYLwd9tz2VtiGe/DBxIP6eqpKYcn5CyqprzAghPgYWSCl3VHvl1jOCZZXuERRKKf0qzS+QUp73a48QYgIwASA0NLTX4sWLq7t5AEpLS/Hy8qpR27qk4qqexhDXScNJFuUuIt2YfsFlZrWYdd7prvHx+Hy5AJycKHxwPMb27dEVSk6uL0df4Uw/38Vo+najMLDmTwk/ta6cPN25/58D3AQfxHrUeL21qTH8O55PY40LHItt6NChu6SUvS+2nD03i4cCE4UQKUAZ1sI0UkrZrQZxZQkhmtvOBppjvfF8XlLKucBcgN69e8vY2NgabA7WrVtHTdvWJRVX9TRkXGaLmUUHFjFrzyy8XbzxxpsSY8k5ywW6BZ4To7RYyP34E3I/nYNb585EzJqJc1gYR7adZO2P+3GRFdzY7nvCHngT/CIdijP/z+XnnV6gk43m31T9fVVffcRmTyK4pha39ztwH/C27edvtbhuRal1aSVpTN00ld3ZuxkWNYxpA6YR4HZ6WIeq/pOaS0o4+fQzlK5bh+9NN9HslZeRzi5s/CaevRvzaO58hBFD0vC86XNwcqwvf06JHldnDTqj5Zx5avA45WLsSQRvSCnvrTxBCPEVcO8Flj+1zHdYbwwHCSHSgJexJoAfhBDjgBPAbTUJWlHqmpSSJceW8M72d9AIDW8OfpPro6+3u7eN/tgx0iY9iiE9ndCXXsT/7rutVcTe20BGmoVuXn8w8O5eaHtOcjzOPem8tuwARpMFJ43AZDl9ecjdWcvTI9o7tA3l8mdPIjjj8UghhBbbzeOqSCnvusAs+565V5QGkluRy6ubX2Vd2jr6NuvL64NeJ8wrzO72xSv/IuP55xHu7rRY8CUevXtz8kgBKz/dgUFn4uqwH2k34Ulo1tWhONMLK5i6JIF1h3Po3cKft2/pxr70otO9hvzceXpEe27qEe7QdpTLX1XDUD8PvIC1cH3xqcmAAdu1e0W53KxOWc1rW16jzFjGM32eYUzHMWiEfT14pNlMzkczyZs7F7fu3YiYOROnkBDiVx7jn1+T8dFkc0PPzQTe8yG4+118hRdgsUi+2ZbC238cQgKvXN+JsQNaotEI2oR4cVOP8EZ9zVtpfKoaffQt4C0hxFtSyufrMSZFqXclhhLe2vYWS48vpVNgJ6YPnk5rv9Z2tzcXFpL+v6cp27QJv9tuI/SlFzFZNKz6ZCtHEypo5bqTYTe44jp0tkOjhibllvHsz3vZnpTPFW2DmD66K5EBjaNHkHLpsufS0HYhhK+UsgjA9pRxrJTy17oNTVHqx7aMbbz4z4vklOcwsdtEJnafeNHaAdJsJn/BQoI/+YTM22+jZNVqTFlZNHv1VfzvuJ3CrHL++OgfCvIF/f1/pef4WxCtaz6yu8lsYf6mJD5YdQRXJw3v3tqN23pFqCeElVphTyJ4WUq55NQHKWWhbURSlQiUS5rOpOOj3R/x9cGvaenTkkXXLKJb8MV7RRuSk0mb8gSGlBQ0FRUU2J4PCHv3HXyvvZbjuzNZ/cVetOZyrm+7jMjxr4Jvza/TH8wo5tmf97I3rYj/dArljZu6EOLj2AB0ilKZPYngfOexqo6dcknbn7uf5zc9T1JREnd1uIsnej2Bu5N93SyT7x6DubAQLJW6alosZLzxFgdKIti9NpcQpyRGXpWK9+hPwKlmT/TqTWY+XpvIJ2uP4efhzMd39+Tars3UWYBS6+w5oO8UQnwAfIy1ROVjWAefU5RLjtFiZF7CPObGzyXAPYDPrv6MgWED7W4vjUY0Xl6Y8/PPmG7QenCo43hy1+bSyfNvrhjTDaeeE2sc554TBTzz016OZpdyc49wXhrVSQ0RodQZexLBY8BLWIePBvgLeLHOIlKUOpJUlMQLG19gX94+rou+juf7Po+vq/2Vucq2bSfrjTc4Jtvh4++Kf8FRAIq9WxDf9SFMTh4Mbf49nR56HELtK0pztgqDmff+OswX/yTRzMeNL//bh6Ed1JBcSt26aCKQUpYBzwkhvKSUpfUQk6LUKou08N2h7/hw14e4Obnx3pD3GNFyhN3tjVlZZL/zLsUrVuAcFkarAU6szxxH5/3zqXAP4XDb20FouCLoBzo9OwPcalb2cXNiLs/9nMCJ/HLu6R/FsyM74O1mX8F7RXGEPaUqB2KtSOYFRAkhugMTpZSP1HVwiuKozLJMXvrnJbZmbGVw+GBeG/gawR7BdrWVBgP5ixaR88mnYDIRNGkSgePHoTEVYZn+KOucJ2HBGQ1GRvjPIPq5hTVKAsU6I2+tOMR320/QMtCDxRP60z86sNrrUZSasufS0IfACKzjBCGljBdCXFmnUSmKg6SULE9azvSt0zFJE9MGTOPWtrfafaO1dNM/ZL35JoakJLyuuorQ55/DJTISXamR7UvS2VcwEYF1uOceXr8RPagzeNtXnrKyCCxlkwAAH0ZJREFUNQezmLpkH9klOiZcGc0Tw9vh7qLKRyr1y67eP1LK1LP+A5nrJhxFcVyBroDXt77OqpRVxATHMH3wdCJ97BvZ05ieTtbb71CyahXOUVFEfjYHryFDMJstxK88yo5lxzEYtbRw202mvgNdPH5lX/lIIsK7EVGNGPNK9by27AC/xZ2kQzNvPru3F90ja/60saI4wp5EkGq7PCSFEC7A48DBug1LUWpmQ9oGXt78MoX6Qib3nMz9ne9Hq7n4N2yLXk/e/Pnkzf0cgOApUwi4/78IFxeS47P559s9FBY5E+mynzaddGxJGcSIfpsIT15MeNdIVn4TyAgPfyLaV11VTErJ0r0ZvPL7fkp0Rp4Y3o6HY1vj4lTzp40VxVH2JIKHgI+AcCAdWAk4NmSiotSycmM5M3bO4KcjP9HWvy1zhs+hfYB9o26WrF1L1vS3MKam4j1iBKHPPoNzWBh56aVsWvQPaSngp83iurZbaHH7OPYcCGbECB8iwtpTOG85EbeMZ8QAF7KTi6tMBJlFOl78dR+rD2bRPdKPd2/pRvtmqnSk0vDs6TWUCzheQVtR6sjurN1M3TSV9NJ07u9yP4/GPIqL9uJ97g0nTpA1/S1K163DJTqaqC/m4zlwIBUlBv6Zv5UDO0pxEWUMDllFl1uHo+36CQhBz0pXmeJ6TCfWO5SI9lwwCUgp+X5HKm+uOIjRbGHqtR15YHArtBr1YJjSONjTayga6xlBf6wPlG0BnpBSHq/j2BSlSgazgY/jPubLfV8S5hXGlyO/pFfoRUdIx1JRQd7nn5M3bz7CyYmQp58m4N57sAgndi87wK4/UzGaNHT1XkOfa6Nwu+L9GheOOZFXzvNL9vLPsTz6Rwfw9s3daBnkWaN1KUpdsefS0LdYnyoebft8J/Ad0K+uglKUizmcf5gXNr3AkYIj3NL2Fp7u8zSezlUfYKWUlKxaRdbbb2M6mYHPqFGEPP00TiHBJO06yT+L91Jc6koL170MuqIM/1FPg2fNunGaLZIFm5N5b+VhtBrBm6O7cFefKDTqLEBphOxJBEJK+VWlz18LIR6tq4AUpSpmi5kF+xcwO242vi6+zL5qNkMih1y0nf54EllvvEHZ5s24tmtH+Ffv4NGnDzknitn0+ipOnnQiwCmL67skEHXHRAhuV+MYj2aV8MzPe9lzopCrOoTw5uguNPdV5SKVxsueRLBWCPEcsBjrpaE7gOVCiAAAKWV+VY0VpaZiv48lT5d3esLC02+HRw3npQEvnVE/+HzMpWXkzfmUvIWL0Li5ETp1Kv533Ul5mZm/P93IwXgDbqKcIeHr6HTnDWja/rfG8RrNFuasS2TW38fwdNXy0Z0x3NA9TA0SpzR69iSCO2w/zx5B6wGsiSG6ViNSFJszksBZPoj9oMoDrJSS4hUryH7nXUzZ2fjefDMhTz4BPv7sWprA7tVZmM0Q47eO3jd2xrXf+2BHN9NTft2Tfrok5Na/ubNvJCsSMjmYUcyobs155YbOBHk5VpBeUeqLPb2GWtX2RoUQyUAJ1gfTTFLK3rW9DeXyVlUS0B05QtYbb1K+fTtunToRMfMj3Lp359iWFLb8tJ2ScjdaucUx8D9O+I2cCq5e1dr2r3vSef6XBCqM1ucq0wsreP+vI3i7apl7by/+07n6TxgrSkOqqmZxHyBVSplp+zwWuAVIAV6phUtCQ21dUxXlDEX6In468lO125lLSsidPZv8r79B6+VFs1dewe+2W8lOKeaPV/4gI8uNQKcMbuydTMRtE8G3Os8CnzZj5eF/k0BlXm7OKgkol6Sqzgg+A4YD2MYWehvrkNQxWIvX31rn0SlNSlpJGl8f/Jpfjv5ChanC7nbSYqHo99/Jfu99zHl5+N1+O8FTJqPDgzUfr+PwAYG7RsfQluvoMOYONJHjahxjVrGO9MLzx5ZZpKvxehWlIVWVCLSVvvXfAcyVUv4M/CyEiHNwuxL4Swghgc+klHMdXJ9yCYvLjmPRgUWsObEGDRquaXUNYzuP5balt120re7AATJff4OKPXtw696NyDlzcGrXgd2/7GLPhgKkRdIzcCO9bu2HS8w7UIMbt1JKdiQXsHBLMiv3ZV5wuTA/1TNIuTQJKeX5ZwixD4iRUpqEEIeACVLKDafmSSm71HijQoRJKU8KIUKAVcBjp9ZdaZkJwASA0NDQXosXL67RtkpLS/Hyqt414PrQ1OOySAt7y/fyd8nfJOmTcNe4M8hrEEO8h+DnZB187YXUFyixlJzT1lvjzVsBU/H67XfcN25EenlSMvpmKvr1ozSpgpx4MxVGL6LdtxHeqYSC1kORFylGfz56k2RLhok1J0yklljwcIIrI5zwcxX8ctSIoVKlShcN/LeLCwPDGkf9gKb+91VdjTUucCy2oUOH7rLnHmxViWAqcC2QC0QBPaWUUgjRBlgopRxUo8jO3c4rQKmU8r0LLdO7d2+5c+fOGq1/3bp1xMbG1iy4OtRU4yo3lrPk2BK+PvA1aaVphHuFc2+nexndZjQezh7nLC/NZvIXLCTzk09oNukR/O+9l6IlS8j54EPMxcX4jxlD8GOPkpOpY+PCbWTneRLsnMjgfrmE3TwRPKruXno+ybllfLU1hR92plKiM9GxuQ/3DWjBjTHh/w4RfUavIT93nh7Rnpt61LxAfW1rqn9fNdVY4wLHYhNC2JUILnhpSEr5phBiDdAc+EuezhgarPcKahqYJ6CRUpbY3v8HeK2m61MuDVllWXx76Ft+PPIjJYYSugd358neT3JV5FUXHB3UkJxM2pQnMKSkoKmoIOejmeTMnIXU6XDv3YtmL72EITCKNZ+v4+hRVzw0eoa130n7MfcgQtpWKz6zRbL+SDYLN6ew/kgOThrBNV2bc9+AFvRq4X9OL6WbeoRzU4/wRn0AURR7Vdl9VEq59TzTjji4zVBgie0/lhPwrZTyTwfXqTRSh/MPs3D/Qv5I/gOLtDAsahhjO40lJiTmom2T7x6DubAQLNZrMFKvB0Dj6UnY51+y+8etxG1JBgm9Q9fR4+5huLR/tVrxFZYb+GFnKl9vPcGJ/HJCvF15Yng77uobSYiPW3V/XUW5JNlVmKY22Qar617f21Xqj0Va2JS+iUUHFrEtYxvuTu7c0f4OxnQcQ6S3fQViTLm5pIRfhaclHv/Co/9Oz/drR1bnUWx+diXlBg/aeu1hwA0t8B78crUeCNuXXsSiLcn8FncSvclC31YBPDuyA//pHIqzVtUGUJqWek8EyuVLb9azLHEZiw4s4njRcULcQ5jScwq3trsVX9eL1/I1l5ZSsno1xcuWU7ZlC57e0ezrNI4uB+bjX3iUE+GxJLa5GSm0hHKUa4YZaHbjY+Bi32ieBpOFP/ZlsGhLCrtSCnB31nJLrwjGDmhBh2Y+jv76inLJUolAcVi+Lp/vD3/P4kOLydfl0yGgA9MHT2dky5E4a6vuRWMxGCjbsIGiZcspXbsWqdfjHB5O4PjxtHCOh8XzSeg8ARdDEeWezXE2FHNl5FLaP/kawte+m7OZRTq+3ZbCt9tTyS3V0zLQg5dGdeLWXhH4ujeOXj6K0pBUIlBqLKkoia8OfMXvib+jN+u5IvwK7ut8H32b9a16HCCzmfIdOylevozilX9hKS5GGxCA36234jPqOly6dOPk/pMcWNeSpG5mTBYPTM4ehDsncF3zGTg/uQu8Q6uMTUrJtqR8Fm1JZuX+LCxSclX7EMYObMkVbYLUcNCKUolKBEq1SCnZmbWThfsXsj5tPS4aF65vfT33drqX1n6tq2yn23+A4mXLKF6xAlN2NhoPD7yvHo7PqFG49OxDWvwJ9q05SPK8DPRmN5yEMyHOaRiNrnRx/5ODuuFktXiYiCqSQJnexJI96SzaksyRrFJ83Z0ZP7gVY/q1ICrw3O6piqKoRKDYyWgx8lfyXyzcv5CD+Qfxd/Xn4e4Pc0f7Owh0v3DxFkNyMkXLl1O8bDmGpCRwdsbriivwHXUdzv2v4MTu48RvOM6J79ZjsjjjKqCVz16iO7qhadGH1cu6cq3PK0Q47SbKYz8rE15hxOGCc8pCHs8p5autKfy0M40SvYnOYT68e0s3ru8e9m/ff0VRzk8lAqVKJYYSfjryE98c/Ias8ixa+rRk2oBpXB99PW5O5+9eaczOpuTPPylaugxdQgIIgUefPgTc/1+c+seSEpfInm0nSf99Kxa0eGpMdAzYQ6uufoRdeSXaZjeAEOxemcKICT5EHO2B3LmHiP49GNG2279F4s0WydpD2SzckszGo7k4awXXdm3O2AEt6Rnlp+oAKIqdVCJQLlgAxk3rhkZoKDeV07dZX6YNmMbg8MFoxLndK80lJZT8tYri5cso27oNLBZcO3Uk5JlnoO8QUvensCsuj4y/9wIafLU6ujc/RnTPUEIHXYUIOHdcoZ4jWvDrnnQeSRjIVMs/vJEwiHHh5QwZ3Jw56xP5emsKaQUVNPNx46mr23Fn3yiCvVUNAEWpLpUIlAsWgNGZdVwXfR1jO42lU2Cnc+Zb9HpK16+neOkyStevRxoMOEdFEThxIpbeQ0g9msGO/WXkbk8GINC5mL4tkoju14KAfiMRF6kHfHrcf3fuYBoY4Kkf4gGJWUL/6ACmXtuR4Z1U339FcYRKBE2Y0WLkSH7VD4q/fcXbZ3yWZjPl27ZRtGw5JX/9haW0FG1QEL6334EhZhDpaUVsPWym6GAe4EQz11wGtk8hekA7fHveYXeff4B3/jx0zrj/ZinxcNGy5JFBtG/mbfe6FEW5MJUImpACXQHxOfHEZccRnxPPvtx96MwXH0NfSolu3z6Kli6l+I8/MOfkovH0xPPqEVR06U9GjoHNSU6ULTOjwY1wj6PEdEul1eBueHZ+AC7yLAGAyWzhSFYpcamFxKUWEJdaSMYFxvevMJhVElCUWqQSwWXKbDGTWJR4xoE/pTgFACfhRMfAjtza7la6h3Tn6fVPIyyS67ZbGL1VsmSAYHkfDc0KIGfmLIqWL8OYcgLh7Iz7lUMp7dSP9GItKWme6DZ54IQTkT7HaB3jTIsre+PW+pGLjvufWaQjLrWAPamFxJ0oJCG9iHKD9du/n4czPSL9yCjSUaIzndNWjfuvKLVLJYLLRImhhIScBOJyrAf9vTl7KTWWAhDgFkD34O7c3PZmugd3p3Ng5zN6/Pw27w+G70wiOv0obka4a51kRHw0Oo8W5KZ9inO/weiH3kamzp0T2f4YE9xwEWW09D9OdBcvomIH4xx2zQVjKzeYSEgrsn3bL2TPiUIyi63f9p21gk5hvtzeO5KYSD9iIv1oEeiBEOKc2sAA7s5anh7Rvo72oqI0TSoRXIKklKQUp1i/7efEEZcdR2JhIhKJRmho69eWa1tdS0xIDDHBMUR4R5zTldJSXo7+eBKGxGPcsTaJI+3H0bxwPm6FRyn1acvhDuMJydrCsZteIK0gGEuyM+6aQtoFHyM6Jojw2KFo/a8/JzaLRZKYU8oe2wE/LrWQI1klmC3WUcyjAjzo2yrAetCP8qNTcx/cnM/fz//U+P6Nedx/RbkcqERwCagwVbAvdx/xOfHEZ8cTnxNPgb4AAG9nb7qFdGNEyxF0D+5O16CueLmcrmZkLi6mIi4OQ2Ii+mOJ6I8nYjiWiPHkyX+XCQHEkW9J6DwBn+JkCvzbIYWGkxFX412cQ9fII0T3jqTZFVeh8fA7I7bcUj1xtgP+ntQC9qYWUaK3Xs7xdnMiJtKP4R1b//ttP9Cret071bj/ilL3VCKoRxfqrx/oFsi6O9YB1m/7GWUZ/17bj8uJ40j+EUzSenBt5duKIZFDiAmOISYkhla+rRAIzPn56I8lYli/lMxjiegTEzEkJmLKybGuF4HBOwRjVDsMba5E19WPcq0nZWZPinU+GIX1unt+YCecDSWEZ26g09BAWk5+AuFinaczmtmfUmC7vGP9mVZgLeSu1Qg6NPPmhpgwYiL96BHlR3SQlxrTR1EuASoR1KML9dfP0+WxcP/Cf7/xZ1dkA+Du5E7XoK7c3+V+YkJi6BrYFa8ivfWAvzUR/bFFnLAd8M2FhViEExVuAej9wjCEtEDfvjMVXX0pw4dSix9mTvfe0eiN+Djl4uNRjI9nFvp9BWSF9KVZ5jayQ3rhU5zI2q7/xX9/3r+XeA5mFGM0Wy/xhPm6ERPlx9gBLYiJ9KdruK8aykFRLlEqEdSDClMFeRXWJHC+3jlSI3hv53uEe4XTu1lvYgK7EWMOo3muBfPxJPTrjqNP/JjsxETSdCYq3IOpcAtC79ccnV83dF1iKdf6US78kJx+sMpZVODjnEegZyHRvgX4Brvi0zwA36gwvFq2xuLmR36ZgXEzNhEbLhnt9wbREXEc1/dgietUfl2VRqrzCTxctHSL8GXc4Oh/v+2HqupdinLZUImghsqN5eRV5JGny7vgz9yKXPIq8ig3lQMw+OhVZ/TOuX2jpP+xaOLbtuS+9lE4JWShSzxK2YlNFGl9yXQPosItGJ1PKDrPEZT3CkCvObOAirumCB/XQsK98/DxL8A7yB1NkB/6gBDyPFqTq9NwvFRPbqme3FIDual68g6Vk1u6i4JyA1JCX50Tv3pYCNEG0VIKNmkD+dXDRHOzhs+fHkTbEG+06hKPoly2LstEEPt9LPnluf9+837e9s07wCPo32vxZ5NSUmYsu/CB/axpFaaKf9sKKXEzgIcOmlm8aSa9aW/xZJDJF39jCD5GLV56DeYNSRzuMI6QooXI8iyyg2PIiLqOjicT2JYp0Hl2pSI4FlNo5X7yFry0+Xi7lRDgmYHWKxOjpytFHt6kuQaw3xxBTlmo9SBfpKck61S/+xzby8rTRUuQtytBXq60CvKkd8sAgrxcCfZy4cNVR8gvtzDTNJp2mjRmmW4mx9mCJdhVVe5SlCagQRKBEGIk8BGgBeZJKd++SJNq6ZjQ/Zxv3r2SolnTqxXfdvv2jAN6fkU++eW5lBfmoS3X4akDT73EQ289sHtVQKDZmy5Gb7xM7ngYPXEx+ONkckYYnZBmLRazM2atC2atG2atGyYnV8xaV8xaN3RaV8q0bpidXDF1sk7f2/3RM+LNDo7BxyUPd9cS3F2PU+EMOVonEnFnl8GLLJMn4AkGIN/2Avw8TAR5lRPo6UKnMB+CvVwJ8nIhyMuVwErvg7xcq7x+7+3mzPO/JJBj9OcOwzRA9ddXlKak3hOBEEILfAxcDaQBO4QQv0spD9TWNm7emERS63EEli7GqMujwLc1ueGjuGb3WnKP/o2v2Y1gixta6Y2QwSC6Yda6YrIdvM1Orpi0braDuSsIDQYXyHf59xh8Xlqpxwk9TkKPRujRaAxoNAaEpgKNxkyWMZfOSRWUe7aiIKADoZnbCE9fxqLBev5wfRWN8CXAw3oAD7Z9e7/e0+Xfb/KVD+wBni64ONXOQGuqv76iNG0NcUbQFzgmpTwOIIRYDNwI1FoiKHQ9SpcD84nvOglLpXFuUqOuPWM5rUWPVurRSh1a9GiEASEqcBIlaDRmjBoLUivRaQV6jUAnBOVoqRAa9Fpn9E6u6J1d0Dt7YHDxRLh44ursi6uTBlcnDW7OWtt7La7OGn7acIhZRYvJDh1My+QVpIddgX/uLjZo7mTXi8Px93BpsO6Wqr++ojRdQkpZvxsU4lZgpJRyvO3zvUA/KeWjZy03AZgAEBoa2mvx4sV2b+OH5ZMYv9JCSsu7OBk2mNDMbYRkr+OX/ka8ox7E4uqCydkdrZMzzhpw1oKzRljfa8BFa3uvPT3NWSNsy4FWUKOiJzNW6ogtdCbI+zfuc1/Egoqx5JXcyDo/I0+PaBy9cEpLS/Hy8rr4gvVMxVU9Kq7qaaxxgWOxDR06dJeUsvfFlmuIM4LzHUHPyUZSyrnAXIDevXvL6nxLfeaw4JatbckJ6n76m3feFja3T2f7xLtqGrfDUpP2svBQBkWiL53kNj7V9MPXx8R97cKIje3WYHFV1ljPCFRc1aPiqp7GGhfUT2wNkQjSgMhKnyOAkxdYtkbCTD3Z0/1WVrX7kpNDjhFWlMjVHuMJM/1Um5uptnvv74b3nkBmrDzMHYXTCPdz51F1LV5RlAbWEIlgB9BWCNEKSAfuBO6uzQ1Ma/M6R1sYcE0Yw7zcGbzhfg+RIzsyzfn12txMjahr8YqiNDb1ngiklCYhxKPASqzdR7+QUu6vzW30HNGCnsAdo9qybl0QS9UBV1EU5YIa5DkCKeUKYEVDbFtRFEU5k6r4rSiK0sSpRKAoitLEqUSgKIrSxKlEoCiK0sTV+5PFNSGEyAFSatg8CMitxXBqi4qrelRc1aPiqp7GGhc4FlsLKWXwxRa6JBKBI4QQO+15xLq+qbiqR8VVPSqu6mmscUH9xKYuDSmKojRxKhEoiqI0cU0hEcxt6AAuQMVVPSqu6lFxVU9jjQvqIbbL/h6BoiiKUrWmcEagKIqiVOGyTQRCiNeFEHuFEHFCiL+EEGG26UIIMVMIccw2v2c9xzVDCHHItu0lQgg/2/SWQogKW7xxQog5jSEu27znbfvrsBBiRD3HdZsQYr8QwiKE6F1pekPvr/PGZZvXYPvrrDheEUKkV9pH1168VZ3GM9K2T44JIZ5ryFgqE0IkCyESbPtoZwPG8YUQIlsIsa/StAAhxCohxFHbT/862biU8rJ8AT6V3j8OzLG9vxb4A2uBnP7AtnqO6z+Ak+39O8A7tvctgX0NuL8uFFcnIB5wBVoBiYC2HuPqCLQH1gG9K01v6P11obgadH+dFeMrwP8aah+dFYvWti+iARfbPurU0HHZYksGghpBHFcCPSv/XQPvAs/Z3j936v9lbb8u2zMCKWVxpY+enK6CdiOwSFptBfyEEM3rMa6/pJQm28etWAvzNLgq4roRWCyl1Espk4BjWOtO11dcB6WUh+tre/aqIq4G3V+N2L+1yqWUBuBUrXLFRkq5Acg/a/KNwELb+4XATXWx7cs2EQAIId4UQqQCY4BptsnhQGqlxdJs0xrCA1jPTk5pJYTYI4RYL4S4ooFigjPjakz762yNZX9V1tj216O2y31f1NllBfs0tv1SmQT+EkLsstVKb0xCpZQZALafIXWxkQapR1BbhBCrgWbnmTVVSvmblHIqMFUI8TzwKPAydtZMrsu4bMtMBUzAN7Z5GUCUlDJPCNEL+FUI0fmsM5uGiKtR7K/zaBT763zNzjOtzrrmVRUj8Cnwum37rwPvY03yDaFe90s1DZJSnhRChACrhBCHbN/Om4xLOhFIKYfbuei3wHKsiaDOayZfLC4hxH3AKGCYtF38k1LqAb3t/S4hRCLQDqi1m1c1iYtGsL8u0KbB99cF1Pn+qszeGIUQnwPL6ioOO9TrfqkOKeVJ289sIcQSrJexGksiyBJCNJdSZtguYWfXxUYu20tDQoi2lT7eAByyvf8dGGvrPdQfKDp16lVPcY0EngVukFKWV5oeLITQ2t5HA22B4w0dF9b9dacQwlVY60y3BbbXV1wX0tD7qwqNZn+dde9rNLDvQsvWg39rlQshXLDWKv+9AeMBQAjhKYTwPvUea6eJhtxPZ/sduM/2/j7gQmeijmnoO+V19QJ+xvoPuhdYCoTbpgvgY6w9GBKo1OOjnuI6hvVaaZztdao30y3Afqy9KXYD1zeGuGzzptr212HgmnqOazTWb5N6IAtY2Uj213njauj9dVaMX9n+xvdiPaA0b6hYbPFcCxyx7ZupDRlLpZiibX9D8ba/pwaLC/gO6yVPo+1vaxwQCKwBjtp+BtTFttWTxYqiKE3cZXtpSFEURbGPSgSKoihNnEoEiqIoTZxKBIqiKE2cSgSKoihNnEoEiqIoTZxKBEqTI4SYahtC+tQw5f2EEOsqD0EshOgthFhnex8rhCiyjWt0SAjx3lnru0kIMY3zEEKU2n5WHjY7XgixWQjR3javqxBiQV39vopyMSoRKE2KEGIA1mE0ekopuwHDOT0YWogQ4poLNN0opewB9ABGCSEGVZr3DPCJHZtPlFLGSCm7Yx1J8gUAKWUCECGEiKr+b6QojlOJQGlqmgO50jpWEVLKXGkbawaYAbxYVWMpZQXWJ6/DAYQQ7QC9lDLX9rmVEGKLEGKHEOL1KlblAxRU+rwU67ALilLvVCJQmpq/gEghxBEhxCdCiCGV5m0B9EKIoRdqbBvKuS2nByUbhHWIi1M+Aj6VUvYBMs9q3tp2aSgReBL4oNK8nUBjGUpbaWJUIlCaFCllKdALmADkAN8LIf5baZE3OP9ZwRVCiL1YD+7LpJSnDvLNbes5ZRDWMWPAOtZPZacuDbUGpgBzK83LBsKq/xspiuNUIlCaHCmlWUq5Tkr5MtY6FbdUmvc34Ia1jGllG233FLoCDwshYmzTK2zLn7EJO8L4HWtpwlPcbOtSlHqnEoHSpAgh2p81RHkMkHLWYm9ivQF8DinlEeAtrEN2AxwE2lRa5B9OX+sfU0Uog7GOwnlKOxrX8MdKE6ISgdLUeAELhRAHbJd6OmEt8v4vKeUKzrzcc7Y5wJW2egMbgB5CiFMVuCYDk4QQOwDfs9qdukcQD0wHxleaNxRr8SRFqXdqGGpFcZAQ4iNgqZRydQ3buwLrgcFSSlOtBqcodlBnBIriuOmAhwPto4DnVBJQGoo6I1AURWni1BmBoihKE6cSgaIoShOnEoGiKEoTpxKBoihKE6cSgaIoShP3/zcwSZZN3PN/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2cfc6c193c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(SNR_dB, sum_rate_CE_PC/5, marker='o', label='gACE - Partially Connected Archi')\n",
    "plt.plot(SNR_dB, sum_rate_CE_FC/2, marker='v',label='gACE - Fully Connected Archi')\n",
    "plt.plot(SNR_dB, sum_digital_rate/5,marker='s', label='Digital Precoder')\n",
    "plt.plot(SNR_dB, sum_rate_CE_OMP/5,marker='p', label='gACE - OMP Init')\n",
    "# plt.plot(SNR_dB, OMP,marker='D', label='OMP Precoder')\n",
    "# plt.plot(SNR_dB, quan_OMP, label='quan_OMP')\n",
    "# plt.plot(SNR_dB, (sum_rate_CE_FC-sum_rate_CE_PC)/5, label='FC')\n",
    "plt.plot(SNR_dB, sum_rate_coordiante_descent,marker='x', label='Coordiante Descent from [6]')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel('SNR(dB)')\n",
    "plt.ylabel('Spectral Efficiency (bits/s/Hz)')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model['Wopt'] = mat['Wopt']\n",
    "model['Fopt'] = mat['Fopt']\n",
    "FRF_omp, FBB_omp = OMP( model['Fopt'][:,:,channel_index], model['Nrf'], model['At'][:,:,channel_index] ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linalg.norm(FRF_omp@FBB_omp,'fro')\n",
    "print(\"Ns\", model['Ns'] , \"   Frobenius norm FRF*FBB=\", x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "SNR_dB = np.arange(-30,15,5)\n",
    "OMP_rate = np.zeros((SNR_dB.shape[0]))\n",
    "sum_OMP_rate = np.zeros((SNR_dB.shape[0]))\n",
    "for channel_index in range(100):\n",
    "#     print('index: ' + str(channel_index))\n",
    "    FRF_omp, FBB_omp = OMP( model['Fopt'][:,:,channel_index], model['Nrf'], model['At'][:,:,channel_index] )\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        OMP_rate[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],FRF_omp@FBB_omp)\n",
    "    sum_OMP_rate = sum_OMP_rate + OMP_rate\n",
    "#     print(OMP_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_OMP_rate/100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "SNR_dB = np.arange(-30,15,5)\n",
    "OMP_quan_rate = np.zeros((SNR_dB.shape[0]))\n",
    "sum_OMP_quan_rate = np.zeros((SNR_dB.shape[0]))\n",
    "for channel_index in range(100):\n",
    "#     print('index: ' + str(channel_index))\n",
    "    FRF_omp, FBB_omp = OMP( model['Fopt'][:,:,channel_index], model['Nrf'], model['At'][:,:,channel_index] )\n",
    "    Frf_quantize, Fbb_quantize = quantize_Frf(model,nb_bit,FRF_omp,channel_index)\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        OMP_quan_rate[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],Frf_quantize@Fbb_quantize)\n",
    "    sum_OMP_quan_rate = sum_OMP_quan_rate + OMP_quan_rate\n",
    "#     print(OMP_quan_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_OMP_quan_rate/100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize condition\n",
    "x = np.linalg.norm(Frf_quantize@Fbb_quantize,'fro')\n",
    "print(\"Ns\", model['Ns'] , \"   Frobenius norm FRF*FBB=\", x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ps_system_proba_omp = initialize_pmf_omp(model,nb_bit,Frf_quantize,phase_value,3)\n",
    "ps_system_proba_omp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_CE_OMP= np.zeros((SNR_dB.shape[0]))\n",
    "sum_rate_CE_OMP = np.zeros((SNR_dB.shape[0]))\n",
    "for channel_index in range(2):\n",
    "    print('index: ' + str(channel_index))\n",
    "    CE_result = CE_beamforming(nb_elite,S,nb_iterate,model,nb_bit,nb_state,channel_index,10**(10/10),'FC')\n",
    "    for i in range(SNR_dB.shape[0]):\n",
    "        rate_CE_FC[i] = data_rate(model['Ns'],10**(SNR_dB[i]/10),model['Wopt'][:,:,channel_index],model['H'][:,:,channel_index],CE_result['Frf_max']@CE_result['Fbb_max'])\n",
    "    sum_rate_CE_FC = sum_rate_CE_FC + rate_CE_FC\n",
    "    print(rate_CE_FC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}