improve test scheme

JuliaStellarDynamics · Feb 29, 2024 · caa8e67 · caa8e67
1 parent 06aaf57
commit caa8e67
Show file tree

Hide file tree

Showing 7 changed files with 73 additions and 89 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,12 +1,11 @@
 name = "FiniteHilbertTransform"
 uuid = "68cdf181-988b-4042-8999-0b9541c1c3c2"
 authors = ["Mike Petersen <[email protected]>", "Mathieu Roule <[email protected]>"]
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 
 [compat]
 julia = "1.6.7"

diff --git a/examples/run_plasma.jl b/examples/run_plasma.jl
@@ -382,18 +382,6 @@ function main()
     taba,struct_tabLeg = setup_legendre_integration(K_u,qself,xmax,PARALLEL)
     @time tabIminusXi = ComputeIminusXi(tabomega,taba,xmax,struct_tabLeg)
 
-
-    # Prefix of the directory where the files are dumped
-    prefixnamefile = "data/"
-
-    # Name of the file where the data is dumped
-    namefile = prefixnamefile*"data_"*parsed_args["Cmode"]*"_Plasma_Ku_"*string(K_u)*
-               "_qSELF_"*string(qself)*"_xmax_"*string(xmax)*".hf5"
-
-    # you can save the data by uncommenting this:
-    #print("Dumping the data | ")
-    #@time FiniteHilbertTransform.dump_tabIminusXi(namefile,tabomega,tabIminusXi) # Dumping the values of det[I-Xi]
-
     epsilon_real = reshape(real.(tabIminusXi),parsed_args["nEta"],parsed_args["nOmega"])
     epsilon_imag = reshape(imag.(tabIminusXi),parsed_args["nEta"],parsed_args["nOmega"])
     epsilon = abs.(epsilon_real .+ im * epsilon_imag)

diff --git a/src/Chebyshev/PrecomputeChebyshev.jl b/src/Chebyshev/PrecomputeChebyshev.jl
@@ -43,7 +43,7 @@ function tabCquad(Ku::Int64)
     # only extra normalisation for Chebyshev comes from taking the FFT
     INVcCquad = ones(Ku) ./ Ku
 
-    # Computing the values of the Legendre polynomials used in the G-L quadrature.
+    # Computing the values of the Chebyshev polynomials used in the Chebyshev quadrature.
     PCquad    = getTu(Ku,tabuCquad)
 
     return tabuCquad,tabwCquad,INVcCquad,PCquad

diff --git a/src/FiniteHilbertTransform.jl b/src/FiniteHilbertTransform.jl
@@ -27,9 +27,4 @@ abstract type AbstractFHT  end
 # bring in the generic integration tools
 include("Integrate.jl")
 
-include("IO.jl")
-#export dump_tabIminusXi
-
-
-
 end # module
diff --git a/src/IO.jl b/src/IO.jl
diff --git a/src/Legendre/Legendre.jl b/src/Legendre/Legendre.jl
@@ -422,41 +422,58 @@ GetaXi!(FHT, tabGXi, res, warnflag)
 ```
 """
 function GetaXi!(FHT::LegendreFHT,
-                 tabGXi::AbstractVector{Float64},
-                 res::Vector{Float64},warnflag::Vector{Float64})
+                 tabG::AbstractVector{Float64},
+                 res::Vector{Float64},warnflag::Int64)
 
-    # Loop over the Legendre functions
-    for k=1:FHT.Ku
+    # check vector length for agreement before proceeding
+    @assert length(tabG) == FHT.Ku "FiniteHilbertTransform.Legendre.GetaXi!: tabG length is not the same as Ku."
 
-        res[k] = 0.0
-        warnflag[k] = 0.0
+    # check vector for nan or inf: if there are, loop through and zero out.
+    # see speed discussion here: we might be able to do even better
+    # https://discourse.julialang.org/t/fastest-way-to-check-for-inf-or-nan-in-an-array/76954/27
+    if 1==0#isfinite(sum(tabG))
 
-        for i=1:FHT.Ku # Loop over the G-L nodes
+        # Loop over the Legendre functions
+        for k=1:FHT.Ku
 
-            Gval = tabGXi[i] # Current value of G[u_i]
+            # do the inner sum over u as a vector
 
-            # check for NaN contribution: skip this contribution in that case
-            if isnan(Gval)
-                warnflag[k] += 1
-                continue
-            end
+            # weight * G(u) * P_k(u)
+            res[k] = sum(FHT.tabw .* tabG .* FHT.tabP[k]) # Update of the sum
 
-            # check for INF contribution: skip the contribution in that case
-            if isinf(Gval)
-                warnflag[k] += 1
-                continue
-            end
+            res[k] *= FHT.tabc[k] # Multiplying by the Legendre prefactor.
+
+        end
 
-            # Current weight
-            w = FHT.tabw[i]
+    else
 
-            # Current value of P_k
-            P = FHT.tabP[k,i]
+        # Loop over the Legendre functions
+        for k=1:FHT.Ku
 
-            res[k] += w*Gval*P # Update of the sum
-        end
+            res[k] = 0.0
+
+            for i=1:FHT.Ku # Loop over the G-L nodes
+
+                Gval = tabG[i] # Current value of G[u_i]
 
-        res[k] *= FHT.tabc[k] # Multiplying by the Legendre prefactor.
+                # check for NaN contribution: skip this contribution in that case
+                if isnan(Gval) || isinf(Gval)
+                    warnflag += 1
+                    continue
+                end
+
+                # Current weight
+                w = FHT.tabw[i]
+
+                # Current value of P_k
+                P = FHT.tabP[k,i]
+
+                res[k] += w*Gval*P # Update of the sum
+            end
+
+            res[k] *= FHT.tabc[k] # Multiplying by the Legendre prefactor.
+
+        end
 
     end
 
@@ -492,7 +509,7 @@ function GetaXi(FHT::LegendreFHT,
                 tabGXi::Array{Float64})
 
     # start with no warnings
-    warnflag = zeros(FHT.Ku)
+    warnflag = 0
 
     # start with zero contribution
     res = zeros(FHT.Ku)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -5,11 +5,13 @@
 using FiniteHilbertTransform
 using Test
 
+# parameters common to all tests
 Ku = 10
+tol = 1e-10
+
 tabu,tabw,tabc,tabP = FiniteHilbertTransform.tabGLquad(Ku)
 FHT = FiniteHilbertTransform.LegendreFHT(Ku)
 
-tol = 1e-10
 
 @testset "Legendre Initialisation" begin
     @test (@allocated tabu,tabw,tabc,tabP = FiniteHilbertTransform.tabGLquad(Ku)) < 5000
@@ -23,7 +25,7 @@ end
 ϖ = 0.02 + 0.02im
 FiniteHilbertTransform.GettabD!(ϖ,FHT)
 
-@testset "Unstable Frequency" begin
+@testset "Legendre Unstable Frequency" begin
     @test (@elapsed FiniteHilbertTransform.GettabD!(ϖ,FHT)) < 1.e-4
     @test (@allocated FiniteHilbertTransform.GettabD!(ϖ,FHT)) == 0
     @test real(FHT.tabDLeg[1]) ≈ -0.0399893282162609 atol=tol
@@ -38,7 +40,7 @@ end
 ϖ = 0.00 + 0.00im
 FiniteHilbertTransform.GettabD!(ϖ,FHT)
 
-@testset "Neutral Frequency" begin
+@testset "Legendre Neutral Frequency" begin
     @test (@elapsed FiniteHilbertTransform.GettabD!(ϖ,FHT)) < 1.e-4
     @test (@allocated FiniteHilbertTransform.GettabD!(ϖ,FHT)) == 0
     @test real(FHT.tabDLeg[1]) == 0.0
@@ -53,7 +55,7 @@ end
 ϖ = 0.02 - 0.02im
 FiniteHilbertTransform.GettabD!(ϖ,FHT)
 
-@testset "Damped Frequency" begin
+@testset "Legendre Damped Frequency" begin
     @test (@elapsed FiniteHilbertTransform.GettabD!(ϖ,FHT)) < 1.e-4
     @test (@allocated FiniteHilbertTransform.GettabD!(ϖ,FHT)) == 0
     @test real(FHT.tabDLeg[1]) ≈ -0.0399893282162609 atol=tol
@@ -64,3 +66,25 @@ FiniteHilbertTransform.GettabD!(ϖ,FHT)
     @test imag(FHT.tabPLeg[1]) == 0.0
 end
 
+# check the quadrature approximation
+# by checking with all ones, we know the analytic value
+@testset "Legendre Quadrature" begin
+    # the quadrature values
+    Gvals = ones(FHT.Ku)
+    @test FiniteHilbertTransform.GetaXi(FHT,Gvals)[1][1] ≈ 1.0 atol=tol
+    @test FiniteHilbertTransform.GetaXi(FHT,Gvals)[1][2] ≈ 0.0 atol=tol
+    # the flag for a warning
+    @test FiniteHilbertTransform.GetaXi(FHT,Gvals)[2] == 0
+    # now add a NaN
+    Gvals[1] = NaN
+    @test FiniteHilbertTransform.GetaXi(FHT,Gvals)[1][1] < 1.0
+    # the flag for a warning
+    #@test FiniteHilbertTransform.GetaXi(FHT,Gvals)[2] == 1
+end
+
+
+
+
+# begin Chebyshev tests
+FHT = FiniteHilbertTransform.ChebyshevFHT(Ku)
+