commit 0df2172bf0eeb5a2a9748d77830a81f7d8420f0d Author: Joanna Kocot Date: Tue Oct 22 13:46:54 2024 +0200 Populated app repository with code and licence diff --git a/LICENCE.txt b/LICENCE.txt new file mode 100644 index 0000000..fe463e0 --- /dev/null +++ b/LICENCE.txt @@ -0,0 +1,408 @@ +Attribution-NonCommercial 4.0 International + +======================================================================= + +Creative Commons Corporation ("Creative Commons") is not a law firm and +does not provide legal services or legal advice. Distribution of +Creative Commons public licenses does not create a lawyer-client or +other relationship. Creative Commons makes its licenses and related +information available on an "as-is" basis. Creative Commons gives no +warranties regarding its licenses, any material licensed under their +terms and conditions, or any related information. Creative Commons +disclaims all liability for damages resulting from their use to the +fullest extent possible. + +Using Creative Commons Public Licenses + +Creative Commons public licenses provide a standard set of terms and +conditions that creators and other rights holders may use to share +original works of authorship and other material subject to copyright +and certain other rights specified in the public license below. The +following considerations are for informational purposes only, are not +exhaustive, and do not form part of our licenses. + + Considerations for licensors: Our public licenses are + intended for use by those authorized to give the public + permission to use material in ways otherwise restricted by + copyright and certain other rights. Our licenses are + irrevocable. Licensors should read and understand the terms + and conditions of the license they choose before applying it. + Licensors should also secure all rights necessary before + applying our licenses so that the public can reuse the + material as expected. Licensors should clearly mark any + material not subject to the license. This includes other CC- + licensed material, or material used under an exception or + limitation to copyright. More considerations for licensors: + wiki.creativecommons.org/Considerations_for_licensors + + Considerations for the public: By using one of our public + licenses, a licensor grants the public permission to use the + licensed material under specified terms and conditions. If + the licensor's permission is not necessary for any reason--for + example, because of any applicable exception or limitation to + copyright--then that use is not regulated by the license. Our + licenses grant only permissions under copyright and certain + other rights that a licensor has authority to grant. Use of + the licensed material may still be restricted for other + reasons, including because others have copyright or other + rights in the material. A licensor may make special requests, + such as asking that all changes be marked or described. + Although not required by our licenses, you are encouraged to + respect those requests where reasonable. More considerations + for the public: + wiki.creativecommons.org/Considerations_for_licensees + +======================================================================= + +Creative Commons Attribution-NonCommercial 4.0 International Public +License + +By exercising the Licensed Rights (defined below), You accept and agree +to be bound by the terms and conditions of this Creative Commons +Attribution-NonCommercial 4.0 International Public License ("Public +License"). To the extent this Public License may be interpreted as a +contract, You are granted the Licensed Rights in consideration of Your +acceptance of these terms and conditions, and the Licensor grants You +such rights in consideration of benefits the Licensor receives from +making the Licensed Material available under these terms and +conditions. + + +Section 1 -- Definitions. + + a. Adapted Material means material subject to Copyright and Similar + Rights that is derived from or based upon the Licensed Material + and in which the Licensed Material is translated, altered, + arranged, transformed, or otherwise modified in a manner requiring + permission under the Copyright and Similar Rights held by the + Licensor. For purposes of this Public License, where the Licensed + Material is a musical work, performance, or sound recording, + Adapted Material is always produced where the Licensed Material is + synched in timed relation with a moving image. + + b. Adapter's License means the license You apply to Your Copyright + and Similar Rights in Your contributions to Adapted Material in + accordance with the terms and conditions of this Public License. + + c. Copyright and Similar Rights means copyright and/or similar rights + closely related to copyright including, without limitation, + performance, broadcast, sound recording, and Sui Generis Database + Rights, without regard to how the rights are labeled or + categorized. For purposes of this Public License, the rights + specified in Section 2(b)(1)-(2) are not Copyright and Similar + Rights. + d. Effective Technological Measures means those measures that, in the + absence of proper authority, may not be circumvented under laws + fulfilling obligations under Article 11 of the WIPO Copyright + Treaty adopted on December 20, 1996, and/or similar international + agreements. + + e. Exceptions and Limitations means fair use, fair dealing, and/or + any other exception or limitation to Copyright and Similar Rights + that applies to Your use of the Licensed Material. + + f. Licensed Material means the artistic or literary work, database, + or other material to which the Licensor applied this Public + License. + + g. Licensed Rights means the rights granted to You subject to the + terms and conditions of this Public License, which are limited to + all Copyright and Similar Rights that apply to Your use of the + Licensed Material and that the Licensor has authority to license. + + h. Licensor means the individual(s) or entity(ies) granting rights + under this Public License. + + i. NonCommercial means not primarily intended for or directed towards + commercial advantage or monetary compensation. For purposes of + this Public License, the exchange of the Licensed Material for + other material subject to Copyright and Similar Rights by digital + file-sharing or similar means is NonCommercial provided there is + no payment of monetary compensation in connection with the + exchange. + + j. Share means to provide material to the public by any means or + process that requires permission under the Licensed Rights, such + as reproduction, public display, public performance, distribution, + dissemination, communication, or importation, and to make material + available to the public including in ways that members of the + public may access the material from a place and at a time + individually chosen by them. + + k. Sui Generis Database Rights means rights other than copyright + resulting from Directive 96/9/EC of the European Parliament and of + the Council of 11 March 1996 on the legal protection of databases, + as amended and/or succeeded, as well as other essentially + equivalent rights anywhere in the world. + + l. You means the individual or entity exercising the Licensed Rights + under this Public License. Your has a corresponding meaning. + + +Section 2 -- Scope. + + a. License grant. + + 1. Subject to the terms and conditions of this Public License, + the Licensor hereby grants You a worldwide, royalty-free, + non-sublicensable, non-exclusive, irrevocable license to + exercise the Licensed Rights in the Licensed Material to: + + a. reproduce and Share the Licensed Material, in whole or + in part, for NonCommercial purposes only; and + + b. produce, reproduce, and Share Adapted Material for + NonCommercial purposes only. + + 2. Exceptions and Limitations. For the avoidance of doubt, where + Exceptions and Limitations apply to Your use, this Public + License does not apply, and You do not need to comply with + its terms and conditions. + + 3. Term. The term of this Public License is specified in Section + 6(a). + + 4. Media and formats; technical modifications allowed. The + Licensor authorizes You to exercise the Licensed Rights in + all media and formats whether now known or hereafter created, + and to make technical modifications necessary to do so. The + Licensor waives and/or agrees not to assert any right or + authority to forbid You from making technical modifications + necessary to exercise the Licensed Rights, including + technical modifications necessary to circumvent Effective + Technological Measures. For purposes of this Public License, + simply making modifications authorized by this Section 2(a) + (4) never produces Adapted Material. + + 5. Downstream recipients. + + a. Offer from the Licensor -- Licensed Material. Every + recipient of the Licensed Material automatically + receives an offer from the Licensor to exercise the + Licensed Rights under the terms and conditions of this + Public License. + + b. No downstream restrictions. You may not offer or impose + any additional or different terms or conditions on, or + apply any Effective Technological Measures to, the + Licensed Material if doing so restricts exercise of the + Licensed Rights by any recipient of the Licensed + Material. + + 6. No endorsement. Nothing in this Public License constitutes or + may be construed as permission to assert or imply that You + are, or that Your use of the Licensed Material is, connected + with, or sponsored, endorsed, or granted official status by, + the Licensor or others designated to receive attribution as + provided in Section 3(a)(1)(A)(i). + + b. Other rights. + + 1. Moral rights, such as the right of integrity, are not + licensed under this Public License, nor are publicity, + privacy, and/or other similar personality rights; however, to + the extent possible, the Licensor waives and/or agrees not to + assert any such rights held by the Licensor to the limited + extent necessary to allow You to exercise the Licensed + Rights, but not otherwise. + + 2. Patent and trademark rights are not licensed under this + Public License. + + 3. To the extent possible, the Licensor waives any right to + collect royalties from You for the exercise of the Licensed + Rights, whether directly or through a collecting society + under any voluntary or waivable statutory or compulsory + licensing scheme. In all other cases the Licensor expressly + reserves any right to collect such royalties, including when + the Licensed Material is used other than for NonCommercial + purposes. + + +Section 3 -- License Conditions. + +Your exercise of the Licensed Rights is expressly made subject to the +following conditions. + + a. Attribution. + + 1. If You Share the Licensed Material (including in modified + form), You must: + + a. retain the following if it is supplied by the Licensor + with the Licensed Material: + + i. identification of the creator(s) of the Licensed + Material and any others designated to receive + attribution, in any reasonable manner requested by + the Licensor (including by pseudonym if + designated); + + ii. a copyright notice; + + iii. a notice that refers to this Public License; + + iv. a notice that refers to the disclaimer of + warranties; + + v. a URI or hyperlink to the Licensed Material to the + extent reasonably practicable; + + b. indicate if You modified the Licensed Material and + retain an indication of any previous modifications; and + + c. indicate the Licensed Material is licensed under this + Public License, and include the text of, or the URI or + hyperlink to, this Public License. + + 2. You may satisfy the conditions in Section 3(a)(1) in any + reasonable manner based on the medium, means, and context in + which You Share the Licensed Material. For example, it may be + reasonable to satisfy the conditions by providing a URI or + hyperlink to a resource that includes the required + information. + + 3. If requested by the Licensor, You must remove any of the + information required by Section 3(a)(1)(A) to the extent + reasonably practicable. + + 4. If You Share Adapted Material You produce, the Adapter's + License You apply must not prevent recipients of the Adapted + Material from complying with this Public License. + + +Section 4 -- Sui Generis Database Rights. + +Where the Licensed Rights include Sui Generis Database Rights that +apply to Your use of the Licensed Material: + + a. for the avoidance of doubt, Section 2(a)(1) grants You the right + to extract, reuse, reproduce, and Share all or a substantial + portion of the contents of the database for NonCommercial purposes + only; + + b. if You include all or a substantial portion of the database + contents in a database in which You have Sui Generis Database + Rights, then the database in which You have Sui Generis Database + Rights (but not its individual contents) is Adapted Material; and + + c. You must comply with the conditions in Section 3(a) if You Share + all or a substantial portion of the contents of the database. + +For the avoidance of doubt, this Section 4 supplements and does not +replace Your obligations under this Public License where the Licensed +Rights include other Copyright and Similar Rights. + + +Section 5 -- Disclaimer of Warranties and Limitation of Liability. + + a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE + EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS + AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF + ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS, + IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION, + WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR + PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS, + ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT + KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT + ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU. + + b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE + TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION, + NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT, + INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES, + COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR + USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN + ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR + DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR + IN PART, THIS LIMITATION MAY NOT APPLY TO YOU. + + c. The disclaimer of warranties and limitation of liability provided + above shall be interpreted in a manner that, to the extent + possible, most closely approximates an absolute disclaimer and + waiver of all liability. + + +Section 6 -- Term and Termination. + + a. This Public License applies for the term of the Copyright and + Similar Rights licensed here. However, if You fail to comply with + this Public License, then Your rights under this Public License + terminate automatically. + + b. Where Your right to use the Licensed Material has terminated under + Section 6(a), it reinstates: + + 1. automatically as of the date the violation is cured, provided + it is cured within 30 days of Your discovery of the + violation; or + + 2. upon express reinstatement by the Licensor. + + For the avoidance of doubt, this Section 6(b) does not affect any + right the Licensor may have to seek remedies for Your violations + of this Public License. + + c. For the avoidance of doubt, the Licensor may also offer the + Licensed Material under separate terms or conditions or stop + distributing the Licensed Material at any time; however, doing so + will not terminate this Public License. + + d. Sections 1, 5, 6, 7, and 8 survive termination of this Public + License. + + +Section 7 -- Other Terms and Conditions. + + a. The Licensor shall not be bound by any additional or different + terms or conditions communicated by You unless expressly agreed. + + b. Any arrangements, understandings, or agreements regarding the + Licensed Material not stated herein are separate from and + independent of the terms and conditions of this Public License. + + +Section 8 -- Interpretation. + + a. For the avoidance of doubt, this Public License does not, and + shall not be interpreted to, reduce, limit, restrict, or impose + conditions on any use of the Licensed Material that could lawfully + be made without permission under this Public License. + + b. To the extent possible, if any provision of this Public License is + deemed unenforceable, it shall be automatically reformed to the + minimum extent necessary to make it enforceable. If the provision + cannot be reformed, it shall be severed from this Public License + without affecting the enforceability of the remaining terms and + conditions. + + c. No term or condition of this Public License will be waived and no + failure to comply consented to unless expressly agreed to by the + Licensor. + + d. Nothing in this Public License constitutes or may be interpreted + as a limitation upon, or waiver of, any privileges and immunities + that apply to the Licensor or You, including from the legal + processes of any jurisdiction or authority. + +======================================================================= + +Creative Commons is not a party to its public +licenses. Notwithstanding, Creative Commons may elect to apply one of +its public licenses to material it publishes and in those instances +will be considered the “Licensor.” The text of the Creative Commons +public licenses is dedicated to the public domain under the CC0 Public +Domain Dedication. Except for the limited purpose of indicating that +material is shared under a Creative Commons public license or as +otherwise permitted by the Creative Commons policies published at +creativecommons.org/policies, Creative Commons does not authorize the +use of the trademark "Creative Commons" or any other trademark or logo +of Creative Commons without its prior written consent including, +without limitation, in connection with any unauthorized modifications +to any of its public licenses or any other arrangements, +understandings, or agreements concerning use of licensed material. For +the avoidance of doubt, this paragraph does not form part of the +public licenses. + +Creative Commons may be contacted at creativecommons.org. + diff --git a/Mmax.py b/Mmax.py new file mode 100644 index 0000000..ba2f8c1 --- /dev/null +++ b/Mmax.py @@ -0,0 +1,713 @@ +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 + +def 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): + """ + 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 = 1 + + # Setting up the logfile-------------------------------- + logging.basicConfig(filename="application.log", + filemode='a', + format='%(asctime)s,%(msecs)d %(name)s %(levelname)s %(message)s', + datefmt='%H:%M:%S', + level=logging.INFO) + + logger = logging.getLogger(__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 ,5 ], + '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 + 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'] + # 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= coloring[itr_loc[0]], ls='--') + 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 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= coloring[itr_loc[j]], ls='--') + 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[j]+1 <= len(itr_loc) 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= coloring[itr_loc[0]], ls='--') + 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= coloring[itr_loc[j]], ls = '--') + 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 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+3*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!') + Times = int(np.floor((End_time*24*3600 - Cat[0,4]) / time_step)) + else: + Times = int(np.floor((np.max(Cat[:,4]) - Cat[0,4]) / time_step)) + 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.") + raise ValueError("No model is generated since no event in all time windows is available due to imporoper values of either time window, minimum limit of events or completeness magnitude!") + + # 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..5c6a8c7 --- /dev/null +++ b/maxmagnitude_wrapper.py @@ -0,0 +1,55 @@ +# -*- 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) + # Parameters with action='append' are array inputs. We can pass array by using a parameters multiple time. Eg.: + # --Inpar=A --Inpar=B --Inpar=C pass ['A', 'B', 'C'] array as script input. + parser.add_argument("--Inpar", help="Input parameters of the model to be saved and plotted.", type=str, action='append', default=None) + parser.add_argument("--time_inj", help="Injection time for plotting (in day).", type=int, default=None) + parser.add_argument("--time_shut_in", help="Shut-in time for plotting (in day).", type=int, default=None) + parser.add_argument("--time_big_ev", help="Big event time for plotting (in day).", type=int, default=None) + 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..af05e68 --- /dev/null +++ b/util/M_max_models.py @@ -0,0 +1,216 @@ +# ----------------- +# Copyright © 2024 ACK Cyfronet AGH, Poland. +# ----------------- + +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'.") + 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() + + 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