From b3be576e99d1199d9c92ded1d48a2c32e7d488f9 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Sun, 8 Dec 2024 23:14:50 +0900 Subject: [PATCH] Add MaxMagnitudeDPModels application code --- Mmax.py | 723 ++++++++++++++++++++++++++++++++++++++++++++++ maxmagnitude_wrapper.py | 49 ++++ util/CandidateEventsTS.py | 43 +++ util/Find_idx4Time.py | 14 + util/M_max_models.py | 217 ++++++++++++++ util/base_logger.py | 62 ++++ 6 files changed, 1108 insertions(+) create mode 100644 Mmax.py create mode 100644 maxmagnitude_wrapper.py create mode 100644 util/CandidateEventsTS.py create mode 100644 util/Find_idx4Time.py create mode 100644 util/M_max_models.py create mode 100644 util/base_logger.py diff --git a/Mmax.py b/Mmax.py new file mode 100644 index 0000000..a3132a4 --- /dev/null +++ b/Mmax.py @@ -0,0 +1,723 @@ +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) diff --git a/maxmagnitude_wrapper.py b/maxmagnitude_wrapper.py new file mode 100644 index 0000000..de02925 --- /dev/null +++ b/maxmagnitude_wrapper.py @@ -0,0 +1,49 @@ +# -*- coding: utf-8 -*- + +# ----------------- +# Copyright © 2024 ACK Cyfronet AGH, Poland. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# This work was partially funded by DT-GEO Project. +# ----------------- + +import sys +import argparse +from Mmax import main as Mmax + +def main(argv): + parser = argparse.ArgumentParser() + parser.add_argument("Input_catalog", help="Input catalog: path to input file of type 'catalog'") + parser.add_argument("Input_injection_rate", help="Input injection rate: path to input file of type 'injection_rate'") + parser.add_argument("--time_win_in_hours", help="Time window length (in hours- backward from the current time).", type=int, default=6) + parser.add_argument("--time_step_in_hour", help="Time interval for computation (in hours).", type=int, default=3) + parser.add_argument("--time_win_type", help="Time window type for computation.", type=int, default=0) + parser.add_argument("--End_time", help="End time of the computations (in day).", type=int, default=None) + parser.add_argument("--ev_limit", help="Minimum events number required for model computation.", type=int, default=20) + parser.add_argument("--Model_index", help="Model index: parameter of type 'INTEGER'", type=int) + parser.add_argument("--Mc", help="Completeness magnitude.", type=float, default=0.8) + parser.add_argument("--Mu", help="Friction coefficient.", type=float, default=0.6, required=False) + parser.add_argument("--G", help="Shear modulus of reservoir (in Pa).", type=float, default=35000000000) + parser.add_argument("--ssd", help="Static stress drop (in Pa).", type=float, default=3000000) + parser.add_argument("--C", help="Geometrical constant.", type=float, default=0.95) + parser.add_argument("--b_value_type", help="b-value type: parameter of type 'TEXT'", action='append') + parser.add_argument("--cl", help="Confidence level in van der Elst model.", type=float, action='append') + parser.add_argument("--mag_name", help="Magnitude column name", type=str) + + args = parser.parse_args() + Mmax(**vars(args)) + return + +if __name__ == '__main__': + main(sys.argv) \ No newline at end of file diff --git a/util/CandidateEventsTS.py b/util/CandidateEventsTS.py new file mode 100644 index 0000000..cf8dd5e --- /dev/null +++ b/util/CandidateEventsTS.py @@ -0,0 +1,43 @@ +import numpy as np + +class CandidateEventsTS: + def __init__(self, data, current_time, Mc, time_win, space_win=None): + assert time_win > 0, f"Time windows is {time_win}, which should be a positive number" + + self.data = data + self.current_time = current_time + self.Mc = Mc + self.time_win = time_win + self.space_win = space_win + + def filter_by_time(self): + indx = np.where((self.data[:, 4] > (self.current_time - self.time_win)) & (self.data[:, 4] <= self.current_time))[0] + if len(indx) > 0: + self.data = self.data[indx, :] + else: + self.data = [] + + def filter_by_magnitude(self): + if self.Mc: + indx = np.where(self.data[:, 5] > self.Mc)[0] + if len(indx) > 0: + self.data = self.data[indx, :] + else: + self.data = [] + + def filter_by_space(self): + dist = np.sqrt(np.sum((self.data[:, 1:4] - self.data[-1, 1:4]) ** 2, axis=1)) + indx = np.where(dist < self.space_win)[0] + if len(indx) > 0: + self.data = self.data[indx, :] + else: + self.data = [] + + def filter_data(self): + self.filter_by_time() + if len(self.data) > 0: + self.filter_by_magnitude() + if len(self.data) > 0 and self.space_win: + self.filter_by_space() + + return self.data diff --git a/util/Find_idx4Time.py b/util/Find_idx4Time.py new file mode 100644 index 0000000..f6cf838 --- /dev/null +++ b/util/Find_idx4Time.py @@ -0,0 +1,14 @@ +import numpy as np +def Find_idx4Time(In_mat, t): + # In_mat: time array + # t = target time + In_mat = np.array(In_mat) + t = np.array(t) + if len(np.shape(t)) == 0: + return np.where(abs(In_mat - t) <= min(abs(In_mat - t)))[0][0] + else: + In_mat = In_mat.reshape((len(In_mat), 1)) + t = t.reshape((1,len(t))) + target_time = np.matmul(np.ones((len(In_mat),1)), t) + diff_mat = target_time - In_mat + return np.where(abs(diff_mat) <= np.min(abs(diff_mat), axis = 0)) \ No newline at end of file diff --git a/util/M_max_models.py b/util/M_max_models.py new file mode 100644 index 0000000..c2fccc5 --- /dev/null +++ b/util/M_max_models.py @@ -0,0 +1,217 @@ +import numpy as np + +class M_max_models: + def __init__(self, data = None, f_name = None, time_win = None, space_win = None, + Mc = None, b_method = None, num_bootstraps = None, + G = None, Mu = None, + dv = None, Mo = None, SER = None, + cl = None, + ssd = None, C = None, + ): + + self.data = data # Candidate data table: 2darray n x m, for n events and m clonums: x, y, z, t, mag + self.f_name = f_name # Feature's name to be calculated, check: def ComputeFeaure(self) + self.time_win = time_win # Time window whihc a feature is computed in + self.space_win = space_win # Space window ... + + self.Mc = Mc # Magnitude of completeness for computing b-positive + self.b_method = b_method # list of b_methods + self.num_bootstraps = num_bootstraps # Num of bootstraps for standard error estimation of b-value + + self.G = G # Shear modulus + self.Mu = Mu # Friction coefficient + self.dv = dv # Injected fluid + self.SER = SER + + self.Mo = Mo # Cumulative moment magnitude + self.cl = cl # Confidence level + + self.ssd = ssd # Static stress drop (Shapiro et al. 2013) + self.C = C # Geometrical constant (Shapiro et al. 2013) + + + + + def b_value(self, b_flag): + if b_flag == '1': + return 1, None + + # maximum-likelihood estimate (MLE) of b (Deemer & Votaw 1955; Aki 1965; Kagan 2002): + elif b_flag == 'b': + X = self.data[np.where(self.data[:,-1]>self.Mc)[0],:] + if X.shape[0] > 0: + b = 1/((np.mean(X[:,-1] - self.Mc))*np.log(10)) + std_error = b/np.sqrt(X.shape[0]) + else: + raise ValueError("All events in the current time window have a magnitude less than 'completeness magnitude'. Use another value either for 'time window', 'minimum number of events' or 'completeness magnitude'. Also check 'time window type'.") + return b, std_error + + # B-positive (van der Elst 2021) + elif b_flag == 'bp': + + # Function to perform bootstrap estimation + def bootstrap_estimate(data, num_bootstraps): + estimates = [] + for _ in range(num_bootstraps): + # Generate bootstrap sample + bootstrap_sample = np.random.choice(data, size=len(data), replace=True) + # Perform maximum likelihood estimation on bootstrap sample + diff_mat = np.diff(bootstrap_sample) + diff_mat = diff_mat[np.where(diff_mat>0)[0]] + estimate = 1/((np.mean(diff_mat - np.min(diff_mat)))*np.log(10)) + estimates.append(estimate) + return np.array(estimates) + + diff_mat = np.diff(self.data[:,-1]) + diff_mat = diff_mat[np.where(diff_mat>0)[0]] + bp = 1/((np.mean(diff_mat - np.min(diff_mat)))*np.log(10)) + bootstrap_estimates = bootstrap_estimate(diff_mat, self.num_bootstraps) + std_error = np.std(bootstrap_estimates, axis=0) + return bp, std_error + + # Tapered Gutenberg_Richter (TGR) distribution (Kagan 2002) + elif b_flag == 'TGR': + from scipy.optimize import minimize + + # The logarithm of the likelihood function for the TGR distribution (Kagan 2002) + def log_likelihood(params, data): + beta, Mcm = params + n = len(data) + Mt = np.min(data) + l = n*beta*np.log(Mt)+1/Mcm*(n*Mt-np.sum(data))-beta*np.sum(np.log(data))+np.sum(np.log([(beta/data[i]+1/Mcm) for i in range(len(data))])) + return -l + X = self.data[np.where(self.data[:,-1]>self.Mc)[0],:] + M = 10**(1.5*X[:,-1]+9.1) + initial_guess = [0.5, np.max(M)] + bounds = [(0.0, None), (np.max(M), None)] + + # Minimize the negative likelihood function for beta and maximum moment + result = minimize(log_likelihood, initial_guess, args=(M,), bounds=bounds, method='L-BFGS-B', + options={'gtol': 1e-12, 'disp': False}) + beta_opt, Mcm_opt = result.x + eta = 1/Mcm_opt + S = M/np.min(M) + dldb2 = -np.sum([1/(beta_opt-eta*S[i])**2 for i in range(len(S))]) + dldbde = -np.sum([S[i]/(beta_opt-eta*S[i])**2 for i in range(len(S))]) + dlde2 = -np.sum([S[i]**2/(beta_opt-eta*S[i])**2 for i in range(len(S))]) + std_error_beta = 1/np.sqrt(dldb2*dlde2-dldbde**2)*np.sqrt(-dlde2) + return beta_opt*1.5, std_error_beta*1.5 + + def McGarr(self): + b_value, b_stderr = self.b_value(self.b_method) + B = 2/3*b_value + if B < 1: + sigma_m = ((1-B)/B)*(2*self.Mu)*(5*self.G)/3*self.dv + Mmax = (np.log10(sigma_m)-9.1)/1.5 + if b_stderr: + Mmax_stderr = b_stderr/np.abs(np.log(10)*(1.5*b_value-b_value**2)) + else: + Mmax_stderr = None + else: + Mmax = None + Mmax_stderr = None + + return b_value, b_stderr, Mmax, Mmax_stderr + + def Hallo(self): + b_value, b_stderr = self.b_value(self.b_method) + B = 2/3*b_value + if b_value < 1.5: + sigma_m = self.SER*((1-B)/B)*(2*self.Mu)*(5*self.G)/3*self.dv + Mmax = (np.log10(sigma_m)-9.1)/1.5 + if b_stderr: + Mmax_stderr = self.SER*b_stderr/np.abs(np.log(10)*(1.5*b_value-b_value**2)) + else: + Mmax_stderr = None + else: + Mmax = None + Mmax_stderr = None + + return b_value, b_stderr, Mmax, Mmax_stderr + + def Li(self): + sigma_m = self.SER*2*self.G*self.dv - self.Mo + Mmax = (np.log10(sigma_m)-9.1)/1.5 + if Mmax < 0: + return None + else: + return Mmax + + def van_der_Elst(self): + b_value, b_stderr = self.b_value(self.b_method) + # Seismogenic_Index + X = self.data + si = np.log10(X.shape[0]) - np.log10(self.dv) + b_value*self.Mc + if b_stderr: + si_stderr = self.Mc*b_stderr + else: + si_stderr = None + + Mmax = (si + np.log10(self.dv))/b_value - np.log10(X.shape[0]*(1-self.cl**(1/X.shape[0])))/b_value + if b_stderr: + Mmax_stderr = (np.log10(X.shape[0]) + np.log10(X.shape[0]*(1-self.cl**(1/X.shape[0]))))*b_stderr + else: + Mmax_stderr = None + + return b_value, b_stderr, si, si_stderr, Mmax, Mmax_stderr + + def L_Shapiro(self): + from scipy.stats import chi2 + + X = self.data[np.isfinite(self.data[:,1]),1:4] + # Parameters + STD = 2.0 # 2 standard deviations + conf = 2 * chi2.cdf(STD, 2) - 1 # covers around 95% of population + scalee = chi2.ppf(conf, 2) # inverse chi-squared with dof=#dimensions + + # Center the data + Mu = np.mean(X, axis=0) + X0 = X - Mu + + # Covariance matrix + Cov = np.cov(X0, rowvar=False) * scalee + + # Eigen decomposition + D, V = np.linalg.eigh(Cov) + order = np.argsort(D)[::-1] + D = D[order] + V = V[:, order] + + # Compute radii + VV = V * np.sqrt(D) + R1 = np.sqrt(VV[0, 0]**2 + VV[1, 0]**2 + VV[2, 0]**2) + R2 = np.sqrt(VV[0, 1]**2 + VV[1, 1]**2 + VV[2, 1]**2) + R3 = np.sqrt(VV[0, 2]**2 + VV[1, 2]**2 + VV[2, 2]**2) + + L = (1/3*(1/R1**3+1/R2**3+1/R3**3))**(-1/3) + + return R1, R2, R3, L + + def Shapiro(self): + R1, R2, R3, L = self.L_Shapiro() + Sh_lmax = np.log10((2*R1)**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5 + Sh_lint = np.log10((2*R2)**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5 + Sh_lmin = np.log10((2*R3)**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5 + Sh_lavg = np.log10((2*L)**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5 + + return Sh_lmax, Sh_lint, Sh_lmin, Sh_lavg + # return R1, R2, R3, L, np.log10(R3**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5 + + def All_models(self): + + return self.McGarr()[2], self.Hallo()[2], self.Li(), self.van_der_Elst()[-2], self.Shapiro()[-1] + + + def ComputeModel(self): + if self.f_name == 'max_mcg': + return self.McGarr() + if self.f_name == 'max_hlo': + return self.Hallo() + if self.f_name == 'max_li': + return self.Li() + if self.f_name == 'max_vde': + return self.van_der_Elst() + if self.f_name == 'max_shp': + return self.Shapiro() + if self.f_name == 'max_all': + return self.All_models() \ No newline at end of file diff --git a/util/base_logger.py b/util/base_logger.py new file mode 100644 index 0000000..188b72f --- /dev/null +++ b/util/base_logger.py @@ -0,0 +1,62 @@ +# +# ----------------- +# Copyright © 2024 ACK Cyfronet AGH, Poland. +# ----------------- +# +import os +import logging + +def getDefaultLogger(name): + """ + Retrieves or creates a logger with the specified name and sets it up with a file handler. + + The logger is configured to write log messages to the file path specified by the + 'APP_LOG_FILE' environment variable. If the environment variable is not set, + the logger will write to the file 'base-logger-log.log' in the current + working directory. The logger uses the 'INFO' level as the default logging level + and writes log entries in the following format: + + 'YYYY-MM-DD HH:MM:SS,ms LEVEL logger_name message' + + If the logger does not already have handlers, a file handler is created, and the + logging output is appended to the file. The log format includes the timestamp with + milliseconds, log level, logger name, and the log message. + + Parameters: + ----------- + name : str + The name of the logger. This can be the name of the module or any identifier + that you want to associate with the logger. + + Returns: + -------- + logger : logging.Logger + A logger instance with the specified name. The logger is configured with a + file handler that writes to the file specified by the 'APP_LOG_FILE' + environment variable, or to 'base-logger-log.log' if the environment + variable is not set. + + Example: + -------- + logger = getDefaultLogger(__name__) + logger.info("This is an info message.") + try: + # some code causing exception + except Exception: + logger.exception('An error occurred') + + Notes: + ------ + - The 'APP_LOG_FILE' environment variable should specify the full path to the log file. + - If 'APP_LOG_FILE' is not set, logs will be written to 'base-logger-log.log'. + + """ + logger = logging.getLogger(name) + if not logger.hasHandlers(): + file_handler = logging.FileHandler(os.environ.get('APP_LOG_FILE', 'base-logger-log.log'), mode='a') + formatter = logging.Formatter('%(asctime)s,%(msecs)d %(levelname)s %(name)s %(message)s') + file_handler.setFormatter(formatter) + logger.addHandler(file_handler) + logger.setLevel(logging.INFO) + + return logger \ No newline at end of file