Newer
Older
R_phipi / DNN.ipynb
@Davide Lancierini Davide Lancierini on 8 Oct 2018 25 KB Set up NN for sgn/bkg classification
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import pickle\n",
    "import math\n",
    "\n",
    "trunc_normal= tf.truncated_normal_initializer(stddev=1)\n",
    "normal = tf.random_normal_initializer(stddev=1)\n",
    "\n",
    "from architectures.utils.toolbox import *\n",
    "from architectures.DNN import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/disk/lhcb_data/davide/Rphipi/NN_test/MC_for_NN.pickle', 'rb') as f:\n",
    "    MC_sig_dict=pickle.load(f, encoding='latin1')\n",
    "with open('/disk/lhcb_data/davide/Rphipi/NN_test/data_for_NN.pickle', 'rb') as f:\n",
    "    data_bkg_dict=pickle.load(f, encoding='latin1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ds_mass= 1968"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m=MC_sig_dict[\"Ds_ConsD_M\"].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "MC_sig_dict[\"Ds_OWNPV_FIT\"]=MC_sig_dict[\"Ds_OWNPV_CHI2\"]/MC_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
    "MC_sig_dict[\"Ds_ENDV_FIT\"]=MC_sig_dict[\"Ds_ENDVERTEX_CHI2\"]/MC_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "MC_sig_dict[\"Ds_ConsD_M_norm\"]=MC_sig_dict[\"Ds_ConsD_M\"]/Ds_mass\n",
    "#MC_sig_dict[\"Ds_Hlt2Phys_TOS_int\"]=MC_sig_dict[\"Ds_Hlt2Phys_TOS\"].astype(np.float32)\n",
    "#MC_sig_dict[\"Ds_Hlt1TrackMVADecision_TOS_int\"]=MC_sig_dict[\"Ds_Hlt1TrackMVADecision_TOS\"].astype(np.float32)\n",
    "#MC_sig_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS_int\"]=MC_sig_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\"].astype(np.float32)\n",
    "\n",
    "del MC_sig_dict[\"Ds_ConsD_M\"], MC_sig_dict[\"Ds_ENDVERTEX_CHI2\"], MC_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "del MC_sig_dict[\"Ds_OWNPV_CHI2\"], MC_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
    "del MC_sig_dict[\"Ds_Hlt2Phys_TOS\"], MC_sig_dict[\"Ds_Hlt1TrackMVADecision_TOS\"], MC_sig_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\"]\n",
    "del MC_sig_dict[\"Ds_IPCHI2_OWNPV\"] #temporary\n",
    "del MC_sig_dict[\"mu_plus_MC15TuneV1_ProbNNmu\"] #temporary\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "data_bkg_dict[\"Ds_OWNPV_FIT\"]=data_bkg_dict[\"Ds_OWNPV_CHI2\"]/data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
    "data_bkg_dict[\"Ds_ENDV_FIT\"]=data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"]/data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "data_bkg_dict[\"Ds_ConsD_M_norm\"]=data_bkg_dict[\"Ds_ConsD_M\"]/Ds_mass\n",
    "#data_bkg_dict[\"Ds_Hlt2Phys_TOS_int\"]=data_bkg_dict[\"Ds_Hlt2Phys_TOS\"].astype(np.float32)\n",
    "#data_bkg_dict[\"Ds_Hlt1TrackMVADecision_TOS_int\"]=data_bkg_dict[\"Ds_Hlt1TrackMVADecision_TOS\"].astype(np.float32)\n",
    "#data_bkg_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS_int\"]=data_bkg_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\"].astype(np.float32)\n",
    "\n",
    "del data_bkg_dict[\"Ds_ConsD_M\"], data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"], data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "del data_bkg_dict[\"Ds_OWNPV_CHI2\"], data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
    "del data_bkg_dict[\"Ds_Hlt2Phys_TOS\"], data_bkg_dict[\"Ds_Hlt1TrackMVADecision_TOS\"], data_bkg_dict[\"Ds_Hlt2RareCharmD2PiMuMuOSDecision_TOS\"]\n",
    "del data_bkg_dict[\"Ds_IPCHI2_OWNPV\"] #temporary\n",
    "del data_bkg_dict[\"mu_plus_MC15TuneV1_ProbNNmu\"] #temporary\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim=len(MC_sig_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ds_DIRA_OWNPV\n",
      "Ds_IP_OWNPV\n",
      "Ds_OWNPV_FIT\n",
      "Ds_ENDV_FIT\n",
      "Ds_ConsD_M_norm\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(23875, 5)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MC_sig = np.array(\n",
    "    \n",
    "        [np.zeros(shape=(dim), dtype=np.float32) for i in range(m)]\n",
    "    \n",
    "        )\n",
    "\n",
    "for event in range(m):\n",
    "    for i, key in enumerate(MC_sig_dict):\n",
    "        MC_sig[event][i]=MC_sig_dict[key][event]\n",
    "        \n",
    "for key in MC_sig_dict:\n",
    "    print(key)\n",
    "\n",
    "MC_sig.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ds_DIRA_OWNPV\n",
      "Ds_IP_OWNPV\n",
      "Ds_OWNPV_FIT\n",
      "Ds_ENDV_FIT\n",
      "Ds_ConsD_M_norm\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(23875, 5)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_bkg = np.array(\n",
    "    \n",
    "        [np.zeros(shape=(dim), dtype=np.float32) for i in range(m)]\n",
    "    \n",
    "        )\n",
    "\n",
    "for event in range(m):\n",
    "    for i, key in enumerate(data_bkg_dict):\n",
    "        data_bkg[event][i]=data_bkg_dict[key][event]\n",
    "        \n",
    "for key in data_bkg_dict:\n",
    "    print(key)\n",
    "data_bkg.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "val_size=1000\n",
    "test_size=3000\n",
    "\n",
    "train_size=m-val_size-test_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(19875, 1000, 3000)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(train_size, val_size, test_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "MC_sig_label=np.concatenate((MC_sig,np.ones(shape=(MC_sig.shape[0],1))),axis=1)\n",
    "data_bkg_label=np.concatenate((data_bkg,np.zeros(shape=(data_bkg.shape[0],1))),axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99996907, 0.03041327, 0.38916221, 0.28149822, 1.00150383,\n",
       "       1.        ])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "data=np.concatenate((MC_sig_label,data_bkg_label), axis =0)\n",
    "data[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "np.random.shuffle(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "X=data[:,0:dim]\n",
    "Y_labels=data[:,dim].astype(int)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_labels=Y_labels.reshape(Y_labels.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp=np.zeros(shape=(m,2))\n",
    "for i in range(m):\n",
    "    if Y_labels[i]==0:\n",
    "        temp[i][0]=1\n",
    "    else:\n",
    "        temp[i][1]=1\n",
    "Y_labels=temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X[0:train_size]\n",
    "Y_train = Y_labels[0:train_size]\n",
    "\n",
    "X_val = X[train_size:train_size+val_size]\n",
    "Y_val = Y_labels[train_size:train_size+val_size]\n",
    "\n",
    "X_test = X[train_size+val_size:train_size+val_size+test_size]\n",
    "Y_test = Y_labels[train_size+val_size:train_size+val_size+test_size]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some constants\n",
    "\n",
    "LEARNING_RATE = 0.0001\n",
    "BETA1 = 0.5\n",
    "BATCH_SIZE = 64\n",
    "EPOCHS = 20000\n",
    "SAVE_SAMPLE_PERIOD = 2000\n",
    "\n",
    "#task='TRAIN'\n",
    "task='TEST'\n",
    "PATH='test_1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "global sizes\n",
    "sizes = {\n",
    "    'dense_layers': [\n",
    "                        #(512, 'bn', 0.8, lrelu, tf.glorot_uniform_initializer()),\n",
    "                        #(128, 'bn', 0.5, tf.nn.relu, tf.glorot_uniform_initializer()),\n",
    "                        (16, 'bn', 0.5, tf.nn.relu, tf.glorot_uniform_initializer()),\n",
    "                        (8, 'bn', 0.5, tf.nn.relu, tf.glorot_uniform_initializer()),\n",
    "                        (4, 'bn', 0.5, tf.nn.relu, tf.glorot_uniform_initializer()),\n",
    "                    ],\n",
    "    'n_classes':2,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bkg():\n",
    "    \n",
    "    tf.reset_default_graph()\n",
    "    nn = DNN(dim, sizes,\n",
    "              lr=LEARNING_RATE, beta1=BETA1,\n",
    "              batch_size=BATCH_SIZE, epochs=EPOCHS,\n",
    "              save_sample=SAVE_SAMPLE_PERIOD, path=PATH)\n",
    "    \n",
    "    vars_to_train= tf.trainable_variables()\n",
    "    \n",
    "    if task == 'TRAIN':\n",
    "        init_op = tf.global_variables_initializer()\n",
    "        \n",
    "    if task == 'TEST':\n",
    "        vars_all = tf.global_variables()\n",
    "        vars_to_init = list(set(vars_all)-set(vars_to_train))\n",
    "        init_op = tf.variables_initializer(vars_to_init)\n",
    "    \n",
    "    # Add ops to save and restore all the variables.\n",
    "    saver = tf.train.Saver()\n",
    "    gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=0.333)\n",
    "    \n",
    "    with tf.Session(config=tf.ConfigProto(gpu_options=gpu_options)) as sess:\n",
    "        \n",
    "        sess.run(init_op)\n",
    "\n",
    "        if task=='TRAIN':\n",
    "            print('\\n Training...')\n",
    "            \n",
    "            if os.path.exists(PATH+'/CNN_model.ckpt.index'):\n",
    "                saver.restore(sess,PATH+'/CNN_model.ckpt')\n",
    "                print('Model restored.')\n",
    "            \n",
    "            nn.set_session(sess)\n",
    "            nn.fit(X_train, Y_train, X_val, Y_val)\n",
    "            \n",
    "            save_path = saver.save(sess, PATH+'/CNN_model.ckpt')\n",
    "            print(\"Model saved in path: %s\" % save_path)\n",
    "        \n",
    "        if task=='TEST':\n",
    "            print('\\n Evaluate model on test set...')\n",
    "            saver.restore(sess,PATH+'/CNN_model.ckpt')\n",
    "            print('Model restored.')\n",
    "            \n",
    "            nn.set_session(sess)\n",
    "            #nn.test(X_test, Y_test)\n",
    "            \n",
    "            output = nn.predict(X_test)\n",
    "            \n",
    "            return output\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Propagation\n",
      "Input for propagation (?, 5)\n",
      "(?, 5)\n",
      "(?, 16)\n",
      "(?, 8)\n",
      "(?, 4)\n",
      "Logits shape (?, 2)\n",
      "Propagation\n",
      "Input for propagation (?, 5)\n",
      "(?, 5)\n",
      "(?, 16)\n",
      "(?, 8)\n",
      "(?, 4)\n",
      "Logits shape (?, 2)\n",
      "\n",
      " Evaluate model on test set...\n",
      "INFO:tensorflow:Restoring parameters from test_1/CNN_model.ckpt\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:tensorflow:Restoring parameters from test_1/CNN_model.ckpt\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model restored.\n"
     ]
    }
   ],
   "source": [
    "if __name__=='__main__':\n",
    "\n",
    "    if task == 'TRAIN':\n",
    "        if not os.path.exists(PATH):\n",
    "            os.mkdir(PATH)\n",
    "    \n",
    "        elif os.path.exists(PATH):\n",
    "            if os.path.exists(PATH+'/checkpoint'):\n",
    "                ans = input('A previous checkpoint already exists, choose the action to perform \\n \\n 1) Overwrite the current model saved at '+PATH+'/checkpoint \\n 2) Start training a new model \\n 3) Restore and continue training the previous model \\n ')\n",
    "                \n",
    "                if ans == '1':\n",
    "                    print('Overwriting existing model in '+PATH)\n",
    "                    for file in os.listdir(PATH):\n",
    "                        file_path = os.path.join(PATH, file)\n",
    "                        try:\n",
    "                            if os.path.isfile(file_path):\n",
    "                                os.unlink(file_path)\n",
    "                            #elif os.path.isdir(file_path): shutil.rmtree(file_path)\n",
    "                        except Exception as e:\n",
    "                            print(e)\n",
    "                            \n",
    "                elif ans == '2':\n",
    "                    PATH = input('Specify the name of the model, a new directory will be created.\\n')\n",
    "                    os.mkdir(PATH)         \n",
    "        bkg()\n",
    "\n",
    "    elif task == 'TEST': \n",
    "        if not os.path.exists(PATH+'/checkpoint'):\n",
    "            print('No checkpoint to test')\n",
    "        else:\n",
    "            output = bkg()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3000, 2)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=np.argmax(output, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "NN_selected=X_test[a.astype(np.bool)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=[MC_sig_dict[\"Ds_ConsD_M_norm\"][i][0]*Ds_mass for i in range(m)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "b=[NN_selected[i][4]*Ds_mass for i in range(NN_selected.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAD8CAYAAAB3qPkTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHSJJREFUeJzt3X/wZXV93/HnS0CSaX4A8tXSXciSZJOKnRHNtyupTUPVAEKbJa1kcDrKKJ1NWuxozLRC0on5RWfTRElto+0aqJhJgjRq2SrWbFHqOCM/FkUEkbAiDRt2YJNF1ElDCr77x/185bp7v/u9d7/3e++59z4fM3fuOZ/7Off7Pp8533PO+5zP+dxUFZIkSZIkddVzph2AJEmSJElHY+IqSZIkSeo0E1dJkiRJUqeZuEqSJEmSOs3EVZIkSZLUaSaukiRJkqROM3GVJEmSJHWaiaskSZIkqdNMXCVJkiRJnXb8tAM4mlNPPbW2bNky7TAkSXPirrvu+vOqWpp2HLPMY7MkaZyGPTZ3OnHdsmULe/funXYYkqQ5keT/TDuGWeexWZI0TsMem+0qLEmSJEnqNBNXSZIkSVKnmbhKkiRJkjrNxFWSJEmS1GkmrpIkSZKkTjNxlSRJkiR1momrJEmSJKnTTFwlSZIkSZ1m4ipJkiRJ6rTjpx2ApNFsufKjR5Q9vPOiKUQiSZI0OZ4DLTbvuEqSJEmSOs3EVZIkSZLUaSaukiRJkqROM3GVJEmSJHWagzNJHTZoEAJJkiRp0XjHVZIkSZLUaSaukiRJkqROM3GVJGnOJDkuyeeSfKTNn5nk9iQPJvlAkue28hPb/L72+ZZpxi1J0mpMXCVJmj9vBu7vm/8N4Jqq2go8AVzeyi8HnqiqHwSuafUkSeocE1dJkuZIks3ARcDvtvkArwD+qFW5Hri4TW9v87TPX9nqS5LUKSaukiTNl98G/g3wzTb/POCrVfV0m98PbGrTm4BHANrnT7b63ybJjiR7k+w9ePDgRsYuSdJA/hyOJElzIsk/Ah6vqruSnLtSPKBqDfHZswVVu4BdAMvLy0d8Lknj5M8BapCh77iOY6CHJFe18geSnD/ulZEkacG9HPjJJA8DN9DrIvzbwElJVi5WbwYebdP7gdMB2uffCxyaZMCSJA1jlK7C6xroIclZwKXAi4ALgHcnOW594UuSpBVVdVVVba6qLfSOuZ+oqn8GfBJ4Tat2GXBTm97d5mmff6KqvKMqSeqcoRLXMQ30sB24oaqeqqqvAPuAbeNYCUmSdFRvA96aZB+9Z1ivbeXXAs9r5W8FrpxSfJIkHdWwz7iuDPTw3W1+6IEekqwM9LAJuK3vO/uXkSRJY1RVtwK3tumHGHCxuKr+CrhkooFJknQM1rzj2j/QQ3/xgKprDfQw1AAQjlwoSZIkSeo3TFfhcQ308K3yAct8S1XtqqrlqlpeWloaeYUkSZIkSfNlzcR1jAM97AYubaMOnwlsBe4Y25pIkiRJkubSen7H9W3ADUl+Hfgc3z7Qw++1gR4O0Ut2qar7ktwIfBF4Griiqp5Zx9+XJEmSJC2AkRLX9Q70UFVXA1ePGqQkSZIkaXGN8juukiRJkiRNnImrJEmSJKnTTFwlSZIkSZ1m4ipJkiRJ6jQTV0mSJElSp5m4SpIkSZI6zcRVkiRJktRpJq6SJEmSpE4zcZUkSZIkdZqJqyRJkiSp00xcJUmaE0m+I8kdST6f5L4kv9LK35fkK0nubq+zW3mSvCvJviT3JHnpdNdAkqTBjp92AJIkaWyeAl5RVd9IcgLw6SQfa5/966r6o8PqvxrY2l4vA97T3iVJ6hTvuEqSNCeq5xtt9oT2qqMssh14f1vuNuCkJKdtdJySJI3KxFWSpDmS5LgkdwOPA3uq6vb20dWtO/A1SU5sZZuAR/oW39/KJEnqFBNXSZLmSFU9U1VnA5uBbUn+DnAV8LeBvwucArytVc+grzi8IMmOJHuT7D148OAGRS5J0upMXCVJmkNV9VXgVuCCqjrQugM/BfxXYFurth84vW+xzcCjA75rV1UtV9Xy0tLSBkcuSdKRTFwlSZoTSZaSnNSmvxN4FfClledWkwS4GLi3LbIbeH0bXfgc4MmqOjCF0CVJOipHFZYkaX6cBlyf5Dh6F6dvrKqPJPlEkiV6XYPvBn621b8ZuBDYB/wl8IYpxCxJ0ppMXCVJmhNVdQ/wkgHlr1ilfgFXbHRckiStl12FJUmSJEmdZuIqSZIkSeo0E1dJkiRJUqeZuEqSJEmSOs3EVZIkSZLUaSaukiRJkqROM3GVJEmSJHWaiaskSZIkqdNMXCVJkiRJnWbiKkmSJEnqNBNXSZIkSVKnmbhKkiRJkjrNxFWSJEmS1GkmrpIkSZKkTjt+2gFIWr8tV350YPnDOy+acCSSJEnS+HnHVZIkSZLUaSaukiTNiSTfkeSOJJ9Pcl+SX2nlZya5PcmDST6Q5Lmt/MQ2v699vmWa8UuStBoTV0mS5sdTwCuq6sXA2cAFSc4BfgO4pqq2Ak8Al7f6lwNPVNUPAte0epIkdY6JqyRJc6J6vtFmT2ivAl4B/FErvx64uE1vb/O0z1+ZJBMKV5KkoTk4k9QBqw2uJEmjSnIccBfwg8DvAF8GvlpVT7cq+4FNbXoT8AhAVT2d5EngecCfH/adO4AdAGecccZGr4IkSUdY847rOJ+XSXJVK38gyfkbtVKSJC2qqnqmqs4GNgPbgBcOqtbeB91drSMKqnZV1XJVLS8tLY0vWEmShjRMV+GxPC+T5CzgUuBFwAXAu9tVYUmSNGZV9VXgVuAc4KQkK72sNgOPtun9wOkA7fPvBQ5NNlJJkta2ZuI6xudltgM3VNVTVfUVYB+9K8GSJGkMkiwlOalNfyfwKuB+4JPAa1q1y4Cb2vTuNk/7/BNVdcQdV0mSpm2owZmSHJfkbuBxYA8jPC8DrDwv863yAcv0/60dSfYm2Xvw4MHR10iSpMV1GvDJJPcAdwJ7quojwNuAtybZR++YfG2rfy3wvFb+VuDKKcQsSdKahhqcqaqeAc5uV3E/zLE9LzP0czTALoDl5WWv+kqSNKSqugd4yYDyhxjQy6mq/gq4ZAKhSZK0LiP9HM46n5f5VvmAZSRJkiRJGmiYUYXH9bzMbuDSNurwmcBW4I5xrYgkSZIkaT4N01X4NOD6NgLwc4Abq+ojSb4I3JDk14HP8e3Py/xee17mEL2RhKmq+5LcCHwReBq4onVBliRJkiRpVWsmruN8XqaqrgauHj1MSZIkSdKiGukZV0mSJEmSJs3EVZIkSZLUaSaukiRJkqROM3GVJEmSJHWaiaskSZIkqdNMXCVJkiRJnWbiKkmSJEnqNBNXSZIkSVKnmbhKkiRJkjrNxFWSJEmS1GkmrpIkSZKkTjNxlSRpTiQ5Pcknk9yf5L4kb27lv5zkz5Lc3V4X9i1zVZJ9SR5Icv70opckaXXHTzsASZI0Nk8DP19Vn03y3cBdSfa0z66pqt/qr5zkLOBS4EXA3wL+V5IfqqpnJhq1JElr8I6rJElzoqoOVNVn2/TXgfuBTUdZZDtwQ1U9VVVfAfYB2zY+UkmSRmPiKknSHEqyBXgJcHsrelOSe5Jcl+TkVrYJeKRvsf0cPdGVJGkqTFwlSZozSb4L+CDwlqr6GvAe4AeAs4EDwDtWqg5YvAZ8344ke5PsPXjw4AZFLUnS6kxcJUmaI0lOoJe0/n5VfQigqh6rqmeq6pvAe3m2O/B+4PS+xTcDjx7+nVW1q6qWq2p5aWlpY1dAkqQBTFwlSZoTSQJcC9xfVe/sKz+tr9pPAfe26d3ApUlOTHImsBW4Y1LxSpI0LEcVliRpfrwceB3whSR3t7JfAF6b5Gx63YAfBn4GoKruS3Ij8EV6IxJf4YjCkqQuMnGVJGlOVNWnGfzc6s1HWeZq4OoNC0qSpDGwq7AkSZIkqdNMXCVJkiRJnWbiKkmSJEnqNBNXSZIkSVKnmbhKkiRJkjrNxFWSJEmS1GkmrpIkSZKkTjNxlSRJkiR1momrJEmSJKnTTFwlSZIkSZ1m4ipJkiRJ6jQTV0mSJElSp5m4SpIkSZI6zcRVkiRJktRpJq6SJEmSpE4zcZUkaU4kOT3JJ5Pcn+S+JG9u5ack2ZPkwfZ+citPkncl2ZfkniQvne4aSJI0mImrJEnz42ng56vqhcA5wBVJzgKuBG6pqq3ALW0e4NXA1vbaAbxn8iFLkrQ2E1dJkuZEVR2oqs+26a8D9wObgO3A9a3a9cDFbXo78P7quQ04KclpEw5bkqQ1rZm4jrPbUZLLWv0Hk1y2caslSdJiS7IFeAlwO/CCqjoAveQWeH6rtgl4pG+x/a1MkqROGeaO61i6HSU5BXg78DJgG/D2lWRXkiSNT5LvAj4IvKWqvna0qgPKasD37UiyN8negwcPjitMSZKGtmbiOsZuR+cDe6rqUFU9AewBLhjr2kiStOCSnEAvaf39qvpQK35spQtwe3+8le8HTu9bfDPw6OHfWVW7qmq5qpaXlpY2LnhJklYx0jOu6+x2ZHckSZI2UJIA1wL3V9U7+z7aDaw8onMZcFNf+evbYz7nAE+uHNslSeqS44eteHi3o96xcXDVAWV1lPLD/84Oel2MOeOMM4YNT5IkwcuB1wFfSHJ3K/sFYCdwY5LLgT8FLmmf3QxcCOwD/hJ4w2TDlSRpOEMlrkfrdlRVB4bsdrQfOPew8lsP/1tVtQvYBbC8vHxEYitJkgarqk8z+EIxwCsH1C/gig0NSpKkMRhmVOFxdTv6OHBekpPboEzntTJJkiRJklY1zB3XsXQ7qqpDSX4NuLPV+9WqOjSWtZAkSZIkza01E9dxdjuqquuA60YJUJIkSZK02IYenEmSJEmSumTLlR89ouzhnRdNIRJtNBNXSZIkSVMxKPGUBhnpd1wlSZIkSZo0E1dJkiRJUqeZuEqSJEmSOs3EVZIkSZLUaSaukiRJkqROc1RhacIcPU+SJEkajYmrNMf8bTNJkiTNA7sKS5IkSZI6zcRVkiRJktRpJq6SJEmSpE4zcZUkSZIkdZqJqyRJcyLJdUkeT3JvX9kvJ/mzJHe314V9n12VZF+SB5KcP52oJUlam4mrJEnz433ABQPKr6mqs9vrZoAkZwGXAi9qy7w7yXETi1SSpBGYuEqSNCeq6lPAoSGrbwduqKqnquorwD5g24YFJ0nSOpi4SpI0/96U5J7WlfjkVrYJeKSvzv5WJklS55i4SpI0394D/ABwNnAAeEcrz4C6NegLkuxIsjfJ3oMHD25MlJIkHYWJqyRJc6yqHquqZ6rqm8B7ebY78H7g9L6qm4FHV/mOXVW1XFXLS0tLGxuwJEkDmLhKkjTHkpzWN/tTwMqIw7uBS5OcmORMYCtwx6TjkyRpGMdPOwBJkjQeSf4QOBc4Ncl+4O3AuUnOptcN+GHgZwCq6r4kNwJfBJ4GrqiqZ6YRtyRJazFxlSRpTlTVawcUX3uU+lcDV29cRJIkjYddhSVJkiRJnWbiKkmSJEnqNBNXSZIkSVKnmbhKkiRJkjrNxFWSJEmS1GkmrpIkSZKkTjNxlSRJkiR1momrJEmSJKnTTFwlSZIkSZ1m4ipJkiRJ6jQTV0mSJElSp5m4SpIkSZI6zcRVkiRJktRpJq6SJEmSpE4zcZUkSZIkdZqJqyRJcyTJdUkeT3JvX9kpSfYkebC9n9zKk+RdSfYluSfJS6cXuSRJqzNxlSRpvrwPuOCwsiuBW6pqK3BLmwd4NbC1vXYA75lQjJIkjWTNxHVcV26TXNbqP5jkso1ZHUmSFltVfQo4dFjxduD6Nn09cHFf+fur5zbgpCSnTSZSSZKGN8wd1/exziu3SU4B3g68DNgGvH0l2ZUkSRvuBVV1AKC9P7+VbwIe6au3v5VJktQpayauY7pyez6wp6oOVdUTwB6OTIYlSdJkZUBZHVEp2ZFkb5K9Bw8enEBYkiR9u2N9xnXUK7dDX9H14ChJ0tg9ttIFuL0/3sr3A6f31dsMPHr4wlW1q6qWq2p5aWlpw4OVJOlw4x6cabUrt0Nd0QUPjpIkbYDdwMr4EpcBN/WVv76NUXEO8OTKhWlJkrrkWBPXUa/cDnVFV5IkrU+SPwQ+A/xwkv1JLgd2Aj+R5EHgJ9o8wM3AQ8A+4L3Av5xCyJIkren4Y1xu5crtTo68cvumJDfQG4jpyao6kOTjwL/rG5DpPOCqYw9bmg1brvzotEOQtGCq6rWrfPTKAXULuGJjI5Ikaf3WTFzbldtzgVOT7Kc3OvBO4MZ2FfdPgUta9ZuBC+lduf1L4A0AVXUoya8Bd7Z6v1pVhw/4JEmSJEnSEdZMXMd15baqrgOuGyk6SZIkSdLCG/fgTJIkSZIkjZWJqyRJkiSp0451cCZJM2rQgFEP77xoCpFIkiRJw/GOqyRJkiSp00xcJUmSJEmdZuIqSZIkSeo0E1dJkiRJUqc5OJMkSZKkueFAlPPJO66SJEmSpE4zcZUkSZIkdZqJqyRJkiSp00xcJUmSJEmdZuIqSZIkSeo0E1dJkiRJUqf5cziSJC2AJA8DXweeAZ6uquUkpwAfALYADwM/XVVPTCtGSfNt0M/USMPyjqskSYvjH1bV2VW13OavBG6pqq3ALW1ekqTOMXGVJGlxbQeub9PXAxdPMRZJklZlV2FpTOz+IqnjCvjjJAX8l6raBbygqg4AVNWBJM8ftGCSHcAOgDPOOGNS8UqS9C0mrpIkLYaXV9WjLTndk+RLwy7YktxdAMvLy7VRAUqStBq7CkuStACq6tH2/jjwYWAb8FiS0wDa++PTi1CSpNV5x1XSwG7OD++8aAqRSNoISf4G8Jyq+nqbPg/4VWA3cBmws73fNL0oJUlanYmrJEnz7wXAh5NA79j/B1X1P5PcCdyY5HLgT4FLphijJEmrMnGVJGnOVdVDwIsHlP8F8MrJRyRJ0mh8xlWSJEmS1GkmrpIkSZKkTjNxlSRJkiR1ms+4SiMaNAKvJEmSpI3jHVdJkiRJUqeZuEqSJEmSOs3EVZIkSZLUaT7jKmmg1Z7lfXjnRROORJIkzZIujgfiec3sM3GVjqKLO15JkiRp0dhVWJIkSZLUad5xlTSSQXeh7WYjSZKkjWTiKkmSJGkheUF+dpi4aiH57KokSfPDgXek+WfiKmndvFopSdoI6z2+eHzaeN4M0KSYuGquuTOVJGm8Rjm2DkoS13t31GP79Nj2mqaJJ65JLgD+A3Ac8LtVtXPSMWj2uePsPq9yS7PDY7Okwy3yuZbnMN000cQ1yXHA7wA/AewH7kyyu6q+OMk4NFnD/vMv8g5yUaz3Kr2k8fPYrKNZ77G5i8f2LiYlGxGT51/j1cXtZtFM+o7rNmBfVT0EkOQGYDvgwXFM1rvzmVQ3HXeSWsuw28hq26wHGGloHpvHZJIDBE0q0Zmkaf/9Ya3n+LRR6zjs985KG8+69V408HxlsFTV5P5Y8hrggqr6523+dcDLqupNg+ovLy/X3r17x/K3/UeVNI/We2K0aAfHJHdV1fK04+iSrh2bu5jkbcSdq/WexG7ExT1Js2Xad9DHtb8e9tg86cT1EuD8ww6O26rqX/XV2QHsaLM/DDwwhj99KvDnY/ieRWF7Dc+2Go3tNTzbajTDttf3VdXSRgczS6Z4bJ6EWfw/msWYYTbjnsWYwbgnaRZjhtmLe6hj86S7Cu8HTu+b3ww82l+hqnYBu8b5R5Ps9Qr78Gyv4dlWo7G9hmdbjcb2WpepHJsnYRa3i1mMGWYz7lmMGYx7kmYxZpjduNfynAn/vTuBrUnOTPJc4FJg94RjkCRJz/LYLEnqvIneca2qp5O8Cfg4vSH3r6uq+yYZgyRJepbHZknSLJj477hW1c3AzRP+szPXvWnKbK/h2Vajsb2GZ1uNxvZahykdmydhFreLWYwZZjPuWYwZjHuSZjFmmN24j2qigzNJkiRJkjSqST/jKkmSJEnSSGY2cU1yXZLHk9zbV3Z2ktuS3J1kb5JtrfzcJE+28ruT/FLfMhckeSDJviRXTmNdNtoqbfXiJJ9J8oUk/yPJ9/R9dlVrjweSnN9XPvdtBaO1V5ItSf5v37b1n/uW+ZFWf1+SdyXJNNZnIyU5Pcknk9yf5L4kb27lpyTZk+TB9n5yK09ri31J7kny0r7vuqzVfzDJZdNap410DO21sPuuo7TVJW3+m0mWD1tmofddi2CV/fMH+v5HHk5yd99nndgmRom7K8eVVWJe7TyrM/v2EePuxD52lZg7f542Stwd2q5n8rzlGOLuxLY9dlU1ky/gHwAvBe7tK/tj4NVt+kLg1jZ9LvCRAd9xHPBl4PuB5wKfB86a9rpNqK3uBH68Tb8R+LU2fVZrhxOBM1v7HLcobXUM7bWlv95h33MH8KNAgI+tbJvz9AJOA17apr8b+JO2Df174MpWfiXwG236wtYWAc4Bbm/lpwAPtfeT2/TJ016/DrTXwu67jtJWL6T3O6K3Ast99Rd+37UIr0H758M+fwfwS13bJkaMuxPHlUExs/p5Vmf27SPG3Yl97Coxd/48bcS4u7Jdz+R5yzHE3Ylte9yvmb3jWlWfAg4dXgysXJH6Xg77HboBtgH7quqhqvpr4AZg+1gD7YBV2uqHgU+16T3AP23T24EbquqpqvoKsI9eOy1EW8HI7TVQktOA76mqz1RvT/F+4OJxxzptVXWgqj7bpr8O3A9sordtXN+qXc+z674deH/13Aac1NrqfGBPVR2qqifotfEFE1yViTiG9lrN3P8/rtZWVXV/VT0wYJGF33ctglX2z0Dvzgjw08AftqLObBMjxj3QpI8rI55ndWbfPovnh7N6njaL50uzet7i+UPPzCauq3gL8JtJHgF+C7iq77MfTfL5JB9L8qJWtgl4pK/O/la2CO4FfrJNX8KzPz6/WpssclvB6u0FcGaSzyX530l+rJVtotdGK+a+vZJsAV4C3A68oKoOQG9nCzy/VXP7aoZsL3DfdXhbrcZtSz8GPFZVD7b5WdkmDo8buntcWe08q+ttPYvnh7N6njYz50uzet6yyOcP85a4/gvg56rqdODngGtb+WeB76uqFwP/EfjvrXxQH/pFGWb5jcAVSe6i1+Xgr1v5am2yyG0Fq7fXAeCMqnoJ8FbgD9rzHAvVXkm+C/gg8Jaq+trRqg4oW7jta4T2Wvh9l9uWRvBavv2u5axsE4fH3eXjymrnWV1v61k8P5zV87SZOF+a1WPLop8/zFviehnwoTb93+jdDqeqvlZV32jTNwMnJDmV3lWG/itBm1m7+8hcqKovVdV5VfUj9A6YX24frdYmC9tWsHp7ta46f9Gm72rlP0SvvTb3fcXctleSE+jtRH+/qlb+/x5rXWlWugE93soXfvsapb0Wfd+1SlutZuG3rUWW5HjgnwAf6Cvu/DYxKO6OH1cGnmfR/baeufPDWT1Pm4XzpVk9b/H8Yf4S10eBH2/TrwAeBEjyN1dGKEtvJLnnAH9B7wHyrUnOTPJc4FJg98SjnoIkz2/vzwH+LbAyuttu4NIkJyY5E9hK76H5hW0rWL29kiwlOa5Nfz+99nqoddf4epJz2rb3euCmqQS/gdq6XQvcX1Xv7PtoN70TBdr7TX3lr0/POcCTra0+DpyX5OQ2It55rWyujNpei7zvOkpbrcZ912J7FfClqurvcjgL28QRcXf8uDLwPIvu79tn7vxwVs/Tun6+NKvnLZ4/NNWBEaKO5UXvKs4B4P/Ru3pwOfD3gbvojZB1O/Ajre6bgPta+W3A3+v7ngvpjcz1ZeAXp71eE2yrN7f1/hNgJ5C++r/Y2uMB+kZ2W4S2GrW96A06sLJtfRb4x33fs0zvWY8vA/+pv43n5dX+5wq4B7i7vS4EngfcQu/k4BbglFY/wO+0NvkC3z4q7BvpDTKxD3jDtNetI+21sPuuo7TVT7X/y6eAx4CP9y2z0PuuRXgN2j+38vcBPzugfie2iVHi7spxZVDMrH6e1Zl9+4hxd2Ifu0rMnT9PGyXuDm3XM3necgxxd2LbHvdrZWOSJEmSJKmT5q2rsCRJkiRpzpi4SpIkSZI6zcRVkiRJktRpJq6SJEmSpE4zcZUkSZIkdZqJqyRJkiSp00xcJUmSJEmdZuIqSZIkSeq0/w+U3YZ1uS8RSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(1,2,1)\n",
    "plt.hist(a,bins=70);\n",
    "plt.subplot(1,2,2)\n",
    "plt.hist(b,bins=70);\n",
    "\n",
    "fig=plt.gcf();\n",
    "fig.set_size_inches(16,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}