Keras_Galaxies.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_005.fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:316: RuntimeWarning: divide by zero encountered in log10\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_030.fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/matplotlib/pyplot.py:537: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  max_open_warning, RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_060.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_090.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_100.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_120.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_150.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_180.fits\n",
      "m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_200.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_005.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_040.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_060.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_100.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_200.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_300.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_400.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_410.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_420.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_600.fits\n",
      "Outside of the box\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_700.fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:8: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:9: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:10: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:10: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:12: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:13: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:14: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outside of the box\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_800.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_900.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_1000.fits\n",
      "q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_1100.fits\n"
     ]
    }
   ],
   "source": [
    "##### q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1\n",
    "##005,010,020,030,040,060,080,100,120,140,160\n",
    "#q0.5_fg0.3_allrx10_sunruns/hires_kin\n",
    "##170,180,185,190,195,205,210,220,225,230,240,250,260\n",
    "\n",
    "#m0.5_fg0.3_allrx10_isolated_sunruns\n",
    "##005,200\n",
    "#m1_fg0.3_allrx10_isolated_sunruns\n",
    "##005,100,200\n",
    "\n",
    "viewpt_list=[5]\n",
    "import pyfits\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy import shape\n",
    "from scipy import ndimage\n",
    "import scipy.ndimage.filters as filters\n",
    "from scipy.ndimage.morphology import generate_binary_structure, binary_erosion\n",
    "from scipy.ndimage import iterate_structure\n",
    "from scipy.ndimage.filters import maximum_filter\n",
    "import matplotlib\n",
    "import scipy.optimize as opt\n",
    "from astropy.cosmology import WMAP9 as cosmo\n",
    "import os\n",
    "from astropy.convolution import convolve, convolve_fft\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.convolution import Gaussian2DKernel\n",
    "from astropy.convolution import Gaussian1DKernel\n",
    "from astropy.convolution import convolve\n",
    "\n",
    "\n",
    "\n",
    "img_list=['m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_030.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_060.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_100.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_200.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_005.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_060.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_150.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_210.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_240.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_270.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_300.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_320.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_340.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_360.fits',\n",
    "          'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_400.fits',\n",
    "         'q0.2_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_420.fits']\n",
    "myr=[30,60,100,200,5,60,150,210,240,270,300,320,340,360,400,420]\n",
    "merger=[0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1]\n",
    "merger_sim=['fg3_m_15_iso_1','fg3_m_15_iso_1','fg3_m_15_iso_1','fg3_m_15_iso_1',\n",
    "        'fg3_m_15','fg3_m_15','fg3_m_15','fg3_m_15',\n",
    "            'fg3_m_15','fg3_m_15','fg3_m_15','fg3_m_15',\n",
    "            'fg3_m_15','fg3_m_15','fg3_m_15','fg3_m_15']\n",
    "\n",
    "\n",
    "img_list=['q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_180.fits',\n",
    "          'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_185.fits','q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_190.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_195.fits','q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_205.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_210.fits','q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_220.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_225.fits','q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_230.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_240.fits','q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_250.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin/broadband_260.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_005.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_010.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_020.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_030.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_040.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_060.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_080.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_100.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_120.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_140.fits',\n",
    "         'q0.5_fg0.3_allrx10_sunruns/hires_kin_early_cen1/broadband_160.fits',\n",
    "         'isolated_galaxies/m0.5_fg0.3/broadband_005.fits','isolated_galaxies/m0.5_fg0.3/broadband_010.fits','isolated_galaxies/m0.5_fg0.3/broadband_030.fits','isolated_galaxies/m0.5_fg0.3/broadband_050.fits',\n",
    "         'isolated_galaxies/m1_fg0.3/broadband_005.fits','isolated_galaxies/m1_fg0.3/broadband_010.fits','isolated_galaxies/m1_fg0.3/broadband_020.fits','isolated_galaxies/m1_fg0.3/broadband_030.fits','isolated_galaxies/m1_fg0.3/broadband_040.fits','isolated_galaxies/m1_fg0.3/broadband_050.fits','isolated_galaxies/m1_fg0.3/broadband_060.fits','isolated_galaxies/m1_fg0.3/broadband_100.fits']\n",
    "myr=[180,185, 190, 195, 205, 210, 220, 225, 230, 240, 250, 260,\n",
    "     5,10,20,30,40,60,80,100,120,140,160,\n",
    "    5,10,30,50,\n",
    "    5,10,20,30,40,50,60,100]\n",
    "merger=[1,1,1,1,1,1,1,1,1,1,1,1,\n",
    "       1,1,1,1,1,1,1,1,1,1,1,\n",
    "       0,0,0,0,\n",
    "       0,0,0,0,0,0,0,0]\n",
    "merger_sim=['fg3_m12','fg3_m12','fg3_m12',\n",
    "           'fg3_m12','fg3_m12','fg3_m12','fg3_m12',\n",
    "           'fg3_m12','fg3_m12','fg3_m12','fg3_m12',\n",
    "            'fg3_m12','fg3_m12','fg3_m12','fg3_m12',\n",
    "           'fg3_m12','fg3_m12','fg3_m12','fg3_m12',\n",
    "           'fg3_m12','fg3_m12','fg3_m12','fg3_m12',\n",
    "           'iso_m0.5_fg0.3','iso_m0.5_fg0.3','iso_m0.5_fg0.3','iso_m0.5_fg0.3',\n",
    "           'iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3','iso_m1_fg0.3']\n",
    "\n",
    "\n",
    "img_list=['m1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_005.fits',\n",
    "    'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_030.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_060.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_090.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_100.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_120.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_150.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_180.fits',\n",
    "          'm1_fg0.3_BT0.2_allrx10_isolated_sunruns/broadband_200.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_005.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_040.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_060.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_100.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_200.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_300.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_400.fits',\n",
    "          'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_410.fits',\n",
    "          'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_early_cen1/broadband_420.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_600.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_700.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_800.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_900.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_1000.fits',\n",
    "         'q0.1_fg0.3_BT0.2_allrx10_sunruns/hires_kin_late/broadband_1100.fits',]\n",
    "myr=[5,30,60,90,100,120,150,180,200,\n",
    "     5,40,60,100,200,300,400,410,420,\n",
    "     600,700,800,900,1000,1100]\n",
    "merger_sim=['fg3_m_10_iso_1','fg3_m_10_iso_1','fg3_m_10_iso_1','fg3_m_10_iso_1',\n",
    "       'fg3_m_10_iso_1','fg3_m_10_iso_1','fg3_m_10_iso_1','fg3_m_10_iso_1','fg3_m_10_iso_1',\n",
    "       'fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10',\n",
    "        'fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10',\n",
    "        'fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10','fg3_m_10']\n",
    "merger=[0,0,0,0,0,0,0,0,0,\n",
    "       1,1,1,1,1,1,1,1,1,\n",
    "       1,1,1,1,1,1]\n",
    "\n",
    "X=[]\n",
    "y=[]\n",
    "for i in range(len(img_list)):\n",
    "    print(img_list[i])\n",
    "    for j in range(len(viewpt_list)):\n",
    "        plott='no'\n",
    "        viewpt=viewpt_list[j]\n",
    "        im=pyfits.open(img_list[i])\n",
    "        \n",
    "        \n",
    "        a=produce_camera(myr[i],img_list[i],viewpt)\n",
    "        \n",
    "        band=65\n",
    "        try:\n",
    "            fill_img=a[1][band]\n",
    "        except IndexError:\n",
    "            band=2\n",
    "        b=determine_coords(a[1][band], plott)#was c, taking the second dimension\n",
    "        if b[0]==0:\n",
    "            continue\n",
    "        \n",
    "        \n",
    "        low_pass=b[4]\n",
    "        num_sol=b[5]\n",
    "\n",
    "        '''Now, fit a couple 2D gaussians if there are 2 brightest pixels, otherwise\n",
    "        fit only one 2D gaussian. The output of fit_2_gaussians will be the positions of these\n",
    "        maxima'''\n",
    "\n",
    "        if num_sol==1:\n",
    "            #this is if there's only really one solution because the bulges are too close together\n",
    "            c=fit_2_gaussian(b[1],shape(a[1][band])[0]-b[0],b[1],shape(a[1][band])[0]-b[0],low_pass, plott)\n",
    "            if c[8]=='no':\n",
    "                 c=fit_2_gaussian(b[2],shape(a[1][band])[0]-b[0],b[3],shape(a[1][band])[0]-b[1],low_pass, plott)#was d[3],300-d[2],d[1],300-d[0]\n",
    "\n",
    "        else:\n",
    "            c=fit_2_gaussian(b[2],shape(a[1][band])[0]-b[0],b[3],shape(a[1][band])[0]-b[1],low_pass, plott)#was d[3],300-d[2],d[1],300-d[0]\n",
    "\n",
    "        \n",
    "        if c[8] == 'no':\n",
    "            print('Did not fit the 2D Gaussian well')\n",
    "        if c[4] > c[5]:\n",
    "            '''this means point 1 is brighter'''\n",
    "            in_x = c[1]\n",
    "            in_y = c[0]\n",
    "            in_2_x = c[3]\n",
    "            in_2_y = c[2]\n",
    "            \n",
    "\n",
    "        if c[5] > c[4]:\n",
    "            '''point 2 is the brighter source'''\n",
    "            in_x = c[3]\n",
    "            in_y = c[2]\n",
    "            in_2_x = c[1]\n",
    "            in_2_y = c[0]\n",
    "\n",
    "            \n",
    "        \n",
    "        \n",
    "        '''Now place a aperture over each center and figure out which is brighter overall'''\n",
    "        c_final=determine_brighter(a[1][band],  in_y,shape(a[1][band])[0]-in_x,  in_2_y, shape(a[1][band])[0]-in_2_x, a[0])\n",
    "        \n",
    "        if c_final[2]=='yes':\n",
    "            continue\n",
    "\n",
    "        if c_final[0] > c_final[1]:\n",
    "            d=clip_image(a[1][band], a[0], 0.03, int(shape(a[1][band])[0]-in_x), int(in_y))\n",
    "            #this means the first aperture (center) is indeed brighter\n",
    "        else:\n",
    "            d=clip_image(a[1][band], a[0], 0.03, int(shape(a[1][band])[0]-in_2_x), int(in_2_y))\n",
    "\n",
    "        if d[3]=='no':\n",
    "            #this means the cutout dis outside the image --> bad\n",
    "            continue\n",
    "        \n",
    "        \n",
    "        e=nanomags(0.03, a[0], d[0],viewpt, myr[i])#was a[1]\n",
    "        \n",
    "        texp=1\n",
    "        prep=convolve_image(myr[i],e[0],e[1],e[2],0.03,a[0],viewpt, e[3], e[4], e[5], e[6], texp)\n",
    "        \n",
    "        X.append(prep[3]/np.max(prep[3]))\n",
    "        y.append(merger[i])\n",
    "        \n",
    "        plt.clf()\n",
    "        fig=plt.figure()\n",
    "        ax1=fig.add_subplot(111)\n",
    "        \n",
    "        im1=ax1.imshow(abs(prep[0]),norm=matplotlib.colors.LogNorm(vmin=10**(-1), vmax=10**1),cmap='afmhot')\n",
    "        #print(mask_log_image)\n",
    "        #vmin=d[0].min(),vmax=max_img,\n",
    "        \n",
    "        ax1.axis('off')\n",
    "        \n",
    "        #ax1.annotate(r't = '+str(myr[i]/100)+' Gyr', xy=(0.5,0.9), color='white', xycoords='axes fraction', size=15)\n",
    "        plt.savefig('../MaNGA_Papers/Paper_I/SDSS_image_plot_contrast_'+str(myr[i])+'_'+str(viewpt)+'_'+str(merger_sim[i])+'.png',bbox_inches='tight')\n",
    "        \n",
    "        \n",
    "        #plt.imshow(prep[3]/np.max(prep[3]))\n",
    "        #plt.colorbar()\n",
    "        #plt.show()\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiled\n"
     ]
    }
   ],
   "source": [
    "def twoD_two_Gaussian(xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset,\n",
    "                     amplitude_2, xo_2, yo_2, sigma_x_2, sigma_y_2, theta_2):\n",
    "    (x, y) = xdata_tuple \n",
    "    xo = float(xo)\n",
    "    yo = float(yo)   \n",
    "    xo_2 = float(xo_2)\n",
    "    yo_2 = float(yo_2)  \n",
    "    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)\n",
    "    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)\n",
    "    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)\n",
    "    \n",
    "    a_2 = (np.cos(theta_2)**2)/(2*sigma_x_2**2) + (np.sin(theta_2)**2)/(2*sigma_y_2**2)\n",
    "    b_2 = -(np.sin(2*theta_2))/(4*sigma_x_2**2) + (np.sin(2*theta_2))/(4*sigma_y_2**2)\n",
    "    c_2 = (np.sin(theta_2)**2)/(2*sigma_x_2**2) + (np.cos(theta_2)**2)/(2*sigma_y_2**2)\n",
    "    \n",
    "    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) \n",
    "                            + c*((y-yo)**2)))+ amplitude_2*np.exp( - (a_2*((x-xo_2)**2) + 2*b_2*(x-xo_2)*(y-yo_2) \n",
    "                            + c_2*((y-yo_2)**2)))\n",
    "    \n",
    "    return g.ravel()\n",
    "\n",
    "def detect_peaks(image):\n",
    "    \"\"\"\n",
    "    Takes an image and detect the peaks using the local maximum filter.\n",
    "    Returns a boolean mask of the peaks (i.e. 1 when\n",
    "    the pixel's value is the neighborhood maximum, 0 otherwise)\n",
    "    \"\"\"\n",
    "\n",
    "    # define an 8-connected neighborhood\n",
    "    struct = generate_binary_structure(2,1)\n",
    "    \n",
    "    neighborhood = iterate_structure(struct, 10).astype(bool)\n",
    "    \n",
    "    #apply the local maximum filter; all pixel of maximal value \n",
    "    #in their neighborhood are set to 1\n",
    "    local_max = maximum_filter(image, footprint=neighborhood)==image\n",
    "    #local_max is a mask that contains the peaks we are \n",
    "    #looking for, but also the background.\n",
    "    #In order to isolate the peaks we must remove the background from the mask.\n",
    "    \n",
    "    \n",
    "    #we create the mask of the background\n",
    "    background = (image==0)\n",
    "\n",
    "    #a little technicality: we must erode the background in order to \n",
    "    #successfully subtract it form local_max, otherwise a line will \n",
    "    #appear along the background border (artifact of the local maximum filter)\n",
    "    eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)\n",
    "\n",
    "    #we obtain the final mask, containing only peaks, \n",
    "    #by removing the background from the local_max mask (xor operation)\n",
    "    detected_peaks = local_max ^ eroded_background\n",
    "\n",
    "    return detected_peaks\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "def produce_camera(myr,img, viewpoint):\n",
    "    im=pyfits.open(img)\n",
    "\n",
    "    plt.clf()\n",
    "    \n",
    "\n",
    "    camera_data=im['CAMERA'+str(viewpoint)+'-BROADBAND'].data\n",
    "    pixelscale =  im['CAMERA'+str(viewpoint)+'-BROADBAND'].header['CD1_1']\n",
    "\n",
    "    \n",
    "\n",
    "    xs = np.linspace(0,shape(camera_data[0])[0]-1,shape(camera_data[0])[0])\n",
    "    xs_kpc = [pixelscale*x for x in xs]\n",
    "\n",
    "    \n",
    "    \n",
    "    return pixelscale, camera_data\n",
    "def determine_coords(img, plot):\n",
    "    \n",
    "    \n",
    "    kernel = np.ones((10,10))\n",
    "\n",
    "    lp_3 = ndimage.convolve(img, kernel)#was result\n",
    "    if plot=='yes':\n",
    "        plt.imshow(img)\n",
    "        plt.title('before filtering')\n",
    "        plt.colorbar()\n",
    "        plt.show()\n",
    "    \n",
    "    '''Okay here is where you can filter out low stuff'''\n",
    "    \n",
    "    max_value=(lp_3.max())\n",
    "    low = np.where(lp_3 < 0.2*max_value)\n",
    "    \n",
    "   \n",
    "    lp_3[low] = 0\n",
    "    \n",
    " \n",
    "    \n",
    "    \n",
    "    kernel_high = -1*np.ones((3,3))\n",
    "    kernel_high[1,1]=8\n",
    "    \n",
    "    hp_lp = ndimage.convolve(lp_3, kernel_high)\n",
    "    \n",
    "    \n",
    "    '''sharpen'''\n",
    "\n",
    "    kernel_sharpen = np.zeros((3,3))\n",
    "    kernel_sharpen[0,1]=-1\n",
    "    kernel_sharpen[1,0]=-1\n",
    "    kernel_sharpen[1,1]=5\n",
    "    kernel_sharpen[1,2]=-1\n",
    "    kernel_sharpen[2,1]=-1\n",
    "    kernel_sharpen=kernel_sharpen\n",
    "    \n",
    "    hp_lp_sharp = ndimage.convolve(hp_lp, kernel_sharpen)\n",
    "    if plot=='yes':\n",
    "        pp.imshow(hp_lp_sharp)\n",
    "        pp.title('Sharpened')\n",
    "        pp.colorbar()\n",
    "        pp.show()\n",
    "    \n",
    "    data_max = filters.maximum_filter(hp_lp, 3)\n",
    "    #data_max = filters.maximum_filter(data, neighborhood_size)\n",
    "    maxima = (hp_lp == data_max)\n",
    "    data_min = filters.minimum_filter(lp_3, 10)\n",
    "    diff = ((data_max - data_min) > 100)\n",
    "    maxima[diff == 0] = 0\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    indices = np.where(detect_peaks(lp_3) == 1)#was hp_lp_sharp\n",
    "    #number_of_sols=len(indices[0])\n",
    "    number_of_sols=len(indices[0])\n",
    "    \n",
    "    \n",
    "    try:\n",
    "        return indices[0][0],indices[0][-1],indices[1][0],indices[1][-1], lp_3, number_of_sols\n",
    "    except IndexError:\n",
    "        return 0,0,0,0,lp_3,number_of_sols\n",
    "def fit_2_gaussian(x_1,y_1,x_2,y_2, data, plot):\n",
    "# Create x and y indices\n",
    "    data=np.flipud(data)\n",
    "    x = np.linspace(0, 299, 300)\n",
    "    y = np.linspace(0, 299, 300)\n",
    "    x, y = np.meshgrid(x, y)\n",
    "    \n",
    "\n",
    "    # add some noise to the data and try to fit the data generated beforehand\n",
    "    initial_guess = (20,x_1,y_1,7,7,0,10,20,x_2,y_2,7,7,0)\n",
    "    data=data.ravel()\n",
    "    \n",
    "   \n",
    "    \n",
    "    #And plot the results:\n",
    "    \n",
    "    data_here = twoD_two_Gaussian((x,y), *initial_guess)\n",
    "    if plot=='yes':\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "        ax.hold(True)\n",
    "        ax.imshow(data.reshape(300, 300), cmap=plt.cm.jet, origin='bottom',\n",
    "            extent=(x.min(), x.max(), y.min(), y.max()))\n",
    "        ax.contour(x, y, data_here.reshape(300, 300), 8, colors='w')\n",
    "        plt.show()\n",
    "    try:\n",
    "        popt, pcov = opt.curve_fit(twoD_two_Gaussian, (x, y), data, p0=initial_guess)\n",
    "        fit='yes'\n",
    "    except RuntimeError:\n",
    "        popt=[0,0,0,0,0,0,0,0,0,0,0,0,0]\n",
    "        fit='no'\n",
    "    data_fitted = twoD_two_Gaussian((x, y), *popt)\n",
    "    \n",
    " \n",
    "    if plot=='yes':\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "        ax.hold(True)\n",
    "        ax.imshow(data.reshape(300, 300), cmap=plt.cm.jet, origin='bottom',\n",
    "            extent=(x.min(), x.max(), y.min(), y.max()))\n",
    "        ax.contour(x, y, data_fitted.reshape(300, 300), 8, colors='w')\n",
    "        plt.show()\n",
    "    \n",
    "    return popt[1], popt[2], popt[8], popt[9], popt[0], popt[7], np.sqrt(popt[3]**2+popt[4]**2), np.sqrt(popt[10]**2+popt[11]**2), fit \n",
    "def determine_brighter(img, x, y, x2, y2, pix):\n",
    "    kpc_arcmin=cosmo.kpc_proper_per_arcmin(0.03)#insert the redshift to get the kpc/arcmin scaling\n",
    "\n",
    "    \n",
    "    ap_size=(3*(kpc_arcmin.value/60))/pix ###arcsec * kpc/arcsec   / kpc/pix\n",
    "    from photutils import CircularAperture,aperture_photometry\n",
    "    '''step 1: define the circular aperture'''\n",
    "    positions = [(x, y), (x2, y2)]\n",
    "    apertures = CircularAperture(positions, ap_size)\n",
    "    phot_table = aperture_photometry(img, apertures)\n",
    "    total_light_1=phot_table['aperture_sum'][0]\n",
    "    total_light_2=phot_table['aperture_sum'][1]\n",
    "\n",
    "    \n",
    "    masks = apertures.to_mask(method='center')\n",
    "    mask = masks[0]\n",
    "\n",
    "    image = mask.to_image(shape=((shape(img)[0], shape(img)[0])))\n",
    "    skip='no'\n",
    "    try:\n",
    "        data_cutout = mask.apply(img)\n",
    "        \n",
    "    except TypeError:\n",
    "        skip='yes'\n",
    "\n",
    "    '''plt.imshow(data_cutout, norm=matplotlib.colors.LogNorm())\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    plt.imshow(img, norm=matplotlib.colors.LogNorm())\n",
    "    plt.scatter(x,y, color='white')\n",
    "    plt.scatter(x2,y2, color='red')\n",
    "    plt.show()'''\n",
    "    return total_light_1, total_light_2, skip\n",
    "\n",
    "def clip_image(ins, pixelscale, redshift, xcen, ycen):\n",
    "    plt.clf()\n",
    "    '''I count 50\" to a side images'''\n",
    "    #print('pixelscale', pixelscale)\n",
    "    kpc_arcmin=cosmo.kpc_proper_per_arcmin(redshift)#insert the redshift  \n",
    "    #print(kpc_arcmin.value/60, 'kpc per arcsec')\n",
    "    '''Divide the pixelscale (kpc) by kpc/arcsec to get arcsec\n",
    "    size of pixels'''\n",
    "    size_a=pixelscale/(kpc_arcmin.value/60)\n",
    "    #print('size', size_a, 'arcsec pixels')\n",
    "    num_pix_half=int(25/size_a)\n",
    "    if xcen-num_pix_half < 0 or ycen-num_pix_half < 0 or xcen+num_pix_half > 300 or ycen+num_pix_half > 300:\n",
    "        print('Outside of the box')\n",
    "        clipped=0\n",
    "        tag='no'\n",
    "    \n",
    "    else:\n",
    "        clipped=(ins[xcen-num_pix_half:xcen+num_pix_half,ycen-num_pix_half:ycen+num_pix_half])\n",
    "        tag='yes'\n",
    "    \n",
    "    return clipped, size_a, num_pix_half, tag, xcen, ycen\n",
    "'''Now I have to convert the units of LAURAS sims into nanomaggies and AB mags (mags of DR7 and DR13)\n",
    "'''\n",
    "def nanomags(z, pixscale, camera_data, view, number):\n",
    "#first, go from specific intensity (because we have per steradian) to \n",
    "#janskies (W/Hz/m^2)\n",
    "    c = 299792.458*1000#to get into m/s\n",
    "\n",
    "#this 1.35e-6 comes from the arcsin(R_sky/Distance to object)\n",
    "#the answer needs to be in radians\n",
    "\n",
    "#J=10^-26 W/m^2/Hz, so units of flux density\n",
    "#reference site: http://www.cv.nrao.edu/course/astr534/Brightness.html\n",
    "#We need to go from a spectral brightness (I_nu) which is in m units\n",
    "#To a flux density (S_nu) which is in units of Janskies (W/m^2/Hz)\n",
    "\n",
    "#So you need to multiply the Fν by c / λ^2 to convert it into Fλ. \n",
    "#But we are not done yet! Recalling from above, the units of Fλ \n",
    "#are not an energy density. You need to get another factor of λ \n",
    "#in there to make the units work out to be energy density: \n",
    "#calculate λFλ to get units of ergs/s/cm^2.\n",
    "\n",
    "#1.35e-6 comes from: sin^-1(pixelscale=0.166666#kpc/pix/ comoving distance in kpc)\n",
    "    pixelscale=pixscale#0.166666#kpc/pix/\n",
    "    d_a = cosmo.angular_diameter_distance(z).value#, **fidcosmo)\n",
    "    d_co = cosmo.comoving_distance(z).value\n",
    "    d_lum = cosmo.luminosity_distance(z).value \n",
    "    #print(\"Angluar-diameter distance to z=0.03 is\", 1000*(d_a), 'kpc')\n",
    "    #print(\"Comoving distance to z=0.03 is\", 1000*((d_co)), 'kpc')\n",
    "    #print(\"Luminosity distance to z=0.03 is\", 1000*((d_lum)), 'kpc')\n",
    "\n",
    "\n",
    "\n",
    "    #here's a good review of all the different distances\n",
    "    #http://www.astro.ufl.edu/~guzman/ast7939/projects/project01.html\n",
    "\n",
    "    #masked=ma.masked_where(camera_data[2] <= 0, camera_data[2])\n",
    "    Janskies=np.array(10**(26)*camera_data*(pixelscale/(1000*d_lum))**2*np.pi*((6185.2*10**(-10))**2/c), dtype='>f4')\n",
    "    #A Jansky is 10−26 watts per square metre per hertz.\n",
    "\n",
    "\n",
    "\n",
    "    nanomaggy=Janskies/(3.631*10**(-6))\n",
    "\n",
    "    \n",
    "\n",
    "    #first, convert to counts\n",
    "    #dn=img/cimg+simg\n",
    "    cimg=0.005005225 #nanomaggies/count\n",
    "    gain=4.735\n",
    "    darkvar=1.1966\n",
    "    simg=121.19590411\n",
    "    counts=(nanomaggy)/cimg+simg\n",
    "    sigma_counts=np.sqrt(counts/gain+darkvar)\n",
    "    sigma_nanomags=sigma_counts*cimg\n",
    "    \n",
    "    nanomag_bg = counts*cimg\n",
    "\n",
    "\n",
    "\n",
    "    '''The sky resids are given by:''' \n",
    "    sky_resids_mine=cimg*np.random.normal(0.331132,5.63218,shape(nanomaggy))\n",
    "    sky_resids_mine_counts=np.random.normal(0.331132,5.63218,shape(nanomaggy))\n",
    "    d_image=(nanomaggy)+sky_resids_mine\n",
    "    degraded_image=d_image\n",
    "    degraded_image_counts=d_image/cimg\n",
    "    #degraded_image=ma.masked_where(d_image <= 0, d_image)#numpy.ma.masked_less(x, value, \n",
    "\n",
    "    #in nanomags\n",
    "    #poisson_plus = np.sqrt((np.flipud(nanomaggy)+ 0.004701377)/1367.29546)+sky_resids_mine\n",
    "    #poisson_plus_nanomags = poisson_plus*1367.29546\n",
    "    \n",
    "  \n",
    "    mag=22.5-2.5*np.log10(nanomaggy)\n",
    "\n",
    "    \n",
    "    #return n, sky resids added (nanomags), poisson noise, resids in nanomags, sky resids added (counts), poisson noise, resids in counts\n",
    "    return nanomaggy, degraded_image, sigma_nanomags, sky_resids_mine, degraded_image_counts, sigma_counts, sky_resids_mine_counts\n",
    "#was degraded image, sigma_nanomags\n",
    "\n",
    "def convolve_image(number,nanomaggy,degraded_image,sigma_nanomags,z,pixscale, view, background, counts, counts_sig, background_counts, t_exp):#all of these are in nanomags\n",
    "    #PSF = 1.61 arcsec\n",
    "    \n",
    "\n",
    "    \n",
    "\n",
    "    kpc_arcmin=cosmo.kpc_proper_per_arcmin(z)#insert the redshift to get the kpc/arcmin scaling\n",
    "\n",
    "    sigma=1.61/2.355#apparently the preimage sigma is not large :)\n",
    "    ##kpc/pix is the pixelscale\n",
    "    ##conversion factor is kpc/\"\n",
    "    #1.61 is the FWHM in arcsec of the psf\n",
    "    #pixelscale is kpc/pix\n",
    "    kernel_sigma_pix=(sigma*(kpc_arcmin.value/60))/pixscale\n",
    "  \n",
    "\n",
    "    gaussian_2D_kernel = Gaussian2DKernel(kernel_sigma_pix)#standard deviation in pixels\n",
    "\n",
    "\n",
    "\n",
    "    #result = convolve(np.sum(CAMERA0.data[:,:,:],axis=0), gaussian_2D_kernel)\n",
    "    result_nano = (convolve(nanomaggy, gaussian_2D_kernel))\n",
    "    result = (convolve(degraded_image, gaussian_2D_kernel))\n",
    "    result_bg = (convolve(background, gaussian_2D_kernel))\n",
    "    result_bg_counts = (convolve(background_counts, gaussian_2D_kernel))\n",
    "    result_error = (convolve(sigma_nanomags, gaussian_2D_kernel))\n",
    "    \n",
    "    result_counts = (convolve(counts, gaussian_2D_kernel))\n",
    "    result_error_counts = (convolve(counts_sig, gaussian_2D_kernel))\n",
    "\n",
    "   \n",
    "    \n",
    "    \n",
    "    \n",
    "   \n",
    "    #def convolve_image(number,nanomaggy,degraded_image,sigma_nanomags,z,pixscale, view, background, counts, counts_sig, background_counts, t_exp):#all of these are in nanomags\n",
    "    #return degraded image, bg, nanomaggy\n",
    "    return result, result_bg, result_nano, result_counts, result_error_counts, result_bg_counts\n",
    "    \n",
    "\n",
    "print('Compiled')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'list'> <class 'list'>\n",
      "(22, 182, 182)\n",
      "(22,)\n",
      "<class 'list'> <class 'list'>\n",
      "(17, 182, 182)\n",
      "(17,)\n",
      "shape bf (17,)\n",
      "[1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1]\n",
      "shape af (17, 2)\n",
      "[[0. 1.]\n",
      " [0. 1.]\n",
      " [0. 1.]\n",
      " [1. 0.]\n",
      " [1. 0.]\n",
      " [0. 1.]\n",
      " [1. 0.]\n",
      " [1. 0.]\n",
      " [1. 0.]\n",
      " [0. 1.]\n",
      " [0. 1.]\n",
      " [1. 0.]\n",
      " [0. 1.]\n",
      " [1. 0.]\n",
      " [1. 0.]\n",
      " [0. 1.]\n",
      " [0. 1.]]\n",
      "(17, 182, 182, 1)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "print(type(X), type(y))\n",
    "print(np.shape(X))\n",
    "print(np.shape(y))\n",
    "\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)\n",
    "    \n",
    "print(type(X_train), type(y_train))\n",
    "print(np.shape(X_train))\n",
    "print(np.shape(y_train))\n",
    "\n",
    "X_train = np.array(X_train)\n",
    "X_test = np.array(X_test)\n",
    "\n",
    "X_train = X_train.reshape(X_train.shape[0],  X_train.shape[1], X_train.shape[2], 1)\n",
    "X_test = X_test.reshape(X_test.shape[0],   X_test.shape[1], X_test.shape[2], 1)\n",
    "X_train = X_train.astype('float32')\n",
    "X_test = X_test.astype('float32')\n",
    "\n",
    "print('shape bf',np.shape(y_train))\n",
    "print(y_train)\n",
    "\n",
    "y_train = np_utils.to_categorical(y_train, 2)\n",
    "y_test = np_utils.to_categorical(y_test, 2)\n",
    "print('shape af',y_train.shape)\n",
    "print(y_train)\n",
    "\n",
    "print(X_train.shape)#print(y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(None, 180, 180, 32)\n"
     ]
    }
   ],
   "source": [
    "model = Sequential()\n",
    " \n",
    "model.add(Convolution2D(32, (3, 3), activation='relu', input_shape=( X_train.shape[1], X_train.shape[2], 1), data_format='channels_last'))\n",
    "print(model.output_shape)\n",
    "\n",
    "\n",
    "model.add(Convolution2D(32, (3, 3), activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2,2)))\n",
    "model.add(Dropout(0.25))\n",
    " \n",
    "model.add(Flatten())\n",
    "model.add(Dense(128, activation='relu'))\n",
    "model.add(Dropout(0.5))\n",
    "model.add(Dense(2, activation='softmax'))\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/keras/models.py:981: UserWarning: The `nb_epoch` argument in `fit` has been renamed `epochs`.\n",
      "  warnings.warn('The `nb_epoch` argument in `fit` '\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/10\n",
      "17/17 [==============================] - 3s 171ms/step - loss: 0.6955 - acc: 0.4118\n",
      "Epoch 2/10\n",
      "17/17 [==============================] - 2s 89ms/step - loss: 0.5233 - acc: 0.7059\n",
      "Epoch 3/10\n",
      "17/17 [==============================] - 2s 89ms/step - loss: 0.4090 - acc: 0.7647\n",
      "Epoch 4/10\n",
      "17/17 [==============================] - 1s 81ms/step - loss: 0.7926 - acc: 0.6471\n",
      "Epoch 5/10\n",
      "17/17 [==============================] - 1s 76ms/step - loss: 0.2459 - acc: 0.9412\n",
      "Epoch 6/10\n",
      "17/17 [==============================] - 1s 74ms/step - loss: 0.2100 - acc: 0.8824\n",
      "Epoch 7/10\n",
      "17/17 [==============================] - 1s 77ms/step - loss: 0.0964 - acc: 1.0000\n",
      "Epoch 8/10\n",
      "17/17 [==============================] - 1s 79ms/step - loss: 0.0412 - acc: 1.0000\n",
      "Epoch 9/10\n",
      "17/17 [==============================] - 1s 76ms/step - loss: 0.0674 - acc: 1.0000\n",
      "Epoch 10/10\n",
      "17/17 [==============================] - 1s 78ms/step - loss: 0.0384 - acc: 1.0000\n",
      "5/5 [==============================] - 0s 46ms/step\n"
     ]
    }
   ],
   "source": [
    "# 8. Compile model\n",
    "model.compile(loss='categorical_crossentropy',\n",
    "              optimizer='adam',\n",
    "              metrics=['accuracy'])\n",
    " \n",
    "# 9. Fit model on training data\n",
    "model.fit(X_train, y_train, \n",
    "          batch_size=32, nb_epoch=10, verbose=1)\n",
    " \n",
    "# 10. Evaluate model on test data\n",
    "score = model.evaluate(X_test, y_test, verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.009745675139129162, 1.0] loss accuracy\n",
      "Saved trained model at path \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/matplotlib/pyplot.py:537: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  max_open_warning, RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.01044631 0.9895537 ]\n",
      "(182, 182)\n",
      "5\n",
      "[0. 1.]\n",
      "[2.1732468e-04 9.9978274e-01]\n",
      "(182, 182)\n",
      "5\n",
      "[0. 1.]\n",
      "[2.570763e-04 9.997429e-01]\n",
      "(182, 182)\n",
      "5\n",
      "[0. 1.]\n",
      "[2.0246756e-04 9.9979752e-01]\n",
      "(182, 182)\n",
      "5\n",
      "[0. 1.]\n",
      "[0.96314615 0.03685388]\n",
      "(182, 182)\n",
      "5\n",
      "[1. 0.]\n"
     ]
    }
   ],
   "source": [
    "print(score, 'loss','accuracy')\n",
    "model_path='path'\n",
    "model.save(model_path)\n",
    "print('Saved trained model at %s ' % model_path)\n",
    "plt.clf()\n",
    "fig, axs = plt.subplots(2,3, figsize=(9, 6), facecolor='w', edgecolor='k')\n",
    "fig.subplots_adjust(hspace = .1, wspace=.1)\n",
    "\n",
    "axs = axs.ravel()\n",
    "for j in range(len(X_test)):\n",
    "    print(model.predict(X_test)[j])\n",
    "    print(shape(X_test[j].reshape(X_test[j].shape[0],X_test[j].shape[1])))\n",
    "    print(len(X_test))\n",
    "    axs[j].imshow(abs(X_test[j].reshape(X_test[j].shape[0],X_test[j].shape[1])),norm=matplotlib.colors.LogNorm(vmin=0.01))#norm=matplotlib.colors.LogNorm(vmax=0.8))\n",
    "    #axs[j].annotate('Prediction [Nonmerging Merging]', color='white', xy=((0.1,0.95)), xycoords='axes fraction')\n",
    "    \n",
    "    #axs[j].annotate('                 ['+str(round(model.predict(X_test)[j][0],6))+' '+str(round(model.predict(X_test)[j][1],6))+']', color='white', xy=((0.1,0.9)), xycoords='axes fraction')\n",
    "    axs[j].annotate('['+str(round(model.predict(X_test)[j][0],6))+' '+str(round(model.predict(X_test)[j][1],6))+']', color='white', xy=((0.1,0.9)), xycoords='axes fraction')\n",
    "    \n",
    "    axs[j].annotate(str(y_test[j]), color='white', xy=((0.1,0.8)), xycoords='axes fraction')\n",
    "    #plt.colorbar()\n",
    "    axs[j].set_xticks([])\n",
    "    axs[j].set_yticks([])\n",
    "    print(y_test[j])\n",
    "    if j == len(X_test):\n",
    "        axs[j].set_visible(False)\n",
    "fig.delaxes(axs[j+1])\n",
    "plt.savefig('../MaNGA_Papers/Paper_I/keras_CNN_q0.1_fg0.3_BT0.2.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADoBJREFUeJzt3X2MXOV1x/HfyXq9jo1JvHHYboiL\nHeMEiGlMOjIgLKCiuA5CMiiKiRVFDiFxmuCktK4EdavGrWjlVgmRQynS0ri2I95CAsJ/0CR0FUGi\nwpbFMeYtvJlNY7PsYjZgQ4i9Xp/+sdfRBnaeWc/cmTu75/uRVjtzz71zj6792zszz8x9zN0FIJ53\nFd0AgGIQfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivADQU1r5M6mW5vP0KxG7hII5bd6U4f9kE1k\n3ZrCb2YrJG2W1CLpP9x9U2r9GZqls+2iWnYJIKHHuye8btVP+82sRdJNkj4h6QxJq83sjGofD0Bj\n1fKaf6mk5919j7sflnSHpJX5tAWg3moJ/8mSfjXm/t5s2e8xs7Vm1mtmvcM6VMPuAOSp7u/2u3uX\nu5fcvdSqtnrvDsAE1RL+fZLmjbn/wWwZgEmglvA/ImmRmS0ws+mSPi1pRz5tAai3qof63P2Ima2T\n9CONDvVtcfcnc+sMQF3VNM7v7vdJui+nXgA0EB/vBYIi/EBQhB8IivADQRF+ICjCDwRF+IGgCD8Q\nFOEHgiL8QFCEHwiK8ANBEX4gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ERfiBoAg/\nEBThB4Ii/EBQhB8IivADQRF+IKiaZuk1sz5JByWNSDri7qU8mkJ+bFr6n7jl/XPruv9n/np+2drI\nzKPJbU9ZOJisz/yKJesv3zC9bG1n6c7ktvtH3kzWz75rfbJ+6l89nKw3g5rCn/kTd9+fw+MAaCCe\n9gNB1Rp+l/RjM3vUzNbm0RCAxqj1af8yd99nZidJut/MfuHuD45dIfujsFaSZmhmjbsDkJeazvzu\nvi/7PSjpHklLx1mny91L7l5qVVstuwOQo6rDb2azzGz2sduSlkt6Iq/GANRXLU/7OyTdY2bHHuc2\nd/9hLl0BqLuqw+/ueyR9LMdepqyW0xcl697Wmqy/dMF7k/W3zik/Jt3+nvR49U8/lh7vLtJ//WZ2\nsv4v/7YiWe8587aytReH30puu2ng4mT9Az/1ZH0yYKgPCIrwA0ERfiAowg8ERfiBoAg/EFQe3+oL\nb+TCjyfrN2y9KVn/cGv5r55OZcM+kqz//Y2fS9anvZkebjv3rnVla7P3HUlu27Y/PRQ4s7cnWZ8M\nOPMDQRF+ICjCDwRF+IGgCD8QFOEHgiL8QFCM8+eg7ZmXkvVHfzsvWf9w60Ce7eRqff85yfqeN9KX\n/t668Ptla68fTY/Td3z7f5L1epr8X9itjDM/EBThB4Ii/EBQhB8IivADQRF+ICjCDwRl7o0b0TzR\n2v1su6hh+2sWQ1eem6wfWJG+vHbL7hOS9ce+cuNx93TM9fv/KFl/5IL0OP7Ia68n635u+au7930t\nuakWrH4svQLeoce7dcCH0nOXZzjzA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQFcf5zWyLpEslDbr7\n4mxZu6Q7Jc2X1Cdplbv/utLOoo7zV9Iy933J+sirQ8n6i7eVH6t/8vwtyW2X/vNXk/WTbiruO/U4\nfnmP82+V9PaJ0K+T1O3uiyR1Z/cBTCIVw+/uD0p6+6lnpaRt2e1tki7LuS8AdVbta/4Od+/Pbr8s\nqSOnfgA0SM1v+PnomwZl3zgws7Vm1mtmvcM6VOvuAOSk2vAPmFmnJGW/B8ut6O5d7l5y91Kr2qrc\nHYC8VRv+HZLWZLfXSLo3n3YANErF8JvZ7ZIekvQRM9trZldJ2iTpYjN7TtKfZvcBTCIVr9vv7qvL\nlBiwz8nI/ldr2n74wPSqt/3oZ55K1l+5uSX9AEdHqt43isUn/ICgCD8QFOEHgiL8QFCEHwiK8ANB\nMUX3FHD6tc+WrV15ZnpE9j9P6U7WL/jU1cn67DsfTtbRvDjzA0ERfiAowg8ERfiBoAg/EBThB4Ii\n/EBQjPNPAalpsl/98unJbf9vx1vJ+nXXb0/W/2bV5cm6//w9ZWvz/umh5LZq4PTxEXHmB4Ii/EBQ\nhB8IivADQRF+ICjCDwRF+IGgKk7RnSem6G4+Q58/N1m/9evfSNYXTJtR9b4/un1dsr7olv5k/cie\nvqr3PVXlPUU3gCmI8ANBEX4gKMIPBEX4gaAIPxAU4QeCqjjOb2ZbJF0qadDdF2fLNkr6oqRXstU2\nuPt9lXbGOP/k4+ctSdZP3LQ3Wb/9Qz+qet+n/eQLyfpH/qH8dQwkaeS5PVXve7LKe5x/q6QV4yz/\nlrsvyX4qBh9Ac6kYfnd/UNJQA3oB0EC1vOZfZ2a7zWyLmc3JrSMADVFt+G+WtFDSEkn9kr5ZbkUz\nW2tmvWbWO6xDVe4OQN6qCr+7D7j7iLsflXSLpKWJdbvcveTupVa1VdsngJxVFX4z6xxz93JJT+TT\nDoBGqXjpbjO7XdKFkuaa2V5JX5d0oZktkeSS+iR9qY49AqgDvs+PmrR0nJSsv3TFqWVrPdduTm77\nrgpPTD/z4vJk/fVlrybrUxHf5wdQEeEHgiL8QFCEHwiK8ANBEX4gKIb6UJjv7U1P0T3Tpifrv/HD\nyfqlX72m/GPf05PcdrJiqA9ARYQfCIrwA0ERfiAowg8ERfiBoAg/EFTF7/MjtqPL0pfufuFT6Sm6\nFy/pK1urNI5fyY1DZyXrM+/trenxpzrO/EBQhB8IivADQRF+ICjCDwRF+IGgCD8QFOP8U5yVFifr\nz34tPdZ+y3nbkvXzZ6S/U1+LQz6crD88tCD9AEf7c+xm6uHMDwRF+IGgCD8QFOEHgiL8QFCEHwiK\n8ANBVRznN7N5krZL6pDkkrrcfbOZtUu6U9J8SX2SVrn7r+vXalzTFpySrL9w5QfK1jZecUdy20+e\nsL+qnvKwYaCUrD+w+Zxkfc629HX/kTaRM/8RSevd/QxJ50i62szOkHSdpG53XySpO7sPYJKoGH53\n73f3ndntg5KelnSypJWSjn38a5uky+rVJID8HddrfjObL+ksST2SOtz92OcnX9boywIAk8SEw29m\nJ0j6gaRr3P3A2JqPTvg37qR/ZrbWzHrNrHdYh2pqFkB+JhR+M2vVaPBvdfe7s8UDZtaZ1TslDY63\nrbt3uXvJ3UutasujZwA5qBh+MzNJ35H0tLvfMKa0Q9Ka7PYaSffm3x6AepnIV3rPk/RZSY+b2a5s\n2QZJmyR9z8yukvRLSavq0+LkN23+Hybrr/9xZ7J+xT/+MFn/8/fenazX0/r+9HDcQ/9efjivfev/\nJredc5ShvHqqGH53/5mkcvN9X5RvOwAahU/4AUERfiAowg8ERfiBoAg/EBThB4Li0t0TNK3zD8rW\nhrbMSm775QUPJOurZw9U1VMe1u1blqzvvDk9Rffc7z+RrLcfZKy+WXHmB4Ii/EBQhB8IivADQRF+\nICjCDwRF+IGgwozzH/6z9GWiD//lULK+4dT7ytaWv/vNqnrKy8DIW2Vr5+9Yn9z2tL/7RbLe/lp6\nnP5osopmxpkfCIrwA0ERfiAowg8ERfiBoAg/EBThB4IKM87fd1n679yzZ95Vt33f9NrCZH3zA8uT\ndRspd+X0Uadd/2LZ2qKBnuS2I8kqpjLO/EBQhB8IivADQRF+ICjCDwRF+IGgCD8QlLl7egWzeZK2\nS+qQ5JK63H2zmW2U9EVJr2SrbnD38l96l3SitfvZxqzeQL30eLcO+FD6gyGZiXzI54ik9e6+08xm\nS3rUzO7Pat9y929U2yiA4lQMv7v3S+rPbh80s6clnVzvxgDU13G95jez+ZLOknTsM6PrzGy3mW0x\nszlltllrZr1m1jusQzU1CyA/Ew6/mZ0g6QeSrnH3A5JulrRQ0hKNPjP45njbuXuXu5fcvdSqthxa\nBpCHCYXfzFo1Gvxb3f1uSXL3AXcfcfejkm6RtLR+bQLIW8Xwm5lJ+o6kp939hjHLO8esdrmk9HSt\nAJrKRN7tP0/SZyU9bma7smUbJK02syUaHf7rk/SlunQIoC4m8m7/zySNN26YHNMH0Nz4hB8QFOEH\ngiL8QFCEHwiK8ANBEX4gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiCoipfuznVnZq9I+uWY\nRXMl7W9YA8enWXtr1r4keqtWnr2d4u7vn8iKDQ3/O3Zu1uvupcIaSGjW3pq1L4neqlVUbzztB4Ii\n/EBQRYe/q+D9pzRrb83al0Rv1Sqkt0Jf8wMoTtFnfgAFKST8ZrbCzJ4xs+fN7LoieijHzPrM7HEz\n22VmvQX3ssXMBs3siTHL2s3sfjN7Lvs97jRpBfW20cz2Zcdul5ldUlBv88zsJ2b2lJk9aWZ/kS0v\n9Ngl+irkuDX8ab+ZtUh6VtLFkvZKekTSand/qqGNlGFmfZJK7l74mLCZnS/pDUnb3X1xtuxfJQ25\n+6bsD+ccd7+2SXrbKOmNomduziaU6Rw7s7SkyyR9TgUeu0Rfq1TAcSvizL9U0vPuvsfdD0u6Q9LK\nAvpoeu7+oKShty1eKWlbdnubRv/zNFyZ3pqCu/e7+87s9kFJx2aWLvTYJfoqRBHhP1nSr8bc36vm\nmvLbJf3YzB41s7VFNzOOjmzadEl6WVJHkc2Mo+LMzY30tpmlm+bYVTPjdd54w++dlrn7xyV9QtLV\n2dPbpuSjr9maabhmQjM3N8o4M0v/TpHHrtoZr/NWRPj3SZo35v4Hs2VNwd33Zb8HJd2j5pt9eODY\nJKnZ78GC+/mdZpq5ebyZpdUEx66ZZrwuIvyPSFpkZgvMbLqkT0vaUUAf72Bms7I3YmRmsyQtV/PN\nPrxD0prs9hpJ9xbYy+9plpmby80srYKPXdPNeO3uDf+RdIlG3/F/QdLfFtFDmb4+JOmx7OfJonuT\ndLtGnwYOa/S9kaskvU9St6TnJP23pPYm6u27kh6XtFujQessqLdlGn1Kv1vSruznkqKPXaKvQo4b\nn/ADguINPyAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivADQf0/sEWOix6VKakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(60000, 1, 28, 28)\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:3: UserWarning: Update your `Conv2D` call to the Keras 2 API: `Conv2D(32, (3, 3), activation=\"relu\", input_shape=(1, 28, 28...)`\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Negative dimension size caused by subtracting 3 from 1 for 'conv2d_2/convolution' (op: 'Conv2D') with input shapes: [?,1,28,28], [3,3,28,32].",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mInvalidArgumentError\u001b[0m                      Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/framework/ops.py\u001b[0m in \u001b[0;36m_create_c_op\u001b[0;34m(graph, node_def, inputs, control_inputs)\u001b[0m\n\u001b[1;32m   1566\u001b[0m   \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1567\u001b[0;31m     \u001b[0mc_op\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mc_api\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTF_FinishOperation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mop_desc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1568\u001b[0m   \u001b[0;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInvalidArgumentError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mInvalidArgumentError\u001b[0m: Negative dimension size caused by subtracting 3 from 1 for 'conv2d_2/convolution' (op: 'Conv2D') with input shapes: [?,1,28,28], [3,3,28,32].",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-ab32396f1ebf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSequential\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mConvolution2D\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m32\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'relu'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minput_shape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m28\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m28\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mConvolution2D\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m32\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'relu'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMaxPooling2D\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpool_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/keras/models.py\u001b[0m in \u001b[0;36madd\u001b[0;34m(self, layer)\u001b[0m\n\u001b[1;32m    495\u001b[0m                 \u001b[0;31m# and create the node connecting the current layer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    496\u001b[0m                 \u001b[0;31m# to the input layer we just created.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 497\u001b[0;31m                 \u001b[0mlayer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    498\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    499\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlayer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_inbound_nodes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutput_tensors\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/keras/engine/topology.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inputs, **kwargs)\u001b[0m\n\u001b[1;32m    617\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    618\u001b[0m             \u001b[0;31m# Actually call the layer, collecting output(s), mask(s), and shape(s).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 619\u001b[0;31m             \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    620\u001b[0m             \u001b[0moutput_mask\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_mask\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprevious_mask\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    621\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/keras/layers/convolutional.py\u001b[0m in \u001b[0;36mcall\u001b[0;34m(self, inputs)\u001b[0m\n\u001b[1;32m    166\u001b[0m                 \u001b[0mpadding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpadding\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    167\u001b[0m                 \u001b[0mdata_format\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata_format\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 168\u001b[0;31m                 dilation_rate=self.dilation_rate)\n\u001b[0m\u001b[1;32m    169\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrank\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    170\u001b[0m             outputs = K.conv3d(\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py\u001b[0m in \u001b[0;36mconv2d\u001b[0;34m(x, kernel, strides, padding, data_format, dilation_rate)\u001b[0m\n\u001b[1;32m   3339\u001b[0m         \u001b[0mstrides\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstrides\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3340\u001b[0m         \u001b[0mpadding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpadding\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3341\u001b[0;31m         data_format=tf_data_format)\n\u001b[0m\u001b[1;32m   3342\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3343\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mdata_format\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'channels_first'\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mtf_data_format\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'NHWC'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/ops/nn_ops.py\u001b[0m in \u001b[0;36mconvolution\u001b[0;34m(input, filter, padding, strides, dilation_rate, name, data_format)\u001b[0m\n\u001b[1;32m    778\u001b[0m         \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    779\u001b[0m         data_format=data_format)\n\u001b[0;32m--> 780\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    781\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    782\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/ops/nn_ops.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inp, filter)\u001b[0m\n\u001b[1;32m    866\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    867\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# pylint: disable=redefined-builtin\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 868\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconv_op\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    869\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    870\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/ops/nn_ops.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inp, filter)\u001b[0m\n\u001b[1;32m    518\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    519\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# pylint: disable=redefined-builtin\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 520\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    521\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    522\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/ops/nn_ops.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inp, filter)\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0mpadding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpadding\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    203\u001b[0m         \u001b[0mdata_format\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata_format\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 204\u001b[0;31m         name=self.name)\n\u001b[0m\u001b[1;32m    205\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    206\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/ops/gen_nn_ops.py\u001b[0m in \u001b[0;36mconv2d\u001b[0;34m(input, filter, strides, padding, use_cudnn_on_gpu, data_format, dilations, name)\u001b[0m\n\u001b[1;32m    954\u001b[0m         \u001b[0;34m\"Conv2D\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfilter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstrides\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstrides\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    955\u001b[0m         \u001b[0mpadding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpadding\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0muse_cudnn_on_gpu\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0muse_cudnn_on_gpu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 956\u001b[0;31m         data_format=data_format, dilations=dilations, name=name)\n\u001b[0m\u001b[1;32m    957\u001b[0m     \u001b[0m_result\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_op\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutputs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    958\u001b[0m     \u001b[0m_inputs_flat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_op\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py\u001b[0m in \u001b[0;36m_apply_op_helper\u001b[0;34m(self, op_type_name, name, **keywords)\u001b[0m\n\u001b[1;32m    785\u001b[0m         op = g.create_op(op_type_name, inputs, output_types, name=scope,\n\u001b[1;32m    786\u001b[0m                          \u001b[0minput_types\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minput_types\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mattrs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mattr_protos\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 787\u001b[0;31m                          op_def=op_def)\n\u001b[0m\u001b[1;32m    788\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0moutput_structure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mop_def\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_stateful\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mop\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    789\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/framework/ops.py\u001b[0m in \u001b[0;36mcreate_op\u001b[0;34m(self, op_type, inputs, dtypes, input_types, name, attrs, op_def, compute_shapes, compute_device)\u001b[0m\n\u001b[1;32m   3390\u001b[0m           \u001b[0minput_types\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minput_types\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3391\u001b[0m           \u001b[0moriginal_op\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_default_original_op\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3392\u001b[0;31m           op_def=op_def)\n\u001b[0m\u001b[1;32m   3393\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3394\u001b[0m       \u001b[0;31m# Note: shapes are lazily computed with the C API enabled.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/framework/ops.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, node_def, g, inputs, output_types, control_inputs, input_types, original_op, op_def)\u001b[0m\n\u001b[1;32m   1732\u001b[0m           op_def, inputs, node_def.attr)\n\u001b[1;32m   1733\u001b[0m       self._c_op = _create_c_op(self._graph, node_def, grouped_inputs,\n\u001b[0;32m-> 1734\u001b[0;31m                                 control_input_ops)\n\u001b[0m\u001b[1;32m   1735\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1736\u001b[0m       \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_c_op\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/tensorflow/python/framework/ops.py\u001b[0m in \u001b[0;36m_create_c_op\u001b[0;34m(graph, node_def, inputs, control_inputs)\u001b[0m\n\u001b[1;32m   1568\u001b[0m   \u001b[0;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInvalidArgumentError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1569\u001b[0m     \u001b[0;31m# Convert to ValueError for backwards compatibility.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1570\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1571\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1572\u001b[0m   \u001b[0;32mreturn\u001b[0m \u001b[0mc_op\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Negative dimension size caused by subtracting 3 from 1 for 'conv2d_2/convolution' (op: 'Conv2D') with input shapes: [?,1,28,28], [3,3,28,32]."
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/beckynevin/anaconda/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAByBJREFUeJzt3V+s13Udx/Hz4xzATpJKpOlmeAIM\npqQVK5gM2hrkRRe1dmLOm2hdpKGraLNc69+o2dbciMgLN0O3LDuulhf9GWuNuSWnzGarljSBNYVO\nHc7AQgr8/X7d2FV93+iP8//1eNy++J7vd2NPvhcfzu/X6na7fUCeBTP9AMDMED+EEj+EEj+EEj+E\nEj+EEj+EEj+EEj+EGpjOm21ZMOy/E8IU298Zab2SP+fND6HED6HED6HED6HED6HED6HED6HED6HE\nD6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6HE\nD6HED6HED6HED6HED6HED6HED6HED6HED6HED6HED6EGZvoBmFqtgfqvuP8Ny6b0/s98+prGrT3Y\nKa9dvuJv5T54e6vc/3rvosbtqXWPlNeOt0+X+7tGdpb7yk8dLPfZwJsfQokfQokfQokfQokfQokf\nQokfQjnnnwb9a1aVe3fxwnI/tvnScj+zvvlMeukl9Xn14zfU590z6ScvLin3r33z5nIfXftw43bk\n3Jny2nvGtpT7VY93y30u8OaHUOKHUOKHUOKHUOKHUOKHUI76JkH73W8v93v37S33axc2/+rpfHau\n2y73z+/5cLkPnK6P2zaM7Gjcljz/Unnt4vH6KHDwydFynwu8+SGU+CGU+CGU+CGU+CGU+CGU+CGU\nc/5JsPiZY+X+m39dXe7XLhybzMeZVDuPry/3w/+sP/p734pHG7dTnfqc/opv/LLcp9Lc/4Xd8/Pm\nh1Dih1Dih1Dih1Dih1Dih1Dih1Ctbnf6TjS3LBhOOD79HxPbN5T7CzfXH6/d/7uLy/3p2/e86mf6\nr13jby33X2+uz/HbJ0+Ve3fDDY3b0TvLS/uGbnm6/gP8X/s7I/V3l7/Mmx9CiR9CiR9CiR9CiR9C\niR9CiR9COeefBfqXvb7c2ycmyv3Iw81n9X/Y9EB57Tu/eke5X7535n6nnt445wdK4odQ4odQ4odQ\n4odQ4odQ4odQPrd/FmiPn7ig68+9sKjna6+79Y/l/vf7+usf0Gn3fG9mljc/hBI/hBI/hBI/hBI/\nhBI/hHLUNw+suetQ47Z97XvKa7+9/Oflvnn44+W+5JGD5c7s5c0PocQPocQPocQPocQPocQPocQP\noZzzzwPV12SfuG1Nee1fHjtT7p/Z9VC5f/ZDHyj37m8vadyu/soT5bV90/ix8om8+SGU+CGU+CGU\n+CGU+CGU+CGU+CGUr+gON/GRDeX+nS98vdyHBi7q+d7XPbSj3Ffdf7zcXzp8tOd7z2e+ohsoiR9C\niR9CiR9CiR9CiR9CiR9COeen1L3pxnJ/3T3Plft33/yznu+9+hcfLfe3fKn5cwz6+vr62n8+3PO9\n5zLn/EBJ/BBK/BBK/BBK/BBK/BBK/BDKOT8XpP+Ky8v92LaVjdvoXbvLaxec591065Gt5X5q44ly\nn6+c8wMl8UMo8UMo8UMo8UMo8UMoR33MmO8/V39F92BrUbm/2D1b7u+74xPNP/uHo+W1c5mjPqAk\nfgglfgglfgglfgglfgglfgg1MNMPwOzW2Vh/dPezw/VXdF9/49HG7Xzn+OezZ+Jt5T74oycv6OfP\nd978EEr8EEr8EEr8EEr8EEr8EEr8EMo5/zzXWnd9uR+6sz5rv/+mB8t900X179RfiH93z5X7wYmh\n+gd0jk/i08w/3vwQSvwQSvwQSvwQSvwQSvwQSvwQyjn/HDAwtLzcn91+VeP2xW3fK6/94MXjPT3T\nZLh7bF25H9i9vtwve7D+3H9q3vwQSvwQSvwQSvwQSvwQSvwQylHfNBi45k3lfuodV5b7ti//tNw/\ndukPXvUzTZadx+vjuCe+1Xyct3Tfr8prL+s4yptK3vwQSvwQSvwQSvwQSvwQSvwQSvwQyjn/KzRw\n5Rsbt4kHXltee9vQgXK/ZclYT880GXY8v7Hcn7qv/oruZY/+vtyX/sNZ/WzlzQ+hxA+hxA+hxA+h\nxA+hxA+hxA+hYs75z763/pjos5+cKPe7V/64cdv6mtM9PdNkGWufadw2PbazvHb15/5U7ktP1uf0\nnXJlNvPmh1Dih1Dih1Dih1Dih1Dih1Dih1Ax5/xH31//O3do7ciU3XvvyRXlvvvA1nJvtVvlvnrX\nkcZt1dhoeW27XJnPvPkhlPghlPghlPghlPghlPghlPghVKvb7U7bzbYsGJ6+m0Go/Z2R+j+GvMyb\nH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJH0KJ\nH0KJH0KJH0JN60d3A7OHNz+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+E\nEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+EEj+E+g86JSIeTPOXvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADuNJREFUeJzt3X+QVfV5x/HPw3bll+hIDBtCSIkK\nUkobiBuMjQlJrA7YTNGZhoTpGEptyUyixWjbOLYzddKZDs2YWNNgUhKJmB+YzqiR6VCjbplaE0JY\nkIiKBkOWCiJEoAV/4S779I89pBvd872Xe8+95+4+79fMzt57nnPueebCZ8+993vO/Zq7C0A8o8pu\nAEA5CD8QFOEHgiL8QFCEHwiK8ANBEX4gKMIPBEX4gaB+o5k7O81G+xiNb+YugVBe08t63Y9bNevW\nFX4zWyDpNkltkr7h7itT64/ReF1ol9SzSwAJm72r6nVrftlvZm2SVklaKGmWpCVmNqvWxwPQXPW8\n558n6Vl33+3ur0u6W9KiYtoC0Gj1hH+KpOcG3d+bLfs1ZrbczLrNrLtXx+vYHYAiNfzTfndf7e6d\n7t7ZrtGN3h2AKtUT/n2Spg66/45sGYBhoJ7wb5E03czeZWanSfqEpPXFtAWg0Woe6nP3PjO7RtIP\nNDDUt8bdnyysMwANVdc4v7tvkLShoF4ANBGn9wJBEX4gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIrw\nA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivADQRF+ICjCDwRF+IGgCD8QVFOn6MbI0/eRC5L1\n/Z/On6LtpxetTW777k1Lk/W3rzotWW/buC1Zj44jPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8EVdc4\nv5n1SDom6YSkPnfvLKIptI7++XOT9S+v+Uqyfl57/n+x/gr7fuyibybrz3SeSNb/atr7KuwhtiJO\n8vmwu79YwOMAaCJe9gNB1Rt+l/SgmW01s+VFNASgOep92X+xu+8zs0mSHjKzp939kcErZH8UlkvS\nGI2rc3cAilLXkd/d92W/D0q6T9K8IdZZ7e6d7t7ZrtH17A5AgWoOv5mNN7MJJ29LukzSE0U1BqCx\n6nnZ3yHpPjM7+TjfdfcHCukKQMPVHH533y3p3QX2ghL0XpY+NeOvb/9Wsj6jPX1NfX9iNH93b29y\n2//tT79NnFvhXeTxhe/NrY3duCO5bf9rr6UffARgqA8IivADQRF+ICjCDwRF+IGgCD8QFF/dPQK0\nnXFGbu3lD85MbvvZW7+brH947EsV9l778ePOI7+XrHfdflGy/sObv5ysP/SNr+XWZn37muS253xu\nU7I+EnDkB4Ii/EBQhB8IivADQRF+ICjCDwRF+IGgGOcfAfbeNSW3tuW9q5rYyan5/KQtyfoDp6fP\nA1jWc1myvnbaw7m1M2YdSm4bAUd+ICjCDwRF+IGgCD8QFOEHgiL8QFCEHwiKcf5hoO8jFyTr6+bk\nT5M9Sumv1q5k2Z5LkvXuh38rWd9xdX5vG18dk9x2UveryfqzR9LfVdD+Dxtza6MsuWkIHPmBoAg/\nEBThB4Ii/EBQhB8IivADQRF+IChz9/QKZmskfVTSQXefnS2bKOl7kqZJ6pG02N2PVNrZGTbRL7T0\nuHFE/fPnJuv/tPb2ZP289tpP1/jDp69M1tv+6OVk/fAfnJ+sH5qdP6A+Y9VzyW37ntubrFfyb/u2\n5tb2n0ifQ/CnS/8iWW/buK2mnhpts3fpqB+u6iyGao78d0pa8IZlN0rqcvfpkrqy+wCGkYrhd/dH\nJB1+w+JFktZmt9dKuqLgvgA0WK3v+TvcfX92+wVJHQX1A6BJ6v7Azwc+NMj94MDMlptZt5l19+p4\nvbsDUJBaw3/AzCZLUvb7YN6K7r7a3TvdvbNdo2vcHYCi1Rr+9ZKWZreXSrq/mHYANEvF8JvZOkmb\nJJ1vZnvN7GpJKyVdama7JP1+dh/AMFJxgNjdl+SUGLCvkl3w28n6i9enx5xntKevyd+a+CjlP16a\nldz20N1Tk/W3HEnPU3/mt3+cridqfcktG6ujLf0W9NB1ryTrk/K/KmDY4Aw/ICjCDwRF+IGgCD8Q\nFOEHgiL8QFB8dXcBRo0bl6z3feFosv7jmfcm67/oez1Zv/6mG3JrZ/3Xfye3nTQ+9+RMSdKJZHXk\nmjd5T7Le05w2GoojPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ExTh/AV6dn75k9wcz01+9Xcmfrfhs\nsj7h+/mX1ZZ52SxaG0d+ICjCDwRF+IGgCD8QFOEHgiL8QFCEHwiKcf4C/O7fb0/WR1X4G7tsT/pb\n0Md+/yen3BOkdmvLrfWmZ6ZXm1VYYQTgyA8ERfiBoAg/EBThB4Ii/EBQhB8IivADQVUc5zezNZI+\nKumgu8/Olt0s6c8l/TJb7SZ339CoJlvB/1x1UW7tbztuSW7brwpTbD+Ynkb7nfpRso6h9Xr+rAP9\n6k9u+8DO9L/JdG2rqadWUs2R/05JC4ZYfqu7z8l+RnTwgZGoYvjd/RFJh5vQC4Amquc9/zVm9riZ\nrTGzswrrCEBT1Br+r0o6V9IcSfslfTFvRTNbbmbdZtbdq+M17g5A0WoKv7sfcPcT7t4v6euS5iXW\nXe3une7e2a7RtfYJoGA1hd/MJg+6e6WkJ4ppB0CzVDPUt07ShySdbWZ7Jf2dpA+Z2RxJroHZij/V\nwB4BNEDF8Lv7kiEW39GAXlpa39j82pmj0uP4m15Lv905567n0/tOVkeuUePGJetP3zK7wiNsza38\n8e6FyS1nrvhFsp5/BsHwwRl+QFCEHwiK8ANBEX4gKMIPBEX4gaD46u4mOHTi9GS9b3dPcxppMZWG\n8p5Z+TvJ+tOLvpKs//srZ+bWnl91XnLbCUfypz0fKTjyA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQ\njPM3wV/+8GPJ+ozEpafDXf/8ubm1g9e/mtx2Z2d6HP+SHR9P1scv2J1bm6CRP45fCUd+ICjCDwRF\n+IGgCD8QFOEHgiL8QFCEHwiKcf5qWX5pVIW/obddvC5ZX6UZtXTUEvZ8Pn/qckm655Nfyq3NaE9/\n5fl7frI0WX/7lU8l60jjyA8ERfiBoAg/EBThB4Ii/EBQhB8IivADQVUc5zezqZLuktQhySWtdvfb\nzGyipO9JmiapR9Jidz/SuFZL5vmlfvUnN50/9lCyft2dFyTr534z/fjtLxzLrR2Y/9bkthM/vjdZ\nv/adXcn6wnHp7yJY/3JHbu2TOxYktz37X8Yn66hPNUf+Pkk3uPssSe+T9BkzmyXpRkld7j5dUld2\nH8AwUTH87r7f3bdlt49J2ilpiqRFktZmq62VdEWjmgRQvFN6z29m0yTNlbRZUoe7789KL2jgbQGA\nYaLq8JvZ6ZLukXSdux8dXHN3V867YjNbbmbdZtbdq+N1NQugOFWF38zaNRD877j7vdniA2Y2OatP\nlnRwqG3dfbW7d7p7Z7tGF9EzgAJUDL+ZmaQ7JO1098GXaK2XdPKyq6WS7i++PQCNUs0lve+XdJWk\nHWa2PVt2k6SVkv7VzK6WtEfS4sa0OPyNsfTTvPPSryXrj35gTLK+6/jbcmvLzuxJbluvFc9/IFl/\n4EdzcmvTV/D12WWqGH53f1T5V7NfUmw7AJqFM/yAoAg/EBThB4Ii/EBQhB8IivADQdnAmbnNcYZN\n9AtteI4Ots04N7c2Y92e5Lb/+LZNde270leDV7qkOOWx4+nHXvKfy5P1GctG7vTiw9Fm79JRP5z4\novn/x5EfCIrwA0ERfiAowg8ERfiBoAg/EBThB4Jiiu4qnfjZz3Nruz42LbntrGuvTdafWvzPtbRU\nlZkbPp2sn3/7K8n6jMcYxx+pOPIDQRF+ICjCDwRF+IGgCD8QFOEHgiL8QFBczw+MIFzPD6Aiwg8E\nRfiBoAg/EBThB4Ii/EBQhB8IqmL4zWyqmW00s6fM7EkzW5Etv9nM9pnZ9uzn8sa3C6Ao1XyZR5+k\nG9x9m5lNkLTVzB7Kare6+y2Naw9Ao1QMv7vvl7Q/u33MzHZKmtLoxgA01im95zezaZLmStqcLbrG\nzB43szVmdlbONsvNrNvMunt1vK5mARSn6vCb2emS7pF0nbsflfRVSedKmqOBVwZfHGo7d1/t7p3u\n3tmu0QW0DKAIVYXfzNo1EPzvuPu9kuTuB9z9hLv3S/q6pHmNaxNA0ar5tN8k3SFpp7t/adDyyYNW\nu1LSE8W3B6BRqvm0//2SrpK0w8y2Z8tukrTEzOZIckk9kj7VkA4BNEQ1n/Y/Kmmo64M3FN8OgGbh\nDD8gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQTZ2i28x+\nKWnPoEVnS3qxaQ2cmlbtrVX7kuitVkX29pvu/tZqVmxq+N+0c7Nud+8srYGEVu2tVfuS6K1WZfXG\ny34gKMIPBFV2+FeXvP+UVu2tVfuS6K1WpfRW6nt+AOUp+8gPoCSlhN/MFpjZM2b2rJndWEYPecys\nx8x2ZDMPd5fcyxozO2hmTwxaNtHMHjKzXdnvIadJK6m3lpi5OTGzdKnPXavNeN30l/1m1ibpZ5Iu\nlbRX0hZJS9z9qaY2ksPMeiR1unvpY8Jm9kFJL0m6y91nZ8u+IOmwu6/M/nCe5e6fa5Hebpb0Utkz\nN2cTykwePLO0pCsk/YlKfO4SfS1WCc9bGUf+eZKedffd7v66pLslLSqhj5bn7o9IOvyGxYskrc1u\nr9XAf56my+mtJbj7fnfflt0+JunkzNKlPneJvkpRRvinSHpu0P29aq0pv13Sg2a21cyWl93MEDqy\nadMl6QVJHWU2M4SKMzc30xtmlm6Z566WGa+Lxgd+b3axu79H0kJJn8le3rYkH3jP1krDNVXN3Nws\nQ8ws/StlPne1znhdtDLCv0/S1EH335Etawnuvi/7fVDSfWq92YcPnJwkNft9sOR+fqWVZm4eamZp\ntcBz10ozXpcR/i2SppvZu8zsNEmfkLS+hD7exMzGZx/EyMzGS7pMrTf78HpJS7PbSyXdX2Ivv6ZV\nZm7Om1laJT93LTfjtbs3/UfS5Rr4xP/nkv6mjB5y+jpH0k+znyfL7k3SOg28DOzVwGcjV0t6i6Qu\nSbskPSxpYgv19i1JOyQ9roGgTS6pt4s18JL+cUnbs5/Ly37uEn2V8rxxhh8QFB/4AUERfiAowg8E\nRfiBoAg/EBThB4Ii/EBQhB8I6v8AG8x2aarNGp8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADV9JREFUeJzt3W+MXXWdx/HPp8O0tVUiU+zsCJWy\nCCaEZAczFlf+LJsiQcKmEE0jiW43IdYHkl0SH8B2d7MYH4hmFYkakhG6lo2Cu1FCHwACEyMhktoB\nKwWLgliW1tKpFtMipX+/PpiDGWDuubf3nnvPnX7fr6SZe8/vnHs+Oelnzr333Lk/R4QA5DOv7gAA\n6kH5gaQoP5AU5QeSovxAUpQfSIryA0lRfiApyg8kdVIvdzbfC2KhFvdyl0Aqr+tPOhQH3cq6HZXf\n9hWSbpM0IOmOiLilbP2FWqwLvLKTXQIosSkmWl637af9tgckfUvSxySdK+la2+e2+3gAequT1/wr\nJD0fES9ExCFJ90haVU0sAN3WSflPk/TSjPs7imVvYnut7Unbk4d1sIPdAahS19/tj4jxiBiLiLFB\nLej27gC0qJPy75S0bMb904tlAOaATsq/WdLZts+0PV/SJyVtrCYWgG5r+1JfRByxfb2kH2n6Ut/6\niHimsmQAuqqj6/wRcb+k+yvKAqCH+HgvkBTlB5Ki/EBSlB9IivIDSVF+ICnKDyRF+YGkKD+QFOUH\nkqL8QFKUH0iK8gNJUX4gKcoPJEX5gaQoP5AU5QeSovxAUpQfSIryA0lRfiApyg8kRfmBpCg/kBTl\nB5Ki/EBSlB9IivIDSXU0S6/t7ZL2Szoq6UhEjFURCqjCnz5xQcOxL3/l9tJtv7j6H0vHY/LptjL1\nk47KX/j7iPh9BY8DoId42g8k1Wn5Q9JDtp+wvbaKQAB6o9On/RdFxE7bSyU9bPvZiHh05grFL4W1\nkrRQizrcHYCqdHTmj4idxc8pSfdKWjHLOuMRMRYRY4Na0MnuAFSo7fLbXmz7XW/clnS5pLn/FiiQ\nRCdP+4cl3Wv7jcf5XkQ8WEkqAF3Xdvkj4gVJf1Nhlq46sOptr0jePL5koHR8aP3jVcZBD0yNNX5i\n+8Xt/9DDJP2JS31AUpQfSIryA0lRfiApyg8kRfmBpKr4q7454XeXlP+eW3TWH8sfYH2FYVCNeeWX\nZ+N9BxqOrVz6bOm2E/5IW5HmEs78QFKUH0iK8gNJUX4gKcoPJEX5gaQoP5BUmuv8X7jq/0rHv7zt\n8h4lQVUGzjqjdPzZv2v84YzRn32qdNv3bt7aVqa5hDM/kBTlB5Ki/EBSlB9IivIDSVF+ICnKDySV\n5jr/oI/UHQEVO+mO19re9sBvTq4wydzEmR9IivIDSVF+ICnKDyRF+YGkKD+QFOUHkmp6nd/2eklX\nSZqKiPOKZUOSvi9puaTtklZHxCvdi9ncsYtGS8cvXvhYj5KgV5Yv/kPb2y575GiFSeamVs7835F0\nxVuW3SRpIiLOljRR3AcwhzQtf0Q8KmnvWxavkrShuL1B0tUV5wLQZe2+5h+OiF3F7ZclDVeUB0CP\ndPyGX0SEpGg0bnut7Unbk4d1sNPdAahIu+XfbXtEkoqfU41WjIjxiBiLiLFBLWhzdwCq1m75N0pa\nU9xeI+m+auIA6JWm5bd9t6THJX3A9g7b10m6RdJHbT8n6bLiPoA5pOl1/oi4tsHQyoqzdOTFq95R\nOr50YFGPkqAqJy1/X+n4J4Y2tv3Y7/ht+cdSMnwKgE/4AUlRfiApyg8kRfmBpCg/kBTlB5I6Yb66\n+6T37+9o+9effXdFSVCVl76+uHT8wgXHSsfv3Hd648E/7msn0gmFMz+QFOUHkqL8QFKUH0iK8gNJ\nUX4gKcoPJHXCXOfv1NLJ8mvGmN3AqUtKx3d//JyGY0Ord5Ru+5Nz7myy94Wlo7d/q/H3yi7d/dMm\nj33i48wPJEX5gaQoP5AU5QeSovxAUpQfSIryA0lxnb9wYKj892D5X5Z35tjF55eOx4BLx1+6rPFM\nSIfee7h023nzy7+k+qGLv1E6PlgeTS8fbZztP164pnTbvcfKP3uxaF559uFNjb/joeH8colw5geS\novxAUpQfSIryA0lRfiApyg8kRfmBpJpe57e9XtJVkqYi4rxi2c2SPiNpT7Hauoi4v1shW3Hw9cHS\n8WNNruz+97pbS8c3Xj963JladeOSO0rH56n8YvqBONRw7HdHy6+Ff3PPpaXjlz1yQ+n4u38+v3R8\n5KHdDcf8Yvnf8+/ZVj7t+vBA+WcYYvPW0vHsWjnzf0fSFbMsvzUiRot/tRYfwPFrWv6IeFTS3h5k\nAdBDnbzmv972U7bX2z6lskQAeqLd8t8u6SxJo5J2SfpqoxVtr7U9aXvysA62uTsAVWur/BGxOyKO\nRsQxSd+WtKJk3fGIGIuIsUE1/iMPAL3VVvltj8y4e42kp6uJA6BXWrnUd7ekSyWdanuHpP+UdKnt\nUU3/ZeR2SZ/tYkYAXeCI3v1l88keigu8smf7m+m3X/rb0vFlH9rZoyTHb88DJfPMS1ryTOPr3fMf\n3Fx1nMrsvPEjpeO/+Odvlo7f8+p7Ssfv+sCy4840122KCe2LvU2+ZWEan/ADkqL8QFKUH0iK8gNJ\nUX4gKcoPJJXmq7vP/NfH647QthH9f90RumLRJXuar1Ti33/88dLxc/Szjh7/RMeZH0iK8gNJUX4g\nKcoPJEX5gaQoP5AU5QeSSnOdHyeeM+5jou1OcOYHkqL8QFKUH0iK8gNJUX4gKcoPJEX5gaQoP5AU\n5QeSovxAUpQfSIryA0lRfiApyg8kRfmBpJr+Pb/tZZLukjQsKSSNR8RttockfV/ScknbJa2OiFe6\nFxXZDLj83PTKOYOl43/1QJVpTjytnPmPSPp8RJwr6cOSPmf7XEk3SZqIiLMlTRT3AcwRTcsfEbsi\n4sni9n5J2ySdJmmVpA3FahskXd2tkACqd1yv+W0vl3S+pE2ShiNiVzH0sqZfFgCYI1ouv+13SvqB\npBsiYt/MsYgITb8fMNt2a21P2p48rIMdhQVQnZbKb3tQ08X/bkT8sFi82/ZIMT4iaWq2bSNiPCLG\nImJsUAuqyAygAk3Lb9uS7pS0LSK+NmNoo6Q1xe01ku6rPh6Abmnlq7svlPRpSVttbymWrZN0i6T/\ntX2dpBclre5ORGR1NI6Vr8CnVDrStPwR8ZgkNxheWW0cAL3C704gKcoPJEX5gaQoP5AU5QeSovxA\nUkzRjTnrtQ+9VneEOY0zP5AU5QeSovxAUpQfSIryA0lRfiApyg8kxXV+9K1mX92NznB0gaQoP5AU\n5QeSovxAUpQfSIryA0lRfiAprvOjNgcfeU/p+NHRJt/bj45w5geSovxAUpQfSIryA0lRfiApyg8k\nRfmBpBwR5SvYyyTdJWlYUkgaj4jbbN8s6TOS9hSrrouI+8se62QPxQVmVm+gWzbFhPbFXreybisf\n8jki6fMR8aTtd0l6wvbDxditEfFf7QYFUJ+m5Y+IXZJ2Fbf3294m6bRuBwPQXcf1mt/2cknnS9pU\nLLre9lO219s+pcE2a21P2p48rIMdhQVQnZbLb/udkn4g6YaI2CfpdklnSRrV9DODr862XUSMR8RY\nRIwNakEFkQFUoaXy2x7UdPG/GxE/lKSI2B0RRyPimKRvS1rRvZgAqta0/LYt6U5J2yLiazOWj8xY\n7RpJT1cfD0C3tPJu/4WSPi1pq+0txbJ1kq61Parpy3/bJX22KwkBdEUr7/Y/Jmm264al1/QB9Dc+\n4QckRfmBpCg/kBTlB5Ki/EBSlB9IivIDSVF+ICnKDyRF+YGkKD+QFOUHkqL8QFKUH0iq6Vd3V7oz\ne4+kF2csOlXS73sW4Pj0a7Z+zSWRrV1VZjsjIsrnPi/0tPxv27k9GRFjtQUo0a/Z+jWXRLZ21ZWN\np/1AUpQfSKru8o/XvP8y/ZqtX3NJZGtXLdlqfc0PoD51n/kB1KSW8tu+wvavbD9v+6Y6MjRie7vt\nrba32J6sOct621O2n56xbMj2w7afK37OOk1aTdlutr2zOHZbbF9ZU7Zltn9s+5e2n7H9L8XyWo9d\nSa5ajlvPn/bbHpD0a0kflbRD0mZJ10bEL3sapAHb2yWNRUTt14RtXyLpVUl3RcR5xbKvSNobEbcU\nvzhPiYgb+yTbzZJerXvm5mJCmZGZM0tLulrSP6nGY1eSa7VqOG51nPlXSHo+Il6IiEOS7pG0qoYc\nfS8iHpW09y2LV0naUNzeoOn/PD3XIFtfiIhdEfFkcXu/pDdmlq712JXkqkUd5T9N0ksz7u9Qf035\nHZIesv2E7bV1h5nFcDFtuiS9LGm4zjCzaDpzcy+9ZWbpvjl27cx4XTXe8Hu7iyLig5I+JulzxdPb\nvhTTr9n66XJNSzM398osM0v/RZ3Hrt0Zr6tWR/l3Slo24/7pxbK+EBE7i59Tku5V/80+vPuNSVKL\nn1M15/mLfpq5ebaZpdUHx66fZryuo/ybJZ1t+0zb8yV9UtLGGnK8je3FxRsxsr1Y0uXqv9mHN0pa\nU9xeI+m+GrO8Sb/M3NxoZmnVfOz6bsbriOj5P0lXavod/99I+rc6MjTI9deSflH8e6bubJLu1vTT\nwMOafm/kOklLJE1Iek7SI5KG+ijb/0jaKukpTRdtpKZsF2n6Kf1TkrYU/66s+9iV5KrluPEJPyAp\n3vADkqL8QFKUH0iK8gNJUX4gKcoPJEX5gaQoP5DUnwER0gZdW5joZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADN1JREFUeJzt3X+s3XV9x/Hn23pbtLiF4qwN1OEI\naBjJirtDnQx1iEHCLPyxSs1MtxCrmWxjcckI+0P+cFmjE0fioimjUjZFFwHhD5xiM0cMjHFhHVC6\nyY8VaVMoBDbBhXKh7/1xv5AL3PM9l/P79v18JDfnnO/7+z3fd77pq9/vOZ9zzicyE0n1vG7cDUga\nD8MvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqmo149yZ8tjRR7BylHuUirlWX7Oc3kwFrNuX+GP\niLOAy4FlwN9l5pa29Y9gJe+OM/rZpaQWt+eORa/b82V/RCwD/hb4CHASsDEiTur1+SSNVj+v+U8F\nHsjMhzLzOeBbwPrBtCVp2PoJ/zHAI/Me722WvUxEbI6ImYiYmeVgH7uTNEhDf7c/M7dm5nRmTk+x\nYti7k7RI/YR/H7B23uNjm2WSloB+wn8HcEJEvD0ilgPnAzcOpi1Jw9bzUF9mPh8RFwLfZ26ob1tm\n7hpYZ5KGqq9x/sy8CbhpQL1IGiE/3isVZfilogy/VJThl4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWi\nDL9UlOGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8VZfil\nogy/VJThl4oy/FJRfc3SGxF7gKeBF4DnM3N6EE1JAA9+8b2t9d0f/0prfSqWdayd/oebW7d9w3f/\nrbV+OOgr/I0PZuYTA3geSSPkZb9UVL/hT+AHEXFnRLRfR0maKP1e9p+Wmfsi4i3AzRHxn5l5y/wV\nmv8UNgMcwRv73J2kQenrzJ+Z+5rbA8D1wKkLrLM1M6czc3qKFf3sTtIA9Rz+iFgZEW968T7wYeDe\nQTUmabj6uexfDVwfES8+zzcz858G0pWkoes5/Jn5EPBrA+xFxTz6p7/ZWv/Rx77QWp/N5b3vPHvf\n9HDhUJ9UlOGXijL8UlGGXyrK8EtFGX6pqEF8q0/qyTNrD7XWV72uj6E8deWZXyrK8EtFGX6pKMMv\nFWX4paIMv1SU4ZeKcpxfQ/XM7767Y+3a8y7vsnW0Vr/2P+9srf9wQ+dfkl/58K7Wbds/gXB48Mwv\nFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0U5zq++PHvOqyZpepnP/dW2jrUTp9rH8bvZfsVZrfW33ndr\nX89/uPPMLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFdR3nj4htwDnAgcw8uVm2Cvg2cBywB9iQmU8N\nr01Nqv2/92xr/YNvaKsva912054Ptdbfernj+P1YzJn/KuCVn6a4GNiRmScAO5rHkpaQruHPzFuA\nJ1+xeD2wvbm/HTh3wH1JGrJeX/Ovzsz9zf1HgdUD6kfSiPT9hl9mJpCd6hGxOSJmImJmloP97k7S\ngPQa/sciYg1Ac3ug04qZuTUzpzNzeooVPe5O0qD1Gv4bgU3N/U3ADYNpR9KodA1/RFwD3Aa8IyL2\nRsQFwBbgzIi4H/hQ81jSEtJ1nD8zN3YonTHgXjSBXn/sMa31Xb/19db6bL7QsbZ7tn3fP73sxNb6\nSm5vfwK18hN+UlGGXyrK8EtFGX6pKMMvFWX4paL86e7ilv3qO1rr09+8d2j7/th1f9xaP/7afx3a\nvuWZXyrL8EtFGX6pKMMvFWX4paIMv1SU4ZeKcpy/uIc/enRr/TtH/3uXZ2j/+e2PP/g7HWsnbnmw\nddvOXwbWIHjml4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWiHOc/zD35B+9trV//6S92eYap1uqnH3l/\na312U+dZml54/Kdd9q1h8swvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0V1HeePiG3AOcCBzDy5WXYp\n8Eng8Wa1SzLzpmE1qXZtv71/6+e/0mXrI/ra9217j2utr90zvN/9V38Wc+a/CjhrgeVfzsx1zZ/B\nl5aYruHPzFuAJ0fQi6QR6uc1/4URcXdEbIuIowbWkaSR6DX8XwWOB9YB+4EvdVoxIjZHxExEzMxy\nsMfdSRq0nsKfmY9l5guZeQi4Aji1Zd2tmTmdmdNTdP6Sh6TR6in8EbFm3sPzAN/SlZaYxQz1XQN8\nAHhzROwFPgd8ICLWAQnsAT41xB4lDUHX8GfmxgUWXzmEXtSjn1zyxo612Rzur9+/bUt7PYe6d/XD\nT/hJRRl+qSjDLxVl+KWiDL9UlOGXivKnu5eAQ+8/pbX++envDm3fZ957fmv9yBk/37VUeeaXijL8\nUlGGXyrK8EtFGX6pKMMvFWX4paIc518C/vKqra31k6d6/+Lsn+0/vbX+ixufaq0P9wvDGibP/FJR\nhl8qyvBLRRl+qSjDLxVl+KWiDL9UlOP8S8Apy9v/j+7n57lv+/q7WutveerWnp9bk80zv1SU4ZeK\nMvxSUYZfKsrwS0UZfqkowy8V1XWcPyLWAlcDq5mbcXlrZl4eEauAbwPHAXuADZnZ/uVvLeiR75zc\nWp+KnUPb95ofPdFa9/v6h6/FnPmfBz6bmScB7wE+ExEnARcDOzLzBGBH81jSEtE1/Jm5PzPvau4/\nDewGjgHWA9ub1bYD5w6rSUmD95pe80fEccApwO3A6szc35QeZe5lgaQlYtHhj4gjgWuBizLzZ/Nr\nmZnMvR+w0HabI2ImImZmOdhXs5IGZ1Hhj4gp5oL/jcy8rln8WESsaeprgAMLbZuZWzNzOjOnp1gx\niJ4lDUDX8EdEAFcCuzPzsnmlG4FNzf1NwA2Db0/SsCzmK73vAz4B3BPx0pjTJcAW4B8j4gLgYWDD\ncFpc+rpNsf036/6htd7tK7v/e+jZjrXf+N5Frdu+8+H7Wus6fHUNf2b+GIgO5TMG246kUfETflJR\nhl8qyvBLRRl+qSjDLxVl+KWi/OnuEXh21fLW+mlH/LzLMyxrrX7//97WsXbi5jtatz3UZc86fHnm\nl4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paL8Pv8I\n/MLOR1vrf7T3t1vrX1v7L4NsRwI880tlGX6pKMMvFWX4paIMv1SU4ZeKMvxSUV3H+SNiLXA1sBpI\nYGtmXh4RlwKfBB5vVr0kM28aVqNL2fP//XBrfe972rc/h18fYDfSnMV8yOd54LOZeVdEvAm4MyJu\nbmpfzsy/Hl57koala/gzcz+wv7n/dETsBo4ZdmOShus1veaPiOOAU4Dbm0UXRsTdEbEtIo7qsM3m\niJiJiJlZDvbVrKTBWXT4I+JI4Frgosz8GfBV4HhgHXNXBl9aaLvM3JqZ05k5PcWKAbQsaRAWFf6I\nmGIu+N/IzOsAMvOxzHwhMw8BVwCnDq9NSYPWNfwREcCVwO7MvGze8jXzVjsPuHfw7UkalsW82/8+\n4BPAPRGxs1l2CbAxItYxN/y3B/jUUDqUNBSLebf/x0AsUHJMX1rC/ISfVJThl4oy/FJRhl8qyvBL\nRRl+qSjDLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFGX6pqMjM0e0s4nFg/u9Yvxl4YmQNvDaT2tuk\n9gX21qtB9vbLmflLi1lxpOF/1c4jZjJzemwNtJjU3ia1L7C3Xo2rNy/7paIMv1TUuMO/dcz7bzOp\nvU1qX2BvvRpLb2N9zS9pfMZ95pc0JmMJf0ScFRH/FREPRMTF4+ihk4jYExH3RMTOiJgZcy/bIuJA\nRNw7b9mqiLg5Iu5vbhecJm1MvV0aEfuaY7czIs4eU29rI+KfI+K+iNgVEX/SLB/rsWvpayzHbeSX\n/RGxDPgJcCawF7gD2JiZ9420kQ4iYg8wnZljHxOOiNOBZ4CrM/PkZtkXgCczc0vzH+dRmfnnE9Lb\npcAz4565uZlQZs38maWBc4HfZ4zHrqWvDYzhuI3jzH8q8EBmPpSZzwHfAtaPoY+Jl5m3AE++YvF6\nYHtzfztz/3hGrkNvEyEz92fmXc39p4EXZ5Ye67Fr6WssxhH+Y4BH5j3ey2RN+Z3ADyLizojYPO5m\nFrC6mTYd4FFg9TibWUDXmZtH6RUzS0/MsetlxutB8w2/VzstM98FfAT4THN5O5Fy7jXbJA3XLGrm\n5lFZYGbpl4zz2PU64/WgjSP8+4C18x4f2yybCJm5r7k9AFzP5M0+/NiLk6Q2twfG3M9LJmnm5oVm\nlmYCjt0kzXg9jvDfAZwQEW+PiOXA+cCNY+jjVSJiZfNGDBGxEvgwkzf78I3Apub+JuCGMfbyMpMy\nc3OnmaUZ87GbuBmvM3Pkf8DZzL3j/yDwF+PooUNfvwL8R/O3a9y9Adcwdxk4y9x7IxcARwM7gPuB\nHwKrJqi3vwfuAe5mLmhrxtTbacxd0t8N7Gz+zh73sWvpayzHzU/4SUX5hp9UlOGXijL8UlGGXyrK\n8EtFGX6pKMMvFWX4paL+H5OL6YVERhITAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADktJREFUeJzt3X+MXHW5x/HP07JbSqn3duF201vK\nD6WAFbnFO2lRiD9CQSRKwR+EemNqUl0hVC83knuxxtg/MCEKkkr8wYJN2xsuoCmERlGEaiDeaGUh\npYAVWsnWti79QYWWy2273T7+sadmhT3fmc6cmTPd5/1KNjtznnPmPJnsZ8/MfM+cr7m7AMQzruwG\nAJSD8ANBEX4gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCOq4Vu6s0yb48ZrUyl0CoezX/+mgH7Ba1m0o\n/GZ2maRlksZLutvdb0mtf7wmaa5d3MguASSs87U1r1v3y34zGy/pu5I+ImmWpAVmNqvexwPQWo28\n558jabO7v+TuByXdJ2l+MW0BaLZGwj9d0tYR97dly/6OmfWYWZ+Z9Q3qQAO7A1Ckpn/a7+697l5x\n90qHJjR7dwBq1Ej4t0uaMeL+KdkyAMeARsL/pKSZZnaGmXVKukbSmmLaAtBsdQ/1ufshM1ss6REN\nD/Utd/fnC+sMQFM1NM7v7g9LerigXgC0EKf3AkERfiAowg8ERfiBoAg/EBThB4Ii/EBQhB8IivAD\nQRF+ICjCDwRF+IGgCD8QFOEHgiL8QFCEHwiK8ANBEX4gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIrw\nA0ERfiAowg8ERfiBoAg/EFRDs/SaWb+kfZKGJB1y90oRTQFovobCn/mQu+8u4HEAtBAv+4GgGg2/\nS/qFmT1lZj1FNASgNRp92X+Ru283s6mSHjWzP7j7EyNXyP4p9EjS8Tqhwd0BKEpDR35335793inp\nQUlzRlmn190r7l7p0IRGdgegQHWH38wmmdnkI7clXSrpuaIaA9Bcjbzs75b0oJkdeZz/cfefF9IV\ngKarO/zu/pKkfymwFwAtxFAfEBThB4Ii/EBQhB8IivADQRF+IKgivtWHNnbww+lvWW/5t8PJ+nXv\neTxZv2HKi0fd0xHvvvuLyfoJA56sv/q+A8n6affkH9s6H+lLbhsBR34gKMIPBEX4gaAIPxAU4QeC\nIvxAUIQfCIpx/jFg17Xvza3d8Z/fTW5bmTCUrI+rcnxY2D8vWT//H/6UW3vmc8uS21ZTrbf3dS3I\nrXU90tCuxwSO/EBQhB8IivADQRF+ICjCDwRF+IGgCD8QFOP8bcA6OpP1/fPSV0hf/ZVv5db++bj0\nLEmLtlySrG+59exkfdJP1yfrvzrh1Nza4w+eldx29cw1yXo1e9eflFvrauiRxwaO/EBQhB8IivAD\nQRF+ICjCDwRF+IGgCD8QVNVxfjNbLumjkna6+7nZsi5J90s6XVK/pKvd/S/Na3NsG1icvrb+726s\n9r33/LH8T23+WHLLQ58YTNZP2L0uWU9fWV/6c8+/5tbWzWzs+/w/e2Nysn7mnVtza4ca2vPYUMuR\nf4Wky9607CZJa919pqS12X0Ax5Cq4Xf3JyTtedPi+ZJWZrdXSrqy4L4ANFm97/m73X0gu/2ypO6C\n+gHQIg1/4OfursRbPzPrMbM+M+sbVHpuNQCtU2/4d5jZNEnKfu/MW9Hde9294u6VjsQHUwBaq97w\nr5G0MLu9UNJDxbQDoFWqht/M7pX0G0lnm9k2M1sk6RZJl5jZJknzsvsAjiFVx/ndPe/i5xcX3MuY\ntemOucn6Cx+/I1k/XOXx3/notbm1c27sT247tPuVKo/emGuva96Lwpu/sTBZn7L1N03b91jAGX5A\nUIQfCIrwA0ERfiAowg8ERfiBoLh0dwH+eNsFyfoLH09Pk/3a4f3J+qf+8Olk/ewvvphbG9q3L7lt\nNeMmTUrWX/nkecn6/BPzLys+ThOT257z4+uT9TNXMJTXCI78QFCEHwiK8ANBEX4gKMIPBEX4gaAI\nPxAU4/w1Gt89Nbe28qrvJbc9XOVLudXG8Tsv2VLl8es3bvasZP3c5RuT9Zu7v1NlD/lXb7pw/TXJ\nLc9emt73UJU9I40jPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ExTh/jez4/PHqyoTGRpwnfqkzve/T\nZiTrm649Jbd26bynk9v+x9TeZP3U49Lfua92jsGQ50/ibfefnN721U1VHh2N4MgPBEX4gaAIPxAU\n4QeCIvxAUIQfCIrwA0FVHec3s+WSPippp7ufmy1bKunzknZlqy1x94eb1WQ78P0HcmvrDnQkt507\nYTBZf+ix+5L1atcDaMRj/58ea980mD9OL0kfmvh6st53MP8chn9cxXX3y1TLkX+FpMtGWX67u8/O\nfsZ08IGxqGr43f0JSXta0AuAFmrkPf9iM9tgZsvNbEphHQFoiXrD/31J75A0W9KApNvyVjSzHjPr\nM7O+QeW/bwbQWnWF3913uPuQux+WdJekOYl1e9294u6VjsTFHAG0Vl3hN7NpI+5eJem5YtoB0Cq1\nDPXdK+mDkk42s22Svi7pg2Y2W5JL6pf0hSb2CKAJzBPfty7a26zL59rFLdtfqxz8cCVZv/UH6ev6\nn9c5PllftXd6sn7z41fk1s5asT+57XE7XkvWp96bHuj5wYxfJuvn/Py63NpZi/qS2+LorfO12ut7\nrJZ1OcMPCIrwA0ERfiAowg8ERfiBoAg/EBSX7i5A5yPpIaslZ+SeAFmIs/S7urfdNz/d209PfShZ\nH/T08WNif/qy5CgPR34gKMIPBEX4gaAIPxAU4QeCIvxAUIQfCIpx/uAOTUz//x/09PTj1S4rfsaK\nP+XvO7klmo0jPxAU4QeCIvxAUIQfCIrwA0ERfiAowg8ExTh/cJPv+216hdyJ2HCs48gPBEX4gaAI\nPxAU4QeCIvxAUIQfCIrwA0FVHec3sxmSVknqluSSet19mZl1Sbpf0umS+iVd7e5/aV6raIZ911xQ\nZY2nWtIHWq+WI/8hSV9291mSLpB0vZnNknSTpLXuPlPS2uw+gGNE1fC7+4C7P53d3idpo6TpkuZL\nWpmttlLSlc1qEkDxjuo9v5mdLul8Seskdbv7QFZ6WcNvCwAcI2oOv5mdKGm1pBvcfe/Imru7hj8P\nGG27HjPrM7O+QR1oqFkAxakp/GbWoeHg3+PuD2SLd5jZtKw+TdLO0bZ19153r7h7pUMTiugZQAGq\nht/MTNIPJW1092+PKK2RtDC7vVBSejpXAG2llq/0XijpM5KeNbP12bIlkm6R9CMzWyRpi6Srm9Mi\nmum1t3OqR1RVw+/uv5ZkOeWLi20HQKvwbx8IivADQRF+ICjCDwRF+IGgCD8QFJfuDm76428k6x2L\nxyfrg6Oe1I1jAUd+ICjCDwRF+IGgCD8QFOEHgiL8QFCEHwiKcf7g7H/XJ+sr9k5N1hdM3p6sv/Gu\nabm1zq3bktuiuTjyA0ERfiAowg8ERfiBoAg/EBThB4Ii/EBQjPMj6fY7P5msL7hxWbI+7Wubc2uv\nvHpeeue/3ZCuoyEc+YGgCD8QFOEHgiL8QFCEHwiK8ANBEX4gKHNPX3jdzGZIWiWpW5JL6nX3ZWa2\nVNLnJe3KVl3i7g+nHutt1uVzjVm9jyXjTz4pWe9cnT5V5P4zf5Jb+8AzC5Lbdn16V7I+9OpryXpE\n63yt9voeq2XdWk7yOSTpy+7+tJlNlvSUmT2a1W5391vrbRRAeaqG390HJA1kt/eZ2UZJ05vdGIDm\nOqr3/GZ2uqTzJa3LFi02sw1mttzMpuRs02NmfWbWN6gDDTULoDg1h9/MTpS0WtIN7r5X0vclvUPS\nbA2/MrhttO3cvdfdK+5e6dCEAloGUISawm9mHRoO/j3u/oAkufsOdx9y98OS7pI0p3ltAiha1fCb\nmUn6oaSN7v7tEctHXpb1KknPFd8egGap5dP+CyV9RtKzZnbkOs9LJC0ws9kaHv7rl/SFpnSIUg3t\nfiVZP/iJ9FDgO2/L/7PYOO/O5LZXnLMoWecrv42p5dP+X0sabdwwOaYPoL1xhh8QFOEHgiL8QFCE\nHwiK8ANBEX4gqKpf6S0SX+kFmutovtLLkR8IivADQRF+ICjCDwRF+IGgCD8QFOEHgmrpOL+Z7ZK0\nZcSikyXtblkDR6dde2vXviR6q1eRvZ3m7v9Uy4otDf9bdm7W5+6V0hpIaNfe2rUvid7qVVZvvOwH\ngiL8QFBlh7+35P2ntGtv7dqXRG/1KqW3Ut/zAyhP2Ud+ACUpJfxmdpmZvWBmm83spjJ6yGNm/Wb2\nrJmtN7O+kntZbmY7zey5Ecu6zOxRM9uU/R51mrSSeltqZtuz5269mV1eUm8zzOxXZvZ7M3vezP49\nW17qc5foq5TnreUv+81svKQXJV0iaZukJyUtcPfft7SRHGbWL6ni7qWPCZvZ+yW9LmmVu5+bLfum\npD3ufkv2j3OKu/9Xm/S2VNLrZc/cnE0oM23kzNKSrpT0WZX43CX6ulolPG9lHPnnSNrs7i+5+0FJ\n90maX0Ifbc/dn5C0502L50tamd1eqeE/npbL6a0tuPuAuz+d3d4n6cjM0qU+d4m+SlFG+KdL2jri\n/ja115TfLukXZvaUmfWU3cwourNp0yXpZUndZTYziqozN7fSm2aWbpvnrp4Zr4vGB35vdZG7v0fS\nRyRdn728bUs+/J6tnYZrapq5uVVGmVn6b8p87uqd8bpoZYR/u6QZI+6fki1rC+6+Pfu9U9KDar/Z\nh3ccmSQ1+72z5H7+pp1mbh5tZmm1wXPXTjNelxH+JyXNNLMzzKxT0jWS1pTQx1uY2aTsgxiZ2SRJ\nl6r9Zh9eI2lhdnuhpIdK7OXvtMvMzXkzS6vk567tZrx295b/SLpcw5/4/1HSV8voIaevt0t6Jvt5\nvuzeJN2r4ZeBgxr+bGSRpJMkrZW0SdJjkrraqLf/lvSspA0aDtq0knq7SMMv6TdIWp/9XF72c5fo\nq5TnjTP8gKD4wA8IivADQRF+ICjCDwRF+IGgCD8QFOEHgiL8QFB/BcMMVHsmbz+8AAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5 0 4 1 9]\n",
      "<class 'numpy.ndarray'>\n",
      "<class 'numpy.ndarray'>\n",
      "(60000, 28, 28)\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'STOP' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-526fe55edef4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_train\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 20\u001b[0;31m \u001b[0mSTOP\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     21\u001b[0m \u001b[0;31m# 5. Preprocess input data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     22\u001b[0m \u001b[0mX_train\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mX_train\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m28\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m28\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'STOP' is not defined"
     ]
    }
   ],
   "source": [
    "# 3. Import libraries and modules\n",
    "import numpy as np\n",
    "np.random.seed(123)  # for reproducibility\n",
    " \n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, Dropout, Activation, Flatten\n",
    "from keras.layers import Convolution2D, MaxPooling2D\n",
    "from keras.utils import np_utils\n",
    "from keras.datasets import mnist\n",
    " \n",
    "# 4. Load pre-shuffled MNIST data into train and test sets\n",
    "(X_train, y_train), (X_test, y_test) = mnist.load_data()\n",
    "for j in range(5):\n",
    "    plt.imshow(X_train[j])\n",
    "    plt.show()\n",
    "print(y_train[0:5])\n",
    "print(type(X_train))\n",
    "print(type(y_train))\n",
    "print(np.shape(X_train))\n",
    "STOP\n",
    "# 5. Preprocess input data\n",
    "X_train = X_train.reshape(X_train.shape[0], 28, 28, 1)\n",
    "X_test = X_test.reshape(X_test.shape[0],  28, 28, 1)\n",
    "X_train = X_train.astype('float32')\n",
    "X_test = X_test.astype('float32')\n",
    "X_train /= 255\n",
    "X_test /= 255\n",
    " \n",
    "# 6. Preprocess class labels\n",
    "print(np.shape(y_train))\n",
    "Y_train = np_utils.to_categorical(y_train, 10)\n",
    "Y_test = np_utils.to_categorical(y_test, 10)\n",
    " \n",
    "print(np.shape(Y_train))\n",
    "    \n",
    "STOP\n",
    "# 7. Define model architecture\n",
    "model = Sequential()\n",
    " \n",
    "model.add(Convolution2D(32, 3, 3, activation='relu', input_shape=(28,28,1)))\n",
    "model.add(Convolution2D(32, 3, 3, activation='relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2,2)))\n",
    "model.add(Dropout(0.25))\n",
    " \n",
    "model.add(Flatten())\n",
    "model.add(Dense(128, activation='relu'))\n",
    "model.add(Dropout(0.5))\n",
    "model.add(Dense(10, activation='softmax'))\n",
    " \n",
    "# 8. Compile model\n",
    "model.compile(loss='categorical_crossentropy',\n",
    "              optimizer='adam',\n",
    "              metrics=['accuracy'])\n",
    " \n",
    "# 9. Fit model on training data\n",
    "model.fit(X_train, Y_train, \n",
    "          batch_size=32, nb_epoch=10, verbose=1)\n",
    " \n",
    "# 10. Evaluate model on test data\n",
    "score = model.evaluate(X_test, Y_test, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}