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)