2D_Strang_SSP_RK3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code used in Chapter 4: DFT Numerics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import some important libraries. Most Python distributions do not come with pyFFTW, might need to pip install pyfftw."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numexpr as ne\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "import time\n",
    "import pyfftw\n",
    "import h5py\n",
    "import os\n",
    "from numba import jit\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define functions for central difference based on the 4th-order isotropic central difference stencil.\n",
    "<img src=\"grid.png\" alt=\"central difference grid\" style=\"width: 500px;\"/>\n",
    "\n",
    "Function `shift` gets the stencil values around each $(m,n)$ and stores them in an intermediate array. These intermediate arrays are multiplied by the corresponding weights in function `centralDiffNumexpr` (`centralDiffShift` is a non-numexpr equivalent)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def shift(A,num):\n",
    "    E = np.empty_like(A)\n",
    "    F = np.empty_like(A)\n",
    "    G = np.empty_like(A)\n",
    "    H = np.empty_like(A)\n",
    "    \n",
    "    E[:num] = A[-num:]\n",
    "    E[num:] = A[:-num]\n",
    "    F[:-num] = A[num:]\n",
    "    F[-num:] = A[:num]\n",
    "    \n",
    "    G[:,:num] = A[:,-num:]\n",
    "    G[:,num:] = A[:,:-num]\n",
    "    H[:,:-num] = A[:,num:]\n",
    "    H[:,-num:] = A[:,:num]\n",
    "    return E + F + G + H\n",
    "\n",
    "@jit(nopython=True)\n",
    "def centralDiffShift(A,h,Nx):\n",
    "    A = A**3\n",
    "#     D = np.empty_like((Nx,Nx))\n",
    "    D = -1.*shift(A,2)+ 16.*shift(A,1) - 60.*A\n",
    "    return D/(12.*h**2)\n",
    "\n",
    "# def centralDiffNumexpr(A,h,Nx):\n",
    "#     A = ne.evaluate('A**3')\n",
    "#     F = shift(A,2)\n",
    "#     E = shift(A,1)\n",
    "#     D = np.empty_like(A)\n",
    "#     D = ne.evaluate('(-1.*F + 16.*E - 60.*A) / (12.*h**2)')\n",
    "#     return D\n",
    "\n",
    "class writer(object):\n",
    "    \n",
    "    def __init__(self,prefix,chiN,eps,sigma):\n",
    "        self.prefix = prefix\n",
    "        self.chiN = chiN\n",
    "        self.eps = eps\n",
    "        self.sigma = sigma\n",
    "        \n",
    "        self.path = self.get_file_path()\n",
    "\n",
    "        self.groups = ['avg_conc', 'energy', 'order_parameter']\n",
    "        \n",
    "        \n",
    "    def get_file_path(self):\n",
    "        base_output_folder = \"./output\"\n",
    "        base_output_name = \"/%s_chiN=%.2f_eps=%.4f_sigma=%.4f\" %(self.prefix,self.chiN,self.eps,self.sigma)\n",
    "        hdf5_format = '.h5'\n",
    "        return base_output_folder + base_output_name + hdf5_format\n",
    "        \n",
    "        \n",
    "    def create_output_file(self, options):\n",
    "        print(self.path)\n",
    "\n",
    "        if os.path.exists(self.path):\n",
    "            os.rename(self.path, self.path+'_old')\n",
    "\n",
    "        file = h5py.File(self.path, 'a')\n",
    "\n",
    "        for group in self.groups:\n",
    "            # check if groups have been created\n",
    "            # if not created, create empty groups\n",
    "            if not (group in file):\n",
    "                file.create_group(group,track_order=True)\n",
    "\n",
    "        for key in options:\n",
    "            if type(options[key]) != str:\n",
    "                file.attrs.create(key,options[key])\n",
    "        file.close()\n",
    "        \n",
    "    def write_data(self, it, data):\n",
    "        file = h5py.File(self.path , 'r+')\n",
    "        \n",
    "        for group, datum in zip(self.groups,data):\n",
    "            name = group + '_%i' %it\n",
    "            if isinstance(datum, (list, tuple, np.ndarray)):\n",
    "                file.create_dataset(group + '/' + name, data=datum ,chunks=True, compression='gzip', compression_opts=4, dtype=np.float32)\n",
    "            else:\n",
    "                file.create_dataset(group + '/' + name, data=datum , dtype=np.float32)\n",
    "        # print(\"writing time = %.1f for arrays %s\" %(name,path))\n",
    "        file.close()\n",
    "        \n",
    "    def get_dataset(self,group,it):\n",
    "        file = h5py.File(self.path, 'r')\n",
    "        \n",
    "        arr = np.copy(file[group][group + '_' + str(it)])\n",
    "        \n",
    "        file.close()\n",
    "        return arr\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialise wisdom to None for pyFFTW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "wisdom = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main looping consisting of the Strang-Splitting involving:\n",
    "\n",
    "    Substep 1: Fourier-PS (eqn 4.5a)\n",
    "    Substep 2: Strong Stability Preserving 3rd-order Runge Kutta + 4th-order central finite difference scheme (eqn 4.5b)\n",
    "    Substep 3: Fourier-PS (eqn 4.5c)\n",
    "\n",
    "Energy is calculated from the free energy functional (eqn 1.1a), and the average concentration from (eqn 1.1b)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(options, wisdom=None):\n",
    "    L = options['L']\n",
    "    Nx = options['Nx']\n",
    "    eps = options['eps']\n",
    "    sigma = options['sigma']\n",
    "    m = options['m']\n",
    "    dt = options['dt']\n",
    "    it = options['it']\n",
    "    chiN = options['chiN']\n",
    "    prefix= options['prefix']\n",
    "    \n",
    "    io = writer(prefix,chiN,eps,sigma)\n",
    "    io.create_output_file(options)\n",
    "    \n",
    "    h = L/Nx\n",
    "    dx = np.linspace(0,L,Nx)\n",
    "    \n",
    "    # lv defines the levels for which sampling of the systems are made, e.g. sampling of energy and average concentration.\n",
    "    lv = int(it/100)\n",
    "\n",
    "    # define seed and uniform distribution.\n",
    "    # seed = 693\n",
    "    # np.random.seed(seed)\n",
    "    r = np.random.uniform(-1.,1.,(Nx,Nx))\n",
    "\n",
    "    # zeroing the mean of the uniform distribution.\n",
    "    if r.mean()<0:\n",
    "        r -= r.mean()\n",
    "    else:\n",
    "        r += r.mean()\n",
    "\n",
    "    # a scaling factor\n",
    "    C = 0.01\n",
    "\n",
    "    # a constant\n",
    "    tpiL = 2.*np.pi/L\n",
    "\n",
    "    # initialise arrays for pyFFTW.\n",
    "    u = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    v = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    vp = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    vpp = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "\n",
    "    # define initial conditions.\n",
    "    u[:,:] = m + C*r\n",
    "\n",
    "    # get wavenumber/frequency grids (k2 == k^2; k4 == k^4)\n",
    "    freqs = np.fft.fftfreq(Nx,(1.0/float(Nx)))\n",
    "    kx, ky = np.meshgrid(freqs,freqs)\n",
    "    k = kx + ky\n",
    "    k2 = kx**2+ky**2\n",
    "    k4 = k2**2\n",
    "    k = np.sqrt(k2)\n",
    "\n",
    "    # initialise arrays to store avg. concentration and energy values.\n",
    "    avgConc = np.zeros(int(it/lv+1))\n",
    "    energy = np.zeros(int(it/lv+1))\n",
    "\n",
    "    # if there is wisdom for FFT plans, import it.\n",
    "    if wisdom != None:\n",
    "        pyfftw.import_wisdom(wisdom)\n",
    "\n",
    "    # define pyFFTW plans for the forward and inverse FFTs.\n",
    "    fft_object = pyfftw.FFTW(u, vpp, threads=4, axes=(0,1))\n",
    "    ifft_object = pyfftw.FFTW(vpp, v, direction='FFTW_BACKWARD', threads=4, axes=(0,1))\n",
    "\n",
    "    fft_object = pyfftw.FFTW(v, vp, threads=4, axes=(0,1))\n",
    "    ifft_object = pyfftw.FFTW(vp, u, direction='FFTW_BACKWARD', threads=4, axes=(0,1))\n",
    "\n",
    "    # initialise some counters.\n",
    "    j = 0\n",
    "    t1 = 0\n",
    "    t2 = 0\n",
    "    t3 = 0\n",
    "\n",
    "    for i in range(it+1):\n",
    "        # the first Strang-splitting substep (eqn 4.5a): Fourier-PS.\n",
    "        tic1 = time.time()\n",
    "        vpp = np.exp(0.5*(-1.*eps**2*k4*tpiL**4+k2*tpiL**2-sigma)*dt)*fft_object(u)\n",
    "        v = ifft_object(vpp).real\n",
    "        toc1 = time.time()\n",
    "        t1 += (toc1-tic1)\n",
    "\n",
    "        # the second Strang-splitting substep (eqn 4.5b): SSP-RK3.\n",
    "        tic2 = time.time()\n",
    "        Fu = centralDiffShift(v,h,Nx)\n",
    "        v1 = ne.evaluate('v + dt * (Fu + sigma*m)')\n",
    "        Fu = centralDiffShift(v1,h,Nx)\n",
    "        v2 = ne.evaluate('3./4. * v + 1./4. * v1 + 1./4. * dt * (Fu + sigma*m)')\n",
    "        Fu = centralDiffShift(v2,h,Nx)\n",
    "        v = ne.evaluate('1./3. * v + 2./3. * v2 + 2./3. * dt * (Fu + sigma*m)')\n",
    "        toc2 = time.time()\n",
    "        t2 += (toc2-tic2)\n",
    "\n",
    "        # the third Strang-splitting substep (eqn 4.5c): Fourier-PS.\n",
    "        tic3 = time.time()\n",
    "        vp = np.exp(0.5*(-1.*eps**2*k4*tpiL**4+k2*tpiL**2-sigma)*dt)*fft_object(v)\n",
    "        u = ifft_object(vp).real\n",
    "        toc3 = time.time()\n",
    "        t3 += (toc3-tic3)\n",
    "\n",
    "        # plot intermediate morphologies.\n",
    "    #     if i%int((it-1)/10)==0:\n",
    "    #         fig1 = plt.figure(figsize=(3,3))\n",
    "    #         plt.imshow(u,origin='left', extent=[0,L,0,L])\n",
    "\n",
    "        # sample energy and average concentration at levels.\n",
    "        if i%lv==0:\n",
    "            dim = 2\n",
    "            u0 = u\n",
    "\n",
    "            # get average concentration.\n",
    "            for d in range(dim):\n",
    "                u0 = integrate.trapz(u0,dx)\n",
    "                \n",
    "            avg_conc = u0 / L**2\n",
    "            avgConc[j] = avg_conc\n",
    "\n",
    "            # get energy.\n",
    "            T1 = -1.j * k * tpiL * np.fft.fft2(u)\n",
    "            T1 = np.real(np.fft.ifft2(T1))\n",
    "            T1 = 0.5*eps**2*abs(T1)**2\n",
    "\n",
    "            T2 = 1./4*(np.real(u)**2-1.)**2\n",
    "\n",
    "            up = np.real(u)-m\n",
    "            for d in range(dim):\n",
    "                up = integrate.trapz(up,dx)\n",
    "            T3 = -1.*up\n",
    "            T3 = 0.5*sigma*abs(T3)**2\n",
    "\n",
    "            en = T1 + T2 + T3\n",
    "            en = np.real(en)\n",
    "            for d in range(dim):\n",
    "                en = integrate.trapz(en,dx)\n",
    "            energy[j] = en\n",
    "            j += 1\n",
    "            \n",
    "#             print(\"writing output at i = %i\" %i)\n",
    "            data = [avg_conc, en, u]\n",
    "            io.write_data(i, data)\n",
    "        \n",
    "        if it % 1000 == 0:\n",
    "            print(\"iteration = %i\" %it)\n",
    "\n",
    "    return pyfftw.export_wisdom()\n",
    "\n",
    "# normalise against initial energy. Possible since energy is unitless.\n",
    "# energy /= energy[0]\n",
    "\n",
    "# save array outputs\n",
    "# np.savetxt('uArray-dt=%s,it=%s,Nx=%s,sigma=%s,m=%s,seed=%s.txt' %(dt,it,Nx,sigma,m,seed), u)\n",
    "# np.savetxt('energy-dt=%s,it=%s,Nx=%s,sigma=%s,m=%s,seed=%s.txt' %(dt,it,Nx,sigma,m,seed), energy)\n",
    "# np.savetxt('avgConc-dt=%s,it=%s,Nx=%s,sigma=%s,m=%s,seed=%s.txt' %(dt,it,Nx,sigma,m,seed), avgConc)\n",
    "\n",
    "# print time taken.\n",
    "# print(\"time taken1: %.2f\" %(t1))\n",
    "# print(\"time taken2: %.2f\" %(t2))\n",
    "# print(\"time taken3: %.2f\" %(t3))\n",
    "# print(\"total time taken: %.2f\" %(t1+t2+t3))\n",
    "\n",
    "# # plot final figure.\n",
    "# fig = plt.figure(figsize=(3,3))\n",
    "# plt.imshow(u,origin='left', extent=[0,L,0,L])\n",
    "# plt.colorbar()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.08\n",
      "5.0\n",
      "20.0\n"
     ]
    }
   ],
   "source": [
    "# L = length of domain; Nx = spatial resolution;\n",
    "# h = spatial stepsize; dx = 1D spatial grid;\n",
    "# dt = temporal stepsize\n",
    "L = 2.*np.pi\n",
    "L = 1.0\n",
    "Nx = 128\n",
    "dt = 0.0005\n",
    "\n",
    "# define small parameters.\n",
    "eps = 0.08\n",
    "sigma = 5.0\n",
    "\n",
    "# L = 6.0\n",
    "chiN = 20.0\n",
    "N = 20\n",
    "l = 1./N\n",
    "f = 1.-0.3\n",
    "chi = chiN / N\n",
    "omega = (L**2)**(2/3)\n",
    "# omega = 1.0\n",
    "\n",
    "# eps2 = l**2 / (3. * f * (1. - f) * chi * omega)\n",
    "# eps = np.sqrt(eps2)\n",
    "\n",
    "# sigma = 36. * omega / (1. * N**2 * l**2 * f**2 * (1.-f)**2 * chi)\n",
    "\n",
    "# print(eps2)\n",
    "print(eps)\n",
    "print(sigma)\n",
    "print(chiN)\n",
    "\n",
    "# u-bar == m.\n",
    "m = 0.0\n",
    "\n",
    "# number of iterations\n",
    "it = 10000\n",
    "it += 1\n",
    "\n",
    "options = {\n",
    "    'L' : L,\n",
    "    'Nx' : Nx,\n",
    "    'eps' : eps,\n",
    "    'sigma' : sigma,\n",
    "    'm' : m,\n",
    "    'dt' : dt,\n",
    "    'it' : it,\n",
    "    'prefix' : 'dft',\n",
    "    'N' : N,\n",
    "    'l' : l,\n",
    "    'f' : f,\n",
    "    'chi' : chi,\n",
    "    'chiN' : chiN,\n",
    "    'omega' : omega\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./output/dft_chiN=20.00_eps=0.0800_sigma=5.0000.h5\n"
     ]
    }
   ],
   "source": [
    "iwisdom = dft(options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array(0.24999997, dtype=float32), array(0.2499999, dtype=float32), array(0.24999776, dtype=float32), array(0.24993956, dtype=float32), array(0.24829315, dtype=float32), array(0.21716946, dtype=float32), array(0.1530988, dtype=float32), array(0.13617572, dtype=float32), array(0.1334275, dtype=float32), array(0.13309368, dtype=float32), array(0.1330296, dtype=float32), array(0.13299419, dtype=float32), array(0.13296665, dtype=float32), array(0.13294499, dtype=float32), array(0.13292821, dtype=float32), array(0.13291538, dtype=float32), array(0.13290565, dtype=float32), array(0.13289836, dtype=float32), array(0.13289295, dtype=float32), array(0.132889, dtype=float32), array(0.13288616, dtype=float32), array(0.13288414, dtype=float32), array(0.13288277, dtype=float32), array(0.13288188, dtype=float32), array(0.13288134, dtype=float32), array(0.13288108, dtype=float32), array(0.13288099, dtype=float32), array(0.13288103, dtype=float32), array(0.13288116, dtype=float32), array(0.13288136, dtype=float32), array(0.13288158, dtype=float32), array(0.13288184, dtype=float32), array(0.1328821, dtype=float32), array(0.13288236, dtype=float32), array(0.13288261, dtype=float32), array(0.13288285, dtype=float32), array(0.13288309, dtype=float32), array(0.1328833, dtype=float32), array(0.13288349, dtype=float32), array(0.13288367, dtype=float32), array(0.13288383, dtype=float32), array(0.13288398, dtype=float32), array(0.13288411, dtype=float32), array(0.13288425, dtype=float32), array(0.13288435, dtype=float32), array(0.13288446, dtype=float32), array(0.13288455, dtype=float32), array(0.13288462, dtype=float32), array(0.1328847, dtype=float32), array(0.13288477, dtype=float32), array(0.13288482, dtype=float32), array(0.13288487, dtype=float32), array(0.13288492, dtype=float32), array(0.13288496, dtype=float32), array(0.132885, dtype=float32), array(0.13288502, dtype=float32), array(0.13288505, dtype=float32), array(0.13288508, dtype=float32), array(0.13288511, dtype=float32), array(0.13288513, dtype=float32), array(0.13288514, dtype=float32), array(0.13288516, dtype=float32), array(0.13288517, dtype=float32), array(0.13288519, dtype=float32), array(0.1328852, dtype=float32), array(0.1328852, dtype=float32), array(0.13288522, dtype=float32), array(0.13288522, dtype=float32), array(0.13288523, dtype=float32), array(0.13288523, dtype=float32), array(0.13288525, dtype=float32), array(0.13288525, dtype=float32), array(0.13288525, dtype=float32), array(0.13288526, dtype=float32), array(0.13288526, dtype=float32), array(0.13288526, dtype=float32), array(0.13288526, dtype=float32), array(0.13288526, dtype=float32), array(0.13288526, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32), array(0.13288528, dtype=float32)]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAURklEQVR4nO3dbYxc133f8e9vl4+RRFKyGMflg8g4amAlaqVmTSUQ7AKJLNNtIOWFgsiFWwVVK8SxgBRCUSgwaqfKGycGGr+oipioVbRpU0WxU5d1laiCLTdwE0lcyrJiSlVM0460pQoxoqwH64FP/76YWXI43iUvyVkO9+z3Awzm3nPPvfO/Evibu2fO3ElVIUlq18S4C5AkLSyDXpIaZ9BLUuMMeklqnEEvSY1bNu4Chl1++eW1ZcuWcZchSYvK7t27/7qq1s+17YIL+i1btjA9PT3uMiRpUUnyV/Ntc+hGkhpn0EtS4wx6SWqcQS9JjTPoJalxnYI+yfYkzybZm+TuObbfleTpJE8l+XKSKwa2HU3yZP+xc5TFS5JO77TTK5NMAvcCHwBmgF1JdlbV0wPdvg5MVdUbST4K/DbwS/1tb1bVNSOuW5LUUZd59NuAvVW1DyDJ/cDNwPGgr6pHBvo/CnxklEV2cuj78LXPzL99YrL/WA6bfwY2bYPk/NUnSWPSJeg3AM8PrM8A152i/+3AHw+sr0oyDRwBPlVVXxzeIckdwB0Amzdv7lDSHA6/CX/66Xk2znHP/Xf9bbjuV+Anb4FlK87uNSVpEegS9HNd9s75ayVJPgJMAX93oHlzVe1P8qPAV5L8RVV9+6SDVe0AdgBMTU2d3S+hXHQ5/Mb35t9+7BjUUTj0Onzzj+DxHfDFj8LB78DPfvysXlKSFoMuH8bOAJsG1jcC+4c7JbkB+DhwU1W9PdteVfv7z/uArwLXnkO9Z29iAiaXw+pL4b23w68+Cu+8Gl54cizlSNL50iXodwFXJtmaZAVwK3DS7Jkk1wKfpRfyLw60X5pkZX/5cuB6Bsb2xyqBy7bA954bdyWStKBOG/RVdQS4E3gIeAZ4oKr2JLknyU39bp8GLgb+cGga5XuA6STfAB6hN0Z/YQQ9wNrNvaD3d3MlNazT3Sur6kHgwaG2Twws3zDPfn8GXH0uBS6odZvh8Bvwxku9MX5JatDS/mbsuv4MH4dvJDVsiQd9/zNmg15Sw5Z20K816CW1b2kH/ep1sGqtQS+paUs76KE3Tv/K86fvJ0mLlEE/O8VSkhpl0K9zLr2kthn06zb37n/z5svjrkSSFoRB71x6SY0z6J1LL6lxBr1X9JIaZ9CvWgcr1zjFUlKzDPrkxMwbSWqQQQ+9WyEY9JIaZdCDV/SSmmbQQy/o334V3jzFb85K0iJl0INTLCU1zaAHp1hKappBD7Duit6zQS+pQQY9wOpLYcXFzqWX1CSDHnpz6S/5EXjthXFXIkkjZ9DPWrUW3np13FVI0sgZ9LNWrulNsZSkxhj0s1athbdeGXcVkjRyBv2sVWscupHUJIN+1qq1Dt1IapJBP2vlWjj8Bhw9PO5KJGmkDPpZq9b0nh2+kdQYg37WqrW957e8sZmkthj0s1b2r+gdp5fUGIN+1vEreqdYSmpLp6BPsj3Js0n2Jrl7ju13JXk6yVNJvpzkioFttyX5Vv9x2yiLHynH6CU16rRBn2QSuBf4EHAV8OEkVw11+zowVVV/C/g88Nv9fS8DPglcB2wDPpnk0tGVP0IO3UhqVJcr+m3A3qraV1WHgPuBmwc7VNUjVfVGf/VRYGN/+YPAw1V1sKpeBh4Gto+m9BFz6EZSo7oE/QZg8P69M/22+dwO/PGZ7JvkjiTTSaYPHDjQoaQFsPKS3rNDN5Ia0yXoM0dbzdkx+QgwBXz6TPatqh1VNVVVU+vXr+9Q0gKYmOwN33hFL6kxXYJ+Btg0sL4R2D/cKckNwMeBm6rq7TPZ94LhHSwlNahL0O8CrkyyNckK4FZg52CHJNcCn6UX8i8ObHoIuDHJpf0PYW/st12YVnlFL6k9y07XoaqOJLmTXkBPAvdV1Z4k9wDTVbWT3lDNxcAfJgF4rqpuqqqDSX6T3psFwD1VdXBBzmQUvFWxpAadNugBqupB4MGhtk8MLN9win3vA+472wLPq5Vr4PX/N+4qJGmk/GbsIK/oJTXIoB/kj49IapBBP2h2emXNOXtUkhYlg37QqrVQR3s/QCJJjTDoB3ljM0kNMugHeb8bSQ0y6Aet7Ae9346V1BCDftDxoRuv6CW1w6Af5NCNpAYZ9INWekUvqT0G/aBVjtFLao9BP2j5aphY5vRKSU0x6Acl/viIpOYY9MNWrXXoRlJTDPph/viIpMYY9MNWrXWMXlJTDPphjtFLaoxBP8wxekmNMeiHOXQjqTEG/bCVa+DQa3Ds6LgrkaSRMOiH+e1YSY0x6Id5B0tJjTHoh630V6YktcWgH+bQjaTGGPTDHLqR1BiDftjxHx/xil5SGwz6YSv9lSlJbTHoh80O3ThGL6kRBv2wyeWw/Ie8opfUDIN+Lt7YTFJDOgV9ku1Jnk2yN8ndc2x/f5InkhxJcsvQtqNJnuw/do6q8AW1eh28cXDcVUjSSCw7XYckk8C9wAeAGWBXkp1V9fRAt+eAXwb++RyHeLOqrhlBrefPD18Fzz8+7iokaSS6XNFvA/ZW1b6qOgTcD9w82KGqvltVTwHHFqDG82/TdfDqDLwyM+5KJOmcdQn6DcDzA+sz/bauViWZTvJokl+Yq0OSO/p9pg8cOHAGh14gm7b1np9/bLx1SNIIdAn6zNFWZ/Aam6tqCvgHwGeSvPsHDla1o6qmqmpq/fr1Z3DoBfIjV/dm3jh8I6kBXYJ+Btg0sL4R2N/1Bapqf/95H/BV4NozqG88JpfDhp/yil5SE7oE/S7gyiRbk6wAbgU6zZ5JcmmSlf3ly4HrgadPvdcFYtM2eOEpOPT9cVciSefktEFfVUeAO4GHgGeAB6pqT5J7ktwEkOS9SWaAXwQ+m2RPf/f3ANNJvgE8AnxqaLbOhWvTT0Mdhf/7xLgrkaRzctrplQBV9SDw4FDbJwaWd9Eb0hne78+Aq8+xxvHYONV7fv4x2Pq+8dYiSefAb8bO54cug8t/3HF6SYueQX8qm6/rzbw51sbXAyQtTQb9qWy6Dt76Hrz0rXFXIklnzaA/lU0/3Xt+7tHx1iFJ58CgP5V3vBtWXwbPPghHj4y7Gkk6Kwb9qSQw9Y/hL/8E/v12OPidcVckSWes0/TKJe3n/iX88HvgS3fB774Ptv0TWLMBLlrf+zWqiWUwsRwmJiGTvTeHTPQf/WVyYp0MPfOD64N3ncjgHSg6tM8nHfqcizqTu2JImtPkCrh49LeBMei7uPqW3jdl/9vH4Gu/M+5qJLVqwxT80y+P/LAGfVfrNsNt/703Vv/GS/D9A/D2a3DsSP9xFOpY/3G0d4U7u06dWIf+1W/94PPxbbMGlju1z+d8XW0v8F8NUusuunxBDmvQn6nJZXDJO3sPSVoE/DBWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN6xT0SbYneTbJ3iR3z7H9/UmeSHIkyS1D225L8q3+47ZRFS5J6ua0QZ9kErgX+BBwFfDhJFcNdXsO+GXg94f2vQz4JHAdsA34ZJJLz71sSVJXXa7otwF7q2pfVR0C7gduHuxQVd+tqqeAY0P7fhB4uKoOVtXLwMPA9hHULUnqqEvQbwCeH1if6bd1cS77SpJGoEvQZ4626nj8TvsmuSPJdJLpAwcOdDy0JKmLLkE/A2waWN8I7O94/E77VtWOqpqqqqn169d3PLQkqYsuQb8LuDLJ1iQrgFuBnR2P/xBwY5JL+x/C3thvkySdJ6cN+qo6AtxJL6CfAR6oqj1J7klyE0CS9yaZAX4R+GySPf19DwK/Se/NYhdwT79NknSepKrrcPv5MTU1VdPT0+MuQ5IWlSS7q2pqrm1+M1aSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMZ1Cvok25M8m2Rvkrvn2L4yyR/0tz+WZEu/fUuSN5M82X/87mjLlySdzrLTdUgyCdwLfACYAXYl2VlVTw90ux14uap+LMmtwG8Bv9Tf9u2qumbEdUuSOupyRb8N2FtV+6rqEHA/cPNQn5uB/9Bf/jzwc0kyujIlSWerS9BvAJ4fWJ/pt83Zp6qOAK8A7+hv25rk60n+V5L3zfUCSe5IMp1k+sCBA2d0ApKkU+sS9HNdmVfHPi8Am6vqWuAu4PeTrPmBjlU7qmqqqqbWr1/foSRJUlddgn4G2DSwvhHYP1+fJMuAtcDBqnq7ql4CqKrdwLeBv3muRUuSuusS9LuAK5NsTbICuBXYOdRnJ3Bbf/kW4CtVVUnW9z/MJcmPAlcC+0ZTuiSpi9POuqmqI0nuBB4CJoH7qmpPknuA6araCXwO+L0ke4GD9N4MAN4P3JPkCHAU+JWqOrgQJyJJmluqhofbx2tqaqqmp6fHXYYkLSpJdlfV1Fzb/GasJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJalynoE+yPcmzSfYmuXuO7SuT/EF/+2NJtgxs+/V++7NJPji60iVJXZw26JNMAvcCHwKuAj6c5KqhbrcDL1fVjwG/A/xWf9+rgFuBnwC2A/+2fzxJ0nmyrEOfbcDeqtoHkOR+4Gbg6YE+NwO/0V/+PPBvkqTffn9VvQ18J8ne/vH+fDTln3DoyDEe3ffSSW3JwDIhgQATE2Gy/1i1bJKLVk6yesUka1YtZ9Xy0b0PVRVVcKyKY/3nE+u9thp+ptdntl8NHAdOHON4+/HXAhjc1m/nxL7VPz4n7TPH8nx9hs5t3vOe97/Hme7RZd/TO4ddLxjncv5aPC5aOclP/I21Iz9ul6DfADw/sD4DXDdfn6o6kuQV4B399keH9t0w/AJJ7gDuANi8eXPX2k/y2luH+Uf3PX5W+w5asWyCNauWs3rFBJPpvRkAx4P66LETj2NVHJldPlYcnQ30YyeCXJK6umbTOr74setHftwuQZ852oYjbL4+XfalqnYAOwCmpqbOKh7XrF7OFz76MwPHPPkFq05cAc8G9pFjxduHj/LGoaN8/9BRXnvrMK+8eZhX3zzMW4ePHQ90AhMJk7PP/b8GJibCsomc1JZw/A0iCRP99Yn+G8bkRK8t9NomBv7KCEB/n4n01tPvy3BbTvyVMitD+8xum22bXTv5L52T9z1xrHmWmWPnOcy3KZl7yykONWcdZyqdXuHCdi7nr8Xh4pVdIvnMdTnqDLBpYH0jsH+ePjNJlgFrgYMd9x2J5ZMT/NQVly3EoSVpUesy62YXcGWSrUlW0PtwdedQn53Abf3lW4CvVG8Qdydwa39WzlbgSuDcx1ckSZ2d9oq+P+Z+J/AQMAncV1V7ktwDTFfVTuBzwO/1P2w9SO/NgH6/B+h9cHsE+FhVHV2gc5EkzSGnmj0xDlNTUzU9PT3uMiRpUUmyu6qm5trmN2MlqXEGvSQ1zqCXpMYZ9JLUOINekhp3wc26SXIA+KtzOMTlwF+PqJzFYqmd81I7X/Ccl4pzOecrqmr9XBsuuKA/V0mm55ti1Kqlds5L7XzBc14qFuqcHbqRpMYZ9JLUuBaDfse4CxiDpXbOS+18wXNeKhbknJsbo5cknazFK3pJ0gCDXpIa10zQJ9me5Nkke5PcPe56FlqS+5K8mOSb467lfEmyKckjSZ5JsifJr427poWWZFWSx5N8o3/O/2rcNZ0PSSaTfD3Jl8Zdy/mS5LtJ/iLJk0lGegvfJsbok0wCfwl8gN6vWu0CPlxVT59yx0UsyfuB14H/WFU/Oe56zock7wLeVVVPJLkE2A38QuP/nwNcVFWvJ1kOfA34tap69DS7LmpJ7gKmgDVV9fPjrud8SPJdYKqqRv4lsVau6LcBe6tqX1UdAu4Hbh5zTQuqqv6U3o+8LBlV9UJVPdFffg14hjl+bL4l1fN6f3V5/7H4r85OIclG4O8D/27ctbSilaDfADw/sD5D4wGw1CXZAlwLPDbeShZefxjjSeBF4OGqav2cPwP8C+DYuAs5zwr4n0l2J7ljlAduJegzR1vTVz1LWZKLgS8A/6yqXh13PQutqo5W1TXARmBbkmaH6pL8PPBiVe0edy1jcH1V/R3gQ8DH+sOzI9FK0M8AmwbWNwL7x1SLFlB/nPoLwH+uqj8adz3nU1V9D/gqsH3MpSyk64Gb+uPV9wM/m+Q/jbek86Oq9vefXwT+K70h6ZFoJeh3AVcm2ZpkBb0fJ9855po0Yv0PJj8HPFNV/3rc9ZwPSdYnWddfXg3cAPyf8Va1cKrq16tqY1Vtoffv+CtV9ZExl7XgklzUn2BAkouAG4GRzahrIuir6ghwJ/AQvQ/oHqiqPeOtamEl+S/AnwM/nmQmye3jruk8uB74h/Su8p7sP/7euItaYO8CHknyFL0LmoeraslMOVxC3gl8Lck3gMeB/1FVfzKqgzcxvVKSNL8mruglSfMz6CWpcQa9JDXOoJekxhn0ktQ4g17qIMm6JL867jqks2HQS92sAwx6LUoGvdTNp4B397+k9elxFyOdCb8wJXXQv1vml5bKvf/VFq/oJalxBr0kNc6gl7p5Dbhk3EVIZ8OglzqoqpeA/53km34Yq8XGD2MlqXFe0UtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Lj/D2pc48KvbYZaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD7CAYAAABqkiE2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df6wmV1nHP89799670FJKRWDpNm5JKloNSkOwiMGGqkAlNCZAikQXrGk0gCga28ofaCIJqEEwGnBDgdaULuWHtkEUsEKMiawsglioleWH7dKFFvnZ0u798T7+cc7MO3c7Z+/73vd5zz0z83ySzZkzM3dmdt73nPeZM+f7fURVcRxnuIx2+wIcx9ldvBNwnIHjnYDjDBzvBBxn4Hgn4DgDxzsBxxk4C+sEROS5InKniBwTkWsWdR7HceZDFjFPQESWgP8Bfh44DnwSeImqft78ZI7jzMWeBR336cAxVf0SgIgcBi4HWjuBFdmre+UMQMGiT5KtC5OqPGzblnUik23V8inbdDTZXxv76OiU/YQ6zlLZWjJqLDe31cfQrece6eRy4rbRSBnF5SUZh3Wi9fJSvJGjqi7KiMl+YZ/xw9aN0Mlp4zEmt0A45bIRJutorJsGbXzYp37s1TYFqh+qap8xsmW5KsdaLYcbWdU3GW1ZBtjUxrrq73RUL2u9blLXcfVhVRcpxNtXf2TVNhk3licX+/D9FGSsD1sXTw7x/z75u8k62tbRsi3y3fH/fUNVf5BTWFQncC5wd6N+HPip5g4ichVwFcBeHsnFq8+DzU20viHx7k4bqdStZISMJssAshRb12iE7In/5aWlUO7ZM9m+vBzXLaHLcb9Y6nLYf7yyxHhlsgywuTpivBK/UI1yc7VaJu4Xv1irzXUajwXjank1/t9XYgNd3WR5ZQOA1Vg+YmWdM1bWADhjOZRn7jnJo5YfAuCsPaE8c+kkAI9aeogzl8K6M0Zx3ejBenmvrIdtss6qbMZ14fzL8XauiLAcG91yvLfLssQoNqxR3LYk0z1lbsbPeIwyjq1pM37e64RrGKvyUNxvPf7dQyqsx173Ad0T1y3z/fFqWFeX4SZ/b/wI7t/cG5Zj+d2NvTywEfb7XizvX1/l/rWw/OB6+C48uBbKtZN72FgLn7eeDKWsjRidDP/n0Vr8v1f1k7AUPhbiR8BoTevlpTWty6WT1XL8vOtyk9FauA+yHsu1DVgP3wHZCOtY30A3wjqqcjNs081xvfyRB274X1pYVCfQ9lOwpTWr6iHgEMBZo3O07gDmaPyhkNbGD4QOoNH4630ajR8IHUBL46/KZuMP62RL4w/bpLXxh322Nv6wTVsbP8DyysaWxg9wxsralsYP8Kjlh1obf6g/tKXxQ+gMmo0fYFU2Wxs/wDKypfEDjBjN1fhDOW5t/AAP6XhL4wdY19GWxg/w/fFqa+MHuH9z75bGD/DAxuqWxg9w/9pqa+MH2Fhb2tL4AUYnpbXxQ+gAmo0/7LO18Yd12tr4q3JL4wdY39jS+IHQAbQ1/ljXuC7FojqB48B5jfp+4J7k3sqkA5hljEJka+MPC+2NH0IH0Gz8EDqAZuMHWN7T2virstn4If7qNxp/WNfe+MM+uqXxh3Lc2vgh/Po3Gz+EX/9m44fw69/W+CE0+Gbjh/Dr32z8Yd24tfFD+PVvNv5QytSNH0IH0Gz8YZ22Nn4Iv/7Nxg/h17/Z+CH8+rc1fgi//s3GD+HXv9n4Ifz6tzV+CL/+zcYP4de/rfFDaPDNxh+26ZbGH9aNWxs/hF//ZuOH+OvfbPwQOoCWxh/qm/VyikW9HfgkcIGInC8iK8AVwK0LOpfjOHOwkEhAVTdE5JXAh4El4B2q+rnT/9EMUcB2z/9tEQDM9PzfFgHAbM//bRFAOFb7839bBACzPf+3RQAw2/N/WwQQysU9/7dFAKE+/fN/WwQAsz3/t0UAkH7+b4sAwj7TP/+3RQDhnNM//7dGANW28enb1aIeB1DVDwEfWtTxHcexYWGdwMxMEwVMOQjYFgHU+0w5CNgWAYT69IOAbRFA2NY+CNgWAcBsg4BtEUBVn3YQsC0CgMUOArZFADDbIGBbBACzDQK2RQCQHgRsiwCq+rSDgK0RAMw2CNgSAQBbB9sT+LRhxxk45UQC2zHLm4C2CABmehPQFgHAbG8C2iKAULa/CWiLAGC2NwFtEQDM9iagLQII5eLeBLRFADDbm4C2CABmexPQFgFA+k1AWwQQtk3/JqAtAoDZ3gS0RgBhYdsou/xOYCeDgG2NH2YaBGxr/FU57SBga+OH5CBgW+OH2QYB2xo/zDYI2Nb4YbGDgG2NP+wz/SBgW+OH2QYB2xo/pAcB2xp/VU47CNja+GN96kHAtsYPUz1m++OA4wycciOBOQYB2yIAmG0QsC0CCNumHwS0mAYMsw0CtkUAMNsgoEUEEMrpBwHbIgCYbRCwLQKA2QYBLaYBh31mGARsiQBgtkHAHU+3xyMBxxk85UUCbRFArM+lBYCZBgHbIoBQz6sFgNkGAdsiAJhtENAiAgjrpx8EbIsAYLZBwLYIAGYbBLTQAoT69IOArREAzDYIuIMIoKKsTiDxBgBmHARsafyhPv0gYFvjr+o5tQAw2yBgW+OH2QYBLd4AwGyDgG2NP5TTDwK2NX6YbRDQQgsQjjXDIGBb44fZBgHn8AXxxwHHGTjlRAIxCuiTF4CFFgDcCyCnF4CFFiCUtl4AVoOAbXgk4DgDp6BIYNQ7QxALLUCouyEI5DEEsdACVKWpIcg8g4DVYHti9zI6AZnMA+iTIYjFNGBwQ5CchiAW04DB3hBkx4OAjcH2FP444DgDp4xIgEYU0CNDEAstALghSE5DEAstQDinsSGIhfFuwmXMIwHHGThFRALCZDJQnwxBLLQAVd0NQfIYgphoAcDeEMTAfTuFRwKOM3CKiAQQmUQBPTIEsdACQH5DEAstAHTTEMRCCwD5DEFqptHcJCikE2jMA+iRIYiFFgDyG4JYaAHCPt0zBDHRAsR6DkOQcNLpNDcp/HHAcQZOGZEAUkcBfTIEsdACQH5DEAstAHTTEMRCCwCZDEFmzb+ZwCMBxxk4ZUQCIpMooEeGIBZaAMhvCGKhBYBuGoKYaAFgsYYgO03Cm2DHnYCInAfcADwBGAOHVPUtInIO8B7gAPAV4MWq+q3TH6wxD6BHhiAWWgDIbwhioQUIZfcMQUy0ALA4Q5A5kvCmmOdxYAP4XVX9UeBi4BUiciFwDXCbql4A3BbrjuMUyo4jAVU9AZyIy98TkTuAc4HLgUvibtcDHweuPu3BROoooE+GIBZaAMg/CGihBYBuGoJYaAFgFwYBp9HcJDAZGBSRA8BTgSPA42MHUXUUj0v8zVUiclREjq5tPmhxGY7j7IC5BwZF5Ezg/cBvq+p3peqxtkFVDwGHAB69+gSto4AeGYJYaAEg/yCghRYAumkIYqIFiPWcg4BTaW4SzNUJiMgyoQO4UVU/EFd/XUT2qeoJEdkH3DvFgeoOoE+GIBbTgCH/IKDFNGDopiGIxTTgUM88CDjFdPsUO34ckPCTfx1wh6q+qbHpVuBgXD4I3LLTcziOs3jmiQSeCfwK8F8i8pm47g+ANwA3i8iVwF3Ai7Y9klBHAX0yBLHQAoQy7yCghRYAumkIYqIFqLblHAScRnOTYJ63A/8KpAYALt3pcR3HyUsxMwarKKBPhiAWWgDIPwhooQWAjhqCGGgBIP8g4DSamxSuHXCcgVNEJKAjqaOAPhmCWGgBQpnXEMRCCwDdNAQx0QKEBVtDEIMkvCmK6ASQyTyAPhmCWGgBIL8hiIUWADpqCGKhBQBzQxCLJLwp/HHAcQZOEZGAitRRQJ8MQUrRAsBsg4AWWgDopiFIMVoAmGkQcBrNTQqPBBxn4BQRCSCTyUB9MgQpRQsA+ZODQkcNQQrRAsCMg4BTaG5SFNEJ6GgyD6BPhiClaAEgf3JQ6KghSClaAJhpEHAazU0KfxxwnIFTRCSATKKAPhmClKIFCPW8yUGhm4YgxWgBwDwJbwqPBBxn4BQRCahMJgP1yRCkFC0A5E8OCh01BClEC1DvY5iEN0URnQCjyTyAPhmCWEwDhvyGIBbTgKGbhiAm04DB3BDEIv9mCn8ccJyBU0QkEB4H+mcIYqEFgPyGIBZaAOioIYiFFgDsDUEMkvCm8EjAcQZOOZFAjAL6ZAhiEQGEMq8hiIUWALppCGIRAYTC1hDEIglvCo8EHGfgFBEJMJpMBuqTIYhFBBDW5zUEsdACQEcNQSwigFi3NASxSMKboohOYDIw2C9DEIvXgJDfEMRCCwAdNQQxeA0I9oYgFkl4U/jjgOMMnCIigSAlVjcEKcQQxEILAG4IYmkIYpGEN4VHAo4zcIqIBFQmk4HcEGT3DUEstADghiCWhiAWSXhTWGQlXgKOAl9V1eeLyPnAYeAc4D+AX1HVtdMeZNToANwQZNcNQSy0AOCGIJaGIBZJeFNYPA68GrijUX8j8OeqegHwLeBKg3M4jrMg5k1Nvh/4ReD1wGtipuJnA78cd7ke+EPgrac/0CQKcEOQ3TcEsdAChH3cEAQwMQSxSMKbYt5I4M3A70P8RsMPAN9W1fjJcBw4t+0PReQqETkqIkc3739gzstwHGen7DgSEJHnA/eq6qdE5JJqdcuurd2Qqh4CDgGsHtivVRTghiC7bwhioQWo6m4IgokhiEkS3gTzPA48E3iBiFwG7AXOIkQGZ4vInhgN7AfumeMcjuMsmB13Aqp6LXAtQIwEfk9VXyoi7wVeSHhDcBC4ZduDjbSOAvrkCmShBYD8rkAWWoCwrYOuQBZaADB3BbJIwptiEfMErgYOi8gfA58GrtvuD0Qm8wD6ZAhioQUI++Q1BLHQAlRl5wxBLLQAYG4IYpGEN4VJJ6CqHwc+Hpe/BDzd4riO4yyeImYMimgdBfTJEMRCCwD5DUEstABhn+4ZglhoAcDeEMQiCW8K1w44zsApIhIYjbSOAvpkCGKhBYD8hiAWWoBQ76AhiIUWAMwNQSyS8KYooxMQrTuAPhmCWGgBQpnXEMRCCxCO1T1DEAstQKjbGoJYJOFN4Y8DjjNwiogElmRcRwF9MgSx0AJAfkMQCy1AKDtoCGKgBQB7QxCLJLwpPBJwnIFTRCQwEq2jgD4ZglhoASC/IYiFFqAqu2YIYqEFqEpLQxCLJLwpiugElmRcdwB9MgSxmAYM+Q1BLKYBQ0cNQQymAVelpSGIRRLeFP444DgDp4xIAK2jgD4ZglhoASC/IYiFFiCcs4OGIAZaALA3BLFIwpvCIwHHGThFRAIjGddRQJ8MQSy0AJDfEMRECwCdNASx0AKEurEhiEES3hQeCTjOwCkiElgSraOAPhmCWGgBIL8hiIUWADpqCGKgBQB7QxCLJLwpiugERozrDqBPhiAWWgDIbwhiogWI9c4ZghhoAarS0hDEIglvCn8ccJyBU0YkIFpHAX0yBLHQAkB+QxALLQB00xDEQgsQttkaglgk4U3hkYDjDJwiIoElxnUU0CdDEAstAOQ3BDHRAkA3DUEMtAChbmsIYpGEN0URncCIcd0B9MkQxEILAPkNQUy0ANBJQxALLUBVtzQEsUjCm8IfBxxn4JQRCYjWUUCfDEEstACQ3xDEQgsA3TQEsdACgL0hiEUS3hQeCTjOwCkjEkDrKKBPhiAWWgDIbwhiogWI9d3WAtT7ZNQChH1sDUEskvCmmKsTEJGzgbcDPw4o8GvAncB7gAPAV4AXq+q3TnscJvMA+mQIYjENGPIbglhMAw717hmCWEwDDvvYGoJY5N9MMe/jwFuAf1TVHwF+ArgDuAa4TVUvAG6LdcdxCmXHkYCInAU8C3gZgKquAWsicjlwSdztekKOwqtPd6wRWkcBfTIEsdACQH5DEBMtQLWtY4YgFlqA8Pe2hiAWSXhTzBMJPAm4D3iniHxaRN4uImcAj1fVEwCxfFzbH4vIVSJyVESOfvOb4zkuw3GceZhnTGAPcBHwKlU9IiJvYYbQX1UPAYcAnvKUZa2igD4ZglhoAWAXDEEMtADQTUMQCy1A2GZrCGKRhDfFPJHAceC4qh6J9fcROoWvi8g+gFjeO8c5HMdZMDuOBFT1ayJyt4g8WVXvBC4FPh//HQTeEMtbtjuWiNRRQJ8MQSy0AJDfEMRECxAWOmcIYqEFCKWtIYhFEt4U884TeBVwo4isAF8CXk6ILm4WkSuBu4AXbXcQYTIPoE+GIBZaANgFQxALLQB00hDERAsA5oYgFkl4U8zVCajqZ4CntWy6dJ7jOo6TjyJmDAqTyUB9MgSx0AJAfkOQYrQAkN0QxEILAPaGIBZJeFO4dsBxBk4hkYDUUUCfDEEstACwC4YghWgBIL8hiIUWAOwNQSyS8KYopBOYzAPokyGIhRYAdsEQpBQtAGQ3BLHQAoC9IYhFEt4U/jjgOAOniEgAdj4PoGRDEAstQNgnryFIMVoAyG4IYqEFAHtDEIskvCk8EnCcgVNEJCBTRgFdMwSx0AJU9SEmB6336VhyULA3BLFIwpvCIwHHGThFRALT0EVXIAstQNiW1xXIRAsAnXQFstACgL0r0CKT8BbfCXTZEMRCC1CVnUsOCp00BLHQAoC9Icgik/D644DjDJxiI4E+GIJYaAHCPnkNQSwigFB0zxDEQgsA9oYgi8q/Ef7GcZxBU1wk0CdDEAstQKhnNgSxiABivWuGIBZaALA3BLHOv9GkqE6gb4YgFlqAcKzMhiAGbwCgo4YgBloAsDcEscy/cSr+OOA4A6eISEDROgrokyGIhRYglG4IUi3DYg1BLLQAYG8IYjUI2IZHAo4zcIqIBGAyGahPhiAWWoCqdEOQPIYgFlqAULc1BLHIv5GiiE5AmcwD6JMhiMU0YHBDkK4lBwV7QxCLJLwp/HHAcQZOEZEATOYB9MkQxEILEM7phiCQxxDEQgsA9oYgFu7bKTwScJyBU0QkoEyigD4ZgphoAcANQTqWHLSqWxqCWLhvp5grEhCR3xGRz4nI7SJyk4jsFZHzReSIiHxBRN4TU5Q5jlMoO44ERORc4LeAC1X1QRG5GbgCuAz4c1U9LCJvA64E3nq6YymTyUB9MgSx0ALALhiCWGgBoJOGIBZaAMhnCFIxjeYmxbyPA3uAR4jIOvBI4ATwbOCX4/brgT9ku06g0QH0yRDERAsQ651LDgqdNASx0AJAPkMQmF5zk2LHjwOq+lXgzwiZh08A3wE+BXxbVeM3lePAuW1/LyJXichRETn6f988/eil4ziLY57HgccAlwPnA98G3gs8r2XX1m5IVQ8BhwCe8pRlraKAPhmCWGgBIL8hiIUWALppCGKhBYA8hiCzJuFNMc/A4M8BX1bV+1R1HfgA8NPA2SJSdS77gXvmOIfjOAtmnjGBu4CLReSRwIPApcBR4GPAC4HDwEHglu0ONEbqKKBPhiAmWgDIbwhioQWAThqCWGgBYLGGIDtNwptix52Aqh4RkfcB/wFsAJ8mhPd/DxwWkT+O667b9lhM5gH0yRDERAsA2Q1BLLQAod49QxALLQAszhBkniS8KeZ6O6CqrwNed8rqLwFPn+e4juPko4gZg2OkjgL6ZAhioQWAXRgENNACQDcNQSy0AJB/EHAazU0K1w44zsApJhKoooA+GYKYaAFivWtagKrsmiGIhRYA8g8CTqO5SVFGJ6BSdwB9MgSxmAYc6pkHAQ2mAVdl1wxBLKYBQ/5BwGmm26fwxwHHGThlRAKM6iigT4YgJlqAalvHkoNCNw1BLLQAocw7CDiN5iaFRwKOM3DKiARU6iigV4YgBloAyD8IaKEFCPXuGYJYaAEg/yDgNJqbFB4JOM7AKSIS2GRURwF9MgQx0QKEhbyGIAZaAOimIYiFFiCUtoYgFkl4UxTRCYxV6g6gV4YgFloAyG8IYqAFqMquGYJYaAHA3hDEIglvCn8ccJyBU0QksMmojgL6ZAhSjBYAZhoEtNAChG3dMwQpRQsAsw0CTqO5SeGRgOMMnDIiAR3VUUCvDEEK0QLAjIOABlqAUO+eIUgpWgCwT8KboohOYKxSdwC9MgQpRQsAMw0CWmgBqnrXDEFK0QKAfRLeFP444DgDp4hIYFOljgL6ZAhSjBYAsicHhW4agpSiBQh12yS8KTwScJyBU0QkMNZRHQX0yhCkEC1AvU9GLUDYp3uGIKVoAcA+CW+KIjqBTZW6A+iTIYjJNGDIbghiMQ047NM9QxCLacBgbwhikX8zhT8OOM7AKSISUJU6CuiVIYiFFgCyG4JYaAHC33fPEMRCCwD2hiAWSXhTeCTgOAOniEhgrFJHAX0yBLGIAEKR1xDEQgsQtnXPEMQiAgilrSGIRRLeFNv+D0XkHSJyr4jc3lh3joh8VES+EMvHxPUiIn8hIsdE5LMictF2x3ccZ3eZJhJ4F/CXwA2NddcAt6nqG0Tkmli/mpCa/IL476eAt8bytKhKHQX0yhDEIgKI9a4lBw1l9wxBLCKAsN7WEMQiCW+KbTsBVf0XETlwyurLgUvi8vXAxwmdwOXADaqqwCdE5GwR2aeqJ057jrHUHUCvDEEMXgNCfkMQEy0AdNIQxOI1INgbglgk4U2x04HBx1cNO5aPi+vPBe5u7Hc8rnsYInKViBwVkaOb33tgh5fhOM68WA8MSsu61p9DVT1ESGXO6oH9WkUBbgjCrhuCWGgBwA1BLA1BLJLwpthpJPB1EdkHEMt74/rjwHmN/fYD9+zwHI7jZGCnkcCtwEHgDbG8pbH+lSJymDAg+J3txgMAUKmjADcE2X1DEAstALghiKUhiEUS3hTbdgIichNhEPCxInIceB2h8d8sIlcCdwEvirt/CLgMOAZ8H3j5dscHYDyZB+CGIOy6IYiFFgDcEMTSEMQiCW+Kad4OvCSx6dKWfRV4xXbHdBynHIqYMSg6mQfghiC7bwhioQUANwSxNASxSMKbwrUDjjNwiogE0MlkIDcEYdcNQSy0AOCGIJaGIBZJeFN4JOA4A6eISEDGk8lAvXIFstACQHZXIAstAHTTFchCCwD2rkAWSXhTFNEJoJMOoFeGIBZaAMhuCGKhBYBuGoJYaAHCPraGIBZJeFP444DjDJwiIgHRSRTQJ0MQCy0A5DcEsdACQDcNQSy0AGBvCGKRhDeFRwKOM3CKiAQYT6KAXhmCWGgBILshiIUWALppCGKhBQB7QxCLJLwpiugEqseBvhmCWGgBQj2zIYiBFgC6aQhioQUIpa0hiEUS3hT+OOA4A6eISCC8ItT+GYIYaAEgvyGIhRYAumkIYqEFAHtDEIskvCk8EnCcgVNEJCBVJNAzQxALLUBVdi05aKh3zxDEQgsA9oYgFkl4U5TRCYwnHUCvDEEMpgFXZdeSg0I3DUEspgGDvSGIRRLeFP444DgDp4hIIAwMjvtnCGKgBYD8hiAWWgDopiGIhRYA7A1BLJLwpvBIwHEGThGRQLAXG/fOEMRCCxDqeQ1BLLQAVb1rhiAWWgCwNwSxSMKbwiMBxxk4RUQCjLWOAnplCGKgBYD8hiAWWgDopiGIhRYA7A1BLJLwpiijE1CtO4BeGYIYaAGqsmvJQaGbhiAWWgCwNwSxSMKbwh8HHGfgFBEJiDYmA/XIEMRCCxC25TUEsdACQDcNQSy0AGBvCGKRhDfFtndaRN4hIveKyO2NdX8qIv8tIp8Vkb8VkbMb264VkWMicqeIPGe74zuOs7tMEwm8C/hL4IbGuo8C16rqhoi8EbgWuFpELgSuAH4MeCLwTyLyw6p6+q5ItY4CemUIYqAFCPW8hiAWWgDopiGIhRYA7A1BLJLwppgmF+G/iMiBU9Z9pFH9BPDCuHw5cFhVTwJfFpFjwNOBf9vmJJMOoEeGIBZagKreteSg0E1DEAstANgbgpgk4U1gMTD4a8A/xOVzgbsb247HdQ9DRK4SkaMicnRt80GDy3AcZyfMNTAoIq8FNoAbq1Utu7X+LKvqIeAQwKNXn6BVFNAnQxALLQDkNwSx0AJANw1BLLQAYG8IYpGEN8WOOwEROQg8H7g0piSH8Mt/XmO3/cA9Oz2H4ziLZ0edgIg8F7ga+FlV/X5j063Au0XkTYSBwQuAf9/2gHFMoG+GIBZagLBPXkMQCy0AdNMQxEILAPaGICZJeBNs2wmIyE3AJcBjReQ48DrC24BV4KMSPvBPqOpvqOrnRORm4POEx4RXbPtmAEB10gH0yBDEYhpw2CevIYjFNGDopiGIxTRgsDcEsci/mWKatwMvaVl93Wn2fz3w+u2O6zhOGRQxYxC0jgL6ZAhioQUIf5/XEMRCCxDK7hmCWGgBwN4QxCQJbwLXDjjOwCkjElAmUUCPDEEstABhW15DEAstAHTTEMRCCwD2hiAmSXgTeCTgOAOnkEhA6yigT4YgFlqAUOY1BLHQAoSye4YgFloAsDcEsUjCm6KMTgAmHUCPDEFMtACQ3RDEQgsA3TQEsdACgL0hiEkS3gT+OOA4A0d0lsk2i7oIkfuAB4Bv7Pa1AI/Fr6OJX8dWunwdP6SqP3jqyiI6AQAROaqqT/Pr8Ovw68h7Hf444DgDxzsBxxk4JXUCh3b7AiJ+HVvx69hK766jmDEBx3F2h5IiAcdxdgHvBBxn4BTRCYjIc2OegmMick2mc54nIh8TkTtE5HMi8uq4/hwR+aiIfCGWj8l0PUsi8mkR+WCsny8iR+J1vEdEVjJcw9ki8r6YU+IOEXnGbtwPEfmd+JncLiI3icjeXPcjkWej9R5I4C/i9/azInLRgq9jMfk+VHVX/wFLwBeBJwErwH8CF2Y47z7gorj8KOB/gAuBPwGuieuvAd6Y6T68Bng38MFYvxm4Ii6/DfjNDNdwPfDrcXkFODv3/SC4U38ZeETjPrws1/0AngVcBNzeWNd6D4DLCE7bAlwMHFnwdfwCsCcuv7FxHRfGdrMKnB/b09LU51r0F2uK/+wzgA836tcSEpvkvo5bgJ8H7gT2xXX7gDsznHs/cBvwbOCD8Uv1jcYHvuUeLegazoqNT05Zn/V+MLGtP4egbfkg8Jyc9wM4cErja70HwF8DL2nbbxHXccq2XwJujMtb2gzwYeAZ056nhMeBqXMVLIqYXOWpwBHg8ap6AiCWj5/g+30AAAIhSURBVMtwCW8Gfh+olB4/AHxbVSsheI578iTgPuCd8bHk7SJyBpnvh6p+Ffgz4C7gBPAd4FPkvx9NUvdgN7+7O8r30UYJncDUuQoWcnKRM4H3A7+tqt/Ndd7G+Z8P3Kuqn2qubtl10fdkDyH8fKuqPpWg5cgyPtMkPm9fTghrnwicATyvZdcS3m3vynd3nnwfbZTQCexargIRWSZ0ADeq6gfi6q+LyL64fR9w74Iv45nAC0TkK8BhwiPBm4GzRaSSeue4J8eB46p6JNbfR+gUct+PnwO+rKr3qeo68AHgp8l/P5qk7kH2724j38dLNcb+815HCZ3AJ4EL4ujvCiGh6a2LPqkEr/TrgDtU9U2NTbcCB+PyQcJYwcJQ1WtVdb+qHiD83/9ZVV8KfIxJjscc1/E14G4ReXJcdSnBOj7r/SA8BlwsIo+Mn1F1HVnvxymk7sGtwK/GtwQXA9+pHhsWQSPfxwv04fk+rhCRVRE5n2nzfVQscpBnhgGQywij818EXpvpnD9DCJk+C3wm/ruM8Dx+G/CFWJ6T8T5cwuTtwJPiB3kMeC+wmuH8Pwkcjffk74DH7Mb9AP4I+G/gduBvCKPeWe4HcBNhLGKd8At7ZeoeEMLwv4rf2/8Cnrbg6zhGePavvq9va+z/2ngddwLPm+VcPm3YcQZOCY8DjuPsIt4JOM7A8U7AcQaOdwKOM3C8E3CcgeOdgOMMHO8EHGfg/D/weEPH8JYGZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "prefix = options['prefix']\n",
    "chiN = options['chiN']\n",
    "eps = options['eps']\n",
    "sigma = options['sigma']\n",
    "it = options['it']\n",
    "dt = options['dt']\n",
    "\n",
    "lv = int(it/100)\n",
    "         \n",
    "io = writer(prefix,chiN,eps,sigma)\n",
    "its = np.arange(0,it,lv)\n",
    "\n",
    "\n",
    "t = np.arange(0,it/lv)*dt*lv\n",
    "\n",
    "# arr = io.get_dataset('energy', 1000)\n",
    "# print(arr)\n",
    "\n",
    "energy = []\n",
    "for it in its:\n",
    "    arr = io.get_dataset('energy',it)\n",
    "    energy.append(arr)\n",
    "    \n",
    "avg_conc = []\n",
    "for it in its:\n",
    "    arr = io.get_dataset('avg_conc',it)\n",
    "    avg_conc.append(arr)\n",
    "    \n",
    "print(energy)\n",
    "    \n",
    "fig2 = plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.plot(t,avg_conc)\n",
    "plt.plot(t,energy)\n",
    "# plt.ylim([0,1])\n",
    "plt.show()\n",
    "\n",
    "order_param = io.get_dataset('order_parameter',options['it']-1)\n",
    "fig = plt.figure()\n",
    "plt.imshow(order_param)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}