Merge branch 'master' into feat-libfasttransforms

Author: MikaelSlevinsky

.travis.yml

@@ -14,6 +14,7 @@ julia:
   - 1.0
   - 1.1
   - 1.2
+  - 1.3
   - nightly
 matrix:
   allow_failures:

Project.toml

@@ -1,11 +1,12 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -14,9 +15,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
-AbstractFFTs = "≥ 0.3.1"
-DSP = "≥ 0.4.0"
-FFTW = "≥ 0.0.4"
-SpecialFunctions = "≥ 0.3.4"
-ToeplitzMatrices = "≥ 0.4.0"
+AbstractFFTs = "≥ 0.4"
+DSP = "≥ 0.6"
+FastGaussQuadrature = "≥ 0.4"
+FFTW = "≥ 0.3"
+SpecialFunctions = "≥ 0.8"
+ToeplitzMatrices = "≥ 0.6"
 julia = "≥ 1.0"

appveyor.yml

@@ -1,9 +1,9 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1.0
   - julia_version: 1.1
-  - julia_version: 1.2    
+  - julia_version: 1.2
+  - julia_version: 1.3
   - julia_version: nightly
 
 platform:

src/FastTransforms.jl

@@ -1,7 +1,6 @@
-__precompile__()
 module FastTransforms
 
-using AbstractFFTs, DSP, FFTW, Libdl, LinearAlgebra, SpecialFunctions, ToeplitzMatrices
+using AbstractFFTs, DSP, FastGaussQuadrature, FFTW, Libdl, LinearAlgebra, SpecialFunctions, ToeplitzMatrices
 
 import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
              inv, size, view
@@ -13,6 +12,8 @@ import AbstractFFTs: Plan
 
 import DSP: conv
 
+import FastGaussQuadrature: unweightedgausshermite
+
 import FFTW: dct, dct!, idct, idct!,
              plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
              plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct,
@@ -61,6 +62,10 @@ export plan_fejer1, fejernodes1, fejerweights1,
 include("clenshawcurtis.jl")
 include("fejer.jl")
 
+export hermitepoints, weightedhermitetransform, iweightedhermitetransform
+
+include("hermite.jl")
+
 include("fftBigFloat.jl")
 
 export gaunt
@@ -72,6 +77,8 @@ export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
        diskones, diskzeros, diskrand, diskrandn,
        triones, trizeros, trirand, trirandn, trievaluate
 
+lgamma(x) = logabsgamma(x)[1]
+
 include("specialfunctions.jl")
 
 end # module

src/hermite.jl

@@ -0,0 +1,38 @@
+
+
+hermitepoints(n) = unweightedgausshermite( n )[1]
+
+
+# exp(-x^2/2) H_n(x) / sqrt(π*prod(1:n))
+struct ForwardWeightedHermitePlan <: FastTransformPlan{FORWARD,Float64}
+    Vtw::Matrix{Float64} # vandermonde
+end
+
+struct BackwardWeightedHermitePlan <: FastTransformPlan{BACKWARD,Float64}
+    V::Matrix{Float64} # vandermonde
+end
+
+
+function _weightedhermite_vandermonde(n)
+    V = Array{Float64}(undef,n,n)
+    x,w = unweightedgausshermite( n )
+    for k=1:n
+        V[k,:] = FastGaussQuadrature.hermpoly_rec(0:n-1, sqrt(2)*x[k])
+    end
+    V,w
+end
+function ForwardWeightedHermitePlan(n::Integer) 
+    V,w = _weightedhermite_vandermonde(n)
+    ForwardWeightedHermitePlan(V' * Diagonal(w / sqrt(π)))
+end
+BackwardWeightedHermitePlan(n::Integer) = BackwardWeightedHermitePlan(_weightedhermite_vandermonde(n)[1])
+
+
+*(P::ForwardWeightedHermitePlan, v::AbstractVector) = P.Vtw*v
+*(P::ForwardWeightedHermitePlan, v::AbstractVector{Float64}) = P.Vtw*v
+*(P::BackwardWeightedHermitePlan, v::AbstractVector) = P.V*v
+*(P::BackwardWeightedHermitePlan, v::AbstractVector{Float64}) = P.V*v
+
+weightedhermitetransform(v) = ForwardWeightedHermitePlan(length(v))*v
+iweightedhermitetransform(v) = BackwardWeightedHermitePlan(length(v))*v
+

test/hermitetests.jl

@@ -0,0 +1,24 @@
+using FastTransforms
+
+@testset "Hermite" begin
+    @test hermitepoints(1) == [0.0]
+    @test hermitepoints(100_000)[end] ≈ 446.9720305443094
+
+    @test weightedhermitetransform([1.0]) == [1.0]
+    @test weightedhermitetransform(exp.(-hermitepoints(2).^2/2)) == [1.0,0.0]
+    @test weightedhermitetransform(exp.(-hermitepoints(3).^2/2)) ≈ [1.0,0.0,0.0]
+    @test weightedhermitetransform(exp.(-hermitepoints(1000).^2/2)) ≈ [1.0; zeros(999)]
+    @test weightedhermitetransform(exp.(-hermitepoints(3000).^2/2)) ≈ [1.0; zeros(2999)]
+
+    for n in (1, 5,100)
+        x = randn(n)
+        @test iweightedhermitetransform(weightedhermitetransform(x)) ≈ x
+        @test weightedhermitetransform(iweightedhermitetransform(x)) ≈ x
+    end
+
+    x = hermitepoints(100)
+    @test iweightedhermitetransform([0.0; 1.0; zeros(98)]) ≈ (exp.(-x.^2 ./ 2) .* 2x/sqrt(2))
+    @test iweightedhermitetransform([0.0; 0; 1.0; zeros(97)]) ≈ (exp.(-x.^2 ./ 2) .* (4x.^2 .- 2)/(sqrt(2)*2^(2/2)))
+    @test iweightedhermitetransform([0.0; 0; 0; 1.0; zeros(96)]) ≈ (exp.(-x.^2 ./ 2) .* (-12x + 8x.^3) / (sqrt(2*3)*2^(3/2)))
+    @test iweightedhermitetransform([0.0; 0; 0; 0; 1.0; zeros(95)]) ≈ (exp.(-x.^2 ./ 2) .* (12 .- 48x.^2 .+ 16x.^4) / (sqrt(2*3*4)*2^(4/2)))
+end

test/runtests.jl

@@ -16,9 +16,10 @@ end
 include("nuffttests.jl")
 
 include("fftBigFloattests.jl")
-
 include("paduatests.jl")
 
 include("gaunttests.jl")
 
+include("hermitetests.jl")
+
 include("toeplitztests.jl")