Update src/seismic_hazard_forecasting.py

log more package version numbers

src/seismic_hazard_forecasting.py

@@ -1,5 +1,189 @@
 # -*- coding: utf-8 -*-
+import numpy as np
+import pandas as pd
+import Rbeast as rb
+import matplotlib.dates as mdates 
+from matplotlib.dates import DateFormatter, AutoDateLocator
+import matplotlib.pyplot as plt
+import pip
 
+pip.main(['install', '--verbose', 'igfash==0.2.2'])
+
+from igfash.rate import datenum_to_datetime, av_rates
+
+def apply_beast(act_rate, **kwargs):
+    """
+    Applies BEAST to the smmothed rate data using different smoothing windows.
+    Input
+    act_rate : The activity rate data array to smooth and apply BEAST.
+    Output
+    out : A list of BEAST results for each smoothed rate array.
+    prob : A list of probabilities and change points extracted from BEAST results.
+    """
+    out = []; prob = []
+    # Define the range of smoothing windows based on 15-20% of the size of rate array
+    ll = int(np.ceil(0.15 * len(act_rate)))
+    ul = int(np.ceil(0.20 * len(act_rate)))
+    
+    # Apply BEAST for each smoothed rate array
+    for ii in range(ll, ul + 1):
+        # Apply Gaussian smoothing with window size ~15-20%
+        bb=pd.Series(act_rate).rolling(window=ii, win_type='gaussian', center=True, min_periods=1).mean(std=ii).to_numpy()  
+        # Apply BEAST on the smoothed array discarding periodicity
+        beast_result = rb.beast(bb, season='none', **kwargs)
+        out.append(beast_result)
+        
+    # Extract probabilities and change points from BEAST results
+    # Gets the median number of changepoints or if median = 0 gets 90th percentile
+    for result in out:
+        if result.trend.ncp_median != 0 and not np.isnan(result.trend.ncp_median):
+            ncp_median = int(result.trend.ncp_median)  # No need for np.max if it's scalar
+            cp_probs = result.trend.cpPr[:ncp_median]
+            cps = result.trend.cp[:ncp_median]
+            sorted_indices = np.argsort(cps)
+            sorted_cps = cps[sorted_indices]
+            sorted_probs = cp_probs[sorted_indices]
+        else:
+            if not np.isnan(result.trend.ncp_pct90):
+                ncp_pct90 = int(result.trend.ncp_pct90)
+                cp_probs = result.trend.cpPr[:ncp_pct90]
+                cps = result.trend.cp[:ncp_pct90]
+                sorted_indices = np.argsort(cps)
+                sorted_cps = cps[sorted_indices]
+                sorted_probs = cp_probs[sorted_indices]
+            else:
+                # Optional fallback in case both are NaN
+                cp_probs = []
+                cps = []
+
+#    for result in out:
+#        if result.trend.ncp_median != 0:
+#            ncp_median = int(np.max(result.trend.ncp_median))  #Ensure it's integer
+#            cp_probs = result.trend.cpPr[:ncp_median]
+#            cps = result.trend.cp[:ncp_median]
+#        else:
+#            ncp_pct90 = int(np.max(result.trend.ncp_pct90))  #Ensure it's integer
+#            cp_probs = result.trend.cpPr[:ncp_pct90]
+#            cps = result.trend.cp[:ncp_pct90]
+            
+        # Sort the change points and corresponding probabilities
+        #sorted_indices = np.argsort(cps)
+        #sorted_cps = cps[sorted_indices]
+        #sorted_probs = cp_probs[sorted_indices]
+        
+        # Store the sorted change points and probabilities
+        prob.append(np.column_stack((sorted_probs, sorted_cps)))
+        
+    return out, prob
+
+def calc_bins2(dates, unit, bin_dur, dates_calc, rate_forecast, rate_unc_high, rate_unc_low, multiplicator, filename = "activity_rate.png", **kwargs):
+
+    start_date = dates.min()
+    end_date = dates.max()
+    end_date_dt=datenum_to_datetime(end_date)
+
+    # Unit and duration check, default to 15 DAYS or 12 HOURS
+    # valid_units = ['hours', 'days']
+    # if unit not in valid_units:
+    #     unit = 'days'
+    #     bin_dur = 15
+
+    # Adjust to '12 hours' if data duration is less than 15 days
+    # if (end_date - start_date) < 15 and unit == 'days':
+    #     unit = 'hours'
+    #     bin_dur = 12
+
+    # -------- bin, bin edges, bin midpoint calc -------------
+    # Determine number of bins
+    bin_edges=[end_date]
+    while bin_edges[-1] > start_date:
+        bin_edges.append(bin_edges[-1] - (bin_dur/multiplicator))
+    bin_edges = bin_edges[::-1]
+    bin_edges_dt = [datenum_to_datetime(d) for d in bin_edges]  
+    # Discretize the time data
+    bin_counts, _ = np.histogram(dates, bins=bin_edges)
+    # Calculate the rate
+    act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator) * multiplicator for i, count in enumerate(bin_counts)]
+    
+    # ---------Apply BEAST (Zhao et al., 2019) 10 times--------
+    ppr_all = []
+    for _ in range(10):
+        out, ppr = apply_beast(act_rate, **kwargs)  # Apply BEAST
+        ppr_iter = np.concatenate([p for p in ppr if p is not None], axis=0)  # Concatenate results
+        ppr_all.append(ppr_iter[:, 1])
+    # Concatenate all probabilities from all BEAST runs
+    pprs = np.concatenate(ppr_all, axis=0)
+    pprs = pprs[~np.isnan(pprs)]
+    
+    # -------Changepoint decision and rate calculations--------
+    rt, idx, u_e = av_rates(act_rate, pprs)
+    
+    # --------------------- Line Plot -------------------------
+    fig, ax = plt.subplots(figsize=(14, 5))
+    
+    # Plot activiry rate for all catalog
+    ax.plot(bin_edges_dt[1:len(bin_edges)], act_rate, '-o', linewidth=2, markersize=6)
+    
+    # Plot forecasted value with uncertainties with red line connected to last value
+    next_date = end_date + (bin_dur/multiplicator) 
+    ax.plot(datenum_to_datetime(next_date), rate_forecast, 'ro', label='Forecasted Rate', markersize=6)
+    ax.plot([bin_edges_dt[-1], datenum_to_datetime(next_date)], [act_rate[-1], rate_forecast], 'r-')
+    ax.vlines(datenum_to_datetime(next_date), rate_unc_low, rate_unc_high, colors='r', linewidth=2, label='Bootstrap uncertainty')  # Vertical line for the uncertainty range
+    
+    # Format the x-axis with dates
+    ax.xaxis.set_major_locator(AutoDateLocator())
+    ax.xaxis.set_major_formatter(DateFormatter('%d-%b-%Y'))
+    plt.xticks(rotation=45)
+    plt.title(f'Activity rate (Time Unit: {unit}, Time Window Duration: {bin_dur} {unit})')
+    plt.xlabel('Time (Window Center Date)')
+    plt.ylabel('Activity rate per selected time window')
+    plt.grid(True)
+    
+    # Add Vertical Line and Shaded Area
+    sel_data = bin_edges[-3]
+    sel_data_dt=datenum_to_datetime(sel_data)
+    ax.axvline(sel_data_dt, color='r', linestyle='--', linewidth=2,  label='Time window used')
+    def convert_dates_to_num(dates):
+        return [mdates.date2num(date) if isinstance(date, datetime) else date for date in dates]
+    
+    # Shade the area to the right of the vertical line
+    ylim = ax.get_ylim()
+    
+    # Add the shaded area using the Polygon function
+    ax.fill_between([sel_data_dt, end_date_dt], 0.0, ylim[1]+0.1, color='r', alpha=0.5)
+    
+    # Add rates and changepoints from BEAST - changepoint belongs to NEXT period
+    if 'rt' in locals():
+        for i in range(len(idx)):
+            # Plot first period
+            if i == 0 and len(idx) > 1:
+                ax.plot(bin_edges_dt[1:len(bin_edges)][:int(u_e[idx[i]])], [rt[i]] * int(u_e[idx[i]]), linewidth=2)
+            # Plot Last and second to last periods when more than 1 changepoints
+            elif i == len(idx) - 1 and len(idx) > 1:
+                ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]])], 
+                        [rt[i]] * (int(u_e[idx[i]]) - int(u_e[idx[i - 1]])), linewidth=2)
+                ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i]]):], 
+                        [rt[i + 1]] * (len(act_rate)  - int(u_e[idx[i]])), linewidth=2)
+            # Plot first and last period if only one changepoint
+            elif i == 0 and len(idx) == 1:
+                ax.plot(bin_edges_dt[1:len(bin_edges)][:int(u_e[idx[i]])], [rt[i]] * int(u_e[idx[i]]), linewidth=2)
+                ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i]]):], 
+                         [rt[i + 1]] * (len(act_rate) - int(u_e[idx[i]])), linewidth=2)
+            # Plot intermediate value if it's only one 
+            elif len(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[1 - 1]]):int(u_e[idx[1]])])==1:
+                nnn = multiplicator / bin_dur
+                ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]]) + 1], 
+                        [rt[i]] * (int(u_e[idx[i]]) + 1 - int(u_e[idx[i - 1]])), linewidth=2)                
+            # Plot intermediate periods    
+            else:
+                ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]])], 
+                        [rt[i]] * (int(u_e[idx[i]]) - int(u_e[idx[i - 1]])), linewidth=2)
+
+    plt.legend(loc='best')
+    plt.savefig(filename)
+    plt.show()
+    
+    return act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e
 
 def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max,
          m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration,
@@ -45,21 +229,18 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     """
 
     import sys
+    from importlib.metadata import version
     import logging
     from base_logger import getDefaultLogger
     from timeit import default_timer as timer
     from math import ceil, floor, isnan
     import numpy as np
-    import scipy
-    import obspy
     import dask
     from dask.diagnostics import ProgressBar  # use Dask progress bar
     import kalepy as kale
     import utm
     from skimage.transform import resize
     import psutil
-    import openquake.engine
-    import igfash
     from igfash.io import read_mat_cat, read_mat_m, read_mat_pdf, read_csv
     from igfash.window import win_CTL, win_CNE
     import igfash.kde as kde
@@ -105,11 +286,13 @@ verbose: {verbose}")
 
     # print key package version numbers
     logger.debug(f"Python version {sys.version}")
-    logger.debug(f"Numpy version {np.__version__}")
+    logger.debug(f"Numpy version {version('numpy')}")
-    logger.debug(f"Scipy version {scipy.__version__}")
+    logger.debug(f"Scipy version {version('scipy')}")
-    logger.debug(f"Obspy version {obspy.__version__}")
+    logger.debug(f"Obspy version {version('obspy')}")
-    logger.debug(f"Openquake version {openquake.engine.__version__}")
+    logger.debug(f"Openquake version {version('openquake.engine')}")
-    logger.debug(f"Igfash version {igfash.__version__}")
+    logger.debug(f"Igfash version {version('igfash')}")
+    logger.debug(f"Rbeast version {version('rbeast')}")
+    logger.debug(f"Dask version {version('dask')}")
 
     # print number of cpu cores available
     ncpu = psutil.cpu_count(logical=False)
@@ -365,7 +548,7 @@ verbose: {verbose}")
         act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit,
                                                                              time_win_duration, dates_calc,
                                                                              rate_forecast, rate_unc_high, rate_unc_low,
-                                                                             multiplicator, quiet=True)
+                                                                             multiplicator, quiet=True, figsize=(14,9))
 
         # Assign probabilities 
         lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges)