import sys import numpy as np import matplotlib.pyplot as plt import pandas as pd # import os import scipy.io import logging from geopy.distance import geodesic from geopy.point import Point from util.base_logger import getDefaultLogger def main(Input_catalog, Input_injection_rate, time_win_in_hours, time_step_in_hour, time_win_type, End_time, ev_limit, Model_index, Mc, Mu, G, ssd, C, b_value_type, cl, mag_name, time_inj=None, time_shut_in=None, time_big_ev=None, Inpar = ["SER", "dv", "d_num"]): """ Main function for application: MAXIMUM_MAGNITUDE_DETERMINISTIC_MODELS Arguments: Input catalog: path to input file of type 'catalog' Input injection rate: path to input file of type 'injection_rate' **kwargs: Model paramters. Returns: 'PLOT_Mmax_param': Plot of all results (also check 'application.log' for more info) 'DATA_Mmax_param': Results with 'csv' format 'application.log': Logging file """ f_indx = Model_index b_method = b_value_type Plot_flag = 0 # Setting up the logfile-------------------------------- logger = getDefaultLogger(__name__) # Importing utilities ---------------------- logger.info("Import utilities") from util.Find_idx4Time import Find_idx4Time from util.CandidateEventsTS import CandidateEventsTS from util.M_max_models import M_max_models def latlon_to_enu(lat, lon, alt): lat0 = np.mean(lat) lon0 = np.mean(lon) alt0 = 0 # Reference point (origin) origin = Point(lat0, lon0, alt0) east = np.zeros_like(lon) north = np.zeros_like(lat) up = np.zeros_like(alt) for i in range(len(lon)): # Target point target = Point(lat[i], lon[i], alt[i]) # Calculate East-North-Up east[i] = geodesic((lat0, lon0), (lat0, lon[i])).meters if lon[i] < lon0: east[i] = -east[i] north[i] = geodesic((lat0, lon0), (lat[i], lon0)).meters if lat[i] < lat0: north[i] = -north[i] up = alt - alt0 return east, north, up # Importing data logger.info("Import data") mat = scipy.io.loadmat(Input_catalog) Cat_structure = mat['Catalog'] Cat_id, Cat_t, Cat_m = [], [], [] Cat_x, Cat_y, Cat_z = [], [], [] Cat_lat, Cat_lon, Cat_elv, Cat_depth = [], [], [], [] for i in range(1,Cat_structure.shape[1]): if Cat_structure.T[i][0][0][0] == 'X': Cat_x = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_x)): raise ValueError("Catalog-X contains infinite value") if Cat_structure.T[i][0][3][0] == 'km': Cat_x *= 1000 if Cat_structure.T[i][0][0][0] == 'Lat': Cat_lat = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_lat)): raise ValueError("Catalog-Lat contains infinite value") if Cat_structure.T[i][0][0][0] == 'Y': Cat_y = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_y)): raise ValueError("Catalog-Y contains infinite value") if Cat_structure.T[i][0][3][0] == 'km': Cat_y *= 1000 if Cat_structure.T[i][0][0][0] == 'Long': Cat_lon = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_lon)): raise ValueError("Catalog-Long contains infinite value") if Cat_structure.T[i][0][0][0] == 'Z': Cat_z = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_z)): raise ValueError("Catalog-Z contains infinite value") if Cat_structure.T[i][0][3][0] == 'km': Cat_z *= 1000 if Cat_structure.T[i][0][0][0] == 'Depth': Cat_depth = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_depth)): raise ValueError("Catalog-Depth contains infinite value") if Cat_structure.T[i][0][3][0] == 'km': Cat_depth *= 1000 if Cat_structure.T[i][0][0][0] == 'Elevation': Cat_elv = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_elv)): raise ValueError("Catalog-Elevation contains infinite value") if Cat_structure.T[i][0][3][0] == 'km': Cat_elv *= 1000 if Cat_structure.T[i][0][0][0] == 'Time': Cat_t = Cat_structure.T[i][0][2] if not all(np.isfinite(Cat_t)): raise ValueError("Catalog-Time contains infinite value") if Cat_structure.T[i][0][0][0] == mag_name: Cat_m = Cat_structure.T[i][0][2] # if not all(np.isfinite(Cat_m)): if np.argwhere(all(np.isfinite(Cat_m))): raise ValueError("Catalog-Magnitude contains infinite value") if any(Cat_x): Cat_id = np.linspace(0,Cat_x.shape[0],Cat_x.shape[0]).reshape((Cat_x.shape[0],1)) arg = (Cat_id, Cat_x, Cat_y, Cat_z, Cat_t, Cat_m) Cat = np.concatenate(arg, axis=1) elif any(Cat_lat): if any(Cat_elv): Cat_x, Cat_y, Cat_z = latlon_to_enu(Cat_lat, Cat_lon, Cat_elv) elif any(Cat_depth): Cat_x, Cat_y, Cat_z = latlon_to_enu(Cat_lat, Cat_lon, Cat_depth) else: raise ValueError("Catalog Depth or Elevation is not available") Cat_id = np.linspace(0,Cat_x.shape[0],Cat_x.shape[0]).reshape((Cat_x.shape[0],1)) arg = (Cat_id, Cat_x, Cat_y, Cat_z, Cat_t, Cat_m) Cat = np.concatenate(arg, axis=1) else: raise ValueError("Catalog data are not available") mat = scipy.io.loadmat(Input_injection_rate) Inj_structure = mat['d'] if 'Date' in Inj_structure.dtype.names: inj_date = Inj_structure['Date'][0,0] inj_rate = Inj_structure['Injection_rate'][0,0] Hyd = np.concatenate((inj_date,inj_rate), axis=1) else: raise ValueError("Injection data are not available") if Cat[0,4]>np.max(Hyd[:,0]) or Hyd[0,0]>np.max(Cat[:,4]): raise ValueError('Catalog and injection data do not have time coverage!') if Hyd[0,0] < Cat[0,4]: same_time_idx = Find_idx4Time(Hyd[:,0], Cat[0,4]) Hyd = Hyd[same_time_idx:,:] Hyd[:,0] = (Hyd[:,0] - Cat[0,4])*24*3600 Cat[:,4] = (Cat[:,4] - Cat[0,4])*24*3600 logger.info('Start of the computations is based on the time of catalog data.') # Model dictionary Feat_dic = { 'indx' :[0 ,1 ,2 ,3 ,4 ,5 ], 'Full_names':['All_M_max' ,'McGarr' ,'Hallo' ,'Li' ,'van-der-Elst','Shapiro' ], 'Short_name':['max_all' , 'max_mcg', 'max_hlo', 'max_li', 'max_vde' , 'max_shp'], 'f_num' :[5 ,4 ,4 ,1 ,6 ,4 ], 'Param' :[{'Mc': 0.8, 'b_method': ['b'], 'cl': [0.37], 'Inpar': ['dv','d_num'], 'Mu':0.6, 'G': 35*10**(9), 'ssd': 3*10**6, 'C': 0.95, 'ev_limit': 20, 'num_bootstraps': 100}] } Feat_dic['Param'][0]['Inpar'] = Inpar Feat_dic['Param'][0]['ev_limit'] = ev_limit num_bootstraps = 100 # Number of bootstraping for standard error computation Feat_dic['Param'][0]['num_bootstraps'] = num_bootstraps if f_indx in [0,1,2,4]: if Mc < np.min(Cat[:,-1]) or Mc > np.max(Cat[:,-1]): raise ValueError("Completeness magnitude (Mc) is out of magnitude range") Feat_dic['Param'][0]['Mc'] = Mc if f_indx in [0,1,2]: if Mu < 0.2 or Mu > 0.8: raise ValueError("Friction coefficient (Mu) must be between [0.2, 0.8]") Feat_dic['Param'][0]['Mu'] = Mu if f_indx in [0,1,2,3]: if G < 1*10**(9) or G > 100*10**(9): raise ValueError("Shear modulus of reservoir rock (G) must be between [1, 100] GPa") Feat_dic['Param'][0]['G'] = G if f_indx in [0,5]: if ssd < 0.1*10**(6) or ssd > 100*10**(6): raise ValueError("Static stress drop (ssd) must be between [0.1, 100] MPa") Feat_dic['Param'][0]['ssd'] = ssd if f_indx in [0,5]: if C < 0.5 or C > 5: raise ValueError("Geometrical constant (C) of Shaprio's model must be between [0.5, 5]") Feat_dic['Param'][0]['C'] = C if f_indx in [0,1,2,4]: Feat_dic['Param'][0]['b_method'] = b_method if f_indx in [0,4]: for cl_i in cl: if cl_i < 0 or cl_i > 1: raise ValueError("Confidence level (cl) of van der Elst model must be between [0, 1]") Feat_dic['Param'][0]['cl'] = cl # Setting up based on the config and model dic -------------- ModelClass = M_max_models() Model_name = Feat_dic['Full_names'][f_indx] ModelClass.f_name = Feat_dic['Short_name'][f_indx] f_num = Feat_dic['f_num'][f_indx] if Feat_dic['Param'][0]['Mc']: ModelClass.Mc = Mc else: Mc = np.min(Cat[:,-1]) ModelClass.Mc = Mc time_win = time_win_in_hours*3600 # in sec ModelClass.time_win = time_win if f_indx == 0: # Only first b_methods is used Feat_dic['Param'][0]['b_method'] = b_method[0] ModelClass.b_method = Feat_dic['Param'][0]['b_method'] logger.info(f"All models are based on b_method: { b_method[0]}") Feat_dic['Param'][0]['cl'] = [cl[0]] ModelClass.cl = Feat_dic['Param'][0]['cl'] logger.info(f"All models are based on cl: { cl[0]}") ModelClass.Mu = Feat_dic['Param'][0]['Mu'] ModelClass.G = Feat_dic['Param'][0]['G'] ModelClass.ssd = Feat_dic['Param'][0]['ssd'] ModelClass.C = Feat_dic['Param'][0]['C'] ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps'] if f_indx == 1 or f_indx == 2: ModelClass.b_method = Feat_dic['Param'][0]['b_method'] ModelClass.Mu = Feat_dic['Param'][0]['Mu'] ModelClass.G = Feat_dic['Param'][0]['G'] ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps'] Feat_dic['Param'][0]['cl'] = [None] if f_indx == 3: Feat_dic['Param'][0]['b_method'] = [None] ModelClass.G = Feat_dic['Param'][0]['G'] Feat_dic['Param'][0]['cl'] = [None] if f_indx == 4: ModelClass.b_method = Feat_dic['Param'][0]['b_method'] ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps'] ModelClass.G = Feat_dic['Param'][0]['G'] ModelClass.cl = Feat_dic['Param'][0]['cl'] if f_indx == 5: ModelClass.ssd = Feat_dic['Param'][0]['ssd'] ModelClass.C = Feat_dic['Param'][0]['C'] Feat_dic['Param'][0]['b_method'] = [None] Feat_dic['Param'][0]['cl'] = [None] # Output dictionary -------------------------------- Output_dict = { 'Type' :['idx', 'Time[day]'], 'label' :[None, None], 'b_method' :[None, None], 'cl' :[None, None], } c_out = 2 if any(Feat_dic['Param'][0]['b_method']) and any(Feat_dic['Param'][0]['cl']): for i in range(len(Feat_dic['Param'][0]['b_method'])): for j in range(len(Feat_dic['Param'][0]['cl'])): if f_indx == 0: # f_index == 0 for i in Feat_dic['Full_names'][1:]: Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(i) Output_dict['b_method'].append(None) Output_dict['cl'].append(None) c_out += 1 elif j == 0: # f_index == 4 Output_dict['Type'].append('b_value') Output_dict['label'].append('b_value') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('b_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Seismogenic Index') Output_dict['label'].append('SI') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('si_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(Model_name) Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(Feat_dic['Param'][0]['cl'][j]) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('M_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) c_out += 6 else: # f_index == 4 Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(Model_name) Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(Feat_dic['Param'][0]['cl'][j]) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('M_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) c_out += 2 elif any(Feat_dic['Param'][0]['b_method']): # f_index == 1, 2 for i in range(len(Feat_dic['Param'][0]['b_method'])): Output_dict['Type'].append('b_value') Output_dict['label'].append('b_value') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('b_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(Model_name) Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) Output_dict['Type'].append('Standard Error') Output_dict['label'].append('M_std_err') Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i]) Output_dict['cl'].append(None) c_out += 4 elif f_indx == 5: # f_index == 5 # for i in ['L(max)','L(int)','L(min)', 'L(avg)']: # Output_dict['Type'].append('Length[m]') # Output_dict['label'].append(i) # Output_dict['b_method'].append(None) # Output_dict['cl'].append(None) # Output_dict['Type'].append('Maximum magnitude') # Output_dict['label'].append(Model_name) # Output_dict['b_method'].append(None) # Output_dict['cl'].append(None) # c_out += 5 for i in ['Shapiro (Lmax)','Shapiro (Lint)','Shapiro (Lmin)', 'Shapiro (Lavg)']: Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(i) Output_dict['b_method'].append(None) Output_dict['cl'].append(None) c_out += 4 else: # f_index == 3 Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append(Model_name) Output_dict['b_method'].append(None) Output_dict['cl'].append(None) c_out += 1 # Add True maximum magnitude to the outputs Output_dict['Type'].append('Maximum magnitude') Output_dict['label'].append('True Max-Mag') Output_dict['b_method'].append(None) Output_dict['cl'].append(None) c_out += 1 # if any(Feat_dic['Param'][0]['Inpar']): # Input parameters if Inpar: for i in Feat_dic['Param'][0]['Inpar']: if i == 'd_num': Output_dict['Type'].append('Number of events') Output_dict['label'].append('Ev_Num') elif i == 'dv': Output_dict['Type'].append('Volume[m3]') Output_dict['label'].append('Vol') elif i == 'Mo': Output_dict['Type'].append('Seismic moment') Output_dict['label'].append('Mo') elif i == 'SER': Output_dict['Type'].append('Seismic efficiency ratio') Output_dict['label'].append('SER') Output_dict['b_method'].append(None) Output_dict['cl'].append(None) c_out += 1 # Functions ------------------------------ # Computing in extending time or time window def Feature_in_Time(In_arr): SER_l = 0 SER_c = 0 i_row = 0 for i in np.arange(1,Times): # Defining time window and data if time_win_type == 0: ModelClass.time_win = i*time_step candidate_events = CandidateEventsTS(Cat, i*time_step, None, i*time_step) else: candidate_events = CandidateEventsTS(Cat, i*time_step, None, time_win) data = candidate_events.filter_data() # USER: Check if cumulative dv is correctly computed !!! end_time_indx = Find_idx4Time(Hyd[:,0], i*time_step) # ModelClass.dv = np.sum(Hyd[:end_time_indx,1]) ModelClass.dv = np.trapz(Hyd[:end_time_indx,1], Hyd[:end_time_indx,0])/60 if len(data) > Feat_dic['Param'][0]['ev_limit'] and ModelClass.dv > 0: ModelClass.Mo = np.sum(10**(1.5*data[:end_time_indx,5]+9.1)) if f_indx != 5: # Parameters have not been assinged for f_index == 5 SER_c = ModelClass.Mo/(2.0*ModelClass.G*ModelClass.dv) SER_l = np.max((SER_c, SER_l)) ModelClass.SER = SER_l ModelClass.data = data In_arr[i_row,0] = i # index In_arr[i_row,1] = i*time_step # Currecnt time c = 2 if any(Feat_dic['Param'][0]['b_method']) and any(Feat_dic['Param'][0]['cl']): for ii in range(len(Feat_dic['Param'][0]['b_method'])): ModelClass.b_method = Feat_dic['Param'][0]['b_method'][ii] for jj in range(len(Feat_dic['Param'][0]['cl'])): # f_index == 4 ModelClass.cl = Feat_dic['Param'][0]['cl'][jj] Out = ModelClass.ComputeModel() if jj == 0: In_arr[i_row,c:c+f_num] = Out c += f_num else: In_arr[i_row,c:c+2] = Out[-2:] c += 2 elif any(Feat_dic['Param'][0]['b_method']): # f_index == 1, 2 for ii in range(len(Feat_dic['Param'][0]['b_method'])): ModelClass.b_method = Feat_dic['Param'][0]['b_method'][ii] In_arr[i_row,c:c+f_num] = ModelClass.ComputeModel() c += f_num else: # f_index = 0, 3, 5 In_arr[i_row,c:c+f_num] = ModelClass.ComputeModel() c += f_num # Compute true maximum magnitude in the data In_arr[i_row,c] = np.max(data[:end_time_indx,5]) c += 1 if Inpar: for ii in range(len(Feat_dic['Param'][0]['Inpar'])): # Add input parameters if Feat_dic['Param'][0]['Inpar'][ii] == 'd_num': In_arr[i_row,c+ii] = data.shape[0] elif Feat_dic['Param'][0]['Inpar'][ii] == 'dv': In_arr[i_row,c+ii] = ModelClass.dv elif Feat_dic['Param'][0]['Inpar'][ii] == 'Mo': In_arr[i_row,c+ii] = ModelClass.Mo elif Feat_dic['Param'][0]['Inpar'][ii] == 'SER': if ModelClass.SER: In_arr[i_row,c+ii] = ModelClass.SER i_row += 1 return In_arr[:i_row,:] # Plotting models and parameters def Plot_feature(Model_Param_array,Output_dict): myVars = locals() # Function for findin min-max of all similar parameters def Extermom4All(Model_Param_array, itr_loc): Mat1D = np.reshape(Model_Param_array[:,itr_loc], -1) NotNone = np.isfinite(Mat1D) if np.min(Mat1D[NotNone])>0: return [np.min(Mat1D[NotNone])*0.95, np.max(Mat1D[NotNone])*1.05] elif np.min(Mat1D[NotNone])<0 and np.max(Mat1D[NotNone])>0: return [np.min(Mat1D[NotNone])*1.05, np.max(Mat1D[NotNone])*1.05] elif np.max(Mat1D[NotNone])<0: return [np.min(Mat1D[NotNone])*1.05, np.max(Mat1D[NotNone])*0.95] # Function for setting relevant lagends in the plot def Legend_label(loc): l = Output_dict_c['label'][loc] if Output_dict_c['b_method'][loc]: if Output_dict_c['cl'][loc]: l+='('+Output_dict_c['b_method'][loc]+', cl='+str(Output_dict_c['cl'][loc])+')' else: l+='('+Output_dict_c['b_method'][loc]+')' return l c_NotNone = [] # Removing all parameters with None or constant value for i in range(Model_Param_array.shape[1]): NotNone = np.isfinite(Model_Param_array[:,i]) Eq_value = np.mean(Model_Param_array[:,i]) if any(NotNone) and Eq_value != Model_Param_array[0,i]: c_NotNone.append(i) else: logger.info(f"No-PLOT: All values of {Output_dict['Type'][i]} are {Model_Param_array[0,i]}!") if len(c_NotNone) > 1: Model_Param_array = Model_Param_array[:,c_NotNone] # New output dictionary based on valid parameters for plotting Output_dict_c = {'Type':[], 'label':[], 'b_method':[], 'cl':[]} for i in range(len(c_NotNone)): Output_dict_c['Type'].append(Output_dict['Type'][c_NotNone[i]]) Output_dict_c['label'].append(Output_dict['label'][c_NotNone[i]]) Output_dict_c['b_method'].append(Output_dict['b_method'][c_NotNone[i]]) Output_dict_c['cl'].append(Output_dict['cl'][c_NotNone[i]]) coloring=['blue','g','r','c','m','y', 'brown', 'darkolivegreen', 'teal', 'steelblue', 'slateblue', 'purple', 'darksalmon', '#c5b0d5', '#c49c94', '#e377c2', '#f7b6d2', '#7f7f7f', '#c7c7c7', '#bcbd22', '#dbdb8d', '#17becf', '#9edae5', 'brown', 'darkolivegreen', 'teal', 'steelblue', 'slateblue',] # All parameters to be plotted All_vars = Output_dict_c['Type'][2:] Uniqe_var = list(dict.fromkeys([s for s in All_vars if 'Standard Error' not in s])) #list(set(All_vars)) # defining handels and labels to make final legend All_handels = ['p0'] for i in range(1,len(All_vars)): All_handels.append('p'+str(i)) handels = [] labels = [] # itr_loc: location of paramteres with similar type itr_loc = np.where(np.array(All_vars) == Uniqe_var[0])[0]+2 fig, myVars[Output_dict_c['Type'][0]] = plt.subplots(1,1,figsize=(8+int(len(All_vars)/3),6)) fig.subplots_adjust(right=1-len(Uniqe_var)*0.09) if Output_dict_c['label'][itr_loc[0]] == 'True Max-Mag': # plot with dash-line myVars[All_handels[itr_loc[0]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= 'k', ls='--', lw = 2) else: myVars[All_handels[itr_loc[0]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]]) handels.append(All_handels[itr_loc[0]-2]) labels.append(Legend_label(itr_loc[0])) myVars[Output_dict_c['Type'][0]].set_ylabel(Output_dict_c['Type'][itr_loc[0]]) myVars[Output_dict_c['Type'][0]].set_ylim(Extermom4All(Model_Param_array, itr_loc)[0], Extermom4All(Model_Param_array, itr_loc)[1]) myVars[Output_dict_c['Type'][0]].set_xlabel('Day (From start of the recording)') if End_time: myVars[Output_dict_c['Type'][0]].set_xlim(0,End_time) # Plotting statndard error (if exists) if itr_loc[0]+1 < len(Output_dict_c['Type']) and Output_dict_c['Type'][itr_loc[0]+1] == 'Standard Error': myVars[Output_dict_c['Type'][0]].fill_between(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]] - Model_Param_array[:,itr_loc[0]+1], Model_Param_array[:,itr_loc[0]] + Model_Param_array[:,itr_loc[0]+1], color= coloring[itr_loc[0]], alpha=0.1) # Plotting similar parameters on one axis for j in range(1,len(itr_loc)): if Output_dict_c['label'][itr_loc[j]] == 'True Max-Mag': # plot with dash-line myVars[All_handels[itr_loc[j]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= 'k', ls='--', lw = 2) else: myVars[All_handels[itr_loc[j]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]]) handels.append(All_handels[itr_loc[j]-2]) labels.append(Legend_label(itr_loc[j])) # Plotting statndard error (if exists) if itr_loc[0]+1 < len(Output_dict_c['Type']) and Output_dict_c['Type'][itr_loc[j]+1] == 'Standard Error': myVars[Output_dict_c['Type'][0]].fill_between(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]] - Model_Param_array[:,itr_loc[j]+1], Model_Param_array[:,itr_loc[j]] + Model_Param_array[:,itr_loc[j]+1], color= coloring[itr_loc[j]], alpha=0.1) first_itr = 0 # Check if there is any more parameter to be plotted in second axes # The procedure is similar to last plots. if len(Uniqe_var) > 1: for i in range(1,len(Uniqe_var)): itr_loc = np.where(np.array(All_vars) == Uniqe_var[i])[0]+2 myVars[Uniqe_var[i]] = myVars[Output_dict_c['Type'][0]].twinx() # if it is third or more axis, make a distance between them if first_itr == 0: first_itr += 1 set_right = 1 else: set_right = 1 + first_itr*0.2 first_itr += 1 myVars[Uniqe_var[i]].spines.right.set_position(("axes", set_right)) if Output_dict_c['label'][itr_loc[0]] == 'True Max-Mag': # plot with dash-line myVars[All_handels[itr_loc[0]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= 'k', ls='--', lw = 2) else: myVars[All_handels[itr_loc[0]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]]) handels.append(All_handels[itr_loc[0]-2]) labels.append(Legend_label(itr_loc[0])) myVars[Uniqe_var[i]].set_ylabel(Output_dict_c['Type'][itr_loc[0]]) myVars[Uniqe_var[i]].yaxis.label.set_color(coloring[itr_loc[0]]) myVars[Uniqe_var[i]].spines["right"].set_edgecolor(coloring[itr_loc[0]]) myVars[Uniqe_var[i]].tick_params(axis='y', colors= coloring[itr_loc[0]]) myVars[Uniqe_var[i]].set_ylim(Extermom4All(Model_Param_array, itr_loc)[0], Extermom4All(Model_Param_array, itr_loc)[1]) if itr_loc[0]+1 < len(Output_dict_c['Type']) and Output_dict_c['Type'][itr_loc[0]+1] == 'Standard Error': myVars[Uniqe_var[i]].fill_between(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]] - Model_Param_array[:,itr_loc[0]+1], Model_Param_array[:,itr_loc[0]] + Model_Param_array[:,itr_loc[0]+1], color= coloring[itr_loc[0]], alpha=0.1) for j in range(1,len(itr_loc)): if Output_dict_c['label'][itr_loc[j]] == 'True Max-Mag': # plot with dash-line myVars[All_handels[itr_loc[j]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= 'k', ls = '--', lw = 2) else: myVars[All_handels[itr_loc[j]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]]) handels.append(All_handels[itr_loc[j]-2]) labels.append(Legend_label(itr_loc[j])) if itr_loc[j]+1 < len(Output_dict_c['Type']) and Output_dict_c['Type'][itr_loc[j]+1] == 'Standard Error': myVars[Uniqe_var[i]].fill_between(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]] - Model_Param_array[:,itr_loc[j]+1], Model_Param_array[:,itr_loc[j]] + Model_Param_array[:,itr_loc[j]+1], color= coloring[itr_loc[j]], alpha=0.1) # If there are timing, plot them as vertical lines if time_inj: myVars['l1'], = plt.plot([time_inj,time_inj], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='--', c='k') handels.append('l1') labels.append('Start-inj') if time_shut_in: myVars['l2'], = plt.plot([time_shut_in,time_shut_in], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='-.', c='k') handels.append('l2') labels.append('Shut-in') if time_big_ev: myVars['l3'], = plt.plot([time_big_ev,time_big_ev], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='dotted', c='k') handels.append('l3') labels.append('Large-Ev') box = myVars[Output_dict_c['Type'][0]].get_position() if len(handels) < 6: myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9]) plt.legend([myVars[ii] for ii in handels], labels, loc='upper center', bbox_to_anchor=(0.5+0.06*first_itr, -0.15), fancybox=True, shadow=True, ncol=len(handels)) elif len(handels) < 13: myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.04*int(len(handels)/2), box.width, box.height * (1 - 0.04*int(len(handels)/2))]) plt.legend([myVars[ii] for ii in handels], labels, loc='upper center', bbox_to_anchor=(0.5+0.1*first_itr, -0.04*int(len(handels)/2)), fancybox=True, shadow=True, ncol=int(len(handels)/2)+1, handleheight=2) else: myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.04*int(len(handels)/2), box.width, box.height * (1 - 0.04*int(len(handels)/2))]) plt.legend([myVars[ii] for ii in handels], labels, loc='upper center', bbox_to_anchor=(0.6+0.1*first_itr, -0.04*int(len(handels)/2)), fancybox=True, shadow=True, ncol=int(len(handels)/2)+1, handleheight=2) plt.title(Model_name) # plt.savefig(cwd+'/Results/'+Model_name+'.png', dpi=300) plt.savefig('PLOT_Mmax_param.png', dpi=300) # plt.show() # Run functions based on the configurations ------------------- # Computing model if time_step_in_hour > time_win_in_hours: raise ValueError('Time steps should be <= time window.') if (time_win_in_hours+time_step_in_hour)*3600 > np.max(Cat[:,4]): raise ValueError('Time window and time steps are like that no more than three computations is possible. Use smaller value for either time window or time step.') time_step = time_step_in_hour*3600 if End_time: if End_time*24*3600 > np.max(Cat[:,4])+24*3600: raise ValueError(f'End_time is longer than the maximum time of catalog, which is {np.ceil(np.max(Cat[:,4])/24/3600)} day!') elif time_step > 0: Times = int(np.floor((End_time*24*3600 - Cat[0,4]) / time_step)) else: Times = 2 time_step = time_win elif time_step > 0: Times = int(np.floor((np.max(Cat[:,4]) - Cat[0,4]) / time_step)) else: Times = 2 time_step = time_win Model_Param_array = np.zeros((Times,c_out)) logger.info("Computing input parameter(s) and the model(s)") Model_Param_array = Feature_in_Time(Model_Param_array) # Plotting if Plot_flag > 0: if Model_Param_array.any(): logger.info("Plotting results") Plot_feature(Model_Param_array, Output_dict) else: logger.info("Model_Param_array is empty or not enough values to plot. Check 'csv' file.") # Saving logger.info("Saving results") # Save models and parameters in def OneLineHeader(loc): l = Output_dict['Type'][loc] if Output_dict['b_method'][loc]: if Output_dict['label'][loc] != 'b_value' and Output_dict['label'][loc] != 'b_std_err' and Output_dict['label'][loc] != 'M_std_err': l += '_'+Output_dict['label'][loc] else: l = Output_dict['label'][loc] if Output_dict['cl'][loc]: l += '(b_method: '+Output_dict['b_method'][loc]+', q='+str(Output_dict['cl'][loc])+')' else: l += '(b_method: '+Output_dict['b_method'][loc]+')' elif Output_dict['label'][loc]: l += '('+ Output_dict['label'][loc] +')' return l db_df = pd.DataFrame(data = Model_Param_array, columns = [OneLineHeader(i) for i in range(len(Output_dict['Type']))]) db_df = db_df.round(4) # db_df.to_excel(cwd+'/Results/'+Full_feat_name+'.xlsx') db_df.to_csv('DATA_Mmax_param.csv', sep=';', index=False) # np.savez_compressed(cwd+'/Results/'+Full_feat_name+'.npz', Model_Param_array = Model_Param_array) # change the name! if __name__=='__main__': # Input_catalog = 'Data/Cooper_Basin_Catalog_HAB_1_2003_Reprocessed.mat' # Input_injection_rate = 'Data/Cooper_Basin_HAB_1_2003_Injection_Rate.mat' # Model_index = 4 # b_value_type = 'TGR' main(Input_catalog, Input_injection_rate, time_win_in_hours, time_step_in_hour, time_win_type, End_time, ev_limit, Inpar, time_inj, time_shut_in, time_big_ev, Model_index, Mc, Mu, G, ssd, C, b_value_type, cl, mag_name)