eqcorrscan / EQcorrscan

Earthquake detection and analysis in Python.
https://eqcorrscan.readthedocs.io/en/latest/
Other
166 stars 86 forks source link

Fix segmentation fault with distance decluster int32 limit #523

Closed flixha closed 1 year ago

flixha commented 1 year ago

What does this PR do?

Fix segmentation fault when using utils.find_peaks.decluster_distance_time with more than 46340 detections (46341^2 --> int32 overflow for distance_index in find_peaks.c/decluster_dist_time and find_peaks.c/decluster_dist_time_ll).

Why was it initiated? Any relevant Issues?

Only distance_index was defined as int, while all other indices / lengths in the c-libraries are at least of type long. So I expect no negative side effects from this fix. Let me know if we need an explicit test for this.

Minimal working example to demonstrate the previous segmentation fault:

import numpy as np
import ctypes
from eqcorrscan.utils.libnames import _load_cdll

# 46341: works   -  46342: fails
limit = 46342

arr = np.zeros([limit, ], dtype='float32')
inds = np.arange(0, limit, dtype='int32')
distance_matrix = np.zeros([limit, limit], dtype='float32')

threshold = 0
trig_int = 20000000 # microseconds?
hypocentral_separation = 200 # km

long_type = ctypes.c_long

arr = np.ascontiguousarray(arr, dtype=np.float32)
inds = np.ascontiguousarray(inds, dtype=long_type)
distance_matrix = np.ascontiguousarray(
    distance_matrix.flatten(order="C"), dtype=np.float32)

out = np.zeros(len(arr), dtype=np.uint32)
length = len(arr)

utilslib = _load_cdll('libutils')

long_type = ctypes.c_long
func = utilslib.decluster_dist_time
func = utilslib.decluster_dist_time_ll

func.argtypes = [
    np.ctypeslib.ndpointer(dtype=np.float32, shape=(length,),
                            flags='C_CONTIGUOUS'),
    np.ctypeslib.ndpointer(dtype=long_type, shape=(length,),
                            flags='C_CONTIGUOUS'),
    np.ctypeslib.ndpointer(dtype=np.float32, shape=(length * length,),
                            flags='C_CONTIGUOUS'),
    long_type, ctypes.c_float, long_type, ctypes.c_float,
    np.ctypeslib.ndpointer(dtype=np.uint32, shape=(length,),
                            flags='C_CONTIGUOUS')]
func.restype = ctypes.c_int

ret = func(
    arr, inds, distance_matrix, long_type(length), np.float32(threshold),
    long_type(trig_int), hypocentral_separation, out)

PR Checklist

flixha commented 1 year ago

Seeing some test failures because NCEDC cannot be reached via FDSN web services...

calum-chamberlain commented 1 year ago

Thanks!