Rejection based on selection mode causes a missing window with respect to Flexwin

[rdno]

Hi,

While selection_mode is set to all_waves, pyflex is not selecting the first window that Flexwin selects. You can see the outputs below.

I am using the default configuration for Flexwin. Reason for discrepancy turned out to be the rejection rules. Flexwin does not reject the windows that fall in the noise area. Yet, pyflex does.

I added a custom option to test this issue. With this option reject_on_selection_mode method does not do any rejection. Last plot shows the output that uses this option.

Rejecting the windows in the noise region seems like a sane default. Although, I think it can be helpful to give the user full control of the rejection mechanisms.

What do you think?

Flexwin:

Pyflex:

Pyflex with no rejection:

[wjlei1990]

Hi Ridvan,

Thanks for the very careful check.

I think there are a few issues as you mentioned:

1. pyflex has the mechanism to reject the windows based on [noise_start_index, noise_end_index]. But it seems the value of the noise_end_index is not calculated in a proper way, which makes the pyflex reject the first phase. I revisited the code that used to generate the noise_start_index. The original code is:

                 self.config.noise_end_index = \
                     int((self.ttimes[0]["time"] - 2 * self.config.min_period) *
                         self.observed.stats.sampling_rate)
   

   It seems in the code, the offset is not counted in, which might cause problems in the future.

   But anyway, I think in the example you mentioned, it is clear we need to change this strategy a little bit to make the noise region correct. What do think about this modification:

         # reject windows which has overlap with the noise region
         self.windows = [win for win in self.windows
                         if (win.right > self.config.noise_end_index)]
   

   Instead of using win.left, we use win.right larger than the noise_end_index?

2. the selection_mode is not quite clear.
   So what is all_waves mean? Should it just include every thing, even in the noise region? or it should still work in the current way? Thus, the user might get confused by the definition of "all_waves" option.
   Currently, the selection mode has ["all_waves", "body_waves", "surface_waves", "body_and_surface_waves", "mantle_waves"].

   I think it is a good suggestion to add another mode called "custom" so reject_on_selection_mode will be skipped.

[krischer]

self.config.noise_end_index = \
    int((self.ttimes[0]["time"] - 2 * self.config.min_period) *
        self.observed.stats.sampling_rate)

This, if I remember correctly, is just used to determine the end of the noise based on the first arrival. What do you mean by offset in this example?

Twice the minimum period is a bit of an ad-hoc value. It assumes that a ancausal zerophase filter has been used and this causes the first arrival to have some energy a bit earlier. Thus the noise should end a bit earlier so pyflex can still pick it up.

Instead of using win.left, we use win.right larger than the noise_end_index?

Hmm...why not I guess - the later cross correlation and other checks should ensure that the picked windows are at least reasonable.

I added a custom option to test this issue. With this option reject_on_selection_mode method does not do any rejection. Last plot shows the output that uses this option.

I have to look at the code again but would this not potentially also pick some very late windows? This would be a bit dangerous as synthetic as well as observed data are both band limited and thus there is a high change for some good fit late in the wave train which has no proper physical reason and the kernels would be meaningless.

[rdno]

the selection_mode is not quite clear.
So what is all_waves mean? Should it just include every thing, even in the noise region? or it should still work in the current way? Thus, the user might get confused by the definition of "all_waves" option.
Currently, the selection mode has ["all_waves", "body_waves", "surface_waves", "body_and_surface_waves", "mantle_waves"].

I think it is quite reasonable for all_waves to exclude the noise region.

I have to look at the code again but would this not potentially also pick some very late windows? This would be a bit dangerous as synthetic as well as observed data are both band limited and thus there is a high change for some good fit late in the wave train which has no proper physical reason and the kernels would be meaningless.

Yes. It can be harmful. But I am suggesting this option for a specific use case. Let's say selection_mode options do not cover what you want. Instead of editing the source code, you can tell pyflex to do no rejection based on selection_mode. After that, you can use waterlevel/signal2noise ratio/dlnA arrays (like you do in Flexwin) to select the parts of the seismogram that you are interested in. Basically, user should adjust the parameters to avoid the false positives while using this option.

[krischer]

Yes. It can be harmful. But I am suggesting this option for a specific use case. Let's say selection_mode options do not cover what you want. Instead of editing the source code, you can tell pyflex to do no rejection based on selection_mode. After that, you can use waterlevel/signal2noise ratio/dlnA arrays (like you do in Flexwin) to select the parts of the seismogram that you are interested in. Basically, user should adjust the parameters to avoid the false positives while using this option.

👍

Another option might be to set the left borders of windows that are overlapping with the noise to the left-most accepted sample (first arrival - 2 * min period). I think this would also solve the particular issue shown here.

[wjlei1990]

Hi Lion,

For the offset issue, I mean here:

                self.config.noise_end_index = \
                    int((self.ttimes[0]["time"] - 2 * self.config.min_period) *
                        self.observed.stats.sampling_rate)

The noise_end_index is calculated based on first arrival, so the timing is related to the event origin time. However, when we reject the windows,

        self.windows = [win for win in self.windows
                        if (win.left > self.config.noise_end_index)]

The win.left timing is related the first point of seismogram. So the timings for these two cases are different.

Correct me if I am wrong.

[wjlei1990]

I think we all agree with the "custom" mode should be added. Ok, so one more thing, if using custom mode, should pyflex still reject windows in noise region? I think the answer is yes, right? So we should split the function reject_on_selection_mode to reject_on_noise_region and reject_on_selection_mode. If using custom mode, then reject_on_selection mode should be skipped.

Second issue, for modifying the noise_end_index. There are few scenarios you guys and I already mentioned about:

1. using the win.right rather than the win.left
2. reset the win.left to noise_end_index if the window overlaps with the noise region
3. using larger values of the ad-hoc variable, for example,
   self.config.noise_end_index = \ int((self.ttimes[0]["time"] - 4 * self.config.min_period) * self.observed.stats.sampling_rate)
   Also, we need to make sure noise_end_index is larger than 0.

[rdno]

using the win.right rather than the win.left

I can confirm using win.right instead of win.left solves this particular issue.

reset the win.left to noise_end_index if the window overlaps with the noise region

I also tried this using the code below:

for window in self.windows:
    if window.left < self.config.noise_end_index:
        window.left = self.config.noise_end_index

# reject windows which has overlap with the noise region
self.windows = [win for win in self.windows
                if (win.left >= self.config.noise_end_index)]

using larger values of the ad-hoc variable, for example, self.config.noise_end_index = \ int((self.ttimes[0]["time"] - 4 * self.config.min_period) * self.observed.stats.sampling_rate) Also, we need to make sure noise_end_index is larger than 0.

I am not sure if this is a right solution or not. How would you decide how much further you need to go? Why not use 6 minimum period offset?

I am including the plots of the outputs. Maybe, it can be helpful to decide which way to go. cc/dt/dA values stay the same for both case (First arrival is at 697 s, min period is 45 s, noise_end is 607 s).

If you reject the window based on win.right instead of win.left, first window starts at 590:

If you set the win_left to noise_end first window starts at 607 (noise_end):

[wjlei1990]

Hi Ridvan,

Thanks for you experiments.

I think both methods has its cons and pros.

1. usinng win.right. If so, we lose control of the win.left, which may cause win.left to be on the very left(window too large)...

2. using noise_end_index as win.left when overlapping. First, the noise_end_index is what we set, so it is subjective.

   self.config.noise_end_index = \
   int((self.times[0]["time"] - 2 * self.config.min_period) *
       self.observed.stats.sampling_rate)
   

   Second, the window might still get rejected at the later stage because of window length being too short.

But personally, I favor the second option. I think it is safer, though windows might still get rejected.

[rdno]

Hi Wenjie,

How about using win.center instead of win.left?

Using center value ensures that majority of the window (at least 50%) in the signal region.

[wjlei1990]

I think the win.center not necessary keep the 50% of windows. The win.center is the a maximum value of STA/LTA in that window...So theoretically, it could be very close to win.left or win.right.

But I think that might be a valid option. It at least ensures the maximum value is not in noise region...

[rdno]

Oh, I see. I thought it was the central point. "Center point" can be calculated locally to achieve this.

I suggested this because it has the advantage of not manipulating window borders manually and noise region inclusion is potentially a little bit better than win.right method.

Doesn't changing 'win.left' value has some potential to result in not up-to-date cc/dt/dlna values? Maybe this has to be checked as well.

[wjlei1990]

Yes. That is true. If we change the win.left, then the measurements needs to be made again. The other way is we call move the win.left to noise_end_idx at the very beginning when we generate those windows. But I think that is probably not so good?

Talking about the win.center method, if one is concerned if the win.left is too far away, then we can also check the win.left, say:

win.left >= noise_end_idx - 2 * min_period / dt

To have a little more constrain on the win.left.

Anyway, I will implement it tomorrow and please let me if you have more suggestions.

:)

[wjlei1990]

Hi Ridvan, I have some updates on the dev branch of pyflex. Please have a look.

There are a few things you may need to pay attention:

1. in config, I added one more parameters, which will be used in the noise region calculation and signal region calculation. So you may need try a few numbers to make it the way you want.

However, I still don't consider the condition for major-arc surface wave. In this case, you need explicitly set the singal_end_index to -1 or npts.

2. I added another selection_mode called custom, so no windows will be rejected in reject_on_selection_mode.

Please check my commit log for more details.

[rdno]

Hi,

However, I still don't consider the condition for major-arc surface wave. In this case, you need explicitly set the singal_end_index to -1 or npts.

max_time_after_last_arrival can be also used for this purpose, right?

I added another selection_mode called custom, so no windows will be rejected in reject_on_selection_mode.

👍

I encountered one minor thing by accident, if the max_time_before_first_arrival is set to high, noise_end becomes zero then if condition at window_selector.py#L487-L489 does not pass. This results in None value for signal_start_index. Hence, window selection does not work. Maybe explicit error stating noise region length is zero can be helpful for this edge case.

[wjlei1990]

Nice catch. I really shouldn't use if 0 since it will catch both None and 0. I fixed it.

Another point, I double check the code. For current version of the code, if you set noise_end_index to 0, explicitly or implicitly(by max_time_before_first_arrival), it won't pass because I have a check in _convert_negative_index that make sure noise_end_index is larger than noise_start_index, not even equal. Because if noise_end_index - noise_start_index == 0, then there are two functions, which evaluate the noise level, will fail.

So I can either put this restriction at the very beginning or put this check in the noise level evaluation functions.

I am not sure if it is too strict. What do you think?