uafgeotools / rtm

A Python package for locating infrasound sources using reverse time migration
https://uaf-rtm.readthedocs.io/
MIT License
38 stars 13 forks source link

Speed up AGC with pandas.rolling #70

Open flixha opened 2 years ago

flixha commented 2 years ago

Hi, I was looking for a an AGC implemented in Python to apply to an obspy.core.Stream and your package came up as the first result. When I tried the AGC for longer streams (1 - 24 hours), I found that it gets quite slow so I looked a bit into how to speed it up. I tried a few things like convolution, but achieved the best speed up with pandas.DataFrame.rolling - for 24 hours of one 100 Hz trace, the speed up was up to x50, e.g., from ~4200 ms to 90 ms for method='walker'.

As I haven't used or thought about RTM much yet I don't know if that is in any way relevant for your implementation; but I thought it's best to share it. For now I kept an option to switch between old and new implementation so that they can be easily compared; and I put some extra code lines in the comments that show other implementations that were not as fast. If you do want to merge the implementation I can of course also clean that up.

Demonstration that implementations give the same result:

import numpy as np
from obspy import read
from rtm.waveform import _agc

st = read()
st.detrend()

# True (though tiny floating point differences, apparently due to mean function)
np.allclose(_agc(st, win_sec=5, method='gismo', method_exec='new')[0].data,
            _agc(st, win_sec=5, method='gismo', method_exec='old')[0].data,
            rtol=1e-12, atol=1e-15)
# True (exact)
(_agc(st, win_sec=5, method='walker', method_exec='new') ==
_agc(st, win_sec=5, method='walker', method_exec='old'))
liamtoney commented 2 years ago

Hi @flixha many thanks for this improvement. I don't recall us noticing issues with the AGC speed in our RTM use cases, but it seems like if this is faster and produces identical results (within machine epsilon) then it'd make sense to include. We'll converse and get back to you.