Skip to content

WIP: Weighted correlations #508

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
calum-chamberlain wants to merge 10 commits into
base: develop
Choose a base branch
from
Open

WIP: Weighted correlations #508

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
695ff58
WIP:Add test for weighted correlations
calum-chamberlain Aug 4, 2022
ff66165
Implement weights
calum-chamberlain Aug 5, 2022
3d62d54
Test array correlations
calum-chamberlain Aug 5, 2022
487f99d
Document weighting
calum-chamberlain Aug 5, 2022
45befc7
Merge branch 'develop' into weighted-correlations
calum-chamberlain Aug 8, 2022
bb9c1c0
Merge branch 'develop' into weighted-correlations
calum-chamberlain Aug 8, 2022
3a38bc6
Merge branch 'develop' into weighted-correlations
calum-chamberlain Aug 9, 2022
20109be
Merge branch 'develop' into weighted-correlations
calum-chamberlain Jan 3, 2023
f5b8aa5
Merge branch 'develop' into weighted-correlations
calum-chamberlain Jan 4, 2023
8410d97
Merge branch 'develop' into weighted-correlations
calum-chamberlain Mar 31, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,8 @@
- Made handling of horizontal_chans and vertical_chans consistent so that user
can freely choose relevant channels.
* utils.correlate
- Weighted correlations supported. Template weights are expected in
individual trace.stats.extra.weight attributes.
- Fast Matched Filter now supported natively for version >= 1.4.0
- Only full correlation stacks are returned now (e.g. where fewer than than
the full number of channels are in the stack at the end of the stack, zeros
Expand Down
16 changes: 14 additions & 2 deletions eqcorrscan/doc/submodules/utils.correlate.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,18 @@ supported natively, and can be used as a backend by setting `xcorr_func="fmf"`.

<a href="https://github.com/beridel/fast_matched_filter" target="_blank">Fast Matched Filter</a>

Using weights
~~~~~~~~~~~~~

New for version 0.5.0: support for weighting cross-correlation stacks. All array
correlation functions have an optional `weights` argument to allow you to weight
different templates. This makes most sense when using the stream correlation functions
where different template channels might have different quality. For the stream
correlation functions the weight for each channel of each template is required to
be stored in `trace.stats.extra.weight`. Individual channel cross-correlations
are multiplied by this weight before stacking. Weights can be any floating point
number, although weights between 0 and 1 make the most sense.


Switching which correlation function is used
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down Expand Up @@ -141,10 +153,10 @@ for example:
... 1, False, xcorr_func='numpy')

>>> # do correlation using a custom function
>>> def custom_normxcorr(templates, stream, pads, *args, **kwargs):
>>> def custom_normxcorr(templates, stream, pads, weights *args, **kwargs):
... # Just to keep example short call other xcorr function
... print('calling custom xcorr function')
... return numpy_normxcorr(templates, stream, pads, *args, **kwargs)
... return numpy_normxcorr(templates, stream, pads, weights, *args, **kwargs)

>>> detections = match_filter(
... ['1'], [template], stream, .5, 'absolute', 1, False,
Expand Down
101 changes: 96 additions & 5 deletions eqcorrscan/tests/correlate_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
import pytest
from multiprocessing import cpu_count
from obspy import Trace, Stream, read
from obspy.core.util import AttribDict
from scipy.fftpack import next_fast_len

import eqcorrscan.utils.correlate as corr
Expand Down Expand Up @@ -98,6 +99,16 @@ def generate_multichannel_templates():
return templates


def generate_weighted_multichannel_templates():
random.seed(42)
templates = generate_multichannel_templates()
for template in templates:
for tr in template:
tr.stats.extra = AttribDict()
tr.stats.extra.update({"weight": random.random()})
return templates


def read_gappy_real_template():
return [read(join(pytest.test_data_path, "DUWZ_template.ms"))]

Expand Down Expand Up @@ -164,6 +175,8 @@ def print_class_name(request):
# array fixtures

starting_index = 500
GLOBAL_WEIGHT = 0.5 # Weight applied to array ccs funcs
# global to allow comparison to unweighted


@pytest.fixture(scope='module')
Expand Down Expand Up @@ -194,6 +207,29 @@ def pads(array_template):
return np.zeros(array_template.shape[0], dtype=int)


@pytest.fixture(scope='module')
def array_ccs_weighted(array_template, array_stream, pads):
""" Use each function stored in the normxcorr cache to correlate the
templates and arrays, return a dict with keys as func names and values
as the cc calculated by said function"""
out = {}
weights = np.ones(array_template.shape[0]) * GLOBAL_WEIGHT
for name in list(corr.XCORR_FUNCS_ORIGINAL.keys()):
func = corr.get_array_xcorr(name)
print("Running %s" % name)
cc, _ = time_func(
func, name, array_template, array_stream, pads, weights)
out[name] = cc
if "fftw" in name:
print("Running fixed len fft")
fft_len = next_fast_len(
max(len(array_stream) // 4, len(array_template)))
cc, _ = time_func(func, name, array_template, array_stream, pads,
weights, fft_len=fft_len)
out[name + "_fixed_len"] = cc
return out


@pytest.fixture(scope='module')
def array_ccs(array_template, array_stream, pads):
""" Use each function stored in the normxcorr cache to correlate the
Expand Down Expand Up @@ -251,6 +287,12 @@ def multichannel_templates():
return generate_multichannel_templates()


@pytest.fixture(scope='module')
def weighted_multichannel_templates():
""" create weighted multichannel templates """
return generate_weighted_multichannel_templates()


@pytest.fixture(scope='module')
def multichannel_stream():
""" create a multichannel stream for tests """
Expand Down Expand Up @@ -290,8 +332,7 @@ def real_multichannel_stream():
if fname != 'default'}


@pytest.fixture(scope='module')
def stream_cc_output_dict(multichannel_templates, multichannel_stream):
def _stream_cc_output_dict(multichannel_templates, multichannel_stream):
""" return a dict of outputs from all stream_xcorr functions """
# corr._get_array_dicts(multichannel_templates, multichannel_stream)
out = {}
Expand Down Expand Up @@ -320,12 +361,35 @@ def stream_cc_output_dict(multichannel_templates, multichannel_stream):
return out


@pytest.fixture(scope='module')
def stream_cc_output_dict(multichannel_templates, multichannel_stream):
""" return a dict of outputs from all stream_xcorr functions """
return _stream_cc_output_dict(
multichannel_templates=multichannel_templates,
multichannel_stream=multichannel_stream
)


@pytest.fixture(scope='module')
def stream_cc_weighted_output_dict(weighted_multichannel_templates, multichannel_stream):
return _stream_cc_output_dict(
multichannel_templates=weighted_multichannel_templates,
multichannel_stream=multichannel_stream
)


@pytest.fixture(scope='module')
def stream_cc_dict(stream_cc_output_dict):
""" return just the cc arrays from the stream_cc functions """
return {name: result[0] for name, result in stream_cc_output_dict.items()}


@pytest.fixture(scope='module')
def stream_cc_weighted_dict(stream_cc_weighted_output_dict):
""" return just the cc arrays from the stream_cc functions """
return {name: result[0] for name, result in stream_cc_weighted_output_dict.items()}


@pytest.fixture(scope='module')
def real_stream_cc_output_dict(real_templates, real_multichannel_stream):
""" return a dict of outputs from all stream_xcorr functions """
Expand Down Expand Up @@ -589,6 +653,21 @@ def test_single_channel_similar(self, array_ccs):
np.save("cc2.npy", cc2)
assert np.allclose(cc1, cc2, atol=self.atol)

def test_known_weight_application(self, array_ccs, array_ccs_weighted):
cc_names = list(array_ccs.keys())
print(cc_names)
for key1, key2 in itertools.combinations(cc_names, 2):
cc1 = array_ccs_weighted[key1]
cc2 = array_ccs_weighted[key2]
print(f"Comparing {key1} and {key2}")
assert np.allclose(cc1, cc2, atol=self.atol)
for cc_name in cc_names:
print(f"Checking for {cc_name}")
assert np.allclose(
array_ccs_weighted[cc_name],
array_ccs[cc_name] * GLOBAL_WEIGHT,
atol=1e-5)

def test_autocorrelation(self, array_ccs):
""" ensure an autocorrelation occurred in each of ccs where it is
expected, defined by starting_index variable """
Expand Down Expand Up @@ -623,6 +702,17 @@ def test_multi_channel_xcorr(self, stream_cc_dict):
for cc_name, cc in zip(cc_names[2:], cc_list[2:]):
assert np.allclose(cc_1, cc, atol=self.atol)

def test_weighted_multi_channel_xcorr(self, stream_cc_weighted_dict):
""" test various correlation methods weighted with multiple channels """
# get correlation results into a list
cc_names = list(stream_cc_weighted_dict.keys())
cc_list = [stream_cc_weighted_dict[cc_name] for cc_name in cc_names]
cc_1 = cc_list[0]
# loop over correlations and compare each with the first in the list
# this will ensure all cc are "close enough"
for cc_name, cc in zip(cc_names[2:], cc_list[2:]):
assert np.allclose(cc_1, cc, atol=self.atol)

def test_real_multi_channel_xcorr(self, real_stream_cc_dict):
""" test various correlation methods with multiple channels """
# get correlation results into a list
Expand Down Expand Up @@ -837,10 +927,11 @@ def test_callable_registered(self, multichannel_templates,
""" ensure a callable can be registered """
small_count = {}

def some_callable(template_array, stream_array, pad_array):
def some_callable(template_array, stream_array, pad_array,
weight_array):
small_count['name'] = 1
return corr.numpy_normxcorr(template_array, stream_array,
pad_array)
pad_array, weight_array)

func = corr.get_stream_xcorr(some_callable)
func(multichannel_templates, multichannel_stream)
Expand Down Expand Up @@ -883,7 +974,7 @@ class TestRegisterNormXcorrs:
# helper functions
def name_func_is_registered(self, func_name):
""" return True if func is registered as a normxcorr func """
# Note: don not remove this fixture or bad things will happen
# Note: do not remove this fixture or bad things will happen
name = func_name.__name__ if callable(func_name) else func_name
return name in corr.XCOR_FUNCS

Expand Down
Loading