search

jpjones76 / SeisIO.jl

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

detrend! fails with NodalData #54

Closed tclements closed 3 years ago

tclements commented 3 years ago

On the dev branch 

(@v1.5) pkg> st
Status `~/.julia/environments/v1.5/Project.toml`
  [b372bb87] SeisIO v1.1.0 `https://github.com/jpjones76/SeisIO.jl.git#dev`

with the test TDMS data 

julia> N = SeisIO.Nodal.read_nodal("Node1_UTC_20200307_170738.006.tdms")

detrend! fails 

julia> detrend!(N)
ERROR: MethodError: no method matching dtr!(::SubArray{Float32,1,Array{Float32,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}, ::Array{Int64,2}, ::Float64, ::Int64)
Closest candidates are:
  dtr!(::Array{T,1}, ::Array{Int64,2}, ::Float64, ::Int64) where T<:AbstractFloat at /home/timclements/.julia/packages/SeisIO/NGyBv/src/Processing/detrend.jl:74
Stacktrace:
 [1] detrend!(::SeisIO.Nodal.NodalData; chans::Array{Int64,1}, n::Int64) at /home/timclements/.julia/packages/SeisIO/NGyBv/src/Processing/detrend.jl:140
 [2] detrend!(::SeisIO.Nodal.NodalData) at /home/timclements/.julia/packages/SeisIO/NGyBv/src/Processing/detrend.jl:135
 [3] top-level scope at REPL[71]:1
 [4] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1088

I believe this is an issue with the definition of dtr!

Using the NodalData.x[i] syntax returns a SubArray type: 

julia> typeof(N.x[1])
SubArray{Float32,1,Array{Float32,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}

julia> isa(N.x[1],Array)
false

julia> isa(N.x[1],AbstractArray)
true

but dtr! only accepts Arrays. Changing dtr! to accept AbstractArray as function dtr!(x::AbstractArray{T,1}, ti::Array{Int64,2}, fs::Float64, n::Int64) where T <: AbstractFloat should fix this: 

julia> detrend!(N)

works! Should be an easy change.

jpjones76 commented 3 years ago

Testing the fix now. I'm adding automated tests for all processing operations on NodalData (except ungap and merge, which are currently inapplicable). Since the structure has significant differences from SeisData and EventTraceData, I want to be sure that you don't run into problems further down the line.

Incidentally, how do views of an Array{Float32, 2} (equivalently, Matrix{Float32}) behave on GPU with your package?

tclements commented 3 years ago

Views work just the same as normal arrays but I believe they have to be contiguous. Views do not interact well with transposes or reshapes at the moment - https://github.com/JuliaGPU/Adapt.jl/issues/21 but apparently a fix could be on the way at some point https://github.com/JuliaLang/julia/pull/31563

jpjones76 commented 3 years ago

Perfect. I'm pretty sure I have no transpose or reshape calls in SeisIO processing functions at the moment. Eventually, as your group gets further into nodal arrays, you'll be able to write NodalData processing methods into the submodule for two-dimensional FFTs on the underlying :data field. So, I think this will need only a minor tweak to allowed data types to be fully GPU compatible.

jpjones76 commented 3 years ago

Fixed on dev

jpjones76 commented 3 years ago

Fixed on master and in the new release. Will close this issue as soon as Julia Registrator merges my PR(s).