Prediction_PowerSpec.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuda:0\n"
     ]
    }
   ],
   "source": [
    "import argparse\n",
    "import os\n",
    "import random\n",
    "import shutil\n",
    "import time\n",
    "import warnings\n",
    "import os\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.backends.cudnn as cudnn\n",
    "import torch.optim\n",
    "import torch.utils.data\n",
    "import torch.nn.functional as F\n",
    "\n",
    "from torch.utils import data\n",
    "import random\n",
    "import numpy as np\n",
    "from itertools import product\n",
    "import argparse\n",
    "\n",
    "from train_f import *\n",
    "from Dataset import Dataset\n",
    "#from Models import *\n",
    "# from Models2 import *\n",
    "# from Models4 import *\n",
    "from args import args\n",
    "from Models import *\n",
    "\n",
    "device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')\n",
    "print(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(filenames, model_name, C_model = None, vel = 0, redshift = 0, label = ''):\n",
    "    dim_out = 1\n",
    "    dim_in = 1\n",
    "        \n",
    "    for i, filename in enumerate(filenames):\n",
    "        if model_name == 'reg_Incep':\n",
    "            model = Inception(dim_in, 6, 8, 10, reg = 1).to(device)\n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "        elif model_name == 'reg_R2Unet':\n",
    "            model = R2Unet(dim_in,dim_out,t=3,reg = 1).to(device)         \n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "        elif model_name == 'subhalo_r2unet':\n",
    "            model = R2Unet_subhalo(dim_in,dim_out,t=3,reg = 1).to(device)         \n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "        elif model_name == 'R2Unet':\n",
    "            state_dict_mask = torch.load('pretrained/' + C_model)\n",
    "            mask_model = R2Unet(dim_in, dim_out, t= 3).to(device)\n",
    "            mask_model.load_state_dict(state_dict_mask)\n",
    "            pred_model = R2Unet(dim_in,dim_out,t=3,reg = 1).to(device)\n",
    "            model = two_phase_conv(mask_model,pred_model)\n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "        elif model_name == 'Incep':\n",
    "            state_dict_mask = torch.load('pretrained/' + C_model)\n",
    "            mask_model = Inception(dim_in, 6, 8, 10).to(device)\n",
    "            mask_model.load_state_dict(state_dict_mask)\n",
    "            pred_model = R2Unet(dim_in,dim_out,t=3,reg = 1).to(device)\n",
    "            model = two_phase_conv(mask_model,pred_model)\n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "        elif model_name == 'subhalo_2phase':\n",
    "            state_dict_mask = torch.load('pretrained/' + C_model)\n",
    "            mask_model = Inception(dim_in, 6, 8, 10).to(device)\n",
    "            mask_model.load_state_dict(state_dict_mask)\n",
    "            pred_model = R2Unet_subhalo(dim_in,dim_out,t=3,reg = 1).to(device)\n",
    "            model = two_phase_conv(mask_model,pred_model)\n",
    "            state_dict = torch.load('pretrained/' + filename)\n",
    "            model.load_state_dict(state_dict)\n",
    "\n",
    "        with torch.no_grad():\n",
    "            for ID in ranges:\n",
    "                if vel == 1:\n",
    "                    inp=dark_velo[:,ID[0]:ID[0]+32,ID[1]:ID[1]+32,ID[2]:ID[2]+32]\n",
    "                    inputs=torch.from_numpy(inp).unsqueeze(dim = 0).to(device).float()\n",
    "                else:\n",
    "                    #inp=dark[ID[0]:ID[0]+32,ID[1]:ID[1]+32,ID[2]:ID[2]+32]\n",
    "                    inp=dark[ID[0]:ID[0]+32,ID[1]:ID[1]+32,ID[2]:ID[2]+32]\n",
    "                    dark_count=torch.from_numpy(inp).unsqueeze(dim = 0).unsqueeze(dim = 1).to(device).float()\n",
    "                    sub=subhalo[ID[0]:ID[0]+32,ID[1]:ID[1]+32,ID[2]:ID[2]+32]\n",
    "                    sub=np.log(sub+1)\n",
    "                    subhallo_mass=torch.from_numpy(sub).unsqueeze(dim = 0).unsqueeze(dim = 1).to(device).float()\n",
    "                output = model(dark_count)\n",
    "                #print(' output.size= ', output.size())\n",
    "#                 outputs = F.softmax(output, dim=1)\n",
    "        #     predicted = outputs.max(1, keepdim=True)[1]\n",
    "                predictions[:,ID[0]:ID[0]+32,ID[1]:ID[1]+32,ID[2]:ID[2]+32]=output.cpu().numpy()[0]\n",
    "                if (ID[1]==0 and ID[2]==0):\n",
    "                    print(ID)\n",
    "                #print(output.cpu().numpy()[0].shape)\n",
    "        predictions[predictions < 0] = 0\n",
    "        #np.save('/scratch/ys3202/predictions_' + 'red' + str(redshift) + '_' + filename + '.npy', predictions)\n",
    "        np.save('/scratch/xz2139/predictions_' + label + '_' + filename + '.npy', predictions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## modify\n",
    "normalize = 0  #modify\n",
    "vel = 0 #modify\n",
    "redshift = 0  #modify \n",
    "tng = 0\n",
    "if tng == 1:\n",
    "    dark=np.load('/oldscratch/xz2139/TNG_100_dark.npy')\n",
    "elif redshift == 0:\n",
    "    dark=np.load('/oldscratch/xz2139/cosmo_dark/Dark_zeros.npy')\n",
    "    subhalo = np.load('/oldscratch/xz2139/Fullmass_zeros.npy')\n",
    "if  normalize:\n",
    "    dm_mean=5.614541471004486\n",
    "    dark=(dark-dm_mean)/dm_mean\n",
    "if vel:\n",
    "    dark_velo = np.load('/oldscratch/xz2139/Dark_velo_zeros_4.npy')\n",
    "\n",
    "predictions=np.zeros((1,1024,1024,1024))\n",
    "data_range = 1024\n",
    "pos=list(np.arange(0,data_range,32))\n",
    "ranges=list(product(pos,repeat=3))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 0)\n",
      "(32, 0, 0)\n",
      "(64, 0, 0)\n",
      "(96, 0, 0)\n",
      "(128, 0, 0)\n",
      "(160, 0, 0)\n",
      "(192, 0, 0)\n",
      "(224, 0, 0)\n",
      "(256, 0, 0)\n",
      "(288, 0, 0)\n",
      "(320, 0, 0)\n",
      "(352, 0, 0)\n",
      "(384, 0, 0)\n",
      "(416, 0, 0)\n",
      "(448, 0, 0)\n",
      "(480, 0, 0)\n",
      "(512, 0, 0)\n",
      "(544, 0, 0)\n",
      "(576, 0, 0)\n",
      "(608, 0, 0)\n",
      "(640, 0, 0)\n",
      "(672, 0, 0)\n",
      "(704, 0, 0)\n",
      "(736, 0, 0)\n",
      "(768, 0, 0)\n",
      "(800, 0, 0)\n",
      "(832, 0, 0)\n",
      "(864, 0, 0)\n",
      "(896, 0, 0)\n",
      "(928, 0, 0)\n",
      "(960, 0, 0)\n",
      "(992, 0, 0)\n"
     ]
    }
   ],
   "source": [
    "model_subhalo = ['regression1.pth']\n",
    "C_model = 'model_2@6@8@10_80_v0_n0.pth'\n",
    "predict(model_subhalo, C_model=C_model, model_name = 'Incep', vel = 0, label = '1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def power_spectrum_np(cube, mean_raw_cube, SubBoxSize):\n",
    "\n",
    "    nc = cube.shape[2] # define how many cells your box has\n",
    "    delta = cube/mean_raw_cube - 1.0\n",
    "\n",
    "    # get P(k) field: explot fft of data that is only real, not complex\n",
    "    delta_k = np.abs(np.fft.rfftn(delta)) \n",
    "    Pk_field =  delta_k**2\n",
    "\n",
    "    # get 3d array of index integer distances to k = (0, 0, 0)\n",
    "    dist = np.minimum(np.arange(nc), np.arange(nc,0,-1))\n",
    "    dist_z = np.arange(nc//2+1)\n",
    "    dist *= dist\n",
    "    dist_z *= dist_z\n",
    "    dist_3d = np.sqrt(dist[:, None, None] + dist[:, None] + dist_z)\n",
    "\n",
    "    \n",
    "    ################ NEW #################\n",
    "    dist_3d  = np.ravel(dist_3d)\n",
    "    Pk_field = np.ravel(Pk_field)\n",
    "    \n",
    "    k_bins = np.arange(nc//2+1)\n",
    "    k      = 0.5*(k_bins[1:] + k_bins[:-1])*2.0*np.pi/SubBoxSize\n",
    "    \n",
    "    Pk     = np.histogram(dist_3d, bins=k_bins, weights=Pk_field)[0]\n",
    "    Nmodes = np.histogram(dist_3d, bins=k_bins)[0]\n",
    "    Pk     = (Pk/Nmodes)*(SubBoxSize/nc**2)**3\n",
    "    \n",
    "    k = k[1:];  Pk = Pk[1:]\n",
    "    \n",
    "    return k, Pk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "target=np.load('/oldscratch/xz2139/Full_zeros.npy')\n",
    "prediction=np.load('/scratch/xz2139/predictions_1_regression1.pth.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred=prediction[0,448:,448:,448:]\n",
    "tar=target[448:,448:,448:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_k,p_pk0=power_spectrum_np(pred,pred.mean(), 42.1875)\n",
    "t_k,t_pk0=power_spectrum_np(tar,tar.mean(), 42.1875)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAF6CAYAAADh86m8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XmYjfX/x/Hnx4wx1sGQ7GPJMox17GRLoUhFkkiyp117318q7au1Eu2SUopSIvsWg7FkTXayxVgGs31+f3wIJQYzc5855/W4rnON+z73fZ/3zPd7rl7XZzXWWkREREQyuyxeFyAiIiKSFhRqRERExC8o1IiIiIhfUKgRERERv6BQIyIiIn5BoUZERET8gkKNiIiI+AWFGhEREfELCjUiIiLiF4K9LiA9FChQwEZERHhdhoiIiKSBJUuW7LPWFrzQdX4ZaiIiIoiJifG6DBEREUkDxpgtqblO3U8iIiLiFxRqRERExC8o1IiIiIhf8MsxNSIiIuktMTGR7du3c/z4ca9L8RuhoaEUK1aMrFmzXtL9CjUiIiKXYPv27eTOnZuIiAiMMV6Xk+lZa9m/fz/bt2+nVKlSl/QMdT+JiIhcguPHjxMeHq5Ak0aMMYSHh19Wy5dCjYiIyCVSoElbl/v3VKgRERERv6BQIyIiIuTKlQuAnTt30r59+/Ne+/bbbxMfH//3cevWrTl48GC61pcaCjUiIiJ+Kjk5+aLvKVKkCOPHjz/vNf8MNZMnTyZv3rwX/VlpTaFGRETkMhmTPq/z2bx5MxUqVKBz585UrFiR9u3bEx8fT0REBI899hg1atTgq6++YuPGjbRs2ZKaNWvSqFEj1q5dC8CmTZuoV68eUVFRPP3002c9t3LlyoALRQMGDKBy5cpUqVKFoUOHMmTIEHbu3EnTpk1p2rQp4LYn2rdvHwBvvvkmlStXpnLlyrz99tt/P7NixYr07NmTSpUqce2113Ls2LG0/p9BoSbVDh2CF16AM5KpiIiIl9atW0e/fv1Ys2YNefLkYcSIEQCEh4ezdOlSbrvtNnr16sXQoUNZsmQJr7/+Ov369QPg/vvvp2/fvqxcuZLChQuf8/kjR45k8+bNxMbGsmLFCjp37sx9991HkSJFmDFjBjNmzDjr+iVLlvDhhx/y66+/snDhQt5//32WLVsGwIYNG7jnnnv47bffyJs3L19//XWa/z0UalJr4EB4+mmoVAkmTfK6GhER8SHWps/rQooXL06DBg0AuOOOO5g7dy4AHTt2BODIkSPMnz+fDh06UK1aNXr37s2uXbsAmDdvHp06dQKgS5cu53z+tGnT6N27N8HBblm7/Pnzn7eeuXPnctNNN5EzZ05y5crFzTffzJw5cwAoVaoU1apVA6BmzZps3rz5wr/gRfL5xfeMMRWB+4FwYIq1dpQnhXToANOnw/Ll0Latew0eDBERnpQjIiLyzynQp45z5swJQEpKCnnz5iU2NjZV96enbNmy/f3voKAg/+l+MsZ8YIzZY4xZ9Y/zLY0x64wxvxtjHgew1q6x1vYBOgLXeVEvwLZi9ehYJoZ9T78NuXPDxIkQGQkvvggnTnhVloiIBLCtW7eyYMECAD7//HMaNmx41vt58uShVKlSfPXVV4BbtXf58uUANGjQgC+++AKAMWPGnPP5LVq04L333iMpKQmAv/76C4DcuXNz+PDhf13fqFEjvv32W+Lj4zl69CgTJkygUaNGafCbpo5X3U8fAS3PPGGMCQKGA62ASKCTMSby5HttgcnAFxlb5mn/+x98+U0wxV67n1e7ryWpQyc4dgyeegqqVoVffvGqNBERCVDly5dn+PDhVKxYkQMHDtC3b99/XTNmzBhGjx5N1apVqVSpEt999x0AgwcPZvjw4URFRbFjx45zPr9Hjx6UKFGCKlWqULVqVT7//HMAevXqRcuWLf8eKHxKjRo16NatG7Vr16ZOnTr06NGD6tWrp/Fv/d+MTU2nXXp8sDERwPfW2sonj+sBA6211508fgLAWvvSGfdMtNa2vdCzo6OjbUxMTJrWu2sXPPYYfPqpOy5RAj676xcafnEPZt06d7JTJ3jjDfiPAVciIuI/1qxZQ8WKFT37/M2bN3PDDTewatWqC1+ciZzr72qMWWKtjb7Qvb40ULgosO2M4+1AUWNME2PMEGPMSGDmf91sjOlljIkxxsTs3bs3zYsrXBg++QTmzHENM1u3wtXPNueGYsvZ++CLkD07jB0L5cu7sTYnm+pEREQkY/hSqDkna+1Ma+191tpe1to3z3PdSGtttLU2umDBgulWT8OGEBMDw4dD3rww+ZdsFB32BC93XU1i67Zw+DA88ABER8PJfk4REZG0FhER4XetNJfLl0LNDqD4GcfFTp7zOcHB0K8frF8PPXq4Rpkn3osgIvY7Zj08ERsR4WZJ1a/vLji5IJGIiIikH18KNYuBq4wxpYwxIcBtwMSLeYAxpo0xZmRcXFy6FPhPBQvC++/DwoWuYWbnTmjyRhtaFvuNPT2fgqxZYfRo1yU1ahSkpGRIXSIiIoHIqyndY4EFQHljzHZjzN3W2iSgPzAFWAN8aa397WKea62dZK3tFRYWlvZFn0ft2vDrry7ghIfDz3NzUOSDQbzYaSWJjZvDX39Bz57QoAGcXFlRRERE0pYnocZa28laW9ham9VaW8xaO/rk+cnW2nLW2jLW2he8qO1SZcnieprWr3ddU9bCU5+Up9jqqczq+wW2cOHTTTr33w8Z1JokIiISKHyp++myZXT307nkz+8GEcfEuCE1e/YamrzTkRbF1rL79gfcDmVDhkCFCm62lEdT6kVEJHM7ePDg33s9paeZM2cyf/78dP+ctOBXocar7qdzqV4d5s5108ALFYJfFueh8Ni3eOHmJSTWrg9//gm33w7Nm8OaNV6XKyIimczFhhprLSmXMLZToUYA1yjTpQusWwcPPui6qJ7+qipFfp/DrK6jseHhMGOGW/jmySfh6FGvSxYRkUthTPq8zuPxxx9n48aNVKtWjQcffJDmzZtTo0YNoqKi/l41ePPmzZQvX56uXbtSuXJltm3bxujRoylXrhy1a9emZ8+e9O/fH4C9e/dyyy23UKtWLWrVqsW8efPYvHkz7777Lm+99RbVqlX7e3NKn2Wt9btXzZo1rS9atcrapk1P77/avNo+u/vGnqdPlChh7bffWpuS4nWpIiJyAatXrz59kF4bdZ/Hpk2bbKVKlay11iYmJtq4uDhrrbV79+61ZcqUsSkpKXbTpk3WGGMXLFhgrbV2x44dtmTJknb//v02ISHBNmzY0N5zzz3WWms7depk58yZY621dsuWLbZChQrWWmufeeYZ+9prr6Xp3+58zvq7ngTE2FT899/nd+m+GMaYNkCbsmXLel3KOVWq5LaI+vJLePhh+CU2nEKxIxl0Q3ce3dKPrCuXQbt2cMMNbtxNqVJelywiIqnh8fhIay1PPvkks2fPJkuWLOzYsYPdu3cDULJkSerWrQvAokWLaNy4Mfnz5wegQ4cOrF+/HoBp06axevXqv5956NAhjhw5ksG/yeXxq+4n60Njav6LMdCxI6xdC48/7payefr7uhTeuojZ7Ydg8+SB7793O4APGqQdwEVE5ILGjBnD3r17WbJkCbGxsRQqVIjjx48DkDNnzlQ9IyUlhYULFxIbG0tsbCw7duwgV65c6Vl2mvOrUJOZ5MoFL70EK1fCtdfC/rhgGo+/l+ZF17GnRWc4ftxtDR4VBVOnel2uiIj4mNy5c3P48GEA4uLiuOKKK8iaNSszZsxgy5Yt57ynVq1azJo1iwMHDpCUlMTXX3/993vXXnstQ4cO/fs4Njb2X5/j6xRqPFa+PPz0E0yYACVLwow1V1Jo6me8eM10Eq+qCBs2uNRz223wH1vDi4hI4AkPD6dBgwZUrlyZ2NhYYmJiiIqK4pNPPqFChQrnvKdo0aI8+eST1K5dmwYNGhAREcGp3o0hQ4YQExNDlSpViIyM5N133wWgTZs2TJgwIVMMFDbWj9ZJOWNMTc8NGzZ4Xc5Fi4+HV1+Fl192vU75cibwbeO3aDTzOUx8vGveee45uPdetwGViIh4Zs2aNVSsWNHrMi7akSNHyJUrF0lJSdx00010796dm266yeuy/nauv6sxZom1NvpC9/pVS01mGFNzPjlywMCBsHo1tG0LB46G0HjyY7Qospo9DdrBkSPw0ENQsybMm+d1uSIikgkNHDiQatWqUblyZUqVKkW7du28LinN+FVLzSnR0dE2JibG6zIu2+TJcN99sHGjO36x4Q88suVegrdtcifuugteecXtrCkiIhkqs7bU+Dq11Pip1q1h1Sp44QXInh2enHs9hfb9xtym/8OGhMCHH7pBOe+9px3ARUQ84I8NA1663L+nQo2PCw11iw2vXQvt28Nfx7LTaMZztLhyFXtrXAsHDkCfPlCvHixZ4nW5IiIBIzQ0lP379yvYpBFrLfv37yc0NPSSn6Hup0xm2jQ3TnjtWgDLK9HjeWj7gwT/ucPtw9CvHzz/POTN63WpIiJ+LTExke3bt/+9HoxcvtDQUIoVK0bWrFnPOp/a7ie/CjWZffZTaiUkwNChblDxkSMQHnKYiTWfpd6itzHJyW4Hzddfh86dL7h3iIiIiK8LyDE1mX32U2qFhLhtFtatc7llf0JuGix4nVaFlrG/YgPYvdvtpNm0KZxcPElERMTf+VWoCTRFisBnn8GsWVClCkzZGUXBNbN5o/KHJOUr4N6oXh2aNIHx4yEx0euSRURE0o1CjR+4+mo3RnjIEMgTloUBq7pR9PA65kY/gM2Vy4WbDh3cBpmDBrmWHBERET+jUOMngoPdAOL166F7d9iTlJ9GMW9RKusOJl47jMSyFdw2C//7HxQvDnfcAQsWeL6zrIiISFpRqPEzV1wBo0e7vFKnDmw5kIcbf76H7H+s5rmrp7G/UTtscjKMGQP160N0tFvv5tgxr0sXERG5LJr95MeshYULXbfU+PGQlOTOt660hZdKvkvUr+9j9u93J/Pnhx49oG9fiIjwrGYREZF/Csgp3af48zo1l2rHDnjnHbf48L597lzR8OMMrj+ONluGErLi5MJ9xkCbNtC/P1xzjaaEi4iI5xRqFGrO6fhx+OIL13qzbJk7FxxkeazpIu4LGkbBGV9iEhLcG+XLwz33wJ13Qp483hUtIiIBTaFGoea8rHUbfQ8ZAt98A8nJ7nyLKrt55apRVFv4DmbHDncyVy7o2tUFnMhI74oWEZGApFCjUJNq27a5rqmRI+HUEJvCBZN4s8l33LRjGNnmzzx9cbNmrmuqTRs35UpERCSdKdQo1Fy0Y8dg7FgYPBhWrHDngoPhoWtX8VC24RT6+VM4etS9Uby420izRw835UpERCSdKNQo1Fwya2HOHBduvv0WUlLc+aY14ni10sfUXDAM8/vJ2WUhIdCxo2u9qV3bu6JFRMRvKdQo1KSJLVtgxAh4/304cMCdK1wohVdaTOPWPcPINvX70wv41arlws2tt8JlbB0vIiJyJoUahZo0FR/v1usbMgRWrXLnsmaF/tdv4pFc71D4h1GnU0+RIvDRR9CihWf1ioiI/wjIXbqNMW2MMSPj4uK8LsXv5MgBPXu6sTbTp0O7dm4xv7e+LUWRz16l6VXb+bXXaGxUFdi5E6691m0lfuKE16WLiEiAUEuNXLJNm1zX1KhRcPCgO1escDJjol6m0S/PYJKToVo1+PxzqFjR22JFRCTTCsiWGslYpUrBa6/B9u1uSnjFirB9VxCNf36KLqXmcaJ4GYiNhRo13AV+GKBFRMR3KNTIZcuZ083u/u03mDzZhZ0xv9eh8K5lLK/RzS1j3K8f3Hgj7N3rdbkiIuKnFGokzRgDrVq5cTd9+sCBpNxUW/ohT5YZR3KevDBpElSpAj//7HWpIiLihxRqJM3lyuV6m6ZMgWLF4KWNt1LxxHK2l74a/vwTrrsOHnzQteCIiIikEYUaSTfXXgsrV7r9MDecKEHJP6YzMuJFbHAwvP021Knj+qxERETSgEKNpKu8ed2SNd9+CwWuCKL35idoHDyfuIJlXT9VdDQMH65BxCIictkUaiRD3Hija5S59VaYc7wWRfcu46ei3V0X1KkNMvfs8bpMERHJxBRqJMMUKADjxsEXX0C2/LlotWM0XXN8xYmc+eCHH9wg4p9+8rpMERHJpBRqJMN17OhabW64AT6Nb0/Zo8tZWaAJ7N7tpk/df78GEYuIyEXzq1CjbRIyjyuvhIkT4cMP4VCe4lTbN41nc7xMSlCw22CqVq3Tm0yJiIikgl+FGmvtJGttr7CwMK9LkVQwBrp1czOkml0TxMD4x6iVvJCducu5QBMdDUOHahCxiIikil+FGsmcSpRw6/GNGAFrc9TkqsNLGZO9h9sM87774PrrXdeUiIjIeSjUiE8wBvr2dbO8azTMyR3H3udmvuZItvzw448QFeX2YBAREfkPCjXiU8qUgZkz4Y03YHK2myl/YgXzsjVze0Zdfz088ggkJnpdpoiI+CCFGvE5QUHw0EOwbBkUrVWURiem8hgvk2yC4PXXoVEj2LzZ6zJFRMTHKNSIz6pYEebPh+cHZeGtrI/RyM5mV9bi8OuvUL26W6ZYRETkJIUa8WnBwfDUU7BwIewvV5/Kicv4PksbOHgQbroJHnjADSgWEZGAp1AjmUKNGrBkCbTrHk6blO94kDdJNFlh8GBo0AA2bvS6RBER8ZhCjWQauXLB6NEwbpzhw7AHaWDnsjUowqWdGjXgq6+8LlFERDykUCOZzq23QmwsZK1fmyrJy/iam+HQIfdGv37aYkFEJEAp1EimFBEBs2bB/f+Xl1vNePozlAQTAu+8A3Xrwvr1XpcoIiIZTKFGMq3gYHj2WZg5yzCxeH/q2gX8YcrA8uVQsyZ8/rnXJYqISAZSqJFMr1Ejl2NK31KDanYpX9ARjhyBzp2hZ0+Ij/e6RBERyQAKNeIX8uVz44TfGJmH7qFj6cV7nDDZYNQoqFMH1qzxukQREUlnCjXiN4xxDTNLlhoWVe1FLbuItZSHVauw0dHw8cdelygiIulIoUb8TsWKbrG+ZvdXIZoYPuUOTHw8dOsGXbq4N5OSvC5TRETSmM+HGmNMO2PM+8aYccaYa72uRzKH0FB4+2348odcPFzgE+7iA+LJDp99BvXquf6q1q3htdcgJgaSk70uWURELpOx1mb8hxrzAXADsMdaW/mM8y2BwUAQMMpa+/IZ7+UDXrfW3n2h50dHR9uYmJi0L1wypV274M47YdvUNdzPYFoEzaBM8j+mfIeFQePG0LSpe0VFQRafz/wiIgHBGLPEWht9wes8CjVXA0eAT06FGmNMELAeaAFsBxYDnay1q0++/wYwxlq79ELPV6iRf0pJgbfegldegb17oQg7aMJMmjKDFsEzKJn0x9k3hIefHXIiI92gHRERyXA+HWoAjDERwPdnhJp6wEBr7XUnj584eenLJ19TrbXTzvO8XkAvgBIlStTcsmVL+hUvmZa18NtvMGOGe82cCQcOQHG20pQZNGM6LYJnUCRp29k3FizoxuMMGgTZs3tSu4hIoMqMoaY90NJa2+PkcRegDq715k5cy02stfbdCz1bLTWSWikpbo2bUyFn9mw4dMhSmj9oygzXkhM0gyuSd7kbKlaEMWOgenVvCxcRCSCpDTU+P2jAWjvEWlvTWtsnNYFG5GJkyeLyyUMPwaRJsH8/LFpk6P1KGbZf14PeOcdQKHkHdVnAGirAmjXYOnVcP5YGF4uI+BRfCjU7gOJnHBc7eU4kwwQHQ61a8Oij8NNPrmtq3jxD40frUjd4CcO4B5OYCI8/Ds2bw9atXpcsIiIn+VKoWQxcZYwpZYwJAW4DJl7MA4wxbYwxI+Pi4tKlQAk8WbNC/fquYWbeshx8UmsYrfmBPykEs2aRElUFxo71ukwREcGjUGOMGQssAMobY7YbY+621iYB/YEpwBrgS2vtbxfzXGvtJGttr7CwsLQvWgJe5cowfz40fbU1tUJWMJE2ZDkUB7ffDnfcAQcPel2iiEhA82ygcHrSQGFJb+vWwd3dLZHz3+ctHiQn8SQXK0HQmE/h6qu9Lk9ExK/4zUDhi6HuJ8ko5cvD7DmGyoN7UT90GYuJJmj7VmyTJtjHn4CEBK9LFBEJOH4VatT9JBkpSxa47z6Y8Fs5nmw8n0E8RYo1mFdeJjG6Hqxd63WJIiIBxa9CjYgXSpeGn2dk5Yr3BtEq+yw2EUHWlUtJqloD+8U4r8sTEQkYCjUiacAY6NULRq9ryCPXxPIpdxCccAzT6TaSn3nOLWUsIiLpyq9CjcbUiNeKF4evfg4jafQnPBr8JikYgp57hoQOneH4ca/LExHxa34VajSmRnyBMXBXd8Mtcx+ka9hEDpOLkK/HcrxeU9i92+vyRET8ll+FGhFfUqcOvLD8BrqVnccWShAau5DjVWvDihVelyYi4pcUakTSUcmS8OGSKjzedBELqEvo7q0k1mkA33/vdWkiIn7Hr0KNxtSIL8qTBz79uRDjes/gczqR9fgRbJu22Dfe1ABiEZE05FehRmNqxFcFB8Pb74ayf/AYnjHPYrCYAQ+T1KO3FuoTEUkjfhVqRHzdvfcZ6nz/f3TNNo5jhBL8wfskNGsJf/3ldWkiIplesNcFiASa1q2h2KJbubVFBO/vacuV82aQVC6S4No1oFSps18REZAvn5tSJSIi56UNLUU8smsX9Gy5jedW3EgNlv3ndUk582BLliL4qlKYLnfALbdkYJUiIt5L7YaWCjUiHoqPh769klkx9jdKpGyiFP9+5eLoWfckDnmHrPf28ahiEZGMF5ChxhjTBmhTtmzZnhs2bPC6HJFUS0pyLTfbtsHWrWf83Go5vGkfQVs3UT9uMgN5FoDjb44g9MG+HlctIpIxAjLUnKKWGvFHK1bAlw0HM+jwAwAcfXU4OR/p53FVIiLpL7WhRrOfRDKJKlWg29L7GZhvMAA5H72HuBeHe1yViIjvUKgRyUTKloUeK+5j0BVDAAh7qj/7nx3mcVUiIr5BoUYkkylWDHqvupdXig0FIHzgvfz51FCPqxIR8Z5CjUgmVLAg9FnVn7fKuFaaK1+8j22PDPG4KhERb/lVqNHeTxJIwsKg94p7GBY5AoDir9/PxvsGe1yViIh3/CrUaO8nCTQ5ckCvZX0ZVfMdAMoMfYDV3V/zuCoREW/4VagRCUQhIXDXr334rOG7AER++CjLWz6mHcBFJOAo1Ij4gaAg6Dy7N+PbfUYiwVSd8ipLavTAJiZ5XZqISIZRqBHxE8ZA+wmdmXbvROLJTs3YD1hatgOJh497XZqISIZQqBHxM62GtGLZq9M4QF5qbv2W1SVbcXjHIa/LEhFJdwo1In6owSP12T5mNn9mKUzVAzPZcVUT9qza43VZIiLpSqFGxE9F3R7F8Wnz2Jy1LBWOLeNo9QZs/GWz12WJiKQbvwo1WqdG5GwRTUuRa9lc1uWoRqmk38nRoj7fN3qFuW/HcOhAstfliYikKe3SLRIA4nfFsTHqRqL2z/r73F/kY0W+Jhyq1ZwCHZtTrWN5cuQ0HlYpInJuqd2lW6FGJECkHE9gzQvfcGzSNIqs/YUiJzaf9f52irK49UDafd8Do2wjIj5EoUahRuS8jqzcxKZRv5D08y+U+H064UluIPEX0a/TfsHDBAd7XKCIyEmpDTV+NaZGRFIvV1Qpogb3oPqasYQn/Mnq/m4PqdtiBjC28gvEx3tcoIjIRVKoEREwhsihfdn49IekYOiy7mm+rPB/HPjL/1pyRcR/KdSIyN/KPN+NXa9+RhJBdNv2PBPKP87OHQo2IpI5KNSIyFmKPnI7B0Z8QSLBdN/3Kj9VfJAN6xVsRMT3KdSIyL8U7NueY59+TYIJofvhwcyr1o+lMSlelyUicl4KNSJyTnnuaEvy19+RkCUb3Y69y8Z6dzBzygmvyxIR+U8KNSLyn7Lf1BIzeTLHgnPRIWksKa2uZ+Jn2hxTRHyTQo2InFfW65qRbeFsDuUoRDP7C8W7NObT13Z5XZaIyL/4VajR3k8i6SNLzerkXrmA/eFXUZ1YGj5an3ceXI8frt0pIpmYX4Uaa+0ka22vsLAwr0sR8TumdCnC18xjd0RtSrGZDm/X5+1Ov5Ki8cMi4iP8KtSISDorWJBCq6azq3prCrCf3uOa8tJVHzBnlpKNiHhPoUZELk7OnBRe9B07W3YnB8d46o+7ydakLg83WMjKlV4XJyKBTKFGRC5ecDBFJo/i2PufcSh3EWqzmDfm1yO2Shfu77CTLVu8LlBEApFCjYhcGmPI3qMzeXau4+gDT5EYlI0ufMYL48sxqsxLjPk4yesKRSTAKNSIyOXJlYucbw0i6/rVHL32JnJxlOeTnyS0W0eGvJ7gdXUiEkAUakQkbZQuTc4p38CUKRzPnpdb+IaIR9rz9CMnNPVbRDKEQo2IpK1rryV07i+cyJmftkyi4es30rfbMZLUGyUi6UyhRkTSXo0aZJs/gxNhBWnJFDp8cgOd2x3l2DGvCxMRf6ZQIyLpo0oVss2fSUL4lTRnOv1+aE2Hloc5ftzrwkTEXynUiEj6iYwkZP4sEgsVpTGzeXx2K7q1P0JioteFiYg/UqgRkfRVrhxZ580i4criNGQevX5oS++ux7S9goikOYUaEUl/ZcoQMvsXEgoUphkzuPWLm3i4v2ZFiUjaUqgRkYxx1VWEzJpGQl43eLjxOx159mn1Q4lI2lGoEZGMExlJyMypJOTKRzu+I/LFzrz2UpJabEQkTSjUiEjGqlqVkBk/k5A9D7fyFUWf7Ervbic0K0pELpvPhxpjTGljzGhjzHivaxGRNBIdTci0H0nMlpPbGcudnzSjbZ3dbN3qdWEikpl5EmqMMR8YY/YYY1b943xLY8w6Y8zvxpjHAay1f1hr7/aiThFJR/Xrk3XBHBKuLE4D5jNqRS3urBrL9OleFyYimZVXLTUfAS3PPGGMCQKGA62ASKCTMSYy40sTkQxTvTohyxaRFF2XEmzj+4MNGHHNN4wY4XVhIpIZeRJqrLWzgb/+cbo28PvJlpkE4AvgxtQ+0xjTyxgTY4yJ2bt3bxpWKyLGY9O7AAAgAElEQVTp6sorCZ4zg5QuXclJPOPtLWy85w1GjfK6MBHJbHxpTE1RYNsZx9uBosaYcGPMu0B1Y8wT/3WztXaktTbaWhtdsGDB9K5VRNJSaChZPv4IXnsNawxvMIBfe45i3DivCxORzCTY6wIuxFq7H+jjdR0iks6MgQEDMDlzQr9+vEtvOt4eTs6cN3HDDV4XJyKZgS+11OwAip9xXOzkuVQzxrQxxoyMi4tL08JEJAP17Yt9ZiBBpPBZSieG3DyTGTO8LkpEMgNfCjWLgauMMaWMMSHAbcDEi3mAtXaStbZXWFhYuhQoIhnDPPN/2L79COUE4xPb8lTrZSxZ4nVVIuLrLhhqjDHVjTGvGWN+NcbsNsbsMsYsPHmu+qV8qDFmLLAAKG+M2W6MudtamwT0B6YAa4AvrbW/XcrzRSSTMwYzdAi2w63k4TATjrekb6vNbN/udWEi4suMPc/65MaYybhZSpOARcBOwACFcbOV2gJ5rbXXp3+pqRcdHW1jYmK8LkNELteJE6S0voEs06exhBr0jZrH9Pmh5MrldWEikpGMMUustdEXvO4CoeYKa+2eC3xQQWutT8yhNsa0AdqULVu254YNG7wuR0TSwoEDJNeIJmjzH4zibr67YRTffgtBQV4XJiIZJbWh5rzdT6cCjTGm1Tk+oM/Ja3wi0IDG1Ij4pXz5CJrwNSnZQunBaAp9P4oBA7wuSkR8UWoHCv/PGNPs1IEx5lEuYmE8EZHLUq0aWUa+B8Aw+jP77SUMG+ZxTSLic1IbatoCLxpjGhljXgDqoFAjIhmpa1fo04dQTvA1tzDw3v2MHet1USLiS1IVaqy1+3DBZjhQBGh/cisDn6J1akT83NtvQ+3aRLCFSdzAvV0OMnmy10WJiK+40EDhw8CZF4QASSfPWWttnvQt79Jo9pOIH9u2DRo2hK1biaEm7UKn8MXUcBo29LowEUkvaTVQOLe1Ns8Zr1Brba5T59OuXBGRVCpeHObMwZYpQzRLmHy8KXe13k1srNeFiYjXzhtqjDGlL/SA1FwjIpKmSpTAzJ6NrVCBKqxk4uEm3NF0B0uXel2YiHjpQmNqXjTGTDLG9DLG1DDGFDbGFDn5717GmB+AFzKi0NTQmBqRAFKkCGbWLFKiqlCRtUw42ITbGu9i/nyvCxMRr5x3TA2AMaYMbh+mhkCJk6e3APOAsdbaP9K1wkugMTUiAeSvv0i5pgVZli1lNRW5PsdMRk+6gmbNLnyriGQOaTKmBsBauxF4E5gGrAfWAVOBN3wx0IhIgMmfnyxTf8ZWqkwka/g2vgWdW/3FTz95XZiIZLTUrlPzMVARGAIMBSKBT9KrKBGRixIejvllGrZ8eaqygokJ19Hz1jg2b/a6MBHJSKkNNZWttT2stTNOvnoCldOzMBGRi1KoEOaXX7ClS1OLGD4+fBNdOiWRlOR1YSKSUVIbapYaY+qeOjDG1AE0aEVEfEvRopjp00m5ohDNmEHrhf/j2We9LkpEMkpqQ01NYL4xZrMxZjOwAKhljFlpjFmRbtVdJM1+EhFKliTLl+OwWYJ4gpdZMWgiM2d6XZSIZIQLzn4CMMaUPN/71totaVZRGtDsJxHhtdfg0Uc5SBitr1jCxN/KUKCA10WJyKVIs9lP4ELL+V6XX66ISBobMICUG9uRlzhG7LmF65sd488/vS5KRNJTarufREQyF2PI8vFHJEaUpRrLeWFlG66rF8fGjV4XJiLpRaFGRPxXWBhZJ00gpcAVXMMvfLa5Ie3rbGPZMq8LE5H0oFAjIv6tcmWyLFpIcrkKRLGKH/bX4d6Gy5gxw+vCRCSt+VWo0ewnETmnUqUIWjiflEaNKcIufoy/muevncM333hdmIikJb8KNdbaSdbaXmFhYV6XIiK+Jl8+skydgu14G7k5wsSkVgxpP5uRI70uTETSil+FGhGR88qWDTPmM2yXLuTiKD/YVnzeeyZDhnhdmIikBYUaEQksQUGYDz+Ebt3ISTyTac03D8xm+nSvCxORy6VQIyKBJygIRo+G7t3JwTG+srfwYPttbN3qdWEicjkUakQkMGXJAiNHYltcS0H28d6BDnRsd4Ljx70uTEQulUKNiASuoCDM2M9JLlaCuvxKl2UP0qMHJCZ6XZiIXAqFGhEJbOHhBH0znpSsIfTjHbKM+YTrroN9+7wuTEQull+FGq1TIyKXpFYtsgwbCsAHdCdqxmCia1piYz2uS0Quil+FGq1TIyKXrGdPeOopgklmMA/w3NY7aV7/GIsXe12YiKSWX4UaEZFLZgwMGgTjxmFz5KArn/LzsYbc1WYfO3Z4XZyIpIZCjYjImW69FbNwIbZUaWqylA92t6bTDYeJj/e6MBG5EIUaEZF/iorCzJ1DcokIarOYZ2Lb0eOO46SkeF2YiJyPQo2IyLkUKULQL1NJCi9Ec6bTfsLtPPFIktdVich5KNSIiPyXsmUJnv4ziTnDuJkJlH+zF6+9ar2uSkT+g0KNiMj5VKlC1ik/kBSSne58iHnsET78QMFGxBcp1IiIXEiDBgR/9w3JWYIZwBv83uNlJk70uigR+SeFGhGR1GjZkqAxn5KC4QX7JItveZk5c7wuSkTOpFAjIpJat92Gee89UjA8n/QEi1s8wYrl6ooS8RUKNSIiF8H06gkff0KyCeKhEy+ztN49rIxN9rosEcHPQo32fhKRjJCl6x0kf/kNCVmy0e3YO2yvdRNzfjzidVkiAc+vQo32fhKRjBLSvi12ys8cyZqPVkmTyHN9QyYO3+Z1WSIBza9CjYhIRsp2zdVkX76QPXmvoqpdTvX+9Rn70mavyxIJWAo1IiKXIahiOQr+vpDtJetTnO3UfrI57z+70+uyRAKSQo2IyGUy4fkptnwyu4vXpAx/UH9gC958cp/XZYkEHIUaEZG0EBZGoaU/caBIJJVYTYuXmvLyvTuwmvEtkmEUakRE0kqBAuRbPJW4IhWIYhW3D6vHS13XaHdvkQyiUCMikpaKFCFsxVz+Kl+PEmyj92cNee3GuSRrKRuRdKdQIyKS1sLDyb90GnvqtiGcv3jg++a81+gzEhO9LkzEvynUiIikhxw5uGLON+y45V6ykUC/BV0YX/kZjhzWIBuR9KJQIyKSXoKDKTp+CJsHDCWZLHRa/xxzSnRm++/Hva5MxC8p1IiIpLOI1/rz58hJHDG5aHVwLH9GNmP9gv1elyXidxRqREQyQNGerUmaOY/d2YoTnbiA5CbN2L92r9dlifgVhRoRkQyS9+oq5F61kM3ZylMxYQVxNZtxYtser8sS8RsKNSIiGShH2SKELpjBhuAKlI5fxd6opiTv3O11WSJ+QaFGRCSDXVm9MMcmz2SNiaRY3Gr+jGxK8o4/vS5LJNNTqBER8UCVFoU4OGEGq0wURePWsK9iI5LmLPC6LJFMzedDjTEmpzHmY2PM+8aYzl7XIyKSVurdeAVxE6az3FSj0OHfyXJ1AxJ63QNxcV6XJpIpeRJqjDEfGGP2GGNW/eN8S2PMOmPM78aYx0+evhkYb63tCbTN8GJFRNJRgxsLcHTqfN4MfZJkggh5fwRJ5SNhwgSvSxPJdLxqqfkIaHnmCWNMEDAcaAVEAp2MMZFAMWDbycu0e4qI+J36zbPTZsULtCu+lAXUJXj3Trj5Zrj3XrTNt0jqeRJqrLWzgb/+cbo28Lu19g9rbQLwBXAjsB0XbOA89RpjehljYowxMXv3au0HEclcrroKPlkWxatt5tKP4RwnGwwbBo8/rmAjkkq+NKamKKdbZMCFmaLAN8Atxph3gEn/dbO1dqS1NtpaG12wYMH0rVREJB2Eh8M33wVRaVg/bs3yNYkEw6uvwoABcFxbK4hciC+FmnOy1h611t5lre1rrR3jdT0iIunJGLjnHmg17Ho6M4ZkssCbb0LVqjBnjtflifg0Xwo1O4DiZxwXO3ku1YwxbYwxI+M0c0BEMrm+faHwfbfSmFmsNRVh/Xq45hr46SevSxPxWb4UahYDVxljShljQoDbgIkX8wBr7SRrba+wsLB0KVBEJCO98QaU7tKQqnYZ79IHEhLgpptg+nSvSxPxSV5N6R4LLADKG2O2G2PuttYmAf2BKcAa4Etr7W9e1Cci4guCg+Hjj+F/z2ejLyN4j15ubE2bNjB3rtflificYC8+1Frb6T/OTwYmZ3A5IiI+yxh4+mkoWNDQt887hJDAXfEfQevWMHUq1KnjdYkiPsOXup8um8bUiIi/6t0bXnwpCz0YxVjTCQ4fhubN4dNPvS5NxGf4VajRmBoR8WePPQYPDQiii/2Ez4K6wtGj0LUrdOsGR454XZ6I5/wq1IiI+DNj3LI13e4OpkvyR/QNGU1SSHY38KZmTYiN9bpEEU8p1IiIZCLGwHvvQdeuhncTulMlIYYd+Sq7Kd916sBbb0FKitdlinjCr0KNxtSISCAICoKPPoKRI2FTaCRlDyzix5Inp3w/9BDccgucOOF1mSIZzq9CjcbUiEigMAZ69oSZMyF7vuy03vIOT0R+R0refPDtt3DjjXDsmNdlimQovwo1IiKBpk4dt3tC0aLw8uq2tMw2k8R8BWHKFKhfH1au9LpEkQyjUCMikslVqgSLF0PdujB1dxWuTp5JQrFSbuBwdDQMH66dviUg+FWo0ZgaEQlUhQu7rqi2bWHhoUiq2uUc6dzLjbPp399N+9ZO3+Ln/CrUaEyNiASybNlg7FjX67R2R24qzXmPrS9/DtmzwyefQMuWcPCg12WKpBu/CjUiIoEuRw6YNMl1RW3dCpHPd2Luq/NdU86sWVCrFnz3nbqjxC8p1IiI+Jn8+WHGDOjc2S063PTBakx8YoEbfPP779CuHdx5JyQne12qSJpSqBER8UOhoW5bqAEDICkJ2t1fkkE3LyVl8BDIlcu92b27e1PET/hVqNFAYRGR005tq/Dqq+74f8+H0GbKvRz68kfImdONs2nYENau9bZQkTTiV6FGA4VFRM5mDDzyCPz4o+uWmjwZqt3TkI3Df3KL2/z6K1SrBm+8oe4oyfT8KtSIiMi5XXcdLFni9r3ctAnqP9qQzd+vgrvuclsqDBgAzZrBtm1elypyyRRqREQCREQEzJ0L11wDe/bAdR3z8tvDH8APP8CVV8Ls2RAVBS+9BPHxXpcrctEUakREAkhoKHz9NVSt6jb2joqCnhNac3TBCrdyX1wcPPmke2P6dK/LFbkoCjUiIgEmTx745Re30HBwMIwaBfVvLMimt7+DqVOhShX44w9o3txtBy6SSSjUiIgEoPBwGDoUli2DcuVgxQpo2hS2lb/GbST11FPuwrvvhmefhQ0bvC1YJBX8KtRoSreIyMWpVAkWLXIrEG/ZAo0bwxffhJA0cJALMykpMHCgSz433+wW7xPxUcb64VLZ0dHRNiYmxusyREQyjYMHXW/T0qXuuG5d+OZrS+Flk2HcOPjqK7chZqFCbufMChU8rVcCizFmibU2+kLX+VVLjYiIXJq8eWHePHjvPShWDBYuhOhahsVXXO8W6Ts1xmb3brdg37BhbgdwER+iUCMiIoCbGdWrl1vPplEj2LnT/fzsM9yGmBMnup2+9++He++FihXdOREfoVAjIiJnueIKmDYN+vRx6/J16eJWJU7OlsMtSTxhgut++uMPuPFGeO45SEz0umwRhRoREfm3kBB45x33Cg6G11+H66+HAweN2+V75Up45RW3D8Mzz7iQM2mS12VLgFOoERGR/9Snj1vTpkABmDLFrUZ86BAu6Tz6qFuNuFw512rTti306AGHD3tdtgQohRoRETmvq6+GmBgoU8bNjmrVCtatO/lmq1bw229uQ8xs2WD0aLca8bhx4Ieza8W3+VWo0To1IiLpo2RJt9hwkSIwf75b32bgwJMbewcHw0MPuRHG1au7BW9uuw1q13ZTqkQyiF+FGmvtJGttr7CwMK9LERHxO6VKudzSs6dbk+/ZZ904m7/3vjy1kt/IkW62VEyMm/7dty8cO+Zp7RIY/CrUiIhI+rrySpdZpk51s6SmTIHWrWH79pMXBAe71LNhAzz9NGTNCu++CzVquLVtFG4kHSnUiIjIRWveHGbNcg0ys2a5scIjRpwxjCZnTnj+ebePVLlysHbt6bVtvv5a420kXSjUiIjIJalQAX79FTp0cA0w99zjGmmSks64qGpVWL4cvvjC/XvLFmjfHlq0cAOMRdKQQo2IiFyy4sXhyy9hzBjInt1NfrrlFjhrvkZoKHTs6MbYjBgB+fO7eeJVq7pdwKdM8ax+8S8KNSIictluvx2mT4d8+dzOCVWrui2jztoeKjjYDRpevx769XNdUB984LZeuPdeWLHCs/rFPyjUiIhImqhb122EWbOm62W68043PvhfWSU8HIYPd6sSP/WUG0w8bJhLQl27ur2lRC6BQo2IiKSZcuVgwQLXDVW2rBs207ChW7TvXyIjYdAgmDHDdU9lzw6ffuoWxXnooTOmVImkjkKNiIikqaxZoXt3iI11Y4IPH3azpf73Pzh69Bw3NGjgBhLHxLiuqKNH4a23oHRpt+3C+vUZ/jtI5qRQIyIi6SJnTjeAuG1bOHjQNco0aQKbNv3HDZGR8OOPrlnn1lvdcsWjR7tpVh06/Edzj8hpCjUiIpJuQkLg229h5ky3InFMjOuW6tXr5MaY51K9uts7au1aN0c8a1YYP94N1qldG156CXbtyshfQzIJvwo12vtJRMT3GAONG7uxNl26QFAQvP++Gxc8c+Z5brzqKrd88R9/wMMPQ65cbjG/J5907z3yiHtP5CRj/XBVx+joaBsTE+N1GSIicg6rV7twc6o3qW1bl10KFbrAjfHx8PPPbhr4pEnunDFuIb9ataBPHyhWLF1rF28YY5ZYa6MvdJ1ftdSIiIjvi4x0U7+fe86Nu5k40U39nj//AjfmyAHt2rkbFi1y079DQlzQeeEFN/bm/vtdakpJyZDfRXyLWmpERMQzO3e62dxz57q1+f7v/1yvUmhoKh+wbx/88AN8840LO6dceSVcd517eKtW6VK7ZJzUttQo1IiIiKcSE+Gxx9wsbnC7f/fvD/fdB2FhF/Gg2FgYPBh++gn+/PP0+erV4eqr3ejkyMg0rV0yhkKNQo2ISKYyfToMGADLlrnjiAgYO9atVHxRrHUDdn78Ed54w80nP6VWLWjTBrp1cxtXSaagUKNQIyKS6VjrZkQ98ggsWeK6pAYMcC05efNewgOPH4fZs9288g8/dMenVK4MN93kWnA0wNinKdQo1IiIZFoJCW7m9htvuOMSJVz3VJs2btmaS3LggBuN/MEHbnDxkSPufFCQCzdNmrip4o0bQ7ZsafFrSBpRqFGoERHJ9BYudGNrFi92x4ULw113uTE3hQtfxoMTElwLzvvvw9dfu9WLT8md2y2i06QJNG0KUVFQsODl/BpymRRqFGpERPxCYiKMGAHvvusWGQaXO557zoWb4ODL/IAdO1y4Wb/e7cD5z23Fg4OhShW3VUNUFFSrBkWLXuaHysVQqFGoERHxK9a6qd+vvgrff+/OVakC77wD9eun4Qdt3uzWwfn2W/fvhQvdh58pIsJ1U9WqBc2auTVyjEnDIuRMCjUKNSIifmvSJNcttXmzO+7WDXr3voSZUqmxfz9Mm+a6q9audfs9HDt29jWlS0OZMm4O+m23uUFAlStD9uzpUFDgUahRqBER8Wvx8fDii67lJjHRnevc2Z0rUSIdPzghASZPdtPGf/8dpk51iwD+U758rruqTh3XpHTddVCggFp0LoFCjUKNiEhAWLcOhg1zk5ri493sqC5d3LTwChUyoIDjx2H5crfg35o1LvDs3u3G6PxTgQKur+zqq13giYhwrTsX3PgqsCnUKNSIiASU9eth4EAYN+701k/t2sGbb0KpUhlcjLVuXM7GjbBhA0yZ4oJPfPy5ry9SxA1ADg93QadUKfcqXdoNSg4KytDyfY1CjUKNiEhA2rABXnsNPv8cjh51eeC229wCflFRHhaWmAjz5rkNN3/7DVauhE2bIC4ODh/+7/vy5HHjdfLlc6GnaFGoXdsFniuvdK+QEL/u1lKoUagREQloO3a4BfzGjDm9DE3Dhm4aeMeO3tZ2lpQUNzZn1Sq3QODmzS7s/PGH+3nmPlb/JTQUoqPdTubh4a6bq3BhF4CKFIFy5dJ5oFH68ptQY4wpDTwFhFlr26fmHoUaERE5ZcsWtzLxqFGnJy1dc40LOAMGQM6c3tZ3Qdu3u2Bz4IAbkLxmDfz6qwtCe/a4MT1JSRd+Tp48bvxOrlyu+SpHDvfL58jhFv4pVswForAw9zp1/al/58kDJ064HUczuDvMJ0KNMeYD4AZgj7W28hnnWwKDgSBglLX25VQ8a7xCjYiIXKoDB9z2T08+6f7bDG5XhK5d3SsTN2S4gLNpk2uS2rXLbeK5c6drrtq5022kdWpbiMtljNtGIm9eF3Ry53Y/U1Lc9hMhIWnzOWd9pG+EmquBI8Anp0KNMSYIWA+0ALYDi4FOuIDz0j8e0d1au+fkfQo1IiJy2fbscYv3vf66a/QA11gxfLibEn7Je0v5suRkOHTIvQ4fdsfHjrmgs3atWzX5r79cGIqLc9ed+fPUvxMS/r1Gz5n27nVdX2nMJ0LNyUIigO/PCDX1gIHW2utOHj8BYK39Z6D553POG2qMMb2AXgAlSpSouWXLljSpX0RE/FNCAnz5pXtNmuTO5csHN9/spoOXK+fXY28vXVKSC0P79rk/4oEDblZXUBA0aJAum4H6cqhpD7S01vY4edwFqGOt7f8f94cDL+BadkZdKPyAWmpERCT1rIWRI2HoUDcp6ZToaGjZEjp1gshI7+qT1IeaLBlRzOWw1u631vax1pZJTaARERG5GMa4LRZWrXKhpmtXN342JgYGDYLq1d36N3v2eF2pXIgXoWYHUPyM42Inz102Y0wbY8zIuLi4tHiciIgEmMhI+PhjN6t62DDXSpOQAM8+69bE69LFbQGVmslGkvG86H4Kxg0Ubo4LM4uB2621v/3XMy6Wup9ERCStzJzppoSf2hkcoHx5eOABaNUKSpb0rLSA4RPdT8aYscACoLwxZrsx5m5rbRLQH5gCrAG+TMtAIyIikpaaNHEDideudasSFyvm9pvq29ct9Nu8uZsqfr5JQZIxfH7xvUuhlhoREUkvJ0646d8//ggzZpxerTg8HHr0cGFHrTdpy2dmP2UkY0wboE3ZsmV7btiwwetyRETEz/31F3z3HYwY4QYWA2TJ4mZO1asHffpk0E7hfi4gQ80paqkREZGMdGpT7qFD3bo3iYmn32vd2r1uv92tgyMXT6FGoUZERDywdy8sX+420hw79vSWDDlywK23Qv36cMcdkD27t3VmJgEZatT9JCIivmTLFjf25ptvYOrU0+eNgWbNXMhp1QqKF//vZ0iAhppT1FIjIiK+ZvVqGDfOzaRatuzs95o1g7vuct1U+fN7U58vU6hRqBERER+1cyeMH+9eixfD8ePufEiIW8G4fXu3uWbhwt7W6SsUahRqREQkEzh4EL74wo2/mTPHDTo+pVw513rTtasLO4FKoUahRkREMpndu2HBAjdFfOFCOHz49HuVKsH117uft9zi9qcKFAEZajRQWERE/EVSEsyfDx99BBMmuBadU3LnhqpV3S7iPXvCFVd4VmaGCMhQc4paakRExJ/Ex7vVixcuhB9++PdA41MBp1kzt/Cfvw02VqhRqBERET+1fTtMnw6ffurG4ZxaC+eURo3cnlUdOrjuqizputNj+lOoUagREZEAcPw4zJ0LU6bA7Nlu4b8zQ07u3HDNNXDttZl3V3GFGoUaEREJQIcOuQX/vvrKhZ3du89+v1w5aNPGteY0apQ5uqoCMtRooLCIiMjZtmyBn35yKxpPngzHjp1+zxg3HqdWLWjaFGrXhjJlvKv1vwRkqDlFLTUiIiL/lpTkxuBMn+5+LlgACQlnX1O1KjRoADVquKBTurQ3tZ5JoUahRkRE5LyOHXPTxqdOdSsbL1589to4ACVKuG6q665zLTlly0JQUMbWqVCjUCMiInJR4uNh0SLXijN3LsybB0ePnn1NjhyuNadOHdeiU726a80xJv3qUqhRqBEREbksSUmwdKlrzZk+HWJjYdu2f18XGgoNG8K776bPmByFGoUaERGRNLd/v2vBmTPHhZxVq+DPP922DTt3Qp48af+ZqQ01wWn/0SIiIuKvwsOhbVv3ArcB5969sHJl+gSai5HJ1xg8mzGmjTFmZFxcnNeliIiIBARj3N5TzZt7XYmfhRpr7SRrba+wsDCvSxEREZEM5lehRkRERAKXQo2IiIj4BYUaERER8QsKNSIiIuIXFGpERETELyjUiIiIiF/wq1CjdWpEREQCl1+FGq1TIyIiErj8KtSIiIhI4FKoEREREb+gUCMiIiJ+wVhrva4hzRlj9gJbLvK2MMBXRhhnRC1p/RmX87xLvfdi7kvrawsA+1L5PH+g74d3z9P3w/fp+5H+zytprS14wSdZa/VywW6k1zVkZC1p/RmX87xLvfdi7kvra4EYr/7/4cVL3w/vnqfvh++/9P3wneep++m0SV4XcIaMqCWtP+Nynnep917Mfel1baDwpb+Jvh9pf5++H5fHl/4mgfb9OItfdj+JpDdjTIy1NtrrOkR8kb4f4hW11IhcmpFeFyDiw/T9EE+opUZERET8glpqRERExC8o1IiIiIhfUKgRERERv6BQIyIiIn5BoUbkMhljchpjPjbGvG+M6ex1PSK+xBhT2hgz2hgz3utaxP8p1IicgzHmA2PMHmPMqn+cb2mMWWeM+d0Y8/jJ0zcD4621PYG2GV6sSAa7mO+HtfYPa+3d3lQqgUahRuTcPgJannnCGBMEDAdaAZFAJ2NMJFAM2Pb/7d3v651zHMfx50sjKWUxrPnxRfqaaFHbLBJSSylJmh93yJ2FSHFDbvAHuGw80CcAAAOhSURBVKGUzR22IptYEbEoP2drC60wW9tqzCahLLM2m7cb55LzPU5z9mXnnJ3v83HrOp/r/bn6XDc+p9f5XFfn05Qd6uMYpUFZRu/zQ+obQ43URVV9CPzc0TwP2Nr88jwArABuAnbSCjbgnNIUcITzQ+obv4Cl3s3i7xUZaIWZWcAq4JYkSxiuPWCkfuo6P5KcmmQpcFmSRwczNE0V0wY9AOlYV1V7gbsHPQ5pGFXVT8DiQY9DU4MrNVLvvgPObvt8VtMmyfmhIWCokXq3AbgwyXlJTgBuA14f8JikYeH80MAZaqQukrwErAXGk+xMck9VHQTuB1YDm4CXq+rLQY5TGgTnh4aVu3RLkqSR4EqNJEkaCYYaSZI0Egw1kiRpJBhqJEnSSDDUSJKkkWCokSRJI8FQI6kvkowl+aJL+6VJnm+On0jy8GGusTTJlUmWJfktyclt555KUklOm8TY5iR57kj7SRouhhpJg/YIsLTH2iuAdc3xVppdoJMcB1zHJP+Wv6o2AhckOWMy/SUNB0ONpL5Lcn6Sz5PMBxYA69tOX5zk/STbkzzQ1mc2sKWqDjVNK4BFzfE1wBrgYFM7luTrJC8m2ZTklSQnNefmJvkkycYk69tWe94Cbj1a9yzp6DPUSOqrJOPAq8BdQIDNNfGvzS8CFgLzgMeTHN+03wC83Va3BZiRZDpwO62Q024ceKaqZgN7gHubPYlWAg9W1RzgemBfU78euPp/uUlJA2GokdRPM4DXgDubRz7nArs7at6sqv1V9SPwA/DXI6GFTAw1AKtobZw4H/io49y3VbWmOX4BuIpW0NldVRsAqmpPs2cRwC5g7D/cm6QBmzboAUiaUn4BvqEVML5q2tJRs7/t+BAwrXl0dEpV7eqoXQl8Ciyvqj+SCZfq3Nju3za6Sw81koaYoUZSPx0AbgZWJ/kV2A6c2UO/a4H3OhurakeSx4B3u/Q5J8mCqloL3AF8DGwGZiaZW1Ubmvdp9jWrNTOBHZO6K0lDwcdPkvqqqvYCNwIPAacD4+lYYumi832a9us9W1XbupzaDNyXZBMwHVhSVQdovVz8dJKNwDvAiU39PP75CEvSMSQT38+TpP5KspxW4Fh3mJrPgPlV9XuP1xwD3qiqS45gHB8Ai6rq+177SBourtRIGrQngcWHK6iqy3sNNJORZA6wzUAjHdtcqZEkSSPBlRpJkjQSDDWSJGkkGGokSdJIMNRIkqSRYKiRJEkj4U8k5MTxt1aldwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b5e0df4fd30>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[9,6])\n",
    "\n",
    "plt.plot(p_k,p_pk0,label='prediction',linewidth=2,color='b')\n",
    "plt.plot(t_k,t_pk0,label='target',linewidth=2,color='r')\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('k(h/Mpc)')\n",
    "plt.ylabel('p(k)')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}