dSCAdemo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "incorporate-provision",
   "metadata": {},
   "source": [
    "# demixed Sparse Component Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "nervous-inside",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "import dill\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bizarre-hierarchy",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "mpl.rcParams['font.size'] = 9\n",
    "mpl.rcParams['lines.linewidth'] = 1\n",
    "mpl.rcParams['savefig.transparent']=True\n",
    "mpl.rcParams['font.family'] = 'sans-serif'\n",
    "mpl.rcParams['font.sans-serif'] = ['Arial']\n",
    "mpl.rcParams['savefig.bbox']='tight'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "jewish-mortgage",
   "metadata": {},
   "outputs": [],
   "source": [
    "kurtosis_centered=lambda x: ((x**4).sum()/x.size)/((x**2).sum()/x.size)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "005e5ee8-e8da-4f8b-9b77-e7af01a231a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "names=['Sample','Memory1','Distractor','Memory2','Test'];\n",
    "names=np.array(names)\n",
    "def DrawLines(kk=[0.51809945, 1.03684093, 2.05390537, 2.57065977, 3.59486604 ,4.5],\n",
    "              scale=4.5,oscale=0,linetype='--w',adjustment=0,axis=['x','y'],\n",
    "              discrete=False,epochs=np.arange(5),colors=np.array(['C0','k','C3','k','gray']),\n",
    "              names=['S','M1','D','M2','Test'],rotationx=0,rotationy=0,lw=1,shade=[0,2,4],fontsize=None):\n",
    "    if not fontsize:\n",
    "        fontsize=int(round(mpl.rcParams['font.size']*3/4))\n",
    "    scale=plt.gca().get_ylim()[1]*0.95\n",
    "    kk=np.array(kk);\n",
    "    if discrete:\n",
    "        kk*=10;\n",
    "    kkc=(kk[1:]+kk[:-1])/2;\n",
    "    epochs=np.array(epochs)\n",
    "    names=np.array(names)\n",
    "    bgcolors=np.array(['none',None,'none','none','none'])#np.array(['C0',None,'C3','none','k'])\n",
    "\n",
    "    ylim=plt.gca().get_ylim()\n",
    "    ylen=0.03*(ylim[1]-ylim[0])\n",
    "    xlim=plt.gca().get_xlim()\n",
    "    xlen=0.03*(xlim[1]-xlim[0])\n",
    "    borders=list(set(np.concatenate((epochs[:-1],epochs[:-1]+1))))\n",
    "    if lw>0:\n",
    "        for x in kk[borders]:\n",
    "            if 'y' in axis: plt.plot([xlim[0]+xlen,xlim[1]-xlen],[x+adjustment,x+adjustment],linetype,lw=lw);\n",
    "            if 'x' in axis: plt.plot([x+adjustment,x+adjustment],[ylim[0]+ylen,ylim[1]-ylen],linetype,lw=lw);\n",
    "    shadecolor=['C0','gray','C3','gray','k']\n",
    "    if shade:\n",
    "        for x in shade:\n",
    "            plt.fill_betweenx(ylim,[kk[x]]*2,[kk[x+1]]*2,facecolor=shadecolor[1],edgecolor=None,alpha=0.15)\n",
    "    \n",
    "    \n",
    "    if 'x' in axis: \n",
    "        plt.xticks(kkc[epochs],names[epochs],rotation=rotationx,fontsize=fontsize)#('sample','delay1','distractor','delay2','test')\n",
    "        for ticklabel, tickcolor , bgcolor in zip(plt.gca().get_xticklabels(), colors[epochs],bgcolors[epochs]):\n",
    "            ticklabel.set_color(tickcolor)\n",
    "            if bgcolor==None: continue;\n",
    "            ticklabel.set_bbox(dict(boxstyle='square, pad=0.001',edgecolor='none'))\n",
    "            #print(ticklabel.get_bbox_patch())\n",
    "            ticklabel.set_backgroundcolor(bgcolor)\n",
    "        plt.tick_params(bottom=False)\n",
    "    if 'y' in axis: \n",
    "        plt.yticks(kkc[epochs],names[epochs],rotation=rotationy,fontsize=fontsize)#('sample','delay1','distractor','delay2','test'))\n",
    "        for ticklabel, tickcolor , bgcolor in zip(plt.gca().get_yticklabels(), colors[epochs],bgcolors[epochs]):\n",
    "            ticklabel.set_color(tickcolor)\n",
    "            if bgcolor==None: continue;\n",
    "            ticklabel.set_bbox(dict(boxstyle='square, pad=0.001',edgecolor='none'))\n",
    "            ticklabel.set_backgroundcolor(bgcolor)\n",
    "        plt.tick_params(left=False)\n",
    "            \n",
    "    \n",
    "    #plt.xlabel('Testing time point')\n",
    "    #plt.ylabel('Training time point')\n",
    "    #\n",
    "    print(kk)\n",
    "    return()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "981f9f83-d720-41ba-998e-1f65da8adf98",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "cmap = mpl.cm.get_cmap('Blues')\n",
    "\n",
    "cmapS=cmap(np.arange(0.35,1,0.2))\n",
    "\n",
    "cmap = mpl.cm.get_cmap('Reds')\n",
    "\n",
    "cmapD=cmap(np.arange(0.35,1,0.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5d8b32fc-6cd2-4fcc-a814-72c80cb1e2cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "def print_sig_ax(ax,Ui,cmapX=cmapS,stim='Sample'):\n",
    "\n",
    "        \n",
    "    n_timepoints=45\n",
    "    #timepoints=np.linspace(0,5.5,n_timepoints)\n",
    "    timepoints=np.linspace(0.,4.5,n_timepoints)\n",
    "    sig=Ui.reshape(4,n_timepoints)\n",
    "    lines=[]\n",
    "    for j in range(4):\n",
    "        lines.append(ax.plot(timepoints,sig[j,:],color=cmapX[j],label='%d'%(j+1))[0])\n",
    "    #if i==0:\n",
    "    #    ax.legend(lines,['%d'%(j+1) for j in range(4)],title=stim+' no.',loc='upper left',bbox_to_anchor=[0.7,0.5])\n",
    "    ylims=ax.get_ylim()\n",
    "    ax.set_ylim(ylims[0]-0.8*(ylims[1]-ylims[0]),ylims[1])\n",
    "    if i==0:\n",
    "        tmpposes=np.linspace(*ylims,4)\n",
    "        for j in range(4):\n",
    "            ax.text(4.55,tmpposes[3-j],'%d'%(j+1),color=cmapX[j])\n",
    "    \n",
    "    tmpax=ax.twinx()\n",
    "    tmpax.fill_between(timepoints,sig.std(0),color=cmapX[-1],zorder=-2)\n",
    "    ylims=tmpax.get_ylim()\n",
    "    tmpax.set_ylim(ylims[0],ylims[1]+1.25*(ylims[1]-ylims[0]))\n",
    "    tmpax.set_yticks([])\n",
    "    if i==1:\n",
    "        tmpax.set_ylabel(stim+' info                   ',color=cmapX[-1])\n",
    "    plt.sca(ax)\n",
    "    #DrawLines(axis=['x'],linetype='--k',lw=0)\n",
    "    for side in ['top','right','left']:\n",
    "        ax.spines[side].set_visible(False)\n",
    "        tmpax.spines[side].set_visible(False)\n",
    "    ax.set_yticks([],[])\n",
    "    tmpax.set_ylabel('s.d.',color='k')#color=cmapX[-1])\n",
    "    tmpax.yaxis.set_label_coords(-0.1,0.2)\n",
    "    #ax[i,0].set_ylabel(DPFC_comp_names[i],fontsize=24)\n",
    "    ax.set_zorder(1)  # default zorder is 0 for ax1 and ax2\n",
    "    ax.patch.set_visible(False)\n",
    "    DrawLines(axis=['x'],lw=0)\n",
    "    ax.tick_params(axis='x', pad=-2)\n",
    "    for side in ['top','right','left']:\n",
    "        ax.spines[side].set_visible(False)\n",
    "        tmpax.spines[side].set_visible(False)\n",
    "    #ax[i,0].tick_params(axis='x', which='major', pad=15)\n",
    "    ax.tick_params(width=10)\n",
    "    ax.set_yticks([],[])\n",
    "    ax.set_ylabel('Activity')\n",
    "    ax.yaxis.set_label_coords(-0.,0.75)\n",
    "    ax.set_title('SC%d  '%(i+1))#+PFC_comp_names_new[i])#,fontsize=24)\n",
    "    ax.set_zorder(1)  # default zorder is 0 for ax1 and ax2\n",
    "    ax.patch.set_visible(False)\n",
    "#ax[0,0].set_title(r'data SC response',color='C9')#$\\vec{u}_i$')\n",
    "    ax.set_zorder(1)  # default zorder is 0 for ax1 and ax2\n",
    "    ax.patch.set_visible(False)\n",
    "    DrawLines(axis=['x'],lw=0)\n",
    "    return(f,ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "harmful-professor",
   "metadata": {},
   "source": [
    "## demixing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "equivalent-hebrew",
   "metadata": {},
   "source": [
    "the data should be organized into the shape of (profile,neuron) to feed into SCA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recreational-bulgarian",
   "metadata": {},
   "source": [
    "below is an example function of demixing data of shape (neuron,Stim1,Stim2,Time) and transform it to shape (profile,neuron)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "streaming-negative",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_XnXT(MF,stim_shape=(4,4),stim_keep=0):\n",
    "    #MF, the original data, takes shape (neuron,Stim1,Stim2,...,Stim_n,Time) or (n_neuron,n_profile), n_profile=Stim1*Stim2*...*Time\n",
    "    #\n",
    "    #stim_shape should be defined according to how many values there is for each stimulus\n",
    "    #\n",
    "    #stim_keep means which stimulus dimension to keep\n",
    "    #The result is the neural data contributed by the kept stimulus and its interaction with Time\n",
    "    tensor=MF.reshape((MF.shape[0],*stim_shape,-1)).copy()\n",
    "    tensor-=tensor.reshape((MF.shape[0],-1)).mean(1).reshape((MF.shape[0],*([1]*(len(stim_shape)+1))))\n",
    "    XT=tensor.copy();\n",
    "    for i in range(len(stim_shape),0,-1):\n",
    "        if i-1==stim_keep: continue;\n",
    "        XT=XT.mean(i);\n",
    "    T=XT.mean(1)\n",
    "    return((XT-T[:,None,:]).reshape((MF.shape[0],-1)).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "molecular-senior",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_X1X2T(MF,stim_shape=(4,4)):\n",
    "    # this is an example of getting the interaction between two stimuli and its interaction with time\n",
    "    # this currently is not generalizable to more than 2 stimulus dimensions\n",
    "    MF-=MF.reshape((MF.shape[0],-1)).mean(1)[:,None,None]\n",
    "    nneuron=MF.shape[0]\n",
    "    Xtensor=MF.reshape((MF.shape[0],*stim_shape,-1))\n",
    "    \n",
    "    Xt=Xtensor.mean(1,keepdims=True).mean(2,keepdims=True)\n",
    "    Xs=Xtensor.mean(2,keepdims=True).mean(3,keepdims=True)\n",
    "    Xd=Xtensor.mean(1,keepdims=True).mean(3,keepdims=True)\n",
    "    Xst=Xtensor.mean(2,keepdims=True)-Xs-Xt\n",
    "    Xdt=Xtensor.mean(1,keepdims=True)-Xd-Xt\n",
    "    \n",
    "    return((Xtensor-Xt-Xs-Xd-Xst-Xdt).reshape((nneuron,-1)).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "jewish-christmas",
   "metadata": {},
   "source": [
    "import demo data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "academic-herald",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('demoOriginalData.pkl','rb') as f:\n",
    "    MF=dill.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "located-difficulty",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(467, 16, 45)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MF.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "optional-deadline",
   "metadata": {},
   "source": [
    "it has the shape of (neuron,stim1*stim2,time), the number of dimensions doesn't matter, i.e. you can either merge adjacent dimensions or not, but the order has to be like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "powerful-newport",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(180, 467)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpX=get_XnXT(MF,stim_shape=(4,4),stim_keep=0)\n",
    "tmpX.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polar-representative",
   "metadata": {},
   "source": [
    "here we keep the first stimulus dimension, each profile has a size of 4*45=180"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bridal-clinton",
   "metadata": {},
   "source": [
    "## storing the data to be fed into SCA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inclusive-samuel",
   "metadata": {},
   "source": [
    "make a folder to run SCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "informational-governor",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "folder='demo/'\n",
    "# this should be the same as what you use in the SCA script\n",
    "if not os.path.exists(folder):\n",
    "    os.mkdir(folder)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "amazing-terrorist",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(folder+'tmpX.pkl','wb') as f:\n",
    "    dill.dump(tmpX,f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "certified-albuquerque",
   "metadata": {},
   "source": [
    "## run SCA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outstanding-electricity",
   "metadata": {},
   "source": [
    "if it crashes, run it in raw python or ipython console"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "built-setup",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dSCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "internal-telling",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function dSCA.Learn(X, X2=None, test=False, n_components=3, alpha=0.5, n_iter=600, lossname='lossEla')>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dSCA.Learn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conceptual-mailman",
   "metadata": {},
   "source": [
    "in this demo we learn the decomposition $$X=UV$$ \n",
    "by minimizing \n",
    "$$Loss=||X-U V^T||_{frobenius} +\\alpha\\sum_i^n||v_i||$$\n",
    "where $$ V=(v_1,v_2,...,v_n)$$ $$U=(u_1,u_2,...,u_n)$$\n",
    "under the constraint $$||u_i||_2^2=1 \\,\\forall u_i \\,in \\,U$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ranking-wireless",
   "metadata": {},
   "source": [
    "by default, a small amount (0.01) of L2 regularization is also added to the loss function, this makes the optimization more stable. Alternatively it could be set as a fixed fraction of $\\alpha$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "latter-railway",
   "metadata": {},
   "source": [
    "in the equation, there are 2 important hyperparameters, n or n_components and alpha, the strength of L1 regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "iraqi-contribution",
   "metadata": {},
   "outputs": [],
   "source": [
    "U,V,losses,errs=dSCA.Learn(tmpX,alpha=0.4)\n",
    "#this uses the default (arbitrary) hyperparameters, \n",
    "#you should choose your hyperparameter for a good reason, or do a grid search for it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "available-chemical",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(180, 3)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olympic-comment",
   "metadata": {},
   "source": [
    "this is the shape of our profile matrix (size_of_each_profile, number_of_components)\n",
    "\n",
    "to plot it, take each individual profile vector and reshape it into (stimulus,time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "marine-distributor",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n",
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n",
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n",
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n",
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n",
      "[0.51809945 1.03684093 2.05390537 2.57065977 3.59486604 4.5       ]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAAHdCAYAAAAHCXsOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7lklEQVR4nO2dd3wUZf7H3zNbk2x6bySUUAOBRLGigIiKoCJYEBE9ORuinj/OwqnHqXd6HtgLdyqCiKKiCAhK7whISYCQhISE9N6T7Tvz+yMSCQnJbrIJiez79eL1IrMzz3xm9rNPf76PIMuyjAsXXYh4oQW4uPhwmc5Fl+MynYsux2U6F12Oy3QuuhyX6Vx0OS7TuehyXKZz0eW4TOeiy3GZ7jeWLFlCQkICnp6e+Pr6MmLECJ5++ulm5xUXF/PUU0/Rt29fNBoNvr6+3HTTTWzYsKHxnIyMDB5++GHi4uJQKBSMHj26C5+k++MyHfDaa68xa9YsbrjhBr7//ns+//xzbr31VtasWdPkvLS0NEaMGMH69euZO3cuGzdu5PPPPyc6OppbbrmFpKQkAJKTk1m/fj39+/enf//+F+KRujWCa+wVwsPDue222/jggw+aHJdlGUEQGv++5JJLMJlM7NmzBy8vrybnHj16FB8fH3r16oUkSYhiw+956tSplJWVsX379k5/jp6C8kIL6A5UVVUREhLS7PjZhtu5cyeHDh1izZo1zQwHMGzYsMb/nzGci5ZxvR0gPj6e9957j6VLl1JeXt7iOTt27EChUDBu3LguVvfHw2U64IMPPkCn03H//fcTGBjIkCFDeOmll6ipqWk8Jz8/n8DAQNzc3C6g0j8GLtPRUDSmpKSwZs0aHnvsMWRZ5pVXXuGSSy6hrq6u8byzi1sX7cdlut/QaDRMmjSJ999/nxMnTvDJJ5+Qnp7Op59+CjQ0NkpLSzEajRdYac/HZbrz8OCDD+Ln50dqaioAo0ePxmq1smXLlgusrOfjMh1QUlLS7FhpaSnV1dUEBwcDMGrUKBISEpg3bx61tbXNzj927Bi5ubmdrvWPgKvLBBg6dCi33nor48ePJygoiOzsbBYsWIC7uzszZ85sPG/58uWMGTOGSy+9lKeeeorBgwdTU1PDhg0b+Pjjj9m/fz+RkZHo9XrWr18PNDRAampqWLlyJQATJkzA3d39gjxnt0F2Ib///vvy9ddfL4eGhsoajUaOioqSp02bJqekpDQ7t7CwUH7iiSfk3r17y2q1Wvbx8ZHHjx8vf/fdd43nZGVlyUCL/7KysrrwyTqHtJJa+eeUYnnFkTx5x6kyh693jUi4cJjcKgMCUFRrQm+xcU0ff4eudxWvLhwm0qehr7LSYEFvsTl8vash4aLLcZnORZfjMp2LLsdlOhddjqsh4cJhJFlGlkH+rSPIJjV0gChE+8amXV0mLhzmWGENx4uajsoE6dRcFxNo1/U9wnQZGRkOX9OvX79OUNIUR3V1hSZo3/uyB2fp7/F1uvIaAz3gd+PiLHqs6SxWiQUrDjL6yW959r+7qNWbL7QkF3bSI01XUqnnwTc2sHxzCneNHcCuo/nc8fcfOZZZdqGlubCDHmm6fyz9hbzSOj579gZemHEZ3/5jIn5eWh5esInMwuoLLc9FG/Q405VVG9h9tIBHboljeEwQABGBnvxv7jiCfN2Z8/ZWqutMF1ili9bocaZbvy8LhULghpFRTY7r3NS89+RYavRm/u+jHVitjg9Eu+gaHDbdwYMHO0OH3azde4rRwyPw9tA0+ywyyJM3Z1/LwdRi3vl67wVQ58IeHDbd+PHjGTp0KG+++SalpaWdoem8nMytJDWnkklX9j3vOZcODGHauIH8a8l2souquk6cC7tx2HSFhYXMmzePjRs3EhUVxeTJk1mzZg02W+cXZz/+komPTsPVQ8NaPe/x24bjrdMy9531na7JheM4bDqNRsO0adP4+eefycjI4LLLLmPu3LmEh4czd+7cTusNt0kSP/6SyU2XRaNSKlo918NNxX+euIkfd6eybk9qp+hx0X461JDIyMggLS2N4uJihgwZQmVlJQkJCbz22mvO0tdIanYFpVUGbhwZbdf5t48ewvjL+vF/76zHaLI4XY+L9uPwLJPs7GyWLl3K0qVLsVqtzJw5k8TERHr37g3AvHnziI+P5/nnn3eq0OTTFShEgcHR/pTUWtiaXkNqsZHCagszR/pzZR/PJucLgsCCJyYQf9/7vPP1Xp6971qn6rmYsdgklv6aT3JhHbUmK77uKiYODuTafvatlXDYdAMHDmTSpEl88MEH3HDDDc1CLfTt25dp06Y5mmybnMgup2+4D6Io8uya05TVWRkY7IaXVsH8nwp45Oogbo/zbXJNTK8AZt9xOW8s28m9N40gPLB5tCUXjmOTZXzcVDw3rg9BOjWnyvS8sS0LP3c1Q8M827zeYdOlpaXRq1evZsfz8vKIiIgAYNGiRY4m2yYnTpczKMqPLw9VUFRjYdFd0UT5aZBkmU9/KeWj3SUUVpt55OqgJvO6np85mq82JPHCRxv57KWpTtd1MaJVKpga93totX6BHgwO1pFWWm+X6Ryu08XGxrZ4/Oz4bM7GYrWRnldFcIA3Kw6Xc3e8P1F+Df10oiDw5yuDeOLaYNYcr2L+T/kYzBI2SeaX9DKyyo3846FxrNh0lK0HT3Waxp6MTZL4ZlsaLy3ey8ncSoevN9skTpXr6eWjtet8h3O6lqYRGY3GTo1olJ5XhdUmcahUIsxLzbQEv2bnTIr1IdhTyasbCnjs29MYrNmU15lRK0Q+/XMC113al/vmf8ueTx4hKsSn07ReKCRZps5kRelgQMa0nAr+sXQfxzLLCPRx44fdGUy6si//d2cCfl5tm0iWZT75JZcQTzWX9PK26552m27EiBEIgoDRaCQ+Pr7JZ4WFhVxzzTX2JuUwJ7LLEQTI0yv5z41BqJUtv9iRUTreur0Xn/xSSlyfYCbEhfL2zyd55LPDfPTgeP78j6+Y9rev2PLhLNw0qk7T62xkWabWZKW4zkx5vRlfNxVRvm6oFCI5lQZOlddTobdgk2UUAvR3E9Ap2p5jmFVYzYx//kx4oI7P591IbO8AFq46xoqNyexIKeOxu69i3EBvgjxbfleyLPPZgXwKa0w8P64vop0Zj90zh5cuXYosyzz66KNN6myiKBIcHMzYsWNRKjtnycXMl75gw+F8Bo6M5+Np0XblqmdmuepNVmYs2k9WaT1/n9CH+1/8kntuiOODZ27tsK6umjl8tLCG5KJaBMBbq6TGZEWWQaUQMNtkQjw1hHpp8dYqOVZYQ43BzGB3I+pWMj2L1ca9//wZvdHC1/NvRhAVLP6llB+OVRHtZmLfrkR8IiPxi4zgk2nRBOpUTfTLssySX/PJKNUzb1wfPDT2f/d2n3kmkMzw4cMZPny43TdwBscyy7Go3LhhkLfDxbi7RsniP49k2gf7eOWnLJ6eMZp/fbqZP91yCQkDwztJsfOoMVo4UVTLwCAdsSGeqBQiJquN7EoDBouN3n7ueGl/z4l83FSsSy4kw6BhoLuJ862V+fCHJE7mVvDR3PGsTKpm9bEqDBaJ2aOCuHWoD+9521i8/jg+/j68v7OYZ8YEozeacdeqAVj6az7pJfU8f31fhwwHduZ069evZ8KECQDNwtyfzS233OLQze3BYrXhN+4V3MMiWfvsVfh52PeA5+YqFXVm7v1oP8VVBoyZKYQHeLLh3QfaVRfNrTSw+ngJ9fW1eKtFor1E+vm0rcvRnE6WZbadKqfeZGXCoGC7V1slpZ0iRa8hUmMhWP378GReaS1H0ks5mF7K99vTiOwXhdUnBI1S4MbB3kyJ8yPEq8HAFqvE/a//zPHMcqTfLPKvx8bzl2lXU1Zn5qkfUlCJAuJZmq7q7cufLotoU59dpouNjeX48eMAjZ3AzRISBDIzM9u8oaMczSjisgc+ZNToeD6c2XLLuSVa+oLPGO9kZj75x0/wzb+mMWnUILvTNFklVh0r4qcTpQTq1KixUmWSqTHLXBep4uowVasmdtR0uVUGdmdVcE0fP8K97Y91nJGRQUq9GrUo09etYTQmLaeCe15Zj9kqodJo0Pn5cNt1wxgc6s6VfXR4aZsPLZbXGNhw4DTrTtRQaZL59pnriesb7NAztIRd2cYZwwFkZWV1+KaO8NP+BiPfflnHi0I/nZqvZl/OY0sOsS4nnzlvrueGy2NQq9p+DRV6Cwu3ZVJQbeL2YSHcPDiQ01mZyLLMjnwLW3ItVJlkrghV4a0RULWSK5XUmvD3ULeac1XqzRzKqyLMS+uQ4c7goZCotDYYSW+y8NdFu+gV4kXgoMHUmAXem9qLQF3rjSl/LzfuGTeIMZeamfXVabanVzvFdA73073xxhsUFBS062ZWm8SSnVnUGOwbC5UkmW92pKF2c+PaGJ923fNcvN1VLH14JNMmXUZxWTUD7v+Ud34+SV6F/rzX5FQamP9zOrUmG/+4KYZbhwajVDS8OkEQGB2h5tY+ao6UWnk/ycA/D+hZnGzAZGteiGxLL+f/Vqfy3q5srDap2eeyLJNaUsvGk6VolAouibSvG+JcdAoJsyxiluCNrw5SWF5H9JCBFNZKvHpzeDPD2WSZonobKRVWDhRZqDL9ri3UW82rE8P585g+7dJyLg6ve73pppvYtm0bo0aN4r777mPKlCl2R5ZMzK7irvd+QaMSuefKXjx4bW8Cz9MXZLLamP3pAb78cgPXXtKH9x690hGZdhVlL3y6g4VLtuAbGYFXRCRP3dCfWaN7NxoKYF92FZ/syyVYp2HumN74uv/+ZZ3beq01S5QZZCpNEj9nm4n2UnB3f01jV8JpyYclB/IZEe7FscJahoV5MmdUFApRQG+2kVNlILNcT63JyoBAD+LCvO2ux51NRkYGZgmS6t0oSk1n3odbGTlyEPl48Y+bwrksWtd4bmG9jcRSK8fLrdSflRe4K+HuAVp6ef5e7Dpr3Wu7FluXlJSwfPlyli1bRnp6OpMnT2bmzJlcd911bV9bbWTxziy+3JuDQhT4z7Q4xsUGYzRZWLn1OH5e7gwbFMn/fZnE1h1HkKoq2PCf2/DR2dfbfQZ7X9CC5bt4cdEmxl47nBS9hhAfN0b09mVgqCe1VonkUj1XRPtybbQ3/92aSZ3Ryn+mDSPQS9tql0l6pZUv00xcGqykt7eClAorR8ts3DgwgMlDg9mZWcGKw4XotEr83FV4aRWEeKrpG+BBTIAHgbrmM6PP5mRpPSdL6hnZy5sgz6bnntGVVKfh5YU/UF5nQR3dn+fHhzEmpmH8ucwgsSXXTEqFDQ+VwFB/BYP8lPhrRQQBvj5pJL9O4ubeaoYFKFGKwoU13dkkJyc3bs7W1kTO4so6lm88RqCXGwqFyJKdWSSdrqCPr5KM9Fzq6hvC5fv17Ye3vw+5hw/zt/tHM/WqtltE52LvC6qoNzP9X6vZvvMYETFReEf2otZgxWxpKF48NEqGhHvxa1YFoT5uWKwSClHgowcS0FlaXvJYqzdzIKWQ74+UcjS7Cv9gf665rC+X9wsi1EtNmd6CQEMMkIwyA2V1ZioNVhSiQHy4F9cPCGBwiI6KGj2+nm5NGieSLLPuRCnfJhb+9jcMDtExNNST3n5uBHtqyMw6jcEqsyKxnKWLf8anTz+GJcQQ6OeOWgSlCAYreKkFxkSqGBqgRHFOA8gqyazNNJNUZsVNCYP8lExO6E1v/47HS2636QoLCxtzu9OnTzNlyhQWL17c6jU/HzrN5KeanyMqlWj9/HEPCsZQXIi+rJTgED8MdQbSv/0LxYV5DutrzXSSLFNSayaxoIbvkopQKUQ8qsr435fbGdovhEFRgZwqrMJikwkM9gc3HVOujuGOyyKprDfz2JJDHM+tIS7CjchAHQqFQEZBLckFenSSgZwTadTVGxFEEbVGjclgRFQoiB3ci1uuGcTNV/SnsKCM1TtOkHK6BG8PLR7uGgySQIVZoqzGgKG0gqLSal740xj+9sAYAOpNVhbtzeVIfg23xAYxcXAQh/Oq2Z1ZSUaZHqO1aR3x8I5DZJ/K56EZY5kYF0C5UcJkazCUu1JgWKCy1QYPQFG9xPHyhuJ3clwY1/UPcPi7OBeHTffFF1/w+eefs3PnTsaMGcOMGTOYPHmyXdsXGS02TpXpKajSczSvipJaM96eWjzUSvZnVWK2SOjclPyy9SA56TlcPnoEHz4+DlVtscMP1pLpcioNrE0u4Uh+DcbfcrIx/fy4a0QoOo2S/cm5PLHwR2yAl5c79XoTqRmFmC02Xnh4PJOvj8Nik6g3Wll1MJ/E05XklOqxSTK9gz2Qyov55Zdk3L109B0ygIfG9iIizJedqeV8s/4IRfklWA2GRj19I/y4IrYXdQYzlbUGKmoMVNYYsMjgHuCHv7uKQ4dP8uObM+ndO5R3d2ZjsNh45KpeDA9vOk1LkmWKakyU11soLirAZDRz17zv8IuMYOVz1+Kn6djYuCzL9O7bD2U76pjn4rDphg8fzowZM5g+fXqLOwe2RUG1kYN5VRgsNoaEeDI42BNREKg1WtmSXo7BYiMh3It9R0+zudhKoKeGPw1UolE49rBnm67GaOWTfbkczqsh0EPNtf386OvvTpSfG17a37tLUkvqOJJfjU6jIMBdjVIhUlil58MvtvPLrxm8+vxUosL8UCtENEoRrbkWH4UVqySzeGMai1b+yvjRQ5k5eSRapUC1TYFSkInUWPAWrLz7Sw0bD+Vw93Bfrh0eRXatxIHMSjKK68gt1zeG3IKGIS61WkHh8WTMej0jrruSK4aE8texfQnQqVt99oyMDP7+ZSLfbzrK44/dym1xfk06idtLt6nTOUJZvYlNJ8sI9tRwaYQPntrW+8cO5Vbzzs7TRHgquae/Gq3SfuOdeUHl9Wb+vSWTOrONe+JDuTzat8Vfa16VgV1ZFQwK0jE8vGk3RZ3BxJUPLsLTQ8P2j/7cuEbjTIX9p/1ZPLNoF3+aEMv9t11CqUWJVRYIUlvxU9oah6J6RffhsSWH2JFSilWSUStERvb1o3+IJ32DPQjw1ODtrsJqkzlZVEtaYS0pueVs/Wk3olpLYEw/7riqL1MujSA+2rdxNMAmyYhCQ/dNZb2ZV77axydf7yY4wJOFz92GRpQaO4k7QpeabtGiRTzyyCMAvPvuu+c974knnmg1HVmWKa4zE6xT2z389PWRAn5OKcNDBVP7aQh0FzFaZdyUAupWcr++ffuSUabn/V3ZCAI8d11fQrxabhGW15vZklFGmJeGq6L9WtR2MCWfMY9+zLXxvXHTqCiuqEMl2PD3duPnA6e54dJo/vXnq9ockTBZbXywKYOoAA+ujw3Gy63t2S57jmYz4+/fUFpZj3+vCCxqHQE+bsT28qewzkpOuR6VQsBfpyEvr5iyjAxUIiz6v3EERoZRYRWJ8zDR0dlnXWq6CRMmNO4AM2bMmJYTEgS2bt3qFFFnk1pSxy9Z5aQW1JBf/3tF2U0JN0SpiQtQNvmizTaZbXlmTlYLlOstRHhreea6Pvi5t/zl5lYZ+CW7El83FWP6BbRaZ/lk9a8sXnOQYD8dQX46CorLKanUEx3ixfwHrmhzlVpHvrR6g5n/fLGTt77ag/msMPqeOjeiwv2xWGyUV9VTVlHDyMGhvP7nqwj0cafOJpKi1xDjZsJH2bwz2hF6ZPHaHvKrDezMrCDWTU9GVcPL1irgaJmVY+U2enuJTOitIdBNxGSTWZ5qpKheYnRMAPERXgwM1rVoJFmWOV5Uy/GiWnr5uHFZlI/DEyAvRFDE6jojReW1VNYayCupIelkIcmZxbhpVYQHejFiQBgJ0e5NBuJP1GsQBZmB7h0Lp+Ys0zk8AW7o0KEcO3as2fHhw4eTmJjoDE1N8Pxt2owFBUMDfn+RMb5K4gKtrMsy82GSgZEhSvLrJEoNEjMGaRkzovWx2tJ6M8eLahka4smQEM8es5ert06L91kd5VPHNp8Ece6PIURt4ZRRQ71NwMOOyZ328PfvjrPpeDG1BiseWgUT4kJ5btKg806wPRuHTXf69OkWj+fk5DialF14qJWAjElqbop+PkpmxynYV2hhZ74FUYD7BmkJ17VezAFklNXjqVH0KMO1F1+lhFqQKDYr6eOEBgXAjKujeHbiQNw1SsrrTDy+9DD/3XqKOeNj2rzWbtM9/fTTAFgslsb/nyEzM5M+fZwzGHwuClFAI8gYWzAdgFIUuDpczfAgFZIs49XadNnfMFps5FYZiAvz+sMbDkAQIFhtJc+kIkKytDqj2F76BTdd9SUKAqdL6+261m7TVVY2rBKSJKnx/9AwXX3w4MG8+eab9iblMBpRxii3/qZ0KgGwz0CnyvUIQG8/j46L6yEEqmwUmFTkmNT01Zo73JIF+GhzBh9szqDeZMPXQ8WzEwfadZ3dpvvss88AiI+PZ86cOe1T2U60okytzTmh9CRZ5lR5Pb183dHYUf/4o6AQoLfWTIZRQ7FFQYiDncUlZgVBRkuTqfGPjuvHo+P6kVFcyw+HCgg8T5fUuTj81gcMGEBaWlqTY6mpqWzatMnRpOxGI0oYJQFntLMLa0zUm230C7h4crkz+KokQlQW8kwqaq32f/UGm0C2SUVRbcsRTvsFezIozIu/fplkV3oOm27OnDl4eDT9wjw8PDo199OKMjICZrnj44cpxbX4uanwP0+/3R+dcI0VD4VEplGFZOePON+sQi3I9PU//w/VapPIKjv/RNizaVd8ujPhI84QGRlJfn6+o0nZjVZseDsttWAdIa/aSGm9mWEXSQOiJUQBojUWzLJIubXtVn69TaDSqiBcY22cUFpvsvLt/lxqDJaGmc4FNby/KYNrBtg3A8XhLpPw8HASExObLENMSkoiLKz1QIUdQS00bEJllATaGwLHJskkFlQT+tsa0YsZN4WMj9JGkVlJgNLWaqMiz6RCK0r4K3+vAwrA6sMF/GtNCmarhL9OzY1xofzlxv523d9h0z366KNMnTqVf/zjH/Tr14+MjAzmz5/PY4895mhSdiMKoBFa7quzl/SyeupNNq7p7djW339UQlRWUg0aqm3ieYfHaq0iNTYFfbVNx23dNUq+ePSydt/bYdOdqbv985//JDs7m+joaB5//HHi4uLaLcIe7Ok2OR9mq0RyUQ19/N3xtmOA/WLAUynhIdooNCvxUbY8PFZsUeImSvh2cMz2XBz+FgVB4IknnuDEiROkp6czffp03nvvPW6++WanCjsXrXj+DuK2OFFci02GoaGu+HRnE6q2UmdTtNgdZZGgyioSoLI6pU/vbBw2nc1m44cffmDSpElER0fz4osvMn369E4bBjuDRpQwtaPbRG+2cbK0joGBOtxUbVecLyZ8lBJuokR2Cy3ZcmtDIRigcn4Ac7tNl5aWxjPPPENYWBgzZszA19eX9evXExgYyOzZs/H379y6Unu7TY4X1aAQRQYG69o++SJDEKCP1oxREsg3/V7TkmUotSjwVdpwYN6s3dhdpxs0aBD+/v78+9//5q677mrsq+uqrocz3SYGSUAj2pfd1RgtZJbrGR7uhVpx8Yw+OIK7QiZcYyXPpMRbKeGllKiTRIySSC9N5wQIt/ubuP/++zEajfztb39j/vz5TUJNdAUaQUYpyNTbORxmtknsy67ETaUgJsCVy7VGiMqKp0Iiw6Am06CiwKRELUh4KZzbgDiD3aZbvHgxRUVFvPLKK+zZs4e4uDji4+Opra2lurrzdx4UhIZQCfaMwVpl2J5RRo3JytW9/dq1Sv5iQhCgr9ZMoNqKXmroJglUtd5/1xEcKnM8PDyYNWsWe/fu5dixY4wZMwYPDw/i4uKYNWtW5yg8C50oUW8TW21MWGVI02uoNVkZ2y8Af4/WV065aEAlQqTGSqyHieEeBkLV1k67V7srOoMHD2bhwoXk5+ezbNkyCgsLnamrRXQKCQkB/Xm6TiQZMgxqTJLA2H6B+Lm7DNceVCKdlsuBE7beVCqVTJkyhXXr1jlDT6t4KCQEZOpaKGJlGU4Z1dTbRGLczU0C3bjoXvSoJp0ogLsotWi6PLOSKqtIXzcznp1UAXbhHHqU6aChiD3XdFYZSsxKQtXWDi+zc+EYBpOFIXe/TchN/7L7mh5nOs/fgv2dPfhfYlYiA8Gqzqv8umiZlz/Z6vD2Vz3OdLrfis4zuZ0kNwxMB6psqHrc0/RsjqQVsGHfSf567yiHrutxX5NKBI0gUfNb10mZRYFVbljt5KLrsFptzH5jNW8/PRGN2rHJSj3OdNAwUF1mUXKsXkO+WYWv0tY4TOaia3j7673E9g3hmhEtR9tvjc7Z4qaTidRY8FHaqLAqqLGKhLlyuS4lM7+C/35/gH2LH23X9T3SdIIAXr8NTrvoevYkZVNWVU/Cfe8DYLbYqKk3EX3rG3z72j1cOrj1cL090nQuLixTr4vl+st+D6az73guD7+2in2LH8Xfu+2YxC7TuXAYN42qyS6Sfl4NwbhD/NveYBh6aEPCRffimhG9Kfppnt3n93jTfb8tmR92JF9oGS4coMcWr1abxPvf/MKWg6dQiCLXjuiDr5fje2i56Hp6ZE5ns0n845MtbD14CmgIF7F+T1obV7noLvRI06XnlZN4spAz3cGSLPPjnlTMFld/XU+gR5rudEFls0h0dQYz2w917bagLtpHzzRdYSWKc1Z3CcCq7a4GRU+gR5ruVF55s71SZSCvtIZDqZ0XPcqFc+iRpssqqGzxuCgIfLnBvsB8Li4cPc50xZX1mCwthzqQZJmTOWUcTS/qYlUuHKHHme50QUWrnzfkdoldI8ZFu+iBpqtEbGVnG0mWSc4qYc/R7C5U5cIRepzpsgoqaSt0kygIzP/fZrr5DlQXLT3OdKfyypHaMJMky+xOymbl1q6Nt+LCPnqU6UxmC8UVdXadKwjwf++sp1bfchh6FxeOHmW67OIa7C0wZRnKq/W8+qnztwN10TF6lOlO57fecj0XSZJ5/9t9nMwp6yRFLtpDjzJdVmElSgeDGwoCri6UbkaPMl1haU2z4a+2sEkyKzYddbVkuxFOn8SZl5fXbEcdZ1FSaV8j4lyyC6tISi9ieP9QJyv6Y5FXUsOepNPsO56LJEm89ZeJTXbIdhZON93gwYOpqalxdrIAlFXbt/fUuShEge+2HneZrhVKKut5YuFabDapsUvqRFYJsX2DnX4vpxevycmdM73IaLZgMLVvkqZNkvl6s6uIbY3Faw9ik+RGwylEkV2JnTM/0emmi4yMdHaSAJRXGTp0fW5xNUdOdn600J5IcmYxe5KykaTf68s2SWJXYjaSvVslOoBDxevnn39u13n33Xdfu8S0Rnk7i9YzKBQi3209TvyAzts4ryciSTKLvj+AKAjNRnpq9SaSs0oY6uQi1iHTvfXWW03+Tk5Oxt3dndDQUAoLC9Hr9cTGxnaK6dpbnzuDzSax/OdEXn5oXLNZxxcz2w9ncbqw5fmJZ4pYZ5vOobd/5MiRxn8333wzzzzzDGVlZaSkpFBaWsqzzz7baXuElVXVoWhldok9FFfUselAhpMU/TFYtycVsdmKkwZsksTuTihi2/0t/ve//2X+/PkolQ2ZpUql4qWXXmLRokVOE3c25dX6Dkf8Vogin6456BxBfwDySmo4mVOG1MrgYq3exPFTxU69b7tN5+bmxuHDh5scS0xMRKvtnA18S6vqHe4YPhebJPHT3jSKymudpKpns+XXDMQ2fskKUWTjgZNOvW+7++meeeYZxo8fz7333ktUVBTZ2dksX76cV155xZn6GimuqHdKOrIMX/ycyNzpjoUs/aMhSTKbDmS0OU3MJknsOpLNA5P09Gv1TPtpd073+OOP8/3332M0Gtm2bRtGo5GVK1fy+OOPO0laUzraej2DJMt88sOvF32f3ZGT+VTXGe0+/+c9zsvtOjQiMXbsWMaOHessLefFbLFSb2h59+X2kF1UxZ6kbK4eHu20NHsam/Y1xICxSW1XWc5EUHj2IQtaTcc3hXF630FeXp6zk6SipmMdw+eiVIis3pni1DR7ErnF1fxyPMcuw52hzmDm683HnHJ/p5tu8ODBzk6SsirnFK1nsNokftiRfNEWse9/s9fhzb8E4N2v9zrl/j1i7LWsyjmNiLPJK6khObPE6el2dzYfyODE6dImQ172INOwPsUZtNt0ixcvbtxoOCkpiWHDhpGQkNApe7+WVetbXXbYHhSiwI+7/zhFrMFkIfV0KZsPZGAyt7wjdU29kU+6QT9luxsSr776Kvv37wfg2Wef5cYbb0Sn0zFnzhy2bdvmNIEA5VV6RAGcGUvdJsms2n6C52aOdmKqnU9KVgkHU/KJCPIiyE/H9sNZfLUhkUOpBY3nzLqhL0/efVWza1fvSMFg7Jwt0h2h3aYrKysjMDAQo9HIL7/8wurVq1Eqlbz99ttOlNdAaXUdtg52DLfE0Ywi8ktrHN7b6kJwNKOIVxdvZe2u1CbHW6qabf71FNfG92k2f3DrwVNt9st1Be0us3x9fUlPT+enn34iISEBjUaDxWJxuK5gDyXl9XavAnMEAVi3O7XN8y402YWVXP3n/7J+b/O+MlluvvZcFATeWbEH41nFbGp2aYcnTTiLdpvuySefZPjw4dx777089thjAOzevZtBgwY5TdwZOutlCYLAFz8ndkou6kxeXbwNSZbt1inJMhXVBr74KbHx2M7DWR2eMOEs2l28Pv3000yaNAmFQkGfPn0A6NWrFx9//LHTxEFDx3BnLZiWZJlfT+TxyuKtzP/zOLuuyS6qYuWWY3y9+Rjm2lJ0bmr8vN15ZMrl9I/0d7rGtOxSlm9IbCuSRjMkZNbsTOGSQREM6xfCjsNZDvXLdSYdGpGIiYlp8nf//v07JKYlThdWOT3Nc/n35zuJHxDOLde0nksfzSji6j//F5skIUky0X5gshiorDHy6qdbeHfuJHx0zo3w/vInWxFFsV25sYDAf5bt5KHbL6WmG0U66B75bSv89/v9bc6E6CgC8MArK0nLLj3vOVarjVmvfof0m+HORkKmus7Ev5futMsc5dV6jma0HUPvSFoB329PbnfxLyFTZzDz1pd7u03RChfAdAaTBZPZvgU2ReW1fLL6YKe3uGTAZLYy9bkvz1uUv7ViD8dPFWM7z4RGSZZJzizmX0u3s2TdIb7ZfLTFAXWbTeL2Z77g8j99xIcr951X0y/Hcpj49FIUHVwCKMkyNknqNkUrXADTzX1nPTFTFrLXjvhxb3+1p8telk2SySyo4KF/rWo2PJaWXcorn25tswUtA4dTC1i9I4XlPyfx0v82NWlBAixadYADJ/KQZZn/e2c9T725lrySaiRJQpZlcour+e+qA4yfs5iqOuN5Td6TEeQuHoC8c95XrN2VgkIUeGPOTTw65TKEForP0sp6YqYuxGS2Eu3n+H1OOxb2pAlP3X0VT99zNQE+7qzYdJSn3vyReoO5mQHa0iUKAvEDQnnhT2NRKERENz/i73u/SU4v0GBWtVKBUimid0LnbXvelz0U1iqo2vr3DqdzwbZpskkNv/Tjp4p45/8moVIqOJiSz+tLt1NdZ6Sksg6rteXYwp3NOyv28O43e4mJ9CctuwxBaDMOY4tIssyh1AJe/O8mIoK82ZVa3WJUeACz1Yb5Aj1vV3PB9wZbsu4wp/IrSBgYztsr9iAKwgUvUmRAlmTSshuiPXWkLJCB5MwSTmSVkl1Jt+8T7AouuOlkGXYeOc3OI6cBsHWDYRpnI/02bHCRZGRt0uWmu35kPzLzm06RKSqvw2K1ERns3eI1AW6Ozxr28FG3S58jOKqrKzRB+96XPVye0Msp6XRaQ0KSJAoKCto+0Q5ycnIcvqZXL+e8oNZwVFdXaIL2vS97OJ/+sLAwh6aedVpOV1BQ0GlxTVx0L3Jzcx0KD9cjcjoX3RtHc7ou76dz4aL7DMg5yDub0xn/1g5ueGsnk97bTW6FfdOfBEFonIoFUFhYiEKhYP78+QD8+OOPxMbGIopi43R8e0gZEkvm5Ns5dfNEsm6fQuU33zj0PG3pWrhwIYMGDWLYsGFMnjy5zcCTkiSxaNEiFi1axIIFC3jzzTdZtGgRn332mV16EhMTqa93/toU6AZdJu3hUHYF+zLLWf/EKJQKkcJqA+4q+x7Fz8+Pffv2YbPZUCgUrFy5kiFDhjR+PmDAAFauXMkjjzzikCaFpyd9Vn0PgCU/n9w5c0AG37vudIquhIQEZs+ejVar5YUXXmDBggW8/PLL501PFMXGZ9i+fTvu7u6MHDnS7udJTEwkLCwMDw8Pu6+xlx5putJaE95uqsZI66He9k8nEgSBUaNGsWPHDsaOHcuqVau4/fbbGz8/d7pWe1CFhxP8zLMUv/663aZrS9fo0aMb/x8fH8+6desc1pWfn8/GjRsxm834+vpy2223oVar2bhxIydPnkShUDBkyBACAwMpKCjgm2++QavVMmvWLIfv1Ro90nRXxwTy1qZ0rlu4nVExgUweEU5cpI/d1995550sW7aMgQMHolarCQgIoKzMuXtNaIcMxpzlWPhUe3UtWbKE6dOnO5S2LMts3LiRu+++Gzc3N/bs2cOBAweIj48nOTmZp556CkEQMBqNaLVa9u/fz4QJEwgKCnLoPvbQI+t0Oo2SdU9czfxbhqBRidz76X52njz/XLhzufLKKzlw4AArVqxg6tSpnSOyHe0ze3SdWfh01113OZx+cXExS5cuZdGiRSQmJlJdXY1Wq0Wr1bJ69WpSU1NRqzu/A7tH5nTQEBpiVEwgo2IC8XNXszmlmGv6B9p1rSAIXHPNNbz++uukpKTw1VdfOV2fMSUV9W/T+O2lLV1r165l2bJl7Nixo12aQkNDmTlzZrPjf/7zn8nMzOTo0aMcPXqUO++0r0rQXnqk6U6V1qEUBaL8PZBlmZPFdfQP1jmUxuzZs4mLi8Pf3/nrGiyFhZS88QZ+0+9x+Nrz6Tp06BBz585ly5Yt6HSOPesZqqurKSoqIiQkBLPZTF1dHTqdDovFQv/+/QkJCWHJkiUAaDQaTKbOmeLeI02nN9l4ac1x6owN89KGhnsz88poh9KIiYlpsdGwYcMGHnzwQUpLSxk3bhxjxoyxKye01daSedtkZIsFUaPBZ9rd+LSj6D6frmeffZaamhomTpwIwFVXXcUHH3xgd7qCIDBlyhTWrVuH2dwwNnv99dejUqlYsWIFNlvDbIRx4xoWKMXFxbF69epOaUi4OodddDk9siHhomfjMp2LLsdlOhddjst0Lrocl+lcdDku07noclymc9HluEznostxmc5Fl+My3W8sWbKEhIQEPD098fX1ZcSIETz99NPNzisuLuapp56ib9++aDQafH19uemmm9iwYUPjOd9++y233HIL4eHh6HQ6EhISOmVSQU/FNQwGvPbaa7z44os888wzjBkzBqPRyKFDh/jiiy/IyPh9q860tDTGjBmDTqfjL3/5C4MHD6ampob169ezePFiDhw4QFxcHFdccQXR0dHcdtttBAYGsn79ehYuXMi7777LnDlzLuCTdg9cpgPCw8O57bbbmg2gy7LcJLjPJZdcgslkYs+ePXh5NQ2OffToUXx8fOjVqxdlZWUEBAQ0+fyee+7hl19+IcvBiZ1/RHrkLBNnU1VVRUhISLPjZxtu586dHDp0iDVr1jQzHMCwYcMa/3+u4QBGjBjB6tWrnaS4Z+Oq09Gw5uC9995j6dKllJe3vCvMjh07UCgUjVN/HGXv3r2dsoVVT8RlOuCDDz5Ap9Nx//33ExgYyJAhQ3jppZeaLPPLz88nMDAQNzfHYwpv2bKF1atXN1lieFEju5BlWZaNRqO8Zs0aefbs2fKgQYNkQI6JiZFra2tlWZblhx9+WA4NDXU43aysLDkoKEi+5ZZbZEmSnC27R+Iy3Xn45JNPZEB+++23ZVmW5ZdffllWKpWywWCwO43y8nJ54MCB8iWXXCLX1dV1ltQeh6t4PQ8PPvggfn5+pKY27KgzevRorFYrW7Zsset6vV7PxIkTMZvNrFu3rlMWLfdUXKYDSkqab8FZWlpKdXU1wcHBAIwaNYqEhATmzZtHbW1ts/OPHTtGbm4uAFarlTvuuKNxG6vOWDvak3H10wHBwcHceuutjB8/nqCgILKzs1mwYAFZWVkkJSXRu3dv4PfOYS8vL5566qnGzuENGzbw8ccfs3//fuLi4njooYf4+OOPeeedd5qFchgxYgQajeZCPGb34UKX792B999/X77++uvl0NBQWaPRyFFRUfK0adPklJSUZucWFhbKTzzxhNy7d29ZrVbLPj4+8vjx4+Xvvvuu8ZyoqCiZ30IXn/svKyurC5+se+LK6Vx0Oa46nYsux2U6F12Oy3QuuhyX6Vx0OS7TuehyXKZz0eW4TOeiy+kRkzjPnjJuL/369esEJU1xVFdXaIL2vS97cJZ+V07nostxma4d5BZX85e31pGeV3mhpfRIekTx2l2QZZmvNiTx9Dvrqak3sWKjmo//ej0De3XSVtJ/UFw5nQO8+N9NPPjP75lw5QBOrHiKiEAdD76xkeQs524H8EfHZTo72Z10mje/3MPLD41j8YtTiA7z5X9zryc62Isn39uOxbWDsN24TNcGBrONL3dncfeL3xIc5ItbcAhnJuZ4uqt5+U9XUlypZ90+13pWe3HV6c5DvcnK57tO88m2U5xOTcdQXc9V113Jwp/SMVgk/nrzQAD6hvswengEn/2UzC1X9kUUhTZSduFwTnfw4MHO0NGtqKgzM2beap5/cxUn9+5DX1LCf+bcyKaXxvO3Wwfx4eZTfLI9s/H8P02IJbOgmh1JeRdQdc/B4Zxu/PjxhIeH88ADDzBjxgwCA+3bpaanUFFnZvyLa0g+eJSEwZHcMeZKrhwWxaWDG3ZunjW6D5X1Zv65OoXckf7ck+DPiJggRsQEsnj9ccaMcO3m3RYO53SFhYXMmzePjRs3EhUVxeTJk1mzZk3j5hc9meO51Yx/aS0nDh3jxisHsv2DP/Hk3Vc1Gu4McycM4C839mfZgXKeX5tHpd7KnybEkphRypH05ot8ehKSDKUWBSap86oJHZquXlBQwOeff87ixYupqanh3nvv5ZFHHnH6cI8zh8H2J+dyMCWfIF8P/L3dySmq4tfUQjYfziavqALJYuG6kTGs+vc9qJSKVu/x7fYkXttYiCgKzL8xjL++s4FBUX688cg1DmlyNu15X5IMZRYFBWYVFllALUgMdDejEX+3h7P0d6ghkZGRQVpaGsXFxVxyySVUVlaSkJDAc889x/PPP+8Ugc6irKqeFxZtYum6wygVIlab1PiZQqPBTefBjVcP5u5r+3PbtYPbNBzAiAgPPrormvnr85m7OpdLh0bxw9ZkyqoNBDiwB213IM+kpNiixF9pI1BlI8uoIk2vZqC7CbWT+zgcNl12djZLly5l6dKlWK1WZs6cSWJiYuMyvXnz5hEfH9+tTFdYVsvlf/oQs8XGgqduRuMXwCdbM8gvqyMhJpBpV0Zz8/AwdFrHf4P+HkoWTI5kwZYitqZaEASBlTvSeeSWYW1f3E0wSgIlFiXhaithmob91ga4m0nVq0nVa+jtZsFTIbWRiv04/JYHDhzIpEmT+OCDD7jhhhuahNMC6Nu3L9OmTXOaQGfwxrKdmK02npt9K5/tL6BaX87EEWH8+ZHLGRzu3eH0NUqReeNDCfJU8XGeH59vTGXWzbGNO293d/JMSpQCBKutjcc0osxAdzOnDGpS9WqCVTaibZJTnsnhOl1OTg69evVqLjwvj4iIiBau6DgdqdPlFFcx9O53COodheQTyNSREcwZH0OEn7vTdcmyzJs/Z7Hkm93cMj6ef06LbVFTZ+PI+6qziaToNfTWmglQNW8MyjIUWxTkm1QMCfEiNrR5bD5Hcdi2sbGxLR4/Oyhgd+L1JTsQFAq8Q0P5+Zlr+PfdcU4xXEsIgsD/3dSHiFBf1u9K4bMdOZ1yH2chy5BjVOIuSvgrW+59EAQIUdsY4mFiYFD79pk9F4dN11LGaDQamxWzncnWkzUs2l3CyRJji3rOkJlfwdL1h1EHhvDmjHhiQjy7RN8bsy7HXSnw5pLtzHl/Fxn5VY2fWa02jqQVcCg1n9TTpdTpO2cjX3uotonUSwoiNBba+vq0ouy06oLddboRI0YgCAJGo5H4+PgmnxUWFnLNNS13Ezib5EIDb2wpRCUKfJdUSS9fNVf10XF5tI4BQVoUvw1DybLMI2+sBYWShydfyjUDu64Te2ifQDb/5zZmfvAru5Ky2H4oiz5hPgzpG8quI1lU1Rkbz/XzcuPLV+7i2njHtl7vKLIM+SYVngobXk5sJNiD3XW6pUuXIssyjz76KIsWLWo8LooiwcHBjB07FqWyc4Zyz9RRKvVWHvsmm2AvFf++JYKj+Qa2ptdwILueGqMNT41IXLg7/QK1rN6VQcrRNGIvjWP3v29Do2q7C6S9us5HjdHGyz/lceBEIYaKcjwUMndcO5A7xgxC566hpt7Ia0u2s/PIad76y83MuvXS86aVX1rDj7tT2ZOUzYO3XEJmHaw5XEBVvZk6k40JcSFcFenGoRO53HxpMO4aVavaKiwip4waBrqZ8FTaZzpn1UkdbkgkJiYyfPhwp9zcXjIyMrBJMs+vzSOr3MRHd0YRoGt4qUXl9fh6ackoM3Mot54jeXqOnSqlJOUEE8cM4/O/3YZW7XzDndFlD9UGKwdy6vnuaC3ZZXr+dG1vEnr7IsugFGS+XLufL386zJ3Xx/HwXaMoqjFzsqiWU8V1VJVXcuJYBpnZRSgUIkG+Okqq6vHuN4Bbr4oh0s8dGVi06gClpzJBlgny8+Af91/O1UPDW9Qjy3Bcr0EtyAxwN9v9vF1quvXr1zNhwgQA1qxZc97zbrnlllbTkSQJSZJR2tHxejZpJ9N5Y3MhO07V8vqkSIZHuCPLMotWH+XD1UkEeLsx6co+9A71Ji2ngp8OnGZgdBA/v/OAXZ287cXRVnVkdG8+3pbJ+5syMFma5i76slJqck6j1GjQhYajlcwYq6uoqqzB3VOHe1AwoqcPgiBQn5UOZiMrXr0Lg8nKpv3p/O+HXxk4MJpypRfV2VmYa2uJv2QQT9wymPgI9yZ17hKzgmyTmkHuRnQK+/OcLjVdbGwsx48fB2jsBG6WkCCQmZnZ4mdnOH6qmEvv/wB3rQovdw23XjuYt/5yc6uNEEmSmf3JHjamVjNvfBjX9vNEkmT+s+IgX2xKYdbNsRhMVn78JZPqejNRwZ7E9g7gvWenEhrQuQ2H9q4GM5htGMw2BKFhClVehYGCSgOVlTUs/HQjWQUVBPl6MGpEb+6bMILrR/ZDkiG/wkBOuZ4+AVqmPrOMIycLAfDy0PC3B8Yw584rMFokNu1P5p+f/8KJU8X4DRlGdKAH0y/x59p+ntRKCtINagJVNqK1lnbp7yhdGiqsvFrPut2p1OhN5BZX8+7Xe/nbA6N54U9jWzz/VEkdr/2QzIbDuUyN9SDYTSY1u4KkU6Wk5Vbyt3sv466xAwCwWG2YLRIebg3Fbk9dgmgyWyksqyUq1KfVH2Od3kTiyUJ6h/sRFuDZ5NyMjAzKqg1Mev4HLhsaiUdkNAey6wn3UdM/0od+gVrG9VLh5+FYHfyCme6NN97g3nvvJSwsrMM3//fnO5j/8RZef2oioRGh7EorpazGiLG6mlOZ+WTllmE1GhoqITT0GfUJ9WZItD83jIzmmrjzd0Y78oIkScYmy6gc7BLo7utev9iYwhsrfmXFSzdjVOv44nAl2SX11BisuKtEFk6OpF+g1u50L5jpbrrpJrZt28aoUaO47777mDJlCu7u9ne2WiWZzOI6FqxP49DpCrKT06gvK0Op1aLVajAb9JiNZjTubgzuF8a9Ywfgp7UR4udOkI876jZaobIsk18vERjS8KMI89Li637+lpzRauPNbafJrjRw06BAxg8IwN3Ohkdnm06WZYxWiTqTDbVSwFvbeov0XF1Wm8Qdf/+RGoOFXpFB+HlqmDYqmugwX+b/VEBJrYXXb+2FWRDJr5eoMcl4qQW81AI2GQxWGbVC4JJgJRqFcGGL15KSEpYvX86yZctIT09n8uTJzJw5k+uuu67V6/ZnVvJ/K5LIL63HV6cmvo8vnhqRXw9noJAsyFYrwX46po0fRu9egaSV1hOgUxMqVVGqlyiolxgeqEStaF7sSLLMiQobu/ItFOt/r6S7qxU8f10fevs3/2EYrTYWbssiq9zAyChvfsmqQqEQ8NQokSSZAJ2amwcHMiLcq8WirjNNl1Jcx4e7s6k0NIyHCgIMC/Xk2n5+aJQi5fUWZBmGhOgI9lQ3K17PsC+rmsXrT2DS6ykrqyGnpJbbru7LnybF8dKGEqpNMkMHBOLjrsRXI1Brlqkxy4gCuCsF9FYZN6XAdZEqplwxENEJgwAdrtMlJyczd+5cNm7c2OZEzqTcSu754AAhQTqmXxlJLz93vjxUgKdWibtKQVm9GYNFQiEK9PFzAwGq9BZqjFZMtgaZl4couTG6aaBoqyTzVZqJU9U2+niLXB2mZnj/aKySzH9/yaWw2siz1/Wlb8DvxtObbby94zSnyvU8M7Y3A4J0VOot7DxVgdEqIQoNX/zJUj29fLTMuSaaUK+m9+0M08myzOaT5XxxMJ8BQTpG9/PDU6OktN7M9oxyMssNQIMJRcAmQ6CHGk+tAkkGb62SoV5m+norsCGwr1KFSrKS4CchAt/tTOftlYdx16i47vax7E4uI8xLxcLbIvD8bZaNJMsINDQOq0wSm3PMHC+3MWVYMJOHNd9DzVHabbrCwsLG3O706dNMmTKFxYsXt3md2WpjbXIpq44VAzA2xp/pCWFolCKyLFNUa+Z4YS3ppfVIMthkmUq9mRgfBRaTic25FmYN0RLh2VAE2mSZlekmTlbauHuAhhifhhd35gvWm238Z1smeVVG7hgeyui+fpTWm3lr+2lqTFaeHh3d6phianEdn+zLRa0QmX9TDOqz6n32mE6WZSpNMnl1EnqlJwU1JgqrTZhtErIMGpVImJeGUC8NlQYLOZVGCmtMjB8QwD0JYSjPWehTVGtCKQj4uqsw2yRSius4UVSHySohCgKZ5XqyKgwEuQkYbAK15oZcXyVCtJeCGB8FGlMdM+avJf6aETw+vh//+imPCB81r98Sift5Js/l1NpIGNgHbzf7ivjWcNh0X3zxBZ9//jk7d+5kzJgxzJgxg8mTJzu8Z1ZyUS1WSSYurO1ZC3tPZJFtUtNXa+TbVD02WeahWDdMNtiYbeZomZW7B2gY4Pt7a+zsXMVgsbH0QD57TlfipVFiskkEuKt5anQ0IZ5th9fPqTTw95/SubafH/eP/L3x0pLpiuptJJVZya2VMFhl6i0yxt8KgCCdmggfLWFeGrQqBYLQ8KPIrzZSVGPCx01FL183YkN1xEe0b8qVLMtsOnySA8UWjKjo5y0S7SGRVW3jVLWNnFoJSYZ9m/Zhqapmwxu3kVlu5q+rc+kXoOHVmyNwO8d4VXVGPvs5hRk3X8blQ5vPMHIUh003fPhwZsyYwfTp01vcrrIzSE/PIN2gpl4SsZjN/HSyFk+NSK2p4Vd8W18NcYFNm/8tFWVFtSbWJZcgyTDjkjC0DgyNbT5ZxpID+dwRF4KbWkGt0UppeQUyYLY11IPKDRJlRhl3JfTzUeChEnBTCoS4i0ToFAwbFNPqPWRZprTOjFmSiOjAzOOMjAyyjCqqrAqGehhRnpVZGq0yGdU20nPKeendDbzxyChuuqw3yYUG5q3NI9pfzb8mRuChUVBdZ2LR2qN8u+0kZpvMXZMu57O5N7Zb1xl6REj/jIwMLBKcNqmRZDhZZiSv1kq0t4qrgsBL3bxy6+zuCVmW+WhPDntPV6EQBTw1CpBsCIBKAd5qES+1QIxPQxGmaGH9a2uacqsMnCiqpcLQ0GF7dW8/In3aZ7xDqZmcNKjppbEQrD5/PfvhBZsorzHy7T8mIggCqcUGnl+bR6iXimkjfFiwdCf5pXWoA4IZHhvF0jnXEtpOTWdjl+kWLVrEI488AsC777573vOeeOKJDgtqiZaKsWqrSIZBjZsoE6mxoFNITabndEaf2JkuDK1SRBAEpzUkMsrq+TW3imCdhoHBOjLL9RTWGLm+fyA+Dtah6s1W1p8owkOUiHEztzpl6UBKEQ++sZH/PHoNN46MbtBSauS5H7LJPJqM1WDAb8Bgpl8VwQOXBTBwQOs5tb3Y1SW9Zs2aRtOtWrWqxXMEQeg007WEt1JioLupYTq1QYNWlAhUWQlQ2ZoUJ85EEATcnDxbJa/KwMHcKvoHehAf7o0gCAR5qNl0spSdmeWM7ReATvP71yTLMrUmK0W1JqoMFqySjFWS8dYqCdJpOF5Ui4hMb23rhgO4dGAw4xJ68fz/dqNWioyN74VGtqArO43CbOSlh65haJ9A+gbY34FsDz2meD0fsgw1NpFSi4IqqwIB8FfZuHZQpMMjDM7U1RKR0b2pNFioNlixSBJWm8zJ0jrCvLVcGe3XpA+szmRl88lSjFaJSB83gj01lNabKKk1o7fYEAXw1qpQKQREQaDSYPmtBQsD3OwfyLdYbTz3v91sOZTDDSOj2fhrNv5eWl576GouHdi0zn7BOoeHDh3KsWPHmh0fPnw4iYmJThF1LvZ+uRYJSi1KCs1K+gXqSIjw6RQ9juiSZKiwKigxK6mXGn4EogBqhYhCFPBzV3NFlG+LdUCrTSKzQk9aSR11Zhu+biqCdGqCPbUE6dRNflSyLFNttCLJMhUFjk2Tt9ok/r54L5sO5fDghFjuu2EwbprmheAFM52np2eLW0/6+flRUVHhFFHn4miOUmhWkm9SccPAIHyd0K90PtrSZZEgWa/FIgt4K2wMjgjAz12FTqN0qGdflmVskv3TxduzkEmWZSxWqdVhxi5fbP30008DYLFYGv9/hszMTPr06drp1q0RrLJSgxuHcqu4LiagS9dvnE2eSYUkQ6yHETdRJrqdC4IEQUDZwtCfMxEEoc1xbWdhd6WnsrKSyspKJElq/H9lZSXV1dUMHjyYb775pjN1OoQoQHyEN6X1ZrIrDRdEQ51NoMyqJEJjwU3s9tXmLsXunO6zzz4DID4+njlz5nSaIGcR6qUlwlvL0cIaevm6OWWg2l5kGbKNatxFicAW1pJe7DjcvBswYABpaWlNjqWmprJp0yaniXIWg4I9qTfbKKzp2mV+ZRYFekmkl7btpX0XIw6bbs6cOXh4eDQ55uHh0S1zP393FX7uKtJL67r0vqUWJd4Km1Pjf/yRaFd8unPDR0RGRpKfn+80Uc5CEARiAjworDVRa7S2fYETMEkC9ZKIv6tYPS8Omy48PLxZf1xSUpJTpq93BlG+7qgVIullXZPbVVoVCMj4nCdMg4t2mO7RRx9l6tSpLF++nP3797N8+XKmTp3aOEzW3VCIAn0D3Mms0GOxdX5xV2FR4K2U6OQejh6Nw0vyz9Td/vnPf5KdnU10dDSPP/44cXFxThfnLGICPDhZWs+GtFIuj/IhwKPtOXTt4UzR2kdt/wLmixGHc7ozA/snTpwgPT2d6dOn895773HzzTd3hj6n4KFWcuOAQNQKgc0nyzheWNMp96lwFa124bDpbDYbP/zwA5MmTSI6OpoXX3yR6dOnk5PTvcNieWlVjOsfyJAQT44V1XKyE1q0la6i1S7sNl1aWhrPPPMMYWFhzJgxA19fX9avX09gYCCzZ8/G39+/M3U6BVEQGBrqxYBAHYfzqsmvdt5ohfG3otXPlcu1id11ukGDBuHv78+///1v7rrrrsa+ugs1rtkRRoR7UW+2sud0JeP7Kx2eKNkSFRYFoqtotQu7c7r7778fo9HI3/72N+bPn98Y26QnIggCV0T74qFWcCCnEskJUworrAp8lDZX0WoHdptu8eLFFBUV8corr7Bnzx7i4uKIj4+ntraW6urqztTYKShFkct6+VCut3CytL5DaeltAgZJxM/VIWwXDjUkPDw8mDVrFnv37uXYsWOMGTMGDw8P4uLimDVrVmdp7DQCPDTEBHhwtLCGOlP7RywqrAoUyHi7hr3sot3zuQcPHszChQvJz89n2bJlFBYWOlNXlxEX5oVGIbI1o4xDeVXkVRtajWN8LrLcUJ/zVTVMIXfRNh1eRKBUKpkyZQrr1q1zhp4uR6UQuaaPH0E6DfnVRnZlVpBRZn9xWy8JmGTxvNHJXTTHtd8r4Ouu5vIoNQB7T1dworiOPv4eLa5bOJdSixKVILtmlDhAz9jSpQsZEuyJ3mIjq0Lf5rl6m0CZRUGI2jVvzhFcpjsHbzcVvXzcOFFc22pXiizLZJtUaEWZIFer1SFcpmuBISENM45by+1yqgzU2RT00lhcDQgHcdXpWsDHTUWkj5YDOVUkF9US4KEmSKch2FODKEBelZETxbX4KG1427kHg4vfcZnuPFwe5UsvHxNl9SbK6s3kVBo4U9iKAoR4agiwtV3vc9Ecl+nOg1IU6eXrRi/fhihFZptEaZ0JmwQhXhrUCpGMjKoLK7KH4jKdnagVIuE9bLfq7oqrIeGiy3GZzkWX4zKdg5zIKqHe4FoD0RFcdToH+N+qAzz55o9o1UpuuWYQd46KZECvgAstq8fhyuns5IufjvDkmz8CYDRb+W7rcV5ctBGDybFN3Vy4TNcmJrOVd1bs4aHXmoa9tUkyRpOVn/edvEDKei6u4rUVlq47zN//t5mSijpaGoWVge+2JjPxyoGouii22x8BV053Hg6m5PPI6z9QfB7DnaG6zsj2I1ldpuuPgMt05+HTNb/aFW5VAL7ZfBRJcgU+tBeX6VqgVm9ixcajWO2IfSIDReV1HErN63xhfxBcpmuBbzcfw2i2f6GOQhQ5lFrQiYr+WLhM1wKfrP7VoXCxNkniYEr3i8/XXXGZ7hyOZRRx5GShwwuwiyvqKK7o2oifPRWX6X5DlmU27U9n5svfomjnTjuJaa4i1h5c/XQ0NBxu++sy9h7NQSEK2NrREhUFgSNphdxwRf9OUPjHwpXTAc+9/zP7jucCtMtw0LAF+eGTBa6uEzu46E23+UAGi9cecopZDCYLGfmds1XVH4mL2nRVtQZm/fN7p21sIgoCiWmuVmxbXLSmM5mtPPDKd5RV1TslVBg0NEYOubpO2uSiNJ3BZGHKc8vZuC+93XW4lpCBtJxyTGbXdKfWuOhMZzBZuHXuMrYdynRaDnc2Nkki5XSZ09P9I+F00+Xldd8xSKvVxr0vfcOepOxOa2UqRJFjGUWdkvYfBaebbvDgwc5O0inIssyTb/7IT7+kdUoOdwabJHH4pKuTuDWcbrrk5GRnJ+kUPvpuP4vXHqIT/dZIZl6Faxp7KzjddJGRkc5O0il8vPrXLruXJMucyCrusvv1NBwaBvv888/tOu++++5rl5jOIr+0htTTpV12P4UocjS9mISBEW2ffBHikOneeuutJn8nJyfj7u5OaGgohYWF6PV6YmNju53pNu5L79L72SSJIycLeICELr1vT8Eh0x05cqTx/y+88AI333wz8+fPR6lUYrFYePnll50u0Bn8/MvJdg/kt5fTBZXUGczo3NRdds+egiA7Ekr8LAIDAyksLESp/N23FouFsLAwSkudW5RlZGQ4fE2/fv0AMFushE54Db3R+RX7aL/WP5977yiuHdG7mabOpj3vyx6cpb/dDQk3NzcOHz7c5FhiYiJarbbDopzJvuO5nWK4thAFge+29NxdhTqTds+ne+aZZxg/fjz33nsvUVFRZGdns3z5cl555RVn6uswG/alo1SIdi2ycSaSLJNVWMmh1HwSBoZ36b07Sk5xFWWV9cR3ku52m+7xxx9n8ODBfPnll2zbto2wsDBWrlzJdddd50x9HebH3aldbrgziILAio1JPcZ0x04Vs2JDEkdPNYyozJs5miuG9XL6fTo0c3js2LGMHTvWWVqcTnpOGSdzLtw4qCTLpGaXkZxZzJA+wRdMhz2UVdXz4qJNjbsFCcDCL3fxVvBEIoO9nXqvP/TY63++2Nnu9Q7OQhREPv7hVypqnLe3bGewdncqsiw3DhHKgMUm88rireiNzg2N9ocde80uquLLDUnYLlDRegZJlsgqqOSR11ex/OdEJKn7RWM3mCys39N8TFqSJIrL63j5k60YnThdq91dJucjNzfX6UNh7ekCeH9NKp+s/rVT++ba6jI5l9MV4OOp5fqR/Rg5JBJvDy1+3m5cP7IfapXz1kg5+r5W7zjBp2sOnjdmiygIDOodyOevPYS7tuP9ju3O6RYvXty40XBSUhLDhg0jISGhy/d+NZgs1J0TGbOyxsDitQe7tDPYXqpqjXy//QTPvvczD722iqnPfcmsV79vc+fF8mo9JZV1WK3O3Z3HZpNYtSO51SBBkiyTklXKPS9+7ZR7tvvn9eqrr7J//34Ann32WW688UZ0Oh1z5sxh27ZtThHXGicyi9mwP51didnYbBJD+gQxckgE5dUGDqfmd+tVWecW+d9uPU7CoHCevPuqZudmF1byz8+28cXPSY3GDAvwZNoNcdw1bhgebmoKympIySply68Z7E7K5tm7hnLTFQPs0rL7aDbl1W3XNyVZZvuhTLvSbIt2F69eXl7U1NRgNBoJDg6mpKQEpVJJYGAgFRXOXRF1bnFxIDmXVxZvQyGK2H6rI4kIyMgoRBGrJHG6CxZltad4PR+CAOvfup/RCX0aj7311W5e/O9moLlRzx3WEwDxt2N9/ARefngccf1D29T05MK1nC6oQmo1r2ugsFZB1da/t3leW7S7ePX19SU9PZ2ffvqJhIQENBoNFoul0yvKNfUm3lmxFwEaDQcgISMD1m5YUbcHAYHbnvmCD1fuQ5IkXvhoI/M+3IjNJrXYGDq36iCfc+xfS7dTVF7b6j2Pnyoms6DSLsM5k3YXr08++STDhw8HYOnSpQDs3r2bQYMGOUXY+fjfqv3UGcxd/Jo6H0mWMZmt/N8763n7qz3klrS/biwhYzTbePnTrbz3f5PO22303dZjiILQqTOpW6Ldpnv66aeZNGkSCoWCPn0aioRevXrx8ccfO03cufxyLIcdR053WvrdhfzSmg6nIUkSucXVbD14iusvi2n2eW5xNQcvUHizDrXTY2KaPkz//p0Xx6Om3sR73/yCAH+4XO5cnJXzCMDn649wzYhoNGpVk89+2H6iSZ24K+kxSxD/t2o/9X/AYrUzkWmIibx2d1qT4+v3prHpQPoFMRxcoKhN5dV6jmUUUVJZz63XDEKjbl3GvuO5F0Wx2hnINMREHto3GB8vd9bvTuX77Rd28VSXm+6FRRtZuHx3498zJ4zgo+duQ2ghnogsy2w6kMF7X++9KIrVzsJktjH33Z8utIxGutx0J3PKm/y9dP0RhsWE8tjUyxuPZRdWsu1QJh9+t59jGUX08RdchusAXd06bYtuERTxr+/+REFZDacLKtl3PLex9SaKDblfd3tpLjpGtzCdjMzC5bsRBJoshu7OQ1ku2k+Xt161LTQazhjNlaF1b9pq8NmL06c2nUGSJAoKmnc+llTUceyUYwFm6qocn/2r8+n8LTEd1dUVmqB978se+vfry6Degc2Oh4WFIYr251+dVrwWFBR02xATLpxLbm4uERH2RzPo8pzOxR8PR3O6TjOdCxfno8cMg53LO5vTGf/WDm54ayeT3ttNboXerusEQeCxxx5r/LuwsBCFQsH8+fMB+PHHH4mNjUUUxcaZ0faQMiSWzMm3c+rmiWTdPoXKb75x6Hna0rVw4UIGDRrEsGHDmDx5MjU1rU8KkCSJRYsWsWjRIhYsWMCbb77JokWL+Oyzz+zSk5iYSH19vUPPYC/dosvEUQ5lV7Avs5z1T4xCqRAprDbgbucaAz8/P/bt24fNZkOhULBy5UqGDBnS+PmAAQNYuXIljzzyiEOaFJ6e9Fn1PQCW/Hxy58wBGXzvutMpuhISEpg9ezZarZYXXniBBQsWtBo7RhTFxmfYvn077u7ujBw50u7nSUxMJCwsDA8PD7uvsZceabrSWhPebqrG/VhDvd3svlYQBEaNGsWOHTsYO3Ysq1at4vbbb2/8/NyZM+1BFR5O8DPPUvz663abri1do0ePbvx/fHw869atc1hXfn4+GzduxGw24+vry2233YZarWbjxo2cPHkShULBkCFDCAwMpKCggG+++QatVsusWbMcvldr9EjTXR0TyFub0rlu4XZGxQQyeUQ4cZE+dl9/5513smzZMgYOHIharSYgIICyMud2M2iHDMac5diO1/bqWrJkCdOnT3cobVmW2bhxI3fffTdubm7s2bOHAwcOEB8fT3JyMk899RSCIGA0GtFqtezfv58JEyYQFBTk0H3soUfW6XQaJeueuJr5twxBoxK599P97Dxpf6SoK6+8kgMHDrBixQqmTp3aOSLb0T6zR9fbb78NwF133eVw+sXFxSxdupRFixaRmJhIdXU1Wq0WrVbL6tWrSU1NRa3u/NBmPTKnA1AqREbFBDIqJhA/dzWbU4q5pn/zjsuWEASBa665htdff52UlBS++uorp+szpqSi7tOn7RMd0LV27VqWLVvGjh072qUpNDSUmTNnNjv+5z//mczMTI4ePcrRo0e58077qgTtpUea7lRpHUpRIMrfA1mWOVlcR/9gnUNpzJ49m7i4OPz9/Z2uz1JYSMkbb+A3/R6Hrz2frkOHDjF37ly2bNmCTufYs56hurqaoqIiQkJCMJvN1NXVodPpsFgs9O/fn5CQEJYsWQKARqPBZDK16z5t0SNNpzfZeGnNceqMDVueDw33ZuaV0Q6lERMT02KjYcOGDTz44IOUlpYybtw4xowZY1dOaKutJfO2ycgWC6JGg8+0u/FpR9F9Pl3PPvssNTU1TJw4EYCrrrqKDz74wO50BUFgypQprFu3DrO5YXH69ddfj0qlYsWKFdhsDYu4x40bB0BcXByrV6/ulIaEq3PYRZfTIxsSLno2LtO56HJcpnPR5bhM56LLcZnORZfjMp2LLsdlOhddjst0Lrocl+lcdDku07noclym+40lS5aQkJCAp6cnvr6+jBgxgqeffrrZecXFxTz11FP07dsXjUaDr68vN910Exs2bGg8Z/369YwaNYqAgAC0Wi39+/fn5Zdf7rQB9J6Ga+wVeO2113jxxRd55plnGDNmDEajkUOHDvHFF180iXeclpbGmDFj0Ol0/OUvf2Hw4MHU1NSwfv16Fi9ezIEDB4iLi+OLL77g+PHjjBw5El9fX3799Vfmz5/PAw884NAg/R8W2YUcFhYmP/bYY82OS5LU5O+EhAQ5NjZWrq6ubnZuUlKSnJ2dfd57zJs3T/by8mqW5sVIj5za5GyqqqoICQlpdvzs8GU7d+7k0KFDrFmzBi8vr2bnDhs2rNV7+Pv7Y7F0/Rag3RGX6WhY6PLee+/Rq1cvJk6c2OLEzh07dqBQKBrnm9mDzWbDbDZz+PBh3n33XR5//PEW4/BddFzorLY7kJSUJPfu3VsGZEEQ5MGDB8svvvhik2L04YcflkNCQhxK19vbW6YhlqM8c+ZMV9H6G67WKw1FY0pKCmvWrOGxxx5DlmVeeeUVLrnkEurq6hrPczSX2r59O7t27eLNN99kzZo1zJ4929nSeyYX2vXdlU8++UQG5LfffluWZVl++eWXZaVSKRsMhnalt2zZMhmQMzIynCmzR+LK6c7Dgw8+iJ+fH6mpqUDDYmer1cqWLVvald6ZjV6yHFwL+0fEZTqgpKSk2bHS0lKqq6sJDm7YkXrUqFEkJCQwb948amubb3907NgxcnNzz3uPPXv2ANC7d28nqe65uDqHgeDgYG699VbGjx9PUFAQ2dnZLFiwgKysLJKSkhqNcqZz2MvLi6eeeqqxc3jDhg18/PHH7N+/n7i4OKZMmcLIkSMZOnQoKpWKPXv2sGDBAiZOnMiKFSsu8NN2Ay50+d4deP/99+Xrr79eDg0NlTUajRwVFSVPmzZNTklJaXZuYWGh/MQTT8i9e/eW1Wq17OPjI48fP17+7rvvGs957bXX5GHDhsk6nU729PSUhw8fLr/zzjuy2WzuysfqtrhyOhddjqtO56LLcZnORZfjMp2LLsdlOhddjst0Lrocl+lcdDku07nocnrEfLqzp4zbS79+/TpBSVMc1dUVmqB978senKXfldP1UKxWG2VVnbPPQ2fjMl0PxGyxcsfzXxJ923948NXvOJHVfMJCd8Zluh6GzSbxp1e/Z+uhTJ6860p2HMki4b73ee6Dn7FabRdanl24TNeDkCSJOQvWsmp7Mp/Pv4N/PjqeEyue4l+Pjee9b37h9meXU1lruNAy28Rluh5CvcHMffO/Zcm6w7z/11u4ZGhvympNqFVK/jLtatYuvI+DKfmMm/0pBpP1QsttlR4xy+Ribb2W1ZpYfSifr/Zk8euew5jr6giMiUHh5QuAUhSYOjKC2df3I8LPnRNZJVz95/9y61V9+NuMyxy6V2foPx8Od5kcPHiQSy65xCk3d9EcWZbZm17Ol3tz2JBUQH1JCYbiAlSiwNMPT6BXRBBebioCPNWkFdXyv62ZfPdrHm9Mi+O2hHD+9dh4/vLWOq4dHsHVQ8Mv9OO0iMM5nZ+fH+Hh4TzwwAPMmDGDwED7dqnpCBdLTidJMnOWHuS7bSfQWPSYqquorTdy/83xzHtgDOGBzRd515usvPDtcdYnFrLi8csZHuXD9Y/9l5N5lXz/yiR8dFqHNHZEv704XKcrLCxk3rx5bNy4kaioKCZPnsyaNWsaN79w0T5kWeaVH5JZ/v1uqrNOEayRuPfGOA4tfZwPnrm1RcMBeGiUvH73UIZGevPw4kMUV5t4+U9XYjBZWfLTiS5+Cvtw2HQajYZp06bx888/k5GRwWWXXcbcuXMJDw9n7ty5ndYb/kdn0ZZTfLT6MKaqSr585S4OLJnNgicmMDC67ZJEo1Tw0Z8SUCoEHl1yCH9vN269qi+rdmVgtnS/zKBDrdeMjAzS0tIoLi5myJAhVFZWkpCQwGuvveYsfd0Wi9VGRY2R3JJa6gzmDqW1PrGQ11cfx1qUy23XDmby6CFtX3QOgZ4a3p8ZT1JOFWuPV3HnmP5U1BrZfCinQ9o6A4cbEtnZ2SxdupSlS5ditVqZOXMmiYmJjSum5s2bR3x8PM8//7zTxXYHauqNvPXVHt79ei96Y0NAHD9PLYufHU/fcB9sksTry3/l19Qi7hjdn1uv7ovO7fzbWSbnVfPk4n14VBdjFgXe+svN7dYWH+3LtMt78dn+PBbf05tLBgTzzfaTTLi8ey17dLgh4ebmxqRJk/jTn/7EDTfc0GKohUceeYRFixY5TWR3aEhk5Jbz7ZZjfPTdfmr1Jh6bejlR/krcNEoWfH2IihojH//1ev675igbf83mithQ9p8oRK1S8MzdlzLl2phmmt5a8Qsv/m8TNosVURT49IUp3H1969Gf2qJab2H0q1sYHuFOgreJvy7axapXb6FfuE+H0gXnvVOHTZeTk0OvXr2aHc/LyyMiIsIpos7lQpjuRFYJP+5KIS2njGMZRRw7VYyHm5o7xsYy74ExRAZ7N+oqrzHw4L83crqoBlEUeOPhUYy7JIriSj0f/ZDIdzszuHvsAP734t2olAoA1uxM4a6/fYV3cBBvPzaGGy/tg4+n/dvCt8ZH6w7xxuYi/jI6iH//bxM3XBrFvHtb77fLKaklJaeS8QmR543ZcsFM5+XlRU1NTbPjfn5+VFRUOEXUuWRkZGCy2PhhVwZXxoYRGeTZ5jXtfUE19UZeXbyND7/bj85NzcDoQAb0CuCGy/tz4xUxuGt/LyrP/jGUVun51xcHuGNMf64cEtYkzW+2pfHa8gMM7RvCQ5NH0i/Cn5uf/hxB58l3/7qH62KD26X1fKSnp/OfrUVsSq0hxFxKSmo2X74wgf6Rvs3O3Z9Vw3urj3P0WCayJBEVFcxbj1xNTIhHs3MvmOk8PT2bhVUwGo2Eh4dTXl7uFFHncuToCea8u41DacWIgsD4S6OYdGUfQvw8CPJ1w0OrRqVs2iZqzwvaezSbe//+DdV1Rp6/fzRP3HkFatX5q72O5MBH0ktYuimdrQczkWUZtYcHj91/A6/dPdxhnW2RkZGBLMv8mFzNB9sLqUw7gU4N3788EW8PDQA2Sebl71NZvfkYNrORS4b1JjLEm9VbjiGICsZfNYDnJg/Cz1OLxSqRXVzDsCH9CQtouevGEew23YgRIxAEgWPHjjF06NAmnxUWFnL55ZezatWqDgs6l6LyWm564hPyy+p5c/a1ZBfV8NnPyeSX1jU5T6kQuGxwKG/OvhZ3jarRdCeySvh601G2HjzFG3Nu4oqhzasGsizzyeqD/N8767lsSASfvjiFXsE+bWpLT0/HYmt4fSqF0GYosX79+rEvtZCH3t2KT1AgG+Zdh5taYeebaBlZlqkz2dBpFI33P/vHkF5q5PX1p/l11yECA3y4ddwwqupN7DyURXFuISHBvrz96FUMifIDIKuolic/2k1WTimCINAryJOC8josVomXHx7HX++9pkN6wQHTLV26FFmWefTRR5s0EkRRJDg4mLFjx6JUOncicr3BzNjHPqGorJpF/zeOmIiG4kGSZEoq9RRV6imt0qM3WqmuN/HBD4kMifbng6euQ+sZwJNv/si2Q5n46LQE+HhQU2/kmwUzEbRa3FQK3FUKvDUiL324gc9+PMQjt1/GG3NubKx3nYskyRzIrGDVwXw2HiuiSv97OFeFCDq1gomxPtyT4If6t5y3Um+lziRhlWSOVSj579ZMAnRqPp51KYPD25dr1JmsbM+o4HBeDXlVBvQWiUsjvXnkql5olGKzHNgqSfzjh3R+WLu/8ZigUDBgxCAGx/fHTysyNUZDuO73595yopx/fZ9CfW09bjo3PHQe/OW24Tx8/YB2aT4bh4vXxMTExrBXnYksy9w3/1t++uUknz9/Q4v1kXM5fLKYR97cQkSgjpySOkL9Pfnno+MZd1kMG04U8fj8rxCUSkbfOhqlSklVRS37N++nvrqW6NgB3HT1YMYMDOT62GA0qt+/gFqjhZX781i6+zTZZXoi/dyYFB+GxlqHStmQuxgtMvlVZlYfqyTUS801/Tw5kF1HeunvYfzVCpGHxvbhsXH92pXD1ZttfJNYyK5TFUgyjIjworefG1qlyNdHigj30fD06N6U5Wc3nG+ROVlpZW+hhVKDTIitljCNjWh/N6KDPdG5KakwynyfYaJIL3FFqIoArYC7SiDATUS22tiTVYfZKmO2yUy8LIaE3n4O6z4Xu0y3fv16JkyYAMCaNWvOe94tt9zSajrHs0q47bkVRAZ60ifUm5GDI7j3xuFo1ErMVondJ8vwclNSJ8GH3+7lp58P8v5ztzOqvwfiOUWXLMuNxUmJXmJjtpkai4S2poIvvtvHvRPiefKeURwv1vPDsSJqjFZClFY+/mQ9nu4aDBaJunojolLJwIRhaL09Ka02YbbYCPLW8q+pscRGevPp9iyW783GZJW4ckAAd4yMYMKwEBRi8xwF4HS5ibd3FJNZZmRklI4re+sI1ClRiBAVFUV+rYmMMj1X9fYlIdK7rVffSGpJHR/tyUFvtjFxcBBjYvzx0v5esmSV63lzexY1JhvuStAqoMwgIwN9vEXG9VIT5tGy0a2SzNZcC4dKLJjOGsBwU0KkTkEfbwV9fRRcHhvjlJjJdpkuNjaW48ePA+ePryYIApmZma2mszUxl5mvraGmzgg2K2a9HqVSQVBUBDYPXwSVGkGW0BoqOJ1yiv5xA4i9LBY/rcCYCDWD/UQOZOtZ9msZGaUmNEoBlVIkMsKHvsHuhHmIpFbYMFolFGJDHUsAroj2YUpcCMGeGt749hBvrfwVUSEyalAw/3zwWvqG+SDLMmml9XyxP48ff82jutaEIIBCFAgO1BEaqEPzW+7kpVUyLNQTrbUOb42IWZLJqrZRUC8R669kVLgKpdDwTiRZJq3Sxs58C4X1EiqFQLBOQ161kXH9/bknPgwEMFsldJrm1RNJkll1rJgfjhfTP9CDR6/sRYCu5c7mSr2Fg7nVnC4oQW+VCXEXGeCrRKe23yhWSUZvlSmul8irk8iutZFbK2GT4Y7hIdzqhJZ2l86ny68y8G1SEakFtaQX1lBZXktxVg6VBYXINomQ8EDKymqwWSw8O/NaBsXHUGWwUVBey5F8Pbn51VTWmokO0BIV5E5urY3KaiM1dWb+cVMYIyLc+e/eUtYcq6JXgDsTRoRxx8gI+gQ0NP+3Jhfz0OJDTBweyt9vH4KvR8tfXnppHQs3ZlBUZeLWhDCu7uuHj5sKUYDSOjNJBbUcL6ylsNqA0QYCEOYhEuAmcqzcio9GICFISZlBJqfWRrlRJtpL5Pb4XsSGeKJSCGxJL2f5oQJskoz02zcwNsafmZeGoxAbTFKpt/DhnmxSS+q5fWgIt8YGIYptG8jZ499mm0x2jY3hA3oT6qXpcHoOm+6NN97g3nvvJSwsrO2TzyGrQs/8n9KJ8NHSL9CdAYEe9PZzR4HExr1pfLUxCVlQkGbVMfXKvtw3Koodp8rYc6KUX9LK8fNU07+XNwqtGneVyOUhSuICFPxncyGHcvWEeqkorLEwJjaQkmozSaer8PFQ8/CYPgzr5c2Dnxzk6v4BfPRAQuMX2xEyMjIwWhten1b5e1H/Y5aJ/DqJIHeREHeR4YFKorwUzbpxCqqNHC+qw10lUmWw8k1iISPCvZieEMaOUxVsPlmOWiHw2NVRDArWOaSrM7hg/XQ33XQT27ZtY9SoUdx3331MmTIFd3d3u66VZRmLTW5s2Z2Pr/fl8NzXxxr/1qhEJsSH8KfhOrRKAZskIwg01vPMNolXNxSSW23l9quiGBCgRufhQWpRDZuTijl8qgKbJDOslzcrZl/R4W6KM7T25Uqy3Kwe2taXlphfw3s7szHZJLRKkdH9/LglNrhJ3a2jujrCBTMdNMToXb58OcuWLSM9PZ3Jkyczc+ZMrrvuOqeIAjiRX02t0YpKIWKoLqFQ9kBExk0hY5YElIKMj9KGTiFRblFSYhYRgGg3K35KGzEx/ZBlmdJ6MweyKvgpqYj4vn4Ee2uI8nVnYJAOlaJjS0Q6Y7p6Vrme9NJ6rurjh0c7fxx/SNOdTXJyMnPnzmXjxo2dNpEzIyMDgySQY1QhACqxwXi1NhEZAQUyQWorwSorqt98dO4LstgkimtNFNQYOV1hQKkQiA3xpF9A85axI7ocwbXCv4F29+YWFhY25nanT59m5syZThF0PtxEmQHuTeetWWWot4noFBKKNnyjUohE+LgR4eNGbIgXRwtrOJRXTXGtiSui/VA6oY7nwj4cLl+++OILxo8fT+/evdmyZQvPPvssRUVFLF68uDP0tYpSAG9l24Y7F3e1gsujfLmmjx+FNSa2ZZRhskqNn9eZrOzPrmTv6QqsktRKSi7ag8M53YIFC5gxYwaff/55izsH9iTCvd0YG6Ngx6ky1iQX4eOmQqsUya82olGKWCQZw6lyrunj3+H6X0/FKoNJEjBJAtVGC95aVYfTdNh0iYmJHb5pdyLAQ80NA4LIrTJQabBQZ7ISF+ZFTKAHlXoLO06Vsy2jjGv7+qM5z5hsW0gy1NkazBzmpekxOyEWmpXkmX43mUeVEe+QLjLdokWLeOSRRwB49913z3veE0880WFBFwKdRsmg4OZz9AJ1GsbGBLD9VDkb0kq5urcffu7nn3pulaHKqqDCosD0WwtbFBrqnTYE0jLL8dYqiQ3xJNLHrVubr8YqkmdSEqSy4q+yohFkBgY73jfbEna1XidMmMD69esBGDNmTMsJCQJbt251iqhzudDT1evNVnZnVVBlsJAQ4UNff3cEQWiiq8IikmVUIyGgU9jwEGWsMthkAXeFhI/SRnhEBMcKaymqNREf7s2AIPs7fB2hpfdlkSDdoEEGdL/p8Va2XF81S5Cs1+IuSvR3M3Pmt9Ftuky6ggttOmiY9Hgor4pT5XqCdGoujfShJC8bWYZ8s5JCswo/pZVIjQX1eap/ZzTtz6kkv9rIpMHBnVJXPPd9STKkGdQYJRFvhY06m4hJFhjgZsbrHOPJMqQa1JgkkSHuxsYuqLP1dxSHn/jcCZxn6IrpThcShSgwspcvY/r6o7fYWJ9SwqFaLQfrtBSalUSoLfTRnt9wZzM0xAurTSK1pK7tk51AjklFvU2kn5uZPm4WhnqY8FRIZBrVWM/JcvLMSupsIn215iaGcyYOJ3v69OkWj+fkdL/1lZ1BiJeWCQODuTTSh3CNhSiNhUHuJkI1VuytormrFcQE6kgtqcPQyYuhKywKSi1KojQWPBUNuZogQB+tGUmGLKOaM2VdpVWkyKwiUmPF8zxFrzOwu/X69NNPA2CxWBr/f4bMzEz69OnjXGXdGIUo0DfAA7mq/YYZHOzJqfJ6kotquSTSx3nizqHEosBTYSNQ3VSrWoTeWjMZRg3Jeg0eColKiwIfpY1gVeeGGrPbdJWVlUBDYL4z/4eG6eqDBw/mzTffdL66PzAapciQYE+SCmro7eeO/3mmWXUEswS1NpForaXFz31VEn0xUW1ToLeJuCkkemvNdufY7cVu03322WcAxMfHM2fOnE4TdDExIEhHdqWB/TmV3DAgyCnTrc6mwqpAAHyV58+R/VQSfqquHXVxuE43YMAA0tLSmhxLTU1l06ZNThN1sSAKApdF+VJrspJc1Hy37I5SYVHirZRQdrPuQIdNN2fOHDw8mi7E9fDwcOV+7cTXTcXgYE9OFNdSqe9YIJ6zMUoC9ZKIXyu53IWiXfHpzg0fERkZSX5+vtNEXWwMDvbEU6vkUF41zuo2rbAoEGmYc9jdcNh04eHhzcZfk5KS2jV93UUDClEgPtyb0nozuVVGp6RZYW1oiTo6A6crcNh0jz76KFOnTmX58uXs37+f5cuXM3Xq1MaxWRftI9RLS5iXhsSCamxSx3I7oyRg6KZFK7RjlsmZuts///lPsrOziY6O5vHHHycuLs7p4i42RoR7sz6lhNSSOoaEtB0k6HxU/la0njvE1V1wOKcTBIEnnniCEydOkJ6ezvTp03nvvfe4+eb2B/Nz0YCXVsWAIB3Hi2oorjW1fcF5qLQq8GrH5NauwmHT2Ww2fvjhByZNmkR0dDQvvvgi06dPv2iGwTqbYaFeBOk07Mosp1Lfcqdua+jNNuolsdW+uQuN3aZLS0vjmWeeISwsjBkzZuDr68v69esJDAxk9uzZ+Pv7d6bOiwaFKHB1bz90GiU7TpVRb3ZsSCq/2oDQTVutZ7C7Tjdo0CD8/f3597//zV133dXYV9edJyL2VFQKkWv7+rMxrZSDudVc29f+H3RutRFPRffrED4bu3O6+++/H6PRyN/+9jfmz5/fGNvERefgplIQH+FNQY2R/Gr7NpkzWyVKak3dumgFB0y3ePFiioqKeOWVV9izZw9xcXHEx8dTW1tLdXV1Z2q8aInw1hLiqeFQXjVWO7pR8muMyNCti1ZwsCHh4eHBrFmz2Lt3L8eOHWPMmDF4eHgQFxfHrFmzOkvjRYsgCFwS4YPBYuNEcdtjszmVegI81HZNJL2QtFve4MGDWbhwIfn5+SxbtozCwkJn6nLxG55aJYOCPEkprqXGeP7WrMkqUVhjIsrXORHaO5MO/yaUSiVTpkxh3bp1ztDjogUGh3jirlZwILfqvGOzuVUN9b5In4vAdC46H6UocGmkL6V1ZjLL9S2ek1NpIMhTg5vKORGpOhOX6XoIIZ4aevu5c6Sgutm6CoPFRnFdzyhawWW6HsWIcC9EQeBwXtPegpxKA6IAkd4u07lwMhqlgvhwb3KqDORXN0yBkmWZ05V6Qr20bQab7C70DJUuGonydSPEU8PB3Cr0Fhu7siqo0Fvo6998W6Xuist0PQxBELg00geTVeLH5GJKak1c08efcG/nbZve2Th3ixsXXYJOo2REhDdZ5fVcHuWLlxPCd3UlLtP1UGICPIgJ6DlF6tm4ilcXXY7LdC66HJfpXHQ5LtP1UKxWG7PfWM3Wg6cutBSHcZmuh7J0/REWrz3E5L9+wc4jWRdajkO4TNcDqTeYmf+/zQBYbRK3zl3GvuM9Z2GUy3Q9kPe+2UtFTcNUJkmWMVtt3Pn8V9hs3XOd67m4TNfDKK2s5z9f7EI6a16dJMmUVtWzPzn3AiqzH5fpehivfLoFUwvLEpUKkXV70lq4ovvhMl0PIvFkIZ+sOdhirBOrTWLV9uQLoMpxXKbrIciyzBML1yKK5//KsgoqOZlT1oWq2ofLdD2ErzYk8euJvFYbC6IgsHZXSheqah8u0/UAauqNPPP+z20GoJZlmdU7TnSNqA7gMl0P4PWlO6isMdBWkE4ZOJiST1WdfREBLhQu03VzMnLLefebX5p0kbSGDBw80b1D8bpM182Z++56h85XKEQS0ws6SY1zcE3i7Mb8tDeNDfvSHbrGZpM4ktq9oy24crpuyomsEmbM/xaxHaHYavQmcoqrnC/KSbhM1w0pLKtl4tNLMZosdtflzkZAIOlkUScocw4u03UzbDaJW/+6jJLK+nZHWRcEOHKy+zYmXKbrZmw5eIpjGUUdmjEiyTJH0zuWRmfidNPl5eU5O8mLis/WHkLphN2uTRYbGXnlTlDkfJxuusGDBzs7yYuG0sp61u5OxeqEHEoUBJLSu2cr1ummS07uGTMduiNfbUxCkpxTJMqyzJG0i8R0kZGRzk7yokCWZT5Z/WubQ112pwekZpditnTuLtXtwaHO4c8//9yu8+677752ibmY+fVEHum5zq2DWW0SJ7JKGd4/1KnpdhSHTPfWW281+Ts5ORl3d3dCQ0MpLCxEr9cTGxvrMl07eGfFXhQK0aktToXYMCTWo0135MiRxv+/8MIL3HzzzcyfPx+lUonFYuHll192usCLgYMp+XzfCbN+bZLE4dQC7r85welpdwRBbueutoGBgRQWFqJU/u5bi8VCWFgYpaWlThMIkJGR4fA1/fr1c6qGlnBUV0uaZFlm/JzF/HI812m5XLTf7/8XgC9evgsvD02H03XWO213Q8LNzY3Dhw83OZaYmIhW23PipHUHNh/IYHdSdqd15MrAsYzuNSTW7lkmzzzzDOPHj+fee+8lKiqK7Oxsli9fziuvvOJMfX9oZFnm+Q82IIoCUgc3Fj4fZ+p1V8VFdUr67aHdOd3jjz/O999/j9FoZNu2bRiNRlauXMnjjz/uTH1/aHYnniY5q6TTDAe/1+taw2qT2J+cy5vLd3Eis7jTtJyhQ/Ppxo4dy9ixY52l5aLjsx8bhrycMQLRGiWV9RRX1BHsp2v2WWJaAW98sYtavQkBOHaqmEXP3YpG3XnRPV1jr+3AZpP467vruefFr7nrb18x/cWv+eWYY6vrq+uMfLctudMNd4Y9R7ObHTOYLLz51W7q9A27aMtARbWBrzYe7VQtrrFXB9Ebzdw57ys+WLmPWr0JvdFCjd7Egi92kl1UaXc632w+htnSdbsVfvlzImVV9U2OrdiYRHWdibMLdwmZVdtPkFVg/7M4imvs1QHKquoZP2cxP+872Wy4ymqTePmTrdQZzBSW1bAr8TTJmcXojeYW0/pk9a9tLil0JhabzIff7Wv8+3RhJT9sT2l5kqgg8O7XezutrtnuOt3ixYsZOXIksbGxJCUlMWPGDFQqFUuXLv1Djr+eLqhkwl+WkFNc1eKXIckyZVV6Zvz9m2ZFZkxkAE/fczVnurmOZhRxtIu7MSRJ4tcT+azfm4ZCEFizK6WhE68FX0mSREZeOT/uTuWWawY5XUu7O4f79OnD/v37CQwM5MYbb2TYsGHodDq2bdvGtm3bnCryQncOJ6UXMvEvS6msM2Cz/f66zu6EbQ1RFFCIIo9Nv5Gi8lr+98Ov1NQZ2z0zuC3Op+tsj4mC0OZUeJVSwUfP3trYAHHWO2236by8vKipqcFoNBIcHExJSQlKpZLAwEAqKiqcIu4MF9J0G/adZNoLKzBbbM1MYq/pznC6AhSi0GlmO4Ojus6HKAoM7RPMq4+OB5z3TttdvPr6+pKens7x48dJSEhAo9FgNBqdNh/sQmC2WFm9M4X0nDJC/D0pq6rnH59sRUZ22pSjzjacM5EkmaSMIh56bRXeHlr69+/Hvx+/qcPpttt0Tz75JMOHDwdg6dKlAOzevZtBg5xfB+hsLFYbry3ZzqJVB6isMXRJ31lPorCslsKyWrYdr7iwpnv66aeZNGkSCoWCPn36ANCrVy8+/vjjDovqSqxWG/f/YyWrdpxo3DXaZbjOpUMjEjExMU3+7t+/f4fEtIYsy3y4ch/ZRVXcfPVArhoWhVIhNhpFaEf/gyRJPPTaKlbtSHZa8emibXpMWIn/rtrPz/vSERFIOb2Lj7T7EAWR+t/6wXw9tYQGePHEXVcQGuDV5No6vYk5C9aSWVDB32ddx9hL+nIso4i/fbSRzQcyWuo1cNGJdHvTybLM5+sPs27PSaChxxyg3mhpcl55jYHKWiNvfL6TBU9OaDyeerqUO55fTmZBJQJw81+WMiAqgJPZZYgK0WW4C0C3XmxdXFHHHc9/ybdbjtt1viTLnMqv4Jstx4CGoZ8rZ31EVkElkiQ3thzTssuQodsuRv6j0y1zOlmWWbn1OE8sWEut3kSkjwPXAis2HOXwaQM//3KysyS66ABdbrrVO0+wOzGbW0YN5MphUSjOWc2enlPGEwvXsv1w1vlGadpGgE37HQux5aLr6HLTLf85ibW7Unj/21/w8dRy85UDGH95fyxWG19vOtqwwdpvLdH21rckWcZVcnZfLmjxWlVrZPmGJJZvSDrnE1f1/o9Ml5vu+pH9yMxvuqi4uKIOAYEgv5a3Bw9wa3l6UGt4+Kjbpc8RHNXVFZqgfe/LHi5P6OWUdNo94N8WkiRRUOCc2Lc5OY7v8Nerl3NeUGs4qqsrNEH73pc9nE9/WFhYq5uqnEun5XQFBQV/yHl1LpqTm5tLRESE3ef3iJzORffG0Zyu00znwsX56NYjEq3xzuZ0xr+1gxve2smk93aTW6G36zpBEHjsscca/y4sLEShUDB//nwAfvzxR2JjYxFFkePH7RsJAUgZEkvm5Ns5dfNEsm6fQuU33zj0PG3pWrhwIYMGDWLYsGFMnjyZmpqaVtOTJIlFixaxaNEiFixYwJtvvsmiRYv47LPP7NKTmJhIfX192ye2g245ItEWh7Ir2JdZzvonRqFUiBRWG3BX2fcofn5+7Nu3D5vNhkKhYOXKlQwZMqTx8wEDBrBy5UoeeeQRhzQpPD3ps+p7ACz5+eTOmQMy+N51p1N0JSQkMHv2bLRaLS+88AILFixoNWCRKIqNz7B9+3bc3d0ZOXKk3c+TmJhIWFgYHh4t9yh0hB5putJaE95uqsbYvKHebnZfKwgCo0aNYseOHYwdO5ZVq1Zx++23N35+7nSt9qAKDyf4mWcpfv11u03Xlq7Ro0c3/j8+Pp5169Y5rCs/P5+NGzdiNpvx9fXltttuQ61Ws3HjRk6ePIlCoWDIkCEEBgZSUFDAN998g1arZdasWQ7fqzV6pOmujgnkrU3pXLdwO6NiApk8Ipw4BwZo77zzTpYtW8bAgQNRq9UEBARQVubcfVK1QwZjzspy6Bp7dS1ZsoTp06c7lLYsy2zcuJG7774bNzc39uzZw4EDB4iPjyc5OZmnnnoKQRAwGo1otVr279/PhAkTCAoKcug+9tAj63Q6jZJ1T1zN/FuGoFGJ3PvpfnaetD882ZVXXsmBAwdYsWIFU6dO7RyR7Wif2aPr7bffBuCuu+5yOP3i4mKWLl3KokWLSExMpLq6Gq1Wi1arZfXq1aSmpqJWd34Hdo/M6aBhz/pRMYGMignEz13N5pRirukfaNe1giBwzTXX8Prrr5OSksJXX33ldH3GlFTUv03jt5e2dK1du5Zly5axY8eOdmkKDQ1l5syZzY7/+c9/JjMzk6NHj3L06FHuvNO+KkF76ZGmO1Vah1IUiPL3QJZlThbX0T+4eXCY1pg9ezZxcXH4+/s7XZ+lsJCSN97Ab/o9Dl97Pl2HDh1i7ty5bNmyBZ3OsWc9Q3V1NUVFRYSEhGA2m6mrq0On02GxWOjfvz8hISEsWbIEAI1Gg8lkatd92qJHmk5vsvHSmuPUGRsihw8N92bmldEOpRETE9Nio2HDhg08+OCDlJaWMm7cOMaMGWNXTmirrSXztsnIFguiRoPPtLvxaUfRfT5dzz77LDU1NUycOBGAq666ig8++MDudAVBYMqUKaxbtw6zuWFs9vrrr0elUrFixQpstoa4KuPGjQMgLi6O1atXd0pDwtU57KLL6ZENCRc9G5fpXHQ5LtO56HJcpnPR5bhM56LLcZnORZfjMp2LLsdlOhddjst0Lrocl+lcdDn/DyPrv9zYVe27AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 200x480 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=plt.subplots(3,1,figsize=(2,3*1.6),dpi=100)\n",
    "fig.tight_layout()\n",
    "for i in range(3):\n",
    "    print_sig_ax(ax[i],U[:,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inclusive-evanescence",
   "metadata": {},
   "source": [
    "variable V contains the loadings of the neurons onto each profile\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;(note that variable \"V\" is $V^T$ in the equation)\n",
    "\n",
    "you could show the scatter plot of loadings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "peaceful-hungary",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f868f928550>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD6CAYAAAC8sMwIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAXX0lEQVR4nO3df4xlZX3H8c93h0sZwDqsbFN3YF1KzdDqVjadgHabpljsEK06Uq1S0jRpE+IfpmjbiUs1BSwGmk21Sf9yE0ltQilYcIulzaoBYyXyY5ZdhGVZi6JlB1PWsqPijjK7++0fMwPn3jn33nPuec6P5973K9lk77kz5zz3zsz3PPf7PM/3MXcXACBeG+puAACgGAI5AESOQA4AkSOQA0DkCOQAEDkCOQBErnAgN7MzzOxhM3vMzA6a2Y0hGgYAyMaKziM3M5N0lru/aGYtSV+XdK27P9jte84991zfunVroesCwKjZt2/fD9x9U+fx04qe2FfuBC+uPmyt/ut5d9i6davm5+eLXhoARoqZfS/teJAcuZmNmdkBSc9L+rK7PxTivACA/oIEcnc/6e4XSzpP0iVm9sbOrzGza8xs3szmjx49GuKyAAAFnrXi7ouSvirpipTndrv7tLtPb9q0LsUDABhQiFkrm8xsYvX/45Iul/RU0fMCALIpPNgp6bWSPmdmY1q5Mdzp7v8e4LwAgAxCzFr5pqTtAdoCABhAiB45gAD27F/Qrr2H9dzikjZPjGtuZkqz2yfrbhYiQCAHGmDP/gVdd/fjWlo+KUlaWFzSdXc/LkkEc/RFrRWgAXbtPfxyEF+ztHxSu/YerqlFiAmBHGiA5xaXch0HkgjkQANsnhjPdRxIIpADDTA3M6Xx1ljbsfHWmOZmpmpqEWLCYCfQAGsDmsxawSAI5EBDzG6fJHBjIKRWACByBHIAiByBHAAiRyAHgMgRyAEgcgRyAIgcgRwAIkcgB4DIEcgBIHIEcgCIHIEcACJHIAeAyBHIASByBHIAiByBHAAiRyAHgMgRyAEgcgRyAIgcgRwAIkcgB4DIEcgBIHIEcgCIHIEcACJHIAeAyBHIASByBHIAiByBHAAiRyAHgMgRyAEgcoUDuZmdb2b3m9khMztoZteGaBgAIJvTApzjhKS/cPdHzexVkvaZ2Zfd/ckA5wYA9FG4R+7u33f3R1f//2NJhyRNFj0vACCboDlyM9sqabukh0KeFwDQXbBAbmZnS7pL0ofd/Ucpz19jZvNmNn/06NFQlwWAkRciRy4za2kliN/m7nenfY2775a0W5Kmp6c9xHWBMuzZv6Bdew/rucUlbZ4Y19zMlGa3ky1EcxUO5GZmkj4r6ZC7f6p4k4D67Nm/oOvuflxLyyclSQuLS7ru7scliWCOxgqRWtkh6Y8kvdXMDqz+e3uA8wKV27X38MtBfM3S8knt2nu4phYB/RXukbv71yVZgLYAtXtucSnXcaAJWNkJJGyeGM91HGgCAjmQMDczpfHWWNux8daY5mamamoR0F+QWSvAsFgb0GTWCmJCIAc6zG6fJHAjKqRWACByBHIAiByBHAAiRyAHgMgx2ImoURcFIJAjYnXWReEGgiYhtYJo1VUXZe0GsrC4JNcrN5A9+xdKvS7QDT1yRKuMuihZetq9biD0ylEHeuSIVui6KFl72hTWQtMQyBGt0HVRsqZqKKyFpiGQI1qz2yd185XbNDkxLpM0OTGum6/cNnB6I2tPm8JaaBpy5IhayLoomyfGtZASzDt72hTWQtMQyIFVczNTbdMZpe49bQproUkI5MAqetqIFYEcSKCnjRgx2AkAkSOQA0DkSK0AFaA2C8pEIAdKVmdxL4wGAjnQIXTvmdosKBuBHEgoo/dMbRaUjcFOIKGM0rjUZkHZCORAQtoS/V7Hs6A2C8pGagVIGDPTSffU44NixSjKRiAHEtKCeK/jWbFiFGUitQIkTHbJW3c7DjQBgRxIuOyiTbmOA01AIAcS7n/qaK7jQBOQIwcSss75Zsk9moQeOZCQZc531k2agaoQyIGEuZkptcbapxq2xqxtzncZi4b27F/Qjlvu0wU779WOW+7jpoBcCORAp86Zhh2PQy+5T+vhf+SOA/r4nscHOh9GD4G8ZvTEmmXX3sNaPtUeuZdPeVtvO/SS+7Qevku67cH/4fcBmQQJ5GZ2q5k9b2ZPhDjfqCDX2jxZetuhl9x3u6ZLhdI1GB2heuT/KOmKQOcaGWXkWlFMlt727PZJ3XzlNk1OjMu0sljo5iu3vTxrJe+nrF49eSokIosg0w/d/WtmtjXEuUYJ5U2bZ25mSnOff6wtvdLaYOt6292W3A9SBnduZkofuePAutS8RIVEZEOOvEaUNw0j+DhDZ32sHPWyBvmUNbt9Ule/ecu6y1AhEVlVFsjN7Bozmzez+aNHWSUnUd40hNDjDLv2HtbyyY7BzpOeOd016Kesm2a36dPvv7hrugbopbKVne6+W9JuSZqeni5WSm5IUN60uNDbqBVNd22eGE+tXZ7lUxYVEjEolujXjD/eYkKPM0yc2dKx48upx7OYm5lqy5FLfMpC+UJNP7xd0jckTZnZETP70xDnBfoJPc7w047efb/jnfrNaAHKEGrWylUhzgPkFboHvLR8KtfxNHzKQtVIrSBqjDMABHI0VJ4ysfSAMeoI5GicQRbVhGJaXzNr7TjQVCwIQuPUWbrg6jdvyXUcaAICORqnztIFN81u044LN7Yd23HhRt00u630awODIpCjceosXbBn/4Ie/u6xtmMPf/cYFSnRaARyNE6dpQtu/OLB1CX6f37ngcbVjKeWPdYw2Il16t5YuM4phWmrOiVprRhiqIHXou9xnQPCaB4COdo0JUA0eUphkVouUpj3OHSNGcSN1ArasNlFNkUGXkO8x9SyRxKBHG0IENkUGXgN8R5Tyx5JBHK0IUD0V3TgNcR7TC17JBHI0abuAFHHTIzkNXsJVc0wxHtMlUUkMdiJNnXOGKljoLXzmr08c8s7glwz1Hvc5AFhVItAjnXqChB1zMRIu2aaczJuLJEVQRghEcjRGFUOtK7N407blq1Ta8x0/TvfELwNQCgEcjRGkf0u88iTTpGkS7aes673XMWiqboXZiEeDHaiMaoaaM2aTlnzwLdf0Mf3PP7y47UbwcLiklyv5PJDDsxWcQ0MDwI5GqOqmRiDpGpuf+jZl/9fxaIpFmYhD1IraJQqBgEnzmyl1lSZ7JLakaST/kohrSpy+SzMQh70yDFS9uxf0Is/PbHueGvMNDczpTFL3wsoebyKRVMszEIeBHKMlF17D2v51PrN3M46/TTNbp/Um3/pnNTvSx6vIpdf98IsxIXUCkZKt9TED5dWUi1Pfv/Hqc8njxdZ0JN1JkqdC7MQHwI5Rkq/KY7d6pF3Hh8kl5935SqLhpAVqRWMlCIpi6K1X5iJgrLQIx9SsS4mydPuQV5jkZRF0dovzERBWQjkQ6js4lNl3STytLvf1/ZqY5GURZHaL1WtXMXoIZAPoTKLT6UF0I/ccUAfvuOAJgfce3It4G4wa5uv3avd/dIUg9zIsqZNBu1Bz81MrSsNwEwUhECOfAiV+RE+LYCuhd68y8g7l6F3BvE1ae3u9RoHzUVnzVUP2oOmhjjKQo98CJX5Eb7fzSBPzz9rzZO0dnd7ja8ebw18I8tyoyvag2YmCspAj3wIlbmYJMvNIGvPv0jgnJuZUmvD+lWYP3nphF49nl47vF/b+z0/Md6iB41GIpAPoTI/wqfdJDpl7fl3+7pkeD6jlf4rOrt9Umefsf4D5fJJl5kGupH1e/5nJ071fB6oi3mXvGSZpqenfX5+vvLrIozkpgymV3Lk0krAzHrTSKsL3tpgkq0E5LRzJgdHu/3mmqRPv//igWbWbO2zb6ekgQZ1gRDMbJ+7T3ceJ0eO3JJ53iJTEdPmdB9/6cS6VZTdZqN0s3livNRcdBV7iQJ50CNHo1yw897Unrap+wBnUp5PBGmy9MjXTE6M64Gdbx3oOsAguvXIyZGjUXqVb+01OFrHdD5WZKIpCORolF4zbroF+cmJcT1zyzv0wM63VprqmDgzfXYMULUggdzMrjCzw2b2tJntDHFOxGvP/gXtuOU+XbDz3tyFpnrNuEkL8ibpsos2hX0BGR07vqztn/gS+2iidoVz5GY2Julbkt4m6YikRyRd5e5PdvueQXPk/QbWRmVn8z37F3TDPQe1uFpD+5wzW7r+nW/ItdN7kdexZ/+CbvziwbZBSTPp6ku3aPp1GzX3r4+1zTqRpNYGaTkxe69ztkvSBpP+8NItuml227rX2o1J+uVfOEv//fxP2o79xoUb9Y1vv6CmTRycGG/JbOVmYCZ1+zOcnBjXZRdt0v1PHR3o59j5fK9zoTyh4ka3HHmIQP4WSTe4+8zq4+skyd1v7vY9gwTytKlqndPSej0fQhXXyNKGuc8/tm6Xm9aYadd739T2B96trdL62R9ZX8ee/QupgXrN6WOml7o8l9eOCzfq4WeOpe7oM8qy/hzTfge6nYtgXp6QcaPMwc5JSc8mHh9ZPRZUv/oZZdR67kwR3HDPwdrrSXfbqmz5pLe1o9f7UeS92rX3cNcgLilYEJekB779AkE8RdafY5YSCNRDL18VdehDzCNP26123V+fmV0j6RpJ2rJlS+6L9KufEbpQVFqVv7xtK0OvayWfG+T9yPI6mKnRDFl+jiFLJWBwVdShD9EjPyLp/MTj8yQ91/lF7r7b3afdfXrTpvyDU/12FQ+963jWgk5FrjGIXtdKPtfr/SjyXlE7uxmy/ByLlkpAGKFjU5oQgfwRSa83swvM7HRJH5B0T4DztulXCCp0oaisd8uq60l3KxbVGrO2dvR6P4q8V3MzU2qNpX0IW7Hjwo3B5rTuuHBj6msddVl/jlnq4lAPvXxlFrFbUzi14u4nzOxDkvZKGpN0q7sfLNyyDv226Aq963i3VYTnnNnSmaefVtuo/9q1+s1ayfJ+DPJerX1Nt1kraTNNzjp9TKfctZSYtjLorJVe3xeTULNWpHx/E8xaqV7o2JSGJfpdNGGGCnpLm0GTnL2z45b7+i7pT5rMUAKg3/dnXbIfchprE6bEohos0c9pVHZzKbJ4p26z2ye1671vavsZJadgdltAlCZLEB4zk2mlN92ZXsrzUblzZ6S8OyuVdS7Ei+qHPQz7bi5lb9JchV4/o26phbv2LQy0b+bf/UH7PP1Be8Eh91Qtc39WxINAPsJGIQikBfrp120sND7Q7bxZhZyOVsXUNjQfgXyEjWoQqPuTVsg9VcvcnxXxIEc+wqqY3xqTXnnlHjMuM503OQ5x2UWbgk1Hq2JqG5qPQD7CCALtei2ZHrTyQNpg5F37FvT7vz4ZZCB9VAbl0RuplRFWxfzWKhWdhldGSqnbOMT9Tx0NtrtQ3aki1I9APuKGJQiEmIGTZSu5tOv2unmM6jgEqkVqBUMhRIW5vCmlLHO4GYdAFQjkGAoher4h53CvYRwCVSC1gqFQ5TS8tXRKtzRM8uYxbOMQaCZ65ChVVSUAqur5JtMp3SRvHtRBQRXokaM0VZYACNXz7dxbNHl87fz9tk5bu3kMQwkExIFAjtJUXQIgxAycs89otZXnTR6XeufcJztuHqNQAgHNQCBHaWKcereYEsSTx7vl4tOqJ8b4+hEncuQoTYxT7/q1OU8ufuLMVuq5uh0HBkWPHOv0GqDLM3g3NzOVujlH0QHIMgcQt74mvce99TUrgTxPLr7bni017OWCIUcgR5teA3SScg3elTH1ruwBxAe/c6zv8ay5+B8upadpuh0HBkUgR5t+i1zyDt6FLgFQ9gDiyS7d5W7He6HELKpCjhxteg3QNWHwruw2jFl6vdpux3thVSeqQiBHm16DfU0YvCy7DVdden6u473UVWI25n1YMRhSK2jTb4CyjMHLkO0r6qbZbZKk2x96VifdNWamqy49/+XjeVVdXZJFSKOJQI42WQYo61xyHnoANW0GzE2z2wYO3HVjEdJoMq9hLtT09LTPz89Xfl0gqbP3Kq307mPeYeeCnfcq7S/aJD1zyzuqbg4CM7N97j7deZweOUZWt97rjV88GG2hK2bKjCYGOzGyus10OXZ8uedmEU1WxUwZBlObh0COkZW1l5p3p6E6lT1TJsuuSKgeqRWMrLQZMN3EVOiqzJkyDKY2E4EcIyttBsxPfnZCiylL6Mkxr2jCojCsRyDHSOvsvXabycJqzBUMpjYTOXIgoa7VmLGg7EAz0SMHOlS9GjMmbCbdTARyALlwo2seUisAEDkCOQBEjkAOAJEjkANA5AoFcjN7n5kdNLNTZrauIhcAoHxFZ608IelKSZ8J0BZgpKTVQmc2CAZRKJC7+yFJsgH2MwRGGTv5ICTmkSOYOnuYZV879PkpPoWQ+gZyM/uKpF9Meepj7v5vWS9kZtdIukaStmzZkrmBiEOdPcyyr13G+Sk+hZD6BnJ3vzzEhdx9t6Td0spWbyHOieaos4dZ9rW7nf+GewbfSYjiUwiJ6YcIos4eZtnX7naexaXBdxKi+BRCKjr98D1mdkTSWyTda2Z7wzQLsenWk6yih1n2tcvYSYgqiwip6KyVL0j6QqC2IGJpu+1U1cMs+9pl7SRE8SmEwqwVBFFnedOyr512/uMvndCx4+wkhGYw9+rHHaenp31+fr7y6wKhdNtJiPQIymRm+9x93Sp6euTAANhgAU1CIAcGRI4bTcH0QwCIHIEcACJHIAeAyBHIASByDHYCHagTjtgQyIEE6oQjRqRWgIRelRSBpiKQAwnUCUeMSK0ACSHrhJNrR1XokQMJoeqEr+XaB61XDuRBIAcSQtUJJ9eOKpFaATqEqKFCrh1VokcOlKDOHZMwegjkQAnYkxNVIrUClIB65agSgRwoCfXKURVSKwAQOQI5AESOQA4AkSOQA0DkCOQAEDkCOQBEzty9+ouaHZX0vYoud66kH1R0rTLE3P6Y2y7R/rrR/vVe5+6bOg/WEsirZGbz7j5ddzsGFXP7Y267RPvrRvuzI7UCAJEjkANA5EYhkO+uuwEFxdz+mNsu0f660f6Mhj5HDgDDbhR65AAw1IY+kJvZ35jZN83sgJl9ycw2192mPMxsl5k9tfoavmBmE3W3KQ8ze5+ZHTSzU2YWzQwEM7vCzA6b2dNmtrPu9uRhZrea2fNm9kTdbcnLzM43s/vN7NDq7821dbcpDzM7w8weNrPHVtt/YyXXHfbUipn9vLv/aPX/fybpV939gzU3KzMz+11J97n7CTP7W0ly94/W3KzMzOxXJJ2S9BlJf+nu8zU3qS8zG5P0LUlvk3RE0iOSrnL3J2ttWEZm9luSXpT0T+7+xrrbk4eZvVbSa939UTN7laR9kmYjeu9N0lnu/qKZtSR9XdK17v5gmdcd+h75WhBfdZakqO5c7v4ldz+x+vBBSefV2Z683P2Qu8e24/Alkp529++4+0uS/kXSu2tuU2bu/jVJL9TdjkG4+/fd/dHV//9Y0iFJ0RR19xUvrj5srf4rPeYMfSCXJDP7pJk9K+lqSX9dd3sK+BNJ/1l3I0bApKRnE4+PKKJgMizMbKuk7ZIeqrcl+ZjZmJkdkPS8pC+7e+ntH4pAbmZfMbMnUv69W5Lc/WPufr6k2yR9qN7Wrtev/atf8zFJJ7TyGholS/sjYynHovokFzszO1vSXZI+3PGpuvHc/aS7X6yVT8+XmFnp6a2h2OrN3S/P+KX/LOleSdeX2Jzc+rXfzP5Y0u9J+h1v4KBGjvc/FkcknZ94fJ6k52pqy8hZzS3fJek2d7+77vYMyt0Xzeyrkq6QVOrA81D0yHsxs9cnHr5L0lN1tWUQZnaFpI9Kepe7H6+7PSPiEUmvN7MLzOx0SR+QdE/NbRoJq4OFn5V0yN0/VXd78jKzTWszy8xsXNLlqiDmjMKslbskTWll5sT3JH3Q3RfqbVV2Zva0pJ+T9H+rhx6MbNbNeyT9g6RNkhYlHXD3mXpb1Z+ZvV3S30sak3Sru3+y5iZlZma3S/ptrVTf+19J17v7Z2ttVEZm9puS/kvS41r5m5Wkv3L3/6ivVdmZ2a9J+pxWfm82SLrT3T9R+nWHPZADwLAb+tQKAAw7AjkARI5ADgCRI5ADQOQI5AAQOQI5AESOQA4AkSOQA0Dk/h8VZeLq8oGGEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(V[0,:],V[1,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "multiple-cooling",
   "metadata": {},
   "source": [
    "or calculate the kurtosis of loadings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "younger-culture",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(2, 80, 'kurtosis=17.14')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAD4CAYAAACuX5+7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdl0lEQVR4nO3de7hddX3n8feHBCgIFizRhgAGLWJBK4wRqxYbLgLx0qCWGkCKYouM0IJlngHsiFCf9FEr2NFRBxyU1AsxiIwg5X6phqvhInItUQJEKASRBmiFIfnOH3udZBNOznWd7HOS9+t59rPXWntdvvsYPJ+zfpeVqkKSJKlNG/W6AEmStP4xYEiSpNYZMCRJUusMGJIkqXUGDEmS1LrJvS4AYJtttqnp06f3ugxJmlBuvvnmx6tqSq/rkPozLgLG9OnTWbRoUa/LkKQJJckDva5BWhubSCRJUusMGJIkqXUGDEmS1DoDhiRJap0BQ5IktW7QgJHkt5LclOSnSe5Mcmqz/ZQkv0xyW/N6Z9cxJyVZnOTeJPuP5ReQJEnjz1CGqT4L7F1VTyfZGFiY5OLmsy9U1ee7d06yCzAH2BXYFrgiyWuqakWbhUuSpPFr0DsY1fF0s7px8xroGe+zgflV9WxV3Q8sBvYYdaWSJGnCGFIfjCSTktwGPAZcXlU3Nh8dk+T2JF9PsnWzbRrwUNfhS5tta57zyCSLkixatmzZKL6CJEkab4Y0k2fTvLFbkq2A85O8Dvgq8Gk6dzM+DZwGHAGkv1P0c84zgTMBZsyYMdAdEamnpp94UU+uu+Qz7+rJdSWpDcMaRVJVTwLXAAdU1aNVtaKqVgJfY3UzyFJg+67DtgMebqFWSZI0QQxlFMmU5s4FSTYD9gXuSTK1a7f3Anc0yxcAc5JsmmRHYCfgpnbLliRJ49lQmkimAvOSTKITSBZU1Q+TfDPJbnSaP5YAHwWoqjuTLADuAp4HjnYEiSRJG5ZBA0ZV3Q7s3s/2wwY4Zi4wd3SlSZKkicqZPCVJUusMGJIkqXUGDEmS1DoDhiRJap0BQ5Iktc6AIUmSWmfAkCRJrTNgSJKk1hkwJElS6wwYkiSpdQYMSZLUOgOGJElqnQFDkiS1zoAhSZJaZ8CQJEmtM2BIkqTWGTAkSVLrDBiSJKl1BgxJktQ6A4YkSWqdAUOSJLVu0ICR5LeS3JTkp0nuTHJqs/1lSS5Pcl/zvnXXMSclWZzk3iT7j+UXkCRJ489Q7mA8C+xdVW8AdgMOSPKHwInAlVW1E3Bls06SXYA5wK7AAcBXkkwai+IlSdL4NGjAqI6nm9WNm1cBs4F5zfZ5wIHN8mxgflU9W1X3A4uBPVqtWpIkjWtD6oORZFKS24DHgMur6kbgFVX1CEDz/vJm92nAQ12HL222rXnOI5MsSrJo2bJlo/kOkiRpnBlSwKiqFVW1G7AdsEeS1w2we/o7RT/nPLOqZlTVjClTpgytWkmSNCEMaxRJVT0JXEOnb8WjSaYCNO+PNbstBbbvOmw74OFRVypJkiaMoYwimZJkq2Z5M2Bf4B7gAuDwZrfDgR80yxcAc5JsmmRHYCfgprYLlyRJ49fkIewzFZjXjATZCFhQVT9Mcj2wIMlHgAeBgwCq6s4kC4C7gOeBo6tqxdiUL0mSxqNBA0ZV3Q7s3s/2XwH7rOWYucDcUVcnSZImJGfylCRJrTNgSJKk1hkwJElS6wwYkiSpdQYMSZLUOgOGJElqnQFDkiS1zoAhSZJaZ8CQJEmtM2BIkqTWGTAkSVLrDBiSJKl1BgxJktQ6A4YkSWqdAUOSJLXOgCFJklpnwJAkSa0zYEiSpNYZMCRJUusMGJIkqXWDBowk2ye5OsndSe5Mcmyz/ZQkv0xyW/N6Z9cxJyVZnOTeJPuP5ReQJEnjz+Qh7PM8cHxV3ZJkS+DmJJc3n32hqj7fvXOSXYA5wK7AtsAVSV5TVSvaLFySJI1fg97BqKpHquqWZvkp4G5g2gCHzAbmV9WzVXU/sBjYo41iJUnSxDCsPhhJpgO7Azc2m45JcnuSryfZutk2DXio67Cl9BNIkhyZZFGSRcuWLRt24ZIkafwacsBIsgVwHnBcVS0Hvgq8GtgNeAQ4rW/Xfg6vF22oOrOqZlTVjClTpgy7cEmSNH4NKWAk2ZhOuPh2VX0foKoeraoVVbUS+Bqrm0GWAtt3Hb4d8HB7JUuSpPFuKKNIApwF3F1Vp3dtn9q123uBO5rlC4A5STZNsiOwE3BTeyVLkqTxbiijSN4GHAb8LMltzbZPAAcn2Y1O88cS4KMAVXVnkgXAXXRGoBztCBJJkjYsgwaMqlpI//0q/nmAY+YCc0dRlyRJmsCcyVOSJLXOgCFJklpnwJAkSa0zYEiSpNYZMCRJUusMGJIkqXUGDEmS1DoDhiRJap0BQ5Iktc6AIUmSWmfAkCRJrTNgSJKk1hkwJElS6wwYkiSpdQYMSZLUOgOGJElqnQFDkiS1zoAhSZJaZ8CQJEmtM2BIkqTWGTAkSVLrBg0YSbZPcnWSu5PcmeTYZvvLklye5L7mfeuuY05KsjjJvUn2H8svIEmSxp+h3MF4Hji+qn4f+EPg6CS7ACcCV1bVTsCVzTrNZ3OAXYEDgK8kmTQWxUuSpPFp0IBRVY9U1S3N8lPA3cA0YDYwr9ltHnBgszwbmF9Vz1bV/cBiYI+2C5ckSePXsPpgJJkO7A7cCLyiqh6BTggBXt7sNg14qOuwpc22Nc91ZJJFSRYtW7Zs+JVLkqRxa8gBI8kWwHnAcVW1fKBd+9lWL9pQdWZVzaiqGVOmTBlqGZIkaQIYUsBIsjGdcPHtqvp+s/nRJFObz6cCjzXblwLbdx2+HfBwO+VKkqSJYCijSAKcBdxdVad3fXQBcHizfDjwg67tc5JsmmRHYCfgpvZKliRJ493kIezzNuAw4GdJbmu2fQL4DLAgyUeAB4GDAKrqziQLgLvojEA5uqpWtF65JEkatwYNGFW1kP77VQDss5Zj5gJzR1GXJEmawJzJU5Iktc6AIUmSWmfAkCRJrTNgSJKk1hkwJElS6wwYkiSpdQYMSZLUOgOGJElqnQFDkiS1zoAhSZJaZ8CQJEmtM2BIkqTWGTAkSVLrDBiSJKl1BgxJktQ6A4YkSWqdAUOSJLXOgCFJklpnwJAkSa0zYEiSpNYZMCRJUusGDRhJvp7ksSR3dG07Jckvk9zWvN7Z9dlJSRYnuTfJ/mNVuCRJGr+GcgfjbOCAfrZ/oap2a17/DJBkF2AOsGtzzFeSTGqrWEmSNDEMGjCq6kfAE0M832xgflU9W1X3A4uBPUZRnyRJmoBG0wfjmCS3N00oWzfbpgEPde2ztNn2IkmOTLIoyaJly5aNogxJkjTejDRgfBV4NbAb8AhwWrM9/exb/Z2gqs6sqhlVNWPKlCkjLEOSJI1HIwoYVfVoVa2oqpXA11jdDLIU2L5r1+2Ah0dXoiRJmmhGFDCSTO1afS/QN8LkAmBOkk2T7AjsBNw0uhIlSdJEM3mwHZKcA8wEtkmyFPgUMDPJbnSaP5YAHwWoqjuTLADuAp4Hjq6qFWNTuiRJGq8GDRhVdXA/m88aYP+5wNzRFCVJkiY2Z/KUJEmtM2BIkqTWGTAkSVLrDBiSJKl1BgxJktQ6A4YkSWqdAUOSJLXOgCFJklpnwJAkSa0zYEiSpNYZMCRJUusMGJIkqXUGDEmS1DoDhqQJ55prriEJSTjllFN6Xc4GI8mHklTz+lCv69H4Nujj2iVJ7VqyZAlnn302ADNnzmTmzJk9rWdDk2QT4PXAG4EZzfvrgY2bXT5cVWcP4TynAJ8aQQnzqupDIziu77qbA29gde0zgNcCk5pd9qqqa0Z6/uYamwC3ALt2bR7WeQ0YkrSOLVmyhFNPPXXVugFjnbse+C89vP4vRnn8Q8DL2ihkACfxwnAxbAYMSdKQNH/Vn93jMtowaY31fwOeBV45zPPMB24bwn6/zeqfWwHzhnmdNa1Z/4PAJsDvjvK8ACTZBfhEs/oM8JKRnMeAIUna0FwF/AC4Gbi5qn45kuaOqroHuGew/ZIc1bV6ZVU9MJzr9OMHwL8Ci+jU/3iSs4HDR3lekmwEnEUnsFwIvBT445Gcy4AhSdqgVNXfrONLHtG1/I3RnqyqRh0kBnAM8Id07lwcA/zTSE/kKBJJ662f//znvPrVryYJG220EV/4whdWfdY3CmUo/R8G23fmzJmr9gFYsWIFZ599Nu94xzuYNm0akydPZvr06atGv+y1116rjj311FNXHdv9WpvLLruMww47jFe96lUAuyd5Ksk9Sf53kjcO5eeSZHaSBUl+keQ/kvwmyS+T/DTJN5McmmTrfo4b0iiSJDOTzEtyb5KnkzyX5N+S3JHke0n+Ikkrt/PHuyS7Am9qVp8Evt/DcgaU5JXA3Gb1k1X14GjON+gdjCRfB94NPFZVr2u2vQz4LjAdWAL8WVX9uvnsJOAjwArgr6vq0tEUKEkjceuttzJr1iweffRRJk+ezDe+8Q0++MEPjvl1n3jiCWbPns3ChQtbPe/TTz/NIYccwoUXXti9eSNgC2Dn5nVkki8BH6+qlWueI8lmwLnAu/q5xLbN6w+ADwIfB/5xODU2t9fPAP6in49f0bx2Bd4P/E/guOGcf4LqvntxTlX9pmeVDO4MOv+ebgG+ONqTDaWJ5Gzgf/HC2yQn0mlH+kySE5v1E5qOIXPo/APaFrgiyWuqasVoC5Wkobr66qs58MADWb58OZtvvjnnnnsu73znO9fJtT/4wQ+ycOFC3vCGN3DwwQczffp0li9fzu23387rXvc6zj//fO644w4++clPAvCBD3yAOXPmDHjOFStWMGvWrFWhZauttuKII47g9NNPvx84Ffgj4M/ptJv/NbAZcGQ/p/p7VoeLR4BvAXcCT9P5xfJ7wFuAt4/w6/8Vq8PFk835bwX+Hdiczh+lbwb26u/g9U2SyXTCWp+v96qWwSQ5DNifzs2BI9v4vT1owKiqHyWZvsbm2cDMZnkecA1wQrN9flU9C9yfZDGwB50hQZI05s477zwOPfRQnn32Wbbeemsuuugi3vKWt6yz61988cUcd9xxnHbaaWy00YtboQ888EC22mqrVeuvfe1rOfDAAwc85+c///lV4WLnnXfmqquuYtttt+X0009/oqrmAfOSfAW4gs7wxb9MckFV/bDvHEkmAR9uVh8A9qiqx/q7XpIpwO8M/Vuv8pfN+78Db66qf13L+V8KvGotnw38wxi6x6uq3dtIw/du4OXN8s+qalEvi1mb5n/vvvbDL1bVzW2cd6SdPF9RVY8AVNUjSfp+gNOAG7r2W9pse5EkR9Ik7B122GGEZUjSameccQYf+9jHWLlyJdOmTePSSy9l111HNZR/2N74xjeuNVyMxHPPPbeq78jkyZM599xz2XbbbV+0X1XdmuSjdJpAoHNn+Yddu0yhM1wS4AdrCxfNuZYBy0ZQ7u817/+ytnDRnH85ax/eef4Irtuff2H1H8K98uGu5XF794JOc8jv0Jlf4+S2Ttp2J8/+eiZVfztW1ZlVNaOqZkyZMqXlMiRtaD796U9z1FFHsXLlSnbeeWeuu+66dR4uAI4++ujWwgXAddddx6OPPgrArFmzeP3rX7/Wfavqe8DiZvVtXX/8Afxn1/IurRX4Qv/RvO+UZOMB91zPJXkF0Ncu9xyd5qJxJ8m76XRtADi6qp5u69wj/a/g0SRTAZr3viS8FNi+a7/tgIdHXp4kDWzlypUcc8wxnHxy5w+vN73pTSxcuLBnd0b33HPPVs930003rVreb7/9hnLI5V3Lb+5bqKp/B/pOtm+S85Ps3XIQ6Lv279Ppg/eepmPpkFVVWnrNbPF7jcRhrG4luLCqHu9lMf1JsiXw1Wb1vKq6cKD9h2ukAeMCVk/ocTidST/6ts9JsmmSHYGdWP0PWpJa98UvfpEvf/nLAOy7775cddVVbLPNNj2rZ9q0fluFR+yRRx5Ztfya17xmKId0N01MXeOzo4HlzfKBwJXAr5NcmeTUJHs2I0FG6gRW/1H5djq/E36dZGGSzybZfwO6szERmkc+S+dGwHI6nYNbNeg/pCTn0OmkuXOSpUk+AnwGeEeS+4B3NOtU1Z3AAuAu4BI6t1scQSJpzDz//POrlp955hmq+m2VXWc222xYf7AP6qmnnlq1/JKXDGnG5u5b3Ft2f9B0MtyNzqjAviaTlwB702l7/xHw8yQjGs9bVUuA3emMPHyy2bwp8Dbgv9P5vbA0yXEZaLKPCS7Jm1ndDPUwMO6ma0iyJ9A3w+hJVdV6a8NQRpEcvJaP9lnL/nNZPVGHJI2pY489lrvvvpvzzz+f66+/nv33359LL72ULbfccvCDh2DlyhdNJ7FOdX+PZ555ZiiHbNG1/NSaH1bV/cDhTYfQtwBvpTPM9Y/pDG+dDnwzyQ5V9ffDrbfpPPpXSf6GzlM++86/N51pp19OZ8TCLvQzlHY9GUXSPffFvHH6h/YRdPpN/iewTZL/sZb9up/PcliSP2qWFwzUkRecKlzSBLfxxhvz3e9+lzlz5vD9739/Vci45JJLeOlLX7rW4zbZZBOee+45nnvuuQHP//jjvW06nzp1dSvHfffdN5R+GDt1La/1r9JmwqermxdJtgA+RueOdICTk5xRVb8aSd1V9f/o3P2+Hjitefz3oXTa/DelM5T2S1X1szUOndCjSJo+Jx/o2jTqqcHHSN8dpM3ozKUyFN3B6Q5e2Bz3Ik4VLmnC6wsZ73//+wG4/vrrOeCAA1i+fPlaj+mbi+Lhhwe+M3zjjTe2V2ije5TJYE06e+yxx6rlyy+/fIA9V3lH1/KQ+8BV1dNV9TngvGbTpqye4nrUquq5qvoG8KWuzW9r6/zjyPtZPRz4x1V1Xy+L6SUDhqT1wuTJk5k/fz5/+qd/CrDqTsbaQsYuu3SayB944AF+8YtfrPW8X/ziqGdMfpEttljdijFYs8db3/pWfvd3O4/tuOiii7jrrrvWum+S97H6DsbCgea6GMCSruWxuMs94PnXg1EkE6FzJ1X1oaH8HOncCeqzV9dn/3ewaxgwJK03Jk+ezDnnnLMqZNxwww1rDRkHHHDAquUTTjih3zsJJ598MldccUXrde64446rlm+55ZYB991kk034+Mc/DnQ6tB500EEvGFnSJ8kf0HmWRJ/PrPH57kk+2czP0K8k2wAHNasF3D5gcS88dmqSzzcjCNe2z+a88JHiPx3q+SeCZtbrvmnQn2L1pGdDPf7srofJndJudeuefTAkrVf6QkYSzj333FUh49JLL31Bn4wjjjiCz33uczzxxBN873vfY8899+TQQw9lm2224cEHH2T+/PksWrSIOXPmMH/+/FZr3Hrrrdl999259dZbufrqqznqqKPYZ599XtChszsAHX/88Vx44YUsXLiQu+66i1133ZUjjjgC4GVJ/pxOU8PhdJo1AL5WVRetcdnfBv4O+FSSa4Hr6LShP0VnevHXA4c0ywDfHubTNDcFjgeOT/IT4MfA3XRGk/w2nYexHcLq2Z1/DPSkE2aS3ek0ZXTrfv7K+5L83hqfn9V0kB3Ih1jdt2FBVQ2pV+5wJdmbTqfZbrt3LX8kyb5rfP75qnqSdciAIWm9M3nyZL7zne+QhAULFnDDDTew3377cdlll60KGVOmTOFb3/oW73vf+/jNb37Dtddey7XXXvuC87z73e/mrLPOaj1gAMydO5f3vOc9rFixgjPOOIMzzjjjBZ9331GZNGkSF198MQcffDA//OEP+fWvf81pp50GsCOd50GtOgz4MnBsP5fsGw4zic4v04EeaLaA/h+WNpDuW0BvYuD+G1cDB1XvxhS/AfjbAT5/T/PqdgWw1oDRDLvtvjszls0jb2fg+vsbZvx/WD10eJ2wiUTSeqkvZHzgA50O/TfeeCP77bffC5pLZs2axW233caHP/xhdthhBzbZZBOmTJnCXnvtxTe/+U0uuOACNt988zGpb9asWVx77bUccsgh7LjjjoPOn7HFFltw4YUXcskll3DIIYfwyle+Ejq/1J+hcyfiTOBNVfVX/T2qvap+ROcuwn+lEyDuoTNnxsrm/S46vxT/uKo+UFX/ueY5BlJVDwA70OmD8E/Az+g89GwFnSnE7wPOAd5TVXuPdHTKOLY3nSG+APdW1XU9rGVcSK8npQGYMWNGLVo0Lh8yJzH9xDXvNK8bSz7zrsF30gYtyc1VNaPXdUj98Q6GJElqnQFDkiS1zoAhSZJaZ8CQJEmtM2BIkqTWGTAkSVLrDBiSJKl1BgxJktQ6A4YkSWqdAUOSJLXOgCFJklrn01Q1IfTqeSCSpJHxDoYkSWrdqO5gJFkCPEXncbzPV9WMJC8DvkvnsbVLgD+rql+PrkxJkjSRtHEHY6+q2q3rkcEnAldW1U7Alc26JEnagIxFE8lsYF6zPA84cAyuIUmSxrHRBowCLktyc5Ijm22vqKpHAJr3l/d3YJIjkyxKsmjZsmWjLEOSJI0nox1F8raqejjJy4HLk9wz1AOr6kzgTIAZM2bUKOuQJEnjyKjuYFTVw837Y8D5wB7Ao0mmAjTvj422SEmSNLGMOGAkeUmSLfuWgf2AO4ALgMOb3Q4HfjDaIiVJ0sQymiaSVwDnJ+k7z3eq6pIkPwEWJPkI8CBw0OjLlCRJE8mIA0ZV/QJ4Qz/bfwXsM5qiJEnSxOZMnpIkqXUGDEmS1DoDhiRJap0BQ5Iktc6AIUmSWmfAkCRJrTNgSJKk1o32WSTawEw/8aJelyBJmgC8gyFJklpnwJAkSa0zYEiSpNYZMCRJUusMGJIkqXUGDEmS1DqHqUrjVC+HBC/5zLt6dm1J6wfvYEiSpNYZMCRJUusMGJIkqXUGDEmS1Do7eU5APg9EkjTeeQdDkiS1bswCRpIDktybZHGSE8fqOpIkafwZkyaSJJOALwPvAJYCP0lyQVXdNRbX61WTgXMFSJLUv7Hqg7EHsLiqfgGQZD4wGxiTgNEr9oXQ+mpDC+1Oaia1b6wCxjTgoa71pcCbu3dIciRwZLP6dJJ7x6iWNW0DPL6OrjUWrL+3JnL94772fHbAj8d9/YPot/5BvvNgXjmqo6UxNFYBI/1sqxesVJ0JnDlG11+rJIuqasa6vm5brL+3JnL9E7l2sH5pohmrTp5Lge271rcDHh6ja0mSpHFmrALGT4CdkuyYZBNgDnDBGF1LkiSNM2PSRFJVzyc5BrgUmAR8varuHItrjcA6b5ZpmfX31kSufyLXDtYvTSipqsH3kiRJGgZn8pQkSa0zYEiSpNZtcAEjyaeT3J7ktiSXJdm21zUNR5J/SHJP8x3OT7JVr2sajiQHJbkzycokE2bI3kSe+j7J15M8luSOXtcyEkm2T3J1krubfzvH9rqm4UjyW0luSvLTpv5Te12TtC5scH0wkry0qpY3y38N7FJVR/W4rCFLsh9wVdOR9rMAVXVCj8sasiS/D6wEzgD+W1Ut6nFJg2qmvv9Xuqa+Bw4eq6nv25bk7cDTwD9V1et6Xc9wJZkKTK2qW5JsCdwMHDiBfv4BXlJVTyfZGFgIHFtVN/S4NGlMbXB3MPrCReMlrDEB2HhXVZdV1fPN6g105hiZMKrq7qpaV7O2tmXV1PdV9RzQN/X9hFBVPwKe6HUdI1VVj1TVLc3yU8DddGYLnhCq4+lmdePmNaH+f0caiQ0uYAAkmZvkIeBQ4ORe1zMKRwAX97qIDUB/U99PmF9w65Mk04HdgRt7W8nwJJmU5DbgMeDyqppQ9UsjsV4GjCRXJLmjn9dsgKr626raHvg2cExvq32xwepv9vlb4Hk632FcGUr9E8ygU99r7CXZAjgPOG6NO5HjXlWtqKrd6Nxx3CPJhGuqkoZrrJ5F0lNVte8Qd/0OcBHwqTEsZ9gGqz/J4cC7gX1qHHaiGcbPf6Jw6vsea/ounAd8u6q+3+t6RqqqnkxyDXAAMCE73UpDtV7ewRhIkp26Vv8EuKdXtYxEkgOAE4A/qar/6HU9Gwinvu+hppPkWcDdVXV6r+sZriRT+kZ7JdkM2JcJ9v870khsiKNIzgN2pjOS4QHgqKr6ZW+rGroki4FNgV81m26YYKNg3gt8CZgCPAncVlX797aqwSV5J/CPrJ76fm6PSxqyJOcAM+k8LvxR4FNVdVZPixqGJH8E/Bj4GZ3/bgE+UVX/3Luqhi7JHwDz6Pzb2QhYUFV/19uqpLG3wQUMSZI09ja4JhJJkjT2DBiSJKl1BgxJktQ6A4YkSWqdAUOSJLXOgCFJklpnwJAkSa37/0Am5Q6/mIk9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(V[0,:])\n",
    "plt.text(2,80,'kurtosis=%.2f'%kurtosis_centered(V[0,:]),fontsize=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sudden-fancy",
   "metadata": {},
   "source": [
    "note that V is the learned best loading that satisfies sparsity constraint, it does not equal the mere projection of data onto the learned profiles, which will fit the data slightly better but have a lot more unecessary non-zero loadings.\n",
    "\n",
    "nevertheless, if you want to compare the kurtosis of SCA and PCA loadings, projecting the data onto activity SCs is better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "polish-mills",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/xiaoxiong/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'kurtosis=5.78')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5gcZZn38e+dyQATcJkg2VVGIIgrKCJEoqK4rmHVrItAJKuooIIKq6K+ZDUKiiQgSjSLeNZF3UUF2ajgCIIGMa6yaNBgglmQrMohOBwMhwmSTGCS3O8fVR1qeqq6q6qrD9X9+1xXX5Oprq5+utJTdz2n+zF3R0REpFFT2l0AERHpDgooIiJSCAUUEREphAKKiIgUQgFFREQKMbXdBWinPffc02fOnNnuYoiIlMpNN930gLvPqN7e0wFl5syZrFq1qt3FEBEpFTO7K267mrxERKQQCigiIlIIBRQRESmEAoqIiBRCAUVERAqhgCIiIoXo6WHDItIcw6tHWLp8HfeMjrHX4AAL5x7AvFlD7S6WNJkCiogUanj1CGdesZax8W0AjIyOceYVawEUVLqcmrxEpFBLl6/bEUwqxsa3sXT5ujaVSFpFAUVECnXP6Fim7dI9FFBEpFB7DQ5k2i7dQwFFRAq1cO4BDPT3Tdg20N/HwrkHtKlE0irqlBeRQlU63jXKq/cooIhI4ebNGlIA6UFq8hIRkUIooIiISCEUUEREpBClDShmdoCZXWpmvzOzjWa22cxuM7NPmdlT210+EZFeU+ZO+acBTwW+B/wJ2AocDJwKvN7MDnX3P7exfCIiPaW0AcXdfwL8pHq7mf0c+DZwEvDJFhdLRKRnlbbJq4a7wp/T21oKEZEeU9oaSoWZ7QLsBuwCPBv4RPjUNW0rlIhID+qGGsrbgQ3A3cByYBA40d2vb2upRER6TOlrKMAwcBtBLWUWcAwwI2lnMzuVoOOeffbZpxXlExHpCaUPKO7+J4JRXgDDZnY58GszG3D382P2vwi4CGD27NneupKKiHS3bmjymsDdfwusBt7V7rKIiPSSrgsooQFgj3YXQkSkl5Q2oJjZUxK2zwGeA6xsbYlERHpbmftQvhSmWFlBMPdkF+Aw4PXAX4D3tbFsIiI9p8wB5TLgLcCbCEZ1OUFg+Xdgqbuvb2PZRER6TmkDirt/myDFioiIdIDS9qGIiEhnUUAREZFCKKCIiEghFFBERKQQCigiIlIIBRQRESmEAoqIiBRCAUVERAqhgCIiIoUo7Ux5kV42vHqEpcvXcc/oGHsNDrBw7gHMmzXU7mJJj1NAESmZ4dUjnHnFWsbGtwEwMjrGmVesBVBQkbZSk5dIySxdvm5HMKkYG9/G0uXr2lQikYACikjJ3DM6lmm7SKsooIiUzF6DA5m2i7SKAopIySycewAD/X0Ttg3097Fw7gFtKpFIQJ3yIiVT6XjXKC/pNAooIiU0b9aQAoh0nEKavMxsqpntVMSxRESknFIHFDP7FzP7jZn93MzmhtsONrOfA2PAmJmtNrNXNauwIiLSuVI1eZnZMcCXgE3hY9jMjgSuDI/xw/DnS4ErzexId7++OUUWEZFOlLaG8v+AW4Ehd38KcBlwOTAKHOjux7j7PwEHAw8DZzSjsCIi0rnSBpRnAxe7+yPh758CngJc4O73V3Zy9zuALwOHF1pKERHpeGkDyl8BD0V+fzj8ORKz793AkxoplIiIlE/aYcPrgedGfq/8+xDgqqp9DwUebLBcIqWhzL8igbQB5SrgNDP7A3AvcD7wJ+CFZvYKd/8xgJn9M/A2gs76pjKzZwInAq8E9gd2Af4IfAf4tLtvanYZRHo986+CqUSlDSifAOYBnw1/3wq8DtgI/NjMHiVoPtsVeAz4WMHljPNW4DSC4HUpMA7MAc4DXmdmh7u7suVJIZIunLUy/3b7hbXXg6lMliqguPuDZnYIcCxBf8oKd/8DgJkdBXwQ2ItgJNjH3f3mJpU36rvA+e6+MbLty2b2e+DDBDWlz7egHNLlal04eznzby8HU4mXOvVKeLf/XzHblwPLiyxUyvKsSnhqGUFAeU4LiyNdrNaFc6/BAUZigkcvZP7t5WAq8box2/DTwp/319xLJKVaF85ezvyrNPpSLVVAMbM3mNkezS5Mo8ysDziboI/nWwn7nGpmq8xs1YYNG1paPimnWhfOebOGOP+4gxkaHMCAocEBzj/u4K5u8hlePcIRS1YwMjqGVT3XK8FU4pm719/JbDtBZ/v3ga+6+3XNLlgeZvY54N3Ah9z9/Hr7z54921etSmo5EwlU96FAcOHs9sARJ+5cGOAEwVSjvHqDmd3k7rOrt2dJX7+FYGTXa81sPfA1gtnzfyqojA0xs48SBJOL0gQTkbTatf5IJw7JjetPqgSTG844sj2Fko6RJaCcBtwOnAK8FjgXWGRm1wJfBa509201Xt80ZrYYOAv4T+Ad7SiDdLdWrz/SqUNy1REvtWTqlHf3le7+NuCpBIFlFfAqgiG8I2b2STNraQOqmS0CFgHfAN7uadrwRDpcrZFl7aSOeKkl1ygvd9/k7l9z9xcBBwGfJmhKfT9wa7hGStOZ2dnAYuCbwMnuvr0V7yvSbJ1aE+jlUW1SX8NLALv774D3mdkHCWbTvw14eaPHrcfMTgPOIcgzdh3wRrMJY07ur6SEESmbTp3fovXspZbC1pR3960ETV/fNbNWfLueH/7cB/h6zPM/AxRQpJQWzj0gdmRZJ9QEtJ69JCksoES5e1xa+6Lf4yTgpGa/j0g7qCYgZZQ2oOwHaBagSAu1sybQiUOWpfOlTQ55V7MLIiKdoVOHLEvna0qTl4gk66S7/7iyKIuw5FVoQDGzvyIY6YW7f6PIY0s+nXTxks66+48ry8Lv3Mz49vipXO0esiydr+gaylOBiwmyMSigtFknXbwk0K67/7Q1kaRgAsUMWdYNTncrOn39AwQpWc4t+LiSQ6fOtu5l7ZiwWLmxGBkdw3nixiJunkuSIoYsx5VjwbI1nDW8tqHjSucoNKC4+4PuvtjdzynyuJJPp8627mXtSF2SdGPRZ9XJ55MVkVk5KbHkpSvXM7y66TMNpAW6cYEtCSnvUudpR+qSpBuIbe6TyhJnKFz3pVnlcFCtuUsooHQx5V3qPO1YkCvpBqLy3pWyDA700983sdZS5Pel1o2Mas3dIVWnfLha41/cfTyy7SkEubSOAvYkWHL3SmCxuz/YhLJKRpptnV2zO42zHr+I8tRK41I9ebKZn3/h3ANYsGwNcd3+qjV3h7QrNm4D3uTu3wp/fypwI8H67XeGj2eEv/8eeL67P9KcIhdHKzZKVLNXZsx6/Fr7Q7YbhU4ZXXXW8FouXbl+QlDp1dUvyyxpxcYsSwCfGAkoFwMnAG9298si+50GfA44390/XFDZm0YBRaIq66RXK2o1wqzHP/ScaxkdG5+0fXCgn8e2bi/tksSdEtwkvyKWAI46FvhKNJgAuPsXzOzlBJMbOz6giERlHRWX9cKY5fjDq0digwkQu71MM9mVrbh7Ze6UN7MnAbsD/52wy08JkkmKlEqWUXFJcztqDX/Ncvw8o57UsS3tliWgVNrGNgPjwNaE/bZG9hUpjSyj4vJMGs1y/FrBYded4of6ZunYHl49whFLVrDfGVdzxJIVTZ0H0sr3kvbK0uR1lpmdEv57K3Bgwn4zUap7KaEso+LyTBrNcvykFRsB+vumMNBPpsW3os1zg9P6eXTL1h1pVpqZkkfpf3pL2oCyHpjGE01ZfwZelLDvPwJrGiyXSFPU6/dI276fd4netMdfOPcATl8W/2e0cWycC48/NHX/TfVF/eHNxfTBpOlDUubi3pJ2PZSZafYzs+nA14BfN1AmkaYo8m652Uv0zps1xOIrb4ntgN8rnLmetsxxF/U4Wfpg0p5Lpf/pLUXn8nrY3T/j7r8o8rhSPp3Ybl5kssxWzHhffMxBhWQ6SHvxztIHk/ZcKv1Pb9ECW1K4Tm03L/puOa6WUOQci6IyHdTqj6nIGqjSnstm1+Sks+QKKGY2FXgxcLO7bzSzfoI+lZvdfWORBZTy6dR287z9HmllCaRpA08RczbiLur9fcauO01l49h4rkCV9lwq/U9vyVtD2Z1gvskrgBXAHlW/Sw/r1HbzZt8tpw2kra7BNeOinuVcaiJj72ikyat6MYX0iytIV2t2TaCWWnf+RV1Yk94jbSDNWoMrohmt6Iu6ah4Sp7R9KGZ2JvA84DCC4cx3pR2NJs3VrnbzNHf+jV5Ya71H2kCaNQVL9fstWLaG05etYajNF3HVPKRamddD+ThwJPBH4OE2l0Ui2rHmBzR3yePKqLXTl61JfI85B86YVE2PC6RZU7DErXII6dK9iLRSaWsowP7ufjuAmf0vsFubyyMR7bh7bVbfTVwa+Wojo2NcftPIhJxDBsw/bPJ5yFKDq1f2Zgx2UDZgyau0NZRKMBGpaNachzQTA/vMYmsSP71tchaiLDW4NGUvcrBDnqSXIhVlrqGITNCsvpt6F+yB/r7EgJP02iJSsFQUOdih1UO+VRvqLqWtoeRlZqea2SozW7Vhg3JYdpNm9d3UumBH12XP+to05s0aYnCgP/H5ogc7tHLIt2pD3afnaijufhFwEQQrNra5OFKwZvTdzDlwBpesXD9p+4mH78N58w7e8XuzRrYtPuag2D4cs4mDDor43FmHfDdSw+jUCbCSX5E1FF2cpVCdkg8srh+kenszR7ZVH3twoJ/+PqOyeneRd/ZZ1mxptIbRqRNgJb8iayia2CiF6aR8YGkvfM0c2RY99hFLVkzKQjw2vo3Tl61h6fJ1Lcsf1mgNo50TYKU5cgUUd3+QSO3G3e+nB/tjpHk6qTmk0y58te7gswbepCarNK9ttIahxJHdR0FAOlK7mkPimtmyNAO1Qr1AlnYyZ1yT1YJla5iZsomx0WHa7ZoAK81TWJOXme0JTHf33xd1zDrv9yZg3/DXGcBOZnZW+Ptd7v7NVpRDmjP0s6haQZayJTWznX/cwZx/3MEdM7w17s6+WprAm2YWPiTXdIqoYSh9S3fJHFDM7M3AS9z91Mi284EPhP9eCfyju/+lsFLGexvw91XbPhr+/BmggNICzerrmHPgDC5duX7CSI+sF6t6ZasONpse25rYzHbDGUcWujxu0vNpAmC0nyNpnZNK4K11vEZn4StBpFQz92yDs8zsBmCdu781/H028Cvg58BtBBf689z9nILLWrjZs2f7qlWr2l2MUjtiyYrYi9rQ4AA3nHFkrmPGpTox4ISqYbqNlC3NXX70ve9YclTq961X/rjnB/r7mH/YEJffNDJpe61moKS0MAP9U5h/2NNqHi/p/DTy2aU3mNlN7j67enuePpRnAL+N/P5a4CHgle7+DuCrwOtylVJKpxl9HUlNMUnDd5PUKlvaddYBpphlGrqcVP5LVq7n0HOu5ZyrbomtCV12492Zk1vOmzXE/MMmB5ux8e1csnJ97PEWX3kLQGwyy2oacSVZ5AkouwPRVRn/AbjO3R8Pf18F7NNowaQcmpE/q6gglVSGSgd0WtvcM82zqFXO0bFxHt48HvvctoTWgnqfO2ugHR0b56zhtZOSWVYzgs/czjlAUi55Asp9wN8CmNkM4FDg+sjzuwHpbv2k9JoxAqpekEo74TGubPVMn9a/Y9RRn02+f08zgipvMI17vzTHy1MbjKsNRctgKE2+ZJcnoKwATjOz9wMXE3zvro48fwCgb16PaMbQz1pBKsvs7GjZ0hjo72PR0QdxwxlHcseSo9ies8Yw58AZqd4vyoA3vHDvXME5TwBLqg1td2docGBSzaWodWWku+UZNnw28GLgk+Hv57n7nQBmNhWYD1xeSOmkFFq5vOwRS1ZkmvBYKdt+Z1yd2LxjEDtCKc/Q5eHVI1x+U/b7KQcuXbmewWn97Dx1ChvHxlOPmlo49wD+9dtr2B7zAXeeOoXHtm6ftL3PLDao7DU4oJQoklvmGoq7/wk4CDgEmOnuZ0eengacyhPBRiSXebOGdtQUokN2817skoLA0OAAFx5/KAALlq2Z0ISWpzkvS2d/NQce3jzO6Ng4g9P6mXPgDJYuX8fMM65m/zOvSZxwOG/WELsnZCQe6O+L/Qy1akPNWldGul+umfLuvs3d17r7+qrtj7j79ys1FpGi5b3YLZx7AP1TJvZR9E8x5hw4I7EJLU9zXlF38Q9vHueSlet31JAqtYmkJr7RhI7+0bHxCc1+lYXAfnDzvURPx+BA/47PltRkl6cpT3qLUq9Ix6nV6d7QIIDqPm+DH9x8b82hukk1pSStuIuP689Iet/KR66ct0pgGh0bZ9PjT3zu0bFxFl95C8OrR/jBzffGHitpu0hF3YBiZtvNbKuZ7RT5fVudx9bmF126Ub1O91q1hlqBaOnydYxvm9hnML7NJ2Xtrcg7XDbPyLI8RkbHJgXauDFiTvDZ0zTFjY6Nc+YVaxPPSdJ2kYq6M+XN7GKC7+Xb3X1b5Pea3P3kIgrYTJop33nSzryvTiky58AZNWeF1+qUr6XeTPU4Zw2v5bIb704cSVW0ysz/eksFF2X6tH4WHX1Q6lQxUEy+Ny0X3DmSZspnTr3STTo9oPTiH1DShT+aAiQptUnc6yqBKClQTZ/Wz5bx7TXv3rOkkUlKhdJsSZ+/Wfr7jOOfv3eqVDFJqWayBOoijiHFKTL1irRAr663nabTvVaW3GqVTvKkvpdFRx9Ud65Klo72RkZ5VRhBoJuSYcm6IoPJ9GnBipC1jG/z1Kliaq1tk1YRx5DmyxxQzOzJZvasqm37mdnnzOxSM5tbXPF6V6/+AaXpdM9yga8EokYmYGbpaG90lNfQ4AB3LDmKaTtNjZ1X0mz9U4xFRx/EC2ZOr7tv2lQxRcxr0dyYcsgzsfEzwDOBFwCY2W4EqVf2Cp8/3syOdPefF1PE3tSrf0BpUqInTTisbvapDkRxEzDrNVH1T7FMaWSSypZWZWhuM/+f+6cY4wnRarddgkvCL/74UO7jVwfgtBNEazXxdtqqmRIvT0B5EXBJ5PfjCYLJPwFrgB8TrI2igNKAsv0BNdLfE/faWn0WSQs7zT9siJ/etiFTGeo1UaXtWK98hpHRsUmBLUv/xmU33s2lK9fX37EBScEEgvksS5evy92EZjApAKdZiKve2jVaLrgc8gSUvwGi3/hXAavc/UewY1TYvzZetM7Uqo7yMv0BNbLIVtJrV931UGJwKHJhp3o1ge1O7GeJfg92H+hn0+NbdwxLdp4IIkPhCLRlv7570rDlOEWNDMvbST84rb/uOZk+rT8xY7IT/3++89QpO/6Po6PEKmo18UZrlr02SKVs8gSUcSB6m/z3BEkiK0aBJzdQpo7VrNUJ45TpD6jexSDPa6OrNcad57T5w+qtjJjmolv9Waq/B3HzM5xgVvo9o2P89LYNHP/8vbn6t/cmXoiLljcsbRnfVrNJ8cLjD625OFf14Ia4JsUt45Nzi6Vp4tVywZ0vzyiv/wPmW+AYYA/gJ5Hn9yZYcKvrtLqjPOss7XZppL8naZ8ist2eNbyWBcvWTBopd9bw2h0j6NKKljPtSK7oOiqXrFzP41u37xgQ8DdP2inTZ2mVsfHtsQMjKitOVr6DaTMWpP2bUf6w7pCnhvIFghrJwwTJIG9nYkB5KbC24ZJ1oF7tKK+nkf6eLJ3YWc7z8OqRSWvSwxMrI2ZtWop+lrz/35VUJ4102rdCmtpx2hp02r+ZMjXxllmzm+wzBxR3/4aZbQdeQ7By48fdfRyCIcUEKzp+sbASdpCydZS3SiMXg7jXJrX/ZznPtZqzsgaTgf4+5hw4gyOWrOCe0TGmJKR+7wbTpwVZi9M0LyWNmotesHbpn8JYTBNXdXbkaIAaGR3bkcQymlNNGtOKJvu82YYvcff57v5Wd/9DZPuD7n6Yu3+tkNJ1mGasTtgNGpnjEffaEw7fZ9J57u8zNj22NdXa7sOrR2rWApJWRkzad/5hQ1x+08iOprNuDSYA7qQ6x3HiJuPGBROAuP+CymiuaBLLXpnQ2wqtaLJvKPVKWCPZL/z1Dnd/sJBStUie1Cu9mA6lHaLneXBaP49u2TphuGsl7QZMbHaJy+lV7Yj99+A36zem6gepLL7V6c1UaSQtqlVL3IisJEkd9XGiqXTSHCNL+huJlyatUVqF5vIys0OAzwIvqXrqeuC97v7bzAdtg07P5SWBpIvM4EA/j23dnqq5LCpuzsrmx7fGjsAaClcw7IY6yRH778ENOSYsRodA17qBypKAMylAFHnRk4mKDNaF5fIys+cA/0OwDPCVwMfDx/eBI4DrzeygrMcVSZLUsTs6Np46p1fU2Pg2fnrbhgkj6BYdfVDmFQyz2nlqe1Pnrbz94Vyvqx6+ndT8lPY81Wom1miv5mlFk32eb/i5BHNRnufur3H3j4SP44BZwLZwn6YysylmtsDMbjOzLWZ2t5ldYGa7Nvu9e1GttUaarRkXk+ogVasfqNE1TvrM+PTxh/KJ+c9tsNSNKaLvp1abe9yqmABTLGg6S9O/pn7K5mmkrzOtzE1eZvYA8CV3/0jC8+cB73D3PQsoX61yfAZ4L/A94IfAs4D3EDS7vdzd43sDI9TklU6jqcMb7XdKev9d+qfknihYq5ofV16Ac666Jdf7VZprsvQxdLJazU+zzr028RzVazKrUD9l50tq8sozD2VX4L4az98b7tM0YZPae4Ar3H1+ZPsdBH07rwe+1cwy9JJGZsIXMVQxac4DUHftkcGqtChQ+443qbznH3cwq89+5YScXWk5MPOMq1Pv3+lq1RiT1raH4FwuWLaGVXc9xHnzDk7cTzPiyytPQLkdeDXBBMc4rw73aaY3ENwofbpq+1eAJcCJKKAUppEJnWmCUb30KPXuVJNWKjRgzaJXZrrjTSrvOVfdMuEYvSouGEfPb705Og5cunI9s/fdQ0GjC+UJKN8AzjezbwEfA24Ltz8LOBN4JXBGMcVL9HxgO/Cr6EZ332Jma8LnpSCNTOisF4xqJYeMDv9NqtnMmzWUWGOYYsbw6pFMd7xJ5X148/iOppy4jMK9onpgQfX/X5p+mso69woo3SdPp/y/Ad8haFb6LbAlfNxMUHP4DnBBUQVMsBfwgLs/FvPcCLCnmcUmSzKzU81slZmt2rBhQ1ML2S0a6SitN2onqUaQdjXApPJBcHHLOikube2jklG414yOjU84p3lXqOz1dEXdKnNAcfdt7n48MBf4MsH6J9cBXwJe6e6vT9Mh3qBpQFwwgSC4VfaZxN0vcvfZ7j57xowZTSlct2lkdEi9YJR0YUm7GmC0fHEz4MfGt/G+b9+cenRalhFdnVxDaWawiwb2vIGhl5sNu1meJi8A3P3HBMGkHTYDf53w3C6RfaQgeTtK6yURTGpOS5rVnXQhmjdriAUJfSnVaTyi5UpT3k2PbY1NUd+ppk/r57HxbWxOSHtShEogyZNFQMOAu1eeiY23h2nrk55/tZk1u1P+HoJmrZ1jnhsiaA57vMllkJRqpeFPqsG84YV7Z25mS3PXmyZ3UbS8ZbvwDfT3cdRzn9rUYAJPnOuFcw/IVBtqxtwH6Rx5+lBmArvVeH5XYN9cpUnv1wRlf0F0o5ntAhwKaHJJSSQ1p5037+DMzWxpm6vSNtNUOpzLUjsZGhxg/mFDXHbj3U19n2hgnzdrKFPTXyev6yONy93kVcPf0PzmpmXAh4DTCSYyVpxC0HdyaZPfXwqU1JyWtZmturkqaQhr2vb7vB3O7TIyOha7BkwacfN14hgw/7CJ/y9DKZu9qldzlO6TKqCY2UuBl0U2HWdmz4jZdQ+C0V/xjdkFcfe1ZvYF4N1mdgVwDcGw5fcCP0NzUHpWNAglzbBP24xVxpFIWYPJiYfvs2OSYZpJmw5cduPdXLpy/Y7+sLg1baqp36Q3pEq9YmaLgEXhr/VGTP4BeKO7N7XZycz6CGoopxI0wz1AUHM5290fTXMMpV5Jr6zpMBopd7ekSqmnkqI/em5qpVCJqrWMQDSbc1m+L5JOQ+nrzWx3YJDgu3c7sAAYrtrNgUfdvTTrySugpNNoLq80x+/EYBX3ubt5QmP0/7SIVPTSvRrK5eXuG4GNYa3gZOAad9eswB7RSC6vJNHmlehFuhnLkmYpT1xQi27v5hrL2Pg2PnTFb5k3ayjTZy1j06A0R9ZRXv3A14CTii+KdKpGcnnFiS4VC5Pv+IteljRLeSpL11Zmg1cPea7XsVz22fObx7dz1vDaTBM8NUlRKjKN8gpzZT0AbGpSeaQDNZLLK06a0VOtvOutVwObUJuqEzG6oTnsshvv3tFRX2955TSd7Z3apCnFyzNs+BqCjMJfLLgs0qHiRvE0MmonTbBo5V1vUtPOPaNjk/pRClijquNVhlrHDdueve8eO4Jrn9mE2mRckChi+QIpjzwTGz8APNXMvm5mB4eTCXteO1c0bLaiV3qrFyxaOcR0ePVIYjPVXoMDpZuLUoS4nGgV0RUsq1PaxH3na9X+pPvkqaH8maBmfwjBuiPY5C+gu3szJk12pF64Cyty0aO4Gk+lYz7tqn5FWbp8XWIz1cK5ByTmB2u1/imQJZtK1v2j3vDCvWs+n3bNmDkHzqhZ+5Puk3c9lB6o+KfXjFFQ3axewshWqndhqzfaqRXDiPv7LHyT9O+UN5hM65/C7H33qLlP2jVjLlm5PvEY6sjvTpkDiruf1IRylFrRo6B6QdHLvObt+K0VMJYuX1dzFngrgkmfGbvuNLVl+cQ2j29PrF1XznGjn1mz5rtXnj4UqVJvESlprlrDfmu9pt5M+HtGxyb0H8ET/QtDgwMtqZlc8LpD2Nji5JRxfRzVQ70boWzD3StzDcXM9kmzn7sn13e7TNGjoCSbrE2OcTPg40wxY78zrk6s8TQzNcv0af0sOvogoHYtqFm1pOradb3BCWbpRsANDQ4omHSxPH0od5LuO5xuVlQX6KQ+gV6Utckx7citWgtzDa8eYfPjWye9ZqC/L/eosGgKk+HVIyy+8pa6TV3NqiVV167rNd+mCSa6yep+eQLKuUz+Hk8F9geOBdYCP2ywXKVTdJ9AGpowFsg68TJP31ZlKeEFy9awe0Kq9/UIryEAABSSSURBVMGBfhYfc1DdjL1xohfbtDWoZqq+8KdNxdJnxnb3XAki9X0uvzyd8ouTnjOzpwO/RAtcNV0vDFVOK2uTY96cXJUaS1KtYdedp+449wu/e3PdtUWGBgcSc4c1EkyiQ7DzfE4DFixbs2NQQmXuSZogt92dO5Yclfk99X3uDoV2yrv77cC/A+cUeVyZTBPGnpB14mWtPFUD/X0MDvTnKkel5jNv1hC77lT7Xq3SvBW3LHIjowOnGFx4/KHcmTL3WBwPH9WDG3aeWv9yUWQ6nl79PpdZMyYfjgDPbsJxJUJDlSfK0uQYrQlUUohsc98xqRLI1eQUvZjWGplVry+hkazG233iHX3amkWSyoTFLePb6x6jGel4evX7XFbNCCjzgIebcFyJKDphY69JE4DqLSUcNdDfx5wDZ3DEkhU1X9NnVrP2lNTZn8URS1bsaKqKGzCy+fGtqRbPqkizb73PVY++z90hz7DhsxOe2gM4EngO8MlGCiX1aahyc9VbSrh/irHbLlMZ3Twem4k3LpjUW5SsqM746v6HyqPS6Z0lmKRRxGJr+j53hzw1lMU1nrsPOAv4RK7SSGoaqtw6ac71EUtWxAaC6Kinev8/RSairJ6HkyZYTZ/WP6lpa6C/j52nTkkciFBU7jV9n7tDnoCyP1CdKciBhypruZvZXwGPNFg2qaMdQ5V7Vb1zndTWn2XUU9H9BdHj1QtWA/19OyZSVl/UYXKfUpFLQFfo+1x+eQLK6e7+/5KeNLMnAT8CXpy7VNKTOmkeQtayFNEHkKUzfted+tj8+Db2Ghxg02NbY2sQ0feuF6zmH/bExTzpc3bK/410rjzDht9jZgvjnjCzXQkmNR7aUKmk58Tl41qwbA0z27C+TJ7cYHFDkbP2AWRZdndw2k47hhwvPuaguu9dL7D99LYNsdsrOc8qafxPODzIvLRg2ZpU/y/dvE6QTJYnoHwEWGJmJ0Q3mtkAcDUwG5hfQNmkh8Q1yVS6tdNc0JtdlnpzIopYhKxyjOnT6s+DidY40rz3nANnpD5eRVxgvWTl+gm/n75sDbPOvTb2/yZPYJZyyzNT/mNmNgR8zczud/frwlUbrwJeBPyzu/dc6hVpTL0mmVauL5N3TkQRfQDVI7KSmsCqaxy13nt49QiX31T7Ih5Xg0k7SODhzeOxs9q1TlDvyTtT/jSCpq3LzezFwPeBlwJvdPeriipcLWb2L2Z2qZndZmbbzEyLfpVYmr6GVk1ya/dyBNH+m+nT+umfMnlJ3s2Pb02Vnn+/M67mfd++uWZg6J9isU1zWc53XA1OkxV7T66A4u4OvIEgEeTPgTnAie5+eYFlq+dM4BiCJYnvaeH7ShOk6T9o1QW9iP6Qiqx9CNXNRA9vHgeDgf6Jf6qVWkGapqZakzIHB/pZ+tpDYmsMWc93daBod2CW1qsbUMzspXEP4AXABcAm4D+B+6qeb7aXAbu7+0uBm1vwftJE1QtZVd+Tt3KSWxH9IZCvDyGumWh8m/P41slBIalfJ21T1dDgAGsWvbJmzrMsqgNFkYFZyiFNH8p/U3+Nn1OAt0d+d5q8Hoq739nM40vrVc9Ob+cw1SL6Q/L0ISQ1ByXVMuL2T9OklObCPm/WEOdcdUuqmfVxx9Nkxd6TJqCc3PRSSM+pFzC6YZJbnj6EpLkofQm5weKaj2odI+2s/YpFRx8UO6lx/mFDqdY6yfr/2O4bCWlM3YDi7l9vRUFaxcxOBU4F2GefVKsZS8F6Ze2LPJMdk3JazT9saEKusMr2uFpG0jHyNNu1spbRK9+LbmaeZu3OZr252SBweoaXfNbdH4o5zg+Ao9x98nCYGmbPnu2rVmktsFZLWos9ugRuuxR5hxyXPyvNhT2pDFnKVsY7/U7+XshEZnaTu8+u3t6M9PVZDAKLMux/CTApoEi5dOpw0qLvkPPe3Sc1E2Vd86XTA0i1Tv1eSHptDShhx3qmWoWUX6N5r5p1992MiXhlvLC3SzPWRCljTa3MCl0CWCSNRoaTNjOdh+6Q26voYcZK/dJ6CijSco3M82jm2uNFTcRTQsR8ipr/U6F16luv3X0ouZnZ0cAh4a/PCLedFf4+6u6fb0vBJJW8TUHNrEUUsWpgo/0wvd5EU2QToWqcrVfagEKQ0fgtVds+Gv68C1BA6ULNXHu8iCGyjfTDaNhssbROfeuVNqC4+0nASW0uhrRYs9ceb/QOuZG7YmXnLZbWqW+90gYU6U2dns6jkbtiNdEUq9O/K91IAUVKp5OH4jZyV6wmmuJ18nelG2mUl0iBGhmppOy8UnaqoYgULO9dsZpopOwUUKTlen1obC1qopEyU0CRltLQ2O6lGwVRQJGWatfQWF3s4hV1XnSjIKBOeWmxdgyNVU6neEWeF6U5EVBAkRYrKl9WFrrYxUs6L6cvW5M5B5nm0AgooEiLtWNorC528Wp9/qy1lXbcKEjnUUCRlio6o2wautjFq/f5s9TiNIdGQJ3y0gatHhqrnE7x4s5LtbS1OM2hEVBAkR6gi1286HmJS/kC2WpxmkMj5u7tLkPbzJ4921etWtXuYoi0XfWwXwhqcc1ujpRyMrOb3H129XbVUEREtTgphAKKiABqspLGaZSXiIgUQgFFREQKoYAiIiKFUEAREZFCKKCIiEghNMpLpEN1Usr9TiqLdC4FFJEO1Enri3RSWaSzqclLpAN1Usr9TiqLdLZSBhQzGzKzM83sZ2Z2r5ltMrNbzGypmT253eUTaVQnpdzvpLJIZytlQAGOBhYDDwJLgdOBX4Q/V5vZU9pXNJHGdVLK/U4qi3S2sgaU64F93f04d/+Uu3/F3U8B3gnsDby/vcUTaUwnrS/SSWWRzlbKTnl3vyXhqWXAV4DntLA4IoXrpGSNnVQW6WylDCg1PC38eX9bSyFSgE5K1thJZZHOVdYmryTnhD+/nrSDmZ1qZqvMbNWGDRtaVCwRke7X1hqKmQ0SdKSn9Vl3fyjhWO8DXgtc5O4rkg7g7hcBF0GwwFaG9xYRkRra3eQ1CCzKsP8lwKSAYmZvJxjtdTXw7mKKJiIiWbQ1oLj7nYA1cgwzeytBjeNaYL67jxdQNBERyajUfShmdjLBqK7rgHnu/libiyQi0rNKG1DM7CTgq8BPgWPdfUt7SyQi0tvMvXz90mZ2DPA94BHgA0B1DohH3X04xXE2AHfV2GVP4IG85ZQddB6Lo3NZHJ3L/PZ19xnVG9vdKZ/X8whqV4OEI7aq3AXUDShxJyTKzFa5++xcJZQddB6Lo3NZHJ3L4pWyycvdF7u71XjMbHcZRUR6TSkDioiIdB4FlNrimtMkO53H4uhcFkfnsmCl7JQXEZHOoxqKiIgUQgFFREQKoYAiIiKFUEDJwcymmNkvzczN7AftLk+ZmNmQmZ1pZj8zs3vNbJOZ3WJmS83sye0uX6cJv2sLzOw2M9tiZneb2QVmtmu7y1YmZvZMMzvXzFaa2QYz+4uZrTGzD+tcFkcBJZ93AQe3uxAldTSwGHiQIEP06cAvwp+rzewp7StaR7oQ+BRwK/Ae4DvAe4GrzEx/v+m9FVgA/BE4F1gIrAPOA35hZgNtLFvX0CivjMzsaQR/3IuBC4Cr3f3VbS1UiZjZQcCD7n5f1fa3EyT6vMDd39+WwnWY8FytBb7n7vMj298DfBY4wd2/1a7ylYmZzQZ+7+4bq7afB3wYeI+7f74thesiusPJ7vPA7cBn2l2QMnL3W6qDSWhZ+PM5rSxPh3sDwfIOn67a/hVgM3Biy0tUUu6+qjqYhPS9K1BZc3m1hZn9M3AM8GJ332bW0FIuMtHTwp/3t7UUneX5wHbgV9GN7r7FzNaEz0tj9L0rkGooKZnZ7gTNDP/u7ivbXZ4udE748+ttLUVn2Qt4IGGdnxFgTzPbqcVl6hpm1gecDWwF1HRYgJ6qoTS4hv0nCQLwmYUXrIQaPJfVx3of8FrgIndfUUT5usQ0IGnRuC2RfR5vTXG6zqeBw4EPufu6dhemG/RUQCHnGvZm9hLgFOBN7j7alJKVT65zWb0x7IxfClwNvLuYonWNzcBfJzy3S2QfycjMPkrwfbvI3c9vd3m6RU8FlAbWsP8CcDNwo5k9o+q5aeG2UXfvmcV6GjiXO5jZWwkS9F0LzHf38QKK1k3uAZ5tZjvHNHsNETSHqXaSkZktBs4C/hN4R3tL013Uh5LOvsChwO+rHgBzwn8vbkvJSsrMTiYYrXQdMC+hn6DX/Zrgb/QF0Y1mtgvB93FVOwpVZma2iKBm/Q3g7a55E4XqqRpKA94MxHV+fge4CVgC/KGlJSoxMzsJ+CrwU+BYd99S+xU9axnwIYK+qusj208h6Du5tB2FKiszO5vgxu+bwMnuvr29Jeo+mtjYADNzNLExEzM7Bvge8AjwAWCsapdH3b3u8s29wsw+R9DW/z3gGuBZBDPlbwCO1EUxHTM7jWAO2XrgIwTDsaPud/cft7xgXUY1FGm15xE04wwSv8DRXYACyhNOB+4ETgWOAh4APgecrWCSSWXOzj7ED03/GaCA0iDVUEREpBDqlBcRkUIooIiISCEUUEREpBAKKCIiUggFFBERKYQCioiIFEIBRURECqGAIh3PzE4yMzezl7W7LM1gZheHWRdESk0BRaSKmQ2a2eJuDWBRZjbNzM42s1vMbMzMHjKzX5rZazIe501mdoOZPWJmj5rZ/5rZR2L229PMPmlmt5nZZjO7z8xWmNmxxX0qaRelXhGZLLrWy3+34P1OoQ1p1M1sOvAT4G8JUrl/CtiVIF/YvhmO8x/AW4DLCRJWbgP2qz6GmU0DfgHsTZBp+rfAHsBJwLCZvcvdv9TQh5K2UkARCZnZANDyNVnCdWDasRbMZ4FnAC9091vzHMDM3gacDLzZ3b9ZZ/djCYLX6e7+mcgxLgL+BPwLoIBSYmryktIysw+HfSufM7Onh/9eHLPf4vC5mZFtF4fbZpjZf5jZ/cAm4ETgjnC3ReE+bmZ3Rl471cw+aGa3mtkWM3vQzL5nZgfHvPebzexXZjZqZpvM7HYzu9TMZlSXpep1e4flusvMHjOzP5vZL8zsLY2dtR3Hnwm8EfiKu99qZn1mtlvGYxjBkti/qQQTM3tSuD3OX4U/76navpHg3G/K8v7SeVRDkdIxsykEqcjfCZzp7kuiwSKjHwP3AR8laO75IbAAuJAgZfwV4X6PRl5zKfC68LVfAp4CnAb80sz+zt1Xh+U8kSCz7fXA2QSp+vcBXkWwtO+GhM83NTz2EPBF4P+A3YHnAn9HJFuume2Z4bNujKyK+Y8EN5S3mtk3w8+zk5mNABe4+4UpjncAsD/w+bC/5HSCJqxHzOwy4P3uHj1vK4CtwPlmtomgyWs6wfkeBD6W4bNIJ3J3PfTo6AdBG7sDLwMGCC70jxM0s1T2mRnuszjm9YvD52ZGtl0cbrskZv9ax3pF+Nwywmzd4fbnElwsr49su4Jg3ZepdT7fxcGf4oRjOfCBFOfGMzxeFnndheG2PwO3EjRbnUAQ/Bw4J8V7Hxs5xijwfmA+QX+MEyygZlWvOQ64u6pc9wFHtPt7pkfjD9VQpEz2ILhzPwQ42t2XF3DMf8u4f2X008c8vEICuPtvzewHwLFmNsPdNxA05UwDjjKzK6P717Ex/DnHzC529z/X2PcVGcp+c+TfTwp/7gT8nbs/CGBm3yYIMB8ws0+7+8M1jlc5xgzgFe5+Xfj75WGz11sIakI/jLxmlKBm8lVgDUEt7H3A983sH9w9WkYpGQUUKZOLgd2Al7r7/xR0zP/LuP9+BKv9/S7muf8luGvfj6A56+PASwkWDHvQzH5GcHFd5u5/SXoDd7/LzD5G0D9xr5mtIRiN9R13/3XVvtfFHSOFykqZP6gEk/B442b2LYImusOZGAySjjESU46vEwSUl1WOYWZzCVadPMrdf1TZ0cyuAG4DvgC8JOfnkQ6gTnkpk2UEF/OzwxFZUbXu/hNvnNx9c8YyJHU4xx3798CzCVZa/DrBMNqvALeZ2f51XnsW4Ygo4I/A24FfmdknJhTG7CkZHjtFXvqn8Od9MW9/b/hzep2PmPUYHwQ2RYNJ+FnvI2hqe1FVGaVkFFCkTC4lGIU1B/hBOK+h4qHw5x4xr3t6xvepFZz+SPB386yY554d/qyMEsPdH3P3a9z9fe4+myC47AX8a91CuN/u7p9z99eFr/k5QVPUX0d2uzfD48WR1/0q/Pm0mLeubKvV1AawlqCWMpTyGEPAlIRRYFMJzquuSSWmJi8pFXf/LzPbCnwL+KGZHeXuj7r7X8zsPuBIM7NKf4WZPR2Yl/FtKiOT4oLTMPAu4Ewze2PkfZ4DHAP8T9h/gpnt6e4PVL3+NzWOTfi63YHN/sSILNx9i5n9jqAJbTpPXKjz9qH8HLgLONrMhtx9JHzvXYE3E/R1/DJSpj2BPYF73X1jWKbNYXPVCWb2Gnf/XuT47wx/XhPZdivwTOC1wLcjx94v/Fxr3X1Lhs8jHUYBRUrH3b9rZuMEF6XlZvYqd3+EYCjxeQSBZpjgrv4dBH0bz89w/AfN7A/A683sj8D9BE01V7n7j8OO69cD08OO+Mqw4S3AeyOHutbMNhJcvO8mGBp7EkENqNYkwDnARWZ2ObCOIMAdRtDsdaO7r4uUNVcfirtvM7N3AlcRDHf+IsHIuZMJZrK/zd2j80LeTZA94GSCvqyKDwEvB75lZp8D7gT+iaAm9g13/0Vk348TdNJfYkFamzUENZl3AruEx5Iya/cwMz30qPcgMmy4avtRBBfxlQQX66nAJwmad7YQ1AaOpsaw4Rrv+QLgBoLJdg7cGXluKkF/wO+Axwia24aBg6uOcQpPzHN5PCzXNcCcqv0mlIWgU//L4fEfCcvwO+BcYPeCz+3fE8wP+QuwOfzMR8fsVzmHJ8U8N5OgOXJD+Dl/RzCEeErMvrOA7xIE6a3Aw8CPqv9v9Sjnw8L/ZBERkYaoA0xERAqhgCIiIoVQQBERkUIooIiISCEUUEREpBAKKCIiUggFFBERKYQCioiIFEIBRURECvH/AYjsYcDPVee4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V_proj=np.linalg.lstsq(U,tmpX)[0]\n",
    "plt.scatter(V_proj[0,:],V_proj[1,:])\n",
    "plt.xlabel('kurtosis=%.2f'%kurtosis_centered(V_proj[0,:]))\n",
    "plt.ylabel('kurtosis=%.2f'%kurtosis_centered(V_proj[1,:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "severe-throw",
   "metadata": {},
   "source": [
    "# grid search"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "generic-arrival",
   "metadata": {},
   "source": [
    "You SHOULD MODIFY it according to your data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "senior-amendment",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current Time = 23:30:11\n"
     ]
    }
   ],
   "source": [
    "from  datetime import datetime\n",
    "now = datetime.now()\n",
    "\n",
    "current_time = now.strftime(\"%H:%M:%S\")\n",
    "print(\"Current Time =\", current_time)\n",
    "\n",
    "# this folder is where you have put your data and where you will store the grid search results\n",
    "folder='demo'\n",
    "\n",
    "#these are specific to my analysis, which you can just leave here or change it to your liking\n",
    "\n",
    "#the optimization method, here I use ADMM followed by normal gradient descent. It's just a name for my log file.\n",
    "version='MixOpt'\n",
    "type_restructure='s100s50'\n",
    "stitching='SDT'\n",
    "\n",
    "#here you may want to iterate through all the stimuli and there interaction\n",
    "for Stim in ['ST']:#['DT+SDT']:\n",
    "    \n",
    "    #here I wanted to compare L1 and L2 regularizer, you could just use L1\n",
    "    for regu in ['L1']:#,'L2']:\n",
    "        \n",
    "        #this the data filename\n",
    "        filename=folder+'/splitX.pkl'\n",
    "        #+stitching+'_'\\\n",
    "        \n",
    "        #here my data contains already the training and test set.\n",
    "        #you could just import the raw data and split them into train and test set here\n",
    "        with open(filename,'rb') as f:\n",
    "            tmpX1=dill.load(f)\n",
    "            tmpX2=dill.load(f)\n",
    "    \n",
    "        # the name for log file\n",
    "        name='GS_'+Stim+'_'+stitching+'_'+regu+'_'+type_restructure+'_'+version;\n",
    "        \n",
    "        if regu=='L1':\n",
    "            lossname='lossEla'\n",
    "        else:\n",
    "            lossname='lossL2'\n",
    "        \n",
    "        #number of runs.      if the data is well-structured, you should change this loop into repeat&fold by e.g. using the method in sklearn\n",
    "        n_run=1#4\n",
    "        \n",
    "        max_components=3#11\n",
    "        \n",
    "        #how fine grained you want to explore the different alpha values in [10^alpha_s,10^alpha_e]\n",
    "        n_alphas=3#80\n",
    "        alpha_s=-1\n",
    "        alpha_e=0.5\n",
    "        \n",
    "        #log training loss\n",
    "        lossall=np.zeros((n_run,max_components,n_alphas));\n",
    "        valossall=np.zeros((n_run,max_components,n_alphas));\n",
    "        \n",
    "        #log testing loss\n",
    "        errall=np.zeros((n_run,max_components,n_alphas));\n",
    "\n",
    "        \n",
    "        \n",
    "\n",
    "\n",
    "\n",
    "        os.mkdir(folder+'/ckpt_'+name)\n",
    "        os.mkdir(folder+'/ckpt_'+name+'/weights')\n",
    "        \n",
    "        for n_run_i in range(n_run):\n",
    "\n",
    "            current_time = datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\")\n",
    "            logfile=open(folder+'/logfile_'+name+'.txt','a')\n",
    "            logfile.write('n_run %d/%d   \\t\\t\\t\\t\\t %s\\n'%(n_run_i,n_run,current_time));\n",
    "            \n",
    "            #here I switched train and test set, which may not suit you (you may have more than 2 folds)\n",
    "            tmpX1,tmpX2=tmpX2,tmpX1\n",
    "\n",
    "            for j,n_components in zip(range(0,max_components),range(1,max_components+1)):\n",
    "                current_time = datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\")\n",
    "                logfile.write('\\t component %d/%d\\t\\t\\t\\t\\t %s\\n'%(n_components,max_components,current_time));\n",
    "                for k,al in zip(range(n_alphas),np.logspace(alpha_s,alpha_e,n_alphas)):\n",
    "                    current_time = datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\")\n",
    "                    logfile.write('\\t \\t alpha %d/%d\\t\\t\\t\\t\\t %s\\n'%(k,n_alphas,current_time));\n",
    "                    \n",
    "                    \n",
    "                    U,V,losses,errs=dSCA.Learn(tmpX1,tmpX2,True,n_components,al,lossname=lossname)\n",
    "                        \n",
    "                    #log training and validation loss\n",
    "                    lossall[n_run_i,j,k]=losses[-1]\n",
    "                    errall[n_run_i,j,k]=errs[-1]\n",
    "                    \n",
    "                    \n",
    "                    \n",
    "                    with open(folder+'/ckpt_'+name+'/weights/nn%02dweights%02d%02d.pkl'%(n_run_i,j,k),'wb') as f:\n",
    "                        dill.dump(U,f)\n",
    "                        dill.dump(V,f)\n",
    "                        dill.dump(losses,f)\n",
    "                        dill.dump(errs,f)\n",
    "\n",
    "            logfile.write('\\n')\n",
    "            logfile.close()\n",
    "\n",
    "\n",
    "        with open(folder+'/ckpt_'+name+'/losses_TrVa.pkl','wb') as f:\n",
    "            dill.dump(lossall,f)\n",
    "            dill.dump(errall,f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expired-hotel",
   "metadata": {},
   "source": [
    "# Best parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "third-conducting",
   "metadata": {},
   "outputs": [],
   "source": [
    "folder='demo'\n",
    "n_iter=1200;\n",
    "type_restructure='s100s50'\n",
    "version='MixOpt'\n",
    "\n",
    "for Stim in ['ST']:#,'DT']:\n",
    "    for regu in ['L1']:#,'L2']:\n",
    "        \n",
    "\n",
    "    \n",
    "        \n",
    "        name='GS_'+Stim+'_'+stitching+'_'+regu+'_'+type_restructure+'_'+version;\n",
    "        \n",
    "        #load the logged losses\n",
    "        with open(folder+'/ckpt_'+name+'/losses_TrVa.pkl','rb') as f:\n",
    "            lossall=dill.load(f)\n",
    "            errall=dill.load(f)\n",
    "            \n",
    "        \n",
    "        #tmp is the matrix containing the loss for each parameter combination, averaged across repeats\n",
    "        tmp=errall.mean(0)[:,:]\n",
    "        \n",
    "        #extract the best parameters which you can use to run dSCA\n",
    "        n_components,ind_alpha=np.unravel_index(tmp.argmin(),tmp.shape)\n",
    "        n_components+=1\n",
    "        al=np.logspace(alpha_s,alpha_e,80)[ind_alpha]\n",
    "        \n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adolescent-resolution",
   "metadata": {},
   "source": [
    "rerun dSCA using the best parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "international-duplicate",
   "metadata": {},
   "outputs": [],
   "source": [
    "U,V,losses,errs=dSCA.Learn(tmpX,n_components=n_components,alpha=al)\n",
    "#this uses the default (arbitrary) hyperparameters, \n",
    "#you should choose your hyperparameter for a good reason, or do a grid search for it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "inside-simulation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(180, 3)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "confidential-asbestos",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(3,1,figsize=(2,3*1.6),dpi=100)\n",
    "fig.tight_layout()\n",
    "for i in range(3):\n",
    "    print_sig_ax(ax[i],U[:,i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "persistent-landscape",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}