jpjones76 / SeisIO.jl

Julia language support for geophysical time series data
http://seisio.readthedocs.org
Other
47 stars 20 forks source link

Possible performance improvement for t_expand and t_collapse #31

Closed tclements closed 4 years ago

tclements commented 4 years ago

Running sync is a slight performance bottleneck. For instance read_data with 20 Hz data takes less than 1 ms and 6 MB

@btime S = read_data("sac",files[1])
  883.761 μs (107 allocations: 6.60 MiB)

while sync on the same file (after merging, ungapping) is more time and memory intensive

@btime sync(S,s=starttime,t=endtime)
  20.251 ms (261 allocations: 76.47 MiB)

Much of this time and memory is spent in t_expand and t_collapse

@btime t = SeisIO.t_expand(S.t[1], S.fs[1])  
4.852 ms (18 allocations: 26.37 MiB)

@btime SeisIO.t_collapse(t, S.fs[1])
  6.903 ms (48 allocations: 39.76 MiB)

For channels without gaps, it could be more efficient to use a range rather than array in t_expand

@btime t_range = range(S[i].t[1,2], stop= S[i].t[1,2] + dt * (S[i].t[2,1]-1),length=S[i].t[2,1])
  5.554 μs (66 allocations: 5.11 KiB)

Here's a check that the two methods return the same representation

dt = round(Int64, 1.0/(S[1].fs*μs))
t_range = range(S[i].t[1,2], stop= S[i].t[1,2] + dt * (S[i].t[2,1]-1),length=S[i].t[2,1])
t = SeisIO.t_expand(S.t[1], S.fs[1])
all(t .== t_range)
true

This change would give different types of output for data with and without gaps, though it looks like t_collapse is only used in sync.jl while t_expand is in sync.jl and SAC.jl (but the Float32.(μs*(t_expand(t, fs) .- ts)) line will create an Array if the input is a range or an array).

This could be a simple change with multiple-dispatch handling the discrepancy between the array and range. Thoughts?

jpjones76 commented 4 years ago

It's definitely worth it to reoptimize these functions if that won't break anything else.

t_expand and t_collapse are intended to resolve more difficult timing issues (secondarily, t_expand is a good fallback for visualization). I've known from the beginning that this could be done with ranges but there were major memory problems with time gaps. I'll have to check whether or not the Julia developers improved that.

I can definitely recode sync for better performance without gaps, but it's not a good idea to return different Types from the same function. That can lead to instabilities.

A workaround might be for me to only call t_expand and its inverse on gapped data. I'll try this if it turns out that I can't feasibly change t_expand itself.

jpjones76 commented 4 years ago

Your idea works.

I finally had time to work on this a few days ago. I'm revising sync! to avoid t_expand and t_collapse -- which I intend to keep in the code for "last resort" procedures, but which shouldn't be needed in such a routine processing operation.

The revised sync, based on your suggestion, produces identical results in my tests and doesn't change any APIs ... and only scales to ~4K memory, for any trace length. This is tremendous improvement. Many thanks for the great suggestion!

I'm going to include the sync change in a minor version release (v1.1.0) in a couple of weeks, but the new version will be live on master "soon" (next week, at the latest; hoping for today, but tomorrow is more realistic). The version number is already incremented to 1.1.0 on master, but I want this change in the "release" version.

My hope is that you can report on the relative performance of the new sync with real data, assuming you aren't too busy. Lately you've had the best data sets for identifying bugs.

Thanks for your patience here.

tclements commented 4 years ago

Sounds great! I can test this on a few TB of data once you commit to master.

jpjones76 commented 4 years ago

Change is now live on master. Please let me know how your testing goes.

jpjones76 commented 4 years ago

Did the new version of sync work for you?

tclements commented 4 years ago

Yes can safely close