text: from typing import Tuple
  import pandas as pd
  import numpy as np
  import os
  import pickle
  import obspy
  from subprocess import call
  import json
  import datetime
  from scipy.spatial import distance_matrix
  import matplotlib.pyplot as plt
  
  class PhaseNet_Analysis (object):
  
      '''
      This class analysis and compare the 'P' picks and 'S' Picks
      with the given catalog.
      '''
  
      def __init__(self,phasenet_direc: 'str', chile_GFZ_online_direc:'str', export_DF_path:'str',
                  export_mseed_path:'str', working_direc:'str', picks_name:'str',
                  start_year_analysis:'int', start_day_analysis:'int', 
                  end_year_analysis:'int', end_day_analysis:'int', analysis:'bool', time_lag_threshold:'int',
                  station_name_list:'str', apply_filter: 'bool', freqmin:'float', freqmax:'float'):
  
          '''
          Parameters initialization:
              - phasenet_direc: The path of PhaseNet package
  
              - chile_GFZ_online_direc: The directory of stored all three components mseed files.
  
              - export_DF_path : The path of all mseed files DataFrame.
                                  This directory must contain 'DF_chile_path_file.pkl'.
                                  This DataFrame is created once using "generate_DF_file_path" function.
                                  This function is not used often.
                                  After run "DF_path" function, another two DataFrames
                                  ("DF_selected_chile_path_file.pkl" and "DF_auxiliary_path_file.pkl") 
                                  will be created and stored in this directory as a feed path to run PhaseNet.
  
                                  It should be noted that after running "events_data_frame", "events.pkl" will be created
                                  and stored in this directory.
  
              - export_mseed_path: The folder path of writing three components mseed.
                  This folder will be used by PhaseNet.
  
              - working_direc: the path of working directory
  
              - picks_name: the name of picks Dataframe existed in 'export_DF_path ' directory
  
              - start_year_analysis: start year of analysis with 4 digit (like 2011)
  
              - start_day_analysis: start day of analysis related to 
                                  the start_year_analysis. This variable should between 1 to 365. 
  
              - end_year_analysis: end year of analysis with 4 digit (like 2011)
  
              - end_day_analysis: end day of analysis related to 
                                  the end_year_analysis. This variable should between 1 to 365. 
  
              - analysis (boolean): Because running PhaseNet is time consumming, analysis variable helps to use existing 
                  "PhaseNet_result_s_picks.pkl", "PhaseNet_result_s_picks.pkl" which are generated by Phasenet without running PhaseNet again.
                  
                  If you want to run PhaseNet for the first time in the given interval, this variable should be set to True
                  If you already have "PhaseNet_result_s_picks.pkl", "PhaseNet_result_s_picks.pkl" for the given interval, this variable should be set to False.
                                  
                                  - True: start running PhaseNet and perform visualization.
                                  - False:  Use existing PhaseNet results and perform visualization.
  
              
              - time_lag_threshold (int) : Proper time lag threshold in millisecond.
                                          This variable has been use to perform qaulity control of PhaseNet and existing catalog.
  
              - station_name_list: the name of selected station stored in the text file
  
              - apply_filter (boolean): 
                                  - True: Apply filtering on raw data before feeding to phasenet
                                  - False:  feeding raw data without filtering.
              
              - freqmin (float): lower bound filtering (if apply_filter ==True)
  
              - freqmax (float): upper bound filtering (if apply_filter ==True)
  
          '''
          os.chdir('{0}'.format(phasenet_direc))
          self.PROJECT_ROOT = os.getcwd()
  
          self.chile_GFZ_online_direc = chile_GFZ_online_direc
          self.export_DF_path = export_DF_path
          self.export_mseed_path = export_mseed_path
          self.working_direc = working_direc
          self.picks_name = picks_name
          self.start_year_analysis = start_year_analysis
          self.start_day_analysis = start_day_analysis
          self.end_year_analysis = end_year_analysis
          self.end_day_analysis = end_day_analysis
          self.analysis = analysis
          self.time_lag_threshold = time_lag_threshold
          self.station_name_list = station_name_list
          self.apply_filter = apply_filter
          self.freqmin = freqmin
          self.freqmax = freqmax
  
      
      def __call__ (self):
  
  
          if self.analysis == True:
  
              self.DF_path ()
  
              # load DF_auxiliary_path_file.pkl
              with open(os.path.join(self.export_DF_path, "DF_auxiliary_path_file.pkl"),'rb') as fp:
                  DF_auxiliary_path_file = pickle.load(fp)
  
              # load DF_selected_chile_path_file.pkl
              with open(os.path.join(self.export_DF_path, "DF_selected_chile_path_file.pkl"),'rb') as fp:
                  DF_selected_chile_path_file = pickle.load(fp)
              
  
              # Creat S picks DataFrame
              df_S_picks = pd.DataFrame(index =[])
  
              # Creat P picks DataFrame
              df_P_picks = pd.DataFrame(index =[])
  
  
              # Iterate over all mseed files which contains 3 components .
              
              for i in range (DF_auxiliary_path_file.shape[0]):
  
                  # Write mseed file to mseed folder path
                  mseed_name = self.three_components_mseed_maker (i,DF_auxiliary_path_file, 
                                          DF_selected_chile_path_file)
                  
                  # Write the name of mseed file in mseed.csv
                  self.write_mseed_names(mseed_name)
  
                  # Run PhaseNet
                  self.run_phasenet()
  
                  # Remove created mseed in mseed folder to free up memory
                  self.remove_mseed (mseed_name)
  
                  # Read the output of PhaseNet and store P picks and S picks in two data frames
                  p_picks, s_picks = self.read_picks()
  
                  # check
                  if ((p_picks.empty == True) and (s_picks.empty == True)):
  
                      continue
  
                  # save data in data frame
                  #df_total = self.save_DF (df_p_picks, df_s_picks, mseed_name, df_picks)
                  #df_picks = df_total
                  df_S_picks = pd.concat([df_S_picks, s_picks], axis=0)
                  df_P_picks = pd.concat([df_P_picks, p_picks], axis=0)
  
              # save PhaseNet p picks result file in the export directory as 'PhaseNet_result_p_picks.pkl'
              df_P_picks.to_pickle(os.path.join(self.export_DF_path, 'PhaseNet_result_p_picks.pkl'))
  
              # save PhaseNet s picks result file in the export directory as 'PhaseNet_result_s_picks.pkl'
              df_S_picks.to_pickle(os.path.join(self.export_DF_path, 'PhaseNet_result_s_picks.pkl'))
  
              # Apply filter on catalog data between start time and end time of analysis
              catalog_DF_P_picks, catalog_DF_S_picks = self.filter_picks_DF ()
  
              # save catalog p picks result file in the export directory as 'catalog_p_picks.pkl'
              catalog_DF_P_picks.to_pickle(os.path.join(self.export_DF_path, 'catalog_p_picks.pkl'))
  
              # save PhaseNet s picks result file in the export directory as 'catalog_s_picks.pkl'
              catalog_DF_S_picks.to_pickle(os.path.join(self.export_DF_path, 'catalog_s_picks.pkl'))
  
              # Perform visulization & qaulity control of P picks
              # The results of visulization will be found at self.export_DF_path directory
              self.compare_PhaseNet_catalog_P_picks()
  
              # Perform visulization & qaulity control of S picks
              # The results of visulization will be found at self.export_DF_path directory
              self.compare_PhaseNet_catalog_S_picks()
          
  
          # Perform visulization without running PhaseNet
          else:
  
              # Perform visulization & qaulity control of P picks
              # The results of visulization will be found at self.export_DF_path directory
              #self.compare_PhaseNet_catalog_P_picks()
  
              # Perform visulization & qaulity control of S picks
              # The results of visulization will be found at self.export_DF_path directory
  
              self.compare_PhaseNet_catalog_S_picks()      
     
      
      def DF_path (self):
  
          '''
          This function loads "DF_chile_path_file.pkl" file (path of all mseed files
          from export_DF_path and filter the data between the given interval.
          After using this function, 'DF_selected_chile_path_file.pkl' and 
          'DF_auxiliary_path_file.pkl' will be created to feed data to PhaseNet.
  
          Important note: This function will not consider mseed file with less than 3- components.
          '''
  
          # Read pickle data (Path of all chile stream data)
          with open(os.path.join(self.export_DF_path, "DF_chile_path_file.pkl"),'rb') as fp:
              chile_path_file = pickle.load(fp)
  
          chile_path_file = chile_path_file[(chile_path_file['year']>= self.start_year_analysis) & 
                      (chile_path_file['year']<= self.end_year_analysis)]
  
          chile_path_file['convert_yeartoday']= 365*chile_path_file['year']+chile_path_file['day']
          
          # creat upper and lower limit to filter
          lower_limit = 365*self.start_year_analysis + self.start_day_analysis 
          upper_limit = 365*self.end_year_analysis   + self.end_day_analysis
  
          # Apply filter
          chile_path_file = chile_path_file[(chile_path_file['convert_yeartoday']>= lower_limit) & 
                  (chile_path_file['convert_yeartoday']<= upper_limit)]   
  
          chile_path_file = chile_path_file.drop_duplicates()
  
          # creat new DataFrame to make sure all 3-components are existed
          df_counter = chile_path_file.groupby(['network','station', 'year', 'day']).size().reset_index(name='count')
          df_counter = df_counter[df_counter['count']==3]
  
          # drop the 'count' column
          df_counter = df_counter.drop(columns=['count'])
  
          # Save selected DataFrame based on given time interval
          chile_path_file.to_pickle(os.path.join(self.export_DF_path , 'DF_selected_chile_path_file.pkl'))
  
          # Save auxiliary DataFrame based on given time interval
          df_counter.to_pickle(os.path.join(self.export_DF_path , 'DF_auxiliary_path_file.pkl'))
  
      
      def three_components_mseed_maker (self,i:'int',DF_auxiliary_path_file:'pd.DataFrame', 
                                      DF_selected_chile_path_file: 'pd.DataFrame'):
          '''
          This function write a single three components mseed file in mseed folder.
          Parameters:
                      - i (int): the ith row of "DF_auxiliary_path_file" data frame
                      - DF_auxiliary_path_file (data frame): A data frame containing mseed file name 
                          with existed three componets.
                      - DF_selected_chile_path_file(data frame): A selected data frame based on the given time interval.
                      
          '''
  
          # Apply filter to determine the three components in "DF_selected_chile_path_file"
          df = DF_selected_chile_path_file[['network', 'station', 'year','day']]==DF_auxiliary_path_file.iloc[i]
          df_components = DF_selected_chile_path_file[(df['network']== True) & 
                          (df['station']== True) & (df['year']== True) &
                          (df['day']== True)]
  
          # Read three components mseed file 
          streamZ = obspy.read(df_components.path.iloc[2])
          streamN = obspy.read(df_components.path.iloc[1])
          streamE = obspy.read(df_components.path.iloc[0])
          streamN += streamE
          streamN += streamZ
          stream = streamN.sort()
  
          # Apply filter if apply_filter ==True
          if self.apply_filter == True:
              stream = self.filter (stream)
  
          # Write the mseed file (three components) in mseed folder
          string = df_components.file_name.iloc[2]
          mseed_name = string.replace("HHE", "").replace("HHN", "").replace("HHZ", "")
          stream.write(os.path.join(self.export_mseed_path, mseed_name), sep="\t", format="MSEED")
  
          return mseed_name
  
      def write_mseed_names(self, mseed_name:'str') -> None:
          '''
          This function write the name of mseed file to mseed.csv.
          This mseed.csv will be used by PhaseNet
          '''
          df = pd.DataFrame ([mseed_name], columns = ['fname'])
  
          df.to_csv((os.path.join(self.working_direc, 'mseed.csv')),index=False)
      
      def run_phasenet (self) -> None:
          '''
          Run the predefined PhaseNet model (model/190703-214543) to pick S and P .
          '''
          #cmd = '/home/javak/miniconda3/envs/phasenet/bin/python phasenet/predict.py --model=model/190703-214543 --data_list=/home/javak/Sample_data_chile/mseed.csv --data_dir=/home/javak/Sample_data_chile/mseed --format=mseed --plot_figure'.split()
          cmd = '/home/javak/miniconda3/envs/phasenet/bin/python phasenet/predict.py --model=model/190703-214543 --data_list=/home/javak/Sample_data_chile/mseed.csv --data_dir=/home/javak/Sample_data_chile/mseed --format=mseed'.split()
          
          call(cmd)
      
      def read_picks (self):
  
          '''
          Read the csv file of PhaseNet output and return the P picks and S picks.
  
              Output:
                      - df_p_picks (dataframe): Phasenet P picks dataframe
                      - df_s_picks (dataframe): PhaseNet S picks dataframe
              
          '''
          picks_csv = pd.read_csv(os.path.join(self.PROJECT_ROOT, "results/picks.csv"), sep="\t")
          picks_csv.loc[:, 'p_idx'] = picks_csv["p_idx"].apply(lambda x: x.strip("[]").split(","))
          picks_csv.loc[:, 'p_prob'] = picks_csv["p_prob"].apply(lambda x: x.strip("[]").split(","))
          picks_csv.loc[:, 's_idx'] = picks_csv["s_idx"].apply(lambda x: x.strip("[]").split(","))
          picks_csv.loc[:, 's_prob'] = picks_csv["s_prob"].apply(lambda x: x.strip("[]").split(","))
  
          with open(os.path.join(self.PROJECT_ROOT, "results/picks.json")) as fp:
              picks_json = json.load(fp)
          
          if len (picks_json) == 0:
              df_p_picks = pd.DataFrame({'A' : []})
              df_s_picks = pd.DataFrame({'B' : []})
          else:
              df = pd.DataFrame.from_dict(pd.json_normalize(picks_json), orient='columns')
              df_p_picks = df[df["type"] == 'p']
              df_s_picks = df[df["type"] == 's']
  
          return df_p_picks, df_s_picks
      
      def remove_mseed (self,mseed_name:'str') -> pd.DataFrame:
  
          '''
          This function removes created mseed in mseed folder to free up memory.
          '''
  
          file_path = os.path.join(self.export_mseed_path, mseed_name)
          os.remove(file_path)
      
      def save_DF (self, df_p_waves, df_s_waves, daily_data, df):
  
          data = {'P_waves':[df_p_waves],
                  'S_waves':[df_s_waves]
                  }
  
          # Creates pandas DataFrame.
          df_new = pd.DataFrame(data, index =[daily_data])
          df_total = df.append(df_new)
          return df_total
      
      def filter_picks_DF (self)-> tuple:
  
          '''
          This function apply filter on existing picks DataFrame according to 
          start_year_analysis, end_year_analysis, start_day_analysis, end_day_analysis.
  
              Output:
                      - catalog_DF_P_picks (dataframe): catalog P picks dataframe
                      - catalog_DF_S_picks (dataframe): catalog S picks dataframe
              '''
          # load picks_2007_2020.pkl
          with open(os.path.join(self.export_DF_path, self.picks_name),'rb') as fp:
              DF_picks = pickle.load(fp)
  
          # convert the Day of Year in Python to Month/Day
          start_date = datetime.datetime.strptime('{} {}'.format(start_day_analysis, start_year_analysis),'%j %Y')
          end_date   = datetime.datetime.strptime('{} {}'.format(end_day_analysis, end_year_analysis),'%j %Y')
  
          start_date_obspy = obspy.UTCDateTime(year=start_year_analysis, month=start_date.month, day=start_date.day, strict=False)
          end_date_obspy = obspy.UTCDateTime(year=end_year_analysis, month=end_date.month, day=end_date.day, hour=24, strict=False)
  
   
          catalog_DF_P_picks = DF_picks[(DF_picks['picks_time']>= start_date_obspy) & (DF_picks['picks_time']<=end_date_obspy) & (DF_picks['phase_hint']=='P')]
          catalog_DF_S_picks = DF_picks[(DF_picks['picks_time']>= start_date_obspy) & (DF_picks['picks_time']<=end_date_obspy) & (DF_picks['phase_hint']=='S')]
  
          return catalog_DF_P_picks, catalog_DF_S_picks
  
  
      def compare_PhaseNet_catalog_P_picks (self):
  
          '''
          This function compares the result of PhaseNet and existing catalog.
          '''
          
          # catalog_DF_P_picks, df_P_picks
          with open(os.path.join(self.export_DF_path, "PhaseNet_result_p_picks.pkl"),'rb') as fp:
              df_P_picks = pickle.load(fp)
  
          with open(os.path.join(self.export_DF_path, "catalog_p_picks.pkl"),'rb') as fp:
              catalog_DF_P_picks = pickle.load(fp)
  
          # creat extra columns
          df_P_picks[['network', 'others']] = df_P_picks['id'].str.split('.', 1, expand=True)
          df_P_picks[['station_code', 'date']] = df_P_picks['others'].str.split('.', 1, expand=True)
          df_P_picks = df_P_picks.drop(['date', 'others'], axis=1)
  
          # find common station_code in catalog and PhaseNet
          boolean_column = catalog_DF_P_picks['station_code'].isin(df_P_picks['station_code'])
          catalog_DF_P_picks = catalog_DF_P_picks[(boolean_column==True)]
          all_dists = np.array([])
          common_stations = catalog_DF_P_picks['station_code'].unique()
  
          # Creat an empty DataFrame file to store all UTC time of PhaseNet in common station
          all_p_picks_exist_in_catalogtory = pd.DataFrame(index =[])
  
          # loop over all common station
          for i in common_stations:
              bo = catalog_DF_P_picks['station_code']==i
              catalog_filter_station = catalog_DF_P_picks[(bo==True)]
              ao = df_P_picks['station_code']==i
              phasenet_filter_station = df_P_picks[(ao==True)]
  
              # Convert UTC time to datetime64[ms] (millisecond)
              a = catalog_filter_station.picks_time.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
              b = phasenet_filter_station.timestamp.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
  
              # Calculate P1 norme of all datetime64[ms]
              dist_mat = distance_matrix(a,b, p=1)
              dists = np.min(dist_mat, axis=1)
              all_dists = np.append(all_dists, dists)
  
              # append phasenet_filter_station
              min_index = np.argmin(dist_mat, axis=1)
  
              #phasenet_filter_station = phasenet_filter_station[min_index]
              phasenet_filter_station = phasenet_filter_station.iloc[min_index,:]
  
              all_p_picks_exist_in_catalogtory = pd.concat([all_p_picks_exist_in_catalogtory, phasenet_filter_station], axis=0)
  
          # Filter the time lag with the given threshold and capture the picks with more than a given time lag
          dists_filter_lag_time_m=all_dists[all_dists >= self.time_lag_threshold]
  
          # Perform P picks Quality control of PhaseNet by using existing P picks catalog with more than a given time lag
          fig_lag_m, ax_lag_m = plt.subplots(figsize=(20,10))
  
          label_more = '{0}{1}{2}'.format('The total number of common P picks with more than 2 seconds (', dists_filter_lag_time_m.shape[0], ')')
  
          n_lag_m, bins_lag_m, patches_lag_m = ax_lag_m.hist(dists_filter_lag_time_m, 20, density=False, facecolor='b', alpha=0.75, label=label_more)
          steps = (max(dists_filter_lag_time_m) - min(dists_filter_lag_time_m))/20
          plt.xticks(np.arange(min(dists_filter_lag_time_m), max(dists_filter_lag_time_m), step=steps))
          plt.xlabel('Time residual between catalog and PhaseNet (ms)', fontsize=20)
          plt.ylabel('Number of P picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('Common P picks distribution with more than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_lag_m), max(bins_lag_m))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('P picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with more than 2s absolute time residual', extention='png')
          fig_lag_m.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          # Filter the time lag with the given threshold and capture the picks with less than a given time lag
          dists_filter_lag_time=all_dists[all_dists < self.time_lag_threshold]
  
          # Perform P picks Quality control of PhaseNet by using existing P picks catalog with less than a given time lag
          fig_lag, ax_lag = plt.subplots(figsize=(20,10))
          label_less = '{0}{1}{2}{3}{4}'.format('The total number of common P picks ', dists_filter_lag_time.shape[0], ' out of all catalog P picks (', all_dists.shape[0],')')
  
          n_lag, bins_lag, patches_lag = ax_lag.hist(dists_filter_lag_time, 20, density=False, facecolor='b', alpha=0.75, label=label_less)
          steps = (max(dists_filter_lag_time) - min(dists_filter_lag_time))/20
          plt.xticks(np.arange(min(dists_filter_lag_time), max(dists_filter_lag_time), step=steps))
          plt.xlabel('Time residual between catalog and PhaseNet (ms)', fontsize=20)
          plt.ylabel('Number of P picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('Common P picks distribution with less than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_lag), max(bins_lag))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('P picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with less than 2s time residual', extention='png')
          fig_lag.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          
          # Plot the histogram probability of all P picks PhaseNet in the common stations 
          
          all_p_picks_exist_in_catalogtory=all_p_picks_exist_in_catalogtory.iloc[all_dists < self.time_lag_threshold]
          
          fig0, ax0 = plt.subplots(figsize=(20,10))
          label_co = '{0}{1}{2}{3}{4}'.format('The total number of common P picks ',all_p_picks_exist_in_catalogtory.shape[0], ' out of all catalog P picks (', all_dists.shape[0],')')
          n_co, bins_co, patches_co = ax0.hist(all_p_picks_exist_in_catalogtory.prob, 21, density=False, alpha=0.75, color = "blue", label=label_co)
          #label_all = '{0}{1}{2}'.format('All PhaseNet P picks (',df_P_picks.shape[0], ' P picks)')
          #ax0.hist(df_P_picks.prob, 21, density=False, color = "skyblue", ec="skyblue", alpha=0.75, label= label_all)
          steps = (max(all_p_picks_exist_in_catalogtory.prob) - min(all_p_picks_exist_in_catalogtory.prob))/21
          #plt.xticks(np.arange(min(all_p_picks_exist_in_catalogtory.prob), max(all_p_picks_exist_in_catalogtory.prob), step=steps))
          plt.xticks(np.around(np.arange(min(all_p_picks_exist_in_catalogtory.prob), max(all_p_picks_exist_in_catalogtory.prob), step=steps), decimals=2))
          plt.xlabel('PhaseNet P picks Probability', fontsize=20)
          plt.ylabel('Number of P picks', fontsize=20)
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
          plt.title('Common P picks with less than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_co), max(bins_co))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('Common P picks in catalog & PhaseNet \n (2012-01-01 to 2012-01-31)', extention='png')
          fig0.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
  
  
  
  
  
          fig11, ax11 = plt.subplots(figsize=(20,10))
          label_pha = '{0}{1}{2}'.format('The total number of PhaseNet P picks (',df_P_picks.shape[0], ')')
          #n_co, bins_co, patches_co = ax0.hist(all_p_picks_exist_in_catalogtory.prob, 21, density=False, alpha=0.75, label=label_co)
          #label_all = '{0}{1}{2}'.format('All PhaseNet P picks (',df_P_picks.shape[0], ' P picks)')
          ax11.hist(df_P_picks.prob, 21, density=False, color = "slateblue", alpha=0.75, label =label_pha)
          steps = (max(df_P_picks.prob) - min(df_P_picks.prob))/21
          plt.xticks(np.around(np.arange(min(df_P_picks.prob), max(df_P_picks.prob), step=steps), decimals=2))
  
          plt.xlabel('PhaseNet P picks Probability', fontsize=20)
          plt.ylabel('Number of P picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('PhaseNet output P picks(2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_co), max(bins_co))
          plt.grid(True)
          plt.legend(loc='upper center',fontsize=15)
          file_name = '{0}.{extention}'.format('PhaseNet output P picks(2012-01-01 to 2012-01-31)', extention='png')
          fig11.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
          
          # Perform P picks Quality control of PhaseNet by using existing P picks catalog with defined time residual 
          lag_time = all_dists[all_dists < self.time_lag_threshold]
          for j in range (1, bins_lag.shape[0]):
              select_p_picks=all_p_picks_exist_in_catalogtory.iloc[(lag_time < bins_lag[j]) & (lag_time >= bins_lag[j-1])]
              fig, ax = plt.subplots(figsize=(20,10))
              ax.hist(select_p_picks.prob, 21, density=False, color = "b", ec="b", alpha=0.75)
  
              steps = (1 - 0.3)/21
              plt.xticks(np.around(np.arange(min(select_p_picks.prob), max(select_p_picks.prob), step=steps), decimals=2))
  
              plt.tick_params(axis='x', labelsize=15)
              plt.tick_params(axis='y', labelsize=15)
              plt.xlabel('PhaseNet P picks Probability', fontsize=20)
              plt.ylabel('Number of P picks', fontsize=20)
              title_name = '{0}{1}{2}{3}{4}'.format('Common P picks with ',round (bins_lag[j-1]),' - ', round (bins_lag[j]), ' time residual (ms)')
              plt.title(title_name, fontsize=24, pad=23)
              file_name = '{0}{1}.{extention}'.format('PhaseNet_result_P_bins: ',round (bins_lag[j]), extention='png')
              fig.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
  
      def compare_PhaseNet_catalog_S_picks (self):
  
          '''
          This function compares the result of PhaseNet and existing catalog.
          '''
  
          # catalog_DF_P_picks, df_P_picks
          with open(os.path.join(self.export_DF_path, "PhaseNet_result_s_picks.pkl"),'rb') as fp:
              df_S_picks = pickle.load(fp)
  
          with open(os.path.join(self.export_DF_path, "catalog_s_picks.pkl"),'rb') as fp:
              catalog_DF_S_picks = pickle.load(fp)
  
          # creat extra columns
          df_S_picks[['network', 'others']] = df_S_picks['id'].str.split('.', 1, expand=True)
          df_S_picks[['station_code', 'date']] = df_S_picks['others'].str.split('.', 1, expand=True)
          df_S_picks = df_S_picks.drop(['date', 'others'], axis=1)
  
          # find common station_code in catalog_DF_P_picks
          boolean_column = catalog_DF_S_picks['station_code'].isin(df_S_picks['station_code'])
          catalog_DF_S_picks = catalog_DF_S_picks[(boolean_column==True)]
          all_dists = np.array([])
          common_stations = catalog_DF_S_picks['station_code'].unique()
  
          # Creat an empty DataFrame file to store all UTC time of PhaseNet in common station
          all_s_picks_exist_in_catalogtory = pd.DataFrame(index =[])
          
          # loop over all common statio
          for i in common_stations:
              bo = catalog_DF_S_picks['station_code']==i
              catalog_filter_station = catalog_DF_S_picks[(bo==True)]
              ao = df_S_picks['station_code']==i
              phasenet_filter_station = df_S_picks[(ao==True)]
  
              # Convert UTC time to datetime64[ms] (millisecond)
              a = catalog_filter_station.picks_time.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
              b = phasenet_filter_station.timestamp.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
  
              # Calculate P1 norme of all datetime64[m
              dist_mat = distance_matrix(a,b, p=1)
              dists = np.min(dist_mat, axis=1)
              all_dists = np.append(all_dists, dists)
  
              # append phasenet_filter_station
              min_index = np.argmin(dist_mat, axis=1)
  
              #phasenet_filter_station = phasenet_filter_station[min_index]
              phasenet_filter_station = phasenet_filter_station.iloc[min_index,:]
  
              all_s_picks_exist_in_catalogtory = pd.concat([all_s_picks_exist_in_catalogtory, phasenet_filter_station], axis=0)
      
          # Filter the time lag with the given threshold and capture the picks with more than 2 second time lag
          dists_filter_lag_time_m=all_dists[all_dists >= self.time_lag_threshold]
  
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with more than a given time lag
          fig_lag_m, ax_lag_m = plt.subplots(figsize=(20,10))
  
          label_more = '{0}{1}{2}'.format('The total number of common S picks with more than 2 seconds (', dists_filter_lag_time_m.shape[0], ')')
  
          n_lag_m, bins_lag_m, patches_lag_m = ax_lag_m.hist(dists_filter_lag_time_m, 20, density=False, facecolor='r', alpha=0.75, label=label_more)
          steps = (max(dists_filter_lag_time_m) - min(dists_filter_lag_time_m))/20
          plt.xticks(np.arange(min(dists_filter_lag_time_m), max(dists_filter_lag_time_m), step=steps))
          plt.xlabel('Time residual between catalog and PhaseNet (ms)', fontsize=20)
          plt.ylabel('Number of S picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('Common S picks distribution with more than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_lag_m), max(bins_lag_m))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with more than 2s absolute time residual', extention='png')
          fig_lag_m.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          # Filter the time lag with the given threshold and capture the picks with less than a given time lag
          dists_filter_lag_time=all_dists[all_dists < self.time_lag_threshold]
  
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with less than a given time lag
          fig_lag, ax_lag = plt.subplots(figsize=(20,10))
          label_less = '{0}{1}{2}{3}{4}'.format('The total number of common S picks ', dists_filter_lag_time.shape[0], ' out of all catalog S picks (', all_dists.shape[0],')')
  
          n_lag, bins_lag, patches_lag = ax_lag.hist(dists_filter_lag_time, 20, density=False, facecolor='r', alpha=0.75, label=label_less)
          steps = (max(dists_filter_lag_time) - min(dists_filter_lag_time))/20
          plt.xticks(np.arange(min(dists_filter_lag_time), max(dists_filter_lag_time), step=steps))
          plt.xlabel('Time residual between catalog and PhaseNet (ms)', fontsize=20)
          plt.ylabel('Number of S picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('Common S picks distribution with less than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_lag), max(bins_lag))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with less than 2s time residual', extention='png')
          fig_lag.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          
          # Plot the histogram probability of all S picks PhaseNet in the common stations 
          
          all_s_picks_exist_in_catalogtory=all_s_picks_exist_in_catalogtory.iloc[all_dists < self.time_lag_threshold]
          
          fig0, ax0 = plt.subplots(figsize=(20,10))
          label_co = '{0}{1}{2}{3}{4}'.format('The total number of common S picks ',all_s_picks_exist_in_catalogtory.shape[0], ' out of all catalog S picks (', all_dists.shape[0],')')
          n_co, bins_co, patches_co = ax0.hist(all_s_picks_exist_in_catalogtory.prob, 21, density=False, alpha=0.75, color = "r", label=label_co)
          #label_all = '{0}{1}{2}'.format('All PhaseNet P picks (',df_P_picks.shape[0], ' P picks)')
          #ax0.hist(df_P_picks.prob, 21, density=False, color = "skyblue", ec="skyblue", alpha=0.75, label= label_all)
          steps = (max(all_s_picks_exist_in_catalogtory.prob) - min(all_s_picks_exist_in_catalogtory.prob))/21
          #plt.xticks(np.arange(min(all_p_picks_exist_in_catalogtory.prob), max(all_p_picks_exist_in_catalogtory.prob), step=steps))
          plt.xticks(np.around(np.arange(min(all_s_picks_exist_in_catalogtory.prob), max(all_s_picks_exist_in_catalogtory.prob), step=steps), decimals=2))
          plt.xlabel('PhaseNet S picks Probability', fontsize=20)
          plt.ylabel('Number of S picks', fontsize=20)
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
          plt.title('Common S picks with less than 2 seconds time residual (2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_co), max(bins_co))
          plt.grid(True)
          plt.legend(loc='upper center', fontsize=15)
          file_name = '{0}.{extention}'.format('Common S picks in catalog & PhaseNet \n (2012-01-01 to 2012-01-31)', extention='png')
          fig0.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
  
  
  
  
  
          fig11, ax11 = plt.subplots(figsize=(20,10))
          label_pha = '{0}{1}{2}'.format('The total number of PhaseNet S picks (',df_S_picks.shape[0], ')')
          #n_co, bins_co, patches_co = ax0.hist(all_p_picks_exist_in_catalogtory.prob, 21, density=False, alpha=0.75, label=label_co)
          #label_all = '{0}{1}{2}'.format('All PhaseNet P picks (',df_P_picks.shape[0], ' P picks)')
          ax11.hist(df_S_picks.prob, 21, density=False, color = "deeppink", alpha=0.75, label =label_pha)
          steps = (max(df_S_picks.prob) - min(df_S_picks.prob))/21
          plt.xticks(np.around(np.arange(min(df_S_picks.prob), max(df_S_picks.prob), step=steps), decimals=2))
  
          plt.xlabel('PhaseNet S picks Probability', fontsize=20)
          plt.ylabel('Number of S picks', fontsize=20)
  
          plt.tick_params(axis='x', labelsize=15)
          plt.tick_params(axis='y', labelsize=15)
  
          plt.title('PhaseNet output S picks(2012-01-01 to 2012-01-31)', fontsize=24, pad=23)
          plt.xlim(min(bins_co), max(bins_co))
          plt.grid(True)
          plt.legend(loc='upper center',fontsize=15)
          file_name = '{0}.{extention}'.format('PhaseNet output S picks(2012-01-01 to 2012-01-31)', extention='png')
          fig11.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
          
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with defined time residual 
          lag_time = all_dists[all_dists < self.time_lag_threshold]
          for j in range (1, bins_lag.shape[0]):
              select_s_picks=all_s_picks_exist_in_catalogtory.iloc[(lag_time < bins_lag[j]) & (lag_time >= bins_lag[j-1])]
              fig, ax = plt.subplots(figsize=(20,10))
              ax.hist(select_s_picks.prob, 21, density=False, color = "r", ec="r", alpha=0.75)
  
              steps = (1 - 0.3)/21
              plt.xticks(np.around(np.arange(min(select_s_picks.prob), max(select_s_picks.prob), step=steps), decimals=2))
  
              plt.tick_params(axis='x', labelsize=15)
              plt.tick_params(axis='y', labelsize=15)
              plt.xlabel('PhaseNet S picks Probability', fontsize=20)
              plt.ylabel('Number of S picks', fontsize=20)
              title_name = '{0}{1}{2}{3}{4}'.format('Common S picks with ',round (bins_lag[j-1]),' - ', round (bins_lag[j]), ' time residual (ms)')
              plt.title(title_name, fontsize=24, pad=23)
              file_name = '{0}{1}.{extention}'.format('PhaseNet_result_S_bins: ',round (bins_lag[j]), extention='png')
              fig.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
  
  
          '''
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with more than 2 second time lag
          fig_lag_m, ax_lag_m = plt.subplots(figsize=(20,10))
          n_lag_m, bins_lag_m, patches_lag_m = ax_lag_m.hist(dists_filter_lag_time_m, 20, density=False, facecolor='r', alpha=0.75)
          steps = (max(dists_filter_lag_time_m) - min(dists_filter_lag_time_m))/20
          plt.xticks(np.arange(min(dists_filter_lag_time_m), max(dists_filter_lag_time_m), step=steps))
          plt.xlabel('Time lag between catalog and PhaseNet (ms)', fontsize=18)
          plt.ylabel('Frequency', fontsize=18)
          plt.title('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with more 2s time lag', fontsize=21)
          plt.xlim(min(bins_lag_m), max(bins_lag_m))
          plt.grid(True)
          file_name = '{0}.{extention}'.format('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31) with more 2s time lag', extention='png')
          fig_lag_m.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
  
  
          # Filter the time lag with the given threshold and capture the picks with less than 2 second time lag
          dists_filter_lag_time=all_dists[all_dists < self.time_lag_threshold]
  
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with less than 2 second time lag
          fig_lag, ax_lag = plt.subplots(figsize=(20,10))
          n_lag, bins_lag, patches_lag = ax_lag.hist(dists_filter_lag_time, 20, density=False, facecolor='r', alpha=0.75)
          steps = (max(dists_filter_lag_time) - min(dists_filter_lag_time))/20
          plt.xticks(np.arange(min(dists_filter_lag_time), max(dists_filter_lag_time), step=steps))
          plt.xlabel('Time lag between catalog and PhaseNet (ms)', fontsize=18)
          plt.ylabel('Frequency', fontsize=18)
          plt.title('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31)', fontsize=21)
          plt.xlim(min(bins_lag), max(bins_lag))
          plt.grid(True)
          file_name = '{0}.{extention}'.format('S picks Quality control of PhaseNet (2012-01-01 to 2012-01-31)', extention='png')
          fig_lag.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          
          # Plot the histogram probability of all S picks PhaseNet in the common stations 
          all_s_picks_exist_in_catalogtory=all_s_picks_exist_in_catalogtory.iloc[all_dists < self.time_lag_threshold]
          
          fig0, ax0 = plt.subplots(figsize=(20,10))
          label_co = '{0}{1}{2}'.format('Common S picks in catalog with less than 2 seconds time-lag (',all_s_picks_exist_in_catalogtory.shape[0], ' S picks)')
          n_co, bins_co, patches_co = ax0.hist(all_s_picks_exist_in_catalogtory.prob, 21, density=False, alpha=0.75, label=label_co)
          label_all = '{0}{1}{2}'.format('All PhaseNet S picks (',df_S_picks.shape[0], ' S picks)')
          ax0.hist(df_S_picks.prob, 21, density=False, color = "skyblue", ec="skyblue", alpha=0.75, label= label_all)
          steps = (max(all_s_picks_exist_in_catalogtory.prob) - min(all_s_picks_exist_in_catalogtory.prob))/21
          plt.xticks(np.arange(min(all_s_picks_exist_in_catalogtory.prob), max(all_s_picks_exist_in_catalogtory.prob), step=steps))
          plt.xlabel('PhaseNet S picks Probability', fontsize=18)
          plt.ylabel('Frequency', fontsize=18)
          plt.title('PhaseNet output S picks(2012-01-01 to 2012-01-31)', fontsize=21)
          plt.xlim(min(bins_co), max(bins_co))
          plt.grid(True)
          plt.legend(loc='upper right')
          file_name = '{0}.{extention}'.format('PhaseNet output S picks(2012-01-01 to 2012-01-31)', extention='png')
          fig0.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          
          # Perform S picks Quality control of PhaseNet by using existing S picks catalog with defined time lag 
          lag_time = all_dists[all_dists < self.time_lag_threshold]
          for j in range (1, bins_lag.shape[0]):
              select_s_picks=all_s_picks_exist_in_catalogtory.iloc[(lag_time < bins_lag[j]) & (lag_time >= bins_lag[j-1])]
              fig, ax = plt.subplots(figsize=(20,10))
              ax.hist(select_s_picks.prob, 21, density=False, color = "r", ec="r", alpha=0.75)
              plt.xlabel('PhaseNet S picks Probability', fontsize=18)
              plt.ylabel('Frequency', fontsize=18)
              title_name = '{0}{1}{2}{3}{4}'.format('Common S picks in catalog with ',round (bins_lag[j-1]),' - ', round (bins_lag[j]), ' time lag (ms)')
              plt.title(title_name, fontsize=21)
              file_name = '{0}{1}.{extention}'.format('PhaseNet_result_S_bins: ',round (bins_lag[j]), extention='png')
              fig.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
      '''
      
      def mismatched_picks (self, start_time:'str', dt:'int')-> pd.DataFrame:
  
          '''
          This function plot the PhaseNet picks which are mismatched with picks in Catalog.
          Note: make sure all following pickle files exist in the "export_DF_path" directory:
                  - PhaseNet_result_p_picks
                  - PhaseNet_result_s_picks
                  - catalog_p_picks.pkl
                  - catalog_s_picks.pkl
  
          Parameters:
                  - dt(int): delta time in second. example: 7000
                  - start_time: start time. example: "2020-12-31T08:19:57.480000Z"
          Output:
                  - png figure in the "export_DF_path" directory
          '''
          # load all picks in Catalog and PhaseNet based on the existed pickle files in the export_DF_path
  
          # load PhaseNet_result_p_picks.pkl
          with open(os.path.join(self.export_DF_path, "PhaseNet_result_p_picks.pkl"),'rb') as fp:
              PhaseNet_result_p_picks = pickle.load(fp)
  
              # add picks_time to PhaseNet_result_p_picks data frame for synchronization
              PhaseNet_result_p_picks['picks_time']= PhaseNet_result_p_picks.timestamp
  
          # load catalog_p_picks.pkl
          with open(os.path.join(self.export_DF_path, "catalog_p_picks.pkl"),'rb') as fp:
              catalog_DF_P_picks = pickle.load(fp)
              
  
          # load PhaseNet_result_s_picks.pkl
          with open(os.path.join(self.export_DF_path, "PhaseNet_result_s_picks.pkl"),'rb') as fp:
              PhaseNet_result_s_picks = pickle.load(fp)
  
              # add picks_time to PhaseNet_result_s_picks data frame for synchronization
              PhaseNet_result_s_picks['picks_time']= PhaseNet_result_s_picks.timestamp
  
          # load catalog_s_picks.pkl
          with open(os.path.join(self.export_DF_path, "catalog_s_picks.pkl"),'rb') as fp:
              catalog_DF_S_picks = pickle.load(fp)
  
  
  
          # filter path of the existed streams for visualization 
          df_stream_traj = self.filter_stream()
  
          # catalog data prepration for visualization
          catalog_DF_P_picks = self.catalog_data_prepration_for_vis(catalog_DF_P_picks,start_time,dt)
          catalog_DF_S_picks = self.catalog_data_prepration_for_vis(catalog_DF_S_picks,start_time,dt)
  
          # PhaseNet data prepration for visualization
          PhaseNet_result_p_picks = self.phasenet_data_prepration_for_vis(PhaseNet_result_p_picks,start_time,dt)
          PhaseNet_result_s_picks = self.phasenet_data_prepration_for_vis(PhaseNet_result_s_picks,start_time,dt)
  
          # extract mis match picks
          catalog_DF_P_picks = self.extract_mismatch_picks (catalog_DF_P_picks , PhaseNet_result_p_picks) 
          
          df_stream_traj = df_stream_traj.iloc[0:23,:]
  
          fig, ax = plt.subplots(df_stream_traj.shape[0],1,figsize=(40,90),constrained_layout = True)
  
          for i in range (0,int (df_stream_traj.shape[0]/3)):
  
              # Read the N element of the stream
              streamN = obspy.read(df_stream_traj.stream_traj.iloc[3*i])
  
              # Perform Slicing based on the start_time and start_time + dt       
              streamN = streamN.slice (obspy.UTCDateTime(start_time),obspy.UTCDateTime(start_time)+dt)
              
              # Perform filtering based on the freqmin and freqmax
              streamN.filter('bandpass', freqmin= 1, freqmax=20)
  
              # Read the E element of the stream
              streamE = obspy.read(df_stream_traj.stream_traj.iloc[3*i+1])
  
              # Perform Slicing based on the start_time and start_time + dt      
              streamE = streamE.slice (obspy.UTCDateTime(start_time),obspy.UTCDateTime(start_time)+dt)
  
              # Perform filtering based on the freqmin and freqmax
              streamE.filter('bandpass', freqmin= 1, freqmax=20)
  
              # Read the Z element of the stream
              streamZ = obspy.read(df_stream_traj.stream_traj.iloc[3*i+2]) 
  
              # Perform Slicing based on the start_time and start_time + dt 
              streamZ = streamE.slice (obspy.UTCDateTime(start_time),obspy.UTCDateTime(start_time)+dt)
  
              # Perform filtering based on the freqmin and freqmax
              streamZ = streamZ.filter('bandpass', freqmin= 1, freqmax=20)
  
              # add all 3 components and sort them
              streamN += streamE
              streamN += streamZ
              st = streamN.sort()
  
              # Perform visualization data stream
              ax[3*i].set_title(fontsize=25,label="Station: {}".format(st[0].stats.station), fontdict=None, loc='center')
              ax[3*i].plot(st[0].times('matplotlib'), st[0].data, 
                          markersize=1, label = 'E Stream', color = 'k')
              ax[3*i+1].plot(st[1].times('matplotlib'), st[1].data,
                          markersize=1, label = 'N Stream', color = 'k')
              ax[3*i+2].plot(st[2].times('matplotlib'), st[2].data,
                          markersize=1, label = 'Z Stream', color = 'k')
  
              plt.setp(ax[3*i].get_xticklabels(), visible=False)
              plt.setp(ax[3*i+1].get_xticklabels(), visible=False)
  
              
              # Visualization of P Picks imported from catalog
              if catalog_DF_P_picks.shape[0] > 0: 
                  ax[3*i].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_P_picks[catalog_DF_P_picks['station_code'] == st[0].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[0].max()),
                      ymax = (st[0].max()),
                      color='green', linestyle='dashdot', label = 'P pickes from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i].xaxis_date()
  
                  ax[3*i+1].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_P_picks[catalog_DF_P_picks['station_code'] == st[1].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[1].max()),
                      ymax = (st[1].max()),
                      color='green', linestyle='dashdot', label = 'P pickes from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i+1].xaxis_date()
  
                  
                  ax[3*i+2].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_P_picks[catalog_DF_P_picks['station_code'] == st[2].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[2].max()),
                      ymax = (st[2].max()),
                      color='green', linestyle='dashdot', label = 'P pickes from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i+2].xaxis_date()  
  
              # Visualization of S Picks imported from catalog
              if catalog_DF_S_picks.shape[0] > 0:           
                  ax[3*i].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_S_picks[catalog_DF_S_picks['station_code'] == st[0].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[0].max()),
                      ymax = (st[0].max()),
                      color='khaki', linestyle='dashdot', label = 'S picks from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i].xaxis_date()
  
                  ax[3*i+1].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_S_picks[catalog_DF_S_picks['station_code'] == st[1].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[1].max()),
                      ymax = (st[1].max()),
                      color='khaki', linestyle='dashdot', label = 'S picks from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i+1].xaxis_date()
  
              
                  ax[3*i+2].vlines([obspy.UTCDateTime(t).matplotlib_date for t in catalog_DF_S_picks[catalog_DF_S_picks['station_code'] == st[2].stats.station]['picks_time'].tolist()], 
                      ymin = (-st[2].max()),
                      ymax = (st[2].max()),
                      color='khaki', linestyle='dashdot', label = 'S picks from catalog', linewidth=7.0, alpha=0.8)
                  ax[3*i+2].xaxis_date() 
              
              # Visualization of P Picks imported from PhaseNet 
              if PhaseNet_result_p_picks.shape[0] > 0:
                  ax[3*i].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[0].stats.station]['timestamp']], 
                      ymin = (-st[0].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[0].stats.station]['prob'])).tolist(),
                      ymax = (st[0].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[0].stats.station]['prob'])).tolist(),
                      color='b', linestyle='solid', label = 'P picks by PhaseNet', alpha=0.6)
                  ax[3*i].xaxis_date()
  
                  ax[3*i+1].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[1].stats.station]['timestamp']], 
                      ymin = (-st[1].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[1].stats.station]['prob'])).tolist(),
                      ymax = ( st[1].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[1].stats.station]['prob'])).tolist(),
                      color='b', linestyle='solid', label = 'P picks by PhaseNet', alpha=0.6)
                  ax[3*i+1].xaxis_date()
  
                  ax[3*i+2].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[2].stats.station]['timestamp']], 
                      ymin = (-st[2].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[2].stats.station]['prob'])).tolist(),
                      ymax = ( st[2].max()*np.array (PhaseNet_result_p_picks[PhaseNet_result_p_picks['station_code'] == st[2].stats.station]['prob'])).tolist(),
                      color='b', linestyle='solid', label = 'P picks by PhaseNet', alpha=0.6)
                  ax[3*i+2].xaxis_date()
  
  
              
              # Visualization of S Picks imported from PhaseNet 
              if PhaseNet_result_s_picks.shape[0] > 0:
                  ax[3*i].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[0].stats.station]['timestamp']], 
                      ymin = (-st[0].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[0].stats.station]['prob'])).tolist(),
                      ymax = ( st[0].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[0].stats.station]['prob'])).tolist(),
                      color='r', linestyle='solid', label = 'S picks by PhaseNet', alpha=0.6)
                  ax[3*i].xaxis_date()
  
                  ax[3*i+1].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[1].stats.station]['timestamp']], 
                      ymin = (-st[1].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[1].stats.station]['prob'])).tolist(),
                      ymax = ( st[1].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[1].stats.station]['prob'])).tolist(),
                      color='r', linestyle='solid', label = 'S picks by PhaseNet', alpha=0.6)
                  ax[3*i+1].xaxis_date()
  
                  ax[3*i+2].vlines([obspy.UTCDateTime(t).matplotlib_date for t in PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[2].stats.station]['timestamp']], 
                      ymin = (-st[2].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[2].stats.station]['prob'])).tolist(),
                      ymax = ( st[2].max()*np.array (PhaseNet_result_s_picks[PhaseNet_result_s_picks['station_code'] == st[2].stats.station]['prob'])).tolist(),
                      color='r', linestyle='solid', label = 'S picks by PhaseNet', alpha=0.6)
                  ax[3*i+2].xaxis_date()
              
  
              ax[3*i].legend(loc='lower right')
              ax[3*i+1].legend(loc='lower right')
              ax[3*i+2].legend(loc='lower right')
          
          # Save png file in the export_DF_path
          file_name = '{0}{1}.{extention}'.format('Comparison_PhaseNet_catalog_',obspy.UTCDateTime(start_time), extention='png')
          fig.savefig(os.path.join(self.export_DF_path, file_name), facecolor = 'w')
          print ('--------------------------------------------------------------------------------')
          print ('Figure was created and save it in the following directory:')
          print(os.path.join(self.export_DF_path, file_name))
          print ('--------------------------------------------------------------------------------')
          print ('end time is')
          print (st[0].stats.endtime)
          
  
      def filter (self, stream:'obspy') -> obspy:
  
          '''
          Filter the data of all traces in the Stream. This can just support "bandpass" filter.
          Parameters:
                      - stream: obspy stream data
                      - freqmin: minimum frequency
                      - freqmax: maximum frequency
          Return:
                      - filtred stream
          '''
          #sliced_stream = stream.filter('bandpass', freqmin= 1, freqmax=2
          
          stream[0] = stream[0].filter('bandpass', freqmin= self.freqmin, freqmax=self.freqmax)
          stream[1] = stream[1].filter('bandpass', freqmin= self.freqmin, freqmax=self.freqmax)
          stream[2] = stream[2].filter('bandpass', freqmin= self.freqmin, freqmax=self.freqmax)
          
          stream = stream.slice (stream[0].stats.starttime+1000, stream[0].stats.endtime - 1000)
          return stream
  
      def data_slicing (self, starttime:'str', dt:'int', daily_data:'str') -> obspy:
          '''
          Perform Slicing on stream.
              Parameters:
                          - starttime (str): start time for slicing (like 2020-12-31T12:30:58.180000Z)
                          - dt (int): time interval for sclicing
                          - daily_data (str): The name of daily mseed file ( like CX.PB06..HHZ.D.2020.366)
              
              return:
                      - sliced_stream: obspy stream data
          '''
          start=obspy.UTCDateTime(starttime)
          sliced_stream = self.read_data (daily_data).slice (start,start+dt)
          return sliced_stream
  
      def read_data (self, daily_data:'str') -> obspy:
          '''
          Read the mseed daily data and return stream.
              Parameters:
                  - daily_data (str): The name of daily mseed file ( like CX.PB06..HHZ.D.2020.366)
              return:
                      - stream: obspy stream data
          '''
          stream = obspy.read(os.path.join(self.working_direc, 'mseed', '{0}'.format(daily_data)), sep="\t")
          return stream
      
      def filter_stream (self) -> pd.DataFrame: 
          
          '''
          This function filter path of the existed streams based on the "DF_selected_chile_path_file.pkl".
          Make sure the following pickle exists in the "export_DF_path" directory:
              - DF_selected_chile_path_file.pkl
          
          Output:
                  - df_stream_traj (dataframe): This dataframe in a filtered path dataframe based on the the following criterias:
                          1- Rows order based on the stations latitude (CXstatlist.txt)
                          2- Filter the data based on the station list in the "CXstatlist.txt"
                          3- Select stream with have 3 components
                  
                  This data frame contains following columns:
                      'stream_traj': Directory of the existed stream (like: /data2/chile/CHILE_GFZ_ONLINE/2012/CX/PB16/HHN.D/CX.PB16..HHN.D.2012.001')
                      'year': Year of study for visualization according to the input
                      'day' : day of study for visualization according to the input
                      'file_name': stream file name (like: 'CXPB16D2012001')
                      'station_code': station code
                      'network_code': network code
          '''
  
          # load DF_selected_chile_path_file.pkl
          with open(os.path.join(self.export_DF_path, "DF_selected_chile_path_file.pkl"),'rb') as fp:
              DF_selected_chile_path_file = pickle.load(fp)
  
          
          # Path of all streams
          stream_traj = DF_selected_chile_path_file["path"].values.tolist()
  
          # convert stream_traj to data frame
          df_stream_traj = pd.DataFrame({'stream_traj': stream_traj})
  
          # creat year and day columns for filtering
          df_stream_traj[['rest','year', 'day']] = df_stream_traj['stream_traj'].str.rsplit('.', 2, expand=True)
  
          # convert columns to int
          df_stream_traj['year'] = df_stream_traj['year'].astype('int')
          df_stream_traj['day'] = df_stream_traj['day'].astype('int')
  
          # drop the column
          df_stream_traj = df_stream_traj.drop(['rest'], axis=1)
  
          # Filter df_stream_traj based on the given year and day
          df_stream_traj = df_stream_traj[(df_stream_traj['year']>= self.start_year_analysis) & (df_stream_traj['year']<= self.end_year_analysis)]
          df_stream_traj = df_stream_traj[(df_stream_traj['day']>= self.start_day_analysis) & (df_stream_traj['day']<= self.end_day_analysis)]  
          df_stream_traj = self.filter_year_day (df_stream_traj)
  
          # read the station list for visualization
          stations = self.get_stations ()
  
          # Sort stations based on the latitude
          stations = self.sort_stations_latitude(stations)
  
          # Filter df_stream_traj based on the stations dataframe
          df_stream_traj = self.filter_station(df_stream_traj, stations)
  
          # Select the stream with three components
          proper_stream = df_stream_traj.groupby('file_name').file_name.count() > 2
          proper_stream = proper_stream[proper_stream==True]
          df_stream_traj = df_stream_traj[df_stream_traj.file_name.isin(proper_stream.index)]
  
          # Filter df_stream_traj based on the stations dataframe
          df_stream_traj = self.filter_station(df_stream_traj, stations)
  
          # Filter stream station order based on the station code of station
          df_stream_traj = self.DF_sort_station (df_stream_traj, stations)
  
  
          return df_stream_traj
  
      def filter_year_day (self, df_stream_traj: 'pd.DataFrame') -> pd.DataFrame:
  
          '''
          This function filter df_stream_traj based on the given year and day
              Parameters:
                         df_stream_traj (datafarme): file directory 
          '''
          # extract file name
          df_stream_traj[['rest','file_name']] = df_stream_traj['stream_traj'].str.rsplit('/', 1, expand=True)
  
          df_stream_traj[['rest','component']] = df_stream_traj['rest'].str.rsplit('/', 1, expand=True)
  
          # extract station
          df_stream_traj[['rest','station_code']] = df_stream_traj['rest'].str.rsplit('/', 1, expand=True)
  
          # extract network
          df_stream_traj[['rest','network_code']] = df_stream_traj['rest'].str.rsplit('/', 1, expand=True)
  
          # drop the "rest" column
          df_stream_traj = df_stream_traj.drop(['rest','component'], axis=1)
          
          DF_file_name = df_stream_traj["file_name"].str.replace(".","")
          DF_file_name = DF_file_name.str.replace("HHN","")
          DF_file_name = DF_file_name.str.replace("HHE","")
          DF_file_name = DF_file_name.str.replace("HHZ","")
  
          df_stream_traj["file_name"] = DF_file_name
          
          return df_stream_traj
      
  
      def get_stations (self) -> pd.DataFrame:
          '''
          Return the list of stations stored in station_name_list
          
          Output:
              - stations (dataframe): stations dataframe contains the following columns:
                  - 'network_code'
                  - 'station_code' 
                  - 'latidude'
                  - 'longitude'
          '''
          stations = pd.read_csv(os.path.join(self.working_direc, self.station_name_list), sep="\t")
          return stations
      
  
      def sort_stations_latitude(self,stations:'pd.DataFrame') -> pd.DataFrame:
          '''
          Sort station based on the latidude.
              Parameters:
                      - stations (DataFrame): This dataframe contains the following columns:
                          - 'network_code'
                          - 'station_code' 
                          - 'latidude'
                          - 'longitude'
              Output:    
                      - stations (dataframe): sorted stations.
          '''
  
          return stations.sort_values("latidude", ascending=False)
      
      
      def sort_stations_longitude(self,stations:'pd.DataFrame') -> pd.DataFrame:
  
          '''
          Sort station based on longitude.
              Parameters:
                          stations (DataFrame): This data frame contains at least two culomns (station name and coresponding longitude)
          '''
  
          return stations.sort_values("longitude", ascending=False)
          
      
      def filter_station(self,df_stream_traj:'pd.DataFrame', stations:'pd.DataFrame') -> pd.DataFrame:
          '''
          This function filter df_stream_traj dataframe based on station code according to the statlist.txt file.
  
                  Parameters:
                          - df_stream_traj (DF): data frame. This data frame contains following columns:
                              - 'stream_traj': Directory of the existed stream (like: /data2/chile/CHILE_GFZ_ONLINE/2012/CX/PB16/HHN.D/CX.PB16..HHN.D.2012.001')
                              - 'year': Year of study for visualization according to the input
                              - 'day' : day of study for visualization according to the input
                              - 'file_name': stream file name (like: 'CXPB16D2012001')
                              - 'station_code': station code
                              - 'network_code': network code
  
                          - stations (DF): stations data frame. This dataframe contains the following columns:
                              - 'network_code'
                              - 'station_code' 
                              - 'latidude'
                              - 'longitude'
                  
                  Output:
                         - df_stream_traj (DF): filtered dataframe based on the station dataframe
          '''
  
          df_stream_traj = df_stream_traj[df_stream_traj.station_code.isin(stations.station_code)]
  
          return df_stream_traj
      
      def filter_date (self, catalog_DF:'pd.DataFrame',start_time:'str',dt:'int') -> pd.DataFrame:
  
          '''
          This function filter events of catalogs based on the given interval for visualization.
  
              Parameters:
                      - catalog_DF (dataframe) : catalog data frame with the following columns:
                              - 'picks_time'
                              - 'picks_uncertainty'
                              - 'network_code'
                              - 'station_code'
                              - 'phase_hint'
                              - 'origins_time'
                              - 'origins_longitude'
                              - 'origins_latitude',
                              - 'magnitudes' 
                      
                      - start_time (str): UTC data time format like '2012-01-01T17:00:01.820000Z'
                      - dt (int): the time interval in second
              
              output:
                      - catalog_DF (dataframe): prepared catalog dataframe for visualization
          '''
          # calculate end time
          end_time = obspy.UTCDateTime(start_time) + dt
  
          # convert end time to proper format
          end_time = end_time.strftime('%Y-%m-%d %H:%M:%S.%fZ')
  
          # creat a 'UTC_time' column
          catalog_DF['UTC_time'] = pd.to_datetime(catalog_DF.picks_time, utc=True)
  
          # Filter catalog_DF based on the start_time and end_time
          catalog_DF = catalog_DF[(catalog_DF['UTC_time'] >= pd.Timestamp(start_time,tz='UTC')) & (catalog_DF['UTC_time'] <= pd.Timestamp(end_time,tz='UTC'))]
  
          return catalog_DF
      
      def DF_sort_station (self,dataframe:'pd.DataFrame', sorted_station_DF:'pd.DataFrame') -> pd.DataFrame:
          '''
          This function sort the order of stations information based on the sorted_station_DF dataframe.
  
              Parameters:         
                      - dataframe (DF): data frame. This data frame contains following columns:
                          - 'stream_traj': Directory of the existed stream (like: /data2/chile/CHILE_GFZ_ONLINE/2012/CX/PB16/HHN.D/CX.PB16..HHN.D.2012.001')
                          - 'year': Year of study for visualization according to the input
                          - 'day' : day of study for visualization according to the input
                          - 'file_name': stream file name (like: 'CXPB16D2012001')
                          - 'station_code': station code
                          - 'network_code': network code
  
                      - sorted_station_DF: sorted stations data frame. This dataframe contains the following columns:
                          - 'network_code'
                          - 'station_code' 
                          - 'latidude'
                          - 'longitude'
              
              Output:
                      - sorted_dataframe (DF): filtered dataframe based on the sorted_station_DF
  
          '''
          sorted_dataframe  = pd.DataFrame()
  
          for i in range(sorted_station_DF.shape[0]):
              frame = [sorted_dataframe, dataframe[(dataframe['station_code']==sorted_station_DF['station_code'].iloc[i])]]
              sorted_dataframe = pd.concat(frame)
          
          return sorted_dataframe
  
      
      def catalog_data_prepration_for_vis (self, catalog_DF:'pd.DataFrame',start_time:'str',dt:'int') -> pd.DataFrame:
  
          '''
          This function prepare catalog data for visualization. This preparation contains the following:
                  - Filter catalog_DF based on the stations dataframe
                  - Filter catalog based on the given interval for visualization
                  - Filter catalog station order based on the station code of station
              
              Parameters:
                      - catalog_DF (dataframe) : catalog data frame with the following columns:
                              - 'picks_time'
                              - 'picks_uncertainty'
                              - 'network_code'
                              - 'station_code'
                              - 'phase_hint'
                              - 'origins_time'
                              - 'origins_longitude'
                              - 'origins_latitude',
                              - 'magnitudes' 
                      
                      - start_time (str): UTC data time format like '2012-01-01T17:00:01.820000Z'
                      - dt (int): the time interval in second
              
              output:
                      - catalog_DF (dataframe): prepared catalog dataframe for visualization
          '''
          # get the stations stored in txt file
          stations = self.get_stations ()
  
          # Sort the stations based on the latitude
          stations = self.sort_stations_latitude(stations)
  
          catalog_DF['network_station']=catalog_DF['network_code'].astype(str)+'.'+catalog_DF['station_code']
  
          # Filter catalog_DF_P_picks based on the stations dataframe
          catalog_DF = self.filter_station(catalog_DF, stations)
  
          # Filter catalog based on the given interval for visualization
          catalog_DF = self.filter_date (catalog_DF,start_time,dt)
  
          # Filter catalog station order based on the station code of station
          catalog_DF = self.DF_sort_station (catalog_DF, stations)
      
          return catalog_DF
  
      
      def phasenet_data_prepration_for_vis (self, phasenet_picks_DF:'pd.DataFrame' ,start_time:'str',dt:'int'):
  
          '''
          This function prepare phasenet data for visualization. This preparation contains the following:
                  - Filter phasenet_picks_DF based on the stations dataframe
                  - Filter phasenet_picks_DF based on the given interval for visualization
                  - Filter phasenet_picks_DF station order based on the station code of station
              
              Parameters:
                      - phasenet_picks_DF (dataframe) : catalog data frame with the following columns:
                              - 'id': the id of stream in data
                              - 'timestamp': the UTC date time data of picks
                              - 'prob': Probability of picks
                              - 'type': Type of picks (S or P)
                              - 'picks_time': the UTC date time data of picks (like timestamp)
                      
                      - start_time (str): UTC data time format like '2012-01-01T17:00:01.820000Z'
                      - dt (int): the time interval in second
              
              output:
                      - phasenet_picks_DF (dataframe): prepared catalog dataframe for visualization
          '''
          # read the station data
          stations = self.get_stations ()
  
          # Sort the station based on the latitude
          stations = self.sort_stations_latitude(stations)
  
          # creat extra columns
          phasenet_picks_DF[['network', 'others']] = phasenet_picks_DF['id'].str.split('.', 1, expand=True)
          phasenet_picks_DF[['station_code', 'date']] = phasenet_picks_DF['others'].str.split('.', 1, expand=True)
          phasenet_picks_DF = phasenet_picks_DF.drop(['date', 'others'], axis=1)
  
  
          # Filter phasenet_picks_DF based on the stations dataframe
          phasenet_picks_DF = self.filter_station(phasenet_picks_DF, stations)
  
          # Filter phasenet_picks_DF based on the given interval for visualization
          phasenet_picks_DF = self.filter_date (phasenet_picks_DF,start_time,dt)
      
          # Filter phasenet_picks_DF station order based on the station code of station
          phasenet_picks_DF = self.DF_sort_station (phasenet_picks_DF, stations)
      
          return phasenet_picks_DF
  
  
      def extract_mismatch_picks (self, catalog_DF:'pd.DataFrame', phasenet_DF:'pd.DataFrame') -> pd.DataFrame:
  
          '''
          This function extract mismatch picks between phasenet and catalog.
                  parameters:
                          - catalog_DF (dataframe): catalog picks Data Frame (s or p)
                          - phasenet_DF (dataframe): phasenet picks Data Frame (s or p)
                  
                  output:
                          - mis_match_catalog_DF (dataframe)
          '''
          
          # Creat an empty DataFrame
          mis_match_catalog_DF = pd.DataFrame(index =[])
  
          # Find common station code in phasenet_DF and catalog_DF
          common_station = np.intersect1d(catalog_DF.station_code.unique(), phasenet_DF.station_code.unique())
          
          # loop over all common station
          for i in common_station:
              bo = catalog_DF['station_code']==i
              catalog_filter_station = catalog_DF[(bo==True)]
              ao = phasenet_DF['station_code']==i
              phasenet_filter_station = phasenet_DF[(ao==True)]
  
              # Convert UTC time to datetime64[ms] (millisecond)
              a = catalog_filter_station.picks_time.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
              b = phasenet_filter_station.timestamp.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
  
              # Calculate P1 norme of all datetime64[m
              dist_mat = distance_matrix(a,b, p=1)
              dists = np.min(dist_mat, axis=1)
  
              # filter the catalog_filter_station based on the given time_lag_threshold
              catalog_filter_station = catalog_filter_station[dists >= self.time_lag_threshold]
              mis_match_catalog_DF = pd.concat([mis_match_catalog_DF, catalog_filter_station], axis=0)
      
          return mis_match_catalog_DF
      
      def proximity_matrix (self, phase_hint:'str') -> Tuple:
  
          '''
          This function measure the proximity matrix of PhaseNet performance assuming that the catalog is the ground-true.
          In order to measure the proximity matrix, the following variables are using:
  
                  - True Positive (Tp): 
                                      1- Peak probabilities of the phase is above 0.5 
                                      2- Arrival-time residuals is less than 0.1 second
  
                  - False Positive (Fp):
                                      1- Peak probabilities of the phase is above 0.5 
                                      2- Arrival-time residuals is more than 0.1 second
  
                  - True Negative (Tn): 
                                      1- Peak probabilities of the phase is below 0.5 
                                      2- Arrival-time residuals is less than 0.1 second
  
                  - False Negative (Fn):
                                      1- Peak probabilities of the phase is below 0.5 
                                      2- Arrival-time residuals is more than 0.1 second            
          
  
                  parameters:
                          - Phase_hint (S or P): phase hint
  
                  Outputs:
                          - Precision
                          - recall
                          - F1 score
          '''
          if phase_hint == 'P':
              # load PhaseNet_result_p_picks.pkl
              with open(os.path.join(self.export_DF_path, "PhaseNet_result_p_picks.pkl"),'rb') as fp:
                  phasenet_DF = pickle.load(fp)
  
                  # add picks_time to PhaseNet_result_p_picks data frame for synchronization
                  phasenet_DF['picks_time']= phasenet_DF.timestamp
  
              # load catalog_p_picks.pkl
              with open(os.path.join(self.export_DF_path, "catalog_p_picks.pkl"),'rb') as fp:
                  catalog_DF = pickle.load(fp)
              
          else:
              # load PhaseNet_result_s_picks.pkl
              with open(os.path.join(self.export_DF_path, "PhaseNet_result_s_picks.pkl"),'rb') as fp:
                  phasenet_DF = pickle.load(fp)
  
                  # add picks_time to PhaseNet_result_s_picks data frame for synchronization
                  phasenet_DF['picks_time']= phasenet_DF.timestamp
  
              # load catalog_s_picks.pkl
              with open(os.path.join(self.export_DF_path, "catalog_s_picks.pkl"),'rb') as fp:
                  catalog_DF = pickle.load(fp)
  
          # creat extra columns
          phasenet_DF[['network', 'others']] = phasenet_DF['id'].str.split('.', 1, expand=True)
          phasenet_DF[['station_code', 'date']] = phasenet_DF['others'].str.split('.', 1, expand=True)
          phasenet_DF = phasenet_DF.drop(['date', 'others'], axis=1)
  
          # Intialize the true postive counter, false postive counter, true negative counter, and false negative counter 
          true_pos_count = 0
          false_pos_count = 0
  
          true_neg_count = 0
          false_neg_count = 0
          all = 0
  
          # Find common station code in phasenet_DF and catalog_DF
          common_station = np.intersect1d(catalog_DF.station_code.unique(), phasenet_DF.station_code.unique())
          
          # loop over all common stations
          for i in common_station:
              bo = catalog_DF['station_code']==i
              catalog_filter_station = catalog_DF[(bo==True)]
              ao = phasenet_DF['station_code']==i
              phasenet_filter_station = phasenet_DF[(ao==True)]
  
              # Convert UTC time to datetime64[ms] (millisecond)
              a = catalog_filter_station.picks_time.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
              b = phasenet_filter_station.timestamp.to_numpy(dtype='datetime64[ms]')[:, np.newaxis].astype("float")
  
              # Calculate P1 norme of all datetime64[m
              dist_mat = distance_matrix(a,b, p=1)
              dists = np.min(dist_mat, axis=1)
  
              all = all + catalog_filter_station.shape[0]
              probability = phasenet_filter_station.prob.iloc[np.argmin(dist_mat, axis=1)].to_numpy()
  
              #calculate True poitive
              Tp = catalog_filter_station[(dists <= self.time_lag_threshold) & (probability >= 0.5)].shape[0]
              true_pos_count = true_pos_count + Tp
  
              # calculate false poitive
              Fp = catalog_filter_station[(dists > self.time_lag_threshold) & (probability >= 0.5)].shape[0]
              false_pos_count = false_pos_count + Fp
  
              # calculate true negative
              Tn = catalog_filter_station[(dists > self.time_lag_threshold) & (probability < 0.5)].shape[0]
              true_neg_count = true_neg_count + Tn
  
              # calculate false negative
  
              Fn = catalog_filter_station[(dists <= self.time_lag_threshold) & (probability < 0.5)].shape[0]
              false_neg_count = false_neg_count + Fn
  
          # calculate Precision
          precision = true_pos_count /(true_pos_count + false_pos_count)
  
          # calculate Recall
          recall = true_pos_count /(true_pos_count + false_neg_count)
  
          # calculate f1 score
          f1_score = (2*precision*recall)/(precision + recall)
  
          # calculate TRP
          TPR = true_pos_count /(true_pos_count + false_neg_count)
  
          # calculate FPR
          FPR = false_pos_count /(false_pos_count + true_neg_count)
  
          return precision, recall, f1_score, TPR, FPR
          
  
  
  if __name__ == "__main__":
  
      phasenet_direc = '/home/javak/phasenet_chile-subduction-zone'
      chile_GFZ_online_direc = '/data2/chile/CHILE_GFZ_ONLINE'
      export_DF_path = '/home/javak/Sample_data_chile/Comparing PhaseNet and Catalog'
      export_mseed_path = '/home/javak/Sample_data_chile/mseed'
      working_direc = '/home/javak/Sample_data_chile'
      picks_name = '2012_events.pkl'
  
      start_year_analysis = 2012
      start_day_analysis = 1
      end_year_analysis = 2012
      end_day_analysis = 1
      analysis = False
  
      apply_filter = False
      freqmin = 0.2
      freqmax = 10
  
      time_lag_threshold = 2000
  
      station_name_list = 'CXstatlist.txt'
  
      obj = PhaseNet_Analysis (phasenet_direc,chile_GFZ_online_direc,export_DF_path, 
                              export_mseed_path, working_direc, picks_name, 
                              start_year_analysis, start_day_analysis,
                              end_year_analysis, end_day_analysis, analysis, 
                              time_lag_threshold, station_name_list,
                              apply_filter, freqmin, freqmax)
      
      
      start_time ="2012-01-01T05:15:38.820000Z"
      dt = 25
      result = obj.mismatched_picks(start_time,dt)
      #result = obj.get_stations()
      #result = obj()
      #precision, recall, f1_score, TPR, FPR = obj.proximity_matrix('P')
      #print('precision is \n', precision)
      #print('recall is \n', recall)
      #print('f1_score is \n', f1_score)
  
  
  