diff --git a/DNN.ipynb b/DNN.ipynb new file mode 100644 index 0000000..9f07f7e --- /dev/null +++ b/DNN.ipynb @@ -0,0 +1,604 @@ +{ + "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": [ + "
" + ] + }, + "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 +}