uafgeotools / mtuq

moment tensor uncertainty quantification
BSD 2-Clause "Simplified" License
68 stars 23 forks source link

Polarity misfit is not computed in Waveforms+Polarities.py #285

Open Rinku-274 opened 1 week ago

Rinku-274 commented 1 week ago

In the example file Waveforms+Polarities.py, the polarity misfit is missing in the final results results = results_bw + results_sw (line number 206). Also, it seems that the weight for polarity and waveform misfit has not been implemented (equation 3, Silwal et al. 2018). This needs to be addressed for consistency.

rmodrak commented 1 week ago

Very helpful, thanks for pointing this out.

Yes, so far the Waveforms+Polarities.py example only calculates predicted polarities.

I can't access the full text for Silwal et al. 2018, but I'd agree that the relative weighting between waveforms and polarities is important. I'd imagine some sort of user choice / tuning capability may be helpful?

rmodrak commented 1 week ago

before going forward on this, I wonder if anyone might be able to respond to issue #204 ?

Rinku-274 commented 1 week ago

Very helpful, thanks for pointing this out.

Yes, so far the Waveforms+Polarities.py example only calculates predicted polarities.

I can't access the full text for Silwal et al. 2018, but I'd agree that the relative weighting between waveforms and polarities is important. I'd imagine some sort of user choice / tuning capability may be helpful?

Thank you, @rmodrak , for the prompt response. I agree that adding a user choice could be very helpful. Regarding the line: polarities = np.array([-1, -1, -1, 1, 1, 0, 1, 1, -1, 1, 1, 1, 0, 1, 1, 1, -1, 1, 1, 0]) (line 91 in Waveforms+Polarities.py), it can be confusing for users as it’s unclear which station corresponds to which polarity value. A potential improvement would be to allow the code to read polarity information from a separate file. For example, this could be an additional column in the weight file that specifies the polarity values for the corresponding stations.

thurinj commented 4 days ago

Hi @Rinku-274 ,

The current implementation calls the following method internally to get the first motion polarity input from "data" :

def get_observed(self, data):
    """ Extracts observed polarities from data
    """

    if isinstance(data, mtuq.Dataset):
        observed = np.array([_get_polarity(stream) for stream in data])

    elif isinstance(data, list):
        observed = np.array(data)

    elif isinstance(data, np.ndarray):
        observed = data

    else:
        raise TypeError

    return observed

with _get_polarity not being implemented, I believe (at least it is not in the PolarityMisfit class). So, the current method requires passing an array or a list -- meaning there is room for improvement. Do you have an idea of how the external file would be structured? Is there a typical file format for first motion polarities that would be readily compatible with how data are passed in MTUQ already?

As for the total misfit results, thanks for pointing out that we should absolutely include the polarity misfit result here, and that will be an easy fix. However, there is no dedicated implementation of the polarity vs. waveform misfit, as the choice of how to combine misfit is entirely open to the user.

This is an example of how it can be done:

from mtuq.misfit.waveform import calculate_norm_data
norm_bw = calculate_norm_data(data_bw, 
    misfit_bw.norm, ['Z', 'R'])
norm_sw = calculate_norm_data(data_sw, 
    misfit_sw.norm, ['Z', 'R', 'T'])

results = (results_bw/norm_bw + results_sw/norm_sw + results_polarity/len(polarities))/3

This would put all three misfit measures at even weight. Then, you can start experimenting and refine how you want to combine the results (and having this being open-ended is a good way to make sure this is as flexible as possible).

@rmodrak I am going to have a look at issue #204.

Thanks a lot for the feedback 👍.