From 65759b86f12ac90b0af8e3b31785b4d16729326d Mon Sep 17 00:00:00 2001 From: ftong Date: Tue, 9 Sep 2025 10:56:35 +0200 Subject: [PATCH 01/11] change search interval for PGV to be different than that for PGA/SA --- src/seismic_hazard_forecasting.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 4ebc6ba..41b3871 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -88,9 +88,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max else: logger.setLevel(logging.INFO) - # temporary hard-coded configuration - # exclude_low_fxy = False - exclude_low_fxy = True + exclude_low_fxy = True # skip low probability areas of the map thresh_fxy = 1e-3 # minimum fxy value (location PDF) needed to do PGA estimation (to skip low probability areas); also should scale according to number of grid points # log user selections @@ -125,10 +123,6 @@ verbose: {verbose}") logger.info("No magnitude label of catalog specified, therefore try Mw by default") mag_label = 'Mw' - # if cat_label == None: - # print("No magnitude label of catalog specified, therefore try 'Catalog' by default") - # cat_label='Catalog' - time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog') # check for null magnitude values @@ -459,15 +453,19 @@ verbose: {verbose}") if use_pp: # use dask parallel computing pbar = ProgressBar() pbar.register() - # iter = range(0,len(distances)) iter = indices iml_grid_raw = [] # raw ground motion grids for imt in products: logger.info(f"Estimating {imt}") + if imt == "PGV": + IMT_max = 200 # search interval max for velocity (cm/s) + else: + IMT_max = 2.0 # search interval max for acceleration (g) + imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, - log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, rx_label=i, + log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i, rtol=0.1, use_cython=True) for i in iter] iml = dask.compute(*imls) iml_grid_raw.append(list(iml)) @@ -476,11 +474,17 @@ verbose: {verbose}") iml_grid_raw = [] iter = indices for imt in products: + + if imt == "PGV": + IMT_max = 200 # search interval max for velocity (cm/s) + else: + IMT_max = 2.0 # search interval max for acceleration (g) + iml = [] for i in iter: iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, - IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=True) + IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True) iml.append(iml_i) logger.info(f"Estimated {imt} at rx {i} is {iml_i}") iml_grid_raw.append(iml) -- 2.16.5 From 63351ceb1027885673242c5d4c8a62d2fa646460 Mon Sep 17 00:00:00 2001 From: ftong Date: Tue, 9 Sep 2025 11:03:05 +0200 Subject: [PATCH 02/11] fix weighting option selection --- src/seismic_hazard_forecasting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 41b3871..0703b50 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -252,7 +252,7 @@ verbose: {verbose}") # %% compute KDE and extract PDF start = timer() - if xy_win_method == "TW": + if xy_win_method: logger.info("Time weighting function selected") x_weights = np.linspace(0, 15, len(t_windowed)) -- 2.16.5 From a1c0ae36bb1a29c7d547b5c315360ae7c3741479 Mon Sep 17 00:00:00 2001 From: ftong Date: Tue, 9 Sep 2025 14:41:02 +0200 Subject: [PATCH 03/11] set a minimum number of computed grid values to trigger upscaling of grid image --- src/seismic_hazard_forecasting.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 0703b50..3cd03c6 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -519,10 +519,14 @@ verbose: {verbose}") dtype=np.float64) # this reduces values to 8 decimal places iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes - # upscale the grid - up_factor = 4 - iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), - mode='reflect', anti_aliasing=False) + # upscale the grid if there are at least 10 grid values + if np.count_nonzero(iml_grid_tmp) >= 10: + up_factor = 4 + iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), + mode='reflect', anti_aliasing=False) + else: + iml_grid_hd = iml_grid_tmp + iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan # trim edges so the grid is not so blocky -- 2.16.5 From 5a1f43d6cd27c48fecded1a47d6e27eb37f84694 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 10 Sep 2025 12:00:50 +0200 Subject: [PATCH 04/11] enforce: user must have "activity rate estimation" unselected for custom rate to be used Previously, user could enter a value enter the custom rate box, enable "activity rate estimation" and the custom rate box would disappear but the program would still see the value previously entered and use it even though it was no longer visible in the user interface --- src/seismic_hazard_forecasting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 3cd03c6..6d314f7 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -313,7 +313,7 @@ verbose: {verbose}") # run activity rate modeling lambdas = [None] - if custom_rate != None and forecast_select: + if custom_rate != None and forecast_select and not rate_select: logger.info(f"Using activity rate specified by user: {custom_rate} per {time_unit}") lambdas = np.array([custom_rate], dtype='d') lambdas_perc = np.array([1], dtype='d') -- 2.16.5 From 00bd39a098501b2d15d34eada8e095814e59a6c5 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 10 Sep 2025 16:33:11 +0200 Subject: [PATCH 05/11] impose requirement of minimum size of range of output data to do image processing --- src/seismic_hazard_forecasting.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 6d314f7..12e9abe 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -519,8 +519,8 @@ verbose: {verbose}") dtype=np.float64) # this reduces values to 8 decimal places iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes - # upscale the grid if there are at least 10 grid values - if np.count_nonzero(iml_grid_tmp) >= 10: + # upscale the grid if there are at least 10 grid values with range greater than 0.1 + if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1: up_factor = 4 iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), mode='reflect', anti_aliasing=False) @@ -529,9 +529,10 @@ verbose: {verbose}") iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan - # trim edges so the grid is not so blocky - vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x)) + #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x)) vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x)) + + # trim edges so the grid is not so blocky trim_thresh = vmin iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan -- 2.16.5 From f7eb39c43cee5f9f4aecf88033bb24ef94fef8f5 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 10 Sep 2025 18:39:43 +0200 Subject: [PATCH 06/11] add final image smoothing through binlinear interpolation --- src/seismic_hazard_forecasting.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 12e9abe..257c0f2 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -519,13 +519,17 @@ verbose: {verbose}") dtype=np.float64) # this reduces values to 8 decimal places iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes - # upscale the grid if there are at least 10 grid values with range greater than 0.1 + # upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1 if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1: up_factor = 4 iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), mode='reflect', anti_aliasing=False) + interp_method = 'bilinear' + trim_thresh = vmin + iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan else: iml_grid_hd = iml_grid_tmp + interp_method = 'none' iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan @@ -533,8 +537,8 @@ verbose: {verbose}") vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x)) # trim edges so the grid is not so blocky - trim_thresh = vmin - iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan + #trim_thresh = vmin + #iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan # generate image overlay north, south = lat.max(), lat.min() # Latitude range @@ -547,7 +551,7 @@ verbose: {verbose}") cmap_name = 'YlOrRd' cmap = plt.get_cmap(cmap_name) fig, ax = plt.subplots(figsize=(6, 6)) - ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax) + ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation=interp_method) ax.axis('off') # Save the figure -- 2.16.5 From fe9d886499e12cb97bc52a9abcea60bac04bf5e8 Mon Sep 17 00:00:00 2001 From: ftong Date: Fri, 12 Sep 2025 10:37:03 +0200 Subject: [PATCH 07/11] interpolation is always used on the final grid --- src/seismic_hazard_forecasting.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 257c0f2..ad58232 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -524,22 +524,16 @@ verbose: {verbose}") up_factor = 4 iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), mode='reflect', anti_aliasing=False) - interp_method = 'bilinear' trim_thresh = vmin iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan else: iml_grid_hd = iml_grid_tmp - interp_method = 'none' iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x)) vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x)) - # trim edges so the grid is not so blocky - #trim_thresh = vmin - #iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan - # generate image overlay north, south = lat.max(), lat.min() # Latitude range east, west = lon.max(), lon.min() # Longitude range @@ -551,7 +545,7 @@ verbose: {verbose}") cmap_name = 'YlOrRd' cmap = plt.get_cmap(cmap_name) fig, ax = plt.subplots(figsize=(6, 6)) - ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation=interp_method) + ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear') ax.axis('off') # Save the figure -- 2.16.5 From deb700560481e6a5749df3230516b615eaa8d7ff Mon Sep 17 00:00:00 2001 From: ftong Date: Tue, 23 Sep 2025 11:41:08 +0200 Subject: [PATCH 08/11] Force use of fork in multiprocessing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit From Tomasz Balawajder: "Since we are using a Java service to launch the Python process, its behavior differs from running the script directly on the cluster. By default, Dask uses fork() to create worker processes. However, when running under the JVM, the start method defaults to spawn, which does not share memory between processes. This caused the slowdown and unexpected behavior. I’ve forced Python to use fork() in the configuration, and now the application completes in the same time as when executed with sbatch." --- src/seismic_hazard_forecasting.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index ad58232..95297ba 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -69,6 +69,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max from matplotlib.contour import ContourSet import xml.etree.ElementTree as ET import json + import multiprocessing as mp logger = getDefaultLogger('igfash') @@ -448,9 +449,10 @@ verbose: {verbose}") start = timer() - use_pp = False + use_pp = True if use_pp: # use dask parallel computing + mp.set_start_method("fork", force=True) pbar = ProgressBar() pbar.register() iter = indices -- 2.16.5 From 3136c4985ddd8a2a87db2ce132357ed889a0d467 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 24 Sep 2025 14:05:22 +0200 Subject: [PATCH 09/11] disable cython --- src/seismic_hazard_forecasting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 95297ba..d653038 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -468,7 +468,7 @@ verbose: {verbose}") imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i, - rtol=0.1, use_cython=True) for i in iter] + rtol=0.1, use_cython=False) for i in iter] iml = dask.compute(*imls) iml_grid_raw.append(list(iml)) -- 2.16.5 From d661cad9916274d38386056e3d1c989940b21b4f Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 24 Sep 2025 14:13:21 +0200 Subject: [PATCH 10/11] disable progress bar --- src/seismic_hazard_forecasting.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index d653038..5e3614c 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -52,7 +52,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max from math import ceil, floor, isnan import numpy as np import dask - from dask.diagnostics import ProgressBar # use Dask progress bar + #from dask.diagnostics import ProgressBar # use Dask progress bar import kalepy as kale import utm from skimage.transform import resize @@ -453,8 +453,8 @@ verbose: {verbose}") if use_pp: # use dask parallel computing mp.set_start_method("fork", force=True) - pbar = ProgressBar() - pbar.register() + #pbar = ProgressBar() + #pbar.register() iter = indices iml_grid_raw = [] # raw ground motion grids for imt in products: @@ -468,7 +468,7 @@ verbose: {verbose}") imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i, - rtol=0.1, use_cython=False) for i in iter] + rtol=0.1, use_cython=True) for i in iter] iml = dask.compute(*imls) iml_grid_raw.append(list(iml)) -- 2.16.5 From a5534212ba9c912bc9abe7c4c12fbf15b74a3fbc Mon Sep 17 00:00:00 2001 From: ftong Date: Thu, 25 Sep 2025 12:07:02 +0200 Subject: [PATCH 11/11] cleanup --- src/seismic_hazard_forecasting.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 5e3614c..9e43458 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -52,7 +52,6 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max from math import ceil, floor, isnan import numpy as np import dask - #from dask.diagnostics import ProgressBar # use Dask progress bar import kalepy as kale import utm from skimage.transform import resize @@ -453,8 +452,6 @@ verbose: {verbose}") if use_pp: # use dask parallel computing mp.set_start_method("fork", force=True) - #pbar = ProgressBar() - #pbar.register() iter = indices iml_grid_raw = [] # raw ground motion grids for imt in products: -- 2.16.5