Differences in PZResp and MulitStageResponse

[dylanmikesell]

I don't think there is an error/bug anymore. Once I started looking at things correctly and paying a little more attention. Here is a minimum working example.

using SeisIO

S1 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=600, msr=true)
print(S1.resp[1].stage[1],"\n")

S2 = get_data("FDSN", "TA.I02A..BHZ", src="IRIS", s=string("2007-04-16"), t=600)
print(S2.resp)

Here is the output:

PZResp64 with fields: a0: 3.48462e17 f0: 0.2 p: Complex{Float64}[-13300.0 + 0.0im, -10530.0 + 10050.0im, -10530.0 - 10050.0im, -520.3 + 0.0im, -374.8 + 0.0im, -97.34 + 400.7im, -97.34 - 400.7im, -15.64 + 0.0im, -0.037 + 0.037im, -0.037 - 0.037im, -255.1 + 0.0im] z: Complex{Float64}[0.0 + 0.0im, 0.0 + 0.0im, -463.1 + 430.5im, -463.1 - 430.5im, -176.6 + 0.0im, -15.15 + 0.0im]

InstrumentResponse[PZResp with fields: a0: 3.48462e17 f0: 0.2 p: Complex{Float32}[-13300.0f0 + 0.0f0im, -10530.0f0 + 10050.0f0im, -10530.0f0 - 10050.0f0im, -520.3f0 + 0.0f0im, -374.8f0 + 0.0f0im, -97.34f0 + 400.7f0im, -97.34f0 - 400.7f0im, -15.64f0 + 0.0f0im, -0.037f0 + 0.037f0im, -0.037f0 - 0.037f0im, -255.1f0 + 0.0f0im] z: Complex{Float32}[0.0f0 + 0.0f0im, 0.0f0 + 0.0f0im, -463.1f0 + 430.5f0im, -463.1f0 - 430.5f0im, -176.6f0 + 0.0f0im, -15.15f0 + 0.0f0im] ]

The difference is that the MSR is a 64 bit float, while the normal PZResp is 32 bit. I see now after digging deeper into the first stage of the MSR that it even says it is PZResp64 and not a PZResp.

Sorry to bother you will this. Not sure if it causes problems down the line in SeisNoise though.

[jpjones76]

Flagging as a consistency issue. Do you need 64-bit precision, or 32-bit ok?

If you're comfortable with 32-bit, an easy way to achieve consistency is for me to change the XML parser so that MultiStageResp uses 32-bit PZResp objects, rather than PZResp64.

64-bit precision is more difficult. I've been developing SeisIO to work with Float32 by default because it speeds up the arithmetic. @mdenolle and @tclements assure me that this becomes much greater than a factor-of-2 speed increase on current GPUs, but I've done no GPU benchmarking in Julia. Hopefully they can explain what they've found.

I'm almost certain that @tclements is developing SeisNoise.jl to also use 32-bit precision by default. However, I want to be certain of that before I make changes that could affect what he's doing downstream.

[dylanmikesell]

32-bit works just fine for my needs right now. Thanks,

Dylan

On Wed, Mar 25, 2020, 4:08 PM Joshua Jones notifications@github.com wrote:

Flagging as a consistency issue. Do you need 64-bit precision, or 32-bit ok?

If you're comfortable with 32-bit, an easy way to achieve consistency is for me to change the XML parser so that MultiStageResp uses 32-bit PZResp objects, rather than PZResp64.

64-bit precision is more difficult. I've been developing SeisIO to work with Float32 by default because it speeds up the arithmetic. @mdenolle https://github.com/mdenolle and @tclements https://github.com/tclements assure me that this becomes much greater than a factor-of-2 speed increase on current GPUs, but I've done no GPU benchmarking in Julia. Hopefully they can explain what they've found.

I'm almost certain that @tclements https://github.com/tclements is developing SeisNoise.jl to also use 32-bit precision by default. However, I want to be certain of that before I make changes that could affect what he's doing downstream.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jpjones76/SeisIO.jl/issues/41#issuecomment-604113019, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVIE4VVTZNW7KCIEZ3MTDTRJJ6F7ANCNFSM4LTXVNLQ .

[tclements]

Almost all processing in SeisNoise uses Float32 by default. I try to make all functions type stable to prevent slowdowns from casting from Float32 to Float64 or vice-versa. I'm in the process of moving all parameters, such as the sampling rate, from Float64 to AbstractFloat. The only thing I do specifically in Float64 is the creation of Butterworth filters on the GPU due to the need for extra precision when dividing by small numbers https://github.com/tclements/SeisNoise.jl/blob/GPU/src/filter.jl#L75-L135

[dylanmikesell]

Can you educate me on the reason for going to AbstractFloat rather than Float32? All of the floats fall under AbstractFloat; will Float32 be automatically chosen by Julia? I still consider myself new to Julia, so not sure why you would not be explicitly with a Float32 type. Thanks.

On Thu, Mar 26, 2020 at 10:27 AM Tim Clements notifications@github.com wrote:

Almost all processing in SeisNoise uses Float32 by default. I try to make all functions type stable to prevent slowdowns from casting from Float32 to Float64 or vice-versa. I'm in the process of moving all parameters, such as the sampling rate, from Float64 to AbstractFloat. The only thing I do specifically in Float64 is the creation of Butterworth filters on the GPU due to the need for extra precision when dividing by small numbers https://github.com/tclements/SeisNoise.jl/blob/GPU/src/filter.jl#L75-L135

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jpjones76/SeisIO.jl/issues/41#issuecomment-604530018, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVIE4QL7UV3XKTL3DX25K3RJN667ANCNFSM4LTXVNLQ .

[tclements]

Defining a function with AbstractFloat inputs allows for inputs with any Float but will throw an error if the input is an Int or Rational. In the long run this will allow users to be more flexible - GPU users can use Float32 while people who need more precision to do inversions can use Float64 using the same code through multiple dispatch. Here's a simple example where square2 is more flexible than square while both disallow Int inputs:

julia> square(x::Float32) = x * x # function for only Float32
square (generic function with 1 method)

julia> square2(x::AbstractFloat) = x * x # function for any type of Float
square2 (generic function with 1 method)

julia> square(Float32(2)) # works 
4.0f0

julia> square(Float64(2)) # fails
ERROR: MethodError: no method matching square(::Float64)
Closest candidates are:
  square(::Float32) at REPL[1]:1
Stacktrace:
 [1] top-level scope at REPL[4]:1

julia> square2(Float32(2)) # returns Float32
4.0f0

julia> square2(Float64(2)) # returns Float64
4.0

julia> square2(2) # fails for Int 
ERROR: MethodError: no method matching square2(::Int64)
Closest candidates are:
  square2(::AbstractFloat) at REPL[2]:1
Stacktrace:
 [1] top-level scope at REPL[7]:1

[dylanmikesell]

Got it. Thanks!

On Thu, Mar 26, 2020 at 11:37 AM Tim Clements notifications@github.com wrote:

Defining a function with AbstractFloat inputs allows for inputs with any Float but will throw an error if the input is an Int or Rational. In the long run this will allow users to be more flexible - GPU users can use Float32 while people who need more precision to do inversions can use Float64 using the same code through multiple dispatch https://docs.julialang.org/en/v1/manual/methods/. Here's a simple example where square2 is more flexible than square while both disallow Int inputs:

julia> square(x::Float32) = x * x # function for only Float32 square (generic function with 1 method)

julia> square2(x::AbstractFloat) = x * x # function for any type of Float square2 (generic function with 1 method)

julia> square(Float32(2)) # works 4.0f0

julia> square(Float64(2)) # fails ERROR: MethodError: no method matching square(::Float64) Closest candidates are: square(::Float32) at REPL[1]:1 Stacktrace:

julia> square2(Float32(2)) # returns Float324.0f0

julia> square2(Float64(2)) # returns Float644.0

julia> square2(2) # fails for Int ERROR: MethodError: no method matching square2(::Int64) Closest candidates are: square2(::AbstractFloat) at REPL[2]:1 Stacktrace:

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jpjones76/SeisIO.jl/issues/41#issuecomment-604572675, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVIE4VO4WAQLPAYZJS26TTRJOHGVANCNFSM4LTXVNLQ .

[jpjones76]

Almost all processing in SeisNoise uses Float32 by default. I try to make all functions type stable to prevent slowdowns from casting from Float32 to Float64 or vice-versa. I'm in the process of moving all parameters, such as the sampling rate, from Float64 to AbstractFloat. The only thing I do specifically in Float64 is the creation of Butterworth filters on the GPU due to the need for extra precision when dividing by small numbers https://github.com/tclements/SeisNoise.jl/blob/GPU/src/filter.jl#L75-L135

This discussion worries me a little. I allow Float32 precision in SeisIO filtering because conversion to Float64 (e.g., what DSP.jl does) leads to 400% memory overhead on Float32 data.

I use a water level correction in the spectral domain as a workaround. That used to be a standard practice for seismic data processing. SAC and ObsPy both use water level corrections in some algorithms; ObsPy doesn't mention one in its bandpass filter documentation, but it's in their instrument response translation.

I haven't yet encountered a case where my approach led to a divide by zero, but it sounds like you've done much more extensive testing.

How often are you encountering divide by zero errors at Float32 precision? Can you give me an example set of filtering parameters that throws an error?

[tclements]

I construct filter coefficients in Float64 on the GPU then convert back to Float32 before filtering like this

b = coefb(tf)
a = coefa(tf)
newb = CuArrays.zeros(T,N)
newa = CuArrays.zeros(T,N)
copyto!(newb,b)
copyto!(newa,a)
bafft = rfft(newb) ./ rfft(newa)

# convert back to Float32 for speed
bafft = map(x->convert(ComplexF32,x), bafft)

Doing rfft(newb) ./ rfft(newa) in Float32 led to floating point errors due to some non-zero coefficients being less than eps(Float32) .. I've been too lazy to write CUDA kernels using blocks and threads for anything yet so I'm doing filtering on the GPU as multiplication in frequency domain

# apply filter
A .= irfft(rfft(A,1) .* bafft, N, 1)
reverse!(A,dims=1)
A .= irfft(rfft(A,1) .* bafft, N, 1)
reverse!(A,dims=1)

[jpjones76]

Has this been resolved, or should I leave it open a bit longer?

[tclements]

Yes this is good to close on my end.