Gpu (#42)

* Change cc_len to number of samples for GPU testing

* Data length check

* Add ungap to RawData creation

* fix error for reference to cc_len in correlate

* Fix clipping routine

* Add zerophase option to clean_up!

* Fix setting freqmin/freqmax in lowpass/highpass

* Detrend improvements

* Add channels to saving FFTData and CorrData with JLD2

* Deprecate compute_fft, compute_cc. Move phase to compute_fft.jl. Add phasecorrelate

* Convert cc_len and cc_step to Float64 but accept either Float or Int input

* update READme with new processing

* Fix RawData slicing

* Add NoiseData + ArrayFuncs tests

* Remove old submodules

* Remove availability, downsample and transfer. Remove InputParams from SeisNoise.jl

* Fix downsample removal

* Add distance testset

* Fix get_loc error where lon < -180

* Remove double distributed from SeiseNoise.jl

Co-authored-by: timclements <timclements@pop-os.localdomain>

Project.toml

@@ -34,5 +34,5 @@ Glob = "1.2"
 JLD2 = "0.1"
 Plots = "0.29"
 SeisIO = "1.0"
-StatsBase = "0.32,0.32, 0.33"
+StatsBase = "0.32,0.33"
 julia = "1"

README.md

@@ -26,7 +26,7 @@ cross-correlating. For example
 ```Julia
 using SeisNoise, SeisIO
 fs = 40. # sampling frequency in Hz
-freqmin,freqmax = 0.1,0.2 # min and max frequencies 
+freqmin,freqmax = 0.1,0.2 # min and max frequencies
 cc_step, cc_len = 450, 1800 # corrleation step and length in S
 maxlag = 60. # maximum lag time in correlation
 s = "2019-02-03"
@@ -39,9 +39,9 @@ R = RawData.([S1,S2],cc_len,cc_step)
 detrend!.(R)
 taper!.(R)
 bandpass!.(R,freqmin,freqmax,zerophase=true)
-FFT = compute_fft.(R)
+FFT = rfft.(R)
 whiten!.(FFT,freqmin,freqmax)
-C = compute_cc(FFT[1],FFT[2],maxlag)
+C = correlate(FFT[1],FFT[2],maxlag)
 clean_up!(C,freqmin,freqmax)
 abs_max!(C)
 corrplot(C)
@@ -54,51 +54,49 @@ will produce this figure:
 
 SeisNoise can process data and compute cross-correlations on the GPU with CUDA. The [JuliaGPU](https://github.com/JuliaGPU) suite provides a high-level interface for CUDA programming through the CuArrays.jl and CUDAnative.jl packages. CuArrays.jl provides an array type for storing data on the GPU. Data in SeisNoise structures (`R.x`, `F.fft`, and `C.corr` fields, for `RawData`, `FFTData`, and `CorrData`, respectively) can move between an `Array` on the CPU to a `CuArray` on the GPU using the `gpu` and `cpu` functions, as shown below.   
 
-> :warning: Only **Nvidia** GPUs are suported at the moment. Hold in there for AMD/OpenCL support... 
+> :warning: Only **Nvidia** GPUs are suported at the moment. Hold in there for AMD/OpenCL support...
 
 ```julia
-# create raw data and send to GPU 
-R = RawData(S1, cc_len, cc_step) |> gpu 
+# create raw data and send to GPU
+R = RawData(S1, cc_len, cc_step) |> gpu
 R.x
 72000×188 CuArrays.CuArray{Float32,2,Nothing}
 
 # send data back to the CPU
-R = R |> cpu 
+R = R |> cpu
 R.x
 72000×188 Array{Float32,2}
 ```
 
-All basic processing remains the same on the GPU. Here is a complete cross-correlation routine on the GPU: 
+All basic processing remains the same on the GPU. Here is a complete cross-correlation routine on the GPU:
 
 ```julia
-# send data to GPU 
-R1 = RawData(S1, cc_len, cc_step) |> gpu 
+# send data to GPU
+R1 = RawData(S1, cc_len, cc_step) |> gpu
 R2 = RawData(S2, cc_len, cc_step) |> gpu
 R = [R1,R2]
 
-# preprocess on the GPU 
+# preprocess on the GPU
 detrend!.(R)
 taper!.(R)
 bandpass!.(R,freqmin,freqmax,zerophase=true)
 
-# FFT on GPU 
-FFT = compute_fft.(R)
+# Real FFT on GPU
+FFT = rfft.(R)
 whiten!.(FFT,freqmin,freqmax)
 
-# compute correlation and send to cpu 
-C = compute_cc(FFT[1],FFT[2],maxlag) |> cpu 
+# compute correlation and send to cpu
+C = correlate(FFT[1],FFT[2],maxlag) |> cpu
 ```
 
 ### Routines Implemented on the GPU
-- Preprocessing: 
+- Preprocessing:
   - `detrend`,`demean`, `taper`, `onebit`, `smooth`
-- Filtering: 
-  - `bandpass`, `bandstop`, `lowpass`, `highpass` 
+- Filtering:
+  - `bandpass`, `bandstop`, `lowpass`, `highpass`
 - Fourier Domain:
-  - `whiten`, `rfft`, `irfft` 
+  - `whiten`, `rfft`, `irfft`
 - Cross-correlation:
-  - `correlate`, `cross-coherence`, `deconvolution` 
+  - `correlate`, `cross-coherence`, `deconvolution`
 - Post-processing:
   - `stack`, `filter`s, etc..
-
-

src/ArrayFuncs.jl

@@ -1,91 +1,40 @@
 import SeisIO: demean, demean!, taper, taper!,detrend, detrend!
-import DSP: hilbert
-export detrend, detrend!, taper, taper!, demean, demean!, phase, phase!, hanningwindow
-export hilbert
+export detrend, detrend!, taper, taper!, demean, demean!, hanningwindow
 # Signal processing functions for arrays (rather than SeisData or SeisChannel)
 
 """
-    detrend!(X::AbstractArray{<:AbstractFloat,1})
+    detrend!(X::AbstractArray{<:AbstractFloat})
 
 Remove linear trend from array `X` using least-squares regression.
 """
-function detrend!(X::AbstractArray{<:AbstractFloat,1})
+function detrend!(X::AbstractArray{<:AbstractFloat})
     T = eltype(X)
-    N = length(X)
-    A = similar(X,T,N,2)
-    A[:,2] .= T(1)
-    A[:,1] .= range(T(1/N),T(1),length=N)
-    factor = prefactor(A)
-    coeff = factor * X
-    X[:] .-= A *coeff
-    return nothing
-end
-detrend(A::AbstractArray{<:AbstractFloat,1}) = (U = deepcopy(A);
-        detrend!(U);return U)
-
-"""
-    detrend!(X::AbstractArray{<:AbstractFloat,2})
-
-Remove linear trend from columns of `X` using least-squares regression.
-"""
-function detrend!(X::AbstractArray{<:AbstractFloat,2})
-    T = eltype(X)
-    Nrows, Ncols = size(X)
+    N = size(X,1)
 
     # create linear trend matrix
-    A = similar(X,T,Nrows,2)
+    A = similar(X,T,N,2)
     A[:,2] .= T(1)
-    A[:,1] .= range(T(1/Nrows),T(1),length=Nrows)
-    factor = prefactor(A)
+    A[:,1] .= range(T(0),T(1),length=N)
+    # create linear trend matrix
+    R = transpose(A) * A
+
+    # do the matrix inverse for 2x2 matrix
+    # this is really slow on GPU
+    Rinv = inv(Array(R)) |> typeof(R)
+    factor = Rinv * transpose(A)
 
-    # solve least-squares
-    for ii = 1:Ncols
-        coeff = factor * X[:,ii]
-        X[:,ii] .-=  A *coeff
-    end
+    # remove trend
+    X .-= A * (factor * X)
     return nothing
 end
-detrend(A::AbstractArray{<:AbstractFloat,2}) = (U = deepcopy(A);
+detrend(A::AbstractArray{<:AbstractFloat}) = (U = deepcopy(A);
         detrend!(U);return U)
 detrend!(R::RawData) = detrend!(R.x)
 detrend(R::RawData) = (U = deepcopy(R); detrend!(U.x); return U)
 detrend!(C::CorrData) = detrend!(C.corr)
 detrend(C::CorrData) = (U = deepcopy(C); detrend!(U.corr); return U)
 
 
-"""
-  prefactor(A)
-
-Computes prefactor (A^T * A) ^ -1 * X^T when solving Y = AX using linear-least squares.
-"""
-function prefactor(A::AbstractArray{<:AbstractFloat})
-    T = eltype(A)
-    N = size(A,1)
-
-    # create linear trend matrix
-    R = transpose(A) * A
-
-    # do the matrix inverse for 2x2 matrix
-    Rinv = inv2x2(R)
-
-    # return
-    return Rinv * transpose(A)
-end
-
-
-"""
-    inv2x2(X::AbstractArray{<:AbstractFloat,2})
-
-Inverse of 2x2 matrix `X`.
-"""
-function inv2x2(X::AbstractArray{<:AbstractFloat})
-    Xinv = -X
-    Xinv[1,1] = X[2,2]
-    Xinv[2,2] = X[1,1]
-    Xinv ./=  (X[1,1] .* X[2,2] .- X[1,2] .* X[2,1])
-    return Xinv
-end
-
 """
     demean!(A::AbstractArray{<:AbstractFloat})
 
@@ -159,38 +108,3 @@ function hanningwindow(A::AbstractArray, n::Int)
    win .= sin.(win).^2
    return win
 end
-
-"""
-    phase!(A::AbstractArray)
-
-Extract instantaneous phase from signal A.
-
-For time series `A`, its analytic representation ``S = A + H(A)``, where
-``H(A)`` is the Hilbert transform of `A`. The instantaneous phase ``e^{iθ}``
-of `A` is given by dividing ``S`` by its modulus: ``e^{iθ} = \\frac{S}{|S|}``
-For more information on Phase Cross-Correlation, see:
-[Ventosa et al., 2019](https://pubs.geoscienceworld.org/ssa/srl/article-standard/570273/towards-the-processing-of-large-data-volumes-with).
-"""
-function phase!(A::AbstractArray)
-    A .= angle.(hilbert(A))
-    return nothing
-end
-phase(A::AbstractArray) = (U = deepcopy(A);phase!(U);return U)
-phase!(R::RawData) = phase!(R.x)
-phase(R::RawData) = (U = deepcopy(R);phase!(U.x);return U)
-
-"""
-hilbert(A)
-
-Computes the analytic representation of x, x_a = x + j hilbert{x}.
-
-Only works for arrays on the GPU!
-"""
-function hilbert(A::AbstractGPUArray{Float32})
-    Nrows = size(A,1)
-    T = eltype(A)
-    f = fft(A,1)
-    f[2:Nrows÷2 + isodd(Nrows),:] .*= T(2.)
-    f[Nrows÷2+1 + isodd(Nrows):end,:] .= complex(T(0.))
-    return ifft(f,1)
-end

src/SeisNoise.jl

@@ -2,7 +2,7 @@ module SeisNoise
 
 using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO
 using Statistics, StatsBase, Plots, Distributed
-using Distributed, CuArrays, Adapt, CUDAnative, GPUArrays
+using CuArrays, Adapt, CUDAnative, GPUArrays
 
 # check use of cuda
 const use_cuda = Ref(false)
@@ -12,12 +12,8 @@ if !CuArrays.functional()
 end
 
 # import types first
-# include("Types/RawData.jl")
-# include("Types/FFTData.jl")
-# include("Types/CorrData.jl")
 include("Types/NoiseData.jl")
 include("Types/show.jl")
-include("Types/InputParams.jl")
 
 # import pre and post processing tools
 include("distance.jl")
@@ -26,8 +22,6 @@ include("ArrayFuncs.jl")
 include("tools.jl")
 include("slicing.jl")
 include("filter.jl")
-include("downsample.jl")
-include("availability.jl")
 include("phase_shift.jl")
 include("stacking.jl")
 
@@ -37,4 +31,7 @@ include("correlation.jl")
 include("rotation.jl")
 include("Plotting/plotting.jl")
 
+# out with the old
+include("deprecated.jl")
+
 end # module