Added support and tests for GPU acceleration

.github/workflows/CI.yml

@@ -4,6 +4,7 @@ on:
   workflow_dispatch: #get a button to trigger it by hand
   push:
     paths: # run CI/CD when the tests change, the dependencies change, or the version changes
+      - src/*
       - test/*
       - Project.toml
     tags: "*"

Project.toml

@@ -12,17 +12,26 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 
+[weakdeps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+
+[extensions]
+CUDAExt = "CUDA"
+
 [compat]
 AbstractFFTs = "1.0"
+CUDA = "4, 5, 6"
 Documenter = "0.27, 1"
 FFTW = "1.4"
 Interpolations = "0.13, 0.14, 0.15, 0.16"
 LinearAlgebra = "1"
 SpecialFunctions = "1.3, 2"
 Wavelets = "0.9, 0.10"
-julia = "1"
+julia = "1.9"
 
 [extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
@@ -33,4 +42,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 
 [targets]
-test = ["Test", "FFTW", "Wavelets", "Interpolations", "LinearAlgebra", "Logging", "Random", "Documenter"]
+test = ["Test", "FFTW", "Wavelets", "Interpolations", "LinearAlgebra", 
+        "Logging", "Random", "Documenter", "CUDA", "BenchmarkTools"]

docs/make.jl

@@ -5,9 +5,10 @@ ENV["LINES"] = "9"
 ENV["COLUMNS"] = "60"
 makedocs(sitename = "ContinuousWavelets.jl",
     format = Documenter.HTML(),
+    modules = [ContinuousWavelets],
     authors = "David Weber",
     clean = true,
-    strict = true,
+    checkdocs = :none,
     pages = [
         "basic usage" => "index.md",
         "Install" => "installation.md",
@@ -19,5 +20,4 @@ makedocs(sitename = "ContinuousWavelets.jl",
             "Inversion" => "inverse.md"
         ]
     ])
-
-deploydocs(repo = "github.com/UCD4IDS/ContinuousWavelets.jl.git")
+deploydocs(repo = "github.com/UCD4IDS/ContinuousWavelets.jl.git")

ext/CUDAExt.jl

@@ -0,0 +1,20 @@
+module CUDAExt
+
+using ContinuousWavelets: ZPBoundary, SymBoundary
+using ContinuousWavelets, CUDA
+
+function ContinuousWavelets.reflect(Y::CuArray, bt)
+    n1 = size(Y, 1)
+    if typeof(bt) <: ZPBoundary
+        base2 = ceil(Int, log2(n1))
+        padding = fill!(similar(Y, 2^(base2) - n1, size(Y)[2:end]...), 0)
+        x = cat(Y, padding, dims = 1)
+    elseif typeof(bt) <: SymBoundary
+        x = cat(Y, reverse(Y, dims = 1), dims = 1)
+    else
+        x = Y
+    end
+    return x
+end
+
+end

src/CWTConstruction.jl

@@ -1,3 +1,7 @@
+"""
+    CWT(wave::ContWaveClass, Q=8, boundary::WaveletBoundary=SymBoundary(),
+    averagingType::Average = Father(), averagingLength::Int = 4, frameBound=1, p::N=Inf, β=4)
+"""
 # the parameters are:
 #   B Boundary condition
 #   S Storage data type
@@ -29,6 +33,7 @@ struct CWT{B,S,W<:ContWaveClass,N,isAn} <: ContWave{B,S}
     # wavelets as the scale changes. The conjugate
     # p-norm is preserved for the  signal. Should be
     # larger than 1
+    fsample::S
 end
 
 # aliased = ((:Q,:s,:scalingFactor), (:β,:decreasing), (:p, :normalization))
@@ -54,19 +59,17 @@ function processKeywordArgs(Q, β, p; kwargs...)
 end
 
 
-@doc """
-    CWT(wave::ContWaveClass, Q=8, boundary::WaveletBoundary=SymBoundary(),
-    averagingType::Average = Father(), averagingLength::Int = 4, frameBound=1, p::N=Inf, β=4)
-"""
+defaultAveragingLength(::ContWaveClass) = 4
 function CWT(wave::WC,
     Q = 8,
     boundary::B = DEFAULT_BOUNDARY,
     averagingType::A = Father(),
-    averagingLength::Real = 0,
+    averagingLength::Real = defaultAveragingLength(wave),
     frameBound = 1,
     p::N = Inf,
     β = 4;
     extraOctaves = 0,
+    fsample=2000,
     kwargs...) where {WC<:ContWaveClass,A<:Average,B<:WaveletBoundary,N<:Real}
     Q, β, p = processKeywordArgs(Q, β, p; kwargs...) # some names are redundant
     @assert β > 0
@@ -93,7 +96,8 @@ function CWT(wave::WC,
         S(averagingLength),
         averagingType,
         S(frameBound),
-        S(p))
+        S(p),
+        S(fsample))
 end
 
 
@@ -151,10 +155,9 @@ function waveletType(::CWT{B,T,W,N}) where {B,T,W,N}
 end
 
 function Base.show(io::IO, cf::CWT{W,S,WT,N}) where {W,S,WT,N}
-    print(io,
-        "CWT{$(cf.waveType), $(cf.averagingType), Q=$(cf.Q), β=$(cf.β)," *
-        "aveLen=$(cf.averagingLength), frame=" *
-        "$(cf.frameBound), norm=$(cf.p), extraOctaves=$(cf.extraOctaves)}")
+    print(io, "CWT{$(cf.waveType), $(cf.averagingType), Q=$(cf.Q), β=$(cf.β)," *
+        "aveLen=$(cf.averagingLength), frame=$(cf.frameBound), norm=$(cf.p), " *
+        "extraOctaves=$(cf.extraOctaves), fsample=$(cf.fsample)}")
 end
 
 
@@ -172,6 +175,7 @@ function wavelet(wave::WC;
     frameBound = 1,
     p::N = Inf,
     β = 4,
+    fsample=2000,
     kwargs...) where {WC<:ContWaveClass,A<:Average,T<:WaveletBoundary,N<:Real}
     return CWT(wave,
         Q,
@@ -181,5 +185,6 @@ function wavelet(wave::WC;
         frameBound,
         p,
         β;
+        fsample=fsample, 
         kwargs...)
 end

src/ContinuousWavelets.jl

@@ -44,6 +44,8 @@ export qmf,
     getDualCoverage,
     caveats
 
+export reflect
+
 @doc """
     ContWave{Boundary,T}
 The abstract type encompassing the various types of wavelets implemented in

src/apply.jl

@@ -46,7 +46,8 @@ function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N
         OutType = T
     end
 
-    wave = zeros(OutType, size(x)..., nScales)  # result array
+    # wave = zeros(OutType, size(x)..., nScales)  # result array (CPU Safe)
+    wave = fill!(similar(x, OutType, size(x)..., nScales), 0) # (GPU Safe)
     # faster if we put the example index on the outside loop through all scales
     # and compute transform
     if isAnalytic(cWav.waveType)
@@ -131,46 +132,49 @@ end
 function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
     outer = axes(x̂)[2:end]
     n1 = size(x̂, 1)
-    isSourceEven = mod(size(wave, 1) + 1, 2)
+    nWave = size(wave, 1)
+    isSourceEven = mod(nWave + 1, 2)
+    negFreqEnd = n1 - isSourceEven
     # the averaging function isn't analytic, so we need to do both positive and
     # negative frequencies
     @views tmpWave = x̂ .* daughters[:, 1]
-    @views wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceEven,
-            outer...]),
-        dims = 1)
-    @views wave[1:n1, outer..., 1] = tmpWave
-    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
+    @views wave[(n1+1):end, outer..., 1] .= reverse(conj.(tmpWave[2:negFreqEnd,
+            outer...]), dims = 1)
+    @views wave[1:n1, outer..., 1] .= tmpWave
+    @views wave[:, outer..., 1] .= fftPlan \ copy(wave[:, outer..., 1])  # averaging
     for j = 2:size(daughters, 2)
-        @views wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
-        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
+        @views wave[1:n1, outer..., j] .= x̂ .* daughters[:, j]
+        @views wave[:, outer..., j] .= fftPlan \ copy(wave[:, outer..., j]) # wavelet transform
     end
 end
 
 # analytic on complex data with an averaging function
 function analyticTransformComplex!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
     outer = axes(x̂)[2:end]
     n1 = size(daughters, 1)
-    isSourceEven = mod(size(wave, 1) + 1, 2)
+    nWave = size(wave, 1)
+    isSourceEven = mod(nWave + 1, 2)
+    negFreqStart = n1 - isSourceEven + 1
     # the averaging function isn't analytic, so we need to do both positive and
     # negative frequencies
     @views positiveFreqs = x̂[1:n1, outer...] .* daughters[:, 1]
-    @views negativeFreqs = x̂[(n1-isSourceEven+1):end, outer...] .*
+    @views negativeFreqs = x̂[negFreqStart:end, outer...] .*
                            reverse(conj.(daughters[2:end, 1]))
-    @views wave[(n1-isSourceEven+1):end, outer..., 1] = negativeFreqs
-    @views wave[1:n1, outer..., 1] = positiveFreqs
-    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
+    @views wave[negFreqStart:end, outer..., 1] .= negativeFreqs
+    @views wave[1:n1, outer..., 1] .= positiveFreqs
+    @views wave[:, outer..., 1] .= fftPlan \ copy(wave[:, outer..., 1])  # averaging
     for j = 2:size(daughters, 2)
-        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
-        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
+        @views wave[1:n1, outer..., j] .= x̂[1:n1, outer...] .* daughters[:, j]
+        @views wave[:, outer..., j] .= fftPlan \ copy(wave[:, outer..., j])  # wavelet transform
     end
 end
 
 function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType)
     outer = axes(x̂)[2:end]
     n1 = size(x̂, 1)
     for j = 1:size(daughters, 2)
-        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
-        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
+        @views wave[1:n1, outer..., j] .= x̂[1:n1, outer...] .* daughters[:, j]
+        @views wave[:, outer..., j] .= fftPlan \ copy(wave[:, outer..., j])  # wavelet transform
     end
 end
 
@@ -180,8 +184,8 @@ function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::NoAve)
     n1 = size(x̂, 1)
     # the no averaging version
     for j = 1:size(daughters, 2)
-        wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
-        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
+        @views wave[1:n1, outer..., j] .= x̂ .* daughters[:, j]
+        @views wave[:, outer..., j] .= fftPlan \ copy(wave[:, outer..., j])  # wavelet transform
     end
 end
 
@@ -196,7 +200,7 @@ function otherwiseTransform!(wave::AbstractArray{<:Real},
     n1 = size(x̂, 1)
     for j = 1:size(daughters, 2)
         @views tmp = x̂ .* daughters[:, j]
-        @views wave[:, outer..., j] = fromPlan \ tmp  # wavelet transform
+        @views wave[:, outer..., j] .= fromPlan \ tmp  # wavelet transform
     end
 end
 
@@ -210,20 +214,21 @@ function otherwiseTransform!(wave::AbstractArray{<:Complex},
     outer = axes(x̂)[2:end]
     n1 = size(daughters, 1)
     isSourceEven = mod(size(fromPlan, 1) + 1, 2)
+    negStart = n1 - isSourceEven + 1
     for j = 1:size(daughters, 2)
-        @views wave[1:n1, outer..., j] = @views x̂[1:n1, outer...] .* daughters[:, j]
-        @views wave[n1-isSourceEven+1:end, outer..., j] = x̂[n1-isSourceEven+1:end,
-            outer...] .* reverse(conj.(daughters[2:end,
-            j]))
-        @views wave[:, outer..., j] = fromPlan \ (wave[:, outer..., j])  # wavelet transform
+        @views wave[1:n1, outer..., j] .= x̂[1:n1, outer...] .* daughters[:, j]
+        @views wave[negStart:end, outer..., j] .= x̂[negStart:end, outer...] .*
+                                            reverse(conj.(daughters[2:end, j]))
+        @views wave[:, outer..., j] .= fromPlan \ copy(wave[:, outer..., j])  # wavelet transform
     end
 end
 
-function reflect(Y, bt)
+function reflect(Y, bt) 
     n1 = size(Y, 1)
     if typeof(bt) <: ZPBoundary
         base2 = ceil(Int, log2(n1))   # power of 2 nearest to N
-        x = cat(Y, zeros(2^(base2) - n1, size(Y)[2:end]...), dims = 1)
+        # x = cat(Y, zeros(2^(base2) - n1, size(Y)[2:end]...), dims = 1)
+        x = cat(Y, zeros(eltype(Y), 2^(base2) - n1, size(Y)[2:end]...), dims = 1)
     elseif typeof(bt) <: SymBoundary
         x = cat(Y, reverse(Y, dims = 1), dims = 1)
     else
@@ -232,7 +237,6 @@ function reflect(Y, bt)
     return x
 end
 
-
 function cwt(Y::AbstractArray{T},
     c::CWT{W};
     varArgs...) where {T<:Number,W<:WaveletBoundary}