detecting _regional_phases_HZ_git.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "# Detecting Regional Seismic  Phases  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the input from the SAC files\n",
    "This step can be for three components record or one component record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "\n",
    "sac_dir = '/localdata/rayan/dataprc_sac_hvsr_data/2012*/*.BHZ.D.SAC'\n",
    "out_sac = open('sacfiles.txt', 'w')\n",
    "for sacfile in (glob.glob(sac_dir)):\n",
    "    if len(glob.glob(sacfile[:-9]+'*.SAC')) == 3:\n",
    "        #out_sac.write('\\n'.join([sacfile[:-9]+\"BHZ.D.SAC  \" +\n",
    "        #                    sacfile[:-9]+\"BHN.D.SAC  \"+sacfile[:-9]+\"BHE.D.SAC  \"+'\\n']))\n",
    "        out_sac.write('\\n'.join([sacfile[:-9]+\"*BHZ.D.SAC\"+'\\n']))\n",
    "    else:\n",
    "        print ('there are no BHE BHN BHZ')\n",
    "out_sac.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/localdata/rayan/generalized-phase-detection\n"
     ]
    }
   ],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the input of the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import string\n",
    "import time\n",
    "import argparse as ap\n",
    "import sys\n",
    "import os\n",
    "import obsplus as op\n",
    "import numpy as np\n",
    "import obspy.core as oc\n",
    "import matplotlib as mpl\n",
    "import pylab as plt\n",
    "import h5py\n",
    "import time\n",
    "from scipy.signal import spectrogram \n",
    "#import tensorflow as tf\n",
    "import tensorflow.compat.v1 as tf\n",
    "tf.disable_v2_behavior()\n",
    "##\n",
    "freq_min=1\n",
    "freq_max=10\n",
    "filter_data = True\n",
    "decimate_data = True\n",
    "\n",
    "#\n",
    "# Number of GPUs to use (if any)\n",
    "n_gpu = 1 \n",
    "#####################\n",
    "\n",
    "\n",
    "half_dur = 0.64 \n",
    "only_dt = 0.01\n",
    "n_win = int(half_dur/only_dt)\n",
    "n_feat = 2*n_win\n",
    "n_feat_spec=65\n",
    "n_shift = n_feat  # Number of samples to shift the sliding window at a time\n",
    "sample_rate=1/only_dt\n",
    "##using gpu here is much faster\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\" \n",
    "os.environ['CUDA_VISIBLE_DEVICES']=\"0\" ##0=gpu,\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.allow_growth=True\n",
    "sess = tf.Session(config=config)\n",
    "#using obsplus makes it faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#functions\n",
    "def sliding_window(data, size, stepsize=128, padded=False, axis=-1, copy=False):\n",
    "    #stepsize =400 so devided with 400 sample not over laping. \n",
    "    shape = list(data[0].shape)\n",
    "    shape[axis] = np.floor(data[0][0].shape[axis] / stepsize - size / stepsize + 1).astype(int)\n",
    "    shape.append(size)\n",
    "\n",
    "    strides = list(data[0].strides)\n",
    "    strides[axis] *= stepsize\n",
    "    strides.append(data.strides[axis])\n",
    "    print((\"shape:{0},strides:{1}\").format(shape,strides))\n",
    "\n",
    "    strided = np.lib.stride_tricks.as_strided(data[0], shape=shape, strides=strides)\n",
    "    \n",
    "    if copy:\n",
    "        return strided.copy()\n",
    "    else:\n",
    "        return strided\n",
    "    \n",
    "def spect(x):\n",
    "    frequencies, times, spectrog = spectrogram(np.hstack(x),\n",
    "                                              sample_rate,nfft=128,noverlap=0, \n",
    "                                              nperseg=128,window='triang',\n",
    "                                              detrend='linear',\n",
    "                                              scaling='spectrum'\n",
    "                                             \n",
    "                                             )\n",
    "    return spectrog\n",
    "\n",
    "#\n",
    "def phases_to_train(st,freq_min,freq_max,op):\n",
    "    st.detrend(type='linear')\n",
    "   \n",
    "    st.filter(type='bandpass', freqmin=freq_min, freqmax=freq_max)\n",
    "    print(\" filtering is done###\")\n",
    "   \n",
    "    st.interpolate(100.0)\n",
    "    print(\" decimateing is done###\")\n",
    "\n",
    "    st_dict = {f'event{x}': st.copy() for x in range(1)}\n",
    "    from_dict = op.obspy_to_array(st_dict)\n",
    "    chan = st[0].stats.channel\n",
    "    seedid=from_dict.coords['seed_id'].values[0]\n",
    "    sr = from_dict.stats['event0'][seedid].sampling_rate\n",
    "    dt = from_dict.stats['event0'][seedid].delta\n",
    "    net = from_dict.stats['event0'][seedid].network\n",
    "    sta = from_dict.stats['event0'][seedid].station\n",
    "\n",
    "    t_p=from_dict.stats['event0'][seedid].sac.t1\n",
    "    t_s1=from_dict.stats['event0'][seedid].sac.t5\n",
    "    t_s2=from_dict.stats['event0'][seedid].sac.t4\n",
    "    t_lg_b=from_dict.stats['event0'][seedid].sac.t2\n",
    "    t_lg_e=from_dict.stats['event0'][seedid].sac.t3\n",
    "\n",
    "    Noise_data=from_dict.sel(time=slice(0, t_p)).data\n",
    "\n",
    "    P_data=from_dict.sel(time=slice(t_p, t_s1-2)).data\n",
    "\n",
    "    S_data=from_dict.sel(time=slice(t_s2, t_lg_b)).data\n",
    "\n",
    "    Lg_data=from_dict.sel(time=slice(t_lg_b, t_lg_e)).data\n",
    "    \n",
    "   \n",
    "    Noise_train1=sliding_window(Noise_data, 128, stepsize=128)\n",
    "    shp=(Noise_train1.shape[0],Noise_train1.shape[1],65)\n",
    "    N_T=[spect(Noise_train1[:,i,:]) for i in range(0,Noise_train1.shape[1])]\n",
    "    Noise_train=np.reshape(N_T, shp)\n",
    "    #\n",
    "    P_train1=sliding_window(P_data, 128, stepsize=128)\n",
    "    shp=(P_train1.shape[0],P_train1.shape[1],65)\n",
    "    P_T=[spect(P_train1[:,i,:]) for i in range(0,P_train1.shape[1])]\n",
    "    P_train=np.reshape(P_T, shp)\n",
    "    #\n",
    "    S_train1=sliding_window(S_data, 128, stepsize=128)\n",
    "    shp=(S_train1.shape[0],S_train1.shape[1],65)\n",
    "    S_T=[spect(S_train1[:,i,:]) for i in range(0,S_train1.shape[1])]\n",
    "    S_train=np.reshape(S_T, shp)\n",
    "    #\n",
    "    LG_train1=sliding_window(Lg_data, 128, stepsize=128)\n",
    "    shp=(LG_train1.shape[0],LG_train1.shape[1],65)\n",
    "    L_T=[spect(LG_train1[:,i,:]) for i in range(0,LG_train1.shape[1])]\n",
    "    LG_train=np.reshape(L_T, shp)\n",
    "    #\n",
    "    tr_win0 = np.zeros((Noise_train.shape[1], n_feat_spec, 1))\n",
    "    tr_type0 = np.zeros((Noise_train.shape[1]))\n",
    "    tr_type0[:] = 0\n",
    "    tr_win0[:, :, 0] = Noise_train[0]\n",
    "    #tr_win0[:, :, 1] = Noise_train[1]\n",
    "    #tr_win0[:, :, 2] = Noise_train[2]\n",
    "    tr_win0 =tr_win0/np.max(np.abs(tr_win0),axis=(1,2))[:,None,None] #normalizine \n",
    "    \n",
    "#\n",
    "    tr_win1 = np.zeros((P_train.shape[1], n_feat_spec, 1))\n",
    "    tr_type1 = np.zeros((P_train.shape[1]))\n",
    "    tr_type1[:] = 1\n",
    "    tr_win1[:, :, 0] = P_train[0]\n",
    "    #tr_win1[:, :, 1] = P_train[1]\n",
    "    #tr_win1[:, :, 2] = P_train[2]\n",
    "    tr_win1 =tr_win1/np.max(np.abs(tr_win1),axis=(1,2))[:,None,None]#normalizine \n",
    "    #\n",
    "    tr_win2 = np.zeros((S_train.shape[1], n_feat_spec, 1))\n",
    "    tr_type2 = np.zeros((S_train.shape[1]))\n",
    "    tr_type2[:] = 2\n",
    "    tr_win2[:, :, 0] = S_train[0]\n",
    "    #tr_win2[:, :, 1] = S_train[1]\n",
    "    #tr_win2[:, :, 2] = S_train[2]\n",
    "    tr_win2 =tr_win2/np.max(np.abs(tr_win2),axis=(1,2))[:,None,None]#normalizine                         \n",
    "                            \n",
    "                            \n",
    "#\n",
    "    tr_win3= np.zeros((LG_train.shape[1], n_feat_spec, 1))\n",
    "    tr_type3 = np.zeros((LG_train.shape[1]))\n",
    "    tr_type3[:] = 3\n",
    "    tr_win3[:, :, 0] = LG_train[0]\n",
    "    #tr_win3[:, :, 1] = LG_train[1]\n",
    "    #tr_win3[:, :, 2] = LG_train[2]\n",
    "    tr_win3 =tr_win3/np.max(np.abs(tr_win3),axis=(1,2))[:,None,None]#normalizine                         \n",
    "   \n",
    "    return Noise_train1,Noise_train ,tr_win0, tr_type0,tr_win1, tr_type1,tr_win2, tr_type2,tr_win3, tr_type3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == \"__main__\":    \n",
    "    \n",
    "    tr_win=np.zeros([1,n_feat_spec,1])\n",
    "    tr_type=np.zeros([1])\n",
    "    time0=time.time()\n",
    "   \n",
    "    with open('sacfiles.txt') as f:\n",
    " \n",
    "       \n",
    "        for count,line in enumerate (f.read().splitlines()):\n",
    "\n",
    "            time1=time.time()\n",
    "            ff=oc.read(str(line))\n",
    "            Noise_train1,Noise_train,tr_win0, tr_type0,tr_win1, tr_type1,tr_win2, tr_type2,tr_win3, tr_type3=phases_to_train(ff,freq_min,freq_max,op)\n",
    "\n",
    "            #np.fromiter(tr_type4,dtype=\"float\")\n",
    "            tr_win=np.concatenate((tr_win,tr_win0,tr_win1,tr_win2,tr_win3))\n",
    "            tr_type=np.concatenate((tr_type,tr_type0,tr_type1,tr_type2,tr_type3))\n",
    "            time2=time.time()\n",
    "            total_time=time2-time1\n",
    "            print('time for line number {} is {} second'.format(count,total_time))\n",
    "          \n",
    "    \n",
    "time3=time.time()\n",
    "time=time3-time0\n",
    "print('all_time is {} second'.format(time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "with h5py.File('data100fs_0.5_10HZ_spec.h5', 'w') as  hf:\n",
    "    hf.create_dataset('X', data=tr_win)\n",
    "    hf.create_dataset('Y', data=tr_type)\n",
    "#read the hf5 data to check the output\n",
    "\n",
    "\"\"\"\n",
    "with h5py.File('data.h5', 'r') as of:\n",
    "    X=of.get('X')\n",
    "    Y=of.get('Y')\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "with h5py.File('data100fs_0.5_10HZ_spec.h5', 'r') as of:\n",
    "    X=of.get('X')\n",
    "    Y=of.get('Y')\n",
    "    print (of['Y'])\n",
    "    print (of['X'])\n",
    "    print (of['X'][700].shape)\n",
    "    #print (np.unique(of['X'][1:]))\n",
    "    sample_count1=of['X'].shape[0]\n",
    "    print(sample_count1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#based on [Ross et al., 2018]\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "import warnings\n",
    "import os\n",
    "import matplotlib as mpl\n",
    "# mpl.use('Agg')\n",
    "import keras\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Dropout, Activation, Flatten\n",
    "from tensorflow.keras.layers import Conv1D, MaxPooling1D, LSTM\n",
    "from tensorflow.keras import losses\n",
    "from tensorflow.keras.models import model_from_json\n",
    "from keras.utils.np_utils import to_categorical\n",
    "#import tensorflow as tf\n",
    "import tensorflow.compat.v1 as tf\n",
    "tf.disable_v2_behavior()\n",
    "\n",
    "# important steps\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.allow_growth=True\n",
    "sess = tf.Session(config=config)\n",
    "#\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "\n",
    "warnings.simplefilter(\"ignore\")\n",
    "\n",
    "#hdf_file = \"scsn_ps_2000_2017_shuf.hdf5\"\n",
    "hdf_file = \"data100fs_0.5_10HZ_spec.h5\"\n",
    "sample_count = sample_count1 #2421000 ##total number of 400 windows we have. \n",
    "show_plot = 1\n",
    "train = 1\n",
    "total_dur = 1.28\n",
    "win_dur = 1.28  \n",
    "dt = 0.01\n",
    "phase_dir = \"phase\"\n",
    "#n_win = int(win_dur / dt)\n",
    "n_win=int(65)\n",
    "#n_full = int(total_dur / dt)\n",
    "n_full=int(65)\n",
    "n_start = (n_full - n_win) // 2\n",
    "gpu_count=1\n",
    "batchsize=1000\n",
    "epoches=15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#functions \n",
    "def make_parallel(model, gpu_count):\n",
    "    \"\"\"\"\n",
    "    This should be only used when we have many GPU\n",
    "    \"\"\"\n",
    "    from tensorflow.keras.layers import add\n",
    "    from tensorflow.keras.layers import Lambda\n",
    "    from tensorflow.keras.models import Model\n",
    "\n",
    "    def get_slice(data, idx, parts):\n",
    "        shape = tf.shape(data)\n",
    "        size = tf.concat([shape[:1] // parts, shape[1:]], axis=0)\n",
    "        stride = tf.concat([shape[:1] // parts, shape[1:]*0], axis=0)\n",
    "        start = stride * idx\n",
    "        return tf.slice(data, start, size)\n",
    "\n",
    "    outputs_all = []\n",
    "    for i in range(len(model.outputs)):\n",
    "        outputs_all.append([])\n",
    "\n",
    "    # Place a copy of the model on each GPU, each getting a slice of the batch\n",
    "    for i in range(gpu_count):\n",
    "        with tf.device('/gpu:%d' % i):\n",
    "            with tf.name_scope('tower_%d' % i) as scope:\n",
    "\n",
    "                inputs = []\n",
    "                # Slice each input into a piece for processing on this GPU\n",
    "                for x in model.inputs:\n",
    "                    input_shape = tuple(x.get_shape().as_list())[1:]\n",
    "                    slice_n = Lambda(get_slice, output_shape=input_shape, arguments={\n",
    "                                     'idx': i, 'parts': gpu_count})(x)\n",
    "                    inputs.append(slice_n)\n",
    "\n",
    "                outputs = model(inputs)\n",
    "\n",
    "                if not isinstance(outputs, list):\n",
    "                    outputs = [outputs]\n",
    "\n",
    "                # Save all the outputs for merging back together later\n",
    "                for l in range(len(outputs)):\n",
    "                    outputs_all[l].append(outputs[l])\n",
    "                print('ff')\n",
    "                print(outputs_all)\n",
    "\n",
    "    # merge outputs on CPU\n",
    "    with tf.device('/cpu:0'):\n",
    "        merged = []\n",
    "        for outputs in outputs_all:\n",
    "            #merged.append(merge(outputs, mode='concat', concat_axis=0))\n",
    "            merged.append(tf.keras.layers.Concatenate(axis=0)(outputs))\n",
    "           \n",
    "\n",
    "        return Model(model.inputs, merged)\n",
    "\n",
    "\n",
    "def train_model(model, seis_data_fn=\"scsn_all.hdf5\",\n",
    "                validation_ratio=0.25, batch_size=32,\n",
    "                sample_count=None, callbacks=None,epochs=100):\n",
    "    \"\"\" Train the classification model:\n",
    "\n",
    "    \"\"\"\n",
    "    with h5py.File(seis_data_fn, \"r\") as seis_data:\n",
    "        if sample_count is None:\n",
    "            sample_count = int(seis_data.attrs[\"sample_count\"])\n",
    "        sample_idxs = range(0, sample_count)\n",
    "        #sample_idxs = np.random.permutation(sample_idxs)\n",
    "        train_samp_idxs = sample_idxs[0:int(\n",
    "            (1 - validation_ratio) * sample_count)]\n",
    "        val_samp_idxs = sample_idxs[int(\n",
    "            (1 - validation_ratio) * sample_count):]\n",
    "        print(\"Using %d samples for training set\" % len(train_samp_idxs))\n",
    "        print(\"Using %d samples for validation set\" % len(val_samp_idxs))\n",
    "        train_seq_gen = gen_train_seq(batch_size=batch_size,\n",
    "                                      seis_data=seis_data,\n",
    "                                      train_samp_idxs=train_samp_idxs)\n",
    "        val_seq_gen = gen_val_seq(batch_size=batch_size,\n",
    "                                  seis_data=seis_data,\n",
    "                                  val_samp_idxs=val_samp_idxs)\n",
    "        hist = model.fit_generator(generator=train_seq_gen,\n",
    "                                   validation_data=val_seq_gen,\n",
    "                                   steps_per_epoch=len(\n",
    "                                       train_samp_idxs)/batch_size,\n",
    "                                   validation_steps=len(\n",
    "                                       val_samp_idxs)/batch_size,\n",
    "                                   epochs=epochs,\n",
    "                                   max_queue_size=1000,\n",
    "                                   verbose=1,\n",
    "                                   workers=1,\n",
    "                                   use_multiprocessing=False,\n",
    "                                   callbacks=callbacks)\n",
    "        return hist\n",
    "\n",
    "\n",
    "def predict_test(model, seis_data_fn=\"scsn_all.hdf5\",\n",
    "                 validation_ratio=0.3, batch_size=32,\n",
    "                 sample_count=None):\n",
    "    \"\"\" predict the classification model\n",
    "    \"\"\"\n",
    "    with h5py.File(seis_data_fn, \"r\") as seis_data:\n",
    "        if sample_count is None:\n",
    "            sample_count = int(seis_data.attrs[\"sample_count\"])\n",
    "        sample_idxs = range(0, sample_count)\n",
    "        #sample_idxs = np.random.permutation(sample_idxs)\n",
    "        val_samp_idxs = sample_idxs[int(\n",
    "            (1 - validation_ratio) * sample_count):]\n",
    "        print(\"Using %d samples for validation set\" % len(val_samp_idxs))\n",
    "        val_seq_gen = gen_val_seq(batch_size=batch_size,\n",
    "                                  seis_data=seis_data,\n",
    "                                  val_samp_idxs=val_samp_idxs)\n",
    "        y_pred = model.predict_generator(generator=val_seq_gen,\n",
    "                                         steps=len(val_samp_idxs)/batch_size,\n",
    "                                         max_queue_size=1000,\n",
    "                                         verbose=1,\n",
    "                                         workers=1,\n",
    "                                         use_multiprocessing=False)\n",
    "        val_seq_gen_y = gen_val_seq_y(batch_size=batch_size,\n",
    "                                      seis_data=seis_data,\n",
    "                                      val_samp_idxs=val_samp_idxs)\n",
    "        Y_test = []\n",
    "        snr = []\n",
    "        dist = []\n",
    "        mag = []\n",
    "        count = 0\n",
    "        for vals in val_seq_gen_y:\n",
    "            count += 1\n",
    "            Y_test += vals.tolist()\n",
    "            #Y_test += vals[0].tolist()\n",
    "            #snr += vals[1].tolist()\n",
    "            #dist += vals[2].tolist()\n",
    "            #mag += vals[3].tolist()\n",
    "            if count >= len(val_samp_idxs)/batch_size:\n",
    "                break\n",
    "\n",
    "    # , np.array(snr), np.array(dist), np.array(mag)\n",
    "    return np.array(y_pred), np.array(Y_test).reshape(len(Y_test), 1)\n",
    "\n",
    "\n",
    "def gen_train_seq(batch_size, seis_data, train_samp_idxs):\n",
    "    \"\"\" Generates training data \n",
    "    \"\"\"\n",
    "    while True:\n",
    "        # generate sequences for training\n",
    "        train_samp_count = len(train_samp_idxs)\n",
    "        batches = int(train_samp_count / batch_size)\n",
    "        remainder_samples = train_samp_count % batch_size\n",
    "        if remainder_samples:\n",
    "            batches = batches + 1\n",
    "        # generate batches of samples\n",
    "        for idx in range(0, batches):\n",
    "            if idx == batches - 1:\n",
    "                batch_idxs = train_samp_idxs[idx * batch_size:]\n",
    "            else:\n",
    "                batch_idxs = train_samp_idxs[idx * batch_size:idx\n",
    "                                             * batch_size + batch_size]\n",
    "            batch_idxs = sorted(batch_idxs)\n",
    "\n",
    "            X = seis_data[\"X\"][batch_idxs, n_start:n_start+n_win]\n",
    "            Y = seis_data[\"Y\"][batch_idxs]\n",
    "            yield (np.array(X).reshape((X.shape[0], X.shape[1], 1)),\n",
    "                   to_categorical(Y, 4))#number of the classes\n",
    "\n",
    "\n",
    "def gen_val_seq(batch_size, seis_data, val_samp_idxs):\n",
    "    \"\"\" Generates validation data\n",
    "    \"\"\"\n",
    "    while True:\n",
    "        # generate sequences for validation\n",
    "        val_samp_count = len(val_samp_idxs)\n",
    "        batches = int(val_samp_count / batch_size)\n",
    "        remainder_samples = val_samp_count % batch_size\n",
    "        if remainder_samples:\n",
    "            batches = batches + 1\n",
    "        # generate batches of samples\n",
    "        for idx in range(0, batches):\n",
    "            if idx == batches - 1:\n",
    "                batch_idxs = val_samp_idxs[idx * batch_size:]\n",
    "            else:\n",
    "                batch_idxs = val_samp_idxs[idx * batch_size:idx\n",
    "                                           * batch_size + batch_size]\n",
    "            batch_idxs = sorted(batch_idxs)\n",
    "\n",
    "            X = seis_data[\"X\"][batch_idxs, n_start:n_start+n_win]\n",
    "            Y = seis_data[\"Y\"][batch_idxs]\n",
    "            yield (np.array(X).reshape((X.shape[0], X.shape[1], 1)),\n",
    "                   to_categorical(Y, 4))#classes\n",
    "\n",
    "\n",
    "def gen_val_seq_y(batch_size, seis_data, val_samp_idxs):\n",
    "    \"\"\" Generates validation data\n",
    "    \"\"\"\n",
    "    while True:\n",
    "        # generate sequences for validation\n",
    "        val_samp_count = len(val_samp_idxs)\n",
    "        batches = int(val_samp_count / batch_size)\n",
    "        remainder_samples = val_samp_count % batch_size\n",
    "        if remainder_samples:\n",
    "            batches = batches + 1\n",
    "        # generate batches of samples\n",
    "        for idx in range(0, batches):\n",
    "            if idx == batches - 1:\n",
    "                batch_idxs = val_samp_idxs[idx * batch_size:]\n",
    "            else:\n",
    "                batch_idxs = val_samp_idxs[idx * batch_size:idx\n",
    "                                           * batch_size + batch_size]\n",
    "            batch_idxs = sorted(batch_idxs)\n",
    "\n",
    "            Y = seis_data[\"Y\"][batch_idxs]\n",
    "\n",
    "            yield np.array(Y) \n",
    "\n",
    "\n",
    "def WritePhaseFile(fdir, evid, fm):\n",
    "    \"\"\" Writes out an STP format phase file using picked first motions\"\"\"\n",
    "    line_count = 0\n",
    "    picks = {}\n",
    "    g = open(\"%s/%s.phase\" % (fdir, evid), 'w')\n",
    "    with open(\"%s/%s.phase\" % (phase_dir, evid), 'r') as f:\n",
    "        for line in f:\n",
    "            if line_count > 0:\n",
    "                temp = line.split()\n",
    "                net = temp[0]\n",
    "                sta = temp[1]\n",
    "                chan = temp[2]\n",
    "                if chan[0] == 'B':\n",
    "                    chan = list(chan)\n",
    "                    chan[0] = 'H'\n",
    "                    chan = \"\".join(chan)\n",
    "                phase = temp[7]\n",
    "                pick = temp[12]\n",
    "                pol = temp[8][0]\n",
    "                temp2 = line\n",
    "                if phase == 'P':\n",
    "                    if chan[2] != 'Z':\n",
    "                        chan = list(chan)\n",
    "                        chan[2] = 'Z'\n",
    "                        chan = \"\".join(chan)\n",
    "                    temp2 = list(temp2)\n",
    "                    if (net, sta, chan) in fm:\n",
    "                        temp2[49] = fm[(net, sta, chan)]\n",
    "                    else:\n",
    "                        temp2[49] = '.'\n",
    "                    temp2 = \"\".join(temp2)\n",
    "                g.write(temp2)\n",
    "            else:\n",
    "                g.write(line)\n",
    "            line_count += 1\n",
    "    g.close()\n",
    "    return picks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if __name__ == \"__main__\":\n",
    "    import tensorflow.compat.v1 as tf\n",
    "    tf.disable_v2_behavior()\n",
    "\n",
    "    from sklearn.externals import joblib\n",
    "    import h5py\n",
    "    from tensorflow.compat.v1.keras.layers import BatchNormalization\n",
    "    from tensorflow.compat.v1.keras.layers import MaxPooling1D\n",
    "    from tensorflow.compat.v1.keras.layers import GaussianNoise\n",
    "   \n",
    "    if train:\n",
    "        model = Sequential()\n",
    "        model.add(Conv1D(32, 21, padding='same', input_shape=(n_win, 1)))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "        model.add(MaxPooling1D(pool_size=2))\n",
    "\n",
    "        model.add(Conv1D(64, 15, padding='same'))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "        model.add(MaxPooling1D(pool_size=2))\n",
    "\n",
    "        model.add(Conv1D(128, 11, padding='same'))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "        model.add(MaxPooling1D(pool_size=2))\n",
    "\n",
    "        model.add(Conv1D(256, 9, padding='same'))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "        model.add(MaxPooling1D(pool_size=2))\n",
    "   \n",
    "        model.add(Flatten())\n",
    "\n",
    "        model.add(Dense(200))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "\n",
    "        model.add(Dense(200))\n",
    "        model.add(BatchNormalization())\n",
    "        model.add(Activation('relu'))\n",
    "\n",
    "        model.add(Dense(4)) #4 is number the output samples N pn sn Lg\n",
    "        model.add(Activation('softmax')) # software for multiple sampels \n",
    "\n",
    "        print(model.summary())\n",
    "        adam = tf.keras.optimizers.Adam()\n",
    "\n",
    "        #model = make_parallel(model, gpu_count)#in case we have many GPUS\n",
    "\n",
    "        model.compile(loss=losses.categorical_crossentropy,\n",
    "                      optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "        early_stop = tf.keras.callbacks.EarlyStopping(patience=5,\n",
    "                                                   monitor='val_acc') #monitoring the accuracy\n",
    "        model_cp = tf.keras.callbacks.ModelCheckpoint(\"model_pol_best.hdf5\",\n",
    "                                                   # model_cp = keras.callbacks.ModelCheckpoint(\"model.{epoch:02d}-{val_loss:.2f}.hdf5\",\n",
    "                                                   save_best_only=True)\n",
    "       \n",
    "        hist = train_model(model, seis_data_fn=hdf_file,\n",
    "                           validation_ratio=0.25, batch_size=batchsize,\n",
    "                           callbacks=[early_stop, model_cp],\n",
    "                           sample_count=sample_count,epochs=epoches)\n",
    "\n",
    "        hist = hist.history\n",
    "        joblib.dump(hist, \"hist_pol.pkl\")\n",
    "\n",
    "        print(\"Saved picking model to disk\")\n",
    "\n",
    "        # serialize model to JSON\n",
    "        model_json = model.to_json()\n",
    "        with open(\"model_pol.json\", \"w\") as json_file:\n",
    "            json_file.write(model_json)\n",
    "\n",
    "        # serialize weights to HDF5\n",
    "        # model.save_weights(\"model_pick_best.hdf5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction code \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    # load json and create model\n",
    "    json_file = open('model_pol.json', 'r')\n",
    "    loaded_model_json = json_file.read()\n",
    "    json_file.close()\n",
    "    model = model_from_json(loaded_model_json, custom_objects={\n",
    "        'tf': tf})\n",
    "    # load weights into new model\n",
    "    model.load_weights(\"model_pol_best.hdf5\")\n",
    "    #model = make_parallel(model, 1)\n",
    "    hist = joblib.load(\"hist_pol.pkl\")\n",
    "    print(\"Loaded model from disk\")\n",
    "\n",
    "    if True:\n",
    "        plt.clf()\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "        ax.plot(np.arange(len(hist['loss'])), hist['loss'], c='k', marker='o',\n",
    "                label='Training')\n",
    "        ax.plot(np.arange(len(hist['val_loss'])), hist['val_loss'], c='b',\n",
    "                marker='o', label='Validation')\n",
    "        ax.set_xlabel(\"Epoch\")\n",
    "        ax.set_ylabel(\"Loss\")\n",
    "        ax.legend(loc=0)\n",
    "        \"\"\"\"\n",
    "        if show_plot:\n",
    "            plt.show()\n",
    "        else:\n",
    "            plt.savefig(\"pick_loss.pdf\")\n",
    "         \"\"\"\n",
    "        plt.savefig(\"pick_loss.pdf\")\n",
    "        plt.show()\n",
    "\n",
    "    print(\"Running forward predictions on test data\")\n",
    "    y_pred_proba, y_test = predict_test(model, seis_data_fn=hdf_file,\n",
    "                                        # y_pred, y_test, snr, dist, mag = predict_test(model, seis_data_fn=hdf_file,\n",
    "                                        validation_ratio=0.8,\n",
    "                                        batch_size=10000,\n",
    "                                        sample_count=sample_count)\n",
    "    y_pred = np.argmax(y_pred_proba, axis=1)\n",
    "    print(\"\")\n",
    "    print(\"Finished predictions\")\n",
    "\n",
    "    from sklearn.metrics import classification_report\n",
    "    target_names = ['Noise', 'P','S','Lg']\n",
    "    print(classification_report(y_test, y_pred, labels=[0,1,2,3]))\n",
    "    y_test = y_test.flatten()\n",
    "    FP = []\n",
    "    FN = []\n",
    "    for pr_min in np.arange(0, 1.0, 0.01):\n",
    "        y_pred0 = np.argmax(y_pred_proba, axis=1)\n",
    "        y_prob0 = np.max(y_pred_proba, axis=1)\n",
    "                       \n",
    "            \n",
    "\n",
    "        idx = np.where(y_prob0 < pr_min)[0]\n",
    "        y_pred0[idx] = 2\n",
    "\n",
    "        idx = np.where(y_pred0 < 2)[0]\n",
    "        precision = np.sum(y_pred0[idx] == y_test[idx]) / float(idx.size)\n",
    "        idx = np.where(y_test < 2)[0]\n",
    "        recall = np.sum(y_pred0[idx] == y_test[idx]) / float(idx.size)\n",
    "        FP.append(precision)\n",
    "        FN.append(recall)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "    pr_min = np.arange(0, 1.0, 0.01)\n",
    "    plt.scatter(FN, FP, c=pr_min, cmap='jet')\n",
    "    plt.xlabel(\"Recall\")\n",
    "    plt.ylabel(\"Precision\")\n",
    "    plt.colorbar()\n",
    "    \n",
    "    \n",
    "    plt.savefig(\"recall_precision.pdf\")\n",
    "    plt.xlim(0.1, 1.0)\n",
    "    plt.ylim(0.1, 1.0)\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:deep-learning] *",
   "language": "python",
   "name": "conda-env-deep-learning-py"
  },
  "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.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "205px",
    "width": "218px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "196px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 560.61666,
   "position": {
    "height": "582.85px",
    "left": "1653px",
    "right": "78.4667px",
    "top": "72px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "block",
   "window_display": true
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}